The Hubble Space Telescope PanCET Program: Exospheric Mg ii and Fe ii in the Near-ultraviolet Transmission Spectrum of WASP-121b Using Jitter Decorrelation

The jitter files include 28 different engineering measurements (see Table 1). To determine the suitability of each of the jitter optical state vectors for decorrelating transit or eclipse light curves, we measured the linear Pearson correlation coefficient, R, between each jitter vector and the out-of-transit white-light-curve photometric data for 61 separate HST visits covering close to 300 spacecraft orbits. Specifically, we used all of the STIS G430L and G750L data from 23 visits of program GO-12473 (P.I. Sing), G430L, and G750L data from 34 visits of GO-14767 (P.I.s Sing & López-Morales), and four G430L visits from GO-14797 (P.I. Crossfield). The data from each HST STIS visit were reduced in the same manner as detailed in Sing et al. (2016), and we summed the counts of each extracted spectra to produce a white-light curve for each visit. Using the out-of-transit data, the percent of variance in common between jitter decorrelation parameters and the STIS white-light-curve photometry, R², was determined for each of the 61 visits. We also measured the quantity for the traditional detrending variables (ϕ_HST, X_psf, Y_psf, and S_λ) with the results given in Table 1. As a point of comparison, the individual R² values for a G430L eclipse of WASP-12b are reported in addition to the two STIS E230M transits of WASP-121b which are visits 97 and 98 of the PanCET program. The raw white-light curve of the STIS G430L eclipse of WASP-12b reaches a photometric precision of 979 ppm, as measured by the standard deviation of the light curve about the mean, which is a factor of 6.4× higher than expected from the photon noise levels.

Table 1.  Jitter Engineering Data Optical State Vectors (Column 1) and the Correlation R² Values (in %) with the White Light Curve Photometry for the Visits Covering a STIS/G430L Eclipse of WASP-12 (Column 1) and the STIS/E230M Data of WASP-121 for Visits 97 and Visit 98 (Columns 3 and 4, Respectively)

Vector	W-12	W-121	W-121	$\overline{{R}^{2}}$	${\sigma }_{\overline{{R}^{2}}}$
(1)	(2)	(3)	(4)	(5)	(6)
ϕHST	36	54	17	34	28
Sλ	86	0	2	30	27
Xpsf	26	30	1	19	22
Ypsf	31	30	1	16	20
V2_dom	16	13	6	7	9
V3_dom	1	1	0	6	8
V2_roll	57	36	7	36	29
V3_roll	54	34	9	35	29
SI_V2_AVE	0	11	19	4	5
SI_V2_RMS	56	8	0	19	20
SI_V2_P2P	58	15	0	19	20
SI_V3_AVE	4	14	15	5	6
SI_V3_RMS	55	29	6	12	15
SI_V3_P2P	55	5	2	13	15
R.A.	55	37	8	31	28
Decl.	57	36	7	29	27
Roll	56	50	13	14	16
LimbAng	8	0	3	10	14
TermAng	0	7	3	8	9
LOS_Zenith	8	0	3	10	14
Lat	13	41	7	20	21
Long	76	6	0	38	29
Mag_V1	1	22	8	13	15
Mag_V2	12	4	0	16	19
Mag_V3	12	29	4	15	16

Note. The average R² value from 61 HST STIS G430L & G750L CCD visits, $\overline{{R}^{2}}$, is also reported in column 5, and the standard deviation of R² in column 6.

Download table as:  ASCIITypeset image

An additional systematic trend is seen in the STIS CCD time-series data, where the first exposure of each spacecraft orbit is consistently found to exhibit significantly lower fluxes than the remaining exposures (Brown et al. 2001; Sing et al. 2011). Attempts were made to mitigate this effect in the observational setup by employing a quick 1 s integration (Sing et al. 2015), which was discarded. However, higher overall correlation values were found between the STIS CCD white-light curves and the optical state vectors if, in addition, the first long science exposure of every orbit was also systematically discarded from the analysis. We suspect the 1 s integration technique is not an overall effective method to mitigate the first-orbit-exposure systematic, and the first long science exposure of any STIS CCD orbit in a time-series analysis should be treated with caution. For the correlation analysis and for all subsequent CCD light curve fitting throughout this paper, we chose to discard the first long science exposure of each HST orbit.

From the correlation analysis, we identify several pairs of jitter vectors that are often found to correlate well with the STIS data (R² ≳ 30%, also see the Appendix). The roll of the telescope along the V2 and V3 axis (V2_roll, V3_roll) and the R.A. and decl. of the aperture reference R.A. and decl. are two pairs of vectors that are both related to the pointing of the telescope and in some data sets are found to highly correlate with the data with R² values as high as 90%. In addition, the HST sub-point latitude and longitude (Lat and Long respectively) are seen to have large correlations. The latitude and longitude are similar in nature (but not identical) to ϕ_HST, which is already included in the traditional model. From HST visit to visit, the correlation values themselves for any particular vector are found to change by large amounts, which is reflected in Table 1 by the large standard deviation values, ${\sigma }_{{R}^{2}}$. However, in retrospect a large spread in correlation values related to the telescope's position is to be expected. Different targets across the sky observed at different telescope orientation values will not necessarily favor any particular telescope axis or R.A. and decl. direction, but rather the actual vector direction(s) the target PSF takes with respect to the instrument during a transit observation. We note a general trend where the R² values in visits with relatively good pointing show lower levels of correlation (e.g., visit 98), while visits with poorer quality pointing show higher levels of correlation (e.g., visit 97).

