Hydrostatic escape

Where the gas is too rarefied for a wind, escape proceeds particle by particle from the exobase, the level where the mean free path first reaches the local scale height. Each species escapes at the smaller of two bottlenecks: how fast it effuses from the exobase, and how fast diffusion resupplies it from below. The branch is natively per-species. Symbols shared across the section are in the notation tables.

The upper structure

The atmosphere profile usually stops well below the exobase, so the branch builds an extended upper structure above the profile top, with composition and mean molecular mass \(\mu_0\) frozen at their values there. The temperature follows the Bates profile in the form Yelle (2024) uses, their Eq. (19) ¹, and the radius follows from hydrostatic balance, both in the log-pressure coordinate \(\zeta = \ln(P_\mathrm{top}/P)\):

with \(r(\zeta = 0) = R_0\), the shape parameter \(\gamma_\mathrm{B}\) (gamma_bates, default 0.75), and the exospheric temperature \(T_\infty\) that the profile approaches (T_exo_value, default 1000 K). Because the Bates profile integrates analytically, so does the radius,

so the structure carries no step error. It ends where the local Jeans parameter reaches 2, since the structure is unbound beyond that point, and the flag extension_unbound is raised. The exobase is the lowest point of this structure where

with \(\sigma_\mathrm{C}\) the mixture cross section of the collisionality switch, taken over the species at the profile top at the local temperature. This defines \(R_\mathrm{exo}\), \(n_\mathrm{exo}\), and \(T_\mathrm{exo} = T(R_\mathrm{exo})\).

The extension is a requirement, not a refinement. The Jeans parameter at the true exobase can differ from its photospheric value by an order of magnitude, and evaluating the escape on photospheric values biases rates toward false retention by up to three orders of magnitude (Johnson et al. 2013 ²).

The exospheric temperature. \(T_\infty\) is prescribed, and it is the branch's dominant sensitivity, because the rate depends on it exponentially. A value below \(T_\mathrm{top}\) would build a thermosphere that cools with height, whose exobase is more strongly bound than its anchor, so it is raised to \(T_\mathrm{top}\) and flagged. An optional estimator (T_exo_mode = 'thermostat') solves the heating and cooling balance of the wind temperature at the profile top instead. A local balance without conduction is biased high by construction, so it is not the default.

Jeans effusion

Each species \(i\), of particle mass \(m_i\), leaves the exobase with the Jeans effusion velocity of Yelle (2024), their Eqs. (20) and (21) ¹,

where \(\lambda_i\) is the species Jeans parameter. Direct simulation Monte Carlo runs find the escape rate above the equilibrium Jeans rate, because collisions above the exobase and a non-Maxwellian distribution there feed the escaping tail: about 1.7 times at \(\lambda = 6\), falling to about 1.4 at \(\lambda = 15\) (Volkov et al. 2011 ³⁴). The code interpolates that factor linearly,

holding the end values outside the simulated range and flagging it (volkov_extrapolated). The effusion flux, per unit area of the sphere of radius \(R_0\) (Yelle 2024, their Eq. (15)), is

with \(\tilde{X}_i(\zeta_\mathrm{exo})\) the mole fraction at the exobase from Eq. (6).

Diffusive supply

Below the exobase each species diffuses through the frozen background with molecular diffusion coefficient \(D_i = b_i / n\) and eddy diffusion coefficient \(K_{zz}\) (the profile's top value, or the setting kzz, default 300 m² s⁻¹). The mixture diffusion parameter follows Blanc's law, \(b_i = (1 - X_i) / \sum_{k \ne i} (X_k / b_{ik})\), over binary parameters \(b_{ik}\) that each carry a provenance class. Thermal diffusion enters through an effective mass, Yelle (2024) Eq. (4), with the thermal diffusion factor \(\alpha_\mathrm{T} = -0.25\) that Yelle adopts for light species. The mole fraction in diffusive equilibrium, their Eq. (9), is

with \(X_i(\zeta = 0)\) the species mole fraction at the profile top. The limiting flux of their Eqs. (10) and (11), the largest flux diffusion can deliver to the exobase, is

The exobase is located by root finding on the closed-form structure, so it does not move with the quadrature grid of these integrals. Both integrals are first order in the log-pressure step. The quadrature grid is therefore doubled from hydrostatic_levels_min until two successive refinements each change the bulk rate, and every element whose rate clears the rate floor, by less than hydrostatic_rtol (default 1%). Two, so that one chance agreement between two grids does not end the refinement, and every element, since a trace whose supply integral is still moving can sit under a bulk rate that the dominant species sets alone. Refinement also stops at hydrostatic_levels_max, which must be at least four times hydrostatic_levels_min for two changes to fit. The levels used and the last change are reported in the diagnostics, and the finest grid's values are returned, never an extrapolation.

The combined rate

The two bottlenecks combine by the harmonic mean of Yelle (2024), their Eq. (14) ¹, which lies below both and tends to the smaller:

The most abundant species at the profile top supplies itself: it takes \(\tilde{X}_i = X_i(\zeta = 0)\) and the Jeans flux alone, \(\Phi_i = \tilde{\Phi}_{\mathrm{J},i}\). So does any species whose Jeans rate is already below the rate floor, to save the cost of its integrals. Its rate is then an estimate, not a bound, since diffusive enrichment can raise a light species' full rate above it (2 to 5 times on the trace hydrogen and helium profiles tested), enough to carry a rate just under the floor above it. Such species are marked pruned in the per-species detail. Every flux on this page is per unit area of one reference sphere, the profile top at radius \(R_0\), following Yelle (2024). The two bottlenecks live at different heights (effusion at the exobase, diffusion between the profile top and the exobase), and Eq. (8) can only combine them once both are per unit area of the same sphere, which is why Eq. (5) rescales the exobase flux by \((R_\mathrm{exo}/R_0)^2\). The choice of sphere does not change the result: in a steady outflow the same number of particles crosses every sphere, so the flux times the sphere's area, which is the rate, is the same at \(R_0\) as at \(R_\mathrm{exo}\). The hydrostatic rate is

with each species rate mapped onto elements by the mass fractions of its formula. Heavy-element rates are lower limits, since the nonthermal channels that dominate heavy-species loss from real exospheres are absent, and every hydrostatic result carries the flag hydrostatic_lower_limit.

The exosphere stability test

Two escape temperatures decide whether the exobase can stay hydrostatic at all:

The neutral escape temperature \(T_\mathrm{esc}\) is where thermal energy rivals binding energy. The plasma escape temperature \(T_\mathrm{esc,+}\) is half of it, because in an ionised exosphere the ambipolar electric field shares each ion's binding energy with its electron (Chatterjee & Pierrehumbert 2026, their Eq. (34) ⁵). An exobase with \(T_\mathrm{exo}\) above half the selected escape temperature is unstable by their Figure 10 criterion.

In such states the hydrodynamic candidate competes with the Jeans rate, the larger one is dispatched with its own label, and the state is flagged gate_rerouted. Taking the larger keeps the dispatched rate from falling as the exosphere heats past the threshold. Neither candidate describes the outflow of an unstable exosphere, as the limitations page explains.

The setting gate chooses the neutral (default) or the plasma escape temperature. Both are always computed, with the ionisation fraction at the exobase under the surface-averaged flux \(F_\mathrm{XUV}/4\). States where the two choices disagree are flagged contested_ion, with both branch rates recorded, because the ion physics that would decide them is not modelled.