Hydrodynamic Escape of Mineral Atmosphere from Hot Rocky Exoplanet. I. Model Description

Following previous models of hydrodynamically escaping atmospheres (e.g., Watson et al., 1981; García-Muñoz, 2007b; Tian et al., 2008), we integrate the equations of continuity, motion and energy for a spherically symmetric flow of an inviscid, multi-component gas:

where $t$ is time and $r$ is the radial distance from the planetary center. All the other quantities are functions of $r$ and $t$: $\rho_{s}$ and $\dot{\rho}_{s}$ are respectively the mass density and net mass production rate per unit volume of gas species $s$, $\rho$ is the total mass density (i.e., $\rho=\sum_{s\in\mathscr{S}}\rho_{s}$ $=\sum_{s\in\mathscr{S}}m_{s}n_{s}$, where $m_{s}$ and $n_{s}$ are the mass and number density of species $s$, respectively, and $\mathscr{S}$ represents the set of gas species, namely, $\mathscr{S}$ = Na, Na⁺, Na²⁺, etc.), $u$ is the bulk velocity, $u_{s}$ is the diffusion velocity of species $s$ (or the component of a velocity difference caused by diffusion), $P$ and $E$ are the pressure and the total energy density (or the sum of the kinetic energy density and the specific internal energy of the bulk gas flow), respectively, $f_{\rm{ext}}$ is the external force, $q$ is the heat flux, and $Q_{\rm{net}}$ is the net energy deposition rate.

We assume that the atmosphere consists of perfect monoatomic gases. The total energy density is given by

and the sound speed is $c_{S}=\sqrt{\gamma{P}/{\rho}}$ with the heat capacity ratio $\gamma=$ 5/3.

The external force $f_{\rm{ext}}$ is given by

where $a$ is the separation between the planetary and stellar centers, $M_{p}$ and $M_{*}$ are the planetary and stellar masses, respectively. The first, second, and third terms on the right-hand side represent the planetary gravity, the stellar gravity, and the centrifugal force, respectively. The last two terms are collectively termed the tidal force.

Finally, the quantities $\dot{\rho_{s}}$, $Q_{\rm{net}}$, $u_{s}$ and $q$ are determined from the models of chemical reaction, heating/cooling processes and diffusion/thermal conduction, which are described below.

Since the atmosphere considered in this study is composed of atoms and their ions, ionization and recombination dominate the atmospheric chemistry. Here, we consider reactions relevant both to radiative and thermal processes. All those processes except inverse ones are listed in Table 3: Such inverse processes are the thermal recombination of each ion (i.e., three-body recombination opposite to the thermal ionization reactions named TI 1–16 in Table 3). For the thermal recombination, we assume an electron as the third body that receives the reaction heat and calculate the reaction rate constant of thermal recombination from that of thermal ionization and the equilibrium constant given by the Saha’s ionization equation (Eq. 9.45 in Rybicki & Lightman, 1986). We ignore direct exchange of electric charges between atoms and ions for simplicity, since those reaction rates are unavailable for almost all of the species considered in this model. Note that, instead of the charge exchange, the thermal ionization and recombination reactions of each atom or ion with electrons lead to removing and adding their electrons, and then these two processes consequently exchange the charge between different atoms and ions in this model.

Consider a photo-ionization process such that an $i$-times charged ion of chemical element $Y$, which is denoted by $Y^{i}$, leads to a further loss of $x$ electrons and produces an ion with an $i+x$ charge (i.e., $Y^{i+x}$). Using the symbol $n(Y^{i})$, instead of $n_{s}$, the rate of photo-ionization from $Y^{i}$ to $Y^{i+x}$ is given by

where $F^{*}_{\lambda}$ is the emergent photon flux from the host star, $\tau_{\lambda}$ is the optical depth at wavelength $\lambda$ measured from infinity, $\sigma_{Y^{i},\lambda}$ is the photo-ionization cross section of species $Y^{i}$, $\eta_{i\rightarrow(i+x),Y,\lambda}$ is the probability that a collision with a photon results in removing $x$ electrons from $Y^{i}$ and yielding $Y^{i+x}$.

