# Theory for planetary exospheres: II. Radiation pressure effect on exospheric density profiles

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: II. Radiation pressure effect on exospheric density profiles. Icarus, 2016, 266, pp.423-432. 10.1016/j.icarus.2015.08.023. insu-03670535

HAL Id: insu-03670535

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

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: II. Radiation pressure effect on exospheric density profiles

<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 27 March 2015

Revised 6 August 2015

Accepted 12 August 2015

Available online 2 September 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 generalization of the study by Bishop and Chamberlain (1989). In this second part of our work, we present here the density profiles of atomic Hydrogen in planetary exospheres subject to the radiation pressure. We first provide the altitude profiles of ballistic particles (the dominant exospheric population in most cases), which exhibit strong asymmetries that explain the known geotail phenomenon at Earth. The radiation pressure strongly enhances the densities compared with the pure gravity case (i.e. the Chamberlain profiles), in particular at noon and midnight. We finally show the existence of an exopause that appears naturally as the external limit for bounded particles, above which all particles are escaping.

© 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).

To describe correctly the exospheric population, we distinguish three types of particles: escaping, ballistic and satellite (Chamberlain, 1963; Banks and Kockarts, 1973).

* 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).

* The ballistic particles also come from the exobase but with 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.

* 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 from ballistic particles undergoing few collisions mainly near the exobase. Thus, they do not exist in a collisionless model of the exosphere.

By definition, their trajectories are conics in the pure gravity case. Chamberlain (1963) proposed an analytical approach to estimate the density of each population via Liouville's theorem which states that the distribution function remains constant along a dynamical trajectory. A Maxwellian distribution function is assumed at the exobase and propagated to the upper layers via

<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>E-mail addresses: arnaud.beth@gmail.com (A. Beth), pgarnier@irap.omp.eu (P. Garnier).</footnote>

<footnote n="1">¹ Corresponding author at: Université de Toulouse; UPS-OMP; IRAP, Toulouse, France.</footnote>
Liouville's theorem. The density for each population is then derived as the product between the barometric law and a partition function $\zeta$.

$$n(r) = n_{\text{bar}} \zeta(\lambda) = n(r_{\text{exo}}) e^{\lambda - \lambda_{\text{exo}}} (\zeta_{\text{bal}} + \zeta_{\text{esc}}) \quad (1)$$

where $\lambda$ is the ratio between the gravitational and thermal energies.

$$\dot{\lambda}(r) = \frac{GMm}{k_B T_{\text{exo}} r} = \frac{v_{\text{esc}}(r)^2}{U^2} \quad (2)$$

with $r$ the distance from the center of the body, $v_{\text{esc}}(r)$ the escaping velocity, $U$ the most probable velocity for the Maxwellian distribution, $G$ the gravitational constant, $M$ the mass of the planet or the satellite and $T_{\text{exo}}$ the temperature at the exobase considered constant in the exosphere.

The radiation pressure disturbs the ellipses or hyperbolas described by these particles. The resonant scattering of solar photons leads to a total momentum transfer from the photon to the atom or molecule. 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 differential of momentum between before and after the scattering each second imparts a force, the radiation pressure. Bishop and Chamberlain (1989) proposed to take into account this effect on the structure of planetary exospheres. In particular, they highlighted analytically the “tail” phenomenon: the density for atomic Hydrogen is higher in the nightside direction than in the dayside direction, as observed for the first time by OGO-5 (Thomas and Bohlin, 1972; Bertaux and Blamont, 1973).

This problem is similar to so-called Stark effect (Stark, 1914): the effect of a constant electric field on the atomic Hydrogen's electron. Its study shows it can be transposed to celestial mechanics in order to describe the orbits of artificial and natural satellites in the perturbed two-body problem. A first but incomplete work was performed by Bishop and Chamberlain (1989). They focused on the density profiles along the Sun-planet axis: in the velocity phase space, the problem is thus only 2D (one component of the angular momentum is null, $p_{\phi}$, and thus the problem takes place on a hyperplane in the 3D-velocity phase space). They determined the density profiles for bounded trajectories (only ballistic particles, neither escaping nor satellite particles) for atomic Hydrogen along the Sun-planet axis, on the dayside and the nightside, for Earth, Venus, Mars or for sodium at Mercury.

In this work, we generalize the formalism developed by Bishop and Chamberlain (1989) to the whole exosphere (3D case) and highlight several phenomena. Our study is based on Beth et al. (2016), where we detailed the dynamical aspects induced by the radiation pressure on the trajectories of exospheric particles. We now present the implications on exospheric density profiles, local time asymmetries as well as a specific study of the particles with satellite orbits. We will briefly describe the formalism used in Section 2, before we derive the neutrals density in Section 3, and present the results in Section 4 and conclude in Section 5.

## 2. Model

For this work, we decide to study the effect of the radiation pressure on atomic Hydrogen in particular. We model the radiation pressure by a constant acceleration $a$ coming from the Sun. As previously defined by Bishop (1991), this acceleration depends on the line center solar Lyman-$\alpha$ flux $f_0$, in $10^{11}$ photons cm$^{-2}$ s$^{-1}$ Å$^{-1}$:

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

This problem is similar to the classical Stark effect (Stark, 1914): a constant electric field (here the radiation pressure) applied on 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) on the electron (the Hydrogen atom), 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) adopting 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} \quad (4)$$

