Alkaline exospheres of exoplanet systems: evaporative transmission spectra

ABSTRACT

Hydrostatic equilibrium is an excellent approximation for the dense layers of planetary atmospheres, where it has been canonically used to interpret transmission spectra of exoplanets. Here, we exploit the ability of high-resolution spectrographs to probe tenuous layers of sodium and potassium gas due to their formidable absorption cross-sections. We present an atmosphere–exosphere degeneracy between optically thick and optically thin mediums, raising the question of whether hydrostatic equilibrium is appropriate for Na i lines observed at exoplanets. To this end we simulate three non-hydrostatic, evaporative, density profiles: (i) escaping, (ii) exomoon, and (iii) torus to examine their imprint on an alkaline exosphere in transmission. By analysing an evaporative curve of growth, we find that equivalent widths of |$W_{\mathrm{Na D2}} \sim 1{\!-\!} 10\, \mathrm{m\mathring{\rm A}}$| are naturally driven by evaporation rates ∼10³−10⁵ kg s⁻¹ of pure atomic Na. To break the degeneracy between atmospheric and exospheric absorption, we find that if the line ratio is D2/D1 ≳ 1.2 the gas is optically thin on average roughly indicating a non-hydrostatic structure of the atmosphere/exosphere. We show this is the case for Na i observations at hot Jupiters WASP-49b and HD189733b and also simulate their K i spectra. Lastly, motivated by the slew of metal detections at ultra-hot Jupiters, we suggest a toroidal atmosphere at WASP-76b and WASP-121b is consistent with the Na i data at present.

line: profiles, radiative transfer, techniques: spectroscopic, planets and satellites: atmospheres

1 INTRODUCTION

The alkali metals, sodium (Na i) and potassium (K i), have long been predicted to be observable spectroscopically during an exoplanet transit (Seager & Sasselov 2000; Brown 2001; Hubbard et al. 2001). While this launched the study of extrasolar atmospheres, we remark here that alkaline observations are also consistent with extrasolar exospheres of diverse geometries in terms of the required source rates (Johnson & Huggins 2006; Oza et al. 2019b). This subtle degeneracy between collisional atmospheres and collisionless exospheres is due to the large absorption cross-sections of the Na and K atoms, enabling small column densities of gas to illuminate optically thin gas.

While transmission spectra can be simulated analytically (Lecavelier Des Etangs et al. 2008; de Wit & Seager 2013; Bétrémieux & Swain 2017; Heng & Kitzmann 2017; Jordán & Espinoza 2018; Fisher & Heng 2019), these approaches rely on the assumption of an atmosphere being in local hydrostatic equilibrium, an assumption that generally holds in dense atmospheric layers where gravity and pressure gradients are the dominant forces. Hydrostatic equilibrium becomes less accurate in more tenuous atmospheric layers, and can break down in certain circumstances such as in an escaping atmospheric wind. In fact, soon after the first detection of Na i at HD209458b [interpreted as an exoplanet atmosphere by Charbonneau et al. (2002), and recently contested by Casasayas-Barris et al. (2020) attributing the signal to the Rossiter-McLaughlin effect and centre-to-limb variations], Vidal-Madjar et al. (2003) detected neutral hydrogen beyond the Hill sphere and interpreted this as an evaporating component of the atmosphere. Close-in exoplanet atmospheres were hence observed to be non-hydrostatic, and shown to be hydrodynamically evaporating ∼10⁷ kg s⁻¹ of gas due to XUV-driven escape (e.g. Murray-Clay, Chiang & Murray 2009). Considerable endogenic modelling of evaporative transmission spectra is indeed underway, largely aiming at exospheric signatures of atomic hydrogen (Bourrier & Lecavelier des Etangs 2013; Christie, Arras & Li 2016; Allan & Vidotto 2019; Murray-Clay & Dijkstra 2019; Wyttenbach et al. 2020), but also by atomic helium (Oklopčić & Hirata 2018; Lampón et al. 2020), atomic magnesium (Bourrier, Lecavelier des Etangs & Vidal-Madjar 2015), and ionized magnesium (Dwivedi et al. 2019). The problem of a hydrostatic assumption is further amplified and fundamentally different if the gas were to be detached from the planet. Such an exogenic¹ alkaline source, such as an outgassing satellite or a thermally desorbing torus, would not be in hydrostatic equilibrium. The impact of exogenic sources on transmission spectra, one of the key questions we seek to answer in this study, has not yet been investigated until present. In the following, we shall use evaporative and non-hydrostatic interchangeably.

The retrieval of atmospheric parameters from observations using various techniques such as χ²-minimization (e.g. Madhusudhan & Seager 2009), Bayesian analysis (e.g. Madhusudhan et al. 2011; Lee, Fletcher & Irwin 2012; Benneke & Seager 2013), or advanced machine-learning methods (e.g. Márquez-Neila et al. 2018; Hayes et al. 2020) is generally carried out within a hydrostatic set-up. Various hydrostatic retrieval codes exist in the community (e.g. Madhudsudhan & Seager 2009; tau-rex i Waldmann et al. 2015; bart Blecic, Dobbs-Dixon & Greene 2017; exo-transmit Kempton et al. 2017; atmo Goyal et al. 2018; ^πη Pino et al. 2018,²; aura Pinhas et al. 2018; helios-t Fisher & Heng 2018; platon Zhang et al. 2019; Brogi & Line 2019; Fisher & Heng 2019), mostly differing in terms of their treatment of atmospheric chemistry (ranging from constant abundances over chemical equilibrium to atmospheres in disequilibrium), assumptions of temperature profiles (using isothermal or parametric profiles, or invoking radiative-convective equilibrium modelling), and retrieval routines. For a comprehensive overview of different codes and retrieval techniques (see Madhusudhan 2018), here we briefly point out some of the codes that pose exceptions to the common assumption of an atmosphere in hydrostatic equilibrium. The nemesis (Irwin et al. 2008) and chimera(Line et al. 2013) codes are flexible in terms of their atmospheric structure such that they can process arbitrary pressure profiles. The more recent tau-rex iii code (Al-Refaie et al. 2019) allows retrievals in chemical equilibrium or disequilibrium for custom (non-hydrostatic) pressure profiles (note that the default retrieval set-up is still hydrostatic).

In this study, we examine evaporative transmission spectra of the Na i doublet (⁠|$\lambda _{\mathrm{D2}}=5889.95\, \mathring{\rm A},\, \lambda _{\mathrm{D1}}=5895.92\, \mathring{\rm A}$|⁠) and the K i doublet (⁠|$\lambda _{\mathrm{D2}}=7664.90\, \mathring{\rm A},\, \lambda _{\mathrm{D1}} =7698.96\, \mathring{\rm A}$|⁠). The lines when viewed in absorption are extremely bright due to resonance scattering (Brown & Yung 1976; Draine 2011) off the sodium atoms. The resonant scattering cross-section is large, producing considerable absorption in highly tenuous columns of gas at pressures where one expects large deviations from hydrostatic equilibrium. Fortunately, astronomers and planetary scientists have had decades of fundamental understanding of the physics of the Na and K alkaline resonance lines by directly observing comets and moons in situ from within our Solar system (Oza et al. 2019b, see table 1 and references therein). The most spectrally conspicuous aforementioned evaporative sources are Jupiter’s and Saturn’s moons Io & Enceladus, motivating the mass-loss model of an evaporating exomoon shown to be roughly consistent with several extrasolar gas giant planets today (Oza et al. 2019b). Coupling a metallic evaporation model (termed dishoom) to a radiative transfer code capable of treating non-hydrostatic profiles (termed prometheus) is the crux of our present study on evaporative transmission spectra. This coupling permits a holistic approach to transit spectra, capable of breaking endogenic–exogenic degeneracies in alkaline exospheres today.

We use a custom-built radiative transfer code to simulate high-resolution transit spectra in the sodium doublet for four scenarios, seeking the precise imprint of an evaporative transmission spectrum of an alkaline exosphere in regards to a canonical hydrostatic atmosphere. We present this code in Section 2.1. The mass-loss model is laid out in Section 2.2. Our four examined scenarios, corresponding to a particular geometry and spatial distribution of the sodium atoms, are presented in detail in the following sections:

• Section 2.3 – Hydrostatic: A spherically symmetric hydrostatic atmosphere.

• Section 2.4 – Escaping: A spherically symmetric envelope undergoing atmospheric escape.

• Section 2.5 – Exomoon: A spherically symmetric cloud sourced by an outgassing satellite.

• Section 2.6 – Torus: An azimuthally symmetric torus sourced by a satellite or debris.

We describe each of these scenarios with a particular number density profile n(r) (since an exospheric collisionless gas cannot be described with a pressure profile) and two free parameters. These three components fully determine our simulated transit spectra. Of course, an exhaustive description of metals at a close-in gas giant system would involve a careful treatment of atmospheric collisional processes (e.g. Huang et al. 2017), exospheric physical processes (Leblanc et al. 2017), an atmospheric escape treatment in 3D (e.g. Debrecht et al. 2019), and also the ability to track ions in the presence of a magnetic field (e.g. Carnielli et al. 2020). Since we reduce every scenario to a number density profile with two free parameters, our model is heuristic at present. Hence, this study is not targeted at providing an exhaustive model of hot Jupiter exospheres with the corresponding transit spectra, but rather encourages a novel approach by elucidating fundamental differences between hydrostatic and non-hydrostatic assumptions.

We use the synergy of our two distinct approaches (radiative transfer and mass-loss) to gain a physical intuition on the four density scenarios described above starting in Section 3. We find that evaporative sodium profiles naturally allow vastly more extended, yet tenuous distributions of the absorbing atoms than hydrostatic profiles. This property of evaporative sodium can lead to absorption in a primarily optically thin regime, resulting in high-resolution transit spectra differing from spectra computed within a hydrostatic framework. In Section 3.4, we provide a simple diagnostic gathered from physics of the interstellar medium (e.g. Draine 2011) to determine if an alkaline gas is optically thin or thick at a transiting exoplanet. Given that several hot Jupiters observed in high-resolution appear to reveal an optically thin regime, we simulate the spectral imprint for each evaporative scenario at HD189733b in Section 4. The forward model in this section serves to demonstrate the differences between hydrostatic and evaporative transmission spectra by coupling mass-loss a priori.

For those more acquainted with inverse modelling, we use the observed transit spectra of the hot Jupiters WASP-49b and HD189733b (system parameters given in Table 1) to portray the behaviour of evaporating alkalis in parameter space (Section   5). In our observational analysis, best-fitting parameters are found within the radiative transfer code using a bounded χ²minimization. We then check the plausibility of these retrieved parameters by comparing to the calculated values from the alkali mass-loss model. We choose to comparatively analyse these two planets, as they were observed by the same spectrograph (HARPS) and reduced by the same authors (Wyttenbach et al. 2015, 2017) enabling a consistent comparison to the data of both hot Jupiters. These planets represent the first high-resolution Na i detections and therefore have been extensively studied by independent groups (Louden & Wheatley 2015; Cubillos et al. 2017; Huang et al. 2017; Pino et al. 2018; Fisher & Heng 2019) compared in Section 6.1. Furthermore, these planets were shown to be able to host exogenic sources of alkali metals (i.e. satellites) in terms of tidal stability and average column densities (Oza et al. 2019b). For our study of alkaline exospheres, we focus on the physics of evaporating Na i acknowledging the behaviour of evaporating K i is nearly identical, simulated in Section 6.2. We discuss the imminence of metal detections at ultra-hot Jupiters (WASP-76b & WASP-121b) in the context of toroidal atmospheres in Section 6.3 and conclude in Section 7.

Table 1.

System parameters of WASP-49b (Wyttenbach et al. 2017) and HD189733b (Wyttenbach et al. 2015), lifetimes of neutral sodium (τ_Na) are estimated from Huebner & Mukherjee (2015). γ denotes the systemic velocity of the exoplanetary systems, which is important for the comparison of the observational data to our simulated spectra.

Parameter	WASP-49b	HD189733b
R* (R⊙)	1.038	0.756
R0 (RJ)	1.198	1.138
Mp (MJ)	0.399	1.138
Teq (K)	1400	1140
γ (km s−1)	41.7	−2.28
tNa (s)	241	1010

2 METHODS

2.1 prometheus: Alkaline transmission spectra of atmospheres and exospheres

We simulate transit spectra using a simple and flexible custom-built python code: prometheus. Our code computes absorption either in the Na i or K i doublet for an exoplanet system in transit geometry. For a spherically symmetric system (our first three scenarios), we divide the atmosphere using a linear grid in x-axis and a logarithmic one in z-direction with adjustable spatial resolutions to adapt to each scenario (see Fig. 1). For the torus scenario exhibiting azimuthal geometry, we also linearly discretize the y-axis that affects the computation of the transit spectra (equations 1 and 3; see Section 2.6 for details). The transit spectrum is computed in the following two steps. First, we calculate the optical depth along the chord (x-axis, Fig. 1) at a certain altitude z and wavelength λ:

$$\begin{eqnarray*} \tau (z,\lambda)=2\int _0^{\infty }\! n(r) \sigma \bigl (\lambda ,T(r)\bigr) \chi _{\mathrm{i}}(r) \, \mathrm{d}x, \end{eqnarray*}$$

(1)

given that |$r=\sqrt{x^2+z^2}$| and where χ_i(r) is the volume mixing ratio³ of the neutral absorber (either Na i or K i), σ(λ, T(r)) the absorption cross-section and n(r) the number density profile that depends on the scenario (equations 8, 9, 12, 15). n(r) is generally calculated assuming hydrostatic equilibrium in this study, we investigate how different spatial distributions of the absorbing atoms – i.e. different number density profiles – compare to the canonical hydrostatic atmosphere in terms of transmission spectra.

Figure 1.

Sketch of an exoplanet in transit geometry. In this scheme, the three coordinate axes build a regular right-handed coordinate system and the host star is located at (−a_P, 0, 0), with a_P being the orbital radius of the exoplanet. Our four scenarios exhibit two different geometries: spherical symmetry (the orange cloud) and azimuthal symmetry (the red torus). For a detailed scheme of the architecture of a sodium exosphere with exogenic sources illustrated (see fig. 4 of Oza et al. 2019b).

