Theory for planetary exospheres: I. Radiation pressure effect on dynamical trajectories

The planetary exospheres are poorly known in their outer parts, since the neutral densities are low compared with the instruments detection capabilities. The exospheric models are thus often the main source of information at such high altitudes. We present a new way to take into account analytically the additional effect of the radiation pressure on planetary exospheres. In a series of papers, we present with an Hamiltonian approach the effect of the radiation pressure on dynamical trajectories, density profiles and escaping thermal flux. Our work is a generalization of the study by Bishop and Chamberlain (1989). In this first paper, we present the complete exact solutions of particles trajectories, which are not conics, under the influence of the solar radiation pressure. This problem was recently partly solved by Lantoine and Russell (2011) and completely by Biscani and Izzo (2014). We give here the full set of solutions, including solutions not previously derived, as well as simpler formulations for previously known cases and comparisons with recent works. The solutions given may also be applied to the classical Stark problem (Stark,1914): we thus provide here for the first time the complete set of solutions for this well-known effect in term of Jacobi elliptic functions.


Introduction
The exosphere is the upper layer of any planetary atmosphere: it is a quasi-collisionless medium where the particle trajectories are more dominated by gravity than by collisions.
Above the exobase, the lower limit of the exosphere, the Knudsen number (Ferziger and Kaper, 1972) becomes large, collisions become scarce, the distribution function cannot be considered as maxwellian anymore and, gradually, the trajectories of particles are essentially determined 1 INTRODUCTION 2 by the gravitation and radiation pressure by the Sun. The trajectories of particles, subject to the gravitational force, are completely solved with the equations of motion, but it is not the case with the radiation pressure (Bishop and Chamberlain, 1989).
In the absence of radiation pressure, we can distinguish three types of trajectories for the exospheric particles : -the escaping particles come from the exobase and have a positive mechanical energy: they can escape from the gravitational influence of the planet with a velocity larger than the escape velocity. These particles are responsible for the Jeans' escape (Jeans, 1916). They can also be defined as crossing only once the exobase, -the ballistic particles also come from the exobase but have a negative mechanical energy, they are gravitationally bound to the planet. They reach a maximum altitude and fall down on the exobase if they do not undergo collisions. They cross the exobase twice, -the satellite particles never cross the exobase. They also have a negative mechanical energy but their periapsis is above the exobase: they orbit along an entire ellipse around the planet without crossing the exobase. The satellite particles result in their major part from ballistic particles undergoing few collisions mainly near the exobase (Beth et al. (2014)). Thus, they do not exist in a collisionless model of the exosphere.
The radiation pressure disturbs the conics (ellipses or hyperbolas) described by the particles under the influence of gravity. The resonant scattering of solar photons leads to a total momentum transfer from the photon to the atom or molecule (Burns et al., 1979). In the nonrelativistic case, assuming an isotropic reemission of the solar photon, this one is absorbed in the Sun direction and scattered with the same probability in all directions. For a sufficient flux of photons in the absorption wavelength range, the reemission in average does not induce any momentum transfer from the atom/molecule to the photon. The momentum variation, each second, between before and after the scattering imparts a force, the radiation pressure. Bishop and Chamberlain (1989) proposed to analyze its effect on the structure of planetary exospheres. In particular, they highlighted analytically the "tail" phenomenon at Earth: the density for atomic Hydrogen, which is sensitive to the Lyman-α photons, is higher in the nightside direction than in the dayside direction in the Earth corona. Nevertheless, their work was limited only to the Sun-planet axis, with a null component assumed for the angular momentum around the Sun-planet axis. We thus generalize here their work to a full 3D calculation, in order to investigate the influence of the radiation pressure on the trajectories (this paper), as well as the density profiles and escape flux (following works).
This problem is similar to the so-called Stark effect (Stark, 1914): the effect of a constant electric field on the atomic Hydrogen's electron. Its study can be transposed to celestial mechanics in order to describe the orbits of artificial and natural satellites in the perturbed (e.g. by the radiation pressure force) Two-Body Problem. A recent description of the Stark effect solutions was already given by Lantoine and Russell (2011). However, they give the analytical solutions of all trajectories only in the 2D case (and for bounded trajectories in the 3D case), and their formulas have some issues as will be discussed later. Another analytical study proposed by Biscani and Izzo (2014) uses the Weierstrassian formulations to solve the motions for bounded and unbounded trajectories and to find periodic motions. Also, the motion can be approached numerically by developing the equations of motion in Taylor series but this leads to some issues for high eccentricities (Pellegrini et al., 2014). Hatten and Russell (2014) compared recently these methods and their computing efficiencies.
In this paper, based on the same formalism as Bishop and Chamberlain (1989), we provide for the first time the complete exact 3D solutions of the Stark effect (and its celestial mechanics analogue) for any initial condition and for both bounded and unbounded trajectories.
The first section describes the formalism used, before the sections 2/3/4 provide the equations of motion and time. We then discuss about circular orbits in section 6, while a comparison with previous works is given in section 6, before we conclude in section 7.