where $x$ is the sunward coordinate and $\theta$ the angle with respect to the Sun-planet axis. Along the Sun-planet axis, $w$ is null in the sunward direction whereas $u$ is null in the nightside direction. 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 (5)$$

independent of $t$, the time and $\phi$, the azimuth about the planet-Sun axis. $p_u$, $p_w$ and $p_\phi$ are the conjugate momenta, $GM$ the standard gravitational parameter of the planet and $m$ the mass of the species. According to canonical relations, we have:

$$\begin{aligned} p_u &= \frac{m(u+w)}{4u} \frac{du}{dt} \\ p_w &= \frac{m(u+w)}{4w} \frac{dw}{dt} \\ p_\phi &= muw \frac{d\phi}{dt} \end{aligned} \quad (6)$$

We do not assume $p_\phi = 0$ as Bishop and Chamberlain (1989) did: their study is restricted to the Sun-planet axis where either $u = 0$ or $w = 0$.

As shown by Bishop and Chamberlain (1989), the problem has three constants of the motion: $E$ the mechanical energy, $p_\phi$ and $A$ defined as

$$E = \mathcal{H} \quad (7)$$

because the forces are conservative,

$$\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 separation constant (Bishop and Chamberlain, 1989) which is similar to the norm of the Laplace–Runge–Lenz vector and $p_\phi$ because

$$\frac{dp_\phi}{dt} = - \frac{\partial \mathcal{H}}{\partial \dot{\phi}} = 0 \quad (9)$$

As function of these three constants, we can rewrite the conjugate momenta:

$$\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 (10)$$

with

$$\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 (11)$$
<table><caption>**Table 1**
Conditions on the roots of the polynomials and on the exobase location to define the ballistic, satellite and escaping trajectories.</caption><thead><tr><th>Type of trajectory</th><th>Ballistic</th><th>Satellite</th><th>Escaping</th></tr></thead><tbody><tr><td>Positive real roots of P<sub>3</sub></td><td>2</td><td>2</td><td>2</td></tr><tr><td>Positive real roots of Q<sub>3</sub></td><td>3</td><td>3</td><td>1</td></tr><tr><td>Condition on R<sub>pressure</sub><sup>reco</sup></td><td>&gt; U<sub>+</sub>+W<sub>-</sub>/2</td><td>&lt; U<sub>+</sub>+W<sub>-</sub>/2</td><td>&gt; U<sub>+</sub>+W<sub>0</sub>/2</td></tr><tr><td>Necessary and sufficient conditions?</td><td>Yes</td><td>Yes</td><td>No, need for tracking</td></tr></tbody></table>

We can remark here that for $p_{\phi} = 0, 0$ is a simple root of $P_3$ and $Q_3$.

Based on the formalism detailed in Beth et al. (2016), we use dimensionless quantities and the same annotations:

$$ \varepsilon = \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) \quad (12) $$

$$ \begin{aligned} A &= 2\varepsilon U - 4UP_{U}^{2} - \frac{P_{\phi}^{2}}{U} + 2\lambda_{a} - \lambda_{a}U^{2} \\ &= -2\varepsilon W + 4WP_{W}^{2} + \frac{P_{\phi}^{2}}{W} - 2\lambda_{a} - \lambda_{a}W^{2} \end{aligned} \quad (13) $$

$$ \begin{aligned} P_3(U) &= \lambda_a U^3 - 2\varepsilon U^2 + (A - 2\lambda_a)U + P_\phi^2 \\ &= \lambda_a(U - U_0)(U - U_-)(U - U_+) \end{aligned} \quad (14) $$

$$ \begin{aligned} Q_3(W) &= \lambda_a W^3 + 2\varepsilon W^2 + (A + 2\lambda_a)W - P_\phi^2 \\ &= \lambda_a(W - W_-)(W - W_+) (W - W_0) \end{aligned} $$

with $U_0, U_-, U_+$ and $W_0$ real roots such as $U_0 < 0 < U_- < U_+; W_-$ and $W_+$ can be real (then $W_- < W_+ < W_0$) or complex conjugates. For each polynomial, one root can be 0 if $P_{\phi} = 0$. We define $\lambda_a$ as:

$$ \lambda_a = \frac{\sqrt{GMam}}{k_B T} = \frac{GMm}{k_B T R_{\text{pressure}}} = \lambda(R_{\text{pressure}}) \quad (15) $$

the Jeans parameter at the distance $R_{\text{pressure}} = \sqrt{GM/a}$.

### 3. Calculation of exospheric densities

The density in a gas is given from the distribution function $f$ by:

$$ n(\vec{r}) = \int f(\vec{r}, \vec{v}) d^3 \vec{v} \quad (16) $$

The bounds of integration depend on the type of particle. As for the so-called Chamberlain (1963) model, we distinguish three types of particles: ballistic, satellite and escaping particles. As detailed in the next section (and summarized in Table 1), the trajectory of these particles are not conics at all but they keep a part of their definition in this problem: the ballistic particles cross the exobase twice, the satellite particles do not cross the exobase and do not escape, and the escaping particles cross the exobase once before they escape. As shown by Beth et al. (2016), each type of particles must respect conditions about the roots of $P_3$ and $Q_3$, summarized in Table 1.

#### 3.1. Reminder about the previous analytical work by Bishop and Chamberlain (1989)