In addition, our code can be coupled to a chemistry code such as fastchem(Stock et al. 2018) to calculate these mixing ratio profiles, but we only examine constant mixing ratio profiles in the scope of this paper. Given that the mixing ratio χ_i for a species i can be described as χ_i(r) = χ_i(r) × (1 − f_ion(r)), the assumption of a constant mixing ratio of the absorber implies a constant ionization fraction f_ion(r) for the absorbing species, given that the total (neutral plus ionized) mixing ratio of the absorber is not expected to vary significantly. At present, we focus solely on the absorption of alkali metals, as they dominate the absorption in the narrow wavelength regions around the doublets. We also include Rayleigh scattering from a background H₂ atmosphere, but find this to be negligible and therefore do not include it for the present application.

We compute the absorption cross-section σ(λ, T) as the sum of the individual D2 and D1 absorption lines. These lines are modelled as Voigt profiles using scipy.special.wofz. Line broadening is included in our model due to the gas temperature T(r) (Doppler broadening). Since the absorbing atoms reach velocities significantly above thermal speeds for our three evaporative, non-hydrostatic scenarios, we must also consider this non-thermal Doppler broadening. We incorporate this broadening by treating the average velocity |$\bar{v}_{\mathrm{i}}$| of the absorbing atoms as an effective line temperature:

$$\begin{eqnarray*} T(r)=\frac{\pi m_{\mathrm{i}}}{8 k_{\rm B}}\bar{v}_{\mathrm{i}}(r)^2, \end{eqnarray*}$$

(2)

where m_i is the atomic mass of the absorber. The line is then broadened using Voigt profiles, with a Doppler broadening parameter given by this line temperature for the three evaporative scenarios. We remark that the alkali D lines have strong sub-Lorentzian wings due to pressure broadening (e.g. Heng et al. 2015; Allard et al. 2019), which is an observable effect in low-resolution transmission spectroscopy (e.g. Nikolov et al. 2018). Since we investigate high-resolution spectra in this analysis, the probed pressures are comparatively lower. Upon verifying the influence of pressure broadening in our models, we find the effect negligible for our simulated spectra. In line with Heng et al. (2015), we neglect pressure broadening throughout this work.

Secondly, we average all chords over the stellar disc to obtain the flux decrease due to a canonical atmosphere and more comprehensively, any exospheric source of absorption above the reference radius R₀, where F_out denotes the flux out of transit and F_in is the flux during transit of the hot Jupiter:

$$\begin{eqnarray*} \Re :=\frac{F_{{\rm in}, \lambda }}{F_{{\rm out}, \lambda }}=\frac{2}{R_{\ast }^2-R_0^2}\int _{R_0}^{R_{\ast }}\!z {\rm e}^{-\tau (z,\lambda)}\, \mathrm{d}z. \end{eqnarray*}$$

(3)

This scheme is a general prescription for the computation of transmission spectra. To compare our calculations to observations WASP-49b (Wyttenbach et al. 2017) and HD189733b (Wyttenbach et al. 2015), we use a convolution, binning, and normalization routine for all simulated transit spectra as in Pino et al. (2018, see Section 2.7 for details). While our code is comparatively simple in terms of atmospheric chemistry, wavelength coverage, line broadening, and absorbing species, it has the stark advantage that the number density profile can be arbitrarily defined without relying on the assumption of the atmosphere being in local hydrostatic equilibrium. This allows the simulation of transit spectra in two endogenic (hydrostatic and escaping) and two exogenic scenarios (exomoon and torus).

2.2 dishoom: Alkaline mass-loss to atmospheres and exospheres

We couple a semi-analytic metal mass-loss model dishoom(Desorbing Interiors via Satellite Heating to Observing Outgassing Model) described comprehensively in Oza et al. (2019b), Section 3 (Jupiter’s Atmospheric Sodium), and Section 4 (Jupiter’s Exospheric Sodium) to prometheus in order to simulate evaporative transmission spectra for the Na i and K i lines.

2.2.1 Observing outgassing during a planetary transit

An outgassing source undergoing mass-loss can produce a prodigious quantity of foreground gas capable of generating a spectral signature during transit. Our escape model generalizes outgassing of planetary surfaces to include a wide-range of phenomena capable of producing spectral signatures as observed in the Solar system. Escape of the so-called supervolatiles N₂, CH₄, and CO is included (see Johnson et al. 2015), as well as water products H₂O, O₂ (see Oza et al. 2019a) where the mechanisms leading to outgassing and eventual escape include solar heating, magnetospheric ion sputtering, and general space weathering. For close-in systems, like the ones studied here, nearly molten temperatures lead to volcanic and magma-products due to tidal heating or direct sublimation of grains (see Section 2.2.2). For a volcanically active system, Na and K are products of the parent molecules NaCl and KCl as observed at Io. The close-in systems are jeopardous in that the atomic lifetimes are strongly limited by photoionization. Following Huebner & Mukherjee (2015), we can estimate the atomic lifetime for a species i as |$t_{\mathrm{i}} \sim \frac{1}{k_{\mathrm{i}}}$| where k_i is the photoionization rate coefficient in s⁻¹:

$$\begin{eqnarray*} k_{\mathrm{i}, \lambda } = \int _{\lambda }^{\lambda + \Delta \lambda } \!\sigma _{\gamma }(\lambda) \Phi (\lambda)\, \mathrm{d}\lambda , \end{eqnarray*}$$

(4)

in the wavelength of interest λ + Δλ, a photoionization cross-section σ_γ(λ) and the spectral photon flux:

$$\begin{eqnarray*} \Phi (\lambda) = \frac{2 \pi c}{4 \lambda ^4 [\exp (hc/ \lambda k_{\rm B} T) -1 ]}, \end{eqnarray*}$$

(5)

for a Blackbody of temperature T (e.g. Rybicki & Lightman 1979). As we describe in Section 2.2.2, photoionization time-scales provide a critical lower limit to the atomic lifetime, as recombination processes due to the ambient plasma could boost these lifetimes. Independent of the mass-loss mechanism, we can write the general mass-loss rate for a source venting a species i depending on the number of atoms |$\mathcal {N}_{\mathrm{i}}$|⁠:

$$\begin{eqnarray*} \dot{M}_{\mathrm{i}}\approx m_{\mathrm{i}} k_{\mathrm{i}} \mathcal {N}_{\mathrm{i}} \approx m_{\mathrm{i}} t_{\mathrm{i}}^{-1} \mathcal {N}_{\mathrm{i}} . \end{eqnarray*}$$

(6)

The number of evaporating atoms |$\mathcal {N}_{\mathrm{i}}$| is then fed into the number density profiles for each evaporative mechanism, enabling a straightforward computation of an evaporative transmission spectrum.

2.2.2 Metallic mass-loss

In addition to the evaporation of Na i and K i, the code is equipped to estimate the destruction of rocky bodies of arbitary composition (e.g. MgSiO₃ and Fe₂SiO₄). This description of atmospheric loss holds observational relevance given the recent explosion of heavy metal detections at gas giants (Hoeijmakers et al. 2018, Hoeijmakers et al. 2019; Sing et al. 2019; Cabot et al. 2020; Gibson et al. 2020; Hoeijmakers et al. 2020). In our study of the alkali metals, we shall use chondritic ratios based on Fegley & Zolotov (2000) for Na and K, and rely on observations by Lellouch et al. (2003) and Postberg et al. (2009) to constrain the Na i abundance when estimating mass-loss rates due to thermal evaporation of silicate grains (i.e. Section 2.6).

In practice, the metal evaporation code computes different regimes of escape (e.g. thermal and non-thermal) due to several heating mechanisms (e.g. XUV and tidal heating) for a close-in irradiated body. The dominant mass-loss rate is then either supplied to prometheus to generate a forward-model transit spectrum, or compared to retrieved mass-loss rates and velocities in the inverse modelling as a plausibility check if the different evaporative scenarios can indeed provide the required absorber source rates. The dominant mass-loss mechanism varies between the different evaporative scenarios. We show the equations we use to estimate |$\dot{M}_{\mathrm{i}}$| in the respective subsection of the scenario (Sections 2.4–2.6). Generally, as data regarding the chemistry and plasma conditions are largely unknown at an exoplanetary system at present, we will find it useful to use scaling relations to Solar system bodies observedin situ (see Gronoff et al. 2020 and references therein for a full review on escape from Solar system and exoplanet bodies).

In our heuristic model, we provide upper limits to the required mass-loss rate based on the assumption that photoionization is the dominant process regulating the alkali lifetime. While to first order this is valid, two of our evaporative scenarios (Exomoon: 2.5 and Torus: 2.6 described) are directly analogous to the radiation environment of a gas giant magnetosphere. As observed and simulated in the Jupiter-Io plasma torus system, recombination and charge exchange could considerably extend the net lifetime of the alkali atoms (see Wilson et al. 2002 and references therein). Radiative recombination is additionally important for a close-in magnetosphere as ions sourced by photo- or electron-impact ionization of the evaporating atoms have the ability to accumulate in the magnetosphere so long as advection is small. We find, similar to Vidal-Madjar et al. (2013) for Mg at HD209458b, that n_e ∼ 10⁸ cm ⁻³ are required for electronic recombination of Na. In a toroidal B-field, we remark that ion-recombination and charge exchange can be effective to source the extended alkali clouds we describe here. In the absence of a toroidal B-field, a conservative ion density assuming charge neutrality based on the simulations by Dwivedi et al. (2019), ion densities of |$\sim 10^7 {\!-\!} 10^8\, \mathrm{cm}^{-3}$| were simulated up to ∼3R₀ validating a viable recombination mechanism at distances corresponding to our exogenic evaporating scenarios (Sections 2.5 and 2.6). Furthermore, we find that for the alkali atoms studied here, loss rates due to radiation pressure are small mostly due to the domineering ionization rates discussed above. Balancing the radiation pressure force with that of gravitation (e.g. Vidal-Madjar et al. 2003; Murray-Clay et al. 2009), upper limits to the mass-loss rate due to radiation pressure |$\dot{M}_{\mathrm{rad.press.}} \sim 10 {\!-\!} 10^4\,$| kg s⁻¹ are found for the range of sodium column densities studied in Section 3. In Tables 2 and 3, we qualitatively present the plasma processes regulating our evaporative scenarios and the average speeds of the alkali atoms |$\bar{v}_{\mathrm{i}}$|⁠.

Table 2.

Sodium plasma processes.

Plasma process	Reaction
Photoionization	Na + γ → Na+ + e−
Recombination	Na+ + e− → Na + γ
Charge exchange	Na+ + Na →Na* + Na+
Dissociative recombination	NaX+ + e− → Na* + X*

Table 3.

Sodium velocity distributions due to plasma-driven escape of Na at an exo-Io or exo-torus. The nominal velocity distribution gives the range of velocities based on thermal and non-thermal flux distributions used to model the analogous physical process in the Solar system. If an example is available, the reference is provided.

Physical process	Nominal velocity distribution
(Analog)	|$\bar{v}$|
Atmospheric sputteringa, b	1–30 km s−1
(Atmospheric escape at Io)	10 km s−1
Charge exchangeb	30–100 km s−1
(Sodium exosphere at Jupiter)	60 km s−1
Pickup ionsb	10–60 km s−1
(Io plasma torus)	74 km s−1
Resonant orbitc	10-42 km s−1
(Io motion)	17 km s−1
Thermal	0.6–2 km s−1
(Volcano temperature)	1.4 km s−1

^aSmyth & Combi (1988); Smyth (1992); ^b Wilson et al. (2002); ^c Cassidy et al. (2009).

The simple coupling via mass-loss rate and equation (6) we perform is to demonstrate the dramatic influence of mass-loss and number densities on an alkaline transmission spectrum. We remark that more robust hydrodynamic codes focusing on individual planets are especially suited for this problem, as has been shown in Ly α at HD189733b (Bourrier & Lecavelier des Etangs 2013; Christie et al. 2016), H α, H β, and H γ at KELT-9b (Wyttenbach et al. 2020), Mg i at HD209458b (Bourrier et al. 2015), Mg ii at WASP-12b (Dwivedi et al. 2019), and He i at HD209458b (Oklopčić & Hirata 2018; Lampón et al. 2020).

2.3 Hydrostatic scenario

This is the canonical scenario for transmission spectra. The number density profile is straightforward to derive using the definition of the pressure scale height: |$H(r) = \frac{k_{\rm B} T(r)}{\mu (r) g(r)}$|⁠, where μ(r) is the mean molecular weight and g(r) the local gravitational acceleration. Solving the equation of hydrostatic equilibrium dP(r)/dr = −n(r)μ(r)g(r) for number density, n(r), yields the hydrostatic number density profile:

$$\begin{eqnarray*} n_{\mathrm{hyd}} (r) = n_0 \frac{T_0}{T(r)} \exp \Bigl [-\!\int _{R_0}^r\!\frac{\mathrm{d}r^{\prime }}{H(r^{\prime })}\Bigr ], \end{eqnarray*}$$

(7)

where R₀ is a fixed reference radius, and quantities with subscript zero are evaluated at this reference radius. Although our code can process arbitrary temperature and mixing ratio profiles, we restrict ourselves to isothermal and vertically-mixed (i.e. with constant mixing ratios of the neutral absorber throughout the atmosphere) models, as we want to isolate the effect of varying the number density profiles within our different scenarios. Therefore, T(r) and μ(r) are constant and the number density profile in equation (7) can be integrated:

$$\begin{eqnarray*} n_{\mathrm{hyd}} (r) = n_0 \exp \bigl (\lambda (r)-\lambda _0\bigr), \end{eqnarray*}$$

(8)

where λ(r) is the Jeans parameter given by |$\lambda (r)=\frac{GM_p \mu }{r k_{\rm B}T}$|⁠. We use the canonical value of |$\mu =2.3\, \mathrm{amu}$| (corresponding to a mass mixing ratio x_He = 0.25) for our hydrostatic model. Using the number density profile of equation (8) with constant temperature and mixing ratios throughout the atmosphere leads to three free parameters in this scenario: temperature T, pressure at reference radius P₀, and absorber mixing ratio χ_i. However, P₀ and χ_i are mutually degenerate in our model since we neglect pressure broadening, collision-induced absorption and other absorbers (see Heng & Kitzmann 2017; Welbanks & Madhusudhan 2019 for a detailed discussion on this degeneracy). Therefore, we have the partial pressure of the absorber at the reference radius (|$P_{ 0,\rm{i}}=P_0\times \chi_{\rm{i}}$|) as a second free parameter. We show how these two parameters, T and P_{0, i}, affect the transmission spectrum for a hot Jupiter with planetary parameters of WASP-49b in Figs 2 and 3. We will find it useful to define auxiliary parameters that are easier to interpret than the free parameters, but do not fundamentally affect the transmission spectrum. For the hydrostatic scenario, we use P₀ as an auxiliary parameter (e.g. in Fig. 2), calculated from P_{0, i} under the assumption of a specific absorber mixing ratio. A summary of free, auxiliary, and fixed parameters for all scenarios can be found in Table 4.

