ATMOSPHERE EXPANSION AND MASS LOSS OF CLOSE-ORBIT GIANT EXOPLANETS HEATED BY STELLAR XUV. I. MODELING OF HYDRODYNAMIC ESCAPE OF UPPER ATMOSPHERIC MATERIAL

In the first studies of the planetary atmosphere expansion aimed at estimating the material escape rate and related mass loss, an isothermal Parker-type solution was imposed with an a priori prescribed boundary condition at the planetary inner atmosphere (or surface) related to the XUV heating of a fixed, predefined thin layer (Watson et al. 1981). However, it was shown in Tian et al. (2005) that in a thin-layer-heating approximation the expanding atmosphere mass loss rate strongly depends on the altitude of the layer where the radiation energy is deposited.

Moreover, as was suggested already by Parker, the realistic solution that corresponds to an outflowing (expanding) wind regime requires an additional volume heating in order to compensate for the adiabatic cooling of the expanding gaseous envelope. This was also confirmed in the follow-up numerical modeling studies. For example, in the case of the solar wind, an additional acceleration factor is called for to appear at a few solar radii. This additional driver is attributed to the dissipation of Alfvén waves (Usmanov et al. 2011). For an outflowing atmosphere of a hot Jupiter, one of the possible candidates for the role of such a distributed accelerator is the absorption of XUV, which takes place everywhere where atomic or molecular hydrogen is present. However, it turns out that at distances of about several planetary radii, where the additional forcing of the expanding planetary atmospheric material is needed to keep it in the outflowing regime (i.e., to escape the gravity tension), the gas is highly ionized, and the absorption of XUV decreases.

Nevertheless, the same fact of the increasing ionization degree, which reduces the efficiency of XUV heating, contributes to the additional acceleration of the wind and appears as an additional booster of the planetary plasma wind. That is because the thermal pressure p = (n_a + n_i + n_e)kT increases (up to two times) at certain intervals of heights (in the ionization region) because of the contribution of the electron and ion pressure, whereas the density of material ≈(n_a + n_i)m_i remains practically unchanged. The created pressure gradient provides an additional driving force for the expansion of ionized atmospheric material. Here we use the standard notations with n_a, n_i, n_e, m_i, k, and T standing for the number density of neutrals, ions, electrons, ion mass, Boltzmann constant, and temperature, respectively. In the course of the ionization, the expanding atmospheric wind experiences in the ionization region an additional acceleration and reaches velocities that sustain the continuous outflow regime. Note that the adiabatic cooling and the plasma pressure gradient mechanism caused by the ionization act opposite to one another, and the simple isothermal model that does not consider them properly produces an erroneously too strong planetary wind. To our knowledge, the role of the ionization layer and acceleration of the atmospheric plasma by the pressure gradient were not elucidated thus far for the problem of a hot Jupiter plasma wind generation. At the same time, this effect is implicitly included in several reported simulations (e.g., García Muñoz 2007; Guo 2011, 2013; Koskinen et al. 2010, 2013).

Properly defined boundary conditions at the inner edge of the simulation box appear to be a crucial issue for the models of atmosphere expansion driven by XUV heating. These boundary conditions are usually associated with a deep layer of inner planetary atmosphere, which we conventionally call a planet surface.

Besides the plasma outflow and the related adiabatic cooling, because of plasma expansion, another efficient energy drain that acts to compensate for the XUV heating is an excitation of the second level of hydrogen atoms by electron impact and fast emission of Lyα photons (Murray-Clay et al. 2009). At the same time, because of the exponential decrease with temperature, this mechanism can fully compensate for XUV heating of the deep and cold layers of the thermosphere if the latter decreases with height fast enough. However, because of the wide spectrum of stellar XUV illuminating a planet, its heating rate decreases with height (i.e., with the increase of integrated column density), not exponentially, but more likely as a power law. Therefore, a model that is based only on the Lyα cooling requires a limiting inner boundary at the layer with a temperature of several thousand Kelvin. To admit this condition, one needs an additional justification (i.e., model) for the inner boundary at the lower thermosphere. An estimation of the overall energy budget shows that the inner "cold" thermosphere has to have a temperature of about 1000–2000 K. The cooling mechanism due to infrared emission that involves the H₂ and ${\rm H}_{3}^ +$ species is usually addressed in that respect (Yelle 2004). Although complex models that include the major photochemistry effects already exist (Yelle 2004; García Muñoz 2007; Koskinen et al 2007), in the present paper we adopt a bit different and more simple approach to take into account the infrared cooling effects.

While infrared cooling due to ${\rm H}_{3}^ +$ may dominate in the deep thermosphere at pressures larger than ∼1 μbar, it has been concluded in Koskinen et al. (2013), based on a comparison of kinetic and hydrodynamic models, that this process does not affect formation of the expanding planetary wind. In that respect, Koskinen et al. (2013) used a hydrodynamic model with an inner boundary fixed at the level of 1 μbar pressure. On the other hand, regardless of Lyα or infrared cooling, the expanding plasma motion alone, because of the adiabatic cooling mechanism, can effectively remove heat from the deep layers of a thermosphere affected by XUV radiation. In some models this effect is known to be so strong that it eventually overcomes the heating and leads to an unrealistic drop in temperature (Watson et al. 1981). However, in a case when radiation (integral over spectrum) attenuates not exponentially but more likely as the inverse power of the density column, the domination by the adiabatic cooling effect is diminished. Below, using the results of an analytical solution that takes into account only adiabatic cooling and a proxy model of XUV heating at large integrated column densities, as well as by means of the numerical simulation, we show that the inner boundary can be extended down to the deeper layers with a much higher pressure on the order of $\tilde{P}_0 = p_{\max } = 10\,{\rm \mu bar}$. García Muñoz (2007) also arrives at a similar conclusion.

The question of the inner boundary pressure appears to be a crucial one also because, according to Tian et al. (2005) and Trammell et al. (2011), the material escape at low pressures is practically proportional to it. This has to be avoided in a self-consistent model. Therefore, taking the inner boundary sufficiently deep in the atmosphere, where the XUV energy deposition per particle becomes practically zero, seems to be a reasonable choice for a more realistic model that is independent of the particular parameter values at the inner boundary, or in other words, independent of the size of the simulation box. Indeed, we found (Section 5.2) that at pressures starting from $\tilde{P}_0 = p_{{\rm sat}} \sim 1\,\mu {\rm bar}$ and higher, the obtained solution shows quite weak dependence on $\tilde{P}_0$.

Another important modeling parameter related to the inner boundary pressure that has to be properly defined is the base (or the inner boundary) temperature $\tilde{T}_o$. General energy budget considerations constrain it to be at around 1000 K. The particular values of 750 K (Yelle 2004), 1200 K (García Muñoz 2007), and 1350 K (Koskinen et al. 2013) were used in the models of a hot Jupiter. In the model presented in this paper the behavior of a hot Jupiter environment insignificantly depends on a particular value of $\tilde{T}_o$. That is due to taking the inner boundary of the simulation box at sufficiently deep layers with high enough pressure that are unaffected by the stellar XUV. At such a deeply located inner boundary, an extremely steep (exponential) density growth takes place such that the condition of a constant mass flux supposes very small values for the velocity. Therefore, taking the velocity to equal zero at the inner boundary does not affect the numerical solution at all if the boundary is sufficiently deep.