Before performing the density profile calculation, we specify here the main differences with the initial work by Bishop and Chamberlain (1989).

Bishop and Chamberlain (1989) performed 1D calculations along the Sun-planet axis only, implying some assumptions on

the trajectories of the particles (described below). We explain here rigorously why their study and ours are different and complementary.

So that a particle crosses the Sun-planet axis, $U$ or $W$ will be null once and thus $P_{\phi}$ too (cf. Eq. (6)) because it is independent of the time. Thus, $P_3$ and $Q_3$ have 0 as common root. For bounded trajectories, $W_-, W_+$ and $W_0$ are real positive and $0 < W_- < W_+ < W_0$ thus $W_- = 0$. The question now is: which root of $P_3$ is null? As explained further, $U_0 + U_+ + U_- = -W_+ - W_0 < 0$ and $U_+$ is strictly positive. Then, between $U_-$ and $U_0$, one is negative, the other one is null thus $U_0 < 0$ and $U_- = 0$.

The $(U, W)$-motion of a particle crossing the Sun planet axis is $(U, W) \in [0; U_+] \times [0; W_+]$. According to the Poincaré recurrence theorem (Poincaré, 1890) as used by Beth et al. (2016), these particles will necessarily pass as close as we want to the $(0, 0)$ position. Consequently, as explained by Bishop and Chamberlain (1989) for the planar motion ($P_{\phi} = 0$), all bounded trajectories cross the exobase. These authors considered the particles orbiting several times around the planet as satellite particles with a finite lifetime. We will show in the following sections for our 3D-case, that pure satellite particles exist which never cross the exobase (see Section 4.3). Finally, our integrals here will be calculated in another way compared with Bishop and Chamberlain (1989) since $P_{\phi}$ can vary. Although our model cannot take into account the singularity $P_{\phi} = 0$, this is not a problem for the integration (for the numerical evaluation, we just avoid this value).

#### 3.2. Densities of ballistic particles

We will first focus on the densities of ballistic particles. These particles are trapped by the potential that includes both the gravity and the radiation pressure. Consequently, the parabolic coordinates $U$ and $W$ must have finite values and belong to a closed interval from $\mathbb{R}^+$. On the first hand, there is no further constraint for $U$-values, since the $U$-motion is always limited by two paraboloids defined by $U_-$ and $U_+$ (see Section 2.4 in Beth et al., 2016). On the other side, the $W$-motion must respect the following condition to allow for ballistic particles: the roots $W_-$ and $W_+$ must be real positive (otherwise, they can be complex or real negative). This implies a condition on the energy $\varepsilon$ according to Eq. (14):

$$ \lambda_a(W_0 + W_- + W_+) = -2\varepsilon > 0 \quad (17) $$

As explained by Beth et al. (2016), the $U$-motion does not have the same period as the $W$-motion or the $\phi$-motion. In the $(U, W)$-plane, the trajectory fills the whole square $[U_-: U_+] \times [W_-; W_+]$ (see Beth et al., 2016, Fig. 1). If a particle belongs to this square, after a certain time, the particle will pass as close as we want to any position chosen in this square (except for specific periodic cases, cf. Biscani and Izzo (2014), occurring only if we are looking for them). Thus, if a part of the exobase surface belongs to this square, the particle in this square will necessarily cross several times the exobase (at the distance $r_{\text{exo}}$). We will thus be able to call it a ballistic particle. By definition $u + w = (U + W)R_{\text{pressure}} = 2r$, and the nearest distance from the
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Representations of areas where ballistic (blue) and satellite (green) particles evolve in the usual coordinate system (left panel) and in the (U, W) plane. To pass from the left grid to the right one, we apply the transformation (x, y) → (u = √x² + y² + x, w = √x² + y² - x). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

center of the planet where the particle can pass is $r_{-} = (U_{-} + W_{-})R_{pressure}/2$. If $W_{-}$ is real positive, then $r_{-}$ allow us to define the type of this particle: if $r_{-} < r_{exo}$ the particle is ballistic, if $r_{-} > r_{exo}$ the particle is satellite.

In summary, the bounds of integration needed to calculate the density of ballistic particles from Eq. (16) is the part of the velocity space where $W_{-}$ and $W_{+}$ are real positive and $r_{-} < r_{exo}$.

A last parameter must be determined: what is the distribution function to be used? The calculation of the ballistic density (detailed below) is based on the Liouville theorem, which uses the conservation of the distribution function along a dynamical trajectory. Since the ballistic particles cross the exobase, i.e. the external limit of the collisional atmosphere, one can reasonably assume (as in Chamberlain, 1963) that the distribution function is Maxwellian at the exobase.

### 3.3. Description of the algorithm for ballistic particles

**First step.** We choose a specific location in the exosphere and we thus fix the values of U and W.

**Second step.** Then, we scan all possibilities for the velocity vector. For each 3-tuple (three components of the velocity vector), we calculate the corresponding three constants of the motion ($\mathcal{E}$, $A$, $P_{\phi}$).

**Third step.** For each tuple ($\mathcal{E}$, $A$, $P_{\phi}$), we calculate the three roots of $P_3$ and $Q_3$.

**Fourth step.** We test if the three roots of $Q_3$ are real and positive.

