ACCESS: Tentative Detection of H₂O in the Ground-based Optical Transmission Spectrum of the Low-density Hot Saturn HATS-5b

We observed HATS-5b over a period of 4 yr, from 2014 to 2018, with the IMACS instrument on the Magellan Baade 6.5 m telescope at Las Campanas Observatory as part of the ACCESS project ¹⁹ (Rackham et al. 2017). We observed seven transits over seven nights in total, but we found two of these transits to be poorly affected by weather after reducing them, and so we excluded them from the final analysis. The five transits analyzed in this work were obtained on 2014 October 29, 2015 December 17, 2017 November 30, 2017 December 20, and 2018 November 27, which we refer to in order as T1–T5, respectively. In addition, we were able to extract one extra transit light curve by making use of the second-order spectrum from the 2015 observing run, as the setup of the observation allowed for this light to be captured without contamination. We refer to this hereafter as T2.2.

For all observations, we used the f/2 camera in the multi-object spectrograph configuration. IMACS has eight chips, arranged in a 4 × 2 configuration, and in this f/2 mode the light is typically split across both vertically aligned chips. This camera has a very large field of view ($27\buildrel{\,\prime}\over{.} 4$), which, combined with the ability for custom slit masks, allows for simultaneous observations of multiple comparison stars in addition to our target star. For each observation of HATS-5, between five and seven comparison stars were also observed in order to track various systematics—including atmospheric conditions—during data collection. These stars are all of similar spectral type and are relatively close by on sky; we list them and their properties in Table 1. All seven comparison stars were used for T1 and T2, although the second order for comp2 fell off the chip and so could not be used for T2.2. For T3–T5, only comp1, comp2, comp4, comp8, and comp9 were observed, and all five of these comparison stars were used for T3 and T4. For T5, comp4 showed odd behavior in the white light curve and so was left out of the analysis. The custom slit masks were designed to be wide enough to minimize the potential for slit losses at typical seeing conditions at the observatory, while also being narrow enough to prevent contamination from nearby stars, or any more than necessary telluric signals. We use 10''-wide and 22''-long slits for T1, 10''-wide and 10''-long slits for T2 and T2.2, and 20''-wide and 40''-long slits for the remainder. The move to wider slits was made to see whether obtaining more sky background would allow for a more complete removal of background noise, but no difference is seen in the final product's result between the two setups in our data. The lengths are chosen to cover enough sky past the spectra to properly subtract the sky background. These slit widths and lengths are used for all of the stars in the given setup.

The first transit observation, T1, used the "Spectroscopic2" setting, which provides unfiltered spectra covering the wavelength range ∼4000–9500 Å, in combination with the grism with 300 lines mm⁻¹ and a blaze angle of 175 (Gri-300-17.5). As this was still during the preliminary stages of ACCESS, the setup was changed in following years to take into account blocking filters to prevent higher-order light contamination. However, during this analysis we looked for second-order contamination in the T1 data set and did not find any evidence for this to any significant degree, consistent with other early ACCESS studies (Rackham et al. 2017). The T2 setup used the WB5600-9200 filter, which allows light at the same throughput in the range 5600–9200 Å with lower throughput tails down to 5350 Å and up to 9350 Å, although we only take advantage of the shorter-wavelength extension. All wavelengths beyond this range are blocked. This was combined with the grism with 150 lines mm⁻¹ and a blaze angle of 188 (Gri-150−18.8). This combination caused there to be a significant separation between the first- and second-order light that was fully contained on the detector, which allowed for the use of T2.2. The final three transits, T3–T5, all used the same setup with the GG455 filter, which has a sharp cutoff at 4550 Å and covers out to approximately 10000 Å, and once again the Gri-300-17.5 disperser. Our exposure times varied from night to night, as they were optimized with respect to atmospheric conditions to maintain counts at or below half-saturation, which is what we have found works best for these precise observations (Bixel et al. 2019) and keeps us within the linearity regime of the CCD. We use 2 × 2 binning for all of the observations to decrease the read-out time. Our observations are summarized in Table 2.

Stellar inhomogeneities like starspots and faculae cause an observable impact on transmission spectroscopy, even if these active regions are not directly transited (e.g., Rackham et al. 2018, 2019). HATS-5 is a typical G-type star, and Zhou et al. (2014) note a lack of chromospheric activity and a slow rotation rate, which points to a star with little predicted stellar activity. Nevertheless, to check for the likelihood of stellar activity impacting our results, we inspected photometric data from the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017). ASAS-SN photometrically monitors the entire sky in the V band with 24 telescopes worldwide. There are 320 observations of HATS-5 taken with cameras at the Cerro Tololo International Observatory in Chile and the Haleakala Observatory in Hawaii from 2013 October 31 to 2018 November 27, which covers our entire observation period, with an average photometric uncertainty of 18,500 ppm per datapoint.

We detrended these data using Gaussian processes (GPs, implemented with george; Ambikasaran et al. 2015), and we find that the resulting data are flat to within uncertainties. Nevertheless, we looked for periodicity using a generalized Lomb−Scargle periodogram, which excels at detecting periodic signals in unevenly sampled data (Lomb 1976; Scargle 1982; VanderPlas 2018). We find two potentially significant period signals in our periodogram that predict a ∼5000 ppm amplitude of variability, but with significant error bars (at around the same 5000 ppm level when phase-folded and binned to 30 minutes) and corresponding unconvincing phase-folded results. Thus, we additionally looked for potential stellar activity signals in data from the Transiting Exoplanet Survey Satellite (TESS; e.g., Ricker et al. 2014), which observed HATS-5 in Sectors 5 and 32, both with 2-minute cadence. We are unable to reproduce any of the signals from the ASAS-SN data in the TESS data, and while we find two more potentially significant (4σ) peaks in the TESS periodogram, neither of them looks substantially periodic when phase-folded either. Any periodicity down to around ∼500 ppm should be easily identifiable in the TESS data, and so we conclude that this star is quiet at the precision of the available measurements and that the signals we were resolving were simply due to unmodeled systematics. We also note that the star was quiet during previous radial velocity (RV) observations (rms scatter less than 5 m s⁻¹ for both RV instruments; Zhou et al. 2014), which fits with this view of a quiet host star. Nonetheless, we still model for stellar contamination as mimicking features in our transmission spectrum (see Section 5). Additionally, we acknowledge that Evans et al. (2016) found a stellar companion to HATS-5, but that this companion is at a separation of 174 ± 006 from HATS-5, meaning that it would be outside even our widest slit (slit width is given in diameter, meaning that a companion would need to be within 10'' to be in the slit for our widest slit configuration, and a visual inspection along the slit length in an image taken before the disperser was put in place did not show any additional sources).

