On the effect of the Sun on Kordylewski clouds

,


Introduction
The presence of Trojan asteroids near the triangular points of the Sun-Jupiter (S-J) system leads to a natural question: Are there similar objects in the vicinity of the Earth-Moon (E-M) triangular points?The first observational work addressing this question appeared in 1961 when the polish astronomer K. Kordylewski published a paper [Kor61].He reported the presence of a region populated with small particles near the triangular point L 5 at the Earth-Moon system.That type of dust region in the E-M system has been called Kordylewski cloud.
During the subsequent years, several authors tried to reproduce the observations of Kordylewski with positive [Roa75] and negative [FV80,FV82,Win89] results.These ambivalent results put the existence of Kordylewski clouds under the shadow of doubt.The letter by R. G. Roosen and C. L. Wolf [RW69] (1969) opened the debate to the scientific community lurching between the existence and non-existence of Kordylewski clouds.It is clear that their existence cannot be followed by analogy with the existence of Trojan asteroids.Indeed, the Restricted Three-Body Problem (RTBP) seems to be a reasonable model for the S-J system as these two bodies represent more than the 99% of the total mass of the solar system and perturbations due to the gravitational fields of other planets are small.On the other hand, the motion in the Earth-Moon system is severely affected by the presence of Sun.
The first paper proposing a model for the study of the libration clouds [Poh64] considered, however, a modified Kepler problem with Earth as central mass.This work, due to Pohle, was published in 1964, 5 years before the letter by Roosen and Wolf.The models got more complicated as the controversy increased.Katz, in 1975 [Kat75], presented a numerical study using a model considering Sun's gravitational field together with Earth's and Moon's.Sun in this model is frozen at infinity but this suffices to capture the unstable nature of the E-M triangular points.It is also remarkable that the author considered the importance of Solar Radiation Pressure (SRP) in the problem.Indeed, as the libration clouds are to be formed by small particles, SRP should be considered as a relevant force.In 1979, Burns, Lamy and Soter [BLS79] published a seminal paper on the effect of different kinds of radiation forces on small particles in the solar system.In this work it is analysed how SRP affects the dynamics of dust.Also the impact of other radiation forces as the solar wind, the Poynting-Robertson effect and the Yarkowsky effect are analysed.Other references concerning the relation between SRP and dust are [KOM02,Vin09,Sri99,LC15,LCG16].
In this work we consider several models that include the gravitational effects of Sun, Earth and Moon on a dust particle with the dissipative PR effect.Although very realistic models have the obvious advantage of being accurate, they are very complex and difficult to understand: as the model usually includes different effects, it is non-trivial to find the relevance of each of them on the final properties of the system.The approach we have used here to deal with these difficulties is to construct an scaffolding of increasingly accurate (and increasingly complicated) models.In each of these models we will focus on some dynamical structures and uncover new features that we study in detail.This detailed study of features of one model (which involves a combination of analytical and numerical methods) can be taken as the starting point of another study in the next model in our hierarchy of accuracy.At each stage new phenomena appear.At the end, therefore, we uncover some phenomena present in the most realistic model but not present (even in a qualitative form) in the most simplified one.

Organisation
The manuscript starts with a description of the proposed model in Section 2; first describing the Hamiltonian part and after the dissipative contribution part.Section 3 describes the dynamical equivalents of the triangular points in the proposed model and it is followed, in Section 4, by a study of effective stability regions nearby.Section 5 discusses the same regions in a restricted N -body problem including the solar radiation pressure and the Poynting-Robertson drag effect.Finally, the paper ends with a conclusion, Section 6, and some future works.Some technical details can be found on Section 7.

The ABCP model
To study the influence of SRP upon a test particle we first make use of a simplified model; called Augmented Bicircular Problem (ABCP).The ABCP is a modification of the well-known Bicircular Problem (BCP) to account for Solar Radiation Pressure (SRP) and Poynting-Robertson effect (PR).
In our first simplified setting, it is assumed that the only forces acting on the particle are the gravitational pull of three primaries: the Moon, the Earth and the Sun.As we assume that the dust particle is spherical, SRP acts as a force which is opposite to Sun's gravity and PR is a force opposite to the direction of the motion, see Figure 1 for a depiction of the model.
Further standard considerations on the motion of the primaries are the following: The Earth and the Moon are assumed to revolve, around its common centre of masses, describing a circular orbit.At the same time, the centre of masses of the Earth-Moon system revolve together with the Sun, around the barycentre of the system, describing again a circular orbit.Moreover, we assume that motion of the three primaries (the Moon, the Earth and the Sun) is contained in a plane.Motion on the orthogonal direction to this plane will be ignored in the present work.
A synodical frame of reference is taken, so that Earth and Moon are fixed on the horizontal x axis.The plane px, yq contains the motion of the three primaries.The unit of distance is the one between Earth and Moon, the unit of mass is the sum of Earth's and Moon's masses and the unit of time is taken so that the period of Earth and Moon around its centre of masses is 2π.In these units the universal gravitation constant equals one.Hence, we will denote by µ the mass of the Moon, m S the mass of the Sun, ω S the mean angular velocity of the Sun in these synodic coordinates and a S the distance between the Sun and the Earth-Moon barycentre.
We notice that the units and the frame of reference are the standard ones for the Earth-Moon circular Restricted Three Body Problem (RTBP).Sun's gravity and SRP appear then as periodic-time dependent perturbations of period T S " 2π{ω S .
The ABCP equations of motion are discovered in two steps.Firstly, incorporating a SRP to the BCP which will still keep the Hamiltonian structure, Section 2.1; and secondly adding a dissipation term due to the PR, Section 2.
Figure 1: Sketch of the ABCP.Primaries located at M " p1 ´µ, 0q, E " pµ, 0q and S " px S , y S q where x S " a S cos θ and y S " ´aS sin θ.Perturbed Lagrangian points L j for j " 1, . . ., 5. PR force opposite to particle motion (whose direction is denoted by the red arrow) and SRP opposite to Sun's gravity.

