Developed by Chamberlain (1963), this model provides insights into the density distribution in planetary exospheres influenced by planetary gravity. TheThe model assumes a Maxwellian velocity distribution at the base of a spherically symmetric exosphere, where particles are injected. he formalism of the model is discussed in detail, e.g., Shen (1963), Chamberlain (1963), Stern and Flynn (1995), Lammer and Bauer (1997), Valeille (2009), and Öpik and Singer (1961). Despite its limitations, the Chamberlain model remains instrumental in studying the evolution of planetary atmospheres.

Assuming that (1) the exosphere is spherically symmetric, and (2) at the location of the exobase, which serves as the lower boundary for the model, the atmosphere is in equilibrium, the density in the entire exosphere n a is in Equations 1, 2. n a r = n 0 e − 1 − α E − 1 − α 2 1 / 2 e − E / 1 + α − n 0 α 2 2 π 1 / 2 1 + E e − E E 1 2 H E − 1 + 1 − 2 3 H E α 1 / 2 , (1) H E = 3 2 E + 1 1 + π e E α E 1 / 2 1 − erf E ≈ 1 + 2 E  for large  E . (2) Here, n a is the density of a collisionless exosphere, n 0 is the density at the exobase, E = G M m R k T , α = r c r , where r c is the location of the exobase, and T is the temperature at the exobase.

A traditional definition of the exobase’s location for an exosphere (Equation 3) is where the collision mean free path, λ = σ n ( r c ) , equals the scale height, H s h , ∫ r c ∞ σ n r d r = σ n r c H s h r c = H s h r c λ r c = 1 , (3) where H s h ( r c ) = k T 0 m g 0 , and σ is a total collision cross-section (Schunk and Nagy, 2000).

Most radiative emissions of gaseous species in cometary comae, observable in visible light, are primarily due to fluorescence induced by sunlight within an optically thin medium. Quantitative analyses of gas species within these comae revealed that the spatial variation in the brightness of certain species (e.g., C 2 , CN) is inconsistent with the inverse distance relationship. Further investigations have recognized that most gas species detected through visible-range spectroscopy are not stable molecules but relatively unstable radicals produced by the photodissociation of parent molecules.

The quantitative framework for analyzing these phenomena, developed by Haser (1957), remains widely utilized in comet research today. Haser’s model characterizes the distribution of secondary species, such as cometary radicals, resulting from the photodissociation of parent species. The model assumes that the coma acts as a spherically symmetric point source from which parent molecules uniformly outflow. The depletion of these species’ densities, n p ( r ) , is modeled using exponential decay: n p r = Q 4 π r 2 v e − r / γ p , (4) where Q is the production rate of the comet, r is the distance from the nucleus, v is the radial outflow speed, γ p = v / ν p h is the parent species scale lengths and ν p h is the photo-destruction frequency. The density, n d , of a daughter species created by the destruction of its parent species is given in Equation 5. n d r = Q 4 π r 2 v γ d γ p − γ d e − r / γ p − e − r / γ d (5) Here, γ p is the daughter species scale lengths.

The processes occurring within the coma of a comet are generally more complex. Numerous daughter species detected do not originate from a singular, direct decomposition of a parent molecule. Instead, a daughter molecule may be produced from different parent molecules via diverse reaction pathways and may involve several intermediate stages. Moreover, if the dissociation processes are exothermic, the excess energy may enhance the reaction products’ kinetic energy, increasing the velocity of the lighter fragments (Shou, 2017; Combi et al., 2004). Despite its simplifications, the Haser model is extensively utilized as it provides a robust estimation of production rates by employing scale lengths readily available in the literature.

The 2005 discovery by the Cassini spacecraft of gas and ice grain plumes emanating from Enceladus’ south pole has ignited significant interest in and research into the mechanisms behind these phenomena (Waite et al., 2006). The multiplume model developed to study the Enceladus’ plume integrates contributions from discrete gas vents, diffuse fissures or crevices along the Tiger Stripes - large, parallel cracks near the south pole that emit water vapor and ice particles and a global spherical source. The gas vents are treated as individual point sources where gas is ejected into the exosphere at high velocities (approximately 1 km/s), significantly exceeding Enceladus’ escape velocity (235 m/s) (Tenishev et al., 2014). Therefore, a number density of water vapor in the exosphere that is due to an individual gas vent can be evaluated as n x = A p r 2 ∫ v 2 ⁡ exp − β 2 v 2 ⁡ sin 2 ⁡ θ + v ⁡ cos ⁡ θ − v p 2 d v if x − r p ⋅ r p ≥ 0 0 otherwise, (6) where x is a point where a number density is calculated, ( r , ϕ , θ ) are coordinates of the point x in a spherical coordinate frame associated with the vent, r p is the vent location, v p is the bulk speed of the gas flow ejected from the vent, β = m H 2 O 2 k T , and A p is a normalization constant defined as A p = F p ∭ v 3 ⁡ exp − β 2 v 2 ⁡ sin 2 ⁡ θ + v ⁡ cos ⁡ θ − v p 2 d v d θ d ϕ (7) Here, F p is the total source rate ( s − 1 ) of the vent. The integration limits in Equation 7 are constrained by the condition ( v ( v , ϕ , θ ) ⋅ r p ) ≥ 0 . Here, v is the velocity vector that corresponds to that defined in the spherical frame by the combinations of ( v , ϕ , θ ) . The integrals in Equations 6, 7 are calculated numerically.

