The TESS-Keck Survey. II. An Ultra-short-period Rocky Planet and Its Siblings Transiting the Galactic Thick-disk Star TOI-561

We obtained a high signal-to-noise spectrum at R = 60,000 (Section 6) of TOI-561 to determine atmospheric parameters and detailed chemical abundances using the line list and forward modeling procedure of Brewer et al. (2016). The modeling uses Spectroscopy Made Easy (SME; Valenti & Piskunov 1996; Piskunov & Valenti 2017) in an iterative scheme that alternates between solving for global stellar properties and a detailed abundance pattern. We begin by estimating T_eff from B–V colors then fitting for T_eff, $\mathrm{log}g$, [M/H], Doppler line broadening, and the abundances of the α elements calcium, silicon, and titanium. All other elements are scaled solar values based on the overall metallicity given by [M/H] and the initial abundances are set to solar. The temperature of the resulting model is perturbed by ±100 K and used as input to refit the spectrum. The χ² weighted average of the global stellar parameters is then fixed and used as the input for the next step of simultaneously fitting for the abundances of 15 elements.

Simultaneous fitting of the elements is critical in obtaining precise abundances due to chemical processes in the stellar photosphere (e.g., Ting et al. 2018). The formation of molecules in cooler stars, even in very low numbers, can alter the atomic number densities and hence measured abundances using only isolated atomic lines.

The global parameters and abundance pattern obtained in the first iteration are then used as an initial guess for a second fitting following the same steps. Finally, the macroturbulence is set using a T_eff relation from Brewer et al. (2016), and we solve for the projected rotational velocity, $v\sin i$, with all other parameters fixed. The resulting gravities have been shown to be consistent with those from asteroseismology to within 0.05 dex and the abundance uncertainties are between 0.01 and 0.04 dex (Brewer et al. 2015). An empirical correction is applied to the abundances as a function of temperature (Brewer et al. 2016), which adds additional uncertainty to the absolute abundance, especially at temperatures between 5000 and 5500 K, and we adopt 0.05 dex uncertainty for most elements. Our analysis yielded a low stellar metallicity and high α abundance ([Fe/H] = −0.41 ± 0.05, [α/Fe] = + 0.23 ± 0.05, see Table 1).

Table 1. Host-star Characteristics

Basic Properties
Tycho ID	243-1528-1
TIC ID	377064495
Gaia DR2 ID	3850421005290172416
R.A.	09:52:44.44
Decl.	+06:12:57.00
Tycho VT Magnitude	10.25
TESS Magnitude	9.49
2MASS K Magnitude	8.39
Gaia DR2 Astrometry
Parallax, π (mas)	11.627 ± 0.067
Radial Velocity (km s–1)	79.54 ± 0.56
Proper Motion in RA (mas yr–1)	−108.432 ± 0.088
Proper Motion in DEC (mas yr–1)	−61.511 ± 0.094
High-resolution Spectroscopy
Effective Temperature, Teff (K)	5326 ± 64
Surface Gravity, $\mathrm{log}g$ (cm s−2)	4.52 ± 0.05
Projected rotation speed, $v\sin i$ (km s−1)	<2.0
log ${R}_{\mathrm{HK}}^{{\prime} }$ (dex)	−5.1
Iron Abundance, [Fe/H] (dex)	−0.41 ± 0.05
Carbon Abundance, [C/H] (dex)	−0.19 ± 0.05
Nitrogen Abundance, [N/H] (dex)	−0.51 ± 0.05
Oxygen Abundance, [O/H] (dex)	+0.09 ± 0.05
Sodium Abundance, [Na/H] (dex)	−0.39 ± 0.05
Magnesium Abundance, [Mg/H] (dex)	−0.20 ± 0.05
Aluminum Abundance, [Al/H] (dex)	−0.19 ± 0.05
Silicon Abundance, [Si/H] (dex)	−0.24 ± 0.05
Calcium Abundance, [Ca/H] (dex)	−0.27 ± 0.05
Titanium Abundance, [Ti/H] (dex)	−0.20 ± 0.05
Vanadium Abundance, [V/H] (dex)	−0.27 ± 0.05
Chromium Abundance, [Cr/H] (dex)	−0.43 ± 0.05
Manganese Abundance, [Mn/H] (dex)	−0.60 ± 0.05
Nickel Abundance, [Ni/H] (dex)	−0.37 ± 0.05
Yttrium Abundance, [Y/H] (dex)	−0.42 ± 0.05
Alpha Abundance, [α/Fe] (dex)	+0.23 ± 0.05
Distance Modulus and Isochrone Modeling
Stellar Luminosity, L⋆ (L⊙)	0.522 ± 0.017
Stellar Mass, M⋆ (M⊙)	0.805 ± 0.030
Stellar Radius, R⋆ (R⊙)	0.832 ± 0.019
Stellar Density, ρ⋆ (ρ⊙)	1.38 ± 0.11
Surface Gravity, $\mathrm{log}g$ (cgs)	4.500 ± 0.030
Age (Gyr)	10 ± 3
Transit Modeling
Limb Darkening (TESS band), q1	${0.2}_{-0.2}^{+0.2}$
Limb Darkening (TESS band), q2	${0.4}_{-0.2}^{+0.3}$

Note. The TESS magnitude is adopted from the TESS Input Catalog (Stassun et al. 2018), and the kinematics are taken from Gaia DR2 (Lindegren et al. 2018). Stellar parameters from isochrone modeling are formal uncertainties only and do not incorporate systematic errors from different model grids. The transit modeling is described in Section 4.

Download table as:  ASCIITypeset image

The effective temperature derived from the SME analysis (5326 ± 25 K) ⁵⁹ is in good agreement with alternative estimates using SpecMatch Synth (5249 ± 110 K; Petigura 2015), SpecMatch-Emp (5302 ± 110 K; Yee et al. 2017), and color–T_eff relations applied in the TESS Input Catalog (5440 ± 110 K; Stassun et al. 2018), and applying a J − K color–T_eff relation (5300 ± 110 K, Casagrande et al. 2010). We adopted the SME-derived solution, with an error bar calculated from the standard deviation of T_eff estimates from different methods: 5326 ± 64 K.

For each spectrum, we measure the Mt. Wilson S value, an indicator of the chromospheric magnetic activity. The Mt. Wilson S value is a measure of the strength of the emission cores in the Ca ii H and K lines relative to the nearby continuum flux. Our procedure for determining the S values is described in Isaacson & Fischer (2010). See Section 6, Table 2, for the full S-value time series. A Lomb-Scargle periodogram of the S values results in peaks near 100 days and 230 days, neither of which is near the expected rotation period or magnetic activity cycle of this old K dwarf (Section 6).