Trends with telescope position are consistent with photometric systematics caused by slit light losses. If correct, the trends appear even when using slit sizes that are much greater than the size of the PSF. Our hypothesis is that the known telescope breathing leads to small changes in the position of the PSF on the detector (usually at the sub-pixel level), which largely manifest themselves as slit light losses. Even if the central PSF is well centered on a much wider slit, diffraction spikes and the wide wings of the PSF could contribute to light losses. If the dominant sources of systematics with STIS are position-related slit light losses, assuming similar guiding performance, then presumably there could be larger systematics when using small slits such as the 02 × 02 used with the STIS E230M compared to the 52'' × 2'' slit used with the CCD. In addition, there would be larger systematics in visits with poorer guiding performance, which appears to be the case. Several STIS visits in the PanCET program had guide star acquisition problems, which resulted in reduced pointing accuracy and light curves with comparably larger systematic trends. We also observe that the jitter vectors (V2_roll, V3_roll) and (R.A., decl.) often have the same general trends as the traditional vectors (S_λ, X_psf, and Y_psf), which are measured directly from the spectra. Similar trends are generally expected, as changes in the PSF position (measured through the telescope position in the jitter files) will also be recorded by changes in the placement of the target stellar spectra on the detector. These vectors are not identical, however, as the detector-measured vectors will also be sensitive to further detector effects such as the pixel-to-pixel flat-fielding, while the jitter vectors are able to record the telescope position regardless of whether the target PSF is contained within the slit.

As a test, we implemented jitter detrending on the STIS G430L eclipse data of WASP-12b, finding good agreement in the eclipse depth with Bell et al. (2017) who used a Gaussian process method to decorrelate the time series data (see the Appendix). However, several improvements can be seen utilizing jitter detrending. The best-fit measured eclipse depth was found to be positive while Bell et al. (2017) found a negative value. In addition, with jitter detrending the first orbit was successfully recovered for use in the analysis, while Bell et al. (2017) had to discard it.

We observed two transits of WASP-121b with the HST STIS E230M echelle grating during 2017 February 23 (visit 97) and 2017 April 10 (visit 98). The observations were conducted with the NUV-MAMA detector using the Echelle E230M STIS grating and a square 02 × 02 entrance aperture. The E230M spectra had resolving power of λ/(2Δλ) = 30,000 and was configured to the 2707 Å setting to cover the wavelength ranges between 2280 and 3070 Å in 23 orders (see Figure 1). The resulting spectra had a dispersion on the detector of approximately 0.049 Å/pixel. Both transits covered five HST spacecraft orbits, with the transit event occurring in the third and fourth orbits. During each HST orbit, the spectra were obtained in TIME-TAG mode, where the position and detection time of every photon was recorded in an event list, which had 125 μs precision.

Zoom In Zoom Out Reset image size

We sub-divided the spectrum of each HST spacecraft orbit into sub-exposures with the Pyraf task inttag, each with a duration of 274.01996 s. This allowed orbits 2 through 5 to be divided into 10 sub-exposures each, while the first orbit was divided into seven sub-exposures. These sub-exposures were then each reduced with CALSTIS version 3.17, which includes the calibration steps of localization of the orders, optimal order extraction, wavelength calibration, and flat-field corrections. The mid-time of each exposure was converted into barycentric Julian dates in barycentric dynamical time (BJD_TBD) for use in the transit light curves (Eastman et al. 2010).

The light curves were modeled using the analytical transit models of Mandel & Agol (2002). For the white-light curves, the central transit time, planet-to-star radius ratio, stellar baseline flux, and instrument systematic trends were fit simultaneously, with flat priors assumed. The inclination and stellar density (or equivalently the semimajor axis to stellar radius ratio, a/R_star) were held fixed to the values found in Evans et al. (2018), as fitting for these parameters found values consistent with the literature, though with an order of magnitude lower precision. For example, we find a value of a/R_star = 3.63 ± 0.16 from the NUV white-light curve of visit 98, while Evans et al. (2018) reports a value of 3.86 ± 0.02. Compared to the NUV data here, the study of Evans et al. (2018) benefits from much higher photometric precision in a larger data set with a wider array of wavelengths, including STIS data with complete phase coverage of the transit, which constrains the planet's orbital system parameters to a much greater degree than the NUV data alone.

