A CRITERION FOR THE VALIDITY OF PARKER'S MODEL IN THERMAL ESCAPE PROBLEMS FOR PLANETARY ATMOSPHERES

Hydrodynamic modeling based on equations of continuum fluid flows is one of the major approaches for simulation of thermally escaping atmospheric flows. Historically, such models were first applied to study a similar phenomenon of solar wind flow. In particular, Parker (1964a, 1964b) developed a one-dimensional hydrodynamic model of thermal escape from solar corona, assuming that the temperature drops to zero far from the star. Later on, Parker's model became a basis for numerous studies of thermal escape from atmospheres of planetary bodies in the outer solar system (e.g., Watson et al. 1981; Chassefière 1996; Krasnopolsky 1999; Strobel 2008; Zhu et al. 2014) and exoplanets (e.g., Tian et al. 2005; Muñoz 2007).

In order to satisfy the condition of zero temperature far from the surface, it is necessary to continue solution of Parker's model above the exobase into exosphere, where the Fourier law for the heat flux fails (e.g., Volkov et al. 2011a). In spite of this basic flaw, Parker's model has never been systematically tested against kinetic models, which can accurately predict the heat flux and structure of non-equilibriums flow in the exosphere. Thus, the usefulness of Parker's model for thermal escape problems remains questionable.

One-dimensional steady-state hydrodynamic equations of inviscid thermally conducting gas flows have a few distinct types of solutions depending on conditions far from the surface (Freeman & Johnson 1972; Roberts & Soward 1972), including the Chamberlain-type solar breeze solutions with decreasing velocity in the far field (Chamberlain 1961, 1965) and Parker-type solutions with monotonously increasing gas velocity (Parker 1964a, 1964b). The kinetic simulations predict the solar-breeze-type thermally escaping solutions only in conditions of the atmospheric blow-off or transition from the blow-off to Jeans-like escape (Volkov 2014). In the present letter, thermal escape in the Jean-like regime is considered, when the kinetic simulations (Volkov et al. 2013) predict monotonously increasing gas velocity, so that only Parker-type hydrodynamic solutions are chosen for comparison.

The present letter reports for the first time results of systematic comparison of thermal escape rates obtained in hydrodynamic and kinetic simulations. Based on this comparison, the range of surface conditions is found, where Parker's model can be used for calculations of the escape rates. It is shown that Parker's model satisfactorily predicts the escape rate only in the so-called high-density (HD) asymptotic regime. The criterion, determining surface conditions for which Parker's model is capable of quantitative prediction of the escape rate is then established theoretically. Although both hydrodynamic and kinetic simulations are performed for atmospheres without heating, the obtained criterion can be also used to assess the applicability of Parker's model in the externally heated atmospheres based on the parameters on the top of heated atmospheric layer as it is discussed in Section 6.

In one-dimensional steady-state atmospheric flows without heating, the distributions of gas density ρ, velocity u, and temperature T can be described by the Navier–Stokes equations

$$\begin{eqnarray}&&{r}^{2}\rho u={{\rm{\Phi }}}_{{\rm{m}}},\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&\rho {uu}^{\prime} =-(\rho {\mathfrak{R}}T)^{\prime} +{\tau ^{\prime} }_{{rr}}+\displaystyle \frac{3{\tau }_{{rr}}}{r}-\displaystyle \frac{{GM}}{{r}^{2}}\rho ,\end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{m}}}\left({c}_{p}T+\displaystyle \frac{{u}^{2}}{2}-\displaystyle \frac{{GM}}{r}\right)+{r}^{2}\left({q}_{r}-u{\tau }_{{rr}}\right)={{\rm{\Phi }}}_{{\rm{e}}},\end{eqnarray} \tag{ 1c }$$

where r is the radial distance, the prime symbol denotes derivatives with respect to r, ${\mathfrak{R}}={k}_{{\rm{B}}}/m,$ ${c}_{p}=(5+\zeta ){\mathfrak{R}}/2,$ k_B is the Boltzmann constant, m and ζ are the mass and number of internal degrees of freedom of a molecule, M is planet's mass, G is the gravitational constant, Φ_m and Φ_e are the mass and energy escape rates per steradian, and τ_rr and q_r are the viscous stress and heat flux. Atmospheric flows are usually sought with constrains on the surface density ρ₀ and temperature T₀:

$$\begin{eqnarray}&&\mathrm{at}\quad \quad r={r}_{0}:\quad \rho \left({r}_{0}\right)={\rho }_{0},\quad T\left({r}_{0}\right)={T}_{0},\end{eqnarray} \tag{ 2 }$$

where r₀ is the surface radius, and the density is assumed to drop to zero far from the planet:

$$\begin{eqnarray}&&\mathrm{at}\quad \quad r\to \infty :\quad \rho \left(r\right)\to 0.\end{eqnarray} \tag{ 3 }$$

Two more boundary conditions are necessary to solve Equations (1), since Φ_m and Φ_e are unknown.

Parker (1964a, 1964b) made three basic assumptions, which helped him to resolve the issue of the missing boundary conditions. First, he considered outflows with a monotonically increasing velocity, which is much smaller than the isothermal sound velocity $\sqrt{{\mathfrak{R}}{T}_{0}}$ at the surface. Second, he assumed that the molecular viscosity plays a marginal role, in which case Equations (1) can be reduced by neglecting τ_rr, so that

$$\begin{eqnarray}&&u^{\prime} =u\;\displaystyle \frac{\frac{2{\mathfrak{R}}T}{r}-{\mathfrak{R}}T^{\prime} -\frac{{GM}}{{r}^{2}}}{{u}^{2}-{\mathfrak{R}}T},\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&T^{\prime} =\displaystyle \frac{{{\rm{\Phi }}}_{{\rm{m}}}\left({c}_{p}T+\displaystyle \frac{{u}^{2}}{2}-\displaystyle \frac{{GM}}{r}\right)-{{\rm{\Phi }}}_{{\rm{e}}}}{{r}^{2}{\kappa }_{00}{\left(T/{T}_{00}\right)}^{\omega }},\end{eqnarray} \tag{ 4b }$$

where the Fourier law, ${q}_{r}=-\kappa {dT}/{dr},$ with thermal conductivity, $\kappa ={\kappa }_{00}{\left(T/{T}_{00}\right)}^{\omega }$ is assumed for the heat flux, and ω is the viscosity exponent. Equations (4) require four boundary conditions, so that Parker (1964a, 1964b) finally assumed that the temperature should drop to zero far from the surface:

$$\begin{eqnarray}&&\mathrm{at}\quad \quad r\to \infty :\quad T\left(r\right)\to 0.\end{eqnarray} \tag{ 5 }$$

In this model, the isothermal Mach number ${Ma}=u/\sqrt{{\mathfrak{R}}T}$ evolves from ${{Ma}}_{0}={u}_{0}/\sqrt{{\mathfrak{R}}{T}_{0}}\ll 1$ at r = r₀ to infinity at $r\to \infty ,$ i.e., the flow passes through the isothermal sonic point, where $u=\sqrt{{\mathfrak{R}}T}$ and the denominator of Equation (4a) is equal to zero. Distributions of gas parameters are continuous only if the numerator becomes zero simultaneously with the denominator, i.e., the isothermal sonic point is the critical point of Equations (4). Using Equation (4b) in order to exclude T' from the consideration at the critical point, one can obtain

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{e}}}={{\rm{\Phi }}}_{{\rm{m}}}\left({c}_{p}{T}_{{\rm{c}}}+\displaystyle \frac{{u}_{{\rm{c}}}^{2}}{2}-\displaystyle \frac{{GM}}{{r}_{{\rm{c}}}}\right)+\displaystyle \frac{{GM}{\kappa }_{{\rm{c}}}}{{\mathfrak{R}}}-2{T}_{{\rm{c}}}{r}_{{\rm{c}}}{\kappa }_{{\rm{c}}},\end{eqnarray} \tag{ 6 }$$

where subscript "c" denotes parameters in the critical point. This equation shows that gas parameters at the critical point and escape rates are not independent of each other. For a solution passing through the critical point with monotonically increasing velocity, the density drops as 1/r² far from the surface and Equation (3) is satisfied automatically, so it should be excluded from consideration and replaced by Equation (6).

Thus, Parker's model implies a solution of Equation (4) with boundary conditions given by Equations (2) and (5) starting at the surface with ${{Ma}}_{0}\lt 1$ and passing through a critical point, where atmospheric parameters satisfy Equation (6).

In the present paper, Parker's model is analyzed in a dimensionless form in order to reduce the number of parameters governing the solutions. Equation (4) is re-written in terms of dimensionless variables $\bar{r}=r/{r}_{{\rm{c}}},$ $\bar{T}=T/{T}_{{\rm{c}}},$ and $\bar{u}=u/\sqrt{{\mathfrak{R}}{T}_{{\rm{c}}}}$ as