The stellar atmosphere and interior models are typically calculated using a solar-scaled α-element abundance mixture and thus assume [Fe/H] = [M/H]. To account for the nonsolar α abundances of TOI-561, we averaged the individual abundance measurements for [Mg/H], [Si/H], [Ca/H], and [Ti/H] to derive [α/Fe] = +0.23 ± 0.05, and then applied the calibration by Salaris et al. (1993) to convert the measured [Fe/H] value into an overall metal abundance, yielding [M/H] = −0.24 ± 0.10:

$$\begin{eqnarray}&&[{\rm{M}}/{\rm{H}}]=[\mathrm{Fe}/{\rm{H}}]+{\mathrm{log}}_{10}(0.694\times {10}^{[\alpha /\mathrm{Fe}]}+0.306).\end{eqnarray} \tag{ 1 }$$

We adopted a conservative uncertainty of 0.1 dex for [M/H] to account for potential systematics in the Salaris et al. (1993) calibration.

Next, we used T_eff, [M/H], $\mathrm{log}g$, Gaia DR2 parallax (adjusted for the 0.082 ± 0.033 mas zero-point offset for nearby stars reported by Stassun & Torres 2018), 2MASS K-band magnitude, a 3D dust map, and bolometric corrections to calculate a luminosity by solving the standard distance modulus, as implemented in the "direct mode" of isoclassify (Huber et al. 2017). We then combined the derived luminosity with T_eff and [M/H] to infer additional stellar parameters (mass, radius, density) using the "grid mode" of isoclassify, which performs a probabilistic inference of stellar parameters using a grid of MIST isochrones (Choi et al. 2016). The isochrone-derived $\mathrm{log}g$ (4.50 ± 0.03 dex) is in excellent agreement with spectroscopy (4.52 ± 0.05 dex), confirming that no additional iteration in the above steps is required for a self-consistent solution. The derived age of the isochrone fit is 10 ± 3 Gyr, consistent with the mean age of a galactic thick-disk star (see the following section).

The full set of stellar parameters is listed in Table 1. The results show that TOI-561 is an early-K dwarf with a radius of R_⋆ = 0.832 ± 0.019 R_⊙ and mass M_⋆ = 0.805 ± 0.030 M_⊙. We note that the quoted uncertainties are formal error bars and do not include potential systematic errors due to the use of different model grids (Tayar et al. 2021). For example, the stellar radius in Table 1 is 3% lower than predicted from an application of the Stefan–Boltzmann law using either the "direct mode" of isoclassify or SED fitting (Stassun et al. 2017), both of which yield 0.86 ± 0.02 R_⊙. However, this 3% (≈1 σ) difference does not significantly affect our main conclusions on the properties of the planets in the TOI-561 system, because the planet density errors are dominated by uncertainties in the planet masses (see Section 7).

We used the stellar evolution model fitting tool kiauhoku (Claytor et al. 2020) to estimate the rotation period of TOI-561. Using the stellar T_eff, [Fe/H], and [α/Fe] from Table 1 as inputs and assuming an age of 10 ± 3 Gyr, we found two different model-dependent estimates of the rotation period. Assuming the magnetic braking law described in van Saders & Pinsonneault (2013), we found P_rot = 38.5 ± 7.3 days, but assuming the stalled-braking law of van Saders et al. (2016), we found P_rot = 35.7 ± 3.4 days. These rotation periods are consistent with the upper limit of $v\sin i$ we determined spectroscopically. However, the estimated rotation periods differ significantly from the periodicity identified in the Mt. Wilson S-value activity indices (see Table 2), suggesting that the rotation period is not detected in the S-value time series. A rotation period of >30 days is likely too long to identify in the single sector of TESS photometry. We checked the archives of several ground-based photometric surveys, but TOI-561 saturates in ASAS-SN and Pan-STARRS, and it is too close to the equator to be included in WASP.

Early studies of star counts in the Milky Way revealed two distinct populations in the galactic disk which dominate at different scale heights, commonly denoted the "thin" and "thick" disk population (Gilmore & Reid 1983). Spectroscopic and photometric surveys have shown that these populations can be approximately separated based on kinematics and chemical abundances, with thick-disk stars being kinematically hotter (e.g., Fuhrmann 1998), older (Bensby et al. 2005), more metal poor, and enriched in α-process elements (e.g., Fuhrmann 1998).

The formation of the thick disk is still debated, with scenarios including external processes such as the accretion of stars from the disruption of a satellite galaxy (e.g., Abadi et al. 2003) and induced star formation from mergers with other galaxies (e.g., Brook et al. 2004), or a natural dynamical evolution of our galaxy including radial migration (Schönrich & Binney 2009a, 2009b). While the mere existence of a distinct thick disk is still in question (Bovy et al. 2012), spectroscopic and asteroseismic surveys have confirmed that chemically identified "thick-disk" stars belong to the old population of our galaxy, with typical ages of ∼11 Gyr (Silva Aguirre et al. 2018).

The detection of exoplanets around different galactic stellar populations can provide powerful insights into their formation and evolution (Adibekyan et al. 2012). For example, the discovery of five sub-Earth-size planets orbiting the thick-disk star Kepler-444 (Campante et al. 2015) demonstrated for the first time that terrestrial planet formation has occurred for at least ∼11 Gyr, and the discovery of a close M-dwarf binary companion demonstrated that this process can even proceed in a truncated protoplanetary disk (Dupuy et al. 2016). While TESS probes nearby stellar populations, it has significant potential to expand this sample. Indeed, Gan et al. (2020) recently presented the first TESS exoplanet orbiting a thick-disk star identified based on kinematics.

Figure 4 compares the chemical properties of TOI-561 with a sample of field stars in the TESS candidate target list (CTL) observed by the GALAH survey (De Silva et al. 2015; Sharma et al. 2018) and a sample of known exoplanet hosts from the Hypatia catalog (Hinkel et al. 2014). We calculated [α/Fe] for stars in the Hypatia catalog in the same manner as for TOI-561 and discarded stars with abundance uncertainties > 0.2 dex (calculated from the scatter between different methods). TOI-561 is consistent with the thick disk in terms of its chemical abundances, in agreement with the high proper motions measured by Gaia (Table 1) and the kinematic classification of TOI-561 by Carrillo et al. (2020). To independently confirm the kinematic classification, we used the UVW velocity vector of TOI-561 via the online velocity calculator of Rodriguez (2016), finding (U, V, W) = ( − 60.0, − 70.9, + 16.7) km s⁻¹. Using the probabilistic framework of Bensby et al. (2004) and Bensby et al. (2014), we find a thick-to-thin disk probability ratio of TD/D = 19, indicating strong evidence that this star is a member of the thick disk.

