An optical transmission spectrum for the ultra-hot Jupiter WASP-121b measured with the Hubble Space Telescope

An optical transmission spectrum for the ultra-hot Jupiter WASP-121b measured with the Hubble Space Telescope

[    David K. Sing    Jayesh Goyal    Nikolay Nikolov    Mark S. Marley    Kevin Zahnle    Gregory W. Henry    Joanna K. Barstow    Munazza K. Alam    Jorge Sanz-Forcada    Tiffany Kataria    Nikole K. Lewis    Panayotis Lavvas    Gilda E. Ballester    Lotfi Ben-Jaffel    Sarah D. Blumenthal    Vincent Bourrier    Benjamin Drummond    Antonio García Muñoz    Mercedes López-Morales    Pascal Tremblin    David Ehrenreich    Hannah R. Wakeford    Lars A. Buchhave    Alain Lecavelier des Etangs    Éric Hébrard    Michael H. Williamson
\par\LTX@newpageABSTRACT

We present an atmospheric transmission spectrum for the ultra-hot Jupiter WASP-121b, measured using the Space Telescope Imaging Spectrograph (STIS) onboard the Hubble Space Telescope (HST). Across the wavelength range, the data imply an atmospheric opacity comparable to – and in some spectroscopic channels exceeding – that previously measured at near-infrared wavelengths (). Wavelength-dependent variations in the opacity rule out a gray cloud deck at a confidence level of and may instead be explained by VO spectral bands. We find a cloud-free model assuming chemical equilibrium for a temperature of  K and metal enrichment of solar matches these data well. Using a free-chemistry retrieval analysis, we estimate a VO abundance of  dex. We find no evidence for TiO and place a upper limit of  dex on its abundance, suggesting TiO may have condensed from the gas phase at the day-night limb. The opacity rises steeply at the shortest wavelengths, increasing by approximately five pressure scale heights from to in wavelength. If this feature is caused by Rayleigh scattering due to uniformly-distributed aerosols, it would imply an unphysically high temperature of  K. One alternative explanation for the short-wavelength rise is absorption due to SH (mercapto radical), which has been predicted as an important product of non-equilibrium chemistry in hot Jupiter atmospheres. Irrespective of the identity of the NUV absorber, it likely captures a significant amount of incident stellar radiation at low pressures, thus playing a significant role in the overall energy budget, thermal structure, and circulation of the atmosphere.

Corresponding author: Thomas M. Evanstmevans@mit.edu

0000-0001-5442-1300]Thomas M. Evans \move@AU\move@AF\@affiliationPhysics and Astronomy, Stocker Road, University of Exeter, Exeter, EX4 3RF, UK \move@AU\move@AF\@affiliationKavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Avenue, 37-241, Cambridge, MA 02139, USA

\move@AU\move@AF\@affiliation

Physics and Astronomy, Stocker Road, University of Exeter, Exeter, EX4 3RF, UK \move@AU\move@AF\@affiliationDepartment of Earth and Planetary Sciences, Johns Hopkins University, Baltimore, MD, USA

\move@AU\move@AF\@affiliation

Physics and Astronomy, Stocker Road, University of Exeter, Exeter, EX4 3RF, UK

\move@AU\move@AF\@affiliation

Physics and Astronomy, Stocker Road, University of Exeter, Exeter, EX4 3RF, UK

\move@AU\move@AF\@affiliation

NASA Ames Research Center, Moffett Field, California, USA

\move@AU\move@AF\@affiliation

NASA Ames Research Center, Moffett Field, California, USA

\move@AU\move@AF\@affiliation

Center of Excellence in Information Systems, Tennessee State University, Nashville, TN 37209, USA

\move@AU\move@AF\@affiliation

Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

\move@AU\move@AF\@affiliation

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

\move@AU\move@AF\@affiliation

Centro de Astrobiologıía (CSIC-INTA), ESAC Campus, Camino Bajo del Castillo, E-28692 Villanueva de la Canada, Madrid, Spain

\move@AU\move@AF\@affiliation

NASA Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

\move@AU\move@AF\@affiliation

Department of Astronomy and Carl Sagan Institute, Cornell University, 122 Sciences Drive, 14853, Ithaca, NY, USA

\move@AU\move@AF\@affiliation

Groupe de Spectroscopie Moléculaire et Atmosphérique, Université de Reims, Champagne-Ardenne, CNRS UMR F-7331, France

\move@AU\move@AF\@affiliation

Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA

\move@AU\move@AF\@affiliation

Sorbonne Universités, UPMC Université Paris 6 and CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98 bis boulevard Arago, F-75014 Paris, France

\move@AU\move@AF\@affiliation

Physics and Astronomy, Stocker Road, University of Exeter, Exeter, EX4 3RF, UK

\move@AU\move@AF\@affiliation

Observatoire de l’Université de Genève, 51 chemim des Maillettes, 1290 Sauverny, Switzerland

\move@AU\move@AF\@affiliation

Physics and Astronomy, Stocker Road, University of Exeter, Exeter, EX4 3RF, UK

\move@AU\move@AF\@affiliation

Zentrum für Astronomie und Astrophysik, Technische Universität Berlin, Hardenbergstrasse 36, D-10623 Berlin, Germany

\move@AU\move@AF\@affiliation

Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

\move@AU\move@AF\@affiliation

Maison de la Simulation, CEA, CNRS, Univ Paris-Sud, UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

\move@AU\move@AF\@affiliation

Observatoire de l’Université de Genève, 51 chemim des Maillettes, 1290 Sauverny, Switzerland

\move@AU\move@AF\@affiliation

Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, Maryland 21218, USA

\move@AU\move@AF\@affiliation

DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 328, DK-2800 Kgs. Lyngby, Denmark

\move@AU\move@AF\@affiliation

Sorbonne Universités, UPMC Université Paris 6 and CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98 bis boulevard Arago, F-75014 Paris, France

\move@AU\move@AF\@affiliation

Physics and Astronomy, Stocker Road, University of Exeter, Exeter, EX4 3RF, UK

\move@AU\move@AF\@affiliation

Center of Excellence in Information Systems, Tennessee State University, Nashville, TN 37209, USA

1 Introduction

Spectroscopic observations made during the primary transit of an exoplanet allow the atmospheric transmission spectrum of the day-night boundary region to be probed (2000ApJ...537..916S), while the same type of observation made during secondary eclipse provides the emission spectrum of the dayside hemisphere (1998ApJ...502L.157S). Much of the transmission and emission spectroscopy work published to date has employed the Hubble Space Telescope (HST), primarily with the Space Telescope Imaging Spectrograph (STIS), covering the UV-optical wavelength range, and Wide Field Camera 3 (WFC3), covering the near-IR wavelength range.

