Multi-Component MHD Model of Hot Jupiter Envelopes

Abstract

A numerical model description of a hot Jupiter extended envelope based on the approximation of multi-component magnetic hydrodynamics is presented. The main attention is focused on the problem of implementing the completed MHD stellar wind model. As a result, the numerical model becomes applicable for calculating the structure of the extended envelope of hot Jupiters not only in the super-Alfvén and sub-Alfvén regimes of the stellar wind flow around and in the trans-Alfvén regime. The multi-component MHD approximation allows the consideration of changes in the chemical composition of hydrogen–helium envelopes of hot Jupiters. The results of calculations show that, in the case of a super-Alfvén flow regime, all the previously discovered types of extended gas-dynamic envelopes are realized in the new numerical model. With an increase in magnitude of the wind magnetic field, the extended envelope tends to become more closed. Under the influence of a strong magnetic field of the stellar wind, the envelope matter does not move along the ballistic trajectory but along the magnetic field lines of the wind toward the host star. This corresponds to an additional (sub-Alfvénic) envelope type of hot Jupiters, which has specific observational features. In the transient (trans-Alfvén) mode, a bow shock wave has a fragmentary nature. In the fully sub-Alfvén regime, the bow shock wave is not formed, and the flow structure is shock-less.

1. Introduction

Hot Jupiters are giant exoplanets with masses on the order of Jupiter’s mass, located in the immediate vicinity of a host star [1]. The first hot Jupiter was discovered in 1995 [2]. Due to the close location to the host star and the relatively large size, gas envelopes of hot Jupiters can overfill their Roche lobes, resulting in intense gas outflows both at the night side (near the Lagrange point $L_2$) and at the day side (near the inner Lagrange point $L_1$) of the planet [3,4]. The presence of such outflows is indirectly indicated by the excessive absorption of radiation in the near ultraviolet range observed in some hot Jupiters during their transit across the disk of their host star [5,6,7,8,9,10,11]. These conclusions are confirmed by direct numerical calculations in the framework of one-dimensional aeronomic models [1,12,13,14,15].

In this case, the matter of expanding the upper atmosphere of hot Jupiter is no longer collision-less, and a gas-dynamic approximation can be used to describe it. Such an exosphere is more correctly called the extended envelope of a hot Jupiter. These are located above the exobase, have a fairly large size, relatively high density, and are characterized by significant deviations from a spherical shape. The structure of an extended envelope and its physical properties are determined by the fact that several other forces act on each element of the flow in the envelope in addition to the planet gravity: the gravitational force of the star, the orbital centrifugal force, the orbital Coriolis force, as well as the forces determined by interaction with the stellar wind, the radiation of the star, and the magnetic field. Therefore, the study of the structure of gas envelopes of such objects is one of the most urgent problems of modern astrophysics [16].

In the three-dimensional approximation, the gas-dynamic structure of the extended envelope of hot Jupiters was studied in [17,18,19,20,21,22]. In these works, it is shown that, depending on the model parameters, structures of three main types can be formed during the interaction of stellar wind with the expanding envelope of a hot Jupiter [17]. If the atmosphere of a planet is completely located inside its Roche lobe, then a closed envelope is formed. If the flow of matter from the inner Lagrange point $L_1$ is stopped by the dynamic pressure of stellar winds, then a quasi-closed envelope is formed. Finally, if the dynamic pressure of the stellar wind is not sufficient to stop outflow from the Lagrange point $L_1$, an open envelope is formed. The type of forming envelope significantly determines the mass loss rate of the hot Jupiter [17]. On the other hand, in the works [23,24,25,26,27], noticeable dependence of the mass loss rate by the planet on the strength of stellar wind wasn’t found, while the type of envelope changed from the open type to the closed one.

Hot Jupiters may have their own magnetic field, which can influence the process of their flow around by stellar wind. However, estimates of the intrinsic magnetic fields of hot Jupiters show that it is most likely quite weak. The characteristic value of the magnetic moment of hot Jupiters, apparently, is 0.1–0.2 ${\mathsf{\mu }}_{\mathrm{J}}$, where ${\mathsf{\mu }}_{\mathrm{J}}=1.53\times {10}^{30}$ G·cm³ is the Jupiter magnetic moment. This value is in a good agreement with both observational [28,29,30,31] and theoretical [32] estimates.

The relatively low value of the dipole moment is explained by the inefficiency of the dynamo process of magnetic field generation in the bowels of these planets. This is due to the fact that hot Jupiters are located close to the host star; therefore, due to strong tidal disturbances, the proper rotation of a typical hot Jupiter should pass into a state of synchronization with its orbital motion over a period of about several million years [33]. In the state of synchronous rotation, the efficiency of dynamo generation of magnetic field drops sharply.

The process of dynamo-generation of magnetic field occurs in the bowels of the planet and is largely determined by its internal structure (see, for example, [34,35]). However, it should be noted that the magnetic field of hot Jupiters can be generated not only in the bowels but also in the upper layers of atmosphere. The estimates made in [21] show that the upper atmosphere of hot Jupiters consists of almost completely ionized gas. This is due to the processes of thermal ionization and hard radiation of the host star.

Therefore, the upper part of hot Jupiter atmosphere can be called the ionospheric envelope. It is shown in [36] that the proper magnetic field of hot Jupiter should influence the formation of large-scale (zonal) currents in its atmosphere. Detailed three-dimensional calculations [37,38] demonstrate a complex picture of the distribution of winds in the upper atmosphere in which magnetic fields play an important role.

In particular, electromagnetic forces can shift the hot spot forming at point facing the host star to the west. This effect can also be manifested on the light curves of hot Jupiters. For example, a comparison of the observed light curves with those calculated for the planet HAT-P-7b allows estimation of the characteristic magnetic field in the atmosphere of 6 G [39]. Most likely, this estimate is greatly overvalued. Recall that, at the level of the Jupiter cloud layer, the magnitude of the magnetic field is 4–5 G.

It is important to note one more circumstance. Due to their close location to the host star, hot Jupiters can have a quite strong magnetic field induced by the stellar wind magnetic field. As the calculations presented in [40] show, the corresponding magnetic moment of such a field can be from 10% to 20% of the Jupiter magnetic moment. Summing up all these comments, we can conclude that the question of magnitude and configuration of magnetic field in hot Jupiters is still open.

Analysis of the influence of stellar wind magnetic field on the process of its flow around the atmosphere of hot Jupiter [41] shows that in the case of hot Jupiters, this effect can be extremely important. This is due to the fact that almost all hot Jupiters are located in the so-called sub-Alfvén zone of stellar wind of the host star, where the wind speed is less than the Alfvénic one. Taking into account the orbital motion of the planet, the flow velocity turns out to be close to the Alfvénic one.

Therefore, depending on the specific situation (the major semi-axis of the planet orbit, the spectral class of host star, the proper rotation period of star, the features of stellar wind), both super-Alfvén and sub-Alfvén flow regimes can be realized. Note that in the super-Alfvén regime, the magnetosphere of hot Jupiter will contain all the main elements (the bow shock wave, the magnetopause, the magnetospheric tail at the night side, and others) present in the magnetospheres of planets of the Solar system [42,43]. In the case of sub-Alfvén flow regime, the bow shock wave will be absent in the magnetosphere structure [44].

The problem of the magnetic field influence on the dynamics of extended envelopes of hot Jupiters is of constant great interest to the scientific community. In recent years, several attempts have been made to account for the magnetic field in one-dimensional [45,46,47,48], in two-dimensional [49] and in three-dimensional [47,50,51] numerical aeronomic models of hot Jupiter atmospheres. However, in these studies, only the immediate vicinity of the planet was considered, and estimates of the mass loss rate were performed without considering the presence of extended envelopes.

An exception is the work [51], in which the authors performed three-dimensional numerical modeling in a wide spatial domain and obtained MHD solutions for exoplanets with open and quasi-closed envelopes. In our recent work, we investigated the influence of both the planet own magnetic field [52,53] and the wind field [41,54,55,56,57] on the envelope dynamics of hot Jupiters. The results of these investigations are presented in the review [16].

In this paper, we consider the further development of our numerical MHD codes. The main attention in the paper is focused on obtaining a numerical model that is self-consistent with the stellar wind and includes additional physical effects due to the complex chemical composition of the hydrogen–helium envelopes of hot Jupiters. Taking into account the self-consistent model with the wind allows, in particular, to more correctly determine the position of the Alfvén point.

In our previous numerical models, this parameter was either estimated approximately from the condition of constancy in the computational domain of the radial wind speed, or was set separately. Another important development direction is the creation of a numerical model of multi-fluid (multi-component) magnetic hydrodynamics for describing the structure of hot Jupiter extended envelope. In the paper, we did not consider specific processes that lead to local changes in the concentration of components (chemical reactions, ionization, etc.), as well as the corresponding heating–cooling processes, since the main goal was to include the stellar wind model.

However, in the future, this will allow for not only consideration of the chemical composition of the hydrogen–helium envelopes of hot Jupiters but also to trace the distributions and dynamics of components of particular interest (e.g., biomarkers). The multi-fluid MHD model described in this paper was developed on the basis of already existing single-fluid MHD model for describing the structure of a hot Jupiter extended envelope [16].

The paper is organized as follows. Section 2 describes the MHD model of stellar wind we used. Section 3 describes the three-dimensional numerical model of a multi-component MHD. Section 4 describes the numerical model of the hot Jupiter extended envelope based on multi-component MHD. Section 5 presents the results of calculations. The main conclusions of the study are summarized in Section 6. Finally, some computation method details are described in the Appendix A.

2. Stellar Wind Model

2.1. Basic Equations

To describe the structure (including the magnetic field) of stellar wind in the vicinity of hot exoplanets in our numerical model, we will rely on the well-studied properties of the solar wind. As shown in numerous ground- and space-based investigations (see, for example, the review of [58]), the picture of magnetic field of the solar wind is fairly complex. Schematically, this structure is illustrated in Figure 1 (see, e.g., our paper [41]). In the corona region, the magnetic field is mainly determined by the Sun intrinsic magnetism, and therefore it is essentially non-radial.

At the boundary of the corona, which is located at a distance of several radii of the Sun, magnetic field becomes purely radial with high accuracy. Beyond this region is the heliospheric area, the magnetic field that is significantly determined by the properties of the solar wind. In this region, the magnetic field lines with a distance from the center gradually twist in the form of a spiral due to rotation of the Sun, and therefore (especially at large distances) the magnetic field of the wind can be described with good accuracy using a simple Parker model [59].

However, the observed magnetic field in the solar wind is not axisymmetric but has a manifested sector structure. This is due to the fact that at different points of the spherical surface of the corona, the field may have different polarity (the direction of field lines relative to the direction of normal vector), for example, due to the inclination of magnetic axis of the Sun to its rotation axis.

As a result, two clearly distinguished sectors with different directions of the magnetic field are formed in the solar wind in the ecliptic plane. In one sector, the magnetic field lines are directed towards the Sun, and in the opposite sector, away from the Sun. These two sectors are separated by the heliospheric current sheet, which rotates together with the Sun and therefore the Earth, during its orbit around the Sun, crosses it many times per year, passing from the solar wind sector with one magnetic field polarity to the neighboring sector with the opposite magnetic field polarity.

In our calculations, we neglect the possible sector structure of the wind magnetic field, as well as the presence of a heliospheric current sheet in it, focusing on the influence of its global parameters. We reasonably assume that the orbit of a hot exoplanet is located in the heliospheric region beyond the boundary of the corona. In Figure 1, it is shown as a dotted circle. To describe the structure of the wind (including its magnetic field) in the heliospheric region, an axisymmetric (or even a spherically-symmetric) model [60] can be used as the first approximation.

We will consider the wind model in an inertial reference frame in spherical coordinates (r, $\theta$, $\phi$). We assume that the center of the spherical coordinate system coincides with the center of the star. At the same time, we can ignore the dependence of the wind parameters on the angle $\theta$, since we are interested in the flow structure only near the orbit plane of a hot exoplanet. Therefore, for simplicity, we will assume that all values depend only on the radial coordinate r.

The stationary structure of the wind under such conditions is determined by the continuity equation

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2\rho v_r\right)=0,$$

(1)

the equation of motion for the radial $v_r$ and azimuthal $v_{\phi }$ components of the velocity vector $\mathit{v}$