As for the stellar spectrum $F^{*}_{\lambda}$, we adopt the solar-type star models that Claire et al. (2012) developed by combining the observed spectra of the Sun and solar analogs of different ages. Three model spectra of different stellar ages 0.1, 1, and 4.56 Gyr are presented in Fig. 2. The photo-ionization cross section $\sigma_{Y^{i},\lambda}$ is given by Verner & Yakovlev (1995). As for $\eta_{i\rightarrow(i+x),Y,\lambda}$, we consider single and multiple electron emission that occurs when X+UV removes an electron from an atom and then such electron emission from an inner-shell of the atom brings about further emission of electrons (which is termed the Auger effect). Since we consider only up to 4$+$ ions in this study, we use the form of $\eta_{i\rightarrow(i+x),Y,\lambda}$ given by Kaastra & Mewe (1993) for $i+x\leq$ 4 and assume $\eta_{i\rightarrow(i+x),Y,\lambda}$ = 0 for $i+x\geq 5$. Note that the function is normalized so that the sum of $\eta_{i\rightarrow(i+x),Y,\lambda}$ is unity.

To integrate Eq. (7), we consider 23 spectral intervals in the wavelength region between 0.1 and 250 nm in such a way

In each bin of width $\Delta\lambda_{k}$ ($\equiv\lambda_{k+1}-\lambda_{k}$), we calculate the averaged value of the photo-ionization cross section, $\bar{\sigma}_{s,k}$, as

where $\langle F_{\lambda}^{\ast}\Delta\lambda\rangle_{k}$ is the photon flux in the $k$th bin given by

Note that we simply use the value of $F^{*}_{\lambda}$ calculated at $\lambda$ = 0.5 nm over wavelength between 0.1 and 0.5 nm, where the data of $F^{*}_{\lambda}$ are unavailable in Claire et al. (2012). Likewise, the averaged photon energy at the $k$th bin is given by

where $c$ is the light velocity and $h$ is the Planck constant. Note that while we use the above averaged values, $\bar{\sigma}_{s,k}$ and $\bar{\lambda}_{k}$, at each bin below, we do not indicate so for the sake of shorthand.

Finally, the net mass production rate, $\dot{\rho_{s}}$ in Eq. (2), is calculated from reactions tabulated in Table 3. For instance, in the case of a reaction of two-body reactants $s_{1}$ and $s_{2}$ and a single product $s_{3}$ with a rate coefficient $k_{r}$, the mass production rate of species $s_{3}$ is given by $\dot{\rho}_{s_{3}}=m_{s_{3}}k_{r}n_{s_{1}}n_{s_{2}}$.

In this model, we consider the absorption of stellar X+UV by atoms and ions via photo-ionization. The absorbed energy is used for removing the outermost electrons (i.e., normal ionization) and also the inner-shell electrons, which induces characteristic X-ray emission, and is also partitioned into the thermal energy. In addition, we consider chemical reactions and radiative transitions as other heating and cooling processes. Thus, the net energy deposition rate, $Q_{\rm{net}}$, is given by

where $Q_{\rm{X+UV}}$ is the total energy of X+UV absorbed per unit time given by

$Q_{\rm{Ei}}$ is the energy loss rate via photo-ionization, which is given by the sum of the photo-ionization rates given in Eq. (7) times corresponding ionization energies; $Q_{\mathrm{X}}$ is the energy loss rate via characteristic X-ray emission, the probability of which is calculated from Kaastra & Mewe (1993) with the assumption that the characteristic X-ray is lost completely from the atmosphere to space, and, thus

where $\eta_{\rm{X,\lambda}}^{k}$ and $E_{\rm{X,\lambda}}^{k}$ are the emission probability and energy of the $k$th characteristic X-ray (see Kaastra & Mewe, 1993); $Q_{\rm{chem}}$ is the endothermic/exothermic rate of thermo-chemical reactions (i.e., reactions except photo-ionization and radiative recombination), whose ionization energy data are taken from the NIST database²²2https://physics.nist.gov/PhysRefData/ASD/ionEnergy.html; $Q_{\rm{rad}}$ is the net radiative energy absorbed/emitted per time by energy level transition.