In most cases of practical interest, studying exospheres of planets, planetary satellites, and comets involves rarefied gas flows under strong non-equilibrium conditions. These conditions can be described only with the Boltzmann equation, which, in its general form, can be written as ∂ f 1 ∂ t + v 1 ∂ f 1 ∂ x + a 1 ∂ f 1 ∂ v 1 = ∫ | v 1 − v 2 | f 1 v 1 ′ f 1 v 2 ′ − f 1 v 1 f 1 v 2 d σ d v 2 , (8) where σ is the scattering cross-section. The integral in the right-hand-side of Equation 8 describes collisions, by which collision partners having velocities v 1 ′ and v 2 ′ get the velocities v 1 and v 2 after the collision. Equation 8 can be generalized for gas mixtures and include the effect of external forces.

The Direct Simulation Monte Carlo (DSMC) method is a method of choice to solve the Boltzmann equation for gas flow in the exospheres of planets, moons, and comets (Combi, 1996; Tenishev et al., 2008; Bird, 1994; Crifo et al., 2002; Tucker et al., 2015). One of the DSMC method’s most important advantages is its ability to incorporate processes that are more complex than elastic collisions without significantly complicating the numerical procedure.

The general scheme of Monte Carlo models can be described using so-called Markov chains. Briefly, a Markov chain is defined as a system, S , consisting of a finite set of states M { s 1 , s 2 , … , s l } . At each discrete sequence of times t = 0,1,2 , … , n , the system is in one of the states, s i , which determines a set of conditional probabilities p i 1 , p i 2 , … , p i l . The quantity p i j is the probability that the system, which is in state s i at the n -th time step, will be in state s j at the ( n + 1 ) -th time step. In other words, p i j is the probability of the transition s i → s j . It is important to note that the probability of a transition depends only on the current state and is not affected by the previous history.

The evolution of the distribution function f ( v , t ) as a Markov process is described by a collision integral f v , t = ∫ f v − Δ v P v − Δ v , Δ v d Δ v , (9) where P ( v , Δ v ) represents the probability that a particle with velocity v at time t will have velocity v + Δ v at time t + Δ t , and it must satisfy the normalization condition ∫ P v , Δ v d Δ v = 1 . (10) A reasonable model for the transition probabilities must be developed to apply Equations 9, 10 to a real gas. This formulation does not require a simultaneous change in the velocity coordinates of both partners during a collision, allowing for the description of more comprehensive relaxation processes compared to the standard Boltzmann collision integral. In most practical cases, models of microscopic processes that define the transition probabilities are available for rarefied gases.

The results of a Monte Carlo numerical simulation are moments of the velocity distribution function for the simulated gas flow, representing measurable parameters such as density, bulk velocity, temperature, and pressure. The major challenge is to develop an approximation of a complex Markov chain representing the dynamics of each atom or molecule of simulated gas flow using a simpler one, ensuring that the mean value of the distribution function moments remains the same in both cases. One technique to reduce the number of possible states in the chain is to decrease the total number of degrees of freedom in the system based on physical considerations. This approach is used in the DSMC method, where a single model particle represents many real gas molecules.

The numerical schemes of the DSMC method are based on physical assumptions that form the basis of the phenomenological derivation of the Boltzmann equation (Bird, 1994). The key concept in developing collision relaxation schemes is the total collision frequency, ν . Using a probability density, ω , of transition ( v i , v j ) → ( v i ′ , v j ′ ) for a pair of particles, the collision frequency can be defined as ν = n N ∑ i < j ∫ ω v i , v j → v i ′ , v j ′ d 3 v i d 3 v j = n N ∑ i < j σ t g i j g i j , (11) where g i j is the relative speed between particles i and j , and σ t ( g i j ) is the total collision cross section. The total collision frequency given by Equation 11 depends on the velocities of all particles. In principle, it should be recalculated after each collision, which is a time-consuming procedure, as the summation is performed over N ( N − 1 ) / 2 possible collision pairs, where N is the number of model particles in a computational cell. To achieve the correct relaxation dynamics in a gas flow, the characteristic size of computational cells must be smaller than the local mean free path.

Due to the statistical nature of the DSMC method, numerical solutions always exhibit some noise. There are two principal sources of error associated with DSMC calculations. One arises from a high real-to-simulated particle number ratio, which becomes especially significant for high-density flows. Another source of statistical noise is notable when the mean flow velocity is much smaller than the mean molecular thermal speed. For low-velocity flows, large statistical fluctuations can obscure some features of the flow structure. Noise filtering techniques have been described by, e.g., Boyd and Stark (1989), Kaplan and Oran (2002), and Garcia et al. (1987).

Significant efforts have been made to establish a theoretical proof of the convergence of numerical solutions obtained using the DSMC method to those of the Boltzmann equation. Notable results from these efforts can be found by, e.g., Babovsky and Illner (1989), Rjasanow and Wagner (1996), Wagner (1992), and Volkov (2011).

Fluid models offer a practical and efficient approach to studying the dynamics of planetary exospheres. Despite kinetic physics dominating the dynamics of exospheres, fluid-type numerical methods are also used to model these environments. One significant advantage of fluid-type models is their ability to be easily coupled with magnetohydrodynamic (MHD) models to simulate complex interactions between the exosphere and the surrounding plasma environment.

