An emission spectrum for WASP-121b measured across the 0.8-1.1 micron wavelength range using the Hubble Space Telescope

TL;DR

This study presents Hubble Space Telescope observations of WASP-121b's emission spectrum from 0.8 to 1.1 microns, revealing atmospheric properties, thermal inversion, and chemical composition through spectral retrieval analysis.

Contribution

It extends the wavelength coverage of WASP-121b's emission spectrum and provides detailed atmospheric retrieval including thermal dissociation and ionization effects.

Findings

Emission due to H- shortward of 1.1 micron

No evidence for the previously reported 1.25 micron VO feature

Estimated atmospheric metallicity and C/O ratio consistent with solar values

Abstract

WASP-121b is a transiting gas giant exoplanet orbiting close to its Roche limit, with an inflated radius nearly double that of Jupiter and a dayside temperature comparable to a late M dwarf photosphere. Secondary eclipse observations covering the 1.1-1.6 micron wavelength range have revealed an atmospheric thermal inversion on the dayside hemisphere, likely caused by high altitude absorption at optical wavelengths. Here we present secondary eclipse observations made with the Hubble Space Telescope Wide Field Camera 3 spectrograph that extend the wavelength coverage from 1.1 micron down to 0.8 micron. To determine the atmospheric properties from the measured eclipse spectrum, we performed a retrieval analysis assuming chemical equilibrium, with the effects of thermal dissociation and ionization included. Our best-fit model provides a good fit to the data with reduced chi^2=1.04. The data diverge from a blackbody spectrum and instead exhibit emission due to H- shortward of 1.1 micron. The best-fit model does not reproduce a previously reported bump in the spectrum at 1.25 micron, possibly indicating this feature is a statistical fluctuation in the data rather than a VO emission band as had been tentatively suggested. We estimate an atmospheric metallicity of [M/H]=1.09(-0.69,+0.57), and fit for the carbon and oxygen abundances separately, obtaining [C/H]=-0.29(-0.48,+0.61) and [O/H]=0.18(-0.60,+0.64). The corresponding carbon-to-oxygen ratio is C/O=0.49(-0.37,+0.65), which encompasses the solar value of 0.54, but has a large uncertainty.

Tables

Table 1. Table 1: MCMC results for the joint fit to the G102v1 and G102v2 white light curves. Quoted values are the posterior medians and uncertainties give the ± 34 plus-or-minus 34 \pm 34 % credible intervals about the median. Values adopted for fixed parameters are reported at the bottom of the table.

Free	G102v1	G102v2
$D$ (ppm)	${682}_{-73}^{+73}$
$T_{\text{mid}}$ (MJD)	${58063.7624}_{-0.0023}^{+0.0053}$	${58096.9047}_{-0.0015}^{+0.0024}$
$\beta$	${1.21}_{-0.08}^{+0.08}$	${1.15}_{-0.09}^{+0.09}$
$\sigma$ (ppm)	${87}_{-6}^{+6}$	${83}_{-6}^{+6}$
Fixed	Value	Reference
$P$ (d)	$1.2749255$	Delrez et al. (2016)
$a/R_{\star }$	$3.86$	Evans et al. (2018)
$i$ (∘)	$89.1$	Evans et al. (2018)
$b$	$0.06$	Evans et al. (2018)

Table 2. Table 2: Similar to Table 1 , but reporting the MCMC results for the joint fit to the G102v1 and G102v2 spectroscopic light curves.

		G102v1			G102v2
$\lambda$ ($\mu \text{m}$)	$D$ (ppm)	$\beta$	$\sigma$ (ppm)		$\beta$	$\sigma$ (ppm)
$0.800$-$0.820$	${367}_{-152}^{+154}$	${1.10}_{-0.08}^{+0.08}$	${476}_{-33}^{+34}$		${1.04}_{-0.08}^{+0.08}$	${449}_{-34}^{+35}$
$0.820$-$0.839$	${330}_{-107}^{+106}$	${1.03}_{-0.08}^{+0.08}$	${391}_{-30}^{+31}$		${1.03}_{-0.07}^{+0.07}$	${389}_{-25}^{+26}$
$0.839$-$0.859$	${487}_{-77}^{+81}$	${0.99}_{-0.07}^{+0.08}$	${342}_{-24}^{+27}$		${1.09}_{-0.07}^{+0.07}$	${377}_{-25}^{+26}$
$0.859$-$0.878$	${657}_{-81}^{+81}$	${1.06}_{-0.07}^{+0.07}$	${350}_{-23}^{+24}$		${1.01}_{-0.08}^{+0.08}$	${330}_{-25}^{+26}$
$0.878$-$0.898$	${595}_{-69}^{+72}$	${1.00}_{-0.07}^{+0.07}$	${305}_{-22}^{+23}$		${0.95}_{-0.08}^{+0.08}$	${288}_{-25}^{+25}$
$0.898$-$0.917$	${574}_{-65}^{+70}$	${1.06}_{-0.07}^{+0.08}$	${310}_{-21}^{+22}$		${1.02}_{-0.07}^{+0.07}$	${299}_{-21}^{+21}$
$0.917$-$0.937$	${667}_{-63}^{+67}$	${1.06}_{-0.07}^{+0.07}$	${298}_{-19}^{+21}$		${0.99}_{-0.07}^{+0.07}$	${277}_{-20}^{+21}$
$0.937$-$0.956$	${671}_{-83}^{+78}$	${1.01}_{-0.08}^{+0.08}$	${277}_{-21}^{+22}$		${1.00}_{-0.08}^{+0.08}$	${274}_{-21}^{+22}$
$0.956$-$0.976$	${743}_{-94}^{+98}$	${1.08}_{-0.10}^{+0.09}$	${291}_{-27}^{+24}$		${1.08}_{-0.08}^{+0.08}$	${289}_{-20}^{+22}$
$0.976$-$0.995$	${733}_{-65}^{+70}$	${1.14}_{-0.07}^{+0.07}$	${303}_{-18}^{+19}$		${1.18}_{-0.07}^{+0.07}$	${314}_{-18}^{+19}$
$0.995$-$1.015$	${797}_{-68}^{+64}$	${1.12}_{-0.07}^{+0.07}$	${298}_{-18}^{+19}$		${1.03}_{-0.07}^{+0.08}$	${276}_{-20}^{+21}$
$1.015$-$1.034$	${795}_{-67}^{+70}$	${1.14}_{-0.07}^{+0.07}$	${298}_{-18}^{+19}$		${0.98}_{-0.08}^{+0.08}$	${258}_{-20}^{+21}$
$1.034$-$1.054$	${736}_{-56}^{+58}$	${1.02}_{-0.07}^{+0.07}$	${269}_{-19}^{+20}$		${1.00}_{-0.07}^{+0.07}$	${265}_{-19}^{+20}$
$1.054$-$1.073$	${852}_{-63}^{+62}$	${1.14}_{-0.07}^{+0.07}$	${306}_{-18}^{+20}$		${1.03}_{-0.07}^{+0.08}$	${277}_{-19}^{+21}$
$1.073$-$1.093$	${832}_{-62}^{+62}$	${1.04}_{-0.07}^{+0.07}$	${283}_{-20}^{+20}$		${1.10}_{-0.07}^{+0.07}$	${300}_{-20}^{+20}$
$1.093$-$1.112$	${791}_{-65}^{+62}$	${1.08}_{-0.07}^{+0.07}$	${299}_{-19}^{+20}$		${1.01}_{-0.07}^{+0.08}$	${279}_{-20}^{+21}$
$1.112$-$1.132$	${895}_{-74}^{+71}$	${1.05}_{-0.08}^{+0.08}$	${292}_{-21}^{+21}$		${1.04}_{-0.07}^{+0.08}$	${289}_{-20}^{+21}$