The total parameterized model of the flux measurements over time, f(t), was modeled as a combination of the theoretical transit model, $T(t,{\boldsymbol{\theta }})$ (which depends upon the transit parameters ${\boldsymbol{\theta }}$), the total baseline flux detected from the star, F₀, and the instrument systematics model $S({\boldsymbol{x}})$, giving

$$\begin{eqnarray}&&f(t)=T(t,\theta )\times {F}_{0}\times S({\boldsymbol{x}}).\end{eqnarray} \tag{ 1 }$$

Based on the results of Section 2 and the Appendix, for our most complex systematics error model tested, we included a linear baseline time trend, ϕ_t, as well as the traditional optical state vectors (ϕ_t, ϕ_HST, ${\phi }_{{HST}}^{2}$, ${\phi }_{{HST}}^{3}$, ${\phi }_{{HST}}^{4}$, X_psf, Y_psf and S_λ) and jitter vectors (V2_roll, V3_roll, R.A., decl., Lat, and Long) resulting in up to 14 terms used to describe $S({\boldsymbol{x}})$ depending on the data set in question as described below.

The errors on each data point were initially set to the pipeline values, which are dominated by photon noise. The best-fitting parameters were determined simultaneously with a Levenberg–Marquardt (L-M) least-squares algorithm (Markwardt 2009) using the unbinned data. After the initial fits, the uncertainties for each data point were rescaled to match the standard deviation of the residuals. A further scaling was also applied to account for any measured systematic errors correlated in time ("red noise"). After rescaling the error bars, the light curves were then refit, thus taking into account any underestimated errors in the data points.

The red noise was measured by checking whether the binned residuals followed an N^−1/2 relation, when binning in time by N points. In the presence of red noise, the variance can be modeled to follow a ${\sigma }^{2}={\sigma }_{{\rm{w}}}^{2}/N+{\sigma }_{{\rm{r}}}^{2}$ relation, where σ_w is the uncorrelated white noise component, and σ_r characterizes the red noise (Pont et al. 2006). For our best-fitting models for $S({\boldsymbol{x}})$, we did not find evidence for substantial red noise.

The uncertainties on the fitted parameters were calculated using the covariance matrix from the Levenberg–Marquardt algorithm, which assumes that the probability space around the best-fit solution is well-described by a multivariate Gaussian distribution. Previous analyses of other HST transit observations (Berta et al. 2012; Line et al. 2013a; Sing et al. 2013; Nikolov et al. 2014) have found this to be a good approximation when fitting HST STIS data. We also computed uncertainties with a Markov chain Monte Carlo (MCMC) analysis (Eastman et al. 2013), which does not assume any functional form for this probability distribution. In each case, we found equivalent results between the MCMC and the Levenberg–Marquardt algorithm for both the fitted parameters and their uncertainties, as the posterior distributions were found to be Gaussian. Inspection of the 2D probability distributions from both methods indicates that there are no significant correlations between the systematic trend parameters and the planet-to-star radius contrast.

The effects of stellar limb darkening are strong at NUV wavelengths. To account for the effects of limb darkening on the NUV transit light curve, we adopted the four-parameter, nonlinear limb-darkening law, calculating the coefficients as described in Sing (2010). For the model, we used the 3D stellar model from the Stagger-grid (Magic et al. 2015) with the model (T_eff = 6500, log g = 4, [Fe/H] = 0.0), which was closest to the measured values of WASP-121A in effective temperature, gravity, and metallicity (T_eff = 6460 ± 140 K, log10 g = 4.242 ± 0.2 cgs, [Fe/H] = +0.13 ± 0.09 dex; Delrez et al. 2016). The Stagger-grid 3D models are calculated at a resolution of R = 20,000 which is close to the native resolution of the E230M data (R = 30,000), and contains a fully line-blanketed NUV region that matches the flux-calibrated data of WASP-121A well (see Figure 2). We included the individual responses from each order in the echelle spectra and converted the stellar model spectral wavelengths from air to vacuum. We also applied a wavelength shift to the data to take into account the systemic radial velocity of the system, 38.36 km s⁻¹ (Gaia Collaboration et al. 2018), shifting the star to the rest frame.

Zoom In Zoom Out Reset image size

In a test to see how well the 3D models performed, we fit visit 98 with a three-parameter limb-darkening law (see Sing 2010) and let the linear coefficient, c₂, fit freely while the other two nonlinear parameters were fixed to the model values. We found good agreement at the 1σ level between the fit coefficient c₂ = 0.881 ± 0.215 and the theoretical value of 1.077. In addition, the fit R_pl/R_star = 0.1415 ± 0.0050 was also consistent (within 1σ) when fitting for c₂ versus fixing the limb darkening (see Section 4.4). For the remainder of the study, we fixed the limb-darkening coefficients to the model values. We note that using the second-order Akaike information criterion (AICc) for model selection (see the Appendix), that the AICc favored fixing the limb-darkening parameters to the model values.