A non-exhaustive list of HST transmission spectroscopy highlights at optical through IR wavelengths include: the detection of Na on HD 209458b (2002ApJ...568..377C); multiple detections of HO (e.g. 2013ApJ...774...95D; 2016ApJ...822L...4E; 2014Natur.513..526F; 2013MNRAS.434.3252H; 2015ApJ...814...66K; 2018AJ....155..156T; 2017Sci...356..628W; 2018AJ....155...29W); widespread evidence for aerosols (e.g. 2014Natur.505...69K; 2014MNRAS.437...46N; 2015MNRAS.447..463N; 2008MNRAS.385..109P; 2015MNRAS.446.2428S; 2016Natur.529...59S); and a detection of He in the extended atmosphere of WASP-107b (2018Natur.557...68S). At UV wavelengths, transit observations made with STIS have probed the hydrogen exospheres of hot Jupiters (e.g. 2003Natur.422..143V) and warm Neptunes (e.g. 2015Natur.522..459E), while heavier elements such as oxygen have been detected using the HST Cosmic Origins Spectrograph (e.g. 2010ApJ...714L.222F; 2013A&A...553A..52B). For emission, a similar list includes: detections of HO absorption (2017AJ....154..158B; 2014Sci...346..838S); evidence for HO emission (2017Natur.548...58E); evidence for TiO emission (2015ApJ...806..146H); constraints on optical reflection spectra (2013ApJ...772L..16E; 2017ApJ...847L...2B); and multiple featureless thermal spectra (e.g. 2018MNRAS.474.1705N; 2018AJ....156...10M).

This paper reports a transmission spectrum measured for the ultra-hot ( K) Jupiter WASP-121b across the wavelength range using STIS. Discovered by 2016MNRAS.tmp..312D, WASP-121b orbits a moderately bright () F6V host star, which has an estimated radius of (2016MNRAS.tmp..312D) and measured parallax of  mas (2018A&A...616A...1G), corresponding to a system distance of  parsec. WASP-121b itself has a mass of , an inflated radius of , and a dayside equilibrium temperature above  K. Together, these properties make WASP-121b an excellent target for atmospheric characterization (2016MNRAS.tmp..312D; 2016ApJ...822L...4E; 2017Natur.548...58E).

We previously published the near-IR transmission spectrum for WASP-121b measured using WFC3 in 2016ApJ...822L...4E. Those data revealed absorption due to the HO band centered at , along with a second bump across the wavelength range, which we suggested could be a signature of FeH or VO. Analyzing the same dataset, 2018AJ....155..156T reproduced the feature and presented a best-fit model including absorption by TiO and VO, although they did not discuss FeH. In 2016ApJ...822L...4E, we also compared the WFC3 transmission spectrum with transits measured at optical wavelengths by 2016MNRAS.tmp..312D using ground-based photometry. This comparison implied significantly deeper transits at optical wavelengths relative to the near-IR, which we speculated could be evidence for a strong opacity source such as TiO and/or VO. Subsequent modeling of these data confirmed such an interpretation to be plausible (e.g. 2017ApJ...845L..20K; 2018A&A...617A.110P).

In 2017Natur.548...58E, we presented a secondary eclipse observation for WASP-121b, also made with WFC3 at near-IR wavelengths. The measured spectrum indicates a mean photosphere temperature of approximately  K and shows the HO band in emission, rather than absorption, implying the dayside hemisphere has a vertical thermal inversion. As for the transmission spectrum, the emission data exhibit a second bump across the wavelength range, which can be fit with VO in emission. To do so, however, requires assuming a VO abundance over higher than expected for solar elemental composition in chemical equilibrium, casting doubt on this interpretation. Models assuming chemical equilibrium and abundances closer to solar do not reproduce the bump (e.g. 2018A&A...617A.110P). For now, we do not have a satisfying explanation for this feature, but the fact that it has been observed in both the transmission spectrum and emission spectrum is intriguing.

Our understanding of the atmosphere of WASP-121b remains a work in progress. For instance, the thermal inversion measured for the dayside hemisphere implies significant heating at low pressures ( mbar), though it is unclear what causes this. One possibility is absorption of incident stellar radiation at optical wavelengths by TiO and VO (e.g. 2003ApJ...594.1011H; 2008ApJ...678.1419F). However, neither of these species have yet been definitively detected in the atmosphere of WASP-121b, despite the hints described above. Furthermore, it has been pointed out that TiO and VO could be removed from the upper atmospheres of even very hot planets by cold-trapping (e.g. 2009ApJ...699.1487S; 2009ApJ...699..564S; 2017AJ....154..158B). Additionally, the dayside temperatures of ultra-hot Jupiters such as WASP-121b are likely high enough for significant thermal dissociation of TiO and VO, along with other molecules such as HO, to occur (2018ApJ...855L..30A; 2018AJ....156...17K; 2018ApJ...866...27L; 2018A&A...617A.110P). Nonetheless, evidence for TiO has been detected on the dayside of WASP-33b (2015ApJ...806..146H; 2017AJ....154..221N), which has a mean photosphere temperature of around  K at near-IR wavelengths, making it even hotter than WASP-121b. An optical transmission spectrum measured for another ultra-hot Jupiter, WASP-19b, also exhibits a prominent TiO band (2017Natur.549..238S), although this may have been the signature of unocculted star spots (2018MNRAS.tmp.2569E). Despite the picture remaining unclear, observations such as these imply TiO, and presumably VO, can perhaps persist at low pressures in ultra-hot Jupiter atmospheres. As will be described in the following sections, the STIS transmission spectrum for WASP-121b provides new evidence for VO absorption at optical wavelengths.

Absorption at UV wavelengths may also play a significant role in heating the upper atmospheres of strongly-irradiated planets such as WASP-121b. For instance, 2009ApJ...701L..20Z examined non-equilibrium sulfur chemistry in the context of hot Jupiter atmospheres and concluded that SH and S could be important absorbers across the wavelength range. These species may be driven to higher-than-equilibrium abundances via reactions involving the photolytic and photochemical destruction of HS. As will be reported below, the measured transmission spectrum for WASP-121b exhibits a strong signal at wavelengths shortward of and absorption by SH appears to provide a viable explanation.

We begin, however, by describing our observations and the steps taken to extract the spectra from the raw data frames in Section 2. We present analyses of the white lightcurves in Section 3 and spectroscopic lightcurves in Section 4. The results are discussed in Section 5, including the implications of the measured transmission for the planetary atmosphere. Our conclusions are given in Section 6.

2 Observations and data reduction

We observed three primary transits of WASP-121b using HST/STIS as part of the Panchromatic Comparative Exoplanet Treasury (PanCET) survey (Program 14767; P.I.s Sing and López-Morales). This was comprised of two visits made on 2016 Oct 24 and 2016 Nov 6 with the G430L grating, and one visit made on 2016 Nov 12 with the G750L grating. In what follows, we shall refer to the first and second G430L visits as the G430Lv1 and G430Lv2 datasets, respectively. For all three STIS visits, the target was observed for 6.8 hours, covering five consecutive HST orbits. Observations were made using the widest available slit ( arcsec) to minimize slit losses and the detector gain was set to 4 /DN. Overheads were reduced by only reading out a pixel subarray containing the target spectrum. Exposure times of 253 s and 161 s were used for the G430L and G750L observations, respectively. We also took a short 1 s exposure at the start of each HST orbit for both gratings, but discarded these exposures in the subsequent analysis. This was done because STIS observations typically suffer from a systematic in which the first exposure of each HST orbit has anomalously lower counts relative to the immediately-following exposures (e.g. 2013ApJ...772L..16E; 2014MNRAS.437...46N; 2015MNRAS.447..463N; 2015MNRAS.446.2428S) and we wanted to minimize the integration time lost to this effect. With this observing setup, we acquired a total of 48 science exposures for each G430L visit and 70 science exposures for the G750L visit.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefExample spectra for the G430L and G750L gratings (top panel), and the G141 grism (bottom panel). Dark and light vertical bands indicate the wavelength channels adopted for the spectroscopic lightcurves.