Table 3. Table 3: Retrieval prior ranges and MCMC results

Parameter	Unit	Allowed range	Result
$\left[\text{M}/\text{H}\right]$	dex	$-1$ to $2$	${1.09}_{-0.69}^{+0.57}$
$\left[\text{C}/\text{H}\right]$	dex	$-1$ to $2$	$-{0.29}_{-0.48}^{+0.61}$
$\left[\text{O}/\text{H}\right]$	dex	$-1$ to $2$	${0.18}_{-0.60}^{+0.64}$
${\log }_{10}\left({\kappa }_{\text{IR}}\right)$	dex cm2 g-1	$-5$ to $0.5$	$-{3.01}_{-0.62}^{+0.56}$
${\log }_{10}\left(\gamma \right)$	dex	$-4$ to $1.5$	${0.64}_{-0.16}^{+0.19}$
$\psi$	—	$0$ to $2$	${0.99}_{-0.09}^{+0.06}$

Full text

††pubyear: 2019††pagerange: An emission spectrum for WASP-121b measured across the $0.8$–$1.1\,\mu\textnormal{m}$ wavelength range using the Hubble Space Telescope–An emission spectrum for WASP-121b measured across the $0.8$–$1.1\,\mu\textnormal{m}$ wavelength range using the Hubble Space Telescope

1 Introduction

Eclipse observations made as an exoplanet passes through superior conjunction allow the emission from the planetary dayside hemisphere to be inferred. Dozens of eclipse measurements have been published for broad photometric passbands, most notably from the Kepler space mission (e.g. Borucki et al., 2009; Kipping & Spiegel, 2011; Faigler et al., 2013; Angerhausen et al., 2015; Esteves et al., 2015), the Spitzer Space Telescope Infrared Array Camera (IRAC) (e.g. Charbonneau et al., 2005; Knutson et al., 2010; Stevenson et al., 2010; Blecic et al., 2014; Garhart et al., 2019), Infrared Spectrograph (IRS) (Stevenson et al., 2012; Cubillos et al., 2014), and Multiband Imaging Photometer (Deming et al., 2005; Charbonneau et al., 2008; Crossfield et al., 2012), as well as ground-based telescopes (e.g. Sing & López-Morales, 2009; Gibson et al., 2010; de Mooij et al., 2011; Mallonn et al., 2019; Kovács & Kovács, 2019). A smaller number of spectroscopic eclipse measurements have also been published, including two with the Spitzer IRS (Richardson et al., 2007; Grillmair et al., 2007, 2008), two with the Hubble Space Telescope (HST) Near-Infrared Camera Multi-Object Spectrometer (Swain et al., 2009a, b), and fourteen with the HST Wide Field Camera 3 (WFC3) infrared spectrograph (Wilkins et al., 2014; Ranjan et al., 2014; Stevenson et al., 2014b, a; Kreidberg et al., 2014, 2018; Crouzet et al., 2014; Haynes et al., 2015; Line et al., 2016; Beatty et al., 2017; Sheppard et al., 2017; Evans et al., 2017; Nikolov et al., 2018b; Arcangeli et al., 2018; Mansfield et al., 2018). Spectroscopic observations are of particular value by allowing opacity bands to be resolved, which in turn encode information about chemical abundances and the vertical temperature profile of the atmosphere.

The WFC3 instrument offers two grisms for infrared spectroscopy: G102 covering wavelengths $0.8$-$1.1\,\mu\textnormal{m}$ and G141 covering wavelengths $1.1$-$1.6\,\mu\textnormal{m}$. All of the exoplanet emission spectra published to date have used the G141 grism for two main reasons: (i) the longer wavelength coverage provides a more favorable planet-to-star brightness ratio; and (ii) the G141 passband provides access to stronger opacity bands than the G102 passband, in particular the $1.4\,\mu\textnormal{m}$ H2O band. The observations made with WFC3 G141 have resulted in detections of H2O absorption for WASP-43b (Stevenson et al., 2014b; Kreidberg et al., 2014), HD 189733b (Crouzet et al., 2014), HD 209458b (Line et al., 2016), and Kepler-13Ab (Beatty et al., 2017), and H2O emission for WASP-121b (Evans et al., 2017). Other spectral features reported include TiO emission for WASP-33b (Haynes et al., 2015) and CO absorption for WASP-18b (Sheppard et al., 2017). The remaining seven spectra do not exhibit significant spectral features and typically appear blackbody-like (Wilkins et al., 2014; Ranjan et al., 2014; Nikolov et al., 2018b; Mansfield et al., 2018; Kreidberg et al., 2018), while the CO absorption feature claimed for WASP-18b has been challenged (Arcangeli et al., 2018). Recent studies (Arcangeli et al., 2018; Lothringer et al., 2018; Parmentier et al., 2018; Kreidberg et al., 2018) have suggested that continuum opacity due to H- and thermal dissociation of H2O itself can explain the lack of spectral features observed for the ultrahot Jupiters, such as WASP-18b, HAT-P-7b, and WASP-103b, which have dayside temperatures well in excess of 2000 K.