Inclusion in the model of deep thermosphere layers with pressures higher than ∼1 μbar requires a certain attention to the infrared cooling effect due to ${\rm H}_{3}^ +$ and an understanding of how it affects the simulation results. To check this, we performed calculations with a proxy model for the infrared cooling with volume rates of about 10⁻⁷ erg cm³ s⁻¹ as reported in Yelle (2004) and García Muñoz (2007). The obtained results enable us to conclude that the temperature profile in the inner atmosphere regions (i.e., close to the planet surface) and the corresponding density distribution are somewhat sensitive to the additional energy-loss mechanism related to the infrared cooling (besides the adiabatic and Lyα cooling). However, because of its local character and the action first of all in the deep layers of the inner thermosphere, the ${\rm H}_{3}^ +$ cooling leads to the reduction of the simulated material escape rate approximately only by a factor of two. Such an insignificant difference gives a reason for an essential simplification of the model of atmosphere expansion and mass loss by ignoring the complex hydrogen chemistry at the base of the cold thermosphere or the use of its simplified modeling approximations.

An appropriate definition of conditions at the outer boundary also plays an important role in the modeling. Similar to the case of the inner boundary conditions discussed above, it would be reasonable to have the outer boundary conditions also defined such that the obtained numerical solutions are quantitatively independent of them in a sufficiently wide range of parameters.

Besides, it is necessary to distinguish between the two different types of outer boundary conditions that are related to different external physical circumstances. One is an open (or free) outer boundary, when the expanding atmosphere material can leave the vicinity of a planet, i.e., when the escape of mass and energy takes place and the escaping plasma can flow through the outer boundary unaffected. Another situation corresponds to the closed (or shut-up) conditions when the material flow is suppressed or restricted at a certain distance from a planet by some external physical factors not included self-consistently in the model. As an example of such restricting factors, an encountering stellar wind or planetary magnetic field can be referred. The choice of a particular type of outer boundary condition finally depends on the goal of the specific modeling and the region of application of the modeling results. In particular, in the case of a hot Jupiter, modeling of different parts of the magnetized planetary extended thermosphere may require different sets of external boundary conditions.

For example, simulation of a so-called wind zone of a hot Jupiter magnetosphere (Adams 2011; Trammell et al. 2011; Khodachenko et al. 2012) at sufficiently high latitudes where plasma leaves the planet along the open magnetic field lines needs the open outer boundary conditions, whereas modeling of the "dead zone" occupied by a stagnated plasma of an inner magnetosphere locked in the region of closed magnetic field lines would require the closed boundary.

The interaction of the magnetized hot Jupiter with a stellar wind results in the formation of a magnetosphere obstacle located at several planetary radii from a planet. Depending on the character of the encountering stellar wind plasma flow (subsonic or supersonic/Alfvénic), this magnetosphere obstacle can have the form of a paraboloidal shock or Alfvén wing system (Ip et al. 2004; Erkaev et al. 2005; Khodachenko et al. 2012) elongated in the direction of the plasma flow, which for the close-orbit planets also has to take into account the Keplerian velocity of orbital motion of the planet (Grießmeier et al. 2007; Cohen et al. 2011; Khodachenko et al. 2012). Because of these factors, the outer boundary conditions in the first numerical modeling of a hydrodynamic atmosphere escape on a hot Jupiter (Yelle 2004) were imposed at a relatively close distance of only 3 R_p from the planet. Note that a question regarding the interaction of an upcoming stellar wind with the expanding planetary plasma wind and the role of a planetary magnetic field in this process was not considered in Yelle (2004) or in the later modeling cases (García Muñoz 2007; Penz et al. 2008; Guo 2011, 2013; Adams 2011; Trammell et al. 2011; Koskinen et al. 2013; Khodachenko et al. 2012; Lammer et al. 2013; Erkaev et al. 2013). Thus it still remains a subject for future investigation.

In this section we describe the basic equations of the presented numerical model. For simplicity's sake we restrict the modeling by assuming a purely hydrogen atmosphere of a hot Jupiter composed in the general case of hydrogen atoms, ions, and electrons. As has been shown in García Muñoz (2007), Trammell et al. (2011), Koskinen et al. (2010), and Guo (2011), the atomic and Coulomb collisions in the modeled regions of a hot Jupiter upper atmosphere are efficient enough to ensure the "nonslippage" approximation for neutral hydrogen atoms, ions, and electrons such that the assumptions of V_a = V_i = V_e = V as well as temperature equilibrium T_a = T_i = T_e = T hold true. Here V_k and T_k with k = a, i, e correspond to the velocity and temperature of components of the partially ionized atmospheric gas, respectively, whereas V and T stand for its velocity and temperature as a whole.

For the typical temperatures ⩽10⁴ K realized in the upper atmosphere of a hot Jupiter (Koskinen et al. 2010), the proton–proton Coulomb collision cross section is well above 10⁻¹³ cm⁻². The proton–hydrogen collisions are dominated by resonant charge exchange with the cross section at the typical velocity of 10 km/ c⁻¹, at approximately 10⁻¹⁴ cm⁻². The location of the exobase for a hot Jupiter that corresponds to the region where the Knudsen number $Kn = \bar l/\Delta$ reaches 1 is determined by the condition when the mean free path $\bar l$ is equal to a typical scale of flow Δ. In the regions far from the planet, Δ may be taken to be on the order of the planetary radius, i.e., Δ ∼ R_p, and the condition Kn ∼ 1 is realized approximately in the layers with density ∼10⁴ cm⁻³. That remains well above the height range (10 R_p) considered in the modeling simulations. Close to the planet surface, the barometric height H, which is smaller than R_p, has to be taken as the typical scale Δ. However, because of much larger densities there, the mean free path $\bar l$ also becomes smaller, and the condition Kn < 1 holds true. Therefore, the fluid approach adopted in the present study is fairly valid.

Within the frame of these assumptions, taking into account the processes of hydrogen ionization (by XUV and electron collision) and recombination, the continuity, momentum, and energy balance equations of the model can be written as follows.

$$\begin{equation} \frac{\partial }{{\partial t}}n + \nabla ({{\bf V}n}) = 0 \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial }{{\partial t}}n_i + \nabla ({{\bf V}n_i }) = n_a \langle {\sigma _{{\rm XUV}} F_{{\rm XUV}} } \rangle - n_i \Lambda _e \sigma _{{\rm rec}} + n_a \Lambda _e \sigma _{{\rm ion}} \end{equation} \tag{ 2 }$$

$$\begin{equation} m\frac{\partial }{{\partial t}}{\bf V} + m({{\bf V}\nabla }){\bf V} = - \frac{1}{n}\nabla nkT - \frac{1}{n}\nabla n_e kT + m\nabla \frac{{{\rm GM}_p }}{r} \end{equation} \tag{ 3 }$$