$$\begin{eqnarray}&&\bar{u}^{\prime} =\bar{u}\displaystyle \frac{\frac{2\bar{T}}{\bar{r}}-\bar{T}^{\prime} -\frac{{\lambda }_{{\rm{c}}}}{{\bar{r}}^{2}}}{{\bar{u}}^{2}-\bar{T}},\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&\bar{T}^{\prime} =\displaystyle \frac{{\varphi }_{{\rm{m}}}}{{\bar{r}}^{2}{\bar{T}}^{\omega }}\left(\displaystyle \frac{\gamma }{\gamma -1}\bar{T}+\displaystyle \frac{{\bar{u}}^{2}}{2}-\displaystyle \frac{{\lambda }_{{\rm{c}}}}{\bar{r}}-{{\rm{\Psi }}}_{{\rm{e}}}\right),\end{eqnarray} \tag{ 7b }$$

where the prime symbol denotes derivatives with respect to $\bar{r},$ $\gamma =(\zeta +5)/(\zeta +3),$ ${\lambda }_{{\rm{c}}}={GM}/({\mathfrak{R}}{T}_{{\rm{c}}}{r}_{{\rm{c}}})$ is the Jeans parameter at the critical point, ${\varphi }_{{\rm{m}}}={\mathfrak{R}}{{\rm{\Phi }}}_{{\rm{m}}}/({\kappa }_{{\rm{c}}}{r}_{{\rm{c}}}),$ and ${{\rm{\Psi }}}_{{\rm{e}}}={{\rm{\Phi }}}_{{\rm{e}}}/({\mathfrak{R}}{T}_{{\rm{c}}}{{\rm{\Phi }}}_{{\rm{m}}}).$ Equation (6) in reduced units transforms into

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{e}}}=\displaystyle \frac{\gamma }{\gamma -1}+\displaystyle \frac{1}{2}+\displaystyle \frac{{\lambda }_{{\rm{c}}}-2}{{\varphi }_{{\rm{m}}}}-{\lambda }_{{\rm{c}}}.\end{eqnarray} \tag{ 8 }$$

Solution of Parker's model in reduced units depends only on two surface parameters, e.g., Jeans parameter ${\lambda }_{0}={GM}/({\mathfrak{R}}{T}_{0}{r}_{0})={\lambda }_{{\rm{c}}}/({\bar{r}}_{0}{\bar{T}}_{0})$ and Knudsen number ${{Kn}}_{0}={l}_{{mfp},0}/{L}_{0},$ characterizing the strength of the gravitational field and degree of atmosphere rarefaction at the surface. Assuming that the mean free path of gas molecules ${l}_{{mfp},0}$ is given by the first approximation of the Chapmen-Enskog theory for the Variable Hard Sphere (VHS) molecular model (Bird 1994), and the atmospheric length scale L₀ is equal to r₀, one can find that ${{Kn}}_{0}={K\kappa }({T}_{0})/\left(\sqrt{{\mathfrak{R}}{T}_{0}}{\rho }_{0}{r}_{0}{\mathfrak{R}}\right)$ = ${{Kn}}_{{\rm{c}}}{\bar{T}}_{0}^{\omega -1/2}{\bar{r}}_{0}{\bar{u}}_{0}$, where ${{Kn}}_{{\rm{c}}}$ is the Knudsen number at the critical point

$$\begin{eqnarray}&&{{Kn}}_{{\rm{c}}}=\displaystyle \frac{K}{{\varphi }_{{\rm{m}}}},\quad \quad K=\displaystyle \frac{2(5-2\omega )(7-2\omega )}{15\sqrt{2\pi }}\displaystyle \frac{\gamma -1}{\gamma }\mathrm{Pr}\end{eqnarray} \tag{ 9 }$$

and Pr is the Prandtl number.

A solution of Equation (7) for the surface conditions given by λ₀ and Kn₀ can be obtained using a two-step approach. At the first downstream step, Equation (7) is solved in the downstream direction from the critical point, $\bar{r}\gt 1,$ for a fixed λ_c with the initial conditions

$$\begin{eqnarray}&&\bar{u}\left(1\right)=\;\;1,\quad \quad \bar{T}\left(1\right)=\;\;1\end{eqnarray} \tag{ 10 }$$