Model
In this work, we decide to study the effect of the radiation pressure on atomic Hydrogen in particular. Nevertheless, this formalism can be applied to any species subject to this force or to the interplanetary dust. We model the radiation pressure by a constant acceleration a coming from the Sun. According to Bishop (1991), this acceleration depends on the line center solar Lyman-α flux f 0 , in 10 11 photons.cm −2 .s −1 .Å −1 : In spherical coordinates, the Hamiltonian of one Hydrogen atom can be written: with r the distance from the planet, θ the solar angle, φ the angle with respect to the ecliptic plane, p r , p θ and p φ the conjugate momenta. −GMm/r represents the gravitational potential and mar cos θ the potential energy from the radiation pressure acceleration a. An example of trajectory of a H atom subject to the radiation pressure is given in the figure 1. This problem is similar to the classical Stark effect (Stark, 1914): a constant electric field (here the radiation pressure) is applied to an electron (here an Hydrogen atom) attached to a proton (here the planet). Both systems are equivalent because the force applied by the proton (the planet) to the electron (the Hydrogen atom), i.e. the electrostatic force, varies in r −2 as the gravitational force from the planet on the Hydrogen atom. Thus, we adopt the same formalism as Sommerfeld (1934) and use the parabolic coordinates. We use the transformation: with u and w always positive. Consequently, the Hamiltonian becomes: independent from t and φ.
According to hamiltonian canonical relations, we have:

Constants of the motion
In this new system of coordinates, we study this Hamiltonian. First, H is independent from t explicitly. H is a conserved quantity along the time and corresponds to the mechanical energy E of the system. Moreover, as H is independent from φ, according to canonical relations: Thus, p φ is another constant of the motion. Once E and p φ defined, the equation 4 can be rewritten: The left hand side is a function dependent only on u and p u , the right hand side depends only on w and p w . As both functions are equal and independent, they are equal to a constant A, a separation constant: The motion possesses three constants: E, A and p φ . The equation 8 allows to express p u (respectively p w ) as a functions of E, A, p φ and u (respectively w):

Effective potentials
We have already introduced the Hamiltonian H of the system. We can extend the approach according to Hamilton-Jacobi equations: where S is the Hamilton's principal function or action. This function depends on initial conditions (as u 0 , w 0 , φ 0 and t 0 ) and the actual position of the particle (as u, w, φ and t). As previously demonstrated, H = E and p φ are constants. Thus, and leads to: withŜ the part of the action independent from t 0 , t, φ 0 and φ.
Moreover, the actionŜ can be separated into two parts: one written as a function of u and p u coordinates, the other one with w and p w . By definition, according to the Hamilton-Jacobi equations, we have: with p u (resp. p w ) a function only of u (resp. w), assuming E, A and p φ values already fixed by initial conditions. Then, we can separate again the action, leading to: According to the equation 9, we have the following relations: V u and V w are effective potentials applied in u and w directions (represented in the figure 2).
These potentials play key roles for the motion because they constrained the motion in u and w directions independently. For the motion of the particle, we must respect two conditions: E > V u (u) and E > V w (w). These conditions are analogous to: These both conditions are more restrictive than the usual E > E p where E p is the potential energy. The different roots of P 3 and Q 3 are displayed and correspond to the intersection of the potentials with the horizontal black line. U − is the dimensionless value referring to u − (cf. section 2.7, table 1). Notice that V W will cross only once the energy level E for too high or too low E values.