$$\begin{eqnarray} &&\frac{3}{2}({n + n_e })k \cdot \left({\frac{\partial }{{\partial t}}T + ({{\bf V}\nabla })T + ({\gamma - 1})T\nabla {\bf V}} \right)\nonumber\\ &&\quad = n_a \langle {({E_{{\rm XUV}} - E_{{\rm ion}} })\sigma _{{\rm XUV}} F_{{\rm XUV}} }\rangle \nonumber\\ &&\qquad - E_{21} \Lambda _e \cdot ({n_a \sigma _{12} - n_{a,n = 2} \cdot \sigma _{21} }) - n_a E_{{\rm ion}} \Lambda _e \sigma _{{\rm ion}} \nonumber\\ && \end{eqnarray} \tag{ 4 }$$

where the following notations are used:

$$\begin{eqnarray} &&n = n_a + n_i,\quad n_e = n_i = x_{{\rm ion}} n\nonumber\\ &&\Lambda _e = n_e V_{Te} \end{eqnarray} \tag{ 5 }$$

In Equations (1)–(4) we use the standard cross sections of major processes, such as ionization by XUV: σ_XUV = 6.3 · 10⁻¹⁸ cm²; ionization by electron impact: $\sigma _{{\rm ion}} = 4\, \cdot\, 10^{ - 16}\, \cdot\, e^{ - {{E_{{\rm ion}} } /T}}\, \cdot\, T^{ - {1}}\, {\rm cm}^{2}$; recombination with an electron: σ_rec = 6.7 · 10⁻²¹ T^−3/2 cm²; hydrogen excitation and deexcitation: $\sigma _{12} = \sigma _{21}\, \cdot\, e^{ - {{E_{21} } /T}}\, {\rm cm}^{2}$ and σ₂₁ = 7 · 10⁻¹⁶ T⁻¹ cm², respectively; Lyα absorption: σ_Lα ∼ (10⁻¹⁵ ⋅⋅⋅ 10⁻¹⁴) T^−1/2 cm² with the temperatures scaled in units of characteristic temperature of the problem ($\tilde{T}_c = 10^4\, {\rm K}$), except for the exponent power indices in the expressions for σ_ion and σ₁₂, where temperature is taken in energy units. In addition, for the wavelength-averaged terms corresponding to XUV ionization and heating rates in Equations (2) and (4) we use the proxy model proposed in Trammell et al. (2011):

$$\begin{equation} \langle {\sigma _{{\rm XUV}} F_{{\rm XUV}} }\rangle = \frac{{t_{{\rm XUV}}^{ - 1} }}{{1 + ({\sigma _{{\rm XUV}} NL})^{1.5} }} \end{equation} \tag{ 6 }$$

$$\begin{equation} \langle {({E_{{\rm XUV}} - E_{{\rm ion}} })\sigma _{{\rm XUV}} F_{{\rm XUV}} }\rangle = \frac{{\varepsilon _{{\rm XUV}} t_{{\rm XUV}}^{ - 1} }}{{({1 + 0.2 \cdot \sigma _{{\rm XUV}}\, NL})^{1.2} }}, \end{equation} \tag{ 7 }$$

where $NL = \int\nolimits_r^\infty {n_a dr}$ is the neutral hydrogen column density.

The typical values for the mean ionization time t_XUV ≈ 6 hr and mean photoelectron energy (or energy per absorption event) ε_XUV ≈ 2.7 eV, which goes to heating, are taken in Equations (6) and (7) to correspond to the integrated stellar XUV energy flux of F_XUV ≈ 2 × 10³ erg cm⁻² s⁻¹ expected at 0.05 AU orbit around a Sun-analog star. Note that the decrease of heating (Equation (7)) is close to an inverse rather than an exponent power of the density column. The reasons to use the proxy model instead of one calculated with the frequency-binned real Solar spectra are as follows. The deviation of the proxy model from the spectrum-based one is sufficiently small to ensure that the error in the computed mass loss rate is less than 25%. Besides that, using the analytical model enables one to make simple analytical estimations of the expected solutions, which help to verify the correctness of the obtained numerical results (Section 5.2) and to understand the physics of the involved processes.

To take into account the fact that the star illuminates only one side of the planet, in our calculations we reduce the radiation flux (and the corresponding values given by Equations (6) and (7)) by a factor of two, i.e., take η_net = 50%. In addition, to be definitive we suppose the heating efficiency η_h of stellar XUV (Chassefière 1996; Yelle 2004; Lammer et al. 2009; Koskinen et al. 2013) to be 50% and therefore divide by two once more the outcome of Equation (7), which is then used as an approximation for the term 〈(E_XUV − E_ion)σ_XUVF_XUV〉 in the energy Equation (4).

Besides the XUV heating term and adiabatic cooling due to plasma expansion (the last term at the left-hand side of Equation (4)), the energy Equation (4) also takes into account energy loss and gain due to excitation and deexcitation of the n = 2 level of hydrogen atom by electron impact as well as cooling due to ionization. We assume that the concentration of excited atoms remains sufficiently low, i.e., n_{a, n = 2} ≪ n_{a, n = 1} ≈ n_a. The validity of this assumption is specially verified in the course of the modeling simulation. Because of a too-low rate at the considered densities, we also neglect in Equation (2) the process of three-body recombination, which is the reverse of the inelastic impact ionization.

To complete the energy Equation (4) one has to define a way to calculate the concentration n_{a, n = 2} of the excited atoms at the level n = 2. In our model it is controlled by the density of the Lyα photons that are emitted by the excited atoms and reabsorbed by the atoms in a nonexcited state (i.e., level n = 1). In the case of an optically thin gas envelope, the concentration of excited atoms remains very small (negligible) because of fast radiative deexcitation. On the other hand, in the optically thick case, the concentration n_{a, n = 2} may remain moderately high because of Lyα photon sequential reabsorptions until the photons leave the system. Such a process of multiple reabsorptions and reemissions of a photon with a chaotic change of its travel direction (after each reemission act) can be described in the first approximation as diffusion with a coefficient determined by the average distance (n_aσ_Lα)⁻¹ and velocity c between reabsorptions. The corresponding diffusion equation for the Lyα photons density N_La and concentration of excited atoms may be written as follows.

$$\begin{eqnarray} \frac{d}{{dt}}n_{a,n = 2} &=& \Lambda _e \cdot ({n_{a,n = 1} \sigma _{12} - n_{a,n = 2} \sigma _{21} })\nonumber\\ &&- \frac{{n_{a,n = 1} }}{{\tau _{21} }} + n_{a,n = 1} \cdot cN_{La} \sigma _{La} \end{eqnarray} \tag{ 8 }$$

$$\begin{equation} \frac{\partial }{{\partial t}}N_{La} = \frac{{n_{a,n = 2} }}{{\tau _{21} }} - n_{a,n = 1} \cdot cN_{La} \sigma _{La} + {div}\frac{c}{{n_a \sigma _{La} }}\nabla N_{La}. \end{equation} \tag{ 9 }$$