**Fifth step.** We calculate $U_{-}/2 + W_{-}/2$. If this value is below $r_{exo}/R_{pressure}$ then this is a bounded trajectory crossing the exobase, that corresponds to ballistic particles. Thus, the Liouville theorem can be applied.

**Sixth step.** We integrate the velocity distribution function and derive the density of ballistic particles.

In order to calculate the ballistic density $n_{bal}$, we propose to define a partition function $\zeta_{bal}(r, \theta)$ that is the generalization, when

the radiation pressure effect is included, of the 1D partition function defined by Chamberlain (1963):

$$n_{bal}(r, \theta) = \underbrace{n_{exo} \exp(\lambda(r) - \lambda_c)}_{\text{barometric law}} \exp\left(-\lambda_a \frac{r - r_{exo}}{R_{pressure}} \cos \theta\right) \zeta_{bal}(r, \theta) \quad (18)$$

with $n_{exo}$ the exobase density, $\theta$ the angle with the respect to the Sun-planet axis and $r$ the distance from the planet. By definition, for a Maxwellian distribution,

$$\zeta_{bal}(r, \theta) = \frac{\int_{bal} \exp\left(-\frac{p^2}{2mk_B T}\right) d^3\vec{p}}{\int \exp\left(-\frac{p^2}{2mk_B T}\right) d^3\vec{p}} \quad (19)$$

We choose to work with dimensionless quantities and the (U, W) coordinates. Operating some transformations and with the correct Jacobian, we obtain:

$$\zeta_{bal}(U, W) = \frac{1}{(2\pi)^{3/2}} \frac{4}{U+W} \int_{bal} \exp\left(-\frac{2UP_U^2 + 2WP_W^2}{U+W} - \frac{P_{\phi}^2}{2UW}\right) dP_U dP_W dP_{\phi} \quad (20)$$

The formula (20) is only available for $U \neq 0$ or $W \neq 0$. Indeed, if $U = 0$ or $W = 0$, then $P_{\phi} = 0$ and this integral cannot be performed. Thus, our formulation is only available for the whole exosphere except along the Sun-planet axis, already studied by Bishop and Chamberlain (1989).

As previously mentioned, $\mathcal{E}$ is negative. It represents the contribution of the potential energy and the kinetic energy. Thus, we have the inequality from the Eq. (12):

$$\frac{2\lambda_a}{U+W} - \frac{\lambda_a}{2}(U-W) > \frac{2UP_U^2 + 2WP_W^2}{U+W} + \frac{P_{\phi}^2}{2UW} > 0 \quad (21)$$

The modulus inside the exponential takes finite values. In this case, we choose the following change of coordinates:
$$
\begin{align}
X &= \sqrt{\frac{2U}{U+W}} \qquad P_U = R \sin \Theta \cos \Phi \notag \\
Y &= \sqrt{\frac{2W}{U+W}} \qquad P_W = R \sin \Theta \sin \Phi \notag \\
Z &= \frac{1}{\sqrt{2UW}} \qquad P_{\Phi} = R \cos \Theta \tag{22}
\end{align}
$$

$$
\begin{align}
\zeta_{\text{type}}(U, W) &= \frac{1}{\pi^{3/2}} \int_{J_{\text{bal}}} \exp(-X^2 - Y^2 - Z^2) \, dX \, dY \, dZ \tag{23} \\
&= \frac{1}{\pi^{3/2}} \int_{J_{\text{bal}}} R^2 \exp(-R^2) \, dR \, d\Theta \, d\Phi \notag
\end{align}
$$

To estimate this integral (23), we use the Gauss–Legendre
quadrature as:

$$
\begin{align}
\zeta_{\text{bal}}(U, W) &= \frac{1}{\pi^{3/2}} \int \d[\text{bal}](R, \Theta, \Phi) R^2 \exp(-R^2) \, dR \, d\Theta \, d\Phi \tag{24} \\
&\approx \frac{1}{\pi^{3/2}} \sum_{i=1}^{N} \sum_{j=1}^{N} \sum_{k=1}^{N} \d[\text{bal}](R_i, \Theta_j, \Phi_k) R_i^2 \exp(-R_i^2) w_i w_j w_k \notag
\end{align}
$$

N is the number of points used for the integration and $\d[\text{bal}]$ is a function taking the value 1 if the particle is bounded and crosses the exobase with these initial conditions or 0 otherwise, as a rejection sampling (cf. Fig. 2).

The Gauss–Legendre method is particularly efficient in this case because all bounds of integration are finite: *R* is finite (see Eq. (21)), *Θ* belongs to [0; *π*) and *Φ* to [0; 2*π*].

3.4. The escaping and the satellite particles

Our approach unfortunately cannot be applied so easily to cal-
culate the densities of particles with escaping or satellite
trajectories.

For escaping particles, the coordinates in the full position-
velocity phase space do not guarantee that the particles come from
the exobase or not. Indeed, for escaping particles, the volume is
opened (no restriction on $\varepsilon$). The Poincaré recurrence theorem
(Poincaré, 1890) thus cannot be applied and a particle can have
the conditions to be escaping without crossing the exobase (the
particle is just passing by the planet): it is necessary to follow
the particle along the time as long as it is inside the exopause. If
a particle reaches the exopause without crossing the exobase, this
is not a escaping particle. We tried to track the particles to know if
they cross or not the exobase but we have some time and precision
issues: compared with the ballistic particles, we shall compute the
trajectory of each particle (problem of time) and integrate the
energies until the infinite (problem of precision). Several attempts
were performed without convincing results. However, as will be
discussed in a future paper, it is possible to calculate the analytical
escape flux at the exobase.