Another such ultrahot Jupiter, WASP-121b, is the subject of this study. Discovered by Delrez et al. (2016), WASP-121b has an exceptionally inflated radius ($1.8\,R_{J}$) and a dayside photospheric temperature of approximately 2700 K (Evans et al., 2017). The high temperature results from WASP-121b orbiting at a distance of only 0.025 AU from its F6V host star, where it is subjected to strong tidal forces and may be undergoing atmospheric escape via Roche lobe overflow (Delrez et al., 2016). Observational support for this picture has recently been provided by near-ultraviolet (NUV) transit measurements made with Swift UVOT that are significantly deeper than those measured at optical wavelengths (Salz et al., 2019). This could be explained by an extended atmosphere filling the Roche lobe, which is relatively opaque at NUV wavelengths due, for instance, to heavy metal absorption lines.

At near-infrared wavelengths, the HST WFC3 spectrograph has been used to measure both the transmission spectrum (Evans et al., 2016) and emission spectrum (Evans et al., 2017) of WASP-121b. Absorption due to the $1.4\,\mu\textnormal{m}$ H2O band is observed in the transmission spectrum, while this same band is seen in emission at secondary eclipse, revealing a thermal inversion on the dayside hemisphere. The latter indicates significant absorption of incident stellar radiation at NUV-optical wavelengths for pressures less than $\sim 100$ mbar on the dayside hemisphere (e.g. Guillot, 2010). This is consistent with the low geometric albedo of $A_{g}=0.16\pm 0.11$ inferred by Mallonn et al. (2019) from the $z^{\prime}$ secondary eclipse measurement of Delrez et al. (2016). Possible absorbers include TiO and VO, both of which have strong opacity bands throughout the optical (Hubeny et al., 2003; Fortney et al., 2008). Indeed, the optical transmission spectrum measured using the HST Space Telescope Imaging Spectrograph (STIS) does show evidence for VO absorption at the day-night limb, although TiO is not seen (Evans et al., 2018). A steep rise toward NUV wavelengths is also recovered in the STIS transmission spectrum, which may be caused by the same absorber/s responsible for the deep Swift UVOT transit. One candidate proposed in Evans et al. (2018) is SH, which has been predicted as a product of non-equilibrium chemistry in hot Jupiter atmospheres by Zahnle et al. (2009), and, if present on the dayside hemisphere, could potentially produce the thermal inversion. As alluded to above, absorption by heavy metals such as Fe and Mg might also simultaneously explain the deep NUV transits and dayside thermal inversion. At optical wavelengths, other candidate absorbers that could play a role in generating the thermal inversion include H- ions and molecules such as NaH, MgH, FeH, SiO, AlO, and CaO (Lothringer et al., 2018; Parmentier et al., 2018; Gandhi & Madhusudhan, 2019), although no compelling evidence has been claimed for any of these species based on the published transmission spectrum (Evans et al., 2016, 2018).

Unlike the transmission spectrum — which probes a region of the atmosphere very different to the ultrahot dayside of WASP-121b — a detection of one or more strong optical absorbers in the emission spectrum would provide a definitive link between the radiatively active species present and the thermal inversion. Motivated by this, we acquired secondary eclipse observations of WASP-121b with the G102 grism of WFC3, extending the wavelength coverage into the red optical where emission bands due to species such as TiO, VO, and FeH may be detectable, as well as H- continuum opacity. We describe these observations and our data reduction in Section 2. Our analyses of the white and spectroscopic light curves are presented in Sections 3 and 4, respectively. The results are discussed in Section 5 and we conclude in Section 6.

2 Observations and data reduction

We observed two secondary eclipses of WASP-121b with HST/WFC3 using the G102 grism, which covers a wavelength range of approximately $0.8$-$1.1\,\mu\textnormal{m}$ with a spectral resolving power of $R\sim 200$ at $\lambda=1\,\mu\textnormal{m}$ (G.O. 15135; P.I. Mikal-Evans). The first visit was made on 2017 November 6 and the second visit was made on 2017 December 9. We refer to these two visits as the G102v1 and G102v2 datasets, respectively. For both visits, the target was observed for 6.9 hours over five consecutive HST orbits with identical observing setups. The relative timing of the two visits was designed to provide full phase coverage of the eclipse, using the previously determined ephemerides of WASP-121b. Observations were made in spectroscopic mode with a forward scanning rate of $0.062$ arcsec s*-1* along the cross-dispersion axis. To reduce overheads, only the $512\times 512$ pixel subarray of the detector containing the target spectrum was read out for each exposure. We adopted the SPARS10 sampling sequence with 15 non-destructive reads per exposure ($\textnormal{NSAMP}=15$) resulting in total integration times of 103 s and scans across approximately 50 pixel rows of the cross-dispersion axis. With this setup, we obtained 14 exposures in the first HST orbit following acquisition and 16 exposures in each subsequent HST orbit. Typical peak frame counts were $\sim 32,000$ electrons per pixel for both visits. This translates to $\sim 13,000$ data numbers (DN) per pixel given the detector gain of 2.5 electrons per DN, which is well within the linear regime of the WFC3 detector (see Figure 1 of Hilbert, 2008).

