Characterizing a World Within the Hot-Neptune Desert: Transit Observations of LTT 9779 b with the Hubble Space Telescope/WFC3

We carry out the analysis of the transit data using Iraclis, our highly specialized software for processing WFC3 spatially scanned spectroscopic images (Tsiaras et al. 2016a, 2016c, 2018), which has been used in a number of studies (e.g., Anisman et al. 2020; Changeat et al. 2020; Edwards et al. 2021; Libby-Roberts et al. 2022; Brande et al. 2022; Garcia et al. 2022). The reduction process includes the following steps: zero-read subtraction, reference-pixel correction, nonlinearity correction, dark-current subtraction, gain conversion, sky-background subtraction, calibration, flat-field correction, and bad-pixel/cosmic-ray correction. Then we extract the white (1.088–1.680 μm) and the spectral light curves from the reduced images, taking into account the geometric distortions caused by the tilted detector of the WFC3 infrared channel. The pointing performance of HST was excellent across both visits and Figure 2 shows the shift in the x- and y-position of the spectrum.

Zoom In Zoom Out Reset image size

In the extraction of the flux from the calibrated images, Iraclis offers two techniques. The first extracts the flux from the image as it is at the end of the scan. The second splits the data into the up-the-ramp reads, subtracting each one sequentially from the next read and extracting the flux from each resulting image. The splitting extraction is useful for ensuring that there are no background stars with overlapping contributions to the target star's spectrum. However, differences between the flux extracted with each method can also occur due to persistence. While the level of persistence is dependent upon the accumulated charge, it is also correlated with the time the charge has spent on the detector as well as the time since the last read (Anderson et al. 2014).

In the context of HST WFC3, this means the persistence is dependent upon the brightness of the host star, the scanning rate, and the readout scheme employed. We show this effect for the second forward scan of the second orbit (the first image used in the analysis) in Figure 3. The differences are clear when comparing the G102 data sets to those taken with the G141 grism: the persistence effect is greater in the latter case, likely due to the higher flux levels as the scan rate and readout schemes were identical. As the extracted flux using the full scan and splitting methods can be different, and to minimize the impact of persistence, we use the splitting extraction in this work.

Zoom In Zoom Out Reset image size

We fit the light curves using our transit model package PyLightcurve (Tsiaras et al. 2016b) with the transit parameters from Table 2. The limb-darkening coefficients are calculated using ExoTETHyS (Morello et al. 2020) and based on the PHOENIX 2018 models from Allard et al. (2012). The stellar parameters used are also given in Table 2.

During our fitting of the white light curve, the planet-to-star radius ratio (R_p /R_s ) and the midtransit time (T₀) were the only free parameters, along with the parameters used to model the systematics (Tsiaras et al. 2016c). It is common for WFC3 exoplanet observations to be affected by two kinds of time-dependent systematics: long-term and short-term "ramps." The first affects each HST visit and is generally modeled as a linear behavior, while the second affects each HST orbit is modeled as having an exponential behavior (see, e.g., Kreidberg et al. 2014; Tsiaras et al. 2016c). The parametric model we use for the white light-curve systematics (R_w ) takes the form:

$$\begin{eqnarray}&&{R}_{w}(t)={n}_{w}^{\mathrm{scan}}(1-{r}_{a}(t-{T}_{0}))(1-{r}_{b1}{e}^{-{r}_{b2}(t-{t}_{o})}),\end{eqnarray} \tag{ 1 }$$

where t is time, ${n}_{w}^{\mathrm{scan}}$ is a normalization factor, T₀ is the midtransit time, t_o is the time when each HST orbit starts, r_a is the slope of a linear systematic trend along each HST visit, and (r_b1, r_b2) are the coefficients of an exponential systematic trend along each HST orbit. The normalization factor we use (${n}_{w}^{\mathrm{scan}}$) is changed to ${n}_{w}^{\mathrm{for}}$ for upward-scanning directions (forward scanning) and to ${n}_{w}^{\mathrm{rev}}$ for downward-scanning directions (reverse scanning). The reason for using separate normalization factors is the slightly different effective exposure time due to the known upstream/downstream effect (McCullough & MacKenty 2012).

We fit the white light curves using the formulae above, considering the uncertainties per pixel as propagated through the data reduction process. However, it is common in HST/WFC3 data to have additional scatter that cannot be captured by the ramp model. For this reason, we apply a scaling factor to each of the individual data point uncertainties to match their median to the standard deviation of the residuals and repeated the fitting (Tsiaras et al. 2018). The resulting fits for both the G102 and G141 observations are shown in Figure 4.

Zoom In Zoom Out Reset image size

Next, we extract spectral light curves using the Iraclis "default_high" binning for both grisms. We fit these spectral light curves with a transit model. In this case, the free parameters are the planet-to-star radius ratio and the parameters for the systematics model. The midtime was fixed to the best-fit value from the white light-curve fit while other orbital parameters were again fixed to those in Table 2. The systematics model, (R_λ ), includes the white light curve and a wavelength-dependent, visit-long slope (Tsiaras et al. 2016c), taking the form:

$$\begin{eqnarray}&&{R}_{\lambda }(t)={n}_{\lambda }^{\mathrm{scan}}(1-{\chi }_{\lambda }(t-{T}_{0}))\displaystyle \frac{{\mathrm{LC}}_{w}}{{M}_{w}},\end{eqnarray} \tag{ 2 }$$

where χ_λ is the slope of a wavelength-dependent linear systematic trend along each HST visit, LC_w is the white light curve, and M_w is the best-fit model for the white light curve. Again, the normalization factor we use (${n}_{\lambda }^{\mathrm{scan}}$) is changed to ${n}_{\lambda }^{\mathrm{for}}$ or ${n}_{\lambda }^{\mathrm{rev}}$ for the upward- or downward-scanning directions, respectively. Also, in the same way as for the white light curves, we perform an initial fit using the pipeline uncertainties and then refit while scaling up these uncertainties, such that their median matches the standard deviation of the residuals. The spectral light curves fits, for both visits, are given in Figure 5.

