*Article* **On the Origin of Persistent Radio and X-ray Emission from Brown Dwarf TVLM 513-46546**

**Alexander Stepanov 1,2,\* and Valery Zaitsev <sup>3</sup>**


**Abstract:** We study the origin of unusually persistent microwave and X-ray radiation from the ultracool dwarf TVLM 513-46546. It is shown that the source of ≈1 keV X-ray emission is not the entire corona of the brown dwarf, but a population of several hundreds of coronal magnetic loops, with 10 MK plasma heated upon dissipation of the electric current generated by the photospheric convection. Unlike models, which assume a large-scale magnetic structure of the microwave source, our model suggests that the microwave radiation comes from hundreds of magnetic loops quasi-uniformly distributed over the dwarf's surface. We propose a long-term operating mechanism of acceleration of electrons generating gyrosynchrotron radio emission caused by oscillations of electric current in the magnetic loops as an equivalent RLC circuit. The second population of magnetic loops—the sources of microwave radiation, is at the same time a source of softer (≈0.2 keV) X-ray emission.

**Keywords:** brown dwarf; X-ray emission; microwave radiation; magnetic loops; particle acceleration

#### the Origin of Persistent Radio and X-ray Emission from Brown Dwarf TVLM 513-46546. *Universe* **2022** , *8* , 77. https://doi.org/10.3390/

**Citation:** Stepanov, A.; Zaitsev, V. On

Academic Editors: Galina L. Klimchitskaya, Vladimir M. Mostepanenko and Nazar R. Ikhsanov

universe8020077

Received: 6 December 2021 Accepted: 25 January 2022 Published: 27 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