The Transiting Exoplanet Survey Satellite (TESS) observed WASP-121b in camera 3 from 2019 January 7 to 2019 February 2, and we obtained the time-series photometry though the Mikulski Archive for Space Telescopes' exo.MAST²⁰ web service. We fit the TESS WASP-121b transits using the Presearch Data Conditioning light curve, which has been corrected for effects such as non-astrophysical variability and crowding (Stumpe et al. 2012; Jenkins et al. 2016). From the time series, we removed all of the points that were flagged with anomalies. The time-series barycentric TESS Julian dates were converted to BJD_TDB by adding 2,457,000 days. The TESS light curve contains data for 16 complete transits of WASP-121b and one partial transit. For each transit in the light curve, we extracted a 0.5 day window centered around the transits and fit each transit event individually. We fit the data using the model as described in Sections 4.1 and 4.2, though only included a linear baseline time trend, ϕ_t, for $S({\boldsymbol{x}})$. We found a weighted-average value of R_pl(TESS)/R_star = 0.12342 ± 0.00015, which is in good agreement with the HST transmission spectrum of Evans et al. (2018). We converted the transit times of Evans et al. (2018) and (Delrez et al. 2016) to BJD_TDB using the tools from Eastman et al. (2010) and fit these along with the TESS transit times and visit 98 from Section 4.4 with a linear function of the period P and transit epoch E:

$$\begin{eqnarray}&&T(E)={T}_{0}+{EP}.\end{eqnarray} \tag{ 2 }$$

We find a period of P = 1.2749247646 ± 0.0000000714 (days) and central transit time of T₀ = 2457599.551478 ± 0.000049 (BJD_TDB). We find a very good fit with a linear ephemeris (see Figure 3), having a χ² value of 16.8 for 21 degrees of freedom (dof). We find no obvious signatures of quasi-periodic photometric variabiltiy due to stellar activity, with the TESS data binned by 62 minutes (excluding the transit and eclipse times) giving a standard deviation about the mean of 0.07%, which is only modestly higher than the expected photometric precision of 0.02%.

Zoom In Zoom Out Reset image size

Visits 97 and 98 show dramatically different levels of instrumental systematic trends (see Figure 4). In particular, visit 97 shows a ∼4% change in flux in the fifth out-of-transit orbit, while for visit 98 there is no such large trend, and the light curve points are within ∼0.4% of each other. Upon inspection of the telescope R.A. and decl. from the jitter files, it is clear that visit 97 suffers from a very large drift during the five orbits covering the transit event, while visit 98 has much more stable pointing (see Figure 5). The large drift likely leads to large slit light losses in visit 97. Most importantly, the R.A. and decl. positions of the telescope during the transit in visit 97 are significantly different from the out-of-transit positions, especially the fifth orbit, which was shifted by about a quarter of a pixel. As discussed in Gibson et al. (2011), detrending with optical state vectors can only be expected to work if we interpolate the vectors for the in-transit orbit(s) with the out-of-transit orbits, which requires the baseline function to be well represented by the linear model over the range of the decorrelation parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each visit, we used the AICc and measured σ_r to determine the optimal optical state vectors to include from the full set without overfitting the data while minimizing the red noise. As found in Section 2, the different position-related vectors from the jitter files typically contained similar trends. For visit 98, we found it was optimal to include ϕ_HST terms up to second order, as well as the jitter detrending vectors V2_roll, V3_roll, and R.A. We rotated the roll vector pair (V2_roll, V3_roll), which are contained in the engineering jitter files relative to the axis of the square slit (Vn_roll, Vt_roll), which is rotated by about 45° relative to the spacecraft orientation reference vector U3, such that $S({\boldsymbol{x}})$ could be written as

$$\begin{eqnarray}\begin{array}{rcl}S({\boldsymbol{x}}) & = & {p}_{1}{\phi }_{t}+{p}_{2}{\phi }_{{HST}}+{p}_{3}{\phi }_{{HST}}^{2}\\ & & +\,{p}_{4}{\rm{R.A.}}+{p}_{5}{{Vn}}_{\mathrm{roll}}+{p}_{6}{{Vt}}_{\mathrm{roll}}+1.\end{array}\end{eqnarray} \tag{ 3 }$$

The fit values of interest, namely R_pl/R_star, did not substantially change for any of the top fitting models (AICc–AICc_min ≲ 6). With nine free parameters (R_pl/R_star, F₀, T₀, p_1..6), the fit achieves a signal-to-noise which is 81% of the theoretical photon noise limit, with no detectable red noise and a ${\chi }_{\nu }^{2}$ of 1.21 for 37 dof. We note that including the jitter corrections and the second-order breathing polynomial not only provides a significantly better fit than the traditional systematics model, but we were also able to make use of the first orbit in the visit. The overall photometric performance is similar to that achieved with the STIS CCD using the G430L or G750L, despite the use of a much narrower slit. With visit 97, we measure the white-light curve R_pl/R_star = 0.1374 ± 0.0026 (see Table 2) and a T₀ = 2457854.536411 ± 0.001073 (BJD_TBD). The central transit time agrees very well with the expected ephemeris, occurring −0.2 ± 1.6 min relative to the expected central transit time using the updated ephemeris from Section 4.3.