An important component of performing precise transit spectroscopy involves obtaining precise and accurate orbital and physical parameters of the system under analysis. In particular, the orbital inclination i and the semimajor axis to stellar radius ratio a/R_*—which in turn define the impact parameter, $b=(a/{R}_{* })\cos i$—are of particular importance, as they have been shown to give rise to spurious features in the transit spectrum if not properly constrained (see, e.g., Alexoudi et al. 2020). Having this picture in mind, we decided to update the parameters of the system using the data presented in the discovery paper by Zhou et al. (2014), combined with additional data recently obtained by the TESS mission in Sectors 5 and 32 and precise Gaia EDR3 (Gaia Collaboration et al. 2016, 2021) parallaxes.

Using the methods outlined in Brahm et al. (2018, 2019), we used all the available photometry for HATS-5, along with its spectra and Gaia parallaxes, in order to derive precise stellar properties. In particular, we obtain a precise estimate on the effective temperature of T_eff = 5315 ± 80 K, a metallicity of [Fe/H] = 0.04 ± 0.05, and a surface gravity of $\mathrm{log}g={4.530}_{-0.022}^{+0.022}$—all consistent with, but more precise than, those obtained by Zhou et al. (2014). The absolute stellar mass and radius we derive from those constraints are ${M}_{* }={0.868}_{-0.028}^{+0.031}\,{M}_{\odot }$ and R_* = 0.838 ±0.009 R_⊙, which give rise to a precise stellar density of ${\rho }_{* }={2080}_{-122}^{+129}$ kg m⁻³.

Given the stellar parameters outlined above, we then perform a full fit to the entire photometric and RV data set for HATS-5b. In particular, we fit the discovery HAT-South light curve and the eight GROND light curves (two nights, four broadband light curves in the $g^{\prime}$, $r^{\prime}$, $i^{\prime}$, and $z^{\prime}$ filters obtained on each night) presented in Zhou et al. (2014), the full TESS data set comprising two sectors worth of data, and the Subaru High Dispersion Spectrograph (HDS) and Magellan Baade Planet Finder Spectrograph (PFS) RVs (also from Zhou et al. 2014).

Our fit is performed with juliet (Espinoza et al. 2019a). For all the photometric data sets except for the discovery HAT-South light curve, we use a GP with a multiplication between an exponential kernel and a Matérn kernel to model instrumental systematics, which we detrend with respect to time. In particular, for the GROND data, we use wide priors on the hyperparameters of this kernel, defining a log-uniform prior on the timescales from 1 minute to 100 days for all bandpasses (which is a fairly wide prior given the 1.4-minute cadence of the GROND data), and a log-uniform prior on the amplitude of this GP from 1 to 1000 ppm. We also add jitter terms to capture any residual white-noise component with a log-uniform prior between 1 and 30,000 ppm. We note that this is a very large upper prior limit, but we intentionally leave the bounds this wide owing to the uncertainties of ground-based data. In addition, the amplitude of the GP does not necessarily correspond back to physical units. For the TESS data, we first constrain the GP hyperparameters using the out-of-transit data, and we use the posteriors on those fits to perform the joint fit, which only includes the in-transit data. We do not include a systematics model for the HAT-South light curve, as these are published with a detrending already applied on the data.

As for the RVs, we assume that these are well explained by a Keplerian orbit. We use a prior on the stellar density obtained in our spectroscopic analysis described above in our fits. However, following Sandford et al. (2019), we multiply the stellar density uncertainty obtained there by two in our prior, in order to account for empirically determined underestimations of the uncertainty in that work. We tried both circular and eccentric fits but found via the log-evidences that both models are indistinguishable by the data—we find the eccentricity of the system to be e < 0.05 with a probability of 99% given the data.

The priors and posteriors for the planetary and stellar parameters constrained by this joint fit are presented in Table 3. The full set of priors and posteriors for all parameters can be found in a dedicated GitHub repository. ²⁰

Table 3. Prior and Posterior Parameters of the Global Fit Performed to HATS-5b

Parameter	Prior	Posterior
Stellar and planetary parameters
P1 (days)	N(4.7, 0.12)	${4.76339146}_{-0.00000043}^{+0.00000046}$
t0,1 (BJD)	N(2459198.5, 0.12)	${2459198.51299}_{-0.00025}^{+0.00026}$
Rp,1/R⋆	U(0.0, 1.0)	${0.10620}_{-0.00051}^{+0.00054}$
${b}_{1}=(a/{R}_{\star })\cos (i)$	U(0.0, 1.0)	${0.143}_{-0.063}^{+0.061}$
K1 (m s−1)	U(0, 100)	${30.37}_{-1.11}^{+1.09}$
e1	fixed	0 (<0.05 at 99% credibility)
ω1	fixed	90
ρ⋆ (kg m−3)	N(2080, 2602)	${1997}_{-58}^{+44}$
TESS limb-darkening coefficients
q1,TESS a	U(0, 1)	${0.47}_{-0.11}^{+0.12}$
q2,TESS a	U(0, 1)	${0.24}_{-0.10}^{+0.12}$
Ground-based photometry limb-darkening coefficients
q1,HS	U(0, 1)	${0.14}_{-0.09}^{+0.13}$
q2,HS	U(0, 1)	${0.25}_{-0.16}^{+0.23}$
${q}_{1,{\rm{GROND}},g^{\prime} }$	U(0, 1)	${0.79}_{-0.10}^{+0.10}$
${q}_{2,{\rm{GROND}},g^{\prime} }$	U(0, 1)	${0.32}_{-0.06}^{+0.07}$
${q}_{1,{\rm{GROND}},r^{\prime} }$	U(0, 1)	${0.63}_{-0.11}^{+0.11}$
${q}_{2,{\rm{GROND}},r^{\prime} }$	U(0, 1)	${0.18}_{-0.06}^{+0.07}$
${q}_{1,{\rm{GROND}},i^{\prime} }$	U(0, 1)	${0.53}_{-0.10}^{+0.11}$
${q}_{2,{\rm{GROND}},i^{\prime} }$	U(0, 1)	${0.20}_{-0.06}^{+0.08}$
${q}_{1,{\rm{GROND}},z^{\prime} }$	U(0, 1)	${0.36}_{-0.07}^{+0.08}$
${q}_{2,{\rm{GROND}},z^{\prime} }$	U(0, 1)	${0.32}_{-0.09}^{+0.10}$

