Modeling of Absorption by Heavy Minor Species for the Hot Jupiter HD 209458b

The absorption of stellar radiation observed by HD 209458b in the resonant lines of O i and C ii has not yet been satisfactorily explained. We apply a 2D hydrodynamic multi-fluid model that self-consistently describes the expanding planetary wind, driven by stellar XUV radiation and influenced by tidal forces and the surrounding stellar wind. According to this model, HD 209458b has a hydrogen-dominated plasmasphere, expanding beyond the Roche lobe, in the form of two supersonic streams that propagate toward and away from the star. The species heavier than hydrogen and helium are dragged in the escaping material streams and accelerated up to 50 km s−1. Our simulations show that, assuming solar abundances, O i and C ii produce absorption due to the Doppler resonance mechanism at the level of 6%–10%, which is consistent with the observations. Most of this absorption takes place in the streams. The transit depth in the O i and C ii lines is unaffected by the stellar wind, unless it is strong enough to form a compact bowshock around the planet and able to redirect all the escaping material to the tail. In this case, the absorption profile becomes asymmetric due to the prominent blueshifted attenuation. Thus, the spectroscopic measurements enable probing of the planetary wind character, as well as the strength of the stellar wind. The computed absorption at wavelengths of the Si iii, Mg i, and Mg ii lines at solar abundances appears to be much stronger, compared to the observations. This possibly indicates that Si and Mg may be under-abundant in the upper atmosphere of HD 209458b.