Zoom In Zoom Out Reset image size

We explore the nature of LTT 9779 b's lower atmosphere using Bayesian retrievals. These atmospheric retrievals are performed on the extracted transmission spectrum using the publicly available retrieval suite TauREx 3.1 (Al-Refaie et al. 2021, 2022). ³⁰ In our retrievals, we assume that LTT 9779 b possesses a primary atmosphere with a solar ratio of helium to hydrogen (He/H₂ = 0.17). The atmosphere of LTT 9779 b is simulated to range from 10⁻⁴ to 10⁶ Pa (10⁻⁹ to 10 bar) and sampled uniformly in log-space by 100 atmospheric layers. The retrieved radius is, therefore, the 10 bar radius. We include CIA from H₂–H₂ (Abel et al. 2011; Fletcher et al. 2018) and H₂–He (Abel et al. 2012) as well as Rayleigh scattering for all molecules. We model clouds as a uniform opaque deck, fitting only the cloud-top pressure (i.e., gray clouds). We then perform two main types of retrievals, those which fit each molecular species independently (i.e., free chemistry) and those which assume the molecules to be in chemical equilibrium.

For the free chemistry retrievals, we include the molecular opacities from the ExoMol (Tennyson et al. 2016), HITRAN (Gordon et al. 2016), and HITEMP (Rothman & Gordon 2014) databases. Based on the expected chemical species of such a hot planet, we include the opacities of H₂O (Polyansky et al. 2018), CO (Li et al. 2015), CO₂ (Rothman et al. 2010), HCN (Barber et al. 2013), TiO (McKemmish et al. 2019), VO (McKemmish et al. 2016), FeH (Wende et al. 2010), and H⁻ (John 1988). To implement the last of these, we use the methodology described in Edwards et al. (2020b), fixing the neutral hydrogen volume mixing ratio and thus only fitting for the e- volume mixing ratio. We note that Himes & Harrington (2022) showed that Equation (3) from John (1988) does not replicate the table of coefficients in that same paper. Himes & Harrington (2022) recommend using Table 1 from John (1988), which is what is implemented here.

In each free chemistry case, all molecular abundances are allowed to vary from log(VMR) = −1 to log(VMR) = −15. Higher mixing is not expected in the hydrogen-dominated atmosphere of LTT 9779 b and, in any case, these would also necessitate accounting for self-broadening of the molecular lines (Anisman et al. 2022a, 2022b). We explore using both an isothermal temperature profile and an NPoint profile. For the former, we fit only the temperature. For the later, we fit the temperature at the "surface" (10 bar) and the temperature at the "top" of the atmosphere (10⁻⁹ bar) as well as two intermediary points. For these intermediary points, we also fit the pressure at which these points occur. Therefore, the isothermal profile adds only a single free parameter to the retrieval while the NPoint adds six. For the free chemistry retrievals, we explore both chemical profiles which were constant with altitude and more complex, varying with altitude, two-layer profiles (Changeat et al. 2019).

Additionally, we conduct equilibrium chemistry retrievals using the code GGchem (Woitke et al. 2018) via the TauREx plug-in system (Al-Refaie et al. 2022). As with the free chemistry retrievals, we include Rayleigh scattering and CIA as well as simple gray clouds. For these retrievals, the free chemical parameters are the atmospheric metallicity and C/O ratio. For the equilibrium chemistry retrievals, the metallicity is allowed to vary from 10^–4 to 10⁴ and the C/O ratio has bounds of 10^–3 and 2.

Recovering absolute transit depths is difficult, particularly given the strong systematics present within HST observations (e.g., Changeat et al. 2020; Guo et al. 2020). Therefore, one can induce offsets in the transit depth between data sets from different instruments ³¹ which, if left uncorrected, could bias the retrieved atmospheric composition. Hence, we also conduct retrievals where we allow the G141 data set to be shifted relative to the G102 data, with this offset being a free parameter in the retrieval (e.g., Luque et al. 2020; Yip et al. 2021). We allow the G141 data to be shifted due to the strong correlations between the transit depth and systematics model (see Figure 4). To allow for this possibility to be fully explored, the bounds for this offset are set to be extremely broad: ±10^–2 (10,000 ppm).

Stellar activity, in the form of photospheric heterogeneities (i.e., spots and faculae), can contaminate the transmission spectrum of an exoplanet (e.g., Ballerini et al. 2012; McCullough et al. 2014; Rackham et al. 2018; Barclay et al. 2021). The contamination introduced is strongly chromatic with the strongest effects seen at shorter, optical wavelengths. As such, the possibility of stellar contamination is particularly important to consider for the G102 observations. Hence, we conduct retrievals using the Active Stellar Retrieval Algorithm (ASteRA) plug-in for TauREx, outlined in Thompson et al. (2023), to explore potential contamination in the case of LTT 9779 b. In our setup, we allow for possible contamination effects due to both unocculted spots and unocculted faculae. ASteRA requires four additional fitting parameters, which are the spot and faculae temperatures, T_Spot and T_Fac, respectively, and their respective filling factors, F_Spot and F_Fac. We again use intentionally unrestrictive priors, with F_Spot and F_Fac bounds of 0 and 0.99, T_Spot bounds of 4000 and 5440 K and T_Fac bounds of 5450–6000 K.

Finally, we explore the parameter space using the nested sampling algorithm Multinest (Feroz et al. 2009; Buchner et al. 2014) with 1000 live points and an evidence tolerance of 0.5. We then utilize the Bayesian evidence derived by Multinest as a means of selecting the preferred atmospheric model and for judging the significance of molecular detections.

Close-in gaseous planets are expected to be undergoing atmospheric escape due to the high levels of stellar irradiation. Indeed, hydrodynamic escape has been proposed as the cause of the dearth of hot Neptunes, with planets either being completely stripped or having a high-enough starting mass to retain the majority of their primordial envelope (e.g., Lundkvist et al. 2016; Mazeh et al. 2016).