$$v_r\frac{d{v}_r}{dr}-\frac{v_{\phi }^2}{r}=-\frac{1}{\rho }\frac{dP}{dr}-\frac{G{M}_{\mathrm{s}}}{r^2}-\frac{B_{\phi }}{4\pi \rho r}\frac{d}{dr}\left(r{B}_{\phi }\right),$$

(2)

$$v_r\frac{d{v}_{\phi }}{dr}+\frac{v_r{v}_{\phi }}{r}=\frac{B_r}{4\pi \rho r}\frac{d}{dr}\left(r{B}_{\phi }\right),$$

(3)

the equation of induction

$$\frac{1}{r}\frac{d}{dr}\left(r{v}_r{B}_{\phi }-r{v}_{\phi }{B}_r\right)=0$$

(4)

and the Maxwell equation ($\nabla ·\mathit{B}=0$)

$$\frac{1}{r^2}\frac{d}{dr}\left(r^2{B}_r\right)=0.$$

(5)

Here, $\rho$ represents the density, P is the pressure, G is the gravitational constant, and $M_{\mathrm{s}}$ is the mass of the central star. Density, pressure, and temperature satisfy the equation of state for an ideal polytropic gas,

$$P=K_{\kappa }{\rho }^{\kappa }=\frac{2k_{\mathrm{B}}}{m_{\mathrm{p}}}\rho T,$$

(6)

where $K_{\kappa }$ is the constant, $\kappa$ is the polytropic index, $k_{\mathrm{B}}$ is the Boltzmann constant, and $m_{\mathrm{p}}$ is the proton mass. The average molecular weight of the wind matter is considered to be equal to $1/2$, which corresponds to a fully ionized hydrogen plasma consisting only of electrons and protons.

From the Maxwell Equation (5), we find

$$r^2{B}_r=B_{\mathrm{s}}{R}_{\mathrm{s}}^2,$$

(7)

where $R_{\mathrm{s}}$ is the radius of the star, and $B_{\mathrm{s}}$ is the magnitude of field at the star surface. From the continuity Equation (1), we can obtain an integral of motion corresponding to the law of mass conservation,

$$4\pi r^2\rho v_r={\dot{M}}_{\mathrm{s}},$$

(8)

where the integration constant ${\dot{M}}_{\mathrm{s}}$ determines the mass loss rate by star due to the outflow of matter in the form of a stellar wind. The Equation (2) can be rewritten as:

$$\frac{d}{dr}\left(\frac{v_r^2}{2}+\frac{c_s^2}{\kappa -1}-\frac{G{M}_{\mathrm{s}}}{r}+\frac{B_{\phi }^2}{4\pi \rho }\right)-\frac{v_{\phi }^2}{r}-r{B}_{\phi }\frac{d}{dr}\left(\frac{B_{\phi }}{4\pi \rho r}\right)=0,$$

(9)

where the square of the speed of sound

$$c_s^2=\frac{\kappa P}{\rho }=\kappa K_{\kappa }{\rho }^{\kappa -1}.$$

(10)

We can multiply the Equation (3) by $v_{\phi }$ and divide by $v_r$. Then, it can be represented as:

$$\frac{d}{dr}\left(\frac{v_{\phi }^2}{2}-\frac{B_r{B}_{\phi }^2{v}_{\phi }}{4\pi \rho v_r}\right)+\frac{v_{\phi }^2}{r}+r{B}_{\phi }\frac{d}{dr}\left(\frac{B_r{v}_{\phi }}{4\pi \rho r{v}_r}\right)=0.$$

(11)

Summing the Equations (9) and (11), we find:

$$\begin{array}{c}\frac{d}{dr}\left(\frac{v_r^2}{2}+\frac{v_{\phi }^2}{2}+\frac{c_s^2}{\kappa -1}-\frac{G{M}_{\mathrm{s}}}{r}+\frac{B_{\phi }^2}{4\pi \rho }-\frac{B_r{B}_{\phi }^2{v}_{\phi }}{4\pi \rho v_r}\right)\hfill \\ +r{B}_{\phi }\frac{d}{dr}\left(\frac{B_r{v}_{\phi }}{4\pi \rho r{v}_r}-\frac{B_{\phi }}{4\pi \rho r}\right)=0.\hfill \end{array}$$

(12)

The last term on the left side of this equation turns to zero, because

$$\frac{d}{dr}\left(\frac{B_r{v}_{\phi }}{4\pi \rho r{v}_r}-\frac{B_{\phi }}{4\pi \rho r}\right)=\frac{1}{{\dot{M}}_{\mathrm{s}}}\frac{d}{dr}\left(r{B}_r{v}_{\phi }-r{B}_{\phi }{v}_r\right)=0.$$

(13)

Here, we used the law of conservation of mass (8) and the equation of induction (4). This circumstance allows us to derive the integral of motion corresponding to the law of energy conservation,

$$\frac{v_r^2}{2}+\frac{v_{\phi }^2}{2}+\frac{c_s^2}{\kappa -1}-\frac{G{M}_{\mathrm{s}}}{r}+\frac{B_{\phi }^2}{4\pi \rho }-\frac{B_r{B}_{\phi }{v}_{\phi }}{4\pi \rho v_r}=Q_{\mathrm{s}},$$

(14)

where the constant $Q_{\mathrm{s}}$ determines the density of energy flux in the stellar wind, the value of which is ${\dot{M}}_{\mathrm{s}}{Q}_{\mathrm{s}}$.

Note that, from (7) and (8) follows

$$\frac{B_r}{4\pi \rho v_r}=\frac{B_{\mathrm{s}}{R}_{\mathrm{s}}^2}{{\dot{M}}_{\mathrm{s}}}=\mathrm{const}.$$

(15)

This circumstance allows from the Equations (3) and (4) to find two another integrals of motion:

$$r{v}_{\phi }-\frac{B_r}{4\pi \rho v_r}r{B}_{\phi }=L,$$

(16)

$$r{v}_r{B}_{\phi }-r{v}_{\phi }{B}_r=F.$$

(17)

The first integral determines the law of conservation of angular momentum. The second integral is related to the vertical component of the electric field $E_{\theta }$, since it is obvious that $F=c{E}_{\theta }r$, where c is the speed of light. For the convenience of further calculations, we introduce the following notation:

$$a_r=\frac{B_r}{\sqrt{4\pi \rho }},a_{\phi }=\frac{B_{\phi }}{\sqrt{4\pi \rho }}.$$

(18)

Then, from Equations (16) and (17), it can be obtained that:

$$\left(v_r^2-a_r^2\right)v_{\phi }=\frac{1}{r}\left(v_r^{2}L+\frac{a_{r}F}{\sqrt{4\pi \rho }}\right),$$

(19)

$$\left(v_r^2-a_r^2\right)a_{\phi }=\frac{v_r}{r}\left(a_{r}L+\frac{F}{\sqrt{4\pi \rho }}\right).$$

(20)

The value of the constant F in the last integral of motion (17) can be found from the boundary conditions at some point $r=r_0$. However, in the reference frame rotating with the Sun, the value of this constant should be equal to zero. This is due to the fact that, in this frame of reference, the vector of the magnetic field ${\mathit{B}}^{\prime }$ must be collinear with the velocity vector ${\mathit{v}}^{\prime }$ of the ideally conducting wind plasma: ${\mathit{v}}^{\prime }\times {\mathit{B}}^{\prime }=0$. For non-relativistic motion, the magnetic field ${\mathit{B}}^{\prime }=\mathit{B}$, and the velocity

$$v_r^{\prime }=v_r,v_{\phi }^{\prime }=v_{\phi }-{\Omega }_{\mathrm{s}}r,$$

(21)

where ${\Omega }_{\mathrm{s}}$ is the angular velocity of the proper rotation of the star. Therefore,

$$F=-{\Omega }_{\mathrm{s}}{r}^2{B}_r=-{\Omega }_{\mathrm{s}}{R}_{\mathrm{s}}^2{B}_{\mathrm{s}}.$$

(22)

Substituting the expression (22) into the Equations (19) and (20), we find:

$$\left(v_r^2-a_r^2\right)v_{\phi }=\frac{1}{r}\left(v_r^{2}L-a_r^2{\Omega }_{\mathrm{s}}{r}^2\right),$$

(23)

$$\left(v_r^2-a_r^2\right)a_{\phi }=\frac{v_r{a}_r}{r}\left(L-{\Omega }_{\mathrm{s}}{r}^2\right).$$

(24)

In the obtained expressions, there is a singularity at some point $r=r_A$, where the radial wind velocity $v_r$ becomes equal to the Alfvén one

$$u_A=|a_r|=\frac{|B_r|}{\sqrt{4\pi \rho }},$$

(25)

that corresponding to the radial component of magnetic field $B_r$. Let us call this special point as a Alfvén point. Near the surface of the star, the radial wind velocity $v_r$ should be less than the Alfvén one $u_A$. At large distances, the radial velocity $v_r$, on the contrary, exceeds the Alfvén one $u_A$. At the Alfvén point $r=r_A$, there is a transition from the sub-Alfvén flow regime to the super-Alfvén one. Therefore, the region $r<r_A$ can be called the sub-Alfvén zone of the stellar wind, and the region $r>r_A$, respectively, the super-Alfvén zone.

The azimuthal components of the velocity $v_{\phi }$ and the magnetic field $B_{\phi }$ in the expressions (23) and (24) must remain continuous at the Alfvén point. Therefore, in order to satisfy this condition, it is necessary to set the value of the integration constant

$$L={\Omega }_{\mathrm{s}}{r}_A^2.$$

(26)

As a result, we find the final solution

$$v_{\phi }=\frac{{\Omega }_{\mathrm{s}}}{r}\frac{v_r^2{r}_A^2-a_r^2{r}^2}{v_r^2-a_r^2},$$

(27)

$$a_{\phi }=\frac{{\Omega }_{\mathrm{s}}}{r}{v}_r{a}_r\frac{r_A^2-r^2}{v_r^2-a_r^2}.$$

(28)

Let us introduce the Alfvén Mach number for the radial components of the velocity and the magnetic field,

$$\lambda =\frac{\sqrt{4\pi \rho }{v}_r}{B_r}.$$

(29)

It is not difficult to make sure that

$${\lambda }^2=\frac{v_r{r}^2}{v_A{r}_A^2}=\frac{{\rho }_A}{\rho },$$

(30)

where $v_A$ and ${\rho }_A$ are the values of the radial velocity and wind density at the Alfvén point. For numerical modeling, a convenient record of expressions for the azimuthal components of the velocity and magnetic field are

$$v_{\phi }={\Omega }_{\mathrm{s}}r\frac{1-{\lambda }^2{r}_A^2/r^2}{1-{\lambda }^2},$$

(31)

$$B_{\phi }=B_r\frac{{\Omega }_{\mathrm{s}}r}{v_r}{\lambda }^2\frac{1-r_A^2/r^2}{1-{\lambda }^2}.$$

(32)

Note that in the sub-Alfvén zone of the stellar wind $\lambda <1$, while in the super-Alfvén zone $\lambda >1$. At the Alfvén point, $\lambda =1$. The obtained relations, considering the integrals of mass (8) and energy (14), allow us to algebraically find the distributions of all magnetohydrodynamic quantities describing the wind structure.

2.2. Method of Solution

From the Equations (1)–(5) using simple manipulations, it can be derived that

$$\left(v_r^2-u_S^2\right)\left(v_r^2-u_F^2\right)\frac{r}{v_r}\frac{d{v}_r}{dr}=\left(v_r^2-a_r^2\right)\left(2c_s^2+v_{\phi }^2-\frac{G{M}_{\mathrm{s}}}{r}\right)+2v_r{v}_{\phi }{a}_r{a}_{\phi },$$

(33)

where

$$u_{S,F}^2=\frac{1}{2}\left[c_s^2+a^2\mp \sqrt{{(c_s^2+a^2)}^2-4c_s^2{u}_A^2}\right],$$

(34)

$$a^2=a_r^2+a_{\phi }^2.$$

(35)

The quantities $u_S$ and $u_F$ determine the values of slow and fast magnetosonic velocities, respectively. It can be seen from the Equation (33) that there are two special points in the solution: slow magnetosonic$r=r_S$, in which the speed $v_r=u_S$, and fast magnetosonic$r=r_F$, in which the speed $v_r=u_F$. At these points, the coefficient for the derivative $d{v}_r/dr$ turns to zero. In order for the solution to remain smooth, the right side of the Equation (33) at these points must also vanish. However, as we saw above, the Alfvén point $r=r_A$, in which the wind speed $v_r=u_A$, is also a solution critical point. This circumstance can be expressed directly by substituting the relations (27) and (33) into the Equation (28). As a result, we find:

$$\begin{array}{c}{\left(v_r^2-u_A^2\right)}^2\left(v_r^2-u_S^2\right)\left(v_r^2-u_F^2\right)\frac{r}{v_r}\frac{d{v}_r}{dr}={\left(v_r^2-u_A^2\right)}^3\left(2c_s^2-\frac{G{M}_{\mathrm{s}}}{r}\right)\hfill \\ +\frac{{\Omega }_{\mathrm{s}}^2}{r^2}\left(v_r^2{r}_A^2-u_A^2{r}^2\right)\left(v_r^4{r}_A^2+u_A^4{r}^2-3v_r^2{u}_A^2{r}^2+v_r^2{u}_A^2{r}_A^2\right).\hfill \end{array}$$

(36)

Further, we have

$$\begin{array}{c}\left(v_r^2-u_S^2\right)\left(v_r^2-u_F^2\right)=v_r^4-(c_s^2+u_A^2)v_r^2+c_s^2{u}_A^2-a_{\phi }^2{v}_r^2\hfill \\ =\left(v_r^2-c_s^2\right)\left(v_r^2-u_A^2\right)-a_{\phi }^2{v}_r^2.\hfill \end{array}$$

(37)

Therefore, considering (28)

$${\left(v_r^2-u_A^2\right)}^2\left(v_r^2-u_S^2\right)\left(v_r^2-u_F^2\right)=\left(v_r^2-c_s^2\right){\left(v_r^2-u_A^2\right)}^3-\frac{{\Omega }_{\mathrm{s}}^2}{r^2}{v}_r^4{u}_A^2{(r_A^2-r^2)}^2.$$

(38)

Now, the equation for the radial velocity (33) can be written as:

$$\begin{array}{c}\left[\left(v_r^2-c_s^2\right){\left(v_r^2-u_A^2\right)}^3-\frac{{\Omega }_{\mathrm{s}}^2}{r^2}{v}_r^4{u}_A^2{(r_A^2-r^2)}^2\right]\frac{r}{v_r}\frac{d{v}_r}{dr}\hfill \\ ={\left(v_r^2-u_A^2\right)}^3\left(2c_s^2-\frac{G{M}_{\mathrm{s}}}{r}\right)\hfill \\ +\frac{{\Omega }_{\mathrm{s}}^2}{r^2}\left(v_r^2{r}_A^2-u_A^2{r}^2\right)\left(v_r^4{r}_A^2+u_A^4{r}^2-3v_r^2{u}_A^2{r}^2+v_r^2{u}_A^2{r}_A^2\right).\hfill \end{array}$$

(39)

It is known from observations that in the solar wind, the radial velocity monotonically increases with the distance from the Sun. Therefore, we are interested in a solution with an accelerating wind, when the radial velocity gradient has a positive value at any distance from the star, $d{v}_r/dr>0$. This means that in the Equation (33), the coefficient for the derivative in the left part must change sign simultaneously with the right part.

Let us study the asymptotic behavior of the solution of interest to us at small distances in the limit at $r\to 0$. Suppose that the radial velocity changes in this case according to the power law $v_r=A{r}^{\sigma }$, where A is a certain coefficient. Consider the case when the radial velocity $v_r$ is much smaller than any of the characteristic velocities $u_S$, $u_F$ and $u_A$. Then, the coefficient for the derivative on the left side of the Equation is (33)

$$\left(v_r^2-u_S^2\right)\left(v_r^2-u_F^2\right)=u_S^2{u}_F^2=c_s^2{u}_A^2.$$

(40)

From the relations (27) and (28) at our limit, we find

$$v_{\phi }={\Omega }_{\mathrm{s}}r,a_{\phi }=-a_r\frac{{\Omega }_{\mathrm{s}}r}{v_A}.$$

(41)

It follows that the right-hand side of the Equation (33) is approximately equal to

$$u_A^2\frac{G{M}_{\mathrm{s}}}{r}-2v_r{u}_A^2\frac{{\Omega }_{\mathrm{s}}^2{r}^2}{v_A}\approx u_A^2\frac{G{M}_{\mathrm{s}}}{r}.$$

(42)

As a result, neglecting non-essential terms, we come to the equation:

$$c_s^2\frac{r}{v_r}\frac{d{v}_r}{dr}=\frac{G{M}_{\mathrm{s}}}{r}.$$

(43)

Taking into account the law of mass conservation (8) and the expression for the speed of sound (10), we find the values of constants

$$\sigma =\frac{3-2\kappa }{\kappa -1},A=\frac{{\dot{M}}_{\mathrm{s}}}{4\pi }{\left[\frac{3-2\kappa }{\kappa -1}\frac{\kappa K_{\kappa }}{G{M}_{\mathrm{s}}}\right]}^{\frac{1}{\kappa -1}}.$$

(44)

Since in such a solution the index of $\sigma$ must be positive, then the index of polytrope $\kappa <3/2$. Thus, the solution with an accelerating wind (a positive value of the radial velocity gradient, $d{v}_r/dr>0$) is realized only in the case of polytropic index $\kappa <3/2$. This value is less than the adiabatic index $\gamma =5/3$. Therefore, effective sources of heating must be present in the stellar wind.

In fact, the entropy of an ideal gas $s=c_{V}ln(P/{\rho }^{\gamma })$, where $c_V$ is the specific heat capacity at a constant volume. Therefore, the full derivative

$${\displaystyle \frac{ds}{dt}}=c_V\left({\displaystyle \frac{d}{dt}}lnP-\gamma {\displaystyle \frac{d}{dt}}ln\rho \right).$$

(45)

Taking into account the polytropic equation of state (6), we can write

$$T{\displaystyle \frac{ds}{dt}}=-c_{V}T(\gamma -\kappa )\frac{1}{\rho }{\displaystyle \frac{d\rho }{dt}}.$$

(46)

The expression on the right-hand side of this equation determines the heating function $\Gamma$. Since the specific internal energy of an ideal gas $\epsilon =c_{V}T$, then using the continuity Equation (1), we obtain:

$$\Gamma =(\gamma -\kappa )\epsilon \frac{1}{r^2}{\displaystyle \frac{d}{dr}}(r^2{v}_r).$$

(47)

We can see that in the case of an accelerating wind ($d{v}_r/dr>0$) and $\kappa <\gamma$ this value is positive. Effective heating in the solar wind is determined by the fundamental driving mechanism. An overview of possible wind driving mechanisms for stars of various types can be found in [61].

As boundary conditions in our wind model, we will use the values of density ${\rho }_0$, radial velocity $v_0$ and temperature $T_0$ at some point $r_0$. Using these values, we can calculate the speed of sound

$$c_0=\sqrt{\frac{2k_{\mathrm{B}}{T}_0}{m_{\mathrm{p}}}}$$

(48)

and the integration constant

$${\dot{M}}_{\mathrm{s}}=4\pi r_0^2{\rho }_0{v}_0.$$

(49)

The last expression follows from the law of mass conservation (8). The value of polytropic index $\kappa$, which lies in the range $1<\kappa <3/2$, is considered as unknown. Its specific value must be such that the boundary conditions are satisfied.

The law of mass conservation (8) also implies the relation

$${\dot{M}}_{\mathrm{s}}=4\pi r_A^2{\rho }_A{v}_A.$$

(50)

We find from here

$${\rho }_A=\frac{1}{4\pi }{\left(\frac{{\dot{M}}_{\mathrm{s}}}{r_0^2{B}_0}\right)}^2,$$

(51)

where by

$$B_0=B_{\mathrm{s}}{\left(\frac{R_{\mathrm{s}}}{r_0}\right)}^2$$

(52)

the value of the radial magnetic field at the point $r_0$ is indicated. We introduce dimensionless quantities $x_A=r_A/r_0$ and $y_A=v_A/v_0$. It is not difficult to verify that they satisfy the relation

$$y_A=\frac{\beta }{x_A^2},$$

(53)

where is the dimensionless parameter

$$\beta =\frac{{\rho }_0}{{\rho }_A}=\frac{B_0}{v_0\sqrt{4\pi {\rho }_A}}.$$

(54)

To solve the equations describing the wind structure, it is convenient to rewrite them in a dimensionless form. To do this, we denote the dimensionless variables $\xi =r/r_A$ and $\eta =v_r/v_A$. In this case, the Alfvén Mach number (29) turns out to be equal to $\lambda =\xi \sqrt{\eta }$. In addition, we denote

$$z={\left(\frac{\rho }{{\rho }_0}\right)}^{\frac{\kappa -1}{2}}={\left(\beta {\lambda }^2\right)}^{\frac{1-\kappa }{2}}.$$

(55)

In dimensionless variables, the Equation (33) takes the form

$$\frac{\xi }{\eta }\frac{d\eta }{d\xi }=\frac{Y}{X},$$

(56)

where

$$X={({\lambda }^2-1)}^3\left(y_A^2{\eta }^2-\alpha z^2\right){\xi }^2-\omega {\lambda }^4{x}_A^2{({\xi }^2-1)}^2,$$

(57)

$$Y={({\lambda }^2-1)}^3\left(2\alpha \xi z^2-\frac{\mu }{x_A}\right)\xi +\omega x_A^2({\lambda }^2-{\xi }^2)\left({\lambda }^4-3{\lambda }^2{\xi }^2+{\lambda }^2+{\xi }^2\right),$$

(58)

$$\alpha =\frac{c_0^2}{v_0^2},\mu =\frac{G{M}_{\mathrm{s}}}{r_0{v}_0^2},\omega =\frac{{\Omega }_{\mathrm{s}}^2{r}_0^2}{v_0^2}.$$

(59)

As was noted above, the slow magnetosonic ($\xi ={\xi }_S$, $\eta ={\eta }_S$) and fast magnetosonic ($\xi ={\xi }_F$, $\eta ={\eta }_F$) points are critical solution points. However, the Alfvén point ($\xi =1$, $\eta =1$) is also critical, since in it the values X and Y also simultaneously turn to zero.

The expression determining the conservation of energy (14) in dimensionless variables can be written as follows:

$$Z=\frac{y_A^2{\eta }^2}{2}+\frac{\alpha z^2}{\kappa -1}-\frac{\mu }{x_A\xi }+\frac{\omega x_A^2}{2{\xi }^2{({\lambda }^2-1)}^2}\left[{\lambda }^4+(2{\lambda }^2-1)({\xi }^2-2){\xi }^2\right]-q=0,$$

(60)

where it is denoted

$$q=\frac{Q_{\mathrm{s}}}{v_0^2}.$$

(61)

At the point $r_0$ we have ${\xi }_0=1/x_A$, ${\eta }_0=1/y_A$, ${\lambda }_0=1/\sqrt{\beta }$, $z_0=1$. Therefore, at this point, Equation (60) takes the form:

$$Z_0=\frac{1}{2}+\frac{\alpha }{\kappa -1}-\mu +\frac{\omega }{2{({\lambda }_0^2-1)}^2}\left[{\lambda }_0^4{x}_A^4+(2{\lambda }_0^2-1)(1-2x_A^2)\right]-q=0.$$

(62)

We will use ${\xi }_S$, ${\eta }_S$, ${\xi }_F$, ${\eta }_F$, $x_A$, q and $\kappa$ as unknown quantities in solving the problem. In total, we have 7 unknowns. The system includes four equations $X=0$, $Y=0$ at the points ${\xi }_S$ and ${\xi }_F$, two Equation (60) written at the points ${\xi }_S$ and ${\xi }_F$, as well as the Equation (62). This system of non-linear algebraic equations can be solved numerically by iterations. The algorithm for constructing the solution is as follows. First, we numerically solve the system of equations for the wind parameters ${\xi }_S$, ${\eta }_S$, ${\xi }_F$, ${\eta }_F$, $x_A$, q and $\kappa$. After these parameters are determined, for each value of $\xi$, we numerically solve the Equation (60) and find the dependence $\eta (\xi )$ corresponding to the curve passing through all three critical points.

2.3. Calculation Example