The theoretical foundation of fluid models is based on solving the Euler equations for each species: ∂ ρ n ∂ t + ∇ ⋅ ρ n u n = δ ρ n δ t ∂ ρ n u n ∂ t + ∇ ⋅ ρ n u n u n + p n I = δ ρ n u n δ t ∂ p n ∂ t + ∇ ⋅ p n u n + γ n − 1 p n ∇ ⋅ u n = δ p n δ t (17) Here, ρ , p , u , γ , and I denote the mass density, pressure, velocity vector, the specific heat ratio, and the identity matrix, respectively. The right-hand side of Equation 17 represents the source terms that describe the interaction of each simulated species with others: δ ρ s δ t = ∑ n = neutrals ν n → s n n m s − ∑ t = other species ν s → t n s m s (18) δ ρ s u s δ t = ∑ n = neutrals ν n → s n n u n − u s m s + δ ρ s δ t u s + ∑ t = other species ν ̄ s , t n s m s u t − u s (19) ν ̄ s , t = m t m t + m s n t σ s , t 8 k π T t m t + T s m s + u t − u n 2 (20) δ p s δ t = γ s − 1 m s m n ν n → s n n Δ E n → s + 1 2 ν n → s p n + 1 3 ν n → s n n u s − u t 2 m s − ∑ t = other species ν s → t p s + ∑ t = other species ν ̄ s , t m t m s m t + m s n s 2 3 u s − u t 2 + 2 k T t − T s m t (21) The source terms are described in Equations 18–21. Those source terms with chemical reaction frequencies ν s → t are related to photochemical reactions. In the pressure source term, the excess energies are partitioned under the condition that the momenta of the daughter species must be conserved and are, therefore, inversely proportional to the mass. The source terms involving momentum transfer coefficients ν ̄ s , t between species s and t account for collisions. The collision frequency is linearly proportional to the density of each involved gas species and the relative speed of the colliding gas molecules or atoms, calculated from both species’ thermal and bulk velocities.

A detailed discussion of the methodological questions regarding the use of fluid methods in modeling exospheres is presented in, e.g., Rubin et al. (2014a,b), Huang et al. (2016), Shou et al. (2016), and Shou (2017).

The gas and ice grain plumes discovered by Cassini in 2005 have made Enceladus an object of increased scientific interest (Waite et al., 2006). The observed activity in the south pole region has raised questions regarding the source of the plumes’ material, the mechanisms of its delivery to the surface, and the distribution of dust, ice grains, and gas in the exosphere (Spencer, 2013; Hedman et al., 2009; Postberg et al., 2011; Ingersoll and Ewald, 2011; Hansen et al., 2011; Spitale and Porco, 2007).

The multiplume model (see Section 2.3) was utilized to investigate the potential impact of diffuse or multiple small gas sources along the Tiger Stripes on the water vapor distribution in Enceladus’ exosphere (Tenishev et al., 2010). The model accounts for gas production from individual vents, a global spherical source, and Tiger Stripes to study their effect on water vapor distribution in Enceladus’ exosphere. In addition, DSMC models have been used to simulate Enceladus’ plume, considering similar physics of particle ejection into the exosphere (e.g., Tucker et al., 2015; Mahieux et al., 2019; Yeoh et al., 2015).

The Tiger Stripes are simulated using vertically directed point sources, described by Equation 6, which are uniformly distributed along the stripes. Additionally, a background gas distribution is incorporated to account for the global exospheric environment, modeled as n = c 0 + c 1 r α , where c 0 and c 1 are constants and α = − 2 . This background population represents molecules sputtered from Enceladus’ surface.

The results of applying the multiplume model to analyze Cassini’s Ion and Mass Spectrometer (INMS) E3 and E5 flybys, as well as Ultraviolet Imaging Spectrograph (UVIS) observations from 2005, 2007, and 2010 (Teolis et al., 2010; Hansen et al., 2006; Hansen et al., 2008; Hansen et al., 2011), are illustrated in Figure 1.

The best fit was achieved with a fissure temperature of 180 K and bulk velocities of 350 m/s and 950 m/s for the regular and extra fissures, respectively. These values are consistent with Cassini/CIRS measurements of the surface temperature and the gas thermal velocity (Spencer et al., 2009; Hansen et al., 2011; Goguen et al., 2013). The best fit parameters for the background gas distribution are c 0 = 2.7 × 1 0 10 and c 1 = 1.6 × 1 0 23 . The number density due to the background population is at least two orders of magnitude below the peak exospheric density measured by Cassini/INMS. According to our simulations, the column density due to the background population does not exceed 1%–2% at the peak of the column density during UVIS observations. Our findings indicate that Tiger Stripes contribute 23%–32% to the total plume source rate, varying from 6.4 × 1 0 27 s − 1 to 29 × 1 0 27 s − 1 , which is crucial for explaining Cassini/UVIS observations from 2007 to 2010.

Ground-based observations of the lunar exosphere performed at different phases indicate that sodium number density at the subsolar point is close to n ≈ 60 c m − 3 and varies in altitude with a scale height of 75–120 km (Morgan and Shemansky, 1991; Potter and Morgan, 1988a; Sprague et al., 1992; Sprague et al., 1998; Cremonese and Verani, 1997). At such low densities near the surface, the exosphere is surface-bound. This means that collisions can be neglected in the entire exosphere, and sodium atoms move ballistically, affected only by gravity, solar radiation pressure, and interaction with the surface. Since lunar gravity dominates a neutral particle’s motion in the vicinity of the surface, the initial energy distribution of ejected particles, to a large degree, defines the structure of the exosphere.

