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

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.


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 1 INTRODUCTION 2 as maxwellian anymore 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).
-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 Liouville's theorem. The density for each population is then derived as the product between the barometric law and a partition function ζ. n(r) = n bar ζ(λ) (1) = n(r exo )e λ−λexo (ζ bal + ζ esc ) where λ is the ratio between the gravitational and thermal energies.
with r the distance from the center of the body, v 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 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 twobody 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 φ , 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.(2015a), 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.

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-α flux f 0 , in 10 11 photons.cm −2 .s −1 .Å −1 : 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: where x is the sunward coordinate and θ 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: independent of t and φ. p u , p w and p φ 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: We do not assume p φ = 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: H, p φ and A defined as As function of these three constants, we can rewrite the conjugate momenta: We can remark here that for p φ = 0, 0 is a simple root of P 3 and Q 3 .
Based on the formalism detailed in Beth et al. (2015a), we use dimensionless quantities and the same annotations: 6 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 polynom, one root can be 0 if P φ = 0. We define λ a as: the Jeans parameter at the distance R pressure = GM/a.

Calculation of exospheric densities
The density in a gas is given from the distribution function f by: The bounds of integration depend on the type of particles. 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 twice the exobase, the satellite particles do not cross the exobase and do not escape, and the escaping particles cross the exobase once before they escape.
Necessary and sufficient Yes Yes No, need for tracking conditions?  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 φ too (cf. eq 8) 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 According to the Poincaré recurrence theorem as used by Beth et al. (2015a), 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 φ = 0), all bounded trajectories cross the exobase. These authors considered the particles orbiting severa 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 φ can vary. Although our model cannot take into account the singularity P φ = 0, this is not a problem for the integration (for the numerical evaluation, we just avoid this value).

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 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. (2015a)). 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 E according to equation 12: As explained by Beth et al. (2015a), the U-motion does not have the same period as the W -motion or the φ-motion. In the (U, W )-plane, the trajectory fills the whole square  Beth et al. (2015a), 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 exo ). We will thus be able to call it a ballistic particle. By definition u + w = (U + W )R pressure = 2r, and the nearest distance from the center of the planet where the particle can pass is 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 equation 14 is the part of the velocity space where W − and W + are real positive and 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 disitribution function is maxwellian at the exobase.

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 (E, A, P φ ).
Third step For each tuple (E, A, P φ ), 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 ζ bal (r, θ) that is the generalization, when the radiation pressure effect is included, of the 1D partition function defined by Chamberlain (1963): with n exo the exobase density, θ the angle with the respect to the Sun-planet axis and r the distance from the planet. By definition, for a Maxwellian distribution, We choose to work with dimensionless quantities and the (U, W ) coordinates. Operating some transformations and with the correct Jacobian, we obtain: The formula 18 is only available for U = 0 or W = 0. Indeed, if U = 0 or W = 0, then P φ = 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, E is negative. It represents the contribution of the potential energy and the kinetic energy. Thus, we have the inequality from the equation 10: The modulus inside the exponential takes finite values. In this case, we choose the following change of coordinates: To estimate this integral 21, we use the Gauss-Legendre quadrature as: N is the number of points used for the integration and 1 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. figure 2).
The Gauss-Legendre method is particularly efficient in this case because all bounds of integration are finite: R is finite (see eq. 19), Θ belongs to [0; π] and Φ to [0; 2π].

The escaping and the satellite particles
Our approach can unfortunately not be applied so easily to calculate 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 E). The Poincaré recurrence theorem 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 particles are.

Results
Bishop and Chamberlain (1989) provided only the ballistic density 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 exosphere) 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 influenced by a radiation pressure force.

Asymmetries and comparison with Chamberlain profiles
First of all, as shown on the figure 3 (lower panels), close to the planet, our densities (which include the radiation pressure effect) overcome the ballistic densities from the Chamberlain (1963) model. In the nightside direction at Earth, the densities are up to 10 0.6 ≈ 6 times higher than predicted by Chamberlain (1963). However, in the dayside direction at Earth (see figure   8), the densities are less high, with an enhancement factor of about 10 0.4 ≈ 2.5 times. This is also the case at Titan and Mars (see figure 8). We will also later see ( fig. 7) that, even for relative small radiation pressures (i.e. r exo < 10 − 20 R pressure ), our density overcomes the one obtained through the Chamberlain (1963) formalism. A physical explanation will be brought in section 4.3.
In the figure 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. (2013) at the same distance (8 R E ), as a function of the Solar Zenithal 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.
Our model thus reproduces quite well the local time variation in the equatorial plane although we do not know the precise conditions 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 discrepancies (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.

The exopause
We also plotted the ballistic density profiles at Earth as a function of the distance for different 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 = 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 approach, the limit for bounded trajectories is the 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 will also disappear at Rpressure 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 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.

The satellite particles
The figure 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 φ = 0) can not be satellite particles: the bounded trajectories see their periapsis decrease with time and they necessarily cross the exobase if P φ = 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 can not 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: The bounds of integration (or the value of 1 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 the figure 7, where the Earth ballistic density profile at 8 R E (as in figure 4) is given using a fictive 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 figures 5 and 7, we can reasonably assume that for small distances compared with R pressure (generally some planetary radii except for Hot Jupiters for where 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.

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 the figures 3 and 8, the radiation pressure indeed strongly affects the density profiles. The figure 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 disturbance 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 the figure 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. λ c and λ 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. figure 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.

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 asymetries 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) 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). We precise that our model provides only density profiles for a planetary 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 populations. 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 ballistic 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 densities 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 radiation pressure). This explanation is supported by the figures 5 and 7: near the planet or for small radiation pressure accelerations, the densities remain between the ballistic and ballistic plus satellites ones provided by the Chamberlain (1963) formalism. Consequently, including or not the Chamberlain satellite particles partition 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 photogravitationnal 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.