and Ψ_e given by Equation (8). These equations are solved iteratively by varying ${\varphi }_{{\rm{m}}},$ until a ${\varphi }_{{\rm{m}}}$ is found that the condition

$$\begin{eqnarray}&&\mathrm{at}\;\;\bar{r}\to \infty :\quad \bar{T}\to 0\end{eqnarray} \tag{ 11 }$$

is satisfied. By repeating this iterative procedure for various λ_c, the dependence ${\varphi }_{{\rm{m}}}={\varphi }_{{\rm{m}}}({\lambda }_{{\rm{c}}})$ can be obtained numerically. At the second upstream step, Equation (7) are solved with the initial conditions given by Equation (10) in the upstream direction, $\bar{r}\lt 1,$ for fixed λ₀ and Kn₀. For Kn₀ < Kn_c, one can iteratively find such λ_c and ${\bar{r}}_{0}$ that Jeans parameter and Knudsen number at $\bar{r}={\bar{r}}_{0}$ are equal to λ₀ and Kn₀. Once the dependences ${\lambda }_{{\rm{c}}}\left({\lambda }_{0},{{Kn}}_{0}\right)$ and ${\bar{r}}_{0}\left({\lambda }_{0},{{Kn}}_{0}\right)$ are established, they can be used in order to find any solution of Parker's model for a given surface condition. Typical dependences of λ_c on λ₀ and Kn₀ are shown in Figure 1. It is remarkable that λ_c varies within a narrow range at small Kn₀ and large λ₀ for arbitrary γ and ω.

Zoom In Zoom Out Reset image size

Numerical calculations with Parker's model were performed for mon- ($\zeta =0,$ γ = 5/3) and diatomic (ζ = 2, γ = 7/5) gases with Pr = 2/3 and for ω = 3/4, and ω = 1. The mass escape rate is calculated in reduced units as

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Phi }}}_{{\rm{m}}}}{{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}}=\sqrt{2\pi }{{Ma}}_{0}\displaystyle \frac{\mathrm{exp}({\lambda }_{0})}{{\lambda }_{0}+1},\end{eqnarray} \tag{ 12 }$$

where Ma₀ and ${{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}$ are the Mach number and Jeans escape rate at the surface:

$$\begin{eqnarray}&&{{Ma}}_{0}=\displaystyle \frac{{{Kn}}_{0}}{{{Kn}}_{{\rm{c}}}}\displaystyle \frac{{\lambda }_{0}}{{\lambda }_{{\rm{c}}}}{\bar{T}}_{0}^{1-\omega },\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}={r}_{0}^{2}{\rho }_{0}\sqrt{\displaystyle \frac{{\mathfrak{R}}{T}_{0}}{2\pi }}\left({\lambda }_{0}+1\right)\mathrm{exp}(-{\lambda }_{0}).\end{eqnarray} \tag{ 14 }$$

The dependences of ${{\rm{\Phi }}}_{{\rm{m}}}/{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}$ on Kn₀ at fixed λ₀ are found to be qualitatively similar for all γ and ω considered and include two asymptotic regimes, which correspond to the case of large and small Kn₀, and intermediate regime, corresponding to transition between these asymptotic regimes (Figure 2(a)). Following Parker (1964b), the case, when ${{\rm{\Phi }}}_{{\rm{m}}}/{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}$ asymptotically does not depend on Kn₀, can be called the low-density (LD) regime and the case, when ${{\rm{\Phi }}}_{{\rm{m}}}/{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}$ asymptotically varies as 1/Kn₀, can be called the HD regime.

Zoom In Zoom Out Reset image size

The escape rates calculated with Parker's model are found to be in quantitative agreement with results of recent kinetic simulations performed by Volkov et al. (2011a, 2013) only in the HD regime (Figure 2(b)), when values of ${{\rm{\Phi }}}_{{\rm{m}}}/{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}$ calculated with both models asymptotically approach each other with decreasing Kn₀ as it is seen for λ₀ = 5 and λ₀ = 10. At λ₀ = 15, the escape rates disagree in the considered range of Kn₀, but they presumably agree at smaller Kn₀, where the results of kinetic simulations are not currently available. In the LD regime, kinetic ${{\rm{\Phi }}}_{{\rm{m}}}$ approaches ${{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}$ at ${{Kn}}_{0}\to \infty ,$ when the effects of intermolecular collisions become progressively less important and escaping flows evolve in the free molecular flow regime. On the contrary, ${{\rm{\Phi }}}_{{\rm{m}}}/{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}$ predicted by Parker's model at ${{Kn}}_{0}\to \infty$ can be both larger and many orders of magnitude smaller than 1. For instance, in the range ${\lambda }_{0}$ = 15–20, which is characteristic for a number of large Kuiper Belt objects (KBOs, Johnson et al. 2015), the maximum values of ${{\rm{\Phi }}}_{{\rm{m}}}/{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}$ are in the order of 10⁻²–10⁻⁴.