Observational studies have been undertaken in an attempt to detect atmospheric escape and constrain the rate of mass loss. Many of these have focused on using the Lyα line, but despite these successes, attenuation by neutral hydrogen in the interstellar medium removes the majority of the Lyα line profile and so other tracers of mass loss have also been proposed. One such alternative line is that of metastable helium which occurs at 1.083 μm (Oklopčić & Hirata 2018). Ground-based observations have had notable success in detecting escape via this observable (e.g., Allart et al. 2018; Nortmann et al. 2018; Zhang et al. 2023), but the first detection of helium came via observations with HST. The G102 grism of WFC3 was used to observed WASP-107 b (Anderson et al. 2017), with an excess absorption being detected around 1.083 μm (Spake et al. 2018). These data suggested a mass-loss rate of 10¹⁰–3 × 10¹¹ g s⁻¹. Similarly, HST WFC3 G102 data of HAT-P-11 b found evidence for a mass-loss rate of 10⁹–10¹¹ g s⁻¹ (Mansfield et al. 2018) although high-resolution ground-based observations suggested a lower mass-loss rate (Allart et al. 2018). Observations of the helium triplet have also been undertaken for planets on the edge of the hot-Neptune desert, with the resulting mass-loss rates suggesting that the upper edge (i.e., higher-mass planets) of the desert is impervious to atmospheric escape (Vissapragada et al. 2022).

As LTT 9779 b lies within the rarely populated hot-Neptune desert it is expected that atmospheric loss should be ongoing, despite the old age of the system. Therefore, as well as using the HST WFC3 G102 data to constrain the atmospheric composition via a low-resolution spectrum, we also extract a higher-resolution spectrum to search for evidence of mass loss via the helium line. Using the excellent wavelength calibration offered by Iraclis, we extract overlapping bins with a bandwidth of 3 nm (30 Å) and a separation of 0.1 Å across the entire G102 range and perform fits to each of these light curves in the same way as discussed in Section 2.1. We then use the transit depths recovered in these overlapping bins to search for excess absorption around 1.083 μm.

Separately, we estimate the expected atmospheric escape rate using a photoevaporation model (A. Allan et al. 2023, submitted). In our model, the high-energy stellar radiation ionizes hydrogen and helium atoms, and the excess energy released after the ionization heats the planet's upper atmosphere, which more easily evaporates. We construct the incoming high-energy stellar spectrum of LTT 9779 using HST Space Telescope Imaging Spectrograph (STIS) near-ultraviolet (NUV)/optical observations (G230L and G430L) that were taken as part of the MUSCLES Extension for Atmospheric Transmission Spectroscopy program (GO-16166; PI: Kevin France; France et al. 2020). As LTT 9779 is a Sun-like star (Jenkins et al. 2020), the far-ultraviolet (FUV) portion of the spectrum (1170–2200 Å) is assumed to be the spectrum of the quiet Sun (Woods et al. 2009), scaled to the HST STIS G230L spectrum of LTT 9779; the adopted scaling factor was 9.0519 × 10⁻¹⁶. Meanwhile, the HST STIS G140M data are utilized to reconstruct Lyα using the methodology described in Youngblood et al. (2022). The Lyα flux is then used to compute the extreme-ultraviolet (EUV; 100–1170 Å) flux of LTT 9779 using the scaling relations from Linsky et al. (2014). Finally, as X-ray observations with Chandra led to a nondetection, we also use a scaled solar X-ray spectrum (5–100 Å).

Our atmospheric escape model requires as input the incident flux in four energy bins, X-ray (0.517–1.24 nm), hard EUV (10–36 nm), soft EUV (36–92 nm) and FUV+NUV (91.2–320 nm). By integrating the spectrum above, the luminosities of each of these bins are 3.35 × 10²⁶, 3.02 × 10²⁷, 2.30 × 10²⁷, and 2.22 × 10³¹ erg s⁻¹, respectively. Our model then solves for the hydrodynamics equations (conservation of mass, momentum, and energy) and ionization balance of hydrogen (Allan & Vidotto 2019). In our most recent model version (A. Allan et al. 2023, submitted), we also self-consistently compute the helium population (neutral helium in the singlet and triplet states, as well as ionized helium), and their corresponding energetics. The helium abundance is an input parameter of our model and we run two models, the first assuming an abundance by number of 2% and a second at 10%.

Maintaining exoplanet ephemerides is crucial for ensuring that atmospheric studies of exoplanets can be undertaken. Therefore, we utilize the midtimes from our HST light-curve fits, as well as literature observations, to refine the period of LTT 9779 b. Data from TESS Sector 2 led to the discovery of LTT 9779 b and it has since reobserved the host star in Sector 29. ³² These TESS data were analyzed by Kokori et al. (2023) as part of the ExoClock project (Kokori et al. 2021, 2022), which aims to refine the ephemerides of planets which will be studied by the Ariel mission (Tinetti et al. 2018, 2021; Edwards et al. 2019). In our work, we took the TESS midtimes from this previous study and incorporated new unpublished TESS eclipses. As with Kokori et al. (2023), we utilize the 2 minute cadence Presearch Data Conditioning (PDC) light curves (Smith et al. 2012; Stumpe et al. 2012, 2014). Due to the poor S/Ns of the TESS eclipses, we fit them in groups of four eclipses to achieve eclipse midtime uncertainties that were of a similar precision to the transit measurements. Furthermore, before the end of operations, Spitzer observed eclipses and phase curves of LTT 9779 b (Crossfield et al. 2020; Dragomir et al. 2020). We combine the midtimes from these studies with our data to extend the number of epochs during which LTT 9779 b has been observed. A full list of the midtimes utilized here are given in the Appendix in Tables 4 and 5.

To constrain the orbital period further, we utilize the radial velocity data taken of LTT 9779 by CORALIE and HARPS, which were presented in the detection paper of the planet (Jenkins et al. 2020). To model the radial velocity profile, we use RadVel (Fulton et al. 2018). The radial velocity data and the transit midtimes are fitted under one likelihood using Multinest (Feroz et al. 2009; Buchner et al. 2014) to explore the parameter space. For the linear ephemeris model, we assume a circular orbit (e = 0) following the conclusions of Jenkins et al. (2020).

