Hα ABSORPTION IN TRANSITING EXOPLANET ATMOSPHERES

Absorption¹ of starlight by a transiting planet provides a probe of the atmospheric composition and structure. The Lyα and Hα transitions of hydrogen may potentially give rise to large transit depths, as the combination of large cross section and hydrogen abundance put the optical depth unity surfaces at much higher altitudes (μbar–nbar pressure levels) than the continuum photosphere (mbar–bar).

The detection of Lyα absorption (Vidal-Madjar et al. 2003, 2004, 2008; Ben-Jaffel 2007, 2008; Ehrenreich et al. 2008, 2012; Lecavelier Des Etangs et al. 2010; Lecavelier des Etangs et al. 2012) due to atoms in the ground state traces the bulk of the atomic population, and hence the density profile of the planet's atmosphere. By contrast, the recent detection of Hα absorption by Jensen et al. (2012) reveals the small population in the n = 2 excited state. As the number density in the n = 2 state depends sensitively on the excitation rates due to collisions and radiative pumping by the stellar photons, it is a much more sensitive probe of the atmosphere, in particular the temperature profile. Hence, the Lyα and Hα lines are complementary probes of the atomic hydrogen layer.

The goal of this paper is to construct a detailed model of the n = 2 state population by which to understand the Jensen et al. (2012) Hα transit depth of HD 189733b. Section 2 reviews the relevant observations of Lyα and Hα observations for HD 209458b and HD 189733b. Section 3 presents the details of the atmospheric model. Section 4 presents numerical results for HD 189733b, and a detailed comparison to the Jensen et al. (2012) data. Section 5 presents the predictions of our model for HD 209458b, and they disagree significantly with the data. Section 6 addresses the possibility of observing Balmer continuum absorption. We summarize our results and compare to previous investigations in Section 7. In the Appendix, we study the effect of the stellar Lyα on the cooling rate.

We first review the observations of Hα absorption that motivate our work, as well as previous attempts to determine excited state abundance, or excitation temperature, from the observations. These previous estimates of density and temperature will be compared to the calculations in this paper in Section 7.

The initial detection of Lyα absorption during the transit of HD 209458b (Vidal-Madjar et al. 2003) indicated a 15% decrease in flux over a 5.1 Å wavelength band around line center. The implied occultation radius of 4.3 R_J, where R_J = 7.149 × 10⁹ cm is the radius of Jupiter, is larger than the Roche lobe radius of 3.6 R_J, leading the authors to suggest that the planet must be in a state of Roche lobe overflow. Due to absorption by the interstellar medium (ISM) and geocoronal emission, the center of the Lyα line cannot be used, and absorption is only reliably observed outside ∼75 km s⁻¹ from line center, indicating either a population of "hot hydrogen" (T ∼ 10⁶ K) formed through charge exchange with stellar wind protons (Holmström et al. 2008; Ekenbäck et al. 2010; Tremblin & Chiang 2013) or large columns of "warm hydrogen" (T  ∼  10⁴ K; Yelle 2004), which yield a saturated line. This paper attempts to explain Hα transit spectra by the latter scenario of a large column of warm hydrogen.

For T  ∼  10⁴ K, the observed Lyα absorption is well out on the damping wing of the line, and hence the cross section is not a strong function of temperature. The Lyα transit depth is, however, dependent on the temperature through the atmospheric scale height. As we will show, the Hα absorption is expected to have a stronger dependence on temperature due to the collisional excitation rate of the n = 1 to n = 2 transition.

The first attempt to detect Hα was by Winn et al. (2004) who observed HD 209458. For the remainder of this section, we summarize their observational results, assumptions, and inferences. Integrating over a 5.1 Å band, they find an equivalent width

$$\begin{equation} W_\lambda \equiv \int d\lambda \left(\frac{F_{\rm out}(\lambda)-F_{\rm in}(\lambda)}{F_{\rm out}(\lambda)} \right) < 1.7 {\rm \, m\mathring{\rm{A}}}, \end{equation} \tag{ 1 }$$

where F_in and F_out are the flux in and out of transit, respectively. If the n = 2 state hydrogen occults an area ΔA, and the area of the star is A, then in the optically thin limit, Equation (1) can be used to relate the column density of the n = 2 state, N₂, to the equivalent width as ²

$$\begin{equation} W_{\rm H_\alpha } = \frac{\pi e^2}{m_ec^2}f_{23}\lambda _{\rm H_\alpha }^2N_2 \frac{\Delta A}{A} \,. \end{equation} \tag{ 2 }$$

Here $\lambda _{\rm H_\alpha }=6562.8\, \mathring{\rm{A}}$ is the wavelength of the Hα transition, and f₂₃ = 0.64 is the oscillator strength. Using a fractional area ΔA/A = 0.15 from Lyα measurements (Vidal-Madjar et al. 2003), and assuming the same occulting area for Lyα and Hα lines, the column density can be constrained to be N₂ ≲ 2.4 × 10¹¹ cm⁻².

Hα absorption was subsequently detected by Jensen et al. (2012) in the transits of HD 209458b and HD 189733b. The line profile of HD 209458b showed an excess to the blue and a deficit to the red of line center, for which we have no simple explanation. The line profile of HD 189733b offered no such consternation, exhibiting a symmetric absorption feature about line center. HD 189733b is the focus of our modeling effort.

To quantify the absorption, Jensen et al. (2012) introduce the absorption measure M_abs, defined as

$$\begin{equation} M_{\rm abs} = \left< S_T\right>_{\rm central}-\frac{\left<S_T\right>_{\rm blue} + \left<S_T\right>_{\rm red}}{2}, \end{equation} \tag{ 3 }$$

where

$$\begin{equation} \left< S_T \right>_{\it i} = \left(\frac{F_{\rm in}}{F_{\rm out}}\right)_{\it i} - 1 \end{equation} \tag{ 4 }$$