We define the primary X+UV heating rate, $Q^{\prime}_{\rm{X+UV}}$, as

and, hereafter, use this quantity, instead of $Q_{\mathrm{X+UV}}$, namely,

We consider multi-component diffusion, following the formula described by García-Muñoz (2007a). This formula is based on the momentum equations for multi-component gas derived by Burgers (1969) and arranged in the classical Stefan–Maxwell form. In this study, we adopt the first order approximation of this method for the mineral atmosphere in ambipolar constraint, which is also used in hydrodynamic simulations for hot Jupiters (García-Muñoz, 2007b). For calculating the diffusion velocity of species $s$, $u_{s}$, from the approximate diffusion matrix (see Appendix A in García-Muñoz, 2007b, for the details), one needs the binary diffusion coefficient between species $i$ and $k$, $\mathscr{D}_{i,k}$. Following García-Muñoz (2007a), for the atom-atom and atom-electron binary diffusions, we calculate $\mathscr{D}_{i,k}$ from the hard sphere model as

where $\bar{m}_{i,k}$ = $m_{i}m_{k}/(m_{i}+m_{k})$. For the atom-ion binary diffusion, the interaction between ions and atoms is assumed to be due to the induced polarization potential and $\mathscr{D}_{i,k}$ is given by

where $\alpha_{n}$ is the polarizability of the atom and $Z_{i}$ is the charge number for the ion. The values of $\alpha_{n}$ are taken from CRC Handbook (2011). For the ion-ion binary diffusion, we use the formula of the Coulomb interaction as

where $\beta$ is given as ${\beta}=1.26\times 10^{4}(T^{3}/n_{e})^{1/2}$ (Burgers, 1969). For the ion-electron binary diffusion, following Schunk & Nagy (2000), we calculate $\mathscr{D}_{i,k}$ as

Finally, we calculate $u_{s}$ from Eqs. (A3), (A4) and (A7)–(A12) of García-Muñoz (2007a) with the binary diffusion coefficients described above.

We assume that heat transport occurs via thermal conduction and heat exchange due to chemical diffusion. The heat flux $q$ is calculated as (García-Muñoz, 2007b)

where $\kappa$ is the thermal conductivity and $h_{s}$ is the specific enthalpy of species $s$ given by $h_{s}$ = $\gamma k_{b}T/[(\gamma-1)m_{s}]$. $\kappa$ is given as a sum of the thermal conductivities of species $s$, $\kappa_{s}$, namely

for neutrals, $\kappa_{s}$ is given from Eqs. (14.30) and (14.42) of Banks & Kockarts (1973) as

where $d_{s}$ is the atomic diameter for which we adopt the van der Waals size; for ions, $\kappa_{s}$ is given from Eq. (22.105) of Banks & Kockarts (1973) as

where $m_{p}$ is the proton mass; for electrons, $\kappa_{s}$ is given from Eq (22.122) of Banks & Kockarts (1973) as

For simulating the hydrodynamic flow, we integrate Eq. (2) for $S$ gas species and Eqs. (3)–(4). Thus, we need $S+2$ boundary conditions.

The lower boundary of the escaping atmosphere is defined as a spherical surface above which the incident stellar X+UV is completely absorbed. To be exact, the optical depths at all the wavelength bins in the X+UV are larger than ten at the lower boundary. Since the hydrostatic equilibrium approximation is valid in such a deep atmosphere, we determine the lower-boundary temperature and mass density from our hydrostatic model of the mineral atmosphere developed in Ito et al. (2015). Figure 3 shows an example of the temperature-pressure profile in the hydrostatic atmosphere at the sub-stellar point of the planet with $T_{\mathrm{eq}}$ of $3000$ K, which corresponds to 0.02 AU around a Sun-like star (see Eq. (16) in Ito et al., 2015, $A_{p}=0$).

Also, the hydrostatic model provides the molar fractions of neutral atomic species at the lower boundary. For instance, in the above case, the molar fractions of Na, O, Si and Mg are 0.375, 0.484, 0.127 and 0.0139. At the lower boundary, the molar fractions of all the ionic species are set to zero. In addition, we set the velocity at the lower boundary, $u_{1}$, to the value extrapolated from the upper atmosphere at each time step in the following way:

where $R_{p}$ and $R_{o}$ the radial distances of the lower and upper boundaries from the planet’s center, respectively, and $\rho_{1}$ is the density at $r$ = $R_{p}$. This extrapolation prevents the mass flux from depending on the lower boundary condition.

The upper boundary conditions are less important in this study, since we focus on an escaping atmosphere with a super-sonic velocity, which means no information of the upper boundary is propagated to regions with subsonic velocities where the escaping flow is driven. Thus, we set simply the quantities such as $n_{s}$, $T$ and $u$ at the upper boundary to values extrapolated from the interior of the computational domain, following García-Muñoz (2007b). We set the upper boundary radius to 10 $R_{p}$.

To find the steady state of the hydrodynamic flow, we perform time integration of the conservative forms of Eqs. (2)–(4) based on spatial- and time-discretization:

where

$\Delta r$ is the size of the non-overlapping spatial cell at $r$, and

Using these expressions, we explain the spatial discretization and the numerical schemes for the advection term and marching time, below.

The computational region of the atmosphere is divided into $N_{r}$ layers. The thickness of each layer is assumed to increase with altitude in such a way that the thickness ratio between two neighboring layers, $l_{r}$, is constant. Following García-Muñoz (2007b), we use $N_{r}=600$ and $l_{r}=1.014$. Also, all the variables in each cell are defined at the center of each cell. For $\Delta F_{d}$, we calculate $F_{d}$ at the cell boundaries using the values of variables at the centers of the neighboring cells. For $\Delta F$, we use the finite volume method with second-order accuracy.

The above spatial discretization sometimes causes a numerical instability during time integration of $\bm{\Delta F}$, because of unattenuated short-wavelength numerical disturbances (see also Hirsch, 1990). To ensure the numerical stability, we add an artificial dissipation flux $\bm{F_{\rm{AD}}}$ to $\bm{F}$ as (Swanson & Turkel, 1992; Swanson et al., 1998)

where $\triangle$ and $\nabla$ are the forward and backward spatial difference operators, respectively, $\zeta$ is a constant (= 0.01 in this study), $\bm{U_{AD}}$ is given by

and $\bm{A^{\prime}}$ is a square matrix of $S+2$ rows and columns which is given by

$\bm{\Re}$ is the eigenmatrix of $\bm{\partial F}/\bm{\partial U}$ defined as

where $\bm{\hat{\lambda}}$ is the eigenvalue matrix, which is a diagonal one with $S+2$ rows and columns given by

The three dots mean that all the elements except the first one ($\hat{\lambda}_{1}^{\prime}$) and the final one ($\hat{\lambda}_{3}^{\prime}$) are $\hat{\lambda}_{2}^{\prime}$.

In this matrix,

with reduction constants $V_{n}$ and $V_{l}$, and $\hat{\lambda}_{{j}}$ is the eigenvalue of $\bm{\partial F}/\bm{\partial U}$. Here we use $V_{n}=0.25$ and $V_{l}=0.025$, which are their typical values (Swanson & Turkel, 1992). From Eqs. (44) and (50), $\hat{\lambda}_{1}=u-c_{S}$, $\hat{\lambda}_{2}=u$, and $\hat{\lambda}_{3}=u+c_{S}$. Note that when $\bm{F_{\rm{AD}}}$ makes a major contribution to $\bm{F^{\prime}}$ throughout the atmosphere, this computation may be invalid. Basically, the steeper the gradients of physical quantities in the atmosphere are, the larger $\bm{F_{\rm{AD}}}$ is, because $\bm{F_{\rm{AD}}}$ includes the third order differential of $\bm{U_{AD}}$. We have, however, made sure that the contribution of $\bm{F_{\mathrm{AD}}}$ to $\bm{F}^{\prime}$ is less than 1 % in the simulation region above 1.1 $R_{p}$ where is the atmospheric region of interest in this study.