Given the short period of LTT 9779 b, one might expect that the planet's orbit should be shrinking gradually over time due to the transfer of angular momentum from the planet to the host star, eventually leading to planetary engulfment. This may occur when the orbital period of a planet is shorter than the rotational period of its host star, as is the case for LTT 9779 b (Jenkins et al. 2020). In this scenario the star's tidal bulge lags behind the tidal bulge of the planet, generating a net torque that spins up the star and causes the planet to spiral inwards (Levrard et al. 2009; Matsumura et al. 2010). For a given arrangement, the rate of the decay of the orbital period is dependent on the efficiency of tidal energy dissipation within the star, parameterized by the modified stellar tidal quality factor ${Q}_{* }^{{\prime} }$ (Penev et al. 2018). Directly measuring the orbital decay rate via long-term transit timing measurements provides an estimate of ${Q}_{* }^{{\prime} }$, with a larger value being associated with less efficient dissipation and therefore a slower decay rate. While many studies have searched for evidence of orbital decay (e.g., Patra et al. 2020; Hagey et al. 2022), the only highly convincing detection so far has been for WASP-12 b (e.g., Maciejewski et al. 2016; Patra et al. 2017; Yee et al. 2019; Turner et al. 2020; Wong et al. 2022).

Therefore, in addition to fitting the timing data with a constant-period model, we fit an orbital decay model to the midtimes and radial velocity data of LTT 9779 b. For a planet on a Keplerian circular orbit with a constant orbital period P, we expect the central time of the transits t_tra and the occultations t_occ to increase linearly. Therefore, these can be described by:

$$\begin{eqnarray}&&{t}_{\mathrm{tra}}={t}_{0}+\mathrm{PE}\,\mathrm{and}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{occ}}={t}_{0}+\displaystyle \frac{P}{2}+\mathrm{PE},\end{eqnarray} \tag{ 4 }$$

where t₀ is the midtransit time of the reference epoch and E is the orbit number. In the case of orbital decay, the gradual shrinking of the orbital period gives rise to nonlinear drift in the observable transit and eclipse midtimes. The simplest way to represent this behavior is to include a quadratic term in the expected transit and eclipse center times. The model is given by:

$$\begin{eqnarray}&&{t}_{\mathrm{tra}}={t}_{0}+\mathrm{PE}+\displaystyle \frac{1}{2}\displaystyle \frac{{dP}}{{dE}}{E}^{2}\,\mathrm{an}{\rm{d}}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{occ}}={t}_{0}+\frac{P}{2}+\mathrm{PE}+\frac{1}{2}\frac{{dP}}{{dE}}{E}^{2},\end{eqnarray} \tag{ 6 }$$

where $\tfrac{{dP}}{{dE}}$ is the period derivative, which we assume to be constant. If negative, it signifies that the orbit of the planet is decaying. We again sample the parameter space using Multinest and compare the Bayesian evidence between each fit to determine the preferred model. The constant-period model has two free parameters: the reference epoch t₀ and the period P. The decay model has three parameters, adding the dP/dE term. As we fit the midtimes and radial velocity data simultaneously, each model has the additional radial velocity model parameters of the velocity semiamplitude (K) and systemic velocities of HARPS (V_H) and Coralie (V_C).

On the other hand, if the orbit of LTT 9779 b is slightly eccentric, one may expect the argument of periapse to precess over time due to a number of effects (see Section 4.2) which would induce sinusoidal variations in the transit and eclipse timing. Hence, we also attempt to fit an apsidal precession model to the timing data. In this case, the transit and eclipse times can be expressed as:

$$\begin{eqnarray}&&{t}_{\mathrm{tra}}={t}_{0}+{P}_{s}E-\displaystyle \frac{{{eP}}_{a}}{\pi }\cos \omega ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{occ}}={t}_{0}+\displaystyle \frac{{P}_{a}}{2}+{P}_{s}E+\displaystyle \frac{{{eP}}_{a}}{\pi }\cos \omega ,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\omega (E)={\omega }_{0}+\displaystyle \frac{d\omega }{{dE}}E,\,\mathrm{and}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{P}_{s}={P}_{a}\left(1-\displaystyle \frac{d\omega /{dE}}{2\pi }\right),\end{eqnarray} \tag{ 10 }$$

for argument of pericenter ω, phase ω₀, and precession rate d ω/dE. In these equations, P_s represents the planet's sidereal period, which is assumed to be fixed, while P_a is the "anomalistic" period. This latter period is also fixed, but accounts for the additional period signal due to precession. We note that this approach is an approximation and is only valid for e ≪ 0.1 (Ragozzine & Wolf 2009), but, given the very short period of LTT 9779 b and an upper limit of e < 0.058 at 95% confidence from Jenkins et al. (2020), this assumption is valid in this scenario. For the precession model there are a total of eight free parameters (t₀, P, ω₀, d ω/dE, e, K, V_H, and V_C).

We conduct a number of atmospheric retrievals using TauREx 3.1 using both free chemistry and chemical equilibrium models. Based on the Bayesian evidence, the preferred model is a free chemistry retrieval which assumed an isothermal atmosphere and for abundances of H₂O and CO₂ that changed with altitude. The model is preferred to that with constant chemical abundances by 2.41σ. The preference between this free model and the chemical equilibrium retrieval with a nonisothermal temperature profile is 2.40σ. The nonisothermal chemical equilibrium retrieval is preferred to its isothermal counterpart by 2.48σ. We show the best-fit spectra for the preferred free and chemical equilibrium models in Figure 6.

Zoom In Zoom Out Reset image size

Our preferred free chemistry retrieval favors the presence of H₂O, CO₂, and FeH while also placing upper limits on the presence of TiO and VO. The retrieved abundances for our preferred free and chemical equilibrium retrievals are compared in Figure 7. While many of the abundances agree across these two models, the CO₂ abundance is far higher in the free chemistry case, with the H₂O abundance in the lower half of the atmosphere also being noticeably higher. While our models preferred the presence of CO₂, WFC3 data only have limited sensitivity to carbon-bearing species and this can mean that the molecules retrieved are very dependent upon how the retrieval is set up (e.g., Changeat et al. 2020).