The satellite particles cannot exist in our model because this is a
collisionless model. We have previously proposed a formalism to
estimate their density (Beth et al., 2014); the trajectories are not
closed here (no periodic motions for all bounded cases) and the
formalism of this paper cannot be adjusted in this way. However,
we will show below (Section 4.3) that satellite particles can exist
in the presence of a radiation pressure force, and where these par-
ticles are.

4. Results

Bishop and Chamberlain (1989) provided only the ballistic den-
sity along the Sun-planet axis, along the dayside and nightside
directions. Here, we generalize this approach with a 2D model
(3D if an axisymmetric symmetry is considered) and provide the

ballistic density (main component in the lower part of the exo-
sphere) in every direction from the planet. We present here the
results for different planetary exospheres such as for Earth, Titan
and Mars, but the main features derived (comparison with pure
gravity case, asymmetries, exopause) are general results that
may be applied to any planet hosting a dense atmosphere influ-
enced by a radiation pressure force.

4.1. Asymmetries and comparison with Chamberlain profiles

First of all, as shown on Fig. 3 (lower panels), close to the planet,
our densities (which include the radiation pressure effect) over-
come the ballistic densities from the Chamberlain (1963) model.
In the nightside direction at Earth, the densities are up to
10<sup>0.6</sup> ≈ 6 times higher than predicted by Chamberlain (1963).
However, in the dayside direction at Earth (see Fig. 8), the densities
are less high, with an enhancement factor of about 10<sup>0.4</sup> ≈ 2.5
times. This is also the case at Titan and Mars (see Fig. 8). We will
also later see (Fig. 7) that, even for relative small radiation pres-
sures (i.e. r<sub>exo</sub> < 10 − 20 R<sub>pressure</sub>), our density overcomes the one
obtained through the Chamberlain (1963) formalism. A physical
explanation will be brought in Section 4.3.

In Fig. 3 that shows the 2D maps of ballistic particles densities at Earth, we see clearly the asymmetries induced by the radiation pressure. As expected and previously mentioned, the ballistic density is clearly more significant in the nightside direction than in the dayside direction. Most of the particles are preferentially blown behind the planet. This observation is in agreement with the work by Bishop and Chamberlain (1989), who also found such an asymmetry corresponding to the well known geotail phenomenon of enhanced nightside densities observed (Thomas and Bohlin, 1972; Bertaux and Blamont, 1973; Zoennchen et al., 2011; Bailey, 2012). Nevertheless, the densities in the dayside direction are not so low: the density profile in the dayside direction remains significant compared with the transverse direction, the Dusk/Dawn/ North Pole/South Pole plane. For example, at 10 Earth radii the densities in the dayside direction are about 25% larger than in the perpendicular plane (and 50–100% larger in the nightside direction). The exosphere seems to have a prolate shape oriented along the Sun-planet axis.

The exospheric density asymmetries were recently observed
thanks to the TWINS mission and published by Zoennchen et al.
(2011, 2013) and Bailey (2012). In Fig. 4, we compare at constant
distance (8 R<sub>E</sub>) the variability of the Earth exospheric density from
our model and from the observations as a function of the solar
zenith angle. As Zoennchen et al. (2011, 2013) and Bailey (2012)
fit their observations with a limit number of spherical harmonics,
we do the same here for a clear comparison. Zoennchen et al.
(2013) observed a decrease of 45% in the dusk and dawn directions
in the equatorial plane and an enhancement by up to 45% in the
nighside direction compared with the dayside. In Fig. 4, the
decrease for our model is about 36% and the increase about 25%.
The lower difference in the nightside direction can be explained
by the lack of escaping particles in our model, which would cer-
tainly increase the density in the nightside direction.

Our model thus reproduces quite well the local time variation in
the equatorial plane although we do not know the precise condi-
tions in the exosphere during the observations by Zoennchen
et al. (2013), i.e. the exact exospheric temperature and radiation
pressure acceleration (which is known to have a strong variability
(Vidal-Madjar, 1975)). Nevertheless, we have significant discrep-
ancies (compared with Zoennchen et al., 2013) in the meridional
plane in particular during the equinox. One explanation can be
the satellite particles as will be discussed in Section 4.3.
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Algorithm to determine the value, 0 or 1, of the function $l_{ballistic}$ according to the initial conditions. Once this value is known, we use Eq. (24).</figcaption></figure>

## 4.2. The exopause

We also plotted the ballistic density profiles at Earth as a function of the distance for different $\theta$ (angle with the Sun-planet axis) in Fig. 5. Beyond the clear local time asymmetries and enhanced densities compared with the Chamberlain (pure gravity) profiles, we observe a cut-off of the ballistic density at 36 Earth radii. This limit corresponds to an exopause located at $R_{pressure} = \sqrt{GM/a} = 36 R_E$, which was not introduced artificially as was previously done by the previous works (e.g. Bishop, 1991). Unfortunately, we could not demonstrate this property analytically but the modeling proves that the bounded trajectories occur only inside the sphere of radius $R_{pressure}$. With the Chamberlain (1963) approach, the limit for bounded trajectories is infinity. Physically,