The brown dwarf TVLM 513-46546 is currently one of the most studied among ultracool stars of spectral types M > 7, L, T, and Y. Radio, X-ray, ultraviolet, and optical spectroscopic observations have revealed rather unusual properties of TVLM 513-46546 radiation (see, for example, [1–3]). In addition to the surprising radio to X-ray luminosity ratio (*L*R/*L*<sup>X</sup> ≈ 3.3 × <sup>10</sup>−<sup>12</sup> Hz<sup>−</sup>1), by a factor of 104 larger than expected [2,4], a persistent (on a multiyear timescale) quiescent nonthermal radio emission from the brown dwarf TVLM 513-46546 is observed, which assumes a continuous source of accelerated electrons. The persistent quiescent radio emission is accompanied by long-term X-ray emission. By now, in low-mass dwarf stars, two thermal components of X-ray emission are observed: the "soft" one, with the temperature *<sup>T</sup>* ≈ (1–6) × 106 K, and the "hard" one, corresponding to the source temperature *<sup>T</sup>* ≈ <sup>10</sup><sup>7</sup> K [5].

The current models of the microwave [1] and multiwavelength radiation sources in TVLM 513-46546 [2] provide no explanation for the nature of the accelerator of electrons producing the quiescent nonthermal gyrosynchrotron emission for several years. Moreover, the microwave radiation source was represented either by the dipole magnetic field of the dwarf [1] or as a large-scale magnetic structure with a "covering factor" of about 50% [2].

Here, we show that the sources of radio and X-ray emission of TVLM 513-46546 consist of a few ×100 coronal magnetic loops quasi-uniformly distributed over the dwarf surface. We propose the heating mechanism of coronal loop plasma driven by dissipation of electric current generated by photospheric convection and the mechanism of long-term particle acceleration in the atmosphere of a brown dwarf supporting the persistent nature of the radio flux. Electric currents required for the heating of the loop up to *<sup>T</sup>* ≈ <sup>10</sup><sup>7</sup> <sup>K</sup> are determined. Therefore, the magnetic loops responsible for the persistent microwave emission can also be the sources of the "soft" X-ray component.

#### **2. Radio Emission from TVLM 513-46546**

The M8.5V ultracool dwarf TVM 513-46546 (hereafter, TVLM 513) has an effective temperature of the photosphere *T*eff = 2200 K, the mass *M*<sup>∗</sup> = 0.07 *M*, the radius *R*<sup>∗</sup> = 0.1 *R*, and is located at a distance of 10.6 pc from the Sun. The gravitational acceleration on the star surface is *<sup>g</sup>* = 2 × 105 cm/s2, and rotation velocity *v sin i* ≈ 60 km/s. In the case of brown dwarfs, energy is transferred from the star core to its surface by convection. On the photosphere level, the convection velocity for the stars of M ≥ 7 spectral type varies from *<sup>V</sup>* ≈ <sup>10</sup>3–104 cm/s [6] to *<sup>V</sup>* ≈ 1.4 × 105 cm/s [7]. The size of granulation cells for M8.5V stars approximately coincides with that of supergranulation cells, *<sup>D</sup>* ≈ 1.4 × 107 cm [8].

TVLM 513 displays two components of microwave radiation in the range of 4.8–8.5 GHz. The first one consists of periodic bursts with high brightness temperature ≥10<sup>11</sup> K, highly beamed emission, and 100% circular polarization [3]. The periodicity presented in burst emission (≈1.96 h) is due to the rotation of the dwarf. The peculiarity of the microwave emission was explained in terms of the Electron Cyclotron Maser [9] or the plasma radiation mechanism [10]. The second component, which is called "quiescent", displays a brightness temperature of about 109 K and a small degree of circular polarization (<15%) [1,2]. The weak modulation of the quiescent component under the dwarf rotation means that the sources are distributed uniformly in the TVLM 513 atmosphere, and their total area is comparable with the area of the dwarf's surface. In addition, the respective sources of the radio emission should be supplied continuously by energetic particles, to compensate the losses caused by particles precipitation into the loss-cone and to provide the observed brightness temperature. Microwave observations of the TVLM 513 indicate that the quiescent component is stable on a multiyear timescale [2].

Observations with Very Large Array revealed the spectral index α in quiescent radio flux *F*ν~ν<sup>α</sup> for the wavelengths 20–3.6 cm [1]. It is positive between 20 and 6 cm and equal to α = 0.1 ± 0.2. For the wavelengths 6–3.6 cm, the spectral index is negative and equal to α = −0.4 ± 0.1. A change in the sign of the spectral index implies that the spectrum of quiescent emission has a maximum and the maximum frequency νmax is in the range of νmax = 1.4–4.8 GHz (i.e., 20–6 cm) [1]. The change in the frequency index's sign in the gyrosynchrotron mechanism means that the radiation mode transits from optically thick to optically thin. In other words, the source of radiation in the range 4.8–8.4 GHz is most probably optically thin, which is important for further estimations.

We will use the information on the flux of the quiescent component at 8.4 GHz, i.e., where the emission source is optically thin. For the power-law distribution of energetic electrons *n*e(ε) ∝ ε−δ, the flux of gyrosynchrotron emission from an optically thin source at the frequency ν is [11]

$$F\_{\nu} = 3.3 \times 10^{-24} \times 10^{-0.52 \delta} (\sin \theta)^{-0.43 - 0.65 \delta} \left(\frac{\nu}{\nu\_{\varepsilon}}\right)^{1.22 - 0.9 \delta} (n\_{\varepsilon} d) B \frac{S}{R^2} \text{ ergs/cm}^2 \text{s Hz} \tag{1}$$

Here, *ϑ* is the angle between the magnetic field *B* and the direction toward the observer, *<sup>ν</sup>*<sup>e</sup> = 2.8 × <sup>10</sup><sup>6</sup> B is electron gyrofrequency, *ne* is the number density of energetic electrons, *<sup>d</sup>* is the source thickness in the projection to the observer, S is the source area, and *R* = 10.6 pc is the distance to the source. For the maximum frequency in the gyrosynchrotron spectrum and the polarization degree, we use the following formulas [11]:

$$\nu\_{\rm peak} = 2.72 \times 10^3 \times 10^{0.27\delta} (\sin \theta)^{0.41 + 0.03\delta} (n\_\ell d)^{0.32 - 0.03\delta} B^{0.68 + 0.03\delta} \tag{2}$$

$$r\_c = 1.26 \times 10^{0.035 \delta} \times 10^{-0.071 \cos \theta} \left(\frac{\nu}{\nu\_\ell}\right)^{-0.782 + 0.545 \cos \theta} \tag{3}$$

Equations (1)–(3) are true in the angle range *ϑ* = 20◦ –80◦ and for the optical thickness *τν* = *κνd* < 1, where

$$\kappa\_{\nu} = 1.4 \times 10^{-9 - 0.22\delta} (\sin \theta)^{-0.09 + 0.72\delta} \left(\frac{\nu}{\nu\_{\varepsilon}}\right)^{-1.30 - 0.98\delta} \frac{n\_{\varepsilon}}{B} \tag{4}$$

For the model of quiescent component, we determine the visible area of the source as *S* ≈ *Nloopld*, where *Nloop* is the number of elementary radiation sources (magnetic loops) on the visible hemisphere (Figure 1), and *d* and *l* are a typical thickness and length of a single loop.

**Figure 1.** Sketch of brown dwarf disc with two populations of magnetic loops, the sources of "hard" X-ray emission (red), and both microwave and "soft" X-ray radiation (blue).

Using the relation *α* = 1.22 − 0.90*δ* from the spectral index of the optically thin radio emission α = −0.4 ± 0.1, one can determine the spectral index of energetic electrons *δ* ≈ 1.8 ± 0.1, indicating a rather hard spectrum of emitting particles. Equation (3) for the degree of circular polarization (15%) was used to estimate the magnetic field in coronal magnetic loops, the sources of gyrosynchrotron radiation of TVLM 513, B ≈ 100 G [2].

It should be noted that one of the problems with multi-loop radio source is the estimation of the magnetic field using the observed degree of circular polarization, since various orientations of magnetic loops with different magnetic polarity should lead to a decrease in the Stokes V parameter. We assume here that despite this some predominant average polarity is retained in set of the loops, which gives a final degree of observed circular polarization (about 15%) that is rather low as compared to the polarization degree of gyrosynchrotron emission from solar loops (about 30%). The predominant average magnetic polarity can be associated with the large-scale magnetic field of the dwarf, which is partially fragmented into loops by the photospheric convection. A good example of the conservation of the magnetic field polarity in the loop set one can find from the observations of solar magnetic loops with the Solar and Heliospheric Observatory (SOHO) [12]. Among 30 loops recorded on 1999 August 30 in AR 7986, in 21 loops the negative magnetic polarity was retained. This example means that our estimates of the magnitude of the magnetic field from the observed polarization of microwave radiation from TVLM 513 give a lower average magnitude of the magnetic field in the loops.

The values of the magnetic field and the spectral index of energetic electrons, together with the condition of optically thin sources at 8.4 GHz, impose a limitation on the value *ned* ≤ <sup>2</sup> × 1016 cm<sup>−</sup>2, which determines the microwave flux. Assuming the loop thickness of the order of the granulation cell, *<sup>d</sup>* ≈ <sup>10</sup><sup>7</sup> cm, *<sup>l</sup>* ≈ (2.5–5) × 109 cm, we obtain the area of all elementary sources *<sup>S</sup>* ≈ (2.5–5) × 1016 *<sup>N</sup>*loop cm2; from the microwave flux (Equation (1))

$$F\_{\nu=8.4} = 228 \mu f y \approx (2.28 \text{--} 5.6) \times 10^{-46} (n\_{\ell} d) N\_{loop} \text{ergs cm}^{-2} \text{s}^{-1} \text{Hz}^{-1} \tag{5}$$

one can estimate the number of sources of quiescent radio emission in the visible hemisphere required to explain the observed microwave flux, *Nloop* ≈ (0.4–1.0) × <sup>10</sup>19/(*ned*) ≈ 200–500. The filling factor of the stellar surface with the sources is about 3–16%. The number density of energetic particles is found from the condition of an optically thin source (Equations (1) and (4)) at 8.4 GHz, *ned* ≤ <sup>2</sup> × <sup>10</sup><sup>16</sup> cm<sup>−</sup>2, which yields *ne* ≤ <sup>2</sup> × <sup>10</sup><sup>9</sup> cm<sup>−</sup>3.

#### **3. Long-Term X-ray Emission of TVLM 513 and the Coronal Heating Problem**

The information collected by Chandra ACIS-S3 shows that the quiescent radio emission from the brown dwarf TVLM 513 is accompanied by soft X-ray radiation in the energy range 0.3−2 keV [2]. To estimate the plasma temperature from 0.3–2 keV spectrum, the authors of [2] have shown that the best-fit temperature of TVLM 513 is 0.9 keV, or *<sup>T</sup>* ≈ <sup>10</sup><sup>7</sup> K. The X-ray flux in the 0.3–2 keV range is *<sup>F</sup>*<sup>x</sup> ≈ 6.3×10−<sup>16</sup> ergs/cm2 s and the total luminosity *L*<sup>x</sup> = 4π*R*<sup>2</sup> <sup>∗</sup>*F*<sup>x</sup> <sup>≈</sup> 8.4 <sup>×</sup> <sup>10</sup><sup>24</sup> erg/s. This value of luminosity corresponds to the emission measure *EM*<sup>∗</sup> = *L*x/*P*(*T*), where [5]

$$P(T) = 2 \times 10^{-27} \sqrt{T} + 5 \times 10^{-25} \exp\sqrt{2.8 + 10^6 \text{K}/T} \text{ ergs cm}^3 \text{s}^{-1}. \tag{6}$$

putting *<sup>T</sup>* = 107 K in Equation (6), we obtain *EM*<sup>∗</sup> <sup>≈</sup> <sup>9</sup> <sup>×</sup> <sup>10</sup><sup>47</sup> cm<sup>−</sup>3.

There is the problem of heating of the entire corona of TVLM 513. Indeed, if the source of the quasi-stationary X-ray radiation is a corona with a temperature of *<sup>T</sup>* ≈ <sup>10</sup><sup>7</sup> K, then the emitting volume is V ≈ <sup>4</sup>π*R*<sup>2</sup> <sup>∗</sup>*H*, where *<sup>H</sup>* <sup>=</sup> *kBTR*<sup>2</sup> <sup>∗</sup>/*miGM*∗, *kB* is the Boltzmann constant, *<sup>G</sup>* = 6.67 × <sup>10</sup><sup>−</sup>8cm3/gs2 is the gravitational constant. The X-ray luminosity from corona is

$$L\_{\mathfrak{X}} = P(T)\overline{n}^2 \mathcal{V}\_{\prime} \tag{7}$$

where *P*(*T*) is determined by Equation (6), *n* is the average number density of electrons (ions) in the volume V. Using Equation (7), we determine the volume of the emitting corona V and the average values of the number density of particles in the corona from the observed values of the radiation luminosity and coronal temperature. For TVLM-513, we find V ≈ 2.5 × 1030 cm3, the height scale *<sup>H</sup>* ≈ <sup>4</sup> × <sup>10</sup><sup>9</sup> cm, and the average plasma density in the corona *<sup>n</sup>* ≈ <sup>6</sup> × <sup>10</sup><sup>8</sup> cm<sup>−</sup>3.

As a result of hot coronae, brown dwarfs noticeably increase their energy losses due to thermal conductivity. The radiation losses for TVLM-513 are *Qr* ≈ <sup>3</sup> × 1025ergs/s [13], while the energy losses due to thermal conductivity [14,15] *QT* <sup>≈</sup> 0.9 <sup>×</sup> <sup>10</sup>−6*<sup>T</sup>* <sup>7</sup> <sup>2</sup> *VH*−<sup>2</sup> ≈ <sup>4</sup> <sup>×</sup> <sup>10</sup>29ergs/s significantly exceed similar losses in the solar corona, *QT* <sup>≈</sup> <sup>10</sup>28ergs/s. Therefore, to maintain the corona of TVLM 513 with *<sup>T</sup>* ≈ <sup>10</sup><sup>7</sup> K, more powerful heating sources are required than those in the Sun. Let us assume that in the atmospheres of brown dwarfs, as is the case in the solar corona, there are hot Type II spicules, which in the last decade have been considered to be a probable source of solar corona heating [15–17]. Then the heat flux from the single spicule into the corona due to electron thermal conductivity along the magnetic field of the spicule is [15]

$$Q\_{\rm Tsp} \approx \frac{0.9 \times 10^{-6} T^{\frac{7}{2}}}{\Delta z} \pi r\_0^2 \text{ ergs/s.} \tag{8}$$

For the radius of the spicule (magnetic flux tube) of brown dwarfs, the size of the granulation cell can be taken to be *<sup>r</sup>*<sup>0</sup> ≈ 107 cm. Then for *<sup>T</sup>* = 10<sup>7</sup> K and a spicule length <sup>Δ</sup>*<sup>z</sup>* = 10<sup>9</sup> cm, from Equation (8) we obtain *QTsp* ≈ 1024 ergs/s. Therefore, to compensate for the energy losses from the TVLM 513 corona, ~4 × <sup>10</sup><sup>5</sup> hot spicules covered persistently the dwarf surface are required. Moreover, the total foot-points area of spicules should be of the order of the dwarf surface area (1.5 × <sup>10</sup><sup>20</sup> cm2), which significantly differs from the Sun, where spicules occupy ≈1% of the surface.

#### **4. Generation of Electric Current and Plasma Heating in a Magnetic Loop: The Role of Photospheric Convection**

Energy is transferred from the TVLM 513 core to the surface by the convection. At the photosphere level, the convection velocity is V ≈ <sup>10</sup>4–105 cm/s [6,7]. The half-thickness of the magnetic loops is of the order of the size of a granulation cell *<sup>r</sup>*<sup>1</sup> ≈ *<sup>D</sup>* ≈ 107 cm [8]. When the radial component of the convection velocity *V*<sup>r</sup> interacts with the azimuthal component of the magnetic field of the loop *Bϕ*, an EMF arises, which generates an electric current flowing from one foot point of the loop through the coronal part to the other, and closing in the photosphere, where the inequality *ωe*/*ωiνeaνia* 1 is satisfied and the conductivity becomes isotropic. Here, *ω<sup>e</sup>* and *ω<sup>i</sup>* are gyrofrequencies of electrons and ions, *νea* and *νia* are electron-atom and ion-atom collision frequencies. Thus, the coronal magnetic loop with the photospheric current channel is an equivalent RLC-circuit [18–20]. In the self-consistent equation of the equivalent electric circuit, the resistance *R* and capacitance *C* are found to be dependent on the electric current [21]:

$$\frac{1}{c^2}L\frac{\partial^2 y}{\partial t^2} + \left[R(I) - \frac{|\overline{V}\_r|l\_1}{c^2 r\_1}\right]\frac{\partial y}{\partial t} + \frac{1}{C(I)}y = 0,\tag{9}$$

where *y*(*t*) = [Δ*I*(*t*) − *I*]/*I*, *L* is a loop inductance. EMF is in the photospheric foot point of the loop and is equal to

$$\Xi = \frac{l\_1}{\pi c r\_1^2} \int\_0^{r\_1} V\_r B\_\theta 2\pi r dr \approx \frac{\left|\overline{V\_r}\right| I l\_1}{c^2 r\_1} \tag{10}$$

where *l*<sup>1</sup> is the height of the RLC-circuit section in the area of EMF action, *I* is the longitudinal electric current, *B<sup>ϕ</sup>* ≈ 2*I*/*cr*1. Estimates show that the main contribution to the loop resistance is made by the region of photospheric EMF, where the Cowling conductivity due to ion-atom collisions plays a decisive role. The resistance value in this case [22]

$$R(I) \approx \frac{1.5 l\_1 I^2 F^2}{\pi r\_1^4 c^4 n m\_i \nu\_{\rm ia}'(2 - F)} \; , \tag{11}$$

where *F* = *na*/(*n* + *na*) is the relative density of neutrals, *n* is the density of electrons, *ν ia* ≈ 1.6 × <sup>10</sup>−11*na* <sup>√</sup>*<sup>T</sup>* Hz is the effective frequency of collisions of ions with neutrals. The steady-state value of the current in the loop is determined from the condition *R*(*I*) = Ξ(*I*), hence at *l*<sup>1</sup> ≈ *r*<sup>1</sup> we obtain

$$I \approx 0.8[ |V\_r| \pi r\_1^3 c^2 \nu m\_i \nu\_{ia}' (2 - F) F^{-2} ]^{1/2} \text{cgs} \tag{12}$$

For the brown dwarf considered, in the height interval Δ*z*~*l*1, the atomic density *na* ≈ <sup>5</sup> × <sup>10</sup><sup>16</sup> − <sup>10</sup><sup>17</sup> cm<sup>−</sup>3, *<sup>n</sup>* ≈ <sup>5</sup> × 109 − 1010cm−3, *<sup>T</sup>* ≈ 2200 K. In this case, for the rate of photospheric convection *Vr* ≈ 104 − 105 cm/s and *<sup>r</sup>*<sup>1</sup> ≈ <sup>10</sup>7cm, we obtain the estimates of the current value: *<sup>I</sup>* ≈ 2.5 × 109 ÷ <sup>7</sup> × 1010 A. At the coronal level for an optically thin medium at temperatures *<sup>T</sup>* > <sup>2</sup> × <sup>10</sup><sup>5</sup> K, hydrogen makes the main contribution to the neutral component, while the relative content of neutrals is determined by the formula [23]

$$F = 0.32 \times 10^{-3} \frac{1 + \frac{T}{6T\_H}}{\left[\sqrt{\frac{T}{T\_1}}\right]^{2-b} \sqrt{T} \left[1 + \sqrt{\frac{T}{T\_1}}\right]^{1+b}} \times p \frac{T\_H}{T} \tag{13}$$

where *TH* = 1.58 × 105 K, *<sup>T</sup>*<sup>1</sup> = 7.036 × 105 K, *<sup>b</sup>* = 0.748. In the temperature range of interest, 106–107 K, the approximation *<sup>F</sup>* ≈ 0.15/*<sup>T</sup>* is valid. The rate of Joule dissipation of the current per unit volume of the loop, taking into account Equation (11), is determined by the formula [22]

$$q\_I = \frac{j\_z^2}{\sigma} + \frac{F^2 B\_{\phi}^2 j\_z^2}{(2 - F)c^2 nm\_i \nu\_{ia}^{\prime}} \approx 2.2 \times 10^{-9} \frac{I^4}{n^2 r\_1^6 T^{2/3}} \text{ ergs/cm}^3 \text{s} \tag{14}$$

In the coronal part of the loop, the main contribution to the current dissipation is made by the second term in Equation (14). Despite the relatively low density of neutral atoms, *F* 1, an effective dissipation channel associated with Cowling conductivity is switched on. This is due to a decrease in the effective plasma conductivity

$$
\sigma\_{eff} = \frac{\sigma}{1 + \frac{F^2 \omega\_c \omega\_i}{(2 - F)\nu\_e' \nu\_{ia}'}} \tag{15}
$$

because in the corona the second term in the denominator of Equation (15) is 1. The heating of the magnetic loop is determined by the balance of Joule dissipation, thermal conductivity, and radiation losses:

$$\frac{d}{ds}\kappa\_e T^{5/2} \frac{dT}{ds} = q\_r - q\_I. \tag{16}$$

Here, *<sup>κ</sup><sup>e</sup>* <sup>=</sup> 0.9 <sup>×</sup> <sup>10</sup>−6, *qr* <sup>=</sup> *<sup>χ</sup>*<sup>0</sup> *<sup>p</sup>*<sup>2</sup> 4*k*<sup>2</sup> *B <sup>T</sup>*<sup>−</sup>5/2, *qJ* <sup>=</sup> <sup>10</sup>−8*k*<sup>2</sup> *B I*4 *p*2*r*<sup>6</sup> 1 *T*1/2, *s* is the coordinate along the loop, *χ*<sup>0</sup> = 10−19, *p* = 2*nkBT*. Equation (15) is solved under the following boundary conditions at the foot point (*s* = 0) and at the loop top (*s* = *l*/2): *T* = *T*0, *dT*/*ds* = 0 at *s* = 0 and *T* = *T*1, *dT*/*ds* = 0 at *s* = *l*/2. It is assumed that *T*<sup>1</sup> *T*0. From Equation (16), we find the distribution of temperature and density along the loop [22]

$$T\_1 \approx 2 \times 10^{-2} \frac{\left(IL\right)^{4/9}}{r\_1^{2/3}} \text{ K, } n\_1 = \frac{1}{3} \left(\frac{2\kappa\_\varepsilon}{\chi\_0}\right)^{1/2} \frac{T\_1^2}{L} \text{ cm}^{-3}, \text{ } L = \text{l/2} \tag{17}$$

It can be seen from Equation (17) that the temperature at the top of the loop reaches its maximum, and the plasma density its minimum. Both values vary slightly over most of the length of the loop, and at the foot points they reach the equilibrium values. Let us determine the number of hot loops in the corona of TVLM 513, which provide the observed measure of X-ray emission *EM*<sup>∗</sup> <sup>≈</sup> <sup>9</sup> <sup>×</sup> 1047 cm−3. By analogy with the solar corona, we assume that the thickness of the loops does not change with height, i.e., *<sup>r</sup>*<sup>1</sup> ≈ <sup>10</sup><sup>7</sup> cm. Then from Equation (17) we find the electric current required to heat the loop plasma, *<sup>I</sup>* ≈ <sup>6</sup> × 1010 A, the plasma density *<sup>n</sup>* ≈ 4.8 × 1010 cm−<sup>3</sup> and the measure of the emission of an individual loop *EMloop* ≈ 4.6 × 1045 cm<sup>−</sup>3. To ensure the observed emission measure of "hard" X-ray radiation, it is necessary that *Nloop* = *EM*∗/*EMloop* ≈ 200 hot loops should be presented in the TVLM 513 corona. As shown in Section 2, to explain the persistent microwave emission from TVLM 513 as well as the "soft" X-ray radiation, about 200–500 loops are needed (Figure 1).

#### **5. Mechanism of Long-Term Acceleration of Electrons in a Magnetic Loop—An Equivalent Electric Circuit**

The persistent microwave emission generated in the TVLM 513 corona by a set of magnetic loops suggests that magnetic loops must be constantly replenished with energetic particles to compensate for the losses associated with the escape of particles into the losscone. In [24,25] a mechanism of long-term acceleration of electrons, caused by oscillations of the electric current in the loop, generated by photospheric convection, was proposed and developed. Let us show that even with a relatively small electric current, the acceleration of electrons by current oscillations in the loop is effective. From Equation (9) it follows that the current oscillations of the loop as an equivalent electric circuit are excited when the negative resistance of the photospheric EMF exceeds the loop resistance, *R*(*I*) ≤ |*Vr*|*l*1/ *r*1*c*<sup>2</sup> . The

eigen-frequency of the equivalent RLC-circuit depends on the loop radius *r*1, density *n*, and length *l* of the coronal part of the loop [21,22]:

$$\upsilon\_{\rm RLC} = \frac{c}{2\pi\sqrt{LC(l\_0)}} \approx \frac{1}{(2\pi)^{3/2}\sqrt{\Lambda}} \frac{I\_0}{cr\_1^2\sqrt{nm\_i}}, \ \Lambda = \ln\left(\frac{4l}{\pi r\_1}\right) - \frac{7}{4} \tag{18}$$

Assuming *<sup>r</sup>*<sup>1</sup> ≈ <sup>10</sup>7cm, *<sup>n</sup>* ≈ 1010 cm<sup>−</sup>3, *<sup>I</sup>*<sup>0</sup> ≈ 108 A, *<sup>l</sup>* ≈ <sup>6</sup> × <sup>10</sup><sup>9</sup> cm, from Equation (18) we obtain an estimate for the eigen-frequency of the electric circuit: *<sup>ν</sup>RLC* ≈ <sup>2</sup> × <sup>10</sup>−<sup>2</sup> Hz (that is, a period of 50 s). Oscillations of the electric current are associated with those of the azimuthal component of the magnetic field in the magnetic loop, *Bϕ*(*r*, *t*) = 2*r Iz*(*t*)/*cr*<sup>2</sup> <sup>1</sup>. These oscillations, in turn, according to the equation *rot*<sup>→</sup> *E* = −(1/*c*)*∂* → *Bϕ*/*∂t*, lead to the generation of an electric field directed along the loop axis. Assuming *I*(*t*) = *I*<sup>0</sup> + Δ*Isin*(2*πνRLCt*), we obtain the electric field averaged over the loop cross section [24]

$$
\overline{E}\_z = \frac{4\nu\_{RLC}I\_0}{3c^2} \frac{\Delta I}{I\_0}.\tag{19}
$$

Then from Equation (19) we obtain at *<sup>I</sup>*<sup>0</sup> = <sup>10</sup><sup>8</sup> <sup>A</sup> = <sup>3</sup> × 1017 cgs, <sup>Δ</sup>*I*/*I*<sup>0</sup> = <sup>10</sup>−<sup>2</sup> the value of the electric field: *Ez* ≈ 8.9 × <sup>10</sup>−<sup>8</sup> cgs ≈ 2.67 × <sup>10</sup>−<sup>5</sup> V/cm. In such an electric field at a distance of *<sup>l</sup>* = 6 × <sup>10</sup><sup>9</sup> cm, electrons can acquire energy *<sup>ε</sup>* ≈ 160 keV. At <sup>Δ</sup>*I*/*I*<sup>0</sup> = 0.02 more energetic electrons, *ε*≈ 320 keV, will be accelerated. For the plasma density *n* = 10<sup>10</sup> cm−<sup>3</sup> and the temperature *T* = 10<sup>6</sup> K, the Dreicer field is [26]

$$E\_D = 6 \times 10^{-8} (n/T) \text{ V/cm} \approx 6 \times 10^{-4} \text{ V/cm} \tag{20}$$

and the ratio of the Dreicer field to the accelerating field is *x* = *ED*/*Ez*. The kinetic theory gives the rate of electron acceleration per unit volume [27]:

$$
\dot{m}\_{\ell} = 0.35 n \nu\_{\rm el} \mathbf{x}^{3/8} \exp\left(-\sqrt{2\mathbf{x}} - \mathbf{x}/4\right) \tag{21}
$$

where *<sup>ν</sup>ei*≈ 60n*T*−3/2 <sup>s</sup>−1. Substituting the values of *<sup>n</sup>* and *<sup>T</sup>* into Equation (21), we find the acceleration rate . *<sup>n</sup>*<sup>e</sup> ≈ <sup>3</sup> × 107 el/cm3s. The condition for the generation of microwave radiation (Section 2) requires *ned* < 2 × 1016 cm−2, which for *<sup>d</sup>* = 2*r*<sup>1</sup> = 2 × 107 cm gives *n*<sup>e</sup> < 10<sup>9</sup> cm−3. The developing small-scale turbulence (whistlers for example) prevents the free escape of particles from the magnetic loop and contributes to the accumulation of accelerated particles in coronal magnetic loops.

The regimes of pitch-angle diffusion of accelerated particles into the loss-cone depend on the ratio of the following time scales [28]: the diffusion time *tD*, which corresponds to the mean time for the change in the pitch-angle of a fast electron at π/2 due to the wave-particle interaction, the time of filling the loss-cone *tD*/*σ*, and the time of flight when escaping into the loss-cone, *t*<sup>0</sup> = *l*/*v* ≈ 0.3 s. Here, *σ* = *B*foot/*B*top is the mirror ratio of a magnetic loop. In the case of intermediate diffusion, when the condition *t*<sup>0</sup> < *tD* < *σt*<sup>0</sup> is satisfied, the loss-cone is filled faster than the particles escape into the loss-cone. In this case, the lifetime of a particle in a magnetic trap is of the order of *σt*<sup>0</sup> and the number density of accelerated electrons is [29]

$$m\_{\varepsilon} \approx \frac{\dot{n}\_{\varepsilon}(\mathbf{x}) \sigma l}{\sqrt{2\varepsilon/m}}\tag{22}$$

Substituting into Equation (22) . *ne* = 3 × 107 el/cm3s, *<sup>σ</sup>* = 3, *<sup>l</sup>* = 6 × <sup>10</sup><sup>9</sup> cm, and *<sup>ε</sup>* = 160 keV, we obtain *ne* ≈ 2.3 × 107 cm−3. In the case of strong diffusion, when the condition *tD* < *t*<sup>0</sup> is satisfied, the density of accelerated electrons *ne* is by about an order of magnitude higher. Therefore, the proposed acceleration mechanism will provide the persistent microwave emission.

#### **6. Discussion**

In the problem of the origin of X-ray radiation and heating of stellar coronae, the question arises: what is emitting? Is it the entire hot corona or its local regions, such as coronal magnetic loops, and what are the constant sources of the heating for the corona or the local regions? We have shown that plasma heating to the temperature T~10 MK provided by electric current in several hundred magnetic loops quasi-uniformly distributed over the surface of the brown dwarf TVLM 513, which explains the X-ray emission luminosity of ~1 keV, is energetically more favorable than heating the entire corona. The second population of colder, with T~1 MK, coronal loops is a source of persistent nonthermal microwave radiation explained by gyrosynchrotron radiation of ~150–300 keV electrons in a magnetic field of ~100 G. The proposed mechanism of long-term acceleration of electrons by induced electric fields is based on small oscillations of electric current in the coronal loops supported by photospheric convection.

In the model proposed by Osten et al. [1], the source of microwave radiation is represented by the dipole magnetic field of the dwarf. Berger et al. [2] suggested that the source of microwave radiation is a large-scale magnetic structure with a "covering factor" of about 50%. Both models cannot provide a long-term supply of sub-relativistic electrons to the radiation sources. In addition, the Berger's model does not identify the X-ray source in the brown dwarf corona.

It follows from Equation (17) that the plasma heating temperature is proportional to *T*~*I* 4/9. For a current *I* = 10<sup>8</sup> A, at which the acceleration of electrons is sufficiently efficient, the plasma in the loop heats up only to a temperature of 4 × <sup>10</sup><sup>5</sup> K, and at *<sup>I</sup>* = 10<sup>9</sup> <sup>A</sup> up to 10<sup>6</sup> K, i.e., the "microwave" loops are also sources of a softer X-ray component with a maximum in the spectrum of ≈0.2 keV. With an increase in the current value, the lumped-loop approach used by us is invalid, since the period of the RLC oscillations becomes shorter than the propagation time of the Alfvén disturbance along the loop *t*<sup>A</sup> = *l*/*V*<sup>A</sup> ≈ 30 s. During this time, the current oscillations change direction many times and the acceleration of electrons is ineffective. This circumstance can explain the absence of noticeable nonthermal radio emission from hot coronal loops with the plasma temperature *T* = 10 MK. In this work, we did not consider the origin of the intense, with the brightness temperature *<sup>T</sup>*<sup>b</sup> ≥ 1011 K, flare component of microwave radiation from TVLM 513, the nature of which was discussed in [9,10]

#### **7. Conclusions**

Photospheric convection, interacting with the magnetic field at the foot points of the magnetic loops, leads to generation of an electric current, *<sup>I</sup>* ≈ <sup>10</sup>9–7 × <sup>10</sup><sup>10</sup> A, the dissipation of which under the conditions of partially ionized plasma causes the loop plasma heating up to a temperature of *<sup>T</sup>* ≈ 106–107 K. Hot plasma of ~200 loops, quasiuniformly distributed over the dwarf's surface, is a source of ~1 keV X-ray emission from TVLM 513.

Photospheric convection maintains continuous oscillations of electric current in the magnetic loops as an equivalent RLC circuit, generating the induced electric fields that accelerate electrons up to 150–300 keV in magnetic loops distributed over the dwarf's surface. This provides the persistent quiescent microwave radiation of a few ×100 loops of TVLM 513 by the gyrosynchrotron mechanism.

Effective acceleration of electrons by the electric current oscillations is possible even at relatively low values of electric current, *<sup>I</sup>* ≥ <sup>10</sup><sup>8</sup> A. At currents *<sup>I</sup>* ≈ <sup>3</sup> × <sup>10</sup>8–109 A, the loop plasma heats up only to temperature *<sup>T</sup>* ≈ <sup>10</sup><sup>6</sup> K, i.e., magnetic loops emitting microwave radiation can simultaneously be the sources of ~0.2 keV X-ray emission.

Thus, photospheric convection not only forms numerous magnetic loops in the corona of TVLM 513, but also generates an electric current leading to heating of the loop plasma and acceleration of electrons. As a result, two populations of loops quasi-distributed over the TVLM 513 surface are formed, one of which is a source of "hard" X-ray radiation, and the other is simultaneously a source of microwave and "soft" X-ray radiation.

**Author Contributions:** Investigation, A.S. and V.Z.; writing, A.S. and V.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by the Russian Science Foundation, grant No. 20-12-00268, the RFBR, grant No. 20-02-00108, and the State Assignment Nos. 0030-2021-0002, 1021032422589-5, and 0040-2019-0025.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All necessary data are contained in this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Dmitry Yakovlev**

Ioffe Institute, Politekhnicheskaya 26, 194021 St. Petersburg, Russia; yak.astro@mail.ioffe.ru

**Abstract:** A simple and well known model for thermal radiation spectra from a magnetized neutron star is further studied. The model assumes that the star is internally isothermal and possesses a dipole magnetic field (*B* - 10<sup>14</sup> G) in the outer heat-insulating layer. The heat transport through this layer makes the surface temperature distribution anisotropic; any local surface element is assumed to emit a blackbody (BB) radiation with a local effective temperature. It is shown that this thermal emission is nearly independent of the chemical composition of insulating envelope (at the same taken averaged effective surface temperature). Adding a slight extra heating of magnetic poles allows one to be qualitatively consistent with observations of some isolated neutron stars.

**Keywords:** neutron stars; radiation transfer; magnetic fields

#### **1. Yury N. Gnedin**

This paper is dedicated to the memory of Yury N. Gnedin (1935–2018) Figure 1. He was my senior colleague at Theoretical Astrophysics Department of Ioffe Insitute from 1971 to 1984, and we often met later, when he moved to Pulkovo Observatory. He was always interested in many scientific fields. For instance, his pet subjects were radiation transfer, neutron stars and magnetic fields (although he never lost scientific interest in a great amount of other things). When neutron stars were discovered in 1967, he was extremely enthusiastic about them and initiated ultrafast spreading of interest to these objects among colleagues. He was the first author of the first publication [1] on neutron stars at the Ioffe Institute after the discovery of these stars. He wrote—with different co-authors—plenty of (now) classical papers on neutron star physics. In particular, they include a basic formulation of the radiation transfer problem in a strongly magnetized plasma [2], a theoretical prediction of electron cyclotron lines in the spectra of radiation from magnetized neutron stars [3], the first realistic studies of ionization in the magnetized hydrogen atmospheres of neutron stars [4]. He was also greatly fond of history.

Although we did not work much together, I, as many others, respected Yury Gnedin very much. When he could help, he really did, and he was a Real Gentleman—friendly and open to any one, from a prominent scientist to a hopeless student. I will always remember the wonderful atmosphere, created by Yury Gnedin, and warm hospitality of his family and home. It is amazing that both his sons have grown up into excellent scientists.

I present a small piece of neutron star physics on thermal radiation from magnetized neutron stars. Hopefully, he would like it.

**Citation:** Yakovlev, D. A Simple Model of Radiation from a Magnetized Neutron Star: Accreted Matter and Polar Hotspots. *Universe* **2021**, *7*, 395. https://doi.org/ 10.3390/universe7110395

Academic Editors: Nazar R. Ikhsanov, Galina L. Klimchitskaya and Vladimir M. Mostepanenko

Received: 30 September 2021 Accepted: 18 October 2021 Published: 21 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** Yury N. Gnedin. (**Left**): A rare moment of leisure on his birthday, 13 August 2007, near his summer cottage. (**Right**): At a conference (June 2004). Courtesy to Oleg Gnedin.

#### **2. Introduction**

I will mostly consider thermal X-ray radiation from the surface of a middle-aged (*t* - 106 yr) isolated magnetized neutron star. The radiation is supposed to emerge from the warm neutron star interior and is emitted from a thin atmosphere or condensed surface (e.g., [5]), where the radiation spectrum is formed. The interiors of the star are thought to be nearly isothermal because of the high thermal conductivity of dense matter. Nevertheless, the interiors are thermally insulated from the surface by a thin heat blanketing envelope, where the thermal conduction is poorer and more sensitive to the chemical composition and magnetic field. The field redirects the emerging heat flow and results in anisotropic surface temperature distribution. The thermal surface emission can also become anisotropic.

The anisotropy of observed surface radiation is used to infer magnetic field strength and geometry, the composition of the surface layers, global parameters of the star (such as mass and radius), as well as important parameters of superdense matter in stellar interiors (see, e.g., [5,6] and references therein).

Here, I study a simple model of thermal emission from magnetized neutron stars with isothermal interiors, assuming the star emits a blackbody (BB) radiation with a local effective surface temperature *T*s, which varies over the stellar surface under the effects of dipole surface magnetic field. This model has been invented long ago by Greenstein and Hartke [7], elaborated in the literature (e.g., [8–13]) and reviewed in [6,14]. Here I point out some properties of the surface emission, which, to the best of my knowledge, have not been studied in the literature.

Logically, this paper continues the previous one [15] (hereafter Paper I), which shows that the simple thermal spectra of magnetized neutron stars can be accurately approximated by two-BB (2BB) models. Section 3 outlines theoretical formalism. Section 4 is devoted to the effects of chemical composition of the blanketing envelopes and Section 5 extends the solution to the case, in which magnetic poles contain additional hot spots. The conclusions are formulated in Section 6.

It is important to mention more complicated neutron star emission models. They include sophisticated magnetic field effects on thermal emission, leading to specific spectral, polarization and angular properties of radiation (see, e.g., [6,16–18]), which are not reproduced by the given model.

#### **3. Simple Model Spectra**

#### *3.1. General Formalism*

The calculation of radiation spectral flux from a spherical neutron star under formulated assumptions is simple. Here is the summary using the notations of Paper I. Any small surface temperature element emits like a BB, and the flux is obtained by intergating over a visible part of the surface, taking into account that emitted quanta propagate in Schwarzschild space-time and demonstrate the gravitational redshift of photon energies and light bending. The effects of General Relativity are specified by the compactness parameter *x*<sup>g</sup> = *r*g/*R*, where *r*<sup>g</sup> = 2*GM*/*c*<sup>2</sup> is the Schwartschield radius, *M* is the gravitational star's mass, *R* is the circumferential radius, *G* is the gravitational constant and *c* is the velocity of light. One deals either with local quantities at the stellar surface (e.g., the local surface temperature *T*s) or with the quantities detected by a distant observer. The letter ones will often be marked by the symbol '∞'. For instance, *T*<sup>∞</sup> <sup>s</sup> = *T*<sup>s</sup> 1 − *x*<sup>g</sup> is the redshifted surface temperature. As an exclusion, we will denote local (non-redshifted) photon energy by *E*0, and the redshifted energy by *E* ≡ *E*<sup>∞</sup> = *E*<sup>0</sup> 1 − *x*g.

Let *F*<sup>∞</sup> *<sup>E</sup>* be a radiative spectral flux density [erg cm−<sup>2</sup> <sup>s</sup>−<sup>1</sup> keV<sup>−</sup>1] detected at a distance *<sup>D</sup> <sup>R</sup>*. It is customary to express *<sup>F</sup>*<sup>∞</sup> *<sup>E</sup>* as:

$$F\_E^{\infty} = \frac{R^2}{D^2} H\_E^{\infty} \tag{1}$$

where *H*<sup>∞</sup> *<sup>E</sup>* is the effective flux, that is formally independent of *D*. It can be calculated as an integral of the radiating surface flux over the visible part of the surface,

$$H\_{\rm E}^{\infty} = \frac{15\sigma\_{\rm SB}}{16\pi^{5}k\_{\rm B}^{4}} \int\_{\rm viz} d\Omega\_{\rm s} \frac{(1-\chi\_{\rm g})^{-1}\mathcal{P}E^{3}}{\exp(E/k\_{\rm B}T\_{\rm s}^{\infty}) - 1},\tag{2}$$

where *σ*SB is the Stefan–Boltzmann constant, *k*<sup>B</sup> is the Boltzmann constant, dΩ<sup>s</sup> is a surface solid angle element, and P is the light bending function (e.g., [5,19–21]).

It is convenient to integrate over the star's surface and calculate the bolometric luminosity of the star, *L*s, as well as the averaged non-redshifted effective surface temperature *T*eff (e.g., [10]),

$$L\_{\rm s} = \sigma\_{\rm SB} R^2 \int\_{4\pi} \, \text{d}\Omega\_{\rm s} \, T\_{\rm s}^4 \equiv 4\pi \sigma\_{\rm SB} R^2 T\_{\rm eff}^4. \tag{3}$$

For a uniform surface temperature (*T*<sup>s</sup> = *T*eff), one immediately gets the standard BB flux,

$$H\_E^{\rm BB\infty} = \frac{15\sigma\_{\rm SB}}{4\pi^4 k\_{\rm B}^4} \frac{(1-\chi\_{\rm g})^{-1} E^3}{\exp(E/k\_{\rm B}T\_{\rm eff}^{\infty}) - 1} \tag{4}$$

and the bolometric effective flux,

$$H\_{\rm bol}^{\rm BB\infty} = \int\_0^\infty H\_E^{\rm BB\infty} \,\mathrm{d}E = \frac{\sigma\_{\rm SB} T\_{\rm eff}^{\infty 4}}{1 - \chi\_{\rm g}}.\tag{5}$$

#### *3.2. Input Parameters*

Equation (2) allows one to compute thermal spectra for any given temperature distribution *T*<sup>s</sup> over the neutron star surface. We focus on the distribution created by a dipole magnetic field (with the field strength *B*pole at magnetic poles) due to anisotropic heat transport in a thin (maximum a few hundreds meters) heat blanketing envelope. The input parameters are *M*, *R*, chemical composition of the blanketing envelope, and a nonredshifted temperature *<sup>T</sup>*<sup>b</sup> at its bottom (at density *<sup>ρ</sup>*b∼10<sup>10</sup> g cm<sup>−</sup>3; see, for example, [14]). The local *T*<sup>s</sup> is usually determined by solving local quasistationary radial heat transport within the envelope mediated by an effective radial thermal conductivity. For studying the thermal surface emission, it is profitable to use *T*eff [see Equation (3)] instead of *T*b.

#### *3.3. 2BB Representation*

According to Paper I, the spectral fluxes *H*<sup>∞</sup> *<sup>E</sup>* , computed from Equation (2) for iron heat blankets, are accurately fitted by a familiar 2BB model,

$$H\_E^{\infty} = s\_\text{c} H\_E^{\text{BB\infty}}(T\_{\text{effc}}) + s\_\text{h} H\_E^{\text{BB\infty}}(T\_{\text{effh}}).\tag{6}$$

Here, *<sup>H</sup>*BB<sup>∞</sup> *<sup>E</sup>* (*T*eff) is given by Equation (4); 'c' and 'h' refer, respectively, to colder and hotter BB components. Any fit contains four parameters, which are two effective temperatures *T*effc and *T*effh, and two fractions of effective radiating surface areas, *s*<sup>c</sup> and *s*h. Instead of *T*effc and *T*effh, it is convenient to introduce two dimensionless parameters *p*<sup>c</sup> = *T*effc/*T*eff and *p*<sup>h</sup> = *T*effh/*T*eff, with *p*<sup>c</sup> < *p*h. In Paper I the fits have been done for a number of representative values of *M*, *R*, log *T*eff[K] (from 5.5 to 6.8), log *B*pole[G] (from 11 to 14), photon energies (0.064 < *E* - 40 keV, removing those *E* at which the fluxes are negligibly small) and inclination angles *i* (between line of sight and the magnetic axis). Typical relative fit errors have not exceeded a few percent, meaning the fits were very good, providing excellent analytic representation of the original computed data. Therefore, thermal X-ray spectral fluxes of magnetized neutron stars within the given model are nearly identical to those of 2BB spectral models.

Moreover, as shown in Paper I, it is sufficient to calculate the fluxes *H*<sup>∞</sup> *<sup>E</sup>* for two observation directions which are (i) pole observation, *i* = 0, to be denoted as *H*<sup>∞</sup> *<sup>E</sup>* , and (ii) equator observations, *i* = 90◦, to be denoted as *H*⊥<sup>∞</sup> *<sup>E</sup>* . If these fluxes are known, the radiation flux in any direction *i* is accurately approximated as:

$$H\_E^{i\infty} = H\_E^{\parallel \infty} \cos^2 i + H\_E^{\perp \infty} \sin^2 i. \tag{7}$$

We will see (Section 5) that this approximation is more restrictive than Equation (6).

This accurate representation of numerically calculated *Hi*<sup>∞</sup> *<sup>E</sup>* by a map of four fit parameters (*p*c, *p*h, *s*c, *s*h) in the space of input parameters has to be taken with the grain of salt. The problem is that the fluxes *Hi*<sup>∞</sup> *<sup>E</sup>* are close to 1BB fluxes (4), which leads to some degeneracy of fit parameters (with respect to a fit procedure and a choice of grid points). On the other hand, the assumed surface temperature distribution model and the local BB emission model are definitely approximate by themselves, so that a too-accurate fitting of *Hi*<sup>∞</sup> *<sup>E</sup>* is actually purely academic. However, these results can be helpful for the interpretation of observations, especially because the 2BB model is often used by observers.

#### **4. Iron and Accreted Heat Blankets**

Paper I considered heat blankets made of iron. Here, we employ a more general model [10,11,22], in which a heat blanket consists of shells (from top to bottom) of hydrogen, helium, carbon and iron. Light elements (H, He, C) are thought to appear because of accretion of H and/or He and further burning into C. Iron is either formed at neutron star birth or is a final result of carbon burning. The mass of lighter elements Δ*M* in the heat blanket is treated as a free parameter. The heat blanketing envelope may fully consist of iron (Paper I) or contain some mass Δ*M* - <sup>10</sup>−<sup>6</sup> *<sup>M</sup>* of lighter elements. Lighter elements affect the surface temperature distribution and thermal surface spectral fluxes.

To simplify the task, we will study fully accreted blankets and conclude on the partly accreted ones in the end of Section 4.2.

#### *4.1. Ts-Distributions*

Our surface temperature distribution *T*s is axially symmetric and symmetric with respect to the magnetic equator. For illustration, Figure 2 presents calculated *T*s for a star with *<sup>M</sup>* <sup>=</sup> 1.4 *<sup>M</sup>* and *<sup>R</sup>* <sup>=</sup> 12 km [*x*g=0.344 and the surface gravity *<sup>g</sup>*<sup>s</sup> <sup>=</sup> *GM*/(*R*2<sup>1</sup> <sup>−</sup> *<sup>x</sup>*g) = 1.59 × <sup>10</sup><sup>14</sup> cm s−2]. The effective surface temperatures are shown versus colatitude *<sup>ϑ</sup>*, which varies from *ϑ* = 0 at the north pole, to *ϑ* = 90◦ at the magnetic equator, and then to *ϑ* = 180◦ at the south pole. Three panels (a–c) correspond to three values of the magnetic

field *B*pole = 1012, 1013 and 1014 G at the pole. Each panel displays *T*s(*ϑ*) for six values of log *T*eff [K]= 5.8, 6.0, . . . 6.8. Solid lines are plotted for the heat blankets made of iron, while dashed lines are for the fully accreted blankets. The iron and accreted matter have different thermal conductivities and, hence, different *T*s(*ϑ*) distributions at the same *T*eff. The dotted lines correspond to a non-magnetic star to guide the eye (they give *T*<sup>s</sup> = *T*eff).

**Figure 2.** Effective surface temperature *T*<sup>s</sup> versus colatitude *ϑ*, from magnetic pole (*ϑ* = 0) to equator (*ϑ* = 90◦) and the opposite magnetic pole (*ϑ* = 180◦) of a neutron star with *M* = 1.4 *M* and *R* = 12 km at three values of *B*pole = 1012, 1013 and 10<sup>14</sup> G [panels (**a**–**c**)]. The curves are plotted for six effective surface temperatures, log *T*eff [K] = 5.8, 6.0 . . . 6.8. Solid curves refer to the heat-blanket model made of iron, dashed cures refer to the heat blanket of fully accreted matter; dotted curves correspond to the field-free star.

Heat conduction is mainly radiative in the outer non-degenerate layer. In deeper layers, where the electrons become degenerate, heat is mostly transported by electrons. The field affects also thermodynamic properties of the matter, for instance, the pressure (e.g., [23]). The magnetic field effects, which regulate the anisotropy of the surface temperature distribution, are twofold.

The effects of the first type are the classical effects of electron rotation about magnetic field lines. They are especially important near the equator, where the heat, emergent from the stellar interiors, propagates mainly across field lines and is thus suppressed. This may strongly enhance thermal insulation of the equatorial part of the heat blanketing envelope, producing rather narrow equatorial dips of *T*s(*ϑ*) (the cold equatorial belt, Figure 2). The dips for the accreted matter are seen to be much stronger than for the iron matter, especially at lower *T*eff.

The effects of the second type are associated with the Landau quantization of electron motion across *B*-lines. They become more efficient at higher *B* near magnetic poles, where

ϑ

they make the heat insulating layers more heat transparent and tend to increase *T*s. These relatively warmer polar 'caps' are wide (Figure 2); *T*s increases inside them, if we approach a pole, but not dramatically.

As long as *B*pole is rather weak (*B*pole - 10<sup>11</sup> G), the classical effects stay more important, producing a pronounced colder equatorial belt, but weakly affecting the polar zones. The surface becomes overall colder, compared with the field-free star with the same internal temperature *T*<sup>b</sup> (e.g., see Figures 22 and 23 in [14]). With increasing *B*pole, the quantum effects become dominant, creating hot polar zones and making the surface overall hotter (at the same *T*b). At *B*pole 10<sup>12</sup> G, the surface becomes hotter than at *B*pole = 0; hot polar zones make the equatorial belt less significant.

It is remarkable (Figure 2) that, outside the equatorial belt, at fixed *T*eff, the surface temperature profiles *T*s(*ϑ*) for iron and accreted heat blankets are mainly close. Moreover, these profiles almost coincide at *<sup>B</sup>*pole ≈ 1013 G. The latter conclusion is independent of specific values of *M* and *R*.

#### *4.2. Spectral Fluxes*

Despite the dramatic interplay of the opposite effects and substantial temperature anisotropy, the surface-integrated spectral fluxes *H*<sup>∞</sup> *<sup>E</sup>* of thermal radiation behave smoothly under variations of *B*pole and/or *T*eff. This is because of strong effects of General Relativity, which smooth out photon propagation from the surface to the distant observer [8,10,11]. As noted in the cited papers, General Relativity 'hides' magnetic effects on thermal surface emission. Very large magnetic fields, with *B*pole 10<sup>14</sup> (typical of magnetars), start to affect the surface emission much more strongly, but such fields are not considered here. They deserve a special study. In particular, their heat-insulating envelopes can become thick and the 1D approach to calculating the heat insulation within them can be questionable.

We have checked that the spectral fluxes, calculated for the accreted envelopes, possess the same properties (Section 3.3) as for the iron envelopes. In particular, the flux remains to be close to that calculated in the 1BB approximation at uniform surface temperature *T*<sup>s</sup> = *T*eff (see, e.g., Figure 1 in Paper I). It is disappointing, meaning that such a flux, taken at *T*eff as an input parameter, carries little information on the magnetic field. Then, one should study deviations from the 1BB approximation. It is important that if one fixes *T*<sup>b</sup> and varies *B*pole, one obtains noticeable *T*<sup>s</sup> and flux variations (as demonstrated, for instance, in Figures 23 and 24 of [14]) but these variations are accompanied by variations of *T*eff. This would be a good theoretical construction, if one knew *T*b. However, the observer can measure *T*eff and wishes to infer *T*b, which is not easy because at fixed *T*eff, the theoretical spectral flux is slightly sensitive to *T*b.

Figure <sup>3</sup> shows thermal spectral flux densities *H*<sup>∞</sup> *<sup>E</sup>* for polar observations of the same star as in Figure 2 (at seven values of *T*eff and three values of *B*pole). The fluxes are plotted in logarithmic scale versus the decimal logarithm of the redshifted photon energy *E*. For each value of *T*eff we present three curves. The solid curves give the fluxes radiated from the star with iron heat blanket. The dashed curves present similar fluxes for accreted heat blankets. The dotted curves demonstrate the 1BB model, Equation (4), with given constant *T*<sup>s</sup> = *T*eff (as if the magnetic field is absent).

All three fluxes for each *T*eff look close in the logarithmic format of Figure 3. However, they do not coincide (as clearly seen from Figure 3 of Paper I). Note that the difference between the fluxes emitted from the star with accreted and iron envelopes is noticeably smaller than the difference between these fluxes and 1BB ones. The difference from 1BB fluxes monotonically increases with growing *B*pole and *E*, which is quite understandable. Higher *B*pole produces stronger anisotropy of *T*<sup>s</sup> distribution. Radiation at higher energies is predominantly emitted from hotter places of the surface.

Let us remark that at *B*pole = 1013 G the spectral fluxes from the star with iron and accreted heat blankets are almost indistinguishable (but sufficiently different from the 1BB flux). This is because the *T*<sup>s</sup> distributions for iron and accreted blankets are very close (Figure 2, Section 4.1). They differ only in cold equatorial belts, but the contribution into the fluxes from cold narrow belts appears almost negligible (also see Paper I). In contrast, the contribution from hotter surface regions can be important (Section 5). At *B*pole = 1014 G the situation remains nearly the same as at 1013 G.

**Figure 3.** Thermal spectral fluxes observed from magnetic poles for a neutron star with *M* = 1.4 *M* and *R* = 12 km at three values of *B*pole = 1011 G (**a**), 10<sup>13</sup> G (**b**), and 1014 G (**c**) at log *T*eff = 5.6, 5.8, . . . , 6.6 assuming iron (solid lines) or accreted (dashed lines) heat blankets. The dotted lines refer to non-magnetic star to guide the eye.

Figure 4 shows similar spectral fluxes but for equator observations (and only for *B*pole = 10<sup>11</sup> and 10<sup>13</sup> G). In this case, the situation is similar to that for pole observations.

The main outcome of Figures 3 and 4 is that the difference of the fluxes emitted from stars with iron and accreted heat blanketing envelopes is rather small and can be ignored in many applications, especially taking into account approximate nature of heat blanketing models. The results for partly accreted heat blankets would be similar. Accordingly, one can neglect chemical composition of heat blankets for calculating thermal emission from magnetized neutron stars in the adopted model of heat blankets. Therefore, the chemical composition of heat blankets does affect cooling of the neutron stars but almost does not affect thermal emission (under formulated assumptions).

**Figure 4.** Same as in Figure 3 but for equator observations at *B*pole = 10<sup>11</sup> G (**a**) and 1013 G (**b**).

#### *4.3. Phase-Space Maps*

As discussed in Section 3.3 (see also Paper I), thermal spectral fluxes can be presented by maps of fit parameters (*p*c, *p*h, *s*c, *s*h) as functions of input parameters. For completeness, Figure 5 shows these 2BB maps versus log*T*eff at *B*pole = 1011 G [panels (a) and (e)], 1012 G [panels (b) and (f)], 1013 G [(c) and (g)] and 1014 G [(d) and (h)]. Panels (a)–(d) refer to pole observations, while (e)–(h) to equator observations. Solid lines correspond to the iron heat blankets, while dashed lines correspond to the accreted ones.

One can see that the dependence of fit parameters on *T*eff is smooth, without any specific features. The parameters of both effective BB components are of the same order of magnitude. Most importantly, the temperature *T*effh of the hotter BB component exceeds the temperature *T*effc of the colder component by less than 20%. The effective surface fraction *s*<sup>c</sup> of the colder component slightly exceeds the fraction *s*<sup>h</sup> of the hotter component.

The maps in Figure 5 represent the same thermal fluxes, which are plotted in Figures 3 and 4. The main conclusion of Section 4.2, based on Figures 3 and 4, is that the fluxes for iron and accreted envelopes are nearly the same. If they were identical, then respective solid and dashed curves in Figure 5 should have been identical as well. However, these curves differ. The difference is especially visible for *B*pole = 1011 and 10<sup>12</sup> G, but less visible at higher *B*pole. The difference is mainly explained by two things. First, weak variations of 2BB fit parameters lead to weaker variations of the fluxes. Secondly, the 2BB approximation is accurate but not exact by itself, as discussed in the previous section.

It is important that we fit *unabsorbed* calculated spectral fluxes using *unabsorbed* 2BB models. Accordingly, as argued in Paper I, one can compare 2BB parameters, inferred from observations (corrected for interstellar absorption and non-thermal radiation component), with the theoretical parameters of unabsorbed 2BB models.

**Figure 5.** The 2BB fit parameters versus log *<sup>T</sup>*eff for a 1.4 *<sup>M</sup>* neutron star with *<sup>R</sup>* <sup>=</sup> 12 km and *<sup>B</sup>*pole <sup>=</sup> 1011, 1012, 1013, or 1014 G [panels (**a**–**h**)] for polar (*i* = 0) and equator (*i* = 90◦) observations. Solid lines refer to the star with iron heat blanket, while dashed ones are for accreted heat blanket.

#### **5. Extra Heating of Magnetic Poles**

Section 4 describes the thermal emission of cooling neutron stars with dipole magnetic field near their surfaces, using the *T*s model calculated in [10,11]. According to the theory, radiative spectral fluxes should be almost insensitive to the chemical composition of the heat blanketing envelopes; 2BB fits to their thermal X-ray spectral fluxes do not give the difference between the temperatures *T*effh and *T*effc of the hotter and colder BB components higher than 20% (for *B*pole -10<sup>14</sup> G).

A brief comparison with observations of isolated middle-aged neutron stars in Paper I shows that there are no reliable candidates for such objects at the moment; the difference between *T*effh and *T*effc is actually larger, at least ∼50% for RX J1856.5−3754 [24,25], as a promising example.

Here, we propose a possible phenomenological extension of the discussed model, which may simplify the interpretation of observations of some sources. Specifically, let us assume the presence of additional heating of magnetic poles, which produces polar hotspots and raises the pole temperature. Such an extension has been used, for instance, in Paper I but assuming heating of both poles. Now we consider heating either of both poles or one of them. For simplicity, this heating is assumed to be axially symmetric. It does not violate the axial symmetry of the *T*<sup>s</sup> distribution. Note that two equivalent hotspots do not destroy the initial symmetry of the north and south hemispheres, whereas one hotspot does destroy it. Possible physical justifications for extra heating are mentioned in Section 6.

Let *T*s0(*ϑ*) be the basic effective surface temperature used in Section 3. We introduce a small phenomenological angle *ϑ*<sup>0</sup> that deteminies the size of a hotspot. Following Paper I we assume that at *ϑ* < *ϑ*<sup>0</sup> (in the nothern hotspot),

$$T\_{\rm s}(\theta) = T\_{\rm s0}(\theta) \left[ 1 + \delta \cos^2 \left( \frac{\pi \theta}{2 \theta\_0} \right) \right] \quad \text{at } \theta \le \theta\_0 \tag{8}$$

and *T*s(*ϑ*) = *T*s0(*ϑ*) outside the hotspot. The parameter *δ* specifies an extra temperature enhancement at the magnetic pole; the enhancement smoothly disappears as *ϑ* → *ϑ*0. This will be our *model with one hotspot*. Assuming a similar temperature enhancement near the second magnetic pole, we will get the *model with two hot spots*. The presence of spots renormalizes the total effective temperature *T*eff, Equation (3).

Otherwise, the computation of spectral fluxes is the same as in Section 3. These fluxes can also be approximated by 2BB fits (6), and can be analyzed via phase-space maps. In some cases presented below the relative the fit accuracy becomes worse (reaching sometimes ∼10 %) but the fit remains robust. We will again take the star with *M* = 1.4 *M* and *R* = 12 km. For certainty, we use the model of iron heat blanket and *ϑ*<sup>0</sup> = 10◦.

#### *5.1. Extra Hotspots on Both Poles*

The left panel of Figure 6 demonstrates calculated spectral fluxes for *B*pole = 10<sup>14</sup> G and seven values of log *T*eff[K]=5.6, 5.8, . . . , 6.8 at *δ* = 0.3. Since these hotspots are identical, observations of the first and second poles (*i* = 0 and 180◦) give the same fluxes. Such fluxes are shown by dashed lines. For comparison, the solid lines are calculated at *δ* = 0, in which case the extra hotspots are absent and the results of Section 3 apply. The dotted lines show the fluxes for non-magnetic star. One can see that the presence of hotspots enhances the spectral fluxes at high photon energies. This is expected: enhanced surface temperature of hotspots intensifies generation of high-energy radiation. The higher the *T*eff, the larger the photon energies are affected. The solid and dotted lines are the same as in Figure 3c.

**Figure 6.** (**Left**): Thermal spectral fluxes observed from magnetic poles of a neutron star with *M* = 1.4 *M* and *R* = 12 km at *B*pole = 10<sup>14</sup> G and log *T*eff [K]= 5.6, 5.8,. . . ,6.6. The solid lines are obtained without extra heating, while the dashed lines are for two extra hotspots of angular size *ϑ*<sup>0</sup> = 10◦ on magnetic poles at the heating parameter *δ* = 0.3. Dots show 1BB model for a non-magnetic star to guide the eye. (**Right**): 2BB fit parameters for the same star with *T*eff = 1 MK possessing two polar hotspots of angular size *ϑ*<sup>0</sup> = 10◦ on magnetic poles [see Equation (8)] versus extra relative surface temperature increase *δ* at the pole.

The right panel of Figure 6 shows the 2BB fit parameters versus *δ* for the same star with *T*eff=1 MK. At *δ* = 0, the results naturally coincide with those in Figures 5g,h. However, with increasing *δ*, the fit parameters become different. The temperature *T*effh of the hotter BB component and the effective emission surface area *s*c of the colder component noticeably increase, whereas the emission surface area of the hotter component dramatically falls down. Even with really small hotspot temperature enhancement *δ* ≤ 0.3 (which gives the fraction of extra luminosity ≤1.4%), one obtains an absolutely new phase-space portrait with *s*<sup>h</sup> *s*c. Such 2BB fits have been often inferred from observations of cooling isolated neutron stars (see, e.g., Paper I); these sources are usually interpreted as neutron stars, which have small hotspots with noticeably enhanced temperature.

Therefore, the theory with hotspots predicts (e.g., Paper I) two types of neutron stars, whose spectra can be approximated by the 2BB models. The first ones are those with smooth surface temperature distributions, created by nonuniform surface magnetic fields (as considered in Section 3) to be called spectral 2BB models *of smooth magnetic atmospheres*. The second sources are those with distinct hotspot BB component to be called 2BB *with hotspots*. Naturally, there is a smooth transition between them (for instance, by increasing *δ* in the right panel of Figure 6). It seems that the observations do not provide good candidates for the sources of the first type (Paper I), but there are some candidates for the sources of the second and intermediate types as we will outline later.

#### *5.2. Phase Resolved Spectroscopy in Case of Two Hotspots*

Let us present a few calculated lightcurves to illustrate an important problem of pulse fraction. For simplicity, we consider the star as an orthogonal rotator with the spin axis perpendicular to the line of sight. Figure 7 presents several lightcurves for a star with *B*pole = 10<sup>14</sup> G versus phase angle *θ*. The lightcurves are normalized by *HE*(0) (the spectral flux for pole observations at given energy *E*). The displayed ratio of redshifted or non-redshifted fluxes is the same at a given *E*, so that we drop the symbol ∞. Note that the normalization flux *HE*(0) strongly depends on energy by itself. Panel (a) corresponds to *δ* = 0, in which case the hotspots are actually absent and we have the emission from the smooth magnetic atmosphere. In case (b) the hotspots with *δ* = 0.3 are available and strongly affect the lightcurves.

**Figure 7.** Theoretical lightcurves, normalized to pole observations, versus phase angle *θ* for a neutron star with *T*eff = 1 MK and *B*pole = 1014 G; (**a**) *δ* = 0 (no extra pole heating) and (**b**) *δ* = 0.3. Solid curves are for photon energies log *E*[keV]= −1, −0.8, . . . , 0.8 (from top to bottom). The two lowest curves on panel (**b**) (with highest *E*) almost coincide. Dotted curves show the same but neglecting gravitational light bending. The star is assumed to be orthogonal rotator with the spin axis perperndicular to the line of sight.

Solid lines on each panel in Figure 7 present the lightcurves at 10 energies (from top to bottom), log *E*[keV]= −1, −0.8, ... , 0.8. The higher *E*, the stronger phase variations. The last two lines on panel (b) almost coincide. The lightcurves on panel (a) are seen to be smooth. The curves on panel (b) at *E* - 1 keV are although smooth and resemble those on panel (a). However, at higher energies the lightcurves (b) have shapes with erased dips at nearly equator observations [and then Equation (7) becomes inaccurate although

Equation (6) works well]. These curves are typical for lightcurves produced by antipodal point sources on the neutron star surface; see, e.g., Figure 4 in [19] or Figure 6 in [21]. Clearly, at high energies the star emits from the poles, the hottest places on the surface. Nevertheless, the pulse fraction remains not too high (≤30%) in both cases (a) and (b) even at the highest taken energy *E* ≈ 6.3 keV, at which thermal spectral flux becomes typically very low by itself. Therefore, the enhancement of thermal emission by extra heating of two magnetic poles does not lead to sizeable pulse fractions.

The dotted lines on Figure 7 show the same lightcurves as the solid lines but calculated neglecting gravitational bending of light rays. One sees that without the light bending the pulse fraction would be much stronger (as is well known; see, for example, [8]). The difference of solid and dotted curves on Figure 7 has simple explanation. Without light bending, the equator observations (*θ* = 90◦ or 270◦) would collect more light emitted from the cold equator. This would lower the observed emission and produce stronger dips of the lightcurves.

#### *5.3. Extra Hotspot on One Pole*

Finally, consider the case of a single hotspot, which we put at the north pole. This case is different from the previous one.

The left panel of Figure 8 presents spectral fluxes calculated for the same values of *T*eff and *B*eff as in Figure 6. The dashed lines are for the north-pole observations (*i* = 0) at *δ* = 0.3. These fluxes are the same as in Figure 6: the observer sees a large fraction of the surface (mainly the northern hemisphere); the southern polar region is unseen.

The solid line in the left panel of Figure 8 refers to the same case of *δ* = 0.3 but the star is observed from the south pole (*i* = 180◦). Then the observer cannot see the hotspot, so the line exactly coincides with the solid line in Figure 3c. The dotted line in Figure 8 is for a non-magnetic star (as in Figures 3c and 6). Now, with only one extra heated pole, the difference between observations of the north and south poles is pronounced.

δ

**Figure 8.** (**Left**): Theoretical spectral fluxes for a neutron star with *<sup>M</sup>* <sup>=</sup> 1.4 *<sup>M</sup>*, *<sup>R</sup>* <sup>=</sup> 12 km, *<sup>B</sup>*pole <sup>=</sup> 1014 G, different *<sup>T</sup>*eff and one extra hotspot on one magnetic pole at the extra heating parameter *δ* = 0.3. The solid and dashed curves are for non-heated (*i* = 180◦) and heated pole (*i* = 0) observations, respectively. (**Right**): Maps of fit parameters versus *δ* for this star at *T*eff = 1 MK. The thick solid and gray curves are for the north-pole and equator observations, respectively. Thin horizontal curves are for the south-pole observations.

The right panel of Figure 8 presents the maps of 2BB fit parameters versus *δ* for the same star at *T*eff = 1 MK for the north-pole, south-pole, and equator observations (to be compared with Figure 6). The maps the for north-pole and equator observations are nearly the same as those in Figure 6, but the maps for the south-pole observations are plain. Corresponding fit parameters are independent of *δ* just because the observer cannot see the emission from the hotspot (which is the only emission which increases with *δ*).

#### *5.4. Phase Resolved Spectroscopy in Case of One Hotspot*

Phase-resolved lightcurves are plotted in Figures 9 and 10. They are naturally different from the case of two extra heated magnetic poles. Figure 9 is analogous to Figure 7 and shows the light curves at the same *T*eff, *B*pole, and extra heating parameter *δ* = 0.3 (a) and *δ* = 0.14 (b). With one hotspot, the fractions of extra luminosities are 0.7% (a) and 0.3% (b). The phase-space variations are similar to those in Figure 7 only at low photon energies, in which case the extra polar heating is almost unnoticeable and the lightcurves have two shallow dips over one rotation period. However with growing *E*, the extra heating becomes important and increases the depth and width of the dips (such shapes are known to be produced by point sources on neutron star surface, e.g., Figure 4 in [19]). At *E* 4 keV in Figure 9a the pulse fraction becomes close to 100%. In that case the observer would see only one, but very pronounced dip over one rotation, because the extra hotspot would be poorly visible in a wide range of *θ* with the center at *θ* = 180◦. In Figure 9b the heating is weaker, the dips are smaller, but nevertheless they are much larger than for two heated spots at *δ* = 0.3 in Figure 7. Evidently, one extra-heated magnetic pole is much more profitable than two, for producing large pulse fractions.

**Figure 9.** Theoretical lightcurves for a neutron star with *M* = 1.4 *M*, *R* = 12 km, *T*eff = 1 MK and *B*pole = 1014 G with an extra hotspot of angular size *ϑ*<sup>0</sup> = 10◦ on one magnetic pole versus phase angle *θ* at extra heating parameter *δ* = 0.3 (**a**) and 0.14 (**b**). The curves are for photon energies log *E* [keV]= −1, −0.8, . . . , 0.8. See the text for details.

Figure 10 is plotted for *B*pole = 1011 G to illustrate the effects of *B*pole on the pulse fraction. Figure 10a presents the normalized lightcurves for the star with *δ* = 0.3. Now the pulse fraction is lower than at *B*pole = 1014 G. Nevertheless, it remains much higher than it would be if the additional heating was switched on at both poles. Figure 10b demonstrates that at lower *δ* = 0.14 at *B*pole = 10<sup>11</sup> G the pulse fraction is weaker affected by the extra heating but nevertheless reaches about 70% at *E*∼6 keV.

**Figure 10.** Same as in Figure 9 but for *B*pole = 10<sup>11</sup> G.

#### **6. Discussion and Conclusions**

Several isolated middle-aged neutron stars which radiation spectra have been (or can be) fitted by 2BB models have been selected in Paper I, based on the recent catalog [26] presented also on the Web: www.ioffe.ru/astro/NSG/thermal/cooldat.html, accessed on 17 October 2021.

These can be cooling neutron stars with dipole magnetic fields and additionally heated polar caps. This selection is illustrative (does not pretend to be complete).


6. RX J1856.5−3754 is a neutron star with nearly thermal spectrum. It was discovered by Walter et al. [39]. It is rotating with the spin period of ∼7 s; the effective magnetic field is *<sup>B</sup>*eff∼1.5 × <sup>10</sup><sup>13</sup> G; magnetic properties are highly debated (e.g., [18,40] and references therein). The spectrum has been measured in a wide energy range, including X-rays, optical and radio, and interpreted with many spectral models, particularly, with the model of thin partially ionized magnetized hydrogen atmosphere on top of solidified iron surface (see, e.g., [5,41]). Note alternative 2BB, 2BB+PL, and 3BB fits constructed by [24] and [25]. The 2BB fits give *k*B*T*<sup>∞</sup> effc∼40 eV (*T*<sup>∞</sup> effc ≈ 0.46 MK), *T*effh/*T*effc∼1.6, *R*effh∼0.5 *R* and *R*effc∼*R*; they seem closer to the 2BB spectral models with smooth magnetic atmosphere, than 2BB fits for other sources. Extra heating of magnetic poles would simplify this interpretation.

This brief analysis does not reveal any good candidate which would belong to the family of neutron stars with dipole magnetic fields and smooth surface temperature distribution (Section 4). This does not mean that such neutron stars do not exist; it can be difficult to identify them because their spectra are close to 1BB spectrum. To interpret the above sources one needs to complicate the model, for instance, by introducing extra heating of magnetic poles.

Now let us summarize our results. Following Paper I we have studied simple models (Section 3) of thermal spectra emitted from surfaces of isolated neutron stars with dipole surface magnetic fields 1011 - *B*pole - 1014 G. Such fields make the surface temperature distribution noticeably non-uniform. The model assumes BB emission with a local temperature from any surface element.

In Section 4 we have shown that the spectral X-ray fluxes emitted from such neutron stars are almost the same for heat blanketing envelopes composed of iron and fully accreted matter (at the same average effective surface temperatures *T*eff). By measuring the spectral fluxes and *T*eff, it would be difficult to infer composition of their heat blanketing envelopes, although this composition can noticeably affect neutron star cooling (e.g., [14] and references therein). In the presence of fully accreted matter, the fluxes are accurately fitted by 2BB models (as in Paper I for the iron heat blankets). At *M*∼1.4 *M* and *R*∼12 km, the ratio *T*effh/*T*effc of effective temperatures of the hotter to colder BB components cannot be essentially larger than ∼1.2. For less compact stars, with smaller *M*/*R*, this value can be somewhat higher (see e.g., [8], Paper I). If this ratio is larger from observations, then the model of Sections 3 and 4 cannot explain the data.

In Section 5, following Paper I, we have extended the model by assuming sufficiently weak extra heating of one or both magnetic poles. Increasing extra relative temperature rise *δ* at the pole(s), does not greatly violate the accuracy of the 2BB spectral approximation, but allows one to obtain larger *T*effh/*T*effc. Also, it affects the pulse fraction, which increases with *B*pole and photon energy *E*. The pulse fraction can be very high if one magnetic pole is additionally heated, while the other is not.

Let us stress, that there have been many other studies of thermal emission from magnetized neutron stars. The emission produced by anisotropic surface temperature distribution of the star with dipole magnetic fields was analyzed by Page [8]. Quadrupole magnetic field components were added in [9]. A possible presence of toroidal field component in the neutron star crust was studied in [12,13]. Also, there were studies of thermal emission during magnetic field evolution in neutron stars (e.g., [18,40]). These works give a variety of magnetic fields configurations, surface temperature distributions, and phase-resolved spectra of neutron stars to be compared with observations and determine simultaneously both, the surface temperature distribution and magnetic field geometry.

Introducing some extra heating of magnetic poles, we are actually doing the same. Note that hot spots can be produced by pulsar mechanism. Also, the toroidal crustal magnetic fields can significantly widen the cold equatorial belt, creating relatively stronger heated poles [12], which can mimic hotspots. At the first step it might be simpler to introduce phenomenological polar heating, that can be different on the two poles, and to determine the parameters of the extra heating (the *T*s distribution) by varying the parameters like *δ* and extra temperature profile [Equation (8)] at both poles. This would be a simple and physically transparent model of surface temperature distribution. At the next step, one could infer (constrain) the magnetic field geometry. This method can be easily extended to different *T*s distributions, not necessarily axially symmetric.

Another very important direction would be to go beyond the approximation of local BB emission from any surface patch, and use more realistic emission models (with spectral features and anisotropic radiation, as well as with account of polarization effects, for example, [42,43]).

**Funding:** This research was partly (in analyzing the case of accreted matter) supported by the Russian Science Foundation, grant 19-12-00133.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data underlying this article will be shared on reasonable request to the corresponding author.

**Acknowledgments:** I am grateful to Oleg Gnedin for providing me photos of his Father.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

BB blackbody

#### **References**


## *Review* **Progress in Understanding the Nature of SS433**

**Anatol Cherepashchuk**

P. K. Sternberg State Astronomical Institute, M. V. Lomonosov Moscow State University, Universitetskij Prospect, 13, 119234 Moscow, Russia; cherepashchuk@gmail.com

**Abstract:** SS433 is the first example of a microquasar discovered in the Galaxy. It is a natural laboratory for studies of extraordinarily interesting physical processes that are very important for the relativistic astrophysics, cosmic gas dynamics and theory of evolution of stars. The object has been studied for over 40 years in the optical, X-ray and radio bands. By now, it is generally accepted that SS433 is a massive eclipsing X-ray binary in an advanced stage of evolution in the supercritical regime of accretion on the relativistic object. Intensive spectral and photometric observations of SS433 at the Caucasian Mountain Observatory of the P. K. Sternberg Astronomical Institute of M. V. Lomonosov Moscow State University made it possible to find the ellipticity of the SS433 orbit and to discover an increase in the system's orbital period. These results shed light on a number of unresolved issues related to SS433. In particular, a refined estimate of the mass ratio *Mx Mv* > 0.8 was obtained (*Mx* and *Mv* are the masses of the relativistic object and optical star). Based on these estimates, the relativistic object in the SS433 system is the black hole; its mass is >8*M*. The ellipticity of the orbit is consistent with the "slaved" accretion disc model. The results obtained made it possible to understand why SS433 evolves as the semi-detached binary instead of the common envelope system.

**Keywords:** close binaries; black holes; evolution of binary stars

#### **1. Introduction**

The unique object SS433 shows in its optical spectrum [1], in addition to the standard emission lines of hydrogen and helium, moving emissions, which move along the spectrum with a period of 162.3 d [2]. The amplitude of these displacements reaches an enormous value of ~1000 Å, which amounts up to tens of thousands of km s−<sup>1</sup> on the scale of velocities. Before the discovery of SS433 in the 1970s [2,3], astronomers had never seen anything like this. Therefore, the nature of object SS433 seemed mysterious; some scientists called SS433 "the enigma of the century". In the earliest press releases, there were even mentions that it was possible that we were observing signals from an extraterrestrial intelligence which was shining on us with a powerful re-configurable laser. Fortunately, this hypothesis did not manage to appear in scientific publications since, in 1979, it turned out [4,5] that the moving lines in the spectrum of SS433 are formed in relativistic (*v* 0.26*c*) collimated jets that are ejected from some central source and precess with a period of ≈164 d. The opening angle of the precession cone is ≈20◦, and the precession axis is inclined with respect to the line of sight by an angle of ≈79◦. Although the nature of the central source remained unclear, speculation arose that all of these phenomena were enacted in a binary system [6,7]. The binary system model was proposed in 1980 by [8], who measured the Doppler shifts of the narrow components of the stationary hydrogen lines and found a period of ≈13.1 d, which testified in favor of the close binary system model. The authors from the analysis of their data, assuming that the radial velocity curve measured by them reflects the orbital motion of the components, came to the conclusion that SS433 is a low-mass binary system consisting of a low-mass optical star (1*M*) and a neutron star. As it turned out later, according to our research, this conclusion turned out to be incorrect.

In 1980, we made photometric observations of SS433 and found that SS433 is an eclipsing binary system [9]. The discovery of optical eclipses in the SS433 system made it

**Citation:** Cherepashchuk, A. Progress in Understanding the Nature of SS433. *Universe* **2021** , *8* , 13. https://doi.org/10.3390/ universe8010013

Academic Editors: Nazar R. Ikhsanov, Galina L. Klimchitskaya and Vladimir M. Mostepanenko

Received: 26 October 2021 Accepted: 24 December 2021 Published: 27 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

possible to establish that the radial velocity curve constructed by [8] does not reflect the orbital motion of the components but describes the motion of gas in a gas stream flowing from the Lagrange point L1 onto a relativistic object. Analysis of optical eclipses in SS433, as well as studies of the system's extra eclipse brightness depending on the phase of the precessional period [9], made it possible to construct a self-consistent model of SS433 as a massive X-ray eclipsing binary system at an advanced stage of evolution with an optically bright supercritical accretion disc predicted in [10]. Relativistic jets are ejected from the central parts of the supercritical disc perpendicular to its plane. The plane of the disc is inclined to the orbital plane of the SS433 system at an angle of 20◦, the disc precesses with a period of 162.3 d and relativistic jets track the disc precession. The moving emissions are formed in the outer parts of the jets.

The question arose about the reasons for the disc precession. Back in 1973, Katz proposed a precessing disc model to explain the 35-day X-ray cycle in the Her X-1 system [11]. In 1974, Roberts proposed a slaved precession model for an accretion disc in the Her X-1 system [12]. Even earlier, the idea of the slaved disc was expressed by Shakura [13]. The slaved disc model for SS433 was proposed in 1980 by van den Heuvel et al. [14]. In this model, the axis of rotation of the optical star precesses and the disc tracks this precession. The non-perpendicularity of the axis of rotation of the optical star to the orbital plane may be associated with an asymmetric supernova explosion, which inclined the orbital plane of the binary system relative to the axis of rotation of the star [12,15]. To explain the high and stable velocity of matter in the SS433 jets, a model of radiative acceleration was proposed (the so-called line-looking effect [16,17]).

Thus, by the early 1980s, the main features of SS433 as a massive close binary X-ray system in the supercritical accretion regime had been clarified. However, many problems remained unresolved. The massive binary system model for SS433 [9] was confirmed by [18], who plotted the radial velocity curve from the stationary HeII 4686 line and measured the mass function of a relativistic object (the lower limit of the mass of an optical star), which turned out to be ≈10*M*.

#### **2. Unsolved Problems**

Over the past 40 years, the object SS433 has been studied in detail in the optical, IR, X-ray and radio ranges (see, for example, reviews [6,7] where all the necessary references are given). Many new and interesting results were obtained on the object SS433, which turned out to be the first example of a microquasar in our Galaxy. However, the main problems associated with this unique object have been unresolved until recently. Let us list some of them.


Recently, new observational data have been obtained on the object SS433, which bring us closer to solving these problems.

#### **3. Discovery of the Ellipticity of SS433's Orbit: Strong Support for a Slaved Accretion Disc Model**

The slaved accretion disc model proposed in [14] seems to be very attractive since a natural physical justification can be formulated for it. Indeed, a small (~1%) asymmetry of a supernova explosion, which accompanied the formation of a relativistic object in a binary system, can turn the orbital plane of the binary system relative to the axis of rotation of the optical star by a significant angle [12,15]. A supernova explosion could also significantly increase the eccentricity of the system's orbit. Under the action of tidal interaction from the side of the resulting relativistic object, the axis of rotation of the optical star begins to precess, which leads to the formation of an accretion disc that does not lie in the plane of the system's orbit and precesses with the precession period of the optical star. Since a precessing star has a large mass and a huge rotational moment, its precession is very stable and, as a consequence, the observed precession of the accretion disc should also be stable on average. Our long-term (24 years) spectral observations of SS433 [19,20] using published data indicate the stability of the parameters of the kinematic model SS433 over 40 years. This stability can be considered as an argument in favor of the model of a slaved accretion disc in the SS433 system.

However, until recently, there was one fundamental difficulty in applying the model of a slaved accretion disc to the SS433 system. The orbit of SS433 was considered circular since the secondary minimum on the eclipsing light curve, within the errors associated with significant physical variability of the system, was located approximately in the middle between the two primary minima.

The presence of a circular orbit in the SS433 system during asynchronous rotation of the optical star (the star's rotation axis is not perpendicular to the orbital plane) contradicts the theory of tidal synchronization of stars in binary systems. According to this theory [21,22], the synchronization of axial and orbital rotation due to the dissipation of the energy of orbital motion in dynamic tides should occur before the circularization of the orbit. Therefore, if the orbit of SS433 is circular, then the optical star must be synchronized with the orbital revolution and its rotation axis must be perpendicular to the orbital plane. In this case, the accretion disc should lie in the plane of the system's orbit. The slaved disc model is not applicable in this case. As noted above, a significant initial eccentricity of the orbit could appear in a binary system after a supernova explosion, so the assumption of an initial circular orbit seems unnatural. In this regard, the problem of finding direct signs of the ellipticity of the SS433 orbit arises. Indirect evidence of the ellipticity of the orbit was obtained in [19,20,23]. In these works, it was found that the velocity of matter in relativistic jets SS433 is modulated by the orbital period *P*orb = 13.082 d. In this regard, a hypothesis was put forward about the possible ellipticity of the SS433 orbit: changes in the distance between the components in an elliptical orbit cause a change in the rate of influx of stellar matter into the accretion disc, which can lead to orbital modulation of the velocity of matter in the jets. Other indirect indications of the ellipticity of the SS433 orbit were obtained in [24,25]. In this regard, it is of interest to search for direct evidence of the ellipticity of the SS433 orbit.

During 2018–2021, at the Caucasian Mountain Observatory (KGO) of the GAISH MSU, we carried out intensive photometric monitoring of SS433. Several hundred brightness estimates of SS433 were obtained in BVRI filters during 276 observation nights. These data allowed us, in combination with published photometric observations, to construct a high-precision (including 909 individual observations) mean SS433 light curve in the *V* band in the phases of the precessional period *P*prec = 163.3 d, corresponding to the maximum opening of the accretion disc with respect to the observer [26]. Figure 1 shows this curve. It can be seen that the secondary minimum has a noticeable shift relative to the middle between the two primary minima. Analysis of this displacement and the study of the difference in the widths of the primary and secondary minima performed in [26] made it possible to establish that the eccentricity of the SS433 orbit is *e* = 0.050 ± 0.007, and the longitude of the periastron for the epoch of ≈1990 is *ω* = 40◦ ± 20◦. The ellipticity of the SS433 orbit discovered by us is in good agreement with the model of a slaved accretion disc. Both observational facts—the long-term stability of the kinematic model parameters and the ellipticity of the SS433 orbit—strongly support the model of a slaved accretion disc in the SS433 system [14], as well as the hypothesis that the reason that the disc does not lie in the plane orbit is an asymmetric supernova explosion in the system [12,15].

**Figure 1.** Average orbital light curve of SS433 derived from 909 individual measurements taken during 1978–2012. The observed points were selected in the phases of maximum opening of the accretion disc with respect to the observer. From the paper [26].

#### **4. The Discovery of an Increase in the Orbital Period of the SS433 System: A Strong Argument for the Presence of a Black Hole**

Despite more than forty years of research, the nature of the relativistic object in the SS433 system remained completely unclear until recently (see discussions on this topic in [6,7,20]). In [19,20], using information on the constancy of the orbital period of SS433 over 28 years [27,28], the authors estimated the ratio of the masses of the components in the SS433 *<sup>q</sup>* <sup>=</sup> *Mx Mv* = 0.6 (*Mx* and *Mv* are the masses of a relativistic object and an optical star, respectively). The constancy of the orbital period for the SS433 system is very surprising; after all, the rate of mass loss from the system in the form of wind from the supercritical accretion disc is very high, ~10−4*M* year−<sup>1</sup> (see, for example, [9]). In [20,29], the balance equations for the loss of angular momentum of the SS433 system were derived by various mechanisms: the transfer of mass and angular momentum during the flow of matter of an optical star through the internal Lagrange point *L*<sup>1</sup> to a relativistic object, Jeans mass loss in a symmetric high-velocity outflow matter in the form of wind from the disc as well as the outflow of matter outside the binary system through the outer Lagrange point *L*2. The first mechanism, when flowing from a more massive optical star to a relativistic object, tends to bring the components closer together (and, accordingly, to decrease the orbital period). The second mechanism seeks to increase the distance between the components, that is, to increase the orbital period. The third mechanism leads to a decrease in the distance between the components, i.e., to a decrease in the orbital period. It turns out that, at a huge rate of mass loss in the form of a wind ~10−4*M* year−1, an equilibrium between these three mechanisms is established at a very large mass ratio *q* 0.6. With an optical star mass of ≈10*M* (see, for example, [18]), it follows that the mass of a relativistic object is 6*M*, which is typical of a black hole. In this regard, it is very important to check with additional observations the nature of the constancy or change in the orbital period of SS433. As noted above, in 2018–2021, on the automated 60-cm reflector of the Caucasian Mount Observatory of P. K. Sternberg Astronomical Institute of M. V. Lomonosov Moscow State University, we obtained dense series of photometric observations of SS433 in *BVRI* filters. With the involvement of all published photometric observational data (since 1979), we discovered an increase in the orbital period of SS433 with time [26]. Figure 2 shows the corresponding plots *O* − *C*—the differences between the observed moments of the primary minimum of the light curve of SS433 and theoretical moments calculated with linear elements [28]. Here, the parabola corresponds to an increase in the period, and the inclined line corresponds to a constant period (the slope of this line to the abscissa axis determines the correction to a

constant period). It can be seen that the parabolic approximation of the points in Figure 2 is much more preferable than the linear approximation. Thus, it is shown that the orbital period of SS433 increases. The rate of this increase is *<sup>P</sup>*˙ = (1.0 <sup>±</sup> 0.3) <sup>×</sup> <sup>10</sup>−<sup>7</sup> s per s. It is also shown that the model of the third body in the SS433 system, which can also lead to an increase in the orbital period, is rejected, since the mass of the third body must be more than 16*M*. The lines of such a massive star would be visible in the total spectrum of the SS433 system, which is not observed.

**Figure 2.** *O* − *C* values for SS433, calculated using the ephemeris with a constant orbital period *P*orb = 13.08223 d. Right and left panels: *O* − *C* calculated by different methods. The parabola corresponds to an increase in the orbital period. The inclined straight line corresponds to a constant period with a Δ*P* correction. From the paper [26].

Using the equation to balance the loss of angular momentum and taking into account that the orbital period of SS433 increases, an improved estimate of the mass ratio *q* > 0.8 is given in [26]. With the mass of the optical star *Mv* = 10*M*, an estimate is obtained for the mass of the black hole in SS433 *Mx* > 8*M*. Thus, it has been shown that the SS433 system contains a black hole with a mass close to the average mass of stellar black holes in X-ray binaries (≈8 ÷ 10*M*).

The increase in the orbital period of SS433 found by us allows us to reliably reject the presence of a neutron star in the system since, in this case, the orbital period of SS433 should not increase, but should decrease with time, which contradicts observations.

#### **5. Other Estimates of the Mass Ratio**

The supercritical optically bright accretion disc in the SS433 system makes spectroscopic investigations of the donor star more difficult due to the disc's powerful radiation. Nevertheless, there is a significant progress in this problem by now.

Careful spectroscopic investigations of SS443 made by [30–32] let us estimate the spectral type of the optical star as A7I and to obtain (using absorption lines) the radial velocity curve of this star with the semi-amplitude *<sup>K</sup>*<sup>v</sup> = 58.2 ± 3.1 km s<sup>−</sup>1. It corresponds to the optical star's mass function *f*v(*m*) = 0.268*M*. Similar results were obtained by [33]. In the recent paper by [34], the axial rotation velocity *<sup>v</sup>*rot = <sup>140</sup> ± 20 km s−<sup>1</sup> was measured using the Doppler widening of absorption lines of the optical star. Using this value, the corresponding masses and mass ratio of the components of SS433 were estimated: *<sup>q</sup>* <sup>=</sup> *<sup>M</sup>*<sup>x</sup> *<sup>M</sup>*<sup>v</sup> <sup>=</sup> 0.37 <sup>±</sup> 0.04, the mass of the relativistic object *<sup>M</sup>*<sup>x</sup> <sup>=</sup> 4.2 <sup>±</sup> 0.4*M* and the mass of the optical star *M*<sup>v</sup> = 11.3 ± 0.6*M*.

There are reasons to suppose that the estimates of the mass of the relativistic object in SS433 based on the spectroscopic data cannot be assumed as reasonably reliable, because the SS433 radiation travel through the intercomponent moving gas medium and also through the rotating circumbinary gas envelope. For the gigantic mass loss rate *<sup>M</sup>*˙ <sup>10</sup>−<sup>4</sup> year−1, the density of matter in these structures is very high. Therefore, the selective absorption of the optical star's light in moving dense gas structures can significantly distort the orbital radial velocity curve of the donor star measured using absorption lines, see, e.g., [35–37].

The situation with emission lines is also uncertain. For example, the mass function of the relativistic object calculated using the He II 4686 emission line changed from 10*M* [18] to 2*M* [38]; it depends on the phase of the precessional period. This can be related with the complicated structure of the wind from the precessing accretion disc whose shape can be asymmetric and twisted.

Because of these problems, great hopes were pinned on studies of X-ray eclipses in the SS433 system. Since the orbital inclination of the system *i* = 79◦ is known from the analysis of displacements of moving emission lines, and the optical star fills or even overfills its inner critical Roche lobe, the analysis of the duration of the X-ray eclipse (when the star hides X-ray structures) gives a principal possibility to estimate the mass ratio of the components.

First such estimate was made by [39–41]. From the analysis of X-ray eclipses in SS433 in the range 2–10 keV in the model of eclipse of the thin relativistic jets (radiating in X-rays) by the optical star, there was obtained a very low mass ratio *<sup>q</sup>* <sup>=</sup> *<sup>M</sup>*<sup>x</sup> *M*v 0.15. In this case, it was assumed that the optical star fills its Roche lobe.

Nevertheless, there are observational data indicating that the optical star in SS433 not only fills it inner critical Roche lobe but also overfills it. At the same time, the star outflows both through the inner Lagrange point *L*<sup>1</sup> and through the outer Lagrange point *L*2. From the analysis of spectral data, arguments in favor of the outflow of the optical star's matter through the outer Lagrange point *L*<sup>2</sup> were suggested [42]. It indicates that the star most likely fills its outer Roche lobe whose size is significantly (by 15–20%) greater than the size of the inner critical Roche lobe. In papers by [43,44], the equatorial outflow of the matter in the plane perpendicular to the direction of the jets was discovered in S433. It also indicates that the star outflows through the *L*<sup>2</sup> point and that the star fills its outer Roche lobe. The conclusion about the filling of the outer Roche lobe by the optical star in SS433 and about the outflow through the *L*<sup>2</sup> point is supported by the optical and infrared spectroscopy of this object that revealed two-humped stationary hydrogen emission lines [45,46]. These lines indicate the presence of the circumbinary envelope which most likely is forming during the outflow of the optical star through the *L*<sup>2</sup> point.

Recently, several theoretical studies [47,48] showed that, due to the limited carrying capacity (for the outflowing material) of the surroundings of the *L*<sup>1</sup> point, the star (with the radiative envelope) during the mass transfer in the thermal time scale can be in a state of a significant overflow of its inner critical Roche lobe for a long time and the matter can outflow through the *L*<sup>2</sup> point.

So, taking into account that the radius of the occulting star in SS433 is significantly greater than the size of the inner critical Roche lobe, the estimate of the component mass ratio obtained by [39–41] should be considered as the lower limit of the mass ratio.

In the review by [7], the results of the interpretation of the eclipsing (orbital) and precessional variability in the hard X-ray range (*kT* ≈ 18–60 keV) were summarized. It was shown that, since the X-ray spectrum hardness does not change with the phase of the orbital and precessional variability (the observed X-ray flux changes by about five times), the hard X-ray radiation in contradiction to the soft X-ray radiation (*kT* ≈ 2–10 keV) is formed, not in narrow jets but in an extended quasi-isothermal hot corona of the accretion disc. According to [49], the temperature of the corona is ≈20 keV, its optical depth (Thompson

scattering) is <sup>≈</sup>0.2 and the mass loss rate in the jets is *<sup>M</sup>*˙ <sup>≈</sup> <sup>3</sup> <sup>×</sup> <sup>10</sup>−7*M* year−<sup>1</sup> (the mass loss rate in the wind from the disc is ~10−4*M* year<sup>−</sup>1).

In the paper by [7], to analyze the X-ray eclipses and precessional variability of SS433 in the hard X-ray range (*kT* = 18–60 keV), the following model was used: the optical star overflows its inner critical Roche lobe including the case of the filling of the outer Roche lobe, and the star outflows through the *L*<sup>2</sup> point. The interpretation of observations was made in the model of precessing accretion disc with hot corona in its central parts. It became clear that *q* cannot be unambiguously found from the X-ray eclipsing light curve alone. Co-interpretation of both the eclipsing X-ray light curve and precessional light curve in the model of the filling of the outer Roche lobe let us obtain the estimate of the mass ratio *<sup>q</sup>* <sup>=</sup> *<sup>M</sup>*<sup>x</sup> *M*v > 0.4–0.8. For the mass of the A7I optical star equal to 10*M*, it gives the estimate of the mass of the black hole *M*<sup>x</sup> ≈ 4–8*M*.

So, the conclusion about the high mass ratio of components in the SS433 system and the presence of the black hole in this system are supported by numerous independent observational data.

#### **6. Explanation of the Evolutionary Status of SS433 as a Semi-Detached Binary System**

According to the predictions of the modern theory of the evolution of massive close binary systems [50], at the stage of secondary mass exchange, due to the very high rate of the matter flow from an optical star to a relativistic object, a common envelope should form. In this case, due to deceleration by dynamic friction, rapid decreasing of the separation of components occurs. As a result, depending on the initial angular momentum of the binary system, either a short-period WR + c system (of the Cyg X-3 type, *P*orb 4.8 h) is formed or the Thorne–Zhitkov object [51], a completely convective supergiant with relativistic object in the center, is formed. According to [52], the life time for the Thorne–Zhitkov object can be very short due to the neutrino emission.

The SS433 system is at the stage of secondary mass exchange but its evolution, for some reason, took a different path. The system evolves as a semi-detached binary. The removal of mass and angular momentum from the system is not carried out into the common envelope but in the form of a radial outflow of a powerful wind (*M*˙ <sup>10</sup>−4*M* year<sup>−</sup>1, *<sup>v</sup>* <sup>10</sup><sup>3</sup> km s<sup>−</sup>1) from the supercritical accretion disc. This feature of SS433 has long remained a mystery.

Recently, studies have appeared that have shown that the ratio of the masses of the components is decisive in the evolutionary fate of a massive X-ray binary system at the stage of secondary mass exchange. In a recent work [53], it was shown that when a massive donor star in an X-ray binary system fills or overfills its Roche lobe and outflows onto a relativistic object, then if the component mass ratio is *<sup>q</sup>* <sup>=</sup> *Mx Mv* 0.29, the system can avoid the formation of a common envelope and remain semi-detached. Such a system evolves as a semi-detached system with a stable overflow of the Roche lobe with an optical star and with a stable flow of stellar matter on a relativistic object and with the formation of a supercritical accretion disc with a powerful outflow of stellar wind from it (the so-called isotropic re-emission mode or SS433-like mode). If the mass ratio *q* is small (*q* < 0.29), then a massive X-ray binary system at the stage of secondary mass exchange inevitably passes through the stage of evolution with a common envelope. The evolutionary scenario that described the formation of the object of SS433 type was considered in the paper by [54].

Thus, the central point in understanding the evolution of SS433 is knowing the real mass ratio of the system's components. Until recently, there was no complete clarity on this issue. Different methods gave estimates from *q* 0.14 to *q* = 0.8 (see reviews by [6,7] for all the necessary references). The new estimate *q* > 0.8 [26], obtained on the basis of our discovered increase in the orbital period of SS433, is very reliable and allows us to understand the reason why the SS433 system remains semi-detached. The SS433 system evolves as a semi-detached because the relativistic object here is a relatively large mass black hole (*Mx* > 8*M*), which is typical of stellar mass black holes in X-ray binaries.

#### **7. Conclusions**

Object SS433 ("Enigma of the century") is the first example of a microquasar discovered in the Galaxy. The object is a natural laboratory where complex and extremely interesting physical processes take place that shed light on many problems: from relativistic astrophysics to the physics and evolution of stars and high energy astrophysics.

For example, recently, the gamma radiation of SS433 was discovered; it can be formed as the result of the interaction of the relativistic jets with the matter of the W50 nebula that surrounds SS433 [55]. In the work by [56], the variability of the gamma flux from SS433, with the period of the precession of the accretion disc, was found. The authors noted that such variability can be related with the formation of the gamma radiation close to the accretion disc of SS433. In this regard, let us note that, in papers by [57,58], in the central part of the supercritical accretion disc in the SS433 system, the hot (*kT* ≈ 20 keV) extended corona radiating in the hard (*kT* = 18–60 keV) X-ray range was discovered. Recently, new observational data were obtained that indicated the similarity of physical processes in SS433 and in (still mysterious) objects such as ultra-luminous X-ray sources (ULXs) [59].

In this article, we have summarized the most fundamental problems associated with the object SS433 and showed how the latest observations and theoretical models allow us to understand the nature of this extremely interesting object.

**Funding:** The work was supported by the Russian Science Foundation grant 17-12-01241 and by the Scientific and Educational School of M. V. Lomonosov Moscow State University "Fundamental and applied space research". The author acknowledges support from the M. V. Lomonosov Moscow State University Program of Development.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


## *Article* **Multi-Component MHD Model of Hot Jupiter Envelopes**

**Andrey Zhilkin † and Dmitri Bisikalo \*,†**

Institute of Astronomy of the Russian Academy of Sciences, 48 Pyatnitskaya St., 119017 Moscow, Russia; zhilkin@inasan.ru

**\*** Correspondence: bisikalo@inasan.ru

† These authors contributed equally to this work.

**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.

**Keywords:** numerical simulation; magnetic hydrodynamics (MHD); hot Jupiters

#### **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–11]. These conclusions are confirmed by direct numerical calculations in the framework of one-dimensional aeronomic models [1,12–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].

**Citation:** Zhilkin, A.; Bisikalo, D. Multi-Component MHD Model of Hot Jupiter Envelopes. *Universe* **2021** , *7* , 422. https://doi.org/10.3390/ universe7110422

Academic Editors: Nazar R. Ikhsanov, Galina L. Klimchitskaya and Vladimir M. Mostepanenko

Received: 8 October 2021 Accepted: 3 November 2021 Published: 5 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In the three-dimensional approximation, the gas-dynamic structure of the extended envelope of hot Jupiters was studied in [17–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*<sup>1</sup> 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–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 <sup>μ</sup>J, where <sup>μ</sup><sup>J</sup> = 1.53 × <sup>10</sup><sup>30</sup> <sup>G</sup>·cm3 is the Jupiter magnetic moment. This value is in a good agreement with both observational [28–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–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–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 selfconsistent 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 multicomponent 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].

**Figure 1.** Schematic representation of the solar wind structure in the ecliptic plane. The Sun corresponds to a small colored circle in the center. The arrow shows the direction of rotation of the Sun. The boundary of the middle circle indicates the region of the corona, at the edge of which the magnetic field becomes purely radial. The shaded gray areas correspond to the zones of the heliospheric current sheet (shown by dotted lines running from the corona to the periphery), which separates the solar wind magnetic field with different directions of magnetic field lines (from the Sun or to the Sun). The orbit of a hot exoplanet is shown by a dotted circle, which is located in the heliospheric region.

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*, *θ*, *ϕ*). 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 *θ*, 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,\tag{1}$$

the equation of motion for the radial *vr* and azimuthal *v<sup>ϕ</sup>* components of the velocity vector *v*

$$v\_r \frac{dv\_r}{dr} - \frac{v\_\phi^2}{r} = -\frac{1}{\rho} \frac{dP}{dr} - \frac{GM\_\mathbf{s}}{r^2} - \frac{B\_\phi}{4\pi\rho r} \frac{d}{dr} (rB\_\phi)\_r \tag{2}$$

$$v\_r \frac{dv\_\varphi}{dr} + \frac{v\_r v\_\varphi}{r} = \frac{B\_r}{4\pi\rho r} \frac{d}{dr} \left(r B\_\varphi\right),\tag{3}$$

the equation of induction

$$\frac{1}{r}\frac{d}{dr}\left(rv\_rB\_\varphi - rv\_\varphi B\_r\right) = 0\tag{4}$$

and the Maxwell equation (∇ · *<sup>B</sup>* <sup>=</sup> 0)

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

Here, *ρ* represents the density, *P* is the pressure, *G* is the gravitational constant, and *M*<sup>s</sup> is the mass of the central star. Density, pressure, and temperature satisfy the equation of state for an ideal polytropic gas,

$$P = K\_{\text{k}} \rho^{\text{K}} = \frac{2k\_{\text{B}}}{m\_{\text{P}}} \rho T\_{\text{\textdegree}} \tag{6}$$

where *K<sup>κ</sup>* is the constant, *κ* is the polytropic index, *k*<sup>B</sup> is the Boltzmann constant, and *m*<sup>p</sup> 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\_\text{s} R\_\text{s}^2 \tag{7}$$

where *R*<sup>s</sup> is the radius of the star, and *B*<sup>s</sup> 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}\_{\rm sv} \tag{8}$$

where the integration constant *M*˙ <sup>s</sup> 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{\upsilon\_r^2}{2} + \frac{c\_s^2}{\kappa - 1} - \frac{GM\_\mathbf{s}}{r} + \frac{B\_\wp^2}{4\pi\wp}\right) - \frac{\upsilon\_\wp^2}{r} - rB\_\wp \frac{d}{dr}\left(\frac{B\_\wp}{4\pi\wp r}\right) = 0,\tag{9}$$

where the square of the speed of sound

$$
\sigma\_s^2 = \frac{\kappa P}{\rho} = \kappa K\_\kappa \rho^{\chi - 1}.\tag{10}
$$

We can multiply the Equation (3) by *vϕ* and divide by *vr*. Then, it can be represented as:

$$\frac{d}{dr}\left(\frac{v\_{\varphi}^{2}}{2} - \frac{B\_{r}B\_{\varphi}^{2}v\_{\varphi}}{4\pi\rho v\_{r}}\right) + \frac{v\_{\varphi}^{2}}{r} + rB\_{\varphi}\frac{d}{dr}\left(\frac{B\_{r}v\_{\varphi}}{4\pi\rho rv\_{r}}\right) = 0.\tag{11}$$

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

$$\begin{split} \frac{d}{dr} \left( \frac{v\_r^2}{2} + \frac{v\_\phi^2}{2} + \frac{c\_s^2}{\kappa - 1} - \frac{GM\_\mathbf{s}}{r} + \frac{B\_\phi^2}{4\pi\rho} - \frac{B\_r B\_\phi^2 v\_\phi}{4\pi\rho v\_r} \right) \\ + rB\_\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. \end{split} \tag{12}$$

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

$$\frac{d}{dr}\left(\frac{B\_r v\_\varphi}{4\pi\rho r v\_r} - \frac{B\_\varphi}{4\pi\rho r}\right) = \frac{1}{\dot{M}\_\text{s}}\frac{d}{dr}\left(r B\_r v\_\varphi - r B\_\varphi v\_r\right) = 0.\tag{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\_\wp^2}{2} + \frac{c\_s^2}{\kappa - 1} - \frac{GM\_\text{s}}{r} + \frac{B\_\wp^2}{4\pi\rho} - \frac{B\_r B\_\wp v\_\wp}{4\pi\rho v\_r} = Q\_{\text{S}\prime} \tag{14}$$

where the constant *Q*s determines the density of energy flux in the stellar wind, the value of which is *M*˙ <sup>s</sup>*Q*s.

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

$$\frac{B\_r}{4\pi\rho v\_r} = \frac{B\_\text{s}R\_\text{s}^2}{\dot{M}\_\text{s}} = \text{const.}\tag{15}$$

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

$$r v\_{\varphi} - \frac{B\_r}{4 \pi \rho v\_r} r B\_{\varphi} = L\_{\prime} \tag{16}$$

$$
\Gamma r \upsilon\_r B\_\varphi - r \upsilon\_\varphi B\_r = F.\tag{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θ*, since it is obvious that *F* = *cE<sup>θ</sup> 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}}, \quad a\_\phi = \frac{B\_\phi}{\sqrt{4\pi\rho}}.\tag{18}$$

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

$$
\left(v\_r^2 - a\_r^2\right)v\_\varphi = \frac{1}{r}\left(v\_r^2 L + \frac{a\_r F}{\sqrt{4\pi\rho}}\right),
\tag{19}
$$

$$a\left(v\_r^2 - a\_r^2\right)a\_\varphi = \frac{v\_r}{r}\left(a\_r L + \frac{F}{\sqrt{4\pi\rho}}\right). \tag{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 *B* must be collinear with the velocity vector *v* of the ideally conducting wind plasma: *v* <sup>×</sup> *<sup>B</sup>* = 0. For non-relativistic motion, the magnetic field *B* = *B*, and the velocity

$$v\_r' = v\_{r\prime} \quad v\_\varphi' = v\_\varphi - \Omega\_\text{s} r \tag{21}$$

where Ω<sup>s</sup> 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}.\tag{22}$$

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

$$\left(v\_r^2 - a\_r^2\right)v\_\varphi = \frac{1}{r}\left(v\_r^2L - a\_r^2\Omega\_8r^2\right),\tag{23}$$

$$\left(v\_r^2 - a\_r^2\right)a\_\varphi = \frac{\upsilon\_r a\_r}{r} \left(L - \Omega\_\mathrm{s} r^2\right). \tag{24}$$

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

$$|u\_A = |a\_r| = \frac{|B\_r|}{\sqrt{4\pi\rho}},\tag{25}$$

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

The azimuthal components of the velocity *v<sup>ϕ</sup>* and the magnetic field *B<sup>ϕ</sup>* 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\_{\rm s} r\_A^2. \tag{26}$$

As a result, we find the final solution

$$v\_{\varphi} = \frac{\Omega\_{\rm s}}{r} \frac{v\_r^2 r\_A^2 - a\_r^2 r^2}{v\_r^2 - a\_r^2} \, \text{} \tag{27}$$

$$a\_{\phi} = \frac{\Omega\_{\rm s}}{r} v\_r a\_r \frac{r\_A^2 - r^2}{v\_r^2 - a\_r^2}. \tag{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}.\tag{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},
\tag{30}
$$

where *vA* and *ρ<sup>A</sup>* 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\_{\varphi} = \Omega\_{\mathfrak{s}} r \frac{1 - \lambda^2 r\_A^2 / r^2}{1 - \lambda^2},\tag{31}$$

$$B\_{\varphi} = B\_{r} \frac{\Omega\_{\rm s} r}{v\_{r}} \lambda^{2} \frac{1 - r\_{A}^{2}/r^{2}}{1 - \lambda^{2}}.\tag{32}$$

Note that in the sub-Alfvén zone of the stellar wind *λ* < 1, while in the super-Alfvén zone *λ* > 1. At the Alfvén point, *λ* = 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\_\mathrm{S}^2\right)\left(v\_r^2 - u\_\mathrm{F}^2\right)\frac{r}{v\_r}\frac{dv\_r}{dr} = \left(v\_r^2 - a\_r^2\right)\left(2a\_s^2 + v\_\wp^2 - \frac{GM\_\mathrm{s}}{r}\right) + 2v\_r v\_\wp a\_r a\_\wp,\tag{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] \tag{34}$$

$$a^2 = a\_r^2 + a\_\varphi^2. \tag{35}$$

The quantities *uS* and *uF* 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* = *rS*, in which the speed *vr* = *uS*, and *fast magnetosonic r* = *rF*, in which the speed *vr* = *uF*. At these points, the coefficient for the derivative *dvr*/*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* = *rA*, in which the wind speed *vr* = *uA*, 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{split} \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{dv\_r}{dr} &= \left(v\_r^2 - u\_A^2\right)^3 \left(2v\_s^2 - \frac{GM\_s}{r}\right) \\ &+ \frac{\Omega\_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). \end{split} \tag{36}$$

Further, we have

$$\begin{split} \left(v\_r^2 - u\_S^2\right) \left(v\_r^2 - u\_F^2\right) &= v\_r^4 - \left(c\_s^2 + u\_A^2\right) v\_r^2 + c\_s^2 u\_A^2 - a\_\varphi^2 v\_r^2 \\ &= \left(v\_r^2 - c\_s^2\right) \left(v\_r^2 - u\_A^2\right) - a\_\varphi^2 v\_r^2. \end{split} \tag{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\_s^2}{r^2} v\_r^4 u\_A^2 (r\_A^2 - r^2)^2. \tag{38}$$

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

$$\begin{split} \left[ \left( v\_r^2 - c\_s^2 \right) \left( v\_r^2 - u\_A^2 \right)^3 - \frac{\Omega\_\mathbf{s}^2}{r^2} v\_r^4 u\_A^2 (r\_A^2 - r^2)^2 \right] \frac{r}{v\_r} \frac{dv\_r}{dr} \\ &= \left( v\_r^2 - u\_A^2 \right)^3 \left( 2v\_s^2 - \frac{GM\_\mathbf{s}}{r} \right) \\ &+ \frac{\Omega\_\mathbf{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 - 3 v\_r^2 u\_A^2 r^2 + v\_r^2 u\_A^2 r\_A^2 \right). \end{split} \tag{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, *dvr*/*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* → 0. Suppose that the radial velocity changes in this case according to the power law *vr* = *Arσ*, where *A* is a certain coefficient. Consider the case when the radial velocity *vr* is much smaller than any of the characteristic velocities *uS*, *uF* and *uA*. 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. \tag{40}
$$

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

$$v\_{\varphi} = \Omega\_{\sf s} r, \quad a\_{\varphi} = -a\_{r} \frac{\Omega\_{\sf s} r}{v\_{A}}.\tag{41}$$

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

$$
\mu\_A^2 \frac{GM\_\text{s}}{r} - 2v\_r u\_A^2 \frac{\Omega\_\text{s}^2 r^2}{v\_A} \approx \mu\_A^2 \frac{GM\_\text{s}}{r}.\tag{42}
$$

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

$$c\_s^2 \frac{r}{v\_r} \frac{dv\_r}{dr} = \frac{GM\_\text{s}}{r}.\tag{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}, \quad A = \frac{\dot{M}\_\text{s}}{4\pi} \left[ \frac{3 - 2\kappa}{\kappa - 1} \frac{\kappa K\_\text{s}}{GM\_\text{s}} \right]^{\frac{1}{\kappa - 1}}.\tag{44}
$$

Since in such a solution the index of *σ* must be positive, then the index of polytrope *κ* < 3/2. Thus, the solution with an accelerating wind (a positive value of the radial velocity gradient, *dvr*/*dr* > 0) is realized only in the case of polytropic index *κ* < 3/2. This value is less than the adiabatic index *γ* = 5/3. Therefore, effective sources of heating must be present in the stellar wind.

In fact, the entropy of an ideal gas *s* = *cV* ln(*P*/*ργ*), where *cV* is the specific heat capacity at a constant volume. Therefore, the full derivative

$$\frac{ds}{dt} = c\_V \left( \frac{d}{dt} \ln P - \gamma \frac{d}{dt} \ln \rho \right). \tag{45}$$

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

$$T\frac{ds}{dt} = -c\_V T(\gamma - \kappa) \frac{1}{\rho} \frac{d\rho}{dt}.\tag{46}$$

The expression on the right-hand side of this equation determines the heating function Γ. Since the specific internal energy of an ideal gas *ε* = *cV T*, then using the continuity Equation (1), we obtain:

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

We can see that in the case of an accelerating wind (*dvr*/*dr* > 0) and *κ* < *γ* 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 *ρ*0, radial velocity *v*<sup>0</sup> and temperature *T*<sup>0</sup> 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}}}\tag{48}$$

