On the Sun-shadow dynamics

We investigate the planar motion of a mass particle in a force field defined by patching Kepler's and Stark's dynamics. This model is called Sun-shadow dynamics, referring to the motion of an Earth satellite perturbed by the solar radiation pressure and considering the Earth shadow effect. The existence of periodic orbits of brake type is proved, and the Sun-shadow dynamics is investigated by means of a Poincare'-like map defined by a quantity that is not conserved along the flow. We also present the results of our numerical investigations on some properties of the map. Moreover, we construct the invariant manifolds of the hyperbolic fixed points related to the periodic orbits of brake type. The global picture of the map shows evidence of regular and chaotic behaviour.


Introduction
This paper deals with the study of a mass particle moving under the alternated action of two different force fields: one of Kepler's problem, the other of Stark's problem, where a constant force is added to the central one. This is a basic model to study the short-period evolution of an Earth satellite which alternately spends some time in the Earth shadow and some time in the region where the solar radiation pressure acts. We call this model Sun-shadow dynamics.
The solar radiation pressure (srp) can become the main perturbation to be added to the monopole term of the Earth, when the area to mass ratio of the satellite is large enough. The importance of the srp for accurately predicting the motion of artificial satellites of the Earth was first shown in [12] and [13]. Musen developed an analytic theory that was applied to the Vanguard I satellite. In [13] the motions of the Echo balloon and the Beacon satellite were numerically propagated. In both works it was found that the srp can seriously affect the lifetime of an Earth satellite, especially when a particular resonance condition is satisfied. Since these pioneering studies, the effects of srp have been investigated by many authors. For example, the long and short-period variations of the orbital elements are discussed in detail in [11].
A relevant aspect related to the srp perturbation is the passage through the Earth shadow. Kozai [8] was among the first authors to treat the eclipse perturbation and understand its effects. He developed a semi-analytic method to obtain the first-order variations of the orbital elements when the srp is switched off inside the Earth shadow. It is worth noting that, in general, the accumulation of these short-period effects produces a long-period drift of the semi-major axis, see [11]. Similarly, in [10] Lidov described a semianalytic method to compute the secular variation of the osculating parameters. His study showed that they all oscillate periodically, except in two limiting cases in which either the argument of periapsis or both the semi-major axis and the eccentricity vary monotonically. Also Ferraz-Mello [4] studied the possible secular effects induced by srp with the Earth shadow, by developing an analytic theory in the Hamiltonian formalism; he found that the angular drifts of the longitude, perigee, and node are quite small. A more recent paper, by Hubaux and Lemaître [6], showed that successive crossings of the shadow for long time spans (in the order of 1000 years) cause significant oscillations of the orbital elements, with amplitudes and frequencies that depend on the area-to-mass ratio. Moreover, numerical experiments carried out in [7] indicate that the passage through the shadow is a source of instability for space debris with high area-to-mass ratio at geostationary altitudes.
The idea which inspired this paper emerges in [2]. Beletsky proposed to apply Kepler's dynamics inside the Earth shadow, and Stark's dynamics outside of it. His qualitative analysis of the problem pointed out that the orbital energy has leaps each time the shadow is crossed and that the argument of periapsis, the semi-major axis and the eccentricity have long-period oscillations. More precisely, the semi-major axis decreases and the eccentricity increases while the apse line moves away from the direction of the srp force. In this work we try to go deeper into the subject. We consider the two-dimensional case where the srp force lies in the plane of motion of the satellites. After reviewing Stark's dynamics following [1] and [2], we describe some features occurring when we alternate it with the dynamics of Kepler's problem, using separable variables for the Hamilton-Jacobi equations of both problems. In particular, we prove the existence of periodic orbits of brake type, which are close to the unstable brake periodic orbits of Stark's dynamics. For a further description we introduce a Poincaré-like map S, that we call Sun-shadow map. For this purpose, we choose a section Σ in the boundary of the shadow region, and fix a quantity that has the same value when the section is crossed with the right orientation (but is not conserved along the flow). Then, we present the results of our numerical investigations: we describe the domain of S, prove that the fixed points related to the periodic orbits of brake type are hyperbolic, and compute the stable and unstable manifolds of these points. These manifolds are made of several connected components because there are orbits that either collide with the Earth or go to infinity, therefore they do not go back to Σ. Finally, a global picture of the map is drawn, showing evidence of regular and chaotic behaviour.
The paper is organized as follows. In Section 2 we introduce Kepler's and Stark's dynamics using separable coordinates for the Hamilton-Jacobi equations of the two problems. In Section 3 there is the review of Stark's problem. We investigate the alternation of the two different dynamics in Section 4, proving the existence of a family of periodic orbits. In Section 5 we define the Sun-shadow map and describe our numerical investigations.

