Magnetized winds and their influence in the escaping upper atmosphere of HD 209458b

ABSTRACT

1 INTRODUCTION

The past two decades have been revolutionary when it comes to the study of planetary science. Since the confirmation of the first pulsar planet (Wolszczan 1994), and the first discovery through the transiting technique (Charbonneau et al. 2000), more than two thousand planets beyond our Solar system have been detected. Until then, planetary atmospheres were thought to be subject to a relatively calm environment. Yet, more than two-thirds of the first discovered exoplanets were observed so close to their host star (a ⩽ 0.5 au), as to imply a much more active interaction with the stellar wind and radiation.

For some of these exoplanetary systems, observations in the Lyman α line provided evidence of neutral atmospheric material loss, as a consequence of the heating produced by the enormous amount of UV flux received from their host star, e.g. HD 209458b (Vidal-Madjar et al. 2003), HD 189733b (Lecavelier Des Etangs et al. 2010; Lecavelier des Etangs et al. 2012; Jensen et al. 2012; Ben-Jaffel & Ballester 2013), Wasp 12b (Fossati et al. 2010; Haswell et al. 2012; Nichols et al. 2015), GJ 436b (Kulow et al. 2014; Ehrenreich et al. 2015) and possibly 55 Cnc b (Ehrenreich et al. 2012).

Of particular interest is the work of Vidal-Madjar et al. (2003) where the first Ly α transit observations were analysed, and the presence of an escaping neutral atmosphere was proposed. The authors found a |$10{{\rm \,per\,cent}}$| absorption at −100 km s⁻¹, an asymmetric line profile with more absorption in the blue than in the red part of the line, and a total absorption of |$[15 \pm 4]\,{{\rm \,per\,cent}}$| in the range of ±300 km s⁻¹.

Later, transit absorption observations of atomic lines from H i, O i, C ii, Si iii and Mg i (Vidal-Madjar et al. 2003, 2004; Ballester, Sing & Herbert 2007; Ehrenreich et al. 2008; Linsky et al. 2010; Jensen et al. 2012; Vidal-Madjar et al. 2013) helped to confirm the presence of such a hydrodynamic planetary wind.

We know now that the features present in the Ly α line during the transit of HD 209458b can be explained by a combination of several physical processes such as the stellar–planetary wind interaction, the stellar radiation pressure, and the charge exchange of planetary neutral atoms with stellar ions, etc.

Numerical and theoretical studies have been carried out with the inclusion of one or several of these processes. For instance, the expansion of the planetary atmosphere and the resulting planetary wind, as a consequence of the incident UV stellar flux, have been studied in the works of Murray-Clay, Chiang & Murray (2009), Koskinen et al. (2010, 2013), Guo (2011), Salz et al. (2016), among others. These works give a good description of the lower layers of planetary atmospheres; however, they do not include the dynamical interaction with the stellar wind.

The stellar–planetary wind interaction is considered using hydrodynamic simulations in Schneiter et al. (2007), Villarreal D’Angelo et al. (2014), Schneiter et al. (2016) and through particle simulation including the charge exchange process in Erkaev et al. (2007), Holmström et al. (2008), Tremblin & Chiang (2013), Bourrier & Lecavelier des Etangs (2013), Kislyakova et al. (2014) and Christie, Arras & Li (2016). These works were able to partially reproduce the observed absorption in the Ly α line even though they did not take into account the stellar or the planetary magnetic fields.

Planetary magnetic fields can shield the atmosphere of the planet from the direct interaction with the stellar wind, as they deflect the stellar winds and prevent their penetration into the lower layers of planetary atmospheres. Several studies indicate that, if present, an exoplanetary magnetic field will ultimately determine the amount of atmospheric material that is lost from the planet (Adams 2011; Trammell, Arras & Li 2011; Owen & Adams 2014; Trammell, Li & Arras 2014; Khodachenko et al. 2015). However, the presence of a planetary magnetic field is still an open question.

From theoretical calculations, the expected magnetic moment for hot-Jupiter like planets should be only a few times lower than the magnetic moment of Jupiter (⁠|$\mathcal {M}_\mathrm{J}=1.56\times 10^{27}$| A m²) (Sánchez-Lavega 2004; Durand-Manterola 2009). Nevertheless, Christensen, Holzwarth & Reiners (2009) predict magnetic field strengths, depending on the internal heat flux, that could be an order of magnitude higher than that of Jupiter. Recently, Rogers & McElwaine (2017) proposed that a small planetary dynamo, due to conductivity variations arising from the strong asymmetric stellar heating, could exist in these type of planets.