a first limit is however imposed by the presence of the Sun (or by nearby planets for satellites): the Hill's sphere, that defines the limit of the gravitational influence of the central body. But here, the radiation pressure induces another limit (below the Hill sphere radius) located at $R_{pressure}$, i.e. where the accelerations due to the gravity and radiation pressure forces are equal. Beyond this limit, all particles are escaping. In future works, we will show its implication on the evolution of planetary atmospheres.

The sharp density drop (shown for ballistic populations) at the distance $R_{pressure}$ should still be seen after adding the satellite and escaping populations, and lead to a sharp drop in the total density profiles. The densities of satellite and escaping populations are indeed decreasing with distance: the escaping component should approximately decrease in $r^{-2}$, whereas the satellite component
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Maps of ballistic particles densities at Earth, Mars and Titan (normalized by the exobase density) with an exospheric temperature of respectively $T_{\text{exo}}$ = 800 K, 200 K and 152 K (from the left to the right). The upper panels represent the output of our model. The middle panels represent with the same colors the ballistic density map provided by the Chamberlain (1963) model. The lower panels represent the ratio between the two models in log-scale: red areas show where our density exceeds the Chamberlain's one, the blue areas represent where the Chamberlain's model overestimates the density (except for Titan, our densities overwhelm in any part of the exosphere). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

will also disappear at $R_{\text{pressure}}$ since no bound motion can exist beyond this limit (see also Beth et al., 2014). Thus, no strong increase for these two populations will be able to counter the sharp drop of the ballistic component, which should be seen in the total density measured (unless an external source of neutrals, e.g. the neutral part of the solar wind, adds a significant density that hides the sharp drop).

As far as we know, the only relevant neutral density measurement allowing to investigate this topic is the study by Brandt et al. (2012). These authors provided energetic neutral atom (ENA) images of the Titan environment, and showed that ENA fluxes may be observed up to 50,000 km, i.e. the Hill sphere radius determined by the gravitational influence of Saturn. The ENAs are produced by charge exchange reactions between the magnetospheric ions and the exospheric neutrals, so that the ENA flux profiles provide information about the neutral density profiles. At Titan, the location of the exopause, where a sharp density drop is expected, is due to the gravitational influence of Saturn rather than

to the radiation pressure that leads to an $R_{\text{pressure}}$ distance much further. The observation of ENA fluxes up to the exopause is thus in agreement with a sharp neutral density drop at the exopause distance expected by our model.

### 4.3. The satellite particles

Fig. 5 also provides a comparison between our modeled ballistic profiles and the ballistic and ballistic + satellite densities from Chamberlain (1963). Near the planet, our ballistic density profile remains for any direction between these last two profiles. Why?

The radiation pressure disturbs the trajectories of particles that are initially conics. As shown by Bishop and Chamberlain (1989), the particles crossing the Sun-planet axis ($P_{\phi} = 0$) cannot be satellite particles: the bounded trajectories see their periapsis decrease with time and they necessarily cross the exobase if $P_{\phi} = 0$. The radiation pressure separates the areas where ballistic and satellite
particles can exist: the ballistic particles are preferentially near the Sun-planet axis whereas the satellite particles are mostly located far from this axis in the perpendicular plane (see below). The radiation pressure will convert satellite trajectories into ballistic trajectories and these last ones into escaping particles. This is why we have our density between ballistic and ballistic plus satellite particles densities from Chamberlain (1963): a part of our ballistic particles were probably initially satellite particles converted into ballistic ones by the radiation pressure. Moreover, such a conversion of trajectories is stronger near the axis than in the perpendicular plane. To support this explanation, we cannot calculate the satellite particles density since the Liouville theorem cannot be applied to particles which do not cross the exobase. Nevertheless, we can estimate the phase space volume dedicated to the satellite particles, and look where it is maximum. The phase space volume is given by:

<figure><img src="image_8.png" /><figcaption>Fig. 4. Comparison at Earth between our ballistic particles density and the observations made by Zoennchen et al. (2013) at the same distance (8 R<sub>E</sub>), as a function of the Solar Zenith Angle (SZA) in degrees. We fitted our model (red) with four spherical harmonics (green) as done by Zoennchen et al. (2013) (blue, with only two sinusoids and the constant). For a better comparison, we scaled the blue and green plots in order to be at 1 for SZA = 0. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. Radial density profiles for ballistic particles at Earth with T<sub>exo</sub> = 800 K. The specific directions are in black for the nightside, in green for the dayside and in cyan for the perpendicular plane, whereas a set of profiles for intermediate directions are plotted in red. The corresponding Chamberlain (1963) ballistic and ballistic plus satellite particles densities are given in solid and dashed blue. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

$$V = \int_{R_E^3} \mathcal{I}_{\text{type}} d^3 \vec{p}$$

The bounds of integration (or the value of $\mathcal{I}_{\text{type}}$) depend on the type of particles considered (type = ballistic or satellite; for escaping particles, V is infinite). This volume has no physical meaning but provides the available volume. However, one can extract a useful information if we compare the phase space volumes for satellite and ballistic particles (given in Fig. 6). The satellite particles are thus preferentially located far from the Sun-planet axis and in the perpendicular plane.