The goal of the present section is to formulate a criterion of validity of Parker's model in terms of λ₀ and Kn₀ for arbitrary γ and ω. For this purpose, accurate expressions for ${{\rm{\Phi }}}_{{\rm{m}}}/{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}$ in both HD and LD regimes are obtained first.

In the LD regime, the Knudsen number is large everywhere, including at the critical point, and ${\varphi }_{{\rm{m}}}$ is small according to Equation (9). Following Parker (1964b), one can consider the only case when ${{\rm{\Phi }}}_{{\rm{e}}}={\varphi }_{{\rm{m}}}{{\rm{\Psi }}}_{{\rm{e}}}{T}_{{\rm{c}}}{\kappa }_{{\rm{c}}}{r}_{{\rm{c}}}\gt 0$ and integrate Equation (7b) neglecting all terms in the right-hand side of this equation with exception of Ψ_e. Solution of the obtained equation satisfies Equation (11) if ${\lambda }_{{\rm{c}}}={\lambda }_{{\rm{c}}({\rm{min}})}=(2\omega +3)/(\omega +1),$ ${\varphi }_{{\rm{m}}}{{\rm{\Psi }}}_{{\rm{e}}}=1/(\omega +1),$ and temperature varies as (Parker 1964b)

$$\begin{eqnarray}&&\bar{T}={\bar{r}}^{-\frac{1}{\omega +1}}.\end{eqnarray} \tag{ 15 }$$

The critical Jeans parameter in the LD regime ${\lambda }_{{\rm{c}}(\mathrm{min})}$ corresponds to the minimum possible λ_c for a gas with given ω. If $\bar{r}\ll 1,$ then the term with ${\bar{u}}^{2}$ in the right-hand side of Equation (7a) can be neglected and gas velocity far upstream of the critical point varies as (A. N. Volkov 2015, in preparation)

$$\begin{eqnarray}&&\bar{u}={V}_{0}\left(\omega \right){\bar{r}}^{-\frac{2\omega +3}{\omega +1}}\mathrm{exp}\left[\displaystyle \frac{2\omega +3}{\omega }\left(1-{\bar{r}}^{-\frac{\omega }{\omega +1}}\right)\right],\end{eqnarray} \tag{ 16 }$$

where the coefficient ${V}_{0}(\omega )$ is given by the condition that Equation (16) matches the accurate solution of Equation (7a) at $\bar{r}\ll 1.$ The analysis of numerical solutions of Parkers model showed that V₀(ω ) varies from ∼0.708 to ∼0.6725 when ω increases from 1/2 to 1 and can be fitted to the power law ${V}_{0}(\omega )=0.672/{\omega }^{0.0742}.$ Since ${\bar{r}}_{0}={\left({\lambda }_{{\rm{c}}({\rm{min}})}/{\lambda }_{0}\right)}^{(\omega +1)/\omega },$ Equations (15) and (16) allow one to calculate the surface Mach number ${{Ma}}_{0(L)}$ and escape rate ${{\rm{\Phi }}}_{{\rm{m}}(L)}$ in the LD regime as

$$\begin{eqnarray}{{Ma}}_{0(L)} & = & {V}_{0}\left(\omega \right){\left(\displaystyle \frac{\omega +1}{2\omega +3}{\lambda }_{0}\right)}^{\frac{2\omega +5/2}{\omega }}\\ & & \times \mathrm{exp}\left(\displaystyle \frac{2\omega +3}{\omega }-\displaystyle \frac{\omega +1}{\omega }{\lambda }_{0}\right),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}\displaystyle \frac{{{\rm{\Phi }}}_{{\rm{m}}(L)}}{{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}} & = & \displaystyle \frac{\sqrt{2\pi }}{{\lambda }_{0}+1}{V}_{0}\left(\omega \right){\left(\displaystyle \frac{\omega +1}{2\omega +3}{\lambda }_{0}\right)}^{\frac{2\omega +5/2}{\omega }}\\ & & \times \mathrm{exp}\left(\displaystyle \frac{2\omega +3}{\omega }-\displaystyle \frac{{\lambda }_{0}}{\omega }\right).\end{eqnarray} \tag{ 18 }$$