Zoom In Zoom Out Reset image size

Our preferred chemical equilibrium retrieval is that which uses a nonisothermal profile. It is likely that this allows for the model to have more flexibility in the dissociation of molecules with altitude compared to when an isothermal profile is used. However, when comparing the retrieved metallicity and C/O ratio we find that the temperature profile did not overly impact the results (see Figure 9). We note again that the HST WFC3 range does not give a strong sensitivity to carbon-bearing species and so the C/O ratio is not only poorly constrained but also not likely to be trustworthy (Rocchetto et al. 2016).

Spectral data from different instruments and epochs cannot be automatically assumed to be compatible. Recovering absolute transit depths is extremely difficult and the corner plots of the white light-curve fit for the HST WFC3 G141 observations (Figure 4) show a clear correlation between the transit depth and the systematics model. Therefore, despite the data sets seemingly being compatible when plotted together, we also perform retrievals which allow a global offset to be applied to the G141 data. For these, we utilize the setups from our preferred free chemistry and chemical equilibrium retrievals, with the only change being the addition of this offset parameter.

We find that, in the free chemistry case, the retrieval converges to a solution with a slight offset. However, given the 1σ uncertainties on this parameter, it is also consistent with there being no offset (−10${}_{-38}^{+33}$ ppm). Meanwhile, the chemical equilibrium retrieval prefers a more significant offset (−60 ± 26 ppm) and we show the preferred solution, including the offset to the HST WFC3 G141 data, in Figure 8. However, in both cases, the Bayesian evidence shows that the models without offsets were preferred by 3.96σ and 3.4σ for the free and chemical equilibrium cases, respectively.

Zoom In Zoom Out Reset image size

Due to the slope of the spectrum, we investigate the possibility that unocculted spots and faculae are the cause of the modulation seen. Unocculted faculae in particular are capable of reproducing the observed negative blueward slope. Our free chemistry retrievals in which this effect was accounted for lead to a fit which has a highly similar Bayesian evidence to that of the free chemistry model alone but requires four additional fitting parameters. The inclusion of the modeling of stellar heterogeneities removes any strong constraints on the atmosphere, with no molecules being conclusively detected. However, the retrieval converges to an extremely high faculae filling factor (65% ± 12% of the unocculted stellar surface). Given that the star is known to be slowly rotating and therefore inactive (Jenkins et al. 2020), and that there have been no indications of occulted faculae in any transit observations taken of the LTT 9779 b, this value is likely unrealistic.

Zoom In Zoom Out Reset image size

We provide an overview of our retrieval setups, as well as the associated evidence and reduced chi-square values, in Table 1.

Table 1. Summary of the Atmospheric Retrievals Conducted in This Work

Chemistry	TP Profile	Molecules Removed a	Two-layer Molecules	Offset	Spots/Fac	ln(E)	${\bar{\chi }}^{2}$
Free	Isothermal	None	H2O and CO2	False	False	299.10	1.57
Free	Isothermal	H2O	CO2	False	False	295.64	1.74
Free	Isothermal	CO2	H2O	False	False	297.51	1.62
Free	Isothermal	FeH	H2O and CO2	False	False	297.93	1.73
Free	Isothermal	None	None	False	False	297.77	1.46
Free	Isothermal	None	H2O and CO2	False	True	299.31	2.43
Free	Isothermal	None	H2O and CO2	True	False	293.76	1.68
Free	NPoint	None	H2O and CO2	False	False	298.32	2.15
Free	NPoint	None	H2O and CO2	True	False	294.98	2.27
Free	NPoint	None	None	False	False	297.61	1.87
Chem Eq	Isothermal	N/A	N/A	False	False	295.74	1.57
Chem Eq	NPoint	N/A	N/A	False	False	297.62	1.77
Chem Eq	NPoint	N/A	N/A	True	False	293.30	1.90

Note.

^a Compared to the standard setup of H₂O, CO₂, CO, HCN, TiO, VO, FeH, and H^–.

Download table as:  ASCIITypeset image

As well as the low-resolution spectrum used for our atmospheric retrievals, we also extract a higher-resolution spectrum using the G102 data with the aim of searching for the helium line at 1.083 μm that would be indicative of ongoing mass loss. However, despite sampling this spectral range with overlapping spectral bins which had a spacing of only 0.1 Å and a bandwidth of 3 nm (30 Å), we find no obvious evidence of excess absorption at, or around, this wavelength (see the top panel of Figure 10). The variations that are seen in this spectrum are likely due to correlated noise.

Zoom In Zoom Out Reset image size

We also model the predicted atmospheric evaporation using the methods of Allan & Vidotto (2019) and the reconstructed high-energy spectrum of LTT 9779 presented in Section 2.3. We find that assuming either a 2/98 or 10/90 helium to hydrogen number abundance, the evaporation rate is expected to be ∼10¹¹ g s⁻¹, reaching velocities of up to ∼50 km s⁻¹. However, even though the planet is anticipated to have a substantial evaporation rate driven by photoionization, the expected signature of the evaporation in the triplet lines is very small, particularly in the 2/98 case. The middle panels of Figure 10 show the anticipated transmission spectrum in the helium triplet, where we derive a peak transit depth of 0.29%. Due to the width of the spectral bins (3 nm/30 Å), the helium triplet would not appear as a sharp peak but as a plateau in the data: multiple bins would completely contain the excess absorption. In the bottom panel of Figure 10 we show the magnitude of this plateau when our 10/90 He/H model is convolved with these spectral bins. Even in the 10/90 case, the expected signature strength is far smaller than the uncertainties on the G102 spectral data (∼20 ppm versus ∼135 ppm). Therefore, we conclude the data are not sufficient to detect this tracer as even though we expect LTT 9779 b has a significant evaporation rate, its evaporation is barely seen in the helium triplet.

