Evidence for He i 10830 Å Absorption during the Transit of a Warm Neptune around the M-dwarf GJ 3470 with the Habitable-zone Planet Finder

The blue curve at the top of Figure 1 shows the weighted average of all the in-transit spectra from three transit epochs divided by the average of the out-of-transit spectra. The vacuum wavelengths of the He 10830 Å triplet in the planet's rest frame are marked by the orange vertical lines. Since the telluric correction as well as the sky emission line subtraction are not perfect, regions of sky/telluric contamination have enlarged error bars as discussed in Section 2. The weighted average of all combinations of transits taken two at a time is also shown in the curves below. The single peak observed at 10833 Å is not present in Transits 2 and 3, it is only seen in Transit 1. However, since 10833 Å was separated from the telluric band only during Transit 1, it did not get averaged out in the final weighted average. Since Transit 1 had only half the S/N compared to the other two transits we suspect this is an artifact due to photon and detector noise (see blue curves in Figure 5 for the individual ratio spectrum from Transit 1). For a null result comparison, one of the out-of-transit epochs is also shown in Figure 1. These data are processed the same way as other in-transit spectra, and we do not see any signatures of absorption inside the He 10830 Å triplet window.

Zoom In Zoom Out Reset image size

We measure the integrated EW of the detected broad absorption in the vacuum wavelength range (10831.9 Å to 10833.6 Å ) in Figure 1 to be 0.012 ± 0.002 Å. We ignored the J = 0 line (10832.06 Å in vacuum) of the triplet since its oscillator strength is only one-eighth of the other two lines combined, so it is below our continuum noise (the J = 0 line is the dashed orange line in Figure 1). The resonant scattering absorption of He i 10830 Å is directly proportional to the number density of He i in the ²_3S metastable triplet state. We use the curve of growth analysis in the optically thin regime to estimate the column density of $\mathrm{He}{\,}_{3}^{2}{\rm{S}}$ metastable atoms,

$$\begin{eqnarray}&&{N}_{\mathrm{He}{}_{3}^{2}{\rm{S}}}=\displaystyle \frac{\mathrm{EW}}{8.85\ast {10}^{-13}}\displaystyle \frac{1}{{\lambda }_{2}^{2}\times {f}_{{{ik}}_{2}}+{\lambda }_{3}^{2}\times {f}_{{{ik}}_{3}}}\,{\mathrm{cm}}^{-2},\end{eqnarray} \tag{ 1 }$$

where ${N}_{\mathrm{He}{}_{3}^{2}{\rm{S}}}$ is the column density of $\mathrm{He}{\,}_{3}^{2}{\rm{S}}$ metastable atoms in ${\mathrm{cm}}^{-2}$, EW is the measured combined EW of J = 1 and J = 2 lines of He 10830 Å triplet in units of Å, λ₂ and λ₃ are the wavelengths of the J = 1 and J = 2 triplet lines in units of microns, and ${f}_{{{ik}}_{2}}$ and ${f}_{{{ik}}_{3}}$ are the oscillator strengths of the lines taken from the NIST Atomic Spectra Database Lines Database (Drake 2006).

Substituting our measured EW, we obtain ${N}_{\mathrm{He}{}_{3}^{2}{\rm{S}}}=2.4\times {10}^{10}\ \,{\mathrm{cm}}^{-2}$. This measured EW is a flux-weighted average across the unresolved stellar disk during the transit. Hence, the column density we measure is also a stellar disk flux-weighted average column density of the exosphere during the transit.

We modeled the steady-state helium distribution in the outflowing exosphere of GJ 3470b using the formalism and state transition coefficients outlined in Oklopčić & Hirata (2018). Using hydrodynamic simulations Salz et al. (2016) show that the exosphere of GJ 3470b is not isothermal. We therefore do not use a Parker wind model for our atmospheric analysis, but solve for the steady-state helium distribution using the velocity and density field of the exosphere from the Salz et al. (2016) PLUTO-CLOUDY hydrodynamic simulations. To be internally consistent, we also used the GJ 3470 stellar SED used by Salz et al. (2016).