Figure 2.

Simulation of high-resolution transit spectra in the sodium doublet for different reference pressures in the hydrostatic scenario. The shape of the transit spectra is not determined by P₀ but by the partial pressure of sodium at the reference radius, P_{0, Na}. For better readability, we fix χ_Na and vary P₀. We set |$\chi _{\mathrm{Na}}=1.7\,$| ppm (the solar mixing ratio) and |$T=3000\,$| K. The planetary parameters and data points are for WASP-49b, taken from Wyttenbach et al. (2017).

Figure 3.

Simulation of high-resolution transit spectra in the sodium doublet, for different temperatures in the hydrostatic scenario. We set |$\chi _{\mathrm{Na}}=1.7\,$| ppm (the solar mixing ratio) and |$P_0=1\,$| mbar. The planetary parameters and data points are for WASP-49b, taken from Wyttenbach et al. (2017).

Table 4.

Summary of free, auxiliary, and fixed parameters. The free parameters fundamentally define the transmission spectrum in our setting. We derive auxiliary parameters from them to improve the physical interpretability. For instance, we calculate the reference pressure P₀ in the hydrostatic scenario as auxiliary parameter from the free parameter P_{0, i} assuming a specific mixing ratio of the absorber.

Scenario	Free	Auxiliary	Fixed
Hydrostatic	|$P_{0,\mathrm{i}},\, T$|	P0	|$R_0,\, \mu$|
Escaping	|$\mathcal {N_{\mathrm{i}}},\, \bar{v}_{\mathrm{i}}$|	|$P_0,\, \dot{M}_{\mathrm{i}}$|	|$R_0,\, q_{\mathrm{esc}}$|
Exomoon	|$\mathcal {N_{\mathrm{i}}},\, \bar{v}_{\mathrm{i}}$|	|$\dot{M}_{\mathrm{i}}$|	|$R_s,\, q_{\mathrm{moon}}$|
Torus	|$\mathcal {N_{\mathrm{i}}},\, \bar{v}_{\mathrm{i}}$|	|$\dot{M}_{\mathrm{i}}$|	|$R_0,\, a_t,\, v_{\mathrm{ej}}$|

2.4 Escaping scenario

Hot Jupiters are observed to have escaping planetary winds (e.g. Vidal-Madjar et al. 2003). Currently, there are significant modelling efforts to simulate escaping atmospheres by either solving the Navier–Stokes equations or more computationally expensive yet robust, by directly solving the Boltzmann equation on a particle to particle basis (see Gronoff et al. 2020 for both approaches to atmospheric escape). While atmospheric escape is not the prime focus of the current study, we emphasize that the latter kinetic model is sorely needed at present for a full description of an exoplanet exosphere. Nevertheless, the hydrodynamical approach of escaping atmospheres has been applied to different hot Jupiters, computing pressure, temperature, and velocity profiles of the wind. Murray-Clay et al. (2009) found that at pressures |$P\gt 1\, \mathrm{nbar}$|⁠, the atmosphere can still be treated as hydrostatic, while at lower pressures the atmospheric structure strongly deviates from a hydrostatic profile. As a first-order approximation for the number density profile we use a power law with q_esc = 6 for our analysis, which is roughly in line with the profiles found for WASP-49b (Cubillos et al. 2017) and HD209458b (Murray-Clay et al. 2009):

$$\begin{eqnarray*} n(r)_{\mathrm{esc}} = n_0 \Bigl (\frac{R_0}{r}\Bigr)^{q_{\mathrm{esc}}}. \end{eqnarray*}$$

(9)

The escape index q_esc = 6 is indicative of ion-neutral scattering for an escaping neutral gas interacting with a plasma (Johnson 1990; Johnson et al. 2006a). The reference radius R₀ corresponds to the base of the wind in this scenario. As we want to isolate the different scenarios we neglect the hydrostatic layer below the base of the wind for the computation of the transit spectra. Since the hydrodynamical simulations generally treat only hydrogen and helium we do not know the mixing ratio profile of the absorber, χ_i(r). As in the hydrostatic scenario we set this profile to a constant value, meaning that we encounter a similar degeneracy between n₀ and χ_i as in the hydrostatic scenario.⁴ Hence, we have the number density of the absorber at the base of the wind, n_{0, i}, as a free parameter. We can convert this quantity into the total number of absorbing atoms as the volume integral of equation (9):

$$\begin{eqnarray*} \begin{aligned} \mathcal {N_{\mathrm{i}}}= \iiint _V \!\chi _{\mathrm{i}}\times n_{\mathrm{esc}}(r)\, \mathrm{d}V &=4\pi \chi _{\mathrm{i}} \int _{R_0}^{\infty }\!r^2 n_{\mathrm{esc}}(r)\, \mathrm{d}r \\ &=\left(\frac{4\pi }{q_{\mathrm{esc}}-3}\right) n_{0,\mathrm{Na}}R_0^3. \end{aligned} \end{eqnarray*}$$

(10)

Given that R₀ and q_esc are fixed we can use |$\mathcal {N}_{\mathrm{i}}$| and n_{0, i} interchangeably as free parameters. Since we use |$\mathcal {N}_{\mathrm{i}}$| to compute the required source rate, comprising the coupling between our codes (equation 6), we choose the total number of absorbing atoms in the systems as a free parameter. Since this quantity is difficult to interpret physically, we use the pressure at the base of the wind as auxiliary parameter. Assuming a certain absorbing mixing ratio and temperature we calculate this pressure as P₀ = n_{0, i}/χ_i × k_BT. We expect photoionization of the alkali metals to be significant in the escaping scenario. We remark here that while photoionization probably affects χ_i(r), it does not change the validity of setting χ_i(r) = const. as long as the ionized fraction is constant in the escaping wind. However, the retrieved value of |$\mathcal {N}_{\mathrm{i}}$| cannot be accurately converted into P₀, as we do not know the mixing ratio of the neutral absorber, χ_i. We show how P₀ affects the transmission spectrum for a hot Jupiter with planetary parameters of WASP-49b in Fig. 4. Apparently, the escaping transit spectrum depends strongly on the value of P₀ (compared to the hydrostatic transit spectrum, Fig. 2). We elucidate this behaviour in Section 4.

Figure 4.

Simulation of high-resolution transit spectra in the sodium doublet for different reference pressures in the escaping scenario. The shape of the transit spectra is determined by |$\mathcal {N}_{\mathrm{Na}}$|⁠, but for better readability we show the auxiliary parameter P₀. We use the planetary equilibrium temperature (Table 1) and a solar mixing ratio of |$\chi _{\mathrm{Na}}=1.7\,$| ppm to calculate the reference pressures. We fix |$\bar{v}_{\mathrm{Na}}=10\,$| km s⁻¹ in this figure. The planetary parameters and data points are for WASP-49b, taken from Wyttenbach et al. (2017).

The speeds of the absorbing atoms in escaping winds greatly exceed the thermal speeds. As described in Section 2.1, we use the wind speed to calculate a line temperature, which is then treated as thermal Doppler broadening of the line. We simplify the velocity profile of the escaping wind to a constant value (as in the hydrostatic model with vertical winds of Seidel et al. 2020). We therefore employ the mean speed of the absorbing atoms in the escaping wind, |$\bar{v}_{\mathrm{i}}$|⁠, as second free parameter. The impact of this parameter on a transmission spectrum for a hot Jupiter with planetary parameters of WASP-49b is shown in Fig. 5. We note that line broadening in an escaping wind consists of both the broadening due to the stochastic, thermal motions and of the wavelength shifts according to the line-of-sight wind velocities (Seidel et al. 2020). Our simplified treatment of line broadening consists only of the latter effect (but treating the wind velocity as thermal speeds), since we expect this to be the dominant source of line broadening. Nevertheless, the inclusion of the broadening due to the inherent thermal motions of the escaping gas would alter the Doppler core of the Voigt profile and thereby also the transmission spectrum, which could become important for lower wind speeds.

Figure 5.

Simulation of high-resolution transit spectra in the sodium doublet for different |$\bar{v}_{\mathrm{Na}}$| in the escaping scenario. We set |$\mathcal {N}_{\mathrm{Na}}=2.16\times 10^{33}$|⁠, corresponding to a pressure at the base of the wind of |$P_0=0.1\,$|nbar (using the planetary equilibrium temperature from Table 1 and |$\chi _{\mathrm{Na}}=1.7\,$| ppm). The planetary parameters and data points are for WASP-49b, taken from Wyttenbach et al. (2017).

We estimate the mass-loss rate in this evaporative scenario using dishoom. Hydrodynamically escaping atmospheres have two escape regimes: radiation-recombination limited (Murray-Clay et al. 2009) and energy-limited escape (Watson, Donahue & Walker 1981), where the latter regime is found to be accurate within a factor of 1.1 for a |$1\, M_J$| planet and within a factor of ∼1/4 for a |$0.3\, M_J$| planet (Allan & Vidotto 2019). Here, we choose to crudely estimate the mass-loss rate of the escaping absorber using the energy-limited approximation to hydrodynamic escape (e.g. Johnson, Volkov & Erwin 2013; described as |$\dot{M}_U$| in Oza et al. 2019b):

$$\begin{eqnarray*} \dot{M}_{\mathrm{esc,i}} \sim x_{\mathrm{i}}m_{\mathrm{i}} \frac{Q}{U}, \end{eqnarray*}$$

(11)

where Q is the XUV-heating rate of the upper atmosphere given a heating efficiency between 0.1 and 0.4, U is the binding energy of the atmosphere, x_i is the mass mixing ratio for the atom with mass m_i. For sodium, we use x_Na = 1.7 × 10⁻⁵, which corresponds to the solar volumetric mixing ratio of sodium of 1.7 ppm, multiplied by m_Na/μ (with |$\mu =2.3\,$| amu). This mass-loss rate, together with the required source rate estimated from equation (6), comprises the coupling between dishoom and prometheusfor the escaping scenario. We note that the upper limits we use here are to examine the extremeties of the evaporative wind scenario. As the escaping atmosphere passes the exobase, the energy-limited approximation used here will overestimate the source rate (Johnson et al. 2013). Furthermore, as hot Jupiters are expected to have strong magnetic fields (Yadav & Thorngren 2017), the mass-loss rates may be far less as pointed out by Christie et al. (2016) based on the MHD simulations of Trammell, Arras & Li (2011); Trammell, Li & Arras (2014); Owen & Adams (2014) as well as charge-exchange and hydrodynamical simulations by Tremblin & Chiang (2013).

2.5 Exomoon scenario

As observed in the Jupiter-Io system, a tidally heated satellite can outgas significant amounts of particles, especially neutral sodium, via volcanism. For an exo-Io orbiting a hot Jupiter, the outgassing is further enhanced by sublimating the silicate surface (see equation 18). At Io the observed volcanic escape is not due to outgassing however, but sublimation and volcanic outgassing (see equation 5 in Oza et al. 2019b and surrounding text) coupled with plasma-driven escape to space. We shall therefore refer to the observed gas near the satellite as atmospheric sputtering (Haff, Watson & Yung 1981; Johnson 2004), although outgassing is often used interchangeably to describe the satellite itself. To first order, we approximate the number density profile of the absorber by scaling n_i(r) = n(r) × χ_i(r) to the sputtered Na i number density profile observed at Io (Burger et al. 2001) with a power-law exponent of q_moon = 3.34:

$$\begin{eqnarray*} n_{\mathrm{moon,i}}(r)=n_{0,\mathrm{i}} \left(\frac{R_{\rm s}}{r}\right)^{q_{\mathrm{moon}}}, \end{eqnarray*}$$

(12)

where R_s, the radius of the satellite, is set to Io’s radius. Since we directly calculate the number density of the absorber in this scenario, we can drop the mixing ratio profile χ_i(r) for the computation of the optical depth along the chord in equation (1). As in the escaping scenario, we use the total number of the absorbing atoms in the system as a free parameter. We compute this quantity analogously to equation (9) in the escaping scenario:

$$\begin{eqnarray*} \mathcal {N_{\mathrm{i}}}=\iiint _V \, n_{\mathrm{i}}(r) \mathrm{d}V=\left(\frac{4\pi }{q_{\mathrm{moon}}-3}\right) n_{0,\mathrm{i}}R_{\rm s}^3 . \end{eqnarray*}$$

(13)

Note that the gas in such a sputtered cloud is not in local thermodynamic equilibrium; therefore, we do not use pressure as an auxiliary parameter. Instead, we convert |$\mathcal {N_{\mathrm{i}}}$| into a source rate using equation (6). Transit spectra for different values of this mass-loss rate are shown in Fig. 6. As we want to isolate the different scenarios, we neglect the planetary atmosphere in this scenario and consider only the sputtered cloud from the satellite for the computation of the transit spectrum. For computational convenience, we place the satellite in the centre of our coordinate system (Fig. 1) to exploit the spherical symmetry of the system. The reference radius in this scenario corresponds to the surface of the satellite, R₀ = R_s = R_Io.

Figure 6.

Simulation of high-resolution transit spectra in the sodium doublet for different mass-loss rates of the sodium atoms in the exomoon scenario. We set |$\bar{v}_{\mathrm{Na}}=10\,$| km s⁻¹. The planetary parameters and data points are for WASP-49b, taken from Wyttenbach et al. (2017).

We will again refer to the speeds of atoms instead of temperatures in this scenario as T is not well-defined, moreover irrelevant for an exosphere/collisionless gas. Similar to the escaping scenario, we describe the particles with a line temperature calculated from the mean velocity of the absorbing atoms (which we assume to be constant throughout the system). Therefore, we have the mean velocity of the absorbing atoms in the system, |$\bar{v}_{\mathrm{i}}$|⁠, as our second free parameter. Transit spectra for different values of this velocity are shown in Fig. 7.