Table 2.  WASP-121b Broadband Transmission Spectral Results and Nonlinear Limb Darkening Coefficients for the STIS E230M

λc	Δλ	RP/R*	${\sigma }_{{R}_{P}/{R}_{* }}$	c1	c2	c3	c4
(Å)	(Å)
2673	799	0.13735	0.00257	0.4011	−0.2625	1.2818	−0.4674
2387	236	0.15300	0.00845	0.4282	−0.8777	1.6772	−0.2513
2600	200	0.13962	0.00465	0.4826	−0.6976	1.9529	−0.7708
2800	200	0.13760	0.00366	0.4090	−0.2445	1.1797	−0.3940
2986	172	0.12734	0.00389	0.3374	0.2045	0.8312	−0.4348

Download table as:  ASCIITypeset image

For visit 97, even the most complex model did not achieve fit residuals as small as visit 98. As noted above, the fifth orbit of the visit displays dramatically increased systematic trends compared to the other orbits, with the position of the telescope during the orbit being relatively far from the position of the in-transit orbits (see Figure 5). As such, we find R_pl/R_star can change dramatically (by about 0.016 R_pl/R_star) between differing systematics models when including the data from the fifth orbit in the light-curve fits. Given these factors, we find that it is preferable to drop the last orbit from the analysis when measuring the absolute transit depths in visit 97. This is similar to many HST STIS and WFC3 transit studies that have dropped the first orbit as it usually displays larger non-repeatable trends due to changes in the thermal breathing or charge trapping. When excluding the fifth orbit, we found the fit R_pl/R_star values did not change substantially between differing systematics models (∼0.008 R_pl/R_star), and using the AICc, we found $S({\boldsymbol{x}})=S({RA}$, Vn_roll, Vt_roll, ϕ_t) to be optimal. The best-fit planet radius for visit 97 was found to be R_pl/R_star = 0.1364 ± 0.0110, which matches the value from visit 98 at well within 1σ, and achieves 43% of the theoretical photon noise limit.

When measuring the transmission spectrum, R_pl(λ)/R_star, we fixed the system parameters as they are not expected to have a wavelength dependence. In addition, we fixed the the limb-darkening coefficients to those determined from the stellar model for each spectroscopic passband, using the same method as described in Section 4.2. The light-curve-fitting methods were otherwise identical to those described in Section 4.1. We fit each spectroscopic light curve with the same functional systematics model as was found to optimally fit the white-light curve, fitting the various spectroscopic bins simultaneously for the wavelength-dependent R_pl(λ)/R_star and systematic parameters p_1..6.

Figure 6 shows our resulting broadband spectra from visit 98 in ∼200 Å bins, as well as the NUV white-light-curve value. R_pl(λ)/R_star is observed to sharply rise toward shorter wavelengths. The overall NUV transmission spectrum shows strong extinction and is substantially higher than the optical and near-infrared, reaching altitudes of ΔR_pl(λ)/R_star = 0.0157 ± 0.0026 higher in the atmosphere. This is more than 18× the size of the pressure scale height of the lower planetary atmosphere, H = kT/μg, that is H/R_star = 0.00084 at a temperature T of 2000 K.

Zoom In Zoom Out Reset image size

While the overall transit depths of visit 97 are dependent upon the choice of $S({\boldsymbol{x}})$, the steep rise in the NUV broadband spectrum can be seen regardless of including or excluding the fifth orbit from the analysis. When including the fifth orbit, a systematics model of $S({\boldsymbol{x}})$ = S(ϕ_HST, S_λ, R.A., decl., Vn_roll) produces a fit radius that is both consistent with visit 98 and independent of excluding the fifth orbit. The resulting broadband transmission spectrum is shown in Figure 6, which matches well with the results of visit 98.

Given the much lower precision and the sizeable systematic uncertainty of including or excluding the fifth orbit in visit 97, we elect to report the transmission spectra at higher resolutions using visit 98 only. We fit the data to various resolutions and report the transmission spectra at a resolution of R ∼ 650, corresponding to 4 Å bins with results from 5 Å bins also reported in Section 4.6. At these bin sizes, each light curve point has a few thousand photons and the ∼1.5% transit depth can typically be detected across the E230M wavelength range. However, in the spectral regions with strong stellar absorption lines, or at the edges of the spectral orders where the efficiency is low, the transit itself is not always detected due to the lack of flux. At this resolution, strong atomic transitions can be resolved and probed efficiently and resolution-linked bias is minimized (Deming & Sheppard 2017). Given the high resolution on the E230M, each 4 Å bin still contains ∼60 pixels. To help probe the central wavelength of the observed features more precisely, we also measured the transmission spectra in multiple 4 Å bin sets, each set shifted by 1 Å, corresponding to about 100 km s⁻¹ velocity shifts (see Figure 7). We removed points where the transit was not detected (R_pl(λ)/R_star < 0.05) or had large errors (${\sigma }_{{R}_{P}/{R}_{* }}\gt 0.09$), which predominantly occurred within strong stellar lines or at the order edges. The R ∼ 650 NUV transmission spectra can be seen in Figures 8–11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find no evidence for absorption by Mg i. In a 5 Å band centered on the ground-state Mg i line at 2852.965 Å, we measure a R_pl(Mg i)/R_star = 0.100 ± 0.056, which is consistent at the 1σ confidence level with the optical and near-infrared value.