The X-ray spectrum of GJ 3470 we used from Salz et al. (2016) is calculated using the plasma emission model from CHIANTI (Dere et al. 1997, 2009) normalized to the X-ray luminosity of GJ 3470b estimated based on the stellar rotation period and the stellar mass (Pizzolato et al. 2003). The EUV flux (100 Å to ∼912 Å), which is critical for these calculations due to the strong dependence of helium ionization, comes from the scaling of the model-dependent Lyα flux (Linsky et al. 2014). The model-dependent Lyα flux of GJ 3470 itself is estimated by the Linsky et al. (2013) model based on the X-ray luminosity. See Section 2.2 of Salz et al. (2016) for more details. This model-dependent irradiance spectrum of GJ 3470 is the source of the biggest uncertainty in our calculations. However, the model-dependent Lyα flux Salz et al. (2016) derive, and the value we adopt here, is consistent with the Lyα flux measured by Bourrier et al. (2018) during three transits of GJ 3470b using the Space Telescope Imaging Spectrograph instrument on board the HST.

As most of the free electrons in the atmospheres of exoplanets come from ionized hydrogen, to obtain the electron density of the GJ 3470 atmosphere, we first solved for the steady-state of hydrogen following Oklopčić & Hirata (2018). To obtain the steady-state distribution of helium atoms, we briefly discuss here the relevant system of integro-differential equations from Oklopčić & Hirata (2018), given by,

$$\begin{eqnarray}\begin{array}{rcl}v\displaystyle \frac{\partial {f}_{1}}{\partial r} & = & (1-{f}_{1}-{f}_{3}){n}_{e}{\alpha }_{1}+{f}_{3}{A}_{31}-{f}_{1}{{\rm{\Phi }}}_{1}{e}^{-{\tau }_{1}}\\ & & -{f}_{1}{n}_{e}{q}_{13a}\,+\,{f}_{3}{n}_{e}{q}_{31a}\,+\,{f}_{3}{n}_{e}{q}_{31b}+{f}_{3}{n}_{{{\rm{H}}}^{0}}{Q}_{31},\end{array}\end{eqnarray} \tag{ 2 }$$

and,

$$\begin{eqnarray}\begin{array}{rcl}v\displaystyle \frac{\partial {f}_{3}}{\partial r} & = & (1-{f}_{1}-{f}_{3}){n}_{e}{\alpha }_{3}-{f}_{3}{A}_{31}-{f}_{3}{{\rm{\Phi }}}_{3}{e}^{-{\tau }_{3}}\\ & & +{f}_{1}{n}_{e}{q}_{13a}\,-\,{f}_{3}{n}_{e}{q}_{31a}\,-\,{f}_{3}{n}_{e}{q}_{31b}-{f}_{3}{n}_{{{\rm{H}}}^{o}}{Q}_{31},\end{array}\end{eqnarray} \tag{ 3 }$$

where f₁ and f₃ are the fractions of neutral helium in the ground state of the singlet and triplet states, respectively, v is the velocity of the outflowing exosphere as a function of radius we adopt from Salz et al. (2016), n_e is the density of free electrons we obtained by solving the steady-state of hydrogen, and α₁ and α₃ are the helium recombination rates to singlet and triplet state from Osterbrock & Ferland (2006), respectively. The collision coefficients q_ijk and Q₃₁ are from Oklopčić & Hirata (2018), calculated using coefficients from Bray et al. (2000) and Roberge & Dalgarno (1982). The radiative decay rate A₃₁ from the helium metastable state to singlet state is adopted from Drake (1971). Φ₁ and Φ₃ are the effective photoionization rate coefficients calculated using,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{i}={\int }_{{\lambda }_{1}}^{{\lambda }_{2}}\displaystyle \frac{\lambda {F}_{\lambda }}{{hc}}{a}_{\lambda }d\lambda ,\end{eqnarray} \tag{ 4 }$$

where F_λ is the irradiated flux on GJ3470b described earlier in this section, a_λ is the photoionization cross section taken from Brown (1971) for the singlet state, and from Norcross (1971) for the triplet state. For the singlet state, the integral is evaluated up to the ionization wavelength of helium (504 Å), while for the triplet state the integral is computed in the interval starting at the Lyman limit¹⁹ (911.6 Å) to the metastable state ionization wavelength (2583 Å).