Study of P 3
P 3 is a polynom of degree 3 with lim u→+∞ P 3 (u) = +∞. This polynom possesses three roots, whose one is real at least. As P 3 (0) = p 2 φ > 0, one of these roots is real negative, according to intermediate value theorem since lim u→−∞ P 3 (u) = −∞. Nevertheless, the motion occurs for positive u values, and we know this motion exists. It implies there is an interval in R + such as P 3 < 0 (otherwise, there is no motion in u-direction, no possible physically, cf. eq. 19). To comply with this last condition, both other roots are real and positive. In summary, P 3 has three real roots: one negative and two positive.
We call each root u 0 , u − and u + such as u 0 < 0 < u − < u + and the u-motion is restricted in term of dimensionless quantities, as can be seen in the figure   2).

Study of Q 3
Q 3 is a polynom of degree 3 with lim w→+∞ Q 3 (w) = +∞. This polynom possesses three roots, whose one is real at least. As Q 3 (0) = −p 2 φ < 0, one of these roots is real positive, according to intermediate value theorem since lim w→+∞ Q 3 (w) = +∞. Nevertheless, the motion occurs for positive w values. We have restrictions on both other roots: they must be both real positive, both real negative or both complex conjugates.
In the case where the three roots are real positive, we call each root w 0 , w − and w + such as 0 < w − < w + < w 0 and the motion is restricted to w ∈ [w − ; w + ] ∪ [w 0 ; +∞[ (as can be seen in the figure 2 with the dimensionless quantities, cf. section 2.6).
In the case with one positive root and both other complex or real negative, we call each root w 0 (the positive one), w − and w + (keep the same order as previously defined if they are real) such as the motion is restricted to w ∈ [w 0 ; +∞[ only ([W 0 ; +∞[ in term of dimensionless quantities).

Restriction on the motion
Each constant value of u or w defines a paraboloid in three dimensions. For each interval, constrained by fixed values of u and w, the motion will be contained between the paraboloids defined by these limit values as shown in the figure 3. For the u-motion, this is always limited by two paraboloids defined by u = u − and u = u + as shown by the blue area in the figure 3 (left panel). Similarly, for the w-motion, there are two cases: the motion is constrained between one paraboloid (w = w 0 ) and the infinity or between two paraboloids (w = w − and w = w + ). Both cases are represented by the red area in the figure 3 (right panel).

Summary of restrictions
As previously demonstrated, the motion is constrained in specific areas of the 3D space. The motion is always between two paraboloids due to restrictions on u but it can be constrained between two other paraboloids or only one regarding w. Thus, the motion is constrained by four paraboloids or three, opened to infinity as shown on the figure 4 by the green area. Nevertheless, there is at this point no other information on the exact motion of the particle.
As shown in the figure 4 (right panel), the particle seems to explore all the area when its motion is restricted by four paraboloids. This is simply an observation. How can we prove that without any information on the exact motion of the particle? According to the Poincaré recurrence theorem, if the dynamical trajectory of the autonomous system evolves in a finite volume of the phase space, then in any domain, as small as it could be, there are at least two points which belong to the same trajectory. Here, the motion is constrained in space with the four paraboloids but also in velocity because p u and p w are finite values and p φ is constant. All the positions in this part of the phase space can belong to the same trajectory. This is also linked to the KAM theory: here, the perturbation (the radiation pressure) affects the periodic motion (the ellipse, bounded trajectories) but it can remain quasi-periodic. As we will see, the global motion is not periodic but u and w motions possess their own period according to another parameter. The global motion can be periodic only if all periods are commensurable. This is an important result because if the particle belongs to an area crossing the exobase and if this area is closed in the phase space, along a finite time, the particle will again cross the exobase. We can now define the ballistic and satellite particles as follows : both populations have no elliptic trajectories due to the radiation pressure, but they evolve in a closed domain.
Depending on their constants of the motion, we can easily determine whether they cross (i.e. the domain crosses) the exobase or not, corresponding to ballistic and satellite particles respectively.
Escaping particles are in the case where the initial value of w is higher than the highest real root of Q 3 and their available area is opened to the infinity. Thus, the theorem cannot be applied here. Even if the restriction area crosses the exobase, the particle may come from infinity, come close to the exobase and go away without crossing the exobase. We need, in order to identify escaping particles (that cross the exobase and go to the infinity), to track along the time the particle to know if these particles come from the exobase or not. Thus, it is necessary to solve their trajectory along the time.

Dimensionless expressions
For convenience, as usual in fluid mechanics, we define characteristic quantities. We decide to write all equations with dimensionless parameters. First, for distance, we define: This characteristic value was introduced by Bishop (1991) and defines the limit distance where the radiation pressure overwhelms the gravitation of the planet. Then, we rewrite: where U and W are the dimensionless quantities associated to u and w. The energy E is adimensionned with the use of the characteristic energy k B T exo . For p u and p w , we express choose mk B T exo GM/a as the unit for A and the unit for the adimensionless time τ is derived from the units of u/w and p u /p w . We summarize these changes in table 1.
We express all the previous equations as a function of these new quantities: with λ a : the Jeans parameter at the distance R pressure .
We can also reduce the equations of the motion. We introduce the dimensionless time τ .

Dynamical trajectories
In this part, we give implicit expressions for the dynamical trajectories of the particles, under the influence of both gravity and radiation pressure. Such expressions were already given for the 2D case with P φ = 0 and partially generalized to the 3D case in Lantoine and Russell (2011) with Jacobi elliptic functions (only bounded) (Jacobi, 1829) and in Biscani and Izzo (2014) with the Weierstrass functions for bounded and unbounded trajectories. This last work dealt with the dynamical trajectories of artificial satellites but it can be applied to exospheric species subject to the radiation pressure. We propose here corrections as well as a better way to give "simple" expressions for dynamical trajectories. We also provide expressions for unbounded trajectories that are missing in the literature. In the same way, we introduce our notations for the next papers to be published, where the influence of the radiation pressure on the density profiles and escape flux will be investigated.
According to the previous part, we have different restrictions on the motion and thus, we must distinguish the cases, that will correspond to the different types of possible trajectories.
We may thus define the ballistic/satellite/escaping populations based on the roots of the P 3 and Q 3 polynoms. The trajectories are not elliptic or hyperbolic at all but one can keep their basic definitions: crossing twice the exobase for ballistic particles, orbiting but not crossing the exobase for satellite particles, coming from the exobase and escaping to the infinity for escaping particles.
The number of real roots is a key parameter for the analytical resolution of the trajectories.
The previous equations of the motion, expressed as a function of time, must be rewritten as a function of a new variable T 0 (the subscript 0 is necessary to distinguish from τ ) defined as: Until we know the expression of U(T 0 ) and W (T 0 ), we cannot express solutions as a function of time. The system of equations leads now to: All solutions, trajectories and time, are parametrized according to the variable T 0 without physical sense.

The U-motion
We solve the differential equation describing the u-motion: with U 0 , U − and U + being real roots.
First, we set Y 2 = U − U 0 . The motion occurs always with U > 0 and U 0 < 0 to justify this change. We obtain: Finally, we have: We define k U as: The final equation is: The solution of this equation is: dn is a Jacobi elliptic function and T U depends on initial conditions. The final expression for U is: with cn and sn other Jacobi elliptic functions. These functions are 4K(k U )-periodic (see Appendix, eq. A.1). To determine T U , we preferably choose the second expression for the sign.
The initial conditions give us U(0), then: where arcsn is defined as the reciprocal function of sn: with F the incomplete elliptic integral of the first kind (see Appendix, eq. A.1).
The sign depends on P U : If P U is positive (resp. negative), U increases (resp. decreases) and −sn 2 increases (resp. decreases) too. Then, T U is positive (resp. negative). This solution corresponds to the solution η given by Lantoine and Russell (2011). The U-motion is always constrained. The characteristics of the motion/particle (ballistic, satellite or escaping) is thus determined by conditions on Q 3 (W ).

The W-motion
For the w-motion, we need to solve the differential equation: W 0 is a positive real root but we must distinguish between various cases for the other roots: W + and W − may be real or complex conjugate roots.

Three real roots
We must consider two cases: W > W 0 and W + > W > W − . The second one occurs mathematically and not only physically when W − and W + are positive.
a) For W − < W < W + . This case occurs for bounded trajectories (ballistic or satellite trajectory). Here, this is the same treatment as for the U-motion. After setting the two following , we obtain a similar solution: The final expression for W is The initial conditions giving us W (0) and then leads to: The sign depends on P W : If P W is positive (resp. negative), W increases (resp. decreases) and sn 2 increases (resp. decreases) too. Then, T W is negative (resp. positive). This solution corresponds to the solution ξI given by Lantoine and Russell (2011).
b) For W > W 0 . First, the motion occurs only for W ∈ [W 0 ; +∞[ (corresponding to escaping particles). Then, we set Y 2 = W − W 0 . We obtain: Finally, we have: We define k W as: The solution of this equation is cs is a Jacobi elliptic function and T W depends on the initial conditions.
The final expression for W is: These functions are 4K(k W )-periodic and are defined on T 0 − T W ∈ R/{4mK(k W )|m ∈ Z}.
These functions diverge at 4mK(k W ) but this is not an issue: the motion can diverge with respect to T 0 , but, according to the time τ and the integration (see Section 4), the w-motion remains continuous for τ ∈ R as Let us assume the initial conditions provide W (0), then: ±arccs The sign depends on P W : If P W is positive (resp. negative), W increases (resp. decreases) and cs 2 increases (resp. decreases) too. Then, T W is positive (resp. negative). This solution corresponds to the solution ξII given by Lantoine and Russell (2011). These expressions are useful only for three real roots.
A last case remains: only one real positive root.

Only one real positive root
For any initial conditions, the particle will be escaping if Q 3 has only one real root. The solution for this case is more complex compared with the previous solutions. As done by Lantoine and Russell (2011), we could apply some transformations and obtain a new expression, that would be a combination of cn and sn. Nevertheless, we propose a simpler way to determine the expression with the knowledge of all roots, even imaginaries. We start from the equation 56 with W + and W − not real: After the separation of variables, we obtain: Now, we apply the procedure proposed in Abramowitz and Stegun (1964) (p. 597) in the case where we have only one real root. First, we define: and also According to Abramowitz and Stegun (1964) (p. 597), the left hand side corresponds to: with F the Elliptic function of first kind (defined in the appendix Appendix A) and According to the equation 57: with am the Jacobi amplitude defined as Thus This function is 4K(k W 0 )-periodic and is defined on Z}. This function diverges at 2K(k W 0 ) + 4mK(k W 0 ) but this does not impact our result: the motion can diverge with respect to T 0 , but, according to the time τ and the integration (see Section 4), the w-motion remains continuous for τ ∈ R as As usual, according to the initial conditions, we define T W 0 as If P W is positive (resp. negative) then T W 0 is negative (resp. positive). This solution corresponds to the solution η given by Lantoine and Russell (2011) but our derived solution is less complex to use (with less time computing).

Time equation
We gave analytical formulations for the different kinds of trajectories, expressed implicitly. Now, since we have all expressions of the trajectories as a function of T 0 , we can express the real time τ (or t): Here, we use MATHEMATICA and MAPLE to derive the primitives. It is necessary to be very careful since these different programs can have different definitions of the elliptic functions for example. To avoid these issues, we remind at each use the definition employed. The first part of the integral gives: with α = √ λ a √ U + − U 0 , E the incomplete elliptic function of the second kind (see Appendix, eq. A.1).
The second part of the integral is more complex because we have different expressions according to the number of roots and the initial conditions. In the case of three real roots, if If W (0) > W 0 with three real roots then Finally, in the case of only one real root, the time equation is given by: For both last cases, T 0 − T W belongs to ]0; 4K(k W )[ (respectively ] − 2K(k W 0 ); 2K(k W 0 )[) for the equation 71 (resp. for the equation 72). Nevertheless, when T 0 − T W tends to 4K(k W 0 ) (resp. 2K(k W 0 )), the integration diverges and τ too. The w-motion occurs on an subset of R with respect to T 0 but on R entirely with respect to τ . The transformation of τ into T 0 is bijective because the integrand U + W is strictly positive.

φ-equation
To complete the description of the motion as a function of time, it is also necessary to solve the evolution of the angle φ, obeying to: As already done in the previous part, we separate into two integrals. The first part still gives: with Π the incomplete elliptic integral of the third kind (see Appendix, eq. A.1).
In the three real roots case for Q 3 , if initially we have W + > W (0) > W − : or if initially W (0) > W 0 : and for the last case with one real root: This last formula is only available for continuous.

Circular orbits
With the solutions previously derived, we can know the exact motion of a bounded or unbounded particle as a function of the time such as given in figure 6. It is clear that even a bounded trajectory the motion has no periodicity at all (especially for the φ motion). Nevertheless, it could be interesting to focus on stable bounded orbits and search for periodic motions (as in Biscani and Izzo (2014)), and thus investigate in particular the circular stable orbits for spacecrafts (Namouni and Guzzo, 2007) or the possible positions for satellite particles produced by collisions in the exosphere (Beth et al., 2014). Thus, we dedicate this section to the conditions to obtain such orbits.
For a specific set of initial conditions, it can be possible to obtain circular orbits. This orbit occurs when, on the one hand, the attraction of the planet projected along the x-axis is equal to acceleration due to the radiation pressure: In dimensionless quantity, this will be expressed by: On the other hand, it is also necessary that the centrifugal force induced by the rotation around the x-axis is equal to the acceleration around the planet in the perpendicular plane to the x-axis. Thus, we obtain the secondary equality: In dimensionless quantity, this will be expressed by: Combining these two equations, we obtain: We need to study the polynom with X = sin 2 θ ∈ [0; 1]. Depending on the P 2 φ values, we have zero, one or two solutions for P (X) = 0. Indeed, P (0) = P (1) > 0 so that according to the Rolle theorem, there is a ∈]0; 1[ with P ′ (a) = 0. Here, a is equal to 8/9. Nevertheless, if P (a) is positive then P (X) = 0 does not have solutions and inversely, if P (a) is negative then we have two solutions. This value is: A critical maximum value of P φ thus exists to allow for circular orbits and is Above this value, we cannot find any bounded trajectories: there is no any equilibrium point and thus no circular orbits (stable or not). For lower values of |P φ |, we have two solutions for P (X) = 0: one stable and one unstable as shown by Namouni and Guzzo (2007) and theirs plots of the equipotentials. These two solutions correspond respectively to the stable point around which the equipotentials are closed and to the saddle point, which is the last limit where one can find closed equipotentials, and is the only point where two equipotentials can cross. As long as |P φ | < |P cφ |, these both specific points exist: they have the same U. Physically, the potential V W has two extrema as plotted in figure 2. When |P φ | reaches |P cφ |, the local minimum goes to the right and the local maximum goes to the left at the same location W crit . For higher |P φ | values, V W have no extremum anymore: the potential is strictly decreasing and the particles are unbounded (escaping). Thus, the bounded particles, satellite and ballistic particles, have The distinction between them thus depend on if they cross or not the exobase.
The critical values are given in (z, ρ) coordinates (z is −x for us, in comparison with Namouni and Guzzo (2007)). In dimensionless unities and in (R, θ) using the equations 81 then 79, the critical orbit is: or in (U, W ) coordinates: The real positive roots of the polynomial 83 combined with the equality 79 give the positions of the circular orbits (two coordinates are necessary) allowed to spacecraft or particles under the influence of both gravity and radiation pressure.

Summary
The knowledge of the exact trajectories of particles or satellites under the influence of gravity and radiation pressure needs the calculation of the spatial coordinates, i.e. the U/W/φ motions, as well as the time evolution. We summarize all needed equations in table 7. Besides, our 3D solutions can be easily applied to the 2D case. Indeed, in the 2D case, P φ = 0 and thus, one of the roots for each polynomial P 3 and Q 3 is null: it could be U 0 or U − for P 3 (if U + = 0, there is no possible motion) and any roots of Q 3 . We precise that in this case the φ-motion is not important because the motion is planar. Compared with Lantoine and Russell (2011), our formulations are first developed for the 3D case and can be used easily for the 2D  Figure 5: Summary of the different solutions for each kind of motion. The notation *NEW* corresponds to the new solutions derived in this paper, not yet derived in previous works. The symbol (s) corresponds to an example of formula with a simpler expression than the one proposed by Lantoine and Russell (2011). L2D labeled solutions were explicitely given by Lantoine and Russell (2011) only in the 2D case, whereas L3D labeled solutions were also given by Lantoine and Russell (2011) in the 3D case. The function sign is noted sg.  Lantoine and Russell (2011) gave only the methodology to obtain the 3D solutions based on 2D ones but not the expressions, which apparently leads to complex expressions. Biscani and Izzo (2014) provided also the exact formulas for bounded and unbounded trajectories using the Weierstrass functions but this formulation is also difficult particularly because of the need to use the Inverse Weierstrass function, not implemented in all computer softwares and the need to work with complex values (e.g. the complex logarithm function). In this paper, we solved the motion for the bounded and unbounded trajectories in the 3D case; we provide the exact formulas for all cases, as well as the definitions used in Appendix A. We also highlight in table 7 which solutions are simpler compared with Lantoine and Russell (2011), which are completely new and which were only provided in the 2D case by Lantoine and Russell (2011).  Biscani and Izzo (2014) and Lantoine and Russell (2011). Even if we do not agree with the number of evaluations of each Jacobi elliptic function mentioned by Hatten and Russell (2014) (e.g. to call cn and sn is similar to call am then cos and sin: cn = cos • am and sn = sin • am), two arguments show our solutions are efficient in terms of CPU time and accuracy : first, the solutions expressed in terms of Jacobi elliptic functions (such as in this paper or by Lantoine and Russell (2011)) are more efficient than Weierstrass elliptic functions (used by Biscani and Izzo (2014)) ; second, several solutions given above are less complex to implement than those by Lantoine and Russell (2011), e.g. equation 65. Also, the analytical formulations are preferable for long duration motions.

Conclusions
We determined analytically the trajectories of the particles or spacecraft under the influence of both planetary gravity and stellar radiation pressure. We thus provide for the first time the complete exact solutions of the well-known Stark effect (effect of a constant electric field on the atomic Hydrogen's electron) with Jacobi elliptic functions, for both bounded and unbounded orbits. These expressions may be implemented for modeling spacecraft or particles trajectories: instead of solving the equation of the motion, based on differential equations, with numerical methods such as the Runge-Kutta method where one cumulates errors along the time, it is here possible to obtain precise expressions of the motion with only periodic errors, due to the precision on the evaluation of the elliptic functions used. In particular, we provide the analytical conditions for stable circular orbits. Moreover, we discuss about the possible issues inherent to the formalism used and the importance of being extremely careful with the routines implemented.
The formalism used here will allow us in a next paper to generalize the work by Bishop and Chamberlain (1989) to derive the exact neutral densities and escape flux in planetary exospheres, under the influence of both gravity and stellar radiation pressure. This is important for understanding the atmospheric structure and escape of planets in the inner solar system, as well as the atmospheric erosion during the early ages where the radiation pressure (and UV flux) of the Sun was extreme.

Acknowledgments
This work was supported by the Centre National d'Études Spatiales (CNES).

Appendix A. Elliptic integrals
In this paper, we use the three incomplete elliptic integrals F , E and Π: (1 − n sin 2 θ) 1 − k 2 sin 2 θ dθ (A.1) Sometimes, other formulas (shown below) are proposed with the change sin θ = t but one