On the Hydrosphere Stability of TESS Targets: Applications to 700 d, 256 b, and 203 b

We implemented an adaptation of the Earth Similarity Index (ESI; Schulze-Makuch et al. 2011) to select the targets of interest. The ESI serves as a good basis for our purposes in this study by preferring Earth-like planets in terms of radius and insolation.

Instead of relying on one value for the ESI, however, we used the Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013) to estimate the likelihood that a target has an ESI of above 0.8. Input parameter distributions were obtained from the data validation files produced by the TESS Science Operations Pipeline (Jenkins et al. 2016) for all TOIs. The MCMC was sampled using 30 chains of 5000 steps each, resulting in 150,000 samples per parameter (Figure 1).

Zoom In Zoom Out Reset image size

The ESI of each TOI was calculated using the standard formulation except for using insolation in place of surface temperature. The inner threshold of insolation was set at the values determined by Leconte et al. (2015) to cause a runaway greenhouse effect, with the outer threshold being the CO₂ maximum greenhouse needed to sustain liquid water (Kopparapu et al. 2014). The three TOIs presented in this study were the top targets for this selection criterion. Additional TOIs have since been discovered, but none have had a higher likelihood of being Earth-like.

Initial VPlanet simulations were performed on each planet by sampling the system parameters available from the Exoplanet Follow-up Observing Program (ExoFOP 2019), assuming a Gaussian distribution, and performing the method described in Section 3 for the primary atmosphere test. These initial simulations indicated that TOI 203 b was almost certainly desiccated; 99.08% of the models that started with 10 Earth oceans of water ended with no water. The simulations of TOI 700 d and 256 b, however, only returned 30.35% and 15.18% desiccated models, respectively. From these results, we decided to move forward in applying the process detailed herein to TOIs 700 d and 256 b only.

The detection and characterization of transiting exoplanets begins with the host star, and inaccurate measurements of stellar properties have effects at every level of observation and modeling. In order to ensure that our targets are accurately modeled, we determine the host stars' radius, mean density, and effective temperature using standard methods in the analysis of M dwarfs (Berta-Thompson et al. 2015; Dittmann et al. 2017; Kostov et al. 2019; Ment et al. 2019) .

We acquire the measured K-band magnitude and parallax from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) and Gaia survey (Gaia Collaboration et al. 2016, 2018; Luri et al. 2018), respectively. The input parameters are assumed to be Gaussian, and the uncertainties on the output parameters are estimated with a simple Monte Carlo simulation using 200,000 samples from the input distributions. The necessary stellar parameters, as well as their estimated uncertainties, are included in Table 1.

The mass of the stellar host is determined using the mass–luminosity relation from Benedict et al. (2016). This mass is then used to estimate the stellar radius using a relationship for single stars (Boyajian et al. 2012), and the mean density is calculated. We estimate the bolometric luminosity as the mean of the corrected K and V magnitudes (Mann et al. 2015, erratum, and Pecaut & Mamajek 2013, respectively) and calculate the effective temperature using the Stefan–Boltzmann law. Finally, we confirm these parameters using the Mann et al. (2019) K-magnitude luminosity-distance–mass relationship.

As with other studies using TESS data (e.g., Gilbert et al. 2020), we decided to apply a reduction scheme that differs from the standard pipeline. These changes consist of creating a custom photometry aperture for the files, applying reduction and decorrelation techniques while masking transits, and moving the Gaussian process (GP) fitting/removal to the data fitting portion. This process necessitates that we begin with the 2 minute cadence TPFs that have undergone an initial cleaning for cosmic rays and flagged observations.

We used the Python program lightkurve (Lightkurve Collaboration et al. 2018) to download the publicly available TPFs and created a custom aperture for each using the threshold method. This method selects pixels that have a flux higher than 3 standard deviations above the median brightness and are contiguous with the center pixel. Each aperture is visually inspected to ensure that it includes only the target star and is devoid of other contamination.

To avoid overcorrecting the data during transit, each TPF is then masked for the duration of each planet's transit(s) with room for error. Expected transit times and lengths, as well as the variance for each, are taken from ExoFOP and included in Table 1.

Using the pixel and time series masks, we perform a pixel-level decorrelation on each TPF using the built-in method in lightkurve. Finally, we extract each light curve from the TPFs, remove NANs and outliers >7σ outside of the expected transits, and normalize the light curves. An example of the final, cleaned light curve is given in Figure 2.