Figure 7.

Simulation of high-resolution transit spectra in the sodium doublet for different mean velocities of the sodium atoms (corresponding to different effective temperatures) in the exomoon scenario. We set |$\dot{M}_{\mathrm{Na}}=4\times 10^8\,$| g s⁻¹, which corresponds to |$\mathcal {N}_{\mathrm{Na}}=2.5\times 10^{33}$| neutral sodium atoms in the system. The planetary parameters and data points are for WASP-49b, taken from Wyttenbach et al. (2017).

We calculate the plasma-driven mass-loss rate of the absorber within dishoom using the scaling parameters |$\mathcal {P} = \frac{P_s}{P_{\mathrm{Io}}}$| = f_mag + f_ram, |$\mathcal {V} = \frac{v_s}{v_{\mathrm{Io, torus}}}$|⁠, |$\mathcal {U} = \frac{U_s}{U_{\mathrm{Io}}}$|⁠, |$\mathcal {R}_x = \frac{\mathcal {R}_{x,s}}{\mathcal {R}_{x,\mathrm{Io}}}$| describing the total plasma pressure, relative ion velocities, binding energies, and exobase altitudes, respectively. We can therefore express |$\dot{M}_P$| (equation 7 of Oza et al. 2019b) more tractably in terms of a scaled plasma-heating rate |$\propto \mathcal {P} \mathcal {V} \mathcal {R}_x^2$| as

$$\begin{eqnarray*} \begin{aligned} \dot{M}_{\mathrm{moon,i}} &\sim x_{\mathrm{i}} \frac{\mathcal {P} \times \mathcal {V}}{\mathcal {U}} \mathcal {R}_x^2 \dot{M}_{\mathrm{Io}} \\ &\sim \frac{x_{\mathrm{i}} }{\mathcal {U}} \left(\frac{B_r^2}{2 \mu _0} P_{\mathrm{mag,Io}}^{-1} + n_{\mathrm{i}} m_{\mathrm{i}} u_{\mathrm{i}}^2 P_{\mathrm{ram,Io}}^{-1} \right) \mathcal {V} \mathcal {R}_x^2 \dot{M}_{\mathrm{Io}}. \end{aligned} \end{eqnarray*}$$

(14)

The dominant components of the total plasma pressure are the magnetic pressure, where |$P_{\mathrm{mag,Io}} = 1.5 \times 10^{-5}\, \mathrm{dynes\, cm}^{-2}$|⁠, dependent on the magnetic field strength at the satellite radius (B_r) given μ₀ the permeability of free space . The ram pressure is also critical, where |$P_{\mathrm{ram,Io}} = 2.6 \times 10^{-6}\, \mathrm{dynes\, cm}^{-2}$|⁠, dependent on n_i, m_i, and u_i the ion number density, mass, and velocity, respectively. We adopt a mass mixing ratio with respect to SO₂, x_Na = 0.1, corresponding to the that used for an exo-Io; for a full description on extrasolar atmospheric sputtering, see section 4.2.1 in Oza et al. (2019b). As dishoom does not explicitly track the ions that drive escape we use the velocity distributions modelled explicitly for Io by Smyth & Combi (1988); Smyth (1992) to constrain the velocities of the sodium gas sputtered from our satellite. The velocities range between 2 and 30 km s⁻¹, whereas speeds approaching 100 km s⁻¹ are possible due to charge exchange. A summary of the plasma processes and mean velocities associated with various processes are presented in Tables 2 and 3, respectively.

2.6 Torus scenario

As observed at Saturn, a moon or debris around a planet can source a circumplanetary torus with neutrals (in our study Na & K) in two ways. Case (1) Direct outgassing from an active satellite: Enceladus outgasses water-products generating an OH torus along its orbital path (Johnson et al. 2006b). Case (2) Desorption of grains from a toroidal atmosphere Johnson & Huggins (2006): UV photoionization of ice grains generates a close-in O₂ torus at Saturn (Johnson et al. 2006a). We can model both sources of an exoplanet torus via the number density of the absorber as

$$\begin{eqnarray*} n_{\mathrm{tor,i}}(a,z)=n_{0,\mathrm{i}}\times \exp \left [- \left(\frac{z}{H_t} \right)^2\right ]\times \exp \left [-\left(\frac{a-a_t}{4H_t}\right)^2\right], \nonumber\\ \end{eqnarray*}$$

(15)

which depends on the torus scale height H_t and the satellite semimajor axis a_t, which approximates the distance between the planet and the circumplanetary torus. a denotes the radial distance in the orbital x − y-plane (see Fig. 1), |$a=\sqrt{x^2+y^2}$|⁠. We set a_t = 2R₀ for our analysis as in Oza et al. (2019b) based on Domingos, Winter & Yokoyama (2006) and Cassidy et al. (2009), although in principle this number is allowed to decrease until the Roche radius for a silicate grain (Roche 1849). The torus scale height can be expressed as H_t = a_t × v_ej/v_o, where |$v_{\mathrm{o}}=\sqrt{GM_p/a_t}$| is the orbital velocity of the debris or venting moon and v_ej denotes the ejection velocity of the atoms, which we fix to 2 km s⁻¹ based on Johnson & Huggins (2006). As in the exomoon scenario, we can drop the mixing ratio profile of the absorber χ_i(r) for the computation as it is already incorporated into the number density profile.

The geometry for the computation of the transit spectra is more complex in this scenario due to the number density profile being azimuthally symmetric instead of spherically symmetric. Hence, we also need to (linearly) discretize our coordinate grid along the y-axis. Equation (1) then needs to be adjusted to the new geometry:

$$\begin{eqnarray*} \tau (y,z,\lambda)=2 \int _0^{\infty }\! n_{\mathrm{tor,i}}\Bigl (a,z\Bigr) \sigma \bigl (\lambda ,T(a,z)\bigr) \mathrm{d}x\, . \end{eqnarray*}$$

(16)

To obtain the flux decrease at a certain wavelength, the optical depth has to be averaged both over the y- and z-coordinate. Equation (3) changes to

$$\begin{eqnarray*} \Re =\frac{4}{\pi (R_{\ast }^2-R_0^2)}\int _0^{R_{\ast }}\!\mathrm{d}z\int _{\operatorname{Re}\bigl (\sqrt{R_0^2-z^2}\bigr)}^{\sqrt{R_{\ast }^2-z^2}}\!\mathrm{d}y\, \Bigl [1-{\rm e}^{-\tau (y,z,\lambda)}\Bigr ]. \end{eqnarray*}$$

(17)

As in the exomoon scenario, we have |$\bar{v}_{\mathrm{i}}$| and |$\mathcal {N_{\mathrm{i}}}$| as our free parameters. Transit spectra for some choices of these parameters are shown in Figs 8 and 9. There is again the following conversion between |$\mathcal {N_{\mathrm{i}}}$| and n_{0, i}, which has to be solved numerically: |$\mathcal {N_{\mathrm{i}}}=\iiint _V \!n_{\mathrm{tor,i}}(a,z)\, \mathrm{d}V$|⁠. Since the circumplanetary torus is not in thermodynamical equilibrium, we again use the absorber mass-loss rate |$\dot{M}_{\mathrm{i}}$| (instead of pressures) as an auxiliary parameter for an easier interpretation of |$\mathcal {N_{\mathrm{i}}}$| and to couple the radiative transfer code to calculations within dishoom.

Figure 8.

Simulation of high-resolution transit spectra in the sodium doublet for different mass-loss rates of the sodium atoms in the torus scenario. We set |$\bar{v}_{\mathrm{Na}}=10\,$| km s⁻¹. The planetary parameters and data points are for WASP-49b, taken from Wyttenbach et al. (2017).

Figure 9.

Simulation of high-resolution transit spectra in the sodium doublet for different mean velocities of the sodium atoms in the torus scenario. We set |$\dot{M}_{\mathrm{Na}}=2\times 10^9\,$| g s⁻¹, which corresponds to |$\mathcal {N}_{\mathrm{Na}}=1.3\times 10^{34}$| neutral sodium atoms in the system. The planetary parameters and data points are for WASP-49b, taken from Wyttenbach et al. (2017).

Case 1: Outgassed Torus: While the dominant escape mechanism that sources a cloud uniquely due to a satellite (exomoon scenario) is atmospheric sputtering, a plasma torus can be fuelled by several of the plasma processes described in Tables 2 and 3. We first focus on a directly outgassed torus: a tenuous Na i torus fuelled by a small outgassing satellite similar to Enceladus (e.g. Porco et al. 2006; Johnson et al. 2006b; Postberg et al. 2009) where the probability of escape, |$\mathcal {R} (1+ \lambda _0) \exp (-\lambda _0) \approx 1$|⁠. Here, λ₀ is the Jeans parameter, and |$\mathcal {R}$| a parameter to account for surface heating (see Johnson et al. 2015). This case of an exo-Enceladus fuelling a torus is therefore in contrast with the exo-Ios described in Oza et al. (2019b) in that the thermal desorption |$\dot{M}_{0, d}$| (see equations 10– 12 in Oza et al. 2019b and Table 5) approaches the thermal evaporation rate of silicate grains so that

$$\begin{eqnarray*} \begin{aligned} \dot{M}_{\mathrm{tor1,i}} &\sim \dot{M}_{0, \mathrm{evap}} \\ &\sim x_{\mathrm{i}} 4 \pi R_{\rm s}^2 P_{\mathrm{vap}} (T_0) \left(\frac{m_{\mathrm{i}}}{2 \pi k_{\rm B} T_0} \right)^{1/2}, \end{aligned} \end{eqnarray*}$$

(18)

The toroidal mass rate for Case 1 is then limited by the grain temperature T₀ and the surface area of the evaporating body, |$4\pi R_s^2$|⁠. A fundamental parameter driving the evaporation rate is the vapour pressure of the rocky mineral in question, P_vap. In principle, this formalism can compute the sublimation of arbitrary grain compositions, given that experiments have constrained the necessary coefficients to estimate the vapour pressure (c.f. Table 3, equation 13; van Lieshout, Min & Dominik 2014, where we use MgSiO₃ and Fe₂SiO₄). x_i is the mass fraction of the absorbing atom outgassing off of grains of exo-Io composition (Oza et al. 2019b), and m_i the mass. As little is known regarding the volatile composition of silicate grains at hot Jupiters, we use a conservative estimate of x_Na = 0.003 assuming the grain has experienced substantial volatile loss, like at Io. The value is consistent with direct observations of NaCl by Lellouch et al. (2003), where outgassed values of NaCl/SO₂ = 0.3–1.3  per cent are consistent with CI chondritic compositions of Na/S = 0.13 and Cl/S = 0.01. The value is also roughly consistent with grains ejected from Enceladus (NaCl/H₂O = 0.005–0.02) observed by the Cosmic Dust Analyser aboard Cassini (Postberg et al. 2009). We comment that due to the magmatic nature of these grains, more recent geophysical modelling of high-temperature rocky bodies is warranted to better predict x_Na (e.g. Noack, Rivoldini & Van Hoolst 2017; Bower et al. 2019; Sossi et al. 2019) and its origin.

Table 5.

Sodium source rates computed within dishoom. The enhanced sodium source rate |$\dot{M}_{\mathrm{tor2,Na}}$| corresponds to a desorbing torus (equation 19). The minimally required source rate |$\dot{M}_{\mathrm{min,Na}}$| is calculated using equation (29) with the observed equivalent widths at the Na D2 line (Wyttenbach et al. 2015; Wyttenbach et al. 2017) for the respective planet.

Source rate	WASP-49b	HD189733b
|$\dot{M}_{\mathrm{esc,Na}}$| (kg s−1)	102.4 ± 0.3	102.3 ± 0.3
|$\dot{M}_{\mathrm{moon,Na}}$| (kg s−1)	104.3 ± 1.5	104 ± 1.2
|$\dot{M}_{\mathrm{tor1,Na}}$| (kg s−1)	103.1 ± 1.3	101.0 ± 1.5
|$\dot{M}_{\mathrm{tor2,Na}}$| (kg s−1)	106.4 ± 1.3	104.1 ± 1.5
|$\dot{M}_{\mathrm{min,Na}}$| (kg s−1)	104.7 ± 0.5	103.9 ± 0.1

Case 2: Desorbing Torus: In the case of a directly outgassed torus, the surface area is quite small yielding a lower limit to the toroidal mass rate. However for a toroidal atmosphere, the desorption of Na from a silicate torus (similar to Saturn’s O₂ ring atmosphere from H₂O grains, Johnson et al. 2006a) can considerably enhance the toroidal mass rate. The enhancement factor ξ is simply the ratio of the surface areas:

$$\begin{eqnarray*} \begin{aligned} \dot{M}_{\mathrm{tor2,i}}&=\xi \times \dot{M}_{\mathrm{tor1,i}} \sim \frac{4 \pi ^2 a_t H_t}{4 \pi R_s^2} \dot{M}_{\mathrm{tor1,i}} \\ &=4\pi \left(\frac{R_p}{R_s} \right)^2 \left(\frac{v_{\mathrm{ej}}}{v_{\mathrm{o}}} \right) \dot{M}_{\mathrm{tor1,i}}. \end{aligned} \end{eqnarray*}$$

(19)

For a range of semimajor axes a_t we find this enhancement ranges from ξ ∼ 10^3 ± 0.5. In this way, distinguishing between a directly outgassed toroidal exosphere and a desorbing toroidal atmosphere may be achieved by high-resolution data sets (see Section 6.3.1).

2.7 Comparison to observations

We compare simulated spectra for all four scenarios to high-resolution observations of the sodium D doublet at WASP-49b and HD189733b. We restrict our retrieval analysis to the data points within the wavelength interval |${}[5888,5898]\, \mathring{\rm A}$|⁠. We note that the Na i detections at these two planets were reduced and analysed by the same algorithm, using the same instrument (HARPS) thereby validating a 1-to-1 comparison. We retrieve the two free parameters in every scenario using the reduced chi-squared statistic (e.g. Ocvirk et al. 2006; Andrae, Schulze-Hartung & Melchior 2010):

$$\begin{eqnarray*} \chi ^2_r=\frac{\chi ^2}{\nu }=\frac{1}{\nu }\sum _i\Bigl (\frac{O_i-C_i}{\sigma _i}\Bigr)^2, \end{eqnarray*}$$