and the integration constant

$$
\dot{M}\_{\rm s} = 4\pi r\_0^2 \rho\_0 v\_0. \tag{49}
$$

The last expression follows from the law of mass conservation (8). The value of polytropic index *κ*, which lies in the range 1 < *κ* < 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}\_{\rm s} = 4\pi r\_A^2 \rho\_A v\_A. \tag{50}
$$

We find from here

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

where by

$$B\_0 = B\_\mathrm{s} \left(\frac{R\_\mathrm{s}}{r\_0}\right)^2\tag{52}$$

the value of the radial magnetic field at the point *r*<sup>0</sup> is indicated. We introduce dimensionless quantities *xA* = *rA*/*r*<sup>0</sup> and *yA* = *vA*/*v*0. It is not difficult to verify that they satisfy the relation

$$y\_A = \frac{\beta}{x\_A^2},$$
 
$$x\_A^2 \text{ }$$

where is the dimensionless parameter

$$
\beta = \frac{\rho\_0}{\rho\_A} = \frac{B\_0}{v\_0 \sqrt{4\pi\rho\_A}}.\tag{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 *ξ* = *r*/*rA* and *η* = *vr*/*vA*. In this case, the Alfvén Mach number (29) turns out to be equal to *λ* = *ξ* √*η*. In addition, we denote

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

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

$$\frac{\xi}{\eta} \frac{d\eta}{d\xi} = \frac{\chi}{X'} \tag{56}$$

where

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

$$\mathcal{Y} = (\lambda^2 - 1)^3 \Big( 2a\xi z^2 - \frac{\mu}{\varkappa\_A} \Big) \mathfrak{F} + \omega \mathfrak{x}\_A^2 (\lambda^2 - \xi^2) \Big( \lambda^4 - 3\lambda^2 \xi^2 + \lambda^2 + \xi^2 \Big), \tag{58}$$

$$
\mu = \frac{c\_0^2}{v\_0^2}, \quad \mu = \frac{GM\_\text{s}}{r\_0 v\_0^2}, \quad \omega = \frac{\Omega\_\text{s}^2 r\_0^2}{v\_0^2}. \tag{59}
$$

As was noted above, the slow magnetosonic (*ξ* = *ξS*, *η* = *ηS*) and fast magnetosonic (*ξ* = *ξF*, *η* = *ηF*) points are critical solution points. However, the Alfvén point (*ξ* = 1, *η* = 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{az^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] - \eta = 0,\tag{60}$$

where it is denoted

$$q = \frac{Q\_s}{v\_0^2}.\tag{61}$$

At the point *r*<sup>0</sup> we have *ξ*<sup>0</sup> = 1/*xA*, *η*<sup>0</sup> = 1/*yA*, *λ*<sup>0</sup> = 1/ *β*, *z*<sup>0</sup> = 1. Therefore, at this point, Equation (60) takes the form:

$$Z\_0 = \frac{1}{2} + \frac{\mathfrak{a}}{\kappa - 1} - \mu + \frac{\omega}{2(\lambda\_0^2 - 1)^2} \left[ \lambda\_0^4 \mathbf{x}\_A^4 + (2\lambda\_0^2 - 1)(1 - 2\mathbf{x}\_A^2) \right] - q = 0. \tag{62}$$

We will use *ξS*, *ηS*, *ξF*, *ηF*, *xA*, *q* and *κ* 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 *ξ<sup>S</sup>* and *ξF*, two Equation (60) written at the points *ξ<sup>S</sup>* and *ξ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 *ξS*, *ηS*, *ξF*, *ηF*, *xA*, *q* and *κ*. After these parameters are determined, for each value of *ξ*, we numerically solve the Equation (60) and find the dependence *η*(*ξ*) 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*<sup>0</sup> = 1400 cm<sup>−</sup>3, velocity *v*<sup>0</sup> = 130 km/s, and temperature *<sup>T</sup>*<sup>0</sup> <sup>=</sup> 7.3 <sup>×</sup> <sup>10</sup><sup>5</sup> <sup>K</sup> at a distance of *<sup>r</sup>*<sup>0</sup> <sup>=</sup> <sup>10</sup>*R* 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*<sup>s</sup> = 1.1*M*, and radius *R*<sup>s</sup> = 1.2*R*. The period of star proper rotation *<sup>P</sup>*rot = 14.4 days, which corresponds to the angular velocity <sup>Ω</sup><sup>s</sup> = 5.05 × <sup>10</sup>−<sup>6</sup> <sup>s</sup>−<sup>1</sup> or the linear velocity at the equator *v*rot = 4.2 km/s.

The major semi-axis of the planet orbit *A* = 10.2*R*, which is close to the selected boundary point *r*<sup>0</sup> and corresponds to the period of rotation around the star *P*orb = 84.6 h. The magnetic field at the stellar surface was set to *B*<sup>s</sup> = 0.5 G, which approximately corresponds to the value of average magnetic field at the surface of the Sun in a quiet state *<sup>B</sup>* <sup>=</sup> 1 G, if we compare the corresponding magnetic fluxes, *<sup>B</sup>*s*R*<sup>2</sup> <sup>s</sup> <sup>≈</sup> *<sup>B</sup>R*<sup>2</sup> .

It should be noted that the parameter *B*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*<sup>s</sup> is obtained from the average value of the radial magnetic field *B*<sup>c</sup> at the corona boundary *r* = *R*c, recalculated by considering the conservation of the magnetic flux by the star radius, *B*s*R*<sup>2</sup> <sup>s</sup> = *B*c*R*<sup>2</sup> c .

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 Ω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 *<sup>M</sup>*˙ <sup>s</sup> <sup>=</sup> 1.85 <sup>×</sup> 1011 g/s = 2.94 <sup>×</sup> <sup>10</sup>−<sup>15</sup> *<sup>M</sup>*/year, density at the Alfvén point *<sup>ρ</sup><sup>A</sup>* <sup>=</sup> 2.25 <sup>×</sup> <sup>10</sup>−<sup>22</sup> g/cm3, dimensionless parameters *α* = 0.713, *β* = 10.424, *μ* = 1.241, *ω* = 0.0725. Note that the magnitude of the mass loss rate of the star *M*˙ <sup>s</sup> 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 *ξ*, *η* 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.

**Figure 2.** Integral curves corresponding to various types of solutions for the magnetohydrodynamic stellar wind. Critical points are shown as circles. The Alfvén critical point has coordinates (1, 1). The bold line corresponds to the solution passing through all three critical points. The right diagram shows the vicinity of the Alfvén point.

The numerical solution of a system of non-linear algebraic equations describing the stellar wind structure gives the following coordinates of magnetosonic singular points: *ξ<sup>S</sup>* = 0.328, *η<sup>S</sup>* = 0.514, *ξ<sup>F</sup>* = 1.067, *η<sup>F</sup>* = 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 *xA* = 2.511 is obtained, which corresponds to the radius *rA* = 25.114*R*. The polytropic index turned out to be *κ* = 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* < 15*R*), the effective adiabatic index *γ* = 1.1 [65,66]. At large distances *r* > 25*R*, the effective adiabatic index can be estimated by the value *γ* = 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* < 20*R*, the wind structure in vicinity of the planet is found to be almost iso-thermal with good accuracy.

The calculation results are shown in Figures 3 and 4. Figure 3 shows the distributions of the radial velocity *vr* (solid line), the speed of sound *cs*, the Alfvén velocity *uA*, as well as the slow *uS* and fast *uF* 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 *vr* = *cs*. 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 *uS* ≈ min(*uA*, *cs*) = *cs*.

Figure 4 shows the radial profiles of the particle number density *n*(*r*), temperature *T*(*r*), azimuthal velocity *vϕ*(*r*), as well as the radial *Br*(*r*) (dotted line) and azimuthal *Bϕ*(*r*) (solid line) components of the magnetic field. The dependence *vϕ*(*r*) is non-monotonic. The maximum values of the azimuthal velocity are reached approximately in the area of the hot Jupiter's location 13*R*. However, the characteristic values of 15 km/s are small compared to the radial velocity of 130 km/s.

**Figure 3.** Distributions of the radial velocity *vr* (solid line), the speed of sound *cs*, the Alfvén velocity *uA*, as well as the slow *uS* and fast *uF* magnetosonic velocities depending on the radius *r*.

**Figure 4.** Particle number density profiles *n*(*r*) (**top left**), temperature profiles *T*(*r*) (**top right**), azimuthal velocity *vϕ*(*r*) (**bottom left**), as well as the radial *Br*(*r*) (dotted line) and azimuthal *Bϕ*(*r*) (solid line) components of the magnetic field (**bottom right**).

#### **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 *v* the average mass velocity of matter. Then, the continuity equations for each component of the kind *s* can be written as:

$$\frac{\partial \rho\_s}{\partial t} + \nabla \cdot (\rho\_s \mathbf{v}) = D\_s + S\_s. \tag{63}$$

where *ρ<sup>s</sup>* is the density of the component *s*, the value

$$D\_s = \nabla \cdot (\rho\_s \mathfrak{w}\_s) \tag{64}$$

determines diffusion, *<sup>w</sup><sup>s</sup>* <sup>=</sup> *<sup>v</sup>* <sup>−</sup> *<sup>v</sup><sup>s</sup>* 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\_{\mathbf{s}} \rho\_{\mathbf{s}} \tag{65}
$$

due to the conditions

$$
\sum\_{\mathbf{s}} \rho\_{\mathbf{s}} \mathfrak{w}\_{\mathbf{s}} = 0, \quad \sum\_{\mathbf{s}} \mathbb{S}\_{\mathbf{s}} = 0,\tag{66}
$$

we have the usual continuity equation,

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0. \tag{67}$$

Let us introduce the values *ξ<sup>s</sup>* that describe the mass fraction of the component *s* and satisfy the relations:

$$
\rho\_s = \zeta\_s \rho\_\prime \quad \sum\_s \zeta\_s = 1. \tag{68}
$$

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

$$\frac{\partial}{\partial t}(\rho\_{\rm s}^{\sf s}) + \nabla \cdot (\rho\_{\rm s}^{\sf s} \boldsymbol{\sigma}) = D\_{\rm s} + S\_{\rm s}.\tag{69}$$

Note that one of the values *ξ<sup>s</sup>* 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 *ξ<sup>s</sup>* can be considered independent. The mass fraction of electrons in this case is equal to

$$
\zeta\_{\rm e} = 1 - \sum\_{s} \zeta\_{s}.\tag{70}
$$

Since the mass of electrons is much smaller than the mass of ions and neutrals, *ξ*<sup>e</sup> 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[ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} \right] + \nabla P + \frac{\mathbf{B} \times \nabla \times \mathbf{B}}{4\pi} = \rho f - \mathbf{D}\_{\upsilon} \tag{71}
$$

$$
\rho \left[ \frac{\partial \varepsilon}{\partial t} + (\mathbf{v} \cdot \nabla) \varepsilon \right] + P \, \nabla \cdot \mathbf{v} = D\_\varepsilon + \rho Q\_\prime \tag{72}
$$

where *P* is the pressure, *B* is the magnetic field, *f* is the specific external force, the values

$$D\_{\upsilon} = \sum\_{\mathbf{s}} \nabla \cdot (\rho\_{\mathbf{s}} \mathbf{w}\_{\mathbf{s}} \mathbf{w}\_{\mathbf{s}})\_{\prime} \tag{73}$$

$$D\_{\varepsilon} = \sum\_{s} [\nabla \cdot (\rho\_s \varepsilon\_s \mathfrak{w}\_s) + P\_s \nabla \cdot \mathfrak{w}\_s] \tag{74}$$

are determined by the diffusion velocities *ws*, and the value of *Q* describes the heating– cooling sources. In the expression (74), the notations *Ps* and *ε<sup>s</sup>* 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,

$$\frac{\partial \mathbf{B}}{\partial t} - \nabla \times (\boldsymbol{\sigma} \times \mathbf{B}) = -\nabla \times (\boldsymbol{\eta} \, \nabla \times \mathbf{B}) + \nabla \times (\boldsymbol{\sigma}\_{\rm D} \times \mathbf{B}).\tag{75}$$

Here, the first term in the right hand side describes the ohmic diffusion of the magnetic field, where *η* is the corresponding magnetic viscosity. The second term determines the effect of ambipolar diffusion occurring in an incompletely ionized plasma. The vector *v*<sup>D</sup> 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 *η* = 0, *v*<sup>D</sup> = 0.

The density *ρ*, pressure *P*, internal energy *ε*, and temperature *T* satisfy the equation of state for an ideal gas

$$P = \frac{k\_{\rm B}}{\mu m\_{\rm P}} \rho T, \quad \varepsilon = \frac{k\_{\rm B} T}{\mu m\_{\rm P} (\gamma - 1)} \tag{76}$$

where *k*<sup>B</sup> is the Boltzmann constant, *m*<sup>p</sup> is the proton mass, *γ* is the adiabatic index. The average molecular weight *μ* is determined by the expression *ρ* = *μm*p*n*, where *n* is the total number density. Let us write down the condition of quasi-neutrality of plasma

$$
\hbar - e\mathfrak{n}\_{\text{e}} + \sum\_{\text{s}} e\_{\text{s}} \mathfrak{n}\_{\text{s}} = \mathfrak{0},
\tag{77}
$$

where *e* is the elementary charge, *n*<sup>e</sup> is the electron number density, *es* is the charge of particles of the kind *s*, and *ns* is the number density of the component *s*. Hence, the number density of electrons

$$
\mathfrak{m}\_{\mathfrak{e}} = \sum\_{\mathbf{s}} Z\_{\mathbf{s}} \mathfrak{n}\_{\mathbf{s}\prime} \tag{78}
$$

where *Zs* = *es*/*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\_{\rm P}} = \sum\_{\rm s} (1 + Z\_{\rm s}) \frac{\mathfrak{f}\_{\rm s}}{m\_{\rm s}} \, \, \, \tag{79}$$

where *ms* 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 *B* is represented as a superposition of the background magnetic field *H* and the magnetic field *b* induced by electrical currents in plasma itself, *B* = *H* + *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, ∇ × *<sup>H</sup>* <sup>=</sup> 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, *<sup>∂</sup>H*/*∂<sup>t</sup>* <sup>=</sup> 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, *v<sup>s</sup>* = *v*. We can also assume that all the components are in thermodynamic equilibrium, and therefore their temperatures are equal to the temperature of matter, *Ts* = *T*.

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

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0,\tag{80}$$

$$\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} = -\frac{\nabla P}{\rho} - \frac{\mathbf{b} \times \nabla \times \mathbf{b}}{4\pi \rho} - \frac{\mathbf{H} \times \nabla \times \mathbf{b}}{4\pi \rho} + f\_{\prime} \tag{81}$$