and i = blue, red, central indicates the domain of the wavelength integration. The purpose of subtracting off the bands just outside the line center is to derive the drop in flux due to the upper atmosphere, where the Hα absorption takes place; absorption of the neighboring continuum takes place much deeper in the planet's atmosphere, near the continuum photospheres at mbar–bar pressures. In their analysis, Jensen et al. (2012) considered three 16 Å bands with the central band covering the Hα line and blue and red being adjacent bands at shorter and longer wavelengths, respectively. The absorption measure M_abs can be related back to equivalent width by W_λ ≈ M_absΔλ. For HD 189733b, the absorption measure is M_abs = (− 8.72 ± 1.48) × 10⁻⁴ (Jensen et al. 2012).

Winn et al. (2004) and Jensen et al. (2012) estimate the excitation temperature, T_exc, starting with the Boltzmann distribution for the number densities n₁ and n₂ of atoms in states n = 1 and n = 2:

$$\begin{eqnarray} \frac{n_2}{n_1} & =& \frac{g_2}{g_1} e^{-10.2\ {\rm eV}/k_B T_{\rm exc}}, \end{eqnarray} \tag{ 5 }$$

where g₂ = 8 and g₁ = 2 are the degeneracies of each state, including both 2s and 2p upper states, and the energy difference is 10.2 eV. The authors then estimate the left hand side by equating n₂/n₁ ≃ N₂/N₁, where N₁ is the column density of ground state hydrogen. This assumption relies on the distribution of n₁ and n₂ being similar in the atmosphere, so that the effective path length in the integration is comparable. Lastly, they assume that N₂/N₁ ≃ W_Hα/W_Lyα = 0.0128 Å/0.32 Å = 0.04, which is valid if both transitions are optically thin and the occulting area ΔA is the same for each. Plugging this into the left hand side and solving then gives T_exc = 2.6 × 10⁴ K. This high excitation temperature is not achieved for thermal gas in any published models to date (e.g., Yelle 2004; Murray-Clay et al. 2009; García Muñoz 2007). We will revisit the estimates of T_exc in Section 7 in the context of our model for the atmospheric structure and level populations.

In this section we present our model for the density and temperature profiles and the hydrogen level populations. Section 3.1 outlines the different layers of the atmosphere of interest. Section 3.2 discusses thermal and ionization equilibrium. Section 3.3 presents the rate-equilibrium equations which determine the ground and excited state populations of hydrogen. Section 3.4 discusses the numerical solution of the coupled equations of hydrostatic balance, thermal equilibrium photoionization equilibrium and detailed balance for the level populations.

3.1. Basic Structure of the Atmosphere

We consider a spherically symmetric atmosphere above the photospheric radius at r = R_p, where r is the spherical radius. The atmosphere is divided into three zones (see Figure 1), according to the state of hydrogen: molecular, atomic and ionized. Inside the photospheric radius, r < R_p, we consider the planet to be opaque to Hα. The shell R_p < r < R_b is composed mainly of molecular gas, due to the low temperature. Assuming a Boltzmann distribution and a mean temperature of T_eq = 1200 K, the equilibrium temperature (for zero albedo) of HD 189733b, the relative abundance of n = 2 hydrogen is n₂/n₁ = 5.65 × 10⁻⁴³, where n₁ and n₂ are the number density of hydrogen in the n = 1 and n = 2 states, respectively (we use n for both number density and radial quantum number. The meaning is clear from the context). Due to the small abundances and low temperatures, we assume that the contribution of the molecular layer to the Hα absorption is negligible. The transition from molecular to atomic hydrogen takes place at a radius R_b, and it is this "atomic" layer at r > R_b and P ≲ 1 μbar that gives rise to the Hα absorption in our model. At pressures P ≲ 1 nbar hydrogen will be mostly photoionized, forming an "ionized" layer that contributes little to the Hα absorption.

Zoom In Zoom Out Reset image size

A parameter that enters our model is the base radius of the atomic layer R_b. It is determined by the photospheric radius R_p and the thickness of the molecular layer R_b − R_p. If the molecular layer is isothermal at temperature T, the thickness is

$$\begin{equation} R_{\rm b}-R_{\rm p} \simeq \frac{k_{\rm B}T R_{\rm p}^2}{\mu m_{\rm H}G M_{\rm p} }\ln \left(\frac{P_{\rm p}}{P_{\rm b}}\right), \end{equation} \tag{ 6 }$$

where μ ≃ 2.3 is the mean molecular weight, M_p is the mass of the planet, P_p is the pressure of the optical photosphere, and P_b is the pressure at the molecular-to-atomic transition. For T = T_eq, P_b = 1 μbar and P_p = 1 bar, the thickness is R_b − R_p = 0.032 R_p. An upper bound to the thickness is found using the dissociation temperature at which the gas transitions from molecular to atomic. Plugging P_b = 1 μbar into the Saha equation, and using the roto-vibrational energies from Borysow et al. (1989), we find T = 1934 K and R_b − R_p = 0.054 R_p. Since the equilibrium and dissociation temperature estimates differ by only 40%, and both give altitudes which are a small fraction of R_p, we conclude the thickness of the molecular layer will not greatly affect the Hα transit depth. Henceforth we use R_b − R_p = 0.054 R_p in numerical calculations for HD 189733b.

3.2. Ionization and Thermal Equilibrium

The electron and ion abundances are set by rate equilibrium between photoionization of ground-state hydrogen and radiative recombination,

$$\begin{equation} \Gamma _{\rm 1s}\left(N_{\rm 1s}\right)n_{\rm 1s} = \alpha _{\rm B}(T) n_{\rm e}^2\,\,. \end{equation} \tag{ 7 }$$

Here Γ_1s is the photoionization rate from the 1s state, defined in Equation (9) below, and α_B is the Case B recombination coefficient defined in Table 2. As discussed below, we assume all free electrons come from hydrogen so that n_e = n_p.

The heating rate is dominated by photoelectrons from the ionization of hydrogen. In the high-UV case we consider here, cooling will be dominated by Lyα emission (Murray-Clay et al. 2009). The temperature is then determined by balancing photoelectric heating (Q_1s) and Lyα cooling (Λ_Lyα),

$$\begin{equation} Q_{\rm 1s}\left(N_{\rm 1s}\right)n_{\rm 1s} = \Lambda _{\rm Ly\alpha }\left(T\right) n_{\rm e}n_{\rm 1s}, \end{equation} \tag{ 8 }$$