The magnetic field value of HD 209458b remains unknown. Moreover, the magnetic field strength for its host star has not yet been detected via spectropolarimetric observations (Mengel et al. 2017). Analytical calculations for the planetary magnetic moment have given values in the range [1–|$8]\,\mathcal {M}_\mathrm{J}$| (Sánchez-Lavega 2004; Durand-Manterola 2009). Kislyakova et al. (2014) proposed a magnetic moment of |$0.1\,\mathcal {M}_\mathrm{J}$| based on the transit observations in the Ly α line. Trammell et al. (2014) also employed the Ly α observations to constrain the magnetic field of HD 209458b. In that work, the authors found that atmospheric material can be supported within closed magnetic field lines in what is known as a wind ‘dead zone’ [also found by Trammell et al. (2011) and Khodachenko et al. (2015)]. When compared to the Ly α transit absorption, this trapped material can reproduce the observations of Ben-Jaffel (2008) for a planetary magnetic field of 10 G (⁠|$2.8\,\mathcal {M}_{\rm J}$|⁠) or less. Others' works that aimed to explore the effects that the magnetic field has on the atmospheric expansion have adopted similar values (see Grießmeier et al. 2004; Khodachenko et al. 2012, 2015; Owen & Adams 2014). Erkaev et al. (2017) used a 3D MHD model to reproduce the interaction of a non-magnetized planet and the stellar wind. In order to reproduce the Ly α observations as in Kislyakova et al. (2014), they concluded that the planet should have an intrinsic magnetic moment of about [0.13–|$0.22]\,\mathcal {M}_\mathrm{J}$|⁠.

In previous works (Schneiter et al. 2007; Villarreal D’Angelo et al. 2014; Schneiter et al. 2016) we have studied how the hydrodynamic wind–wind interaction shapes the Ly α profile during the planetary transit. Using 3D hydrodynamics simulations, we correlated the planetary mass-loss rate (⁠|$\dot{M}_\mathrm{p}$|⁠) with the observed total absorption in this line. Schneiter et al. (2007) assumed solar wind conditions for a fully ionized stellar wind and a neutral planetary wind, proposing an upper limit for |$\dot{M}_\mathrm{p}$|⁠. Villarreal D’Angelo et al. (2014) explored the influence of different stellar and planetary wind conditions and structures for the planetary wind and, finally, Schneiter et al. (2016) included explicitly the photoionization process in the calculations, thus making the proposed model closer to self-consistent. When contrasted with the Ly α transit observations, all of them successfully reproduced the main features mentioned in Vidal-Madjar et al. (2003), but they fail to reproduce the absorption in the red part of the line profile, which we have successfully reproduced in the present work by including magnetic fields.

Until now, there are only a few 3D MHD simulations that consistently model the stellar–planetary wind interaction (see Cohen et al. 2011; Matsakos, Uribe & Königl 2015; Tilley, Harnett & Winglee 2016) but none of them is based on HD 209458b. This work is an effort to explore how the magnetic field affects the Ly α line during transit, considering both the star and the planetary magnetic fields with different magnetic moments. Due to the collisional coupling between ions and neutrals in the upper planetary atmosphere, we expect a measurable influence of the magnetic field on the observed Ly α profile, as suggested by Trammell et al. (2014) and Erkaev et al. (2017). Our 3D MHD models have the star and the planet in the same physical domain. They include a proper treatment of the stellar photoionization and the effects of cooling and heating of the neutral material.

Section 2 introduces the new MHD version of the guacho code. Appendix  A presents some of the 1D/2D/3D MHD tests that were performed in order to validate the MHD implementation. The numerical setup is presented in Section 3. Results, analysis and comparison with observations are treated in Section 4. Finally, the conclusions are given in Section 5.

2 THE MHD CODE: GUACHO

The set of equations (1)–(4) is advanced in time with a second-order Godunov method. The fluxes are calculated with the HLLD Riemmann solver of Miyoshi & Kusano (2005), and a linear reconstruction of the primitive variables is applied using the minmod slope limiter to ensure stability.

As in every MHD code, the magnetic field divergence must be maintained close to zero. For this purpose, the flux-interpolated Central Difference scheme developed by Tóth (2000) was implemented.

The right-hand sides of equations (2) and (3) show the source terms included in our simulations. We will explain them separately in the following subsections.

2.1 Gravity and the stellar radiation pressure

In the current models, we ignore the dependence of μ with the relative velocities of the H atoms with respect to the star. A more careful treatment of the stellar radiation pressure force that includes such dependency is currently in consideration (Schneiter et al. in preparation), but the preliminary results show that it is not very relevant for a stellar wind with the characteristics of HD 209458.

2.2 Heating, cooling and radiative transfer

The radiative field of the star is included with the ray tracing method used in Schneiter et al. (2016). The total stellar EUV flux is divided in 10⁶ photon packages that are launched from the stellar surface at random positions in random directions. As these packages travel through the grid, they are attenuated by a factor e^−Δτ, with |$\Delta \tau =a_0\,n_\mathrm{{H_I}}\,\Delta l$|⁠. This factor depends on the neutral material within each cell, the path-length (Δ_l) and the photo-ionization cross-section that we have assumed to be a₀ = 6.3 × 10⁻¹⁸ cm², considering that all the photons are at the Lyman limit. The photo-ionization rate ϕ is then calculated by equating the ionizing photon-rate S_⋆ with the ionization per unit time within each cell (⁠|$S_\star =n_{\mathrm{H_I}}\,\phi \,\mathrm{dV}$|⁠). The contribution of the photo-ionization rate is then added to equation (6) and also used to calculate the heating term in equation (3). The heating per unit time and volume is |$\psi =\phi \,E_0$|⁠, with an energy gain per ionization of E₀ = 0.86 eV.