We fit the transit and eclipse midtimes with three models: linear ephemeris, orbital decay, and apsidal precession. We find no statistical evidence to prefer orbital decay over a linear ephemeris: the linear model has a Bayesian log evidence of ln(E) = −225.1 while the decay yielded ln(E) = −224.9. The apsidal precession model, however, is highly favored with a Bayesian evidence of ln(E) = −201.2. Examination of the data shows that the Spitzer transits have significant observed minus calculated (O − C) residuals compared to the uncertainties on the midtime (see Figure 11). These were taken from fits conducted by Crossfield et al. (2020) to the phase curves of LTT 9779 b. The eclipse midtimes from these phase curves (which were taken from the same data but independently fitted by Dragomir et al. 2020) do not show, in contrast, the same magnitude of deviation, thereby suggesting these points may be outliers rather than an actual detection of a significant deviation from the expected time.

Zoom In Zoom Out Reset image size

Due to the shape of the decay and precession models within an O − C plot (i.e., a parabola and a sine curve), we surmise that these Spitzer transit midtimes could be driving the model toward a solution that favors precession, with the large gaps in the transit coverage also enabling such a solution. Therefore, we also explore fits to the data without the Spitzer transit midtimes. In this case, the linear model is preferred with ln(E) = −181.0 versus ln(E) = −185.8 and ln(E) = −181.7 for the orbital decay and apsidal precession models, respectively. Though the evidence against apsidal precession versus a linear ephemeris is only marginal, the best-fit solution for the precession model is linear within the uncertainties (see Table 6). The inclusion or exclusion of the Spitzer transit midtimes has a minimal impact on the best-fit linear period but, given the potentially spurious nature of the midtimes, we report the recovered period without these points in Table 2. For completeness, the best-fit parameters for all fits are given in the Appendix.

Table 2. Stellar and Planetary Parameters Used or Derived in This Study

Parameter	Unit	Value	Source
R*	R⊙	0.949 ± 0.006	J20
M*	M⊙	${1.02}_{-0.03}^{+0.02}$	J20
T*	K	5480 ± 42	J20
log(g)	log10(cm s−2)	4.47 ± 0.11	J20
Fe/H	dex	0.25 ± 0.04	J20
Age	Gyr	2.0 ${}_{-0.9}^{+1.3}$
Prot	days	<45	J20
Rp	R⊕	4.72 ± 0.23	J20
Mp	M⊕	${29.32}_{-0.81}^{+0.78}$	J20
a/R*		${3.877}_{-0.091}^{+0.090}$	J20
i	degrees	76.39 ± 0.43	J20
e		<0.058	J20
Teq	K	1978 ± 19	J20
T0	BJDTDB	2,459,043.310602 ± 0.000090	TW
P	days	0.79206410 ± 0.00000014	TW

Note. J20: Jenkins et al. (2020); TW: this work.

Download table as:  ASCIITypeset image

We have presented transmission observations of the hot Neptune LTT 9779 b, extracting a spectrum from 0.8–1.6 μm. Our best-fit free chemistry model shows evidence for high abundances of H₂O and CO₂. To test the strength of each of these detections, we run additional retrievals. Without H₂O in the model, the best-fit model provides a worse fit to the data with the Bayesian evidence suggesting that the presence of H₂O is preferred by 3.12σ. In this model without H₂O, the retrieved CO₂ abundance is even higher than in our preferred model as shown in Figure 12. We also run a retrieval without CO₂, finding that the Bayesian evidence was reduced in this case: the inclusion of CO₂ was preferred by 2.01σ. The lower significance of this detection is likely due to the fact that WFC3 data only have limited sensitivity to carbon-bearing species. Furthermore, it has been shown that carbon-bearing species identified within a preferred model can be dependent upon how the retrieval is set up (e.g., Changeat et al. 2020). Interestingly, the removal of CO₂ from the model does not lead to the detection of CO or HCN. Instead, the H₂O abundance is affected with a slightly higher volume mixing ratio being preferred (see Figure 12). The other molecular species that was constrained in the free chemistry case was FeH and we found that our retrievals preferred its presence by 2.13σ. Hence, for all these species we find a relatively weak detection significance.

Zoom In Zoom Out Reset image size

Our retrievals with ASteRA prefer models with a high spot coverage. Such coverages are unrealistic, but could potentially be caused by assuming the wrong star temperature in our models. As the spots in the model are hotter than the rest of the stellar surface, the unrealistic coverage might indicate that our assumed stellar temperature is too low. We explore this possibility by independently determining the basic stellar parameters. We perform an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia DR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016) and Stassun et al. (2017, 2018). We pull the B_T V_T magnitudes from Tycho-2, the JHK_S magnitudes from the Two Micron All Sky Survey (2MASS), the W1–W4 magnitudes from the Wide-field Infrared Survey Explorer (WISE), the G_BP G_RP magnitudes from Gaia, as well as the NUV flux from the Galaxy Evolution Explorer (GALEX). Together, the available photometry spans the full stellar SED over the wavelength range 0.2–20 μm (see Figure 13).

Zoom In Zoom Out Reset image size

We perform a fit using PHOENIX stellar atmosphere models, with the free parameters being the effective temperature (T_eff) and metallicity ([Fe/H]); we set the extinction A_V ≡ 0, due to the close proximity of the system, and we assume a surface gravity $\mathrm{log}g\approx 4.5$ as appropriate for a main-sequence dwarf. Finally, we use the Gaia spectrum as a consistency check for the overall absolute flux calibration.

The resulting fit (Figure 13) has a best-fit T_eff = 5450 ± 75 K and [Fe/H] = 0.0 ± 0.3, with a reduced χ² of 1.4. Integrating the model SED gives the bolometric flux at Earth, F_bol = 3.503 ± 0.041 × 10⁻⁹ erg s⁻¹ cm⁻². Taking the F_bol together with the Gaia parallax directly gives the bolometric luminosity, L_bol = 0.7173 ± 0.0084 L_⊙, which with the T_eff gives the stellar radius, R_⋆ = 0.951 ± 0.018 R_⊙. In addition, we can estimate the stellar mass from the empirical eclipsing-binary-based relations of Torres et al. (2010), giving M_⋆ = 0.96 ± 0.06 M_⊙. Finally, the mass and radius together give the mean stellar density, ρ_⋆ = 1.58 ± 0.13 g cm⁻³.