Sodium atoms can be released into the exosphere of the Moon through various source processes, with the most significant being thermal desorption, photo-stimulated desorption, sputtering by solar wind, and micrometeorite vaporization (e.g., Morgan and Killen, 1997; Stern, 1999; Morgan and Shemansky, 1991; Sprague, 1992; Smyth and Marconi, 1995; Sarantos et al., 2011; Cintala, 1992).

Micrometeorite vaporization, initially suggested by Potter and Morgan (1988b) as a sodium source in the Moon’s exosphere, has been confirmed by observational studies during meteor stream passages (Hunten et al., 1991; 1998). Numerical modeling by Cintala (1992); Mangano et al. (2007); Sarantos et al. (2008); Lee et al. (2011) investigated its impact on exospheric dynamics. The sodium source rate, ranging from approximately 2 × 1 0 4 to 4.9 × 1 0 4 atoms c m − 2 s − 1 (Smyth and Marconi, 1995; Morgan and Killen, 1997; Bruno et al., 2006; Bruno et al., 2007), is usually assumed to be uniformly distributed over the lunar surface and fits a Maxwellian distribution with temperatures of 3,000–5000 K (Sarantos et al., 2010; Sarantos et al., 2012b). The deviation from the uniform distribution was investigated by Cremonese et al. (2013), suggesting that the distribution peaks at the equator.

The sputtering of the lunar surface by solar wind and magnetospheric ions significantly contributes to sodium mobilization from the lunar regolith into the exosphere, enhancing sodium atom diffusion from deeper layers for exospheric injection via desorption processes. The sodium flux from sputtering is quantified as Φ s = Y f Na Φ i , where f Na = 0.0053 is the surface fraction of sodium, Y is the total sputtering yield, and Φ s and Φ i are sputtering and incident ion fluxes, respectively (Killen et al., 2007; Mura et al., 2009; Mura et al., 2010; Burger et al., 2010a). The sputtering flux can be reformulated as Φ s = σ s σ Φ i , where σ is the surface density of sodium, and σ s is a sputtering cross-section. Yakshinskiy and Madey (1999) and Killen et al. (2012) suggest that the surface sodium number density is σ = ( 3 − 3.9 ) × 1 0 12 c m − 2 . Laboratory studies by Dukes et al. (2011) show the cross-section for sputtering of Na on mineral surfaces is σ s = 1 0 − 15 c m 2 a t o m s − 1 , with most sodium sputtered as ions. With a sputtering yield Y = 0.05 and precipitating flux Φ SW = 1 − 4 × 1 0 8 c m − 2 s − 1 , the flux of sputtered sodium is about 2.65 × 1 0 4 Na c m − 2 s − 1 at the sub-solar point (Wurz et al., 2007; Sarantos et al., 2010; Sarantos et al., 2012b). The energy distribution of sputtered sodium atoms, independent of surface temperature, is given in Equation 22. f E e ∼ E e E e + E b 3 1 − E e + E b E i m H + m Na 4 m H m Na 2 (22) Here, E e is the ejected energy, E b = 2 eV is the binding energy, E i is the kinetic energy of the incident ion, and m H and m Na are the masses of solar wind protons and sodium atoms, respectively (McGrath et al., 1986; Mura et al., 2007; Sarantos et al., 2012a). The exit angle with respect to the local surface normal, θ , follows the Knudsen cosine law f ( θ ) = cos ⁡ θ (Cassidy and Johnson, 2005).

Photon stimulated desorption (PSD) has been recognized as a significant sodium source in the Moon’s exosphere, initially suggested by McGrath et al. (1986). The process involves ejecting particles influenced by solar photons, with the total flux of Φ psd = ( Φ ph ⁡ cos ⁡ ψ σ psd σ ) / r h 2 , where Φ ph = 2 × 1 0 14 photons c m − 2 s − 1 is the flux of solar photons with h ν ≥ 5 eV, σ is the surface sodium number density, σ psd is the cross-section of photon stimulated desorption, r h is the heliocentric distance, and ϕ is the solar angle (Killen et al., 2012). The typical value of the flux has been evaluated by Sarantos et al. (2012b) to be 2 × 1 0 6 atoms c m − 2 s − 1 . Yakshinskiy and Madey, (1999), Yakshinskiy and Madey, (2004) discuss the experimental investigations of sodium ejection from lunar samples. The energy distribution of ejected sodium atoms is f psd ( E ) ∼ ( E U β ) / ( ( E + U ) 2 + β ) , where U = 0.052 eV is characteristic energy related to the surface binding energy and β = 0.7 is a shape parameter (Mouawad et al., 2011; Sprague et al., 2012; Burger et al., 2010b).