The cooling is computed using the prescription described in Biro, Raga & Canto (1995), which includes contributions from recombination and collisional ionization of hydrogen, collisional excitation of the H Ly α line and [Oi], [O ii] forbidden lines. The latter are obtained assuming that the excitation is in the low-density regime (for which we solve the five level atom problems with the atomic parameters from Pradhan & Nahar 2011), and assuming that the ionization state of oxygen follows that of hydrogen (which can be justified by the efficient charge exchange between them). The oxygen cooling rate is multiplied by a factor of ∼7 in order to account for other important coolants (such as C, N and S; see Biro et al. 1995). At temperatures above 5 × 10⁴K (where oxygen is expected to be more than singly ionized) we switch to a coronal equilibrium cooling curve (in fact the parametrization of the cooling calculated by Raymond, Cox & Smith 1976 is given in Raga & Cantó 1989).

3 THE NUMERICAL SETUP

Our models are based on the planetary system HD 209458b, composed by a solar like star and a hot Jupiter like planet. The whole system is simulated in a Cartesian mesh of 460 × 230 × 460 cells with a resolution of 3 × 10⁻⁴ au.

All the models have the same geometrical configuration, with a non-rotating star located at the centre of the coordinate system, and a non-rotating planet in a circular orbit with a radius 0.047 au in the xz-plane.

To initialize our models, we fill the entire grid (except inside the planet wind boundary) with the stellar wind using the isothermal Parker wind solution (Parker 1958) corresponding to a coronal temperature of T_⋆ = 1.5 × 10⁶ K and the dipolar magnetic field of the star. In this medium we overwrite a region around the planet position with the planetary wind structure (see below). For subsequent times, the stellar and planetary wind configurations are re-imposed at every time-step within spheres of radii |$1\,R_\star$| for the stellar wind and |$3\,R_\mathrm{p}$| for the planetary wind. The position of the planetary boundary changes according to the orbit. We also include the orbital speed to the velocity field in the planet wind. All external boundaries of the simulation are treated with outflow boundary conditions.

The simulations evolve until they reach a steady state, which is usually after the planet completes ∼3/4 of its orbit.

The physical parameters of the star/planet system (Torres, Winn & Holman 2008), together with the adopted values for initializing the winds, are presented in Table 1.

3.1 The stellar wind

The initial conditions set inside the star will give rise to a continuously blowing wind. Using the stellar parameters for the solar like star HD 209458, we calculate the values of density, velocity and temperature at the surface of the star (R_⋆) and extrapolate them as constant values inside the stellar radius to fill the entire stellar volume. Since the values inside this boundary are reimposed at every time-step the flow within this region does not evolve withtime.

Fixing the stellar mass-loss rate to the solar value of 2.0 × 10⁻¹⁴ M_⊙ yr⁻¹, we can derive the value of the density at the base of the wind (⁠|$\dot{M}_{\star }=4 \pi \rho _{\star }R_{\star }^2v_\star$|⁠).

The radiative field of the star is simulated with the emission of photons at random positions from the stellar surface. The total amount of photons emitted in random directions is calculated using the stellar luminosity in the EUV, log₁₀(L_EUV) < 27.74 erg s⁻¹, derived in the work of Sanz-Forcada et al. (2011), which is then divided in 10⁶ photon-packages. The adopted stellar photon rate is therefore |${S_{0}=2.5\times 10^{38}\,\mathrm{s^{-1}}}$|⁠, which, assuming that each photon is at the Lyman limit, corresponds to a flux of |${F_0=884\,\mathrm{erg\,cm^{-2}\,s^{-1}}}$| at the orbital distance of HD 209458b.

3.2 The planetary wind

In simulating the wind of HD 209458b we are making the assumption that the upper part of the planetary atmosphere is under an hydrodynamic escape as discussed in Murray-Clay et al. (2009), Guo (2011), Salz et al. (2016). The approach is the same as the one implemented on the star, specifying the initial values of the hydro-dynamical variables for the planetary wind inside a radius, which in this case is |$3\,R_\mathrm{p}$|⁠. At smaller radii the values are kept constant. Launching the planetary wind at this distance allows us to have a lower resolution than the one that would be required to model the generation of the wind from its base (⁠|$1\,R_\mathrm{p}$|⁠). The adopted wind parameters at |$3\,R_\mathrm{p}$| for an |$80{{\rm \,per\,cent}}$| ionized hydrogen atmosphere are presented in Table 1 and are in agreement with the standard planetary wind model derived in Murray-Clay et al. (2009) for the case of an energy limited driven wind. A velocity of 10 km s⁻¹ and a temperature of 10⁴ K are then imposed every time step at this boundary.

The planetary density is obtained from its mass-loss rate, taken to be |$\dot{M}_\mathrm{p}=2\times 10^{10}$| g s⁻¹. This value is consistent with our previous works (Villarreal D’Angelo et al. 2014; Schneiter et al. 2016) and it is also in agreement with recent estimations derived from observations (Salz et al. 2015).