(20)

where O_i are observations with the corresponding errors σ_i, C_i are computed values, and ν are the degrees of freedom of the model given by ν = N_O − N_P, the difference between the number of data points and the number of free parameters (in our analysis, N_P = 2 for all scenarios). The standard deviation of this distribution is given by |$\sigma =\sqrt{2/\nu }$|⁠. Lower values of |$\chi ^2_r$| indicate that the model is a better fit to data. If |$\chi ^2_r\lt 1$|⁠, the model overfits the data.

To compare our simulated transit spectra to the observations, we apply a convolution with the instrumental line spread function (LSF), normalization, and binning routine to the raw simulated spectra. For the convolution of the raw spectrum with the LSF of the instrument, we use a Gaussian with full width at half-maximum of 0.048 |$\mathring{\rm A}$|⁠. We then bin the convolved spectrum to 0.2 |$\mathring{\rm A}$| wide bins centred on the D2 and D1 absorption lines, which maximizes the signal-to-noise ratio according to Wyttenbach et al. (2017). Finally, we normalize the spectrum by the average transit spectrum ℜ in two reference bands, |$B=[5874.94;5886.94]\, \mathring{\rm A}$| and |$R=[5898.94;5910.94]\, \mathring{\rm A}$|⁠.

3 OPTICALLY THIN GAS IN HIGH-RESOLUTION TRANSMISSION SPECTRA OF EXOPLANETS

Consider a foreground gas of atomic species i, illuminated by a background radiation field F_{out, λ}. If we treat the gas as a slab the observed intensity can be written according to Beer’s law or Lambert’s law (e.g. Rybicki & Lightman 1979):

$$\begin{eqnarray*} F_{{\rm in}, \lambda } = F_{{\rm out}, \lambda } {\rm e}^{-\tau (\lambda)}, \end{eqnarray*}$$

(21)

appropriate for a single line of sight. Using equation (1) and assuming that the absorption cross-section is solely a function of wavelength (which holds for our scenarios as we do not vary the Doppler broadening parameter and neglect pressure broadening), the slant optical depth for the slab can be written as

$$\begin{eqnarray*} \tau (\lambda) = N_{\mathrm{i}} \sigma (\lambda), \end{eqnarray*}$$

(22)

with the line-of-sight column density of species i:

$$\begin{eqnarray*} N_{\mathrm{i}}=\int _{-\infty }^{+\infty }\!n_{\mathrm{i}}(x)\, \mathrm{d}x. \end{eqnarray*}$$

(23)

For transmission spectroscopy geometry, equation (21) needs to be averaged over infinitely many line of sights. In the context of a hydrostatic atmosphere, this can be done analytically. We briefly review this formalism in the next section.

3.1 Canonical hydrostatic gas: an effective column density at τ ∼ 1

For a hydrostatic planetary atmosphere with a constant pressure scale height H, the line-of-sight column is given by (Fortney 2005)

$$\begin{eqnarray*} N_{\mathrm{i}}(z) = N_{v,\mathrm{i}}(z)\sqrt{\frac{2 \pi R_0}{ H}}= n_{\mathrm{i}}(z)\sqrt{2 \pi R_0 H}, \end{eqnarray*}$$

(24)

where N_{v, i}(z) = n_i(z)H denotes the vertical column density above altitude z of species i and |$\sqrt{(2 \pi R_0)/ H}$| is the Chapman enhancement factor. We have assumed that the mixing ratio of atomic species i does not change throughout the line of sight. Combining equations (22) and (24), one can define a reference optical depth:

$$\begin{eqnarray*} \tau (R_0,\lambda)\equiv \tau _0(\lambda) = n_{0,\mathrm{i}}\sigma (\lambda) \sqrt{2 \pi R_0H}, \end{eqnarray*}$$

(25)

which can also be written in terms of a reference pressure. Continuing with a hydrostatic profile decaying over an altitude z:

$$\begin{eqnarray*} \tau (z,\lambda) = \tau _0(\lambda) {\rm e}^{- z/H}. \end{eqnarray*}$$

(26)

Equation (26) is an analytical expression for the general optical depth profile (equation 1) needed for the calculation of transmission spectra. Integrating over all lines of sight (equation 3) using an identity (Chandrasekhar 1960), a closed form expression for the atmospheric transit radius R_λ can be obtained (de Wit & Seager 2013; Bétrémieux & Swain 2017; Heng & Kitzmann 2017; Jordán & Espinoza 2018) given reasonable assumptions for an isothermal atmosphere and τ₀(λ) → ∞. The transit radius and transit depth are related via (c.f. equation 3)

$$\begin{eqnarray*} 1 - \Re (\lambda) = \frac{R_{\lambda }^2 - R_{0}^2 }{R_{*}^2}. \end{eqnarray*}$$

(27)

Remarkably, by assuming an optically thick reference pressure (τ₀(λ) → ∞), it turns out that the optical depth at the transit radius τ(R_λ, λ) ≡ τ_eff(λ) (termed effective optical depth) is ≈0.56, independent of wavelength. This analytical result is elegant in that the notion of an atmospheric transit radius can be interpreted as an approximate boundary between opaque and transparent layers. In this sense, it is the (wavelength-dependent) location z at which τ(λ, z) ≈ 0.56, which determines the transit radius and simultaneously the transit depth. This formalism implies that the effective atmospheric column density, at the atmospheric transit radius R_λ, converges for all wavelengths as |$N (R_{\lambda }) \approx \frac{0.56}{\sigma (\lambda)}$|⁠. At T = 10³−10⁵ K, the atmospheric column density probed at Na D2 line centre is |$N(R_{\mathrm{NaD2}}) \sim 5 \times 10^{11} {\text {--}} 1.5 \times 10^{12}\, \mathrm{cm}^{-2}$|⁠.

3.2 Non-hydrostatic Gas: evaporative column densities at τ ≪ 1

In solving the same problem for a collisionless exosphere of arbitrary number density profile n(r), one cannot pinpoint the location of such a boundary, rendering the concept of an atmospheric transit radius inappropriate. On the other hand, as τ ≪ 1, the transit depth reveals the total number of absorbing atoms|$\mathcal {N}_{\mathrm{i}}$| spread over the stellar disc. The remarkable ability to constrain the total number of absorbing atoms |$\mathcal {N_{\mathrm{i}}}$| is a fundamental property of foreground gas well known from studies of the interstellar medium (Spitzer 1978; Draine 2011) using measurements of equivalent widths W_λ (in units of |$\mathring{\rm A}$|⁠):

$$\begin{eqnarray*} W_{\lambda } = \int _{\lambda - \Delta \lambda }^{\lambda + \Delta \lambda } \!\Bigl (1 - \Re (\lambda) \Bigr)\, \mathrm{d} \lambda . \end{eqnarray*}$$

(28)

We can then give a lower bound to the column density of the absorbing species i (equation 23) using the approximation of the curve of growth in an optically thin regime (Draine 2011; 〈N〉_obs in Oza et al. 2019b):

$$\begin{eqnarray*} N _{\mathrm{min,i}}= 1.13\times 10^{-12}\, \mathrm{cm}^{-1} \frac{W_{\lambda }}{f_{ik} \lambda _{ik}^2}, \end{eqnarray*}$$

(29)

where f_ik is the oscillator strength of the line and λ_ik the wavelength of the transition. Unlike the interstellar medium where a single line-of-sight is considered, the background radiation field (star) is far larger than the foreground gas (the planetary atmosphere/exosphere). Therefore, a transmission spectrum requires averaging an arbitrary number density distribution n(r) over infinitely many line-of-sights (or chords) since the gas is not necessarily homogeneous.⁵ In order to apply equation (29), we reduce the various line-of-sight column densities in an exosphere to a disc-averaged column:

$$\begin{eqnarray*} \langle N_{\mathrm{i}} \rangle \equiv \frac{\mathcal {N}_{\mathrm{i}}}{\pi R_*^2}. \end{eqnarray*}$$

(30)

Provided that an outgassing or evaporative source is present, |$\mathcal {N}_{\mathrm{i}}$|⁠, the number of evaporating atoms can be described by an arbitrary mass-loss rate (e.g. equations 11, 14, 18). The mass-loss rates then directly supply the above disc-averaged column density. In other words, an evaporative column density 〈N_i〉 is equivalent to the observed column density N_{min, i} (Johnson & Huggins 2006; Oza et al. 2019b), given that the absorption occurs in an optically thin regime. As we show in the remainder of this paper, the evaporative scenarios do indeed absorb in a primarily optically thin regime. In Table 5, we show calculations from dishoom demonstrating that predicted mass-loss rates for hot Jupiters roughly align with the required equivalent widths (converted into minimal mass-loss rates using equations 6, 29, and 30) found by high-resolution spectroscopy.

3.3 Evaporative curve of growth for sodium at an exoplanet

Using our radiative transfer code, we calculate equivalent widths at Na D2 line centre for our three evaporative scenarios (Fig. 10). The equivalent width depends on the free parameters |$\bar{v}_{\mathrm{Na}}$| and |$\mathcal {N}_{\mathrm{Na}}$| (equivalently: |$\mathcal {N}_{\mathrm{Na}}$|⁠, 〈N_Na〉 or |$\dot{M}_{\mathrm{Na}}$|⁠). The dependence of the equivalent width on a (disc-averaged) column density is reminiscent of the classical curve of growth for a single line of sight. If the equivalent width is known to high precision, one can constrain an evaporative column density of occulting atoms depending on the scenario in Fig. 10. These values coincide with the values from Oza et al. (2019b, their table 5) and indicate how effective each evaporative scenario is in generating the observed absorption.

Figure 10.

Equivalent width at the Na D2 line versus evaporating atoms during transit, for planetary parameters of HD189733b (Table 1). The evaporative column density is essentially a disc-averaged column (equation 30). We also show the classical curve of growth (‘homogeneous’), which is the equivalent width of a single line of sight as a function of column density. We convert 〈N_i〉 into sodium source rates using equations (6), (29), and (30). Bold lines are with |$\bar{v}_{\mathrm{Na}}=10\,$| km s⁻¹, dashed lines have |$\bar{v}_{\mathrm{Na}}=1\,$| km s⁻¹.

To generate typical observed equivalent widths of |$W_{\mathrm{NaD2}}\sim 1{\!-\!}10\, \mathrm{m}\mathring{\rm A}$| at hot Jupiters, the required evaporative column densities range from ∼10¹⁰ to |$\sim 10^{12}\, \mathrm{Na\, cm}^{-2}$|⁠. For a torus, whose number density is a Gaussian distribution (equation 15), the column densities can be far larger. We note that absorption along a single line of sight (or, equivalently, in a homogeneous cloud as in Hoeijmakers et al. 2020) is always absorbing more efficiently than our evaporative sodium distributions.

Independent of the source driving mass-loss, we confirm that extreme mass-loss rates are required for observation of extrasolar evaporating sodium. To observe an equivalent width of 10 m|$\mathring{\rm A}$|⁠, for instance, roughly ∼10⁵ kg s⁻¹ of sodium gas at 10 km s⁻¹ is required for an escaping atmosphere or exo-Io as estimated by Oza et al. (2019b). In comparison, Io’s volcanism outgasses ∼7 × 10⁶ kg s⁻¹ of SO₂ (Lellouch et al. 2015). A thermally desorbing torus, however, requires |$\dot{M}_{\mathrm{Na}}\sim 5 \times 10^{5}$| kg s⁻¹ based on the scenario described. We note, based on the discussion in Sections 2.6 (equation 19) and 6.3, that these rates can be achieved if the grains are thermally desorbing over a large surface area, while being trapped in a toroidal magnetic field for instance. This analysis provides an evaporative curve of growth for atomic sodium viewed at an exoplanet. The quantity of evaporating gas ∼10³¹−10³⁴ Na atoms is able to govern the detection and non-detection of optically thin absorption during transit. The evaporative curve of growth is drastically influenced by the spatial distribution of the Na atoms across the star during transit.

3.4 The D2-to-D1 line ratio

At the Jupiter-Io system, the D2/D1 ratio is able to provide information on the velocity of the Na i atoms, which is observed to be variable over decades of observations. Therefore in our application to extrasolar systems, the ratio may be able to inform predictions on the average velocity distributions of the Na atoms and their ongoing behaviour. However, the ratio first and foremost provides a very simple result.

Fortuitously, the oscillator strength of the Na D2 transition is twice as large as the D1 transition, meaning that the absorption cross-section and optical depth at D2 line centre is double the respective values at D1 line centre (Draine 2011). In this way, the optical depth at sodium D2 line centre is τ(λ_NaD2) ≡ τ_NaD2 = τ(λ_NaD1)/2. This reveals an easy diagnostic:

$$\begin{eqnarray*} f_{\mathrm{D2/D1}}=\frac{1-{\rm e}^{-\tau _{\mathrm{NaD2}}}}{1-{\rm e}^{-\tau _{\mathrm{NaD2}}/2}}. \end{eqnarray*}$$

(31)

This relation, together with the measured line ratios for HD189733b and WASP-49b, is shown in Fig. 11. Throughout this paper, we calculate line ratios by using the transit depths averaged over bandwidths of |$0.2\, \mathring{\rm A}$|⁠, centred on the absorption lines. We note that the exact value of the line ratio depends on the choice of the bandwidth, a negligible effect for modeled transmission spectra but not for the observations. For small bandwidths, the measurement error associated with the binned transit depth is large, while for large bandwidths, the observations at wavelengths which are more than |$\sim 0.5\, \mathring{\rm A}$| away from the line centre mostly contain noise. Hence, we choose an intermediate bandwidth of |$0.2\, \mathring{\rm A}$| (which maximizes the signal-to-noise ratio according to Wyttenbach et al. 2017) for the calculation of the line ratios. While optically thick chords with τ_NaD2 > 10 lead to a line ratio of one, f_D2/D1 transitions over two orders of magnitude in τ_NaD2 to f_D2/D1 = 2 for optically thin chords, τ_NaD2 < 0.1.

Figure 11.