Thermal desorption is a mechanism for transferring sodium atoms from the lunar regolith to the exosphere, primarily through the sublimation of adsorbed sodium atoms with the total flux of Φ td = ν σ ⁡ exp − U td / k T , where ν is the vibrational frequency of the adsorbed atom with a binding energy of U td = 0.259 eV (Stern, 1999; Mura, 2012). The vibrational frequency, typically set at ν = 1 0 13 s − 1 (Killen et al., 2004; Killen et al., 2007; Yakshinskiy et al., 2000; Milillo et al., 2011; Leblanc and Johnson, 2010). The surface temperature, T , varies from 300 ⁡ cos 1 / 4 ⁡ ψ + 100   K on the dayside to 100   K on the nightside, where ψ is the subsolar angle (Killen et al., 2012). A more detailed temperature model based on LRO DIVINER data was developed by Hurley et al. (2015).

When a sodium atom collides with the lunar surface, it may be scattered, adsorbed, or chemically bound (Sprague, 1992). Particles that are directly scattered quickly equilibrate to the local surface temperature. The proportion of particles that become adsorbed or bound is governed by the sticking coefficient, which is highly dependent on the local surface temperature (Yakshinskiy and Madey, 2005). Killen et al. (2012) note that only about 50% of these adsorbed particles can be re-emitted into the exosphere through various desorption processes.

The rate of sodium photoionization by solar UV radiation has been a subject of considerable debate. Huebner et al. (1992) provided both theoretical and empirical photoionization rates of 5.92 × 1 0 − 6 s − 1 and 6.52 × 1 0 − 6 s − 1 at 1 AU for the quiet and active Sun, respectively, and higher laboratory-derived rates of 1.62 × 1 0 − 5 s − 1 and 1.72 × 1 0 − 5 s − 1 for similar conditions. Historically, estimates of sodium photoionization lifetime in the Moon’s exosphere ranged from 12 to 17 h (Potter and Morgan, 1988a; Killen and Ip, 1999; Flynn and Mendillo, 1995; Sprague et al., 1992; Hunten, 1992; Ip, 1991). In more recent work the sodium photoionization lifetime of 36–47 h was assumed (Wilson et al., 1999; Wilson et al., 2003; Line et al., 2012), which is consistent with that by Huebner et al. (1992).

Figure 2 illustrates an example of applying AMPS to model the distribution of sodium in the Moon’s exosphere. Our simulations indicate that photon-stimulated desorption is the dominant source of sodium in the lunar exosphere, with its rate exceeding that of meteoritic impact vaporization by a factor of approximately 8–9. The total source rate of sodium is estimated to be 1.68 × 1 0 23 atoms per second. Surface interactions play a crucial role in the behavior of sodium atoms, with the majority being reabsorbed by the lunar surface rather than escaping into space. Specifically, about 70% of sodium atoms produced by meteoritic impact vaporization and 25% of those produced by photon-stimulated desorption are reabsorbed, leading to an estimated escape rate from the exosphere of 5.3 × 1 0 22 atoms per second.

Both kinetic and fluid-type methods have been successfully used to model cometary comae. This section illustrates the applications of these methods in studying ESA’s Rosetta mission target, comet 67P/Churyumov-Gerasimenko.

Unlike planets’ dense, collisional atmospheres, cometary comae are characterized by minimal interactions between particles. Cometary comae are a distinctive phenomenon within the solar system, functioning as a planetary atmosphere with minimal influence from gravity. As a comet nears the Sun, water vapor and other gases sublimate, forming a cloud of gas, ice, and refractory materials (such as rocky and organic dust) expelled from the nucleus’s surface. The sublimated gas molecules may experience frequent collisions and participate in photochemical reactions near the nucleus. Due to the comet’s negligible gravity, it generates a large and highly variable dusty coma extending far beyond the cometary nucleus’s size (Hässig et al., 2015).

The sublimation of volatiles from a cometary nucleus is the primary source of volatiles in a cometary coma. The nucleus is covered by a porous layer of ice and solid grains, which experiences periodic solar illumination due to the rotation of the nucleus. This leads to the sublimation of volatiles, thereby contributing to the formation of the coma. Thermal re-radiation, solid-state heat conduction, and mass and energy subsurface transport also play significant roles. These processes collectively form the foundation for the thermophysical modeling of the nucleus’s gas production, which is used to define the boundary conditions on the nucleus’s surface for subsequent modeling of gas dynamics in cometary comae (Davidsson, 2024; Davidsson and Gutiérrez, 2005; Davidsson et al., 2010a).

With a typical density near the surface of n ≈ 1 0 19   m − 3 and a water collisional cross-section of σ ≈ 1 0 − 19   m − 2 , the mean free path in the coma is λ = 1 2 n σ < 1   m , which allows for the application of hydrodynamic methods near the nucleus. As the distance from the nucleus increases, the collision rate in the outflowing gas rapidly decreases, making momentum exchange within the gas phase negligible beyond approximately 1 0 3 km from the nucleus. Therefore, accurately modeling a cometary coma from the near nucleus region to the far coma requires a kinetic approach, where gas thermalization is described at the level of individual particle collisions. The cross-section values for collisions between major neutral components used in such modeling are summarized in Table 1 based on the works of Crifo (1989) and Combi (1996).

As the primary species, water dominates the thermodynamic balance of cometary comae through its photodissociation and radiation cooling. Its rotational transitions may allow radiation cooling or heating, which could be essential for controlling H 2 O velocity and temperature in the intermediate and outer coma of active comets (Bockelée-Morvan and Crovisier, 1987; Xie and Mumma, 1996; Marconi and Mendis, 1982). For radiation cooling through the emission of rotational lines to become efficient, the rotational degrees of freedom must be coupled with translational ones. This coupling is possible only when the coma is in the collision-dominated regime. An empirical formula for energy loss by radiation was proposed by Shimizu (1976); Crovisier (1984); Combi (1996).

