The escape of heavy atoms from the ionosphere of HD209458b. I. A photochemical-dynamical model of the thermosphere

The detections of atomic hydrogen, heavy atoms and ions surrounding the extrasolar giant planet (EGP) HD209458b constrain the composition, temperature and density profiles in its upper atmosphere. Thus the observations provide guidance for models that have so far predicted a range of possible conditions. We present the first hydrodynamic escape model for the upper atmosphere that includes all of the detected species in order to explain their presence at high altitudes, and to further constrain the temperature and velocity profiles. This model calculates the stellar heating rates based on recent estimates of photoelectron heating efficiencies, and includes the photochemistry of heavy atoms and ions in addition to hydrogen and helium. The composition at the lower boundary of the escape model is constrained by a full photochemical model of the lower atmosphere. We confirm that molecules dissociate near the 1 microbar level, and find that complex molecular chemistry does not need to be included above this level. We also confirm that diffusive separation of the detected species does not occur because the heavy atoms and ions collide frequently with the rapidly escaping H and H+. This means that the abundance of the heavy atoms and ions in the thermosphere simply depends on the elemental abundances and ionization rates. We show that, as expected, H and O remain mostly neutral up to at least 3 Rp, whereas both C and Si are mostly ionized at significantly lower altitudes. We also explore the temperature and velocity profiles, and find that the outflow speed and the temperature gradients depend strongly on the assumed heating efficiencies...


Introduction
speaking, the equations are valid below the exobase, and terms due to heat 180 conduction and viscosity gain significance as Kn → 1. We note that the 181 exobase on HD209458b is typically located at a very high altitude (see Sec-182 tion 3.1), and viscosity and heat conduction are not particularly important. 183 We included species such as H, H + , He, He + , C, C + , O, O + , N, N + , 184 Si, Si + , Si 2+ , and electrons in the hydrodynamic model. We also generated 185 simulations that included Mg, Mg + , Na, Na + , K, K + , S, and S + , but the 186 presence of these species did not affect the density profiles of H, O, C + , or Si 2+ 187 significantly. The model includes photoionization, thermal ionization, and 188 charge exchange between atoms and ions. The reaction rate coefficients for 189 these processes are listed in Table 1. Multiply charged ions were included only 190 if the ionization potential of their parent ion was sufficiently low compared to 191 the thermal energy and radiation field in the upper atmosphere. We note that 192 our model also includes impact ionization by thermal electrons. In general,193 this can be important for species with low ionization potential such as alkali 194 metals (e.g., Batygin and Stevenson, 2010), although we find photoionization 195 to be more significant in the thermosphere (see Section 3.2).

360
The results of the DSMC calculations can be incorporated into hydro- increases the escape velocity of the electrons.

383
In order to incorporate the ambipolar electric field in the modified Jeans 384 conditions we expressed the Jeans parameters for ions and electrons as: where φ e is the ambipolar electric potential, q i,e is the electric charge and where r 1 = 1.08 R p , δz 0 = 10 km, and f c = 1.014. We solved the equations of  temperature gradient above the peak. As a result, the temperature profile is 461 almost isothermal when η net = 1.

462
It is also interesting that the temperature profiles in the models that are For η net ranging from 0.1 to 1, the pressure averaged temperature below 3 R p 484 varies from 6,200 K to 7,800 K for the average solar flux. In the H50 model 485 the pressure averaged temperature is 7,000 K. We note that T p is a fairly    Table 2 for the input parameters). Model C1 assumes a constant where F E is the flux of energy, Φ pe is the photoelectron flux (cm −2 s −1 ), n e 572 is the density of thermal electrons (cm −3 ) and
Assuming that all of the energy is deposited by electrons that are thermalized 575 within a path element dr, we can estimate the e-folding length scale for 576 thermalization of photoelectrons with different energies as follows: We calculated Λ pe for different photoelectrons based on the C1 model, 578 and compared the result with the vertical length scale H of the atmosphere.