In order to find an equation for the escape rate in the HD regime the definition of λ₀ can be re-written to represent the surface temperature as a function of λ₀ and Kn₀:

$$\begin{eqnarray}&&{\bar{T}}_{0}=\displaystyle \frac{{\lambda }_{{\rm{c}}}\left({\lambda }_{0},{{Kn}}_{0}\right)/{\lambda }_{0}}{\bar{r}{}_{0}\left({\lambda }_{0},{{Kn}}_{0}\right)}.\end{eqnarray} \tag{ 19 }$$

With exception of the case, when λ₀ is close to $\gamma /(\gamma -1),$ solutions of Parker's model at various λ₀ correspond to ${\lambda }_{{\rm{c}}}\approx {\lambda }_{{\rm{c}}({\rm{min}})}$ (Figure 1). Such "upstream" solutions start at large Knudsen numbers at the critical point and then, with decreasing $\bar{r},$ first follow Parker's LD approximation. Then they gradually approach another limit solution corresponding to the case of small surface Knudsen number, i.e., Parker's HD approximation. Then surface temperature and its gradient can be found as

$$\begin{eqnarray}&&{\bar{T}}_{0}=\displaystyle \frac{{\lambda }_{{\rm{c}}({\rm{min}})}/{\lambda }_{0}}{{\bar{r}}_{0}\left({\lambda }_{0},{{Kn}}_{0}\right)},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{d\bar{T}}{d\bar{r}}\right|}_{\bar{r}={\bar{r}}_{0}}\approx -\displaystyle \frac{{\lambda }_{{\rm{c}}\left({\rm{min}}\right)}/{\lambda }_{0}}{{\bar{r}}_{0}^{2}}.\end{eqnarray} \tag{ 21 }$$

Equation (7b) can be re-written at $\bar{r}={\bar{r}}_{0}$ as

$$\begin{eqnarray}&&{\bar{r}}_{0}^{2}{\bar{T}}_{0}^{\omega -1}{\left.\displaystyle \frac{d\bar{T}}{d\bar{r}}\right|}_{\bar{r}={\bar{r}}_{0}}={\varphi }_{{\rm{m}}}\left(\displaystyle \frac{\gamma }{\gamma -1}+\displaystyle \frac{{{Ma}}_{0}^{2}}{2}-{\lambda }_{0}-\displaystyle \frac{{{\rm{\Psi }}}_{{\rm{e}}}}{{\bar{T}}_{0}}\right).\end{eqnarray} \tag{ 22 }$$

By excluding the derivative from Equations (21) and (22) in the HD regime, when ${{Ma}}_{0}\ll 1$ and ${\bar{T}}_{0}\gg 1,$ one can obtain

$$\begin{eqnarray}&&{\varphi }_{{\rm{m}}}{\bar{T}}_{0}^{1-\omega }=\displaystyle \frac{{\lambda }_{{\rm{c}}({\rm{min}})}/{\lambda }_{0}}{{\lambda }_{0}-\displaystyle \frac{\gamma }{\gamma -1}}.\end{eqnarray} \tag{ 23 }$$

This equation can be used in order to calculate the surface Mach number and escape rate in the HD regime based on Equation (12) as

$$\begin{eqnarray}&&{{Ma}}_{0(H)}=\displaystyle \frac{{{Kn}}_{0}}{K}\displaystyle \frac{{\lambda }_{0}}{{\lambda }_{{\rm{c}}({\rm{min}})}}{\varphi }_{{\rm{m}}}{\bar{T}}_{0}^{1-\omega }=\displaystyle \frac{{{Kn}}_{0}}{K\left({\lambda }_{0}-\displaystyle \frac{\gamma }{\gamma -1}\right)},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Phi }}}_{{\rm{m}}(H)}}{{{\rm{\Phi }}}_{{\rm{m}}\left({\rm{J}}0\right)}}=\sqrt{2\pi }{{Ma}}_{0(H)}\displaystyle \frac{\mathrm{exp}({\lambda }_{0})}{{\lambda }_{0}+1}.\end{eqnarray} \tag{ 25 }$$