Equation (40) includes chemical reaction terms, $\mathcal{L_{R}}$, which differ greatly in magnitude (i.e., mathematically stiff for stable time-integration). To time-integrate such a system of stiff differential equations, we use the implicit-explicit Runge-Kutta method with six-step and fourth-order accuracy named ARK4(3)6L[2]SA, which was developed by Kennedy & Carpenter (2003). Treating $\mathcal{L_{R}}$ implicitly and $\mathcal{L_{H}}$ explicitly, we calculate the time marching for $\mathbf{U}$ by the integration of Eq. (40). Then, for the user specific tolerance error, we set that the relative error for all elements of $\mathbf{U}$ is $10^{-4}$ and the absolute error for $\rho_{s}$ is $10^{-10}\rho$. Also, we use the variable time step as (Kennedy & Carpenter, 2003)

where $\Delta t_{\rm CFL}$ is the time step limited by the Courant-Friedrichs-Lewy (hereafter, CFL) condition, and $\Delta t_{\rm PID}$ is the time step controlled by Proportional-Integral-Differential controllers. $\Delta t_{\rm CFL}$ is given by

where $\zeta_{\rm CFL}$ is the CFL number that equals 2.01. And, $\Delta t_{\rm PID}$ is given from Eq. (30) of Kennedy & Carpenter (2003) as

where $\jmath$ is the user specific tolerance error, $p$ is the value that equals the order of accuracy minus one, thus three in this model, and $\delta^{\left(n+1\right)}$ is the numerical error of $\bf{U}(t+\Delta t^{\left(n\right)})$. The error is estimated by the difference between the fourth-order accuracy values and the third-order accuracy values (see Kennedy & Carpenter, 2003, for the details).

This numerical method is similar to that used in García-Muñoz (2007b). The main difference between their and our numerical schemes is whether the artificial dissipative terms of the discrete equations are scaled by the eigenvalues of the Jacobian matrix or by the spectral radius (see Swanson & Turkel, 1992, for the details). The time integration based on a Runge-Kutta method is also used in García-Muñoz (2007b). Following García-Muñoz (2007b), we decide that the system has converged to a steady state, when the relative changes in $\rho$, $\rho u$ and $\rho E$ averaged over the $N_{r}$ layers are below $10^{-6}$ for 1000 seconds. It is also required that the sum of the mass flux terms, including those of the artificial dissipation, from the continuity equations over all species fluctuate within 1 % throughout the calculated region.

For an isothermal transonic atmosphere, the velocity is analytically given by (Parker, 1964)

where

which is termed the critical radius, and $\bar{m}$ is the mean mass per gaseous particle. Our numerical solution is compared with Eq. (58) in Fig. 4 for a super-Earth of mass 10 $M_{\oplus}$ and radius 2 $R_{\oplus}$ having an atmosphere composed only of atomic hydrogen whose temperature is 5000 K. Here we have ignored the tidal force, namely, the last two terms of Eq. (6). Figure 4 shows that the numerically calculated velocity reproduces the analytical solution with enough accuracy (at most one percent uncertainty of sound velocity).

García-Muñoz (2007b) performed numerical simulations of hydrodynamic escape of a highly UV-irradiated hydrogen-dominated atmosphere, supposing the hot Jupiter orbiting at 0.047 AU, HD 209458 b (see also Yelle, 2004). With the same elementary processes and conditions as those considered in García-Muñoz (2007b), we calculate the atmospheric structure. Note that we use the fitting polynomial function of temperature (Koskinen et al., 2007) for the radiative cooling rate of H³⁺ shown in Neale et al. (1996), although García-Muñoz (2007b) gives no clear description of how to use the cooling rate.

In Figs. 5 and 6, we reproduce some of the results of García-Muñoz (2007b): The former shows the velocity, pressure and temperature profiles, which can be compared directly with the previous study’s results (dashed lines) that are shown in Figs. 3 and 4 of García-Muñoz (2007b). Here we consider two cases without and with the tidal effect; the former and latter are referred to as the "SP" case and the "SP-TF" case, respectively. Figure 6 shows the radial distributions of the number densities of five gaseous species (top panel) and the heating and cooling rate (bottom panel) calculated in the SP case; The latter panel is shown in the same form as Fig. 2 of García-Muñoz (2007b).