where the Lyα cooling rate is Λ_Lyα(T) ≃ 10.2 eV(C_{1s → 2s}(T) + C_{1s → 2p}(T)) and C_{1s → 2s} and C_{1s → 2p} are the collisional excitation rates from the 1s state to the 2s and 2p states, respectively (see Table 2). Excitation to the 2s state also contributes to Lyα cooling since the ℓ-mixing reactions rapidly turn 2s into 2p, which then emits a Lyα photon (see Figures 12 and 13). In Equation (8), Q_1s is the heating rate per particle, defined in Equation (10) below.

In cases where there is a non-zero Lyα excitation rate J_Lyα due to the stellar radiation field, there exists the possibility, for sufficiently low gas temperatures, of Lyα heating when the excitation temperature of the radiation field approaches or exceeds the gas temperature. In the Appendix, we estimate the contribution of true Lyα absorption to the overall heating rate as a function of J_Lyα and T. We find that for the solar value of J_Lyα, the cooling rate is not significantly altered by the inclusion Lyα heating for the bulk of the layer that is responsible for the Hα absorption; however, the heating could be important at lower temperatures near the base of the atomic layer (see Figure 16). To address this issue fully it is necessary to have a detailed treatment of Lyα radiation transfer which is beyond the scope of this paper.

While some parameters of the planet are well known (see Table 1), such as planetary mass (M_p), radius (R_p) and semi-major axis (a), other parameters are less certain, and may even vary with time, such as the ionizing and Lyα flux of the star HD 189733. Also uncertain is the role of day–night transport of heat in the upper atmosphere (Koskinen et al. 2007). Given these uncertainties, we parameterize the ionizing flux of HD 189733, and hence the heating and ionization rates, in an attempt to explain the observed transit depth and line profile. Specifically, the ionization rates will be multiplied by a factor ξ/4. The flux averaged over 4π steradians corresponds to ξ = 1. The substellar flux corresponds to ξ = 4. Ionization rates larger than the nominal substellar value have ξ ⩾ 4. For the Lyα rate, we use the solar spectrum and evaluate the flux at the orbital separation of HD 189733b. We take the value of J_Lyα to be constant throughout the atmosphere.

Table 1. HD 189733b Parameters^a

Mass (MJ)	1.138
Radius (RJ)	1.138
a (AU)	0.031

Note. ^aSource: exoplanet.eu

Download table as:  ASCIITypeset image

The ground-state photoionization and photoelectric heating rates are defined as

$$\begin{equation} \Gamma _{\rm 1s}\left(N_{\rm 1s}\right) = \frac{\xi }{4}\int _{\nu _0}^{\infty } \frac{4\pi J_{\nu }}{h\nu }\sigma _{\rm 1s} e^{-\sigma _{\rm 1s}N_{\rm 1s}}d\nu \end{equation} \tag{ 9 }$$

and

$$\begin{equation} Q_{\rm 1s}\left(N_{\rm 1s}\right) = \frac{\xi }{4}\int _{\nu _0}^{\infty }\frac{4\pi J_{\nu }}{h\nu } h\left(\nu - \nu _0 \right) \sigma _{\rm 1s}e^{-\sigma _{\rm 1s}N_{\rm 1s}} d\nu ,\,\, \end{equation} \tag{ 10 }$$

where J_ν is the mean stellar intensity, $N_{\rm 1s}(r) = \int _r^\infty dr^{\prime } n_{\rm 1s}(r^{\prime })$ is the 1s state column density, ν₀ = 13.6 eV/h is the ionization threshold frequency, and σ_1s is the photoionization cross sections for the 1s state (Osterbrock & Ferland 2006). The ultraviolet spectrum is a synthetic spectra for HD 189733b downloaded from the X-exoplanets Archive at the CAB (Sanz-Forcada et al. 2011). The integration over frequency allows higher energy photons to dominate the integral at larger depth (Trammell et al. 2011), leading to a larger rate than if a single frequency was used (Murray-Clay et al. 2009). Figures 2 and 3 show the photoionization and heating rates due to the 1s state as a function of N_1s, evaluated for HD 189733b. At small N_1s, the optically thin limit, the rates are constant, while at large columns they decrease roughly as a power law.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although ionization of the n = 2 state hydrogen contributes negligibly to the electron density, ionization has an important effect on the number densities n_2s and n_2p of the 2s and 2p states. The rates are given by

$$\begin{equation} \Gamma _{\rm 2s,2p} = \int _{\nu _1}^{\nu _0}\frac{4\pi J_\nu }{h\nu } \sigma _{\rm 2s,2p}\left(\nu \right) d\nu , \end{equation} \tag{ 11 }$$

where ν₁ = 10.2 eV/h is the frequency threshold for n = 2 state ionization. The upper frequency bound of ν₀ is included since photoionization of the ground state will both quickly attenuate ionizing radiation above this frequency and will contribute negligibly to the n = 2 photoionization rate due to the small cross-section. The spectrum used is taken from the Castelli and Kurucz Atlas (Castelli & Kurucz 2003). We omit the factor of one quarter found in Equation (9) due to the fact that Balmer continuum photons are optically thin throughout the region of interest, making the photoionization rate insensitive to the specific geometry considered. Additionally, we have ignored attenuation of the radiation field due to the small column densities involved.

3.3. Level Populations

In the upper atmosphere, several different physical processes are important in setting the level populations, including collisional excitation and de-excitation, and bound–bound and bound-free transitions due to the stellar radiation field. The 2s and 2p states must be considered separately due to the lack of a fast radiative transition between 2s and the ground state.³ The equation expressing rate equilibrium for the 2p state is

$$\begin{eqnarray} && B_{\rm 1s \rightarrow 2p} J_{\rm Ly\alpha } n_{\rm 1s} + C_{\rm 1s\rightarrow 2p}n_{\rm 1s}n_{\rm e} + C_{\rm 2s \rightarrow 2p} n_{\rm 2s} n_{\rm e} + \alpha _{\rm 2p}n_{\rm e}^2 \nonumber \\ &&\quad = A_{\rm 2p\rightarrow 1s}n_{\rm 2p} + B_{\rm 2p \rightarrow 1s} J_{\rm Ly\alpha } n_{\rm 2p} \nonumber \\ & &\qquad + C_{\rm 2p\rightarrow 1s}{n_{\rm 2p}} n_{\rm e} + C_{\rm 2p \rightarrow 2s} n_{\rm 2p} n_{\rm e} + \Gamma _{\rm 2p} n_{\rm 2p} \end{eqnarray} \tag{ 12 }$$