Let us consider the results of a demonstration calculation of the stellar wind structure. The model parameters were the values of density $n_0=1400{\mathrm{cm}}^{-3}$, velocity $v_0=130\mathrm{km}/\mathrm{s}$, and temperature $T_0=7.3\times {10}^5\mathrm{K}$ at a distance of $r_0=10R_{\odot }$ from the star [62]. For example, let us take a typical hot Jupiter HD 209458b, which is used in our numerical simulations described below. The host star is characterized by the following parameters: spectral class G0, mass $M_{\mathrm{s}}=1.1M_{\odot }$, and radius $R_{\mathrm{s}}=1.2R_{\odot }$. The period of star proper rotation $P_{\mathrm{rot}}=14.4$ days, which corresponds to the angular velocity ${\Omega }_{\mathrm{s}}=5.05\times {10}^{-6}{\mathrm{s}}^{-1}$ or the linear velocity at the equator $v_{\mathrm{rot}}=4.2\mathrm{km}/\mathrm{s}$.

The major semi-axis of the planet orbit $A=10.2R_{\odot }$, which is close to the selected boundary point $r_0$ and corresponds to the period of rotation around the star $P_{\mathrm{orb}}=84.6$ h. The magnetic field at the stellar surface was set to $B_{\mathrm{s}}=0.5\mathrm{G}$, which approximately corresponds to the value of average magnetic field at the surface of the Sun in a quiet state $B_{\odot }=1\mathrm{G}$, if we compare the corresponding magnetic fluxes, $B_{\mathrm{s}}{R}_{\mathrm{s}}^2\approx B_{\odot }{R}_{\odot }^2$.

It should be noted that the parameter $B_{\mathrm{s}}$ corresponds to the average value of magnetic field on the star surface only formally. The solution under consideration is valid only in the heliospheric area (see Figure 1). In the corona region, the star intrinsic magnetic field plays a decisive role. Therefore, strictly speaking, the value of $B_{\mathrm{s}}$ is obtained from the average value of the radial magnetic field $B_{\mathrm{c}}$ at the corona boundary $r=R_{\mathrm{c}}$, recalculated by considering the conservation of the magnetic flux by the star radius, $B_{\mathrm{s}}{R}_{\mathrm{s}}^2=B_{\mathrm{c}}{R}_{\mathrm{c}}^2$.

On the other hand, it is known that the magnetic fields of solar-type stars can lie in the range from about 0.1 to several Gauss [63,64]. In addition, the host stars for hot Jupiters can be not only of the solar type, since their spectral classes lie in the wide range from class F to class M. In addition to the radial component, the azimuthal component of the magnetic field is also present in the stellar wind, which is due to the star proper rotation. The angular velocity of the star proper rotation ${\Omega }_{\mathrm{s}}$, in turn, also depends on the spectral class [64]. These remarks significantly expand the set of possible models of the stellar wind in vicinity of hot Jupiters.

Hence, the following parameter values are obtained: mass loss rate ${\dot{M}}_{\mathrm{s}}=1.85\times {10}^{11}\mathrm{g}/\mathrm{s}=2.94\times {10}^{-15}{M}_{\odot }/\mathrm{year}$, density at the Alfvén point ${\rho }_A=2.25\times {10}^{-22}\mathrm{g}/{\mathrm{cm}}^3$, dimensionless parameters $\alpha =0.713$, $\beta =10.424$, $\mu =1.241$, $\omega =0.0725$. Note that the magnitude of the mass loss rate of the star ${\dot{M}}_{\mathrm{s}}$ in our stellar wind model is only one of its parameters. It does not coincide with the real mass loss rate, since we solved the problem only in the plane of planet orbit. The real wind does not have spherical symmetry, and therefore these values can differ several times.

The results of calculations of the wind structure are shown in Figure 2. A family of integral curves corresponding to different types of solutions is shown. Some of these curves pass through the critical points, the positions of which are indicated by circles. In this case, the Alfvén critical point in the plane of the variables $\xi$, $\eta$ has coordinates (1, 1). As a solution corresponding to the stellar wind, it is necessary to use an integral curve passing through all three critical points. The corresponding curve in Figure 2 is shown as a bold line.

The numerical solution of a system of non-linear algebraic equations describing the stellar wind structure gives the following coordinates of magnetosonic singular points: ${\xi }_S=0.328$, ${\eta }_S=0.514$, ${\xi }_F=1.067$, ${\eta }_F=1.025$. On the left panel of Figure 2, it can be seen that a slow magnetosonic point is a critical point of the “saddle” type. The fast magnetosonic point is located near the Alfvén point, the vicinity of which is shown on an enlarged scale on the right panel of Figure 2.

It can be seen that a fast magnetosonic point is also a critical point of the “saddle” type. The Alfvén critical point has a more complex topology, since it has a higher singularity order. For the Alfvén point, the value $x_A=2.511$ is obtained, which corresponds to the radius $r_A=25.114R_{\odot }$. The polytropic index turned out to be $\kappa =1.055$, and the parameter $q=12.774$. The polytropic index is close to unity. Therefore, the stellar wind in vicinity of hot Jupiter can be considered almost iso-thermal. This result is in good agreement with the observations.

It is known that, at short distances from the Sun ($r<15R_{\odot }$), the effective adiabatic index $\gamma =1.1$ [65,66]. At large distances $r>25R_{\odot }$, the effective adiabatic index can be estimated by the value $\gamma =1.46$ [67]. Since the orbits of hot Jupiters are located at close distances from the host star in the region of wind acceleration $r<20R_{\odot }$, the wind structure in vicinity of the planet is found to be almost iso-thermal with good accuracy.

The calculation results are shown in Figure 3 and Figure 4. Figure 3 shows the distributions of the radial velocity $v_r$ (solid line), the speed of sound $c_s$, the Alfvén velocity $u_A$, as well as the slow $u_S$ and fast $u_F$ magnetosonic velocities depending on the radius r. The fast magnetosonic point is located close to the Alfvén point, slightly exceeding it. The slow magnetosonic point coincided with the sound point at which the radial wind speed $v_r=c_s$. The latter circumstance is due to the fact that in this region, because of the small contribution of the azimuthal component of the magnetic field, the slow magnetosonic velocity $u_S\approx min(u_A,c_s)=c_s$.

Figure 4 shows the radial profiles of the particle number density $n(r)$, temperature $T(r)$, azimuthal velocity $v_{\phi }(r)$, as well as the radial $B_r(r)$ (dotted line) and azimuthal $B_{\phi }(r)$ (solid line) components of the magnetic field. The dependence $v_{\phi }(r)$ is non-monotonic. The maximum values of the azimuthal velocity are reached approximately in the area of the hot Jupiter’s location $13R_{\odot }$. However, the characteristic values of $15\mathrm{km}/\mathrm{s}$ are small compared to the radial velocity of $130\mathrm{km}/\mathrm{s}$.

3. Multi-Component Magnetic Hydrodynamics

Let us consider an approximation of multi-component (multi-fluid) magnetic hydrodynamics, which uses equations for mass quantities, the induction equation, as well as the continuity equations for individual plasma components. The individual components (electrons, ions, and neutrals of various kinds) of plasma will be marked with the index s. Denote by $\mathit{v}$ the average mass velocity of matter. Then, the continuity equations for each component of the kind s can be written as:

$${\displaystyle \frac{\partial {\rho }_s}{\partial t}}+\nabla ·\left({\rho }_s\mathit{v}\right)=D_s+S_s.$$

(63)

where ${\rho }_s$ is the density of the component s, the value

$$D_s=\nabla ·\left({\rho }_s{\mathit{w}}_s\right)$$

(64)

determines diffusion, ${\mathit{w}}_s=\mathit{v}-{\mathit{v}}_s$ is the diffusion velocities. The last term in the right hand side of (63) takes into account the source function, which describes changes in the number of particles of the kind s due to chemical reactions, as well as the processes of dissociation, ionization, and recombination. For the total density

$$\rho =\sum _s{\rho }_s$$

(65)

due to the conditions

$$\sum _s{\rho }_s{\mathit{w}}_s=0,\sum _s{S}_s=0,$$

(66)

we have the usual continuity equation,

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho \mathit{v}\right)=0.$$

(67)

Let us introduce the values ${\xi }_s$ that describe the mass fraction of the component s and satisfy the relations:

$${\rho }_s={\xi }_s\rho ,\sum _s{\xi }_s=1.$$

(68)

Then, from the Equation (63), we can find

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho {\xi }_s\right)+\nabla ·\left(\rho {\xi }_{s}v\right)=D_s+S_s.$$

(69)

Note that one of the values ${\xi }_s$ can be excluded, since it can always be expressed in terms of the others using conditions (68). It is convenient to consider electrons as such a component. We will use a separate index for it, and, for the other components, we assume that s runs through the values from 1 to N. Therefore, all values ${\xi }_s$ can be considered independent. The mass fraction of electrons in this case is equal to

$${\xi }_{\mathrm{e}}=1-\sum _s{\xi }_s.$$

(70)

Since the mass of electrons is much smaller than the mass of ions and neutrals, ${\xi }_{\mathrm{e}}$ turns out to be a value close to zero.

The equations of motion and energy for mass quantities can be written in the following form:

$$\rho \left[{\displaystyle \frac{\partial \mathit{v}}{\partial t}}+\left(\mathit{v}·\nabla \right)\mathit{v}\right]+\nabla P+\frac{\mathit{B}\times \nabla \times \mathit{B}}{4\pi }=\rho \mathit{f}-{\mathit{D}}_v,$$

(71)

$$\rho \left[{\displaystyle \frac{\partial \epsilon }{\partial t}}+(\mathit{v}·\nabla )\epsilon \right]+P\nabla ·\mathit{v}=D_{\epsilon }+\rho Q,$$

(72)

where P is the pressure, $\mathit{B}$ is the magnetic field, $\mathit{f}$ is the specific external force, the values

$${\mathit{D}}_v=\sum _s\nabla ·\left({\rho }_s{\mathit{w}}_s{\mathit{w}}_s\right),$$

(73)

$$D_{\epsilon }=\sum _s\left[\nabla ·\left({\rho }_s{\epsilon }_s{\mathit{w}}_s\right)+P_s\nabla ·{\mathit{w}}_s\right]$$

(74)

are determined by the diffusion velocities ${\mathit{w}}_s$, and the value of Q describes the heating–cooling sources. In the expression (74), the notations $P_s$ and ${\epsilon }_s$ are used for the partial pressures and specific internal energies of the components. To these equations, we should add the induction equation describing the evolution of the magnetic field,

$${\displaystyle \frac{\partial \mathit{B}}{\partial t}}-\nabla \times \left(\mathit{v}\times \mathit{B}\right)=-\nabla \times \left(\eta \nabla \times \mathit{B}\right)+\nabla \times \left(v_{\mathrm{D}}\times \mathit{B}\right).$$

(75)

Here, the first term in the right hand side describes the ohmic diffusion of the magnetic field, where $\eta$ is the corresponding magnetic viscosity. The second term determines the effect of ambipolar diffusion occurring in an incompletely ionized plasma. The vector ${\mathit{v}}_{\mathrm{D}}$ represents the speed of ambipolar diffusion. Our numerical code takes into account the effects caused by both magnetic viscosity and ambipolar diffusion. However, in this paper we do not consider them, and we will assume that the values $\eta =0$, ${\mathit{v}}_{\mathrm{D}}=0$.

The density $\rho$, pressure P, internal energy $\epsilon$, and temperature T satisfy the equation of state for an ideal gas

$$P=\frac{k_{\mathrm{B}}}{\mu m_{\mathrm{p}}}\rho T,\epsilon =\frac{k_{\mathrm{B}}T}{\mu m_{\mathrm{p}}(\gamma -1)},$$

(76)

where $k_{\mathrm{B}}$ is the Boltzmann constant, $m_{\mathrm{p}}$ is the proton mass, $\gamma$ is the adiabatic index. The average molecular weight $\mu$ is determined by the expression $\rho =\mu m_{\mathrm{p}}n$, where n is the total number density. Let us write down the condition of quasi-neutrality of plasma

$$-e{n}_{\mathrm{e}}+\sum _s{e}_s{n}_s=0,$$

(77)

where e is the elementary charge, $n_{\mathrm{e}}$ is the electron number density, $e_s$ is the charge of particles of the kind s, and $n_s$ is the number density of the component s. Hence, the number density of electrons

$$n_{\mathrm{e}}=\sum _s{Z}_s{n}_s,$$

(78)

where $Z_s=e_s/e$ is the charge number of particles of the kind s (for neutrals it is zero). Taking into account this expression, we find