Zoom In Zoom Out Reset image size

TOI-561 is the first chemically and kinematically confirmed thick-disk exoplanetary system discovered by TESS, the fifth known thick-disk star known to host multiple planets, and the first thick-disk star known to host an ultra-period short planet. This further demonstrates that (1) small, rocky planets can form in metal-poor environments (consistent with Buchhave et al. 2012), (2) USPs are not tidally destroyed around old stars (consistent with Hamer & Schlaufman 2020), and (3) rocky planets have been forming for nearly the age of the universe.

We observed TOI-561 using the Las Cumbres Observatory Global Telescope (LCOGT) 1 m networks (Brown et al. 2013) in the Pan-STARRS z-short (zs) band. The telescopes are equipped with 4096 × 4096 SINISTRO cameras having an image scale of 0389 pixel⁻¹ resulting in a $26^{\prime} \times 26^{\prime}$ field of view. The images were calibrated using the standard LCOGT BANZAI pipeline (McCully et al. 2018), and the photometric data were extracted using the AstroImageJ (AIJ) software package (Collins et al. 2017). A full transit window of TOI-561 b was observed continuously for 205 minutes on 2019 April 19 UT from the LCOGT Siding Spring Observatory (SSO) node. TOI-561 c was observed continuously for 381 minutes on 2020 February 3 UT from the LCOGT McDonald Observatory node and again on 2020 March 17 UT from the LCOGT Cerro Tololo Inter-American Observatory (CTIO) node for 230 minutes and then later on the same epoch from the LCOGT SSO node for 269 minutes. TOI-561 d was observed continuously for 300 minutes on 2020 April 24 UT from the LCOGT SSO node.

The Next Generation Transit Survey (NGTS; Wheatley et al. 2018), located at ESO's Paranal Observatory, is a photometric facility dedicated to hunting exoplanets. NGTS consists of 12 independently operated 20 cm diameter robotic telescopes, each with an 8 deg² field of view and a plate scale of 5'' pixel⁻¹. The NGTS telescopes also benefit from subpixel guiding afforded by the DONUTS autoguiding algorithm (McCormac et al. 2013). By using multiple NGTS telescopes to simultaneously observe the same star, NGTS can achieve ultra-high-precision light curves of exoplanet transits (Bryant et al. 2020; Smith et al. 2020).

TOI-561 was observed on two nights using NGTS multi-telescope observations. On UT 2020 February 2, a full transit of TOI-561 c was observed using three NGTS telescopes. A predicted transit ingress of TOI-561 d was observed on the night of UT 2020 March 5 using four NGTS telescopes. A total of 5179 images were obtained during the first observation and 7791 during the second. Both sets of observations were performed using an exposure time of 10 s and the custom NGTS filter (520–890 nm). The airmass of the target was kept below 2, and the sky conditions were good for all the observations.

We reduced the NGTS images using the custom aperture photometry pipeline detailed in Bryant et al. (2020). This pipeline performs source extraction and photometry using the SEP library (Bertin & Arnouts 1996; Barbary 2016). The pipeline also uses GAIA DR2 (Gaia Collaboration et al. 2016, 2018) to automatically identify a selection of comparison stars that are similar to TOI-561 in terms of brightness, color, and CCD position.

We observed full transit windows of TOI-561 b continuously for 120 minutes on 2019 April 23 UT and 2020 May 24 UT simultaneously in the g, r, i, and z_s bands with the MuSCAT2 multicolor imager (Narita et al. 2019) installed at the 1.52 m Telescopio Carlos Sanchez in the Teide Observatory, Spain. The photometry was carried out using standard aperture photometry calibration and reduction steps with a dedicated MuSCAT2 photometry pipeline, as described in Parviainen et al. (2020).

We observed a full transit window of TOI-561 b continuously for 205 minutes on 2019 April 22 UT in the R_c band from the Perth Exoplanet Survey Telescope (PEST) near Perth, Australia. The 0.3 m telescope is equipped with a 1530 × 1020 SBIG ST-8XME camera with an image scale of 12 pixel⁻¹ resulting in a $31^{\prime} \times 21^{\prime}$ field of view. A custom pipeline based on C-Munipack ⁶¹ was used to calibrate the images and extract the differential photometry.

We observed a full transit window of TOI-561 b continuously for 206 minutes on 2019 April 23 UT in the R_c band from El Sauce Observatory in Coquimbo Province, Chile. The 0.36 m Evans telescope is equipped with a 1536 × 1024 SBIG STT-1603-3 camera with an image scale of 147 pixel⁻¹ resulting in an 188 × 125 field of view. The photometric data were extracted using AIJ.