Computation of the D2/D1-ratio f_D2/D1 as function of the optical depth at sodium D2 line centre (equation 31). In the optically thin regime this ratio approaches two, while in the optically thick regime f_D2/D1 goes to one. For the conversion of τ_NaD2 into a column, we assume a fixed absorption cross-section of σ_NaD2 = 4.62 × 10⁻¹² corresponding to a gas at T = 2000 K (equation 6.39 in Draine 2011).

Unfortunately, for the setting of transmission spectroscopy, this relation is not applicable in a straightforward way since one observes a flux decrease averaged over infinitely many chords. For example, a line ratio of 1.5 can be achieved in different ways: By having a constant line-of-sight column density throughout the stellar disc with τ_NaD2 ≈ 1, or by having 10 per cent of the area of the stellar disc blocked with τ_NaD2 ≈ 10 and the remaining 90 per cent having τ_NaD2 ≈ 0.1. Both of these models (and infinitely many other spatial distributions of the absorber) would lead to f_D2/D1 ≈ 1.5. For f_D2/D1 = 1.5, given that the column density is a smooth function of impact parameter, we can therefore only state that the majority of the absorption occurs along chords with τ_NaD2 roughly around 0.5, but chords with vastly different values of τ_NaD2 are probably also present in the system. Hence, the D2-to-D1 ratio tells us in which regime (optically thin/thick) the majority of the absorption occurs.

In the following, we shall test how the predicted quantities of evaporating gas fare against canonical hydrostatic assumptions for high-resolution Na i spectra at hot Jupiters.

4 FORWARD MODELLING AND SCENARIO COMPARISON

We elucidate the general features of evaporative transmission spectra in this section by constructing a forward model for the hot Jupiter HD189733b, using dishoom to calculate sodium source rates for the evaporative scenarios (listed in Table 5). These rates are converted into |$\mathcal {N}_{\mathrm{Na}}$| (using equation 6), a parameter which is then fed into prometheus to compute the transmission spectra. We emphasize that these values of |$\mathcal {N}_{\mathrm{Na}}$| should be regarded as lower limits, since we use a lower limit to the lifetime of neutral sodium in equation (6). We uniformly set |$\bar{v}_{\mathrm{Na}}=10\,$| km s⁻¹ for the three evaporative scenarios in this forward model, a simplification that does not affect the following analysis as the dominant effect of |$\bar{v}_{\mathrm{Na}}$| is only the broadening of the lines. We also examine two hydrostatic scenarios: our first hydrostatic model has |$T=T_{\mathrm{eq}}=1^{\prime }140\,$| K and |$P_0=1\,$| bar, the second model⁶ has |$T=10 T_{\mathrm{eq}}=11^{\prime }400\,$| K and |$P_0=0.1\, \mu$|bar (for the conversion of the reference pressures into P_0,Na we assume |$\chi _{\mathrm{Na}}=1.7\,$| ppm for both scenarios). This choice of parameters for the hydrostatic scenario is arbitrary, and intentionally spans a large range within the free parameters. The forward-model transmission spectra are shown in Fig. 12.

Figure 12.

Forward model transmission spectra of each scenario studied, for the hot Jupiter HD189733b. The dashed blue line represents an extremely heated hydrostatic atmosphere with |$T=10\, T_{\mathrm{eq}}$|⁠. The dashed red line has an enhanced sodium source rate corresponding to a desorbing torus (equation 19). Planetary parameters and data points are taken from Wyttenbach et al. (2015).

The four evaporative transmission spectra in Fig. 12 share some features such as negligible absorption between the lines and a D2-to-D1 line ratio significantly larger than one. The transit depth on the line cores is much larger for the exomoon scenario and for the desorbing torus (equation 19) compared to the two other evaporative transmission spectra, stemming from the almost three orders of magnitude larger sodium source rates (Table 5). The two hydrostatic scenarios have, in contrast to the evaporative scenarios, a D2-to-D1 line ratio only slightly larger than one. Furthermore, the hydrostatic scenario with T = T_eq exhibits significant absorption between the line cores due to the reference pressure of |$P_0=1\,$| bar leading to a very thick atmosphere.

The origin of the variety of transit spectra becomes apparent when examining the spatial structures of the different forward models. We show the corresponding sodium number density profiles in Fig. 13. While the sodium number density in the hydrostatic scenarios drops fast (especially in the lower-temperature case), the three evaporative scenarios lead to more extended and tenuous exospheres. The sodium number density of the torus scenario peaks at r = 2 × R₀ due to the fixed orbital separation of the torus of a_t = 2 × R₀ (see Section 2.6). These number density profiles can be converted into optical depth profiles at a fixed wavelength using equation (1). We fix the wavelength to the Na D2 line centre and show |$\tau _{\mathrm{Na\, D2}}(r)$| in Fig. 14. Since the optical depth profiles⁷ are equivalent to the line-of-sight integrals of the sodium number density profiles times a constant depending on the Doppler broadening (equations 22 and 23), the curves are similar in Figs 13 and 14, but decay slower in the latter plot (note that both plots cover approximately sixteen orders of magnitude). From the optical depth profiles, we can derive the properties of the different transmission spectra in Fig. 12.

• Hydrostatic, T = T_eq: The hydrostatic scenario with T = T_eq exhibits an optical depth that drops very fast with increasing altitude. This behaviour leads to a small and very optically thick layer, which absorbs all incoming starlight nearly uniformly up to a certain altitude, followed by negligible absorption. This leads to f_D2/D1 being very close to one in this model. We have a large reference pressure of |$P_0=1\,$|bar in this scenario, which leads to the slant optical depth at the reference radius |$\tau _{0,\mathrm{Na\, D2}}$| being of the order 10¹⁰ for the Na D2 line centre. Since the absorption cross-section does not drop by ten orders of magnitude between the line centres, we have significant absorption also between the lines.

• Hydrostatic, |$T=10\, T_{\mathrm{eq}}$|: At a larger temperature the hydrostatic scenario has a more puffed up atmosphere. However, the optical depth still drops very fast, leading to f_D2/D1 ≈ 1.1. Since the slant optical depth at the reference radius is only of the order 10³, the atmosphere does not produce significant absorption between the lines.

• Evaporative: On the other hand, the three evaporative scenarios have extended and mostly optically thin exospheres. Both the escaping and the exomoon scenario have long tails in their optical depth profiles. Since |$\tau _{\mathrm{Na\, D2}}$| is lower than unity in these two scenarios for the largest part of the exosphere, the resulting transmission spectrum is optically thin. For optically thin chords, the proportion of absorbed light is directly proportional to the optical depth, leading to f_D2/D1 ≈ 2. Since the optical depth profile in the escaping scenario is offset by two orders of magnitude due to the lower sodium source rate, the flux decrease in this scenario is accordingly lower. The torus scenario sourced by direct outgassing (equation 18) has a plateau with |$\tau _{\mathrm{Na\, D2}}$| slightly lower than one percent, extending over multiple radii, before the optical depth drops. Hence, the exosphere of the torus in this forward model does not generate significant absorption as seen in Fig. 12. On the other hand, the desorbing torus (equation 19) with a source rate three orders of magnitudes larger has its optical depth profile shifted up by the same factor, leading to significant absorption.

We conclude from this analysis that the three evaporative scenarios generally produce optically thin, extended exospheres (especially the escaping and exomoon scenarios), while hydrostatic models lead to a small, optically thick atmosphere. We find here that despite an atmosphere with a temperature ten times larger than the planetary equilibrium temperature, effectively enhancing the atmospheric scale height by a factor of ten, the hydrostatic atmosphere still drops very fast in number density and optical depth, such that the optically thin layer is negligibly small. This confirms our analysis in Section 3.2, in that the evaporative scenarios absorb mostly in an optically thin regime, unlike the hydrostatic scenario.

Figure 13.

Sodium number density profiles of each scenario studied. The number density profiles correspond to our forward models for the hot Jupiter HD189733b. The dashed blue line represents an extremely heated hydrostatic atmosphere with |$T=10\, T_{\mathrm{eq}}$|⁠, the dashed red line has an enhanced sodium source rate corresponding to a desorbing torus (equation 19). The exomoon profile starts in our model already at r = R_Io, we shift the curve to the planetary radius in this plot to enhance readability. For the torus scenario, we show the number density profiles through the orbital plane, n_{tor, Na}(r, 0).

Figure 14.

Optical depth (at Na D2 line centre) profiles of each scenario studied. The optical depth profiles correspond to our forward models for the hot Jupiter HD189733b. The dashed blue line represents an extremely heated hydrostatic atmosphere with |$T=10\, T_{\mathrm{eq}}$|⁠, the dashed red line has an enhanced sodium source rate corresponding to a desorbing torus (equation 19). The exomoon profile starts in our model already at r = R_Io, we shift the curve to the planetary radius in this plot to enhance readability. For the torus scenario, we show the optical depth profiles through the orbital plane, |$\tau (r,0,\lambda _{\mathrm{Na\, D2}})$|⁠. The line-of-sight column density is shown for a Doppler broadening parameter corresponding to T = 10T_eq, allowing a conversion from |$\tau _{\mathrm{Na\, D2}}(r)$| to N_Na(r) using equation (22).

5 OBSERVATIONAL ANALYSIS OF EVAPORATIVE SODIUM TRANSIT SPECTRA

In the following, we will perform a simple retrieval on two hot Jupiters (WASP-49b and HD189733b) applying the reduced χ²-statistic to determine best fits for all four scenarios. Using high-resolution observations of the sodium doublet, we determine the two free parameters in all scenarios. Next, we scrutinize the physical validity of the retrieved parameters based on our sodium source rate calculations from dishoom(Table 5), which comprises the coupling between our codes in this inverse modelling. The retrieval results for both planets are summarized in Table 6. Error bars for the retrieved parameters correspond to the 1σ-interval around the best-fitting model, where σ is the standard deviation of the |$\chi _r^2$|-statistic (see Section 2.7). We have |$\sigma =\sqrt{2/\nu }\approx 0.2$|⁠, where ν are the degrees of freedom of the model. We also evaluate f_D2/D1 for our best-fitting models. The error bars for the line ratios are calculated such that f_D2/D1 lies within the error bars with a likelihood⁸ of |$\approx 68{{\ \rm per\ cent}}$|⁠.

Table 6.

Retrieval results for WASP-49b and HD189733b. Note that due to the irregularity of the |$\chi _r^2$|-maps of WASP-49b (Fig. 16), the error bars for the retrieved parameters of this planet sometimes extend in only one direction. We also evaluate the D2-D1 line ratio f_D2/D1 for our models. From the measurements, we estimate f_{D2/D1, W49} = 1.28 ± 0.62 (Wyttenbach et al. 2017) and f_{D2/D1, HD189} = 1.74 ± 0.45 (Wyttenbach et al. 2015).

			WASP-49b			HD189733b
Scenario	Parameter	Retrieved	fD2/D1	|$\chi _r^2$|	Retrieved	fD2/D1	|$\chi _r^2$|
Hydrostatic	P0, Na (bar)	|$10^{-12.2^{+1.4}}$|	|$1.17^{+0.12}_{-0.02}$|	1.49	10−10.7 ± 0.5	|$1.086^{+0.002}_{-0.002}$|	1.49
	T (K)	6100−3400			5200 ± 1500
Escaping	|$\mathcal {N}_{\mathrm{Na}}$| (Na atoms)	033.1 ± 1	|$1.68^{+0.04}_{-0.06}$|	1.44	1032.5 ± 0.2	|$1.930^{+0.005}_{-0.005}$|	1.51
	|$\bar{v}_{\mathrm{Na}}$| (km s−1)	9+13			19 ± 9
Exomoon	|$\mathcal {N}_{\mathrm{Na}}$| (Na atoms)	1033.4 ± 0.6	|$1.71^{+0.01}_{-0.02}$|	1.45	1032.7 ± 0.2	|$1.737^{+0.002}_{-0.001}$|	1.50
	|$\bar{v}_{\mathrm{Na}}$| (km s−1)	9+13			18 ± 8
Torus	|$\mathcal {N}_{\mathrm{Na}}$| (Na atoms)	|$10^{33.5^{+2.3}_{-0.8}}$|	|$1.24^{+0.07}_{-0.04}$|	1.45	1032.6 ± 0.2	|$1.64^{+0.02}_{-0.04}$|	1.53
	|$\bar{v}_{\mathrm{Na}}$| (km s−1)	7+14			18 ± 9

5.1 Inverse modelling of WASP-49b

Wyttenbach et al. (2017) obtained a high-resolution spectrum of the hot Jupiter WASP-49b, observing significant absorption in the sodium doublet (more than two percent of the D2 line centre). The spectrum shows negligible absorption between the line cores and a D2-to-D1-line ratio which is f_D2/D1 = 1.28 ± 0.62 (calculated in bands of |$0.2\, \mathring{\rm A}$| centred on the sodium D lines). The observation and the best-fitting models for each scenario are shown in Fig. 15, we summarize the corresponding parameters in Table 6. Given the large uncertainty in f_D2/D1, all four scenarios have line ratios within the observational error bars. We remark that for the case of WASP-49b, this ratio converges to two for larger bandwidths since the D2 line is broader than the D1 line (Wyttenbach et al. 2017). If this line ratio was confirmed in a more precise measurement, the hydrostatic and the torus scenario would significantly underestimate f_D2/D1. The achieved goodness of fit is very similar in all four scenarios |$1.44\lt \chi _r^2\lt 1.49$|⁠, with slightly better fits for the evaporative scenarios. From the |$\chi _r^2$|-statistic, a model which correctly describes data with proper error bars should achieve |$\chi _r^2=1$|⁠. The fact that our best-fitting models have values of |$\chi _r^2$| which are ≳ 2σ away from |$\chi _r^2=1$| implies that our models either poorly describe the data or that the observational error bars are underestimated. One can rescale the |$\chi _r^2$|-values such that the best-fitting model achieves |$\chi _r^2=1$|⁠, which is essentially a rescaling of the error bars. However, this procedure is statistically incorrect (Andrae 2010), hence we refrain from such a rescaling and keep our measured values of |$\chi _r^2$|⁠. We remark that over the spectral range investigated here, our models seem to fit the data well when roughly checked by eye, which leads us to suggest that the oddly large values of |$\chi _r^2$| can primarily be attributed to an underestimation of the observational error bars.

Figure 15.

Transit spectrum for WASP-49b with our best-fitting models for each scenario. Planetary parameters and data points are taken from Wyttenbach et al. (2017). Note the break in the x-axis.