Notes. For the priors, N(μ, σ²) stands for a normal distribution with mean μ and variance σ², and U(a, b) stands for a uniform distribution between a and b, respectively.

^a These parameterize the quadratic limb-darkening law using the transformations in Kipping (2013).

Download table as:  ASCIITypeset image

This fit (Figure 1), benefiting from the inclusion of the TESS data and our stellar density prior, significantly improves the precision of the system parameters (Table 3). First, we see an improvement of a factor of 20 on the uncertainty in the planetary period of this exoplanet, reaching subsecond accuracy (with an error of 0.04 s). We have been careful in checking that this precision is indeed warranted in principle, as all the HAT-South data we have used are in BJD-TDB (Barycentric Julian Date−Barycentric Dynamical Time), as well as the newly included TESS data, which should pose a precision floor of 10⁻⁶ s from the time stamps alone (Eastman et al. 2010). In addition, the transit parameters that define the physical properties of the system are also greatly improved. We reach a planet-to-star radius ratio precision of 0.5%, implying a transit depth constraint based on this joint data set of 11,279 ± 115 ppm. Interestingly, our light-curve fits provide a constraint on the stellar density that is even better than the one we use as a prior by a factor of 4, highlighting the quality of the high-precision light curves used in our joint fit. This constraint, combined with the precisely constrained period, gives rise to a value of $a/{R}_{* }={13.381}_{-0.132}^{+0.098}$ and i = 89.39 ±0.27 deg. Both are fully consistent with the values presented in Zhou et al. (2014) but are more precise by a factor of 2.5 and 1.1, respectively.

Zoom In Zoom Out Reset image size

3.1.1. Stellar Metallicity Analysis

We reduced the raw Magellan/IMACS data using the data reduction pipeline consistently used in ACCESS publications (Jordan et al. 2013; Rackham et al. 2017; Bixel et al. 2019; Espinoza et al. 2019b; McGruder et al. 2020; Weaver et al. 2020; Kirk et al. 2021; Weaver et al. 2021; McGruder et al. 2022), described in detail in Espinoza (2017), and which will be briefly summarized here. This pipeline first applies bias calibration and identifies the slit positions for each of the targets in the detector. Next, the pipeline usually corrects for bad pixels and cosmic rays, but we have omitted that process here, as we observed that it was giving suboptimal results in our particular data set. We thus decided to handle these outliers via optimal extraction (see Section 3.2.1), and we skipped this step in the reduction of our data sets. After this step, the spectra for the target and comparison stars are extracted for each exposure, and then the pipeline applies a wavelength calibration with the use of HeNeAr lamps that were taken on each night of observations with a narrower slit (05) than what is used for our science frames in order to obtain higher spectral resolution. Each line profile is fit with a Lorentzian and then matched with identified lamp spectral features. Then, the profile positions are fit with a polynomial to determine the wavelength solution. This is done for the first science frame of the night, after which each following science frame is cross-correlated with the first frame's identified template to determine and correct for any wavelength offsets. It should be noted that we do not apply a flat-field correction to the data—as with past ACCESS studies, we found that it introduced additional noise into the data rather than making a marked improvement (e.g., Bixel et al. 2019). The result is a cleaned spectrum on a known wavelength grid for the target star and each of the comparison stars.

3.2.1. Optimal Extraction

3.2.2. Wavelength Bins

As with previous ACCESS data sets, some of the spectra of our targets fell on two different chips, which creates a wavelength gap in the spectra. For HATS-5b, these gaps are approximately 5990–6120 Å, no gap, 9230–9330 Å, 9230–9330 Å, and 9280–9380 Å for Transits 1–5, respectively. All but the gap for Transit 1 fell redder than the wavelength coverage for Transits 1 and 2, so they would have been ignored regardless. However, we are careful to make our wavelength bins for the transmission spectra so that they do not cover the 5990–6120 Å region. To determine the wavelength bins themselves, we use ExoCTK's Forward Modeling tool (Bourque et al. 2021) with the Fortney grid (Fortney et al. 2010) to predict the expected atmosphere of HATS-5b. ²¹ Using this prediction, we position bins around the most likely features—Na i and K i—and space the rest evenly around them in bins of width 150 Å. These bins are listed in Table 4. The usable second-order spectrum covers 5300–7200 Å, so the bluest 11 bins have contributions from all six transits. The spectrum actually extends all the way to bin 16, but the uncertainties increase rapidly after bin 11, and so we exclude those. The rest are composed of the other five transits, besides bin 12, which had a complex, moving, bad-pixel contamination inside of an absorption feature in T2 that we found to be uncorrectable, and thus it only has four transits contributing to the final value.

Table 4. Wavelength Bins, with Corresponding Transit Depths in ppm