Hamiltonian part: BCP and SRP
Let H be a 2 1 2 degrees of freedom Hamiltonian of the form: Hpx, θ, βq " H RT BP pxq `HCSG px, θq `HSRP px, θ, βq, where x " px, y, p x , p y q.The first part H RT BP is the Hamiltonian of the Earth-Moon RTBP describing a circular orbit; the second part H CSG is a time-periodic perturbation that incorporates the Sun's gravitational field which also describes a circular orbit w.r.t. the E-M barycentre.These two first parts conform the so-called BCP Hamiltonian, see [Hua60, CRB68, SGJM95] for more details.Finally, the last term incorporates the effect of SRP and depends on a parameter β, called lightness number.
Therefore, the equation of motion of a dust particle in positions-momenta coordinates x " px, y, p x , p y q is described by the Hamiltonian in (1): where p x " 9 x ´y, p y " 9 y `x are the momenta, θ " ω S t is the angle of the E-M frame w.r.t. the Sun, the r's denotes the distance from the particle to each of the primaries respectively, i.e. r 2 P E " px ´µq 2 `y2 , r 2 P M " px ´µ `1q 2 `y2 , r 2 P S " px ´xS q 2 `py ´yS q 2 and β is a parameter.The adimensional parameter β is the ratio between the radiation pressure force and the force due to solar gravity.This means that when β " 1 the effects of SRP and Sun's gravity cancel each other.In some situations [LC15] (such as the two body problem and the RTBP with the Sun being one of the primaries) it is frequent to say that the particle feels a reduced gravity, as if the Sun were lighter.This is not the situation here.Indeed, while it is true that the magnitude of Sun's gravitational force is reduced due to the SRP, the centre of mass of the Earth-Moon system keeps moving according to the real mass of the Sun and therefore, the Coriolis force due to the rotation of the Earth-Moon system around the Sun is not cancelled out.That is, while the direct gravitational effect of Sun's gravity on the particle is cancelled when β " 1, the indirect effect is not.
In [JCFJ18] it is shown that the Coriolis acceleration due to the motion of the Earth-Moon barycentre, cancels out the linear term of Sun's gravitational potential in equation (2).Therefore, the term of lowest order in the perturbation, once SRP is turned on, is This means that SRP is the leading effect of the perturbation.For β " Op1{a S q " 10 ´3, the magnitude of SRP is already comparable to the one of Sun's gravity.For larger values, the SRP is the most important term in the perturbation.As β increases, the role of the Coriolis term takes importance.When the SRP and Sun's gravity cancel out (at β " 1), the only remaining term in the perturbation is the Coriolis acceleration.The parameter β depends on the area, the mass and the reflectivity of the particle.For a perfectly absorbing spherical body this coefficient can be approximated [BFM94] by 0.2{s, where s is the radius of the particle in micrometers.Reflectivity increases the values of β.In this work, we focus on extremely small values of β.
When β " 0 the model reduces to the well-known BCP: Near the triangular points, the perturbation due to Sun's gravity is large enough to produce a bifurcation on L 4 .It is replaced by three periodic orbits, one small and mildly unstable and two larger and stable.Due to the presence of the unstable periodic orbit, the neighbourhood of L 4 is no longer stable.This is a feature that BCP shears with the real System (by real systems we mean the model that takes into account the whole solar system, where the trajectories of the masses are given by the JPL ephemeris), so it is a more realistic model than the restricted three body problem (RTBP).
The third term in (1) of the Hamiltonian is a new periodic perturbation with the same period as Sun, so, if we set β ą 0 and small, the periodic orbits corresponding to the BCP are changed with β.
Notice that, the most easily detectable particles are the bigger ones.These particles, in principle, shall have larger mass and lower lightness number.In this regard, it would seem that the inclusion of SRP could be irrelevant for practical purposes but, as we will show, even small values of the parameter β play an important role.