The evolution of daughter species in the coma primarily depends on the absorption of solar radiation and interaction with the solar wind, which consists mainly of protons, He + ions, and electrons of solar origin. Absorption of solar radiation leads to the excitation of an atom or a molecule, followed by photoionization or photodissociation. Charge exchange and impact ionization due to interaction with the solar wind are the primary channels for the decay of daughter species. The relative density of ions is typically negligible in the collision zone of the coma at moderate to large heliocentric distances. As a result, the effect of ion interactions with neutral species can be neglected when modeling that region.

Starting a few hundred kilometers from the nucleus, the gas dynamics in a cometary coma are primarily influenced by the formation of energetic daughter species (Combi and Smyth, 1988; Xie and Mumma, 1996; Festou, 1999; Gunnarsson et al., 2002). The dominant photolytic process in a coma is the photodissociation reaction H 2 O + h ν → OH + H + Δ E , which occurs at a rate of β = 1.2 × 1 0 − 5   s − 1 at a heliocentric distance of 1 AU. This reaction produces a mean energy excess of Δ E = 1.78   eV , corresponding to a mean ejection velocity of 18.5   km   s − 1 for H atoms in the rest frame of the parent molecule. Other photodissociation reactions critical in modeling a cometary environment are summarized in Table 2. Photolytic heating, caused by momentum exchange between highly energetic daughter species and other components of the coma, is efficient only in the near-nucleus region, where the dissociation products are thermalized through collisions.

The Direct Simulation Monte Carlo (DSMC) method was employed to model the coma of comet 67P/Churyumov-Gerasimenko in both full 3D and axially symmetric 2D, capturing the distribution and dynamics of major volatile species (Combi, 1996; Tenishev et al., 2008). This approach is crucial for understanding the complex interactions within the coma, especially given the non-equilibrium conditions present due to low particle densities and varying illumination conditions. The DSMC model was used to analyze data from the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) and the Visible and Infrared Thermal Imaging Spectrometer (VIRTIS) onboard the ESA Rosetta mission.

In that modeling, the nucleus was represented with a high-resolution shape model, and the activity was distributed over the surface based on local illumination and empirical data derived from observations (Tenishev et al., 2016). The major species needed to be considered are H 2 O, C O 2 , CO, and O 2 .

The surface activity can be described using a 25-term (order 4) spherical harmonic expansion, capturing the complex activity patterns observed in different regions of the comet. The coefficients of this expansion were determined through a least-squares optimization method using Rosetta/ROSINA data (Fougere et al., 2016a; Fougere et al., 2016b). The gas flux at the nucleus’ surface was defined by the local surface temperature derived from the thermophysical model of the comet’s nucleus (Davidsson and Gutiérrez, 2004; Davidsson and Gutiérrez, 2006; Davidsson et al., 2010b). The gas flux at the nucleus surface is in Equation 23 F k = G Θ SZA f k / R AU β (23) Here, G ( Θ SZA ) represents the variation of flux with solar zenith angle (SZA), f k is the local surface activity factor, and R AU is the heliocentric distance in astronomical units, and β is the exponent corresponding to a power law governing the comet’s source rate evolution with heliocentric distance (Fougere et al., 2016a; Fougere et al., 2016b).

AMPS successfully replicated the temporal and spatial variations observed in the coma of comet 67P/Churyumov-Gerasiminko. It accurately captured the strong seasonal variations in outgassing patterns driven by the comet’s axial tilt. The model demonstrated a strong correlation between water vapor and molecular oxygen throughout the observation period, while C O 2 and CO showed varying correlations pre- and post-equinox. The simulated column densities for H 2 O and C O 2 matched well with observational data, particularly post-equinox. Some discrepancies were noted, especially for C O 2 densities early in the mission, which could be attributed to the limitations of a single power-law approximation for the extended period.

Figure 3 compares the model with Rosetta’s VIRTIS data, further supporting the model’s accuracy in capturing both the large-scale coma structure and the finer details of local outgassing features.

Mars’ hot oxygen corona, first predicted by McElroy (1972) as a result of dissociative recombination of ionospheric O 2 + with electrons in the exosphere, was confirmed by Rosetta/ALICE during its 2007 Mars flyby. The spacecraft detected the atomic oxygen resonance line at 130.4 nm, providing direct evidence of the corona’s existence (Feldman et al., 2011; Gröller et al., 2014).

Later, the MAVEN mission provided detailed observations of Mars’ upper atmosphere (Deighan et al., 2015). While direct escape rate measurements are not feasible (Leblanc et al., 2017), MAVEN’s instruments have allowed indirect estimates of hot oxygen density and escape rates, ranging from 1.2 to 5.5 × 1 0 25   s − 1 , depending on season and solar activity (Jakosky et al., 2015; Lillis et al., 2017). Additionally, MAVEN’s Imaging Ultraviolet Spectrograph (IUVS) has systematically observed 130.4 nm coronal emissions since 2014, revealing exobase densities from 6.3 × 1 0 3 to 8.0 × 1 0 3   c m − 3 during low solar activity and up to 1.2 × 1 0 4   c m − 3 during moderate activity, with effective temperatures between 4,100 and 4500 K (Chaufray et al., 2015; Deighan et al., 2015; McClintock et al., 2015; Qin et al., 2024).