Our initial conditions result in a subsonic planetary wind (a ∼ 12.3 km s⁻¹) launched from a region inside the corresponding sonic radius. The planetary wind will then adjust itself to the conditions outside the launching radius. In all our models, the planetary wind is never entirely suppressed by the stellarwind pressure. However, in the models with large planetary magnetic fields, we see evidence of wind suppression at the equator near the wind launching region. This could correspond to regions where the magnetic pressure dominates (Trammell et al. 2011), but our resolution is insufficient to make any conclusions.

4 RESULTS AND DISCUSSION

To investigate the role of magnetic fields on the neutral material carried with planetary wind, we run a total of five models. Each model has a different value for the polar magnetic field at the surface of the star (B_⋆) and the planet (B_p); the remaining parameters for the initial conditions for both winds (stellar and planetary) are shown in Table 2.

A common feature of the interaction between the planetary and stellar wind is the presence of a shock-like region where both winds meet. The location where the planetary wind effectively stops the stellar wind, the stand-off distance, is the key in controlling the extent of the planetary magnetosphere, and it will depend on the pressure balance at this point. The stand-off distance (R₀) for all the models, measured along the line that joins the star and the planet, is presented in Table 2. Trailing the planet, a cometary tail composed of material from the stellar and planetary winds is formed, a known feature of the wind–wind interaction found in previous works (Schneiter et al. 2007; Villarreal D’Angelo et al. 2014; Schneiter et al. 2016). A 3D rendering of model B1.1 is shown in Fig. 1. Magnetic field lines coloured according to the value of the total magnetic field are shown. The planet is surrounded by contours of the number density that shows the extent of the cometary tail.

Figure 1.

3D view of the magnetic field lines coloured by the intensity of the magnetic field for model B1.1. The contours around the planet position show the number density.

Open in new tabDownload slide

To analyse the characteristics of the different models we present in Fig. 2 the density, temperature and magnetic field values for four of the five models in the xy-plane (i.e. perpendicular to the orbitalplane) centred at the planet position. In this figure, the star is at the right edge of each panel magenta (also in the xy-plane) but out of the field of view. Models B5.1 and B5.5 result in very similar stratification and for this reason we only show model B5.5. This is also true for models B1.0 and B1.1 but, in this case, model B1.0 shows the interaction of the stellar wind with a non-magnetized body.

Figure 2.

Contour plots of density (top row), temperature (middle row) and magnetic field value (bottom row), for a cut in the xy-plane (i.e. perpendicular to the orbital plane) centred at the position of planet. In the density plots, the ionization fraction for values of 0.85, 0.95 and 0.999 is denoted by black contours around the planet.

Open in new tabDownload slide

The first three columns of Fig. 2 correspond to models with B_⋆ = 1 G characterized by a slower and less dense stellar wind in the vicinity of the planet. In the equatorial plane (y = 0) a pronounced current sheet develops as a consequence of the stretching of the stellar magnetic field lines, dragged by the stellar wind. The magnetic field strength (bottom row) there becomes almost negligible indicating that the phenomenon of re-connection is taking place. Possible re-connection sites are also visible at the nose of the planetary magnetosphere, at both sides of the equatorial current sheet, for models B1.1 and B1.5. At these points, the planetary magnetic field lines will connect with the magnetic field of the star. In these cases, no Alfvén wings (Neubauer 1980) would develop because the stellar wind at the planet position is already super-alfvénic (as well as super-sonic).

The last column in Fig. 2 shows one of the models with a higher stellar magnetic field value (B_⋆ = 5 G). Increasing B_⋆ leads to an increase in the temperature and velocity of the stellar wind (Vidotto et al. 2009). With a higher magnetic tension, the magnetic loops are stretched but not open as in the case of B_⋆ = 1 G, remaining closed at the planet position. The planet, orbiting within a denser and hotter environment, does not suffer from reconnection events between its magnetic field and the stellar one in the direction of the star.

As the stellar wind flows around the planet a cavity forms and is filled with the planet’s atmospheric material and magnetic field (when it is present). The properties of the material inside this cavity are governed by the planetary magnetic field, the shock with the stellar wind and the radiative field of the star. This last one will only have a significant influence if the neutral density inside this cavity is optically thick to the stellar photons.

In the case of a non-magnetized planet (model B1.0), the planetary wind fills a region around the planet of lower temperature and higher density than its surrounding. This region’s volume and its temperature and density stratification are almost identical to what is observed in model B1.1, where B_p = 1 G, suggesting that a planetary magnetic field value less than 1 G has the same influence in stopping the stellar wind than in the case of a non-magnetized planet. In both cases, a negligible magnetic pressure contribution is present in the planetary side, thus the sum of thermal and ram pressures are responsible for deflecting the stellar wind. When the planetary magnetic field is increased, as in the model B1.5, the resulting cavity has a larger size, in comparison with B1.0 and B1.1 models. The planetary wind gets more easily accelerated along open magnetic field lines (mainly at the poles), and at the same time, the higher magnetic tension at the equator creates a wind ‘dead-zone’ that gets filled with lower temperature material. Both effects result in a larger cavity.