Introduction
The study of planets beyond the solar system (exoplanets) is one of the fastest growing and intriguing fields in astrophysics. Just after 20 years of the first detection of an exoplanet orbiting a solar-like star (Mayor & Queloz 1995) more than 3800 exoplanets, including about 600 multiple systems, have been discovered (e.g.,http://exoplanet.eu/catalog/). Among those more than 2900 are transiting exoplanets (e.g., NASA exoplanet archive https://exoplanetarchive.ipac.caltech.edu/docs/counts_ detail.html), which allow one to detect and characterize planetary atmospheres using spectral analysis methods (Brown et al. 2001;Charbonneau et al. 2002). Also, because of observational biases, the majority of the known exoplanets orbit rather close to the host star, some of them at a distance shorter than 0.05 au. Lammer et al. (2003) were the first to show that the hydrogen-dominated atmosphere of such planets is heated to several thousand Kelvin and expands dynamically. One of the most studied exoplanets is the transiting hot Jupiter HD 209458b, which presents an optical transit depth of about 1.5% (Charbonneau et al. 1999). Vidal-Madjar et al. (2003) observed primary transits of this planet in Lyα with the STIS spectrograph on board HST and reported a 15% absorption signal, located mostly in the high velocity blue wing of the line. To explain this effect, the Doppler resonant absorption caused by the presence of Energetic Neutral Atoms (ENAs) of planetary origin moving away from the star with velocities up to 150 km s −1 (Holmström et al. 2008;Kislyakova et al. 2014) was proposed. However, subsequent re-analyses of the same data led to a slightly smaller and more symmetric absorption of 6%-9% (Ben-Jaffel 2007; Vidal- Ben-Jaffel & Sona Hosseini 2010), which corresponds to the presence of a planetary escaping atmospheric material up to 2-2.5 planetary radii, hence close to the Roche lobe. This reconsidered result can be explained mostly by the natural line broadening. That was also recently confirmed with a self-consistent numerical modeling of the HD 209458b plasma environment and related Lyα in-transit absorption performed by Khodachenko et al. (2017). Soon after the discovery of hot Jupiters, 1D gas-dynamic simulations began to model the dynamics of planetary upper atmospheres to get deeper insights into the mechanisms of the measured transiting spectra (see e.g., Shaikhislamov et al. 2014, and references therein). The results of these simulations showed that the extreme heating of the planetary upper atmosphere, due to absorption of the ionizing stellar XUV radiation, leads to a supersonic outflow of planetary gas in the form of an escaping planetary wind (hereafter PW), that overcomes planetary gravitation and expands beyond the Roche lobe. All of these simulations are in agreement regarding an average exospheric temperature of about 10 4 K, an outflow velocity of about 10 km s −1 , and a mass-loss rate of the order of 10 9 -10 11 g s −1 .
The first explanation of the Lyα transit observations of HD 209458b was based on the assumption that the hydrogen atoms in the escaping PW are accelerated to high velocities by the radiation pressure (Vidal-Madjar et al. 2003;Lecavelier des Etangs et al. 2004, 2008. However, the studies of García Muñoz (2007), Ben-Jaffel (2007, 2008, and Koskinen et al. (2010), based on 1D hydrodynamic models, led to the conclusion that the absorbing planetary material more likely lies within the Roche lobe, and that the absorption over a broad emission profile is caused by the natural line broadening. The same conclusion was also made by Shaikhislamov et al. (2016) and Khodachenko et al. (2017), who used 2D multi-fluid simulations to treat self-consistently the generation of the PW and its interaction with the stellar wind (hereafter SW) plasma. They showed that if the stellar XUV irradiation and the parameters of the SW are similar to those of the Sun, the thermal line broadening or resonant absorption caused by ENAs is small. However, it has also been shown that there might be conditions when the role of ENAs becomes significant. In particular, for XUV irradiation, typical for solar minimum, or in the case of a significantly stronger SW, the number of ENAs increases, and this may lead to an additional 5%-10% absorption in the blue wing of the Lyα line (Khodachenko et al. 2017).
Besides Lyα, the primary transit of HD 209458b has been also observed with HST/STIS at other far-UV (FUV) wavelengths. In particular, Vidal-Madjar et al. (2004) reported absorption depths of 13%±4.5% and 7.5%±3.5% in O I (2p 2 2P-2p 2 2D) and C II (2p 4 3P-2p 4 3S) FUV resonance lines, respectively. Further analysis of the same data set by Ben-Jaffel & Sona Hosseini (2010) led to similar results (10.5%±4.4% and 7.4%±4.7% for O I and C II, respectively). According to Linsky et al. (2010), the analysis of four additional FUV primary transit observations carried out with the COS spectrograph, on board HST, confirmed the previously measured C II transit depth of 7.8%±1.3%. 7 Linsky et al. (2010) also reported the detection of an 8.2%±1.4% absorption signal at the position of the Si III (1206.5 Å) resonance line in contrast with Vidal-Madjar et al. (2004), who reported a non-detection. It should be noted that due to possible strong stellar variability, Ballester & Ben-Jaffel (2015) called into question the detections of C II and Si III absorption features with COS. Besides that, Vidal-Madjar et al. (2013) reported results obtained from three near-UV (NUV) transit observations of HD 209458b carried out with HST/STIS. The data revealed the presence of an absorption signature of 6.2%±2.9% in the Mg I resonance line, while no absorption was detected in the Mg II h&k lines. This is a surprising result, because Mg II absorption has been previously detected for another hot Jupiter, WASP-12b (Fossati et al. 2010;Haswell et al. 2012), and Mg is known to be easily ionized by the stellar FUV flux. These observational results suggest the presence of heavy species in the upper atmosphere of HD 209458b. If the Lyα absorption is dominated by natural line broadening and caused by exospheric gas lying inside the Roche lobe, the absorption by heavier species, with abundances three-to-five orders of magnitude lower than hydrogen, offers the possibility of tracing the planetary gas beyond the Roche lobe. Vidal-Madjar et al. (2004) estimated that, assuming the solar system abundances, 8 the observed O I absorption can be explained by the presence in the planetary exosphere of O I atoms with a density of at least 10 6 cm −3 . García Muñoz (2007), using a 1D aeronomic model, showed that such density is in fact a realistic one, and that the optical depth for the O I and C II line cores averaged over the Roche lobe is larger than ] V 15 km s 1 ∼3.5). However, the 1D model could not explain the required line broadening width up to 40 km s −1 . For example, Ben-Jaffel & Sona Hosseini (2010), using the density distributions provided by García Muñoz (2007), obtained absorption signatures of just 3.9% and 3.3% for the O I and C II multiplets, respectively. A number of processes have been suggested to account for this discrepancy. Ben-Jaffel & Sona Hosseini (2010) suggested preferential heating of minor heavier species inside the Roche lobe to explain the absorption and concluded that the effective temperatures should be about 10 times larger for oxygen and 5 times larger for carbon particles than for hydrogen. However, the existence of such a large temperature difference is hardly possible for the entire population of heavier species, mostly because the hydrogen density inside of the Roche lobe is higher than 10 7 cm −3 and the partially ionized plasma is strongly collisional. Tian et al. (2005) proposed that sufficient broadening of the oxygen lines could be achieved if the bulk velocity of about 10 km s −1 is added to the thermal motion of the atoms. Koskinen et al. (2010) calculated the absorption by O I atoms based on a 1D hydrostatic model, obtaining an absorption signature of 4.3%. The subsequent addition of the radial velocity of 10-30 km s −1 to the static solution then increased the absorption up to 5.8%-7.2%. A number of other factors, which could increase the O I absorption, were further investigated by Koskinen et al. (2010), including sporadic hot spots on the stellar disk, limb darkening or limb brightening, super-solar abundances, and the position of the H 2 /H dissociation front. However, none of the considered effects could provide an increase of the O I transit depth, comparable to, or higher than that of H I. Additional computations, based on the same model, led to the absorption of 5.2% and 4.6% for C II and Si III, respectively (Koskinen et al. 2013b). They also showed that two times higher stellar XUV flux leads to a significantly higher outflow velocity that increases the C II transit depth up to 6%.
Altogether, the previous works, based on 1D hydrodynamic models, hitherto were not able to find reasonable conditions within which the simulated O I, C II, and Si III absorptions are comparable to (or even exceed) those observed in Lyα. The main point here is that, while Lyα absorption can be explained by natural line broadening with a sufficiently large number of hydrogen particles inside the Roche lobe, the density of heavier species is too small for that, because the typical temperature of 10 4 K is too low to generate thermal broadening over typical line widths of 45-60 km s −1 . However, in the present paper we show that the necessary broadening may be generated by a global flow of planetary material outside the Roche lobe, as was qualitatively proposed by Tian et al. (2005).
The study reported below is based on a fully self-consistent 2D multi-fluid numerical model (Shaikhislamov et al. 2014(Shaikhislamov et al. , 2016Khodachenko et al. 2015Khodachenko et al. , 2017) that does not rely on any empirical simplifications/assumptions regarding the motion of the planetary atmospheric material, and properly combines the micro-physics of the PW generation with its large-scale expansion and interaction with the SW far beyond the Roche lobe. As shown in Shaikhislamov et al. (2016), depending on 7 We note that due to contamination of the COS spectra by the geocoronal emission, no measurements for O I were possible. 8 For the solar system abundances, we suppose (and use) in this paper the element constituents relative to hydrogen, which are presumably typical for the solar system protoplanetary disk and confirmed with aeronomic modeling of HD 209458b (e.g., Burrows & Sharp 1999): O/H=8.5×10 −4 , C/H=3.6×10 −4 , Si/H=3.7×10 −5 , Mg/H=3.7×10 −5 . These values are also suggested by Anders & Grevesse (1989), but at the same time, they admit certain deviations (Lodders 2003;Asplund et al. 2009). Note that these solar system abundances' values are not the same as those estimated for the Sun itself, which are called solar abundances. the parameters of a stellar-planetary system and the corresponding intensity of the SW and the escaping PW, hot Jupiters can have two essentially different regimes of material outflow. The first one, called "blown by the wind," occurs in the presence of a sufficiently strong SW, that stops the escaping planetary material on the day-side and channels it away from the star, confined within a kind of paraboloidal structure with a bowshock, ionopause, and a tail elongated along the SW. The second regime, called "captured by the star," occurs when the action of the tidal force overcomes the SW pressure, forming a double stream structure of the escaping planetary material directed toward and away from the star. One of the main conclusions following from the simulations in Shaikhislamov et al. (2016) is that the flow of the escaping planetary material remains sufficiently dense to be strongly collisional even very far (up to 30 planetary radii, R p ) from the planet. This model has been further applied in Khodachenko et al. (2017) to interpret the Lyα transit observations of HD 209458b, while in the present paper a similar approach is followed to simulate the measured spectral absorption features of heavier species. As it will be shown later on in the paper, the "blown by the wind" and the "captured by the star" regimes, together with the account of the collisional character of the escaping PW and its interaction with the entire SW plasma, have important implications for the blue/red shape of the transit spectra.
The paper is organized as follows. In Section 2, we briefly describe the model used for the simulation of the PW and SW interaction, paying attention to the basic processes, important details regarding the calculation of absorption, and key parameters of the O I, C II, Si III, Mg I, and Mg II spectral lines of interest. In Section 3 the absorption in O I, C II, Si III, Mg I, and Mg II lines is calculated for the modeled interaction of the PW and SW by HD 209458b under different conditions and for different abundances, then the results are compared to the available measurements. Section 4 presents a discussion of the obtained results and our conclusions.

The Numerical Model and Absorption Parameters
The details of our numerical model have been described in Shaikhislamov et al. (2014Shaikhislamov et al. ( , 2016 and Khodachenko et al. (2015). Here, we provide just a short description of the most important characteristics of the code. We employ a 2D axially symmetric hydrodynamic numerical model in the cylindrical coordinate system with the symmetry axis, z, taken along the planet-star line, and the reference frame attached to the planet. Such geometry is well-suited for the simulation of tidally locked non-magnetized planets, with the stellar illumination coming from only one (dayside) direction. Within this approximation, we disregard the Coriolis force and ignore the lateral component of the SW plasma due to the planetary orbital velocity. The plasma is treated as a quasi-neutral fluid in thermal equilibrium with equal electron and ion temperatures. The modeled HD 209458 system consists of a planet with mass M p =0.71 M J and radius R p =1.38 R J orbiting at a distance of D=0.047 au around a solar-type main-sequence G-star with a mass M st =1.148 M Sun and a radius R st =1.13 R Sun . The planetary and stellar radii and masses are given in the units of radius and mass of the solar system Jupiter (R J =71.5×10 3 km, M J =1.89×10 27 kg) and the Sun (R Sun =69.6×10 4 km, M Sun =1.99×10 30 kg), respectively. As the initial state, we consider a fully neutral atmosphere in a barometric equilibrium composed of molecular hydrogen and neutral helium. At the inner boundary of the simulation domain r=R p we fix the temperature 1000 K, the pressure p=0.03 bar, and the mixing ratio x He /x H2 =1/5. We employ a multi-fluid approach, which includes H, H + , He, He + . To treat the thermosphere more accurately, the model considers also the + H 2 , + H 3 molecular ions. The ENAs, formed by charge-exchange between the planetary atomic hydrogen and SW protons, are also treated as a separate fluid. Any minor species, as well as the SW (composed of H + ) are added to these components as separate fluids. However, we do not consider chemical reactions among the different minor species. As the characteristic values of the numerical problem, we use those according to our previous simulations (Shaikhislamov et al. 2014(Shaikhislamov et al. , 2016Khodachenko et al. 2015Khodachenko et al. , 2017, the planet radius R p , the temperature 10 4 K, and the corresponding thermal velocity of ions V o =9.07 km s −1 . The obtained numerical solutions are correspondingly scaled in these units. The processes affecting the gas particles are photo-ionization, excitation, electron impact ionization, dielectronic recombination, and charge-exchange. Photo-ionization also results in a strong heating of the planetary material, which is the main driver of the PW. The stellar XUV flux is assumed to be equal to that of the Sun (Tobiska 1993) and covers the 10-912 Å spectral range, binned with 1 Å steps. The spectrum is based on measurements of the solar radiation under the conditions of moderate activity with index P 10.7 =148. The integrated XUV flux at 1 au of this spectrum is F XUV =4.466 erg s −1 cm −2 . The XUV photons ionize hydrogen atoms, as well as H 2 and He. The model assumes that the energy released in the form of photo-electrons is rapidly and equally redistributed between all locally present particles with an efficiency of η h =0.5. The resulting total efficiency calculated as the heating per absorbed radiation flux is about 8%. The attenuation of the XUV flux inside the planetary atmosphere is calculated for each spectral bin according to the wavelength dependent cross-sections.
The minor species are added to the initial atmosphere assuming solar system abundances (Anders & Grevesse 1989) and evolve according to the dynamic equations. The population of different ionization states for each element is calculated assuming the specific photo-ionization (Verner et al. 1996) and recombination rates (Nahar & Pradhan 1997;Le Teuff et al. 2000). Besides photo-ionization and recombination, the population of atomic oxygen is affected by the resonant charge-exchange reaction with hydrogen ( + « + The corresponding crosssection, of the order of 10 −15 cm 2 (Le Teuff et al. 2000) is sufficiently large, so that oxygen ionization follows the ionization of hydrogen (+ Carbon atoms are ionized in the lower thermosphere, close to the conventional planetary surface in the model, at about 1.1 R p (García Muñoz 2007), because the photo-ionization threshold for carbon is significantly smaller than that for hydrogen. Therefore, we add already ionized carbon atoms to the atmosphere and follow the dynamics of C II and C III, taking into account their further photo-ionization and recombination.
We also follow the same approach for silicon. In the upper atmosphere, it is present mostly as SiO (Visscher et al. 2010). The dissociation of SiO and further ionization of silicon take place relatively close to the planetary surface (Koskinen et al. 2013a) and therefore we follow the dynamics of only Si II and Si III ions. Similarly to oxygen, these components are in charge-transfer balance with hydrogen, meaning Because of the charge-exchange, the Si III abundance is relatively high. We allow also for ionization of Si III, but do not consider recombination of SiIV, which therefore leads to a possible minor underestimation of the Si III density.
Magnesium is instead treated in a slightly different way as compared to other elements. While the Mg I atom is rapidly photoionized, the photo-ionization of a hydrogen-like Mg II is slow (Verner et al. 1996) whereas the recombination is fast. Therefore, both components appear in the ionization-recombination balance. At typical temperatures of less than 10 4 K and electron densities of less than 10 9 cm −3 , Mg II prevails (Vidal-Madjar et al. 2013). However, there is also the resonant charge-exchange plus ionization reaction of Mg I with ionized Helium, which double ionizes Mg I, while skipping the Mg II stage: Mg+He + → Mg 2+ +He+e (DuBois 1986). At relatively high impact energies ∼40 keV it has a maximum cross-section of about 10 −16 cm 2 (DuBois 1986). Depending on the cross-section of this reaction at very low energies, which is currently unknown, Mg I and Mg II might rapidly deplete when the helium in the thermosphere becomes photo-ionized. We study quantitatively this effect below in Section 3.2.2.
For the typical parameters of the planetary exosphere, Coulomb collisions with protons effectively couple all ions. For example, at T<10 4 K and n H+ >10 6 cm −3 the collisional equalization time for temperature and momentum t » p is less than 1 s for protons and less than 10 s for oxygen and carbon ions (Braginskii 1965). This is several orders of magnitude less than the typical gas-dynamic timescale of the problem treated here, which is of the order of 10 4 s. There is also another physical reason for the strong coupling of charged particles in the PW on the considered typical spatial scale of the problem, that is about R p (∼10 10 cm). The chaotic/sporadic magnetic field, which is present in the PW, affects the relative motion of the ions so that they become coupled via the magnetic field due to the Lorentz force, and exchange with the momentum on a timescale of the Larmor period. Even at negligibly small fields of the order of 10 −3 G, the gyroperiod of an oxygen ion is about 1 s. For the same reason, the charged particles can be treated as strongly coupled ones in the hot and rarefied SW as well, even in spite of the fact that Coulomb collisions are negligible there. Therefore, there is no need to calculate the dynamics of every charged component of the plasma fluid species, and we assume in the simulations that all of them (i.e., H + , + H 2 , + H 3 , He + , O n+ , C n+ , Si n+ , Mg n+ ) have the same temperature and velocity. On the other hand, the temperature and velocity of each neutral component are calculated individually by solving the corresponding energy and momentum equations. The neutral hydrogen atoms are also more or less coupled to the main flow by elastic collisions. With a typical cross-section of >10 −16 cm 2 , the mean-free path at a density of 10 6 cm −3 is comparable to R p . Besides elastic collisions, charge-exchange ensures more efficient coupling between hydrogen atoms and protons (Shaikhislamov et al. 2016;Khodachenko et al. 2017). With respect to oxygen atoms the equalization proceeds in three steps. First, atoms are ionized through charge-exchange with protons. Second, ions are rapidly equalized with protons via Coulomb collisions. Third, by charge-exchange, oxygen ions are transformed back to neutral atoms. In the region of interest outside of the Roche lobe, protons dominate over atoms, so the third step is the one that takes the longest time. With a crosssection of about 10 −15 cm 2 , the corresponding timescale is an order of magnitude smaller than the typical gas-dynamic time. So, the oxygen atoms remain strongly coupled to the main flow. We therefore conclude that the minor species considered here are expected to be efficiently dragged by the planetary flow.
The spectral absorption is described by the so-called Voigt convolution of a Lorentz line shape with a natural width and Maxwellian distribution, characterized by the mean temperature T, component of velocity along the line of sight V z , and particle density n: is the Voigt profile, which takes into account the convolution of two broadening mechanisms, one of which alone would produce a Gaussian profile (as a result of Doppler broadening), and another would produce a Lorentzian profile. Here Instead of wavelength, the absorption is expressed in terms of Doppler-shifted velocity in the reference frame of the star. The averaging in Equation (1) covers the stellar disk of radius R St . The fluid values for n, T, V z are calculated by the numerical model, and for the length L along the line of sight, we take an empirical value of L=10 R p . This length is limited by the spiraling of planetary material flow due to the Coriolis force (clockwise in the planet-based reference frame), which is not included in our model. In other words, L corresponds approximately to the size of our modeling box, where the spiraling might still be neglected for the considered system. The part of the stream that is not significantly affected by the Coriolis force, has a length roughly equal to the stream width, which is also about 10 R p . This assumed size, at which the spiraling of the PW stream can be neglected, is also consistent with the large-scale 3D simulations by Kislyakova et al. (2014), Bourrier & Lecavelier des Etangs (2013).
At temperatures above 2 K, which is definitely true in our case, the Voigt integral can be fitted, with an accuracy better than 1%, by the following analytical expression (Tasitsiomi 2006): This analytical fit explicitly shows that the absorption consists of two components. The first one is the resonant absorption by atoms matching the Doppler-shifted velocity of the line profile. This process is characterized by a significant cross-section within the line core. When the average thermal velocity of atoms is comparable with velocities corresponding to a particular part of the line profile » m V 2 kT i then the number of resonant atoms with V z ≈V is statistically large due to the specifics of the usually considered thermal equilibrium Maxwellian distribution function. In this case, the produced absorption is related to the so-called thermal line broadening. Note that the heavier particles at the same temperature have smaller thermal dispersion of the distribution function, and therefore a smaller velocity range for the thermal line broadening. That is, in particular, the reason why it is difficult to explain the absorption in O, C lines only by thermal line broadening. Besides that, a similar resonant absorption will take place when a relatively high bulk velocity of gas motion appears in Doppler resonance with a part of the line profile. This produces absorption that is particularly relevant for this paper. Later on, for simplicity, we call the absorption produced by the first term and proportional to σ res "Doppler resonance absorption." The second term in (2) is due to the non-resonant process related to the natural line broadening caused by the finite lifetime of the absorbing agents affecting the far wings of the Lorentzian profile, which has a much smaller cross-section. The particular values for σ res and σ nat for the considered elements, obtained after substitution of the analytic approximation of the Voigt integral, H, into the general expression for the absorption cross-section, are given in Table 1. The similar formula, albeit differently normalized and without the correcting factor q(x 2 ), can be found in the literature (e.g., Bourrier & Lecavelier des Etangs 2013). Table 1 lists the properties and relevant parameters of the spectral lines considered in this paper. Among them are: the photo-ionization time for each element, depending on the ionization energy threshold and stellar XUV intensity; the measured absorption according to Ben-Jaffel & Sona Hosseini (2010) and Linsky et al. (2010), the statistical weights of the ground states for each line in the multiplet; the relative intensities within each multiplet measured out of transit; and the coefficients σ res , σ nat appearing in (2), in the terms related to Doppler resonance, natural line broadening absorption mechanisms, respectively. The intensity of the unresolved transition without orbital momentum change (3/2-3/2) in the C II multiplet is taken to be zero because its oscillator strength is an order of magnitude lower, as compared to that of the nearby transition 3/2-5/2. This results in a double decrease of the overall absorption of the multiplet. The C II and Si III lines are sufficiently well resolved in the observed spectra (Linsky et al. 2010). For these features, we employ the same analytical fits as those used by Koskinen et al. (2013b). Since the 1334.5 Å C II line corresponds to a ground-state transition, this feature is affected by the interstellar medium (ISM) absorption. Therefore, we do not consider the interval of about ±0.05 Å, i.e., ±10 km s −1 (Ben-Jaffel & Sona Hosseini 2010). One should note that this significantly reduces the overall absorption by the planetary exosphere. For the Si III resonant line the interstellar absorption is likely to be negligible (e.g., Fossati et al. 2015). One can see in Table 1 that in velocity units the half-width-at-half-maximum (HWHM) of the C II and Si III lines reaches a significant value of 35 km s −1 . For the O I, the intrinsic stellar emission line shapes are taken from the solar spectrum. This assumption is justified by the strong similarity between HD 209458 and the Sun (Shkolnik et al. 2005). The double peaked lines of the O I multiplet can be analytically approximated by a double Gaussian profile, with the corresponding parameters from Link et al. (1988). The line cores are not considered, because they are affected by the ISM absorption. Altogether, the shapes used for C II, O I, and Si III lines are shown in Figure 1. They are the same as those in Ben-Jaffel & Sona Hosseini (2010). In particular, according to the plot, the O I absorption supposes the presence of particles with velocities ranging from ±10 to ±25 km s −1 .

Absorption in O I and C II Lines Under Different
Conditions of the Interacting PW and SW Figure 2 shows the simulations of the PW escaping from HD 209458b and passing through the SW. It displays the density of hydrogen atoms n H , as well as the densities n O I and n C II of O I and C II particles, respectively. The assumed abundances of these elements at r=R p are O/H=8.5×10 −4 and C/H=3.6×10 −4 . Note that at r>R p the abundances evolve according to the solution of the model equations. The XUV flux taken for this simulation is 4 erg cm −2 s −1 at 1 au, while the SW was assumed to be very weak, with n SW =1 cm −3 , V SW =45 km s −1 , T SW =10 5 K, so that it does not affect the numerical solution for the outflowing PW. 9 The simulation run under this set of parameters will be referred to later on as the model M1. Shaikhislamov et al. (2016) and Khodachenko et al. (2017) have shown that in the case of neglecting the Coriolis force, the expanding planetary material forms a nearly symmetric double stream structure. Beyond the supersonic transition point located outside the Roche lobe at about 5 R p , the streams of the escaping PW are further accelerated by the tidal force. The width of the stream of outflowing planetary atmospheric material is controlled by the dynamical and thermal pressure balance, whereas the speed of the stream constantly increases beyond the Lagrange point due to the tidal acceleration. Because of that, the width of the streams tends to decrease with the distance from the planet.
The left and middle panels in Figure 3 show the profiles of densities across and along the planet-star line at a fixed position z=0 and at the radial position r=2.5 R p , respectively. One can see that, because of resonant charge-exchange, O I and O II densities have similar (but not identical) profiles as those of H and H + . The density of C II generally follows that of O I. These features are similar to other aeronomic simulations (e.g., García Muñoz 2007). The right panel of Figure 3 shows the profiles of velocity and temperature for oxygen atoms and protons along the planet-star line at the radial position r=2.5 R p , which appear to be rather close to each other. This fact indicates that oxygen is indeed strongly coupled to hydrogen due to charge-exchange, and that the slippage of species is very small. Thus, the oxygen is dragged with the main flow and attains velocities of tens of km s −1 .

Doppler Resonance Versus Natural Line Broadening Absorption Mechanisms in O I and C II Lines, and the Role of the Escaping PW Streams
The absorption produced by O I, carried in the PW flow, is shown in Figure 4 with the line O I 1302.17 Å, taken as an example. This line produces the largest contribution to the whole absorption of the oxygen multiplet because of its highest statistical weight. The calculations of the absorption were performed under the conditions of the simulation run M1. The decomposition of the O I 1302.17 Å absorption shows that its resonant part constitutes a dominant portion of the total (5÷15)% absorption and covers almost the whole line from −30 km s −1 in the blue wing up to 15 km s −1 in the red wing. The asymmetry of the absorption profile is caused by the smaller ionization of oxygen and hydrogen at the night-side, where the stream moves toward the observer (see the middle panel of Figure 3). Note that the resonant part of the absorption has a cutoff at about ±35 km s −1 , caused by the assumed integration limit within L=±10 R p . However, it is outside the actual line width (see Figure 4, right panel) and does not affect the overall absorption of the stellar emission. The same is also true for other lines considered in the paper. For the sake of further comparison, we will call such calculations of the resonant part of the absorption, performed under the same conditions as those considered in M1, model M2. The transit depth of the whole O I multiplet, calculated with the modeled distributions of species, the abundance O/H=8.5×10 −4 , ground level population 5/9, and emission intensities typical for HD 209458b (e.g., F XUV at 1 au 4 erg cm −2 s −1 ), gives a value of 8.5%.
To evaluate the role of the motion of the escaping PW material in the streams and to make a direct comparison with Ben-Jaffel & Sona Hosseini (2010) and Koskinen et al. (2010), who considered the absorption due to only natural line broadening, we calculated the absorption of the O I multiplet, assuming zero velocity for the planetary material flow in Equation (2). The simulation run under these conditions later on will be referred as model M3. This modeling case for the O I multiplet gave the absorption value of 4.3%. This result coincides with that of Koskinen et al. (2010). Its comparison with the previously considered case, which included the effect of escaping PW streams, allows one to conclude that the material motion in Doppler resonance with the corresponding parts of the emission line profile gives a significant contribution of more than 4% to the absorption. Another question concerns the contribution of the region outside the Roche lobe. The calculation of the O I absorption within r<3 R p and L=±3 R p , i.e., inside the Roche lobe, which will be further referred to as model M4, yields 5.7%. This includes the above-mentioned non-resonant part (4.3%) of the absorption, due to natural line broadening, and 1.4% of the resonant part of the absorption by O I within the Roche lobe. Therefore, the major part of the resonant absorption is due to the O I dragged in the PW streams outside the Roche lobe. It should be noted that the calculation results referred to as models M2-M4 are in fact the post-simulation analysis, performed for the same modeling run, M1. They differ only by the conditions and the way in which the absorption is calculated.
The left panel of Figure 5 shows the absorption profile for the C II 1335.71 Å line. It is more symmetric than that for the O I 1302.17 Å line, because the population of C II is much less affected by photo-ionization than the oxygen. The total C II multiplet absorption of 8.7%, calculated for the M1 model case, an abundance C/H=3.6×10 −4 , a ground level population 4/10, and stellar irradiation level typical for HD 209458b (e.g., F XUV at 1 au 4 erg cm −2 s −1 ), is practically equal to that of O I despite the smaller abundance. This is because the C II lines are broader (see the right panel of Figure 5) and the resonant absorption in the streams covers a larger range. It is confirmed by a dedicated calculation of the absorption with the assumption of no PW flow (i.e., for the model M3) that gives 3.8%, which is smaller than that for the O I multiplet, and comparable to the result of Ben-Jaffel & Sona Hosseini (2010), though smaller than that of Koskinen et al. (2013b). The C II absorption within the Roche lobe for r<3 R p and L=±3 R p , (i.e., for model M4) is equal to 4.4%. It appears also lower than that obtained for the O I multiplet. Therefore, the remaining part of the absorption for the C II multiplet (about 4.3%), which is mainly caused by the resonant absorption (i.e., the case of model M2) in the PW streams, is also rather pronounced.

The Influence of Abundances, XUV Flux, and SW Parameters on the Absorption in O I and C II Lines
Next, we study the dependence of the O I and C II absorption on the abundances, stellar XUV flux, and SW parameters. For the sake of comparison, we use the same abundances χ at r=R p as in García Muñoz (2007; O/H=8.5×10 −4 , C/H=3.6×10 −4 , Si/H=3.7×10 −5 , Mg/H=3.7×10 −5 ), which are based on the calculation of major atmospheric molecular constituents at the pressure of 1 atm and temperature of 1200 K (Burrows & Sharp 1999). These values are equal to the solar system abundances (relative to hydrogen), suggested by Anders & Grevesse (1989). However, later they have been revised toward lower values (Lodders 2003;Asplund et al. 2009). Since the minor elements practically do not affect the upper atmospheric material outflow, the absorption at different χ can be calculated using the same simulation run M1. At the same time, for the study of the effect of different XUV and SW conditions, new specific simulations were performed. Figure 6 shows the dependence of absorption in O I and C II multiplets, as well as in Si III, Mg I, and Mg II resonance lines on the corresponding abundances χ of these elements relative to hydrogen. It is slightly weaker for oxygen ( c O 0.4 I ) than for carbon ( c II C 0.5 ). One can see, in particular, that the modeled absorptions match with the values derived from the observations (Ben-Jaffel & Sona Hosseini 2010) for oxygen-at the abundance of 2÷2.5 times higher (O/H=1.5×10 −3 ) than the solar system value, whereas for carbon the corresponding abundance is practically the same as that in the solar system (C/H=2.6×10 −4 ). Therefore, these figures give an essentially lower C/O abundance ratio than the solar one (of about 0.2). However, taking into account the measurement error intransit depth of about ±4%, a conservative conclusion from our modeling could be a ratio C/O<1. This is in agreement with the retrieval atmospheric models, based on infrared transit spectroscopy (Line et al. 2016), though the question remains controversial (MacDonald & Madhusudhan 2017).
As to the influence of XUV flux, the C II absorption changes with F XUV as expected. In particular, for the simulation run under the same conditions as considered for the model M1, but with twice as large XUV flux at 1 au, F XUV =8 erg cm −2 s −1 , which will be referred to later on as model M5, the C II absorption increases by 1.29 times up to 11.2%. Two times decreased XUV flux, F XUV =2 erg cm −2 s −1 at 1 au (later on referred to as model M6) results in a 0.82 times decrease of the transit depth down to 7.2%. This is because the formation of the planetary material outflow and its velocity are driven by the XUV heating, which directly affects the absorption in C II lines. That has also been noted by Koskinen et al. (2013b). In contrast, the absorption by O I has much less variation. For the model M5, it increases only 1.06 times. This is because the O I density is strongly dependent on the density of atomic hydrogen, which varies inversely with respect to XUV flux, due to photo-ionization.
To study the effect of SW on the absorption features of the considered minor species, we include in the simulation two cases with a slow and fast SW, which may better correspond to the reality of HD 209458b. The interaction between the expanding PW and SW plasma flows, and the corresponding Lyα transit effects in the case of HD 209458b have been studied in detail in our previous works (Shaikhislamov et al. 2016;Khodachenko et al. 2017). The addition of minor species considered in the present paper does not change the overall picture of the interacting PW and SW. At the conditions of a slow SW with the total pressure p SW =5×10 −6 μbar, realized for n SW =4600 cm −3 , V SW =230 km s −1 , T SW =1.27 MK, that will be referred to later on as model M7, the transit depth for the O I multiplet increases, as compared to the case of a weak SW (model M1), from 8.5% to 9.2%, while that for C II does similarly increase from 10% to 10.7%. This increase can be explained by the lateral compression of the escaping PW streams by the increased thermal pressure of the SW, as compared to the case of the model M1, which results in the increase of material density in the streams. When the total pressure of the SW is sufficiently high, it stops and compresses the day-side planetary stream and forms a bowshock around the planet. In particular, at a higher pressure, p sw =1.3×10 −4 μbar, realized under the fast SW conditions with n sw =2.5×10 4 cm −3 , V sw =500 km s −1 , T sw =2.90 MK (will be referred to later on as model M8), the sub-stellar position of the bowshock, obtained in the simulations (for details see Khodachenko et al. 2017), is at about 4.3 R p , while the ionopause appears inside the Roche lobe at about 3.3 R p . For all the absorption lines considered in this paper, the transit depths slightly decrease under the fast SW conditions (model M8) as compared to the slow SW case (model M7), because the resonant absorption in the red part of the lines disappears. However, the absorption in the blue part of the lines in the case of model M8 significantly increases. For the O I multiplet the total absorption is 7.9%, while the absorption in the blue part is 12.1% (with 14%, 13%, and 9.9% for each line in the mulitplet). Correspondingly, the transit depth in the red part is only 4.2%. It is produced by oxygen inside the Roche lobe, and natural line broadening is the dominating absorption mechanism in this case. For the C II multiplet the total absorption in the case of model M8 is 7.6%, while the averaged absorption in the blue part is 11.9%. Correspondingly, the transit depth in the red part is only 3.8%, and it is generated by carbon ions inside the Roche lobe.
To summarize the obtained results, we can say that the major contribution to the total absorption in all the considered lines is due to the dominating Doppler resonance mechanism. The Doppler shift in absorption is connected with the material bulk motion in the PW streams escaping beyond the Roche lobe. The increased XUV flux that powers the material escape basically increases the absorption, whereas the stronger SW, which changes the "captured by the star" regime of the material escape to the "blown by the wind" regime decreases the absorption mainly due to the suppression of the PW stream in the stellar direction (for details see Table 2).

Absorption in Si III and Mg I, Mg II Lines under Different
Conditions of the Interacting PW and SW For the absorption in Si and Mg lines we obtained significantly different results, as compared to the absorption in O I and C II lines (see also in Table 2). In spite of Si and Mg abundances being an order of magnitude smaller than the abundances of oxygen and carbon, very large resonant absorption cross-sections for these elements (see σ res in Table 1) result in rather high absorption values.

On the Si III Absorption Features
The silicon ion density profiles along (at r=2.5 R p ) and across (at z=0) the planet-star line, obtained for the simulation run M1, are presented in Figure 7. They show, in particular, that the typical density of Si III varies within the range 10÷100 cm −3 . At these values the corresponding optical depth is rather high: n Si III σ res R p >0.1.
Moreover, the core of the Si III line, being unaffected by the ISM, should also give a significant input (of several percent) to the total absorption in the whole line. Under the conditions of the simulation run M1, presented in Figures 2 and 3, the transit depth reaches values of about 40%, and it is practically all due to resonant absorption by Si III ions moving in the PW streams (see in Figure 8). The specific plateau in the absorption profile indicates that the density of silicon is indeed sufficient to make the streams optically thick. The shape and value of the absorption profile correspond to the stream width of (7÷8)R p . Most of the absorption takes place outside of the Roche lobe, while inside it the absorption is only 5.3%. This effect is caused by the delayed appearance of Si III ions, produced only when the ionization degree of hydrogen becomes sufficiently high. For the same reason, the transit depth of Si III is also strongly dependent on XUV flux. The effect of the PW streams is even more pronounced here than in the case of carbon. In particular, with the zero PW velocity (the case of model M3) the Si III transit depth drops from 40% to 13.4%. Figure 8 also shows the absorption profile, calculated for the silicon abundance 20 times less than that in the solar system. For such a reduced abundance, the simulated transit depth equals the values measured by Linsky et al. (2010). Under the conditions of a slow SW with a total pressure p SW =5×10 −6 μbar (the case of M7), the simulated Si III absorption decreases, as compared to the very weak SW case (the model run M1), down to about 28%. This is because the width of the PW stream decreases to (5÷6)R p .

Calculation of the Magnesium Lines
The observational results regarding Mg I and Mg II resonant lines impose a certain challenge for interpretation. The only available measurements of the absorption of the Mg I line (Vidal-Madjar et al. 2013) were made at low resolution (∼43 km s −1 ). They show about 7.5% difference between the in-transit and out-of-transit Mg I absorption values. At the same time, the absence of any difference in the intensity of Mg II k&h lines in-and out-of-transit has been definitely measured, with an accuracy of better than 2%. However, if Mg I is present in the planetary exosphere in an abundance that yields to a moderate absorption in its line, then Mg II should be present as well in a sufficient quantity to produce a significant (i.e., measurable) absorption. This is because of the high efficiency of the photo-ionization of Mg I, which cannot be easily quenched by recombination. As estimated by Vidal-Madjar et al. (2013), the recombination quenching of the Mg II stage requires an unrealistically high electron density, which itself is proportional to the ionizing XUV flux. Unlike Mg I, Mg II ions should be present over the whole planetary exosphere, producing a large resonant absorption in the ±50 km s −1 interval of the emission line that is comparable in magnitude to the Si III line. The calculated absorption profiles of the Mg lines are presented in Figure 9. One should note that no radiation pressure force acting on Mg was included in our model.
The Mg I atoms mostly exist in the relatively dense upper atmosphere of the planet, close to the considered in the model conventional planetary surface at R p , so the corresponding absorption line is produced mainly due to the natural line broadening mechanism. As can be seen in Figure 6, this absorption is weakly dependent on Mg abundance. At the same time, the Mg II absorption line is produced to a significant degree by particles in Doppler resonance with emission profile. Their number is strongly affected by the magnesium abundance. In particular, at typical magnesium abundance in the solar system (Mg/H=3.7×10 −5 ) the resonant part of the Mg II line absorption comprises practically the whole absorption, whereas at a hundred times reduced abundance (Mg/H=3.7×10 −7 ) it comprises only a half of the total absorption, and at a thousand times reduced abundance the resonant part of absorption becomes negligible. However, at an abundance Mg/H=3.7×10 −8 , the absorption in the core of Mg II line in the ±50 km s −1 interval due to natural line broadening mechanism would still be observable, remaining at the level of about 3%.
There is another specific process that can result in a significant depletion of Mg I and Mg II. This is the resonant double-chargeexchange with a Helium ion: Mg+He + →Mg 2+ +He+e studied in DuBois (1986) at relatively high impact energies ∼40 keV, for which it had maximum rate with a cross-section of about 10 −16 cm 2 . The rate of this reaction at the temperatures of interest (>1000 K) is unfortunately unknown. Therefore, we checked its importance in the context of our model and performed a series of dedicated simulation runs, where the assumed cross-section of the double-charge-exchange reaction with a helium ion was changing in the range from 8×10 −15 cm 2 to 2×10 −24 cm 2 . As shown in Figure 10, if the cross-section at temperatures of ∼10 4 K exceeds 10 −18 cm 2 , the absorption by the magnesium lines becomes negligible, even at the abundances typical for the solar system.

Discussion and Conclusions
We computed and analyzed the absorption in spectral lines of heavy minor elements occurring during the primary transit of hot Jupiter HD 209458b. The calculations were based on the self-consistent 2D hydrodynamic multi-fluid modeling of the expanding planetary upper atmosphere, taking into account the tidal forces, under different SW and XUV radiation conditions, for various abundances of the considered elements.   A summary of all the results presented above on the simulation of absorption by various minor species at HD 209458b under different conditions and model parameter sets (referred as the models M1 KM8) is given in Table 2.
It has been shown for the first time, that transit depths observed for oxygen and carbon resonant multiplets can be reproduced, assuming the typical for solar system abundances' values: O/H=8.5×10 −4 , C/H=3.6×10 −4 . These transit depth scan be even larger than that for the Lyα line, reaching values of about 10%. The performed simulations confirm the previously proposed idea, that the bulk motion of the planetary escaping material with a velocity of several tens of km s −1 , dragging the heavy elements, can sufficiently increase the absorption due to Doppler resonance, which, in contrast to previous studies, includes also the bulk velocity of material in the streams and not only the thermal motion of absorbing particles. This mechanism will dominate in the absorption of the corresponding stellar lines.
The motion of hot Jupiters' upper atmospheric material, that extends far outside the Roche lobe, and its interaction with SW has been studied in our previous works (Shaikhislamov et al. 2016;Khodachenko et al. 2017). At close orbits, the tidal force and XUV heating are strong enough to drive the expanding planetary atmosphere beyond the Roche lobe and to accelerate its material in the form of a double stream PW at the day-and night-sides of the planet. The material velocity inside the streams reaches ±(25÷35) km s −1 at distances of (8÷10)R p along the planet-star line. This fits well with the measured width of many stellar FUV lines, including O I, C II, and Si III. We showed that without such bulk stream motion, the transit depths of the considered elements are below that of the Lyα line and appear approximately the same as those obtained in earlier works based on 1D hydrostatic models, which ignored the motion of escaping material, and especially the presence of the accelerated PW streams outside of the Roche lobe. The 2D axisymmetric model employed in this work allows one to study the flow of the planetary upper atmospheric material beyond the Roche lobe. However, its applicability is restricted within the region r<10 R p around the planet, or even smaller than that. At such distances the PW streams driven by the tidal force should bend clockwise due to the Coriolis force, resulting in a decrease of the acceleration due to gravity along the star-planet line. This is a consequence of momentum conservation, which does not allow material with an orbital velocity to fall directly toward a star, forcing instead its spiraling. Thus, the increase of PW velocity with the distance along the tidally pulled streams, obtained in the present simulation, is most likely an overestimation. A more precise treatment requires 3D modeling. The existing attempts of 3D modeling of hot Jupiters so far (Bisikalo et al. 2013;Bourrier & Lecavelier des Etangs 2013;Kislyakova et al. 2014;Matsakos et al. 2015) still do not include self-consistently the physics of planetary exosphere formation, production and dynamics of ENAs, as well as minor species. In spite of the limitations of the 2D modeling approach, our work clearly shows that the transit depth above 4% in the resonant lines of heavy species comes from the planetary upper atmospheric material streaming beyond the Roche lobe, being accelerated by the tidal force.
While the ions, such as C II, are efficiently dragged by the PW flow via Coulomb collisions with protons, the situation with atoms is more complicated. At typical densities, modeled at the exosphere beyond the Roche lobe, the elastic collisions are only marginally effective for the dragging of heavy neutral particles. However, in the case of O I, the exchange of momentum and temperature with hydrogen atoms and protons proceeds via resonant charge-exchange, that couples the interacting particles and has an order of magnitude larger cross-section than the cross-section of elastic collisions.
Our simulations have shown that the transit depth of the absorption profile of O I and C II lines is not strongly affected by changes in the intensity of the XUV flux or SW parameters. However, if the SW total pressure becomes high enough to stop the PW stream on the day-side and to redirect the material flow toward the tail, the absorption profiles become strongly asymmetric with respect to their red and blue parts. Therefore, spectrally resolved measurements can be used to reveal the presence of the SW plasma flows with pressures exceeding 10 −4 μbar. Note that, as shown in our previous study (Khodachenko et al. 2017), in such a case the Lyα line should also have a strong asymmetry and significant increase in the transit depth. Thus, the detection of significant asymmetry simultaneously in Lyα, O I, and C II lines during the same transit event will unequivocally indicate the presence of a strong SW, as no other mechanism, such as e.g., radiation pressure, can explain this effect.
The results of the present modeling qualitatively agree with the observations of transits assuming the abundances of oxygen and carbon close to the solar system values. However, the simulation of absorption by the resonant lines of silicon and magnesium ions hints that the abundances of these species in the upper atmosphere of HD 209458b should be much less than the typical solar system values. The corresponding lines of Si III and Mg II have an order of magnitude larger cross-section values for the resonant absorption than O I and C II lines. These ions are rapidly produced from atomic components via photoionization and, in case of Si III, via the resonant chargeexchange of Si II with protons. Because of that, the PW streams produce significant transit depths in the resonant lines of Si III and Mg II within the Doppler-shifted velocity range of ±50 km s −1 . The agreement of the simulations with the available measurements of the Si III line by Linsky et al. (2010) is achieved for the silicon abundance 20 times less than that in the solar system. Moreover, re-analysis of the same data of COS and STIS instruments by Ballester & Ben-Jaffel (2015) shows that the absorption by Si III was not detected at all, whereas the previously reported detections may be just the effects of a strong stellar variability. In this case the results of our model indicate that abundance of Mg and Si in the upper atmosphere should be very low indeed, possibly due to silicate condensation deeper down in the lower atmosphere. The absence of absorption in the Mg II lines appears to be an enigma, which is more difficult to explain. To have the Mg II transit depth less than 2%, as reported in Vidal-Madjar et al. (2013), magnesium should practically be absent in the upper atmosphere of HD 209458b, although all thermochemical simulations performed so far do not predict any reduced abundances for silicon or magnesium (Visscher et al. 2010). This paradox can be resolved if the resonant double-chargeexchange reaction of a magnesium atom with a Helium ion (Mg+He + →Mg 2+ +He+e) has a cross-section larger than 10 −18 cm 2 at temperatures of 10 4 K. In this case, the reduction of abundances of silicon and magnesium by an order of magnitude, as compared to the solar system values, would be