Kepler's and Stark's dynamics
We consider Kepler's dynamics, defined bÿ with µ > 0 the Earth gravitation parameter and x = (x, y) ∈ R 2 . Moreover, we take into account Stark's dynamics, given bÿ Both equations (1) and (2) can be written in Hamiltonian form, with Hamilton's functions where p x , p y are the moments conjugated to x, y. Hereafter, the labels k, s will stand for Kepler and Stark, respectively. Besides H k , the angular momentum and the Laplace-Lenz vector are first integrals of Kepler's dynamics. Note that relation holds among these integrals. We denote by L k the opposite of the x-component of A k . On the other hand, besides H s , Stark's dynamics has the first integral which is a generalisation of L k (see [14]), but there are no other integrals independent from L s and H s .

Hamilton-Jacobi equations and separation of variables
The coordinate change (x, y) → (u, v) defined by separates the variables in the Hamilton-Jacobi equations of both Kepler's and Stark's problems. Relations (7) can be completed to a canonical transformation leading to new variables (p u , p v , u, v): Also a transformation of the time variable t can be performed by introducing the fictitious time τ through the differential relation Hereafter, we shall denote with a prime the derivative with respect to τ . Hamilton's functions for the two dynamics in these coordinates are Set U = (p u , p v , u, v) T , where T stands for vector transposition. Stark's and Kepler's dynamical systems can be written as where h k and h s are the values of H k and H s for some given initial conditions. The angular momentum in the new coordinates is while the integrals L k , L s become Hamilton-Jacobi equations for the two problems are where W k , W s are the unknown generating functions. In (15) the variables are separated, so that we obtain where α k is an integration constant. Let k be the value of the integral L k . Substituting p 2 u , p 2 v given by (17) into L k (p u , p v , u, v) = k and simplifying we get In a similar way, from equation (16) we get where s is the value of the integral L s .

Trajectories in Stark's dynamics
As explained in [1], all the possible trajectories of Stark's dynamics can be divided into four categories, depending on the values s , h s of L s , H s . Relations where A 1 and A 2 are integration constants and correspond to the expressions of p 2 restrict the possible configurations on the basis of the values s , h s . Let us set The polynomials U (u) and V (v) can be written as Setting the roots of U (u) are ±u 1 , ±u 2 , and those of V (v) are ±v 1 , ±v 2 . It is convenient to study the problem in the ( s , h s / √ f ) plane. Beletsky showed that this plane can be divided into four regions as shown in Figure 1. These regions do not cover completely the ( s , h s / √ f ) plane: in the remaining part (the brighter one in the figure) the motion is not possible. Each region is characterised by different types of trajectories listed in Table 1, see [1]. The admissible subsets of the configuration space are shown in Figure 2: these are delimited by straight lines in the (u, v) plane, and by parabolas in the (x, y) plane. For completeness we added in Appendix A Tables 2 and 3, describing the features of the trajectories at boundaries of the regions.
, +∞) trajectories type unbounded, self-intersecting, not encircling the origin , u ∈ (−∞, +∞) trajectories type unbounded, self-intersecting, encircling the origin In regions II and III there cannot be zero velocity points.

Remark 2.
For each region of the ( s , h s / √ f ) plane, the v-component of an orbit is periodic (and bounded). The u-component is periodic (and bounded) only if ( s , h s / √ f ) belongs to region IV and u 2 ≤ ξ 2 . On the contrary, the u-component is unbounded if ( s , h s / √ f ) belongs to one among regions I, II, III, or it belongs to region IV , and u 2 ≥ ξ 1 .
Proposition 1. In case u(τ ) and v(τ ) are periodic solutions of (12), relation Proof. The periods T u and T v of the u and v variables can be written as elliptic integrals: Let a u , b u , a v , b v be the positive square roots of the previous quantities. We note that Denoting by M (a, b) the arithmetic-geometric mean of two real numbers a, b (see [3]), we have that corresponds to (26).