We find strong evidence for absorption by Mg ii. We detect and resolve absorption by both the k and h features of the Mg ii ground-state doublet located at 2796.35 and 2803.53 Å respectively (see Figure 11). In 5 Å passbands, the transmission spectra in the Mg ii doublet are found to have radii of R_pl(Mg ii,k)/R_star = 0.284 ± 0.037 and R_pl(Mg ii,h)/R_star = 0.242 ± 0.0431, substantially larger than the values at optical and near-infrared (OIR) wavelengths (R_pl(OIR)/R_star = 0.1217) and larger than the average NUV value of R_pl(NUV)/R_star = 0.1374 ± 0.0026. A 10 Å passband split to cover the Mg ii doublet simultaneously is found to have a radius value of R_pl(Mg ii, h, k)/R_star = 0.271 ± 0.024 (see Figure 12), which is 5.4σ above the white-light-curve R_pl(NUV)/R_star value. In the continuum region surrounding the Mg ii doublet, no other substantial absorption features are observed and almost all of the spectral bins in a 100 Å region around the doublet are consistent with the R_pl(OIR)/R_star value (see Figure 11).

Zoom In Zoom Out Reset image size

The strong Mg ii absorption can also be seen in the transit light curves, with Figure 12 showing transit depths of about 8% within the Mg line, while the surrounding continuum is consistent with the OIR transit depth of 1.5%. The stellar Mg ii double line cores exhibit narrow emission lines, associated with the stellar corona (see Figure 2) with both peaks found just redward of the line centers at velocities near 20–40 km s⁻¹. We performed several checks to ensure that the presence of the emission line did not adversely affect the planetary Mg ii signal. Inspecting the photometric time series of just the stellar Mg ii emission lines, we do not observe any substantial variability differences from that of the surrounding continuum. In addition, when scanning the transmission spectrum over the Mg ii region, we find the peak of the transmission spectrum is found 50–100 km s⁻¹ blueward of the Mg ii line cores, rather than redward where the emission lines are located.

The transmission spectrum contains a very strong absorption feature at 2381.5 Å (R_pl(2381.5)/R_star = 0.332 ± 0.074), which can be identified as a ground-state resonant transition of Fe ii. Fe ii transitions occur in a set of multiplets, with notable transitions at 2600, 2382, and 2344 Å for the UV1, UV2, and UV3 multiplets respectively. For the aforementioned UV1 and UV3 ground-state transitions, the stellar line is strong enough such that almost no flux is measured in the line cores (see Figure 2), so the transmission spectrum of the planet cannot be measured with sufficient precision at those wavelengths. However, other candidate Fe ii features are also seen in the stellar spectrum (see Figures 8–10).

For Fe i, there are strong ground-state transitions in the NUV region at wavelengths of 2484, 2523, 2719, 2913, and 2937 Å. While a significant absorption feature does appear near 2484 Å (see Figure 8), the other transitions of Fe i do not obviously appear in the data.

To help further identify absorption features, we fit a 1D analytic transmission spectral model to the data (Lecavelier des Etangs et al. 2008), following the procedures as detailed in Sing et al. (2015) and Sing (2018). Formally, the analytic model is isothermal and assumes a hydrostatic atmosphere with a constant surface gravity with altitude. These assumptions are not expected to be valid over the large altitude ranges probed by the NUV data, especially at very high altitudes. However, the model is useful for identifying spectral features in the transmission spectra, ruling out absorption by different species, and getting a zeroth-order handle on the atmospheric properties, including the velocity. A more comprehensive and physically motivated modeling effort will be presented in a future work (P. Lavvas et al. 2019, in preparation).