The STIS datasets were reduced following the methodology described in 2014MNRAS.437...46N; 2015MNRAS.447..463N. Raw data frames were bias-, dark-, and flat-corrected using the CALSTIS pipeline (v3.4) with relevant calibration frames. Cosmic ray events and pixels flagged as ‘bad’ by CALSTIS were removed and interpolated over. Overall, we found % of pixels were affected by cosmic rays for all visits with a further % flagged as bad by CALSTIS. To extract spectra from the cleaned 2D frames, we used the IRAF procedure apall with aperture radii of 4.5, 6.5, 8.5, and 10.5 pixels for both the G430L and G750L datasets. The dispersion axis was mapped to a wavelength solution using the x1d files produced by CALSTIS.

In addition to the STIS data, a single primary transit of WASP-121b was observed on 2016 Feb 6 with the G141 grism (Program 14468; P.I. Evans). This dataset was originally published in 2016ApJ...822L...4E, to which the reader is referred for further details.

Example G430L, G750L, and G141 spectra are shown in Figure 2.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHref (Top row) Raw white lightcurves for the G430Lv1, G430Lv2, and G750L datasets. Gray lines show the best-fit transit signals with linear baseline trends. (Middle row) Dispersion drift variable for each dataset. (Bottom row) Cross-dispersion drift variable for each dataset. In all panels, colored symbols indicate data points that were included in the analysis and gray crosses indicate those that were excluded for reasons explained in the main text. The two drift variables are unitless as they have been standardized, i.e. mean subtracted and normalized by their standard deviations.

3 White Lightcurve Analyses

White lightcurves were constructed for each dataset by summing the flux of each spectrum across the full dispersion axis. The resulting lightcurves are shown in the top row of Figure 2. As in our previous work (2013ApJ...772L..16E; 2016ApJ...822L...4E; 2017Natur.548...58E), we followed the methodology outlined by 2012MNRAS.419.2683G and treated each lightcurve as a Gaussian process (GP). Under this approach, the posterior likelihood is described by a multivariate normal distribution of the form , where: is an -length vector containing the flux measurements; is a vector containing the deterministic mean function; is an matrix describing the correlations between data points; and is an diagonal matrix containing the squared white noise uncertainties, , for each data point .

For the mean function, we adopted a 2002ApJ...580L.171M transit model multiplied by a linear trend in time () of the form . We assumed a circular orbit with a period () of 1.2749255 days (2016MNRAS.tmp..312D). We allowed the normalized planet radius () and transit mid-time () to vary as free parameters with uniform priors. As described in Section 3.1, we first performed fits with the normalized semimajor axis () and impact parameter () allowed to vary as free parameters, both with uniform priors. Then, as described in Section 3.2, we fixed and to their weighted-mean values and repeated the fitting.

In all fits, we assumed a quadratic limb darkening law and treated both coefficients (, ) as free parameters. We first estimated values for and by fitting to the limb darkening profile of a stellar model over the appropriate bandpass. Specificially, we used a 3D stellar model from the STAGGER grid (2013A&A...557A..26M) with  K,  cgs, and  dex, as this was the grid point closest to the properties of the WASP-121 host star ( K,  cgs,  dex; 2016MNRAS.tmp..312D). We then applied broad normal priors to and in the model fitting, with means set to these estimated values and standard deviations of 0.6, providing plenty of flexibility for the model to be optimized.

For the GP covariance matrix , we adopted a squared-exponential kernel333We refer the reader to previous studies such as 2012MNRAS.419.2683G, 2013ApJ...772L..16E, and 2014MNRAS.445.3401G for further details of the squared-exponential kernel. with three input variables that it is reasonable to assume could correlate with the instrumental systematics: namely, HST orbital phase (), dispersion drift (), and cross-dispersion drift (). This resulted in four free parameters for each dataset: namely, the covariance amplitude () and correlation length-scales () for each input variable, . For the white noise matrix, , we adopted the formal photon noise values multiplied by a rescaling factor () which was allowed to vary as a free parameter. The latter affords some flexibility to handle high-frequency systematics that are pseudo-white-noise in nature, which would otherwise bias the model toward impractically small values.

For the GP covariance amplitude , we adopted Gamma priors of the form , to favor smaller correlation amplitudes. This can help prevent a small number of outliers having a disproportionate influence on the inferred covariance amplitude. For the correlation length scales , we followed previous studies (e.g. 2017Natur.548...58E; 2017MNRAS.467.4591G) and fit for the natural logarithm of the inverse correlation length scales , adopting uniform priors for each. In practice, this favors longer correlation length scales, with the intention of capturing the lower-frequency systematics present in the data, as these are most degenerate with the planet signal. Higher-frequency systematics can be accounted for through the parameter, for which we adopted a normal prior with mean of 1 and standard deviation of 0.2, to favor values close to the formal photon noise.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefPosterior distributions for , , and obtained from the white lightcurve analyses described in Section 3.1. Blue, pink, and orange regions indicate smoothed contours containing 68% of MCMC samples for the G430L, G750L, and G141 analyses, respectively. Green region indicates the weighted-mean of the HST posterior distributions and gray region indicates the range reported by 2016MNRAS.tmp..312D.

We modeled the white lightcurves for G430L, G750L, and G141 separately. For the G430L lightcurves, we assumed , , and were the same for both visits, while allowing , , , , , and to vary separately for each visit. The posterior distributions were marginalized using affine-invariant Markov chain Monte Carlo (MCMC), as implemented by the emcee Python package (2013PASP..125..306F). In all fits, we randomly distributed five groups of 150 walkers throughout the parameter space and allowed them to run for 100 steps to locate the peak of the posterior distribution. We then re-initialized the five groups of 150 walkers in a tighter ball around this peak and allowed them to run for 500 steps, of which we discarded the first 250 steps as burn-in and combined the remaining 250 steps into a single chain for each walker group. At this point, a comparison of the chains from each walker group confirmed that they appeared well-mixed and converged, with Gelman-Rubin statistic values within 2% of unity for each free parameter (GelmanRubin92). Table 6 summarizes the resulting posterior distributions. For each of the STIS lightcurves produced using the different trial apertures (see Section 2), we obtained results consistent to within for the planet parameters (e.g. ) and report only those for the  pixel aperture.

3.1 and allowed to vary

The purpose of the model fits in which and were allowed to vary as free parameters was to use the HST data to refine our estimates of these system properties. Previously, the only published measurements were those provided in the original discovery paper by 2016MNRAS.tmp..312D, which reported and . Figure 3 shows the posterior distributions obtained from our analyses for comparison, with values reported in Table 6. We find good agreement for both and across our fits to the G430L, G750L, and G141 white lightcurve datasets. Taking the arithmetic weighted-mean of these results, we estimate and , implying  deg. We note that our HST results differ from those of Delrez et al. by for and for . The reason for this disagreement is unclear and will likely be resolved by additional transit observations that are currently planned or in the process of being analyzed (Evans et al., in prep.). For the present study, we note that the primary consequence of assuming slightly different values for and will be to perturb the inferred values for . Importantly, this will be a wavelength-independent effect and thus should not affect our interpretation of the atmospheric transmission spectrum. For this reason, and given the mutual agreement between the G430L, G750L, and G141 datasets, we adopt the HST weighted-mean values for and in all subsequent lightcurve fits.