These values are consistent with those from Jenkins et al. (2020), which we use for all the analyses here. In particular, the recovered temperature is not higher than assumed: it is, in fact, slightly lower (T_eff = 5450 ± 75 K versus 5480 ± 42 K). Therefore, it seems that this is not the cause of the high spot coverage preferred in our retrievals with ASteRA.

Given its high irradiation, it was anticipated that LTT 9779 b would be undergoing significant mass loss. Hence, it might be expected that this escape could be observed using tracers of atmospheric escape such as the helium triplet at 1.083 μm and thus that the nondetection within the HST WFC3 G102 is surprising. However, it is important to consider that the observability of the atmospheric escape via this particular tracer is not only dependent on the level of atmospheric escape. The population of helium in the triplet state also plays a critical role and this has been shown to depend on the nature of the stellar irradiance. The triplet state is more readily populated for K-dwarf hosts with their relatively low mid-UV flux (which can photoionize helium out of the triplet state) to EUV flux (which is both a significant contributor to driving the escape itself and populating the triplet via recombinations) compared to warmer G-type stars (Oklopčić 2019). Our work shows that LTT 9779 b indeed receives a high mid-UV flux, with our model finding that mid-UV photoionizations are the main depopulating source for triplet-state helium and hence largely responsible for the lack of a clear detection of helium at 1.083 μm. To cement this nondetection further, we apply the same fitting procedure to the HST WFC3 G102 data of WASP-107 b, which is shown in Figure 14 and highlight the strength of the signal seen in that case.

Zoom In Zoom Out Reset image size

Alternative tracers, such as the commonly used hydrogen Lyα line may allow the predicted strong escape of LTT 9779 b to be observed. Again utilizing the model of Allan & Vidotto (2019), we predict significant absorption in the uncontaminated Lyα line wings and show the resulting profile in Figure 15. However, we note that three-dimensional modeling considering the interaction with the stellar wind is required for a realistic Lyα prediction, given that high-velocity blueshifted absorption originates due to interaction with the stellar wind (Carolan et al. 2021). We leave such a detailed analysis for further work.

Zoom In Zoom Out Reset image size

In their study of helium escape for planets close to the Neptune desert, Vissapragada et al. (2022) derived a boundary below which a planet could not have lost more than 50% of its initial mass to photoevaporation. As LTT 9779 b is situated below this region, they noted that is was unlikely to be the photoevaporated core of a giant gaseous planet, with high-eccentricity orbital migration being the more likely explanation for the planet's current characteristics as the planet could have been delivered to its current location relatively late. Further evidence for this could be present in the planet's atmospheric composition (e.g., a high C/O ratio) or via the detection of a remaining eccentricity i the orbit of LTT 9779 b. However, the data analyzed here, and in previous studies (e.g., Crossfield et al. 2020; Dragomir et al. 2020), are not able to provide such evidence.

Hence, it is hard to conclude a great deal about the nature of the atmosphere of LTT 9779 b. We expect that far better constraints on the nature of LTT 9779 b's atmosphere from JWST NIRISS (Doyon et al. 2012) phase-curve observations, which have already been taken as part of Cycle 1. The transit portion of the JWST phase curve may also shed light on spectral variability when compared to the data presented here. Given the wider spectral coverage and higher S/N, the NIRISS transmission data should also be able to disentangle atmospheric features more easily from stellar contamination. An analysis of the global temperature and chemistry will test the implications made in this study, and those from the Spitzer data, and provide a more accurate interpretation of both. While we do not detect the signature of atmospheric escape from the helium triplet in this work, the JWST NIRISS observations may provide a detection. Additionally, observations with HST of the Lyα line may provide better constrains on the current escape rate pf LTT 9779 b's atmosphere and thus on the nature of the hot-Neptune desert.

Furthermore, eclipse observations of LTT 9779 b with CHEOPS have found evidence for variability (AO2-013; PI: James Jenkins), with additional observations being taken to confirm this (AO3-022; PI: James Jenkins). HST WFC3 UVIS is also being used to study the eclipse spectrum of LTT 9779 b (GO-16915; PI: Michael Radica; Radica et al. 2022) and these data may further constrain the eclipse-depth variability as well as helping to determine the albedo. Given the HST observations used here were taken roughly a year apart, the conclusions of these eclipse studies may indicate that, while the G102 and G141 data have similar transit depths in their overlapping spectral regions, dynamic changes to the planet may mean the spectra are not compatible.

While we do not find compelling evidence in the data to suggest that LTT 9779 b is undergoing orbital decay, we explore the magnitude of tidal dissipation of energy within the star necessary to yield an observable signal. Assuming the constant-phase lag model for tidal evolution suggested by Goldreich & Soter (1966), the expected decay rate is given by:

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dt}}=-\displaystyle \frac{27\pi }{2Q{{\prime} }_{* }}\left(\displaystyle \frac{{M}_{p}}{{M}_{* }}\right){\left(\displaystyle \frac{{R}_{* }}{a}\right)}^{5},\end{eqnarray} \tag{ 11 }$$

where M_p is the planet's mass, M_* is the star's mass, R_* is the star's radius, a is the semimajor axis of the planet's orbit, and ${Q}_{* }^{{\prime} }$ is the modified tidal quality factor. From our fit of the orbital decay model this yields a ${Q}_{* }^{{\prime} }$ of 1.5 × 10⁴ (without Spitzer transits) or 2.1 × 10³ (with Spitzer transits), both far outside the range of ${10}^{5}\lt {Q}_{* }^{{\prime} }\lt {10}^{8}$ expected given the current theory of energy dissipation within the convective envelopes of cool stars (Penev et al. 2018) and two orders of magnitude smaller than the value of 1.5 × 10⁵ derived for WASP-12 (Patra et al. 2017; Yee et al. 2019; Wong et al. 2022). Based on this analysis and the strong preference of the Bayes factor for a linear ephemeris, we conclude that the best-fit orbital decay rates from the data do not warrant further consideration.