To model the transmission spectra, a continuum slope was included, which was assumed to have cross-section opacity with a power law of index α, (σ(λ) = σ₀(λ/λ_o)^−α). We also included the spectral lines of Mg i, Mg ii, Fe i, and Fe ii with the wavelength-dependent cross-sections calculated using the NIST Atomic Spectra Database (Kramida et al. 2018). Voigt profiles were used to broaden the spectral lines. Given our transmission spectra are at very high altitudes, we did not include effects such as collisional broadening and chose instead to fit for the damping parameter governing the line broadening. We calculated the scale height assuming a mean molecular weight corresponding to atomic hydrogen and fit the model using an isothermal temperature. We also assumed that the number of atoms existing in the various atomic levels for a given species followed a Boltzmann distribution with the statistical weights given in the NIST database, though this assumption may not be accurate at very low pressures. To simplify the calculation, we only included transitions where the lower energy level was beneath a threshold. High energy levels will in practice not be sufficiently populated to produce significant absorption signatures, and for Fe ii we found transitions above 0.43 eV (corresponding to 5000 K) were not prevalent in the data. However, the excited states of Fe ii up to 0.43 eV were needed, as including only ground-state transitions resulted in a much worse fit to the data. The presence of these excited states indicates high temperatures, as a significant population of Fe atoms are found above the ground state.

We fit the model to the NUV data with an L-M least-squares algorithm and MCMC (Eastman et al. 2013). The fit contained five free parameters: the isothermal temperature, α, the Doppler velocity of the planetary atmosphere, the reference planet radius, and the Mg and Fe abundances. We fit for the data between 2300 and 2900 Å, finding a good fit with a χ² of 130.4 for 126 dof. We estimate the detection significance of Fe ii by excluding its opacity from the model and refitting, finding the χ² increases by a Δχ² = 61.8, which corresponds to a 7.9σ detection significance. The velocity of the planetary atmospheric Mg and Fe lines is measured to be $-{56}_{-63}^{+43}\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Given the velocity uncertainties are larger than the spectrograph resolution, the uncertainties can likely be improved using different data reduction methods (e.g., smaller bins or using cross-correlation techniques), which will be presented in a future work (P. Lavvas et al. 2019, in preparation).

A relatively cold planet interior can lead to the condensing of refractory species deep in the atmosphere, thereby removing the gaseous species from the upper layers of the atmosphere, trapping them within condensate clouds at depth (Hubeny et al. 2003; Fortney et al. 2008; Powell et al. 2018), a process that is dependent upon particle settling in the presence of turbulent and molecular diffusion (Spiegel et al. 2009; Koskinen et al. 2013). Al, Ti, VO, Fe, and Mg are among the first refractory species to condense out of a hot-Jupiter atmosphere (Visscher et al. 2010; Wakeford et al. 2017). As such, the presence or absence of these species can then provide constraints to the temperatures of the deeper layers for hot Jupiters, as the presence of these elements in the gas phase of the upper layers limits the global temperature profile to be hotter than the condensation of these species. For WASP-121b, the NUV transmission spectrum provides evidence for Mg and Fe in the exosphere, while the optical transmission spectrum shows signatures of VO but a lack of TiO (Evans et al. 2018). The presence of Fe can provide contraints on the interior temperature, T_int, as the element condenses at the highest temperatures for pressures above about 100 bar (see Wakeford et al. 2017).

We explore the constraints on the deep interior temperature (pressures ⪆10 bar) using the best-fit retrieved temperature–pressure (T–P) profile for WASP-121b, derived from fitting the dayside emission spectra of Evans et al. (2017), which is sensitive to approximately bar to mbar pressure levels and has recently been updated to include HST WFC3 G102 and Spitzer data (Mikal-Evans et al. 2019). The T–P profile was generated using the analytical model of Guillot (2010), which assumes radiative equilibrium, and is parameterized with the Planck mean thermal infrared opacity, κ_IR, the ratio of the optical-to-infrared (OIR) opacity, γ, an irradiation efficiency factor β, and the interior temperature of the planet, T_int. Emission spectra are not directly sensitive to T_int (Line et al. 2013b), so retrievals often do not probe the deep interior pressure layers and fix T_int, with Evans et al. (2017; see also Mikal-Evans et al. 2019) setting T_int to Jupiter-like 100 K values. However, with both Mg and Fe present in the upper layers of the atmosphere, temperatures of T_int = 100 K are clearly too low, as the T–P profile intersects Fe and Mg condensation curves at high pressures (see Figure 13). As the condensation curves are metallicity dependent (Visscher et al. 2010; Wakeford et al. 2017), we illustrate the constraints in Figure 13 using 10× solar metallicity, which is close to the best-fit retrieved abundances for WASP-121b in both the transmission and emission spectra (Evans et al. 2017, 2018; see also Mikal-Evans et al. 2019). Varying T_int in steps of 100 K and using the best-fit T–P profile parameters (κ_IR = 0.0049 cm² g⁻¹, γ = 4.08, β = 1.03), we find internal temperatures near T_int = 500 K are needed to prevent both Fe and Mg from condensing deep in the atmosphere, but this also allows TiO to condense at pressures well below the photosphere. As discussed in Fortney et al. (2008), for adiabatic interiors higher values of T_int lead to warmer interiors and larger planet radii. Given WASP-121b is one of the largest exoplanets found (R_pl = 1.865R_Jup; Delrez et al. 2016), a high T_int would be expected. While there are computational limitations, we note that self-consistent 3D modeling efforts that include both the deep interior and the exosphere in comparison to transmission, emission, and phase curve data can provide better constraints on the irradiation efficiency, atmospheric abundances, atmospheric mixing, and global T–P profiles, which in turn would provide the best condensation constraints on T_int.