$$\frac{1}{\mu m_{\mathrm{p}}}=\sum _s\left(1+Z_s\right)\frac{{\xi }_s}{m_s},$$

(79)

where $m_s$ is the mass of particles of the kind s.

The numerical method we use to solve the equations of multi-component magnetic hydrodynamics is described in the Appendix A.

4. Model for Envelope of Hot Jupiter

4.1. Model Description

The structure and dynamics of plasma flow in the vicinity of hot Jupiter will be described in the approximation of multi-component (multi-fluid) magnetic hydrodynamics, which was discussed in the previous Section 3. In this approximation, it is possible, in particular, to take into account the complex chemical composition of the extended envelope of hot Jupiter. It is convenient to explicitly distinguish the background magnetic field [16,41,52,68,69], when the total magnetic field $\mathit{B}$ is represented as a superposition of the background magnetic field $\mathit{H}$ and the magnetic field $\mathit{b}$ induced by electrical currents in plasma itself, $\mathit{B}=\mathit{H}+\mathit{b}$. This approach allows us to significantly minimize the numerical errors that occur during arithmetic operations with large numbers, provides greater stability of the scheme, and improves the quality of calculations in rarefied regions with a strong magnetic field—for example, in the magnetosphere of hot Jupiter.

In our formulation of the problem, the background field is created by sources distributed inside the star (or, more precisely, inside the corona), as well as in the bowels of the planet. Therefore, these sources are absent in the computational domain, and hence the background field should satisfy the potentiality condition, $\nabla \times \mathit{H}=0$. This allows it to partially exclude from the equations of magnetic hydrodynamics [70,71].

In general case, the background magnetic field is not stationary, $\partial \mathit{H}/\partial t\ne 0$. However, our model assumes that this property is possessed by the magnetic field of the planet. The proper rotation of hot Jupiter, due to strong tidal interactions from a closely located host star, turns out to be synchronized with the orbital rotation. As a result, the period of planet rotation around its proper axis will be equal to its orbital period, and, consequently, a hot Jupiter will always face the same side to its host star.

Therefore, in a rotating frame of reference associated with the orbital motion of the planet, the orientation of its own magnetic field will not change over time. On the other hand, the magnitude of field (for example, the magnetic moment) changes on much larger time scales compared to the orbital period and, in our calculations, this change can be ignored. The same comments can be made about the induced magnetic field of the planet.

As theoretical estimates show (see, e.g., [1,17]), a hydrodynamic approach is applicable to describe the flow of matter near a typical hot Jupiter. This is due to the fact that the extended envelope of hot Jupiter has a sufficiently high density and therefore, its matter is not collision-less. In addition, due to the processes of thermal ionization and hard radiation of the host star, the upper atmosphere and the extended envelope of hot Jupiters consists of almost fully ionized gas [21].

Taking into account the complex chemical composition of the extended envelope, this makes it possible to use the multi-component magnetic hydrodynamics approximation. On the other hand, the analysis of the characteristic collision frequencies [13,46,47,72] of the components in the hydrogen–helium envelope of hot Jupiter allows us to neglect the diffusion effects [48] with fairly good accuracy and assume that the average velocities of all components are equal to the average mass velocity of matter, ${\mathit{v}}_s=\mathit{v}$. We can also assume that all the components are in thermodynamic equilibrium, and therefore their temperatures are equal to the temperature of matter, $T_s=T$.

Considering the above circumstances, the equations of multi-component magnetic hydrodynamics can be written as

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho \mathit{v}\right)=0,$$

(80)

$${\displaystyle \frac{\partial \mathit{v}}{\partial t}}+(\mathit{v}·\nabla )\mathit{v}=-\frac{\nabla P}{\rho }-\frac{\mathit{b}\times \nabla \times \mathit{b}}{4\pi \rho }-\frac{\mathit{H}\times \nabla \times \mathit{b}}{4\pi \rho }+\mathit{f},$$

(81)

$${\displaystyle \frac{\partial \mathit{b}}{\partial t}}=\nabla \times \left(\mathit{v}\times b+\mathit{v}\times \mathit{H}\right)-{\displaystyle \frac{\partial \mathit{H}}{\partial t}},$$

(82)

$${\displaystyle \frac{\partial \epsilon }{\partial t}}+(\mathit{v}·\nabla )\epsilon +\frac{P}{\rho }\nabla ·\mathit{v}=Q,$$

(83)

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho {\xi }_s\right)+\nabla ·\left(\rho {\xi }_s\mathit{v}\right)=S_s,s=1,\dots ,N.$$

(84)

To close this system of equations, the equation of state for an ideal gas (76) with the adiabatic index $\gamma =5/3$ is used. In this paper, we focus on a more correct accounting of the stellar wind. Therefore, in the simulations described below, we assumed that the source functions $S_s$ due to chemical reactions, processes of ionization, recombination and dissociation of molecules, as well as the corresponding contributions to the heating function Q are equal to zero. This implies that various plasma components that we considered as passive admixtures transported together with the matter. In addition, for the same reason, in this work, we did not take into account the effects of magnetic viscosity and ambipolar diffusion.

We assume that the planet orbit is circular. Therefore, in a non-inertial reference frame rotating together with a binary system with an angular velocity $\mathbf{\Omega }$, the locations of the star and the planet centers do not change. In this case, the specific external force is determined by the expression

$$\mathit{f}=-\nabla \Phi -2\left(\mathbf{\Omega }\times \mathit{v}\right).$$

(85)

Here, the first term on the right hand side describes the force due to the gradient of the Roche potential

$$\Phi =-\frac{G{M}_{\mathrm{s}}}{|\mathit{r}-{\mathit{r}}_{\mathrm{s}}|}-\frac{G{M}_{\mathrm{p}}}{|\mathit{r}-{\mathit{r}}_{\mathrm{p}}|}-\frac{1}{2}{\left[\mathbf{\Omega }\times (\mathit{r}-{\mathit{r}}_{\mathrm{c}})\right]}^2,$$

(86)

where $M_{\mathrm{s}}$ is the mass of the star, $M_{\mathrm{p}}$ is the mass of the planet, ${\mathit{r}}_{\mathrm{s}}$ is the radius vector of the star center, ${\mathit{r}}_{\mathrm{p}}$ is the radius vector of the planet center, ${\mathit{r}}_{\mathrm{c}}$ is the radius vector of the center of mass of the system. The second term in the right hand side of (85) describes the Coriolis force.

In this numerical model, the stellar wind is taken into account according to the same scheme as for the purely gas-dynamic case [17]. However, this does not use a constant value of the radial wind velocity but the profile $v_r(r)$ obtained from the wind model (see Section 2). Profiles for the remaining values are calculated using this profile: $\rho (r)$, $v_{\phi }(r)$, $B_{\phi }(r)$. This allows finding the wind parameters at any point $\mathit{r}=(x,y,z)$ of the computational domain. At the same time, in a non-inertial reference frame, the magnetic field of the wind does not change, and the wind velocity is recalculated using the formula

$${\mathit{v}}_{\mathrm{w}}={\mathit{u}}_{\mathrm{w}}-\mathbf{\Omega }\times \left(\mathit{r}-{\mathit{r}}_{\mathrm{c}}\right),$$

(87)

where ${\mathit{u}}_{\mathrm{w}}$ is the wind velocity in the inertial reference frame. These distributions are used to set initial conditions in the region occupied by the stellar wind, as well as to implement boundary conditions. The wind model used assumes a polytropic equation of state (6). The Equations (80)–(84) use the model of adiabatic magnetic hydrodynamics with the adiabatic index $\gamma$. This means that, from the point of view of the adiabatic model, there are effective heating processes in the stellar wind. In order to reconcile these two models, the corresponding function (47) [25] is added to the energy Equation (83)

$$Q_{\mathrm{w}}=(\gamma -\kappa ){\epsilon }_{\mathrm{w}}\nabla ·{\mathit{v}}_{\mathrm{w}},$$

(88)

where ${\epsilon }_{\mathrm{w}}$ is the specific internal energy of the wind matter. This takes into account the expression (87) wind velocity divergence

$$\nabla ·{\mathit{v}}_{\mathrm{w}}=\nabla ·{\mathit{u}}_{\mathrm{w}}=\frac{1}{r^2}{\displaystyle \frac{d}{dr}}\left(r^2{v}_r\right).$$

(89)

Since the stellar wind accelerates (the velocity $v_r$ increases with the distance r from the star), the divergence of its velocity turns out to be positive. Consequently, the function $Q_{\mathrm{w}}>0$.

4.2. Upper Atmosphere

At the initial moment of time, a spherically symmetric iso-thermal atmosphere was set around the planet, the density distribution in which was determined from the hydrostatic equilibrium condition:

$$\rho ={\rho }_{\mathrm{a}}exp\left[-\frac{G{M}_{\mathrm{p}}}{R_{\mathrm{p}}{A}_{\mathrm{gas}}{T}_{\mathrm{a}}}\left(1-\frac{R_{\mathrm{p}}}{|\mathit{r}-{\mathit{r}}_{\mathrm{p}}|}\right)\right].$$

(90)

Here, ${\rho }_{\mathrm{a}}$ is the density at the photometric radius of hot Jupiter $R_{\mathrm{p}}$, $T_{\mathrm{a}}$ is the atmospheric temperature, $A_{\mathrm{gas}}=k_{\mathrm{B}}/(\mu m_{\mathrm{p}})$ is the gas constant, $\mu$ is the average molecular weight. A similar expression can be written for the pressure P. For the hot Jupiter HD 209458b dimensionless parameter,

$$\eta =\frac{G{M}_{\mathrm{p}}}{R_{\mathrm{p}}{A}_{\mathrm{gas}}{T}_{\mathrm{a}}}=10.3\mu {\left(\frac{T_{\mathrm{a}}}{{10}^4\mathrm{K}}\right)}^{-1}.$$

(91)

The initial thickness of the atmosphere was determined from the pressure equilibrium condition with the matter of the stellar wind. The total wind pressure

$$P_{\mathrm{tot}}=P_{\mathrm{w}}+{\rho }_{\mathrm{w}}{v}_{\mathrm{w}}^2+\frac{B_{\mathrm{w}}^2}{8\pi }.$$

(92)

Therefore, from the equality of pressure at the point of a frontal impact, we find the initial radius of the atmosphere

$$R_{\mathrm{a}}={\left(1-\frac{1}{\eta }ln\frac{P_{\mathrm{a}}}{P_{\mathrm{tot}}}\right)}^{-1}{R}_{\mathrm{p}}.$$

(93)

The value $R_{\mathrm{a}}$ is determined to a greater extent by the atmosphere parameters ${\rho }_{\mathrm{a}}$ and $T_{\mathrm{a}}$ and to a lesser extent by the magnetic field of the wind $B_{\mathrm{w}}$. In particular, the ratio between the initial radius of the atmosphere $R_{\mathrm{a}}$ and the size of the Roche lobe of the planet determines the type of extended envelope of the hot Jupiter (closed or open).

We assume that the hydrogen–helium atmosphere of hot Jupiter has a homogeneous chemical composition. In any allocated volume of the atmosphere, the parameter $\chi =[\mathrm{He}/\mathrm{H}]$, equal to the ratio of the helium nuclei number to the hydrogen one number, remains constant. The following significant components of the atmosphere are taken into account [16,73]: molecular hydrogen H₂, atomic hydrogen H, ionized hydrogen H⁺, atomic helium He, once ionized helium He⁺. The mass fraction of hydrogen X and helium Y are determined by the expressions:

$$X=\frac{\rho ({\mathrm{H}}_2)+\rho (\mathrm{H})+\rho ({\mathrm{H}}^{+})}{\rho }=\xi ({\mathrm{H}}_2)+\xi (\mathrm{H})+\xi ({\mathrm{H}}^{+}),$$

(94)

$$Y=\frac{\rho (\mathrm{He})+\rho ({\mathrm{He}}^{+})}{\rho }=\xi (\mathrm{He})+\xi ({\mathrm{He}}^{+}).$$

(95)

From here, we can find: $Y/X=4\chi$. Taking into account the equality $X+Y=1$, we find

$$X=\frac{1}{1+4\chi },Y=\frac{4\chi }{1+4\chi }.$$

(96)

Let us consider the degree of hydrogen ionization $x_1$, the degree of dissociation of molecular hydrogen $x_2$, and the degree of helium ionization $x_3$,