The numerical values of the escape rate closely approach the solutions given by Equations (18) and (25) as shown in Figure 3 for γ = 7/5 and ω = 1. Calculations for γ = 5/3 and/or various ω from 1/2 to 1 demonstrate similar degree of agreement between obtained equations and numerical ${{\rm{\Phi }}}_{{\rm{m}}}/{{\rm{\Phi }}}_{{\rm{m}}({\rm{J}}0)}.$ Both equations, however, become inaccurate when λ₀ closely approaches $\gamma /(\gamma -1).$

Zoom In Zoom Out Reset image size

As it was shown in Section 4, Parker's model enables quantitative prediction of Φ_m only in the HD regime. Based on Equations (17) and (24), this criterion of validity of Parker's model can be now formulated as

$$\begin{eqnarray}&&{{Ma}}_{0(H)}\ll {{Ma}}_{0(L)}.\end{eqnarray} \tag{ 26 }$$

or

$$\begin{eqnarray}&&{{Kn}}_{0}\ll {{Kn}}_{P},\end{eqnarray} \tag{ 27 }$$

where the "turning" Knudsen number Kn_P

$$\begin{eqnarray}{{Kn}}_{p} & = & K\;{V}_{0}(\omega ){\left(\displaystyle \frac{\omega +1}{2\omega +3}{\lambda }_{0}\right)}^{\displaystyle \frac{2\omega +5/2}{\omega }}\\ & & \times \mathrm{exp}\left(\displaystyle \frac{2\omega +3}{\omega }-\displaystyle \frac{\omega +1}{\omega }{\lambda }_{0}\right)\left({\lambda }_{0}-\displaystyle \frac{\gamma }{\gamma -1}\right)\end{eqnarray} \tag{ 28 }$$

corresponds to the nominal point where the dependence of the escape rate Kn₀ turns from the LD to HD regime. A comparison of Parker's solutions with results of kinetic simulations at λ₀ = 5 and λ₀ = 10 shows that Parker's escape rate is within a factor of two from the kinetic rate if ${{Kn}}_{0}/{{Kn}}_{P}\lt {10}^{-2}-{10}^{-1}.$ The turning Knudsen number, Kn_P, is a strong function of ω, but only marginally depends on γ (Figure 4). At fixed λ₀, Parker's model is valid in a substantially broader range of Kn₀ in the pseudo-Maxwellian gas (ω = 1) than in the hard sphere gas (ω = 1/2).

Zoom In Zoom Out Reset image size

The criterion of validity of Parker's model can be also formulated in terms of integral atmospheric parameters, which can be directly deduced from the analysis of absorption spectra of planetary bodies or probing the atmospheres during stellar occultation, e.g., the column density N

$$\begin{eqnarray}&&N={\int }_{{r}_{0}}^{{r}_{x}}\;\displaystyle \frac{\rho }{m}{dr}.\end{eqnarray} \tag{ 29 }$$

Here r_x is the radial position of the nominal exobase, where the mean free path of gas molecules ${l}_{{mfp}}\left({r}_{x}\right)$ is equal to the atmospheric scale height $H\left({r}_{x}\right)={r}_{x}/{\lambda }_{x}$ (${\lambda }_{x}={GM}/({\mathfrak{R}}{T}_{x}{r}_{x})$ is the exobase Jeans parameter). Assuming the barometric distribution of molecules with the scale height ${H}_{0}={r}_{0}/{\lambda }_{0},$ $\rho ={\rho }_{0}\mathrm{exp}(({r}_{0}-r)/{H}_{0}),$ the first approximation to the column density at large λ₀ is

$$\begin{eqnarray}&&N\sim {N}_{0}=\displaystyle \frac{{\rho }_{0}{H}_{0}}{m}=\displaystyle \frac{1}{{{\rm{\Sigma }}}_{T}({T}_{0})}\displaystyle \frac{1}{{\lambda }_{0}{{Kn}}_{0}},\end{eqnarray} \tag{ 30 }$$