Dissipative part: Poynting-Robertson effect
The Poynting-Robertson (PR) drag [BLS79] is an effect related to the solar radiation pressure, specially important for small particles, like interstellar or interplanetary dust.This effect is typically studied for idealised spherical dust orbiting a star.Such that, the dust particle absorbs radiation coming from the star in the radial direction, and then it re-emits the radiation making the particle to loss momentum and eventually spiralling into the star.Notice that, the re-emission would be isotropic seen from the reference frame of the particle, but it is an anisotropic effect seen from the star.
Being (X, Y ) the positions of the particle seen from the star, and ( 9 X, 9 Y ) its velocities, the Poynting-Robertson dissipative force is typically written in an inertial system with origin at its centre of masses, [MD99], as where r is the distance between particle and star, and κ is a parameter that is directly proportional to the radiation parameter, β, and to the mass of the star, M , and inversely proportional to the semi-major axis of the orbit of the particle, a, that is, being G the universal gravitational constant.In our case, the reference frame is the Earth-Moon synodic one, with the Sun orbiting their centre of mass in circular motion.Then, for a particle in this system, the solar radiation pressure is affecting the particle in the radial direction, and the re-emision of radiation due to the Poynting-Robertson drag is again anisotropic.
In order to include the drag force (3) in the model, we need to express the magnitudes involved in the reference frame and units of the model introduced in previous section.Then, the distance between particle and Sun r " r P S " ?X 2 `Y 2 , and the positions and velocities with respect to the Sun (that now is moving) correspond to Assuming that the dust particle is in the neighbourhood of the Earth-Moon system, the semi-major axis of its orbit with respect to the Sun can be taken as the distance between Earth-Moon barycentre and Sun, a S .So the parameter κ in our model, where G equals 1, takes the form With all above considered, the equations of motion for the dissipative model involving Earth-Moon Sunperturbed Bicircular problem under solar radiation pressure and Poynting-Robertson effect in the planar case are written as x " 2 9 y `x ´1´µ x `y ´1´µ Notice that, if we neglect these terms (i.e.formally κ " 0), the equations of motion in Equation ( 7) are the same ones as in the previous model, including only the gravitational effect of the three massive bodies and the solar radiation pressure.

Discussion on the dissipative terms
The dissipative parameter κ in Equation ( 6) linearly depends on the adimensional radiation parameter β, that is the ratio between radiation pressure force and solar gravitational force acting on a dust particle.Making some analysis and assumptions like the particle being a perfect sphere [MD99], one can express β in terms of physical quantities: where L and M are the luminosity and the mass of the star, in this case, the Sun, G is again the universal gravitational constant, c is the speed of light, Q pr is the radiation pressure coefficient, ρ is the density of the particle and s its radius.
Using the values in centimetre-gram-second system of units (CGS) for the physical known magnitudes (L, M , G and c), Equation (8) turns into β " 5.7 ˆ10 ´5 Qpr ρs .It is also standard to take Q pr " 1 and ρ " 3 g{cm 3 , what leads to the relation β « 0.2{s, with s in micrometers, mentioned in previous section.
In spite of the analysis of β is typically made in terms of CGS units in astronomy contributions, it is important to remind that it is a dimensionless parameter, what means that its value is the same regardless the model and units we use.Then, β admits values among 0 and 1.Now, examining Equations ( 3) and ( 6), we see that the dissipative force is proportional to the product β m S , where in our ABCP model m S « 3.289 ˆ10 5 and β P r0, 1s.Also, the first terms of the dissipative force are divided by r 2 P S a 2 S , being r P S and a S both distances close to one astronomical unit.With a S « 388.811 in adimensional units of ABCP, and r P S depending on the relative position to the Sun of the particle, that is in the Earth-Moon system.Then, since these first terms of the PR drag are divided by large quantities, it is expected that they do not affect much to the motion of the particle when the dissipative force is introduced in our model.Regarding to the last term of the dissipative force, it is expected to have even less relevance since it is divided by r 4 P S a 2 S .We note that the ABCP model is a Hamiltonian system (that depends on time in a periodic way) with a very small dissipative perturbation and a single parameter β (from which depend the magnitude of both the SRP and PR).Hamiltonian systems cannot have attractors and, in regions with dynamics close to integrable (like the neighbourhood of elliptic points, periodic orbits or invariant tori) they present a large set of quasi-periodic motions [JV97a,JV97b].However, any small dissipative perturbation can destroy these quasi-periodic motions and create attractors.As we will see, the attractors created by the Poynting-Robertson effect are very small to be relevant in this situation.

Dynamical equivalents of the triangular points
The ABCP model can be seen as a perturbative system of an autonomous part corresponding to the equation of motion of the RTBP with primaries E-M and another perturbative part involving periodic and dissipative effects.
To introduce dynamical phenomena in such a setting, let us consider a general framework given by an ODE of the form # 9 x " f pxq `εgpt, xq, xp0q " x 0 .
Here, f and g are assumed to be smooth enough, g to be T periodic with respect to the first component, that is gpt `T, xq " gpt, xq for each pt, xq and some minimal T ą 0 and ε a perturbing parameter.
Let us suppose that x is an equilibrium point of the unperturbed system (i.e.ε " 0, f pxq " 0).If ε is different from zero, then x is no longer an equilibrium point of (9).However, under generic conditions it can be continued to a periodic orbit of period T if |ε| is small enough.In this context, continued means that, there exists a family of periodic orbits X ε (with X ε pt `T q " X ε ptq), such that X ε Ñ x as ε Ñ 0. The members of the family X ε are called dynamical equivalents of x.This naming is set to denote that periodic orbits play a similar role as the equilibria in the perturbed system.In particular, periodic orbits are the simplest invariant objects there exist for ε ‰ 0 and, moreover, provided that ε is close enough to zero, X ε has the same stability character as x.
In the ABCP model, the autonomous part (f ) correspond to the equations of motion of the Earth-Moon RTBP in the synodical frame and the periodic perturbation (g) encapsulates all the effects due to the Sun (namely, the gravity, SRP and PR).Therefore, in the ABCP, the Lagrangian points existing in the RTBP are replaced by a family of periodic orbits that depend on the parameters β and κ.Note that, physical interpretations of the PR and SRP the two parameters β and κ depend on each other.
Systems such as the ABCP or, in general, (9), i.e, non-autonomous with periodic time dependence, are better handled by means of the so called stroboscopic map P ε .This map is defined for any initial condition within the definition domain of the ODE and is obtained by evaluating the solution of the equation at this initial condition at time T (the period of the perturbation).By definition, an initial condition on a periodic orbit is a fixed point of P ε , i.e.P ε pX ε p0qq " X ε p0q.The stability of a periodic orbit, therefore, is studied as the stability of a fixed point of the stroboscopic map P ε (this fact can be recovered from Floquet Theory).The linear normal behaviour around a fixed point X ε p0q is obtained by studying the spectrum of the Jacobian at the fixed point, which is usually called monodromy matrix of the orbit X ε .An eigendirection of M is said to be stable if its related eigenvalue has modulus less or equal than one.If all the eigenvalues of M are less or equal than one, then the periodic orbit X ε is said to be stable.Otherwise, is unstable.
When the ODE (9) is Hamiltonian (in our case when PR is not included), the stroboscopic map has symplectic structure and the eigenvalues of the monodromy matrix are constrained to certain form.In particular, the eigenvalues are split in pairs of opposite by the product, namely, pλ i , λ ´1 i q.These pairs are further classified according to their modulus.An eigendirection with eigenvalue λ is labelled as: Elliptic (or centre) if λ " exppiαq; Hyperbolic (or saddle) if λ is real with |λ| ą 1 and Complex Hyperbolic (or complex saddle) if λ " ρ exppiαq with ρ ą 1.In this case, the complex conjugate ρ expp´iαq is also an eigenvalue and therefore complex saddles are determined by a quartet of eigenvalues.The linear normal behaviour of fixed point for a symplectic stroboscopic map can be classified by these pairs of eigenvalues.In the case of a 4-dimensional map, there are, generically, four types of fixed points according its stability type: CentreˆCentre, CentreˆSaddle, SaddleˆSaddle and Complex Saddle.
The dynamics of the stroboscopic map near the fixed points is determined by its stability character.In particular, unstable directions (Saddles and Complex Saddles) have some unstable direction attached to the fixed point.
The effect of Sun's gravitational effect on the triangular points is a well studied problem.Indeed, if we set β " 0 in the ABCP, we obtain the standard BCP formulated in [Hua60].It is common knowledge that Sun's gravity destabilises the triangular points in the Earth-Moon system [GLMS87,DJS91].In Figure 2 we display how L 4 is affected by Sun's gravity in the context of the BCP.Panel (a) displays the characteristic curve of fixed points that replace L 4 .The vertical axis displays ε, an artificial parameter that multiplies m S in equation (7) (we remind that β " 0).So that the model reduces to the RTBP for ε " 0 and to the BCP for ε " 1.Notice that the curve crosses one time with ε " 0. This crossing point corresponds to L 4 .The three crossing points with ε " 1 (the top of the figure) corresponds to the three periodic orbits that replace L 4 in the BCP.There are three dynamical replacements of L 4 due to a broken pitchfork bifurcation.The points of the characteristic curve are coloured after its stability character.Points in purple are totally elliptic and points in green are CentreˆSaddle.In Figure 2, panel Notice the small periodic orbit close to the coordinates of L 4 is CentreˆSaddle and, hence, unstable.
(b), we display the projection of the thee orbits onto the configuration space x ´y.The colour of the trajectories, again, correspond to their stability character.
As we have already pointed out, the effect of Sun's gravity is well-known but, in this work, we are interested also on the effect of SRP.In [JCFJ20], it is shown that SRP has a remarkable impact on the stability and the location of the dynamical equivalents of the triangular points.Notice that these two aspects are crucial when looking for regions of effective stability in the Earth-moon system.Indeed, these regions are expected to exist around totally elliptic periodic orbits.Therefore, how the dynamical equivalents change with respect to β provides a way to discard large portions of the phase space for each value of β.
In Figure 3 panels (a) and (b) we show the evolution of the three dynamical equivalents of L 4 (described before, for the standard BCP) as they change with respect to the parameter β.The vertical axis displays the y-component of the orbits as fixed points of the stroboscopic map and the horizontal axis shows β.Panel (a) focuses on small values of β.Notice that the characteristic curve crosses the vertical line β " 0 three times, which corresponds to the three periodic orbits shown in Figure 2. The purple points indicate totally elliptic character while green colour is used to display fixed points with a of character saddleˆcentre.A number of bifurcations take place as the value of β changes.Those values are located in Figure 3 with labels β i .
The characteristic curve undergoes to saddle-centre bifurcations at β 1 « 6.3 ˆ10 ´6 and β 2 « ´5.4 10 ´5, those are the values for which the curve turns from purple to green and vice-versa.The negative values of β are physically meaningless but mathematically significant as they allow us to understand that the three periodic orbits of the BCP are connected by continuation with respect to β.For β 3 « 2.7ˆ10 ´3, the curve of fixed points undergoes a period-doubling bifurcation.This type of bifurcation occurs when a pair of eigenvalues collide at ´1.A family of totally elliptic 2-periodic fixed points branches out.The main family becomes partially hyperbolic and, therefore, unstable.Although the main family loses its stability, the period-doubled family is stable.At β 4 « 4.3 ˆ10 ´3, both families merge again in a period halving bifurcation and the main family recovers its stability character.The stable, period-doubled, family is the In Figure 3, panel (b), it is provided a wider perspective on the behaviour of the characteristic curve.The value of the y-component decreases from its initial value for β " 0 (near 1) until it reaches the value y " 0. Before reaching the horizontal axis, the characteristic curve in Figure 3, panel (b), changes its colour from purple to magenta at a value of β labelled as β 9 " 2.1 ˆ10 ´2.This indicates a Hamiltonian-Hopf bifurcation of the family.In this type of bifurcation two non-conjugate eigenvalues collide at some point of the unit circle (a Krein collision) and get expelled from it.The aftermath of this collision is the fixed point becoming complex saddle.Similarly to the case of the period-doubling bifurcation, the characteristic curve recovers its totally elliptic character in another Krein collision at β 10 « 3.4 ˆ10 ´2.
The characteristic curve reaches, as said before, the horizontal axis at β 11 « 3.7 ˆ10 ´2.Notice that, due to the symmetries of the vector field, the characteristic curve corresponding to the dynamical equivalents of L 5 meets with the one of L 4 at y " 0.Moreover, the curve corresponding to L 3 also meets with them.At β 11 , therefore, the dynamical equivalents corresponding to L 3 , L 4 and L 5 merge in a a single dynamical equivalent, similarly to a pitchfork bifurcation.However, the bifurcation orbit (the orbit corresponding to β 11 ) is an orbit that collides with the Earth.Figure 4: Lyapunov families of invariant tori related to the dynamical equivalent of L 4 for β " 4.05492 10 ´1.The quasi-periodic motion is displayed for the stroboscopic map P β and each invariant torus is seen as a curve surrounding the fixed point that corresponds to the dynamical equivalent.
The stability type of the dynamical equivalents is key to the existence of stability regions.Indeed, totally elliptic points are generically surrounded by quasi-periodic motion (see [JV97a,JV97b]).This behaviour is illustrated in Figure 4.The word generically means that quasi-periodic behaviour surround the totally elliptic fixed points for almost all the cases.However, when an elliptic eigenvalue has the form exppi2πp{qq with p{q a rational number, the motion is resonant and, therefore, the quasi-periodic motion is destroyed.As β changes, the eigenvalues of the dynamical equivalents do so and, hence, their arguments take infinitely many rational values.In Figure 3 panel (c) we display the arguments of the eigenvalues of the dynamical equivalents of L 4 as a function of β.Notice that there are two curves (one for each pair of eigenvalue).The merging of the curves correspond to the first Krein collision at β 9 described before.The two curves split at β 10 .The coloured horizontal lines correspond to low order resonances.Low order resonances have a remarkable impact on the phase space near a fixed point.In fact, low order resonances are, sometimes, called strong resonances.The two crossing points green line (the 1 : 2 resonance) corresponds to β 3 (period doubling) and β 4 (period halving) which have also been described.Notice that there are three coloured lines corresponding to the 1 : 3 resonance (which takes place at β 5 « 1.2 ˆ10 ´2), the 1 : 4 resonance (at β 7 « 1.79 ˆ10 ´2) and a 1 : 5 resonance (at β 9 « 2 ˆ10 ´2).The dashed line corresponds an inner resonance (also displayed in Figure 3, panel (d)) where α 1 " 2α 2 and labelled as β 6 .
In Table 1 we have included a summary of bifurcation the characteristic curves undergoes with respect to β and their respective labels.
On the Poynting-Robertson effect: The influence of the PR effect on the dynamical equivalents is extremely small.The characteristic curve remain close to the one shown in Figure 3.The dissipation due to PR turns elliptic fixed points into attractors.However, this attracting character is remarkably mild.The attracting eigenvalues (i.e.eigenvalues whose modulus is smaller than one) differ from one by a quantity (depending on β) that is at most 10 ´6.This means that a particle attracted by the periodic orbit needs about fifty thousand years to half the distance to it.On the other hand, the PR breaks the symmetry along the horizontal axis.Therefore, the Pitchfork bifurcation that takes place at β 11 breaks when the PR drag is included.

Stability regions in the ABCP
This section is devoted to analyse effective stability regions around the dynamical equivalent of L 4 in the Bicircular Problem under solar radiation pressure for different values of β, with and without dissipative forces due to Poynting-Robertson effect.
In a first analysis, Section 4.1, Poynting-Robertson drag is not included.Later we will discuss its effect on the stability regions in Section 4.2.To finish the discussion about effective stability regions in the ABCP, we compare the dynamics found around triangular points L 4 and L 5 .

Stability regions in the BCP with SRP
For a given value of β, the dynamical equivalent of L 4 is seen as a point in the stroboscopic map P β (of integration time 2π{ω S ).Then, we define a mesh of initial conditions centred at the position of that point, varying x and y coordinates, in order to study their fate for a long period of time.This mesh constitute a cloud of test particles in the neighbourhood of L 4 that aims to reproduce the dynamics of a cloud of dust, like Kordylewski dust clouds, in said a region.
The procedure for defining the mesh is as follows.We vary the positions of L 4 in P β keeping its momenta, such that the initial conditions for, a priori, 10 6 test particles are generated, equispaced by a distance of 10 ´4 and centered at the point occupied by L 4 in P β .Since these points belong to P β we proceed with the integration of each test particle during 1300 applications of the stroboscopic map, that corresponds to a time close to a century.Notice that longer integration times in a simple model like this, would no affect much to the qualitative results.
At each integration step we check whether the particle collides with one of the primaries or if it escapes the Earth-Moon system.Since the primaries are point masses in our model, the collision with them is defined as being at a distance to those points lower than their radius.The escape is set as being at more than 10 Earth-Moon length units from their barycentre, since at this distance a particle would get free of the Earth-Moon gravity and would begin to describe a heliocentric orbit.In the fate colourmaps we have assigned purple colour to the collision with the Earth, red to the collision with the Moon, yellow to the escape and black otherwise.Therefore, black regions in our maps correspond to effective stability regions.Figures 6, 7 and 8 show some of the colourmaps for the fate of a cloud of particles around L 4 in P β for different values of β obtained in the ABCP without PR effect.
We recall that for β " 0, the model corresponds to the Earth-Moon Bicircular Problem, where the existence of effective stability regions around triangular points is already known, [CJ00].Then, as expected, for β close to zero, β " 5.3987 ˆ10 ´10 , first map in Figure 6, we find an effective stability region.Surrounding this region we find purple, red and yellow colours, what means that there are particle trajectories, close to the stable ones, colliding with the Earth, Moon or leaving the Earth-Moon system.
As β grows, this region gets pretty vast and enlarged.For this reason, the meshes needed to be wider for values of β among 10 ´4 and 3.5 ˆ10 ´3 so that the regions are contained in the maps.Following the discussion of Section 2.2.1, if we consider valid the expression relating magnitude of the SRP to the size of the particles, β « 0.2{s with s in µm units, then this big stability region in Figure 6 for the given range of values of β will correspond to cloud of dust particles of radius approximated among 60 and 2000 µm, or among 0.06 and 2 mm.
If we continue increasing β, the effective stability region gets smaller up to a values close to 1.2 ˆ10 ´2 for which it has almost disappeared, see second map in Figure 7.If we look at the continuation of L 4 in Figure 3 (b), we check that L 4 is elliptic for these values of β, then the disappearance of the stability region can not be related to a loss of stability character of the periodic orbit.It must be related to a resonance.Observing Figure 3 (c), we can relate resonances to the values of β.In this case, the resonance affecting the stability region is 1 : 3 denoted by β 5 in Table 1.Moreover, at these values of β we appreciate that the number of trajectories colliding with the Earth increases.This is easily understood looking at Figure 5, since for these values, the periodic orbits start approaching the Earth position.
After the closeness to disappearance of the effective stability region, it gains size again, with shape more or less rounded and enlarged, until it separates in three small regions for β " 1.58510 ˆ10 ´2.We observe these three separated regions for β up to 1.60558 ˆ10 ´2.This behaviour can be explained studying the ratio between the two normal frequencies of the periodic orbit.Observing Figure 3 (d) we see that for β 6 « 1.6 ˆ10 ´2, the ratio of the normal frequencies is 2.
After this division of the stability region, the three parts join again, resulting in one rounded black area.For a value of β 7 « 1.7 ˆ10 ´2, the dynamical equivalent of L 4 crosses resonance 1 : 4, see Figure 3 (c), that corresponds to the second colour map in Figure 8.And for β 8 « 2 ˆ10 ´2 the resonance 1 : 5 is crossed, making the stability region getting small and star-shaped, fourth map in the figure.After this resonance, at β 9 « 2.1 ˆ10 ´2 a Krein Collision is crossed and the dynamical equivalent of L 4 becomes a complex saddle, see Figure 3 (b).Consequently, no effective stability region is found until for β larger than 3.4 ˆ10 ´2, corresponding to β 10 , another Krein Collision is crossed.Then, the dynamical equivalent is elliptic again and a small stability region is found, surrounded by many trajectories colliding with the Earth (purple colour).Again, this is explained due to the closeness of the dynamical equivalent of L 4 to the Earth for high values of β, as shown in Figure 5.

Effect of Poynting-Robertson drag
Now we analyse the same initial conditions for different intensities of the radiation parameter in the ABCP including the Poynting-Robertson drag.In Figure 9, six colourmaps are shown for some of the same β values as in previous section.
As expected and discussed in Section 2.2.1, adding this dissipative force does not have a significant impact on the results.Indeed, the present initial conditions are integrated at most 100 years and recall for the drag to half the distance of the periodic orbit requires more than 50000 years (See comments on PR drag at the end of Section 3).
In the analysis of the colourmaps including the Poynting-Robertson, the only effect that we observe now is a stretching of the structures in the maps.Regardless this small difference, as β increases we observe the same evolution of effective stability regions than the one observed without this effect, including the disappearance or split of the effective stability region when a resonance is crossed or when the stability character of the dynamical substitute changes.

Comparatives between L 4 and L 5
A similar study done for L 4 can be reproduced for L 5 .In this case, we perform colourmaps corresponding to initial conditions taken around the dynamical substitute of L 5 in P β .To better compare the regions nearby L 4 and L 5 , we consider the ACBP without the PR since, as discussed, the PR produces stretching in the effective stability regions.
In order to study the neighbourhood of L 5 , we have performed parameter continuation of the dynamical substitute of L 5 for different values of β as described in Section 3 for L 4 .We found the same resonances and bifurcations as the ones shown in Figure 3 and Table 1 for L 4 .The meshes of initial conditions around L 5 were generated as in Section 4.1.Hence, Figure 10 shows some colourmaps of the neighbourhood of L 5 , for almost the same values of β as in previous section.The difference in some decimals is due to the continuation of the dynamical equivalent.
Notice that there is a symmetry between the triangular points in the ABCP: the dynamics forward in time for L 4 is equal to the dynamics backward in time for L 5 .Then, as one can expect, the effective stability regions found in the neighbourhood of L 5 are the same ones as the ones of L 4 for the same value of β.But, there are differences in the trajectories surrounding the stability regions.Indeed, we observe a higher density of purple and red colours in the maps.What means that a large quantity of the trajectories around the stability regions end up colliding with the Earth and with the Moon respectively.Since these maps are generated forward in time in the vicinity of L 5 , this means that backward in time we would observe several trajectories coming from the Earth and the Moon vicinities to the neighbourhood of L 4 .

Stability regions in a restricted N -body problem with SRP and PR
In previous section we analysed effective stability regions in the neighbourhood of triangular points for different intensities of the solar radiation pressure in the ABCP with and without Poynting-Robertson effect.Now, we aim to reproduce the colourmaps for a restricted N -body problem including solar radiation pressure with and without PR effect, where the initial conditions of the massive bodies involved in the integrations are taken from the JPL ephemeris file DE405 with origin at the barycentre of the Solar system, at the corresponding time.
Since we want to translate positions and velocities from an adimensional synodic reference frame, with origin at the Earth-Moon barycentre, to a dimensional inertial reference frame, with origin at the barycentre of the solar system, we need to apply a meticulous non-autonomous change of coordinates.This change, from adimensional coordinates (a) to ecliptical ones (e) is based in a rotation plus a translation, as the one in [GLMS85], e " kCa `b, (10 where k is a scale factor corresponding to the instantaneous distance between Earth and Moon in astronomical units (AU ), b is the distance to the Earth-Moon barycentre taken from the origin of the ecliptical  where R " R M ´RE and V " V M ´VE , are the position and the velocity of the Moon (R M and V M ) with respect to the Earth (R E and V E ).
It is worth to mention that, for translating velocities, we need to derive Equation (10) with respect to dimensional time, taken in days.
Then we apply this change of coordinates, for positions and velocities, to each initial condition of the maps in Section 4. Since these maps are defined at the temporal Poincaré section P β , and at this time, the configuration of the BCP corresponds to a lunar eclipse, we use the first lunar eclipse of year 2000, which corresponds to the modified Julian day 20.1978749133, for taking the ephemeris of the massive bodies, that are used in the change of coordinates and as starting conditions for their propagation in time.
Once we have initial conditions of both, test particles and massive bodies, we propagate them in time in a restricted N -body problem that accounts for the gravitation interaction among massive bodies and their effect on the massless particle, and also includes the solar radiation pressure and Poynting-Robertson drag on the latter.
The force due to solar radiation pressure in this spatial ecliptical model is written as: where M is the mass of the Sun, now in kilograms, r S its position, in AU and r " r S ´e the distance between particle and Sun in the ecliptical system.Subscripts x, y and z denote the component of position vectors r S and e.Also, we need to write spatial drag force due to Poynting-Robertson for this model that we shortly discussed in Section 2 for the planar case.Following [LZJ95] we write the spatial force as: being (X, Y , Z) the positions of the particle seen from the Sun, ( 9 X, 9 Y , 9 Z) its velocities, r is again the distance between particle and Sun, r " a X 2 `Y 2 `Z2 , κ is the same as in Equation (4), κ " βGM a 2 and a is de semi-major axis of the particle with respect to the Sun.Note that in ABCP we took this parameter as a S , that in this model is 0.999 AU , then, very close to 1.
In the vicinity of triangular points of the Earth-Moon system, the massive bodies playing the most important roles are the Earth, Moon and Sun, and in a fourth place, Jupiter.Then, our N -body problem is defined for these four bodies.
The presence of effective stability regions in a realistic model as the one we treat in present section, is very sensitive.The stability regions tend to disappear or dissipate for long integration times.Therefore, we decided to reduce the time span for this first analysis to 20 years to later study longer periods of time.Under all these considerations, we have found some effective stability regions.Even so, they are only found for low values of β, as seen in Figures 11 and 12.For β larger than, approximately, 1.2 ˆ10 ´3 no trace of stable regions is found.Note that for convenience, the horizontal and vertical axes correspond to x and y coordinates in ABCP units.Colours are assigned again as purple (collision with the Earth), red (collision with the Moon), yellow (escape) and black (otherwise).
The Figures 11 and 12 include the colourmaps obtained for different values of the radiation parameter at each row.The three columns correspond, from left to right, to the propagations in time performed in the neighbourhood of L 4 without including the Poynting-Robertson drag, the second column to the same initial conditions including this dissipative effect due to PR, and the third one to the neighbourhood of L 5 without this effect.They are shown together for easier comparisons.In these maps we observe that in spite of the fact that these regions are more sensitive to be dissipated in the realistic model than they were in the ABCP, the area they cover is pretty large.In fact, it was necessary to enlarge the meshes.
Regarding to the Poynting-Robertson drag, looking at the first two columns of maps in the figures, we appreciate no effect on the dynamics when this dissipative force is included.We do appreciate some differences in the comparison of the maps in the neighbourhood of L 4 and L 5 , first and third columns.The size of the effective stability regions is more or less the same, however, the ones around L 5 seem to be more dense than the ones around L 4 .Again, we also find more red colour, trajectories colliding with the Moon, in the case of L 5 than in the case of L 4 .
A remarkable difference among the results obtained in ABCP and in realistic models, is almost absence of trajectories colliding with the Earth (purple colour) in the latter, whereas in the former, their presence were important.It is worth to mention that these same maps of fate give the same qualitative results be Jupiter included or not.

Time analysis
In previous section, some effective stability regions in the neighbourhoods of triangular points composed by test particles that remain in the Earth-Moon system for 20 years were shown.Now, we aim to analyse the fate of these trajectories in the neighbourhoods of both L 4 and L 5 , for longer periods of time, 30 and 100 years, in order to understand what could be happening to the dust that is accumulated in these regions.In addition, propagations backward in time of the same initial conditions are included for the same reasons.
Figures 13 and 15 show some colourmaps for values of β of 2.084 ˆ10 ´4 and 8.049 ˆ10 ´4, respectively.Now, only SPR is included since PR was shown not to have relevance.The three rows in the figures correspond to integration time of 20, 30 and 100 years.The first two columns correspond to the neighbourhood of L 4 for propagations forward and backward in time, and the third column to the neighbourhood of L 5 for propagations forward in time.
Examining rows from up to down, we observe that as time evolves, be it forwards or backwards, the black regions, composed by stable trajectories, blurs.This does not mean that there are no stable trajectories left for time spans of 100 years.If we only plot the trajectories that have not collided with a massive body nor have left the E-M system, Figure 14, we can observe that there are still stable trajectories after 100 years of time propagation.However, it is clear that their presence has diminished.
Table 2 includes the percentages of the stable trajectories for the colourmaps after 20, 30 and 100 years of time propagation, both around L 4 and L 5 , forward and backward in time, for different values of β.Meanwhile the presence of stable trajectories in the maps after 20 years of propagation is around 1.5 ´2%, it is reduced up to lower than 0.1% after 100 years.
So it is clear that most of trajectories that were stable after 20 years end up colliding or leaving the system as time evolves.In the same table, last four columns, are devoted to understand what happened to those trajectories for a propagation of 100 years.In general, we observe that a low percentage of those trajectories remain stable.Around 6% of them collide with the Earth.A higher percentage is found for collisions with the Moon.But what is remarkable is that the vast majority of them, around a 80% of the trajectories that were stable for 20 years of propagation, end up leaving the Earth-Moon system.If we think about the propagation performed backward in time, this means that the vast majority of the particles that formed the effective stability region 20 years before come from the outside E-M system.
All this information allows us to conclude that, the effective stability regions seem to be mainly conformed by particles that enter our Earth-Moon system, remain in it for some decades, and then leaves it again.This could give an explanation for the presence of Kordylewski dust clouds as an accumulation of dust in slow passage close to triangular points.That is, they would not be conformed by the same particles trapped in those vicinities, but by particles that are in a continuous renewal.

Figure 2 :
Figure2: Numerical continuation of the Earth-Moon L 4 with respect to ε, an artificial parameter multiplying m S , so that, when ε " 0, the model is reduced to the RTBP and, when ε " 1, the model is the standard BCP.Panel (a) displays the x coordinate of the fixed point replacing L 4 as ε changes.The bottom of the plot corresponds to the RTBP and the top to the BCP.The points are coloured according to its stability character.Purple means centreˆcentre and green CentreˆSaddle.Panel (b) shows the periodic orbits of the BCP that corresponds to the three crossing points of the curve in (a) at ε " 1.Notice the small periodic orbit close to the coordinates of L 4 is CentreˆSaddle and, hence, unstable.

Figure 3 :
Figure 3: Behaviour of the dynamical equivalents of L 4 with respect to β. Panels (a) and (b) display the y-component of the fixed point as a function of β.Panel (c) display the normal frequencies (in purple) α 1 and α 2 as a function of β.Panel (d) displays the ratio between the two normal frequencies as a function of β.See text for more details.

Figure 8 :
Figure 8: Fate colourmaps for a cloud of particles around the position of L 4 in P β without Poynting-Robertson effect for values of β: 1.64210 ˆ10 ´2, 1.76828 ˆ10 ´2, 1.97527 ˆ10 ´2, 2.08566 ˆ10 ´2, 2.37762 10 ´2 and 3.40342 ˆ10 ´2.Horizontal: x coordinate.Vertical: y coordinate.Colour: purple to collision with the Earth, red to collision with the Moon, yellow to escape and black otherwise.

Figure 9 :
Figure 9: Fate colourmaps for a cloud of particles around the position of L 4 in P β including Poynting-Robertson effect for values of β: 8.10371 ˆ10 ´5, 7.76812 ˆ10 ´3, 1.18087 ˆ10 ´2, 1.60150 ˆ10 ´2, 1.76828 10 ´2 and 3.40342 ˆ10 ´2.Horizontal: x coordinate.Vertical: y coordinate.Colour: purple to collision with the Earth, red to collision with the Moon, yellow to escape and black otherwise.

Figure 10 :
Figure 10: Fate maps for a cloud of particles around the position of L 5 in P β without Poynting-Robertson effect for values of β: 8.25156 ˆ10 ´5, 7.80816 ˆ10 ´3, 1.18966 ˆ10 ´2, 1.60216 ˆ10 ´2, 1.76884 ˆ10 ´2 and 3.40368 ˆ10 ´2.Horizontal: x coordinate.Vertical: y coordinate.Colour: purple to collision with the Earth, red to collision with the Moon, yellow to escape and black otherwise.

Figure 13 :
Figure 13: Fate colourmaps obtained in a N -body model including SPR for β " 2.084 ˆ10 ´4.Columns: L 4 forward in time, L 4 backward in time, L 5 forward in time.Rows: for 20 years, 30 years, 100 years.

Figure 14 :
Figure 14: Stable trajectories for the analysis of the neighbourhood of L 4 backward in time during 100 years, obtained in a N -Body model including SRP for β " 2.084 ˆ10 ´4.Horizontal: x coordinate.Vertical: y coordinate.

Figure 15 :
Figure 15: Fate colourmaps obtained in a N -body model including SPR for β " 8.049 ˆ10 ´4.Columns: L 4 forward in time, L 4 backward in time, L 5 forward in time.Rows: for 20 years, 30 years, 100 years.See text for more details.

Table 1 :
Labels of the bifurcations of the dynamical equivalents of L 4 with respect to β. See text for more details.