Our results are quite similar to those of García-Muñoz (2007b), although small differences have arisen due to differences in the numerical scheme and the fitting function for the radiative cooling rate of H³⁺. For instance, in the SP case, the sonic point is located at $r$ = 4.6 ${R_{p}}$ in our model, but 4.3 $R_{p}$ in García-Muñoz (2007b). Also, our model estimates the escaping mass flux to be 1.42 $\times 10^{11}$ g/s and 1.47 $\times 10^{11}$ g/s in the SP and SP-TF cases, respectively, while the estimates by García-Muñoz (2007b) are 1.48 $\times 10^{11}$ g/s and 1.52 $\times 10^{11}$ g/s, respectively; that is, the differences are less than 5 %. Thus, it would be fair to say that our model is in good agreement with García-Muñoz (2007b).

Figure 7$a$ shows the radial profiles of mass density $\rho$ (black), temperature $T$ (orange) and velocity $u$ (cyan) of the atmospheric gas in the case of the 0.1-Gyr-old host star. Also, Fig. 8 shows the mass/number density profiles of singly- to quadruply-charged ions for ($a$) all the species, ($b$) Na, ($c$) O, ($d$) Mg and ($e$) Si. First, as found in Fig. 7$a$, the velocity increases with altitude and then exceeds the local speed of sound at $r\simeq$ 3.2 $R_{p}$, the point of which is termed the sonic point. Thus, a transonic, hydrodynamic escape of the atmosphere turns out to occur.

The mass density decreases with altitude rapidly below $r$ = 1.6 $R_{p}$, but slowly above that altitude, indicating that the pressure scale height increases rapidly from that altitude on. This can be readily understood from the fact that temperature increases greatly beyond $r$ = 1.6 $R_{p}$ to $1.8\times 10^{4}$ K at $r$ = 10 $R_{p}$. Also, as shown in Fig. 8$a$, all the neutrals are photo-ionized almost completely at $r\simeq 1.6R_{p}$ by the X+UV irradiation and, then, electrons dominate in number for $r>$ 1.6 $R_{p}$. This reduces the mean mass of a gas particle $\bar{m}$, leading to an increase in the pressure scale height. For instance, the scale height is larger approximately by a factor of fifteen at $r=1.6R_{p}$ than that at $r=1R_{p}$, since the mean mass decreases to one-third of the value at $r=1R_{p}$ by ionization, the temperature doubles and the gravity approximately decreases to $1/1.6^{2}$ at $r$ $\simeq$ $1.6R_{p}$.

Above $r$ = 1.6 $R_{p}$, gas density is too low for recombination to occur efficiently. Thus, regardless of element, all the atoms are completely ionized in the upper atmosphere, as shown in Fig. 8$b$–$e$. Also, the similarity in those profiles indicates that gravitational separation is ineffective. For instance, the fraction of the lightest element O increases by at most 1 % throughout the atmosphere. This is because the bulk velocity of gas is much higher than the diffusion velocity of each species and these ionic gases diffuse along with electrons due to the ambipolar electric force.

Note that a hydrodynamic description is invalid for such a low-density gas that collisions between gaseous particles occur on timescales longer than hydrodynamic one. In other words, results shown in this study are invalid above the exobase of the atmosphere, the altitude $r_{\mathrm{exo}}$ of which is defined as

In the present case, $r_{\rm{exo}}$ is estimated to be $\sim$ 9.8 $R_{p}$ (grey diamond in Fig. 7$a$), using the typical collision cross section $\sigma=3\times 10^{-15}$ cm² and the calculated number densities of ions and neutrals. Since the sonic point ($r\sim$ 3.2 $R_{p}$) is located well below the exobase, the hydrodynamic description for the lower atmosphere is valid. Although the hydrodynamic description is invalid for $r\gtrsim 9.8R_{p}$, the structure of the upper atmosphere has no influence on that of the lower atmosphere, because the fluid motion is already supersonic above the exobase. We also have confirmed that hydrodynamic description is valid in the other simulations of this study.