3.2 and held fixed

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefWhite lightcurves for G430L, G750L, and G141 datasets analyzed in this study. (Top row) Relative flux variation after removing the systematics contribution inferred from the GP analyses (see Figure 3.2), with best-fit transit signals plotted as solid lines. (Middle row) Corresponding model residuals, with photon noise errorbars. (Bottom row) Normalized histograms of residuals obtained by subtracting from the data a random subset of GP mean functions obtained in the MCMC sampling. Solid black lines correspond to normal distributions with standard deviations equal to photon noise (i.e. prior to rescaling by the factor described in the main text).

Inferred values for can be biased by differences in the assumed values for and across datasets. For this reason, we held the latter parameters fixed to the HST weighted-mean values determined in the previous section and repeated the white lightcurve analyses. This is physically motivated by the fact that the true values of and should be constant across our datasets, and we are primarily interested in wavelength-dependent variations of arising due to the planetary atmosphere.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSystematics in the white lightcurves for G430L and G750L datasets. Effectively, these are the residuals after dividing the raw flux time series by the transit signals with linear baseline trends shown in Figure 2. Yellow lines and gray shaded regions, respectively, show the means and , , and ranges of the best-fit GP distributions. Note that in practice the transit signal, linear baseline trend, and GP are fit simultaneously. The purpose of this figure is only to highlight the structure of the systematics.

Figure 3.2 shows the best-fit transit models compared with the data after removing the systematics contribution inferred by the GP. The latter are shown separately in Figure 3.2 and Table 6 summarizes the posterior distributions. The resulting estimates for , , and are all within of those obtained for the fits in which and were allowed to vary. Unsurprisingly, we obtain similar estimates for , as this parameter is sensitive to high-frequency noise in the data that is unlikely to be significantly correlated with and . The inferred values imply scatters that are % and % above the photon noise floor for the STIS and WFC3 datasets, respectively. This is illustrated in Figure 3.2, which shows the model residuals. For , we find the inferred values shift by  sec, but remain within of those obtained for the fits in which and were allowed to vary.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefPosterior distributions obtained from the white lightcurve analyses described in Section 3.2. Top right panels show results for the joint analysis of both G430L visits and bottom left panels show results for the G750L analysis. Plotted contours contain 68% and 95% of the MCMC samples. Panels along the diagonal show marginalized posterior distributions. Note that , , , and have been median-subtracted to allow both G430L visits to be plotted on the same axes. The purpose of this figure is to visually illustrate correlations between model parameters. Numerical values for all parameter distributions are summarized in Table 6.

4 Spectroscopic Lightcurve Analyses

Spectroscopic lightcurves were constructed by first summing the spectra of each dataset within the wavelength channels shown in Figure 2. Median channel widths were: 20 pixels ( Å) for both G430L datasets; 20 pixels ( Å) for the G750L dataset; and 4 pixels ( Å) for the G141 dataset. Care was taken to avoid the edges of prominent stellar lines and to maintain similar levels of flux within each channel. Thus, subsets of the G430L and G750L channels were broader than these nominal widths. The resulting raw lightcurves for the STIS datasets are shown in Figures A-A.

We next generated common-mode (i.e. wavelength-independent) signals for each dataset by dividing the raw white lightcurves by the corresponding best-fit transit signals obtained in Section 3 and shown in Figure 3.2. Each of the raw spectroscopic lightcurves were then divided by the resulting common-mode signals. Note that in addition to removing common-mode systematics, this latter step also has the effect of dividing each spectroscopic lightcurve by the intrinsic scatter of the white lightcurve. However, this is acceptable, as the spectroscopic lightcurves have a larger intrinsic scatter than the white lightcurves: dividing white noise by lower-amplitude white noise should on average have zero net effect on the scatter of the resulting corrected lightcurves. Meanwhile, applying a common-mode correction of this nature – as opposed to dividing through by the best-fit systematics model from the white lightcurve fits – has the potential advantage of removing systematics in the spectroscopic lightcurves that may not be captured by our white lightcurve systematics model. The common-mode corrected lightcurves for the STIS datasets are shown in Figures A-A.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSpectroscopic lightcurves for the G430Lv1 and G430Lv2 datasets after removing the systematics contributions inferred from the GP analyses, with best-fit transit signals plotted as solid lines. Green triangles and purple diamonds correspond to the G430Lv1 and G430Lv2 datasets, respectively.

To fit the spectroscopic lightcurves, we used the same approach as described in Section 3. The only exception was that we fixed to the best-fit values listed in Table 6. Thus, for the spectroscopic transit signals, the free paramaters were the radius ratio () and quadratic limb darkening coefficients (, ). For the G430L analysis, we fit both visits jointly with shared values for , , and , as was done for the white lightcurve analysis. In all fits, we again accounted for systematics by fitting for a linear trend in and a GP with as inputs to a squared-exponential covariance kernel. White noise levels were allowed to vary for each individual lightcurve via rescaling parameters. Marginalization of the posterior distributions was performed in the manner described above, using affine-invariant MCMC.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSimilar to Figure 4, but for the G750L spectroscopic lightcurves.

The best-fit transit signals and model residuals are shown in Figure 4 for G430L and Figure 4 for G750L. Figure 4 shows the systematics and GP fits for each spectroscopic lightcurve. Histograms of residuals are shown in Figures A-A. For the G141 spectroscopic lightcurve fits, the results were essentially identical to those presented in 2016ApJ...822L...4E, so we do not duplicate them here. The only difference for the latter is a wavelength-uniform shift of by , in line with the revised white lightcurve analysis which gives (Table 6), compared with the previous estimate of (2016ApJ...822L...4E).444The revised value for within the G141 bandpass can be attribued to the updated values for and adopted in the present study (Section 3.1).

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSimilar to Figure 3.2, but showing the systematics and GP fits for the spectroscopic lightcurves. In all columns, wavelength increases from top to bottom.

As shown in Figure 4, we obtain means and standard deviations for the inferred values across spectroscopic channels of: for the G430lv1 dataset; for the G430Lv2 dataset; for the G750L dataset; and for the G141 dataset. The consistency of these results with indicate that the GP models are broadly successful at marginalizing over the correlations in the lightcurves, implying in turn that degeneracies between the systematics and planet signal are properly accounted for in our estimates of parameters such as , which we are primarily interested in.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHref(Top row) Inferred white noise rescaling parameters for the GP analyses adopting a squared-exponential covariance kernel. (Bottom row) The same, but for the GP analyses adopting a Matérn covariance kernel.

To investigate the sensitivity of our results to the choice of covariance kernel, we repeated the spectroscopic lightcurve fitting using the Matérn kernel, which can be more suitable for modeling high-frequency signals than the squared-exponential kernel (e.g. see 2012MNRAS.419.2683G). For all channels, we found the inferred values remained unchanged to well-within , regardless of which covariance kernel was used. However, the values inferred using the Matérn kernel were on average slightly closer to unity, as illustrated in Figure 4. This suggests some of the channels may contain high-frequency noise that can be suitably accounted for either by inflating the white-noise level above the photon noise floor via or by employing a covariance kernel with enough flexibility to marginalize over signals of this nature, such as the Matérn . Given the results for are found to be insensitive to the choice of covariance kernel, we adopt those obtained using the squared-exponential for the remainder of this paper.