while that for the 2s state is

$$\begin{eqnarray} && C_{\rm 1s\rightarrow 2s}n_{\rm 1s}n_{\rm e} + C_{\rm 2p\rightarrow 2s}n_{\rm 2p}n_{\rm e} + \alpha _{\rm 2s}n_{\rm e}^2 \nonumber \\ &&\quad = C_{\rm 2s\rightarrow 1s}n_{\rm e}{n_{\rm 2s}} + C_{\rm 2s\rightarrow 2p}n_{\rm e}{n_{\rm 2s}} + \Gamma _{\rm 2s} n_{\rm 2s} + A_{\rm 2s\rightarrow 1s}n_{\rm 2s}. \qquad \end{eqnarray} \tag{ 13 }$$

Here C_{i → j} denotes the collisional transition from state i to state j, and both the ${\rm }1{\rm s} \leftrightharpoons 2{\rm s}$, ${\rm }1{\rm s} \leftrightharpoons 2{\rm p}$ and the ℓ-changing reaction ${\rm }2{\rm s} \leftrightharpoons 2{\rm p}$ are included. Collisional transitions induced by atomic and molecular hydrogen, and helium, are negligible. The terms $\alpha _{\rm 2p} n_{\rm e}^2$ and $\alpha _{\rm 2s} n_{\rm e}^2$ represent the effective recombination coefficients to that level, and include radiative cascade from higher levels (Draine 2011). We ignore ionization of helium, and assume the electrons contributed by other elements, such as sodium and potassium, are negligible in comparison to that from hydrogen, hence n_e = n_p. The A's and B's are the Einstein coefficients for radiative transitions, and J_Lyα is 1s → 2p transition rate associated with Lyα radiation. The rate coefficients are listed in Table 2.

Table 2. Relevant Reactions

Number	Reaction	Symbol	Rate	Reference
R1	H1s + γ → e− + H+	Γ1s	See Equation (9)
R2	e− + p → H + γ	αB	$2.54\times 10^{-13} (T/10^4 {\rm \, K})^{-0.8164-0.0208\log (T/10^4\,{\rm K})} \,{\rm cm^3\, s^{-1}}$	Draine (2011)
R3	H1s + e− → H2s + e−	C1s → 2s	1.21 × 10−8(104 K/T)0.455e−118400/T cm3 s−1	Janev et al. (2003)
R4	H1s + e− → H2p + e−	C1s → 2p	1.71 × 10−8(104 K/T)0.077e−118400/T cm3 s−1	Janev et al. (2003)
R5	H2s + e− → H2p + e−	C2s → 2p	$6.21\times 10^{-5}\left(\log \left(T/T_0\right)-\gamma \right)/\sqrt{T}\,{\rm cm^3\, s^{-1}}$	Janev et al. (2003)
R6	H2s + γ → e− + H+	Γ2s	See Equation (11)
R7	H2p + γ → e− + H+	Γ2p	See Equation (11)
R8	e− + H+ → H2s + γ	α2s	(0.282 + 0.047(T/104 K) − 0.006(T/104 K)2)αB	Draine (2011)
R9	e− + H+ → H2p + γ	α2p	αB − α2s	Draine (2011)
R10	H2s → H1s + 2γ	A2s → 1s	8.26  s−1	Osterbrock & Ferland (2006)
R11	H2p → H1s + γ	A2p → 1s	6.3 × 108  s−1	Osterbrock & Ferland (2006)
R12	Lyα cooling	ΛLyα(T)	10.2 eV(C1s → 2s + C1s → 2p)

Download table as:  ASCIITypeset image

Due to the negligible contribution from collisions, the 2p occupation is set by radiative excitation and de-excitation,

$$\begin{eqnarray} \frac{n_{\rm 2p}}{n_{\rm 1s}} & \approx & \frac{B_{\rm 1s \rightarrow 2p}J_{\rm Ly\alpha }}{A_{\rm 2p \rightarrow 1s}} \nonumber \\ & \simeq & 10^{-9}\ \left(\frac{5\ R_\star }{a} \right)^2 e^{16.9-\left(10.2\ {\rm eV}/k_b T_{\rm Ly\alpha , \star }\right) }. \end{eqnarray} \tag{ 14 }$$

For the analytic estimate of J_Lyα, we have included a dilution factor (R_⋆/2a)²  ∼  10⁻², for stellar radius R_⋆ and orbital separation a, and the line intensity at the stellar surface is ≃ (2hν³/c²)exp (− 10.2 eV/k_bT_{Lyα, ⋆}), where T_{Lyα, ⋆} ≃ 7000 K is the approximate excitation temperature for the solar Lyα. The 2p density is small mainly due to the small stellar excitation temperature, and additionally because of the dilution factor. We show the profile of n_2p in Figure 4 including all the physical effects in Equation (12), verifying that radiative rates dominate, and that the 2p population is far smaller than 2s.

Zoom In Zoom Out Reset image size

At sufficiently large n_e, the 2s state is primarily populated by collisional excitation from the ground state while de-population is primarily due to the ℓ-mixing transition 2s→2p. The 2s abundance can be estimated as

$$\begin{eqnarray} \frac{n_{\rm 2s}}{n_{\rm 1s}} & = & \frac{C_{\rm 1s\rightarrow 2s}}{C_{\rm 2s\rightarrow 2p}} = \frac{C_{\rm 2s\rightarrow 1s}}{C_{\rm 2s\rightarrow 2p}}\exp \left(-\frac{118400\,{\rm K}}{T}\right)\qquad\ \ \ \ \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} & = & 1.627\times 10^{-8}\left(\frac{T}{10^4\,{\rm K}}\right)^{0.045}e^{11.84-118400\, {\rm K}/T} \nonumber \\ & & \times \frac{8.633}{\log \left(T/T_0\right)-\gamma },\,\, \end{eqnarray} \tag{ 16 }$$

where T₀ = 1.02 K and γ = 0.57721... is the Euler–Mascheroni constant.