$$x_1=\frac{\rho ({\mathrm{H}}^{+})}{\rho ({\mathrm{H}}^{+})+\rho (\mathrm{H})},x_2=\frac{\rho (\mathrm{H})}{\rho (\mathrm{H})+\rho ({\mathrm{H}}_2)},x_3=\frac{\rho ({\mathrm{He}}^{+})}{\rho ({\mathrm{He}}^{+})+\rho (\mathrm{He})}.$$

(97)

Then, for the mass fractions of the components, it can be written:

$$\xi ({\mathrm{H}}_2)=\frac{(1-x_1)(1-x_2)}{1-x_1+x_2{x}_2}X,$$

(98)

$$\xi (\mathrm{H})=\frac{(1-x_1)x_2}{1-x_1+x_2{x}_2}X,$$

(99)

$$\xi ({\mathrm{H}}^{+})=\frac{x_1{x}_2}{1-x_1+x_2{x}_2}X,$$

(100)

$$\xi (\mathrm{He})=(1-x_3)Y,$$

(101)

$$\xi ({\mathrm{He}}^{+})=x_{3}Y.$$

(102)

Using these definitions, from the expression (79) for the average molecular weight, we find

$$\frac{1}{\mu }=\frac{1}{2}X\left(1+\frac{x_2+2x_1{x}_2}{1+x_1-x_1{x}_2}\right)+\frac{1}{4}Y(1+x_3).$$

(103)

The total degree of ionization of atmospheric matter

$$\xi =\xi ({\mathrm{H}}^{+})+\xi ({\mathrm{He}}^{+})=\frac{x_1{x}_2}{1-x_1+x_1{x}_2}X+x_{3}Y.$$

(104)

In the simulations of the flow structure in the vicinity of the hot Jupiter HD 209458b given below, we set the following parameters of the chemical composition of the atmosphere: $\chi =[\mathrm{He}/\mathrm{H}]=0.05$ [25], $x_1=x_2=x_3=0.9$. From here, we find the mass fraction of hydrogen $X=0.83$ and helium $Y=0.17$. These parameters correspond to the mass fractions of the components $\xi ({\mathrm{H}}_2)=0.009$, $\xi (\mathrm{H})=0.08$, $\xi ({\mathrm{H}}^{+})=0.74$, $\xi (\mathrm{He})=0.02$, $\xi ({\mathrm{He}}^{+})=0.15$. At the same time, the average molecular weight of the atmospheric matter is $\mu =0.69$, and the total degree of ionization is $\xi =0.89$.

The boundary conditions (density ${\rho }_{\mathrm{a}}$, temperature $T_{\mathrm{a}}$ and velocity $v_{\mathrm{a}}$) at the photometric radius were set based on the results of calculations carried out within the framework of one-dimensional aeronomic models considering supra-thermal particles [16,73]. Therefore, in a certain sense, our numerical model is hybrid. It follows from aeronomic calculations that a planetary wind is formed in the upper atmosphere of hot Jupiter under the influence of star hard radiation, which determines the mass loss rate ${\dot{M}}_{\mathrm{a}}\approx {10}^9$–${10}^{10}$ g/s. Taking into account the expression ${\dot{M}}_{\mathrm{a}}=4\pi R_{\mathrm{p}}^2{\rho }_{\mathrm{a}}{v}_{\mathrm{a}}$ from here, we can find the speed of atmospheric outflow $v_{\mathrm{a}}$ at the photometric radius.

4.3. Magnetic Field

The background magnetic field can be set as a superposition of several separate fields:

$$\mathit{H}={\mathit{H}}_{\mathrm{s}}+{\mathit{H}}_{\mathrm{w}}+{\mathit{H}}_{\mathrm{p}}+{\mathit{H}}_{\mathrm{a}}.$$

(105)

Here, ${\mathit{H}}_{\mathrm{s}}$ describes the proper field of the star, ${\mathit{H}}_{\mathrm{w}}$ is determined by the wind field, ${\mathit{H}}_{\mathrm{p}}$ is the proper field of the planet, and ${\mathit{H}}_{\mathrm{a}}$ is determined by the field of atmosphere. Let us characterize the contribution of each term.

As noted in Section 2, the proper magnetic field of the star ${\mathit{H}}_{\mathrm{s}}$ plays a significant role in the corona region. Therefore, if the planet orbit is located in the heliospheric region, then this field can be neglected. The main role in this zone is played by the wind field. The magnetic field of the host star should be taken into account when the planet orbit is located near the outer boundary of the corona or even inside it. Currently, hot Jupiters of this type are observed and in some papers (see, for example, [51]) the field of the star was taken into account explicitly. In our work, we consider the hot Jupiter HD 209458b, whose orbit is located in the heliospheric region. This allows us to neglect the proper field of the host star in the calculations.

The wind magnetic field, which is determined from the model described in Section 2, does not make sense to include (105) entirely in the background field. The fact is that the total magnetic field of the wind ${\mathit{B}}_{\mathrm{w}}$ does not satisfy the condition of potentiality, $\nabla \times {\mathit{B}}_{\mathrm{w}}\ne 0$, since the rotor of the magnetic field just determines the electromagnetic force affecting the dynamics of the wind plasma. However, the radial magnetic field of the wind $B_R$ can be included in the background field. On the one hand, this field, as it is easy to see, satisfies the condition of potentiality. On the other hand, in the vicinity of hot Jupiter, located close to the host star, the radial component of the magnetic field dominates in the wind plasma, since the ratio $B_{\phi }/b_r$ is small in absolute value (see Figure 4). Taking into account the assumed spherical symmetry of the stellar wind, we obtain (see Section 2)

$${\mathit{H}}_{\mathrm{w}}=\frac{B_{\mathrm{s}}{R}_{\mathrm{s}}^2}{{|\mathit{r}-{\mathit{r}}_{\mathrm{s}}|}^2}{\mathit{n}}_{\mathrm{s}},$$

(106)

where the unit vector ${\mathit{n}}_{\mathrm{s}}=(\mathit{r}-{\mathit{r}}_{\mathrm{s}})/|\mathit{r}-{\mathit{r}}_{\mathrm{s}}|$ determines the direction from the center of star ${\mathit{r}}_{\mathrm{s}}$ to the observation point $\mathit{r}$.

In our numerical model, we assume that the intrinsic magnetic field of hot Jupiter is purely dipole,

$${\mathit{H}}_{\mathrm{p}}=\frac{{\mu }_{\mathrm{p}}}{{|\mathit{r}-{\mathit{r}}_{\mathrm{p}}|}^3}\left[3({\mathit{d}}_{\mathrm{p}}·{\mathit{n}}_{\mathrm{p}}){\mathit{n}}_{\mathrm{p}}-{\mathit{d}}_{\mathrm{p}}\right],$$

(107)

where ${\mu }_{\mathrm{p}}$ is the magnetic moment of the planet, ${\mathit{n}}_{\mathrm{p}}=(\mathit{r}-{\mathit{r}}_{\mathrm{p}})/|\mathit{r}-{\mathit{r}}_{\mathrm{p}}|$, ${\mathit{d}}_{\mathrm{p}}$ is a unit vector directed along the magnetic axis and determining the vector of magnetic moment ${\mathit{\mu }}_{\mathrm{p}}={\mu }_{\mathrm{p}}{\mathit{d}}_{\mathrm{p}}$. In our calculations, we assumed the value of magnetic moment of the hot Jupiter HD 209458b to be ${\mu }_{\mathrm{p}}=0.1{\mu }_{\mathrm{J}}$. The orientation of magnetic dipole was determined by the angles $\theta$ and $\varphi$, which were used as model parameters. The components of the unit vector ${\mathit{d}}_{\mathrm{p}}$ directed along the magnetic axis in the Cartesian coordinate system are described by the expressions:

$${\mathit{d}}_{\mathrm{p}}=\left(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta \right).$$

(108)

At the same time, we assumed that the proper rotation of the planet is synchronized with the orbital one, and the axes of proper and orbital rotation are collinear.

For hot Jupiters, the magnetic field induced in the upper atmosphere and the extended envelope seems to play an important role. Since plasma is almost completely ionized in this region, any movements in it lead to the appearance of electric currents and, consequently, to the generation of the magnetic field. This field, in particular, can arise due to zonal currents in the upper atmosphere caused by uneven heating from the radiation of the host star [39]. On the other hand, due to its close location to the host star, a quite strong magnetic field can arise in the upper atmosphere of hot Jupiter, induced by currents shielding the magnetic field of the stellar wind [40]. The configuration of induced magnetic field is such that inside the atmosphere it completely compensates for the magnetic field of the wind, and outside it is dipole. Therefore, in the region beyond the initial atmosphere, we can write

$${\mathit{H}}_{\mathrm{a}}=\frac{{\mu }_{\mathrm{a}}}{{|\mathit{r}-{\mathit{r}}_{\mathrm{p}}|}^3}\left[3({\mathit{d}}_{\mathrm{a}}·{\mathit{n}}_{\mathrm{p}}){\mathit{n}}_{\mathrm{p}}-{\mathit{d}}_{\mathrm{a}}\right],$$

(109)

where the corresponding magnetic moment ${\mu }_{\mathrm{a}}=R_{\mathrm{a}}^3{B}_{\mathrm{w}}/2$, and the unit vector ${\mathit{d}}_{\mathrm{a}}=-{\mathit{B}}_{\mathrm{w}}/B_{\mathrm{w}}$ is directed against the vector ${\mathit{B}}_{\mathrm{w}}$. The magnitude of induced magnetic moment ${\mu }_{\mathrm{a}}$ depends on the initial radius of the atmosphere $R_{\mathrm{a}}$ and the wind field $B_{\mathrm{w}}$ in the orbit of the planet. If we take an average field of the star $B_{\mathrm{s}}=0.5$ Gs, then, for the characteristic dimensions of the atmosphere $R_{\mathrm{a}}$ = (4–5)$R_{\mathrm{p}}$, we can obtain the values of magnetic moment ${\mu }_{\mathrm{a}}$ = (0.05–0.2)${\mu }_{\mathrm{J}}$. In other words, in order of magnitude, the induced magnetic moment ${\mu }_{\mathrm{a}}$ turns out to be equal to the intrinsic magnetic moment ${\mu }_{\mathrm{p}}$ of the hot Jupiter. Note that the induced magnetic field ${\mathit{H}}_{\mathrm{a}}$ is completely determined by the wind field ${\mathit{B}}_{\mathrm{w}}$. Therefore, the structure of induced field (in particular, the vector of induced magnetic moment ${\mathit{\mu }}_{\mathrm{a}}={\mu }_{\mathrm{a}}{\mathit{d}}_{\mathrm{a}}$) will track the direction to the star when the planet moves along its orbit.

As noted above, in a non-inertial reference frame rotating together with a binary system consisting of a star and a planet, the magnetic fields ${\mathit{H}}_{\mathrm{p}}$ and ${\mathit{H}}_{\mathrm{a}}$ are stationary, $\partial {\mathit{H}}_{\mathrm{p}}/\partial t=0$, $\partial {\mathit{H}}_{\mathrm{a}}/\partial t=0$. The proper field of a star, on the contrary, is non-stationary, $\partial {\mathit{H}}_{\mathrm{s}}/\partial t\ne 0$, since it rotates together with the star at the angular velocity ${\mathbf{\Omega }}_{\mathrm{s}}-\mathbf{\Omega }$. The radial component of the wind magnetic field ${\mathit{H}}_{\mathrm{w}}$ is also rigidly associated with the star. In this case, it can be assumed that the field lines velocity of the radial magnetic field in a non-inertial reference frame is equal to

$${\mathit{v}}_{\mathrm{s}}=-\mathbf{\Omega }\times \left(\mathit{r}-{\mathit{r}}_{\mathrm{s}}\right),$$

(110)

since the azimuthal component of the field, due to the proper rotation of the star ${\mathbf{\Omega }}_{\mathrm{s}}$, has already been taken into account in the vector $\mathit{b}$. The change in the magnetic field ${\mathit{H}}_{\mathrm{w}}$ in time is determined by the equation

$${\displaystyle \frac{\partial {\mathit{H}}_{\mathrm{w}}}{\partial t}}=\nabla \times \left({\mathit{v}}_{\mathrm{s}}\times {\mathit{H}}_{\mathrm{w}}\right).$$

(111)