Zoom In Zoom Out Reset image size

Light curves for each sector's observations are cleaned separately. The entire catalog of observations of TOI 256 was used, but only the first 11 sectors' observations were used for TOI 700 to be consistent with the data used in Gilbert et al. (2020) and Suissa et al. (2020).

Finally, the light curves are stitched together based on the average of the error in flux. The light curves for TOI 256 had similar errors and are stitched together before fitting. As TOI 700 has some differences in errors between sectors, the sectors were grouped together into four groups: sectors 1, 3, 4, and 5 form group 1; sectors 6 and 7 form group 2; sectors 8 and 9 form group 3; and sectors 10, 11, and 13 form group 4.

The reduced and extracted light curves are fit for transits using the exoplanet and its dependencies (Foreman-Mackey et al. 2021; Salvatier et al. 2016; The Theano Development Team, Al-Rfou, R., Alain, G., et al 2016; Kumar et al. 2019; Astropy Collaboration et al. 2018; Luger et al. 2019; Foreman-Mackey et al. 2017). The stellar mass and radius are parameterized via a normal distribution bounded above zero. Limb-darkening parameters are estimated in-model following Kipping (2013a) and were sampled uniformly. The impact parameters were sampled uniformly between zero (the center of the star) and one plus the planet-to-star radius ratio for each planet (the edge of transit overlap). We used a beta distribution as a prior on the eccentricities, as suggested by Kipping (2013b), bounded between zero and one. The periastron angle was sampled randomly. Transit depths and periods were parameterized using their natural logs. All other system parameters are modeled using a Gaussian distribution.

The list of parameters included in the fit, as well as the values, is available in Table 1. Parameters described above that are not included in the table had random initial values. Variance for some parameters in the prior distributions is intentionally left high so as to not unduly influence the fitting procedure.

To model noise, we fit a GP that is a mixture of two simple harmonic oscillators comprised of two hyperparameters, ${\rm{ln}}({\rho }_{\mathrm{GP}})$ and ${\rm{ln}}({\sigma }_{\mathrm{GP}})$, as well as a mean flux offset, μ_F , and the natural log of the variance, ${\rm{ln}}({\sigma }_{F})$. The initial value for μ_F is zero, and ${\rm{ln}}({\sigma }_{F})$ is calculated from the variance of the flux with transits masked. The initial value for ${\rm{ln}}({\sigma }_{\mathrm{GP}})$ is set to be ${\rm{ln}}({\sigma }_{F})$, and ${\rm{ln}}({\rho }_{\mathrm{GP}})$ is initially zero. The variance for the noise parameters is set at 10 for each value to allow adequate flexibility in the modeling process.

We first run the parameters through the version of the PyMC3 (Salvatier et al. 2016) optimization routine that was made for use with exoplanet. The initial transit light curves for each planet are plotted along with the GP model and residual (Figure 2) and visually inspected to confirm that the model transits line up with the observed light curve. In addition, the optimization solution is inspected for significant deviations from the expected parameters, as this can indicate that the model is incorrect and must be changed. In practice, changes are made if the median values differ from the expected by more than a standard deviation, the errors are abnormally large, or the model chains do not converge.

After the initial model is inspected, a mask is created for data outside 5σ of the residual mean. This secondary data clipping enables more accuracy in the GP modeling, resulting in a more accurate light-curve fit.

The resulting light curve is used in a second optimization of the parameters that yields a final model to be sampled using the No U-turn Sampler (Hoffman & Gelman 2011). We use four chains with 6000 tuning steps and 5000 draw steps each, yielding a trace of 20,000 steps total. For all traces, the Gelman–Rubin diagnostic (Gelman & Rubin 1992) was near 1.0, and the number of effective samples was over 1000 for all parameters. The sampled traces are saved as fits files and available upon request. Section 4 shows the deterministic parameters computed as part of this sampling, as well as their median and variance.

Before setting up the input files, we estimate the planet mass using the nonparametric method from Ning et al. (2018). We use the mrexo package for python (Kanodia et al. 2019) with the option to generate a complimentary MCMC trace for a posterior sample of the target planets' radii (Figure 3). A practical example can be seen in Figure 4 for TOI 256 b. This naive approximation can be replaced with the resultant probability distribution from simultaneously fitting transit and radial velocity measurements for more accurate simulations. We instead utilize this naive approach to highlight the potential insights gained from photometric transit data alone.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The atmesc model uses a simple, energy-limited atmospheric escape model (Watson et al. 1981) based on investigations of the solar wind by Parker (1964),