The stand-off distance at the sub-stellar point, presented inTable 2, confirms that models B1.0 and B1.1 share the properties of the resulting cavity since they have the same value of |$R_0\sim 7\,R_{\rm p}$|⁠. Similar values have been found by Weber et al. (2017) and Khodachenko et al. (2012) when the contribution of a magneto-disc in the shape of the magnetosphere is taken into account. Smaller values have also been found for HD 209458b in the works of Grießmeier et al. (2004) and Kislyakova et al. (2014), although the value of the stellar wind’s velocity and density employed in these works were much higher than in our models. A higher stand-off distance is reached for model B1.5, where the magnetic field of the planet is now playing a more important role in the pressure balance. Models B5.1 and B5.5 present a more compressed planetary magnetosphere, with stand-off distances near |$4\,R_\mathrm{p}$|⁠. The shorter standoff distance is a result of the increased ram pressure on the stellar side (due to the strong stellar wind). For these cases, the value of the planetary magnetic field has little effect in the properties inside the magnetospheric cavity (the stratification of models B5.1 and 5.5 is remarkably similar.). This is due to the fact that the shock with the stellar wind determines the temperature and the velocity of the material of this region.

The ionization fraction of the planetary wind (χ) is plotted in the first row of Fig. 2. The three black contour levels around the planet represent the values of χ = [0.85, 0.95, 0.999]. For B_⋆ = 5 G, a very compressed planetary magnetosphere occurs, and these contours are found very close to the planet, showing that almost all planetary material becomes ionized down to very near the base of the wind. When the planetary magnetosphere is more expanded, the planetary wind can remain partially neutral (up to |$5\,{{\rm \,per\,cent}}$|⁠) further from the wind base (⁠|$3\,R_\mathrm{p}$|⁠), and more so in the tail direction, as can be seen in models with B_⋆ = 1 G. As we mentioned above, a higher value of the planetary magnetic field increases the velocity and temperature in the wind.Hence,ions in the planetary atmosphere are accelerated more and collisional ionization becomes more effective. This can be seen in model B1.5, where the material inside the planetary magnetosphere has a temperature higher than 10⁵ K, and neutral material is more concentrated around the planet than in models B1.0 and B1.1. In this case, the shape of the partially neutral region, even when the magnetosphere is larger, is small compared to models with B_p = 1 G.

4.1 Lyman α profile

_Hwhere φ(Δ|$v$|_los) is a Gaussian line profile,¹a₀ = 0.01105 cm² is the Ly α cross-section at the line centre (Osterbrock 1989), n_i is the neutral density and ds is the length measured along the los. In our calculations, the velocities are projected along the los-direction (x-direction) and tilted |$3{^{\circ}_{.}}29$| around the z-axis in order to resemble the orbital inclination (i = |$86{^{\circ}_{.}}71$|⁠) of the planet as seen from the Earth. The integration of |$\tau _{v_{\rm los}}$| is made from the surface of the star to the end of the computational domain. We let the planet complete almost 5/9 of its orbit to compute the transit absorption. The Lyman α absorption fraction 1 − e^−τ at mid-transit, showed in Fig. 3, is then calculated by integrating the optical depth within the velocity range [−300,300] km s⁻¹ over 250 channels.

Figure 3.

Lyman α absorption fraction (1 − e^−τ), integrated between [−300, 300] km s⁻¹, for the four models presented in Fig. 2. The plot is a zoom at the position of the star (white circle) at the moment of mid-transit. In order to calculate the absorption in the Ly α line we have taken into account the orbital inclination of the planet, represented with the black filled circle.

Open in new tabDownload slide

For every velocity channel we have also calculated the normalized stellar transmission during transit integrated over the stellar disc (using the stellar photospheric radius of |$1.2\,R_\star$|⁠). We do not include the stellar limb-darkening and we assume a constant stellar emission over the adopted velocity range. The results are shown in Fig. 4. Solid lines are used to represent models with B_⋆ = 1 G and dashed lines are used for models with B_⋆ = 5 G, while colours represent different values for the planetary magnetic field: green, red and blue show models with B_p = 0, 1 and 5 G, respectively. We compare our results with the observations from Vidal-Madjar et al. (2003) in the same figure. The yellow stripe in both plots represents the part of the Lyman α line that is contaminated with the geo-coronal emission, which is removed from the analysis in order to compare with the observations.

Figure 4.

(a) Normalized stellar transmission as a function of the los velocity in the Ly α line averaged over the stellar disc as seen by an observer. (b) Comparison with the observational line profile during and off transit, taken from the work of Vidal-Madjar et al. (2003). The yellow stripe in both graphics corresponds to part of the line contaminated with the geo-coronal glow. In both panels, solid lines are used to represent models with B_⋆ = 1 G and dashed lines are used for models with B_⋆ = 5 G, while colours represent different values for the planetary magnetic field: green, red and blue show models with B_p = 0, 1 and 5 G, respectively.