The corresponding posterior distributions for , , and are summarized for each STIS dataset in Tables 6 and 6. The median uncertainties on are 800 ppm for G430L, 900 ppm for G750L, and 500 ppm for G141, which translates to uncertainties on the transit depth of approximately 200 ppm, 220 ppm and 125 ppm, respectively. For comparison, a change in the effective planetary radius of one atmospheric pressure scale height corresponds to a transit depth variation of  ppm for WASP-121b, assuming average limb temperatures in the range of  K, a planetary surface gravity of 940 cm s, and an atmospheric mean molecular weight of atomic mass units (i.e. equal to that of Jupiter).

5 Discussion

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefTransmission spectrum for WASP-121b obtained using STIS and WFC3 (colored circles) and ground-based photometry from 2016MNRAS.tmp..312D (unfilled squares). Note that the latter are taken from the re-analysis of 2016ApJ...822L...4E, although very similar results were obtained by Delrez et al. Light blue halos indicate the subset of G430L data that we refer to as the blue data in the main text. Two forward models assuming chemical equilibrium are also shown, both with a temperature of  K and solar metallicity. One model includes TiO/VO opacity (light purple line) and the other does not include TiO/VO opacity (dark purple line).

The measured transmission spectrum is shown in Figure 5 and has a number of notable features. In particular, the G430L data exhibit a steep rise toward shorter wavelengths from , where , to , where . This corresponds to a change in effective planetary radius of approximately five pressure scale heights. At longer optical wavelengths covered by the G430L and G750L gratings (), is measured to vary across spectroscopic channels, implying a wavelength-dependent atmospheric opacity. Within some of these optical channels, the atmospheric opacity is found to be even higher than that measured within the HO band at , which is detected in the G141 bandpass (2016ApJ...822L...4E).

Figure 5 also shows the values measured using ground-based photometry in the , , and bandpasses. The latter were originally reported by 2016MNRAS.tmp..312D and an independent analysis of the same lightcurves was presented by 2016ApJ...822L...4E. Both studies obtained similar estimates for in each bandpass that are larger than those obtained from the HST data. It is unclear what is responsible for this tension. As noted above, updated values for and were used in the lightcurve fits of the present study. However, for the G141 dataset, this had the effect of shifting the mean value to a higher value from (2016ApJ...822L...4E) to (Table 6). A similar upward shift for the ground-based photometry would make those data more discrepant relative to the HST data. Alternatively, during the photometry data reduction, effects such as aperture light-losses or an over-estimated background may have artificially deepened the transit signals, resulting in estimates above the true values. Another more speculative possibility is intrinsic variability of the atmosphere from epoch to epoch. For example, 2013A&A...558A..91P report a 3D GCM study showing that significant variations in passive tracer abundances over day timescales are possible at the planetary limb of hot Jupiters. As those authors note, if this occurs for strongly-absorbing species such as TiO and VO, it could have significant implications for transmission spectra measured at different epochs. Indeed, the ground-based photometry and HST STIS observations were separated by over days. However, the current data are insufficient to test this theory, and we consider it more likely that the difference is due to some unaccounted-for systematic in the ground-based photometry.

To evaluate the robustness of the HST transmission spectrum, we performed a number of tests, full details of which are reported in Appendix B. First, we find that the measured transmission spectrum is insensitive to our treatment of limb darkening. Second, we investigated the inclusion of time as an additional GP input variable in the lightcurve fits and obtain very similar results to those reported here. Third, for the G430L data, we find that the measured transmission spectrum is repeatable when each of the two visits are analyzed separately. Fourth, we conclude that stellar activity is unlikely to have significantly affected the measured transmission spectrum, based on: (1) the lack of photometric variability and modest X-ray flux of the WASP-121 host star; (2) the epoch-to-epoch repeatability of the G430L datasets; (3) the good level of agreement obtained across the overlapping wavelength range of the G430L and G750L datasets; and (4) the inability of unocculted spots to explain the shape of the measured spectrum under reasonable assumptions. In the following sections, we therefore seek to interpret the measurements shown in Figure 5 as the signal of the planetary atmosphere.

5.1 Rayleigh scattering and a gray cloud-deck

The signature of aerosol scattering is ubiquitous in observations of exoplanet atmospheres (e.g. 2008MNRAS.385..109P; 2014Natur.505...69K; 2014MNRAS.437...46N; 2015MNRAS.447..463N; 2015MNRAS.446.2428S). For hot Jupiter transmission spectra, this is unsurprising given the large number of refractory species expected to condense at the temperatures and pressures characteristic of these atmospheres (e.g. 2018A&A...614A...1W), as well as the highly-sensitive nature of the grazing geometry to even trace opacity sources (2005MNRAS.364..649F). Indeed, the rise in opacity toward shorter wavelengths that we measure for WASP-121b is somewhat reminiscent of transmission spectra previously obtained for other hot Jupiters, which can be explained by Rayleigh scattering due to high-altitude layers of submicron aerosols (2016Natur.529...59S). In addition, an optically-thick cloud deck could act as a gray opacity source, if present at low pressures.

We investigated how well the WASP-121b transmission spectrum can be explained by aerosols by first fitting simple Rayleigh scattering and cloud deck models to the STIS data spanning the G430L and G750L gratings. For the Rayleigh component, we followed the methodology outlined by 2008A&A...481L..83L (L08), who provide relations between the slope of the transmission spectrum and the atmospheric temperature, under the assumption of scattering particles distributed uniformly with pressure. For the cloud deck component we assumed a wavelength-independent opacity, implemented as a horizontal flat line in that was allowed to float vertically relative to other spectral features in the transmission spectrum. For this initial analysis, we excluded the G141 dataset, as it exhibits a clear spectral feature due to HO, which would add additional complexity to the model. This is addressed in Section 5.4, where we perform a free-chemistry fit to the combined STIS+WFC3 dataset that includes opacity due to both gas-phase species and aerosols.

Our best-fit model combining a Rayleigh slope with cloud-deck is shown in Figure 5.1. It provides a poor fit to the data, with a reduced of 1.8 for 57 degrees of freedom, allowing us to exclude it at confidence. This is due to the inability of a featureless cloud deck to explain the optical data across the wavelength range. Furthermore, the temperature inferred from the Rayleigh slope is  K, which is improbably high for the atmospheric pressures probed in transmission. For instance, if WASP-121b absorbs all incident radiation on its dayside hemisphere (i.e. the Bond albedo is zero), then the substellar point would have a temperature of  K and the day-night boundary probed by the transmission spectrum should be considerably cooler. Furthermore, at such high temperatures, no condensates are expected to exist and molecules should be thermally dissociated, including H.

To be conservative, we also tried dividing the NUV-optical data into different wavelength sections and fitting them one at a time. For convenience, we will refer to these subsets as the blue () and red () data. In principle, a good fit to one or both of these datasets separately should be easier to achieve than a good joint fit, as the models need not be self-consistent.