To calculate the optical depths τ₁ and τ₃ as a function of radius for the singlet and metastable states of helium, we compute the following integrals,

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{1} & = & {a}_{\mathrm{oH}}\displaystyle \frac{0.9}{1.297\ast {m}_{p}}{\int }_{r}^{\infty }(1-{f}_{{\rm{H}}+})\rho (r){dr}\\ & & +{a}_{\mathrm{oHe}}\displaystyle \frac{0.1}{1.297\ast {m}_{p}}{\int }_{r}^{\infty }{f}_{1}\rho (r){dr},\end{array}\end{eqnarray} \tag{ 5 }$$

and,

$$\begin{eqnarray}&&{\tau }_{3}={a}_{\mathrm{oHe}2{\rm{S}}3}\displaystyle \frac{0.1}{1.297\ast {m}_{p}}{\int }_{r}^{\infty }{f}_{3}\rho (r){dr},\end{eqnarray} \tag{ 6 }$$

where m_p is the proton mass, ρ(r) is the density of the GJ3470b exosphere from Salz et al. (2016), f_H+ is the fraction of ionized hydrogen we obtained by solving the steady-state of hydrogen. a_oH and a_oHe are GJ3470b's irradiation flux-averaged hydrogen and helium ionization cross sections in the wavelength range up to helium ionization (504 Å). a_oHe2S3 is the flux-averaged cross section of $\mathrm{He}{\,}_{3}^{2}{\rm{S}}$ in the wavelength range 911.6 Å (hydrogen ionization) to 2583 Å ($\mathrm{He}{\,}_{3}^{2}{\rm{S}}$ ionization).

We solve Equations (2) and (3) iteratively as an initial boundary value partial differential equation by starting with an initial estimate for the optical depth integral term (Equations (5) and (6)). Using the resulting solution, we updated the integral term in Equations (5) and (6). The final solutions converged within one or two iterations since the integral term has impact only in the high opacity base region of the outflowing atmosphere. The differential equation of $\mathrm{He}{\,}_{3}^{2}{\rm{S}}$ metastable atoms (Equation (3)) is very stiff for high densities of GJ 3470b, and we therefore used a Radau solver with adaptive dense gridding (Hairer & Wanner 1996).

The one-dimensional radial density distribution of $\mathrm{He}{\,}_{3}^{2}{\rm{S}}$ metastable atoms obtained above is then integrated assuming a spherical exosphere to obtain column densities at different impact parameters from the planet (Figure 2(b)). We expect the model to fail beyond the Roche-lobe due to complex stellar wind interactions. We estimate the volume equivalent Roche-lobe radius of the GJ 3470b system to be 3.12 R_p following Eggleton (1983). The teardrop shaped Roche-lobe's extent on the star–planet axis is 5.96 R_p (Salz et al. 2016). Both of these points are explicitly highlighted in Figure 2. Assuming an impact parameter of b = 0.47 (Dragomir et al. 2015), Figure 2(a) shows the 2D projection of the system during the midpoint of the transit. In Figure 2(a) we also illustrate the stellar disk on top of the column density map of $\mathrm{He}{\,}_{3}^{2}{\rm{S}}$ metastable atoms. To estimate the flux-averaged column density, we used a quadratic limb-darkening model where we calculated the limb-darkening coefficients using the EXOFAST web-applet²⁰ using the stellar effective temperature, metallicity, and surface gravity from Table 1, in the J band. This flux-weighted average column density is dependent on how far the exosphere extends beyond the Roche-lobe; the predicted background stellar disk flux-weighted average column density of $\mathrm{He}{\,}_{3}^{2}{\rm{S}}$ metastable atoms as a function of the extent of GJ 3470b's exosphere is shown by the dashed line in Figure 2(b)).

Zoom In Zoom Out Reset image size

3.3.1. Line Profile Model