Mars experiences significant atmospheric loss, primarily of hydrogen and oxygen, with total escape rates ranging from 2 to 3 kg/s, primarily from water and carbon dioxide. Oxygen loss occurs through several processes: solar wind-driven electric field acceleration, which removes oxygen ions at a rate of 5 × 1 0 24 atoms per second (130 g/s); photochemical processes, particularly dissociative recombination, which contribute a neutral oxygen loss rate of 5 × 1 0 25 atoms per second (1,300 g/s); and sputtering caused by solar wind ions, leading to an additional loss of 3 × 1 0 24 oxygen atoms per second (80 g/s) (Jakosky et al., 2018).

Space weather significantly affects hot oxygen production in Mars’ upper atmosphere. For example, the X8.2-class solar flare on 10 September 2017, increased hot oxygen production by up to 45% and photochemical escape rates by 20% due to enhanced ultraviolet flux (Fox, 2004; Fox and Hac, 2009; Lee et al., 2018; Cravens et al., 2017).

Both empirical and physics-based modeling are used to analyze MAVEN’s data. Such, Ramstad et al. (2023) developed an empirical method using MAVEN’s in situ measurements to infer Mars’ hot oxygen density near the exobase, finding densities of 6 – 7 × 1 0 3   cm − 3 and effective temperatures of 3700 – 3800   K . Various physics-based models were developed to investigate the production mechanisms and distribution of those energetic atoms (e.g., Chaufray et al., 2007; Fox and Hac, 2009; Gröller et al., 2014; Hodges, 2000; Lee et al., 2015a; b; Lee, 2014; Valeille, 2009; Valeille et al., 2010).

The most commonly used numerical method for analyzing MAVEN/IUVS observations of Mars’ hot oxygen corona is Monte Carlo modeling, which simulates the transport of hot oxygen from its production in the thermosphere through the exosphere and into the corona to compare with observations (Leblanc et al., 2017; Lee et al., 2015a).

As hot oxygen (O) atoms propagate through the thermosphere, they may collide with ambient atoms and molecules, losing energy and becoming thermalized before reaching the corona. However, collisions between a hot O atom and a thermal O atom in the thermosphere can energize the latter enough to enter the corona. The energy transferred during collisions highly depends on the total and angular differential scattering cross-sections. Forward scattering is crucial in modeling this process because it affects how hot O atoms retain energy after collisions with C O 2 , reducing thermalization and allowing more hot O atoms to escape Mars’ gravity. Accurate forward scattering modeling is essential for estimating the density of the hot O corona and photochemical loss rates, which are critical for understanding Mars’ atmospheric evolution (Gacesa et al., 2020; Fox and Hać, 2014; Gacesa et al., 2017; Lee et al., 2020; Kharchenko et al., 2000; Valeille, 2009; Valeille et al., 2009b).

The results of modeling the hot O population in Mars’ corona using AMPS are illustrated in Figure 8. The simulation was initialized using outputs from the Mars Global Ionosphere Thermosphere Model (M-GITM) (Bougher et al., 2014), which provides a detailed description of the background thermosphere, including temperature, wind, and density profiles of significant species such as O, C O 2 , and N 2 . M-GITM accounts for diurnal, seasonal, and solar cycle variations, ensuring that the initial conditions accurately reflect the current Martian atmospheric state.

Below the exobase, the atmosphere is assumed to be in collisional equilibrium with a Maxwellian velocity distribution. This approach was used to model the transport of hot O atoms through the thermosphere, where they originated and interacted with the background population of major thermospheric species. The simulation incorporated key model parameters, including the distribution of hot O source strength and the density of background species, derived from M-GITM outputs for specific solar cycle and seasonal conditions, ensuring that the thermospheric state was accurately represented.

Energy-dependent forward scattering cross-sections for O-O and O- C O 2 interactions from Kharchenko et al. (2000) were used in the modeling. The escape energy for oxygen atoms was set at 1.97 eV, corresponding to the gravitational potential energy required to overcome Mars’ gravity. To ensure statistical convergence and accurate representation of the hot oxygen population, the simulation utilized a large number of particles, typically on the order of 1 0 6 to 1 0 7 . The effect of producing secondary hot O is discussed in previous studies (e.g., Lee et al., 2015b; Lee et al., 2018; Lee et al., 2020).

The parameters for the dissociative recombination reaction used in this study are summarized in Table 3. The reaction rate constant, as detailed in Equation 27, is adapted from Mehr and Biondi (1969) and has been employed in our previous research on Mars’ hot oxygen corona (e.g., Lee et al., 2015b; Valeille et al., 2010). α = 1.95 × 1 0 − 7 300 T e 0.7 c m 3 s − 1 , 300 < T e < 1200 K 7.39 × 1 0 − 8 1200 T e 0.56 c m 3 s − 1 , 1200 < T e < 5000 K (27)

The Martian hot oxygen corona simulation using the AMPS model, cuand ITM, rvealed that the hot oxygen population is highly variable and strongly influenced by both solar cycle and seasonal changes. The study found that solar maximum periods produce significantly higher densities and escape rates of hot oxygen, with global escape rates ranging from 1.14 × 1 0 25 s − 1 during solar minimum to 5.18 × 1 0 25 s − 1 during solar maximum (Lee, 2014). Seasonal variations were also prominent, with increased densities observed during perihelion due to the associated changes in thermospheric conditions (Lee et al., 2015b).