Figure 9 shows the energy budget in the mineral atmosphere: Panel ($a$) shows the heating and cooling rates (represented by solid and dashed lines, respectively), which include $Q_{\rm{net}}$, $Q_{\rm{X+UV}}$, $Q^{\prime}_{\rm{X+UV}}$, $Q_{\rm{chem}}$ and $Q_{\rm{rad}}$; panel ($b$) shows the contribution of the effective coolants, Na, Mg⁺, Si²⁺, Na³⁺ and Si³⁺to $Q_{\rm{rad}}$. The radial profiles of those heating/cooling rates are found to be qualitatively different from each other, while all the rates decrease with altitude because of the decrease in gas density (see Fig. 7$a$). The stellar X+UV is completely absorbed in the atmosphere without reaching the lower boundary ($r$ = 1 $R_{p}$) and energy deposition (i.e., heating) occurs mainly below $r$ = 1.6 $R_{p}$ (see the pink line, $Q_{\rm{X+UV}}$, and the red line, $Q^{\prime}_{\rm{X+UV}}$). The net effect of the thermo-chemical reactions ($Q_{\rm{chem}}$; green dashed line) is cooling everywhere in the atmosphere and the net cooling rate decreases rapidly with altitude. Because thermal ionization, which converts thermal energy to ionization energy, dominates over thermal recombination, it acts as a net cooling. Transition between energy levels occurs by absorption of external radiation (excitation) and radiative emission (de-excitation). Its net effect turns out to be heating in the lower atmosphere below $r\sim 1.08R_{p}$ and cooling for $r>1.08R_{p}$ (see $Q_{\rm{rad}}$ represented by the cyan solid and dashed lines, respectively, in the inset of Fig. 9$a$).

Atmospheric escape is controlled by energy budget. As shown in Fig. 9$a$, the primary X+UV heating, $Q^{\prime}_{\rm{X+UV}}$, is almost in balance with the radiative cooling, $Q_{\rm{rad}}$, which is due mainly to line emission, in the region between $r$ $\sim 1.2R_{p}$ and $\sim 7.5R_{p}$. Consequently, the net energy deposition rate, $Q_{\mathrm{net}}$, is much smaller than $Q_{\mathrm{X+UV}}^{\prime}$ and $Q_{\mathrm{rad}}$. (Note that the radiative heating due mainly to external-radiation absorption is almost balanced with the energy consumption via chemical reactions, $Q_{\rm{chem}}$, below $\sim 1.08R_{p}$.) This indicates that advective energy transfer is ineffective below $7.5R_{p}$. To be more specific, the Na-D emission makes a major contribution to the radiative cooling below 1.6 $R_{p}$ (see the navy line in Fig. 9$b$). For $r>$ 1.6 $R_{p}$, Na is rapidly photo-ionized and, then, the energy budget becomes dominated by emission due to energy-level transition of outermost electrons of multiply-charged ions such as Mg⁺, Si²⁺, Na³⁺ and Si³⁺, as shown in Fig. 9$b$. Above $7.5R_{p}$, instead, $Q_{\rm{net}}$ is larger than $Q_{\rm{rad}}$, indicating that the absorbed X+UV energy is used to drive the advection (i.e., escape) of the atmosphere.

We define the net X+UV heating efficiency, $\epsilon_{\mathrm{X+UV}}$, as

where $r_{s}$ is the radial distance of the sonic point. Here we have considered only the region below $r=r_{s}$ since information cannot propagate from supersonic regions down to subsonic ones. Integrating Eq. (62) over the range where the X+UV heating rate is higher than the line absorption rate (i.e., $Q_{\rm{X+UV}}\geq Q_{\rm{rad}}$), we estimate that $\epsilon_{\mathrm{X+UV}}$ $\simeq 1.6\times 10^{-3}$. The mass loss rate comes out to be 0.3 $M_{\oplus}\,\mathrm{Gyr}^{-1}$ or $5.7\times 10^{10}$ g s^-1.