The broad absorption we present evidence for during the transit of GJ 3470b is possibly the broadest He 10830 Å absorption reported in the literature to date. It spans a velocity range of −36 to +9 km s⁻¹. Figure 3(a) shows the outflowing exosphere velocity from Salz et al. (2016). The wind reaches only velocities up to ∼10 km s⁻¹ inside 3.12 R_p—the volume equivalent Roche-lobe radius. Both the volume equivalent Roche lobe radius as well as the Roche lobe extent on the star–planet axis (5.96 R_p) is marked by vertical lines in the plot.

Zoom In Zoom Out Reset image size

The line-of-sight velocity of the helium atoms we see during the transit is the sum of the hydrodynamic driven wind velocity plus the stellar radiative acceleration, and the planet's orbital radial velocity at the instant of their escape from the planet's gravitational potential. For example, helium atoms escaping the planet potential ∼2.5 hr before the transit will have a line-of-sight velocity of ∼−30 km s⁻¹ during the transit due to change in the line-of-sight component of planet's orbital velocity. Figure 3(b) shows the time a helium atom will take to travel the distance based on the wind model after it leaves the Roche lobe. Atoms that escape the potential 2.5 hr before the transit will be at ∼8 R_p, which is well within the projected stellar disk.

For calculating the line profile predicted by our model, we need to calculate the line-of-sight velocity field of the GJ 3470b exosphere at all impact parameter positions and radial distances. For the region inside the volume equivalent Roche Lobe, we consider only the line-of-sight projection of the radial hydrodynamic wind velocity. For the region outside the Roche lobe, we added the line-of-sight component of the hydrodynamic wind to the difference in the planet's radial velocity between present and the time at which gas at a certain altitude escaped the Roche lobe. Figure 3(c) shows this velocity field. The observer on Earth is toward the left side of the axis and the host star is toward the right side. The negative line-of-sight velocity implies gas is moving toward the observer. In reality, there will be additional blueshift due to the interaction with stellar wind. We do not include any stellar wind induced blueshift of the exosphere in our line model.

To model the line profile, we used the Voigt profile. The modeling methodology is mostly similar to Oklopčić & Hirata (2018), the only difference is in not assuming a spherically symmetric velocity field, and using the GJ 3470b system's impact parameter and stellar disk limb-darkening model while integrating the absorption line profile across the stellar disk. The temperature of the gas for the thermal Gaussian component of the modeled Voigt profile was taken to be 7000 K (Salz et al. 2016; Bourrier et al. 2018). The He 10830 line's Einstein coefficient for the Lorentzian component of the Voigt profile was taken from NIST Atomic Spectra Database Lines Database (Drake 2006). The oscillator strengths for the cross section calculation of each triplet line is also taken from the NIST Database. These cross sections, along with the radial density distribution of the metastable helium we obtained in our model (Section 3.3) were integrated to obtain the optical depth at different impact parameter distances from the planet. Density of metastable helium outside 13 R_pl (radial limit of our 1D model) was set to zero in our integrals. Using a quadratic limb-darkening model, we averaged the transmission curve across the disk to obtain the net transmission profile of the He 10830 triplet lines in GJ 3470b during the transit. Since the measured column density is less than the model predicted column density by a factor of 10 (see Section 3.4), we scaled down the density distribution from our model also by a factor of 10. This predicted line profile is overplotted over all the transit and out of transit ratio spectra in Figure 1.

3.3.2. Limits of the Model

The simplified one-dimensional model we presented here is to check whether the strength of our claimed helium 10830 Å absorption in the exoplanet atmosphere is within the limits of expected physical conditions. A true three-dimensional wind model where the stellar irradence is only on one side of the planet is known to have lesser density on average in all angles by a factor of 4 than the one-dimensional symmetric wind model we used (Stone & Proga 2009).