$$\frac{\partial \mathbf{b}}{\partial t} = \nabla \times (\mathbf{v} \times \mathbf{b} + \mathbf{v} \times \mathbf{H}) - \frac{\partial H}{\partial t},\tag{82}$$

$$\frac{\partial \varepsilon}{\partial t} + (\boldsymbol{\sigma} \cdot \nabla)\varepsilon + \frac{P}{\rho} \nabla \cdot \boldsymbol{\sigma} = Q\_{\prime} \tag{83}$$

$$\frac{\partial}{\partial t}(\rho \mathfrak{f}\_s^\star) + \nabla \cdot (\rho \mathfrak{f}\_s^\star \mathfrak{v}) = \mathbb{S}\_\mathbf{s}, \quad \mathbf{s} = 1, \dots, \mathbf{N}. \tag{84}$$

To close this system of equations, the equation of state for an ideal gas (76) with the adiabatic index *γ* = 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 *Ss* 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 **Ω**, the locations of the star and the planet centers do not change. In this case, the specific external force is determined by the expression

$$f = -\nabla\Phi - 2(\Omega \times \mathbf{v}).\tag{85}$$

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

$$\Phi = -\frac{GM\_\text{s}}{|\mathbf{r} - \mathbf{r\_s}|} - \frac{GM\_\text{p}}{|\mathbf{r} - \mathbf{r\_P}|} - \frac{1}{2} [\mathbf{\Omega} \times (\mathbf{r} - \mathbf{r\_c})]^2 \tag{86}$$