The terms in the right-hand side of Equation (9) stand for the photon emission by the excited atoms with mean emission time τ₂₁, absorption by nonexcited atoms, and spatial diffusion, respectively. For the typical scale of a system, which approximately equals R_p ∼ 10¹⁰ cm and σ_La ∼ 10⁻¹⁵ cm², the gas becomes optically thick for Lyα photons at densities above 10⁵ cm⁻³.

Introducing a normalized number density of Lyα photons $\tilde{N}_{L\alpha } = \tau _{21} c \cdot N_{L\alpha } \sigma _{L\alpha }$, the stationary concentration of the excited atoms may be expressed from Equation (8) as

$$\begin{equation} n_{a,n = 2} = n_{a,n = 1} \frac{{\tilde{N}_{L\alpha } + \Lambda _e \cdot \sigma _{12} \tau _{21} }}{{1 + \Lambda _e \cdot \sigma _{21} \tau _{21} }}. \end{equation} \tag{ 10 }$$

By substitution of this solution in Equations (4) and (9) and taking into account that n_{a, n = 1} ≈ n_a, we obtain a modified energy-balance equation and the equation for normalized number density of Lyα photons:

$$\begin{eqnarray} &&\frac{3}{2}({n + n_e })k \cdot \left(\frac{\partial }{{\partial t}}T + ({{\bf V}\nabla })T + ({\gamma - 1})T\nabla {\bf V} \right)\nonumber\\ &&\quad = n_a \langle ({E_{{\rm XUV}} - E_{{\rm ion}} })\sigma _{{\rm XUV}} F_{{\rm XUV}} \rangle - n_a E_{21} \Lambda _e \sigma _{21}\nonumber\\ &&\qquad \times ({e^{{{ - E_{21} } /T}} - \tilde{N}_\alpha }) - n_a E_{{\rm ion}} \Lambda _e {\rm }\sigma _{{\rm ion}} \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \frac{\partial }{{\partial t}}\tilde{N}_{L\alpha } &=& \tau _{21} c \cdot n_a \sigma _{L\alpha } \cdot \Lambda _e \cdot \sigma _{21} ({e^{ - {{E_{21} } /T}} - \tilde{N}_{L\alpha } })\nonumber\\ &&+ \nabla \cdot \left({\frac{c}{{n_a \sigma _{L\alpha } }}\nabla \tilde{N}_{L\alpha } } \right) = 0. \end{eqnarray} \tag{ 12 }$$

Because the emission time τ₂₁ ∼ 1 ns is very small, we neglect the terms proportional to Λ_e · σ₁₂τ₂₁ and Λ_e · σ₂₁τ₂₁ in the final form of the Equations (11) and (12), as compared with one.

When the system is essentially optically thick, the diffusion term in Equation (12) tends to zero, and the number density of Lyα photons saturates at a level where collisional excitation and deexcitation balance each other, resulting in $\tilde{N}_{L\alpha } = e^{{{ - E_{21} } /T}}$. Note that for temperatures T ≪ E₂₁ the population of excited atoms is much smaller than that of nonexcited, ground-state atoms: ${{n_{a,n = 2} } /{n_a }} \approx \tilde{N}_{L\alpha } \ll 1$. In that respect, one of the reasons for inclusion in Equation (11) of a relatively unimportant electron impact ionization term is to ensure a nonzero cooling in the case of total saturation of the n = 2 level by Lyα photons.

Further on, we represent the temperature and velocity of the simulated upper atmospheric material in units that correspond to the characteristic temperature of the problem taken as $\tilde{T}_c = 10^4\, {\rm K}$ and the corresponding $\tilde{V}_c = V_{Ti} = 9.8\, {\rm km/s}$, whereas the distances will be normalized to a typical scale of the system, defined by planetary radius R_p.

Below we consider numerical solutions of a spherical 1D HD model that describe different behavior regimes of a hot Jupiter's upper atmosphere, corresponding to the two basic types (open and closed) of outer boundary conditions.

The position of the outer boundary is taken such that the results of the calculations, such as the mass loss rate, do not depend on the distance of its particular location. The question discussed in the literature is the position of the exobase relative to the model outer boundary. For HD modeling of some particular hot Jupiters (e.g., HD 209458b), it was shown in several papers that the exobase is farther than 10 R_p. Thus in our modeling we also take the location of the outer boundary at 10 R_p, where the numerical solution practically does not depend on the boundary's exact position and the use of hydrodynamic modeling approach is justified. More detailed (García Muñoz 2007) or partially kinetic conditions (Koskinen et al. 2013) also lead to similar results and conclusions.

When the cooling mechanism due to Lyα radiation is taken into consideration, the model admits two general solutions related to the global energy balance. One of the solutions represents the case already mentioned above with an atmospheric material free outflow when the XUV heating is balanced by the advection and adiabatic cooling. Another solution corresponds to a restricted-flow regime, in which the XUV heating is balanced mostly by Lyα cooling. The particular values of physical quantities are also different for these two types of solutions. The restricted-flow solution is characterized by a hotter and about an order of magnitude denser stagnated plasma envelope. An important issue regarding the free-outflow and restricted-flow solutions is that their realization in our model is completely controlled only by the specific outer boundary conditions. Note that the free-outflow and restricted-flow regimes are physically relevant (and may be attributed) to the "wind" and "dead" zones of a hot Jupiter magnetosphere, respectively (Trammell et al. 2011; Khodachenko et al. 2012).

As mentioned above (Section 3.2), the Lyα cooling is inefficient at temperatures below several thousand Kelvin. Therefore, a small outflow of material is needed to compensate for the excess XUV heating in such a dense deep thermosphere. In that respect, instead of a totally closed (shut-up) condition at the outer boundary, we allow a small flow of material at a speed on the order of 0.1 km s⁻¹. This in fact constitutes the matter of the restricted-flow regime. Such an assumption is relevant for the conditions in a "dead zone" of a hot Jupiter's inner magnetosphere or in the front region of interaction with a stellar wind because a small portion of the planetary plasma might still leave the planet there by convection to higher latitudes or to the tail.

To deal with the open outer boundary in practice, we assume the spatial derivatives of physical values to be small, i.e., asymptotically approaching zero with a decreasing mesh size and increasing distance.

A similar approach to the study of atmosphere expansion regimes was adopted for the first time in García Muñoz (2007), where instead of velocity a fixed pressure was imposed at the outer boundary, and solutions with an escaping material flow varying from a supersonic up to a slow subsonic one were obtained.

Figure 1 shows the result of a modeling calculation with the typical model parameters summarized in Table 1. The case of an open outer boundary condition was considered. In line with the discussion in Section 3.2, the Lyα reabsorption effects are not important for the studied regime of material free outflow and therefore were not included in the model.

Zoom In Zoom Out Reset image size

Table 1. Model Parameters