Open in new tabDownload slide

Our calculated line profiles show that models with a higher value of the stellar magnetic field (B_⋆ = 5 G) present a wider and shallower line profile than models with B_⋆ = 1 G. When B_⋆ is 1 G, the absorption profiles expand a narrow velocity range, from [−150, 120] km s⁻¹, with a deeper absorption value at the smallest velocities. Models B1.0 and B1.1 have the same absorption profile since they present almost the same physical response to the interaction with the stellar wind. It is for these models that we found |$10\,{{\rm \,per\,cent}}$| absorption at −100 km s⁻¹ and |$5\,{{\rm \,per\,cent}}$| absorption at 70 km s⁻¹, in agreement with the observations of Vidal-Madjar et al. (2003, 2008). The distribution of these neutrals, with an average temperature of 8 × 10⁴ K, is asymmetric with a more extended region oriented towards the observer, as we can see from the line profile. This asymmetric distribution is also visible inFig. 3.

When B_p = 5 G (model B1.5), the calculated line profile shows absorption at larger velocities towards the star (reaching |$5\,{{\rm \,per\,cent}}$| at 100 km s⁻¹). The shape of the profile is more symmetric when comparing with that obtained for models B1.0 and B1.1 (see also Fig. 3), as a result of a more important contribution from the planetary magnetic field. The material that remains neutral inside the magnetospheric cavity reaches temperatures of roughly2 × 10⁵ K.

For the models with B_⋆ = 5 G, the planetary wind is embedded inside a hotter environment, and with higher stellar ram pressure. As a result of the shock with the stellar wind, the planetary material gets heated to temperatures of up to 1 × 10⁶ K. Nearly |$0.1{{\rm \,per\,cent}}$| of the planetary wind (∼2 × 10⁻²¹ g cm⁻³) remains neutral under these conditions and acquires velocities around 150 km s⁻¹ in both wings of the profile. Higher velocities are reached towards the blue, driven by the pressure of the stellar wind. For these models, the different values of the planetary magnetic field that we used make only a marginal difference in the absorption profile as shown in Fig. 4.

The absorption profiles presented in Fig. 4 can be compared with those presented in fig. 5 of Schneiter et al. (2016) for models B2b and C2b. These models have the same initial conditions for launching the planetary wind and comparable stellar wind properties at the position of the planet, but they do not include a magnetic field. Model B2b has stellar wind conditions comparable with models that have B_⋆ = 1 G in our work, and model C2b has conditions similar to models with B_⋆ = 5 G, except for a lower density of the stellar wind near the planet. The main difference between the hydrodynamic models of Schneiter et al. (2016) and our MHD models is that the new computed Ly α profiles have a higher absorption in the red part of the line. For the models with B_⋆ = 1 G the profiles of these MHD models have a lower radial velocity (towards the observer), which is due to the fact that the field lines are not entirely opened away from the star and the material does not flow as freely in this direction. The absorption profiles of the models with B_⋆ = 5 G are broader compared with model C2b of Schneiter et al. (2016), but this can be attributed to the difference in the stellar wind ram pressure, given the differences in density of the stellar winds.

The total Ly α absorption during transit is computed by integrating the line profile over the velocity range ±300 km s⁻¹, excluding the part of the line contaminated with the geo-coronal emission. The results for all the models are presented in the last column of Table 2. These values can be compared with the one obtained in the work of Vidal-Madjar et al. (2003) where |$[15 \pm 4]\,{{\rm \,per\,cent}}$| of absorption during transit was found. Our results show that the total absorption in the Ly α line is not strictly correlated with the stand-off distance (the magnetosphere size). For instance, model B1.5 has a larger stand-off distance but the smallest value for the total absorption in Ly α during transit. A more accelerated planetary wind leads to an increase of collisions between ions and neutrals, producing a larger degree of ionization close to the planet, with more ionization at the poles (see Fig. 3). Models B1.0 and B1.1 have a total absorption during transit of |$12\,{{\rm \,per\,cent}}$| that agrees, within the errors, with the value from Vidal-Madjar et al. (2003). The highest values (⁠|${\sim } 13\,{{\rm \,per\,cent}}$|⁠) for the integrated line profile are found for models B5.1 and B5.5. When compared with the observations, these models produced a more pronounced absorption during transit as shown in Fig. 4(b). In these cases, the stand-off distance is very close to the planet, indicating a more compressed planetary magnetosphere than in the previous models, due to the stronger stellar winds. This wind also produces a stronger shock, raising the temperature of the material behind up to ∼10⁶ K, producing a broad absorption profile from the few surviving neutrals in the wake.

5 CONCLUSIONS

Our aim was to study the effects that the presence of dipolar magnetic fields, in the star and the planet, have on the escaping neutral material from the atmosphere of HD 209458b. In doing this, we have assumed that the star could have a dipole magnetic field of magnitude 1 and 5 G (at the poles). For the planet we choose to model a non-magnetized planet, and two other examples with a dipole magnetic with magnitudes of 1 and 5 G (at the poles). These values correspond to a planetary magnetic moment of 4.8 × 10²⁶ and 2.4 × 10²⁷ A m², respectively.