First, we fit a Rayleigh profile to the blue data, as this is where the transmission spectrum exhibits a strong slope. Although we obtain a better statistical fit with a reduced of 1.6 for 10 degrees of freedom (Figure 5.1), the inferred temperature remains implausibly high at  K. Given this, we conclude that the rise in the measured transmission spectrum toward NUV wavelengths is too steep to be explained by scattering out of the transmission beam. Instead, it would suggest the presence of one or more significant NUV absorbers in the upper atmosphere of WASP-121b, assuming the slope is indeed a feature of the planetary spectrum and not caused by an uncorrected systematic effect in the data.

Second, we tried fitting a gray cloud deck to the G430L and G750L data, with the blue G430L subset excluded. For this scenario (not shown in Figure 5.1), we obtain a reduced of 1.9 for 51 degrees of freedom, which formally rules it out at confidence. Alternatively, if the uncertainties for these optical data have been uniformly underestimated by %, this gray cloud scenario would only be excluded at confidence. However, lacking any reason to doubt our inferred uncertainties, we propose instead that the optical data exhibit significant spectral variations that cannot be explained by a gray cloud deck.

5.2 Forward model comparison with optical-NIR data

The results of the previous section imply the transmission spectrum of WASP-121b exhibits significant wavelength-dependent opacity variations across the wavelength range. To explore this further, we used the ATMO code (2014A&A...564A..59A; 2016A&A...594A..69D; 2018MNRAS.474.5158G; 2015ApJ...804L..17T; 2016ApJ...817L..19T) to generate a small grid of aerosol-free atmosphere models spanning temperature and metallicity, assuming isothermal pressure-temperature (PT) profiles and chemical equilibrium abundances. Specifically, our grid consisted of temperatures ranging from  K to  K in  K increments, each evaluated for metallicities of , , , , , , and solar. ATMO solves for the gas-phase and condensed-phase chemical equilibrium mole fractions for a given pressure, temperature, and set of elemental abundances (2016A&A...594A..69D). For the results presented here we consider local condensation, such that the chemistry calculation in each model pressure level is entirely independent of all other pressure levels. We do not account for rainout chemistry, under which condensation deeper within the atmosphere depletes elemental abundances at lower pressures levels (1999ApJ...512..843B; 2011ApJ...737...34M; 2016ApJ...827..121M). Rainout could be important in the atmosphere of WASP-121b, but we defer investigation of this effect to future work that includes a more realistic treatment of the PT profile than the isothermal assumption made here. Finally, we applied uniform vertical offsets to for each model in order to optimize the match to the data. No further tuning of the models was performed.

None of these equilibrium models are able to explain the absorption at wavelengths shortward of , nor the G141 bump between wavelengths of . We discuss these latter two components of the transmission spectrum further in Sections 5.3 and 5.4, respectively. For the remaining data – namely, the STIS data spanning the wavelength range and the WFC3 data covering the HO band centered at – we find a good match is obtained for the model with a temperature of  K and metallicity of solar (Figure 5), which has a reduced of 1.0 for 69 degrees of freedom. Similarly good matches to the data are obtained for the  K models with metallicities of and solar. These metallicities are broadly consistent with predictions for a planet such as WASP-121b (2016ApJ...831...64T).

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHref(Left panel) Individual contributions to the transmission spectrum due to the major radiatively-active species in the best-match forward model shown in Figure 5, i.e. chemical equilibrium for  K and solar metallicity. Note that continuum opacity due to gas-phase species such as H and He is not shown. (Right panel) Corresponding pressure-dependent abundances.

Aside from collision-induced asorption and gas-phase Rayleigh scattering, the primary opacity sources of these models are Na and VO at optical wavelengths and HO at NIR wavelengths. This is illustrated in Figure 5.2, which shows a break-down of the opacity sources in the best-matching chemical equilibrium model. Interestingly, opacity due to TiO is not as significant as VO in the optical, despite Ti being approximately an order of magnitude more abundant than V for solar elemental composition (2009ARA&A..47..481A). This occurs because for a given pressure, the condensation of Ti species commences at higher temperatures than for V species (e.g. 1999ApJ...512..843B; 2018A&A...614A...1W). The isothermal temperature of the best-match model (i.e. 1500 K) is less than the condensation temperature of both TiO(s) and VO(s), meaning that these are the dominant forms of Ti and V in the model, respectively. However, since the isothermal temperature is closer to the VO(g)/VO(s) condensation temperature than the TiO(g)/TiO(s) condensation temperature, the abundance of VO(g) is larger than for TiO(g).

In contrast, 2002ApJ...577..974L found that calcium titanates (e.g. CaTiO) – which are not currently included in ATMO – are likely to be the first Ti-bearing condensates to form. Furthermore, arguing from trends in solar system meteorite data and M/L dwarf spectra, Lodders notes that V will likely condense in solid solution with the calcium titanates, resulting in VO gas-phase depletion commencing at the same temperature as TiO gas-phase depletion. However, for hot Jupiters, Ti and V condensation may depend on condensation and mixing timescales, both vertical and horizontal, that are very different to the protostellar nebula and M/L dwarfs. Such details are complex and beyond the scope of the present study. At this stage, we simply note that VO absorption is favored by these HST data for WASP-121b, with no evidence for significant TiO absorption, an interpretation that is corroborated by the free-chemistry retrieval presented in Section 5.4 below.

We also note that the best-matching forward model temperature of  K is substantially cooler than that of the dayside photosphere, which is inferred to be  K from secondary eclipse measurements (2017Natur.548...58E). Such a large temperature difference between the dayside photosphere probed during secondary eclipse and the upper atmosphere of the day-night limb probed during primary transit is in fact broadly in line with predictions of 3D general circulation models of ultra-hot Jupiters (e.g. 2016ApJ...821....9K). Furthermore, the best-match temperature of  K is likely to be at the lower end of the plausible range, because, as noted above, the forward models we consider here do not include rainout chemistry. Rainout chemistry will likely result in VO condensing at higher temperatures, as the abundance of VO in the upper atmosphere would be determined by the atmospheric temperature profile at higher pressures where the condensation temperature is also higher. Since the appearance or disappearance of VO spectral bands is primarily what determines the ability of our forward models to match the data (Figure 5), forward models with rainout chemistry would consequently tend to favor higher temperatures. As noted above, we do not consider models with rainout here, as the details will be highly sensitive to the atmospheric PT profile at pressures  bar, which is unconstrained by the current data.

5.3 Absorption at NUV wavelengths

We now consider the steep rise in the transmission spectrum at wavelengths shortward of . As explained in Section 5.1, we consider it unlikely that this feature can be explained by Rayleigh scattering due to gas-phase species such as H or high-altitude aerosols. In addition, our chemical equilibrium models presented in Section 5.2 do not predict significant absorption above the H continuum at these wavelengths. Nonetheless, we find the rise of the transmission spectrum at NUV wavelengths is empirically repeatable. It is recovered by our analysis when the spectroscopic lightcurves for the two G430L visits are fit jointly and also when they are each fit individually (see Section B.4).

One candidate absorber is the mercapto radical, SH, comprised of a sulfur atom and a hydrogen atom. Indeed, SH was predicted by 2009ApJ...701L..20Z (Z09) to be a strong NUV absorber in hot Jupiter atmospheres. Using a 1D photochemical kinetics code, Z09 found the abundance of SH may peak at pressures around mbar in typical hot Jupiter atmospheres, with a mixing ratio of  ppm (see Figure 2 of Z09). At these pressures, HS is the most abundant sulfur-bearing phase under chemical equilibrium (2006ApJ...648.1181V), while atomic H and S are also available due to photodissociation of molecules such as H and HO. The production of SH then proceeds through numerous chemical pathways involving HS, H, and S (Z09; see also 2016ApJ...824..137Z).

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefAbsorption cross-sections for SH. Electronic transitions are from 2009ApJ...701L..20Z and rotational-vibrational transitions are from ExoMol (2018MNRAS.tmp..913Y). Cross-sections for HO, VO, and Na are also shown, weighted by the relative abundances implied by the model shown in Figure 5.2.