We observe in the |$\chi _r^2$|-maps for all four scenarios that the shapes of the best-fitting regions are irregular and extend very far in a particular direction (Fig. 16). Hence, the retrieved parameters in Table 6 can contain error bars in only one direction for the case of WASP-49b. In the hydrostatic scenario, the parameters are driven to values such that the resulting atmosphere is more optically thin. We retrieve a large temperature of |$T=6100_{-3400}\,$| K, significantly above the planetary equilibrium temperature of WASP-49b of 1400 K (Wyttenbach et al. 2017). We retrieve similar velocities and source rates in the three evaporative scenarios, with the most prominent difference being the line ratio in the torus scenario of |$f_{\mathrm{D2/D1}}=1.24^{+0.07}_{-0.04}$| (indicating absorption mostly in an optically thick regime). Although the torus scenario is evaporative, we see in Fig. 14 that the optical depth profile is constant over a large range of radii and in the transition region between optically thick and optically thin chords. Since our best-fitting torus scenario requires a sodium mass-loss rate comparable to the desorbing torus forward model of HD189733b, the optical depth profile of the best-fitting torus model for WASP-49b resembles the dashed red line in Fig. 14, indicating rather optically thick absorption.

Figure 16.

Parameter retrieval for WASP-49b. The colour-coding represents the deviation from |$\chi _r^2=1$|⁠, in standard deviations of the |$\chi _r^2$|-distribution (⁠|$\sigma =\sqrt{2/\nu }$|⁠, where ν corresponds to the degrees of freedom of the model, see Section 2.7). The contours are calculated on a 50 × 50 parameter grid and interpolated linearly using pyplot.contourf. The auxiliary parameter P₀ for the hydrostatic and escaping scenarios is computed using |$\chi _{\mathrm{Na}}=1.7\,$|ppm and (only for the escaping scenario) T = T_eq at the base of the wind.

The comparison of the retrieved mass-loss rates from prometheus with the calculated rates within dishoom is shown in Table 7. Since prometheus does not directly retrieve |$\dot{M}_{\mathrm{Na}}$| but rather |$\mathcal {N}_{\mathrm{Na}}$|⁠, we use equation (6) to obtain an upper limit to the mass-loss rate from the retrieved |$\mathcal {N}_{\mathrm{Na}}$|⁠. The retrieved source rate in the escaping scenario is nearly three orders of magnitude larger than the one we calculate using dishoom, indicating that an escaping wind can probably not provide enough sodium to generate the observed transit depth. The retrieved source rate in the exomoon scenario is also larger than the one we calculate from dishoom, but still within the error bars. For the torus scenarios, direct outgassing (equation 18) within dishoom significantly underestimates the retrieved source rate, while a desorbing torus (equation 19) overestimates it (but still within the error bars).

Table 7.

Mass-loss rate comparison between prometheus and dishoom. We emphasize that the prometheus mass-loss rates are upper limits due to our usage of equation (6) and a minimum lifetime of neutral sodium. The sodium source rates within the escaping scenarios fall short of the required rates. While the torus source rates due to direct outgassing (Torus 1, equation 18) are also significantly lower than the retrieved ones, we remark that the enhanced rates due to a desorbing torus (Torus 2, equation 19) are comparable to the retrieved values. Note that the retrieved mass-loss rates within prometheus are larger than the minimal mass-loss rates (Table 5).

Planet	Scenario	|$\dot{M}_{\mathrm{Na}}\,$|(kg s−1)	|$\dot{M}_{\mathrm{Na}}\,$|(kg s−1)
		prometheus	dishoom
	Escaping	105.3 ± 1	102.4 ± 0.3
WASP-49b	Exomoon	105.6 ± 0.6	104.3 ± 1.5
	Torus 1	|$10^{5.7^{+2.3}_{-0.8}}$|	103.1 ± 1.3
	Torus 2	|$10^{5.7^{+2.3}_{-0.8}}$|	106.4 ± 1.3
	Escaping	104 ± 0.2	102.3 ± 0.3
HD189733b	Exomoon	104.3 ± 0.2	104 ± 1.2
	Torus 1	104.2 ± 0.2	101.0 ± 1.5
	Torus 2	104.2 ± 0.2	104.1 ± 1.5

5.2 Inverse modelling of HD189733b

Wyttenbach et al. (2015) detected sodium at HD189733b in the Na i doublet. Compared to WASP-49b, this transit spectrum has weaker absorption features (less than one percent absorption at D2 line centre), but a larger D2-to-D1 line ratio of f_D2/D1 = 1.74 ± 0.45 (again calculated in bands of |$0.2\, \mathring{\rm A}$| centred on the sodium D lines). The line cores for the transit spectrum of HD189733b are broader than the ones of WASP-49b, which leads to more data points lying on the line cores and enabling a more precise retrieval. Therefore, and due to the comparatively smaller error bars for this observation, the best-fitting regions are more regular and significantly smaller for HD189733b (note that we adjusted the colour map scale for the parameter retrievals of HD189733b). The observation and the best-fitting models for each scenario are shown in Fig. 17. The retrieved parameters are summarized in Table 6, we note that the hydrostatic scenario with f_D2/D1 = 1.086 ± 0.001 significantly underpredicts the observed line ratio. All scenarios achieve a very similar goodness of fit with |$1.49\lt \chi _r^2\lt 1.53$|⁠, significantly larger than |$\chi _r^2=1$| (see Section 5.1 for an interpretation of these high |$\chi _r^2$| values).

Figure 17.

Transit spectrum for HD189733b with our best-fitting models for each scenario. Planetary parameters and data points are taken from Wyttenbach et al. (2015). Note the break in the x-axis.

We show the |$\chi _r^2$|-maps for all four scenarios in Fig. 18. As in the case of WASP-49b, the parameters in the hydrostatic scenario are driven to values such that the resulting atmosphere becomes optically thin. Since there is less absorption compared to WASP-49b, the resulting temperature is also lower: |$T=5200\pm 1500\,$|K, which is clearly above the planetary equilibrium temperature of 1140 K (Wyttenbach et al. 2015). The velocities and mass-loss rates are again very similar between the three evaporative scenarios. Compared to WASP-49b, we retrieve larger velocities and smaller source rates for HD189733b, which is due to the broader lines and smaller transit depth in the observation.

Figure 18.

Parameter retrieval for HD189733b. The colour-coding represents the deviation from |$\chi _r^2=1$|⁠, in standard deviations of the |$\chi _r^2$|-distribution (⁠|$\sigma =\sqrt{2/\nu }$|⁠, where ν corresponds to the degrees of freedom of the model, see Section 2.7). The contours are calculated on a 50 × 50 parameter grid and interpolated linearly using pyplot.contourf. The auxiliary parameter P₀ for the hydrostatic and escaping scenarios is computed using |$\chi _{\mathrm{Na}}=1.7\,$|ppm and (only for the escaping scenario) T = T_eq at the base of the wind.

We compare the retrieved mass-loss rates of prometheus to the rates computed within dishoom in Table 7. Again, the escaping wind does not generate enough sodium (by two orders of magnitude) to reproduce the observed transit depth. While the calculated dishoom source rates in the exomoon and desorbing torus (equation 19) scenarios are very comparable to the retrieved ones, a torus sourced by direct outgassing (equation 18) falls short of the retrieved source rate by three orders of magnitude.

6 DISCUSSION

6.1 Comparison to other retrievals

6.1.1 WASP-49b studies

Different authors conducted a parameter retrieval using the same data from Wyttenbach et al. (2017), as we used. In the following, we compare our findings to these studies. Wyttenbach et al. (2017) could fit the line cores and the line wings of the spectrum separately with hydrostatic, isothermal, and vertically mixed (constant mixing ratios throughout the atmosphere) models (equivalent to our hydrostatic scenario), but these authors could not reproduce the entire spectrum with a hydrostatic model. Their best-fitting model for the line cores has |$T=2950^{+400}_{-500}\,$| K. Cubillos et al. (2017) attempted to fit the spectrum with a more sophisticated (but still hydrostatic) model (with a T-P-profile from a hydrodynamic simulation and variable mixing ratios computed with equilibrium chemistry). However, they run into the same difficulty as Wyttenbach et al. (2017): The observed spectrum at WASP-49b has large absorption on the line cores (approximately two percent on the D2 line), but very little absorption between the line cores. As seen in Section 4, this property along with f_D2/D1 > 1.1 indicates that absorption occurs mainly in an extended, optically thin region. Since the hydrostatic models of Wyttenbach et al. (2017) and Pino et al. (2018) both retrieve temperatures smaller than |$3000\,$| K, they model a small and dense atmosphere leading to absorption in the optically thick regime, which is not able to reproduce the full observed transit spectrum.

Fisher & Heng (2019) used yet another hydrostatic model which is isothermal and vertically mixed, but incorporates NLTE-effects, clouds, and less restrictive priors for the temperature range. Since these authors allow for higher temperatures and lower reference pressures (due to cloud decks at high altitudes), their atmosphere is much more extended and optically thin and they were able to reproduce the observed spectrum. These authors retrieve a temperature of |$T=8415^{+1020}_{-1526}\,$| K (in their LTE scenario |$T=7209^{+1763}_{-1892}\,$| K), which is in line with our hydrostatic fit which has very high temperature (⁠|$T=6100_{-3400}\,$| K).

6.1.2 HD189733b studies

As in the case of WASP-49b, different authors conducted a parameter retrieval using the same data from Wyttenbach et al. (2015) as we used. Wyttenbach et al. (2015) could fit different parts of the transit spectrum of HD189733b with isothermal, vertically-mixed, hydrostatic models and interpreted this as a temperature gradient in the atmosphere. Their hottest model to fit the D2 line core has a temperature of |$3270\,$| K. Pino et al. (2018) combined the high-resolution data with low-resolution transmission spectra for HD189733b (Pont et al. 2008; Sing et al. 2011; Sing et al. 2016) and fit a vertically mixed (⁠|$\chi _{\mathrm{Na}}=1\, \mathrm{ppm}$|⁠), hydrostatic model with a customized T-P-profile to the low-resolution data. However, with the parameters retrieved from the low-resolution spectra these authors could not fit the high-resolution data from Wyttenbach et al. (2015), unless the T-P-profile has significantly larger temperatures at low pressures (reaching |$\approx 8000\,$| K at |$\sim 0.1\,$| nbar). Still, this high-temperature model with f_D2/D1 ≈ 1.1 does not reproduce the large line ratio observed at HD189733b.

Huang et al. (2017) performed a very sophisticated hydrostatic simulation of HD189733b’s atmosphere incorporating many different forms of atmospheric chemistry and NLTE effects. They find that the temperature rises steeply from |$2000\,$| K at |$10\, \mu$|bar to |$12^{\prime }000\,$| K at |$0.1\,$| nbar (their fig. 3). These authors furthermore find that at pressures below |$1\, \mu$|bar most of the sodium atoms are ionized (their fig. 5). Huang et al. (2017) used very different parameters (namely LyC boost and atomic layer base pressure) to fit the full observed spectrum, making a comparison to our scenarios difficult. However, we note that the found atomic layer base pressure of |$10\, \mu$|bar compares well to our retrieved reference pressure (⁠|$P_0=11\, \mu$|bar assuming |$\chi _{\mathrm{Na}}=1.7\,$| ppm) in the hydrostatic scenario, and our temperatures (⁠|$T=5200\pm 1500\,$| K), which are significantly larger than the ones retrieved by Wyttenbach et al. (2015), are in line with the atmospheric model from Huang et al. (2017). We also remark that these authors bin the observations in a slightly different way, using 0.05 |$\mathring{\rm A}$| bins. For this particular choice of the bin size, the D2-to-D1 line ratio of HD189733b is approximately 1.1, in line with the hydrostatic models.

6.2 Potassium in the atmosphere/exosphere of hot Jupiters

We have so far focused on the analysis of the Na i doublet. Our codes and the analysis of the D2/D1 line ratio can be readily applied to the K i doublet at |$7667\, \mathring{\rm A}$| and |$7701\, \mathring{\rm A}$|⁠. Here, we do not compare our models to observational data and refrain from a normalization and binning routine applied to the spectra as at Na i. We employ a HARPS-like instrumental LSF convolution procedure for our model spectrum. The model serves as a prediction for expected high-resolution K i detections embarked by ESPRESSO (Chen et al. 2020) and PEPSI (e.g. Keles et al. 2019).

Identical to the forward model presented in Section 4, we use dishoom to calculate K i mass-loss rates (Table 8), which are then converted into |$\mathcal {N}_{\mathrm{K}}$| (the number of neutral K i atoms in the system) using equation (6). The average lifetime of neutral K i is larger by a factor of 3.75 compared to that of Na i as determined by photoionization (Huebner & Mukherjee 2015). We set the velocities of the K I atoms in the evaporative scenarios uniformly to that of sputtered atoms |$\sim 10\,$| km s⁻¹. For the hydrostatic scenarios, we use the retrieved temperature and reference pressure (Table 6). For the endogenic, planetary scenarios we use a volumetric sodium-to-potassium ratio of Na/K =15.9 corresponding to the solar value (Asplund et al. 2009). For the exogenic, exomoon, and torus scenarios, we use lunar Na/K = 6 (Potter & Morgan 1988) that also corresponds to the maximum Na/K ratios for chondrites as studied by Fegley & Zolotov (2000). We note that the evaporative scenarios are strongly dependent on the Na/K ratios due to their optically thin nature. The simulated transmission spectra are shown in Fig. 19 (WASP-49b) and Fig. 20 (HD189733b). In comparison with the |$0.18{{\ \rm per\ cent}}$| K D1 absorption depth detection by Keles et al. (2019) for HD189733b, we find for the same 0.8 |$\mathring{\rm A}$| bandpass a K D1 absorption depth of |$\approx 0.12{{\ \rm per\ cent}}$| for the hydrostatic scenario, |$\approx 0.05{{\ \rm per\ cent}}$| for a the exomoon scenario, and |$\approx 0.04{{\ \rm per\ cent}}$| for a desorbing torus (with an enhanced source rate according to equation 19). At present, we believe more observations are needed to make interpretations of the Na/K ratio at exoplanets.

Figure 19.

Potassium transmission spectrum for WASP-49b. The dashed red line corresponds to an enhanced potassium source rate for a desorbing torus (equation 19). Note the break in the x-axis.