Mp	Rp	FXUV (0.05 AU orbit, Sun-analog star)	$\tilde{T}_0$	$\tilde{P}_0$	Outer Boundary Condition
0.7 MJ	1010 cm	2 × 103 erg cm−2 s−1	1000 K	10−5 bar	Open

Download table as:  ASCIITypeset image

As can be seen in Figure 1, the maximum temperature realized in the free-outflow regime is about 8500 K and is achieved at a height of about 1.75 · R_p. Then, because of adiabatic cooling, the temperature quickly decreases with height because the XUV heating is not enough to compensate for it once the ionization of the expanding gas takes place. The ionization degree x_ion = n_i/n increases sharply near the temperature maximum, and the velocity of the expanding atmospheric plasma wind driven by the increase in the pressure gradient also increases, mostly at the ionization front. Because the temperature and density quickly decrease with distance, the escaping material reaches supersonic velocity at a distance of about 6 · R_p.

A comparison of the obtained results with the outcomes of the closest (in the sense of model parameters and simulation conditions) model by Koskinen et al. (2013) reveals that, in spite of general (phenomenological) agreement, the models show some nonprincipal difference in specific physical quantities. These may be attributed to the slightly different values of XUV flux used in the modeling simulations, which result in different heating rates q_XUV = n_a〈(E_XUV − E_ion)σ_XUVF_XUV〉 in Equations (4) (and (11) realized in the course of modeling, and to different net and heating efficiencies. In particular, our model gives a slightly higher density and a bit smaller velocity with a slightly smaller (than by Koskinen et al. 2013) maximum XUV heating rate: $q_{{\rm XUV}}^{\max } \approx 4 \cdot 10^{ - 8}\, {\rm erg} \cdot {\rm cm}^{ - 3} \cdot {\rm s}^{ - 1}$.

The calculation of mass loss rate $\dot M = 4\pi mnr^2 V_r$ yields the value 7 · 10¹⁰ g s⁻¹, which is by a factor of almost two larger than that predicted with the model by Koskinen et al. (2013). However, the mass loss rates become practically equal if we take the same net efficiency, $\eta _{^{{\rm net}} } = 25\%$, as that used by Koskinen et al. (2013).

The most significant difference between our model results and those of Koskinen et al. (2013) concerns the predicted maximum temperature, which in our case is 3500 K less. One of the reasons is that Koskinen et al. (2013) use a Lyα cooling reduced by a factor of 0.1. Switching off the Lyα cooling in our model reduces the difference down to 2500 K.

Regardless of the maximum temperature, the value of pressure-averaged temperature

$$\begin{equation} \bar T_p = \frac{{\int\nolimits_{p_1 }^{p_2 } {T(p)\,{\rm d(ln }p{\rm)}} }}{{{\rm ln(}p_2 /p_1 {\rm)}}} \end{equation} \tag{ 14 }$$

over the distance interval from R_p to 3R_p, which is used in Koskinen et al. (2013) as a measure of a possible transit effect of a planet, in the case of our model is 6000 K. That is rather close to the temperature range of 6200–7800 K reported in Koskinen et al. (2013).

The behavior of the temperature and density of the expanding planetary atmosphere as parametric functions of pressure in the deep region of dense thermosphere layers near the inner boundary, obtained by solving numerically the model equations, is shown in Figure 2. The position of the inner boundary itself in this study was taken slightly (a few percent) below the level of r = 1 (in units of R_p). Note that temperature smoothly approaches the base (inner boundary) temperature T_o = 0.1 (in units of $\tilde{T}_c = 10^4\, \rm K$), and velocity decreases down to negligibly small values.

Zoom In Zoom Out Reset image size

To understand the nature of the obtained numerical result, let us compare it with a simple analytical solution. Close to the inner boundary, in the dense layers of the thermosphere where the velocity is rather small, we can ignore inertial terms and ionization in the momentum Equation (3). Then, the barometric density profile n(r) = n₀ exp(− r/H(T, r)) with the local scale height of H(T, r) = r²kT/(GmM_p) is an approximate solution of Equation (3), which under the above-mentioned conditions has the following form:

$$\begin{equation} \frac{k}{n}\nabla nT \approx - \nabla \frac{{GmM_p }}{r}. \end{equation} \tag{ 15 }$$

We assume that XUV absorption takes place in a thin layer near r ≈ R_p with very large column densities NL so that σ_XUVNL ≈ σ_XUVnH ≫ 1, and it is balanced by adiabatic expansion.

Then, from the energy balance Equation (4) (or (11)) in a quasi-stationary state, using the approximation Equation (7) taken at NL ≫ 1, assuming as a zero-order approximation the barometric density distribution with the scale height H(T, r), and expressing r²V_r in terms of mass loss rate $4\pi mnr^2 V_r = \dot M = {\rm const}$, we obtain an approximate relation between temperature and density in the vicinity of the inner boundary:

$$\begin{equation} \left({\frac{T}{{T_o }}} \right)^6 \left({\frac{{R_p }}{r}} \right)^8 \frac{n}{{n_o }} = ({\sigma _{{\rm XUV}} H_o n_o })^{ - 1} \left({\frac{{4\pi R_p^2 mEF}}{{\sigma _{{\rm XUV}} k\tilde{T}_o \dot M}}} \right)^5\! =\! {\rm const,} \end{equation} \tag{ 16 }$$

where EF = ε_XUV/(0.2^1.2 · t_XUV), T_o is the inner boundary temperature measured in units of $\tilde{T}_c = 10^4\, \rm K$, and $H_o = H(\tilde{T}_0,R_p)$ is the local scale height in the vicinity of the inner boundary. In these calculations we continue using t_XUV ≈ 6 hr and ε_XUV ≈ 2.7eV, which correspond to the stellar XUV flux expected at 0.05 AU orbit around a Sun-analog star.

Solution of the barometric Equation (15), with Equation (16) taken into account, yields

$$\begin{equation} \frac{T}{{T_o }} = ({{r /{r_{\min } }}})^{{8 /5}} + \frac{1}{{13\, H_o r}}[ {({{r /{r_{\min } }}})^{{{13} /5}} - 1}]. \end{equation} \tag{ 17 }$$

The starting point r_min of the solution is determined from Equation (16) by taking T = T_o and n = n_o. The maximum pressure (or density, i.e., $p_{\max } /k\tilde{T}_0 ^{}$) up to which solution Equation (17) is valid follows from the right-hand side of Equation (16) with const = 1:

$$\begin{equation} p_{\max } = n_0 k\tilde{T}_0 \approx \frac{{k\tilde{T}_o }}{{\sigma _{{\rm XUV}} H_o }}\left({\frac{{4\pi R_p^2 {mEF}}}{{\sigma _{{\rm XUV}} k\tilde{T}_o \dot M}}} \right)^5 \end{equation} \tag{ 18 }$$

The term in brackets in Equation (18) is approximately a ratio of energy delivered by XUV to energy removed by the expanding material flow, if each particle has energy $k\tilde{T}_o$. Because the material escape is almost energy limited and the particle energy is 10 times higher than $k\tilde{T}_o$, the term in brackets is about 20–30. This means that the maximum pressure is much larger than that at the layer with a density of 1/σ_XUVH_o where a significant portion of XUV is absorbed. Typical hot Jupiter conditions in this case yield an estimate for pressure of p_max ∼ 10 μbar.

The analytic solution Equation (17) expressed by means of Equation (16) as a function of pressure is shown in Figure 2. The value p_max = 10μbar corresponds well to the point (marked by an ellipse) where the numerical solution, approaching the inner boundary level, begins to deviate significantly from the analytical solution of Equation (16). Note that in this region the velocity is already very small. At larger pressures, corresponding to deeper atmospheric layers, the considered model solution is governed more by numerical factors, such as the small viscosity (η = 10⁻⁵ in dimensionless units) introduced into Equation (3) for stability reasons. At such very large pressures and so deep in the thermosphere, the hydrogen chemistry, molecular and eddy diffusion are probably more important than the adiabatic cooling considered in this example.

To investigate how the inner boundary location may influence the obtained numerical results regarding the escape of hot Jupiter atmospheric material and to verify the possibility of neglecting the effects of infrared cooling by ${\rm H}_{\rm 3}^ +$ at the base of the thermosphere (as presently assumed), we performed a series of dedicated calculations and present their outcome below.

Figure 3 shows a dependence of the simulated mass loss rate on a value of fixed pressure $\tilde{P}_o$ at the inner boundary, or rather on the position of the inner boundary relative to the fixed distribution Equation (13). One can see that for$\tilde{P}_o$ less than p_sat ∼ 1 μbar the escaping flow is strongly affected by a particular value of the inner boundary pressure$\tilde{P}_o$. Such a dependence of the numerical solution on the inner boundary condition (i.e., on the location of the inner boundary of the simulation box) has been referred to in several previous studies (e.g., Trammell et al. 2011; Adams 2011), where proportionality of the mass loss rate to the boundary pressure was assumed as an input condition. One can see from Figure 3 that, for $\tilde{P}_o$ higher than p_sat ∼ 1μbar, the influence of the inner boundary pressure on the mass loss rate becomes less significant.

Zoom In Zoom Out Reset image size

As to the behavior of other physical quantities in the model, we find that between the cases of $\tilde{P}_0 = 100{\rm }\, \mu {\rm bar}$ and $\tilde{P}_0 = 0.1\, \mu {\rm bar}$ the material density at distances of a few R_p varies by a factor of two, whereas the location of the temperature maximum for $\tilde{P}_0 = 100\,\mu {\rm bar}$ shifts by ∼0.3 R_p to smaller heights.

The choice of a particular maximum pressure at the inner boundary is also related to the depth of absorption of XUV photons. The absorption parameter σ(ε) · NL for a photon of given energy ε = hν can be estimated as

$$\begin{eqnarray*} \sigma (\varepsilon) \cdot NL &\approx &\sigma _{{\rm XUV}} ({{{E_{{\rm ion}} } /\varepsilon }})^3 n_o H_0 \nonumber\\ &\approx & \sigma _{{\rm XUV}} ({{{E_{{\rm ion}} } /\varepsilon }})^3 {{R_p^2 } /{({GmM_p })}} \cdot \tilde{P}_o\nonumber\\ &\approx & 0.5 \cdot 10^{10} ({{{E_{{\rm ion}} } /\varepsilon }})^3 \tilde{P}_o. \end{eqnarray*}$$

Therefore, fixing the thermosphere inner boundary (i.e., depth) at the level of $\tilde{P}_o = 1\, \mu {\rm bar}$ means that photons with energy higher than ε > 20 · E_ion, which correspond approximately to wavelengths shorter than 5 nm, will be excluded from the consideration. Such a restriction would limit the generality of the model.

The discovered mass loss "saturation" effect for sufficiently large inner boundary pressures ($\tilde{P}_o > 1{\rm }\, \mu {\rm bar}$) is related to a rather weak (for large $\tilde{P}_o$) dependence of the obtained numerical solution on the base temperature $\tilde{T}_o$. In particular, the variation of $\tilde{T}_o$ from 750 K to 1500 K results in a change in the mass loss rate $\dot M$ by only 14%.

Figure 4 shows the dependence of the mass loss rate on the intensity of XUV radiation of a Sun-like star illuminating the planet (full analog of HD 209458b), expressed in terms of the orbital distance to the star. It has a power-law character. The fact that the increase of the mass loss rate with decreasing distance is very close to R⁻² enables one to draw a conclusion on the proportionality between the mass loss rate and the XUV flux. Similar behavior has been reported in a number of papers (e.g., Lammer et al. 2003; Garcia Muñoz 2007).

Zoom In Zoom Out Reset image size

For a Jupiter-type planet at an orbital distance below 0.02 AU near a Sun-like star, the tidal forces and Roche lobe effects, which are not included in our model, become important (Garcia Muñoz 2007; Penz et al. 2008; Erkaev et al. 2007). Therefore, we restrict the considered range of orbital distances by 0.02 AU. To justify the omission of a tidal force, a test calculation was performed for a tidally locked system with the gravity potential in the equatorial plane of mGM_p/r(1 + 1/2(r/r_L)³), where r_L is a distance to the Lagrange point, which for the considered hot Jupiter (HD 209458b) system is r_L ≈ 4.6. It was found that the tidal force increases the escape velocity (on the day side of the planet). However, the mass loss rate increases only by 10% and the maximum temperature decreases by only 10%. Therefore, it may be concluded that within the above assumptions the main results of our study are not strongly affected by the inclusion of a tidal force, though it should be taken into account for a complete model of the HD 209458b system (after all of the major effects are properly studied and understood). At the same time, because the main goal of this paper is to build a comprehensive (1D) self-consistent model of a hot Jupiter's upper atmosphere escape for its further generalization to the case of a magnetized planet, without the inclusion of tidal force the role of other material escape drivers and their mutual relation becomes more clear and easy to understand.

To study the role of the infrared cooling effect in our computations, we followed the approach of Penz et al. (2008) and included the term −n_aC_ir in the right-hand side of Equation (4). The taken coefficient C_ir = 0.05 K/s approximately corresponds to the ${\rm H}_{\rm 3}^ +$ cooling rate reported in Yelle (2004) and García Muñoz (2007).

The ${\rm H}_{\rm 3}^ +$ is produced in reactions that involve ionized molecular hydrogen ${\rm H}_{\rm 2}^ +$ in the region where H₂ is the dominant species instead of atomic hydrogen. To take this fact into account and restrict the infrared cooling by layers populated by molecular hydrogen, we included reactions of molecular hydrogen dissociation and association, H₂ + M↔2H + M, with the rates provided in Yelle (2004) and García Muñoz (2007). Assuming a local balance between H₂ and H as determined by the rates of association and dissociation, the relation between species is given by

$$\begin{equation} \displaystyle n_{\rm H_2} = \frac{{8 \cdot 10^{ - 33} ({{{300} /T}})^{0.6} }}{{1.5 \cdot 10^{ - 9} e^{ - {{48000} /T}} }} \cdot n_H \cdot n_H. \end{equation} \tag{ 19 }$$

In the presence of H₂ the energy Equation (4) should be modified. At short wavelengths the cross section of XUV absorption by H₂ is about 2.5 times larger than that for H. Without making a big error we can put it to be twice larger, that is, the same cross section as for H but for two hydrogen atoms comprising the hydrogen molecule. Finally, the heating term in Equation (4) may be modified corresponding to $n_H \sigma _{{\rm XUV},H} + n_{\rm H_2} \sigma _{{\rm XUV},\rm H_2} \approx ({n_{H\Sigma } + n_{\rm H_2} }) \cdot \sigma _{{\rm XUV},H}$. Therefore, in the infrared cooling term −n_aC_ir, instead of n_a we use the sum density, $n_{H\Sigma } = n_{\rm H} + n_{{\rm H_2}}$. After combining the momentum equations for H₂ and H, the momentum equation can also be expressed in terms of the sum density n_HΣ and average mass $\bar m = m({{{\rm n}_{\rm H} + 2{\rm n}_{{\rm H_2}} }})/({{{\rm n}_{\rm H} + {\rm n}_{{\rm H_2}} }})$:

$$\begin{eqnarray} \frac{\partial }{{\partial t}}{\bf V} + ({{\bf V}\nabla }){\bf V} &=& {-} \frac{1}{{\bar m \cdot n_{H\Sigma } }}\nabla n_{H\Sigma } kT\nonumber \\ && - \frac{1}{{\bar m \cdot n_{H\Sigma } }}\nabla n_e kT + \nabla \frac{{{\rm GM}_p }}{r}. \end{eqnarray} \tag{ 20 }$$

This equation replaces Equation (3) in the considered set of model equations. The average mass $\bar m$ can be expressed through sum density n_HΣ using Equation (19).

The corresponding energy density rates for XUV absorption, advection, as well as adiabatic and radiation cooling obtained by solving the modified equations of the model are shown in Figure 5. With the inclusion of infrared cooling C_ir, the peak XUV absorption rate is about two times larger than in the calculation presented in Figure 1. This happens because of the deeper penetration of XUV that is due to the generally smaller column density.

Zoom In Zoom Out Reset image size

Note that radiative cooling consists of an infrared part at low temperatures and a Lyα part at sufficiently high temperatures. The maximum infrared cooling rate reaches the maximum XUV heating rate, and both have values slightly larger than 10⁻⁷ erg · cm^-3 · s⁻¹, which agrees well with the results of Yelle (2004). The value of pressure at this maximum is about 0.07 μbar. The transition from molecular to atomic hydrogen occurs at a slightly larger height, where a pressure of 0.015 μbar and temperature of about 2000 K are realized. One can see that close to the conventional planet surface, i.e., r = 1 (in units of R_p), infrared cooling clearly dominates. Then, after some short space interval (∼0.05R_p), the Lyα radiation becomes an efficient coolant. Finally, at heights above 1.5 R_p only the adiabatic and advection cooling counteract the XUV heating.

The spatial distribution profiles of the main physical quantities obtained in the course of the calculations including the effects of H₂ chemistry and ${\rm H}_{\rm 3}^ +$ cooling are shown in Figure 6. A comparison of these results with those presented in Figure 1 reveals that accounting for the infrared cooling effects in the deep layers does not lead to any principal difference in the system behavior. Accounting for the ${\rm H}_{\rm 3}^ +$ cooling effect as ∝n_HΣC_ir leads to a reduction of the material escape rate of about two times. Variation of the C_ir value by an order of magnitude results only in a moderate (about a factor of 1.7) change in the mass loss rate. Thus, the effect of ${\rm H}_{\rm 3}^ +$ cooling has been estimated more or less robustly. This result differs from that of Penz et al. (2008), where no influence of infrared cooling was reported. In that respect we note that, in spite of the peak cooling rate in our case and in Penz et al. (2008) being the same, the XUV heating rate in Penz et al. (2008) is about an order of magnitude larger.

Zoom In Zoom Out Reset image size

Here we would like to demonstrate the performance of the above-discussed (Section 3.1) effect of the "ionization boost" (due to the contribution of the electron and ion pressure gradients) in the increase of the atmospheric material outflow speed.

It follows from the integration of Equation (3) over r (from the planet surface to infinity) that the term with the electron pressure gives a positive input to the kinetic energy of the expanding flow. The electron density increases with the radial distance over the integration interval from n_e = 0 at r = R_p up to some maximum value and then again falls to zero, and the total number density of hydrogen n monotonically decreases with distance. Because of that, two areas can be distinguished depending on the radial distance r: (1) the area of increasing electron pressure (closer to the planet) where the pressure gradient acts as a decelerating factor on the plasma flow, and (2) the outer area of decreasing pressure where the electron pressure gradient accelerates the expanding plasma flow. Because the total density n decreases exponentially fast with distance, and the location of maximum of electron pressure is spatially close to that of the maximum temperature, the overall contribution of the electron pressure gradient term in the flow kinetic energy can be estimated as

$$\begin{equation} - \int\nolimits_{R_p }^\infty {\frac{1}{n}\nabla n_e kT} \approx kT_{\max }. \end{equation} \tag{ 21 }$$

The effect of acceleration of the escaping wind by the additional pressure gradient of the ionized gas is illustrated in Figure 7, where besides the typical free-outflow solution (the same as in Figure 1), the calculations without ionization and with zero electron pressure are shown. In the limiting case without ionization, the XUV radiation equally heats the whole volume of gas, resulting in a very large escape velocity and high temperature that continuously increase with distance. The speed of the expanding wind exceeds the adiabatic sonic point $\sqrt {{{\gamma kT} /m}}$ at about 7.3 R_p. In another limiting case, shown in Figure 7, the effect of material ionization is included. This results in the spatial dependence of XUV heating and its strong decrease at heights beyond several R_p. At the same time, the electron pressure in this limiting case (i.e., the ionization boost) is set to zero, so the effect of the additional pressure gradient of the ionized gas is artificially ignored. The resulting flow velocity in this case is essentially smaller than that in the true case, and it never reaches the sound speed values, while the overall mass loss is 1.5 times smaller. Contrary to these limiting cases, the true case, with the effects of ionization and the additional pressure gradient of the ionized gas (electron pressure) properly taken into account, is characterized by an increase in outflow velocity, which reaches a sonic point at the distance of 6.25 R_p. Note that because of the different spatial distributions of temperature in the cases considered here, the corresponding density profiles are also essentially different.

Zoom In Zoom Out Reset image size

The comparison of the limiting cases presented in Figure 7 enables one to conclude that the effect of the ionization boost is essentially important in the case of a hot Jupiter's expanding atmosphere. It appears to be a major mechanism that compensates to a certain degree for the decelerating action of adiabatic cooling and provides the necessary additional propulsion force to support the free-outflow regime of atmosphere expansion.

Inclusion of the tidal force effects in the comparative study of the "ionization boost" slightly but not principally affects the results of the modeling. In particular, whereas in the case without the tidal force the effect of the "ionization boost" at r = 8 Rp causes an increase in the escape velocity of ΔV ≈ 0.5 (see Figure 7), the inclusion of the tidal force under the same other conditions results in ΔV ≈ 0.4.

It should be noted that electron pressure acts primarily on ions. In the outer, highly ionized layers of the expanding atmosphere, where the electron pressure gradient acts as an accelerating force, the ions naturally appear as a major mass-contributing component. As to the inner regions populated by partially ionized plasma, it was shown by Guo (2011) that because of fast charge exchange the added momentum of ions is effectively transferred to atoms.

When the closed condition is applied at the outer boundary of the simulation box, the expanding atmospheric material flow gradually decays until a steady state with an almost zero flux through the outer boundary is achieved. This situation corresponds to the suppressed- or restricted-flow regime. Because of the absence of adiabatic cooling in this case, the energy-balance condition is provided for due to Lyα emission. This is true everywhere except in a cold, nonionized, and dense thermosphere with a temperature of several thousand K and less. In reality, in the case of a complete absence of material escape, a solution with a temperature at the inner boundary on the order of 1000 K, provided only by Lyα cooling, is not possible. Therefore, in our numerical model we allow a small outflow velocity at the outer boundary such that the corresponding small flow in the dense thermosphere can sustain an energy balance with the XUV heating. Under these circumstances, the material escape rate is proportional to the allowed fixed-outflow velocity at the outer boundary. It is worth, in that respect, comparing the different temperature profiles realized for different outflow velocities.

Figure 8 shows the temperature profiles obtained for different values of the outer boundary outflow velocity. With a series of subsequent simulation runs we found that for V_out < 0.1 km s^-1 the solution does not depend on the outflow velocity V_out at the outer boundary any more. The maximum value of the velocity achieved inside the simulation box is about two to three times larger than the outer boundary outflow velocity, but it still remains rather small, V_max < 0.25 km s^-1. At a fixed outer boundary velocity of V_out = 0.1 km s^-1, the mass loss rate is 3.5 times smaller than in the case of free outflow (i.e., open outer boundary).

Zoom In Zoom Out Reset image size

It can be seen in Figure 8 that the material temperature grows much faster near the planet surface as compared to the case of an open outer boundary condition. The temperature grows continuously over large distances because the adiabatic cooling does not work in this case, except in only a restricted region at r ⩽ 1.15 R_p close to the planet surface. The pressure-averaged temperature (see Equation (14)) in the interval of heights R_p ⩽ r ⩽ 3 R_p is about 9000 K. This leads to the formation around a planet of a gas envelope that is several times (up to an order of magnitude) denser and somewhat hotter, as compared to the case of a free-outflow regime.

As already pointed out, the temperature in the case of a restricted-flow regime depends on the efficiency of Lyα cooling. In view of the fact that deep near the planet surface the atmospheric gas is optically thick for Lyα, it is important to also check for the possible influence of Lyα reabsorption (included in Equation (11)). It becomes important when the number of Lyα photons is sufficiently large: $\tilde{N}_{L\alpha } e^{{{E_{21} } /T}} \sim \tau _{21} {\rm }\Lambda _e \cdot \sigma _{21} ({n_a \sigma _{L\alpha } R_p })^2 > 1$.

This estimation follows from Equation (12) after an approximation of spatial derivatives as $\nabla \nabla \sim R_p^{ - 2}$. For the conditions of restricted flow, the reabsorption affects the Lyα cooling rate at heights below 2 R_p (Figure 9). Figure 9 shows the result of a model calculation with a relatively high Lyα absorption cross section of σ_Lα ∼ 10⁻¹⁴T^-1/2 cm².

Zoom In Zoom Out Reset image size

As expected, because of strong reabsorption and an established balance between hydrogen excitation and deexcitation that diminishes the efficiency of cooling, the temperature visibly increases at distances below 2 R_p. The maximum increase in temperature as compared to the case without Lyα reabsorption is about 1500 K. Farther from the planet, the effect of Lyα reabsorption decreases, and for r > 2.5 R_p it becomes negligible. As to the influence of Lyα reabsorption on the gas density at distances up to several R_p, i.e., inside the dead-zone region that might be formed in the presence of a magnetic field, as can be seen in Figure 9, such an influence is rather small and practically undistinguishable in the numerical calculation result. Therefore, depending on the purpose of the particular numerical simulation, in some cases the effect of Lyα reabsorption can be neglected, despite its obvious influence on the temperature of gas close to the planet.

Figure 9 also allows one to compare the density distributions realized in the cases of open and closed outer-boundary conditions, i.e., in the free-outflow and restricted-flow regimes, which may correspond, for example, to the "wind" and "dead" zones of a hot Jupiter magnetosphere (Trammell et al. 2011; Khodachenko et al. 2012). In particular, the plasma density in the case of a restricted-flow regime is up to an order of magnitude higher than in the case of free outflow.

Based on the obtained numerical solutions, it is possible to estimate the strength of a hot Jupiter's magnetic field needed to affect and confine the hot gaseous envelope around the planet, energized by the stellar XUV. For this purpose we check the spatial behavior of the plasma beta 8πnT/B² in the equatorial plane for different values of magnetic field at r = R_p, assuming for the sake of being definitive an undisturbed magnetic dipole field geometry and a quasi-stationary restricted-flow regime. The result is presented in Figure 10, which in particular shows that surface fields larger than 1 G should strongly affect the plasma outflow. The area of strongest magnetic control is located near 2 R_p, whereas the overall size of the magnetic field–dominated region may reach, depending on the field strength, up to several R_p. Note that because of the steep density increase close to the planet surface, thermal pressure exceeds magnetic pressure even for very large surface fields.

Zoom In Zoom Out Reset image size

Another interesting issue concerns the flow above the dipole poles. As is well known, if plasma moves along the magnetic field lines in a strong dipole field, the divergence of material flux (guided by the magnetic field lines) at the pole regions is larger than that in the case of a spherically isotropic expansion. Because adiabatic cooling is crucial for sustaining the temperature of the expanding material, a certain reduction of the temperature and escape velocity in the polar regions, as compared to the isotropic expansion case, may be expected. This effect is demonstrated in Figure 11, which shows the outcome of our 1D HD simulation with an artificially increased ("cubic") divergence term as ∂V/∂r + 3V/r (analogous to the divergence of a dipolar magnetic field in polar regions), taken instead of the true divergence term divV = ∂V/∂r + 2V/r, which corresponds to an isotropic spherically symmetric expansion. Justification for such a substitution based on the dipole field geometry can be found, for example, in Trammell et al. (2011). One can see in Figure 11 that the outflow velocity and density in the case of increased divergence are significantly smaller. The maximum temperature realized in the case of inflated flow is ∼1000 K lower than that for the spherically symmetric isotropic flow. The sonic point is reached at about 8 R_p and farther than in the equatorial case, contrary to expectations.

Zoom In Zoom Out Reset image size