Stellar age affects the structure and escape of the atmosphere, since the stellar emission spectrum including X-ray and UV ($\lesssim$ 250 nm) differs with age (see Fig. 2). Figure 7$b$ is the same as Fig. 7$a$, but for the stellar age of 1 Gyr. Comparing both panels, one finds that each quantity undergoes a qualitatively similar change, although being slightly smaller as a whole in the old-star case than in the young-star case. Such similarity and small difference can be understood from the energy budget as follows.

Figure 10 shows ($a$) the mass density profiles of all the neutrals and ions and ($b$) the energy budget in the atmosphere. Comparing Fig. 10$a$ with Fig. 8$a$, one finds that the ionization degree is slightly lower in the old-star case than in the young-star case. For example, at $r=3R_{p}$, the major species is secondary-charged ions in the former case, while being thirdly-charged ones in the latter case. Such a difference in ionization degree is due simply to that in the incident X+UV flux, which can ionize those atoms. The older star emits $\sim$10 times lower flux of X+EUV than the younger star, whereas the fluxes of FUV and visible light scarcely differ between the two stars.

The mass loss rate $\dot{m}$ is estimated to be $3.7\times 10^{-2}$ $M_{\oplus}$ Gyr^-1 in the old-star case, which is approximately ten times lower than that in the young-star case. Also, the net X+UV heating efficiency, $\epsilon_{\mathrm{{X+UV}}}$, decreases from $1.6\times 10^{-3}$ in the young-star case to $2.2\times 10^{-4}$ in the old-star case. At first glance, it may sound strange that the X+UV heating efficiency is greatly changed despite of similarity in the dominant heating/cooling processes (compare Figs. 9$a$ and 10$b$) and in the X+UV irradiation energy integrated over $\lambda\lesssim$ 250 nm ($\simeq 4.9\times 10^{6}$ erg cm^-2 s^-1 for the younger case and $\simeq 3.4\times 10^{6}$ erg cm^-2 s^-1 for the older case, Fig. 2). Such an interpretation is, however, incorrect. Figure 11 shows the ratio of the X+UV heating rate to the net energy deposition rate ($Q_{\rm{net}}$), equivalent to the inverse of the local heating efficiency, where we include the contributions of respective wavelength regions $\lambda$ = 0.1–25 nm (“X-ray”; dashed light-green), 25–125 nm (“EUV”; dashed brown) and 125–250 nm (“FUV”; dashed dark-green) for the young-star ($a$) and old-star ($b$) cases. Note that we show the inverse of the local heating efficiency in order to clarify the contributions of respective wavelength regions. The inverse of the local heating efficiencies for X-ray and EUV wavelengths change with stellar age only slightly, whereas the peak value for FUV wavelengths is lower by a factor of 10 in the old-star case than in the young-star case. In addition, the X-ray contribution to the inverse of the local heating efficiency is lower than for both EUV and FUV for both stellar ages. Therefore, the X-ray and EUV make the dominant contribution to the atmospheric escape, while FUV accounts for a large proportion of $F_{\mathrm{X+UV}}$; for example, in the young-star case, $\int F_{\mathrm{X+UV}}d\lambda\simeq 3.9\times 10^{5}$, $1.4\times 10^{6}$, and $3.1\times 10^{6}$ erg cm^-2 s^-1 for X-ray, EUV, and FUV, respectively.

To gain a deeper understanding of the sensitivity of mass loss rate to stellar spectrum, we perform the simulations by artificially raising or lowering the young-star emission flux by 10-fold in specific wavelength regions (X-ray, EUV, or FUV). Table 2 shows the calculated mass loss rates for the different cases. The mass loss rate is found to be relatively insensitive to difference in FUV. This means that the X-ray and EUV effectively drive the atmospheric motion rather than FUV. Thus, the mass loss rate in the old-star case is approximately ten times lower than that in the young-star case because the X-ray and EUV fluxes are lower. Note that although FUV is less important for the mass loss rates, the difference in the net X+UV heating efficiencies, $\epsilon_{\mathrm{{X+UV}}}$, between the young and old cases is due to that in FUV heating efficiency because $F_{\mathrm{X+UV}}$ accounts for a large proportion of FUV. In Section 5.1.2, we discuss why such a short-wavelength radiation dominates the escape of the mineral atmosphere.