Figure 20.

Potassium transmission spectrum for HD189733b. The dashed red line corresponds to an enhanced potassium source rate for a desorbing torus (equation 19). Note the break in the x-axis.

Table 8.

Potassium source rates computed within dishoom. |$\dot{M}_{\mathrm{tor2,K}}$| corresponds to a desorbing torus (equation 19).

Source rate	WASP-49b	HD189733b
|$\dot{M}_{\mathrm{esc,K}}$| (kg s−1)	101.4 ± 0.3	101.3 ± 0.3
|$\dot{M}_{\mathrm{moon,K}}$| (kg s−1)	103.2 ± 1.6	102.9 ± 1.3
|$\dot{M}_{\mathrm{tor1,K}}$| (kg s−1)	102.9 ± 1.3	10−0.6 ± 1.5
|$\dot{M}_{\mathrm{tor2,K}}$| (kg s−1)	106.2 ± 1.3	102.6 ± 1.5

6.3 Toroidal atmospheres/exospheres at ultra-hot Jupiters

Ultra-hot Jupiters (⁠|$T\gtrsim 2000\, \mathrm{K}$|⁠) have recently been under spectral scrutiny after a remarkable detection of atomic iron at an exoplanet system (Hoeijmakers et al. 2018). Since then, several authors have detected not only Fe i but also Ti, V (Hoeijmakers et al. 2020), and other heavy metals (Hoeijmakers et al. 2018, 2019; Sing et al. 2019; Cabot et al. 2020; Gibson et al. 2020). However, the detection diaspora are vastly different at these bodies and the ‘ultra-hotness’ of the hot Jupiters does not appear to be a sufficient criterion for Fe i detections. For instance, Cauley et al. (2020) was unable to detect metals at WASP-189b using PEPSI (R ∼50 000) on the Large Binocular Telescope despite a planetary equilibrium temperature exceeding |$2640\, \mathrm{K}$|⁠. Therefore atmospheric heating, leading to heavy atom escape as proposed by Cubillos et al. (2020), should also lead to a Fe i signature at WASP-189b. To better understand the temperature-independent metallic detections, we directly compare the latest Na i detections at two ultra-hot Jupiters by the HARPS spectrograph: WASP-76b (Seidel et al. 2019) and WASP-121b (Cabot et al. 2020; Hoeijmakers et al. 2020). In Table 9, we provide system parameters of the two systems along with the measured D2-to-D1 ratios of the Na i doublet. Building on our ability to differentiate optically thin and optically thick gases from f_D2/D1 in Fig. 11 the planets, despite similar equilibrium temperatures (2360 K and 2190 K) and Fe i detections, appear to be in dramatically different gas regimes based on the Na I observations.

Table 9.

System parameters for the ultra-hot Jupiters orbiting F stars: WASP-121b (Hoeijmakers et al. 2020) and WASP-76b (Seidel et al. 2019). Average sodium lifetimes against photoionization are estimated from Huebner & Mukherjee (2015).

Parameter	WASP-121b	WASP-76b
R* (R⊙)	1.46	1.3
R0 (RJ)	1.87	1.83
Mp (MJ)	1.18	0.92
Teq (K)	2360	2190
γ (km s−1)	38	−1.07
tNa (s)	50	85
fD2/D1	2	1

6.3.1 Toroidal atmospheres

In light of the stark detection of f_D2/D1 ≈ 2 at WASP-121b in a companion paper by Hoeijmakers et al. (2020), we decide to model an identical exogenic, toroidal geometry to WASP-76b to better understand the physical mechanism fueling the (presumably temperature-independent) metal detections at ultra-hot Jupiters.

So far, the Na i and Fe i signatures at WASP-76b have been interpreted as endogenic atmospheric winds (Seidel et al. 2019; Ehrenreich et al. 2020). The scenario of day-night migration coupled with rotation leading to an evening/morning or dusk/dawn asymmetry has surprisingly also been observed in an exospheric regime on other tidally-locked bodies (Ganymede: Leblanc et al. 2017; Europa: Oza et al. 2019a). In Oza et al. (2018) it was shown that the phenomena of asymmetric ’atmospheric bulges’ is degenerate with a collisionless gas on a tidally-locked body. Therefore, while we find an asymmetric Fe I atmospheric wind quite reasonable, we test the degeneracy by simulating Na i rotating in a toroidal atmosphere (red line: WASP-76b f_D2/D1 ≈ 1) or exosphere (green line: WASP-121b f_D2/D1 ≈ 2) in Fig. 21.

Figure 21.

Transit spectrum for WASP 76-b (red squares) and WASP 121-b (green dots) comparing optically thick (red model) and optically thin (green model) toroidal atmospheres and exospheres.

We find the WASP-121b observations (green) are roughly reproduced by fixing a Na source orbiting at ∼ 28 km s⁻¹ ejecting ∼10⁹ g s⁻¹ extending to |$\sim 1.9\, R_0$|⁠, equivalent to the planet’s Hill radius. This corresponds to a toroidal scale height of H_{t, W121} ∼ R_J/4 as described in Section 2.6 (Case 1) and equation (18) appropriate for an exo-Enceladus sourcing a toroidal exosphere. In stark contrast, WASP-76b (red) is twice as thick with H_{t, W76} ∼ R_J/2. At WASP-76b, the torus is confined close to the planetary atmosphere with a much slower orbital speed of v_orb ∼ 7 km s⁻¹ at |$1.125\, R_0$|⁠. For grain densities |$\sim 3\, \mathrm{g}\, \mathrm{ cm} ^{-3}$|⁠, the Roche radius for WASP-76b is roughly |$\sim 0.97\, R_0$|⁠. On the other hand, the required Na rate here for a silicate exoring is alarmingly large |$\sim 10^{9}\,$|kg s⁻¹ of pure atomic Na. This rate is 20× as large as the total energy-limited escape rate of the H/He envelope, assuming an XUV-efficiency of η_XUV = 0.3 (c.f. Section 2.4). To our knowledge, the only physical process capable of generating such a high quantity of Na i in a toroidal distribution is the thermal desorption of silicate grains as shown in table 5 of Oza et al. (2019b) and described in Section 2.6 (Case 2) and equation (19). As indicated, an exo-Io at WASP-76b would be catastrophically destroyed due to the mass-loss as also implied by the radiative hydrodynamic simulations of Perez-Becker & Chiang (2013). The evaporation of such a body would imply a large quantity of metallic gas, as observed.

7 CONCLUSIONS

The advent of high-resolution transmission spectroscopy of transiting exoplanets has enabled a closer look at the environment immediately surrounding the exoplanet. By modelling the sodium D doublet, whose resonance lines are exceptionally bright given even sparse column densities ≳ 10¹⁰ cm⁻² (c.f. Fig. 10), we investigate three non-hydrostatic, evaporative scenarios in addition to the canonically used hydrostatic model atmosphere. To this end we use a custom-built radiative transfer code (prometheus), coupled via sodium mass-loss rates from a metal evaporation model (dishoom) to perform both forward and inverse modelling of high-resolution transmission spectra in the sodium doublet.

The three additional scenarios describe (i) an endogenically escaping medium (equation 9), (ii) an exogenic outgassed cloud sourced by an exomoon (equation 12), and an (iii) exogenic torus representing circumplanetary material (equation 15). Profile (i) of an escaping atmosphere has been simulated in detail and observed across several species over the past decade. Should the exoplanet be close-in, it is highly likely that it is experiencing extreme atmospheric escape on the order of ∼10⁷–10¹⁰ kg s⁻¹ (Wyttenbach et al. 2020) corresponding to ∼10–10⁴ kg s⁻¹ of pure sodium escape given solar abundance. The exogenic profiles are based on a recent exomoon study by Oza et al. (2019b) as well as comparative Solar system studies of Jupiter and Saturn’s (cryo-)volcanically active moons Io and Enceladus (Johnson et al. 2006b; Johnson & Huggins 2006). At first glance, we find that both hydrostatic and non-hydrostatic scenarios can fit HARPS transit spectra of WASP-49b and HD189733b.

In mitigating the apparent sodium degeneracy, we find that first determining whether the absorption occurs in a primarily optically thick or optically thin regime is critical (Section 3). A diagnostic based on the ratio of transit depths at the D2 and D1 line centres, f_D2/D1, is shown to be indicative of an optically thin or optically thick regime in Section 3.4. Hydrostatic models we find, despite arbitrary heating at |$T \sim 10^4\, \mathrm{K}$|⁠, cannot achieve line ratios larger than ≈1.2. This is a consequence of the exponentially decaying number density profile prescribed by hydrostatic equilibrium in contrast to the tenuous and extended profiles of the non-hydrostatic scenarios primarily leading to f_D2/D1 ≫ 1. Given multiple observations with line ratios greater than one (WASP-121b, HD189733b, WASP-49b) an inclusion of evaporative sources of the sodium absorption seems warranted. We find that the evaporative torus scenario is nevertheless able to fit planets with f_D2/D1 ≈ 1 (e.g. WASP-76b) in Section 6.3.1. Upon analyzing the Na I source rates required to fit toroidal atmospheres to ultra-hot Jupiter spectra, we remark that a common source for Fe I and Na I due to the evaporation of a rocky body may be reasonable.

Based on an evaporative curve of growth (Fig. 10), we find that for optically thin exospheres (as produced by our three evaporative scenarios), the number of neutral sodium atoms in the system, |$\mathcal {N}_{\mathrm{Na}}$|⁠, is directly proportional to the transit depth or the equivalent widths of the Na D lines. Equivalently, observed equivalent widths can be used to constrain sodium source rates. For optically thick atmospheres (as in a hydrostatic setting), however, it is not the total amount of sodium which sets the transit depth but rather the (wavelength-dependent) location at which the atmosphere becomes transparent (the transit radius R_λ). While the transit radius can always be calculated from an observed transit depth (equation 27), the analytical framework presented in Section 3.1 is uniquely applicable to hydrostatic atmospheres and breaks down for non-hydrostatic number density profiles. Hence, the transit radius loses its physical meaning as border between optically thick and optically thin chords for the evaporative scenarios, rendering R_λ obsolete.

For observed high-resolution transmission spectra of WASP-49b and HD189733b, we find that all four scenarios can be fit to the data, from a statistical point of view. From a physical point of view, however, we remark that the escaping scenario and the torus scenario with direct outgassing (equation 18) are not able to supply the retrieved source rates, according to our mass-loss calculations. Caution should be used in adjusting the solar abundance here in that an atmospheric metal enrichment enhancement ≳100χ_⊙ is unlikely (Thorngren & Fortney 2019). Furthermore, we note that while the hydrostatic scenario can fit the high-resolution observations in terms of |$\chi _r^2$|⁠, a few critical issues remain. (1) The atmospheres require significant thermospheric heating (T ≈ 5 × T_eq) and (2) an observed line ratio of f_D2/D1 ≥ 1.2 (significant for HD189733b), whereas our best-fitting models yield f_D2/D1 ≤ 1.2 despite any aforementioned heating sources at HD189733b and WASP-49b. Evaporative scenarios naturally explain (1) as line broadening due to fast sodium atoms (c.f. Table 3) and due to their inherently tenuous mediums also yield (2) f_D2/D1 ≫ 1.2 (Fig. 11). Physically, therefore, we find an exogenic cloud outgassed from an exomoon or a desorbing torus to be most favorable for WASP-49b and HD189733b .

At present, acknowledging that hydrostatic profiles have been excellent approximations to low-resolution transmission spectra, we suggest that radiative transfer and atmospheric escape models would benefit from the non-hydrostatic framework we have outlined in this paper for high-resolution transmission spectra. Further transmission spectra observations, at higher SNR, could validate non-hydrostatic line ratios of the alkaline resonance lines, in significant excess of ≈1.2. Time-resolved observations such as phase curves could provide insight regarding the spatial distribution of alkali atoms.

Looking forward, the recent identification of heavy metals beyond the Hill sphere of gas giant exoplanets by high-resolution transmission spectroscopy (Hoeijmakers et al. 2019; Cubillos et al. 2020) is reminiscent of the remarkable identification of exocomets or falling evaporating bodies at the β Pictoris system (e.g. Roberge et al. 2000; Lecavelier des Etangs et al. 2001). In this light, based on the evaporative transmission spectra modelling for gas giant exoplanets carried out here, we specifically suggest further modelling of escaping metals subject to ambient plasma fields, along with the likelihood of satellites and tori as at Jupiter and Saturn.

ACKNOWLEDGEMENTS

We sincerely extend our gratitude to A. Wyttenbach, J.V. Seidel, and H.J. Hoeijmakers for sharing the high-resolution Na i data analysed here. We thank R.E. Johnson for insight on the exomoon and tori models. We appreciate the radiative transfer discussions with E. Lellouch, C. Huang, and D. Kitzmann during the preparation of the manuscript. We also thank R. Ottersberg for computational help within prometheus, and acknowledge insightful discussions on mass-loss with R. Murray-Clay. AO acknowledges support from SNSF grant ‘Planets in Time.’ We thank our reviewer I. Waldmann for helpful comments, which improved the quality of this work. We have benefited from the public available programming language python, including the numpy(van der Walt, Colbert & Varoquaux 2011), matplotlib(Hunter 2007), and scipy(Virtanen et al. 2020) packages.

Data availability: No new data were generated or analysed in support of this research.

Footnotes

1

Exogenic refers to an external source whereas endogenic refers to a planetary source.

2

Based on the earlier η code from Ehrenreich et al. (2006).

3

Relative abundance by number. We denote mass mixing ratios by x_i.

4

We encountered this degeneracy in the hydrostatic setting between P₀ and χ_i, as our knowledge of the temperature allowed for a direct conversion between pressures and number densities. As we do not know the temperature in the escaping wind we formulate the degeneracy in this scenario as one between n₀ and χ_i.

5

For an analogous derivation of optically thin gas, please see appendix B of Hoeijmakers et al. (2020) for a derivation of a homogenous, optically thin slab of gas. Here, we present a heterogenous, dynamic gas.

6

The purpose of this model is to demonstrate the spectral imprint of extreme thermospheric heating.

7

If not specified otherwise, by optical depth we mean the optical depth at sodium D2 line centre for the following discussion.

8

We use that the probability of a model is proportional to exp (− χ²/2).