Bin	Wav Start	Wav End	Transits	T1	T2	T2.2	T3	T4	T5	Weighted Mean
1	5365	5515	6	${10583}_{-248}^{+259}$	${10929}_{-502}^{+510}$	${11093}_{-1166}^{+1311}$	${14363}_{-1059}^{+1142}$	${12324}_{-803}^{+816}$	${9280}_{-1079}^{+1047}$	11,429 ± 366
2	5515	5665	6	${11194}_{-213}^{+221}$	${11939}_{-531}^{+588}$	${12301}_{-612}^{+608}$	${11292}_{-1130}^{+914}$	${10312}_{-776}^{+770}$	${10121}_{-928}^{+1130}$	11,190 ± 307
3	5665	5815	6	${11401}_{-477}^{+421}$	${11764}_{-232}^{+241}$	${10891}_{-366}^{+412}$	${12718}_{-1449}^{+1733}$	${11827}_{-753}^{+785}$	${10935}_{-954}^{+1070}$	11,585 ± 364
4	5815	5965	6	${11329}_{-230}^{+215}$	${12447}_{-455}^{+458}$	${11738}_{-340}^{+338}$	${10783}_{-892}^{+965}$	${13609}_{-817}^{+935}$	${11982}_{-1119}^{+1349}$	11,986 ± 311
5	6120	6260	6	${11631}_{-213}^{+224}$	${12453}_{-464}^{+446}$	${11075}_{-734}^{+811}$	${12178}_{-866}^{+939}$	${10771}_{-701}^{+726}$	${10873}_{-1173}^{+1386}$	11,500 ± 320
6	6260	6410	6	${10543}_{-196}^{+219}$	${11979}_{-293}^{+286}$	${11884}_{-565}^{+561}$	${11659}_{-884}^{+955}$	${11702}_{-663}^{+731}$	${11230}_{-1362}^{+1707}$	11,494 ± 336
7	6410	6560	6	${10755}_{-246}^{+225}$	${12534}_{-535}^{+669}$	${12015}_{-740}^{+705}$	${11140}_{-1005}^{+1392}$	${11646}_{-712}^{+714}$	${11982}_{-1279}^{+1465}$	11,684 ± 367
8	6560	6710	6	${10609}_{-209}^{+216}$	${12618}_{-588}^{+597}$	${11236}_{-573}^{+546}$	${12662}_{-893}^{+975}$	${12042}_{-600}^{+644}$	${10765}_{-1059}^{+1170}$	11,653 ± 295
9	6710	6860	6	${11388}_{-1073}^{+1100}$	${12109}_{-469}^{+393}$	${11055}_{-565}^{+577}$	${10198}_{-1389}^{+1293}$	${12005}_{-760}^{+805}$	${10411}_{-1147}^{+1414}$	11,202 ± 406
10	6860	7010	6	${11052}_{-173}^{+179}$	${12637}_{-620}^{+619}$	${11322}_{-1182}^{+1400}$	${11786}_{-974}^{+1055}$	${11886}_{-795}^{+825}$	${10827}_{-906}^{+1034}$	11,585 ± 357
11	7010	7160	6	${11571}_{-222}^{+212}$	${12269}_{-605}^{+519}$	${11495}_{-969}^{+964}$	${12169}_{-888}^{+966}$	${12365}_{-847}^{+867}$	${12679}_{-1488}^{+1613}$	12,083 ± 383
12	7160	7310	4	${11898}_{-226}^{+202}$	⋯	⋯	${11483}_{-804}^{+891}$	${13464}_{-1363}^{+1365}$	${12336}_{-1107}^{+1110}$	12,292 ± 492
13	7310	7460	5	${12105}_{-426}^{+345}$	${12387}_{-789}^{+716}$	⋯	${11869}_{-1087}^{+1098}$	${12064}_{-946}^{+991}$	${12098}_{-1297}^{+1351}$	12,106 ± 434
14	7460	7610	5	${12178}_{-332}^{+351}$	${10677}_{-1187}^{+1198}$	⋯	${12591}_{-1235}^{+1241}$	${12497}_{-844}^{+974}$	${11579}_{-1248}^{+1469}$	11,908 ± 478
15	7610	7760	5	${12134}_{-384}^{+361}$	${12426}_{-609}^{+598}$	⋯	${12659}_{-974}^{+993}$	${11588}_{-966}^{+1039}$	${11349}_{-868}^{+801}$	12,028 ± 358
16	7760	7910	5	${12771}_{-482}^{+434}$	${12256}_{-696}^{+649}$	⋯	${12948}_{-989}^{+1011}$	${10833}_{-675}^{+702}$	${12180}_{-1353}^{+1472}$	12,203 ± 407
17	7910	8060	5	${11960}_{-561}^{+495}$	${13485}_{-690}^{+591}$	⋯	${11939}_{-918}^{+941}$	${12589}_{-930}^{+1014}$	${10614}_{-1302}^{+1661}$	12,113 ± 434
18	8060	8210	5	${12512}_{-272}^{+267}$	${13108}_{-519}^{+536}$	⋯	${11119}_{-745}^{+663}$	${11467}_{-650}^{+695}$	${11816}_{-1039}^{+1100}$	12,008 ± 313
19	8210	8360	5	${12476}_{-308}^{+277}$	${12602}_{-477}^{+463}$	⋯	${11352}_{-468}^{+432}$	${12465}_{-991}^{+949}$	${13588}_{-931}^{+870}$	12,499 ± 298
20	8360	8510	5	${12270}_{-368}^{+318}$	${12170}_{-416}^{+393}$	⋯	${11741}_{-638}^{+680}$	${11398}_{-717}^{+776}$	${11943}_{-1571}^{+1543}$	11,901 ± 389
21	8510	8660	5	${12096}_{-386}^{+358}$	${11828}_{-198}^{+192}$	⋯	${12083}_{-995}^{+1051}$	${11977}_{-947}^{+975}$	${12579}_{-1223}^{+1239}$	12,113 ± 384
22	8660	8810	5	${12518}_{-560}^{+479}$	${12239}_{-579}^{+454}$	⋯	${11684}_{-672}^{+690}$	${10830}_{-1063}^{+1134}$	${12004}_{-1016}^{+1030}$	11,858 ± 362
23	8810	8960	5	${11783}_{-228}^{+219}$	${12908}_{-476}^{+470}$	⋯	${10781}_{-723}^{+798}$	${13641}_{-1344}^{+1409}$	${10568}_{-989}^{+1011}$	11,940 ± 382
24	8960	9110	5	${11813}_{-251}^{+261}$	${11972}_{-241}^{+254}$	⋯	${13195}_{-611}^{+692}$	${11430}_{-1015}^{+1091}$	${9382}_{-1804}^{+1938}$	11,565 ± 456

Note. Note that bin 5 only spans 140 Å rather than 150 Å like the rest of the bins owing to the spacing requirements around the large predicted features and the chip gap. Bins 1–11 have contributions from all six transits, bin 12 from T1 and T3–5, and bins 13–24 from T1–2 and T3–5.

Download table as:  ASCIITypeset image

Following Espinoza et al. (2019b), we model our target star white light curve in magnitude space as a combination of multiple elements:

$$\begin{eqnarray}&&{m}_{T}(t)=-2.51{\mathrm{log}}_{10}T(t)+{m}_{c}(t)+\alpha ,\end{eqnarray} \tag{ 1 }$$

where m_T (t) is the magnitude of the target star, T(t) is a transit model, m_c (t) is the magnitude of (a combination of) the comparison stars, and α is a systematics model.

For the transit model T(t), we use batman (Kreidberg 2015), which assumes a uniform spherical planet. For our white-light light-curve fits, we fit for the period (P), a/R_⋆, R_p /R_⋆, impact parameter (b), and time-of-transit center (t₀). The priors for these parameters are set so as to be consistent with the posteriors found in Section 3.1, except for the planet-to-star radius ratio and the time-of-transit center, whose uncertainties are inflated to allow for possible wavelength-dependent depth changes, as well as transit-timing variations (which we look for and do not find evidence for, see Appendix D). We also define the eccentricity and associated argument of periastron passage (ω), but these are held constant (to zero and 90°, respectively). For limb darkening, we use the uninformative sampling scheme proposed by Kipping (2013) for a quadratic limb-darkening law; however, following Espinoza & Jordan (2016), we use a square-root law as the limb-darkening law to use in our analysis. Our white-light priors are shown in Table 5.

Table 5. The Priors and Results for Our White-light Light-curve Fits

	P	a/R⋆	Rp /R⋆	b	t0 a	q1	q2
Prior	${4.7633915}_{+0.0000004}^{+0.0000004}$	${13.479}_{-0.073}^{+0.059}$	${0.10628}_{-0.1}^{+0.1}$ b	${0.081}_{-0.048}^{+0.056}$	${9198.51302}_{-0.2}^{+0.2}$ b	0.0–1.0	0.0–1.0
Type	normal	normal	truncated normal c	truncated normal c	normal	uniform	uniform
T1	${4.7633915}_{-0.0000003}^{+0.0000003}$	${13.478}_{-0.037}^{+0.033}$	${0.1071}_{-0.0006}^{+0.0006}$	${0.0953}_{-0.0326}^{+0.0278}$	${6959.720301}_{-0.000097}^{+0.000096}$	${0.482}_{-0.086}^{+0.062}$	${0.342}_{-0.089}^{+0.061}$
T2	${4.7633915}_{-0.0000004}^{+0.0000004}$	${13.447}_{-0.045}^{+0.043}$	${0.1122}_{-0.0015}^{+0.0015}$	${0.0679}_{-0.0349}^{+0.0354}$	${7383.661216}_{-0.000135}^{+0.000134}$	${0.687}_{-0.159}^{+0.181}$	${0.283}_{-0.157}^{+0.129}$
T2.2	${4.7633915}_{-0.0000004}^{+0.0000004}$	${13.487}_{-0.053}^{+0.054}$	${0.1072}_{-0.0096}^{+0.0098}$	${0.0867}_{-0.0407}^{+0.0401}$	${7383.661072}_{-0.000222}^{+0.000227}$	${0.605}_{-0.136}^{+0.183}$	${0.211}_{-0.124}^{+0.136}$
T3	${4.7633915}_{-0.0000004}^{+0.0000004}$	${13.484}_{-0.047}^{+0.045}$	${0.1094}_{-0.0023}^{+0.0025}$	${0.0842}_{-0.0378}^{+0.0390}$	${8088.642899}_{-0.000177}^{+0.000173}$	${0.763}_{-0.166}^{+0.142}$	${0.298}_{-0.134}^{+0.099}$
T4	${4.7633914}_{-0.0000004}^{+0.0000004}$	${13.465}_{-0.044}^{+0.042}$	${0.1090}_{-0.0024}^{+0.0024}$	${0.0830}_{-0.0382}^{+0.0354}$	${8107.698916}_{-0.000160}^{+0.000160}$	${0.590}_{-0.165}^{+0.220}$	${0.341}_{-0.196}^{+0.171}$
T5	${4.7633915}_{-0.0000004}^{+0.0000004}$	${13.455}_{-0.041}^{+0.040}$	${0.1081}_{-0.0035}^{+0.0039}$	${0.0719}_{-0.0329}^{+0.0339}$	${8450.660493}_{-0.000142}^{+0.000128}$	${0.565}_{-0.106}^{+0.167}$	${0.197}_{-0.115}^{+0.160}$
Combined	4.7633915 ± 4 × 10−7	13.469 ± 0.0474	0.1088 ± 0.005	0.0815 ± 0.0378	⋯	⋯	⋯

Notes. The combined final planetary fit parameters are also given, which are used for the wavelength-dependent light-curve fit.

^a As t₀ changes for each transit, we use the ephemerides given here to predict when the time of each transit center should in principle be, and we leave the prior broad to account for any small shifts. The values given are actually t₀ − 2,450,000 for brevity. ^b Note that the range for the priors for most variables are the errors as derived in our fit in Section 3.1, but the ranges are inflated on these parameters for the purposes of our fit. The actual errors are ±0.0004 and ±0.00025 for p and t₀, respectively. ^c For our truncated normal priors, the Gaussian distribution is truncated at ±1σ.

Download table as:  ASCIITypeset image

For m_c (t), we use principal component analysis (PCA) on our comparison stars to extract a set of uncorrelated signals that best explain their variations in time. Following the approach of Jordan et al. (2013), this allows us to extract an ordered set of signals that best explain their light curves (the so-called "principal components"), which act as linear regressors in our fit. In order to account for the fact that we do not know a priori the optimal number of principal components to use in our fit, we follow Espinoza et al. (2019b) and perform Bayesian model averaging: we perform fits using 1, 2, ..., N components, where N is the number of comparison stars used on each data set, and then weight the posterior distributions of each of those model fits based on their Bayesian evidences to get the final, weighted posterior distribution, which now accounts for this uncertainty.

Finally, the α term is a systematics model used to account for any effects not corrected for by the comparison stars. This is modeled with GPs, and we use six regressors: the variations of air mass, wavelength shift of our target spectrum, FWHM, sky flux, trace center position, and time. We again use george for our GP implementation, and we use a multidimensional squared-exponential kernel as the kernel of choice for our fits, which has given excellent results on previous ACCESS analyses (e.g., Weaver et al. 2020). The GP regressors were standardized before the fitting (i.e., they were mean-subtracted and divided by their respective standard deviations).