The TOI-561 b SPOC pipeline transit depth is generally too shallow (290 ppm) for ground-based detection, so we checked all three stars within $2.5^{\prime}$ that are bright enough to have caused the SPOC detection (i.e., TESS magnitude <18.1) for a possible NEB that could be contaminating the SPOC photometric aperture. We estimate the expected NEB depth in each neighboring star by taking into account both the difference in magnitude relative to TOI-561 and the distance to TOI-561 (to estimate the fraction of the star's flux that would be contaminating the TESS aperture for TOI-561). If the rms of the 5 minute binned light curve of a neighboring star is more than a factor of 3 smaller than the expected NEB depth, we consider an NEB to be ruled out. We also visually inspect each star's light curve to ensure that there is no obvious eclipse-like signal, even though the rms to the estimated NEB depth threshold is met. Using a combination of the LCOGT, MuSCAT2, PEST, and El Sauce TOI-561 b follow-up observations, we rule out the possibility of a contaminating NEB at the SPOC pipeline ephemeris.

In the LCOGT observation of TOI-561 c on 2020 February 3 UT, we detected a 142 minute, early (0.3σ) ∼1100 ppm egress, relative to the nominal SPOC ephemeris, in a 97 radius aperture around the target star, which is not contaminated by any known Gaia DR2 stars. NGTS observed and detected the same transit on time (Figure 5). As a result, we revised the follow-up orbital period to 10.778325 days for further scheduling. The 2020 March 17 UT LCOGT CTIO and SSO observations then detected an on-time ingress and egress, respectively, at the revised ephemeris.

Zoom In Zoom Out Reset image size

The LCOGT observation of TOI-561 d on 2020 April 24 UT covered an egress ±150 minutes relative to the nominal SPOC pipeline ephemeris and provides ∼1σ coverage of the original SPOC ephemeris uncertainty on the epoch of observation, but is 2 days away from the transit midpoint predicted by our revised 16.29 day ephemeris. The LCOGT light curve does not show an obvious 923 ppm egress during the limited observation window. However, ingress- or egress-only coverage of transits with depths less than 1000 ppm from the LCOGT 1 m network can be difficult to interpret and reliably model due to potential trends in the data that may inject or mask a shallow ingress or egress. The NGTS observation on 2020 March 5 shows what might be the ingress of a transit of planet d. However, given the low S/N, this is not a high-confidence detection. Furthermore, if our revised 16.29 day ephemeris is correct, this observation was 2 days away from the transit midpoint.

Here we perform a joint analysis of TESS light curve and ground-based follow-up to refine the planetary parameters. We downloaded the TESS light curve from the Mikulski Archive for Space Telescopes (MAST; Figure 1). We isolated the transits of each planet with a window of three times the transit duration. We removed long-term stellar variability/instrumental effect by fitting a cubic spline with 1.5 day knot length to the light curve after removing the transits. We also downloaded the ground-based follow-up observations from the ExoFOP website ⁶² . We used the BATMAN (Kreidberg 2015) package for transit modeling, using the transit ephemerides reported by the TESS team as an initial guess for our model. Our model uses the mean stellar density ρ_⋆ as a global parameter (with a Gaussian prior as derived in Section 3 and Table 1) on all three planets. For each planet, we allowed the radius ratio R_p /R_⋆, the impact parameter b, the orbital period P, and the midtransit time T_c to vary freely. We assumed circular orbits for all three planets. The mean stellar density ρ_⋆ and the orbital period P together constrain the scaled semimajor axes a/R_⋆ of each planet. We adopted a quadratic limb-darkening law as parameterized by Kipping (2013). We allowed the coefficients q₁ and q₂ to vary in different photometric bands. We then performed a Monte Carlo Markov Chain analyses with the Python package emcee (Foreman-Mackey et al. 2013) to sample the posterior distribution of the various transit parameters. The results are summarized in Table 3, and Figure 6 shows the best-fit transit models.

Zoom In Zoom Out Reset image size

Table 3. Planet Parameters

Parameter	Median ± 1σ
Planet b
Orbital Period, Pb (days)	${0.446573}_{-0.000021}^{+0.000032}$
Midtransit Time, Tc (BJD)	2458517.4973 ± 0.0018
Radius Ratio, Rp /R⋆	0.016 ± 0.001
Impact Parameter, b	0.3 ± 0.2
Duration, T14 (hours)	1.42 ± 0.10
Orbital Eccentricity, e	0 (fixed)
RV Semiamplitude, Kb (m s−1)	3.1 ± 0.8
Semimajor axis, ab (au)	0.01064 ± 0.00013
Radius, Rb (R⊕)	1.45 ± 0.11
Mass, Mb (M⊕)	3.2 ± 0.8
Density, ρb (g cm−3)	${5.5}_{-1.6}^{+2.0}$
Equilibrium Temperature, Teq,b (K)	2480 ± 200
Planet c
Orbital Period, Pc (days)	10.77892 ± 0.00015
Midtransit Time, Tc (BJD)	2458527.05825 ± 0.00053
Radius Ratio, Rp /R⋆	0.032 ± 0.001
Impact Parameter, b	0.2 ± 0.2
Duration T14 (hours)	4.04 ± 0.26
Orbital Eccentricity, e	0 (fixed)
RV Semiamplitude, Kc (m s−1)	2.4 ± 0.8
Semimajor axis, ac (au)	0.0888 ± 0.0011
Radius, Rc (R⊕)	2.90 ± 0.13
Mass, Mc (M⊕)	7.0 ± 2.3
Density, ρc (g cm−3)	1.6 ± 0.6
Equilibrium Temperature, Teq,c (K)	860 ± 70
Planet d
Orbital Period, Pd a (days)	16.287 ± 0.005
Midtransit Time, Tc (BJD)	2458521.8828 ± 0.0035
Radius Ratio, Rp /R⋆	0.0256 ± 0.0016
Impact Parameter, b	0.1 ± 0.1
Duration, T14 (hours)	4.45 ± 0.46
Orbital Eccentricity, e	0 (fixed)
RV Semiamplitude, Kd (m s−1)	0.9 ± 0.6
Semimajor axis, ad a (au)	0.1174 ± 0.0015
Radius, Rd (R⊕)	2.32 ± 0.16
Mass, Md a (M⊕)	${3.0}_{-1.9}^{+2.4}$
Density, ρd a (g cm−3)	${1.3}_{-0.8}^{+1.1}$
Equilibrium Temperature, Teq,d a (K)	750 ± 60
Other
RV Zero Point, γ(m s−1)	−0.8291
RV Jitter, σj (m s−1)	${4.09}_{-0.42}^{+0.49}$

Note.

^a The orbital period of planet d was incorrectly identified as P = 16.37 days in the original SPOC pipeline because of a partially masked transit. An alias of the orbital period of planet d, P_d = 8 days, is also consistent with the data. Parameters marked with a dagger would be affected by an incorrect assumption of the orbital period of planet d. See Section 8 for assumptions regarding planet equilibrium temperatures.

Download table as:  ASCIITypeset image

We found that the transit duration for d is long compared to what is expected for a planet in a circular, 16 day orbit, given the well-characterized stellar density. ⁶³ The long transit duration can be resolved with a moderate eccentricity for planet d. We applied the photoeccentric method of Dawson et al. (2012) to estimate the eccentricity of planet d, finding ${e}_{d}={0.24}_{-0.13}^{+0.27}$. However, this eccentricity would likely result in the instability of the 10 and 16 day planets, in contrast to the stable orbits we find assuming low eccentricities (see Section 8). The too-long duration of planet d only worsens if we assume the P_d = 8 day solution. Another resolution of the long-duration transits of planet d is suggested in a contemporaneous paper by Lacedelli et al. (2020): the two apparent transits of planet d might be single-transit events from two distinct planets, each with P > 16 days. The similar depths of the transits could result from the observed "peas in a pod" pattern, wherein planets in the same system often have similar sizes (Weiss et al. 2018). We do not detect a secure RV signal at P > 16 days that would be consistent with the orbit of either such planet (Section 6). Yet another possibility is that the apparent 2σ tension between the transit duration and orbital period is the result of systematic or random errors in the photometry. Ultimately, additional high-precision photometry is needed to test these various possible explanations for the long durations of planet d. For the rest of this paper, we will assume that the transits are due to a 16 day planet (except where stated otherwise).

We utilized both Gemini-North with NIRI (Hodapp et al. 2003) and Palomar Observatory with PHARO (Hayward et al. 2001) to obtain near-infrared adaptive optics imaging of TOI-561, on 2019 May 24 and 2020 January 8 respectively. Observations were made in the Brγ filter (λ_o = 2.1686; Δλ = 0.0326 μm). For the Gemini data, 9 dithered images with an exposure time of 2.5 s each were obtained; at Palomar, 15 dithered frames with an exposure of 2.8 s each were obtained. In both cases, the telescope was dithered by a few arcseconds between each exposure, and the dithered science frames were used to create a sky background. Data were reduced using a custom pipeline: we removed bad pixels, performed a sky background subtraction and a flat correction, aligned the stellar position between images, and coadded. The final resolution of the combined dithers was determined from the FWHM of the point-spread function: 013 and 010 for the Gemini and Palomar data, respectively.

The sensitivities of the final combined AO images were determined by injecting simulated sources azimuthally around the primary target every 20° at separations of integer multiples of the central source's FWHM (Furlan et al. 2017). The brightness of each injected source was scaled until standard aperture photometry detected it with 5σ significance. The resulting brightness of the injected sources relative to the target sets the contrast limits at that injection location. The final 5σ limit at each separation was determined from the average of all of the determined limits at that separation, and the uncertainty on the limit was set by the rms dispersion of the azimuthal slices at a given radial distance.

The sensitivity curves are shown in Figure 7 along with an inset image zoomed to the primary target showing no other companion stars. Both the Gemini and Palomar data reach a Δmag ≈ 2 at 015 with an ultimate sensitivity of 7.7 mag and 8.7 mag for the Gemini and Palomar imaging, respectively. To within the limits and sensitivity of the data, no additional companions were detected.

Zoom In Zoom Out Reset image size

We also searched for stellar companions with speckle imaging on the 4.1 m Southern Astrophysical Research (SOAR) telescope (Tokovinin et al. 2018) on 2019 May 18 UT. The speckle observations complementing the NIR AO as the I-band observations are similar to the TESS bandpass. More details of the observations are available in Ziegler et al. (2020). The observations have a sensitivity of ∼1 mag at a resolution of 006 and an ultimate sensitivity of ∼7 mag at a radius of 3''. The 5σ detection sensitivity and speckle autocorrelation functions from the observations are shown in Figure 7. As with the NIR AO data, no nearby stars were detected within 3'' of TOI-561 in the SOAR observations.

High-resolution speckle interferometric images of TOI-561 were obtained on 2020 March 15 UT using the Zorro ⁶⁴ instrument mounted on the 8 m Gemini-South telescope located on the summit of Cerro Pachon in Chile. Zorro simultaneously observes in two bands, i.e., 832/40 nm and 562/54 nm, obtaining diffraction-limited images with inner working angles 0017 and 0026, respectively. Our data set consisted of 5 minutes of total integration time taken as sets of 1000 × 0.06 s images. All the images were combined and subjected to Fourier analysis leading to the production of final data products including speckle reconstructed imagery (see Howell et al. 2011). Figure 7 shows the 5σ contrast curves in both filters for the Zorro observation and includes an inset showing the 832 nm reconstructed image. The speckle imaging results reveal TOI-561 to be a single star to contrast limits of 5 to 8 mag, ruling out main-sequence stars brighter than late M as possible companions to TOI-561 within the spatial limits of ∼2 to 103 au (at d = 86 pc).

We modeled the RVs with the publicly available Python package radvel (Fulton et al. 2018). We used the default basis of a five-parameter Keplerian orbit, in which the RV component of each planet is described by its orbital period (P), time of conjunction (T_c ), eccentricity (e), the argument of periastron passage (ω), and RV semiamplitude (K). Because the orbital ephemerides from TESS are more precise than what we can constrain with 60 RVs, we fixed P and T_c for each transiting planet at the best-fit values from TESS photometry plus additional ground-based photometry from the TESS Science Group 1. ⁶⁵ The eccentricities of all three planets are expected to be small for dynamical reasons. At P < 1 day, the USP is almost certainly tidally circularized. The other planets have compact orbits (P_d /P_c ≈ 1.5), such that large eccentricities would result in orbit crossings, and even modest eccentricities would likely result in Lagrange instability (Deck et al. 2013). Furthermore, the majority of small exoplanets in compact configurations have low eccentricities (Mills et al. 2019; Van Eylen et al. 2019). For these reasons, and also because modeling eccentricities introduces two free parameters per planet, we only explored circular fits for all three transiting planets.

Thus, of the five Keplerian parameters that describe each transiting planet, only the semiamplitude (K) was allowed to vary along two global terms: an RV zero-point offset (γ) and an RV jitter (σ_j ), which is added to the individually determined RV errors in quadrature to account for non-Gaussian, correlated noise in the RVs from stellar processes and instrumental systematics. Our full likelihood model was

$$\begin{eqnarray}&&-2\mathrm{ln}{ \mathcal L }=\displaystyle \sum _{i}\displaystyle \frac{{\left({x}_{\mathrm{meas},i}-{x}_{\mathrm{mod},i}\right)}^{2}}{{\sigma }_{i}^{{\prime} 2}}+\displaystyle \sum _{i}\mathrm{ln}(2\pi {\sigma }_{i}^{{\prime} 2}),\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&{\sigma }_{i}^{{\prime} 2}={\sigma }_{i}^{2}+{\sigma }_{j}^{2}\end{eqnarray} \tag{ 3 }$$

is the quadrature sum of the internal RV error and the jitter.

We optimized the likelihood function with the Powell method (Powell 1964) and used a Markov Chain Monte Carlo (MCMC) analysis ⁶⁶ to determine parameter uncertainties. We explored the optimization of several models. In Model A, we did not enforce any priors on planet semiamplitudes (thus allowing values of K, and hence planet mass, to be negative). Although negative planet masses are unphysical, their consideration offsets the bias toward high planet masses that occurs when planet masses are forced to be positive (Weiss & Marcy 2014). The best-fit values with Model A were K_b = 3.1 ± 0.8 m s⁻¹, K_c = 2.2 ± 0.8 m s⁻¹, and K_d = 0.3 ± 0.8 m s⁻¹. The rms of the RV residuals was 4.2 m s⁻¹.

In Model B, we restricted K > 0. The advantages of restricting planet masses to be larger than zero are (1) the planet masses are physically motivated, and (2) the residuals are more likely to be useful in searching for additional planets. Model B yielded K_b = 3.1 ± 0.8 m s⁻¹, K_c = 2.4 ± 0.8 m s⁻¹, and K_d = 0.9 ± 0.6 m s⁻¹. The rms of the RV residuals was 4.2 m s⁻¹.

Model C was the same as Model B, except we allowed a linear trend in the RVs, $\dot{\gamma }$, which could be caused by acceleration from a long-period companion. However, the RVs do not strongly favor a trend: $\dot{\gamma }$ = 0.009 ± 0.004 m s⁻¹ m s⁻¹ day⁻¹, and the best-fit K values changed by less than 1σ with the inclusion of a trend: K_b = 3.2 ± 0.8 m s⁻¹, K_c = 2.1 ± 0.8 m s⁻¹, and ${K}_{d}\,={1.0}_{-0.6}^{+0.7}\,{\rm{m}}\,{{\rm{s}}}^{-1}$. The rms of the RV residuals was 4.0 m s⁻¹.

In Model D, we considered the hypothesis that planet d has half of the presumed orbital period, which was possible given the gap in the photometry. This model is the same as Model B, except P_d = 8 days. This model did not affect the amplitudes of planets b or c, but resulted in K_d < 2.04 m s⁻¹ (2σ confidence). The rms of the RV residuals was 4.1 m s⁻¹.

In Models A to C, the choice of model makes very little impact on the best-fitting RV semiamplitudes for each planet, and hence our planet mass determinations are robust with respect to our choice of model. For the rest of this paper, we consider Model B as our default model unless stated otherwise. The fitted and fixed parameters of Model B are provided in Table 3.

The simple Keplerian model fit to observed RVs of TOI-561 displayed an rms of 4.2 m s⁻¹ (Figure 8), which is uncharacteristically high for precision RVs of such a bright star (Howard & Fulton 2016). A Lomb-Scargle periodogram of the RV residuals reveals several peaks that might correspond to correlated noise and/or additional planets in the system. However, none of the peaks were significant at the 1% false-alarm probability (FAP) level, which we determined with bootstrap resampling (Figure 9).

Zoom In Zoom Out Reset image size

Nonetheless, we considered several correlated-noise models in an attempt to model and remove a putative red-noise component responsible for the high rms of the RV residuals. In Model E, we employed Gaussian process regression (GP), which has been previously applied in analyzing RVs of many exoplanets (e.g., Haywood et al. 2014; Grunblatt et al. 2015). For details of the GP model, see Dai et al. (2017). In principle, the light curve and RVs are both affected by stellar activity rotating in and out of view of the observer (e.g., Aigrain et al. 2012). In an attempt to model the correlated noise in the RVs, we trained the various "hyperparameters" of our GP quasi-periodic kernel using the out-of-transit TESS PDCSAP light curve. We did not detect the stellar rotation period in the TESS PDCSAP or SAP light curves, possibly because the star is very inactive (log ${R}_{\mathrm{HK}}^{{\prime} }=-5.1$), or the expected rotation period of the star is longer than the photometric baseline. Our MCMC fit to the light curve produced broad posteriors on the rotation period and the other hyperparameters. Therefore, we did not impose Gaussian priors on the hyperparameters in our GP analysis of the RV data, although we did limit the stellar rotation period <500 days for quicker convergence. We used emcee to constrain the posterior distribution of the GP hyperparameters simultaneously with the orbital parameters of the three planets. We detect the RV signal of planet b: K_b = 2.9 ± 0.7 m s⁻¹ and planet c: K_b = 1.7 ± 1.0 and an upper limit for planet d: K_d < 1.6 m s⁻¹ (95% confidence), all of which are within 1σ of the planet semiamplitudes we determined without the correlated-noise component of the model.

We also attempted to train a correlated-noise model on our spectroscopically determined Mt. Wilson S values. A Lomb-Scargle periodogram of the S values shows two significant peaks: one at 100 days and the other at 230 days, both of which cross the 1% FAP threshold (which we computed with bootstrap resampling, Figure 9). (The forest of regular peaks with P < 1 day are likely aliases produced by our window function, which had poor sampling below 1 day). Note that these periods differ from the most prominent peaks in the RV residuals, which are a doublet near 25 days. Furthermore, the S values are not sampled as frequently as the light curve, but they are sampled simultaneously with the RVs, and they are a direct indicator of the chromospheric magnetic activity during the observations. We tried a GP with a quasi-periodic kernel trained on the S values (Model F), but this did not produce significant changes in the semiamplitudes of the planets or the rms of the RV residuals, possibly because the S values had small variability or were sparsely sampled. We also tried a model in which we decorrelated the RVs with respect to the S values (Model G), which did not reduce the rms (and thus did not remove the correlated noise).

None of our attempts to model correlated noise in the RVs changed the amplitudes of the planets or reduced the rms of the RV residuals, and so we prefer the simpler Keplerian models (without correlated noise) described above. Perhaps TOI-561 is too inactive for models trained on stellar activity to be effective, given the current quality of the data. For comparison, Kepler-10, another system with time-correlated RV residuals, has log ${R}_{\mathrm{HK}}^{{\prime} }=-4.89$, which is more active than TOI-561 log ${R}_{\mathrm{HK}}^{{\prime} }=-5.1$. The use of GPs affected the mass determination of Kepler-10 c, lowering it from 14 M_⊕ (no GP) to 7 M_⊕ (with GP). Perhaps there is a minimum stellar activity for which attempts to decorrelate the stellar activity signal can be successful, given RVs with precision of 2 m s⁻¹.

Nonetheless, there are substantial correlated residuals in the RVs of TOI-561 which are uncharacteristic of the HIRES instrument performance (typically 2 m s⁻¹ for V < 11; Howard & Fulton 2016). The residual RVs of TOI-561 are not well explained by any of our models of stellar activity, and so perhaps additional planets contribute to the RV residuals. More RVs are needed to identify the orbital periods of any such planets and model their Keplerian signals.

The USP TOI-561 b has a typical mass and density for its size. At 1.5 R_⊕ and ${5.5}_{-1.6}^{+2.0}$ g cm⁻³, it is 1σ below the peak of the density–radius diagram identified in Weiss & Marcy (2014), consistent with a rocky composition that is either Earth like or iron poor. Nearly all USPs are smaller than 2 R_⊕ and are expected to have rocky compositions, given their small sizes and extreme stellar irradiation (Sanchis-Ojeda et al. 2015), and TOI-561 b is consistent with this expectation. In a homogeneous analysis of USPs with masses determined from RVs, Dai et al. (2019) found that most USPs with <10 M_⊕ are consistent with having Earth-like compositions, whereas the few USPs with >10 M_⊕ likely have H/He envelopes. TOI-561 b (3.2 ± 0.8 M_⊕) is consistent with the rocky group of that study.

The minimum density of a USP can be determined from its orbital period and the requirement that it orbits outside the Roche limiting distance (Rappaport et al. 2013). We investigated the minimum density of TOI-561 b, with the hope that it would provide additional constraints on the mass of the planet. Using the approximation from Sanchis-Ojeda et al. (2014), we find that the minimum density of the USP is

$$\begin{eqnarray}&&{\rho }_{{\rm{p}}}\ [{\rm{g}}\,{\mathrm{cm}}^{-3}]\geqslant {\left(11.3\mathrm{hr}/{P}_{\mathrm{orb}}\right)}^{2}=1.15{\rm{g}}\,{\mathrm{cm}}^{-3}\end{eqnarray} \tag{ 4 }$$

for TOI-561 b, which is below our measured density and corresponds to a minimum mass of 0.64 M_⊕. Such a low mass is ruled out by the data at nearly 3σ confidence. Thus, TOI-561 b is not close enough to its star for the Roche stability criterion to provide additional information about the density.

An interpretation of the USP composition that involves a H/He layer is disfavored by physical models. At a sufficiently low stellar irradiation level, the low density of TOI-561 b might have been consistent with a composition of an Earth-like core overlaid with a thin H/He envelope. However, at P = 0.44 days (with equilibrium temperature >2000 K) and a mass of ∼3 M_⊕, the USP is too irradiated and too low mass to hold onto a H/He envelope (Lopez 2017).

At R_p > 1.5 R_⊕ and with low densities, TOI-561 c and d have substantial gaseous envelopes by volume, although the gas envelopes likely only constitute ∼1% of the planet masses (Lopez & Fortney 2014). TOI-561 c has a radius and mass consistent with the Weiss & Marcy (2014) empirical mass–radius relationship.

The ambiguity of the orbital period for planet d poses a challenge to accurate mass determination. Our RVs are consistent with a nondetection of planet d at both the P_d = 16 and P_d = 8 day orbits. Assuming P_d = 16 days, the RVs provide an upper limit of K_d < 2.1 m s⁻¹, which corresponds to M_p < 7.0 M_⊕ (2σ confidence; see Section 6). Assuming P_d = 8 days, the RVs provide an upper limit of K < 2.0 m s⁻¹, which corresponds to M_p < 5.6 M_⊕ (2σ confidence). In either scenario, the mass of TOI-561 d is approximately 2σ below the Weiss & Marcy (2014) mass–radius relationship and is too low to be consistent with a rocky composition, given the planet's radius. However, if planet d is actually the transits of two distinct planets as suggested in Lacedelli et al. (2020), the low mass presented here does not apply.

With our 1.5 year baseline of RVs, we were able to place constraints on possible long-period giant planets. In Model C, we found a 3σ upper limit to an RV trend of $\dot{\gamma }\lt 0.02{\rm{m}}\,{{\rm{s}}}^{-1}\ {\mathrm{day}}^{-1}$. This observational limit on the acceleration of the star can be converted to a limit on the mass and orbital distance of a possible perturber by setting it equal to the gravitational acceleration from a planet at distance r: ${{GM}}_{\star }{M}_{{\rm{p}}}\ {r}^{-2}={M}_{\star }\dot{\gamma }$, where G is the gravitational constant and r is the distance between the perturber and the star at the time the RVs were measured. Assuming a circular orbit for the putative giant planet (r = a), we find

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{{\rm{p}}}\sin i}{{M}_{{\rm{J}}}}{\left(\displaystyle \frac{a}{5\mathrm{au}}\right)}^{-2}\lt 1.0\end{eqnarray} \tag{ 5 }$$