Spectra were extracted from the raw data frames using a custom-built Python pipeline, which has been described previously (Evans et al., 2016, 2017). We started with the IMA files produced by the calwf3 pipeline version 3.4.1, which already have basic calibrations such as flat fielding, bias subtraction, and nonlinearity correction applied. The target flux was extracted from each exposure by taking the difference between successive non-destructive reads. To do this, we first estimated and subtracted the background flux for each read, by taking the median pixel count within a $10\times 170$ pixel box which was chosen to be as large as possible while avoiding sources within the field and the detector edges. Typical background levels integrated over the full $103\,$s exposures were approximately $70$–$80$ electrons pixel*-1*, rising to over $100$ electrons pixel*-1* at the end of each HST orbit (Figure 1). For each read-difference frame, we then determined the flux-weighted centre of the scanned spectrum along the cross-dispersion axis. All pixel values located more than 30 pixels above and below this row were set to zero, effectively removing flux contributions from nearby contaminant stars and cosmic ray strikes outside a rectangular aperture. Final reconstructed frames were produced by adding together the read-differences produced in this manner. During this process, we also estimated how the spectrum drifted across detector over the course of the observations. For both visits, we measure a drift of $\sim 0.1$-$0.2$ pixel along the dispersion axis and $\sim 0.6$ pixel along the cross-dispersion axis (Figure 1).

The target spectrum was then extracted from each frame by summing the flux within a rectangular aperture spanning the full dispersion axis and 80 pixels along the cross-dispersion axis, centered on the central cross-dispersion row of the scan. The wavelength solution was determined by cross correlating each of these extracted spectra against a model spectrum for the WASP-121 host star modulated by the throughput of the G102 grism, as in Evans et al. (2016, 2017). In addition to the G102 data, a single secondary eclipse of WASP-121b was observed on 2016 Nov 10 with the G141 grism (G.O. 14767; P.I.s Sing and López-Morales). This dataset was originally published in Evans et al. (2017), to which the reader is referred for further details. Example G102 and G141 spectra are shown in Figure 2. Together, both grisms provide continuous wavelength coverage between $\sim 0.8$-$1.6\,\mu\textnormal{m}$.

3 White Light Curve Analyses

White light curves were produced for both visits by integrating each spectrum across the full dispersion axis (Figure 1). The light curves are affected by the well-known hook systematic that correlates with HST orbital phase and is understood to be caused by charge trapping in the WFC3 detector (Zhou et al., 2017). The baseline flux level also exhibits a longer term drift, which is approximately linear in time for both light curves.

Prior to light curve fitting, we chose to discard the first HST orbit of each visit, as these exhibit much stronger hooks than subsequent orbits. Although methods exist to correct the WFC3 first-orbit systematics (e.g. Zhou et al., 2017; de Wit et al., 2018), we opted for this approach to be consistent with our previous analyses (Evans et al., 2016, 2017) and to avoid modeling the baseline trend over the full five-orbit visits, which is less likely to be well approximated as linear. The resulting four-orbit white light curves were fit using a similar methodology to that described in Evans et al. (2017), in which the systematics are treated as a Gaussian process (GP). In the present study, we modeled the eclipse signal using the batman software package (Kreidberg, 2015), allowing the eclipse depth ($D$) and eclipse mid-time ($T_{\textnormal{mid}}$) to vary as free parameters, while holding the remaining parameters fixed to previously determined values (Table 1). We varied the eclipse depths jointly for both visits and allowed the mid-times to vary separately for each visit. For the GP, we employed a Matérn $\nu=3/2$ kernel with HST orbital phase ($\phi$), dispersion drift ($x$), and cross-dispersion drift ($y$) as input variables. We chose not to use time ($t$) as an input variable, as from past experience we have found its inclusion to make little difference to the final result (Evans et al., 2018). The GP free parameters were the covariance amplitude ($A$) and the correlation length scales for each input variable ($L_{\phi},L_{x},L_{y}$). In practice, as in Evans et al. (2017, 2018), we fit for the natural log of the inverse correlation length scale, $\ln\eta_{k}=\ln L_{k}^{-1}$, where $k=\{\phi,x,y\}$. We also parameterized the white noise as $\sigma=\beta\sigma_{0}$, where $\sigma_{0}$ is the formal photon noise floor and $\beta$ is a rescaling factor that was allowed to vary in the fits. We adopted uniform priors for all eclipse parameters, and adopted the same priors as in Evans et al. (2018) for the remaining parameters. Marginalization of the posterior distribution was performed using affine-invariant Markov chain Monte Carlo (MCMC) as implemented by the emcee software package (Foreman-Mackey et al., 2013).

The resulting posterior distributions are summarized in Table 1, and the best-fit light curve models are shown in Figure 3. We obtain eclipse depth measurements of $682\pm 73$ ppm, with inferred $\beta$ values of $1.21\pm 0.08$ and $1.15\pm 0.09$ for the G102v1 and G102v2 light curves, respectively. The latter imply high frequency scatter approximately $20$ % above the photon noise floor for both light curves, which is not accounted for by the Matérn kernel and is evident in the model residuals shown in Figure 3. As a check, we also repeated the light curve fitting using the squared exponential kernel (see e.g. Gibson et al., 2012) and obtained results for the eclipse depth and mid-times that were fully consistent with those reported in Table 1 to well within $1\sigma$. However, the squared exponential fit gave uncertainties that were approximately 5% smaller for the eclipse depth, 40% smaller for the G102v1 mid-time, and 10% smaller for the G102v2 mid-time. For this reason, we adopt the results obtained with the Matérn kernel to be conservative.

4 Spectroscopic Light Curve Analyses

Spectroscopic light curves were constructed using a similar method to that described by Deming et al. (2013), which we have also used previously in Evans et al. (2016, 2017). This involved cross correlating each spectrum against a master spectrum constructed by taking the median of out-of-eclipse exposures, in order to remove wavelength-independent systematics, including those arising due to pointing drift across the detector dispersion axis over the course of each visit. The flux was then binned into the 17 wavelength channels shown in Figure 2, each spanning 8 pixel columns on the detector ($\Delta\lambda=0.02\,\mu\textnormal{m}$). The resulting light curves are shown in Figure 4.

