The Hubble PanCET Program: A Metal-rich Atmosphere for the Inflated Hot Jupiter HAT-P-41b

Our data analysis procedures follow the general methodology detailed in Nikolov et al. (2014, 2015). We commenced analysis from the flt.fits files, which were reduced (bias-, dark-, and flat-field-corrected) using the latest version of the CALSTIS pipeline and the latest calibration frames. We used median combined difference images to identify and correct for cosmic ray events in the images as described by Nikolov et al. (2014). We found that ∼4% of the detector pixels were affected by cosmic ray events. We also corrected pixels identified by CALSTIS as bad with the same procedure. Together with the cosmic ray identified pixels, this resulted in a total of ∼14% interpolated pixels.

We performed spectral extraction with the IRAF procedure APALL using aperture sizes in the range from 6 to 18 pixels with a step of 0.5. The best aperture for each grating was selected based on the resulting lowest light-curve residual scatter after fitting the white light curves. We found that aperture sizes of 13.5, 13.5, and 10.5 pixels satisfy this criterion for visits 83, 84, and 85, respectively.

We cross-correlated and interpolated all spectra with respect to the first spectrum, to prevent subpixel wavelength shifts in the dispersion direction. The STIS spectra were then used to extract both white light spectrophotometric time series and custom wavelength bands after summing the appropriate flux from each bandpass.

The raw STIS light curves exhibit instrumental systematics on the spacecraft orbital timescale, which are attributed to thermal contraction/expansion (referred to as the "breathing effect") as the spacecraft warms up during its orbital day and cools down during orbital night. We take into account the systematics associated with the telescope temperature variations in the transit light-curve fits by fitting a baseline function depending on various parameters.

White and spectroscopic light curves were created from the time series of each visit by summing the flux of each stellar spectrum along the dispersion axis. We fit each transit light curve using a two-component function that simultaneously models the transit and systematic effects. The transit model was computed using the analytical formulae given in Mandel & Agol (2002), which are parameterized with the mid-transit times (T_mid), orbital period (P) and inclination (i), normalized planet semimajor axis (a/R_s), and planet-to-star radius ratio (R_p/R_*).

Stellar limb-darkening was accounted for by adopting the four-parameter nonlinear limb-darkening law with coefficients c1, c2, c3, and c4, computed using a three-dimensional stellar atmosphere model grid (Magic et al. 2015). We adopted the closest match to the effective temperature, surface gravity, and metallicity values for HAT-P-41 determined by Hartman et al. (2012).

As in our past STIS studies, we applied orbit-to-orbit flux corrections by fitting for a fourth-order polynomial to the spectrophotometric time series phased on the HST orbital period and a linear time term. We also used a low-order polynomial (up to a third degree with no cross terms) of the spectral displacement in the dispersion and cross dispersion directions. The first exposures of each HST orbit exhibit lower fluxes and have been discarded in the analysis. Similar to our past HST STIS analyses (Nikolov et al. 2014; Sing et al. 2016), we intended to discard the entire first orbit to minimize the space craft thermal breathing trend, but found that, for two of the three HST visits, a few of the exposures taken toward the end of the first orbit can be used in the analysis.

We then generated systematics models that spanned all possible combinations of detrending variables and performed separate fits using each systematics model included in the two-component function. The Akaike information criterion (AIC; Akaike 1974) was calculated for each attempted function and used to marginalize over the entire set of functions following Gibson (2014). Our choice to rely on the AIC instead of the Bayesian information criterion (BIC; Schwarz 1978) was determined by the fact that the BIC is more biased toward simple models than the AIC. The AIC therefore provides a more conservative model for the systematics and typically results in larger or more conservative error estimates, as demonstrated by Gibson (2014). Marginalization over multiple systematics models assumes equal prior weights for each model tested.

For the white light curves, we fixed the orbital period, inclination, and a/Rs to the values reported in Table 5 and fit for the transit mid-time and planet-to-star radius ratio. We find central transit times Tc[MJD] = 58000.6958 ± 0.0029 (visit 83), Tc[MJD] = 58245.85414 ± 0.00039 (visit 84), and Tc[MJD] = 58280.87484 ± 0.00036 (visit 85). We derive the white light transit depths to be 10200 ± 104 ppm and 10320 ± 85 ppm for G430L and G750L, respectively.

Table 5.  Transit Parameters for HAT-P-41b

Parameter	Value
Rp/Rs	0.1028 ± 0.0016
a/Rs	${5.44}_{-0.15}^{+0.09}$
i [Degrees]	87.7 ± 1.0
Tc [BJD]	2454983.86167 ± 0.00107
$\mathrm{log}{g}_{p}$ [cgs units]	2.84 ± 0.06
P [days]	2.694047 ± 4 × 10−6

Note. All values are from the discovery paper by Hartman et al. (2012).

Download table as:  ASCIITypeset image

For the spectroscopic light curves, a common-mode systematics model was established by simply dividing the white light curve by a transit model (Berta et al. 2012; Deming et al. 2013). We computed the transit model using the orbital period, inclination, and a/Rs from Table 5, and the central times for each orbit from the white light analysis. The common-mode factors from each night were then removed from the corresponding spectroscopic light curves before model fitting.

We then performed fits to the spectroscopic light curves using the same set of systematics models as in the white light curve analysis and marginalized over them as described above. For these fits, Rp/Rs was allowed to vary for each spectroscopic channel, while the central transit time and system parameters were fixed. We assumed the nonlinear limb-darkening law with coefficients fixed to their theoretical values, determined in the same way as for the white light curve. The detrended spectrophotometric light curves are shown in Figures 19–21. The derived STIS transit spectrum is shown along with the entire transit spectrum in Table 6.

Table 6.  Transit Spectrum