where *M*<sup>s</sup> is the mass of the star, *M*<sup>p</sup> is the mass of the planet, *r*<sup>s</sup> is the radius vector of the star center, *r*<sup>p</sup> is the radius vector of the planet center, *r*<sup>c</sup> 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 *vr*(*r*) obtained from the wind model (see Section 2). Profiles for the remaining values are calculated using this profile: *ρ*(*r*), *vϕ*(*r*), *Bϕ*(*r*). This allows finding the wind parameters at any point *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

$$
\boldsymbol{\sigma}\_{\rm W} = \boldsymbol{\mu}\_{\rm W} - \boldsymbol{\Omega} \times (\boldsymbol{r} - \boldsymbol{r}\_{\rm c}) \tag{87}
$$

where *u*<sup>w</sup> 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 *γ*. 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\_{\rm W} = (\gamma - \kappa) \varepsilon\_{\rm W} \nabla \cdot \boldsymbol{\mathfrak{v}}\_{\rm W} \tag{88}$$

where *ε*<sup>w</sup> is the specific internal energy of the wind matter. This takes into account the expression (87) wind velocity divergence

$$
\nabla \cdot \boldsymbol{v}\_{\rm W} = \nabla \cdot \boldsymbol{u}\_{\rm W} = \frac{1}{r^2} \frac{d}{dr} \left( r^2 v\_r \right). \tag{89}
$$