To explore whether or not SH can explain the observed NUV absorption, we performed a simple fit to the 13 shortest-wavelength data points of the transmission spectrum, spanning the wavelength range (i.e. the blue G430L data subset indicated by light blue halos in Figures 5 and 5.1). As in Section 5.1, we followed the methodology outlined in L08. We computed the change in relative planetary radius due to SH absorption, adopting a planetary surface gravity  cm s and stellar radius (2016MNRAS.tmp..312D). We also assumed atomic mass units (see Section 4) and set as the altitude where H becomes optically thick at grazing geometry for a wavelength of  nm (Figure 5), corresponding to a planetary radius of . This in turn translates to an atmospheric pressure of  mbar, assuming a temperature of K and an H scattering cross-section of  cm molecule for  nm (see Section 4.1 of L08; also, 2016Natur.529...59S). Having thus established the pressure scale, we took the temperature-dependent absorption cross-sections for SH and varied the mixing ratio to optimize the match to the NUV transmission spectrum using Equation 1 of L08. For the SH cross-sections, we combined those derived by Z09 with those recently published by the ExoMol project (2018MNRAS.tmp..913Y). Specifically, the Z09 cross-sections were generated from transitions of the lowest five vibrational levels of the ground electronic state to the lowest three vibrational levels of the upper electronic state (without predissociation), and exhibit a strong NUV signature. These transitions are not considered in the ExoMol cross-sections, which only account for rotational-vibrational transitions. Both the Z09 and ExoMol cross-sections are shown in Figure 5.3.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefAbundance predictions for important sulfur species assuming solar metallicity and  cm s. Dashed green line indicates the adopted PT profile, based on the limb average of a 3D GCM for WASP-121b (Kataria et al., in prep). Calculations were performed using the photochemical kinetics code of 2009ApJ...701L..20Z, assuming a planet with a hydrogen-dominated atmosphere orbiting an F6V host star at the same distance as WASP-121b.

The results of this process are shown in Figure 5.1. We obtain respectable matches to the data with mixing ratios of  ppm and  ppm, respectively, for the  K and  K absorption cross-sections of Z09. For comparison, Figure 5.3 shows predicted abundances from the photochemical kinetics code of Z09 for a planet similar to WASP-121b with solar metallicity and vertical mixing parameter  cm s. We find abundances of  ppm are plausible for SH across the bar to mbar pressure range probed by the transmission spectrum, lending some credibility to the hypothesis that it could be the mystery NUV absorber. However, as stressed by Z09, the SH cross-sections remain subject to considerable uncertainty, due to the paucity of available experimental data. This, combined with the low spectral resolution of the G430L data, prevent us from conclusively confirming or ruling out SH at the present time. Other sulfur-bearing compounds that are likely to be abundant, such as SiS, have strong features at NUV wavelengths but remain poorly modeled. 2018ApJ...866...27L have also flagged gas-phase Fe as an important NUV absorber in ultra-hot Jupiter atmospheres, although we find it is unable to account for the measured signal in the present dataset – at least under assumptions of equilibrium chemistry for pressures  bar – as it was included in the ATMO forward models described in Section 5.2 (see Figure 5.2).

Regardless of the identity of the putative NUV absorber, it likely provides significant heating of the upper atmosphere. For instance, across the wavelength range, the mean SH absorption cross-section varies from to  cm molecule (Figure 5.3). Assuming a mixing ratio of  ppm, in line with our above estimates, this implies a mean atmospheric opacity (i.e. absorption cross-section mixing ratio) of  cm molecule at the pressures probed in transmission. We find incorporating such an absorber into the 1D radiative-convective atmosphere model of Marley and collaborators (e.g. 1999Icar..138..268M; 2002ApJ...568..335M; 2008ApJ...678.1419F; 2008ApJ...689.1327S; 2012ApJ...754..135M) would likely heat the atmosphere of WASP-121b by  K at mbar pressures. Such heating could, for example, help maintain optical absorbers such as VO and TiO in the gas phase, which in turn would provide further heating of the upper atmosphere. Properly accounting for effects such as these may be important for accurate modeling of planetary circulation and energy budgets.

Absorption of incident UV flux exceeding that predicted by models has also been observed in solar system atmospheres. Two well-known examples are Venus and Jupiter. On Jupiter, a broad reflectivity dip near has been attributed to a high-altitude dust or haze (1972Icar...16..557O; 1972ApJ...173..451A). The composition of this chromophore is still not known and is generally attributed to some disequilibrium combination of S, N, C, and P species (for a fuller discussion see 2004jpsm.book...79W). Likewise on Venus, dark markings in the atmosphere at UV wavelengths remain poorly understood, well over four decades after their discovery (e.g., 1997veii.conf..415E). These features have also been attributed to some disequilibrium – perhaps S-bearing – absorber (but see 1980JGR....85.8141P).

5.4 Retrieval analysis of optical-NIR data

In addition to comparing the data with predictions of forward models that assume chemical equilibrium (Section 5.2), we performed a free-chemistry retrieval analysis. For these calculations, we treat the abundances of the radiatively-active chemical species as free parameters in the model, rather than solving for the chemical equilibrium abundances at a given temperature. As for the forward models, this was done using ATMO, which can compute transmission spectra for any given atmospheric composition and PT profile. ATMO has previously been used for retrieval analyses of transmission spectra (2017Sci...356..628W; 2018AJ....155...29W) and thermal emission spectra (2017Natur.548...58E).

Since ATMO does not currently include any opacity sources that can explain the steep rise observed at NUV wavelengths (Figure 5), we restricted the retrieval to optical-NIR wavelengths longward of . We assumed an isothermal PT profile and allowed the limb-averaged temperature () to vary as a free parameter, as well as the reference planet radius corresponding to the 1 mbar pressure level () effectively providing a floating offset between the model and data. The abundances of HO, TiO, VO, Na, and FeH were allowed to vary relative to a background atmosphere composition dominated by H and He, assuming uniform mixing ratios with pressure. Other gas-phase absorbers such as K and CO were fixed to equilibrium abundances for the final analysis, as these were found to be unconstrained by the current data. Opacity due to aerosol Rayleigh scattering and optically-thick gray cloud was treated using the approach of 2016Natur.529...59S. Fitting was performed using differential evolution MCMC (2013PASP..125...83E), as described in our previous work (2017Natur.548...58E; 2017Sci...356..628W; 2018AJ....155...29W).

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefResults of the free-chemistry retrieval analysis. (Off-diagonal panels) Heat maps showing density of samples drawn from the MCMC analysis for different pairs of parameters. (Diagonal panels) Marginalized density distributions for individual parameters. Solid orange lines indicate parameter median values and dashed orange lines indicate ranges spanning 68% of samples.