Instrument	Wavelength (μm)	Depth (ppm)a	Instrument	Wavelength (μm)	Depth (ppm)
STIS G430Lb	0.290–0.350	10091 ± 230	STIS G750L	0.711–0.731	10499 ± 364
	0.350–0.370	10006 ± 249		0.731–0.750	10038 ± 302
	0.370–0.387	10397 ± 198		0.750–0.770	10142 ± 224
	0.387–0.404	10273 ± 163		0.770–0.799	10124 ± 281
	0.404–0.415	9999 ± 137		0.799–0.819	10618 ± 317
	0.415–0.426	9980 ± 185		0.819–0.838	9910 ± 239
	0.426–0.437	10324 ± 145		0.838–0.884	10252 ± 238
	0.437–0.443	10428 ± 287		0.884–0.930	10217 ± 298
	0.443–0.448	10245 ± 168	WFC3c	1.122–1.141	10297 ± 107
	0.448–0.454	10411 ± 165		1.141–1.159	10620 ± 115
	0.454–0.459	10367 ± 192		1.159–1.178	10347 ± 130
	0.459–0.465	10227 ± 145		1.178–1.196	10479 ± 113
	0.465–0.470	10640 ± 164		1.196–1.215	10497 ± 111
	0.470–0.476	10429 ± 165		1.215–1.233	10644 ± 111
	0.476–0.481	10453 ± 144		1.233–1.252	10289 ± 100
	0.481–0.492	10584 ± 125		1.252–1.271	10360 ± 107
	0.492–0.498	10224 ± 192		1.271–1.289	10488 ± 122
	0.498–0.503	10289 ± 166		1.289–1.308	10405 ± 97
	0.503–0.509	10459 ± 149		1.308–1.326	10396 ± 94
	0.509–0.514	10519 ± 183		1.326–1.345	10581 ± 92
	0.514–0.520	10506 ± 168		1.345–1.364	10684 ± 105
	0.520–0.525	10429 ± 186		1.364–1.382	10622 ± 113
	0.525–0.531	10558 ± 128		1.382–1.401	10477 ± 112
	0.531–0.536	10247 ± 174		1.401–1.419	10689 ± 98
	0.536–0.542	10451 ± 180		1.419–1.438	10535 ± 108
	0.542–0.547	10476 ± 177		1.438–1.456	10626 ± 113
	0.547–0.552	10422 ± 206		1.456–1.475	10564 ± 121
	0.552–0.558	10794 ± 153		1.475–1.494	10686 ± 113
	0.558–0.563d	9221 ± 279		1.494–1.512	10799 ± 129
	0.563–0.569	10444 ± 188		1.512–1.531	10551 ± 124
STIS G750Le	0.526–0.555	10356 ± 220		1.531–1.549	10566 ± 111
	0.555–0.575	10497 ± 278		1.549–1.568	10515 ± 134
	0.575–0.594	10317 ± 202		1.568–1.587	10436 ± 110
	0.594–0.614	10061 ± 236		1.587–1.605	10492 ± 145
	0.614–0.633	10386 ± 142		1.605–1.624	10340 ± 122
	0.633–0.653	10542 ± 219		1.624–1.642	10338 ± 142
	0.653–0.672	10418 ± 276		1.642–1.661	10331 ± 138
	0.672–0.692	10182 ± 225	Spitzer IRAC1	3.2–4.0	10191 ± 102
	0.692–0.711	10168 ± 192	Spitzer IRAC2	4.0–5.0	10679 ± 145

Notes.

^aTransit depth $=\,{R}_{\mathrm{planet}}^{2}/{R}_{\mathrm{star}}^{2}$. ^bTypical STIS G430L bin size = 0.0055 μm (median resolution ∼350). ^cWFC3 bin size = 0.0186 μm (median resolution ∼75). ^dOutlier bin strongly affected by systematics and ignored in retrieval analyses. ^eTypical STIS G750L bin size = 0.0196 μm (median resolution ∼130).

Download table as:  ASCIITypeset image

The WFC3 data reduction generally follows the methodology of Sheppard et al. (2017), using programs adapted from IDL for Python. We download the "ima" data files from the HST archive, and remove background contamination following the "difference reads" methods of Deming et al. (2013), which allows us to easily resolve and remove potential contamination from the nearby companion (Evans et al. 2016). We determine a wavelength solution by taking the zero point from the F140W photometric observation and fitting for the wavelength coefficients that allow an out-of-transit spectrum to match the appropriate ATLAS stellar spectrum (Castelli & Kurucz 2003).

We then divide the background-subtracted "ima" science frame by the WFC3 flat-field calibration file, and return the dark-, bias-, and flat-field-corrected flux array, in units of electrons. The uncertainty of the flux at each pixel is taken from the "ima" file's error frame, which accounts for gain-adjusted Poisson noise, read noise, and noise from dark current subtraction. This is further adjusted via error propagation for the added uncertainty from background removal and flat-field correction. We use the "ima" file's data quality frame to mask (i.e., give zero weight to) pixels that are flagged as bad in every exposure in the time series. We then correct for cosmic rays using a conservative time series sigma cut of 8σ, while accounting for changes in flux that occur due to uneven scan rates and the transit itself, and set the affected pixels to the median value of that pixel in the time series. The average amount of pixels either impacted by cosmic ray events or flagged as bad pixels is 1.8% of all pixels in an exposure. The reduced exposures are summed over the spatial scan direction to give a 1D spectrum at each observation time.

We use a marginalization light-curve analysis similar to that of Sheppard et al. (2017), applied to transit curves. This is a Bayesian model averaging method, first described by Gibson (2014) and applied by Wakeford et al. (2016), with further detrending by use of band-integrated (white light) residuals in spectral light-curve fitting (Mandell et al. 2013; Haynes et al. 2015). We first analyzed the band-integrated light curve to simplify the spectral light-curve analysis, then we fit the spectral light curves to derive the WFC3 transit spectrum.

We model the observed light curve as a BATMAN¹⁵ transit model (Kreidberg 2015) combined with an instrumental systematic model. For the transit model, we assume nonlinear limb darkening and derive the coefficients by interpolating the 3D values from Magic et al. (2015) to the central wavelength of WFC3 (1.4 μm). The limb-darkening derivation is consistent with that used in the STIS analysis. We only fit for transit depth and central transit time, since the incomplete coverage of HST makes it difficult to improve constraints on other transit parameters, such as a/R_s or inclination. We fix these values in the light curve analyses of each instrument (STIS, WFC3, and Spitzer), which ensures consistent orbital parameters are used when analyzing different data sets. The transit and system parameters are shown in Table 5.

The systematic model grid comprises a set of 50 polynomial models, which account for a visit-long linear slope, HST orbital phase-dependent systematics (i.e., "breathing," ramp), and subpixel wavelength shifts (for full model grid, see Wakeford et al. (2016)). It is computationally difficult to fully sample the parameter space of all 50 models using Markov Chain Monte Carlo (MCMC) samplers, so we instead fit each model using KMPFIT¹⁶ (Terlouw & Vogelaar 2015), a Python implementation of the Levenberg–Markwardt least squares minimization algorithm, to more quickly determine parameter values and uncertainties. Wakeford et al. (2016) found that uncertainties derived from these two methods typically agree within 10%. We then weight each model by its Bayesian evidence—approximated from AIC—and marginalize over the model grid (assuming a prior that each model is equally likely) to derive the light-curve parameters and uncertainties while inherently accounting for uncertainty in model choice.

As is common practice, we ignore the systematic-dominated first orbit in the white light analysis; however, the use of common-mode detrending allows us to include that orbit in the spectral light-curve analysis. The raw light curve, the light curve with instrumental systematics removed, and the residuals from the highest-weight systematic model are shown in Figure 1. We derive the white light depth to be 10490 ± 51 ppm.

Zoom In Zoom Out Reset image size