The existence and production/loss of satellite particles is however not well-known. A recent work by Beth et al. (2014) showed that the satellite particles are produced by scarce collisions just above the exosphere. In the Earth case, they showed the satellites particles do not contribute significantly to exospheric densities. Nevertheless, we mentioned above a discrepancy between our model and the observations by Zoennchen et al. (2013) in the perpendicular plane during equinox, that could be related to the presence of satellite particles. During this period, the polar cusps are in the perpendicular plane. The production of satellite particles might be in the polar cusps where the densities are large and where the satellite trajectories are more stable.

A last evidence for the conversion between ballistic and satellite particles is given by Fig. 7, where the Earth ballistic density profile at 8 R<sub>E</sub> (as in Fig. 4) is given using a fictitious radiation pressure value divided by 1000. We would expect a density profile similar to the pure gravity case, i.e. the Chamberlain profile of ballistic particles. This is however not the case. The radiation pressure force disturbs the trajectories regardless of the value considered: the particles that are quasi-satellite see their periapsis altitude decreases more slowly if the radiation pressure is smaller, so that they can turn around the planet during a longer time without crossing the exobase than with a larger radiation pressure, but they will anyway cross it after some time. In this work, there is no distinction between these quasi-satellite particles and ballistic particles since no simple physical parameter allows to distinguish them.

According to the densities profiles, given in Figs. 5 and 7, we can reasonably assume that for small distances compared with R<sub>pressure</sub> (generally some planetary radii except for Hot Jupiters for where

<figure><img src="image_10.png" /><figcaption>Fig. 6. Ratio between the phase space volumes for satellite particles and for ballistic particles at Earth with T<sub>exo</sub> = 800 K. In red, the ratio is around 1 and in blue the ratio is around 0. Redder areas are most probable areas where satellite particles can exist. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
the exopause may be located below the exobase), the densities are between $n_{bal}$ and $n_{bal} + n_{sat}$ from the Chamberlain formalism. Even if the collisions are not included in the model (which are the source for satellite particles), assuming Chamberlain density profiles with and without the satellite particles contribution can give a range for exospheric densities at a given distance for the disturbed case by the radiation pressure.

## 4.4. Comparative planetary science

With our studies of different cases, Earth, Mars and Titan, corresponding to different conditions (i.e. exospheric temperatures and radiation pressures), we are able to explore a wide range of configurations and to derive general conclusions.

All our case studies showed, with an exopause (always located at $R_{pressure}$) above the exobase, an enhancement of the ballistic densities compared with those provided by the Chamberlain (1963) model. As shown by Figs. 3 and 8, the radiation pressure indeed strongly affects the density profiles. Fig. 8 reveals a ratio between the disturbed and non disturbed densities that is very similar for

Earth and Mars (up to 3), but larger at Titan (up to 7), with a peak at the same normalized distance. The variability of the induced turbance between the various planets may be explained as follows. The first effect of the radiation pressure is to break the spherical symmetry: the physical consequence is the conversion of satellite particles into ballistic ones as shown by Fig. 7. The secondary consequence depends on the intensity of the radiation pressure: the stronger it is (i.e. the closer to the Sun the planet is), the more ballistic particles become escaping ones. We should thus expect less ballistic particles at Earth than at Mars. However, the equations driving the density profiles mostly depend on the parameters $r_{exo}$ and $R_{pressure}$ (i.e. $\lambda_c$ and $\lambda_a$), in particular on the ratio $r_{exo}/R_{pressure}$. Earth and Mars actually have similar ratios $r_{exo}/R_{pressure}$ (i.e. about 35) and thus have a similar enhancement as a function of $r/R_{pressure}$ (cf. Fig. 8). For the Titan case, the radiation is hundred times weaker, which would suggest a weak density enhancement, but the ratio $r_{exo}/R_{pressure}$ is completely different: $R_{pressure}$ is 100 times larger than $r_{exo}$. Nevertheless, we should take precautions the Titan case because another external force could strongly affect the dynamic of atmospheric species: the gravitational attraction by Saturn.

<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Same plot as for Fig. 4 with a radiation pressure thousand times smaller. We superimposed the densities for ballistic and ballistic plus satellite particles provided by Chamberlain (1963).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Ratio between our ballistic density profile in the dayside direction and the ballistic density profile from Chamberlain (1963) at Earth (blue), Mars (red) and Titan (orange). The two lines for Earth and Mars correspond to two exospheric temperatures: 800 K and 1200 K at Earth, 200 K and 350 K at Mars. For Titan, the exospheric temperature is 152 K. The distance *r* is normalized with the value of the exopause distance (i.e. $R_{pressure}$) for each case. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 5. Conclusions

In this paper, we generalize the initial work by Bishop and Chamberlain (1989) by developing a 2D model (3D if we do not assume the axisymmetry) for the density profiles of ballistic exospheric neutral particles, in order to study the impact of the radiation pressure on the structure of planetary exospheres such as at Earth, Mars or any planet with a dense atmosphere.

We reproduce quite well with our simulations the different exospheric asymmetries observed at Earth:

* the “tail phenomenon”: the Earth exosphere has higher densities for atomic Hydrogen in the nightside direction than in the dayside direction. This is already known (Thomas and Bohlin, 1972; Bertaux and Blamont, 1973) and directly attributed to the radiation pressure.

* dusk/dawn/North Pole/South Pole asymmetries (Bailey and Gruntman, 2011; Zoennchen et al., 2011, 2013): the radiation pressure induces a depletion of particles in the perpendicular plane, observed in the equatorial plane.