$$\begin{eqnarray}&&{F}_{\mathrm{EL}}=\displaystyle \frac{{\epsilon }_{\mathrm{XUV}}{{ \mathcal F }}_{\mathrm{XUV}}{R}_{p}}{4{{GM}}_{p}{K}_{\mathrm{tide}}{m}_{H}},\end{eqnarray} \tag{ 1 }$$

where ${{ \mathcal F }}_{\mathrm{XUV}}$ is the XUV energy flux, M_p is the mass of the planet, R_p is the planet radius, _XUV ≈ 0.1 is the initial XUV absorption efficiency, and K_tide is a tidal correction term of order unity (Erkaev et al. 2007). This model is an inferred trend in atmospheric escape as a function of planet size and irradiation based on statistical studies of all known exoplanets (Owen & Wu 2013; Lopez & Rice 2018). It is applied to these targets as an estimate in the absence of specific knowledge of their individual atmospheric escape processes, especially over such a large mass range. We allow _XUV to vary over the lifetime of the star as in Bolmont et al. (2016). This quantity is integrated over the planet's effective radius to incoming XUV energy,

$$\begin{eqnarray}&&{R}_{\mathrm{XUV}}=\displaystyle \frac{{R}_{p}^{2}}{H\mathrm{log}({p}_{\mathrm{XUV}}/{p}_{s})+{R}_{p}}\end{eqnarray} \tag{ 2 }$$

(Lehmer & Catling 2017), where H is the atmospheric scale height, p_XUV is the pressure at the effective XUV absorption level, and p_s is the pressure at the surface.

As fluctuation in the stellar luminosity can greatly affect the model output, we chose to use the accepted values of 0.0233 L_⊙ for TOI 700 and 0.00441 L_⊙ for TOI 256 instead of the trace values. The current stellar age is also taken to be the accepted values of 1.5 and 5 Gyr, respectively. The radius and mass (and thus the density) are fixed to the average values of the trace. The stellar module interpolates over the mass and time of the Baraffe et al. (2015) evolutionary tracks, and the XUV luminosity is tracked using the broken power-law model of Ribas et al. (2005):