To derive the transit spectrum, we bin the 1D spectra from each exposure between the steep edges of the grism response curve (1.12–1.66 μm), deriving a flux time series for each spectral bin. We use bins of width 0.0186 μm (4 pixels) to maximize resolution without drowning the signal in noise. We note that the atmospheric retrieval is not sensitive to the choice of bin size. We fit each spectral light curve using a model similar to that used for the band-integrated light curve, to derive the transit depth at each wavelength.

The spectral light-curve analysis mimics white light analysis with a few exceptions: all system parameters except transit depth are fixed to the white light values, limb-darkening coefficients are interpolated to the bin's central wavelength, and residuals from the band-integrated light curve are incorporated as a scalable parameter. The band-integrated residuals encode possible instrumental systematics that were not captured by the model grid, and so they can be used to further detrend the binned spectral light curves. The shapes of the residuals are assumed to be constant with wavelength, though the amplitude is allowed to vary. This allows for removal of any wavelength-independent red noise from spectral bin curves, at the penalty of slightly increasing the white noise. Note that the band-integrated uncertainty is sufficiently small, relative to spectral light-curve uncertainty, that the added noise has only a minor effect. The wavelength range, transit depth, and depth uncertainty for each WFC3 bin are shown in Table 6.

The shape of our derived spectrum is in excellent agreement with the literature spectrum (Tsiaras et al. 2018), though it is shifted to higher depths by ∼90 ppm, indicating that we derive a larger white light depth. This difference persists even if we derive the spectrum without using white light residuals. This could plausibly be due to different limb-darkening treatments or different systematic modeling choices. This difference emphasizes the importance of considering offsets between instruments in retrieval analyses (Section 5.1.2). Further, the white light depth is subject to the choice of orbital parameters, which are typically fixed in spectroscopic light-curve fits. We run sensitivity tests to determine that accounting for orbital parameter uncertainties increases the scatter between HST STIS and HST WFC3 by roughly 60 ppm for HAT-P-41b. We capture this scatter by including it in WFC3's depth uncertainty, increasing it from 50 to 80 ppm.

As a check, we performed retrievals using the published spectrum from Tsiaras et al. (2018) in combination with the derived STIS and Spitzer depths, and found differences well within the 1σ uncertainties for the retrieved parameters. The major results and conclusions in this paper are not sensitive to the WFC3 spectrum choice.

Marginalization is only reliable if at least one model is a good representation of the data (Gibson 2014; Wakeford et al. 2016). We therefore checked the goodness-of-fit of the highest-weighted systematic model for each light curve using both reduced χ² and residual normality tests. Further, we explored whether red noise is present in the light-curve residuals, as that can bias inferred depth accuracy and precision (Cubillos et al. 2017).

Though χ² cannot prove that a model is correct, it can demonstrate that the fit of a particular model is consistent with that of the "true" model with "true" parameter values (Andrae et al. 2010). Therefore, it is an informative goodness-of-fit diagnostic, and it is particularly useful due to its familiarity and simplicity. The "true" model with "true" parameters will have a reduced χ² of one, with uncertainty defined by the χ² distribution. For both the band-integrated and spectral light curves (66 and ∼88 degrees of freedom, respectively), this results in an acceptable reduced χ² range of roughly 0.7–1.3.

The band-integrated analysis (${\chi }_{\nu }^{2}=1.17$) and all but one of the spectral bins (median ${\chi }_{\nu }^{2}=0.9)$ fall within this range. The exception is the 1.299 μm light curve, for which the highest-weighted model fit has a reduced χ² of 0.56. This low value indicates that the uncertainties in this light curve are overestimated. This is likely due to incorporating white light residuals, which both inflate uncertainties and can potentially interpret random white noise as structure. However, it is not flagged by the normality or correlated noise analyses (described below), and fitting the light curve without incorporating white light residuals finds a consistent depth with a more reasonable ${\chi }_{\nu }^{2}=0.71$. Further, the derived depth and uncertainty exhibit good agreement with the published (Tsiaras et al. 2018) transit depth at this wavelength (accounting for the white light offset). Therefore, we include it in the transit spectrum. For the other 28 spectral bins, the reduced χ² values provide no evidence against the validity of the derived transit depths and uncertainties.

A residual normality test checks whether the residuals for a model are Gaussian-distributed in order to determine goodness of fit, as this is expected for the "correct" model. Like reduced χ², a normality test cannot prove that a model is correct; rather, it can only diagnose incorrect models. We use the scipy implementation of the common Shapiro–Wilk test for normality (Shapiro & Wilk 1965), and determine for which light curves the highest-evidence model has normality ruled out at the 5% significance level. At a sample size of around 90, this is by no means rigorous, but it is still a useful heuristic for flagging potentially problematic light-curve models.

Normality is rejected at the 5% significance level only for the band-integrated residuals and the 1.243 μm spectral bin residuals. In both cases, normality is ruled out due to a single outlier in the time-series. Normality is rejected at 3% significance for the band-integrated residuals, due entirely to the first exposure in the first orbit. When this exposure is ignored, we recover a consistent depth and uncertainty, and the residuals are consistent with normality. We therefore keep this exposure in the analysis.

A possible cause of the spectral bin's outlier is a minor cosmic ray event or bad pixel that was small enough to both avoid the detection by the sigma cut and not affect the band-integrated curve, but large enough to impact the much smaller bin flux. Removing this spectral bin from the retrieval had no noticeable effect on the results. Further, the derived depth and uncertainty are consistent with literature values (Tsiaras et al. 2018). We therefore decide to include this bin in the retrieval.

Finally, we test for correlated noise in the residuals following the methodology of Cubillos et al. (2017) (also see Pont et al. 2006) and using MC³.¹⁷ Though this method is not necessarily rigorous for HST, due to the incomplete phase coverage and relatively small number of exposures, it is still a practical diagnostic. We find no evidence of correlated noise up to the timescale of half an HST orbit (Figure 2). Beyond this timescale (nine exposures per bin) the standard deviation is based on a light curve with fewer than eight points, and so the relationship is less informative and more subject to small number statistics.

Zoom In Zoom Out Reset image size

With the caveats noted above, marginalization does an excellent job in fitting the spectral light curves. We emphasize that removing any of the flagged bin spectra has no effect on the retrieval. Together, these tests support the validity of the derived transit depths and uncertainties in the WFC3 bandpass.

The Spitzer data consist of cubes of 64 subarray frames in each band, each of size 32 × 32 pixels. We extracted aperture photometry for each frame, totaling 21,632 frames at both 3.6 and 4.5 μm. To extract photometry, we used 11 numerical apertures with radii ranging from 1.6 to 3.5 pixels, and we centered those apertures on the position of the host star determined using both a 2D Gaussian fit to the stellar point-spread function and also a center-of-light calculation. Since HAT-P-41 has a companion star 36 distant (Hartman et al. 2012; Evans et al. 2016), we measured the flux of the companion scattered into each of our numerical apertures, using the method described by Garhart et al. (2020). We adopted the magnitude difference in the Spitzer bands as deduced by Garhart et al. (2020). Accounting for our having used different aperture radii than Garhart et al. (2020), we derive dilution correction factors of 1.0171 and 1.0106 at 3.6 and 4.5 μm, respectively. The transit depths are then multiplied by those factors in order to correct for the presence of the companion star.