Our fits are performed via nested sampling using the pymultinest (Buchner et al. 2014) Python library, which makes use of the MultiNest algorithm (Feroz & Hobson 2008; Feroz et al. 2009, 2019). We fit each transit once, flag any flux value deviating by ≥5σ, and fit the light curve again not considering those flagged time stamps. The results of these fits are given in Table 5 and are presented in Figure 2. We obtain precisions of 174–333 ppm for our first-order white-light light-curve fits and a precision of 839 ppm for our second-order white-light light-curve fit obtained from transit T2.2. A corner plot showing the posteriors of our white-light fits are given in Appendix A. The photon noise levels for our first-order white light curves are 75, 223, 72, 75, and 73 ppm, meaning that our residual standard deviation corresponds to roughly 3×, 1.5×, 4×, 2.5×, and 2.5× the photon noise for T1–T5, respectively. The photon noise level for our second-order white light curve is 342 ppm, meaning that our residual standard deviation corresponds to roughly 2× the photon noise for T2.2. This is broadly consistent with the typical capabilities from the ground (e.g., Rackham et al. 2017; Bixel et al. 2019; Espinoza et al. 2019b).

Zoom In Zoom Out Reset image size

To fit the wavelength-dependent light curves, we perform the same fitting process as described in Section 4.1, except we hold all of the system parameters constant to the combined white-light light-curve fit value, which is given in Table 5. The exception to this is t₀, which varies for each transit, and so we use the white-light fit value for each transit individually for the corresponding wavelength-dependent fit. We only allow the radius ratio R_p /R_⋆ to vary (up to ±0.1 from the fit value), along with the limb-darkening and GP terms. We complete this step for each of the 24 bins for each transit. Just as for the white light curves, we do one fit and then perform a 5σ clip before refitting to obtain our final result. An example of the results of one transit's wavelength-dependent fits is shown in Figure 3. The rest are given in Appendix B, and all of the transit depths for each wavelength bin are given in Table 4. We note that the uncertainties on some of the second-order (T2.2) wavelength-dependent light curves are smaller than that for the corresponding white light curve (839 ppm), and we attribute this to the systematics-dominated nature of the data.

Zoom In Zoom Out Reset image size

4.2.1. Transmission Spectrum

After we complete our wavelength-dependent light-curve fits, we are able to combine these fits into a final transmission spectrum through a weighted average, where for each individual spectrum we calculate 1000 samples normally distributed within the errors, and then find the mean and standard deviation of the summation of each bin's total samples for our final bin value and error. The combined spectrum of all our transits is shown in Figure 4, and the combined transit depths used are given in Table 4. We note that the error bars on the second-order transmission spectrum are not notably larger than those on the first-order ones, which is understandable since we are still systematics limited rather than photon limited, and also that the second-order spectrum is overall consistent with our other spectra. If we take the second-order transmission spectrum out of our final combined spectrum, the overall shape of the weighted average spectrum remains the same. However, the addition of the T2.2 spectrum decreases our error on the bluest 11 points by as much as 18%. Additionally, we note that unlike other ACCESS studies in which we needed to offset each transmission spectrum in order to align them, the spectra of each night were aligned without the need for an offset. We did attempt to mean-subtract the data to align them, but we obtained the same combined transmission spectrum as without the offset. We note that this lack of offset is consistent with no significant stellar activity from the host star, as found in Section 2.2.

Zoom In Zoom Out Reset image size

To search for atomic and molecular signatures in our final transmission spectrum, we employ the atmospheric retrieval technique using the open-source library petitRADTRANS (Molliere et al. 2019). Because this library performs its sampling using Nested Sampling, in particular using the MultiNest (Feroz & Hobson 2008; Feroz et al. 2009, 2019) algorithm via pymultinest (Buchner et al. 2014), it calculates the marginal likelihood—also referred to as the Bayesian evidence—$Z={ \mathcal P }({ \mathcal D }| {{ \mathcal M }}_{i})$, i.e., the probability of the data ${ \mathcal D }$ given model ${{ \mathcal M }}_{i}$, for each atmospheric fit. Combined with the prior probability of each model, ${ \mathcal P }({M}_{i})$, this allows us to compute the probability of each model given the data ${ \mathcal P }({{ \mathcal M }}_{i}| { \mathcal D })\propto Z\times { \mathcal P }({M}_{i})$, which enables a metric for performing model comparison between competing models attempting to explain the data. In this work, we assume that all competing models are equally likely, and so we perform comparisons between models by directly comparing their Bayesian evidences. In particular, following Trotta (2008), we consider models with Δ ln Z ≤ 2 to be indistinguishable (which roughly corresponds to a difference of 2.5σ). Before looking into more complex atmospheric models, we test the simple fit of a flat line to our data, which is consistent with the null hypothesis of no atmospheric signals detected. We calculate the Bayesian evidence for this "Flat" model using the Exoretrievals package (Espinoza et al. 2019b).

As many of the spectral signatures that would allow for a determination of abundances through a free chemistry retrieval lie at redder wavelengths in the near-IR, we decide to first perform retrievals using a chemically consistent framework. The abundances of various atoms and molecules can be determined with respect to each other from just the metallicity, [Fe/H], and the carbon-to-oxygen ratio, C/O, for some given temperature−pressure profile if chemical equilibrium is assumed, which allows for a significant reduction in the number of free parameters in our fit. We implement this in petitRADTRANS using its subpackage poor_mans_nonequ_chem, which interpolates abundances of defined species from a precalculated chemical equilibrium grid as a function of pressure, temperature, [Fe/H], and C/O. This grid was calculated with easyCHEM, described in Molliere et al. (2017). We implement this into our retrieval to obtain our atmospheric mass fractions and corresponding mean molecular weights for each retrieval call. All of the abundances assume that the atmosphere is almost entirely composed of hydrogen and helium before these trace elements are taken into account—76.6% H₂ and 23.4% He. We include H₂–H₂ pair collisionally induced absorption (CIA) for all of the fits, which has been argued to be important for retrieving accurate atmospheric abundances (Welbanks & Madhusudhan 2019).