We have found that the interaction between both winds leads to the development of a magnetospheric cavity around the planet inside which a partially ionized planetary wind is able to expand. For values of B_⋆ = 1 G, this cavity still forms even when the planet magnetic field is set to zero. Moreover, virtually the same region is developed when B_p is increased from zero to 1 G, indicating that for lower planetary magnetic field values, the magnetic pressure of the planet is not a major contribution in stopping the stellar wind. These two models also share the same stand-off distance (⁠|${\sim } 7\,R_\mathrm{p}$|⁠) and the same total absorption of |$12\,{{\rm \,per\,cent}}$| when we integrated the Lyman α profile in the velocity range of ±300 km s⁻¹ during a transit. When comparing with the observational data, these models are the ones that best reproduce the absorption profile during transit from the work of Vidal-Madjar et al. (2003). From these observations, the authors also found a total absorption during transit of |$[15\pm 4]\,{{\rm \,per\,cent}}$| when integrating the Ly α line within the same velocity range. For a larger value of the planetary magnetic field, model B1.5 shows that the degree of ionization of the planetary neutral material is higher as a consequence of the increase in the planetary wind velocity and temperature. Hence, for this model, the total absorption (⁠|${\sim } 10\,{{\rm \,per\,cent}}$|⁠) from the Ly α line is less than what observations show.

When B_⋆ is higher, a more compressed planetary magnetosphere is formed, and the neutral material piles up close to the planet. Our models show that for these cases, the total absorption in the Ly α line is about |$13\,{{\rm \,per\,cent}}$| during transit, but absorption of |$30\,{{\rm \,per\,cent}}$| is found for velocities around −100 km s⁻¹ and beyond, contradicting the observations.

We confirm that magnetic fields are an important aspect of the stellar and planetary wind interaction. Comparing with our previous work (Villarreal D’Angelo et al. 2014; Schneiter et al. 2016), they help to better reproduce the observed absorption in the red wing of the Lyman α line. The planetary magnetic field can have an influence in stopping the stellar wind, but only for values higher than or equal to 5 G will it also influence the amount of neutral material present in the planetary wind. It is encouraging to see that for our initial set of parameters chosen the model that best reproduces the Ly α observations of Vidal-Madjar et al. (2003) suggests that the magnetic field of HD 209458b could be less than or equal to 1 G (in agreement with Kislyakova et al. 2014). We must note that this estimation depends on a large number of model parameters, and a more self-consistent modelling of the generation of the planetary wind is necessary in order to confirm such aconclusion.

Footnotes

ACKNOWLEDGEMENTS

The authors thank the referee for the comments and suggestions that helped to improve this work. CVD acknowledges STFC (ST/M001296/1) and the PhD fellowship from CONICET. MAS acknowledges support from the CONICET via an Assistant Research Fellowship. The authors want to thank the Centro de Simulación Computacional para Aplicaciones Tecnológicas - CONICET, where some of the test of the code was run. AE acknowledges support from CONACYT DGAPA-PAPIIT (UNAM) grants IN 109715 and IG 100516.

REFERENCES

APPENDIX A: MHD IMPLEMENTATION ON GUACHO

We present in this work the new version of the code guacho, which includes the equations for solving the ideal MHD problem in a 3D Cartesian grid. The guacho code is a free access code and it is maintained in http://github.com/esquivas/guacho, from where it can be downloaded.

As we state in Section 2, the set of equations (1)–(4) is advanced in time with a second-order Godunov scheme. The code allows the user to choose between two approximate Riemann solvers, the HLLE solver (Einfeldt et al. 1991, quite robust, but somewhat diffusive), and the less diffusive HLLD of Miyoshi & Kusano (2005). To achieve second-order accuracy, the method performs a linear reconstruction of the primitive variables to solve the Riemann problem at each cell interface. To ensure monotonicity a number of slope limiters for the linear reconstruction step are available; the work presented here makes use of the minmod (the most robust/diffusive).

In order to keep the divergence of the magnetic field equal to zero, two methods have been implemented in the code: the eight-wave method proposed by Powell et al. (1999) and the method of Constrained Transport/Central Difference proposed by Tóth (2000). The former considers the MHD equations keeping the terms with |$\nabla \cdot \mathbf {B}$|⁠, which is zero theoretically, but slightly non-zero when calculated numerically. Such terms are treated as sources, allowing the spurious magnetic field divergence to behave as a scalar which is transported out of the grid without amplification. The second method is more elaborated, including an additional step after advancing the equations to correct the magnetic field, and keeps divergence down to machine precision.

To verify that the numerical scheme implemented in guacho behaves correctly, we reproduce several 1D, 2D and 3D tests that are commonly used in the literature to validate MHD codes. All of them can be compared with results from other MHD codes that are widely used, such as athena (Stone et al. 2008) or pluto (Mignone et al. 2007).

Brio-Wu 1D