We fit transit curves to the 22 sets of photometry at each wavelength (eleven apertures, each with two centering methods). Our default fitting procedure fixes the orbital parameters at the values in the discovery paper by Hartman et al. (2012), fitting only for the central time and depth of the transit. The shapes of the Spitzer transits are well-matched when fixing the orbital parameters to those values. However, we also explored including the orbital inclination and a/R_s in the fit (see below), those being the orbital parameters that most strongly affect the shape of the transit. We adopt quadratic limb-darkening coefficients calculated by least squares for the Spitzer bands by Claret et al. (2013), using 2 km s⁻¹ microturbulence. We choose the values for T_eff = 6400 K and log g = 4.0, without interpolation. We fix those coefficients in the fitting process, and we deem these choices to be appropriate given that the limb darkening is minimal at these infrared wavelengths. Our fits to the transit account for the intrapixel sensitivity variations of the Spitzer photometry using pixel-level decorrelation (PLD; Deming et al. 2015), including a linear baseline (ramp) in time. We use a Bayesian information criterion to decide between a linear versus quadratic ramp. The details of the PLD fit are the same as described for secondary eclipses by Garhart et al. (2020), except that we are fitting transits, so we include limb darkening. Briefly, the fitting code bins the photometry and pixel basis vectors to various degrees, and selects the optimal bin size, aperture radius, and centering method, based on the smallest difference from an ideal Allan deviation relation (Allan 1966). The Allan deviation relation expresses that the standard deviation of the residuals should decrease as the square root of the bin time. Operating on binned data allows the PLD algorithm to concentrate on the longer timescales that characterize the red noise (and also the transit duration), as opposed to the 0.4 s cadence time of the raw photometry.

We determine the errors on our transit depths and times using an MCMC procedure, with a burn-in phase of 10,000 steps, followed by 800,000 steps to explore parameter space. We calculate multiple chains for each transit, and verify convergence using a Gelman–Rubin (GR) statistic (Gelman & Rubin 1992). Our GR values for our PLD fits are very close to unity, being 1.0027 at 3.6 μm and 1.0004 at 4.5 μm, indicating good convergence. Our transit depths and times are listed in Table 7.

Table 7.  Transit Times and Depths for HAT-P-41b in the Spitzer Bands

Wavelength	BJD(TDB)	${R}_{p}^{2}/{R}_{s}^{2}$ (ppm)
3.6 μm	2457772.20440 ± 0.00033	10020 ± 100
4.5 μm	2457788.36860 ± 0.00032	10568 ± 135

Note. These are "as observed" transit depths, not corrected for dilution by the companion star. To correct for dilution, multiply the values of ${R}_{p}^{2}/{R}_{s}^{2}$ by 1.0171 at 3.6 μm, and by 1.0106 at 4.5 μm. The corrected values are shown with the rest of the transit spectrum in Table 6.

Download table as:  ASCIITypeset image

Our derived transit times are in excellent agreement with measurements of the same Spitzer transits by Wakeford et al. (2020). Specifically, using our uncertainties, our fitted times differ by 1.1σ and 0.6σ at 3.6 and 4.5 μm, respectively. Wakeford et al. (2020) use eight HST transits, as well as the discovery results and the Spitzer transits, to derive a new ephemeris. Our fitted times (Table 7) differ from that ephemeris by insignificant amounts (0.2 and 7.8 s).

We explored the effect of uncertainties in the orbital parameters, since those can affect the derived transit depth (Alexoudi et al. 2018). Adopting uniform priors for a/R_s and inclination, we find that they are degenerate when fitting only the Spitzer transits. That degeneracy is illustrated in Figure 3, where it is apparent that inclination and a/R_s can trade off to maintain the observed transit duration and the sharp ingress/egress that characterizes the Spitzer transits. Changing the inclination changes the chord length across the stellar disk, and (when limb-darkening is minimal) that can be compensated by changing a/R_s to maintain the same transit duration. The orbital solution from Hartman et al. (2012) is entirely consistent with our likelihood distribution for those parameters, as shown in Figure 3 for 3.6 μm (4.5 μm is similar). We therefore freeze the orbital parameters at the Hartman et al. (2012) values when fitting our Spitzer transits.

Zoom In Zoom Out Reset image size

Figure 4 illustrates the transits at 3.6 and 4.5 μm. The residuals from the best-fit model are included in the figure, and the right panel shows the residuals binned over increasing timescales, a so-called Allan deviation relation (Allan 1966). The slopes of those relations are close to the −0.5 value expected for photon noise.

Zoom In Zoom Out Reset image size

The partial cloud parameter is motivated by work by Line & Parmentier (2016) and MacDonald & Madhusudhan (2017), which showed that if the gray cloud deck were inhomogeneous, then the spectrum we observe (D) would be a weighted average of the clear atmosphere transit spectrum (D_clear) and the cloudy atmosphere spectrum (D_cloudy), with weights given by the cloud fraction (f_c). We implement this as D = f_c ∗ D_cloudy + (1 − f_c) ∗ D_clear. Since high-altitude gray clouds are seen in spectra as flat lines, averaging a spectrum with features with this line will shrink the features and can mimic the effect of an atmosphere with high mean-molecular mass and small scale height. Thus, including this parameter allows us to account for this possible degeneracy and prevents us from overconfidently claiming a high-metallicity atmosphere.

The instrumental offsets are nuisance parameters that can capture the extent to which transit depths from STIS G430L, STIS G750L, or WFC3 are biased relative to the depths from the other instruments. This is motivated by the use of common-mode corrections in the light-curve analysis of each instrument, which can potentially introduce a uniform bias for the depth at each spectral bin for that instrument. Offsets are also able to account for inter-instrumental transit depth scatter introduced by uncertainty in orbital parameters a/R_s and inclination (Section 4.2.2).

We explore three offset scenarios. The first is physically motivated. In this scenario, we use Gaussian priors with sigmas determined by the uncertainties in the band-integrated transit depths to try to reflect the correlated uncertainty that exists between spectral bins for each instrument, whether due to common mode corrections or orbital parameter uncertainties. This essentially propagates the white light depth uncertainty into the retrieval. Derived in Section 4, these uncertainties are 105 ppm, 85 ppm, and 80 ppm for STIS G430L, STIS G750L, and WCF3, respectively. The second scenario investigates the potential impact of unknown sources of bias by setting a large, uniform offset prior for each instrument's offset. The third scenario extends this by setting two large, uniform priors: a WFC3 offset and a single offset for both STIS instruments. The third scenario allows the absolute depths at STIS to vary while preserving the optical spectral shape.