Besides [Fe/H] and C/O, we allow R_pl (reference radius), $\mathrm{log}{\kappa }_{\mathrm{IR}}$ (the ratio between the infrared and optical opacities), and the Guillot γ parameter to vary. These last two parameters are required because we use the temperature−pressure profile described in Guillot (2010) to define the structure of our atmosphere. We fix P₀, the pressure that corresponds to R_pl, the surface gravity $\mathrm{log}g$, the equilibrium temperature T_eq, the internal temperature T_int, and the stellar radius R_⋆. These are given in Table 7.

We first test the most basic model and include all of the elements and molecules that are included in the default abundance output: VO, CO, H₂O, C₂H₂, TiO, Na, K, HCN, CH₄, PH₃, CO₂, NH₃, H₂S, and SiO. We refer to this fit as "Consistent, Full, Clear" from here on. However, we quickly realized that the inclusion of Na and K seems to create large features not present in our data at the equilibrium parameters that fit the rest of the data. Thus, we leave them out for the remainder of our fits, assuming that these are somehow depleted in the gas form in the terminator region of HATS-5b (potential explanations for this are discussed in further detail in Section 5.3). With the remaining molecules, we obtain the atmospheric model we label "Disequilibrium, No Na, K, Clear," as some amount of disequilibrium must be acting for these species to not be seen (such as condensation or rainout). While we include all of these molecules for completeness and for potential future predictions, we note that some of the molecules do not seem to be contributing to the final predicted atmospheric spectrum, as their opacities are simply flat in this region. However, this additional complexity does not affect our corresponding $\mathrm{ln}Z$, as the fit values are still only [Fe/H] and C/O, which we confirmed via testing the exclusion of various flat molecules and seeing no change in the $\mathrm{ln}Z$.

Thus far, we have ignored the potential effect of clouds, so we additionally test a fit including $\mathrm{log}{P}_{\mathrm{cloud}}$ (the pressure at which an opaque cloud deck is added to the absorption opacity), referred to as "Disequilibrium, no Na, K, Cloudy." These retrievals that fit for [Fe/H] and C/O we refer to collectively as our chemically consistent fits.

We note that in our fits thus far, it seems that the major spectral contribution is due to H₂O. Thus, to test the detection of H₂O, we perform a series of free retrievals on our data using only the parameters that may be potentially important in our spectrum—the alkalis (Na and K), H₂O, and clouds. As above, we also include CIA. As suggested in Welbanks & Madhusudhan (2019), the inclusion of CIA has the potential to guard against biases related to determining abundances without including all the molecules in a fit. Besides the abundances, we leave R_pl free, as before, as well as the temperature. Rather than the Guillot profile from the consistent model, the free retrievals use an isothermal temperature profile, so only one parameter is required. We fix $\mathrm{log}g$ and R_⋆, just as in the chemically consistent case. These models are labeled as "Alkalis + H₂O + Cloud," "H₂O + Cloud," "Cloud," "Alkalis + H₂O," and "H₂O."

We also consider the possibility that any spectral features in our transmission spectrum could be generated by stellar contamination (Rackham et al. 2018, 2019), i.e., by heterogeneities on the star rather than by only the planetary atmosphere. We note, however, that this model seems unlikely a priori to explain the data given that Rackham et al. (2019) predict stellar contamination signatures to be relatively small in amplitude for planets orbiting Sun-like stars (typical G stellar type with approximately solar metallicity) like HATS-5. This prior is further strengthened in this case, as HATS-5 seems to be a relatively inactive star at least from spectroscopic indicators (Zhou et al. 2014) and photometric monitoring of the star (see Section 2.2), as well as a lack of activity indicators in our own analysis of high-resolution observations (PFS, FEROS; see Section 3.1.1). Nevertheless, we test the potential effects of stellar heterogeneities mimicking an atmospheric signal using the Exoretrievals package. The stellar contamination within the framework of that work is modeled following Rackham et al. (2018) and essentially determines the fraction of any stellar heterogeneities and the resulting effect on a transmission spectrum of choice. We allow for temperatures both less than and greater than the stellar photosphere, corresponding to potential spots and faculae, respectively. We model it with an underlying flat line model, to test the idea of this spectrum being purely mimicked by stellar activity. We refer to this model as "Flat + Stellar Contamination."

The results of all our retrievals are given in Table 6. Out of all the models tested, our best-fit (i.e., the one for which we attain the largest Bayesian evidence) chemically consistent model is the "Disequilibrium, No Na, K, Clear" case, and our best free retrieval is the "H₂O" case. This best-fit chemically consistent retrieval is plotted against the data in Figure 5, and the resulting retrieved parameters for both this and the "H₂O" case are given in Table 7. The "Flat," "Flat + Stellar Contamination," "Consistent, Full," and "Disequilibrium, No Na, K, Cloudy" models all have ${\rm{\Delta }}\mathrm{ln}Z\gtrsim 4$ when compared against our chemically consistent best-fit model. That is, given the data, those models are more than $\exp (\mathrm{ln}Z)\approx 50$ times less likely than our best-fit model (assuming that the models themselves are a priori equally likely). This lends moderately strong evidence to our best-fit model over those models. Notably, the higher log-evidence of the "Disequilibrium, No Na, K, Clear" retrieval against the "Disequilibrium, No Na, K, Cloudy" retrieval suggests a preference for a lack of clouds at the altitudes probed in transmission for HATS-5b, which is consistent with the theory for a "clear" window around 1000 K proposed by Gao et al. (2021). A plot showing all of the retrievals is given in Figure 13 in Appendix C.

Zoom In Zoom Out Reset image size

Table 7. The Parameters Used in Our Best-fit Atmospheric Models, Given with Their Description and Prior, as Well as Their Retrieved Fit Value