Venus’s thermosphere, ionosphere, and exosphere have been extensively studied over several decades through observations from a range of spacecraft missions. These investigations began during the Soviet Venera and US Mariner missions, continued through the nearly 14-year Pioneer Venus mission, and have extended into more recent times with the Venus Express mission (Bougher et al., 1997; Schubert et al., 2007; Gérard et al., 2017).

Most of our current knowledge of Venus’ upper atmosphere and corona comes from measurements made by in-situ and remote sensing experiments on the Pioneer Venus Orbiter (PVO) from December 1978 to October 1992. The “hot” O corona was observed with an ultraviolet spectrometer (UVS) onboard PVO by measuring the OI resonance triplet near 1304 Å (Nagy et al., 1981).

Since 2006, systematic monitoring by Venus Express (VEx) instruments has enhanced our understanding of Venus’s atmosphere. VIRTIS observations have measured 3D temperatures and derived thermal wind fields at 40–90 km (Piccialli et al., 2012; Piccialli et al., 2008) and mapped the highly variable O 2 IR nightglow distribution (Drossart et al., 2007; Gérard et al., 2008; Gérard et al., 2009; Gérard et al., 2017; Hueso et al., 2008; Piccioni et al., 2009; Soret et al., 2012; Soret et al., 2014) to trace winds at 90–130 km. SPICAV’s airglow observations of NO emissions (Collet et al., 2010; Gérard et al., 2008; Gérard et al., 2009; Stiepen et al., 2013) and vertical profiles of atmospheric density (Bertaux et al., 2007) have confirmed patterns observed with previous missions (Stewart et al., 1980). These VEx datasets provide critical insights into Venusian atmospheric processes (Brecht et al., 2011).

The BepiColombo spacecraft’s second fly-by of Venus on 10 August 2021, using the Mass Spectrum Analyzer (MSA) on Mio, BepiColombo’s magnetospheric orbiter, revealed cold oxygen ( O + ) and carbon ( C + ) ions at six planetary radii with a flux of 4 ± 1 × 1 0 4   cm − 2   s − 1 . The C + to O + ion ratio, at most 0.31 ± 0.2 , suggests contributions from CO or water group ions as oxygen sources. These findings, observed near the magnetic pileup boundary, have significant implications for understanding Venus’s atmospheric evolution and the planet’s water history, offering insights into the climate and habitability of terrestrial planets and exoplanetary systems (Hadid et al., 2024).

In the upper atmosphere of Venus, suprathermal O atoms are primarily produced by exothermic reactions, such as the electron dissociative recombination of O 2 + ions, similar to the processes observed on Mars (Krestyanikova and Shematovich, 2006; Lichtenegger et al., 2006; Gröller et al., 2010; Gröller et al., 2012; Kella et al., 1997). Unlike the lower thermosphere, where these hot atoms are rapidly thermalized by frequent collisions, in the upper thermosphere and exosphere, where collisions are rare, these hot atoms remain unthermalized, forming extended coronae of H and O atoms above the exobase, located near 200 km. Early estimates of “hot” O density derived from PVO/UVS data suggested a significant population, but later observations did not confirm these findings. Instruments such as SPICAV and ASPERA-4 on board Venus Express, which were sensitive enough to detect “hot” O if present at densities similar to those suggested by PVO data, have not detected this population above their detection thresholds (Bertaux et al., 2007; Galli et al., 2008; Lichtenegger et al., 2009). Despite this, earlier studies using data from the PVO/UVS instrument and the Venera 11 mission had established the presence of a corona of hot O atoms, with densities around 1 0 4   cm − 3 at altitudes of 1000 km (Nagy et al., 1981; Nagy and Cravens, 1988; Bertaux et al., 1981).

The Mars oxygen corona model detailed in Section 4.4 was modified to study Venus’s exosphere and corona (Tenishev et al., 2022). This adaptation takes advantage of the fact that the production and dynamics of hot oxygen on Mars and Venus are similar (Valeille et al., 2009b; Valeille et al., 2009a; Lee et al., 2015a). Similarly, we use the output of the Venus Thermosphere General Circulation Model (VTGCM) of Venus’ thermosphere/ionosphere composition to determine the source of hot O and characterize its interaction with major thermospheric species (Bougher et al., 1988).

Our findings indicate that the altitude distribution of hot oxygen during solar maximum aligns closely with observations from the Pioneer Venus Orbiter. Conversely, during solar minimum, we observe a significant decrease in the oxygen density of the corona, consistent with Venus Express’s non-detection of the oxygen corona. The conditions during moderate solar activity naturally lie between these extremes. Our results indicate variability in the density of the extended oxygen corona around Venus by a factor of six over a solar cycle, aligning with observations suggesting a significant reduction in density during low solar activity periods (Gérard et al., 2017). The lack of corona detection by ASPERA-4 and SPICAV onboard Venus Express during solar minimum further supports our findings, highlighting the effect of solar conditions on the visibility of Venus’ oxygen corona (Lichtenegger et al., 2009). The modeling results of Venus’s exosphere and corona are in Figure 9.