$$\begin{eqnarray}\displaystyle \frac{{L}_{\mathrm{XUV}}}{{L}_{\mathrm{bol}}}=\left\{\begin{array}{ll}{f}_{\mathrm{sat}} & t\leqslant \tau \\ {f}_{\mathrm{sat}}{\left(\displaystyle \frac{t}{{t}_{\mathrm{sat}}}\right)}^{-{\beta }_{\mathrm{XUV}}} & t\gt \tau \end{array}\right..\end{eqnarray} \tag{ 3 }$$

As is standard in these simulations for M dwarf–type stars, we use f_sat = 1e − 3 for the saturation phase ratio, τ = 1 Gyr, and β_XUV = −1.23 for the power-law exponent, which is, in this case, set for Sun-like (smaller) stars (Ribas et al. 2005).

We assume that the planet has no hydrogen/helium envelope, as this is eroded in the model before water loss begins. In this case, the model assumes that atmospheric escape only takes place if the planet enters the runaway greenhouse threshold (Kopparapu et al. 2013). In this case, the surface temperature exceeds 647 K, fully evaporating surface oceans and leading to a water vapor–dominated upper atmosphere (e.g., Kasting 1988).

Efficient oxygen sinks are turned on in the model we use. This instantly removes photolytically produced oxygen from the atmosphere, preventing it from building up. This raises the threshold for the diffusion-limited flux of oxygen atoms through the background atmosphere of hydrogen,

$$\begin{eqnarray}&&{F}_{{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}=\displaystyle \frac{({m}_{{\rm{O}}}-{m}_{{\rm{H}}})(1-{X}_{{\rm{O}}}){b}_{{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}{gm}_{{\rm{H}}}}{{k}_{{\rm{B}}{\rm{o}}{\rm{l}}{\rm{t}}{\rm{z}}}{T}_{{\rm{f}}{\rm{l}}{\rm{o}}{\rm{w}}}},\end{eqnarray} \tag{ 4 }$$

which is the limit above which the model adjusts the energy-limited flux equation as follows:

$$\begin{eqnarray}&&{F}_{\mathrm{EL}}{\left(1+\displaystyle \frac{{X}_{{\rm{O}}}}{1-{X}_{{\rm{O}}}}\displaystyle \frac{{m}_{{\rm{O}}}}{{m}_{{\rm{H}}}}\displaystyle \frac{{m}_{c}-{m}_{{\rm{O}}}}{{m}_{c}-{m}_{{\rm{H}}}}\right)}^{-1},\end{eqnarray} \tag{ 5 }$$

where m_c is the crossover mass, the largest particle mass that can be dragged through the outflow:

$$\begin{eqnarray}&&{m}_{c}=\displaystyle \frac{1+{\textstyle \tfrac{{m}_{{\rm{O}}}^{2}}{{m}_{{\rm{H}}}^{2}}}{\textstyle \tfrac{{X}_{{\rm{O}}}}{1-{X}_{{\rm{O}}}}}}{1+{\textstyle \tfrac{{m}_{{\rm{O}}}}{{m}_{{\rm{H}}}}}{\textstyle \tfrac{{X}_{{\rm{O}}}}{1-{X}_{{\rm{O}}}}}}{m}_{{\rm{H}}}+\displaystyle \frac{{k}_{{\rm{B}}{\rm{o}}{\rm{l}}{\rm{t}}{\rm{z}}}{T}_{{\rm{f}}{\rm{l}}{\rm{o}}{\rm{w}}}{F}_{{\rm{H}}}^{{\rm{r}}{\rm{e}}{\rm{f}}}}{\left(1+{X}_{{\rm{O}}}\left({\textstyle \tfrac{{m}_{{\rm{O}}}}{{m}_{{\rm{H}}}}}-1\right)\right){b}_{{\rm{d}}{\rm{i}}{\rm{f}}{\rm{f}}}g}.\end{eqnarray} \tag{ 6 }$$

In the equations above, m_H is the mass of the hydrogen atoms, m_O is the mass of the oxygen atoms, k_Boltz is the Boltzmann constant, T_flow = 400 K is the temperature of the hydrodynamic flow (Hunten et al. 1987; Chassefière 1996), ${b}_{\mathrm{diff}}=4.8\times {10}^{17}{({T}_{\mathrm{flow}}/K)}^{0.75}\,{\mathrm{cm}}^{-1}\,{{\rm{s}}}^{-1}$ is the binary diffusion ratio for the two species (Zahnle & Kasting 1986), g is the planet's acceleration due to gravity, and X_O is the molar mixing ratio of oxygen at the base of the flow.

This model for determining hydrosphere degradation is a very approximate description and neglects to include many key factors in the loss of an atmosphere, such as the wavelength-dependent nature of heating the upper atmosphere and any nonthermal escape processes (e.g., flares). In addition, other photochemical and 3D hydrodynamic processes within the atmosphere can alter the rate of outflow. However, several statistical studies have shown that the escape rate scales with the stellar XUV flux and inversely with the gravitational potential energy of the gas (e.g., Lopez et al. 2012; Lammer et al. 2013; Lopez & Rice 2018; Owen & Wu 2013, 2017). The uncertainties regarding physics and the specifics of each planet are folded into the XUV escape efficiency _XUV ≈ 0.1, which has predicted the correct escape fluxes within a reasonable margin in these studies. The lack of specific information necessary to accurately model hydrodynamic escape necessitates the use of this kind of model for the purposes of our method, which aims to prioritize targets for follow-up observations to obtain this information.

In order to simulate the hydrosphere degradation during the active phase of the star, we assign the initial water on the simulated planet to be 10 Earth oceans as in Luger & Barnes (2014) for their most water-rich case. As we lack planetary formation and migration information of the target systems, their initial water inventory is difficult to ascertain. However, Earth is thought to have formed and accreted up to 70 oceans worth of water (Morbidelli et al. 2000; Raymond et al. 2006; Chassefière et al. 2012), making this a conservative estimate for initial water inventory.

The stability of a secondary hydrosphere is tested by setting the initial water inventory in Earth oceans to be proportional to the planet's mass in Earth masses. This tests the possibility that the planets can currently support a proportional amount of stable water for a long time and thus retain any primary or secondary hydrosphere for long enough that life could develop. Neither method is intended to be a definitive judgment of the planets' habitability.

The simulation results are given in Section 4. For statistics relating to TOI 700 d, the parameters are weighted according to the number of data points in their respective light curves. The output and aggregate files of the final parameters are available upon request.

We use four desiccation metrics to analyze the VPlanet outputs: desiccation percentage, remaining water inventory, percent of static models, and nonstatic desiccation rate. The necessary components of each metric are calculated during the data aggregation phase along with the results presented in Section 4.

3.2.1. Desiccation Percentage

3.2.2. Remaining Water Inventory

3.2.3. Static Cases

3.2.4. Nonstatic Desiccation Rate

In this section, we present the deterministic parameters computed by sampling the resultant trace of fitting TESS light curves for TOIs 700 and 256, as well as their median and variance. In addition, the naive mass prediction for each planet is given.

The TOIs 700 d and 256 b have posterior predicted periods of ${P}_{700{\rm{d}}}={37.406}_{-0.061}^{+0.067}$ and ${P}_{256{\rm{b}}}={24.737}_{-6.12e-5}^{+6.33e-5}$ days and posterior predicted radii of ${R}_{700{\rm{d}}}={1.123}_{-0.114}^{+0.131}$ and ${R}_{256{\rm{b}}}\,={1.650}_{-0.078}^{+0.087}$ R_⊕, respectively. The results of the fit are well in line with other studies using TESS data (Gilbert et al. 2020; Lillo-Box et al. 2020; ExoFOP 2019). Gilbert et al. (2020) and Lillo-Box et al. (2020) used other observations to model stellar variability from the All-Sky Automated Survey for Supernovae and HARPS + ESPRESSO, respectively, leading to small discrepancies between their studies and this work.

From radial velocity observations, the mass of TOI 256 b is known to be M_256b = 6.48 ± 0.46 M_⊕ (Lillo-Box et al. 2020), well within the standard error on our naive prediction, presented in Table 2, which also considers the potential for the target to be gaseous. The effect of the bimodal mass predictions is evident in Figures 5 and 6, where the predicted density has a degeneracy with the predicted water loss. The values given here and used in the simulations serve as examples for the use of this method when radial velocity measurements are not available.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We present the desiccation metrics for TOIs 700 d and 256 b in Table 3 and the correlation between water loss and density in Figures 5 and 6. In many cases, both targets appear to be able to maintain a stable hydrosphere, which is a strong predictor for their habitability; all cases with terrestrial densities retain some water.

For all planets in the target systems, modern simulations using proportional water inventories resulted in no water loss. This is due in majority to the low activity levels of the host stars, resulting in low XUV flux that did not trigger runaway greenhouse effects and thus water loss in the model. Thus, these results are not presented in Table 3. Unfortunately, no other planet besides the two targets shown in Table 3 had more than 1% of their simulations conclude with remaining water. In addition, of those that retained water, none had any static cases, meaning their hydrospheres would be predicted to degrade until fully desiccated. They are excluded from the table for this reason.

Desiccated models are included when calculating the remaining water inventory statistics, resulting in a peak in the distribution at zero. Because of this, the desiccation percent is a more useful metric for assessing the stability of a primary/formation hydrosphere. As expected, the planets that are safely in the habitable zone have the fewest desiccated cases and retain the most water.

4.2.1. TOI 700 d

4.2.2. TOI 256 b

As the model for water loss is an approximation, and many other factors that influence the hydrodynamic outflow of a planet's hydrosphere have been neglected, the results of this study are no final judgment of a target's habitability. In addition, further and significant loss of water has been shown to occur in a nitrogen-poor Earth atmosphere following the early runaway greenhouse period (Wordsworth & Pierrehumbert 2013). Figure 7 clearly shows this, as the models cease losing water the moment the early runaway greenhouse period ends.

Zoom In Zoom Out Reset image size

To combat these factors, we made an effort to induce a more conservative modeling procedure by turning on efficient oxygen sinks, beginning the model at t = 1 Myr, and removing any initial hydrogen envelope that might shield the water from XUV radiation. However, the results conveyed are likely accurate within a few factors and constitute a sufficiently robust prediction of water retention on the targets to inform follow-up observations that might shed further light on their specific scenarios.

4.3.1. Comparison with Initial Simulations

4.3.2. Mass Dependence for Rocky Bodies

4.3.3. Density Dependence for Intermediate Bodies

4.3.4. Comparison of TOI 700 d to Other Studies

4.3.5. Future Observation Recommendations