(3σ confidence). Thus, our nondetection of an RV acceleration rules out a 1.0 M_J planet at 5 au with 3σ confidence (assuming sin i ≈ 1, e = 0, and that we did not primarily sample the orbit while the planet was moving parallel to the sky plane). In systems with compact configurations of transiting planets, giant planets often need to be approximately coplanar with the transiting planets in order for the inner planets to remain mutually continuously transiting (Becker & Adams 2017), and so we do not expect nontransiting companions to have face-on orbits. Thus, the nondetection of an RV trend rules out a variety of scenarios of a coplanar giant planet near the snow line.

We investigated whether we could constrain the architecture of the TOI-561 system and in particular the orbital period of planet d, with stability arguments. To assess stability, we used Spock, a machine-learning-based approach to inferring orbital stability (Tamayo et al. 2020). Spock incorporates several analytic indicators (including MEGNO, AMD, and mutual Hill radius) with an N-body integration of 10⁴ orbits to compute a probabilistic assessment of the system stability. We tested the architectures of TOI-561 with planet d at P_d = 16 days (Model B) and at P_d = 8 days (Model D). Because stability depends sensitively on the eccentricities and inclinations of planets in compact configurations, we varied the eccentricities, arguments of periastron passage, and inclinations. The eccentricities and inclinations in multiplanet systems are well fit with Rayleigh distributions with characteristic values σ_e = 0.035 and σ_i = 2.45° (Mills et al. 2019), and so we sampled 1000 trials from these distributions (sampling ω uniformly on [0,π]). For Model B, the probability of stability is 95% ± 3%, whereas for Model D, the probability of stability is 84% ± 30%. Thus, although more configurations of Model B are stable, we cannot rule out Model D based on a stability argument.