This 1D test, proposed by Brio & Wu (1988), is a MHD version of the Sod shock tube hydro-dynamical problem. It is a suitable problem to check if MHD codes can accurately resolve shocks, rarefaction, contact discontinuities and other simple MHD structures.

The initial conditions are set at t = 0. The entire domain x ∈ (−0.5, 0.5) is filled with a gas with γ = 2 in two uniform states, separated by a discontinuity at x = 0. The left side of the domain is initialized to (ρ, |$v$|_x, |$v$|_y, |$v$|_z, B_x, B_y, B_z, p) = (1, 0, 0, 0.75, 1, 0, 1) while in the right side, these variables are chosen to be (0.125, 0, 0, 0.75, −1, 0, 0.1). The boundary condition at the outer limits of the domain is reflective.

Fig. A1 shows the final estate of the simulation at t = 0.1 using the HLLD solver and the eight-wave method to constrain |$\nabla \cdot \mathbf {B}$|⁠. The blue lines represent the results for a grid size of 10 000 cells, while the red points show the results for a grid of 400 cells. Our results can be directly compared with fig. 2 of the original work by Brio & Wu (1988). They can also be compared with the results shown in fig. 13 from the work of Stone et al. (2008), where the authors implemented this test with the same two resolutions used here into the athena code. Our code yields basically the sameresults.

Figure A1.

Simulation variables as a function of the position for the Brio–Wu test at t = 0.1. Blue lines represent the results for a grid of 10 000 cells while red dots show the corresponding results for a grid of 400 cells.

Open in new tabDownload slide

Kelvin-Helmholtz 2D

A second test is done in a setup with a shear flow that excites the Kelvin–Helmholtz (Helmholtz 1868; Lord Kelvin 1871) hydro-dynamical instability (and a version modified by magnetic fields). The test consists of a square domain, x ∈ [−0.5, 0.5] and y ∈ [−0.5, 0.5], filled with two fluids of different densities moving in opposite directions. In the region |y| < 0.25 we set ρ = 2 and |$v$|_x = 0.5, while in the region |y| ≥ 0.25 we set ρ = 1.0 and |$v$|_x = −0.5. In the entire domain the pressure has a value of 2.5 and γ = 1.4. The boundary conditions are periodic everywhere.

To initialize the instabilities a random component to the |$v$|_x and |$v$|_y variables was added, with a peak-to-peak amplitude of 0.01. For the MHD case, the x component of the magnetic field is set to 0.5.

Fig. A2 shows the density contours obtained at t = 1, 3 and 5 for a 512 × 512 mesh. The left column shows the hydrodynamic case while the right column corresponds to the magneto-hydrodynamic case. The test was also run using the HLLD scheme and the eight-wave method. Our results can be compared with the results presented in the athena web page for both cases (HD and MHD) https://www.astro.princeton.edu/jstone/Athena/tests/kh/kh.html.

Figure A2.

Density contours from the Kelvin–Helmholtz test. The left column shows the density for the HD case, while the right column corresponds to the MHD case. The rows show the outputs at t = 1, 3 and 5, respectively.

Open in new tabDownload slide

Orszag–Tang 2D

This vortex system was originally proposed and studied by Orszag & Tang (1979). Given that the flow develops a complex structure, it is now used extensively as a test to verify the ability of MHD codes to solve turbulence and shocks (see Tóth 2000; Stone et al. 2008, and references therein).

Fig. A3 shows the density distribution at time t = 0.5 for a 512 × 512 grid. The result can be compared with those presented in Tóth (2000) and Stone et al. (2008). In this test, to keep the |$\nabla \cdot \mathbf {B} = 0$|⁠, we used the method of Constrained Transport with the HLLDsolver.

Figure A3.

Density contours at time t = 0.5 from the Orszag–Tang test for a 512×512 grid.

Open in new tabDownload slide

Ryu–Jones 3D

The aim of this test, proposed by Ryu & Jones (1995), is to illustrate the ability of the numerical scheme to resolve the presence of fast and slow shocks, as well as rotational discontinuities. It consists in a MHD shock tube with a left state |$(\rho , v_x, v_y, v_z, B_x, B_y, Bz,E) = [1.08, 1.2, 0.01, 0.5, 2/\sqrt{4\pi }, 3.6/\sqrt{4\pi }, 2/\sqrt{4\pi }, 0.95)]$| and a right state |$[1.0, 0, 0, 0, 2/\sqrt{4\pi }, 4/\sqrt{4\pi }, 2/\sqrt{4\pi }, 1)]$|⁠. Outflow boundary conditions are used in all the domain and we set γ = 5/3.

The test was run using the HLLD scheme and the eight-wave method to maintain the |$\nabla \cdot \mathbf {B}$| close to zero. The results for t = 2.2 are shown in Fig. A4 for a mesh size of 768 × 8 × 8. The size of the mesh was chosen to be able to directly compare our results with those presented in fig. 14 of Stone et al. (2008).

Figure A4.

Simulation variables as a function of the position for the Ryu–Jones 2a test. All the variables are shown for t = 2.2 in a grid of 768 × 8 × 8 cells.

Open in new tabDownload slide