To fit the spectroscopic light curves, we used the same approach as described in Section 3. The only exception was that we fixed $T_{\textnormal{mid}}$ to the best-fit values listed in Table 1. Thus, for the spectroscopic eclipse signals, the only free parameter was the eclipse depth $D$, which we varied jointly across both the G102v1 and G102v2 light curves. Systematics were again accounted for using GPs with Matérn $\nu=3/2$ kernels and white noise rescaling factors (i.e. $\beta$ parameters). The inferred eclipse depths and $\beta$ values are reported in Table 2.

The resulting emission spectrum is shown in Figure 5, as measured eclipse depths in the top panel and corresponding brightness temperatures in the bottom panel. Also shown are the results obtained when the G102v1 and G102v2 datasets are analyzed separately, and when the data are rebinned into 34 channels. Good agreement is obtained for all spectroscopic channels, verifying the repeatability of the measurement.

We also re-analyzed the G141 eclipse data published in Evans et al. (2017). To be fully consistent with the G102 analysis, we fixed the values of $a/R_{\star}$ and $i$ to the values listed in Table 1, whereas the original G141 analysis had adopted the values reported in Delrez et al. (2016). This gave statistically identical results to those reported in Evans et al. (2017), which is to be expected as $a/R_{\star}$ and $i$ primarily affect the eclipse duration, rather than the eclipse depth.

5 Discussion

The secondary eclipse spectrum measured to date for WASP-121b is shown in Figure 6, along with the corresponding brightness temperatures. In addition to the G102 data spanning $0.8$-$1.1\,\mu\textnormal{m}$ presented in the current study, this includes the G141 data spanning $1.1$-$1.6\,\mu\textnormal{m}$ from Evans et al. (2017), ground-based photometric measurements in the $z^{\prime}$ (Delrez et al., 2016) and $K_{s}$ (Kovács & Kovács, 2019) passbands, and the IRAC data at $3.6\,\mu\textnormal{m}$ and $4.5\,\mu\textnormal{m}$ from Garhart et al. (2019). The G102 and G141 measurements agree extremely well at the point of overlap between the two passbands, without any adjustment to the level of either dataset. Similarly, the G102 data are fully consistent with the $z^{\prime}$ measurement without any adjustment.

5.1 Blackbody fits and heat redistribution

To interpret the data, we first consider the simple case in which the planet is assumed to radiate as an isothermal blackbody. Fitting such a model to the full dataset gives a best-fit temperature of $2720\pm 8$ K, with predictions for the wavelength-dependent secondary eclipse depth indicated by the dark yellow line in Figure 6. The data approximately follow the shape of this curve, however, the reduced $\chi^{2}$ is 2.92 for 47 degrees of freedom, allowing it to be ruled out at $6.5\sigma$ confidence. For comparison, in Evans et al. (2017) we ruled out a blackbody model at $5\,\sigma$ confidence by fitting to the data available at that time; namely, the G141, $z^{\prime}$, and IRAC $3.6\,\mu\textnormal{m}$ measurements. The addition of the G102 and IRAC $4.5\,\mu\textnormal{m}$ measurements has therefore increased the discrepancy between the data and a blackbody model.

We also experimented with fitting a blackbody model to different subsets of the data. For instance, if we repeat the fit to the full dataset with the only exception being that we use the G141 white eclipse depth rather than the G141 spectroscopic eclipse depths, we obtain a best-fit temperature of $2754\pm 11$ K with a reduced $\chi^{2}$ of 1.97 for 21 degrees of freedom, reducing the confidence with which such a model can be ruled out to $2.8\sigma$. Alternatively, if we repeat the fit to the full dataset but exclude the G102 spectroscopic eclipse depths — similar to Evans et al. (2017) but with the addition of the $K_{s}$ and IRAC $4.5\,\mu\textnormal{m}$ points — we obtain a best-fit temperature of $2691\pm 9$ K with a reduced $\chi^{2}$ of 2.88 for 30 degrees of freedom, ruling it out at $5.2\sigma$ confidence.

These results imply that, statistically, the departure from a blackbody spectrum is largely driven by the spectroscopic information contained in the G141 data. In Evans et al. (2017), we attributed this departure to a muted H2O emission band at $1.4\,\mu\textnormal{m}$, as well as a tentative VO emission band at $1.25\,\mu\textnormal{m}$. We revisit this interpretation in Section 5.2 with a retrieval analysis of the updated dataset. Here we note that even without the G141 spectroscopic information, considerable tension remains between the data and a blackbody model, mainly due to the mismatch between the overall slope of the data and a blackbody spectrum. This can be appreciated in Figure 6, which shows a systematic decrease in brightness temperature over the near-infrared wavelength range covered by the G102 and G141 passbands. As we discuss in Section 5.2, much of this variation can be explained by the $1.4\,\mu\textnormal{m}$ H2O band, as well as continuum opacity due to H- ions.

Figure 7 shows the approximate range of Bond albedos ($A_{B}$) and heat redistribution efficiencies ($\varepsilon$) for WASP-121b allowed by the emission data. The $\varepsilon$ parameter is defined following Cowan & Agol (2011), with $\varepsilon=0$ corresponding to zero heat redistribution and $\varepsilon=1$ corresponding to uniform heat redistribution. Measurements of low Bond albedos for numerous hot Jupiters (e.g. Kipping & Spiegel, 2011; Heng & Demory, 2013; Bell et al., 2017; Močnik et al., 2018; Mallonn et al., 2019) suggest a low Bond albedo is also likely for WASP-121b, which would be in line with the evidence for significant optical absorption in the transmission spectrum (Evans et al., 2018). Under this scenario ($A_{B}\lesssim 0.1$), the heat redistribution efficiency would be $\varepsilon\approx 0.4$. However, if heat redistribution is inefficient ($\varepsilon\approx 0$), the maximum allowable Bond albedo is $A_{B}\approx 0.3$. A phase curve measurement would allow this degeneracy to be broken, by providing a direct constraint on the heat flux from the nightside hemisphere, and hence on $\varepsilon$.

5.2 Atmosphere modeling of the dayside hemisphere