Substituting the expressions (106) and (110) here, it is not difficult to verify by direct calculations that the right hand side of this equation vanishes due to the spherical symmetry of the field ${\mathit{H}}_{\mathrm{w}}$. Consequently, the left hand side should also be equal to zero, $\partial {\mathit{H}}_{\mathrm{w}}/\partial t=0$. Thus, in our model, the total background magnetic field (105) turns out to be stationary, $\partial \mathit{H}/\partial t=0$.

4.4. Numerical Method

To numerically solve the equations of multi-component magnetic hydrodynamics (80)–(84), we use a combination of difference schemes of Roe (see Section 3) and Lax–Friedrichs [74,75]. The solution algorithm consists of several successive stages resulting from the application of the splitting method by physical processes. Suppose that we know the distribution of all values on the computational mesh at the moment of time t. Then, to obtain the values at the next time moment $t+\Delta t$, we decompose the complete system of Equations (80)–(84) into two subsystems.

The first subsystem corresponds to the equations of ideal multi-component magnetic hydrodynamics with intrinsic magnetic field $\mathit{b}$ of plasma without considering the background magnetic field $\mathit{H}$:

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho \mathit{v}\right)=0,$$

(112)

$${\displaystyle \frac{\partial \mathit{v}}{\partial t}}+\left(\mathit{v}·\nabla \right)\mathit{v}=-\frac{\nabla P}{\rho }-\frac{\mathit{b}\times \nabla \times \mathit{b}}{4\pi \rho }+\mathit{f},$$

(113)

$${\displaystyle \frac{\partial \mathit{b}}{\partial t}}=\nabla \times \left(\mathit{v}\times \mathit{b}\right),$$

(114)

$${\displaystyle \frac{\partial \epsilon }{\partial t}}+\left(\mathit{v}·\nabla \right)\epsilon +\frac{P}{\rho }\nabla ·\mathit{v}=Q,$$

(115)

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho {\xi }_s\right)+\nabla ·\left(\rho {\xi }_s\mathit{v}\right)=0,s=1,\dots ,N.$$

(116)

In our numerical model, a method based on Roe’s scheme was used to solve this system (see Appendix A).

The second subsystem corresponds to considering the influence of the background field:

$${\displaystyle \frac{\partial \mathit{v}}{\partial t}}=-\frac{\mathit{H}\times \nabla \times \mathit{b}}{4\pi \rho },$$

(117)

$${\displaystyle \frac{\partial \mathit{b}}{\partial t}}=\nabla \times \left(\mathit{v}\times \mathit{H}\right).$$

(118)

The first equation in this subsystem describes the influence of the electromagnetic force caused by the background field, and the second equation does the generation of a magnetic field. At the same time, it is assumed that at this stage of the algorithm, the density $\rho$ and the specific internal energy $\epsilon$ do not change. To solve the second subsystem, the Lax–Friedrichs scheme was used with increasing TVD (total variation diminishing) corrections [68].

To clear the divergence of the magnetic field $\mathit{b}$, we used the method of the generalized Lagrange multiplier [76]. The choice of this method is due to the fact that the flow in vicinity of hot Jupiter is essentially non-stationary, especially in the flow on the night side forming the magnetospheric tail.

5. Results of Simulations

5.1. Model Parameters

As an object of study, we use a typical hot Jupiter HD 209458b, which has the mass $M_{\mathrm{p}}=0.71M_{\mathrm{J}}$ and the photometric radius $R_{\mathrm{p}}=1.38R_{\mathrm{J}}$, where $M_{\mathrm{J}}$ and $R_{\mathrm{J}}$ is the mass and radius of Jupiter. The host star is characterized by the following parameters: spectral class G0, mass $M_{\mathrm{s}}=1.1M_{\odot }$, radius $R_{\mathrm{s}}=1.2R_{\odot }$. The period of proper rotation of the star $P_{\mathrm{rot}}=14.4$ days, which corresponds to the angular velocity ${\Omega }_{\mathrm{s}}=5.05\times {10}^{-6}$ s⁻¹ or linear velocity at the equator $v_{\mathrm{rot}}=4.2$ km/s. The major semi-axis of the planet orbit $A=10.2R_{\odot }$, which corresponds to the period of revolution around the star $P_{\mathrm{orb}}=84.6$ h.

In the simulations, the temperature of atmosphere $T_{\mathrm{a}}$ is varied, while the particle number density at the photometric radius was set equal to a fixed value $n_{\mathrm{a}}={10}^{10}$ cm⁻³. The atmospheric outflow rate ${\dot{M}}_{\mathrm{a}}={10}^9$ g/s corresponds to the velocity of the planetary wind at the photometric radius $v_{\mathrm{a}}=78$ cm/s. The values of these parameters correspond to the results obtained in aeronomic models for atmospheres of hot Jupiters [12,14,15,16,48,73]. We assumed that the magnitude of the magnetic moment ${\mu }_{\mathrm{p}}$ of the planet is 0.1 of the magnetic moment of Jupiter, and the orientation of magnetic dipole axis (108) is determined by the angle values $\theta$ = 90° and $\phi$ = 60°.

At the same time, as already mentioned above, we believed that the proper rotation of the planet is synchronized with the orbital one, and the axes of proper and orbital rotations are collinear. The average magnetic field at stellar surface $B_{\mathrm{s}}$ in various models was set to 0.01 G (weak field), 0.2 G (medium field), and 0.5 G (strong field). These values are in the physically acceptable range (see, e.g., [63]).

Figure 5 shows the distributions of the radial wind velocity $v_r$ (left panel) and the heating function $Q_{\mathrm{w}}$, which is determined by the Equation (88) (right panel), for all three cases of values $B_{\mathrm{s}}$. The characteristic value of the heating function in the vicinity of the planet $Q_{\mathrm{w}}=5\times {10}^9$ erg·g⁻¹·s⁻¹. In this case, the maximum values of the heating function are reached at distances of approximately $4R_{\odot }$ from the star center.

For numerical simulations, we used the Cartesian coordinate system (x, y, z) in a non-inertial reference frame rotating together with the binary system “star–planet” around the center of masses. The origin of the coordinate system was chosen at the planet center ${\mathit{r}}_{\mathrm{p}}=(0,0,0)$. The x axis was located along the line connecting the centers of the star and the planet, while the center of the star was at the point ${\mathit{r}}_{\mathrm{s}}=(-A,0,0)$. The y axis was chosen so that its direction coincided with the direction of the orbital motion of the planet. Finally, considering the selected orientations of the axes x and y, the third axis z will coincide in the direction with the vector of the orbital angular velocity $\mathbf{\Omega }$.

In order to check the correctness of considering the complete wind model, we performed separate numerical calculations of the flow structure in the region where the planet is absent. To do this, it was enough to choose a computational domain symmetrically located relative to the center of the star (i.e., in the vicinity of the point $x=-2A$). Calculation results for cases of weak ($B_{\mathrm{s}}=0.01$ G) and strong ($B_{\mathrm{s}}=0.5$ G) magnetic field of the wind are shown in Figure 6 and Figure 7, respectively.

These figures show the distributions of density (color and iso-lines), velocity (arrows) and magnetic field (solid lines) at the initial moment of time (left panels) and after one orbital period (right panels). The density values are given in units of magnitude ${\rho }_{\mathrm{w}}=2.3\times {10}^{-21}$ g/cm³, corresponding to the wind density at a distance of $10R_{\odot }$ from the center of the star. The dotted line shows the boundary of the Roche lobe. The star is located to the right of the computational domain.

The analysis of these figures allows us to conclude that during the time of order of the orbital period, in the numerical model, the distributions of wind parameters for both the case of weak and strong fields practically do not change. Small perturbations introduced into the solution by boundary conditions are not essential for our purposes. Therefore, we can assume that the accounting of the complete wind model is carried out correctly.

Figure 8 shows the initial distributions of density (color) and magnetic field (arrow lines) in the vicinity of the planet for the case, when the temperature of the atmosphere $T_{\mathrm{a}}=6000$ K. The left panel corresponds to the magnitude of the magnetic field at the surface of the star $B_{\mathrm{s}}=0.01$ G (weak field), and the right panel does to the case for the strong field $B_{\mathrm{s}}=0.5$ G. The dotted line again shows the boundary of the Roche lobe and the white circle corresponds to the photometric radius of the planet. The star is located on the left side. The radius of the atmosphere in the case of the strong field (right panel) is smaller compared to the case of the weak field (left panel). This is due to the fact that an increase in the field $B_{\mathrm{s}}$ leads to an increase in the total wind pressure (92) and, consequently, to a decrease in the initial radius of the atmosphere (93).

It can be noticed that the magnetic field is clearly divided into four zones. In the first zone (above and below the planet in the figures), the magnetic field is characterized by open force lines of the star, which begin at its surface and go to infinity. In the second zone (to the right of the planet), the magnetic field is determined by the open force lines of the planet. In the third zone (to the left of the planet), magnetic lines are common for the star and the planet. They start at the surface of the star and finish at the surface of the planet, forming a kind of magnetic “bridge” connecting the star with the planet and performing orbital rotation with them together. Finally, the last zone consists of the closed lines of the planet forming the inner part of the magnetosphere.

There are two neutral points in the orbital plane, in which, due to the superposition of individual fields, the total induction vector $\mathit{B}=0$ and therefore the direction of the magnetic field becomes indeterminate. In space, the set of these neutral points forms a certain line, the shape of which, in particular, is determined by the orientation parameters of the magnetic axis of the planet. In case of a strong wind field (Figure 8, right panel) neutral points and closed dipole lines are located in a more compact region around the planet. When the stellar wind flows around the planet, a more complex flow pattern is formed and the structure of the magnetic field is significantly distorted. However, in general, the described topology (separation into magnetic zones and the presence of neutral points) is preserved.

Below, we present the results of two blocks of numerical simulations. In the first block, we set a weak field of the wind corresponding to the super-Alfvén regime of the stellar wind flowing around hot Jupiter, and varied the temperature of the atmosphere. This led to the formation of various types of super-Alfvén envelopes. In the second block, with fixed atmospheric parameters, we varied the magnitude of the magnetic field of the wind in order to trace how the envelope structure changes in this case. The parameters of the models, as well as the flow characteristics, are presented in

Table 1. Parameters and characteristics of the models: $T_{\mathrm{a}}$ is the temperature of the upper atmosphere, $B_{\mathrm{s}}$ is the magnetic field of the star, and ${\dot{M}}_{\mathrm{p}}$ is the mass loss rate of the planet.

Model	${\mathit{T}}_{\mathbf{a}}$, K	${\mathit{B}}_{\mathbf{s}}$, G	Flow Regime	Envelope Type	${\dot{\mathit{M}}}_{\mathbf{p}}$, 109 g/s
1	5000	0.01	Super-Alfvénic	Closed	1
2	5500	0.01	Super-Alfvénic	Quasi-closed	1
3	6000	0.01	Super-Alfvénic	Quasi-Opened	4
4	6500	0.01	Super-Alfvénic	Opened	10
5	6500	0.2	Trans-Alfvénic	Quasi-Opened	9
6	6500	0.5	Sub-Alfvénic	Quasi-Opened	7

.

5.2. Super-Alfvén Flow Regime

This section presents the results of numerical simulation of the flow structure in the vicinity of hot Jupiter for the case, when the magnetic field at the surface of the star $B_{\mathrm{s}}=0.01$ G. The Alfvén Mach number (29) in the orbit of the planet is $\lambda =15.5$. This means that the planet is located in the super-Alfvén zone of the stellar wind. The total Alfvén Mach number, considering the velocity of the orbital motion of the planet, $\lambda =23.2$. Therefore, in this case, the flow around of hot Jupiter by the stellar wind should occur in the super-Alfvén regime.

Calculations were performed for four models differing in atmospheric temperature values $T_{\mathrm{a}}$ (see Table 1). Namely, the temperature of atmosphere $T_{\mathrm{a}}$ was set equal to 5000 K (model 1), 5500 K (model 2), 6000 K (model 3), 6500 K (model 4). The simulation was carried out in the computational domain $-30\le x/R_{\mathrm{p}}\le 30$, $-30\le y/R_{\mathrm{p}}\le 30$, $-15\le z/R_{\mathrm{p}}\le 15$ with the number of cells $192\times 192\times 96$.

The simulation results are demonstrated by Figure 9. The density distributions (color, iso-lines), velocities (arrows) and magnetic field (lines) in the orbital plane are represented on various panels of this figure. The density is expressed in units of wind density in the vicinity of the planet ${\rho }_{\mathrm{w}}$. The boundary of the Roche lobe is shown by a dotted line. The planet is located in the center of the computational domain and is represented by a light circle, the radius of which corresponds to the photometric one. All the solutions obtained correspond to the time moment of the order of one third of the orbital period from the beginning of the computation. During this period of time, a stable quasi-stationary flow pattern is formed.