We caution that offsets—especially penalty-free, uniform prior offsets—can cloak missing physics in a model. We do not think they should act as a safety net to achieve a good fit to a spectrum, and the inferred atmospheric properties should be understood in context. However, offsets offer a way to both investigate potential instrumental biases and account for absolute depth uncertainty for each instrument. It is valuable to include offsets as model parameters and marginalize over these possible values in order to understand how the uncertainty in the absolute transit depth for each instrument affects the marginalized posterior distributions of the other model parameters.

Line & Parmentier (2016) showed that partial cloud coverage (i.e., clouds at the same height but not uniformly covering the limb azimuthally) could mimic the effect of a high mean molecular mass atmosphere for WFC3 spectra. When partial clouds are present, the observed spectrum would be the weighted average of the cloudy and clear spectra. The transit depth of an atmosphere dominated by gray clouds does not vary with wavelength, and so the cloudy spectrum is a straight line. Averaging a clear spectrum with molecular features and a cloudy, flat spectrum reduces the size of the features by an amount proportional to the cloud fraction. Given that we find a significantly supersolar metallicity, we investigate whether this possible mean molecular mass–cloud fraction degeneracy affects the results from the fiducial case.

When fit independently and allowing the cloud fraction to vary, both WFC3 and STIS spectra retrievals no longer constrain the metallicity to be supersolar. However, fitting the infrared and optical data jointly breaks this degeneracy, as predicted by Line & Parmentier (2016). Effectively, a low mean molecular mass and nonuniform cloud solution should be impacted by Rayleigh scattering in the optical data, especially the bluest six wavelength bins. The dominance of gas absorption opacity over Rayleigh scattering opacity in the STIS data disallows this solution, breaking the degeneracy in favor of the high mean molecular mass solution.

However, it is possible that removing assumptions made in the fiducial model—such as fixed Rayleigh scattering or no instrumental offsets—could muddle this decisive degeneracy break and allow for a low-metallicity solution. We investigate this below.

The fiducial model assumes Rayleigh scattering. In lieu of complicated microphysics, PLATON allows parametric scattering, in which the slope and the magnitude of Rayleigh scattering vary in order to capture the possible signature of many hazes. For a more detailed explanation, see Zhang et al. (2019). Though there is no obvious signature of haze in the optical data (i.e., no linear slope decreasing with increasing wavelength), it is worth exploring whether loosening the assumption of exact Rayleigh scattering affects the results.

Allowing the full scattering parameter space (see Table 10) has little effect: the clear lack of slope in the STIS data conclusively leads to a haze-free atmosphere. Further, the median scattering factor is 0.01, implying that the data is easiest to fit when opacity from Rayleigh scattering is muted. This complicates the mean molecular weight–cloud fraction (μ − f_c) degeneracy. Lower values of μ are now possible, since the model no longer expects scattering opacity to be important at optical wavelengths. The lower the magnitude of Rayleigh scattering—and the shallower the scattering slope—the lower μ can be. This is because decreases in Rayleigh opacity allow for gas absorption to still be dominant at larger scale heights. As a result, a patchy cloud and low-metallicity solution is viable. Though possible, the low μ solution requires a specific combination of cloud-top pressure, scattering slope, scattering factor, and cloud fraction, and does not improve the fit. Therefore, it is much less likely than the high-metallicity solution. The marginalized posterior probability distribution for metallicity has the same maximum likelihood value as the fiducial model. The difference is that the distribution has a tail extending to lower metallicities (Figure 7). The resulting median log metallicity and 1σ range (as determined by the 16th and 84th percentile values) is ${\mathrm{log}}_{10}Z/{Z}_{\odot }={2.34}_{-0.64}^{+0.27}$.

Zoom In Zoom Out Reset image size

Though all cloud fractions and cloud pressures are allowed, the posterior is consistent with a clear atmosphere due to the likelihood desert in the upper left corner of the cloud fraction–cloud-top pressure pairs plot: clouds are only seen above the altitude corresponding to the ∼10 Pa pressure level at fractions below 0.50.

We model an instrumental offset as a constant value added to the forward model's binned transit depth in the wavelength range of the instrument of interest.

For the physically motivated scenario (Scenario 1 from Section 5.1.2), we set the priors for STIS G430L, STIS 750L, and WFC3to be Gaussians centered on zero ppm, with widths set to the uncertainty on the transit depth from their white light curves (105, 85, and 80 ppm, respectively; see Table 10).

The retrieved median WFC3 offset is nontrivial, with a median of about 1.5× the white light uncertainty (130 ± 50 ppm). The offsets in the STIS G430L and G750L are less significant, at 58 ± 58 ppm and −65 ± 55 ppm, both well within their white light uncertainties. However, there is a significant median offset of ∼120 ppm between the two instruments. This is driven primarily by the retrieval attempting to align transit depths in the overlapping wavelength region between the instruments.

The Spitzer 3.6 μm point drives the WFC3 offset: shifting the WFC3 depths down necessitates a smaller radius ratio, which better captures the relatively low transit depth at 3.6 μm. The ability to better capture the Spitzer 3.6 μm results in a greater Bayesian evidence, indicating that this model is strongly preferred over the fiducial model (for a more detailed discussion, see Section 6.3).

When combined with partial clouds, the instrumental offsets cause a small decrease in the retrieved median metallicity (${\mathrm{log}}_{10}Z/{Z}_{\odot }={2.33}_{-0.25}^{+0.23}$). This is for reasons similar to those explained in the parametric scattering section: whereas parametric scattering justified the absence of an optical scattering slope in the low-metallicity solution by effectively removing Rayleigh scattering opacity, the instrumental offset model can decrease the WFC3 depths relative to the STIS depths in order to artificially allow for it.

Note that increasing STIS depths (instead of decreasing WFC3 depths) has the same effect on Rayleigh opacity and thus metallicity. However, it is not a viable solution: this is because, unlike decreasing WFC3 depths, it does not improve the forward model's ability to capture the low Spitzer 3.6 μm point.

Allowing Gaussian-prior instrumental offsets had no significant effect on the results. However, it is not impossible that there is some unknown wavelength-independent systematic that biases the absolute transit depths of the instruments relative to one another. Though it is unlikely, we chose to explore this by allowing offsets in the STIS G430L, STIS G750L, and WFC3 data to vary by about 5% (500 ppm) in either direction (Scenario 2 from Section 5.1.2). Due to the model preferring a lower radius ratio to best explain the Spitzer 3.6 μm point, the median WFC3 offset is a 250 ppm decrease, about 3× the white light uncertainty. Surprisingly, this large offset does not significantly change the 1σ ranges for metallicity (${\mathrm{log}}_{10}Z/{Z}_{\odot }={2.26}_{-0.40}^{+0.24}$). The size of the molecular features and the large differential between Spitzer photometric points drive the supersolar metallicity in this case. Regardless of the magnitude of the offset, we retrieve a high metallicity. The two large, uniform offsets case, where both STIS instruments are offset by the same amount (Scenario 3 from Section 5.1.2), retrieves posterior distributions effectively identical to those of Scenario 2.

