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

A. Beth, P. Garnier, D. Toublanc, I. Dandouras, C. Mazelle

► To cite this version:

A. Beth, P. Garnier, D. Toublanc, I. Dandouras, C. Mazelle. Theory for planetary exospheres: I. Radiation pressure effect on dynamical trajectories. Icarus, 2016, 266, pp.410-422. 10.1016/j.icarus.2015.10.018. insu-03670536

HAL Id: insu-03670536

https://insu.hal.science/insu-03670536

Submitted on 17 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

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

<figure><img src="image_4.png" /></figure>

CrossMark

A. Beth<sup>a,b,c,*</sup>, P. Garnier<sup>a,b,1</sup>, D. Toublanc<sup>a,b</sup>, I. Dandouras<sup>a,b</sup>, C. Mazelle<sup>a,b</sup>

## ARTICLE INFO

**Article history:**
Received 5 February 2015
Revised 6 August 2015
Accepted 19 October 2015
Available online 28 October 2015

**Keywords:**
Atmospheres, structure
Solar radiation
Terrestrial planets

## ABSTRACT

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 generalisation of the study by Bishop and Chamberlain (Bishop, J., Chamberlian, J.W. [1989]. Icarus 81, 145–163). In this first paper, we present the complete solutions of particles trajectories, which are not conics, under the influence of the solar radiation pressure with some assumptions. This problem is similar to the classical Stark problem (Stark, J. [1914]. Ann. Phys. 348, 965–982). This problem was largely tackled in the literature and more specifically, recently by Lantoine and Russell (Lantoine, G., Russell, R.P. [2011]. Celest. Mech. Dynam. Astron. 109, 333–366) and by Biscani and Izzo (Biscani, F., Izzo, D. [2014]. Mon. Not. R. Astron. Soc. 439, 810–822) as we will discuss in this paper. We give here the full set of solutions for the motion of a particle (in our case for an atom or a molecule), i.e. the space coordinates and the time solution for bounded and unbounded trajectories in terms of Jacobi elliptic functions. We thus provide here the complete set of solutions for this so-call Stark effect (Stark, J. [1914]. Ann. Phys. 348, 965–982) in terms of Jacobi elliptic functions (Jacobi, C.G. J. [1829]. Fundamenta nova theoriae functionum ellipticarum. Sumtibus fratrum), which may be used to model the trajectories of particles in planetary exospheres.

© 2015 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

## 1. 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 any more and, gradually, the trajectories of particles are essentially determined 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, the exospheric particles can be distinguished into three categories according to their trajectories:

* the escaping particles come from the exobase and have a positive mechanical energy such as escape from the gravitational well of the planet because their velocity is larger than the escape velocity. These particles are responsible for the Jeans' escape (Jeans, 1916) or thermal escape: supposing a Maxwellian distribution for atoms or molecules at the exobase, it always exists a piece of its distribution exceeding the escape velocity. Theirs trajectories describe an hyperbola crossing only once the exobase,

* the ballistic particles come from the exobase as well but have a negative mechanical energy. They are gravitationally bounded to the planet. They reach a maximum altitude and fall down

<footnote><sup>a</sup>Université de Toulouse, UPS-OMP, IRAP, Toulouse, France</footnote>

<footnote><sup>b</sup>CNRS; IRAP, 9 Av. Colonel Roche, BP 44346, F-31028 Toulouse Cedex 4, France</footnote>

<footnote><sup>c</sup>Department of Physics/SPAT, Imperial College London, United Kingdom</footnote>

<footnote>* Principal corresponding author at: Department of Physics/SPAT, Imperial College London, United Kingdom.</footnote>

<footnote n="1">¹ Corresponding author at: Université de Toulouse; UPS-OMP; IRAP, Toulouse, France.</footnote>

<footnote>E-mail addresses: arnaud.beth@gmail.com (A. Beth), pgarnier@irap.omp.eu (P. Garnier).</footnote>
on the exobase if they do not undergo collisions. Their trajectories describe an ellipse crossing twice the exobase: their periapsis is below the exobase but their apoapsis is above,