The other major source of uncertainty in our model is the irradiated flux in the wavelengths shorter than helium ionization (λ < 512Å). As briefly outlined in the section above, it is estimated by a chain of empirical models. We expect the helium ionization flux to have a significant influence on the population of metastable He atoms. To study the impact, we simulated how the predicted column density of $\mathrm{He}{\,}_{3}^{2}{\rm{S}}$ metastable atoms changes when we suppress the irradience at wavelength shorter than 504 Å by different factors.²¹ Results are shown in Figure 4. Since we wanted to probe only the impact on the steady state of helium atoms, the underlying velocity field and density used in the differential equation was kept fixed, and is the same as the Salz et al. (2016) hydrodynamic model for the original irradiance. It is very instructive to see that the observed column density of metastable helium is reduced proportionally to the reduction in the helium ionization flux.

Zoom In Zoom Out Reset image size

In real systems, since the mass outflow is proportional to the energy absorbed in the lower atmosphere from irradiance, there will be an additional proportional reduction in density when the EUV flux irradiation is reduced. Hence, the combined effect of reduction in helium ionizing radiation is quadratic on the column density of meta-stable helium atoms.

A few more ways the underlying density of the wind we used in our simulation from Salz et al. (2016) can be impacted is summarized in Table 2 of Salz et al. (2016).

As described in Section 3.2, the measured flux-weighted column density of He ${}_{3}^{2}{\rm{S}}$ is ${N}_{\mathrm{He}{}_{3}^{2}{\rm{S}}}=2.4\times {10}^{10}\,{\mathrm{cm}}^{-2}$. This is a factor of 10 less than predicted by our one-dimensional hydrodynamic model described in Section 3.3. However, as discussed in Section 3.3.2, if we correct our simulation for the 3D model simplification to the 1D model by the factor of 4 (Stone & Proga 2009), and assume the helium ionization radiation is lesser by a factor of 1.6 (which is well within the error of the empirically derived EUV SED of GJ3470 in Salz et al. 2016), we reduce the predicted metastable helium by a factor of $4\times {1.6}^{2}=10$. Hence, our observations are consistent with our simplified model under the caveats outlined in Section 3.3.2.

We discuss a few possible reasons for differences in the shape of the simple model and observed line profile below.

Stellar Activity: Some M dwarfs can be active, and exhibit time variability in the He 10830 absorption line (we note that not all M dwarfs are active). In comparison to other transmission spectroscopy lines like Ca ii K, Na i D, and Hα, active regions typically have He 10830 in absorption (Cauley et al. 2018). Thus, when a planet transits an active spot region on the host star, He 10830 is more likely to result in an emission signature, than an absorption signature. Thus He 10830 absorption signatures in-transit spectroscopy are more likely to originate in the exosphere than from the stellar chromosphere. However, under high density chromosphere conditions (Andretta & Giampapa 1995), or active flares, it is possible to obtain He 10830 in emission from the recombination cascade of ionized helium (Huenemoerder & Ramsey 1987). We do not see any indication of flares during transit in the Ca II IR triplet lines. If the absorption we detected was due to GJ 3470b transiting a very active emission spot, we would expect the absorption signal to not be present in all the in-transit spectra taken (see Figure 5). We have a total of nine in-transit spectra, (three spectra from three transits each). Since the phase of each spectrum is slightly different, it is improbable that there was an unusually strong emission spot of He 10830 during each of our transits. Hence, stellar activity cannot explain the absorption signal we detected during transit of GJ3470b. Stellar activity, however, could probably explain some of the emission structures, and there by the net absorption strength variations we measured across the three transits.

Zoom In Zoom Out Reset image size

Telluric Absorption and Sky Emission: Systematics in telluric and sky emission line modeling could potentially cause profile variations in the regions with high error bars in Figure 1. The regions with small error bars are free of any telluric or sky emission line contamination in at least one epoch out of the three transit observations. Hence, imperfect telluric line modeling errors cannot explain the absorption signature.

Stellar Wind: Even though our simple model, which includes the orbital dynamics and the hydrodynamic outflow velocity, could explain the large blueshifted profile seen in the data, it does not fully explain the shape of the detected line. Some component of this discrepancy could be variability; however, we suspect the stellar wind to play a major role. Since stellar wind accelerates the exosphere toward Earth like a comet tail, it could further blueshift the exosphere at the base of the exosphere skewing the line profile to shorter wavelengths. Further monitoring of the GJ 3470b's He 10830 signature during transits will reveal greater insights on the system.