While it is worthwhile to understand the effect on the retrieval, there is no reason to expect such large instrumental offsets for HAT-P-41b. A transit depth offset can be caused by the necessity of analyzing each instrument differently. For example, not handling limb darkening consistently and not using consistent orbital parameters (i.e., inclination) for each analysis might cause an offset, but this is easily fixed and is not an issue for our data set. Since the instruments' observations are from different dates, it is also possible that stellar variability could cause an offset. However, we have long-term photometry (Section 2.3) that shows no such variability. Additionally, the STIS depths are in good agreement with HST UVIS observations (Wakeford et al. 2020). There is no indication that this particular observation is biased in any way, and unresolved companions are confidently ruled out (Evans et al. 2018). The most plausible source is unaccounted-for uncertainties or bias from the spectral analysis, as the WFC3 spectrum derived in this paper is shifted up ∼90 ppm relative to the literature spectrum (Tsiaras et al. 2018), as noted in Section 4.2.2. However, this is still well below the 250 ppm value preferred by the large, uniform offset models. We determine that offsets beyond the physically motivated values are unlikely.

Section 3.1 demonstrated that HAT-P-41 is consistent with a quiet star and stellar activity is not expected to impact the transit spectrum. However, to be conservative, we investigated whether allowing for greater stellar variability impacted our conclusions.

The typical signature of unocculted, cool starspots is to mimic a haze-like slope in the transit spectrum, and such a signature is clearly absent in the derived transit spectrum of HAT-P-41b. On the other hand, the signature of hot faculae is a steep optical drop-off toward shorter wavelengths (Rackham et al. 2019). Given that we see a drop in transit depths in the optical, the retrieval could plausibly be affected if faculae dominate over star spots, and so we focus on a faculae overabundance.

We assumed that the temperature of the stellar photosphere equals the stellar effective temperature. We modeled the faculae following the prescription from Rackham et al. (2019), and accordingly fixed the faculae temperature to T_phot + 100 K. PLATON weights the contributions from the different temperature regimes via the fractional coverage parameter, which represents the overabundance of faculae in the unocculted regions. Rackham et al. (2019) states that moderately active F5V-dwarfs will have around 1% faculae coverage, and up to about 7% on the more active end. This is the faculae fraction, which is much higher than the faculae overabundance. However, we set a conservative uniform prior on the fractional coverage of 0%–10% in order to determine whether high activity would significant alter our conclusions.

We find that including stellar activity has no effect on the posterior probability distribution. It may seem that the STIS data could be explained by a featureless flat line and stellar activity instead of a TiO feature. However, the overabundance of faculae necessary to explain the drop in the bluest six points (0.32–0.42 μm) produces a poor fit to the rest of the STIS data. Therefore, even when including stellar variability, TiO is necessary to explain the STIS depths. Allowing a wider range of faculae temperatures also had no effect.

In summary: there is no evidence of stellar variability from prior observations, and allowing for activity does not affect the retrieval.

Benneke et al. (2019) recently invoked Mie scattering to explain anomalously low Spitzer transit depths. Given the relatively low value of HAT-P-41b's Spitzer 3.6 μm depth relative to the rest of the spectrum—the fiducial model's predicted depth at 3.6 μm is about 3.3σ away from the observed depth—we included Mie scattering in our analysis. In PLATON, Mie scattering can be used in lieu of parametric Rayleigh scattering.

Each condensate is described by a wavelength-dependent complex refractive index, n-ik, where n is the real part and k is the imaginary part of the index. This index explains how that particular condensate interacts with light with wavelengths similar to the particle size. PLATON assumes log-normal distribution in particle size, with geometric standard deviation 0.5, to determine the abundance of different radii condensates for a given mean particle size (Zhang et al. 2019). The other relevant factors are cloud height (condensates are only relevant above that pressure; below it, gray cloud opacity dominates), particle density at the cloud-top pressure, and condensate scale height. The condensate scale height is parameterized as a fraction of the gas scale height, and it describes how the abundance of Mie scattering particles decreases with height. The refractive index is fixed for a given condensate, and the other four parameters are fit for in the retrieval (see Table 10).

PLATON tests one Mie scattering species at a time. Only a few species expected to form clouds in hot Jupiter atmospheres have condensation temperatures above HAT-P-41b's limb temperature (T ∼ 1600 K; Wakeford & Sing 2015). Those five (SiO₂, Al₂O₃, CaTiO₃, FeO, and Fe₂O₃) fall into two phenotypes: "low-n" with real refractive index n ≈ 1.5, and "high-n" with n ≈ 2.5. Though the k values vary more significantly, we find that they do not have a significant impact on the absorption cross section of the condensates. We use n and k values from Kitzmann & Heng (2018), which we average over our wavelength range (0.3–5 μm). The n values are flat over this range, and so the average is an excellent approximation. We tested retrievals with both the low-n (corundum; Al₂O₃) and high-n (hematite; Fe₂O₃) phenotypes.

The priors for the fittable parameters are shown in Table 10. The prior for cloud-top pressure is the same as the fiducial model. Since the condensate radii must be such that they cause a relative drop in opacity around 4 μm (i.e., they increase opacity in the optical by more than in the mid-infrared), we can constrain the mean particle size reasonably well. We set the prior to be log-uniform with a range that contains all plausible values. The number density is not known ahead of time, so we set an uninformative log-uniform prior; widening the prior further did not affect the retrieval. Finally, it is unclear what physical constraints there are on condensate scale height. Fortney (2005) finds that condensate scale heights can be one-third of the gaseous scale height for hot Jupiters, and Benneke et al. (2019) found H_part/H_gas ≈ 3 for a sub-Neptune. Using these values as guides, we set a conservative uniform prior on the fractional scale height and constrain it to be in the range 0.1–10.

Including Mie scattering opacity does not noticeably affect the results of the retrieval, and the Mie scattering parameters are not constrained by the retrieval. The inferred small gaseous scale height—which dampens features and is necessary to explain STIS and WFC3 feature sizes—makes it difficult to explain the large variations in the radius ratio. Combining Mie scattering with partial clouds—physically, an atmosphere with patchy clouds and Mie scattering particles distributed only above those clouds—alleviates the issue of explaining the large transit depth variation. Since partial clouds allow for higher scale heights, Mie scattering could then cause a larger drop in transit depth near Spitzer without needing to invoke an unreasonably high fractional scale height.

Figure 8 shows the corner plot for this model. Though the Mie scattering parameters are not constrained, at number densities above ∼10⁸ m⁻³, particle radii around 0.15 μm, and condensate scale heights greater than the gaseous scale height, lower metallicity and temperature values are possible. This is because added Mie opacity tends to mute features near its peak opacity. This provides a physical reason to expect smaller spectral features, and so a less small scale height is necessary to fit the features. The net impact is a decreased—but still supersolar—median metallicity of ${\mathrm{log}}_{10}Z/{Z}_{\odot }={2.27}_{-0.55}^{+0.30}.$