Moreover, the radiation pressure entails an increase of ballistic particles densities which are in the lower corona (up to several planetary radii), the main exospheric component compared with satellite and escaping particles. Compared with the Chamberlain (1963) model (i.e. without the radiation pressure), the densities are several times higher (up to 2.5 in the dayside direction and 4 in the nightside one at Earth). Only the ballistic (i.e. not the escaping) particle density calculation is performed for numerical reasons (time and precision issues).

We highlight also the appearance of a characteristic distance of the exosphere: the “exopause”. This concept was introduced by Bishop (1991), who included it artificially in their model: this boundary is the limit where the intensity of the radiation pressure is equal to the gravitational attraction. As shown in this paper, this limit appears naturally in our simulations with a break in our density profiles. Physically, the exopause divides the exosphere into two regions: below the exopause, we can find bounded (satellite and ballistic particles) and unbounded (escaping) trajectories; above the exopause, we find only unbounded trajectories. The exopause will lead to a local sharp drop for the total density (including ballistic, satellite and escaping populations), that is in agreement with the observation of energetic neutral atom fluxes at Titan up
the Hill sphere radius only (i.e. the exopause for Titan) by Brandt
et al. (2012). Our model provides only density profiles for a plane-
tary exosphere where the exopause is located above the exobase (a
future paper will investigate the extreme case of Hot Jupiters
where the strong radiation pressure pushes the exopause down
to the exobase). The exopause boundary also induces a constraint
for modeling the exospheres: the size of the simulations box must
be large enough to contain the exopause in order to take into
account all the asymmetries induced by this force.

Moreover, in our study, we have also shown the influence of the
radiation pressure on the repartition of ballistic and satellite pop-
ulations. On the one hand, near the Sun-planet axis, the periapsis of
the bounded particles decreases slowly because of the radiation
pressure until it crosses the exobase. Thus, we find essentially bal-
listic particles (called also satellite particles with finite lifetime)
near the Sun-planet axis. On the other hand, the available regions
to find bounded trajectories which do not cross the exobase (i.e.
satellite particles) are essentially in the dusk/dawn/North Pole/
South Pole plane. Thus, the radiation pressure separates clearly
the ballistic and satellite particles regions even if these particles
are both bounded.

Finally, the radiation pressure induces an important effect on
the velocity phase space: the radiation pressure converts a part
of the satellite particles into ballistic particles and these ones into
escaping particles. This explains the enhancement of ballistic den-
sities compared with the Chamberlain (1963) model: the radiation
pressure converts efficiently the satellite particles of his model
(without radiation pressure) into ballistic ones for ours (with radi-
ation pressure). This explanation is supported by Figs. 5 and 7: near
the planet or for small radiation pressure accelerations, the densi-
ties remain between the ballistic and ballistic plus satellites ones
provided by the Chamberlain (1963) formalism. Consequently,
including or not the Chamberlain (1963) satellite particles parti-
tion function provides an appropriate range of densities to include
the radiation pressure effect on exospheric density profiles, even in
the absence of collisions that are the source of satellite particles.

In future works, we will study the photogravitational Circular Three-Body Problem and the implication on the stability of planetary exospheres, in particular for Hot Jupiters. Moreover, the escaping particles density could not be calculated here because of numerical issues but the escaping flux will be investigated in details in a future work.

Acknowledgment

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

**References**

Bailey, J.J., 2012. Three-Dimensional Exospheric Hydrogen Atom Distributions Obtained From Observations of the Geocorona in Lyman-α, Ph.D. thesis. University of Southern California.

Bailey, J., Gruntman, M., 2011. Experimental study of exospheric hydrogen atom distributions by Lyman-α detectors on the TWINS mission. J. Geophys. Res.: Space Phys. 116, 9302.

Banks, P.M., Kockarts, G., 1973. Aeronomy. Academic Press, New York.

Bertaux, J.L., Blamont, J.E., 1973. Interpretation of OGO-5 Lyman-α measurements in the upper geocorona. J. Geophys. Res. 78, 80-91.

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 et al., 2016. Theory for planetary exospheres: I. Radiation pressure effect on dynamical trajectories. Icarus 266, 410-422.

Biscani, F., Izzo, D., 2014. The stark problem in the Weierstrassian formalism. Month. Not. Roy. 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.

Brandt, P.C. et al., 2012. The distribution of Titan's high-altitude (out to ~50,000 km) exosphere from energetic neutral atom (ENA) measurements by Cassini/INCA. Planet. Space Sci. 60, 107-114.

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.

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

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

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.

Thomas, G.E., Bohlin, R.C., 1972. Lyman-α measurements of neutral hydrogen in the outer geocorona and in interplanetary space. J. Geophys. Res. 77, 2752-2761.

Vidal-Madjar, A., 1975. Evolution of the solar Lyman-α flux during four consecutive years. Sol. Phys. 40, 69-86.

Zoennchen, J.H. et al., 2011. The TWINS exospheric neutral H-density distribution under solar minimum conditions. Ann. Geophys. 29, 2211-2217.

Zoennchen, J.H., Nass, U., Fahr, H.J., 2013. Exospheric hydrogen density distributions for equinox and summer solstice observed with TWINS1/2 during solar minimum. Ann. Geophys. 31, 513-527.