Unstable periodic orbits of brake type
There exists a family of unstable periodic orbits of brake type, x * = x * (t; s ), parametrised by s ∈ (−µ, µ). It is possible to analyse the behaviour of the u and v-components of the trajectory in the reduced phase spaces with coordinates (u, p u ) and (v, p v ). For this purpose we can take into account the two Hamiltonian dynamics defined by obtained from system (18). H su has the two critical points (p * u , ±u * ), where We can show that they are two unstable equilibrium points for the reduced dynamics in the (u, p u ) plane. In fact, the Jacobian of the Hamiltonian vector field induced by H su , evaluated in both critical points, is At these critical points, the value of H su is Figure 3). Thus, the v-component is periodic. This implies that, if s ∈ (−µ, µ) and Stark's Hamiltonian H s has the value h * s , we have two unstable periodic orbits in the (u, v) plane. They have a constant value of the u-component, equal to ±u * , and a constant value of p u , equal to zero. Moreover, they are of brake type because each of them develops between two zero velocity points. These are given by (u * , ±v * 1 ), with v * 1 = v 1 (h * s ), in one case, and by (−u * , ±v * 1 ) in the other. In the (x, y) plane, they correspond to the same periodic orbit x * , which is a parabolic arc and develops between the zero velocity points (ξ Its trajectory is shown in Figure 4. The family of brake orbits x * exists for values s , h s corresponding to the boundary between regions II and IV . On this boundary, the values of ξ 1 and ξ 2 coincide and are equal to ξ * .

Other periodic orbits
There exists another family of periodic orbits of brake type at the boundary of region I, i.e. for s = −µ: they pass through the origin and have a constant value of the ucoordinate, equal to zero, in the (u, v) plane, therefore they lie on the y axis in the (x, y) plane. Moreover, there are periodic orbits of brake type in correspondence of s = µ and h s < −2 √ f µ. They also pass through the origin and lie on the y axis in the (x, y) plane, but in the (u, v) plane they have a constant value of the v-coordinate, equal to zero. If h s = −2 √ f µ, we obtain the two fixed points (u, v) = (± |h s |/f , 0) for Stark's dynamics, corresponding to a single fixed point in the (x, y) plane. Additional periodic orbits can exist for ( s , h s / √ f ) belonging to region IV , as a consequence of Remark 2. Their peculiarity is that the periods T u and T v are commensurable, that is their quotient T v /T u is rational. In Figure 5 some examples are shown. Note that the orbits in Figures 5(c) and 5(d) are of brake type. In this case, if T u is an odd integer multiple of T v , the trajectory passes through the origin as shown in Figure 5(d).