Parameter	Description	Prior (uniform)	Fit Value
Best-fit Chemically Consistent: Disequilibrium, No Na, K, Clear
Rpl	Reference radius	[0.8, 1.0] RJ	0.92 ± 0.01 RJ
C/O	Atmospheric carbon-to-oxygen ratio	[0.2, 1.7]	0.45 ± 0.27
[Fe/H]	Atmospheric metallicity relative to solar	[−2.0, 1.0]	−0.5 ± 1.0
$\mathrm{log}{\kappa }_{{IR}}$	Ratio between infrared and optical opacities	[−5,0]	−2.3 ± 1.6
γ	Guillot gamma parameter	[0.0, 0.8]	0.46 ± 0.23
P0	Pressure corresponding to Rpl	[0.01]	⋯
$\mathrm{log}g$	Planetary surface gravity	[2.85]	⋯
Teq	Planetary equilibrium temperature	[1025] K	⋯
Tint	Planetary internal temperature	[200] K	⋯
R⋆	Stellar radius	[0.838] R⊙	⋯
Best-fit Free Abundance: H2O
Rpl	Reference radius	[0.8, 1.0] RJ	0.92 ± 0.02 RJ
Temperature	Isothermal atmospheric temperature	[300, 1900] K	985 ± 313
H2O	Atmospheric water abundance	[−6.0, 6.0]	−3.0 ± 1.5
$\mathrm{log}g$	Planetary surface gravity	[2.85]	⋯
R⋆	Stellar radius	[0.838] R⊙	⋯

Note. The fixed parameters are also included with their fixed value.

Download table as:  ASCIITypeset image

Our best-fit chemically consistent model parameters are broadly consistent with expectations. The C/O ratio is subsolar, and less than 1 at ∼3σ, and the metallicity is approximately solar, although relatively unconstrained. The temperature−pressure profile, which is given solely by the Guillot profile parameters considering that we fix the temperature parameters, is also somewhat unconstrained but consistent with expectations.

From our best-fit free model ("H₂O"), which is almost identical in appearance to the chemically consistent best-fit model, we are able to determine the molecular contributions to the features in the consistent model. We derive the abundances of the fit H₂O to be $\mathrm{log}({X}_{{{\rm{H}}}_{2}{\rm{O}}})=-0.78\pm 0.56$. Thus, we constrain H₂O to a nonzero value as our best-fit model, even more likely than the best-fit chemically consistent model. The "H₂O" is a better fit than the "Flat" by ${\rm{\Delta }}\mathrm{ln}Z=5.5$, roughly corresponding to 3.6σ. We note that several of the free retrievals are essentially indistinguishable from each other given the evidence, but that all those with ${\rm{\Delta }}\mathrm{ln}Z\lt 3$ include the presence of H₂O. The retrieved temperature from this free retrieval is quite uncertain (985 ± 313 K) but consistent with the equilibrium temperature of the planet (1000 K). Corner plots of the retrieved parameters from both our best chemically consistent and free chemistry retrievals are given in Figures 14 and 15 in Appendix C, respectively.

From the results of our best fits, we report a tentative detection of H₂O in HATS-5b. H₂O is already known to be common in the atmospheres of hot giant planets and is a common detection in the near-IR (see, e.g., Sing et al. 2015), but not in the optical. Notably, Stevenson et al. (2016) found evidence for H₂O in the optical spectrum of HAT-P-26b with LDSS-3C at Magellan/Clay. However, we note that the spectrum is consistent enough with a flat line (${\rm{\Delta }}\mathrm{ln}Z\lesssim$ 4–6, depending on whether it is compared against our best-fit free or chemically consistent retrieval) that we consider this evidence for water alone to not be significant enough and thus in need of confirmation.

The lack of Na and K atomic features in the transmission spectrum of HATS-5b is, however, a puzzle. Our disfavoring of all of the retrievals including clouds suggests that they are not masked by clouds, but instead likely depleted in the atmosphere as suggested by the indistinguishability of the "Alkalis + H₂O" free abundance retrieval, but strong disfavoring of the "Consistent, Full, Clear" retrieval, which includes the alkalis at their predicted abundance from equilibrium chemistry. From our analysis of the stellar abundances (Section 3.1.1), the abundances of Na and K in the host star seem consistent with solar, so we can rule out an inherent depletion of these elements in the system. One possible pathway to explain it could be for these elements to be trapped in molecular form. It is predicted that Na and K condense at ∼900–1000 K at pressures of 10⁻² to 10⁻³ (which is where we probe with transmission spectroscopy) into Na₂S and KCl, at which point they are sequestered away out of the optical transmission spectrum (Lodders 1999; Madhusudhan et al. 2016; Line et al. 2017). However, we note that these condensation relationships are complicated and somewhat uncertain and that there may be another explanation for the lack of these features. Further chemical condensation and atmospheric structure models could dive deeper into these possibilities to explain the lack of Na and K features in the HATS-5b transmission spectrum. It is also possible that the resolution level probed by our observations is not high enough to identify sharp alkali features. Higher-resolution observations would be able to uncover whether this is the case.

5.3.1. Potential Follow-up Observations

As demonstrated in Figure 5, our chemically consistent best-fit model for the optical data predicts large water features in the near-IR, which would be observable with HST, as has been done often in the literature. For example, Wakeford et al. (2017) followed up on the previously described optical H₂O evidence in HAT-P-26 b with HST observations using STIS 750L, WFC3 G102, and WFC3 G141 (0.56–1.62 μm) and were able to obtain a precise H₂O abundance for the planet. While our visible water features, combined with the potential near-IR features, would confirm the result and allow for a more precise absolute water abundance determination, our simulations show that the result of HST observations would not constrain the general atmospheric parameters of the planet much more. We simulated one transit observation with WFC3/UVIS G280 (0.2–0.8 μm) and one with WFC3/IR G141 (1.1–1.7 μm) using PandExo (Batalha et al. 2017), based on our best-fit model result shown in Figure 5, and found that we would not constrain C/O significantly more and would only improve our estimate of [Fe/H] by a factor of 1.5.

However, observations with JWST would further extend our potential detections beyond H₂O to other molecules in the near-IR, which would significantly improve our understanding of this planet. In particular, observations with NIRSpec PRISM would cover near-IR molecular features other than water, such as CO₂ and CH₄, that would better constrain atmospheric parameters like C/O, as well as being able to confirm H₂O better than WFC3/IR G141 could. We simulated NIRSpec PRISM observations with PandExo (Batalha et al. 2017) and found that we could constrain these features with observations of just one transit. We show this in Figure 6. We improved the precision of our C/O and [Fe/H] determinations from our chemically consistent retrievals by 2× and 3×, respectively, with these additional data (to 0.49 ± 0.12 and −0.58 ± 0.33, respectively), and we increase our abundance retrieval of H₂O by a factor of 2 (in log space) with the corresponding free abundance retrievals.

Zoom In Zoom Out Reset image size