579
The latter is either the scale height or R p , depending on which is shorter.   Figure 6: Volume heating rate as a function of pressure in the C2 model due to the absorption of stellar XUV radiation between 1 and 1000Å in 200Å bins.
temperature at higher altitudes in the C2 model decreases much more rapidly 607 with altitude than in the H50 model. 608 Figure 7 shows the terms in the energy equation based on the C2 model.

609
In line with previous studies, stellar heating is mainly balanced by adiabatic 610 cooling. Advection cools the atmosphere at low altitudes below the temper-611 ature peak, whereas at higher altitudes it acts as a heating mechanism. In 612 fact, above 2 R p the adiabatic cooling rate is higher than the stellar heating 613 rate because thermal energy is transported to high altitudes by advection 614 from the temperature peak. The radiative cooling term that is centered near . Advection acts as a cooling mechanism below 1.5 R p and a heating mechanism above this level.
is not included in the output. Conduction is not significant at any altitude 618 in the model. We note that the rates displayed in Figure 7 balance to high 619 accuracy, thus implying that the simulation has reached steady state.

620
The differences between the C1 and C2 models are subtle. The peak tem- tion. This should not be confused with the assumption of a constant η net , which leads to a different temperature profile when compared with either C1 and C2 (see Figure 3). In general, the maximum and mean temperatures 630 in models C1-C4 are relatively similar. Thus we conclude that the mean 631 temperature in the thermosphere of HD209458b is approximately 7,000 K 632 and the maximum temperature is 10,000-12,000 K.

633
The substellar tide is included in the C3 model. We included it mainly with thermal electrons and other species, decays radiatively. We included this cooling mechanism in the C4 model by using the rate coefficient from Glover and Jappsen (2007) that includes a temperature-dependent correction 653 to the rate coefficient given by Black (1981). We also included an additional where ξ = r/r 0 , c = kT /m is the isothermal speed of sound, W = GM p /r 0 , 669 and m is the mean atomic weight. It is often assumed that the vertical The mean atomic weight can be less than 1 because electrons contribute to the number density but not significantly to the mass density. R p where c(ξ c ) = 4.1 km s −1 . This is because the temperature gradient 691 of the model is steeper than the corresponding gradient in the C1 model.

692
The volume averaged mean temperature below 15 R p in the C2 model is  Another problem is that the atmosphere is not isothermal. In fact, the 703 temperature gradient above the heating peak in models C1-C4 (Table 2) is 704 relatively steep, and in some cases it approaches the static adiabatic gradient 705 (T ∝ r −1 ) as defined by Chamberlain (1961). Assuming that the temperature 706 profile can be fitted with c 2 = c 2 0 /ξ β above the heating peak, the estimated 707 values of β for the C1 and C2 models are 0.7 and 0.9, respectively. We note  Here we evaluate the mass loss rates based on our models. We define the 732 mass loss rate simply as: We note that the solar spectrum that we used in this study contains the total flux of 4 × 10 −3 W m −2 at wavelengths shorter then 912Å (the ionization 735 limit of H) when normalized to 1 AU. This value is close to the average solar 736 flux of 3.9 × 10 −3 W m −2 at the same wavelengths (e.g., Ribas et al., 2005).  The gray assumption also fails to include the fact that the net heating effi-800 3 By chance the incident flux is equal to the mean solar flux divided by 4 that we used as a 'globally averaged' value. Here, however, it is taken to be the substellar value.

815
As an example, we calculated Kn 0 and X 0 (see Section 2.1.2) based on the 816 C1 and C2 models. The Knudsen number Kn depends on the mean collision 817 frequency, and it is much smaller than unity at all altitudes below 16 R p .
is approximately constant. This criteria is consistent with the equations of 823 motion, and it means that r 0 that should be used to calculate X 0 is above