The Sun-shadow dynamics
The Sun-shadow problem arises by switching dynamics each time the satellite passes through the boundary of the Earth shadow. The flow develops by alternating Kepler's regime, corresponding to the shadow region, and Stark's regime, in the out-of-shadow region. As shown in Figure 6, the shadow region is defined as the set {(x, y) : At the time t = 0 of entrance into the shadow, let us consider initial conditions belonging to the set The last condition is needed to have the velocity vector pointing inside the shadow region. We search for the solutions (p u , p v , u, v) of the polynomial system where c k is the value of the angular momentum C k . From system (32), it is possible to obtain the state at the exit point of the shadow region. By eliminating the variables p u , p v , v we obtain an eight degree polynomial equation in u: The roots of (33) come in pairs ±u, which give the same values of x. We can select the right value of x using the y-component of the Laplace-Lenz integral A k . Let us call In this regime, the angular momentum can change not only in value, but even in sign. Indeed, the satellite can re-enter Kepler's regime either in the first or third quadrant of the (u, v) plane.
Proposition 2. Each time the satellite crosses the boundary of the shadow region, we have a leap in energy from h s to h k , or vice versa: the variation is equal to ±f (ū 2 − R 2 /ū 2 )/2, whereū is the value taken by u at the crossing point. A similar leap occurs from s to k , or vice versa. In this case, the variation is equal to ±f R 2 /2. When the satellite goes back to Stark's regime, the value of L s is the same as before crossing the shadow; on the other hand, the energy usually changes unless the orbit is symmetric with respect to the u axis.
Proof. Assume that the body enters Stark's regime at the point ( s be the value of the energy and i,1 s the value of the Laplace-Lenz integral. When the body returns to the shadow region, the integrals vary in the following way: with u o,1 the u coordinate of the point on the shadow boundary where the satellite exits from Stark's regime. When the body enters Stark's regime again, by passing through the point (p u i,2 , p u i,2 , u i,2 , v i,2 ), similar variations occur: Thus, we get

Periodic orbits of brake type
We prove the existence of a family of periodic orbits of brake type, x = x(t; s ), parametrised by s , which are close to the brake periodic orbits x * = x * (t; s ) of Stark's problem, described in Section 3.1. For this purpose, we consider values s in (−µ, µ), for which the periodic orbits x * exist. In the following we shall restrict the interval (−µ, µ) for reasons related to the proof. The idea of the proof is to search for an initial point (x, y) = (x 0 , 0) with x 0 > R, i.e. in Kepler's regime, allowing to arrive at a zero velocity point (x B , y B ) in Stark's regime after passing through an exit point (x E , R) from the shadow region. We look for an orbit that is symmetric with respect to the x axis, like x * . Because of the symmetry, x oscillates between the points (x B , y B ) and (x B , −y B ). We search for an initial value x 0 fulfilling where (ξ * /2, 0) belongs to x * , see Figure 7(a). The idea behind this choice is that, passing through the shadow, the pushing effect of the solar radiation pressure is lacking. Moreover, we require that x 0 is such that at the exit point (x E , R) the energy fulfils with h * s given in (29). If h s > h * s we cannot have zero velocity points, see Remark 1. For the proof, we use the variables u, v. The initial point (x 0 , 0) corresponds to two possible points (± √ 2x 0 , 0) in the (u, v) plane. Similarly, the point (ξ * /2, 0) corresponds to (± √ ξ * , 0). By symmetry, we can focus only on the {u > 0} half-plane of the (u, v) plane, see Figure 7(b).
Proposition 3. If x 0 is selected so that relations (34), (35) hold, then the exit point (x E , R) from the shadow region, with belongs to the unbounded component of the possible configurations corresponding to region IV of Stark's regime. Proof. Because of the symmetry of the orbit with respect to the x axis, in Kepler's regime the initial state has the form By system (32) we have and Note from relations (37), (38) that we need to restrict the interval of s to (−µ, µ − f 2 R 2 ). The exit point from the shadow is obtained by solving where p v > 0.
With the last equation we set to zero the y-component of the Laplace-Lenz vector, which is a necessary condition for the symmetry with respect to the x axis. There are four possible real solutions of (39), and two of them have positive values of u. Between these two, only one corresponds to exiting from the shadow, i.e. fulfils This solution gives the u coordinate at the exit point, whose square is with From ξ E , we obtain x E through relation (36), applying the coordinates change (7). Using relation (34), we can write We have real values of ξ E if and only if The last condition is fulfilled for each choice of ∆x ≥ 0 if which holds if we further restrict the interval of s to [ s − , s + ], with The energy in Stark's regime is For the proof we need this result: The energy h s is a decreasing function of x 0 in the interval [ ξ * 2 , ∞). Moreover, we have Proof. From equation (43), the derivative of h s with respect to x 0 is The denominator is always positive, being x 0 , x T > 0. We prove that the numerator is negative. Because of relations (29), (30) and (38), this corresponds to showing that This follows from (34) and x 0 > x T . We conclude that dhs dx 0 < 0. Next we prove the second statement of the lemma. We have

From relations (34) and (41), it holds
Thus, we have The polynomials g + and g − have three roots, but only one is larger than ξ * /2. Denoting the latter with x + 0 and x − 0 respectively, we get This shows the existence of x * 0 > ξ * /2, solution of h s (x 0 ) = h * s . Next we show that the exit point belongs to the unbounded component of the configuration set, i.e. that ξ E > ξ 1 holds for each x 0 > x * 0 . Since x − 0 < x * 0 < x + 0 and ξ E , given in (40), is an increasing function of x 0 , we have where Relation (42) implies ξ E (x − 0 ) > ξ * . Thus ξ E (x * 0 ) > ξ * . Since in region IV we have ξ 2 < ξ * < ξ 1 , and ξ E < ξ 2 or ξ E > ξ 1 , the latter relation holds.

Remark 5. Proposition 3 yields
To search for a zero velocity point (x B , y B ) we use the coordinates u, v and the fictitious time τ . The maps u(τ ), v(τ ) become stationary at τ = τ u , τ v respectively, where We search for value of x 0 such that which corresponds to reach a zero velocity point.
From now on, we use h s ∈ J = (−∞, h * s ) as independent variable, in place of x 0 . Following [1] we can write the integrals τ u , τ v as where We use the result proved below. Proof. i) The derivative of τ u with respect to h s can be written as Therefore, to prove that τ u is strictly increasing, it is sufficient to show that From relations where the latter follows from (28), (38), (40) and x 0 > x T , we see that in fact, relations (28), (38) and (40) yield that (52) is equivalent to which follows from x 0 > x T , ξ 2 = ξ * 2 /ξ 1 (see (24) and (28) This concludes the proof of i). ii) We have Using the previous result, to show that a solution of (46) exists, we only need to find a value of the energy h s ∈ J such that Lemma 3. There existsh s ∈ J and two functionsτ u (h s ),τ v (h s ) such that Proof. Using relation in the interval J, where the function ξ E η 1 is decreasing. Indeed, we have Moreover, we have ∆η 2 Then, using we get which follows from 0 < r ξ < 1 for h s ∈ J and Remark 5. From relation This concludes the proof of the existence of a family of brake periodic orbits x(t; s ) parametrised by s .