TOI-561 is the second transiting multiplanet system discovered around a galactic thick-disk star (Kepler-444, Campante et al. 2015, is the first), and is the first galactic thick-disk star with a USP. The iron-poor, α-enhanced stellar abundances observed today are likely representative of the nebular environment in which the planets formed. Thus, TOI-561 provides an opportunity to study the outcome of planet formation in an environment that is chemically distinct from most of the planetary systems known to date. Furthermore, its membership in the galactic thick disk indicates old age, making the rocky planet TOI-561 b one of the oldest rocky planets known.

The confirmed old age of the system is relevant for dynamical studies of the planets. For instance, the formation of USPs is still poorly understood. The majority of USPs are the only detected transiting planet in their systems (Sanchis-Ojeda et al. 2014). This is partially due to a bias of geometry (as planet detection probability scales with R_⋆/a). However, Dai et al. (2018) found that the mutual inclinations in multiplanet systems with USPs are significantly larger than those in multiplanet systems without USPs, suggesting that systems with USPs have had more dynamically "hot" histories. One mechanism for generating large mutual inclinations in systems with USPs is if the host star is oblate and misaligned with respect to the planets (Li et al. 2020).

Furthermore, dynamical interactions in a multiplanet system can move short-period planets to ultra-short periods in a manner that may excite large eccentricities. Petrovich et al. (2019) proposed a mechanism of secular chaos in which the innermost planet is kicked to a high-eccentricity orbit, which is then circularized. In contrast, Pu & Lai (2019) proposed a scheme in which the innermost planet's neighbors consistently force it to a low-eccentricity orbit, which results in inward migration and eventual circularization.