We performed an atmospheric retrieval analysis for the secondary eclipse data shown in Figure 8 using the ATMO code, which has been described extensively elsewhere (Amundsen et al., 2014; Tremblin et al., 2015, 2016, 2017; Drummond et al., 2016; Goyal et al., 2018). In brief, ATMO was originally developed to solve the one-dimensional (1D) plane-parallel radiative transfer equation assuming hydrostatic balance and radiative-convective equilibrium, with a two-dimensional implementation subsequently introduced by Tremblin et al. (2017). For a given atmospheric composition and pressure-temperature (PT) profile, chemical equilibrium gas phase abundances can be determined by Gibbs energy minimization, or alternatively, arbitrary mixing ratios can be specified. The PT profile can also be provided as an input, or otherwise computed self-consistently given the atmospheric opacity sources, an internal heat flux, and the irradiation from the host star. Condensation can be treated either locally or using a rainout approach (see Goyal et al., 2019). The planetary spectrum viewed by an external observer is produced as output, allowing ATMO to be used for inferring atmospheric properties from primary transit (Evans et al., 2016, 2018; Wakeford et al., 2017, 2018; Nikolov et al., 2018a; Alam et al., 2018) and secondary eclipse (Evans et al., 2017; Nikolov et al., 2018b) measurements.

We previously used ATMO in Evans et al. (2017) to analyze the G141, $z^{\prime}$, and IRAC $3.6\,\mu\textnormal{m}$ emission data for WASP-121b. In that study, we assumed uniform mixing ratios with pressure and allowed the abundances of H2O and VO to vary as free parameters. In the present study, we add the G102, $K_{s}$, and IRAC $4.5\,\mu\textnormal{m}$ data to our analysis and make a number of changes to our retrieval methodology, motivated by recent work highlighting the importance of thermal dissociation and ionization in ultrahot Jupiter atmospheres (Arcangeli et al., 2018; Parmentier et al., 2018; Lothringer et al., 2018; Kreidberg et al., 2018; Mansfield et al., 2018). First, rather than fitting for unconstrained abundances of pre-defined opacity sources, we assumed chemical equilibrium with relative elemental abundances set to solar values and varied metallicity ($[\textnormal{M}/\textnormal{H}]$), as well as the carbon and oxygen elemental abundances ($[\textnormal{C}/\textnormal{H}]$, $[\textnormal{O}/\textnormal{H}]$) separately.111Under this formulation, $[\textnormal{M}/\textnormal{H}]$ controls the number density of all elements heavier than helium except for carbon and oxygen. Second, we accounted for condensation using the rainout scheme described in Goyal et al. (2019) and included the effect of gas phase scattering. Third, we allowed for thermal dissociation and ionization of atoms and molecules, which can result in strongly pressure-dependent abundances for many chemical species throughout layers of the atmosphere probed in emission. We note that ATMO has always accounted for these effects in determining chemical equilibrium abundances of neutral species, but this was not included in our previous retrieval, as we assumed a uniform distribution of H2O and VO in pressure and allowed the abundances to vary freely. Ionic species were, however, added to ATMO for the present study following Gordon & McBride (1994), and the ion-neutral composition was benchmarked against the open source GGChem code (Woitke et al., 2018). The resulting chemical system consisted of 175 neutral gas phase species, 93 condensates, and the ionized species e-, H+, H-, He+, Na+, K+, C+, Ca+, and Si+. The most important radiatively active neutral gas phase species were H2O, CO2, CO, CH4, NH3, Na, K, Li, R, Cs, TiO, VO, FeH, PH3, H2S, HCN, C2H2, SO2, and Fe(g). Collision-induced absorption due to H2-H2 and H2-He was also included.

As in Evans et al. (2017), we fit for the PT profile using the analytic solution derived by Guillot (2010), which is parameterized in terms of the infrared opacity ($\kappa_{\textnormal{IR}}$), the ratio of the visible-to-infrared opacity ($\gamma=\kappa_{\textnormal{V}}/\kappa_{\textnormal{IR}}$), and an irradiation efficiency factor ($\psi$). The latter is identical to the $\beta$ term defined by Line et al. (2013), but we denote it as $\psi$ to avoid confusion with the white noise rescaling factors used in Sections 3 and 4. It effectively allows for nonzero albedo values and varying degrees of heat recirculation from the dayside to nightside hemisphere.

We adopted uniform priors for all six retrieval parameters (i.e. $[\textnormal{M}/\textnormal{H}]$, $[\textnormal{C}/\textnormal{H}]$, $[\textnormal{O}/\textnormal{H}]$, $\log_{10}\kappa_{\textnormal{IR}}$, $\log_{10}\gamma$, $\psi$), with the allowed ranges listed in Table 3. Differential-evolution MCMC was used to marginalize the posterior distribution, using the publicly available software of Eastman et al. (2013). We ran twelve chains each for 30,000 steps and discarded the first 20% of each chain as burn-in before combining them into a single chain. The resulting parameter distributions are reported in Table 3 and shown in Figure 9. The best-fit emission spectrum is shown in Figure 8 and has a $\chi^{2}$ of 43.61 for $\nu=42$ degrees of freedom, indicating an excellent fit to the data with reduced $\chi^{2}_{\nu}=1.04$.

We infer a metallicity of $[\textnormal{M}/\textnormal{H}]={1.09}_{-0.69}^{+0.57}$, translating to a 68% credible interval of $\sim 5$-$50\times$ solar. This is consistent with the plausible metallicity range of $\sim 10$-$30\times$ solar we found for the transmission spectrum (Evans et al., 2018). Note that the WASP-121 host star itself shows only mild evidence for heavy metal enrichment relative to solar, with $[\textnormal{Fe}/\textnormal{H}]=0.13\pm 0.09$ (Delrez et al., 2016). We also obtain 68% credible intervals of $-0.77$ to $0.33$ for $[\textnormal{C}/\textnormal{H}]$ and $-0.42$ to $0.82$ for $[\textnormal{O}/\textnormal{H}]$. These latter ranges are somewhat lower than that inferred for $[\textnormal{M}/\textnormal{H}]$, but are within an order of magnitude. Figure 10 shows the corresponding carbon-to-oxygen ratio, obtained by combining the $[\textnormal{C}/\textnormal{H}]$ and $[\textnormal{O}/\textnormal{H}]$ MCMC chains. We estimate $\textnormal{C/O}=0.49_{-0.37}^{+0.65}$, which is consistent with the solar value of $0.54$ (Asplund et al., 2009). Together, the updated emission spectrum presented here and the transmission spectrum presented in Evans et al. (2018) suggest a metallicity of $\sim 5$-$50\times$ solar for the atmosphere of WASP-121b at pressures below $\sim 100$ mbar, while the C/O ratio is weakly constrained but fully consistent with that measured for the Sun.