where for the VHS model of intermolecular collisions ${{\rm{\Sigma }}}_{T}\left({T}_{0}\right)={k}_{{\rm{B}}}\sqrt{{\mathfrak{R}}{T}_{0}}/[K\kappa ({T}_{0})]=\sqrt{2}{\sigma }_{T}({T}_{0})$ and ${\sigma }_{T}({T}_{0})$ is the total collision cross section of molecules at temperature T₀ (Johnson et al. 2015) with typical value in the order of 10⁻¹⁴ cm⁻² for atmospheric species like N₂ or CH₄ (Bird 1994). The examples of dependences of N/N₀ on Kn₀ at various λ₀ obtained for a monatomic gas with ω = 1 are shown in the inset of Figure 5. It is seen that N/N₀ is independent of Kn₀ if Kn₀ < 10⁻⁵. Calculated values of N/N₀ in the limit of small Kn₀ are shown in Figure 5 by symbols. These values marginally depend on γ and ω, so that the ratio N/N₀ can be considered as a function of only λ₀ and fitted at λ₀ ≥ 4 by equation

$$\begin{eqnarray}&&\displaystyle \frac{N}{{N}_{0}}=\displaystyle \frac{1+{\lambda }_{0}+\mathrm{exp}(4/{\lambda }_{0})}{{\lambda }_{0}},\end{eqnarray} \tag{ 31 }$$

shown in Figure 5 by the solid curve. Then the criterion given by Equations (27) and (28) can be reformulated in terms of N as

$$\begin{eqnarray}&&N\gg {N}_{P}=\displaystyle \frac{1}{\sqrt{2}{\sigma }_{T}\left({T}_{0}\right)}\displaystyle \frac{1+{\lambda }_{0}+\mathrm{exp}\left(\displaystyle \frac{4}{{\lambda }_{0}}\right)}{{\lambda }_{0}^{2}}\displaystyle \frac{1}{{{Kn}}_{P}}.\end{eqnarray} \tag{ 32 }$$

Zoom In Zoom Out Reset image size

The developed criterion for the validity of Parker's model for thermal escape problems can be directly applied for atmospheres without substantial external heating. Certain KBOs, including dwarf planets in the outer solar system, may have such atmospheres. The degree of absorption of solar radiation in KBO atmospheres depends on atmospheric column (Johnson et al. 2015) and, in particular, on the abundance of methane, which is the primary absorbing agent in presumably nitrogen-rich KBO atmospheres (Strobel et al. 1996; Krasnopolsky & Cruikshank 1999). Currently the abundance of methane is not well-constrained for the majority of KBOs with exception of Pluto and Charon (see review of literature data, e.g., in Johnson et al. 2015). For Sedna and Quaoar, which have relatively thin and thick atmospheres (N ≈ 10¹⁷ cm⁻² and N ≈ 10²¹ cm⁻²) at λ ≈ 15 in solar-maximum conditions (Johnson et al. 2015), Kn₀ ≈ 7 · 10⁻⁵ and Kn₀ ≈ 7 · 10⁻⁹, while Kn_p ≈ 5 · 10⁻⁸. For KBOs like Pluto with larger λ₀, the condition given by Equation (28) is even more restrictive. Thus, Parker's model is not applicable for the majority of KBOs, which atmospheres do not contain sufficient amount of CH₄ to be heated by solar radiation.

For planetary atmospheres with external heating, the criterion given by Equation (28) can be used with certain approximations, assuming that λ₀ and Kn₀ characterize flow conditions at the top of the heated layer. These conditions are more favorable for application of Parker's model above the heated layer because of increase in the local Jeans parameter. With increasing amount of absorbed energy and decreasing Jeans parameter at the upper boundary of the heated layer, the critical point approaches the exobase, and Parker's and kinetic model can produce similar escape rates as it is apparent from the comparison of curves for λ₀ = 5 in Figure 2(b). This is in agreement with findings by Zhu et al. (2014), who concluded that Pluto's exobase is locked at λ_x = 4–6, and Parker's model escape rate is different from the enhanced Jeans rate predicted in kinetic simulations by Volkov et al. (2011a, 2011b) by the factor of 1–0.22. For less intensive heating, which can be specific to other dwarf planets, Parker's model applied to the part of the atmosphere above the heated layer can substantially underestimate the escape rate.

In conclusion, the numerical study of Parker's model reveals that two asymptotic regimes can occur depending on the surface conditions. The escape rates in these regimes are calculated theoretically in the form of Equations (18) and (25). Escape rates predicted by Parker's and kinetic models agree only in the HD asymptotic regime. The criterion of applicability of Parker's model is then given by Equations (27) and (28) or (32).

This research was supported by the NASA Planetary Atmospheres Program. The author thanks Prof. R. E. Johnson for helpful comments on the results discussed in the present paper.