* the satellite particles have elliptic trajectories as well but they 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 cannot exist in collisionless models 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 non-relativistic 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 analyse 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 generalise here their work to a full 3D calculation, in order to investigate the influence of the radiation pressure on the trajectories (this paper) by presenting results from Celestial Mechanics in order to be transposed to this study, as well as the density profiles, the escape flux and the planetary atmosphere stability (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 with the assumption that the Sun is fixed in the planetary frame (i.e. no stellar gravity, no centrifugal force and no Coriolis force). A recent description of the Stark effect solutions was already given by Lantoine and Russell (2014) in terms of Jacobi elliptic functions, by Biscani and Izzo (2014) in terms of Weierstrassian formulations and by Pellegrini et al. (2014) in terms of Taylor series.

Lantoine and Russell (2011) provided the analytical solutions of all trajectories in the planar case (2D coordinates and time solution). Moreover, they proposed the necessary transformations to pass from the 2D case to the 3D one. However, they omitted the analytical solution for the motion around the Sun-planet axis for unbounded trajectories, even if the methodology is given. The knowledge of this solution is essential to us: in order to estimate the escaping density and flux from the exobase or take into account the asymmetries of the exobase for future works. Indeed, at a given position in space, the density of escaping or unbounded particles depends on the number of particles (1) coming from the exobase and (2) passing through the point of interest. Thus, we shall identify which trajectory fills both (1) and (2) conditions. Finally, the Jacobi elliptic functions are often defined by two ways (see Appendix A) in the literature and computational software. One should be careful about the definition used. In this paper, we use the second definition and explain our choice (see Appendix A).

The study proposed by Biscani and Izzo (2014) used the Weierstrassian formulations to solve the motions for bounded and unbounded trajectories and to find periodic motions. This formulation is pretty inconvenient for numerical transpositions: the Weierstrass function, its derivative and its inverse are defined in the complex plane. This leads to time consuming operations to describe the trajectory of a particle according to the initial conditions.

On another hand, the motion can be approached numerically by developing the equations of motion in Taylor series (Pellegrini et al., 2014) but this leads to some issues for high eccentricities (Hatten and Russell, 2015): as pointed out, some mean anomaly values may cause the divergence of the series. By the way, Hatten and Russell (2015) compared recently the three methods and their computing efficiencies as well.

In this paper, based on the same formalism as Bishop and Chamberlain (1989), we provide the keys to introduce our future works whose purpose is to describe the influence of the radiation pressure on planetary exospheres by a Hamiltonian way. In this first paper, we provide the complete exact 3D and time solutions of the Stark effect (and its celestial mechanics analogue) for any initial condition and for both bounded and unbounded trajectories, whose some of them are not explicitly given in the literature up to now and simplify some formulations in order to describe the motion of the exospheric atoms and molecules.

Section 2 describes the formalism used, before the Sections 3/4/5 provide the equations of motion and time. We then discuss about circular orbits in Section 6, whilst a comparison with previous works is given in Section 7, before we conclude in Section 8.

## 2. 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 centre solar Lyman-α flux $f_0$, in $10^{11}$ photons cm$^{-2}$ s$^{-1}$ Å$^{-1}$:

$$a = 0.1774 f_0 (\text{cm s}^{-2}) \qquad (1)$$

In spherical coordinates, the Hamiltonian of one Hydrogen atom can be written in the spherical coordinate system:

$$\mathcal{H}(r, \theta, \phi, p_r, p_\theta, p_\phi, t) = \frac{p_r^2}{2m} + \frac{p_\theta^2}{2mr^2} + \frac{p_\phi^2}{2mr^2 \sin^2 \theta} - \frac{GMm}{r} + m \cos \theta \qquad (2)$$

with $r$ the distance from the planet, $\theta$ the angle between the zenith and the Sun direction, $\phi$ the angle with respect to the ecliptic plane, $p_r$, $p_\theta$ and $p_\phi$ the conjugate momenta. $-GMm/r$ represents the gravitational potential and $m \cos \theta$ 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 Fig. 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:

$$\begin{aligned} u &= r + x = r(1 + \cos \theta) \\ w &= r - x = r(1 - \cos \theta) \end{aligned} \qquad (3)$$
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Example of the trajectory of a Hydrogen atom in the Earth (central sphere) exosphere assuming that the Sun is placed at the infinity and so is fixed in the planetary frame. Left panel: view from the Sun. Middle panel: view from the side. Right panel: view from the top. The position of the Sun is given by the opposite direction of the arrow. The initial position is provided by the black cross and the final one by the blue cross. Without the radiation pressure, the trajectory must be the red ellipse. The trajectory with radiation pressure is clearly far from the one without. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

with $u$ and $w$ always positive. Consequently, the Hamiltonian becomes:

$$ \mathcal{H}(u, w, p_u, p_w, p_\phi) = \frac{2up_u^2 + 2wp_w^2}{m(u+w)} + \frac{p_\phi^2}{2muw} - \frac{2GMm}{u+w} + ma\frac{u-w}{2} \quad (4) $$

independent from $t$ and $\phi$.

According to hamiltonian canonical relations, we have:

$$ \begin{aligned} p_u &= \frac{m(u+w)}{4u} \frac{du}{dt} = \frac{p_r}{2} - \sqrt{\frac{w}{u}} \frac{p_\theta}{u+w} \\ p_w &= \frac{m(u+w)}{4w} \frac{dw}{dt} = \frac{p_r}{2} + \sqrt{\frac{u}{w}} \frac{p_\theta}{w} \\ p_\phi &= muw \frac{d\phi}{dt} \end{aligned} \quad (5) $$

## 2.1. Constants of the motion

In this new system of coordinates, we study this Hamiltonian. First, $\mathcal{H}$ is independent from $t$ explicitly. $\mathcal{H}$ is a conserved quantity along the time. In this case study, the gravity and the radiation pressure are conservative forces, i.e. they depend only on the position of the particle. Thus, the mechanical energy $E$ is conserved and corresponds to the Hamiltonian $\mathcal{H}$. Moreover, as $\mathcal{H}$ is independent from $\phi$, according to canonical relations:

$$ \frac{dp_{\phi}}{dt} = -\frac{d\mathcal{H}}{d\phi} = 0 \quad (6) $$

Thus, $p_\phi$ is another constant of the motion. This corresponds to the component of the angular momentum along the $x$-axis. Indeed, the gravitation does not affect the evolution of this one because this is a central force and the radiation pressure acts only along the $x$-axis and thus do not affect either in this direction. Once $E$ and $p_\phi$ defined, the Eq. (4) can be rewritten:

$$ \begin{aligned} 2muE - 4up_u^2 - \frac{p_\phi^2}{u} - m^2au^2 + 2GMm^2 \\ = -2mwE + 4wp_w^2 + \frac{p_\phi^2}{w} - m^2aw^2 - 2GMm^2 \end{aligned} \quad (7) $$

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:

$$ \begin{aligned} A &= 2muE - 4up_u^2 - \frac{p_\phi^2}{u} - m^2au^2 + 2GMm^2 \\ &= -2mwE + 4wp_w^2 + \frac{p_\phi^2}{w} - m^2aw^2 - 2GMm^2 \end{aligned} \quad (8) $$

$A$ is similar to the norm of the Laplace–Runge–Lenz vector (Redmond, 1964) which is proportional to the eccentricity, constant in the Keplerian problem.

The motion possesses three constants: $E, A$ and $p_\phi$. The Eq. (8) allows to express $p_u$ (respectively $p_w$) as a functions of $E, A, p_\phi$ and $u$ (respectively $w$):

$$ \begin{aligned} p_u &= \pm \sqrt{\frac{-P_3(u)}{4u^2}} \\ p_w &= \pm \sqrt{\frac{Q_3(w)}{4w^2}} \end{aligned} \quad (9) $$

with (see Bishop and Chamberlain (1989) Eq. (6) or Lantoinne and Russell (2011), Eq. (89a) with $P_3(u^2) = -P_\eta(\eta)$ and Eq. (89b) with $Q_3(w^2) = P_\phi(\xi)$)

$$ \begin{aligned} P_3(u) &= mau^3 - 2mEu^2 - (2GMm^2 - A)u + p_\phi^2 \\ Q_3(w) &= maw^3 + 2mEw^2 + (2GMm^2 + A)w - p_\phi^2 \end{aligned} \quad (10) $$

## 2.2. Effective potentials

We have already introduced the Hamiltonian $\mathcal{H}$ of the system. We can extend the approach according to Hamilton–Jacobi equations:

$$ \begin{aligned} \frac{\partial S}{\partial q_i} &= p_i \\ \frac{\partial S}{\partial t} &= -\mathcal{H} \end{aligned} \quad (11) $$

where $S$ is the Hamilton’s principal function or action. This function depends on initial conditions (as $u_0, w_0, \phi_0$ and $t_0$) and the actual
position of the particle (as u, w, φ and t). As previously demonstrated, $\mathcal{H} = E$ and $p_\phi$ are constants. Thus,

$$ 
\begin{aligned}
\frac{\partial S}{\partial \phi} &= p_{\phi} \\
\frac{\partial S}{\partial t} &= -E
\end{aligned}
\qquad (12) $$

and leads to:

$$ S = -E(t - t_0) + p_{\phi}(\phi - \phi_0) + \hat{S}[u_0, w_0, u, w, E, p_{\phi}] \qquad (13) $$

with $\hat{S}$ the part of the action independent from $t_0$, $t$, $\phi_0$ and $\phi$.

Moreover, the action $\hat{S}$ 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:

$$ 
\begin{aligned}
\frac{\partial S}{\partial u} &= \frac{\partial \hat{S}}{\partial u} = p_u \\
\frac{\partial S}{\partial t} &= \frac{\partial \hat{S}}{\partial w} = p_w
\end{aligned}
\qquad (14) $$

with $p_u$ (resp. $p_w$) a function only of $u$ (resp. $w$), assuming $E, A$ and $p_\phi$ values already fixed by initial conditions. Then, we can separate again the action, leading to:

$$ \hat{S}[u, w, u_0, w_0, E, p_{\phi}] = S_u[u, E, A, p_{\phi}] + S_w[w, E, A, p_{\phi}] \qquad (15) $$

$$ 
\begin{aligned}
\frac{\partial S_u}{\partial u} &= p_u \\
\frac{\partial S_w}{\partial w} &= p_w
\end{aligned}
\qquad (16) $$

According to the Eq. (9), we have the following relations:

$$ 
\begin{aligned}
\left(\frac{dS_u}{du}\right)^2 &= \frac{m}{2}\left(E - \frac{p_\phi^2}{2mu^2} + \frac{GMm}{u} - \frac{A}{2mu} - \frac{mau}{2}\right) \\
&= \frac{m}{2}(E - V_u(u)) > 0 \\
\left(\frac{dS_w}{dw}\right)^2 &= \frac{m}{2}\left(E - \frac{p_\phi^2}{2mw^2} + \frac{GMm}{w} + \frac{A}{2mw} + \frac{maw}{2}\right) \\
&= \frac{m}{2}(E - V_w(w)) > 0
\end{aligned}
\qquad (17) $$

with

$$ 
\begin{aligned}
V_u(u) &= \frac{p_\phi^2}{2mu^2} - \frac{GMm}{u} + \frac{A}{2mu} - \frac{mau}{2} \\
V_w(w) &= \frac{p_\phi^2}{2mw^2} - \frac{GMm}{w} - \frac{A}{2mw} - \frac{maw}{2}
\end{aligned}
\qquad (18) $$

$V_u$ and $V_w$ are effective potentials applied in u and w directions (represented in the Fig. 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:

$$ P_3(u) < 0 \quad \text{and} \quad Q_3(w) > 0 \qquad (19) $$

These both conditions are more restrictive than the usual $E > E_p$ where $E_p$ is the potential energy.

### 2.3. Study of $P_3$

$P_3$ is a polynomial of degree 3 with $\lim_{u \to +\infty} P_3(u) = +\infty$. This polynomial possesses three roots, whose one is real at least. As $P_3(0) = p_\phi^2 > 0$, one of these roots is real negative, according to intermediate value theorem since $\lim_{u \to -\infty} P_3(u) = -\infty$.

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Representation in dimensionless units of the shape of both potentials $V_u$ (corresponding to $V_u$, blue line) and $V_w$ (corresponding to $V_w$, red line) for a set of $E, A$ and $p_\phi$ values. The motion is possible only in the area where the potential is below the mechanical energy $E$, represented by the black horizontal line. 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 $\epsilon$ for too high or too low $\epsilon$ values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Nevertheless, the motion occurs for positive u values, and we know this motion exists. It implies there is an interval in $\mathbb{R}^+$ such as $P_3 < 0$ (otherwise, there is no motion in u-direction, not 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 to $u \in [u_-; u_+]$ ($U \in [U_-; U_+]$) in terms of dimensionless quantities, as can be seen in the Fig. 2).

### 2.4. Study of $Q_3$

$Q_3$ is a polynomial of degree 3 with $\lim_{w \to +\infty} Q_3(w) = +\infty$. This polynomial possesses three roots, whose one is real at least. As $Q_3(0) = -p_\phi^2 < 0$, one of these roots is real positive, according to intermediate value theorem since $\lim_{w \to +\infty} Q_3(w) = +\infty$. 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 \in [w_-; w_+] \cup [w_0; +\infty[$ (as can be seen in the Fig. 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 \in [w_0; +\infty[$ only ($|w_0; +\infty[$ in terms of dimensionless quantities).

### 2.5. Restriction on the motion

Each constant value of u or w defines a parabolid 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 Fig. 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 Fig. 3 (left panel). Similarly, for the w-motion, there are two cases: the motion is constrained between one parabolid ($w = w_0$) and infinity or between two paraboloids ($w = w_-$ and $w = w_+$). Both cases are represented by the red area in the Fig. 3 (right panel).
<figure><img src="image_7.png" /><figcaption>Fig. 3. Representation of the available space for the motion. Left panel: the blue area corresponds to the allowed region for the u-motion. Right panel: the red area corresponds to the allowed region for the w-motion. The regions are limited by paraboloids (parabolas here by projection). The Sun location is shown with a yellow circle. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 2.6. 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 Fig. 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 Fig. 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 (Poincaré, 1890), 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_\phi$ is constant. All the positions in this part of the phase space can belong to the same trajectory. This is also linked to the Kolmogorov–Arnold–Moser theorem (Kolmogorov, 1954; Moser, 1962; Arnold, 1963): 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. As an interesting fact, the solution for the u-motion and w-motion are the exactly the same for governing the motion of the simple gravity pendulum.

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 extend the definition of the ballistic and satellite particles as presented in the introduction 1: 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.

## 2.7. 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:

$$R_{pressure} = \sqrt{\frac{GM}{a}} \qquad (20)$$

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:

<figure><img src="image_8.png" /><figcaption>Fig. 4. Left panel: Combination of both panels in the Fig. 3. The final allowed region is in green. Right panel: representation in green of the trajectory in the Fig. 1 in the $(x, \rho)$ frame at Earth. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
$$ 
\begin{aligned}
u &= UR_{\text{pressure}} \\
w &= WR_{\text{pressure}}
\end{aligned} 
\qquad (21) 
$$

where *U* and *W* are the dimensionless quantities associated to *u* and *w*. The energy *E* is a dimensionless quantity with the use of the characteristic energy $k_B T_{\text{exo}}$. For $p_u$ and $p_w$, we express them in units of $\sqrt{mk_B T_{\text{exo}}}$, whereas $p_\phi$ is expressed in units of $\sqrt{mk_B T_{\text{exo}} R_{\text{pressure}}}$.
Finally, we choose $mk_B T_{\text{exo}} \sqrt{GM/a}$ as the unit for *A* and the unit for a dimensionless time $\tau$ is derived from the units of *u/w* and $p_u/p_w$. We summarise these changes in Table 1.

We express all the previous equations as a function of these new quantities:

* the constants of the motion

$$ 
\begin{aligned}
\mathcal{E} &= \frac{2UP_U^2 + 2WP_W^2}{U+W} + \frac{P_\phi^2}{2UW} - \frac{2\lambda_a}{U+W} + \frac{\lambda_a}{2}(U-W)
\end{aligned} 
\qquad (22) 
$$

$$ 
\begin{aligned}
\mathcal{A} &= 2\mathcal{E}U - 4UP_U^2 - \frac{P_\phi^2}{U} + 2\lambda_a - \lambda_a U^2 \\
&= -2\mathcal{E}W + 4WP_W^2 + \frac{P_\phi^2}{W} - 2\lambda_a - \lambda_a W^2
\end{aligned} 
\qquad (23) 
$$

* the polynomials $P_3$ and $Q_3$:

$$ 
\begin{aligned}
P_3(U) &= \lambda_a U^3 - 2\mathcal{E}U^2 + (A - 2\lambda_a)U + P_\phi^2 \\
&= \lambda_a(U - U_0)(U - U_-)(U - U_+) \\
Q_3(W) &= \lambda_a W^3 + 2\mathcal{E}W^2 + (A + 2\lambda_a)W - P_\phi^2 \\
&= \lambda_a(W - W_-)(W - W_+)(W - W_0)
\end{aligned} 
\qquad (24) 
$$

* the effective potentials

$$ 
\begin{aligned}
V_U(U) &= \frac{P_\phi^2}{2U^2} - \frac{2\lambda_a - A}{2U} + \frac{\lambda_a U}{2} \\
V_W(W) &= \frac{P_\phi^2}{2W^2} - \frac{2\lambda_a + A}{2W} - \frac{\lambda_a W}{2}
\end{aligned} 
\qquad (25) 
$$

with $\lambda_a$:

<table><caption>**Table 1**

Compilation of the transformations of the parameters into dimensionless ones.</caption><thead><tr><th>Dimensional parameters</th><th>Unit</th><th>Dimensionless parameters</th></tr></thead><tbody><tr><td>u and w</td><td>√GM/a = R<sub>pressure</sub></td><td>U and W</td></tr><tr><td>p<sub>u</sub> and p<sub>w</sub></td><td>√mk<sub>B</sub>T<sub>exo</sub></td><td>P<sub>u</sub> and P<sub>w</sub></td></tr><tr><td>p<sub>φ</sub></td><td>√mk<sub>B</sub>T<sub>exo</sub>GM/a</td><td>P<sub>φ</sub></td></tr><tr><td>E</td><td>k<sub>B</sub>T<sub>exo</sub></td><td>ε</td></tr><tr><td>A</td><td>mk<sub>B</sub>T<sub>exo</sub>√GM/a</td><td>A</td></tr><tr><td>V<sub>u</sub> and V<sub>w</sub></td><td>k<sub>B</sub>T<sub>exo</sub></td><td>V<sub>u</sub> and V<sub>w</sub></td></tr><tr><td>t</td><td>√GM/mk<sub>B</sub>T<sub>exo</sub></td><td>τ</td></tr></tbody></table>

$$ \lambda_a = \frac{\sqrt{GMm}}{k_B T} = \frac{GMm}{k_B T R_{\text{pressure}}} = \lambda(R_{\text{pressure}}) \qquad (26) $$

the Jeans parameter at the distance $R_{\text{pressure}}$.

We can also reduce the equations of the motion. We introduce the dimensionless time $\tau$.

$$ 
\left\{ 
\begin{aligned}
\left(\frac{dU}{d\tau}\right)^2 &= -\frac{4P_3(U)}{(U+W)^2} \\
\left(\frac{dW}{d\tau}\right)^2 &= \frac{4Q_3(W)}{(U+W)^2} \\
\frac{d\phi}{d\tau} &= \frac{P_\phi}{UW}
\end{aligned}
\right. 
\qquad (27) 
$$

### 3. 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_\phi = 0$ and generalised by mathematical transformations to the 3D case in Lantoine and Russell (2011) with Jacobi elliptic functions (Jacobi, 1829). The analytical expression for the $\phi$-motion is missing in Lantoine and Russell (2011) for unbounded trajectories in 3D, although a computational precision comparison was performed between analytical and numerical approximated solutions (cf. Lantoine and Russell (2011), Table 2). Furthermore, the trajectories were completely solved by Biscani and Izzo (2014) with the Weierstrass functions for bounded and unbounded trajectories. These last works dealt with the dynamical trajectories of artificial satellites but they can apply 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$ polynomials. As presented in the introduction of this paper and detailed by Chamberlain (1963), in planetary exospheres with the collisionless hypothesis and only the gravity as external force, the trajectories of atoms and molecules can be divided into three kinds (cf. Fig. 5). First, the escaping particles come from the exobase and go to infinity, these are unbounded particles describing an hyperbola. Then, the ballistic particles are bounded coming from the exobase and describing an ellipse in order to cross twice the exobase. The trajectories are not elliptic or hyperbolic at all when the radiation pressure is included, but one can keep their basic definitions: crossing twice the exobase

**Table 2**

Summary of the different solutions for each kind of motion. The lack of notations means the equation are not explicitly derived in previous works in terms of Jacobi elliptic function. The symbol (s) corresponds to an example of formula with a less complex expression than the one proposed by Lantoine and Russell (2011). L2D3D labelled solutions mean the solutions were explicitly given by Lantoine and Russell (2011) in the 2D case and can be used in the 3D case by some transformations, whereas L3D labelled solutions were directly derived by Lantoine and Russell (2011) for the 3D case. We specify for each case the most time consuming functions and the number of evaluations at each T step. The function giving the sign is noted sg.

<table><thead><tr><th rowspan="2">Cases<br>Trajectory</th><th>Three real roots for P<sub>3</sub> (always)</th><th colspan="2">Three real roots for Q<sub>3</sub></th><th>One real root for Q<sub>3</sub></th></tr><tr><th>Bounded</th><th>Bounded</th><th>Unbounded</th><th>Unbounded</th></tr></thead><tbody><tr><td>u-motion</td><td>(31) (L2D3D)</td><td></td><td></td><td></td></tr><tr><td>w-motion</td><td></td><td>(36) (L2D3D)</td><td>(39) (L2D3D)</td><td>(44) (L2D3D,s)</td></tr><tr><td>Time equation</td><td>(48) (L2D3D)</td><td>(49) (L2D3D)</td><td>(50) (L2D3D)</td><td>(51) (L2D3D,s)</td></tr><tr><td>φ-motion</td><td>(53) (L3D)</td><td>(54) (L3D)</td><td>(55)</td><td>(56)</td></tr><tr><td>T range</td><td>T ∈ R</td><td>T ∈ R</td><td>T ∈ [0; (1 − sg(T<sub>W</sub>))2K(k<sub>W</sub>) + T<sub>W</sub>]</td><td>T ∈ [0; 2K(k<sub>W</sub>) + T<sub>W</sub>]</td></tr><tr><td>Evaluated functions</td><td>am(1), E(1), Π(1)</td><td>am(1), E(1), Π(1)</td><td>am(1), E(1), Π(1)</td><td>am(1), E(1), Π(1)</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Fig. 5. Definition of the three kinds of particle trajectories in a collisionless exosphere and without radiation pressure: ballistic (dot), satellite (cross) and escaping (dash) particles (from Beth et al. (2014)).</figcaption></figure>

for ballistic particles, orbiting but not crossing the exobase for satellite particles (see Fig. 1), 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 defined as:

$$ (U + W) dT = d\tau \qquad (28) $$

Until we know the expression of $U(T)$ and $W(T)$, we cannot express solutions as a function of time. The system of equations leads now to:

$$ \left\{ \begin{aligned} \left(\frac{dU}{dT}\right)^2 &= -4P_3(U) \\ \left(\frac{dW}{dT}\right)^2 &= 4Q_3(W) \\ \frac{d\phi}{dT} &= P_\phi\left(\frac{1}{U} + \frac{1}{W}\right) \\ \frac{d\tau}{dT} &= U + W \end{aligned} \right. \qquad (29) $$

All solutions, trajectories and time, are parametrized according to the variable T without physical sense.

### 3.1. The u-motion

We solve the differential equation describing the u-motion:

$$ \left(\frac{dU}{dT}\right)^2 = -4P_3(U) = -4\lambda_a(U - U_0)(U - U_-)(U - U_+) \qquad (30) $$

with $U_0$, $U_-$ and $U_+$ being real roots.

The solution for U is (cf. Appendix C):

$$ \begin{aligned} U(T) &= U_+ - (U_+ - U_-)\mathrm{sn}^2[\alpha(T - T_U), k_U] \\ k_U &= \sqrt{\frac{U_+ - U_-}{U_+ - U_0}} \\ \alpha &= \sqrt{\lambda_a} \sqrt{U_+ - U_0} \end{aligned} \qquad (31) $$

with cn and sn other Jacobi elliptic functions (cf. Appendix B). These functions are $4K(k_U)$-periodic (see Appendix A, Eq. (A.1)).

The initial conditions give us $U(0)$, then:

$$ \pm \operatorname{arcsn} \left( \sqrt{\frac{U_{+} - U(0)}{U_{+} - U_{-}}}, k_{U} \right) = \alpha T_{U} \qquad (32) $$

where $\operatorname{arcsn}$ is defined as the reciprocal function of sn:

$$ \operatorname{arcsn}(x, k) = \int_{0}^{x} \frac{1}{\sqrt{1 - t^2} \sqrt{1 - k^2 t^2}} dt = F(\operatorname{arcsin}(x), k) \qquad (33) $$

with F the incomplete elliptic integral of the first kind (see Appendix A, Eq. (A.1)).

The sign depends on $P_U$:

$$ \left(\frac{dU}{dT}\right) = 4UP_U \qquad (34) $$

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 $\eta$ given by Lantoinne 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)$.

### 3.2. The w-motion

For the w-motion, we need to solve the differential equation:

$$ \left(\frac{dW}{dT}\right)^2 = 4Q_3(W) = 4\lambda_a(W - W_0)(W - W_-)(W - W_+) \qquad (35) $$

$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.

#### 3.2.1. 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.

3.2.1.1. 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.

The solution for W is (cf. Appendix C):

$$ \begin{aligned} W(T) &= W_- + (W_+ - W_-)\mathrm{sn}^2[\beta(T - T_W), k_W] \\ k_W &= \sqrt{\frac{W_+ - W_-}{W_0 - W_-}} \\ \beta &= \sqrt{\lambda_a} \sqrt{W_0 - W_-} \end{aligned} \qquad (36) $$

The initial conditions giving us $W(0)$ and then leads to:

$$ \pm \operatorname{arcsn} \left( \sqrt{\frac{W(0) - W_-}{W_+ - W_-}}, k_W \right) = \beta T_W \qquad (37) $$

The sign depends on $P_W$:

$$ \left(\frac{dW}{dT}\right) = 4WP_W \qquad (38) $$

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 $\xi I$ given by Lantoinne and Russell (2011).

3.2.1.2. For $W > W_0$. First, the motion occurs only for $W \in [W_0; +\infty[$ (corresponding to escaping particles).
The final expression for W is (cf. Appendix C):

$$W(T) = W_0 + (W_0 - W_-) \mathrm{cs}^2[\beta(T - \mathcal{T}_W), k_W] \quad (39)$$

These functions are 4K($k_W$)-periodic and are defined on $T - \mathcal{T}_W \in \mathbb{R}/\{4mK(k_W)|m \in \mathbb{Z}\}$. These functions diverge at $4mK(k_W)$ but this is not an issue: the motion can diverge with respect to $T$, but, according to the time $\tau$ and the integration (see Section 4), the w-motion remains continuous for $\tau \in \mathbb{R}$ as

$$T - \mathcal{T}_W \in ]0; 4K(k_W)[ \Leftrightarrow \tau \in \mathbb{R}$$

Let us assume the initial conditions provide $W(0)$, then:

$$\pm \arccos \left( \sqrt{\frac{W(0) - W_0}{W_0 - W_{-}}}, k_W \right) = \beta \mathcal{T}_W \quad (40)$$

with

$$\arccos(x, k) = \int_{0}^{x} \frac{1}{\sqrt{1 + t^2 \sqrt{t^2 - 1 + k^2}}} dt = F\left(\arctan\left(\frac{1}{x}\right), k\right) \quad (41)$$

The sign depends on $P_W$:

$$\left(\frac{dW}{dT}\right) = 4WP_W \quad (42)$$

If $P_W$ is positive (resp. negative), $W$ increases (resp. decreases) and $\mathrm{cs}^2$ increases (resp. decreases) too. Then, $\mathcal{T}_W$ is positive (resp. negative). This solution corresponds to the solution $\xi_{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.

### 3.2.2. Only one real positive root

For any initial conditions, the particle will be escaping if $Q_3$ has only one real root $W_0$. 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 more direct way to determine the expression with the knowledge of all roots, even imaginaries. We provide the demonstration in Appendix C.

The final solution is:

$$\begin{align}
W(T) &= W_0 + \lambda^2 \frac{1 - \mathrm{cn}[2\gamma(T - \mathcal{T}_W_0), k_{W_0}]}{1 + \mathrm{cn}[2\gamma(T - \mathcal{T}_W_0), k_{W_0}]} \nonumber \\
\gamma &= \lambda \sqrt{\lambda_a} \nonumber \\
\lambda^2 &= \sqrt{(W_0 - W_+)(W_0 - W_-)} \tag{43} \\
k_{W_0} &= \sqrt{\frac{1}{2} - \frac{1}{4} \frac{W_0 - W_+ + W_0 - W_-}{\lambda^2}} \nonumber
\end{align}$$

Finally, one can simplify the expression based on Olver et al. (2010) (Eq. (22.6.18)):

$$W(T) = W_0 + \lambda^2 \frac{\mathrm{dn}^2[\gamma(T - \mathcal{T}_{W_0}), k_{W_0}]}{\mathrm{cs}^2[\gamma(T - \mathcal{T}_{W_0}), k_{W_0}]} \quad (44)$$

This function is 4K($k_{W_0}$)-periodic and is defined on $T - \mathcal{T}_W \in \mathbb{R}/\{2K(k_{W_0}) + 4mK(k_{W_0})|m \in \mathbb{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$, but, according to the time $\tau$ and the integration (see Section 4), the w-motion remains continuous for $\tau \in \mathbb{R}$ as

$$T - \mathcal{T}_{W_0} \in ] -2K(k_{W_0}); 2K(k_{W_0})[ \Leftrightarrow \tau \in \mathbb{R}$$

As usual, according to the initial conditions, we define $\mathcal{T}_{W_0}$ as

$$\pm \arccos \left( \frac{\lambda^2 - (W(0) - W_0)}{\lambda^2 + (W(0) - W_0)}, k_{W_0} \right) = 2\gamma \mathcal{T}_{W_0} \quad (45)$$

with

$$\begin{align}
\operatorname{arccn}(x, k) &= \int_x^1 \frac{1}{\sqrt{1 + t^2 \sqrt{1 - k^2 + k^2 t^2}}} \mathrm{d}t \nonumber \\
&= F(\operatorname{arccos}(x), k) \quad (46)
\end{align}$$

If $P_W$ is positive (resp. negative) then $\mathcal{T}_{W_0}$ is negative (resp. positive). This solution corresponds to the solution $\xi_{III}$ given by Lantoine and Russell (2011) but our derived solution is less complex to use. We notice here that to derive the *phi*-motion equation with the help of MAPLE, the form (43) was used as input rather than the form (44).

## 4. 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$, we can express the real time $\tau$ (or $t$):

$$\begin{align}
t &= \sqrt{\frac{GMm}{ak_B T_{\mathrm{exo}}}} \tau \nonumber \\
\tau(T) &= \int_0^T U(T') + W(T') \mathrm{d}T' \tag{47}
\end{align}$$

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:

$$\int_{0}^{T} U(T') \mathrm{d}T' = U_{0}T + \frac{\alpha}{\lambda_{a}}[E(\mathrm{am}[\alpha(T - \mathcal{T}_{U}), k_{U}], k_{U}) - E(\mathrm{am}[-\alpha T_{U}, k_{U}], k_{U})] \quad (48)$$

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 $W_+ > W(0) > W_-$ then

$$\begin{align}
\int_0^T W(T') \mathrm{d}T' &= W_0 T - \frac{\beta}{\lambda_a} [E(\mathrm{am}[\beta(T - \mathcal{T}_W), k_W], k_W) \nonumber \\
&\quad - E(\mathrm{am}[-\beta T_W, k_W], k_W)] \tag{49}
\end{align}$$

If $W(0) > W_0$ with three real roots then

$$\begin{align}
\int_0^T W(T') \mathrm{d}T' &= W_0 T - \frac{\beta}{\lambda_a} [E(\mathrm{am}[\beta(T - \mathcal{T}_W), k_W], k_W) \nonumber \\
&\quad - E(\mathrm{am}[-\beta T_W, k_W], k_W)] \nonumber \\
&\quad - \frac{\beta}{\lambda_a} \left[ \frac{\mathrm{cn}[\beta(T - \mathcal{T}_W), k_W]}{\mathrm{sd}[\beta(T - \mathcal{T}_W), k_W]} - \frac{\mathrm{cn}[-\beta T_W, k_W]}{\mathrm{sd}[-\beta T_W, k_W]} \right] \tag{50}
\end{align}$$

Finally, in the case of only one real root, the time equation is given by:

$$\begin{align}
\int_0^T W(T') \mathrm{d}T' &= (W_0 + \lambda^2)T - \frac{\gamma}{\lambda_a} [E(\mathrm{am}[2\gamma(T - \mathcal{T}_{W_0}), k_{W_0}], k_{W_0}) \nonumber \\
&\quad - E(\mathrm{am}[-2\gamma T_{W_0}, k_{W_0}], k_{W_0})] \nonumber \\
&\quad + \frac{\gamma}{\lambda_a} \left[ \frac{\mathrm{sn}[2\gamma(T - \mathcal{T}_{W_0}), k_{W_0}] \mathrm{dn}[2\gamma(T - \mathcal{T}_{W_0}), k_{W_0}]}{1 + \mathrm{cn}[2\gamma(T - \mathcal{T}_{W_0}), k_{W_0}]} - \frac{\mathrm{sn}[-2\gamma T_{W_0}, k_{W_0}] \mathrm{dn}[-2\gamma T_{W_0}, k_{W_0}]}{1 + \mathrm{cn}[-2\gamma T_{W_0}, k_{W_0}]} \right] \tag{51}
\end{align}$$
This equation should be exactly the same as the Eq. (77) from Lantoine and Russell (2011) in a shorter form.

For the Eq. (50),

$$T - \mathcal{T}_W \in ]0; 4K(k_W)[$$

and for the Eq. (51),

$$T - \mathcal{T}_{W_0} \in ] -2K(k_{W_0}); 2K(k_{W_0})[$$

Nevertheless, when $T - \mathcal{T}_W$ tends to $4K(k_{W_0})$ (resp. $2K(k_{W_0})$), the integration diverges and $\tau$ too. The w-motion occurs on an subset of $\mathbb{R}$ with respect to $T$ but on $\mathbb{R}$ entirely with respect to $\tau$. The transformation of $\tau$ into $T$ is bijective, i.e. only one $\tau$ can be associated to each $T$, because the integrand $U+W$ is strictly positive.

## 5. $\phi$-equation

To complete the description of the motion as a function of time, it is also necessary to solve the evolution of the angle $\phi$, obeying:

$$\begin{align} \frac{d\phi}{dT} &= P_{\phi}\left(\frac{1}{U(T)} + \frac{1}{W(T)}\right) \tag{52} \\ \phi(T) - \phi(0) &= \int_{0}^{T} P_{\phi}\left(\frac{1}{U(T')} + \frac{1}{W(T')}\right) dT' \nonumber \end{align}$$

As already done in the previous part, we separate into two integrals. The first part still gives:

$$\int_0^T \frac{1}{U(T')} dT' = \frac{\Pi\left(1 - \frac{U}{U_+}; \text{am}[\alpha(T - \mathcal{T}_U), k_U], k_U\right)}{\alpha U_+} - \frac{\Pi\left(1 - \frac{U}{U_+}; \text{am}[-\alpha T_U, k_U], k_U\right)}{\alpha U_+} \tag{53}$$

with $\Pi$ 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_-$:

$$\int_0^T \frac{1}{W(T')} dT' = \frac{\Pi\left(1 - \frac{W_+}{W}; \text{am}[\beta(T - \mathcal{T}_W), k_W], k_W\right)}{\beta W_-} - \frac{\Pi\left(1 - \frac{W_-}{W}; \text{am}[-\beta T_W, k_W], k_W\right)}{\beta W_-} \tag{54}$$

or if initially $W(0) > W_0$:

$$\int_0^T \frac{1}{W(T')} dT' = \frac{\beta T - \Pi\left(-\frac{W}{W_0-W_-}; \text{am}[\beta(T - \mathcal{T}_W), k_W], k_W\right)}{\beta W_-} + \frac{\Pi\left(-\frac{W}{W_0-W_-}; \text{am}[-\beta T_W, k_W], k_W\right)}{\beta W_-} \tag{55}$$

and for the last case with one real root:

$$\begin{align} \int_0^T \frac{1}{W(T')} dT' &= \frac{T}{W_0 - \lambda^2} \nonumber \\ &\quad - \frac{(W_0 + \lambda^2)}{4W_0\gamma(W_0 - \lambda^2)} \Pi\left(-\frac{(W_0 - \lambda^2)^2}{4\lambda^2 W_0}; \text{am}[2\gamma(T - \mathcal{T}_{W_0}), k_{W_0}], k_{W_0}\right) \nonumber \\ &\quad + \frac{1}{2P_\phi}\arctan\left(\frac{P_\phi}{2\gamma W_0} \frac{\text{sn}[2\gamma(T - \mathcal{T}_{W_0}), k_{W_0}]}{\text{dn}[2\gamma(T - \mathcal{T}_{W_0}), k_{W_0}]}\right) \nonumber \\ &\quad + \frac{(W_0 + \lambda^2)}{4W_0\gamma(W_0 - \lambda^2)} \Pi\left(-\frac{(W_0 - \lambda^2)^2}{4\lambda^2 W_0}; \text{am}[-2\gamma T_{W_0}, k_{W_0}], k_{W_0}\right) \nonumber \\ &\quad - \frac{1}{2P_\phi}\arctan\left(\frac{P_\phi}{2\gamma W_0} \frac{\text{sn}[-2\gamma T_{W_0}, k_{W_0}]}{\text{dn}[-2\gamma T_{W_0}, k_{W_0}]}\right) \tag{56} \end{align}$$

This last formula is only available for $T - \mathcal{T}_{W_0} \in ] -2K(k_{W_0}); 2K(k_{W_0})[$, range where it is continuous. The Eqs. (55) and (56) are not provided in Lantoine and Russell (2011) because it cannot be derived from the 2D case where $\phi$ is constant (modulo $\pi$). However, a computational comparison is performed between numerical and analytical approaches in Lantoine and Russell (2011), then we can supposed they derived the formula but the current authors of this paper do not know why this is not written explicitly mentioned.

## 6. 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 Fig. 6. It is clear that even a bounded trajectory the motion has no periodicity at all (especially for the $\phi$ 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 spacecraft (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:

$$-\frac{GMm\chi}{r^3} - ma = 0 \tag{57}$$

In dimensionless quantity, this will be expressed by:

$$R^2 + \cos \theta = 0 \tag{58}$$

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:

$$-\frac{GMm\rho}{r^3} + \frac{mv_{\phi}^2}{\rho^2} = 0 \tag{59}$$

In dimensionless quantity, this will be expressed by:

$$\lambda_a R \sin^4 \theta - P_{\phi}^2 = 0 \tag{60}$$

Combining these two equations, we obtain:

$$(\sin^2 \theta)^9 - (\sin^2 \theta)^8 + \frac{P_{\phi}^2}{\lambda_a} = 0 \tag{61}$$

We need to study the polynomial

$$P(X) = X^9 - X^8 + \frac{P_{\phi}^2}{\lambda_a} \tag{62}$$

with $X = \sin^2 \theta \in [0; 1]$. Depending on the $P_{\phi}^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 \in ]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 alternatively, if $P(a)$ is negative then we have two solutions. This value is:

$$P\left(\frac{8}{9}\right) = \frac{P_{\phi}^8}{\lambda_a^4} - \frac{1}{9}\left(\frac{8}{9}\right)^8 \tag{63}$$

A critical maximum value of $P_{\phi}$ thus exists to allow for circular orbits and is:

$$|P_{c\phi}| = \frac{8\sqrt{\lambda_a}}{9\sqrt{3}} \tag{64}$$
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Plots of a bounded particle (or spacecraft) motion in the (X, ρ) plane (upper left panel), of the time as a function of T (upper right panel), of the motion in polar coordinates (φ, ρ) (lower left panel) and the U-W coordinates as a function of the time (red for U, blue for W, green for R, lower right panel). U and W do not show any periodicity because their periods are not commensurable (i.e. the ratio is a rational number) with T and thus with the time τ. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Above this value, we cannot find any bounded trajectories: there is no equilibrium point (i.e. the Eqs. (58) and (60) do not have real solutions) 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). 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φ}|$, both these specific points exist: they have the same U. Physically, the potential $V_W$ has two extrema as plotted in Fig. 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$ has no extremum any more: the potential is strictly decreasing and the particles are unbounded (escaping). Thus, the bounded particles, satellite and ballistic particles, have $|P_φ| ≤ |P_{cφ}|$. 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 Eqs. (60) then (58), the critical orbit is:

$$ (R_{crit}, \theta_{crit}) = \left( \frac{1}{\sqrt{3}}, \pi - \arcsin \left( \frac{2\sqrt{2}}{3} \right) \right) \qquad (65) $$

or in (U, W) coordinates:

$$ (U_{crit}, W_{crit}) = \left( \frac{2}{3\sqrt{3}}, \frac{4}{3\sqrt{3}} \right) \qquad (66) $$

The real positive roots of the polynomial (62) combined with the equality (58) 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.

## 7. 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 summarise all needed equations in Table 2.

The u-motion is provided by the Eq. (31). The w-motion is provided by the Eqs. (36), (39) or (44). The φ-motion is provided by $P_φ \times [(53) + ((54) \text{ or } (55) \text{ or } (56))]$. The time equation is provided by [(48) + ((49) or (50) or (51))]. All the expressions are functions of T that is not the real time. We thus have implicit expressions as a function of time. The function $τ(T)$ is bijective but cannot be reversed analytically, a numerical inversion is needed to derive the real time. The problem is similar for the Keplerian case: the time equation $E - e \sin E = \frac{2\pi}{T}(t - t_0)$ where E is the true anomaly, T is the orbital period, e the eccentricity and t the time cannot allow to pass from t to E without a numerical approach.

According to the initial conditions, we can reconstruct the motion of the particle according to:

$$ \left\{ \begin{array}{l} X(T) = (U(T) - W(T))/2 \\ Y(T) = \sqrt{U(T)W(T)} \cos \phi(T) \\ Z(T) = \sqrt{U(T)W(T)} \sin \phi(T) \\ \tau(T) \text{ the true time} \end{array} \right. \qquad (67) $$

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 note 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 easily
transposed to the 2D case, whereas Lantoine and Russell (2011) provided first the equations for 2D case and then gave the necessary transformations to pass to the 3D case for u and w-motions. The time equation in the 3D case does not present any significant difference with the 2D case (addition of term depending linearly on T). However, the φ-motion expression for unbounded trajectories is missing from their paper. 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 of which the complexity prevent its implementation in all computer software and the need to work with complex values (e.g. the complex logarithm function). In this paper, we present 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 highlight in Table 2 which solutions are already derived by Lantoine and Russell (2011), in which case we propose less complex formulations to implement (e.g. Eq. (44) and (51)) and finally which solutions are not explicitly provided (e.g. Eq. (55) and (56)).

Moreover, beyond the new exact solutions given in this paper, the derivation of our solutions based on Jacobi elliptic functions allows a good computing time and accuracy. Hatten and Russell (2015) compared three types of solutions for the Stark effect: two exact ones, proposed by Lantoine and Russell (2011) (Jacobi elliptic functions) and Biscani and Izzo (2014) (Weierstrass elliptic functions), and a numerical one by Pellegrini et al. (2014) (based on Taylor series). They compared the CPU time, the number of calls for each analytic elliptic function and the accuracy between Biscani and Izzo (2014) and Lantoine and Russell (2011) in Python language. Even if we do not agree with the number of evaluations of each Jacobi elliptic function mentioned by Hatten and Russell (2015) (e.g. as am(x, k) is repeated between the different components of the motion, for each T step, only two different evaluations of am(x, k) are required, cf. Table 2), the use of these functions is always more valuable: 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)). Also, the analytical formulations are preferable for long duration motions.

## 8. Conclusions

We determined analytically the trajectories of the particles or spacecraft under the influence of both planetary gravity and stellar radiation pressure. This work may be used as an alternative method to previous works by Lantoine and Russell (2011), Biscani and Izzo (2014) and Pellegrini et al. (2014). We provide the complete exact solutions (t, x, y, z) 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 modelling 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 and critical value for the existence of bounded motions as well. Moreover, we discuss about the possible issues inherent to the formalism used and the importance of being extremely careful with the routines implemented.

Finally, we want to mention that this work is motivated by further purposes (Beth et al., 2016) to determine the distribution

function for exospheric particles subject to gravity and radiation pressure in the quasi-collisionless case. In the stationary case and if we assume the trajectory of particles is mainly determined by external forces, the Boltzmann equation can be reduced to:

$$ \frac{df}{dt} = v_s \frac{df_s}{ds} = P - L \qquad (68) $$

with s the curvilinear abscissa of the particle along its trajectory, v_s the velocity along, P the production term and L the loss term. Knowing the exact trajectory of the particles, the production and loss terms, one can then derive the distribution at any place according to initial or boundaries conditions. We thus focused on this efficient way to implement the follow-up of each particle for any initial conditions. Considering now the radiation pressure, the Eq. (68) becomes:

$$ \frac{df}{dt} = \frac{1}{U(T) + W(T)} \sqrt{\frac{k_B T_{\text{exo}} a}{G M m}} \frac{df}{dT} = P - L \qquad (69) $$

so that the integration along the trajectory can be performed according to the variable T.

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 (Beth et al., 2016) and later the 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.

## Acknowledgment

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 Π:

$$ \begin{aligned} F(\phi, k) &= \int_{0}^{\phi} \frac{1}{\sqrt{1 - k^2 \sin^2 \theta}} d\theta \\ E(\phi, k) &= \int_{0}^{\phi} \sqrt{1 - k^2 \sin^2 \theta} d\theta \\ \Pi(n; \phi, k) &= \int_{0}^{\phi} \frac{1}{(1 - n \sin^2 \theta) \sqrt{1 - k^2 \sin^2 \theta}} d\theta \end{aligned} \qquad (A.1) $$

Sometimes, other formulas (shown below) are proposed with the change sin θ = t but one needs to be very careful: this change is bijective only for θ ∈ [−π/2; π/2]

$$ \begin{aligned} \int_{0}^{x=\sin\phi} \frac{\sqrt{1-k^2 t^2}}{\sqrt{1-t^2}} dt &= \int_{0}^{\phi} \frac{\cos\theta}{|\cos\theta|} \sqrt{1-k^2 \sin^2\theta} d\theta \\ \int_{0}^{x=\sin\phi} \frac{1}{\sqrt{1-k^2 t^2 \sqrt{1-t^2}}} dt &= \int_{0}^{\phi} \frac{\cos\theta}{|\cos\theta|} \frac{1}{\sqrt{1-k^2 \sin^2\theta}} d\theta \\ \int_{0}^{x=\sin\phi} \frac{1}{(1-nt^2)\sqrt{1-k^2 t^2 \sqrt{1-t^2}}} dt &= \int_{0}^{\phi} \frac{\cos\theta}{|\cos\theta|} \frac{1}{(1-n\sin^2\theta)\sqrt{1-k^2 \sin^2\theta}} d\theta \end{aligned} \qquad (A.2) $$

These expressions are not exactly E, F and Π. They agree with the previous formulas (A.1) in the range φ ∈ [−π/2; π/2]. Lantoine and Russell (2011) did not state which formulations they used.
According to their formulas and results, they used the left-hand side of the Eq. (A.2). This may be a problem for bounded trajectories: for $\phi = am(T) \Rightarrow x = sn(T)$, the integrals (A.2) are not continuous contrary to (A.1). Depending on the computer software, the routines and the definitions used for these functions, the results can show some issues (e.g. no continuous motion).

## Appendix B. Jacobi elliptic functions

Most of the readers of this paper are not familiar with the so-called Jacobi elliptic functions (Jacobi, 1829). We proposed here to give a quick overview of these functions and their utility.

In Appendix A, we have already introduced the incomplete elliptic functions E, F and K.

### B.1. Jacobi amplitude function am(x, k)

By definition, this function is defined as the inverse function of F on $\mathbb{R}$ for $0 < k < 1$ as:

$$F(\text{am}(x, k), k) = \text{am}(F(\text{am}(x, k))) = x \quad (B.1)$$

The Jacobi amplitude function am is used in the case of the simple gravity pendulum: the evolution of the angle according to the time can be expressed with this function for any initial conditions.

### B.2. Jacobi elliptic functions cn(x, k), sn(x, k) and dn(x, k)

By definition, these functions are derived from the previous ones by the relations:

$$\begin{align}
\operatorname{sn}(x, k) &= \sin(\operatorname{am}(x, k)) \nonumber \\
\operatorname{cn}(x, k) &= \cos(\operatorname{am}(x, k)) \tag{B.2} \\
\operatorname{dn}(x, k) &= \sqrt{1 - \operatorname{sn}^2(x, k)} \nonumber
\end{align}$$

and any fractions of two of them define another one by the rules:

$$pq(x, k) = \frac{\operatorname{pr}(x, k)}{\operatorname{qr}(x, k)} = \frac{1}{\operatorname{qp}(x, k)} \quad (B.3)$$

and pp(x, k) = 1 with p, q or r being one of these four letters s, c, d or n. From the Eq. (B.1), the derivative of am is linked to dn by:

$$\frac{\mathrm{d} \operatorname{am}(x, k)}{\mathrm{d} x} = \mathrm{dn}(x, k) \quad (B.4)$$

## Appendix C. Demonstrations

We propose in this appendix the demonstrations of specific sections of the paper for interested readers.

### C.1. Section 3.1

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:

$$\begin{align}
\left(\frac{\mathrm{d}U}{\mathrm{d}T}\right)^2 &= -4\lambda_a(U - U_0)(U - U_)(U - U_+) \nonumber \\
4Y^2 \left(\frac{\mathrm{d}Y}{\mathrm{d}T}\right)^2 &= -4\lambda_a Y^2 \left(Y^2 \underbrace{-U_+ + U_0}_{<0} - U_+(U_0 - U_+)\right) \tag{C.1} \\
\left(\frac{\mathrm{d}Y}{\mathrm{d}T}\right)^2 &= -\lambda_a(Y^2 - U_+ + U_0)(Y^2 - U_+ + U_0) \nonumber \\
\left(\frac{\mathrm{d}Y}{\sqrt{\lambda_a} \mathrm{d}T}\right)^2 &= -(Y^2 - U_+ + U_0)(Y^2 - U_+ + U_0) \nonumber
\end{align}$$

Now, we set $Z = \frac{Y}{\sqrt{U_+ - U_0}}$.

$$\left(\frac{\mathrm{d}Z}{\sqrt{\lambda_a} \mathrm{d}T}\right)^2 = (1 - Z^2)(Z^2(U_+ - U_0) - U_+ + U_0) \quad (C.2)$$

Finally, we have:

$$\left(\frac{\mathrm{d}Z}{\sqrt{\lambda_a} \sqrt{U_+ - U_0} \mathrm{d}T}\right)^2 = (1 - Z^2) \left( Z^2 - \underbrace{\frac{U_+ - U_0}{U_+ - U_0}}_{<1 \text{ and } >0} \right) \quad (C.3)$$

We define $k_U$ as:

$$\frac{U_- - U_0}{U_+ - U_0} = 1 - \underbrace{\frac{U_+ - U_-}{U_+ - U_0}}_{<1 \text{ and } >0} = 1 - k_U^2 \quad (C.4)$$

$$k_U = \sqrt{\frac{U_+ - U_-}{U_+ - U_0}} \quad (C.5)$$

The final equation is:

$$\left(\frac{\mathrm{d}Z}{\sqrt{\lambda_a} \sqrt{U_+ - U_0} \mathrm{d}T}\right)^2 = (1 - Z^2)(Z^2 - (1 - k_U^2)) \quad (C.6)$$

The solution of this equation is:

$$Z = \mathrm{dn}\left[\sqrt{\lambda_a} \sqrt{U_+ - U_0}(T - T_U), k_U\right] \quad (C.7)$$

dn is a Jacobi elliptic function and $T_U$ depends on initial conditions.

### C.2. Section 3.2.1

After setting the two following substitutions $Y^2 = W_0 - W$ and $Z = \frac{Y}{\sqrt{W_0 - W}}$, we obtain:

$$Z = \mathrm{dn}\left[\sqrt{\lambda_a} \sqrt{W_0 - W_-}(T - T_W), k_W\right] \quad (C.8)$$

with

$$k_W = \sqrt{\frac{W_+ - W_-}{W_0 - W_-}} \quad (C.9)$$

### C.3. Section 3.2.1

We set $Y^2 = W - W_0$. We obtain:

$$\begin{align}
\left(\frac{\mathrm{d}W}{\mathrm{d}T}\right)^2 &= 4\lambda_a(W - W_0)(W - W_)(W - W_+) \nonumber \\
4Y^2 \left(\frac{\mathrm{d}Y}{\mathrm{d}T}\right)^2 &= 4\lambda_a Y^2 \left(Y^2 \underbrace{-W_+ + W_0}_{>0} - W_+(W_0 - W_+)\right) \tag{C.10} \\
\left(\frac{\mathrm{d}Y}{\mathrm{d}T}\right)^2 &= \lambda_a(Y^2 - W_+ + W_0)(Y^2 - W_+ + W_0) \nonumber \\
\left(\frac{\mathrm{d}Y}{\sqrt{\lambda_a} \mathrm{d}T}\right)^2 &= (Y^2 - W_+ + W_0)(Y^2 - W_+ + W_0) \nonumber
\end{align}$$

Now, we set $Z = \frac{Y}{\sqrt{W_0 - W}}$.

Finally, we have:

$$\left(\frac{\mathrm{d}Z}{\sqrt{\lambda_a} \sqrt{W_0 - W_-} \mathrm{d}T}\right)^2 = (Z^2 + 1) \left( Z^2 + \underbrace{\frac{W_0 - W_+}{W_0 - W_-}}_{<1 \text{ and } >0} \right) \quad (C.11)$$
$$ \frac{W_0 - W_+}{W_0 - W_-} = 1 - \underbrace{\frac{W_+ - W_-}{W_0 - W_-}}_{<1 \text{ and } >0} = 1 - k_W^2 \quad (C.12) $$

The solution of this equation is

$$ Z = \operatorname{cs} \left[ \sqrt{\lambda_a} \sqrt{W_0 - W_-} (T - T_W), k_W \right] \quad (C.13) $$

cs is a Jacobi elliptic function and $T_W$ depends on the initial conditions.

### C.4. Section 3.2.2

We start from the Eq. (C.14) with $W_+$ and $W_-$ not real:

$$ \left(\frac{dW}{dT}\right)^2 = 4\lambda_a(W - W_0)(W - W_+)(W - W_-) \quad (C.14) $$

After the separation of variables, we obtain:

$$ \int^{W} \frac{dW'}{\sqrt{(W' - W_0)(W' - W_+)(W' - W_-)}} = 2\sqrt{\lambda_a}(T - T_W) \quad (C.15) $$

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:

$$ \begin{aligned} \lambda^2 &= \sqrt{(W_0 - W_+)(W_0 - W_-)} \\ &= \sqrt{(W_0 - \operatorname{Re}(W_+))^2 + \operatorname{Im}(W_+)^2} = \sqrt{\frac{Q_3'(W_0)}{\lambda_a}} \end{aligned} \quad (C.16) $$

and also

$$ k_{W_0} = \sqrt{\frac{1}{2} - \frac{1}{4} \frac{W_0 - W_+ + W_0 - W_-}{\lambda^2}} = \sqrt{\frac{1}{2} - \frac{1}{2} \frac{W_0 - \operatorname{Re}(W_+)}{\lambda^2}} \quad (C.17) $$

According to Abramowitz and Stegun (1964, p. 597), the left hand side corresponds to:

$$ \int_{W_0}^{W} \frac{dW'}{\sqrt{(W' - W_0)(W' - W_+)(W' - W_-)}} = \frac{F(\theta, k_{W_0})}{\lambda} \quad (C.18) $$

with $F$ the Elliptic function of first kind (defined in Appendix A) and

$$ \cos \theta = \frac{\lambda^2 - (W - W_0)}{\lambda^2 + (W - W_0)} \Rightarrow W = W_0 + \lambda^2 \frac{1 - \cos \theta}{1 + \cos \theta} \quad (C.19) $$

According to the Eq. (C.15):

$$ \theta = \operatorname{am} \left[ 2\lambda \sqrt{\lambda_a} (T - T_{W_0}) k_{W_0} \right] \quad (C.20) $$

Thus,

$$ W(T) = W_0 + \lambda^2 \frac{1 - \operatorname{cn}[2\lambda\sqrt{\lambda_a}(T - T_{W_0}), k_{W_0}]}{1 + \operatorname{cn}[2\lambda\sqrt{\lambda_a}(T - T_{W_0}), k_{W_0}]} \quad (C.21) $$

### References

* Abramowitz, M., Stegun, I.A., 1964. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York (ninth dover printing, tenth gpo printing edition).

* Arnol'd, V.I., 1963. Proof of a theorem of A.N. Kolmogorov on the invariance of quasi-periodic motions under small perturbations of the Hamiltonian. Russ. Math. Surv. 18 (5), 9–36.

* Beth, A. et al., 2014. Modeling the satellite particle population in the planetary exospheres: Application to Earth, Titan and Mars. Icarus 227, 21–36.

* Beth, A. et al., 2016. Theory for planetary exospheres: II. Radiation pressure effect on exospheric density profiles. Icarus 266, 423–432.

* Biscani, F., Izzo, D., 2014. The Stark problem in the Weierstrassian formalism. Mon. Not. R. Astron. Soc. 439, 810–822.

* Bishop, J., 1991. Analytic exosphere models for geocoronal applications. Planet. Space Sci. 39, 885–893.

* Bishop, J., Chamberlain, J.W., 1989. Radiation pressure dynamics in planetary exospheres: A natural framework. Icarus 81, 145–163.

* Burns, J.A., Lamy, P.L., Soter, S., 1979. Radiation forces on small particles in the Solar System. Icarus 40, 1–48.

* Chamberlain, J.W., 1963. Planetary coronae and atmospheric evaporation. Planet. Space Sci. 11, 901–960.

* Ferziger, J., Kaper, H., 1972. Mathematical Theory of Transport Processes in Gases. North-Holland Publ.

* Hatten, N., Russell, R.P., 2015. Comparison of three stark problem solution techniques for the bounded case. Celest. Mech. Dynam. Astron. 121, 39–60.

* Jacobi, C.G.J., 1829. Fundamenta nova theoriae functionum ellipticarum. Sumtibus fratrum.

* Jeans, J.H., 1916. The Dynamical Theory of Gases. Cambridge University Press (4th ed., 1925).

* Kolmogorov, A.N., 1954. On conservation of conditionally periodic motions for a small change in Hamilton's function. Dokl. Akad. Nauk SSSR 98, 527–530.

* Lantoine, G., Russell, R.P., 2011. Complete closed-form solutions of the stark problem. Celest. Mech. Dynam. Astron. 109, 333–366.

* Moser, J., 1962. On Invariant Curves of Area-Preserving Mappings of An Annulus. Nachrichten der Akademie der Wissenschaften in Göttingen. Vandenhoeck & Ruprecht.

* Namouni, F., Guzzo, M., 2007. The accelerated Kepler problem. Celest. Mech. Dynam. Astron. 99, 31–44.

* Olver, F.W. et al., 2010. NIST Handbook of Mathematical Functions, first ed. Cambridge University Press, New York, NY, USA.

* Pellegrini, E., Russell, R.P., Vittaldev, V., 2014. F and g Taylor series solutions to the Stark and Kepler problems with Sundman transformations. Celest. Mech. Dynam. Astron. 118, 355–378.

* Poincaré, H., 1890. Théorie des invariants intégraux. Acta Math. 13, 46–88.

* Redmond, P.J., 1964. Generalization of the Runge-Lenz vector in the presence of an electric field. Phys. Rev. 133, B1352-B1353.

* Sommerfeld, A., 1934. Atomic Structure and Spectral Lines, third ed., vol. 1. Methuen & Co.

* Stark, J., 1914. Beobachtungen ber den effekt des elektrischen feldes auf spektrallinien. i. quereffekt. Ann. Phys. 348, 965–982.