These results are broadly in line with theoretical predictions for the heavy element content of gas giant atmospheres. For example, the interior structure models of Thorngren & Fortney (2019) predict a metallicity of $21\pm 4\ \times$ solar for the upper atmosphere of WASP-121b, well within the ranges favored by the transmission and emission spectra. However, our revised metallicity estimate is in contrast to the anomalously high abundances for H2O and VO reported in Evans et al. (2017) for the dayside atmosphere of WASP-121b. In that study we obtained a 95% credible lower limit of $1,000\times$ solar for the VO abundance, which was driven by the apparent flux excess measured across the $1.20$-$1.25\,\mu\textnormal{m}$ wavelength range. Although the equilibrium chemistry model presented here fails to replicate this feature (Figure 8), the overall goodness of fit it provides to the full dataset (i.e. $\chi^{2}_{\nu}=1.04$) suggests it is more likely a statistical fluctuation than a VO emission band.

The ability of the chemical equilibrium model to fit the data with abundances closer to solar values is primarily due to the effects of thermal dissociation and ionization, and the opacities of the resulting products such as H-, which as noted above were not treated in our previous retrieval analysis. Parmentier et al. (2018) were the first to include these effects in a detailed analysis of the WASP-121b dayside emission data and made a similar observation. Those authors employed a three-dimensional (3D) general circulation model (GCM), which has the advantage of treating the radiative transfer and dynamics of the atmosphere self-consistently. However, a disadvantage of 3D GCMs is their computational expense, which typically prevents the free parameters (e.g. metallicity, frictional drag, etc.) being optimized to match a given dataset. This practical limitation may explain the somewhat poorer match to the data provided by the GCM of Parmentier et al. compared to that obtained by our retrieval analysis shown in Figure 8.

While we stress that GCMs represent the state-of-the-art for modelling the 3D interplay between circulation, radiation, and chemistry, by implementing a simpler model, our retrieval analysis is able to more fully explore its parameter space and optimize the match to the data. This comes at the costs of approximating the atmosphere as 1D, not solving for the PT profile self-consistently given the atmospheric opacity, and ignoring dynamical effects. For each MCMC sample, the retrieval instead simply takes a step in the six-dimensional parameter space and uses the resulting values for $\kappa_{\textnormal{IR}}$, $\gamma$, and $\psi$ to evaluate the PT profile, and the values of $[\textnormal{M}/\textnormal{H}]$, $[\textnormal{C}/\textnormal{H}]$, and $[\textnormal{O}/\textnormal{H}]$ to determine the elemental abundances. These two inputs — the PT profile and elemental abundances — are then used to solve for the chemical equilibrium abundances, with the effects of thermal dissociation and ionization accounted for appropriately.

Figure 11 shows the PT profile, contribution functions, and abundances obtained in this way from our retrieval analysis. Specifically, the results shown are for the best-fit model (i.e. the MCMC sample with the lowest $\chi^{2}$ value), with panel (a) also displaying the distribution of PT profiles across all MCMC samples. As in Evans et al. (2017), we find the PT profile exhibits an unambiguous thermal inversion. However, our new retrieval analysis puts the $\sim 2700$ K photosphere much deeper within the atmosphere than implied by the PT distribution presented in Evans et al. (2017), which is indicated by dashed lines in panel (a) of Figure 11. For the latter, $2700$ K coincides with pressures $\sim 10\,\mu$bar, versus $\sim 10$ mbar for the updated PT profile. This is a consequence of the different modeling assumptions made in the two studies, as detailed above. In particular, the Evans et al. (2017) analysis did not enforce chemical equilibrium, and consequently dialled the VO abundance up to unrealistically high values in order to fit the $1.25\mu\textnormal{m}$ bump in the measured spectrum. Since the VO abundance was assumed to be uniform across all pressure levels, this in turn increased the opacity throughout the entire atmosphere column, explaining why the photosphere occurs at such low pressures for that model. In the updated retrieval analysis, the abundances of important absorbers such as H2O, VO, and TiO are greatly reduced at low pressures due to thermal dissociation (Figure 11), moving the photosphere to higher pressures.

We find that the temperature increases from approximately $2500$ K to $2800$ K across the pressures probed by the data (i.e. $\sim 30$ to $5$ mbar for the best-fit model as seen in Figure 11), which is close to the range spanned by the brightness temperatures shown in Figure 6. The inference of a thermal inversion is driven by various spectral features in the data appearing in emission rather than absorption, notably: H- within the G102 passband; H2O within the G141 passband; and unresolved CO and H2O bands within the IRAC $4.5\,\mu\textnormal{m}$ channel (Figures 6 and 8). As others have highlighted (Arcangeli et al., 2018; Lothringer et al., 2018; Parmentier et al., 2018; Kreidberg et al., 2018; Mansfield et al., 2018), the H2O bands are significantly muted due to thermal dissociation in the upper layers of the atmosphere, owing to the intense irradiation by the nearby host star. This can be appreciated by examining panel (c) of Figure 11, which shows the effect of thermal dissociation and ionization for various species in our best-fit model. For H2O, the volume mixing ratio decreases by a factor of $\sim 300$ between pressures of 100 mbar and 1 mbar.

