Portal de Periódicos da CAPES

Novel Subharmonic Resonance Periodic Orbits of a Solar Sail in Earth–Moon System

2019; American Institute of Aeronautics and Astronautics; Volume: 42; Issue: 11 Linguagem: Inglês

10.2514/1.g004377

1533-3884

Ying-Jing Qian, Zixiao Liu, Xiao-Dong Yang, Inseok Hwang, Wei Zhang,

Vibration Control and Rheological Fluids

Open AccessEngineering NotesNovel Subharmonic Resonance Periodic Orbits of a Solar Sail in Earth–Moon SystemYing-Jing Qian, Zi-Xiao Liu, Xiao-Dong Yang, Inseok Hwang and Wei ZhangYing-Jing QianBeijing University of Technology, 100124 Beijing, People's Republic of China, Zi-Xiao LiuBeijing University of Technology, 100124 Beijing, People's Republic of China, Xiao-Dong YangBeijing University of Technology, 100124 Beijing, People's Republic of China, Inseok HwangPurdue University, West Lafayette, Indiana 47907-2023 and Wei ZhangBeijing University of Technology, 100124 Beijing, People's Republic of ChinaPublished Online:12 Aug 2019https://doi.org/10.2514/1.G004377SectionsRead Now ToolsAdd to favoritesDownload citationTrack citations ShareShare onFacebookTwitterLinked InRedditEmail AboutNomenclatureAnamplitudes for free-oscillation solutionsB, Camplitudes for forced-oscillation solutions, respectivelya0solar sail's characteristic acceleration in nondimensional unitscccomplex conjugationDout-of-plane amplitudef, gsolar radiation accelerations in plane and out of plane, respectivelyT0, T2fast time scale and slow time scaleαijk, βijk, γijkknown coefficients for polynomials x2, xy, y2,…, in the Earth–moon system when the partial derivatives of the pseudopotential function Ω are expanded in polynomial formsαn, βnreal and imaginary parts of amplitude Anγpitch angle of the sail with respect to the Earth–moon orbital planeεbookmark symbol used to indicate the order of a small parameter, which equals one in the simulationsμmass parameter of the systemσdetuning parameterΩpseudopotential function of the Earth–moon–solar-sail three-body problemΩij(i, j) partial derivatives of the three-body potential evaluated at the Lagrange pointωsfrequency of the sun line in the synodic coordinate system that equals 0.9252ωnsystem natural frequenciesI. IntroductionAs a special type of spacecraft, a solar-sail spacecraft uses radiation pressure on a reflective sail for propulsion without any fuel consumption [1]. In recent years, solar-sail technology has been successfully applied to various missions, such as the Interplanetary Kite-craft Accelerated by Radiation Of the Sun mission by the Japan Aerospace Exploration Agency [2], the LightSail-1 mission by the Planetary Society [3], and NanoSail-D by NASA [4]. Due to advantageous applications in space mission design [5], solar-sail orbits around the libration points become more attractive: especially in novel orbit construction [6–13].From a mathematical point of view, for the Earth–moon–solar-sail system, the equations of motion become a system of inhomogeneous nonlinear ordinary differential equations due to the periodic solar radiation pressure (i.e., SRP). This situation is contrary to the homogeneous nonlinear case for the sun-Earth–solar-sail system. The method of obtaining periodic solutions for the sun-Earth–solar-sail problem around artificial libration points is similar to the one used in the sun–Earth restricted three-body problem [9,14–19]. However, a first-order approximation of the general solution for a nonlinear inhomogeneous system can be considered to consist of two parts: the general solution of the corresponding homogeneous part of the system, and the particular solution. Generally, the particular solutions are determined by the SRP acceleration, which shows up as external excitation frequencies and amplitudes. So, they are also called forced-oscillation solutions in nonlinear dynamics [20]. On the other hand, the general solution for the homogeneous part is relevant to the system's natural frequencies, and Nayfeh and Mook [20] called them free-oscillation solutions. Note that the initial conditions determine which periodic motion represents the actual response.The most representative studies about the particular solutions (forced-oscillation solutions) for periodic motions around the Earth–moon libration points with a solar sail were demonstrated by Gong et al. [21], Simo and McInnes [22], and Tamakoshi and Kojima [23]. Gong et al. [21] introduced quasi-periodic and periodic orbits of the solar sail around the L2 libration point. In their study, there are two time-varying angles in the expression for the SRP acceleration components, namely, the angle of the moon rotating around the Earth θ and the angle of the sun in the Earth–moon rotating frame λ. Therefore, the planar particular solution in Gong et al.'s study consists of two frequencies corresponding to λ and θ. If the inclination of the moon's orbital plane is set to zero, the result will degenerate to the particular solutions discussed by Simo and McInnes [22] and McInnes [24], which have only one frequency corresponding to λ. Likewise, for the solar-sail problem in binary asteroid systems, Heiligers and Scheeres [25] investigated the particular solution for a periodic orbit for which the frequency equals the sun's synodic frequency.Recently, by means of numerical methods, Heiligers et al. [26,27] investigated solar-sail Lyapunov and halo orbits in the Earth–moon three-body problem. It is a good demonstration of applying the general solutions for the homogeneous part of the system (free-oscillation solutions) to construct the solar-sail periodic orbit. Taking Lyapunov orbits as examples, the in-plane linear natural frequencies ωL1 and ωL2 for the Earth–moon L1 and L2 are 2.3343 and 1.8626, respectively. The synodic lunar frequency ωs equals 0.9252. Because the nonlinear frequencies of periodic orbits will decrease with the increase of amplitudes, the possible resonant ratios of ωL1=2ωs and ωL2=2ωs can be obtained by tuning the amplitudes. Therefore, the Lyapunov orbit with a period of (1/2)(2π/ωs) was chosen as an initial condition to construct periodic orbits in the works of Heiligers et al. [26,27]. Furthermore, Heiligers et al. [27] summarized the corresponding orbital periods of all classical libration point periodic orbits for the solar-sail motions in the Earth–moon system. Their study gave an insight in the connection between the possible resonances ratio and the natural frequencies and synodic lunar period. Mathematically speaking, the particular solutions always exist, whereas the free-oscillation solutions can be artificially designed by choosing the initial conditions. Thus, trajectories in the works of Heiligers et al. [26,27] are restrained to follow only the free-oscillation solutions by differential corrections.In this Note, we focus on a special one-third subharmonic resonance period orbit in the vicinity of the triangular libration point in the Earth–moon–solar-sail system. We show that, by explicitly considering both free-oscillation and forced-oscillation terms related to the solar radiation pressure, new orbits corresponding to subharmonic resonance within a certain domain can be designed. The domain in which the novel orbits exist can be analytically found using the perturbation method. The proposed orbit has exactly three times the period of the traditional one in the Earth–moon system. The amplitude of the steady-state subharmonic resonance orbit is computed using a frequency-response equation. Lastly, the stability analysis and some discussions on the proposed orbits are presented.II. Solar Sail in Earth–Moon Restricted Three-Body ProblemIn the circular restricted three-body problem (CRTBP), two primaries revolve around their barycenter in circular orbits under their mutual gravitational attractions. Another small particle (i.e., solar-sail spacecraft) moves in the plane defined by the two revolving primaries. The particle is considered as "massless," which is attracted by the two primaries, but the motions of the primary bodies are assumed not to be affected by the particle.The geometry of the Earth–moon–solar-sail three-body problem is conveniently described in a synodic coordinate system (O–XYZ), which is centered at the barycenter of the two primaries O, as shown in Fig. 1. In this synodic coordinate system, m1 represents the larger primary (Earth) and m2 is the smaller primary (moon). We define μ=m2/(m1+m2)=0.01215 as the mass parameter of the system and consider the Earth–moon distance as the unit. The time unit t is defined as the period time of m2 around m1. The Earth, with a mass of 1−μ, is located at (−μ,0,0); and the moon, with a mass of μ, is located at (1−μ,0,0) [28].Fig. 1 Solar sail in the Earth–moon restricted three-body problem.Note that the Earth–moon CRTBP is autonomous and conserved. However, when we introduce the solar radiation pressure into the system, the system becomes nonautonomous with external periodic excitation. A schematic of the nonautonomous Earth–moon-sail three-body problem is shown in Fig. 2.Fig. 2 Schematic of the nonautonomous Earth–moon-sail three-body problem.The sun line direction S in the rotating coordinate system can be expressed as S=[cos(ωst)−sin(ωst)0]T(1)where ωs=0.9252 is the angular rate of the sun line in the synodic coordinate system. Note that Eq. (1) ignores the small inclination difference between the sun–Earth and the Earth–moon orbital planes, and the distance between the sun and the Earth–moon system maintains a constant. At a time of t=0, the sun is assumed to be on the negative x axis.Due to the periodic motion of the sun, the solar acceleration as appears as an inhomogeneous term which can be written as as=a0(S⋅n)2n(2)where a0 is the solar sail's characteristic acceleration exerted on the sail in the nondimensional unit. We consider a flat and perfectly reflecting sail, and the force due to the SRP is normal to its surface [27]. The unit normal to the sail n is given by n=[cosγcos(ωst)−cosγsin(ωst)sinγ]T(3)where γ is the pitch angle of the sail with respect to the Earth–moon orbital plane.With the preceding definition, the dynamical equation for the solar sail in O–xyz is described by r¨+2[−Y˙X˙0]T=∂Ω∂r+as(4)where Ω is the pseudopotential function of the three-body problem and is given by Ω=12(X2+Y2)+1−μR1+μR2+12μ(1−μ)(5)in which R1 and R2 are the distances of the particle from P1 and P2, respectively.To investigate the motion around the libration points, it is necessary to move the origin of the coordinate system from the barycenter of the system to the libration point. In the new coordinate system, the directions of the x, y, and z axes are in parallel with those of the barycenter synodic system (Fig. 1). Note that the L4– centered synodic coordinate system is denoted as L4–xyz. In the following discussions, we will take the L4 point as an example because, for the L5 case, the problem can be considered in a similar way.The pseudopotential function of the motion around the L4 point with Legendre polynomials [28] can be expressed as x¨−2y˙−Ωxxx−Ωxyy=∑n=2∞∑i+j+k=nαijkxiyjzk+F1cos(ωst),y¨+2x˙−Ωxyx−Ωyyy=∑n=2∞∑i+j+k=nβijkxiyjzk−F1sin(ωst),z¨+z=∑n=2∞∑i+j+k=nγijkxiyjzk+F2(6)and F1=a0cos3γ,F2=a0sinγcos2γ(7)where Ωxx=3/4, Ωyy=9/4, and Ωxy=(33/4)(1−2μ). The expressions for αijk, βijk, and γijk up to the third order can be found in the Appendix. The expressions for SRP acceleration, F1 and F2, are similar to equation 5.137 in the work of McInnes [1] but with different notations. In particular, the sail acceleration is assumed to be constant under the small displacement for the L4 point. Therefore, the SRP is independent of the position of the solar sail.By inspecting Eq. (6), the radiation pressure acts as a pure external excitation source. Therefore, the system has the same natural frequencies as the Earth–moon CRTBP system. The planar natural frequencies for periodic motion in the vicinity of L4 are ω1=0.2982 and ω2=0.9545, and the possible resonance is ωs≈3ω1. The periodic orbits in the vicinity of the triangular libration points in the Earth–moon–solar-sail system are still stable.III. One-Third Subharmonic Orbit AnalysisIn this section, a special one-third subharmonic resonance period orbit will be obtained by superpositioning the free-oscillation and forced-oscillation solutions.To obtain such an orbit, a perturbation solution for a finite but small amplitude is derived using the method of multiple scales [20]. According to this method, x, y, and z are assumed to be functions of two independent time scales; a fast time of T0=t describing the fast linear oscillation; and a slow time of T2=ε2t describing the nonlinear correction in a secular sense, where t is the dimensionless time for the system in Eq. (6). The ε, x, y, and z are assumed to possess the following uniformly valid expansions for all times: x=εx1+ε2x2+ε3x3,…,y=εy1+ε2y2+ε3y3,…,z=εz1+ε2z2+ε3z3,…,F1=εf,F2=εg(8)Then, the time derivative becomes d/dt=D0+ε2D2+…, Dn=∂/∂Tn.Substituting Eq. (8) into Eq. (6), expanding, and equating the coefficients of equal powers of ε on both sides lead to D02x1−2D0y1−Ωxxx1−Ωxyy1=fcos(ωst),D02y1+2D0x1−Ωxyx1−Ωyyy1=−fsin(ωst),D02z1+z1=g(9)The solution to Eq. (9) can be easily obtained as x1=A1(T2)eiω1T0+A2(T2)eiω2T0+BeiωsT0+cc,y1=Λ1A1(T2)eiω1T0+Λ2A2(T2)eiω2T0+CeiωsT0+cc,z1=DeiT0+g+cc(10)where cc denotes the complex conjugation of all preceding terms of the right-hand side. A1 and A2 are considered as functions of the slow time T2 instead of being constants that are affected by the nonlinear terms. The other symbols can be expressed as Λn=2iωn−Ωxyωn2+Ωyy,B=−(ωs2+Ωyy)f−(2ωs−iΩxy)f2[(ωs2+Ωyy)(ωs2+Ωxx)−4ωs2−Ωxy2],C=−i(ωs2+Ωxx)f−(2iωs−Ωxy)f2[(ωs2+Ωyy)(ωs2+Ωxx)−4ωs2−Ωxy2](11)It is clear that the z component is linearly uncoupled from both x and y. Once the solar sail is pitched (sinγ≠0) and D is set to be zero, the resulting out-of-plane motion is in the form of a periodic oscillation displaced to an out-of-plane distance g. Because the motion in the z direction is uncoupled, we focus on the planar motion in the rest of the Note.A. Planar Forced-Oscillation SolutionsFor the planar motion, the forced-oscillation solutions that are the particular solutions of Eq. (6) can be easily derived. The linear steady-state forced solution is independent of the initial conditions, and it can be expressed as xp=BeiωsT0+cc,yp=CeiωsT0+cc(12)The complex amplitudes B and C of the forced-oscillation solutions are dependent on the system's excited frequency and the solar-sail acceleration f. Based on Euler's formula (ex=cosx+isinx), the imaginary parts of amplitudes B and C will show up as the phase angles in the x and y directions. As we can see from Eq. (11), the phase angles in the x and y directions are different, which can make the semimajor axis of the ellipse no longer in parallel with the x axis of the L4 synodic coordinate system. This phenomenon is similar to that observed for the first-order approximation solutions to the long- and short-period orbits in the classical CRTBP model [28].B. Planar Free-Oscillation SolutionFor the planar motion, the free-oscillation solution involves the terms containing A1, A2, Λ1A1, and Λ2A2 in Eq. (10). For general mechanical systems, the free-oscillation solution always decays to zero because of the presence of damping [29]. However, due to the absence of damping, it is reasonable to consider that both the free-oscillation term and forced-oscillation term appear at the same time.We focus on the free-oscillation term. Substituting Eq. (8) into Eq. (6) and equating the coefficients of equal powers of ε2 and ε3 on both sides lead to D02x2−2D0y2−Ωxxx2−Ωxyy2=α200x12+α110x1y1+α020y12,D02y2+2D0x2−Ωxyx2−Ωyyy2=β200x12+β110x1y1+β020y12(13)D02x3−2D0y3−Ωxxx3−Ωxyy3=−2D0D2x1+2D2y1+α300x13+α210x12y1+α120x1y12+α030y13+x2(2α200x1+α110y1)+y2(α110x1+2α020y1),D02y3+2D0x3−Ωxyx3−Ωyyy3=−2D0D2y1−2D2x1+β300x13+β210x12y1+β120x1y12+β030y13+x2(2β200x1+β110y1)+y2(β110x1+2β020y1)(14)Substituting x1 and y1 of Eq. (10) into Eq. (13), the solution to Eq. (13) is obtained. Then, substituting Eq. (10) and the solution of Eq. (13) into Eq. (14) yields D02x3−2D0y3−Ωxxx3−Ωxyy3=−2D0D2x1+2D2y1+S1e3iω1T0+S2e3iω2T0+S3e3iωsT0+S4eiω1T0+S5eiω2T0+…+S13ei(ωs−2ω1)T0+…D02y3+2D0x3−Ωxyx3−Ωyyy3=−2D0D2y1−2D2x1+W1e3iω1T0+W2e3iω2T0+W3e3iωsT0+W4eiω1T0+W5eiω2T0+…+W13ei(ωs−2ω1)T0+…(15)where S4, S5, S13, W4, W5, and W13 can be found in the Appendix. Note that the other coefficients (Si and Wi) will be found as not useful in a later analysis.Assuming one-third subharmonic resonance, we introduce the detuning parameter σ defined by ωs=3ω1+εσ(16)and express (ωs−2ω1)T0 in Eq. (15) as (ωs−2ω1)T0=ω1T0+σT2.By the solvability conditions of Eq. (15), we seek a particular solution of x3 and y3 as x3=U1eiω1T0+Q1eiω2T0,y3=U2eiω1T0+Q2eiω2T0(17)Equation (17) contains the secular terms that make x3 and y3 become unbounded as T0 goes to infinite, unless the secular terms are eliminated. By substituting Eq (17) into Eq. (15), using (ωs−2ω1)T0=ω1T0+σT2, and equating the coefficients of exp(iω1T0) and exp(iω2T0) on both sides, we obtain (−ω12−Ωxx)U1+(−2iω1−Ωxy)U2=V11,(2iω1−Ωxy)U1+(−Ωyy−ω12)U2=V21(−ω12−Ωxx)Q1+(−2iω1−Ωxy)Q2=V12,(2iω1−Ωxy)Q1+(−Ωyy−ω12)Q2=V22(18)where V11=(−2iω1+2Λ1)D2A1+S4+S13ei(σT2),V12=(−2iω2+2Λ2)D2A2+S5V21=(−2iω1Λ1−2)D2A1+W4+W13ei(σT2),V22=(−2iω2Λ2−2)D2A2+W5(19)Thus, the problem of determining the solvability conditions of Eq. (15) is reduced to that of determining the solvability conditions of Eq. (18): |−2iωn−ΩxyV1n−Ωyy−ωn2V2n|=0or V1n=−Λ¯nV2n,n=1,2(20)Substituting Eq. (19) into Eq. (20), we obtain D2A1=−iΓ1(Λ¯1W4+S4+(Λ¯1W13+S13)ei(σT2)),D2A2=−12iΓ2(Λ¯2W5+S5)(21)where Γn=1ωnΩyy+ωn2Ωyy+Ωxx−4+2ωn2(22)Because S4, S5, S13, W4, W5, and W13 are all complex, Eq. (21) is rearranged as follows: D2A1=−iΓ1((R11+iI11)A12A¯1+(R12+iI12)A1A2A¯2+(R13+iI13)A1+(R14+iI14)A¯12ei(σT2))D2A2=−12iΓ2((R21+iI21)A¯2A22+(R22+iI22)A2A1A1¯+(R23+iI23)A2)(23)where R11,…,R23 and I11,…,I23 can be found in the Appendix. Note that Rij and Iij are all real functions of the solar radiation acceleration f and the coefficients of the nonlinear terms α200,…,α030, β200,…,β030.To solve Eq. (23), we write An (n=1,2) in polar form as An=(1/2)αneiβn, where αn and βn (n=1,2) are real. Then, separating the real and imaginary parts yields α1′=14Γ1I11α13+14Γ1I12α1α22+Γ1I13α1+12α12(Γ1I14cosγ1+sinγ1Γ1R14),γ1′=34Γ1R11α12+34Γ1R12α22+3Γ1R13+32Γ1R14cosγ1α1−32sinγ1Γ1I14α1+σ,α2′=18Γ2I21α23+18Γ2I22α12α2+12Γ2I23α2,γ2′=−18Γ2R21α22−18Γ2R22α12−12Γ2R23(24)where γ1=σT2−3β1 and γ1=β2.Furthermore, we can present the first-order approximation for the general solution [Eq. (10)] as x=0.5α1ei(ω1T0+(1/3)(σT2−γ1))+0.5α2ei(ω2T0+γ2)+BeiωsT0+cc,y=0.5Λ1α1ei(ω1T0+(1/3)(σT2−γ1))+0.5Λ2α2ei(ω2T0+γ2)+CeiωsT0+cc(25)where αn and γn (n=1,2) are given by Eq. (24). Note that Eq. (25) is quite different from the linear solution of the system. The amplitude of the linear solution can be chosen as any value. However, in Eq. (25), the variations of αn and γn (n=1,2) strictly follow Eq. (24), which includes the transient motions and steady-state motions, and is affected by the system's nonlinearity.C. Steady-State Motion for Planar Free Oscillation and its Existence DomainThe steady-state motions occur when α1′=α2′=0 and γ1′=γ2′=0, which correspond to the singular points of Eq. (24), where the amplitude and phase do not change [20]. The second steady-state modal motion demands α2=0. Then, based on Eq. (24), the steady-state motions corresponding to the solutions of the first modal motion (i.e., motion relevant to ω1) must be Γ1I11α12+4Γ1I13+2α1Γ1I14cosγ1+2α1sinγ1Γ1R14=0,34Γ1R11α12+3Γ1R13+32Γ1R14cosγ1α1−32sinγ1Γ1I14α1+σ=0(26)Rearranging cosγ1 and sinγ1, Eq. (26) leads either to α1=0 or (12α1(I142+R142))2=(6I14(I11α12+4I13)+2R14(3R11α12+12R13)+8R14σΓ1)2+(−6R14(I11α12+4I13)+2I14(3R11α12+12R13)+2I144σΓ1)2(27)which is an implicit equation for the amplitude of the response α1 as a function of the detuning parameter σ, which is determined by the frequency of the excitation, and the amplitude of the excitation f. It is called the frequency-response equation [20]. Once the frequency of the excitation ωs and solar radiation acceleration f are given, the amplitude of the response α1 is determined.With the relation that γ1=σT2−3β1, we can present the steady-state solution as x=0.5α1ei(1/3)(ωsT0−γ1)+BeiωsT0+cc,y=0.5Λ1α1ei(1/3)(ωsT0−γ1)+CeiωsT0+cc(28)where α1 and γ1 are given by Eq. (26). In the steady state for the free-oscillation term, the nonlinearity adjusts the frequency of the free-oscillation term to one-third of the excitation frequency. It is a typical subharmonic resonance, or frequency demultiplication [20]. In the resonance case, the value of the amplitude is decided by the frequency-response equation.In the following study, we try to find the domain in which the subharmonic resonance can exist. It is noticed that Eq. (27) is quadratic in α12. The solution of Eq. (27) can be expressed as α12=p±(p2−q)(29)where p=−(2(6I14I11+6R14R11)(24I14I13+24R14R13+8R14(σ/Γ1))+2(6I14R11−6R14I11)(−24R14I13+24I14R13+2I14(4σ/Γ1))−144(I142+R142)2)(6I14I11+6R14R11)2+(6I14R11−6R14I11)2,q=4(24I14I13+24R14R13+8R14(σ/Γ1))2+(−24R14I13+24I14R13+2I14(4σ/Γ1))2(6I14I11+6R14R11)2+(6I14R11−6R14I11)2(30)We note that q is always positive, and thus nontrivial free-oscillation amplitudes occur only when p>0and p2≥q(31)Because R11,…,R23 and I11,…,I23 are all real functions of excitation amplitude f, the nontrivial solution region can be located in the f−σ plane by applying the criteria in Eq. (31).IV. Numerical SimulationsIn this section, the Earth–moon system parameters are used in the numerical simulations to demonstrate the existence of novel subharmonic resonance periodic orbits. The boundary of the region where free-oscillation solution can exist is computed in the f–σ plane. A certain amplitude of steady-state motion of the free-oscillation term is determined by the system parameters.A. Existence Domain of Subharmonic Resonance Period OrbitBased on the criteria in Eq. (31), with variations of f and σ, the conditions p>0 and p2>q are given. In Fig. 3, the red area denotes the region that satisfies the criteria in Eq. (31). In other words, the novel subharmonic resonance periodic orbits can exist in such a domain.Fig. 3 Region where subharmonic responses exist.For a given σ, the nontrivial solution of the free-oscillation term exists within a certain range of values of the solar-sail acceleration f. For the solar sail in the Earth–moon system, the detuning parameter σ is determined by the first natural frequency ω1 (ω1=0.2982) and the sun's excitation frequency ωs (ωs=0.9252). Based on Fig. 3, when σ=60.0336 (σ≈ωs−3ω1), subharmonic responses exist under the condition that the corresponding excitation amplitude f belongs to [0,1.954×10−4].B. Frequency-Response CurvesBased on Eq. (27), several analytical frequency-response curves can be obtained with different values of f and are shown in Fig. 4. Note that, although the excitation f is relatively small, the response can be quite large. For example, when the excitation solar radiation acceleration f is chosen as 0.00005 and σ=0.0336, the response amplitude of the periodic orbit α1 has two values: 0.22 and 0.26, which are more than 4000 times larger than the value of the excitation f. It is a prototypical property of the resonance that a small input results in a large response due to energy accumulation. We would like to mention that one of the two values of α1 is generally unstable. The unstable value of α1 cannot be chosen to construct the periodic orbit.Fig. 4 Amplitude of the free-oscillation term.For the solar sail in the Earth–moon system, the amplitude of the free-oscillation term is plotted as a function of the amplitude of the excitation f when σ equals 0.0336 in Fig. 5. Similar to Fig. 4, based on Eq. (27), one input of f can result in two values of α1. We mark them as P1 and P2 in Fig. 5; and α1(P1)=0.0242 and α1(P2)=0.0092. Besides, in Sec. III.B, we also notice that α1=0 is another steady-state amplitude, which is marked as P3 in Fig. 5. The local Lyapunov stabilities of these values are evaluated by the state plane portrait [20]. The state plane portrait shows that a small change in the initial conditions can produces a large change in the response of the system. The center point and stable focus in the state plane indicate the stable value of α1, and the saddle point and unstable focus indicate the unstable value of α1.Fig. 5 Amplitude of the free-oscillation term versus amplitude of excitation with σ=0.0336.To illustrate stability, the trajectories corresponding to f=0.00015 and σ=0.0336 are computed with Eq. (24). The state plane for the free-oscillation solution is shown in Fig. 6. It is clear that the center point in Fig. 6 reflects the case for α1=0.024238 and γ1=3.325, which is the exact value of P1 in Fig. 5. Similarly, point P2 is a saddle point corresponding to α1=0.0092 and γ1=0.183, and P3 is a neutral equilibrium corresponding to α1=0 and γ1=any value (P3 can be any point on the straight line with α1=0). These points correspond to those labeled in Fig. 5. The arrows indicate the direction of the motion of the representative point in Fig. 6. Thus, all the shaded area constitutes the domain of attraction of point P1. The initial conditions must reside in the shaded area if the subharmonic appears in the first-order approximation of the steady-state response.Fig. 6 State plane for the subharmonic response when f=0.00015 and σ=0.0336.Therefore, when more than one stable steady-state solution exists for the free-oscillation solution, the initial conditions determine which steady-state solution is physically realizable by the system. For the Earth–moon–solar-sail system, when f=0.00015, both α1=0.024238 and α1=0 are stable steady-state motion amplitudes. The desired steady-state solution can be obtained through selecting the appropriate initial condition.C. Periodic and Quasi-Periodic Solar-Sail Trajectories Around L4 PointIn this subsection, we present possible periodic and quasi-periodic orbits that can be used as nominal orbits for practical space missions. Based on the analysis and results in Sec. IV.B, we can conclude that, when the solar acceleration f and detuning parameter σ are within the red region in Fig. 3, both the free-oscillation solution with α1≠0 or α1=0 and the forced-oscillation solution can be considered for planar motion, and the subharmonic resonance period orbit exists. When the solar acceleration and detuning parameter are within the blue region in Fig. 3, only the forced-oscillation solution should be considered for planar motion. The illustrative three cases are studied in this subsection. For case 1, σ=0.0336, f=0.00015, and α1=0.024238. In this case, f and σ are within the red region in Fig. 3. We plot the periodic motion to show the synthesis of the response in Figs. 7a–7d. Although there are two frequencies in the planar motion, due to the contribution of nonlinearity, the frequency of the free oscillation becomes exactly one-third of the frequency of the forced oscillation, which results in the final frequency for the actual response appearing as 0.3084. In Figs. 7–10, the red star denotes the location of the Earth–moon L4 point. It is a new type of periodic solution in the solar-sail Earth–moon restricted three-body problem. It is worth noting that, with different values of f, based on Fig. 5, we can always find the corresponding stable steady-state motion amplitudes α1. Then, a family of novel subharmonic resonance period orbits can be found.Considering the linear motion in the z direction in Eq. (10), two types of spatial trajectories are obtained. With parameters D=0 and g=0.01, the displaced subharmonic resonance periodic orbit is shown in Figs. 8a. When D=0.01 and g=0 are chosen, the orbit becomes subharmonic quasi periodic, and is shown in Fig. 8b, because the out-of-plane frequency equals one and the planar frequency equals 0.3084. We can also set D≠0 and g≠0 at the same time, which results in a displaced subharmonic quasi-periodic orbit sharing the same configuration as the orbit in Fig. 8b, which then hovers above or below the L4 point with distance g.For case 2, σ=0.0336 and f=0.00025. In this case, f and σ are within the blue region in Fig. 3. Only the forced-oscillation solution exists in this case. We plot the response for case 2 in Figs. 9a and 9b. By comparing Figs. 9a and 7b, we can conclude that, with the increase of the solar radiation pressure acceleration f, the amplitude of the forced-oscillation solution increases. The trajectory of the forced-oscillation solution is an ellipse. Due to the phase differences in x and y directions, the semimajor axis of the ellipse is no longer in parallel to the x axis of the L4 synodic coordinate system. The simulation result is consistent with the theoretical analysis in Sec. III.A.Adding the motion in the z direction in Eq. (10), two types of spatial trajectories are obtained. When D=0 and g=0.01, the periodic trajectory that hovers above the L4 point with a distance of 0.01 is shown in Fig. 10a. When we choose D=0.01 and g=0, the quasi-periodic orbit is obtained as shown in Fig. 10b. We can also set D≠0 and g≠0 at the same time, which results in a quasi-periodic orbit that hovers above or below the L4 point with distance g.For case 3, σ=0.0336, f=0.00015, and α1=0. In this case, f and σ are within the red region in Fig. 3. Because α1=0, the resulting orbits will be the same types as the orbits in case 2. We do not provide the results of case 3 for brevity.In addition to the periodic and quasi-periodic orbits we presented previously, the pure free-oscillation periodic solution can also be generated as a nominal orbit for space missions with additional control force to balance the forced-oscillation solution, just like the method proposed by Heiligers et al. [26].Fig. 7 Synthesis of the response of the solar sail in the Earth–moon system to subharmonic resonance.Fig. 8 Two types of spatial subharmonic resonance trajectories.Fig. 9 Free-oscillation trajectory of the solar sail in the Earth–moon system.Fig. 10 Two types of spatial forced-oscillation trajectories.V. ConclusionsIn this study, the periodic and quasi-periodic motions around the L4 point of the Earth–moon–