The inferred distributions for each model parameter are summarized in Table 6 and shown in Figure 5.4. We find no evidence for opacity contribution due to aerosols, irrespective of whether they are treated as an enhanced Rayleigh scattering component or an optically-thick gray cloud. For this reason, we only present the results for the case including gray cloud, as the specific treatment of aerosols has negligible impact on the values inferred for the other model parameters.

We obtain a limb-averaged temperature of  K, in agreement with the best-matching chemical equilibrium model presented in Section 5.2. The inferred abundances for HO ( dex), VO ( dex), and TiO ( dex) are also in good agreement with those predicted by the best-matching equilibrium model (Figure 5.2). The inferred abundance for Na ( dex; lower limit of  dex) is somewhat higher than the solar value of  dex (Figure 5.2). One possibility is that the core of the Na line is probing the planetary thermosphere, where temperatures are higher and the pressure scale height is larger. This would produce a strong Na feature that the retrieval may misinterpret as indicating a high abundance. For instance, 2012MNRAS.422.2477H detected a strong Na line in the STIS transmission spectrum for HD 189733b, which high-resolution spectroscopy showed is caused by a thermosphere (2015A&A...577A..62W).

In addition, the inferred abundance for FeH ( dex) is orders of magnitude greater than expected for solar metallicity (Figure 5.2). Such a high FeH abundance – which we consider implausible – is driven by the bump in the measured transmission spectrum across the wavelength range, where FeH has a significant absorption signature (e.g. see Figure 7 of 2007ApJS..168..140S). This can be seen clearly in Figure 5.4, which shows the distribution of spectra implied by the retrieval analysis, compared with the best-matching chemical equilibrium model. The inability of our model to simultaneously explain the bump and the rest of the data results in a moderately-poor overall fit, with a reduced of 1.5 for 67 degrees of freedom.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSimilar to Figure 5, but showing the distribution of model spectra inferred by the retrieval analysis as well as a hypothetical signal due to SH. The dark green line shows the sample mean at each wavelength and the shaded green areas progressively encompass 68%, 95.5%, and 99.7% of samples about the mean. A significant departure from the chemical equilibrium model occurs between wavelengths of . This is due to the retrieval inferring a high FeH abundance to explain the bump in the transmission spectrum measured over the short-wavelength half of the G141 bandpass.

The bump in the transmission spectrum remains puzzling. It has been recovered by multiple independent analyses of the data performed within our own group, as well as those published by others (e.g. 2018AJ....155..156T). We note that it coincides with a possible spectral feature identified in the dayside thermal spectrum, which was tentatively attributed to VO emission (2017Natur.548...58E). However, it is difficult to reconcile VO with the feature seen in the transmission spectrum, as it would require increasing the abundance to a level that would be incompatible with the data at optical wavelengths, where VO has a higher opacity. On the other hand, although the host star is photometrically quiet and care has been taken to precisely measure the absolute transit depths for each bandpass (G430L, G750L, G141), it is conceivable that small offsets in remain, which, if accounted for, could allow VO to simultaneously explain the transmission spectrum at optical wavelengths along with the feature. For multi-epoch observations that do not overlap in wavelength such as those considered here, it is impossible to rule out such a scenario with absolute confidence. Upcoming G141 observations should allow a determinination of whether or not the bump is repeatable. It is also worth noting that a strong thermal gradient over the pressures probed in transmission – which has not been considered in the present study – could potentially affect the shape of the transmission spectrum by altering the pressure-dependent scale height and chemical abundances.

In summary, the retrieval analysis reveals no evidence for aerosols in the optical-NIR transmission spectrum of WASP-121b. The inferred limb-averaged temperature and gas-phase abundances are overall in good agreement with the best-matching forward models of Section 5.2, which assume chemical equilibrium and solar metallicity. The primary exception is the inferred FeH abundance, which as described above is far higher than expected for chemical equilibrium and solar metallicity. We thus conclude that it is unlikely FeH opacity is the true cause of the spectral bump at wavelengths , the provenance of which remains uncertain.

6 Conclusion

We have presented a STIS transmission spectrum for WASP-121b, spanning the wavelength range, adding to the wavelength coverage of published WFC3 data. The new optical data show an increase in atmospheric opacity for wavelengths shortward of , with a slope that is too steep to be explained by Rayleigh scattering. Instead, assuming the NUV rise is a bona fide signature of the planetary atmosphere, it must be caused by one or more absorbers. We propose SH as a possible candidate, with a mixing ratio of approximately  ppm. Although the identity of the NUV absorber remains uncertain, it should cause substantial heating of the upper atmosphere and therefore could be an important component missing from existing models of highly-irradiated atmospheres. At longer wavelengths between , we measure significant opacity variations that can be well-explained by VO absorption. Analyzing the STIS and WFC3 data with both free-chemistry retrievals and comparisons to chemical equilibrium forward models, we estimate abundances of HO and VO approximately solar. We find no significant evidence for TiO, suggesting it may have condensed from the gas-phase. Our chemical equilibrium forward models are unable to simultaneously reproduce the optical data and the WFC3 bump spanning the wavelength range. Free-chemistry retrievals are able to do so, but only by invoking an unrealistically high FeH abundance.

Overall, the evidence uncovered here for significant NUV and optical absorption implies a substantial fraction of incident stellar radiation is likely deposited at low pressures in the atmosphere of WASP-121b. Heating via this mechanism could be responsible for the thermal inversion detected on the dayside hemisphere. The broad coherence of this picture is tantalizing, but many unknowns remain. Although we consider the evidence for VO in the existing transmission spectrum to be reasonably strong, further observations are required to confirm or rule it out at high confidence. Similarly, additional observations, along with a more extensive exploration of candidates other than SH, are required to identify the NUV absorber. The possible explanation provided by SH, however, flags the potential importance of non-equilibrium sulfur chemistry in highly-irradiated atmospheres, which until now has received relatively little attention.

\bibliography@latex

wasp121

\H@refstepcounter table \hyper@makecurrenttable \@setminipage Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefResults of white lightcurve MCMC analyses. Quoted values give sample medians with uncertainties corresponding to ranges encompassing of samples about the median. Note that values were not fit directly as part of the MCMC analysis but obtained by multiplying the values by the formal photon noise values for each lightcurve. and allowed to vary G430Lv1 G430Lv2 G750L G141 () (MJD) (ppm) (ppm) and held fixed G430Lv1 G430Lv2 G750L G141 (fixed) (fixed) () (fixed) (MJD) (ppm) (ppm)
\H@refstepcounter table \hyper@makecurrenttable \@setminipage Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefResults of G430L spectroscopic lightcurve fits for selected parameters. (Å) (ppm) (ppm) 2898-3499 3499-3700 3700-3868 3868-4041 4041-4151 4151-4261 4261-4371 4371-4426 4426-4481 4481-4536 4536-4591 4591-4646 4646-4701 4701-4756 4756-4811 4811-4921 4921-4976 4976-5030 5030-5085 5085-5140 5140-5195 5195-5250 5250-5305 5305-5360 5360-5415 5415-5469 5469-5524 5524-5579 5579-5634 5634-5688
\H@refstepcounter table \hyper@makecurrenttable \@setminipage Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSimilar to Table 6, but for the G750L spectroscopic lightcurve fits. (Å) (ppm) 5263-5550 5550-5648 5648-5745 5745-5843 5843-5940 5940-6038 6038-6135 6135-6233