The Sun-shadow map
To study the Sun-shadow dynamics, it is useful to construct a Poincaré-like map. Traditionally, the Poincaré maps of autonomous two-dimensional Hamiltonian systems are defined by fixing the value of the Hamiltonian. This is not possible for the Sun-shadow dynamics, since the Hamiltonian is not conserved along the flow, see Proposition 2. However, here we can fix the value of L s . Indeed, even though L s is not a constant of motion as well, it assumes the same value in Stark's regime before and after the satellite crosses the shadow. To define a Poincaré map we also need to introduce a section. The selected section corresponds to the upper boundary of the shadow region in the (x, y) configuration space, i.e. to y = R, x ≥ 0. We consider trajectories leaving the section with p y > 0, in the outward direction with respect to the shadow region. Since the dependence on the coordinates is decoupled in the variables u, v (see (17), (18)), we decided to use them to define the map. Thus, we define the Sun-shadow map as where the domain D is discussed below, and (u, p u ), (u , p u ) belong to the section Σ defined as The conditions |u| ≥ √ R, uv = R are necessary to select the desired section in the (x, y) configuration space. The condition up v > −p u v is equivalent to p y > 0 (see (8)). The additional condition, up v > 0, assures that every point (u, p u ) ∈ Σ corresponds to only one trajectory. Indeed, there are points (u, p u ) for which p y > 0 in both the cases p v > 0 and p v < 0. Proposition 4. The map S is not defined in the points (u, p u ) with Moreover, in the second and fourth quadrant of the (u, it is not defined only in the points (u, p u ) with Proof. For each point (u, p u ) in the domain of the map, we have v = R u , and p v is defined by (18) 2 . Condition (54) corresponds to p 2 v ≤ 0, that is not possible. In the second and fourth quadrant of the (u, p u ) plane, the condition  Figures 11 and 12 respectively.
meaning that the map is not defined. On the other hand, if s < µ − f R 2 2 , the previous condition is fulfilled for The domain D ⊂ R 2 does not include the points defined in Proposition 4, nor the points corresponding to trajectories which go to infinity or collide with the Earth before going back to Σ. In Figure 9, for a specific choice of s and f , the domain D is drawn as the white area in a portion of the (u, p u ) plane. The light grey region represents the set D F of forbidden points in Proposition 4. The other two grey areas contain part of the sets D ∞ (darker) and D C (lighter) corresponding to the trajectories which go to infinity and collide with the Earth, respectively. Figure 10 shows the image of the domain D under S in the same portion of the (u, p u ) plane represented in Figure 9.
Remark 6. In the image of the map we may have points belonging to D ∞ or D C , so that we cannot iterate the map again.
In Figure 11, the magnification of the larger rectangular region appearing in Figure 9 highlights the complexity of the structure of D. We have selected five points, labelled with b, c, d, e, f, in the white 'corridors', and one point, labelled with a, in the larger white region on the left. In the same figure we show the portion of the trajectories corresponding Proof. Let Φ s (τ ; U , τ 0 ) and Φ k (τ ; U , τ 0 ) be the integral flow of Stark's and Kepler's dynamical systems (11), (12). Consider (u, p u ) ∈ D and the corresponding orbit in the Sun-shadow dynamics. Before it goes back to the section Σ, the dynamical regime will change n times, with n depending on the shape of the trajectory. The first regime will be always Stark's, the last will be Kepler's. Let us introduce a finite sequence of sections Σ i , i = 0, . . . , n where the dynamics changes, with Σ 0 = Σ n = Σ. Each of them is given by s i (U ) = 0: for the section Σ we have s i = uv − R, while for the intermediate sections, with i = 1, . . . , n − 1, we have s i = uv + R or s i = uv − R depending on the boundary of the shadow region that is crossed when the dynamics changes. It is possible to define the maps S i where D i is made of the points The Sun-shadow map will be given by Thus, we have The integral flow Φ i+1 (resp. X i+1 ) is equal to either Φ s or Φ k (resp. X s or X k ) depending on the regime between the two sections. The term ∂τ /∂U i can be computed as with I the identity matrix and τ i the value of the fictitious time at the section Σ i . Proposition 6. The map S is not area-preserving.
Proof. We give a numerical proof by showing that a circular region of Σ is mapped into a region with a different area. In the section Σ, we can consider a closed curve γ 0 symmetric with respect to the u axis, defined by We chose c = 1 s 2 /km 2 so that γ 0 becomes a circumference of radius r C centred at (u, p u ) = (u C , 0). We sample it with m points. Each point defines a trajectory in the phase space that we propagate with a Runge-Kutta method of Gauss type, by properly switching dynamics at the boundary between Stark's and Kepler's regimes, until its next intersection with the section Σ. In this way, we obtain the image of the initial points under the Sun-shadow map. The resulting points belong to the closed curve γ 1 = S(γ 0 ). To compute the area A 1 of the region enclosed by γ 1 , first we parametrise it by a variable θ. In particular, for each point on γ 1 we compute the values of the parameter θ: where (u j , p j u ) are the coordinates of the points, and (u m , p m u ) = (u 1 , p 1 u ). Then, we interpolate the points (θ j , u j ), (θ j , p j u ) by cubic splines. Finally, we compute the area by applying the Gauss-Green formula and we get The resulting area A 1 is different from the area A 0 = πr 2 C /c of the region enclosed by γ 0 . Indeed, assuming that s = 348600 km 3 /s 2 , f = 9.12 × 10 −9 km/s 2 , u C = 1250 km 1/2 , r C = 250 km 1/2 , and sampling the initial circumference with 2 × 10 5 points, in double precision we get A 0 ≈ 1.9635 × 10 5 km 2 /s, A 1 ≈ 1.9588 × 10 5 km 2 /s.