Since the stellar wind accelerates (the velocity *vr* increases with the distance *r* from the star), the divergence of its velocity turns out to be positive. Consequently, the function *Q*<sup>w</sup> > 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\_\text{a} \exp\left[-\frac{GM\_\text{p}}{R\_\text{p}A\_\text{gas}T\_\text{a}}\left(1 - \frac{R\_\text{p}}{|r - r\_\text{P}|}\right)\right].\tag{90}$$

Here, *ρ*<sup>a</sup> is the density at the photometric radius of hot Jupiter *R*p, *T*<sup>a</sup> is the atmospheric temperature, *A*gas = *k*B/(*μm*p) is the gas constant, *μ* 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{GM\_{\rm P}}{R\_{\rm P}A\_{\rm gas}T\_{\rm a}} = 10.3 \,\mu \left(\frac{T\_{\rm a}}{10^4 \,\text{K}}\right)^{-1}.\tag{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\_{\rm tot} = P\_{\rm W} + \rho\_{\rm W} v\_{\rm W}^2 + \frac{B\_{\rm W}^2}{8\pi}. \tag{92}$$

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

$$R\_{\rm a} = \left(1 - \frac{1}{\eta} \ln \frac{P\_{\rm a}}{P\_{\rm tot}}\right)^{-1} R\_{\rm P}.\tag{93}$$

The value *R*<sup>a</sup> is determined to a greater extent by the atmosphere parameters *ρ*<sup>a</sup> and *T*<sup>a</sup> and to a lesser extent by the magnetic field of the wind *B*w. In particular, the ratio between the initial radius of the atmosphere *R*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 *χ* = [He/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 H2, 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:

$$\mathbf{X} = \frac{\rho(\mathbf{H}\_2) + \rho(\mathbf{H}) + \rho(\mathbf{H}^+)}{\rho} = \mathfrak{J}(\mathbf{H}\_2) + \mathfrak{J}(\mathbf{H}) + \mathfrak{J}(\mathbf{H}^+), \tag{94}$$

$$Y = \frac{\rho(\text{He}) + \rho(\text{He}^+)}{\rho} = \tilde{\mathfrak{J}}(\text{He}) + \tilde{\mathfrak{J}}(\text{He}^+). \tag{95}$$

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

$$X = \frac{1}{1 + 4\chi'} \quad Y = \frac{4\chi}{1 + 4\chi}.\tag{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,

$$\mathbf{x}\_1 = \frac{\rho(\mathbf{H}^+)}{\rho(\mathbf{H}^+) + \rho(\mathbf{H})}, \quad \mathbf{x}\_2 = \frac{\rho(\mathbf{H})}{\rho(\mathbf{H}) + \rho(\mathbf{H}\_2)}, \quad \mathbf{x}\_3 = \frac{\rho(\mathbf{H} \mathbf{e}^+)}{\rho(\mathbf{H} \mathbf{e}^+) + \rho(\mathbf{H} \mathbf{e})}. \tag{97}$$

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

$$
\zeta(\text{H}\_2) = \frac{(1 - \chi\_1)(1 - \chi\_2)}{1 - \chi\_1 + \chi\_2 \chi\_2} \text{X}\_{\prime} \tag{98}
$$

$$\zeta(\text{H}) = \frac{(1 - \chi\_1)\chi\_2}{1 - \chi\_1 + \chi\_2\chi\_2}X,\tag{99}$$

$$
\xi(\mathcal{H}^+) = \frac{\mathbb{x}\_1 \mathbb{x}\_2}{1 - \mathbb{x}\_1 + \mathbb{x}\_2 \mathbb{x}\_2} \mathcal{X}\_\prime \tag{100}
$$

$$\mathfrak{F}(\text{He}) = (1 - \mathfrak{x}\_{\mathfrak{Z}})\mathcal{Y}\_{\prime} \tag{101}$$

$$
\xi(\mathrm{He}^+) = \mathfrak{x}\_3 \mathrm{Y}.\tag{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{\varkappa\_2 + 2\varkappa\_1\varkappa\_2}{1 + \varkappa\_1 - \varkappa\_1\varkappa\_2}\right) + \frac{1}{4}Y(1 + \varkappa\_3). \tag{103}$$

The total degree of ionization of atmospheric matter

$$\mathfrak{F} = \mathfrak{f}(\mathrm{H}^{+}) + \mathfrak{f}(\mathrm{He}^{+}) = \frac{\mathfrak{x}\_{1}\mathfrak{x}\_{2}}{1 - \mathfrak{x}\_{1} + \mathfrak{x}\_{1}\mathfrak{x}\_{2}}\mathrm{X} + \mathfrak{x}\_{3}\mathrm{Y}.\tag{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: *χ* = [He/H] = 0.05 [25], *x*<sup>1</sup> = *x*<sup>2</sup> = *x*<sup>3</sup> = 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 *ξ*(H2) = 0.009, *ξ*(H) = 0.08, *ξ*(H+) = 0.74, *ξ*(He) = 0.02, *ξ*(He+) = 0.15. At the same time, the average molecular weight of the atmospheric matter is *μ* = 0.69, and the total degree of ionization is *ξ* = 0.89.

The boundary conditions (density *ρ*a, temperature *T*<sup>a</sup> and velocity *v*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 *<sup>M</sup>*˙ <sup>a</sup> <sup>≈</sup> 109–1010 g/s. Taking into account the expression *M*˙ <sup>a</sup> = 4*πR*<sup>2</sup> <sup>p</sup>*ρ*a*v*<sup>a</sup> from here, we can find the speed of atmospheric outflow *v*a at the photometric radius.

#### *4.3. Magnetic Field*

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

$$H = H\_\text{s} + H\_\text{W} + H\_\text{P} + H\_\text{a.} \tag{105}$$

Here, *H*<sup>s</sup> describes the proper field of the star, *H*<sup>w</sup> is determined by the wind field, *H*<sup>p</sup> is the proper field of the planet, and *H*<sup>a</sup> 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 *H*<sup>s</sup> 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 *B*<sup>w</sup> does not satisfy the condition of potentiality, ∇ × *<sup>B</sup>*<sup>w</sup> <sup>=</sup> 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 *BR* 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ϕ*/*br* is small in absolute value (see Figure 4). Taking into account the assumed spherical symmetry of the stellar wind, we obtain (see Section 2)

$$H\_{\rm W} = \frac{B\_{\rm s} R\_{\rm s}^2}{|r - r\_{\rm s}|^2} \mathfrak{n}\_{\rm s} \tag{106}$$

where the unit vector *<sup>n</sup>*<sup>s</sup> = (*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*s)/|*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*s<sup>|</sup> determines the direction from the center of star *r*<sup>s</sup> to the observation point *r*.

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

$$H\_{\rm P} = \frac{\mu\_{\rm P}}{|r - r\_{\rm P}|^3} \left[ \Im(d\_{\rm P} \cdot n\_{\rm P}) n\_{\rm P} - d\_{\rm P} \right],\tag{107}$$

where *<sup>μ</sup>*<sup>p</sup> is the magnetic moment of the planet, *<sup>n</sup>*<sup>p</sup> = (*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*p)/|*<sup>r</sup>* <sup>−</sup> *<sup>r</sup>*p|, *<sup>d</sup>*<sup>p</sup> is a unit vector directed along the magnetic axis and determining the vector of magnetic moment *<sup>μ</sup>*<sup>p</sup> <sup>=</sup> *<sup>μ</sup>*p*d*p. In our calculations, we assumed the value of magnetic moment of the hot Jupiter HD 209458b to be *μ*<sup>p</sup> = 0.1*μ*J. The orientation of magnetic dipole was determined by the angles *θ* and *φ*, which were used as model parameters. The components of the unit vector *d*<sup>p</sup> directed along the magnetic axis in the Cartesian coordinate system are described by the expressions:

$$d\_{\rm P} = (\sin \theta \cos \phi, \quad \sin \theta \sin \phi, \quad \cos \theta). \tag{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

$$H\_{\rm a} = \frac{\mu\_{\rm a}}{|r - r\_{\rm P}|^3} \left[ \Im(\mathbf{d}\_{\rm a} \cdot \boldsymbol{\nu}\_{\rm P}) \boldsymbol{\nu}\_{\rm P} - d\_{\rm a} \right],\tag{109}$$

where the corresponding magnetic moment *μ*<sup>a</sup> = *R*<sup>3</sup> <sup>a</sup>*B*w/2, and the unit vector *<sup>d</sup>*<sup>a</sup> <sup>=</sup> <sup>−</sup>*B*w/*B*<sup>w</sup> is directed against the vector *B*w. The magnitude of induced magnetic moment *μ*<sup>a</sup> depends on the initial radius of the atmosphere *R*<sup>a</sup> and the wind field *B*<sup>w</sup> in the orbit of the planet. If we take an average field of the star *B*<sup>s</sup> = 0.5 Gs, then, for the characteristic dimensions of the atmosphere *R*<sup>a</sup> = (4–5)*R*p, we can obtain the values of magnetic moment *μ*<sup>a</sup> = (0.05–0.2)*μ*J. In other words, in order of magnitude, the induced magnetic moment *μ*<sup>a</sup> turns out to be equal to the intrinsic magnetic moment *μ*<sup>p</sup> of the hot Jupiter. Note that the induced magnetic field *H*<sup>a</sup> is completely determined by the wind field *B*w. Therefore, the structure of induced field (in particular, the vector of induced magnetic moment *<sup>μ</sup>*<sup>a</sup> <sup>=</sup> *<sup>μ</sup>*a*d*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 *H*<sup>p</sup> and *H*<sup>a</sup> are stationary, *∂H*p/*∂t* = 0, *∂H*a/*∂t* = 0. The proper field of a star, on the contrary, is non-stationary, *<sup>∂</sup>H*s/*∂<sup>t</sup>* <sup>=</sup> 0, since it rotates together with the star at the angular velocity **<sup>Ω</sup>**<sup>s</sup> <sup>−</sup> **<sup>Ω</sup>**. The radial component of the wind magnetic field *H*<sup>w</sup> 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 noninertial reference frame is equal to

$$
\sigma\_{\rm s} = -\boldsymbol{\Omega} \times (\boldsymbol{r} - \boldsymbol{r}\_{\rm s}) \tag{110}
$$

since the azimuthal component of the field, due to the proper rotation of the star **Ω**s, has already been taken into account in the vector *b*. The change in the magnetic field *H*<sup>w</sup> in time is determined by the equation

$$\frac{\partial H\_{\rm W}}{\partial t} = \nabla \times (\boldsymbol{v}\_{\rm s} \times \boldsymbol{H}\_{\rm W}).\tag{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 *H*w. Consequently, the left hand side should also be equal to zero, *∂H*w/*∂t* = 0. Thus, in our model, the total background magnetic field (105) turns out to be stationary, *∂H*/*∂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* + Δ*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 *b* of plasma without considering the background magnetic field *H*:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0,\tag{112}$$

$$\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\frac{\nabla P}{\rho} - \frac{\mathbf{b} \times \nabla \times \mathbf{b}}{4\pi\rho} + f\_{\prime} \tag{113}$$

$$\frac{\partial \mathbf{b}}{\partial t} = \nabla \times (\mathbf{v} \times \mathbf{b}),\tag{114}$$

$$\frac{\partial \varepsilon}{\partial t} + (\boldsymbol{\sigma} \cdot \nabla)\varepsilon + \frac{P}{\rho} \nabla \cdot \boldsymbol{\sigma} = Q\_{\prime} \tag{115}$$

$$\frac{\partial}{\partial t}(\rho \mathfrak{f}\_s) + \nabla \cdot (\rho \mathfrak{f}\_s \mathfrak{v}) = 0, \quad s = 1, \ldots, N. \tag{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:

$$\frac{\partial \mathbf{v}}{\partial t} = -\frac{H \times \nabla \times \mathbf{b}}{4\pi\rho},\tag{117}$$

$$\frac{\partial \mathbf{b}}{\partial t} = \nabla \times (\mathbf{v} \times \mathbf{H}).\tag{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 *ρ* and the specific internal energy *ε* 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 *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*<sup>p</sup> = 0.71*M*<sup>J</sup> and the photometric radius *R*<sup>p</sup> = 1.38*R*J, where *M*<sup>J</sup> and *R*<sup>J</sup> is the mass and radius of Jupiter. The host star is characterized by the following parameters: spectral class G0, mass *M*<sup>s</sup> = 1.1*M*, radius *R*<sup>s</sup> = 1.2*R*. The period of proper rotation of the star *<sup>P</sup>*rot = 14.4 days, which corresponds to the angular velocity <sup>Ω</sup><sup>s</sup> = 5.05 × <sup>10</sup>−<sup>6</sup> <sup>s</sup>−<sup>1</sup> or linear velocity at the equator *v*rot = 4.2 km/s. The major semi-axis of the planet orbit *A* = 10.2*R*, which corresponds to the period of revolution around the star *P*orb = 84.6 h.

In the simulations, the temperature of atmosphere *T*<sup>a</sup> is varied, while the particle number density at the photometric radius was set equal to a fixed value *n*<sup>a</sup> = 1010 cm−3. The atmospheric outflow rate *M*˙ <sup>a</sup> = 109 g/s corresponds to the velocity of the planetary wind at the photometric radius *v*<sup>a</sup> = 78 cm/s. The values of these parameters correspond to the results obtained in aeronomic models for atmospheres of hot Jupiters [12,14–16,48,73]. We assumed that the magnitude of the magnetic moment *μ*<sup>p</sup> 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 *θ* = 90◦ and *ϕ* = 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*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 *vr* (left panel) and the heating function *Q*w, which is determined by the Equation (88) (right panel), for all three cases of values *B*s. The characteristic value of the heating function in the vicinity of the planet *<sup>Q</sup>*<sup>w</sup> = <sup>5</sup> × 109 erg·g−1·s<sup>−</sup>1. In this case, the maximum values of the heating function are reached at distances of approximately 4*R* from the star center.

**Figure 5.** Distributions of the radial wind velocity *vr* (**left**) and the heating function *Q*w (88) (**right**) for the magnetic field values at the stellar surface *B*s used in simulations.

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 *r*<sup>p</sup> = (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 *<sup>r</sup>*<sup>s</sup> = (−*A*, 0, 0). The *<sup>y</sup>* 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 **Ω**.

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* = −2*A*). Calculation results for cases of weak (*B*<sup>s</sup> = 0.01 G) and strong (*B*<sup>s</sup> = 0.5 G) magnetic field of the wind are shown in Figures 6 and 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 *<sup>ρ</sup>*<sup>w</sup> <sup>=</sup> 2.3 <sup>×</sup> <sup>10</sup>−<sup>21</sup> g/cm3, corresponding to the wind density at a distance of 10*R* 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.

**Figure 6.** Distributions of density (color), velocity (arrows) and magnetic field (lines) in the stellar wind for the case when the magnetic field at the surface of the star *B*<sup>s</sup> = 0.01 G at the initial time (**left**) and after one orbital period (**right**). The dotted line shows the boundary of the Roche lobe.

**Figure 7.** The same as in Figure 6, but for the case, when the magnetic field at the stellar surface *B*<sup>s</sup> = 0.5 G.

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*<sup>a</sup> = 6000 K. The left panel corresponds to the magnitude of the magnetic field at the surface of the star *B*<sup>s</sup> = 0.01 G (weak field), and the right panel does to the case for the strong field *B*<sup>s</sup> = 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*s leads to an increase in the total wind pressure (92) and, consequently, to a decrease in the initial radius of the atmosphere (93).

**Figure 8.** The initial distributions of density (color) and magnetic field (lines with arrows) for the case, when the temperature of the atmosphere *T*<sup>a</sup> = 6000 K, and the magnetic field at the stellar surface *B*s is 0.01 G (left) and 0.5 G (right). The dotted line shows the boundary of the Roche lobe. The white circle corresponds to the photometric radius of the planet.

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 *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.


**Table 1.** Parameters and characteristics of the models: *T*a is the temperature of the upper atmosphere, *B*<sup>s</sup> is the magnetic field of the star, and *M*˙ <sup>p</sup> is the mass loss rate of the planet.

#### *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*<sup>s</sup> = 0.01 G. The Alfvén Mach number (29) in the orbit of the planet is *λ* = 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, *λ* = 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*a (see Table 1). Namely, the temperature of atmosphere *T*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 ≤ *x*/*R*<sup>p</sup> ≤ 30, −30 ≤ *y*/*R*<sup>p</sup> ≤ 30, −15 ≤ *z*/*R*<sup>p</sup> ≤ 15 with the number of cells 192 × 192 × 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 *ρ*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–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.

**Figure 9.** Distributions of density (color, iso-lines), velocity (arrows) and magnetic field (lines) in the orbital plane at a time moment approximately equal to a one third of the orbital period for models 1 (*T*<sup>a</sup> = 5000 K, top left), 2 (*T*<sup>a</sup> = 5500 K, top right), 3 (*T*<sup>a</sup> = 6000 K, bottom left) and 4 (*T*<sup>a</sup> = 6500 K, bottom right). The density is expressed in units of *ρ*w. The dotted line shows the border of the Roche lobe. The white circle corresponds to the photometric radius of the planet.

In model 1 (*T*<sup>a</sup> = 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*<sup>a</sup> = 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 *<sup>M</sup>*˙ <sup>p</sup> <sup>≈</sup> 109 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*<sup>a</sup> = 6000 K, lower left panel Figure 9) and model 4 (*T*<sup>a</sup> = 6500 K, lower right panel Figure 9). In these models, two powerful flows are formed from the vicinity of Lagrange points *L*<sup>1</sup> 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 *<sup>M</sup>*˙ <sup>p</sup> <sup>≈</sup> <sup>4</sup> <sup>×</sup> 109 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 *<sup>M</sup>*˙ <sup>p</sup> <sup>≈</sup> 1010 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*<sup>a</sup> = 6500 K, which corresponds in the case of a weak field (*B*<sup>s</sup> = 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*<sup>s</sup> = 0.2 G (model 5) and *B*<sup>s</sup> = 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.

**Figure 10.** Distributions of density (color, iso-lines), velocity (arrows), and magnetic field (lines) in the orbital plane at a time moment approximately equal to one third of the orbital period for Model 5 (*B*<sup>s</sup> = 0.2 G, **left**) and Model 6 (*B*<sup>s</sup> = 0.5 G, **right**). The density is expressed in units of *ρ*w. The dotted line shows the border of the Roche lobe. The white circle corresponds to the photometric radius of the planet.

In both models, the planet is located in the sub-Alfvén wind zone, but between the points *rS* (slow magnetosonic) and *rA* (Alfvén). Consequently, in the orbit of the planet, the wind velocity *vw* exceeds the slow magnetosonic velocity *uS*, but becomes less than the Alfvén velocity *uA*. In model 5, the Alfvén Mach number at the planet orbit is *λ* = 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 *λ* = 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 *λ* = 0.31, and the total Alfvén Mach number is *λ* = 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*<sup>1</sup> at the day side weakened and the envelope became closer to the quasi-open type. We obtained that the mass loss rate *<sup>M</sup>*˙ <sup>p</sup> <sup>≈</sup> <sup>9</sup> <sup>×</sup> <sup>10</sup><sup>9</sup> g/s for the model 5 and *<sup>M</sup>*˙ <sup>p</sup> <sup>≈</sup> <sup>7</sup> <sup>×</sup> 109 g/s for the model 6. It follows from these results that, at fixed parameters of the atmosphere (density *ρ*<sup>a</sup> and temperature *T*a), the rate of mass loss *M*˙ <sup>p</sup> for hot Jupiter decreases with an increase in the magnetic field of the wind (the parameter *B*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 quasiopen 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*<sup>1</sup> 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.

**Author Contributions:** A.Z. and D.B. developed the analytical and numerical model, obtained the results and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors are grateful to the Government of the Russian Federation and the Ministry of Higher Education and Science of the Russian Federation for the support (grant 075-15-2020-780 (no. 13.1902.21.0039)).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The study was carried out using computing facilities of the Interdepartmental Supercomputer Center of the Russian Academy of Sciences, www.jscc.ru, accessed on 20 October 2021.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Difference Scheme for the Equations of Multi-Component Magnetic Hydrodynamics**

#### *Appendix A.1. Roe Matrix*

In this section, we describe the adaptation of the Roe scheme [77] to the case of multi-component magnetic hydrodynamics equations. We will discard all source terms, since they can be accounted separately. The hyperbolic part of the system consists of the equations of magnetic hydrodynamics

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0,\tag{A1}$$

$$\frac{\partial}{\partial t}(\rho v) + \nabla \cdot (\rho v v + P\_{\Gamma} + \mathcal{B}\mathcal{B}) = 0,\tag{A2}$$

$$\frac{\partial \mathbf{B}}{\partial t} - \nabla \times (v \times \mathbf{B}) = 0,\tag{A3}$$

$$\frac{\partial E\_{\rm T}}{\partial t} + \nabla \cdot \left[ v(E\_{\rm T} + P\_{\rm T}) - \mathcal{B}(v \cdot \mathcal{B}) \right] = 0 \tag{A4}$$

and equations for the mass fractions of components

$$\frac{\partial}{\partial t}(\rho\_{\varsigma s}^{\sf x}) + \nabla \cdot (\rho\_{\varsigma s}^{\sf x} \boldsymbol{\sigma}) = 0, \quad s = 1, \ldots, N. \tag{A5}$$

Here, all the equations are written in a conservative form and notations for total pressure and total energy density are used

$$P\_{\Gamma} = P + \frac{\mathbf{B}^2}{2}, \quad E\_{\Gamma} = \rho \varepsilon + \rho \frac{\mathbf{v}^2}{2} + \frac{\mathbf{B}^2}{2}. \tag{A6}$$

For the convenience of numerical simulation, a system of units is used in these equations, in which the multiplier 4*π* does not occur in the expression for the electromagnetic force.

For the numerical solving of three-dimensional equations of multi-component magnetic hydrodynamics (A1)–(A5), the splitting technique by spatial directions can be used. As a result, the solving of a three-dimensional problem is reduced to solving a series of one-dimensional problems. In the case of using Godunov-type schemes [78], numerical fluxes in each spatial direction are calculated based on the corresponding one-dimensional Riemann problem on the decay of an arbitrary discontinuity.

Consider the case of a plane flow corresponding to the coordinate direction *x*. Since in the plane flow the magnetic field component *BX* = const, the equations of one-dimensional multi-component magnetic hydrodynamics in a conservative form can be written as

$$\frac{\partial \mathbf{u}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} = \mathbf{0},\tag{A7}$$

where the vectors of conservative variables and fluxes are defined by expressions:

$$\begin{array}{c} \begin{pmatrix} \rho\\ \rho v\_{\mathcal{X}}\\ \rho v\_{\mathcal{Y}}\\ \rho v\_{z}\\ B\_{\mathcal{Y}}\\ B\_{z}\\ E\_{\mathcal{T}}\\ \rho\_{\mathcal{S}s}^{\mathcal{X}} \end{pmatrix}, \quad \mathbf{F} = \begin{pmatrix} \rho v\_{\mathcal{X}}\\ \rho v\_{\mathcal{X}}^{2} + P\_{\mathcal{T}}\\ \rho v\_{\mathcal{X}} v\_{\mathcal{Y}} - B\_{\mathcal{X}} B\_{\mathcal{Y}}\\ \rho v\_{\mathcal{X}} v\_{\mathcal{Z}} - B\_{\mathcal{X}} B\_{\mathcal{Z}}\\ v\_{\mathcal{X}} B\_{\mathcal{Y}} - v\_{\mathcal{Y}} B\_{\mathcal{X}}\\ v\_{\mathcal{X}} B\_{\mathcal{Z}} - v\_{\mathcal{Z}} B\_{\mathcal{X}}\\ \rho h v\_{\mathcal{X}} - B\_{\mathcal{X}} (\mathbf{v} \cdot \mathbf{B}),\\ \rho v\_{\mathcal{X}} \mathfrak{F}\_{\mathcal{S}s} \end{pmatrix}. \end{array} \tag{A8}$$

Here, it is assumed that the index *s* runs through values from 1 to *N* and a notation is used for the total enthalpy density *h*, determined by the relation

$$
\rho \mathbf{h} = E\_{\mathbb{T}} + P\_{\mathbb{T}}.\tag{A9}
$$

Note that the equation for the *Bx* component is excluded from this system. In fact, this value is a flow parameter.

The Roe scheme refers to Godunov-type schemes and is based on an approximate solving of the Riemann problem on the decay of an arbitrary discontinuity. In this method, instead of solving the Riemann problem for the original system of non-linear Equation (A7), a linearized problem is solved

$$
\frac{
\partial \mathbf{u}
}{
\partial t
} + \hat{A}(\mathbf{u}\_L, \mathbf{u}\_R) \cdot \frac{
\partial \mathbf{u}
}{
\partial \mathbf{x}
} = 0
\tag{A10}
$$

with initial conditions: *u*(*x*, 0) = *u<sup>L</sup>* for *x* < 0 and *u*(*x*, 0) = *u<sup>R</sup>* for *x* > 0.

In order for the solutions of the original (A7) and the linearized (A10) problems to be consistent, the matrix *A*ˆ(*uL*, *uR*) must satisfy the following three conditions.


$$
\hat{A}(\mathbf{u}\_{L\prime}\mathbf{u}\_{R}) \cdot \Delta \mathbf{u} = \Delta \mathbf{F}\_{\prime} \tag{A11}
$$

where is denoted *<sup>δ</sup><sup>u</sup>* <sup>=</sup> *<sup>u</sup><sup>R</sup>* <sup>−</sup> *<sup>u</sup>L*, *<sup>δ</sup><sup>F</sup>* <sup>=</sup> *<sup>F</sup><sup>R</sup>* <sup>−</sup> *<sup>F</sup>L*. In this case, solving of the linearized discontinuity decay problem (A10) will satisfy the same integral conservation laws as the solving of the original non-linear problem (A7).

As is known, the solution of the Riemann problem for a linear hyperbolic system of Equation (A10) is a set of strong discontinuities whose velocities are equal to the eigenvalues *λα* of the Roe matrix *A*ˆ(*uL*, *uR*), where *α* is the index of the characteristic. Corresponding jumps of values at discontinuities

[**u**]*<sup>α</sup>* = **r***α*Δ*Sα*, (A12)

where square brackets denote differencies between the right and the left values when crossing the discontinuity, Δ*S<sup>α</sup>* = *l <sup>α</sup>* · *<sup>δ</sup><sup>u</sup>* is the characteristic amplitudes, and *<sup>r</sup>α*, *<sup>l</sup> <sup>α</sup>* is the right and left eigenvectors of Roe matrix. At each discontinuity, the corresponding Hugoniot conditions must be satisfied:

$$
\lambda\_{\mathfrak{a}}[\mathbf{u}]\_{\mathfrak{a}} = [\mathbf{F}]\_{\mathfrak{a}}.\tag{A13}
$$

The Roe matrix for the equations of magnetic hydrodynamics (A1)–(A4), as well as all the characteristic parameters necessary for constructing the scheme for the special case of the adiabatic index *γ* = 2, are described in [79]. In [80], the corresponding Roe matrix with the adiabatic index 1 < *γ* ≤ 2 was constructed for the first time. A detailed description of this scheme is given in our monograph [68]. We will not dwell on this here, considering all these properties to be known. Recall, however, that this matrix of dimension 7 × 7 has the following set of eigenvalues:

$$
\lambda\_{\pm \mathbf{F}} = v\_{\mathbf{x}} \pm u\_{\mathbf{F}\prime} \quad \lambda\_{\pm \mathbf{S}} = v\_{\mathbf{x}} \pm u\_{\mathbf{S}\prime} \quad \lambda\_{\pm A} = v\_{\mathbf{x}} \pm u\_{A\prime} \quad \lambda\_{\pm} = v\_{\mathbf{x}\prime} \tag{A14}
$$

where the indices *F*, *S*, *A*, and *E* correspond to fast, slow, Alfvén, and entropy characteristics. The values *uF* and *uS* describe fast and slow magnetosonic velocities, and *uA* is the Alfvén velocity. Following the paper of [80], we introduce intermediate values (Roe's averages) for density, velocity, magnetic field and total enthalpy:

$$\rho = \sqrt{\rho\_L \rho\_R} \quad \text{ } \mathbf{v} = \frac{\sqrt{\rho\_L} \mathbf{v}\_L + \sqrt{\rho\_R} \mathbf{v}\_R}{\sqrt{\rho\_L} + \sqrt{\rho\_R}},\tag{A15}$$

$$\mathbf{B} = \frac{\sqrt{\rho\_{\rm R}} \mathbf{B}\_{L} + \sqrt{\rho\_{\rm L}} \mathbf{B}\_{\rm R}}{\sqrt{\rho\_{\rm L}} + \sqrt{\rho\_{\rm R}}}, \quad h = \frac{\sqrt{\rho\_{\rm L}} h\_{\rm L} + \sqrt{\rho\_{\rm R}} h\_{\rm R}}{\sqrt{\rho\_{\rm L}} + \sqrt{\rho\_{\rm R}}}. \tag{A16}$$

The Roe matrix for the equations of one-dimensional multi-component magnetic hydrodynamics (A7), (A8) has the following structure

$$
\hat{A} = \begin{pmatrix} A\_{ik} & & & \begin{bmatrix} 0 & \dots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \dots & 0 \end{bmatrix} \\ \hline -\tilde{\xi}\_1 v\_\chi & \tilde{\xi}\_1 & 0 & \dots & 0 & v\_\chi & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ -\tilde{\xi}\_N v\_\chi & \tilde{\xi}\_N & 0 & \dots & 0 & 0 & \dots & v\_\chi \end{pmatrix} . \tag{A17}
$$

As it can be seen, this matrix consists of four blocks. In the upper left block of dimension 7 × 7 there are elements of the *Aik* Roe matrix for magnetic hydrodynamics. The indices *i* and *k* run through values from 1 to 7. All elements of the upper right block of dimension *N* × 7 are equal to zero. The lower left block has dimension 7 × *N*. In the line with the number *s* of this block, the first element is equal to −*ξsvx*, the second element is equal to *ξs*, and the remaining elements are zero. The lower right block of dimension *N* × *N* is diagonal, and all its diagonal elements are equal to *vx*.

Let us write out the Roe's conditions (A11) for the matrix (A17). For the matrix components corresponding to the equations of magnetic hydrodynamics, they do not give anything new:

$$
\Delta F^i = \sum\_{k=1}^7 A\_{ik} \Delta u^k. \tag{A18}
$$

Consequently, the structure of magnetohydrodynamic part of the matrix *Aik* is not affected and remains the same. For the selected component *ξ<sup>s</sup>* = *ξ* we have the relation:

$$
\Delta(\rho \pounds v\_{\text{x}}) = -\pounds v\_{\text{x}} \Delta \rho + \pounds \Delta(\rho v\_{\text{x}}) + v\_{\text{x}} \Delta(\rho \pounds). \tag{A19}
$$

From here, we find the Roe's average for the value *ξ*,

$$\tilde{\xi} = \frac{\tilde{\xi}\_L \sqrt{\rho\_L} + \tilde{\xi}\_R \sqrt{\rho\_R}}{\sqrt{\rho\_L} + \sqrt{\rho\_R}}.\tag{A20}$$

This expression is valid for any *ξ<sup>s</sup>* component.

Taking into account the block structure of the Roe matrix *<sup>A</sup>*<sup>ˆ</sup> <sup>∗</sup> it can be written

$$\det(\hat{A} - \lambda \hat{I}) = \det(A\_{ik} - \lambda \delta\_{ik}) (v\_{\ge} - \lambda)^N,\tag{A21}$$

where *<sup>δ</sup>ik* is a unit matrix of dimension 7 <sup>×</sup> 7, and <sup>ˆ</sup>*<sup>I</sup>* is a unit matrix of full size. It follows that the eigenvalues of the Roe matrix *A*ˆ are the eigenvalues (A14) of the matrix *Aik*, as well as *N*-multiple eigenvalues

*λ<sup>s</sup>* = *vx*, *s* = 1, . . . , *N*, (A22)

corresponding to the Equation (A5) for the mass fractions of *ξs*.

#### *Appendix A.2. Eigenvectors*

The right eigenvectors of the Roe matrix *A*ˆ are denoted as follows:

$$\mathbf{r} = \left(r^1, \dots, r^{\mathsf{T}}, \vec{r}^1, \dots, \vec{r}^N\right)^T,\tag{A23}$$

where *T* denotes transposition, the first 7 components of the vector correspond to the magnetohydrodynamic part, the remaining *N* components, marked with tildes, correspond to the Equation (A5) for the mass fractions *ξs*. Let us write out the equations for the right eigenvectors

$$
\hat{A} \cdot \mathbf{r} = \lambda \mathbf{r}.\tag{A24}
$$

We have:

$$\sum\_{k=1}^{7} A\_{ik} r^k = \lambda r^i, \quad i = 1, \dots, 7,\tag{A25}$$

$$-\mathfrak{J}\_s v\_{\mathbf{x}} r^1 + \mathfrak{J}\_s r^2 + v\_{\mathbf{x}} \widetilde{r}^s = \lambda \widetilde{r}^s, \quad s = 1, \ldots, N. \tag{A26}$$

Let us first consider the MHD characteristics when *λ* are determined by the eigenvalues of (A14) of the matrix *Aik*. It follows from the first Equation (A25) that in this case the components of *r<sup>i</sup>* coincide with the corresponding components of the right eigenvectors in magnetic hydrodynamics. The remaining components of the right eigenvectors can be found from the additional Equation (A26). If *<sup>λ</sup>* = *vx*, then for each component *<sup>r</sup>*˜*<sup>s</sup>* of the right vector, it can be written:

$$\mathcal{F}^s = \frac{\mathfrak{f}\_s}{\lambda - v\_\chi} \left(r^2 - v\_\chi r^1\right). \tag{A27}$$

Hence, for fast and slow characteristics we find *r*˜*<sup>s</sup>* = *α<sup>f</sup> ξ<sup>s</sup>* and *r*˜*<sup>s</sup>* = *αsξs*, respectively, and for Alfvén *r*˜*<sup>s</sup>* = 0. The coefficients *α<sup>f</sup>* and *α<sup>s</sup>* determine the normalization of Roe matrix eigenvectors in magnetic hydrodynamics [79,80] (see also [68]). For the entropy characteristic *λ* = *vx*, *r*<sup>1</sup> = 1, *r*<sup>2</sup> = *vx* and, consequently, all additional Equation (A26) are satisfied automatically. This means that in this case we can choose *r*˜*<sup>s</sup>* in an arbitrary way. Without limiting generality, they can be put equal to zero. As a result, we obtain the following right eigenvectors

$$\begin{aligned} \mathbf{r}\_{\pm F} &= \left(r\_{\pm F\_{\prime}}^{1}, \dots, r\_{\pm F\_{F}}^{\top} \boldsymbol{\alpha}\_{F} \boldsymbol{\xi}\_{\mathbb{S}\_{1}}, \dots, \boldsymbol{\alpha}\_{F} \boldsymbol{\xi}\_{\mathbb{S}}\right)^{T}, \\ \mathbf{r}\_{\pm S} &= \left(r\_{\pm S'}^{1}, \dots, r\_{\pm S'}^{\top} \boldsymbol{\alpha}\_{S} \boldsymbol{\xi}\_{\mathbb{S}\_{1}}, \dots, \boldsymbol{\alpha}\_{S} \boldsymbol{\xi}\_{\mathbb{S}}\right)^{T}, \\ \mathbf{r}\_{\pm A} &= \left(r\_{\pm A'}^{1}, \dots, r\_{\pm A'}^{\top} \boldsymbol{0}, \dots, \boldsymbol{0}\right)^{T}, \\ \mathbf{r}\_{E} &= \left(r\_{E'}^{1}, \dots, r\_{E'}^{\top} \boldsymbol{0}, \dots, \boldsymbol{0}\right)^{T}. \end{aligned} \tag{A28}$$

For additional characteristics corresponding to the eigenvalues of *λs*, we come to the following equations:

$$\sum\_{k=1}^{7} A\_{ik} r^k = v\_{\mathcal{X}} r^i, \quad i = 1, \dots, 7,\tag{A29}$$

$$r^2 - v\_{\mathfrak{X}}r^1 = 0, \quad \text{s} = 1, \dots, \text{N}.\tag{A30}$$

The equations for the *r<sup>i</sup>* components coincide with the corresponding equations for the entropy characteristic. Therefore, there should be *r<sup>i</sup>* = *r<sup>i</sup> <sup>E</sup>*. The relation (A30) is also satisfied identically. The values of *r*˜*<sup>s</sup>* remain arbitrary. They must be selected in such a way as to obtain a linearly independent set of vectors. To do this, it is enough to set *r*˜*<sup>s</sup>* = 1 for the characteristic corresponding to the index *s*, and to set the remaining components equal to zero. As a result, we obtain the following set of additional right eigenvectors

$$\mathbf{r}\_{\mathbf{s}} = \left(r\_{E'}^1 \dots \,\_{r} r\_{E'}^{\mathcal{T}} \,\_{0} \mathbf{0}, \dots, \mathbf{1}, \dots, \mathbf{0}\right)^{\mathcal{T}}.\tag{A31}$$

Here, the unit is located in place of the component numbered 7 + *s*. Since the vectors *r<sup>s</sup>* and *r<sup>E</sup>* are linearly independent, and the eigenvalues for these characteristics are the same (*λ* = *vx*), then, without violating linear independence, the vector *r<sup>s</sup>* can be replaced by the difference *<sup>r</sup><sup>s</sup>* <sup>−</sup> *<sup>r</sup>E*. In this case, we have a set of unit vectors

$$\mathbf{r}\_{\mathbf{s}} = (0, \dots, 0, 0, \dots, 1, \dots, 0)^T. \tag{A32}$$

Let us now consider the left eigenvectors of the Roe matrix *A*ˆ:

$$1 = (l\_1, \dots, l\_7, \tilde{l}\_1, \dots, \tilde{l}\_N)\_\prime \tag{A33}$$

where the first seven components of the vector again correspond to the magnetohydrodynamic part, and the remaining *N* components, marked with tildes, correspond to the Equation (A5) for the mass fractions *ξs*. Describing relations for left eigenvectors

$$\mathbf{1} \cdot \hat{A} = \lambda \mathbf{1},\tag{A34}$$

we come to the following equations:

$$\sum\_{k=1}^{7} l\_k A\_{k1} - \sum\_{r=1}^{N} \tilde{l}\_r \tilde{\zeta}\_r v\_x = \lambda l\_1 \tag{A35}$$

$$\sum\_{k=1}^{7} l\_k A\_{k2} + \sum\_{r=1}^{N} \tilde{l}\_r \mathfrak{f}\_r = \lambda l\_{2\prime} \tag{A36}$$

$$\sum\_{k=1}^{7} l\_k A\_{kn} = \lambda l\_{n\nu} \quad n = 3, \dots, 7,\tag{A37}$$

$$
v\_{\rm x} \tilde{l}\_{\rm s} = \lambda \tilde{l}\_{\rm s} \quad \text{s} = 1, \ldots, \text{N}.\tag{A38}$$

For fast, slow, and Alfvén characteristics, the eigenvalue is *λ* = *vx*. Therefore, for them all the additional components are ˜ *ls* = 0. For the entropy characteristic, the values of ˜ *ls* are arbitrary. It can also be chosen them ˜ *ls* = 0, so as not to change anything in the magnetohydrodynamic part of the left eigenvector. As a result, for all MHD characteristics, the left eigenvectors will have the following form:

$$\begin{aligned} 1^{\pm F} &= \left( l\_1^{\pm F}, \dots, l\_7^{\pm F}, 0, \dots, 0 \right), \\ 1^{\pm S} &= \left( l\_1^{\pm S}, \dots, l\_7^{\pm S}, 0, \dots, 0 \right), \\ 1^{\pm A} &= \left( l\_1^{\pm A}, \dots, l\_7^{\pm A}, 0, \dots, 0 \right), \\ 1^E &= \left( l\_1^E, \dots, l\_7^E, 0, \dots, 0 \right). \end{aligned} \tag{A39}$$

For additional characteristics *λs*, the last Equation (A38) is satisfied automatically. The remaining equations are reduced to the form:

$$\sum\_{k=1}^{7} l\_k A\_{k1} - \sum\_{r=1}^{N} l\_r \mathfrak{F}\_r v\_x = v\_x l\_{1\prime} \tag{A40}$$

$$\sum\_{k=1}^{7} l\_k A\_{k2} + \sum\_{r=1}^{N} l\_r \mathfrak{f}\_r = v\_{\ge} l\_{2\prime} \tag{A41}$$

$$\sum\_{k=1}^{7} l\_k A\_{kn} = v\_{\mathbf{x}} l\_{\mathbf{n}} \quad \text{ } n = \mathbf{3}, \dots, 7. \tag{A42}$$

We can see that in the case of ˜ *lr* = 0, *li* = *l E <sup>i</sup>* follows from these equations. However, then we find a linearly dependent set of vectors. To obtain a linearly independent set of vectors, for an additional characteristic with the index *s*, we choose ˜ *ls* = 1, and set the remaining additional components equal to zero. For the rest of the components, we will look for a solution in a more general form:

$$l\_i = \chi l\_i^E + \mathfrak{x}\_{i\nu} \tag{A43}$$

where *χ* is the normalizing factor, and *xi* is the unknown quantities. Since

$$\sum\_{k=1}^{7} l\_k^E A\_{ki} = v\_x l\_i^E \quad ; \quad i = 1, \dots, 7,\tag{A44}$$

then, for *xi*, we come to the equations:

$$\sum\_{k=1}^{7} \mathbf{x}\_k A\_{k1} - \underline{\mathfrak{F}}\_s \mathbf{v}\_{\mathbf{x}} = \mathbf{v}\_{\mathbf{x}} \mathbf{x}\_{1\prime} \tag{A45}$$

$$\sum\_{k=1}^{7} \mathbf{x}\_k A\_{k2} + \mathfrak{f}\_s = v\_x \mathbf{x}\_{2\prime} \tag{A46}$$

$$\sum\_{k=1}^{7} \mathbf{x}\_k A\_{kn} = \mathbf{v}\_{\mathbf{x}} \mathbf{x}\_{n\prime} \quad n = \mathbf{3}, \dots, 7. \tag{A47}$$

Without limiting generality, we can assume that only the value *x*<sup>1</sup> is non-zero. Recall [80] that the first row of the matrix *Aik* contains a single non-zero element *A*<sup>12</sup> = 1. Therefore, the Equation (A47) is satisfied automatically, and the Equations (A45) and (A46) coincide and give the solution *x*<sup>1</sup> = −*ξs*. As a result, the left eigenvectors for additional characteristics will have the following form:

$$\mathbf{1}^s = \left( \chi l\_1^E - \mathfrak{f}\_{s\prime} \chi l\_2^E, \dots, \chi l\_7^E, 0, \dots, 1, \dots, 0 \right). \tag{A48}$$

Here, the unit is located again in place of the component with the index 7 + *s*.

The coefficient *χ* should be found from the condition of orthonormality of eigenvectors

$$\mathbf{1}^{\alpha} \cdot \mathbf{r}\_{\beta} = \delta^{\alpha}\_{\beta}. \tag{A49}$$

For MHD characteristics, when *α* and *β* correspond to ±*F*, ±*S*, ±*A* or *E*, all components of ˜ *ls* = 0. Therefore

$$\mathbf{1}^{\alpha} \cdot \mathbf{r}\_{\beta} = \sum\_{k=1}^{7} l\_k^{\alpha} r\_{\beta}^k = \delta\_{\beta}^{\alpha} \tag{A50}$$

due to the orthonormality of the MHD vectors. If *α* corresponds to MHD characteristics, and *β* is the additional characteristics, then it's obvious,

$$
\mathbf{^ar} \cdot \mathbf{r}\_s = \mathbf{0}.\tag{A51}
$$

If we consider vectors for additional characteristics *α* = *r*, *β* = *s*, then we also find the necessary relation,

**l**

$$\mathbf{l}^r \cdot \mathbf{r}\_s = \delta^r\_s. \tag{A52}$$

Finally, let us consider the case, when *α* = *s*, and *β* corresponds to the MHD characteristics. We have:

$$\mathbf{1}^{\mathbf{s}} \cdot \mathbf{r}\_{\beta} = \chi \sum\_{k=1}^{7} l\_k^E r\_{\beta}^k - \mathfrak{zeta}\_s r\_{\beta}^1 + \mathfrak{r}\_{\beta}^s. \tag{A53}$$

If *α* = ±*F*, then the first term in the right part vanishes due to the orthonormality of the MHD vectors. Then, *r*<sup>1</sup> <sup>±</sup>*<sup>F</sup>* <sup>=</sup> *<sup>α</sup><sup>f</sup>* , *<sup>r</sup>*˜*<sup>s</sup>* <sup>±</sup>*<sup>F</sup>* <sup>=</sup> *<sup>α</sup><sup>f</sup> <sup>ξ</sup><sup>s</sup>* and, consequently, the entire right part turns out to be zero. The situation is similar in cases where *α* corresponds to ±*S* and ±*A*. For the entropy characteristic *α* = *E* we have:

$$\mathbf{1}^{\rm s} \cdot \mathbf{r}\_E = \chi \sum\_{k=1}^{7} l\_k^E r\_E^k - \mathfrak{f}\_s r\_E^1 = 0. \tag{A54}$$

Hence, using the normalization condition of entropy MHD vectors, as well as the equality *r*1 *<sup>E</sup>* = 1, we find *χ* = *ξs*. Thus, we finally find

$$\mathbf{1}^{s} = \left(\mathfrak{f}\_{s}^{\mathfrak{r}} l\_{1}^{\mathbb{E}} - \mathfrak{f}\_{s\prime}^{\mathfrak{r}} \mathfrak{f}\_{s}^{\mathfrak{r}} l\_{2}^{\mathbb{E}}, \dots, \mathfrak{f}\_{s}^{\mathfrak{r}} l\_{7}^{\mathbb{E}}, 0, \dots, 1, \dots, 0\right). \tag{A55}$$

Here, as before, the unit is located again in place of the component with the index 7 + *s*.

Using the expressions obtained for the left eigenvectors, it is easy to calculate the characteristic amplitudes of Δ*S<sup>α</sup>* = **l** *<sup>α</sup>* · <sup>Δ</sup>**u**. Since for all MHD characteristics the additional components of the left eigenvectors ˜ *ls* = 0, the corresponding expressions for characteristic amplitudes do not change. For additional characteristics of *λ<sup>s</sup>* we have:

$$
\Delta \mathbf{S}^{\mathbf{s}} = \mathfrak{J}\_{\mathbf{s}} \sum\_{k=1}^{7} l\_{\mathbf{k}}^{E} \Delta u^{k} - \mathfrak{J}\_{\mathbf{s}} \Delta \rho + \Delta(\mathfrak{J}\_{\mathbf{s}} \rho). \tag{A56}
$$

Taking into account (A20), this expression can be simplified,

$$
\Delta S^s = \mathfrak{f}\_s \Delta S^E + \rho \Delta \mathfrak{f}\_{s\nu} \tag{A57}
$$

where *ρ* is the Roe's average for density (A15).

#### *Appendix A.3. Test Calculations*

To numerically solve the equations of multi-component magnetic hydrodynamics (A1)–(A5), we use the finite-difference Roe scheme [77], some details of which are described above. To improve the accuracy, we apply the Osher increasing correction [81]. The resulting difference scheme belongs to the class of TVD (total variation diminishing) schemes [82] and for the case of single-fluid magnetic hydrodynamics is described in detail in [83]. The scheme has the first order of approximation in time and the third order in space. Note that the magnetohydrodynamic version of the Roe scheme in our scheme is presented in such a way that in the absence of a magnetic field (*B* = 0) this scheme exactly passes into the Roe-Einfeldt-Osher scheme we used in purely gas-dynamic simulations [17].

The main disadvantage of the Roe method should, apparently, be considered that the linear system (A10), qualitatively repeating the solution of the original problem (A7), does not reproduce centered rarefaction waves. Instead, the solution consists of a jumps system propagating at speeds corresponding to the eigenvalues of the Roe matrix *A*ˆ(*uL*, *uR*), and some of these jumps may not satisfy the condition of evolutionarity. However, for most cases, the Roe method works well, even if the exact solution of the original problem involves rarefaction waves.

The exception is solutions with trans-sonic rarefaction waves, for the correct accounting of which special evolutionary corrections are used. In the case of magnetic hydrodynamics equations (single-fluid or multi-fluid), such corrections should be used for fast and slow magnetosonic characteristics, in which rarefaction waves with corresponding properties may sometimes occur. For fast magnetosonic characteristics, our difference scheme uses the Einfeldt [84] correction. In the case of slow magnetosonic characteristics, the initial Einfeldt correction does not agree with a purely gas-dynamic version of the difference Roe scheme. Therefore, for such characteristics, we use a modified entropy correction [83].

In multi-dimensional problems, when using Godunov's methods for solving equations of magnetic hydrodynamics, a separate procedure is required to clear the divergence of the magnetic field, since the difference scheme operates with values averaged over the volume of a cell and, therefore, does not conserve the magnetic flux. The method we use to clear the magnetic field divergence is described in Section 4.4. In the test calculations, the results of which are given in this section, divergence cleaning was not applied.

In the first test calculation, numerical simulation of the decay of an arbitrary MHD discontinuity was carried out. At the initial moment of time, the resting matter (velocity *<sup>v</sup>* <sup>=</sup> 0) was considered in a homogeneous magnetic field *Bx* <sup>=</sup> 1/4, *By* <sup>=</sup> <sup>√</sup>3/4, *Bz* <sup>=</sup> 0. An arbitrary discontinuity was located at the point *x* = 0.

In the region to the left of the discontinuity *x* < 0, the following values were set: density *ρ* = 1, pressure *P* = 1. In the region to the right of the discontinuity *x* > 0, the parameters were set: density *ρ* = 0.125, pressure *P* = 0.1. The adiabatic index *γ* = 5/3. A calculated grid containing 512 cells was used. A two-component mixture consisting of hydrogen ions H<sup>+</sup> and helium He<sup>+</sup> was considered. We assume that at the initial moment of time, the matter to the left of the discontinuity consists only of helium ions, and the matter to the right of the discontinuity consists only of hydrogen ions.

The results of the numerical solution of this problem obtained at the time *t* = 0.15 are shown in Figure A1. The decay of the initial arbitrary discontinuity leads to the formation of two rarefaction waves (fast 1 and slow 2) propagating into the region of helium matter. The contact discontinuity 3 propagates to the right. Two shock waves are formed in the region of hydrogen plasma (slow 4 and fast 5). The panel on the bottom right shows the resulting distribution of the mass fractions of helium ions *ξ*<sup>1</sup> and hydrogen ions *ξ*2. The contact boundary between the matters shifts to the right, is clearly localized in space and does not contain any non-physical oscillations. Analysis of the figure shows that the scheme approximates all types of running waves quite well.

As a second demonstration example, let us consider the results of a test calculation of the problem on a volume-distributed explosion in a medium with a homogeneous magnetic field. The initial conditions correspond to the situation when, in the entire volume of an infinitely long cylinder with a radius of *R* = 0.2 with a density of *ρ* = 1, the pressure instantly increases by 10 times compared to the pressure in the external environment. As a result of such an explosion, the matter of the cloud begins to expand into the external environment.

The density and pressure in the external environment were set as follows: *ρ*ext = 0.125, *P*ext = 0.1. The magnetic field at the initial time in the entire volume of the computational domain was homogeneous: *Bx* <sup>=</sup> 1/4, *Bx* <sup>=</sup> <sup>√</sup>3/4, *Bx* <sup>=</sup> 0. The cloud was filled with helium ions He+, and the external environment was filled with hydrogen ions H+. Due to the formulation of the problem in each plane *z* = const, the flow pattern will be the same.

**Figure A1.** The results of the test calculation of the Riemann problem on the decay of an arbitrary discontinuity (see the description of the initial conditions in the text). Density distributions are shown (**top left**), velocity components *vx* and *vy* (**top right**), magnetic field components *by* (**bottom left**), as well as mass fractions *ξ*<sup>1</sup> and *ξ*<sup>2</sup> (**bottom right**) at time *t* = 0.15. The figures correspond to the digits: 1—fast rarefaction wave, 2—slow rarefaction wave, 3—contact gap, 4—slow shock wave, and 5—fast shock wave.

Therefore, the problem is actually two-dimensional. Calculations were carried out in the computational domain −0.5 ≤ *x* ≤ 0.5, −0.5 ≤ *y* ≤ 0.5 on a mesh with the number of cells 512 × 512. It should be noted that according to its formulation, the previous test problem corresponds to the decay of an arbitrary discontinuity on the right edge of the cloud along the *x* coordinate.

The results of the calculation of the second test problem at time moment *t* = 0.1 are shown in Figure A2. The left panel shows the distribution of density (color) and velocity (arrows). The boundary of the cloud that separates helium and hydrogen corresponds to a solid line. The right panel shows the distribution of density (iso-lines) and magnetic field (lines with arrows). As it can be seen from the figure, the flow structure is determined by a complex system of strong MHD discontinuities propagating outward (fast and slow shock waves, contact discontinuity), as well as fast and slow rarefaction MHD waves propagating to the center.

The presence of the magnetic field leads to the fact that the flow pattern is anisotropic. In particular, the contact boundary along the magnetic field propagates faster than in the direction across the field. As a result, the initially symmetrical helium cloud becomes elongated along the magnetic force lines. The analysis of the figure allows us to conclude that the computational qualities of the difference scheme we used, noted above for a one-dimensional problem, are preserved in the multi-dimensional case.

**Figure A2.** The results of the test calculation of the problem on a volume-distributed explosion in a homogeneous magnetic field. The density (color) and velocity distributions (arrows) are shown on the left. The solid line corresponds to the border of the cloud. The density distributions (iso-lines) and magnetic field distributions (lines with arrows) are shown on the right. The distributions of the quantities are presented at time moment *t* = 0.1.

#### **References**