Zoom In Zoom Out Reset image size

The ionized Mg and Fe lines are seen up to extremely high altitudes, corresponding to high planetary radii of R_pl/R_star ∼ 0.3. We compare the altitudes of this absorption material to the theoretical distances of the L1 Lagrange point, D_L1, and the L1' Lagrange point, D_L1', which is the distance between the planet's center and the closest point of the equipotential surface including L1, which is the Roche lobe. Using the formalism in Gu et al. (2003) and using a mass ratio of M_pl/(M_pl + M_star) = 8.3 × 10⁻⁴ as well as a/R_star = 3.76, we calculate D_L1/R_star = 0.24 (or 1.96R_pl) and the L1' Lagrange size relative to that of the star to be ${D}_{{\rm{L}}1^{\prime} }/{R}_{\mathrm{star}}=0.209$ (or 1.72R_pl). In a transit configuration, the observed limit of the Roche lobe is perpendicular to the planet–star direction and the Roche lobe extends to about 2/3 of the extension to the L1 Lagrange point (Vidal-Madjar et al. 2008). For WASP-121b, we numerically calculated the equivalent radius²¹ of a transiting Roche lobe, R_eqRL, finding R_eqRL/R_star = 0.158 (or 1.3 R_pl).

The core of the Mg ii k line reaches up to R_pl/R_star = 0.309 ± 0.036 (or 2.52 ± 0.29R_pl), which is in excess of R_eqRL/R_star at 4σ confidence. In addition, the ground-state Fe ii resonance line at 2382 Å reaches up to R_pl/R_star = 0.331 ± 0.074 (or 2.7 ± 0.6R_pl), which is in excess of R_eqRL/R_star at 2.7σ confidence. In a transit configuration, these absorption features indicate that both Mg ii and Fe ii reach altitudes such that they are no longer gravitationally bound to the planet. While large hydrogen tails extending well beyond the Roche lobe have been seen at Lyα (Ehrenreich et al. 2015; Bourrier et al. 2018), elements heavier than hydrogen in exoplanets have not previously been found at distances in excess of the Roche lobe.

The existence of heavy atmospheric material reaching beyond the Roche lobe agrees with a geometric blow-off scenario (Lecavelier des Etangs et al. 2004), where the Roche lobe is in close spatial proximity to the planet and the exobase can extend beyond it, letting atmospheric gas freely stream away from the planet. WASP-121b is an ideal planet to potentially observe this phenomenon, as the planet is on the verge of tidal disruption (Delrez et al. 2016). At the terminator, the planet fills a significant portion of its Roche lobe, with the OIR radii reaching R_pl(OIR)/R_eqRL = 77%. The NUV value reaches R_pl(NUV)/R_eqRL = 87 ± 2%, indicating the NUV continuum contains opaque material nearly up to the transit-projected Roche lobe, while the cores of the Mg ii and Fe ii exceed the Roche lobe. A similar process is likely happening in comparably hot and tidally distorted planets such as WASP-12b, which also shows evidence of Mg ii and Fe ii absorption (Fossati et al. 2010; Haswell et al. 2012). However WASP-121b has a more favorable transmission spectral signal and orbits a significantly brighter star than WASP-12A. These properties allow the WASP-121b NUV data to be fit at relatively high resolution with a full limb-darkened, instrument-systematic-corrected transit model, which allows the transmission spectral line profiles to be well resolved and absolute transit depths preserved. Our transit light curves cover ingress, though no early ingress is observed as the NUV transit time is in excellent agreement with the expected ephemeris (see Section 4.3). Like WASP-12b (Nichols et al. 2015; Turner et al. 2016), there is no evidence in WASP-121b for a bow shock or material overflowing the L1 Lagrange point (Lai et al. 2010; Vidotto et al. 2010; Llama et al. 2011; Bisikalo et al. 2013).

As ionized species, the Mg ii and Fe ii atoms are potentially sensitive to the planet's magnetic field, which if strong enough could lead to magnetically controlled outflows (Adams 2011). As such, departures in the transit light curves from spherical symmetry and the velocity profile of the ionized lines could give insights into the nature of the outflow. Our data lack complete phase coverage, so only have a limited ability to constrain a non-spherical absorption profile and we cannot constrain a possible post-transit cometary-like evaporation tail. However, we note that all of our transit light curve fits were well fit to near the photon-noise limit assuming the planet was a sphere. In addition, our model fit found velocities of $-{56}_{-63}^{+43}\,\mathrm{km}\,{{\rm{s}}}^{-1}$, which does not have sufficient precision to confidently distinguish between outflowing blueshifted gas or zero-velocity gas.