854
We compared the temperature and velocity profiles from the C1 and C2 855 models with results from similar models C5 and C6 that use the modified 856 Jeans conditions. Note again that our version of the modified Jeans condi-857 tions includes the polarization electrostatic field that is required in strongly 858 ionized media. Figure 8 shows the temperature and velocity profiles from the 859 models. There is no difference between the C5 model and the C1 model as altitude. In summary, we have shown that the C1 and C2 models are both 868 consistent with kinetic theory. 869 We note that extending the models to 16 R p or higher is not necessarily   Based on the gas giants in the solar system it might be expected that heavy  Table 2 for the input parameters). prove that diffusive separation does not take place.

896
In order to illustrate the results, Figure 9 shows the density profiles of  Once again, the differences between the earlier models and our work arise 908 from different boundary conditions, and assumptions regarding heating rates 909 Figure 9: Density profiles in the upper atmosphere of HD209458b based on the C2 model (see Table 2 for the input parameters).  secondary in a sense that they require the ions to be produced by some other 952 mechanism. In fact, H + is mainly produced by photoionization (P1), al-953 though thermal ionization (R3) is also important near the temperature peak.

954
The production rates are mainly balanced by loss to radiative recombination . Compared with our models (see Figure 9), the H/H + transition occurs at a significantly lower altitude. The difference arises from the lower boundary conditions and gray approximation to heating and ionization in the MC09 model. peak. The production is balanced by loss to radiative recombination (R10).

967
The chemical loss timescale for C is shorter than the timescale for advection 968 below 1.8 R p . Thus advection is unable to move the C/C + transition to 969 altitudes higher than 1.2 R p . collisions between a singly ionized species and H + . Thus diffusion is even 1051 less significant for a species like Si 2+ . We verified these results from our sim-1052 ulations by switching diffusion off in the model and rerunning the C2 model.

1053
As a result the density of the heavy atoms and ions increased slightly at high 1054 altitudes, but the differences are not significant -the results were nearly 1055 identical to the density profiles of the original simulation. 1056 We note that the atmosphere can also be mixed by vertical motions asso-1057 ciated with circulation that are sometimes parameterized in one-dimensional 1058 models by the eddy diffusion coefficient K zz (e.g., Moses et al., 2011). This 1059 mechanism is efficient in bringing the heavy elements to the lower thermo- 1082 where q R is the volume heating rate, and ν ei and ν in are the electron-ion and ion-neutral momentum transfer collision frequencies, respectively.
the detection of heavy atoms and ions in the thermosphere is not without controversy, and the detection of Si 2+ is particularly intriguing. Thus these observations present several interesting challenges to modelers.

1107
The observed transit depths are large, and substantial abundances of the 1108 relevant species are required to explain the observations. However, on every 1109 planet in the solar system heavier species are removed from the thermo-

1136
For net heating efficiencies between 0.1 and 1, the mean temperature below 1137 3 R p varies from 6,000 K to 8,000 K. This means that 8,000 K is a relatively

1148
This would also lead to higher outflow velocity and mass loss rate. However, 1149 the uncertainty in the observed transit depths is also large (Ben-Jaffel, 2008; atures. Therefore our reference model C2 with a mean temperature of 7,200 tude range above the temperature peak. This conclusion is supported both by the hydrodynamic model and new constraints from kinetic theory (Volkov it should be treated with caution. We used parameterized heating efficiencies 1178 for low energy photons, and the location of the sonic point is very sensitive to 1179 the temperature profile. Also, the stellar tide can enhance the escape rates at 1180 the substellar and antistellar points. We did not include this effect because 1181 it may produce horizontal flows that cannot be modeled in 1D.

1182
As long as the upper boundary is at a sufficiently high altitude, we found of Si 2+ is also significant. We found that neutral heavy atoms are dragged that results such as these cannot be freely generalized to other extrasolar 1241 planets. As in the solar system, each planet should be studied separately.

1242
For instance, the dissociation of molecules depends on the temperature pro-file that is shaped by the composition through radiative cooling and stellar 1244 heating. The mass loss rate and escape velocity, that determine whether dif-1245 fusive separation takes place or not, depends on the escape mechanism that 1246 again depends on the temperature and composition of the upper atmosphere.

1247
The results from different models can only be verified by observations that