The strong temperature dependence in the exponential can be eliminated by using Equations (7) and (8),

$$\begin{eqnarray} n_{\rm 2s} & = & \frac{C_{\rm 2s\rightarrow 1s}Q_{\rm 1s} n_{\rm 1s}}{C_{\rm 2s\rightarrow 2p}\Lambda _{\rm Ly\alpha , 0}\left(T\right)n_{\rm e}}\qquad\ \ \ \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} & \approx & \frac{C_{\rm 2s\rightarrow 1s}\alpha _{\rm 1s}^{1/2}(T) n_{\rm 1s}^{1/2}}{C_{\rm 2s\rightarrow 2p}\Lambda _{\rm Ly\alpha , 0}\left(T\right)}\frac{Q_{\rm 1s}}{\Gamma _{\rm 1s}^{1/2}} \,\,, \end{eqnarray} \tag{ 18 }$$

where we have defined Λ_Lyα(T) ≡ Λ_{Lyα, 0}(T)exp (− 118400 K/T). The heating rate scales with N_1s as $Q(N_{\rm 1s})\approx Q_0/(0.22\sigma _{\rm pi}N_{\rm 1s})^{k_1}$ and the photoionization rate scales as $\Gamma _{\rm 1s} \approx \Gamma _{\rm 1s,0}/(\sigma _{\rm pi}N_{\rm 1s})^{k_2}$ (see Figures 2 and 3; Trammell et al. 2011). Along a radial path the column density can be approximated by N_1s ≈ H(T)n_1s, the 2s volume density scales as $n_{\rm 2s} \propto n_{\rm 1s}^{(1+k_2)/2-k_1`}$. The parameters k₁ and k₂ depend on the EUV spectrum, k₁ = 1.25 and k₂ = 1.09 for the synthetic HD 189733 spectrum giving $n_{\rm 2s} \propto n_{\rm 1s}^{0.035}$.⁴ As a result, we find that n_2s is fairly constant over a large range of pressures in the atomic layer.

3.4. Numerical Method

We consider a one-dimensional hydrostatic profile, with uniform spherical irradiation⁵ from the outside. We integrate the equation of hydrostatic balance, and the definition of the column,

$$\begin{equation} \frac{\partial P}{\partial r} = - \frac{G M_p \rho }{r^2} \end{equation} \tag{ 19 }$$

$$\begin{equation} \frac{\partial N_{\rm 1s}}{\partial r} = -n_{\rm 1s}, \quad\ \end{equation} \tag{ 20 }$$

where the pressure is given by P = ((1 + f_He)n_1s + (2 + f_He)n_e)k_BT and the density is ρ = m_Hn_H + m_en_e + m_pn_p + f_Hem_He(n_H + n_p). f_He = 0.1 is the fraction, by number, of helium per hydrogen nucleus, and although it is considered non-reactive in our model, it affects the profile through its contribution to the density and pressure. The relative helium abundance is considered to be constant throughout the atmosphere and is taken to have the solar value (Osterbrock & Ferland 2006).

To generate the atmospheric profile, the pressure is integrated from the outer boundary using Equation (19) through a fourth-order Runge–Kutta scheme. Given the pressure, temperature and ionization state can be solved for using Equations (7) and (8). The column density is then updated using Equation (20) and the process repeated until the base radius R_p is reached. To achieve the desired pressure P_b = 1 μbar at the base, the pressure at the outer boundary is varied using the secant method and new atmospheric profiles are generated until the target base pressure is reached.

Once the final profile is determined, the profiles for n_2s and n_2p can be calculated for the entire atmosphere using Equations (12) and (13).

With the radial temperature and the n_2s profile now calculated, the 2s state column density in the slant geometry, N_2s, and the optical depth, τ_ν, can be calculated as a function of impact parameter b,

$$\begin{equation} N_{\rm 2}\left(b\right) = \int _{-\infty }^{^\infty } n_{\rm 2s}\left(r\right)d\ell \qquad\quad\ \end{equation} \tag{ 21 }$$

$$\begin{equation} \quad \tau _\nu \left(b\right) = \int _{-\infty }^{\infty } n_{\rm 2s}\left(r\right)\sigma _{{\rm H\alpha },\nu }\left(r\right)d\ell \,\,, \end{equation} \tag{ 22 }$$

where ℓ is the line-of-sight coordinate and is related to the radial coordinate and impact parameter b by $r=\sqrt{b^2+\ell ^2}$. We do not include n_2p in the integrals due to its negligible contribution to the overall n = 2 abundance for HD 189733b. The proportional reduction in flux, ignoring the opaque disk at r < R_p, is then given by

$$\begin{equation} \frac{F^{\rm (out)}_{\nu }-F^{\rm (in)}_\nu }{F^{\rm (out)}_\nu } = \frac{2}{R_\star ^2}\int _{R_{\rm p}}^{R_\star } (1-e^{-\tau _\nu (b)})b{\rm d}b. \end{equation} \tag{ 23 }$$

Integrating Equation (23) with respect to wavelength recovers the equivalent width in Equation (1). We do not use this definition in practice, instead we choose to work with the absorption measure M_abs,

$$\begin{eqnarray} M_{\rm abs} & = & \frac{\int e^{-\tau _{\rm ISM}} \left(F^{\rm (out)}_\lambda -F^{\rm (in)}_\lambda \right)d\lambda }{\int e^{-\tau _{\rm ISM}}F^{\rm (out)}_\lambda d\lambda }\qquad\ \ \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} &&\qquad = \frac{\int d\lambda \int _{R_{\rm p}}^{R_\star } e^{-\tau _{\rm ISM}} (1-e^{-\tau _\lambda (b)})I^\star _\lambda b\,db}{\int d\lambda \int _{R_{\rm p}}^{R_\star }e^{-\tau _{\rm ISM}}I^\star _\lambda bdb}, \end{eqnarray} \tag{ 25 }$$

where τ_ISM is the optical depth of the ISM and $I^\star _\lambda$ is the unabsorbed stellar intensity. For Lyα, we use a Voigt profile with a temperature of 8000 K and assume an interstellar hydrogen column density N_H = 10^18.3 cm⁻² (Wood et al. 2005; Lecavelier Des Etangs et al. 2010). We assume no attenuation of Hα due to the ISM. For a constant $I^\star _\lambda$, Equation (25) is equivalent to M_abs = W_λΔλ. In calculating M_abs for Lyα, we assume that $I^\star _\lambda$ is proportional to

$$\begin{equation} I^\star _\lambda \propto \exp \left(-\frac{1}{2}\left(\frac{\Delta v}{64\,{\rm km\, s^{-1}}}\right)^2\right)\frac{c}{\lambda ^2} \end{equation} \tag{ 26 }$$

from Lecavelier Des Etangs et al. (2010). For Hα transits we assume that $I^\star _\lambda$ is constant with λ.

To study the effect of varied heating and ionization rates, we vary ξ between 1 and 20. For each model atmosphere, we calculate M_abs for both Hα and Lyα. For each ξ, we consider heating and ionization rates as defined by Equations (9) and (10), base radius R_b = R_p + 0.054R_p, (see Equation (6)) and base pressure P_b = 1 μbar.

Figure 4 shows the number density profile for the ξ = 1 case. Throughout both the atomic and ionized layers n_2s is seen to roughly trace n_e while n_2p follows the neutral hydrogen abundance, as expected from Equation (14). The ξ = 14 profiles, shown in Figure 5, exhibit the same behavior. Within the atomic layer, n_2s maintains a roughly constant abundance despite n_1s varying by three orders of magnitude, although it is not explicitly constant, as can be seen in Figures 10 and 11. The 50% ionization point is seen to move inward in pressure approximately linearly in with the increase in ionization rate (see Figures 4 and 5).

Zoom In Zoom Out Reset image size

The simulated transits for the abundance profiles in Figures 4 and 5 are shown in Figure 6, compared to the data for HD 189733b from Jensen et al. (2012). The model profiles have been computed using Equation (23). It is clear that ξ = 1 underestimates the transit depth, while the higher ionization rate ξ = 14 curve has roughly the correct width and depth. We note that there are significant oscillatory features in the measured data, which are much larger than the error bars shown in our Figure 6 and their Figure 3. These features persist in the wavelength regions outside the line. Jensen et al. (2012) attribute these features to systematic errors not included in the error bars shown.

Zoom In Zoom Out Reset image size

Next we vary ξ in order to derive the best fit to the corrected width M_abs = 8.8 × 10⁻⁴ reported by Jensen et al. (2012), with the corresponding values of M_abs shown in Figure 7. Additionally, we calculate M_abs for the unresolved Lyα. We find the Hα data are best fit by ξ = 14 (see Table 3). Lecavelier Des Etangs et al. (2010) determined the unresolved Lyα transit depth to be 5.05% ± 0.75%, including a 2.4% contribution from the photospheric disk. For the ξ = 14 model, we find M_abs = 5.94%, including the opaque disk contribution, within two standard deviations of the observed value.

Zoom In Zoom Out Reset image size

It is of interest to compare the location of the regions which contribute to M_abs for the Hα and Lyα. Figure 8 shows the contribution to the integral in Equation (25) as a function of b for both Hα and Lyα. The primary contribution for Lyα occurs farther out than that of Hα since Lyα becomes optically thick throughout the integrated band at lower column densities than Hα. The linear decrease for decreasing b is a geometric effect arising from the fact that annuli at larger impact parameters contribute more. For Hα, the largest contribution occurs at the base of the atomic layer with a noticeable decrease for impact parameters inside that layer. For large b, the contributions for both Hα and Lyα are small. The constancy at large b is due to the large scale height found at these radii. If a planetary wind is present, the densities at these radii could be much lower than found in our hydrostatic model (e.g., Yelle 2004), further decreasing their contribution to M_abs.

Zoom In Zoom Out Reset image size

Although ξ = 14 constitutes an order-of-magnitude increase in the heating rate, the change in temperature in the atomic layer is relatively modest. If we take the volume-weighted average temperature,

$$\begin{equation} \left<T\right> = \frac{3}{R_{\rm top}^3-R_{\rm b}^3}\int _{R_{\rm b}}^{R_{\rm top}} T\left(r\right) r^2 dr, \end{equation} \tag{ 27 }$$

where R_top is the radius where the electron density equals the neutral hydrogen density (see Figure 1), n_e = n_H. For the ξ = 1 case, we find <T > =8637 K, and for ξ = 14 we find <T > =8902 K. Previous investigators have found the need for additional UV heating (Lecavelier Des Etangs et al. 2010). Although we enforce this increase in temperature through an increase in the photoelectric heating rate for hydrogen, a more complete modeling of the heating and cooling processes in the atmosphere could explain the required temperature difference (Koskinen et al. 2012a, 2012b).

Figure 9 shows the temperature profiles as a function of gas pressure for both ξ = 1 and ξ = 14. Within the atomic layer (n_1s > n_e), the ξ = 14 case is approximately 1000 K above the temperatures in the ξ = 1 case at comparable pressures; however, due to the proportionately smaller ionization rate, the atomic layer extends to lower pressures for ξ = 1 which accounts for the similar values of <T >.

Zoom In Zoom Out Reset image size

We note that the temperatures at R_b are 6000–7000 K, much higher than the temperatures required for the formation of molecules. This is a limitation of our model due to the simplified prescriptions for heating and cooling.

The abundance profiles are shown for ξ = 1 and ξ = 14 in Figures 4 and 5, respectively. In both cases, n_2p traces the abundance of neutral hydrogen, as expected from Equation (14). Within the atomic layer, n_2s maintains a roughly constant abundance despite n_1s varying by three orders of magnitude. To better exhibit the small changes in n_1s, Figures 10 and 11 show a linear scale.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 shows the reaction rates for all reactions involved in the formation of 2s hydrogen for the ξ = 14 model. Within the hydrogen layer, the abundance is set by the balance of collisional excitation from 1s and the ℓ-mixing reaction. As pressures decrease, and the atmosphere becomes ionized, the contribution of radiative recombination increases. Once the atmosphere becomes predominantly ionized, radiative recombination becomes the dominant formation pathway.

Zoom In Zoom Out Reset image size

The reaction rates for 2p hydrogen are shown in Figure 13. Throughout the atmosphere, n_2p is set by the radiative transition between 1s and 2p, with little contribution from other rates. Attenuation of the Lyα radiation could play a role in decreasing n_2p; however, at μbar pressures we find that the radiative transition rates exceed the collisional rates by a factor of 10⁵.

Zoom In Zoom Out Reset image size

A similar analysis has been performed for HD 209458b. Since the line profiles for model and data are discrepant, we only produce a parameter study of M_abs for the same range of ξ used above.

Using parameters taken from exoplanet.eu, we take M_p = 0.714 M_J, R_p = 1.38 R_J, a = 0.04747 au, and use the simulated UV spectrum for HD 209458 from the X-exoplanets Archive at the CAB (Sanz-Forcada et al. 2011).

Figure 14 shows the absorption profile for ξ = 1. It is immediately obvious that there is no match between the model prediction and the observers profile. The model profile is symmetric about line center while and the observed profile is antisymmetric. The observed profile has a width characteristic of the orbital velocity (150 km s⁻¹) while the model is only a couple Doppler widths (10 km s⁻¹) wide. Figure 15 shows the dependence of M_abs on ξ. The curve is qualitatively similar to that of the HD 189733b case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It should be noted that for HD 209458b we find that 2p becomes more important than for HD 189733b. Accurately accounting for 2p requires modeling the radiative transfer of the Lyα which is beyond the scope of this paper.

Given the observation of Balmer continuum absorption by HD 209458b (Ballester et al. 2007), we estimate the transit depth for HD 189733b with our model. At the Balmer edge, the bound-free absorption cross-section is σ_bf ∼ 10⁻¹⁷ cm² with the cross section decreasing for smaller wavelengths. For ξ = 14, the number density is bounded n₂ ≲ 10² cm⁻³ and a characteristic length is L ∼ 10⁹–10¹⁰ cm. This yields an upper limit on the optical depth τ_bf ≲ σ_bfn₂L = 10⁻⁵–10⁻⁶ with smaller values of τ_bf for both larger impact parameters and shorter wavelengths. The reduction in flux due to bound-free absorption relative to Hα absorption should be proportional to the ratio of cross sections,

$$\begin{eqnarray} \left(\frac{\Delta F_\nu }{F_\nu }\right)_{\rm bf} & = & \frac{\sigma _{\rm bf}(\nu)}{\sigma _{{\rm H\alpha }}(\nu)} \left(\frac{\Delta F_\nu }{F_\nu }\right)_{\rm H\alpha }\qquad\quad\ \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} &&\qquad\ \ \approx 10^{-4} \left(\frac{\Delta F_\nu }{F_\nu }\right)_{\rm H\alpha } \approx 10^{-6}. \end{eqnarray} \tag{ 29 }$$

Within the context of our model, the transit depth due to the Balmer bound-free continuum absorption is too small to be observed for both HD 189733b and HD 209458b.

We have modeled the abundance of n = 2 hydrogen in hydrostatic atmospheres. We find that the Hα absorption can be explained by metastable 2s, similar to what is found in the ISM (e.g., Townes 1957). The dominant mechanism for the creation of 2s hydrogen is collisional excitation from the 1s state where it subsequently collisionally transitions to the 2p state and is finally radiatively de-excited. The 2s population dominates 2p throughout the atmosphere by two orders of magnitude for the parameters used in our study, although the specifics depend on the chosen value of J_Lyα. We do not model the spatial variation in the intensity of Lyα, instead choosing a constant value. This assumption allows us to estimate an upper limit on the 2p abundance. Since 2p remains negligible compared to 2s for our chosen J_Lyα, we can assume that it should similarly be negligible in the case where resonance scattering has depleted the available Lyα photons. Unlike n_2p, n_2s has limited dependence on the radiation field, instead depending strongly on the gas temperature through the exponential dependence found in C_{1s → 2s}. We find that the data are best fit by an atomic hydrogen layer approximately 500–1000 K hotter than our ξ = 1 case, which corresponds to ξ = 14. Caution should be used in the physical interpretation of this value due to the simplified heating and cooling present in our model.

Because of the strong dependence on temperature, Hα should be considered a complementary probe to Lyα which is relatively insensitive to the gas temperature, probing instead all the atomic gas.

This model differs from the calculations of Jensen et al. (2012). Assuming that both Hα and Lyα are optically thin, they used Equation (2) to derive an excitation temperature of T_exc = 2.6 × 10⁴ K. Arguing that the levels will approach a Boltzmann distribution deeper in the atmosphere and the excitation temperatures found are not compatible with the expected gas temperature at these densities, they conclude that the absorption must occur in the planetary wind. For our ξ = 14 model, we find n₁/n₂ ∼ 10⁷–10¹⁰ which using Equation (5) gives

$$\begin{equation} T_{\rm exc} = \frac{118400\,{\rm K}}{\log \left(\frac{4n_1}{n_2}\right)} \simeq 4800\hbox{--}7000\,{\rm K}\,\,. \end{equation} \tag{ 30 }$$

This value for T_exc is lower than the temperature in the neutral layer and is significantly lower than the value quoted by Jensen et al. (2012) due to the significantly smaller n = 2 abundance in our model. The former is due to the collisional excitation being balanced by the ℓ-mixing reactions, not collisional de-excitation. The latter is due to the overestimation of the gas temperature due to the assumption of optically thin Lyα.

Tremblin & Chiang (2013) have proposed that the Hα absorption can be explained by the same mechanism they use to explain Lyα absorption, colliding planetary and stellar winds. This mechanism allows for Lyα absorption 100 km s⁻¹ from line center to be explained by the formation of hot (T ∼ 10⁶ K) neutral hydrogen generated through charge exchange with solar wind protons. If Hα is in fact probing the interface between the two winds, the Lyα and Hα lines should have comparable widths; however, the observed width for HD 189733b is ≃ 40 km s⁻¹. Tremblin & Chiang (2013) propose to explain this population of cool hydrogen in its n = 2 state through cooling of the initially hot population formed through charge exchange.

We find that observed absorption of Hα can be explained by our model with the primary signal coming from the neutral atomic layer. Within this layer, the abundance of n = 2 hydrogen is roughly constant, even though the overall abundance of hydrogen is increased by three orders of magnitude from the top of the layer to the bottom because the density increase is offset by a decrease in temperature. The transition to molecular hydrogen will lead to a downturn in the overall n = 2 abundance. As a result, there should not be significant contribution to Hα absorption within this layer.

Although many simplifications have been made, we have included the relevant physics and reproduced the transit signal observed by Jensen et al. (2012) for HD 189733b. There are, however, many avenues for improving the calculation. We do not include the cooling required to cause the transition to the molecular state. This results in our model having too high a temperature at the base radius. The inclusion of the necessary physics will eliminate the need for the intermediate boundary at R_b.

Our model also required a higher rate of heating than expected from UV models from the X-exoplanets Archive at the CAB (Sanz-Forcada et al. 2011). The inclusion of heavier atomic species allows for far-UV absorption higher in the atmosphere (Koskinen et al. 2012a, 2012b) which could reduce the need for the large heating rate.

Our model has assumed spherical symmetry. Due to the strong dependence of the n = 2 abundances on the temperature, temperature variations between the day and night side could induce strong day–night variation in n_2s and n_2p, which will be observable since the transit probes the day–night terminator. This can partially be negated by the redistribution of thermal energy by zonal winds; however, this provides more reason to model the signal in three-dimensions, not less.

We thank Adam Jensen and Seth Redfield for useful discussions, and for kindly providing the data for HD 189733b and HD 209458b, as well as Joshua N. Winn, Remy Indebetouw, and Mark Whittle for helpful discussions regarding Hα observations. We also thank the referee for providing constructive comments and suggestions. The authors acknowledge support from NSF AST-0908079 and NASA Origins NNX10AH29G grants.

In this section we investigate the effect of the radiation field on the Lyα cooling rate. The cooling rate Λ(T), as used in the text, ignores the possibility of heating due to collisional de-excitation. Allowing for this possibility, the cooling rate takes the form

$$\begin{eqnarray} \Lambda (T,J_{\rm Ly\alpha })n_{\rm e}n_{\rm 1s} &=& 10.2\,{\rm eV}(C_{\rm 1s\rightarrow 2s}n_{\rm 1s} + C_{\rm 1s\rightarrow 2p}n_{\rm 1s}\nonumber\\ &&-C_{\rm 2s\rightarrow 1s}n_{\rm 2s} - C_{\rm 2p\rightarrow 1s}n_{\rm 2p})n_{\rm e}. \end{eqnarray} \tag{ A1 }$$

This cooling rate depends on radiative excitation and de-excitation as well as recombination and photoionization implicitly through their effect on the level populations, n_2s and n_2p (see Equations (12) and (13)). For sufficiently low temperature or large J_Lyα, the collisional de-excitation rates can become comparable to the excitation rates, reducing the net cooling rate and, in extreme cases, result in net heating.

To this end, we solve Equations (13) and (12) for n_2s and n_2p given fixed values of n_e and T. In the absence of recombination, the solution for Λ(T, J_Lyα) is independent of n_1s, the dependence having been explicitly factored out in Equation (A1). Even if recombination is included, it can safely be ignored so long as its contribution to the 2s and 2p abundance remains small. For the 2s state, where collisional excitation is the dominant formation pathway, this results in the condition

$$\begin{equation} \frac{n_{\rm e}}{n_{\rm 1s}} \ll \frac{C_{\rm 1s\rightarrow 2s}}{\alpha _{\rm 2s}} = 5.17\times 10^5\left(\!\frac{T}{10^4\, {\rm K}}\!\right)^{0.075}\exp \left(-\frac{118400\,{\rm K}}{T}\right). \end{equation} \tag{ A2 }$$

For the 2p state, radiative excitation sets the occupation, resulting in the requirement,

$$\begin{eqnarray} n_{\rm e} &\ll & 5.19\times 10^{10}\,{\rm cm^{-3}} \left(\frac{T}{10^4\,{\rm K}}\right)^{0.32}\nonumber\\ &&\times \left(\frac{J_{\rm Ly\alpha }}{J_{\rm Ly\alpha ,0}}\right)^{1/2}\left(\frac{n_{\rm 1s}}{10^{10}\,{\rm cm^{-3}}}\right)^{1/2}. \end{eqnarray} \tag{ A3 }$$

For the model of HD 189733b, both of these conditions are satisfied, as is evidenced by Figures 12 and 13. In cases where the contribution of recombination is non-negligible, the effect would be to increase the overall abundance of 2s and 2p hydrogen, thus suppressing Lyα cooling.

We take the electron density to be n_e = 10⁸ cm⁻³, the typical abundance found in the atomic layer in our models. Although the solutions have dependence on n_e beyond that shown in Equation (A1), we have found that varying the abundance changes the results minimally.

Figure 16 shows the cooling function for four values of J_Lyα, with J_{Lyα, 0} = 3.42 × 10⁻¹¹ s⁻¹, the value used in our model above. For the case of J_Lyα = 0, we recover the standard cooling rate.⁶ For the three cases where J_Lyα is non-zero, we find that the cooling rate diverges from the standard case as the temperature decreases with the cooling rate becoming zero as the gas temperature approaches the excitation temperature of the radiation field, T_ex ∼ 118400 K/log ((2hν³/c²)J_Lyαg₁/g₂). For temperatures below this threshold, collisional de-excitations return energy to the gas faster than it is removed through collisional excitation, resulting in net heating. The weak dependence of the C_{2s → 1s} and C_{2p → 1s} results in the heating rate being roughly constant as the temperature decreases.

Zoom In Zoom Out Reset image size

We note that for our model the temperature did not drop below T = 6000 K so the inclusion of Lyα heating would not have changed our results.

The implications of this result, however, are that at low temperatures the Lyα cooling rate is possibly over-estimated and the existence of a transition to heating creates a temperature floor in the limit where Lyα is the only cooling mechanism. The specific location of this temperature floor depend on the value of J_Lyα, potentially making the details of the radiative transfer problem important. Lyα undergoes resonant scattering in the atmosphere, and due to the low absorption probability, is unlikely to attenuate exponentially, leaving the possibility that the intensity is non-negligible within the atomic hydrogen layer and could play a role as the transition to molecular hydrogen is approached.