The contribution of different species to the atmospheric opacity is illustrated in Figure 12. It shows absorption cross-sections for radiatively active species weighted by the corresponding abundances indicated in Figure 11 at a pressure of 10 mbar, coincident with the near-infrared photosphere. According to the model, the primary opacity source is H*-* across the G102 and $z^{\prime}$ passbands, as well as the short-wavelength half of the G141 passband. At longer wavelengths, H2O dominates across the remainder of the G141 passband, as well as the $K_{s}$ and IRAC $3.6\mu\textnormal{m}$ passbands, while CO dominates within the IRAC $4.5\mu\textnormal{m}$ passband.

Although we do not detect any spectral features due to TiO or VO in our emission data, the presence of the muted H2O band in the G141 passband could be hinting at the presence of these strong optical absorbers. Both Lothringer et al. (2018) and Parmentier et al. (2018) find that when TiO and VO are excluded as opacity sources in models of ultrahot Jupiter atmospheres, the $1.4\,\mu\textnormal{m}$ H2O band is entirely absent from the emission spectrum. This is because TiO and VO absorb a significant amount of incident stellar radiation at pressure levels above the near-infrared photosphere, even when their abundances have been depleted by thermal dissociation. When TiO and VO are removed as opacity sources, stellar radiation is able to penetrate deeper into the atmosphere, raising temperatures at the near-infrared photosphere by over 100 K (Lothringer et al., 2018). This in turn should increase the thermal dissociation rates for H2O enough to completely nullify the $1.4\,\mu\textnormal{m}$ spectral band. Therefore, the fact that we observe a muted $1.4\,\mu\textnormal{m}$ H2O emission band is consistent with the presence of one or more optical absorbers, such as TiO and VO.

Direct observational confirmation of the optical absorber/s responsible for the thermal inversion on the dayside hemisphere of WASP-121b will be challenging. Evidence for significant optical absorption at the day-night limb has been uncovered in the transmission spectrum, which may be due to VO, although no strong evidence for TiO has been found (Evans et al., 2018). However, even if one or both of these species are present in significant quantities on the dayside, their broad emission bandheads are likely weakened by thermal dissociation and overlapping H- continuum opacity. One possible workaround could be to use high resolution spectroscopy with large aperture ground-based telescopes to doppler-resolve the narrow cores of the strongest TiO and VO lines just prior to eclipse, as these will be less affected (e.g. Nugroho et al., 2017). Using Figure 12 as a guide, TiO could supersede H*-* as the dominant opacity source regulating the emission for wavelengths shortward of $\sim 0.7\mu\textnormal{m}$. Measurements made by the Transiting Exoplanet Survey Satellite (TESS) – which observed WASP-121 throughout January 2019 – will also help inform this picture. For example, if TiO and/or VO are present, a deeper eclipse depth would be expected within the TESS passband compared to the G102 passband, as the former extends across the $\sim 0.6$-$0.95\mu\textnormal{m}$ wavelength range. Looking further ahead to the James Webb Space Telescope, the second order of the Near Infrared Imager and Slitless Spectrograph (NIRISS) single-object spectrograph (SOSS) provides wavelength coverage across $\sim 0.6$-$0.8\,\mu\textnormal{m}$, which encompasses a number of significant TiO and VO bands (Figure 12). At even shorter wavelengths, the high resolution spectroscopy approach could in principle be used to doppler-resolve emission lines due to metals such as iron and titanium (Hoeijmakers et al., 2018) or photochemical products such as SH (Zahnle et al., 2009; Evans et al., 2018). Although these latter species are not included in the model shown in Figure 12, they have strong absorption lines at wavelengths shortward of $\sim 0.5\,\mu\textnormal{m}$ and could potentially heat the upper atmosphere enough to produce the observed thermal inversion (Lothringer et al., 2018).

Finally, we note that the G102 data presented here provide the most direct evidence yet for H- emission in an exoplanet atmosphere. This is due to the clear departure from a blackbody at wavelengths shortward of $1.1\,\mu\textnormal{m}$ (Figure 8), whereas previous claims have instead relied on model-dependent interpretations of G141 spectra that are indistinguishable from blackbodies (Arcangeli et al., 2018; Kreidberg et al., 2018; Mansfield et al., 2018). The extensive ionization implied by this result for pressures below $\sim 100$ mbar (Figure 11) could have important implications for atmospheric dynamics and energy transfer, including increased day-night heat redistribution due to H2 recombination on the nightside hemisphere (Bell & Cowan, 2018; Komacek & Tan, 2018) and increased magnetic drag due to Lorentz forces (Lothringer et al., 2018). As Kreidberg et al. (2018) have shown for WASP-103b — an ultrahot Jupiter that is in many respects similar to WASP-121b — spectroscopic phase curves offer the most promising means of further constraining these fundamentally 3D phenomena.

6 Conclusion

We have presented new secondary eclipse observations for the ultrahot Jupiter WASP-121b acquired with the G102 grism of HST/WFC3. These data extend the wavelength coverage of the measured emission spectrum from $1.1\,\mu\textnormal{m}$ down to $0.8\,\mu\textnormal{m}$. We performed a retrieval analysis of the combined emission dataset, improving upon our previous efforts by incorporating the effects of thermal dissociation and ionization. We confirm the detection of a thermal inversion and our best-fit model indicates that the temperature profile increases from approximately $2500$ K to $2800$ K across the $\sim 30$ mbar to $5$ mbar pressure range. The spectrum is well explained by H- emission for wavelengths shortward of $1.1\,\mu\textnormal{m}$, a muted H2O emission band at $1.4\,\mu\textnormal{m}$, and overlapping CO and H2O emission bands within the $4.5\,\mu\textnormal{m}$ channel. Under the assumption of chemical equilibrium, we find the dayside atmospheric metallicity is likely enriched by at least a factor of a few relative to solar and uncover no evidence for anomalous carbon and oxygen abundances.

Acknowledgements

The authors are grateful to the anonymous referee for constructive feedback that improved the quality of this manuscript. Support for program GO-15135 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. JMG acknowledges funding from a Leverhulme Trust Research Project Grant and University of Exeter PhD Studentship. ALC is funded by a UK Science and Technology Facilities Council (STFC) studentship. BD acknowledges support from an STFC Consolidated Grant (ST/R000395/1). PT acknowledges support by the European Research Council under Grant Agreement ATMO 757858.