Hyperbolic fixed points
Given s ∈ [ − s , + s ], the periodic orbit of brake type x(t; s ) gives rise to two fixed points υ 1,2 of the Sun-shadow map: where ξ E is given by equation (40) and with h s the energy of x in Stark's regime, h s ∈ [h s , h * s ). The point υ 1 lies in the region included into the smaller rectangle appearing in Figure 9.

Invariant manifolds
Here, we describe the numerical technique used for the computation of the invariant manifolds of the fixed points of the Sun-shadow map. The discussion will be focused on υ 1 , but the procedure is the same also for υ 2 .
We took inspiration from the method in [5], thought specifically for planar maps. This algorithm can be applied only to two-dimensional maps which have saddle-type fixed points and whose Jacobian matrix, evaluated at these points, has two real eigenvalues λ 1 , λ 2 with 0 < λ 1 < 1 < λ 2 . As previously shown, the Sun-shadow map S and its fixed point υ 1 fulfil these requirements. We describe the algorithm for the case of one branch B of the unstable manifold. The stable manifold can be constructed in a similar manner using the inverse map S −1 . The method consists in dividing the branch of the manifold into a sequence of primary segments. A primary segment V is a connected subset of B whose last point is the image of its first point under the map. Given an initial primary V 0 in a neighbourhood of the fixed point, all the following primaries can be obtained by iterating the map m times: The branch will be given by the union of the computed primaries: The initial primary is approximated with a segment very close to υ 1 along an unstable eigenvector ψ of DS( υ 1 ). Then, it is corrected by using the technique described in [9], based on the Modified Fast Lyapunov Indicators (MFLI). The lower is the distance of a point from the manifold, the larger is its MFLI. Thus, for each point υ in the sample of the primary the correction is done as follows: 1. consider a small neighbourhood of υ in the direction orthogonal to the corresponding primary curve, and sample it uniformly; 2. compute the MFLI of υ and of each point of the sample; 3. select the point with the larger MFLI.
There is an issue concerning the iterations of the primaries. We observed that, after a few iterations, portions of the primaries are lost because the corresponding trajectories never return to Σ: by consequence the primaries lose their nature of connected sets. We decided to relax the definition of primaries given by Hobson by admitting primaries with several connected components, that we still denote by V i . We summarize below our algorithm: 1. approximate the initial primary V 0 with a segment aligned with the eigenvector ψ in a small neighbourhood of υ 1 : this segment is sampled with n points distributed according to the exponential law with υ i = (u i , p u i ) and a ∈ R. The last one, υ n−1 , is the image of υ 0 under the map. We callṼ 0 the finite set of points approximating V 0 ; 2. correctṼ 0 by using the MFLI, as previously described; 3. iterate the correctedṼ 0 once, and obtain the set  Figure 9. Here, s = 348600 km 3 /s 2 , f = 9.12 × 10 −9 km/s 2 .
6. iterate the correctedṼ 1 under the map m − 1 times. Figure 12 shows the first steps of the construction of a branch of the unstable manifold. The initial primary V 0 (top left) is first iterated once under the map (top right). Its image V 1 intersects the forbidden region D C (middle left), and at the second iteration of S three components are left (middle right), two of which lie in the {u < 0} half-plane. Again, their image V 2 intersects the forbidden regions D C , D ∞ (bottom left) and the number of components increases at the successive iteration (bottom right).
In Figure 13 we draw the four (disconnected) branches of the stable and unstable manifold of the fixed point υ 1 .