Penev et al. (2018) derived an empirical model for ${Q}_{* }^{{\prime} }$ given the tidal forcing period of the star–planet system P_tide from an analysis of all known exoplanet systems with M_p > 0.1 M_Jup, P < 3.5 days, and T_eff,* < 6100 K, the parameter space in which the LTT 9779 system resides. They found:

$$\begin{eqnarray}&&{Q}_{* }^{{\prime} }=\max \left[\displaystyle \frac{{10}^{6}}{{({P}_{\mathrm{tide}}/\mathrm{days})}^{3.1}},{10}^{5}\right]\,\mathrm{wher}{\rm{e}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{tide}}=1/2\left(\displaystyle \frac{1}{{P}_{\mathrm{orb}}}-\displaystyle \frac{1}{{P}_{\mathrm{spin}}}\right).\end{eqnarray} \tag{ 13 }$$

Based on this result, we predict ${Q}_{* }^{{\prime} }$ for LTT 9779 by adopting a stellar rotation period of 45 days, as found in Jenkins et al. (2020). Taking this rotation period at face value requires the assumption that the orbit of LTT 9779 b is aligned with the stellar rotation axis and so this is, in effect, an upper limit. With this model and assumption, we predict that ${Q}_{* }^{{\prime} }=1.7\times {10}^{7}$ for LTT 9779, which corresponds to an orbital decay rate of −0.0085 ms yr⁻¹ for LTT 9779 b. Figure 16 shows that if orbital decay were occurring at this rate, the O − C expected within the current observational baseline would be far below the average uncertainty of the midtransit times. In the near future, more data will be available from CHEOPS, HST, JWST, and TESS observations. By the time the next TESS data for this planet are taken in Sector 69, LTT 9779 b will have orbited its host star around 2300 times since the first detection of a transit. If we assume the ${Q}_{* }^{{\prime} }$ predicted from the tidal forcing period, the O − C would be a mere −0.05 s by that time. Thus, despite the very short orbital period of LTT 9779 b it is unlikely that, if occurring, tidal orbital decay could be detected.

Zoom In Zoom Out Reset image size

Likewise, our fit to the timing data suggests that there is only evidence for orbital precession when the Spitzer data are included and that the best-fit apsidal precession model seems to be exploiting the large gaps in the coverage of the transit and eclipse data. Nevertheless, we explore the potential magnitude of precession we might expect for this world, comparing it to preferred value from our fit. We estimate a maximum theoretical apsidal precession rate by summing the contributions of general relativistic effects and tidal and rotational bulges of both the planet and star, using Equations (6), (10), and (12) from Ragozzine & Wolf (2009). The magnitude of apsidal precession raised by tidal and rotational bulges depends on the Love number k₂ of the perturbing body, thus we assume a value of k_2* = 0.03 for LTT 9779, typical of Sun-like main-sequence stars (Ragozzine & Wolf 2009), and k_2p = 0.3 for LTT 9779 b. We chose the latter as a generous value for a hot Neptune with a core mass fraction of ∼0.9 (Jenkins et al. 2020) given the results of the case study of the hot Neptune GJ 436 b by Kramm et al. (2011), noting that lower values would reduce the predicted rates of precession. Given that we additionally assume the upper limit of orbital eccentricity (e < 0.058) from Jenkins et al. (2020), we effectively calculate the maximum possible precession rate from the aforementioned effects. Our calculations yield:

$$\begin{eqnarray}\begin{array}{rcl}{\dot{\omega }}_{{\rm{tot}}} & = & {\dot{\omega }}_{{\rm{tide}},p}+{\dot{\omega }}_{{\rm{tide}},\ast }+{\dot{\omega }}_{{\rm{rot}},p}+{\dot{\omega }}_{{\rm{rot}},\ast }+{\dot{\omega }}_{{\rm{GR}}}\\ & = & (4.0\times {10}^{-5}+1.5\times {10}^{-7}+2.6\times {10}^{-6}\\ & & +\,3.6\times {10}^{-8}+1.1\times {10}^{-5})\,{\rm{rad}}/E\\ & \approx & 5.4\times {10}^{-5}\,{\rm{rad}}/E,\end{array}\end{eqnarray} \tag{ 14 }$$

where ${\dot{\omega }}_{\mathrm{tide}}$ and ${\dot{\omega }}_{\mathrm{rot}}$ refer to the contribution from the tidal and rotational bulges, respectively, and ${\dot{\omega }}_{\mathrm{GR}}$ is from general relativistic effects. The cumulative precession rate of 5.4 × 10⁻⁵ rad/E is two orders of magnitude below our best-fit model (d ω/dE = 3.1 × 10⁻³ rad/E). Given that we have assumed characteristics that would maximize the theoretical apisdal precession rate, we conclude that none of the above effects could be the source of the timing variation seen in Figure 11.

We also consider the possibility that apsidal precession may be induced by an additional body within the system. Before the confirmation of LTT 9779 b by Jenkins et al. (2020), there were suggestions of other planets within the LTT 9779 system (Pearson 2019), the best solutions being a 39.4 ± 9.5 M_⊕ planet on a 1.516 day orbit and a 73.4 ± 28.1 M_⊕ planet on a 1.65 day orbit. Using Equation (8) from Heyl & Gladman (2007), we calculate the expected apsidal precession rate raised by these proposed companions. However, we again find that the predicted precession rate(s) would be far below the best-fit value from our timing models (see Figure 17), and note that any companion inducing such a strong precession rate would be easily detectable in the radial velocity data. Hence, we find no reasonable way to explain the precession rate fit when including the Spitzer transits, concluding that it is not a viable orbital solution and that it is likely due to poor sampling of the transit and eclipse epochs since the discovery of LTT 9779 b. The combination of current and future transit and eclipse data from CHEOPS, HST WFC3 UVIS, JWST, and TESS should be capable of providing definitive evidence against the precession model, as highlighted in Figure 18.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size