Zoom In Zoom Out Reset image size

Instead of model selection, it is possible to take a weighted average of the results from each model and therefore automatically take their respective evidences into account (Gibson 2014; Wakeford et al. 2016, 2018). The benefit of Bayesian model averaging is the ability to quantify uncertainty in model selection, as well as avoiding having to arbitrarily choose between models with slightly different evidences. However, it requires a few assumptions: it is only valid if the set of models comprises the full model space, i.e., at least one model is a good description of the data. The weight-averaged uncertainties assume Gaussian-distributed posteriors, which is not strictly correct. However, it is useful in combining information from every model.

Here, we show the assumptions we make to use Bayesian model averaging. The ${\chi }_{\nu }^{2}$ values for the models we tested are clustered around one, so it is fair to assume that a "correct" model is contained in the set. Figure 9 shows that, although the posteriors are not perfectly Gaussian, they have sharp, unimodal peaks, and so the uncertainty derived from marginalization is informative.

The model weights are defined by Equation 2 (adapted from Gibson (2014)). Here, W_q is the weight assigned to model q, $P({M}_{q}| D)$ is the the likelihood of model q given the data, and $P(D| {M}_{q})$ is the likelihood of the data given model q, which is equivalent to the Bayesian evidence of model q, E_q. The denominator is a normalization term, summed over N models. We assume a conservative prior that each model is equally likely (P(M_i) = 1 for all i).

$$\begin{eqnarray}\begin{array}{rcl}{W}_{q} & = & \displaystyle \frac{P({M}_{q}| D)}{{\displaystyle \sum }_{i=1}^{N}P({M}_{i}| D)}=\displaystyle \frac{P(D| {M}_{q})P({M}_{q})}{{\displaystyle \sum }_{i=1}^{N}P(D| {M}_{i})P({M}_{i})}\\ & = & \,\displaystyle \frac{{E}_{q}P({M}_{q})}{{\displaystyle \sum }_{i=1}^{N}{E}_{i}P({M}_{i})}.\end{array}\end{eqnarray} \tag{ 2 }$$

The marginalized log metallicity with 1σ uncertainties is calculated from Equations (15) and (16) from Wakeford et al. (2016). The result is ${\mathrm{log}}_{10}Z/{Z}_{\odot }={2.29}_{-0.36}^{+0.24}$ (${194}_{-109}^{+144}$ × Z_⊙). As expected, the highest-weighted models are the offset models. Bayesian model averaging demonstrates that, in PLATON's chemical equilibrium framework, a supersolar metallicity is the most likely result even after accounting for uncertainty in model selection.

The marginalized metallicity is useful as a reference, but it is valuable to give the metallicity distribution for each specific model assumption. Marginalization is most appropriate when the specific model parameters are unimportant; however, we are interested in the impact that modeling assumptions have on the atmospheric parameters. We emphasize that, even for apparently "data-defined" methods, many assumptions have to be made, and those should be explicitly stated for an appropriate interpretation.

In this section, we compare the results from the preferred PLATON and AURA models. These include the fiducial model with partial clouds and Gaussian instrumental offsets for PLATON (Section 6.2.3) and the Gaussian instrumental offset model for AURA (Model 8; Section 8.2).

The similarities reveal the most robust conclusions of our analysis, since they are retrieved despite the many different assumptions that went into each method. Notably, both retrievals robustly find a metal-rich atmosphere with metallicity (defined as O/H) inconsistent with the solar metallicity at >2σ. Both methods find a decisive (>4.8σ) water vapor detection, and at least a moderate detection (>2.7σ) of a non-haze gas absorption feature in the optical. Further, both PLATON and AURA retrievals are consistent with a mostly clear atmosphere, with neither finding strong evidence of haze or uniform, high-altitude gray clouds.

Though the atmospheric properties derived from PLATON and AURA are similar, there are noteworthy differences. AURA infers a cooler limb temperature at 100 mbar (${1320}_{-200}^{+270}\,{\rm{K}}$ compared to ${1710}_{-80}^{+100}\,{\rm{K}}$ for PLATON) as well as a lower metallicity of ${\mathrm{log}}_{10}Z/{Z}_{\odot }={1.46}_{-0.68}^{+0.53}$ compared to ${\mathrm{log}}_{10}Z/{Z}_{\odot }={2.33}_{-0.25}^{+0.23}$ for PLATON, a difference of 1.3σ. This translates to ${29}_{-23}^{+69}$ × Z_⊙ for AURA and ${214}_{-88}^{+149}$ × Z_⊙. The optical absorber also differs: AURA determines the best description of the STIS feature to be absorption from sodium and AlO, whereas PLATON prefers some combination of TiO and VO absorption. Finally, AURA makes no claim on the C/O ratio, as it is a free retrieval framework and no C-bearing species are detected or meaningfully constrained. On the other hand, the chemical equilibrium assumption allows PLATON to find a 3σ upper limit on the C/O ratio of C/O < 0.83.

To further contextualize the results, we added the functionality to retrieve the abundance profiles of relevant molecules in PLATON. We show abundance profiles for six spectroscopically relevant species from AURA, which are also included in PLATON—H₂O, CO, CO₂, Na, TiO, and VO—in Figure 17. We emphasize that enforcing chemical equilibrium narrows the abundance constraints, and we are not reporting these abundances. Instead, they should be interpreted as the expected abundance profiles under the conditions of stable chemical equilibrium for the reported temperature, metallicity, and C/O ratio. As an example, we find no observational constraint on CO, but its abundance is well-defined under chemical equilibrium for the temperatures and metallicities that we do observationally constrain via the water feature. Still, these profiles provide a useful baseline for comparison to free-chemistry retrieval abundances.

Zoom In Zoom Out Reset image size

The abundance profiles for the optical-wavelength absorbers reflect the disagreement on the primary gas absorber: AURA prefers Na, and so it retrieves ∼10× more Na than PLATON and significantly less TiO and VO. Note that the decreasing TiO and VO abundance with increasing pressure for PLATON is due to those molecules condensing out of the atmosphere. PLATON's inferred water abundance is typically a few times greater than AURA's, reflecting the difference in inferred metallicities. CO₂ and CO are unconstrained by AURA, while PLATON finds a high abundance of CO₂ is consistent with the Spitzer observations. This difference is expected, given that PLATON finds weak evidence of CO₂ while AURA found none. Interestingly, this may relate to AURA's only chemical constraint, which is that CO2 must be less abundant than CO and H₂O, due to the inferred temperatures (Section 5.2).

In total, AURA finds a cooler atmosphere with less oxygen but a large sodium enrichment to explain the optical absorption, while PLATON finds a hotter atmosphere with more oxygen and TiO/VO absorption in the optical.