In all four calculation variants, as a result of the stellar wind flowing around, a wide (on the order of several radii of the planet) hydrogen–helium turbulent plume is formed at the night side of the planet. The interaction of the stellar wind with the envelope of hot Jupiter leads to the appearance of bow shock wave. The position and shape of the bow shock in this numerical model are slightly different from those we received in our previous models [16,17,41,52,54,55,56,57,73]. This is due to a more correct consideration of wind parameters. In the old model, wind velocity and temperature were considered as constants, but in the new model they change in space.

In addition, we previously accelerated the wind by disabling the gravitational forces of the star and the planet in the area occupied by the wind plasma. As a result, the flow velocity increased and the shock wave pressed closer to the planet. In the new numerical model, the flow velocity is calculated from the wind structure and as a result, the shock wave moves further away from the planet. Another effect is due to the fact that the sound point is located inside the computational domain at a distance of about 20 radii of the planet towards the star (see Figure 3). Therefore, the lower left edge of the shock wave, reaching this point, breaks off. This is clearly visible on the upper panels of Figure 9.

In model 1 ($T_{\mathrm{a}}=5000$ K, upper left panel Figure 9) a compact envelope of hot Jupiter is formed, except a relatively weak plume located inside the Roche lobe. The bow shock has an almost spherical shape. According to the classification proposed in [17], such a flow configuration corresponds to the closed envelope of hot Jupiter. In model 2 ($T_{\mathrm{a}}=5500$ K, upper right panel Figure 9) a compact envelope of hot Jupiter is also forming. However, in this case, there is a small cusp in the surface of contact discontinuity (envelope boundary) directed towards the inner Lagrange point $L_1$. Therefore, the outflow of matter from the envelope occurs not only due to a turbulent plume forming at the night side, but also due to a weak outflow from the day side of the planet. The shape of bow shock is also close to spherical. An extended envelope of this type can be called quasi-closed. The mass loss rate for these models ${\dot{M}}_{\mathrm{p}}\approx {10}^9$ g/s and determines by the outflow from the night side of the planet (see the last column in Table 1).

A rather complex flow pattern is observed in model 3 ($T_{\mathrm{a}}=6000$ K, lower left panel Figure 9) and model 4 ($T_{\mathrm{a}}=6500$ K, lower right panel Figure 9). In these models, two powerful flows are formed from the vicinity of Lagrange points $L_1$ and $L_2$. The first stream, as in the previous two models, begins at the night side and forms a wide turbulent plume behind the planet. The second stream is formed at the day side, directed towards the star and, therefore, moves against the wind due to the stellar gravity.

A stream of hydrogen–helium matter from the inner Lagrange point significantly distorts the shape of the bow shock, while pushing it further away from the planet. We can say that the shock wave consists of two separate parts, one of which arises around the atmosphere of the planet, and the other arises around the stream from the inner Lagrange point $L_1$. The flow of matter in the stream is directed against the wind and therefore the Kelvin–Helmholtz instability develops along its surface.

In model 3, the flow has a lower intensity and it is largely deflected and blocked by the stellar wind. Such an extended envelope can be called quasi-open. The corresponding mass loss rate ${\dot{M}}_{\mathrm{p}}\approx 4\times {10}^9$ g/s and determines mainly by the outflow from inner Lagrange point $L_1$. Finally, in model 4, the intensity of the flow is so high that it does not stop under impact of the stellar wind and continues to move away towards the star. This leads to a significant loss of matter from the hot Jupiter envelope. We find that in this model the mass loss rate ${\dot{M}}_{\mathrm{p}}\approx {10}^{10}$ g/s. An extended envelope of this type can be called open.

Thus, in the new numerical model, all the previously discovered [17] types of extended envelopes remain. However, their boundaries in the parameter space have shifted to the region of lower densities and temperatures. This is mainly due to the fact that, in this model, we take into account the chemical composition of the atmosphere.

5.3. Sub-Alfvén Flow Regime

In this section, we present the results of the second block of numerical simulations. The temperature of atmosphere was set equal to $T_{\mathrm{a}}=6500$ K, which corresponds in the case of a weak field ($B_{\mathrm{s}}=0.01$ G) to an open extended envelope under conditions of super-Alfvén flow by the stellar wind (model 4 from the previous Section 5.2). In the second block, calculations were performed for the cases $B_{\mathrm{s}}=0.2$ G (model 5) and $B_{\mathrm{s}}=0.5$ G (model 6). The parameters of the computational domain and the mesh were used the same as in the previous models 1–4. The results of the simulations are shown in Figure 10, which shows the flow structure for these two variants at a time moment of about one third of the orbital period from the beginning of computation.

In both models, the planet is located in the sub-Alfvén wind zone, but between the points $r_S$ (slow magnetosonic) and $r_A$ (Alfvén). Consequently, in the orbit of the planet, the wind velocity $v_w$ exceeds the slow magnetosonic velocity $u_S$, but becomes less than the Alfvén velocity $u_A$. In model 5, the Alfvén Mach number at the planet orbit is $\lambda =0.77$ and therefore the planet is located in the sub-Alfvén zone of the stellar wind. The total Alfvén Mach number, considering the velocity of the orbital motion of the planet $\lambda =1.16$.

We can say that this situation corresponds to the trans-Alfvén regime of the stellar wind flowing around the planet, since hot Jupiter is located in the sub-Alfvén wind zone, but the flow occurs in the super-Alfvén regime. In model 6, the Alfvén Mach number on the orbit of the planet is $\lambda =0.31$, and the total Alfvén Mach number is $\lambda =0.46$. Therefore, in this case, the conditions for the implementation of the sub-Alfvén regime of the hot Jupiter flow around by the stellar wind are completely satisfied.

In the super-Alfvén flow around regime, an extended open-type envelope was formed for these atmospheric parameters (see the lower right panel in Figure 9). Under the conditions of a strong magnetic field of the wind, the envelope structure has undergone significant changes. The intensity of the matter flow from the inner Lagrange point $L_1$ at the day side weakened and the envelope became closer to the quasi-open type. We obtained that the mass loss rate ${\dot{M}}_{\mathrm{p}}\approx 9\times {10}^9$ g/s for the model 5 and ${\dot{M}}_{\mathrm{p}}\approx 7\times {10}^9$ g/s for the model 6. It follows from these results that, at fixed parameters of the atmosphere (density ${\rho }_{\mathrm{a}}$ and temperature $T_{\mathrm{a}}$), the rate of mass loss ${\dot{M}}_{\mathrm{p}}$ for hot Jupiter decreases with an increase in the magnetic field of the wind (the parameter $B_{\mathrm{s}}$ in our numerical model).

In addition, the flow direction changed, since, in a strong field, plasma tends to move mainly along magnetic force lines. This is especially noticeable in the case of model 6 (right panel in Figure 10). Since the wind field is almost radial in the vicinity of the planet, the envelope matter moves directly to the star and, continuing to move along such a trajectory, falls immediately onto the star.

Thus, in the case of a strong magnetic wind field (sub-Alfvén flow around regime), we have some new type of extended envelope, complementing the previous classification [17]. In order to distinguish these envelopes types, we can call them super-Alfvén and sub-Alfvén, respectively. For example, in this case, sub-Alfvén extended envelopes of open or quasi-open type are formed. It is obvious that the observational manifestations of such envelopes should have important differences in comparison with ones formed in the super-Alfvén flow around regime [56].

It is also interesting to note that a magnetic barrier is formed in front of the planet in the direction of its orbital motion. It manifests itself as an elongated area of condensation of magnetic field lines (see Figure 10). Within this region, the magnetic field induction increases.

The process of the stellar wind flowing around the planet is basically shock-less. In model 5 (the left panel in Figure 10) a weak shock wave is formed around the planet, but towards the end of the stream it disappears, since the conditions of the sub-Alfvén flow regime are already being realized in this region. In model 6 (the right panel in Figure 10) in the entire computational domain, the flow occurs in the sub-Alfvén mode. Therefore, shock waves are not formed either around the atmosphere of the planet, or around the flow of matter flowing from the inner Lagrange point $L_1$.

In model 6 (the right panel in Figure 10), a non-physical spherical structure appeared around the planet. Its appearance is due to the fact that to describe the induced magnetic field of the atmosphere, we used the formula (109), which follows from the exact solution for the problem of the magnetic field of an ideally conducting sphere placed in a homogeneous external magnetic field. In our case, the magnetic field of the wind is not uniform. Therefore, this solution does not allow to accurately compensate the external field in the entire volume of the atmosphere and parasitic currents will be induced all the time near the surface of the initial contact discontinuity. This is noticeable in models with a strong wind field. In future works, we will attempt to overcome this problem.

6. Conclusions

To investigate the process of the stellar wind matter flowing around hot Jupiters, considering both the planet own magnetic field and the wind magnetic field, we developed a three-dimensional numerical model based on the approximation of multi-component magnetic hydrodynamics. Our numerical model is based on the Roe–Einfeldt–Osher difference scheme of high resolution for the equations of multi-component MHD.

The total magnetic field is represented as a superposition of the external magnetic field and the magnetic field induced by electric currents in the plasma itself. The superposition of the planet proper magnetic field, the induced magnetic field of the atmosphere, and the radial component of the wind magnetic field was used as the external field. In the numerical algorithm, the factors associated with the presence of an external magnetic field were taken into account at a separate stage using an appropriate Godunov-type difference scheme.

The main attention was focused on the inclusion of a complete MHD model of the stellar wind. This, in particular, makes it possible to more correctly calculate the location of the Alfvén point. As a result, the numerical model is applicable for calculating the structure of the extended envelope of hot Jupiters not only in the super-Alfvén [55] and sub-Alfvén [56] regimes of stellar wind flow around but also in the trans-Alfvén regime. The multi-component MHD approximation in the future will be used by us to account for changes in the chemical composition of the hydrogen–helium envelopes of hot Jupiters. In this paper, we did not consider these processes, assuming all the individual components as passive admixtures moving together with the matter.

As an example of the object to study, we considered a typical hot Jupiter HD 209458b. The numerical simulations presented in the paper can be divided into two blocks. The first block includes calculations with a weak wind field (super-Alfvén flow regime), in which we changed the parameters of atmosphere (temperature). This made it possible to simulate different types of super-Alfvén envelopes: closed, quasi-closed, quasi-open and open. In the second block, with fixed atmospheric parameters (open envelope), we changed the magnetic field of the wind and analyzed how the envelope structure changes.

In the case of the super-Alfvén flow regime, all previously discovered types of extended envelopes are also implemented in the new numerical model. However, their boundaries in the parameter space have shifted to the region of lower densities and temperatures. This is due to the fact that in this model we take into account the chemical composition of the atmosphere. Position and shape of the bow shock in the new numerical model differ from those we obtained in the previous models. This can be explained by a more correct consideration of wind parameters. In the old model, wind velocity and temperatures were considered constant throughout the computational domain, and in the new model they change in space according to the analytical solution.

In addition, in the old model, the gravitational forces of the star and planet were turned off to accelerate the wind. As a result, the flow velocity naturally increased and the shock wave pressed closer to the planet. In the new model, the flow velocity is obtained directly from the wind model (considering all forces) and as a result, the shock wave moves further away from the planet.

With the increase in the magnitude of wind magnetic field, the total wind pressure enlarges. As a result, the stream of matter from the inner Lagrange point $L_1$ is stopped by the stellar wind earlier. Therefore, for example, an open envelope tends to become quasi-open with the growth of field. In addition, in the strong magnetic field of the wind, the direction of movement of stream changes. The stream plasma will attempt to move along the magnetic force lines. As a result, an additional type of envelopes is realized—sub-Alfvén ones, which have their own specific observational features.

In the trans-Alfvén regime, the bow shock wave has a fragmentary nature. In the completed sub-Alfvén flow around regime, the bow shock wave is not formed at all. It should also be noted that with an increase in the magnetic field of the wind, the induced field of the atmosphere growths. The corresponding magnetic moment becomes greater than the planet intrinsic magnetic moment. As a result, the magnetic pole shifts. This may affect the overall configuration of the magnetosphere and, in particular, the position of the dead zones and the auroral zone.