In a study of the Gaia-DR2 kinematics of USP host stars in the Kepler field, Hamer & Schlaufman (2020) found that their motions are similar to those of matched field stars (rather than young stars). The broad range of ages of USPs suggests that USPs do not undergo rapid tidal inspiral during the host star's main-sequence lifetime. The existence of a USP around an 10 Gyr old star is consistent with this finding.

Because planets c and d have an orbital period ratio of ∼1.5 or less, the planets likely perturb each other's orbits, producing transit timing variations (TTVs). The TESS sector 8 baseline was too short to detect TTVs (only two transits of each planet were detected). The approximate amplitude of the TTV signal can be computed analytically using the expressions from Lithwick et al. (2012):

$$\begin{eqnarray}&&| V| \approx P\mu ^{\prime} /{\rm{\Delta }}+{ \mathcal O }(e)/{\rm{\Delta }},\end{eqnarray} \tag{ 6 }$$

where P is the orbital period of the inner planet, $\mu ^{\prime}$ is the mass of the perturbing, outer planet (in units of stellar mass), ${ \mathcal O }(e)$ is a first-order dependency on the free eccentricities of the planets, and Δ is the nondimensional distance from mean motion resonance:

$$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{P^{\prime} }{P}\displaystyle \frac{j-1}{j}-1,\end{eqnarray} \tag{ 7 }$$