Conclusions and open questions
In this paper we have investigated the Sun-shadow dynamics, which is defined by patching together Kepler's and Stark's dynamics. After reviewing some relevant features of Stark's problem, we prove the existence of a family of periodic orbits of brake type. Then, we introduce the Sun-shadow map, by fixing a quantity which is not conserved along the flow. This map is differentiable but is not area preserving. Its domain shows fine structures that underlie interesting phenomena when the map is iterated many times. There is numerical evidence that the fixed points υ 1 , υ 2 related to the periodic orbits of brake type are hyperbolic; their invariant manifolds are constructed by an algorithm specifically created for this purpose. A global picture of this map is drawn in Figure 14, where an enhancement of the region close to the point υ 1 is also displayed. We observe evidence of regular and chaotic behaviour, with the presence of several islands: we checked that some of them surround periodic points. In the central region, where we have smaller values of |u|, the plotted points show a regular structure, similar to the one of the phase portrait in Figure 3. On the other hand, the regular behaviour seems to be lost in a neighbourhood of υ 1 , along the stable and unstable branches of its invariant manifold.
This study opens some interesting questions about the Sun-shadow dynamics, which deserve to be further investigated. First, we may wonder whether the winding number associated to the trajectories corresponding to one iteration of the map is bounded from below, see Figure 11. Another interesting question is whether we can show that the islands appearing in Figure 14 correspond to invariant curves around fixed or periodic points. Moreover, we can ask ourselves whether Melnikov's method can be adapted to prove the existence of chaotic dynamics in this case, where the invariant manifolds of the fixed points υ 1 , υ 2 are made of several connected components (maybe infinitely many) due to the presence of the forbidden regions D ∞ , D C , see Figure 12. Finally, possible future developments of this work are the extension of the Sun-shadow dynamics to the three-dimensional case, the inclusion of other perturbations (e.g. the Earth oblateness), and the study of the effect of the penumbra. Figure 14: On the top we show a global picture of the Sun-shadow map, with s = 348600 km 3 /s 2 , f = 9.12 × 10 −9 km/s 2 . At the bottom we zoom in the region close to υ 1 . Table 3: Stark's problem: boundary points in the ( s , h s / √ f ) plane. The features of the trajectories in the (x, y) plane are qualitatively described. i is the imaginary unit.