The differences between a free-chemistry retrieval (AURA) and one constrained by chemical equilibrium (PLATON) are the natural result of the different assumptions made by each method. We therefore consider the PLATON and AURA retrievals to be two orthogonal analyses. We examine the impact of the differences by first explicitly listing the notable assumptions in each method, and then by providing the rationale for which assumptions are driving the differences.

The relevant methodological differences for PLATON as compared to AURA are: (1) the assumption of chemical equilibrium, (2) fixing the elemental ratios between all metals other than carbon to their solar values, (3) assuming an isothermal profile for the atmosphere, (4) not including opacity from AlO, and (5) not including the Allard et al. (2019) H₂-broadened Na line profile.

We find that the Allard et al. (2019) H₂-broadened Na line profile is the key driver in the differences between the retrievals, and flexible element abundances and chemical equilibrium also play roles. AURA is the more flexible retrieval, so we first describe its solution before addressing why PLATON differs.

AURA's lower-temperature solution is preferred for being able to explain the H₂O feature in the WFC3 spectrum, while also explaining the STIS data with H₂-broadened Na absorption and capturing the Spitzer data. AURA is able to provide a fit to the STIS data by independently increasing the Na abundance and by also invoking AlO at relatively low temperatures. At this lower temperature (T ∼ 1300 K), the amount of oxygen necessary for the water abundance and scale height to explain the observed water feature is about 29 × Z_⊙, with a mean molecular weight of about 2.7 AMU and a scale height of about 440 km.

Since PLATON has not yet incorporated the H₂-broadened Na line profile, the low-temperature solution is a relatively poor fit to the STIS data. Instead, TiO/VO are needed to explain the STIS absorption feature, and these are only abundant enough in chemical equilibrium (with fixed metal ratios) at around 1650 K. At this higher temperature, a higher mean molecular weight is required for the same scale height, which must be small enough to explain the molecular feature sizes as well as the dominance of TiO/VO absorption over Rayleigh scattering. The atmospheric metallicity necessary to achieve the higher mean molecular weight is the much higher ∼200 × Z_⊙. Therefore, the differences make sense in light of the stricter assumptions.

To provide more support to this idea, we compare our results with those of a third retrieval method, ATMO (Amundsen et al. 2014; Tremblin et al. 2015, 2016, 2017; Sing et al. 2016), which acts as a middle ground between PLATON and AURA. ATMO's spectral retrievals can further help us to gain insight into the effect of retrieval assumptions, as it includes the Allard et al. (2007) pressure-broadened sodium line but also has the added flexibility of performing a free-element equilibrium-chemistry retrieval. With this assumption for the chemistry, the elemental abundances for each model are freely fit and calculated in equilibrium on the fly. Four elements were selected to vary independently, as they are major species that are also likely to be sensitive to spectral features in the data, while the rest were varied by a trace metallicity parameter ([Z_trace/Z_⊙]). By separately varying the carbon, oxygen, sodium, and vanadium elemental abundances ([C/C_⊙], [O/O_⊙], [Na/Na_⊙], [V/V_⊙]), we allow for nonsolar compositions—but with chemical equilibrium imposed such that each model fit has a chemically plausible mix of molecules, given the retrieved temperatures, pressures, and underlying elemental abundances.

The resulting retrieved atmospheric parameters describe an atmosphere that is most consistent with the one described by AURA. ATMO prefers a temperature of ${1190}_{-120}^{+170}\,{\rm{K}}$ and a metallicity (as defined by the oxygen abundance) of ${\mathrm{log}}_{10}{\rm{O}}/{\rm{O}}\odot ={\mathrm{log}}_{10}Z/{Z}_{\odot }={1.53}_{-0.67}^{+0.55}$, in excellent agreement with AURA's values, and consistent with PLATON's metallicity to 1.3σ, though the retrieved temperatures differ significantly. Like AURA, ATMO finds an enhanced sodium abundance, though uncertainties are large (${\mathrm{log}}_{10}\mathrm{Na}/\mathrm{Na}\odot \,={1.40}_{-1.80}^{+0.75}$). This supports the idea that the inclusion of H₂-broadened sodium line profiles and the flexibility of nonsolar metal ratios—and not necessarily the equilibrium chemistry constraint—allow for the low-temperature, lower oxygen abundance solution found by AURA. The metallicities on all three retrievals indicate a metal-rich atmosphere and agree at the ∼1.3σ level.

Like PLATON (Section 7), ATMO also finds a subsolar C/O ratio (C/O = ${0.17}_{-0.16}^{+0.53}$ consistent with stellar (C/O = 0.19), though carbon is not well-constrained, so the uncertainties are large. The 3σ upper limit of 0.94 is in good agreement with PLATON's 0.83 upper limit. However, unlike PLATON or AURA, ATMO finds no evidence of optical absorbers beyond Na; instead, it prefers a haze and Na to explain the STIS optical data.

Figure 18 elucidates the differences in retrievals by showing the median retrieved fiducial model for PLATON (red), AURA (black), and ATMO (green) from 0.3 to 10 μm. The 1 and 2σ uncertainty contours are shown for PLATON and AURA, both of which are smoothed with a Gaussian filter with σ = 15 for clarity. The AURA predictions are only shown up to 5 μm—as a free-chemistry retrieval, AURA retrieving on the 0.3–5 μm data does not place meaningful constraints on multiple molecules with significant opacity in the 5–10 μm range. Therefore, a prediction is not warranted.

Zoom In Zoom Out Reset image size

While there are subtle differences, such as PLATON and ATMO's preference for CO₂ at 4.5 μm and sodium's prominence at 0.6 μm in the AURA and ATMO retrievals, the most obvious difference is below 0.5 μm, where ATMO prefers a haze instead of a metallic oxide feature. Though ATMO does not include AlO as an opacity source, this difference is likely due to different condensation schemes. PLATON uses GGchem's prescription (Woitke et al. 2018) such that species condense out when it is energetically favorable. AURA is a free-chemistry retrieval, so there are no restrictions on oxides being in the gaseous phase. ATMO, however, includes rainout chemistry (Goyal et al. 2019), such that if a species condenses at a higher pressure, that then depletes the element above that layer. It is plausible that, although PLATON's condensation scheme allows TiO/VO to be in the gas phase around 1700 K, ATMO's scheme does not, making the metallic oxide feature difficult to capture.

In total, we tentatively favor AURA's derived atmospheric parameters over PLATON's, for two main reasons. First, the inclusion of the most up-to-date sodium line profiles and AlO opacity impact the retrieval. Second, constraints from interior modeling (Section 9.2), though not necessarily decisive, are consistent with AURA and in tension with PLATON. Overall, this paints a picture of an atmosphere with a supersolar—but not necessarily superstellar—metallicity, sodium enrichment, possible disequilibrium metallic oxides (e.g., circulated from dayside, dredged up due to vertical mixing), and a planet with a well-mixed interior and a limb temperature lower than the equilibrium temperature.