where P and $P^{\prime}$ are the inner and outer orbital periods and j is an integer. Equation (6) can be used to approximate the TTV amplitude of the outer planet by setting $P\longrightarrow P^{\prime}$ and $\mu ^{\prime} \longrightarrow \mu$ (for a thorough derivation and caveats, see Lithwick et al. 2012).

Using the values for the planet masses determined in Section 7 and assuming circular orbits, we find that the TTV amplitude of planet c (P = 10.78 days) is about 6 minutes, whereas the TTV amplitude of planet d is either 30 minutes (if P_d = 16 days; j = 3) or 15 minutes (if P_d = 8.14 days, j = 4).

Using the system parameters tabulated in this paper, we calculated the transmission and emission spectroscopy metrics (TSM and ESM, respectively) of Kempton et al. (2018) to determine whether these newly characterized planets are compelling targets for future atmospheric or surface characterization via transit or eclipse spectroscopy. Assuming a flat prior on the planets' Bond albedos from 0 to 0.4 and on their efficiency of day-to-night energy recirculation from 1/4 to 2/3, and propagating the uncertainties on all relevant parameters, we calculate the planets' equilibrium temperatures to be 2480 ± 200 K, 860 ± 70 K, and 750 ± 60 K for planets b, c, and d, respectively. Because of the relatively shallow transit depths (≲1000 ppm for all three planets), space-based spectroscopy is likely to be the only feasible avenue for such studies. In contrast to some recent reports of these quantities for other newly discovered TESS planets, here we propagate all parameter uncertainties in order to report how well the TSM and ESM are constrained, as well as how promising the median values are; this is especially essential when planetary properties are not yet measured to high precision. We report the transmission and emission metrics and their uncertainties in Table 4.

Table 4. Atmospheric Prospects

Planet	TSM	ESM	Notes
b	${11}_{-4}^{+73}$	8.7 ± 1.1	Good eclipse target
c	${93}_{-37}^{+55}$	4.87 ± 0.39	Promising transmission target
d a	${82}_{-33}^{+109}$	2.15 ± 0.29	Promising transmission target

Note.

^a Assuming P_d = 16 days. If P_d = 8 days then TSM${}_{d}={103}_{-43}^{+137}$ (68% confidence interval).

Download table as:  ASCIITypeset image

Our analysis shows that TOI-561 b is among the best TESS targets discovered to date for thermal emission measurements (see Table 3 of Astudillo-Defru et al. 2020). With ESM = 7.1 ± 1.1, planet b is clearly a promising target for observations of its secondary eclipse and/or its full-orbit phase curves, as has previously been done for other irradiated terrestrial planets such as 55 Cnc e (Demory et al. 2012, 2016) and LHS 3844b (Kreidberg et al. 2019). The uncertainty on planet b's ESM is dominated by the uncertainty on its transit depth, but regardless, the planet has a reliably high metric in this category. Because of their cooler temperatures, the lower ESM values for planets b and c mark them as less attractive targets for secondary eclipse studies.

As for transit spectroscopy, the TSM values for the sub-Neptunes TOI-561 c and d listed in Table 4 (${93}_{-37}^{+55}$ and ${82}_{-33}^{+109}$, respectively) indicate that they may be particularly amenable to transmission studies. However, because the TSM scales inversely with planetary surface gravity, this result depends on determining more precise values of these planets' masses. These planets clearly warrant additional precise RV follow-up: if the expectation values of their TSMs do not change as the uncertainties shrink, these two planets would be among the top 20 confirmed warm Neptunes for transmission spectroscopy (see Table 11 of Guo et al. 2020). Due to the small size of the highly irradiated planet b, and because it is unlikely to have retained much of an atmosphere, it is not an appealing target for transmission measurements.

Better ephemerides for planets c and d are necessary in preparation for atmospheric studies. The alias of the orbit of planet d should be resolved prior to interpretation of the planetary atmospheres, as the factor of 2 change in the orbital period produces a factor of ∼1.3 change in the equilibrium temperature. Also, planet d may have significant TTVs with amplitudes of ∼30 minutes (assuming the orbital period is 16 days).