TOI-1235 b: A Keystone Super-Earth for Testing Radius Valley Emergence Models around Early M Dwarfs

TOI-1235 was observed in three nonconsecutive TESS sectors between UT 2019 July 18 and 2020 February 18. TOI-1235 is a member of the Cool Dwarf target list (Muirhead et al. 2018) and was included in the TESS Input Catalog (TIC; Stassun et al. 2018b), the TESS Candidate Target List (CTL), and the Guest Investigator program 22198,⁴⁷ such that its light curve was sampled at 2-minute cadence. TESS observations occurred in CCD 3 on Camera 4 in sector 14 (UT 2019 July 18–August 14), in CCD 1 on Camera 2 in sector 20 (UT 2019 December 24–2020 January 20), and in CCD 2 on Camera 2 in sector 21 (UT 2020 January 21–February 18). Sector 14 was the first pointing of the spacecraft in the northern ecliptic hemisphere. As indicated in the data release notes,⁴⁸ to avoid significant contamination in cameras 1 and 2 due to scattered light by Earth and the Moon, the sector 14 field was pointed to +85° in ecliptic latitude, 31° north of its intended pointing from the nominal mission strategy. Despite this, all cameras in sector 14 continued to be affected by scattered light for longer periods of time compared to most other sectors owing to Earth's position above the sunshade throughout the orbit. Camera 2 during sectors 20 and 21 was largely unaffected by scattered light except during data downloads and at the beginning of the second orbit in sector 21 owing to excess Moon glint.

The TESS images were processed by the NASA Ames Science Processing Operations Center (SPOC; Jenkins et al. 2016), which produce two light curves per sector called Simple Aperture Photometry (SAP) and Presearch Data Conditioning Simple Aperture Photometry (PDCSAP; Smith et al. 2012; Stumpe et al. 2012, 2014). The light curves are corrected for dilution during the SPOC processing, with TOI-1235 suffering only marginal contamination with a dilution correction factor of 0.9991. Throughout, we only consider reliable TESS measurements for which the measurement's quality flag QUALITY is equal to zero. The PDCSAP light curve is constructed by detrending the SAP light curve using a linear combination of cotrending basis vectors (CBVs), which are derived from a principal component decomposition of the light curves on a per-sector, per-camera, per-CCD basis. TOI-1235's PDCSAP light curve is depicted in Figure 2 and shows no compelling signs of coherent photometric variability from rotation. However, the set of CBVs (not shown) exhibits sufficient temporal structure that a linear combination of CBVs can effectively mask stellar rotation signatures greater than a few days. Thus, inferring ${P}_{\mathrm{rot}}$ for TOI-1235 from TESS would be challenging and is addressed more effectively with ground-based photometric monitoring in Section 3.2.

Zoom In Zoom Out Reset image size

Following light-curve construction, the SPOC conducts a subsequent transit search on each sector's PDCSAP light curve using the Transiting Planet Search (TPS) module (Jenkins 2002; Jenkins et al. 2010). The TOI-1235.01 transit-like signal was detected in all three sectors independently and passed a set of internal data validation tests (Twicken et al. 2018; Li et al. 2019). The reported period of the planet candidate was 3.44 days in sectors 14 and 20 and three times that value (i.e., 10.33 days) in sector 21 owing to the low signal-to-noise ratio (S/N) of the individual transits. At 3.44 days, there are eight, six, and eight transits observed in each of the three sectors. The transit events are highlighted in Figure 2. The SPOC reported a preliminary transit depth of 841 ± 72 ppm, which corresponded to a planetary radius of 2.0 ± 0.1 ${R}_{\oplus }$ using our stellar radius (Table 1).

Inactive early M dwarfs have typical rotation periods of 10–50 days (Newton et al. 2017). In Section 3.1 we described how measuring ${P}_{\mathrm{rot}}$ for TOI-1235 with TESS is intractable owing to the flexibility in the systematics model. Fortunately, MEarth-North has archival images of the field surrounding TOI-1235 that span 7.1 yr (UT 2008 October 2–2015 November 10), from which ${P}_{\mathrm{rot}}$ may be measured. MEarth-North is a telescope array located at the Fred Lawrence Whipple Observatory (FLWO) on Mount Hopkins, AZ. The facility consists of eight 40 cm telescopes, each equipped with a 256 × 256 field-of-view Apogee U42 camera, with a custom passband centered in the red optical (i.e., RG715). MEarth-North has been photometrically monitoring nearby mid- to late M dwarfs (<0.33 ${R}_{\odot }$) since 2008, in search of transiting planets (Berta et al. 2012; Irwin et al. 2015) and to conduct detailed studies of stellar variability (Newton et al. 2016). Although TOI-1235 was too large to be included in the initial target list (Nutzman & Charbonneau 2008), its position happens to be within 14' of an intentional target (GJ 1131) such that we are able to construct and analyze its light curve here for the first time.

To search for photometric signatures of rotation, we first retrieved the archival image sequence and computed the differential light curve of TOI-1235 as shown in Figure 3. We then investigated the Lomb–Scargle periodogram (LSP) of the light curve, which reveals a significant peak around 45 days that is not visible in the LSP of the window function (Figure 3). Using this value as an initial guess, we proceeded with fitting the light curve following the methods outlined in Irwin et al. (2006, 2011). The model includes systematics terms, predominantly from variations in the precipitable water vapor (PWV) column above the telescope, plus a sinusoidal term to model rotational modulation. As outlined in Newton et al. (2016), a "common mode" vector is constructed as a low-cadence comparison light curve that tracks variations in the PWV and is included in our systematics model as a linear term along with the FWHM of the MEarth point-spread function. With this full model, we measure ${P}_{\mathrm{rot}}$ = 44.7 ± 4.5 days and a variability semiamplitude of 1.33 ppt. The detrended light curve, phase-folded to ${P}_{\mathrm{rot}}$, is included in Figure 3. Figure 3 also reveals that the subtraction of our systematics plus rotation model from the light curve mitigates the 45-day signal in the LSP with no significant residual periodicities. The shallow variability amplitude is unsurprising for relatively warm early M dwarfs like TOI-1235, whose spot-to-photosphere temperature contrasts are small (Newton et al. 2016). We note that knowledge of ${P}_{\mathrm{rot}}$ can be critical for the interpretation of RV signals, as even active regions with small temperature contrasts can induce large RV variations due to the suppression of convective blueshift (Dumusque et al. 2014).

Zoom In Zoom Out Reset image size

We began to pursue the confirmation of the planet candidate TOI-1235.01 by obtaining reconnaissance spectra with the Tillinghast Reflector Échelle Spectrograph (TRES) through coordination with the TESS Follow-up Observing Program (TFOP). TRES is a fiber-fed R = 44,000 optical échelle spectrograph (310–910 nm), mounted on the 1.5 m Tillinghast Reflector telescope at FLWO. Multiple spectra were obtained to search for RV variations indicative of a spectroscopic binary and to assess the level of surface rotation and chromospheric activity. We obtained two spectra at opposite quadrature phases of TOI-1235.01 on UT 2019 December 1 and 13 with exposure times of 2100 and 1200 s, which resulted in an S/N per resolution element of 31.4 and 26.0, respectively, at 519 nm in the order containing the information-rich Mg b lines.

The TRES RVs phase-folded to the TOI-1235.01 ephemeris are depicted in Figure 4 and show no significant variation, thus ruling out a spectroscopic binary. The cross-correlation function of the median spectrum with a rotating template of Barnard's star is also shown in Figure 4 and reveals a single-lined spectrum with no significant rotational broadening ($v\sin i$ < 3.4 km s⁻¹). Lastly, the Hα feature shown is seen in absorption, which is indicative of a chromospherically inactive star and is consistent with ${P}_{\mathrm{rot}}$ ≳ 10 days (Newton et al. 2017). Taken together, our reconnaissance spectra maintain that TOI-1235.01 is a planetary candidate around a relatively inactive star.

Zoom In Zoom Out Reset image size

TESS's large pixels (21'') can result in significant blending of target light curves with nearby sources. To confirm that the transit event occurs on-target, and to rule out nearby eclipsing binaries (EBs), we targeted a transit of TOI-1235.01 with seeing-limited photometric follow-up on UT 2019 December 31. This observation was scheduled after the planet candidate was detected in TESS sector 14 only and occurred during sector 20. The transit observation was scheduled using the TESS Transit Finder, a customized version of the Tapir software package (Jensen 2013). We obtained a z_s-band light curve from the McDonald Observatory with the 1 m telescope as part of the Las Cumbres Observatory Global Telescope network (LCOGT; Brown et al. 2013). The telescope is equipped with a 4096 × 4096 Sinistro camera whose pixel scale is 54 times finer than that of TESS: 0389 pixel⁻¹. We calibrated the full image sequence using the standard LCOGT BANZAI pipeline (McCully et al. 2018). The differential photometric light curve of TOI-1235, along with seven sources within 25, was derived from 7'' uncontaminated apertures using the AstroImageJ software package (AIJ; Collins et al. 2017). The field was cleared of nearby EBs down to Δz_s = 7.15, as we did not detect eclipses from neighboring sources close to the expected transit time.

A full transit event was detected on-target and is included in Figure 4. We fit the light curve with a Mandel & Agol (2002) transit model calculated using the batman software package (Kreidberg 2015). The shallow transit depth of TOI-1235.01 produces a low-S/N transit that does not provide strong constraints on most model parameters relative to what can be recovered from 22 transits in TESS. Consequently, we fix the model to a circular orbit with an orbital period, scaled semimajor axis, and impact parameter of P = 3.44471 days, a/R_s = 13.2, and b = 0.45, respectively. Furthermore, we set the quadratic limb-darkening parameters in the z_s band to u₁ = 0.25 and u₂ = 0.33 as interpolated from the Claret & Bloemen (2011) tables using the EXOFAST tool (Eastman et al. 2013). We fit the baseline flux, time of midtransit, and planet-to-star radius ratio via nonlinear least-squares optimization using the scipy.curve_fit function and find that f₀ = 1.000, T₀ = 2,458,848.962 BJD, and r_p/R_s = 0.0295. The transit is seen to arrive 63 minutes late relative to the linear ephemeris reported by the SPOC from sector 14 only. The transit depth of 0.867 ppt is 4.5σ deeper than the TESS transit measured in our fiducial analysis (0.645 ppt; Section 4.1). Due to the similar wavelength coverage between the z_s and TESS passbands, and because of the large residual systematics often suffered by ground-based light curves of shallow transits, we attribute this discrepancy to unmodeled systematics rather than to a bona fide chromatic transit depth variation.

TESS's large pixels also make the TESS light curves susceptible to contamination by very nearby sources that are not detected in Gaia DR2, nor in the seeing-limited image sequences. To clear the field of very nearby sources and a possible false positive in the form of a blended EB (Ciardi et al. 2015), we obtained two independent sets of high-resolution follow-up imaging sequences as described in the following sections.

3.5.1. Adaptive Optics Imaging with Gemini/NIRI

3.5.2. Speckle Imaging with Gemini/'Alopeke

3.6.1. HARPS-N

We obtained 27 spectra of TOI-1235 with the HARPS-N optical échelle spectrograph at the 3.6 m Telescopio Nazionale Galileo on La Palma in the Canary Islands. The HARPS-N optical spectrograph, with a resolving power of R = 115,000, is stabilized in pressure and temperature, which enable it to achieve sub–meter-per-second accuracy under ideal observing conditions when sufficient S/N is attainable (Cosentino et al. 2012). The spectra were taken as part of the HARPS-N Collaboration Guaranteed Time Observations program between UT 2019 December 24 and 2020 March 12. The exposure time was set to 1800 s. In orders redward of order 18 (440–687 nm), we achieved a median S/N of 45.2 and a median measurement uncertainty of 1.22 m s⁻¹. TOI-1235 did not exhibit any rotational broadening in the HARPS-N spectra, leading to $v\sin i$ ≤ 2.6 km s⁻¹, a result that is consistent with its measured rotation period ${P}_{\mathrm{rot}}$ = 44.7 ± 4.5 days.

We extracted the HARPS-N RVs using the TERRA pipeline (Anglada-Escudé & Butler 2012). TERRA employs a template-matching scheme that is known to achieve improved RV measurement uncertainties on M dwarfs relative to the cross-correlation function (CCF) technique (Anglada-Escudé & Butler 2012). M dwarfs are particularly well suited to RV extraction via template matching because the line lists used to define the binary mask for the CCF technique are incomplete and often produce a CCF template that is a poor match for cool M dwarfs. A master template spectrum is constructed by first shifting the individual spectra to the barycentric frame using the barycentric corrections calculated by the HARPS-N Data Reduction Software (DRS; Lovis & Pepe 2007), after masking portions of the wavelength-calibrated spectra wherein telluric absorption is ≥1%. A high-S/N template spectrum is then built by co-adding the individual spectra. TERRA then computes the RV of each spectrum relative to the template via least-squares matching the spectrum in velocity space. Throughout the extraction process, we only consider orders redward of order 18 such that the bluest orders at low S/N are ignored. The resulting RV time series is provided in Table 2.

Table 2.  Radial Velocity Time Series of TOI-1235 from HARPS-N and HIRES

Time	RV	σRV	Instrument
(BJD −2,457,000)	(m s−1)	$({\rm{m}}\ {{\rm{s}}}^{-1})$
1890.653258	−0.119	0.975	HARPS-N
1905.851683	−7.358	1.281	HIRES
1906.724763	1.803	1.470	HARPS-N

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

3.6.2. HIRES

We begin our fiducial analysis by modeling the TESS PDCSAP light curve (Figure 2) in which the planet candidate TOI-1235.01 was originally detected. The PDCSAP light curve has already undergone systematics corrections via a linear combination of CBVs; however, some low-amplitude temporally correlated signals that are unrelated to planetary transits are seen to persist. We elect to model these signals as an untrained semiparametric Gaussian process (GP) simultaneously with the transit model of TOI-1235 b. We employ the exoplanet software package (Foreman-Mackey et al. 2019) to construct the GP and transit model in each step in our Markov Chain Monte Carlo (MCMC) simulation. Within exoplanet, analytical transit models are computed using the STARRY package (Luger et al. 2019), while celerite (Foreman-Mackey et al. 2017) is used to evaluate the marginalized likelihood of the GP model.

We adopt a covariance kernel of the form of a stochastically driven simple harmonic oscillator in Fourier space. The power spectral density of the kernel is

$$\begin{eqnarray}&&S(\omega )=\sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{{S}_{0}{\omega }_{0}^{4}}{{\left({\omega }^{2}-{\omega }_{0}^{2}\right)}^{2}+{\omega }_{0}^{2}{\omega }^{2}/{Q}^{2}},\end{eqnarray} \tag{ 1 }$$

which is parameterized by the frequency of the undamped oscillator ω₀; the factor S₀, which is proportional to the spectral power at ω₀; and the fixed quality factor $Q=\sqrt{0.5}$. We also include the baseline flux f₀ and an additive scalar jitter s_TESS in our noise model that we parameterize by $\{\mathrm{ln}{\omega }_{0},\mathrm{ln}{S}_{0}{\omega }_{0}^{4},{f}_{0},\mathrm{ln}{s}_{\mathrm{TESS}}^{2}\}$. Our noise model is jointly fit with a transit model for TOI-1235 b with the following free parameters: the stellar mass M_s, stellar radius R_s, quadratic limb-darkening coefficients {u₁, u₂}, orbital period P, time of midtransit T₀, planet radius r_p, impact parameter b, eccentricity e, and argument of periastron ω. Our full TESS model therefore contains 13 free parameters that are parameterized by $\{{f}_{0},\mathrm{ln}{\omega }_{0},\mathrm{ln}{S}_{0}{\omega }_{0}^{4},\mathrm{ln}{s}_{\mathrm{TESS}}^{2},{M}_{{\rm{s}}},{R}_{{\rm{s}}},{u}_{1},{u}_{2},\mathrm{ln}P,{T}_{0},\mathrm{ln}{r}_{{\rm{p}}},$ b, e, ω}. Our adopted model parameter priors are listed in Table 3.

Table 3.  TESS Light Curve and RV Model Parameter Priors

Parameter	Fiducial Model Priors	EXOFASTv2Model Priors
Stellar Parameters
${T}_{\mathrm{eff}}$ (K)	${ \mathcal N }(3872,70)$	${ \mathcal N }(3872,70)$
Ms (${M}_{\odot }$)	${ \mathcal N }(0.640,0.016)$	${ \mathcal N }(0.640,0.016)$
Rs (${R}_{\odot }$)	${ \mathcal N }(0.630,0.015)$	${ \mathcal N }(0.630,0.015)$
Light-curve Hyperparameters
f0	${ \mathcal N }(0,10)$	${ \mathcal U }(-\inf ,\inf )$
$\mathrm{ln}{\omega }_{0}$ (day−1)	${ \mathcal N }(0,10)$	⋯
$\mathrm{ln}{S}_{0}{\omega }_{0}^{4}$	${ \mathcal N }(\mathrm{lnvar}({{\bf{f}}}_{\mathrm{TESS}}),10)$	⋯
$\mathrm{ln}{s}_{\mathrm{TESS}}^{2}$	${ \mathcal N }(\mathrm{lnvar}({{\bf{f}}}_{\mathrm{TESS}}),10)$	⋯
u1	${ \mathcal U }(0,1)$	${ \mathcal U }(0,1)$
u2	${ \mathcal U }(0,1)$	${ \mathcal U }(0,1)$
Dilution	⋯	${ \mathcal N }(0,0.1\ \delta )$ a
RV Parameters
$\mathrm{ln}\lambda$ (days)	${ \mathcal U }(\mathrm{ln}1,\mathrm{ln}1000)$	⋯
$\mathrm{ln}{\rm{\Gamma }}$	${ \mathcal U }(-3,3)$	⋯
${P}_{\mathrm{rot}}$ (days)	${ \mathcal N }(46.1,4.6)$	⋯
$\mathrm{ln}{a}_{\mathrm{HARPS} \mbox{-} {\rm{N}}}$ (m s−1)	${ \mathcal U }(-5,5)$	⋯
$\mathrm{ln}{a}_{\mathrm{HIRES}}$ (m s−1)	${ \mathcal U }(-5,5)$	⋯
$\mathrm{ln}{s}_{\mathrm{HARPS} \mbox{-} {\rm{N}}}$ (m s−1)	${ \mathcal U }(-5,5)$	${ \mathcal U }(-\inf ,\inf )$
$\mathrm{ln}{s}_{\mathrm{HIRES}}$ (m s−1)	${ \mathcal U }(-5,5)$	${ \mathcal U }(-\inf ,\inf )$
${\gamma }_{\mathrm{HARPS} \mbox{-} {\rm{N}}}$ (m s−1)	${ \mathcal U }(-10,10)$	${ \mathcal U }(-\inf ,\inf )$
${\gamma }_{\mathrm{HIRES}}$ (m s−1)	${ \mathcal U }(-10,10)$	${ \mathcal U }(-\inf ,\inf )$
TOI-1235 b Parameters
P (days)	${ \mathcal U }(-\inf ,\inf )$ b	${ \mathcal U }(-\inf ,\inf )$
T0 (BJD −2,457,000)	${ \mathcal U }(-\inf ,\inf )$ b	${ \mathcal U }(-\inf ,\inf )$
$\mathrm{ln}{r}_{{\rm{p}}}$ (${R}_{\oplus }$)	${ \mathcal N }(0.5\cdot \mathrm{ln}(Z)+\mathrm{ln}{R}_{{\rm{s}}},1)$ c	⋯
${r}_{{\rm{p}}}/{R}_{{\rm{s}}}$	⋯	${ \mathcal U }(-\inf ,\inf )$
b	${ \mathcal U }(0,1+{r}_{{\rm{p}},b}/{R}_{{\rm{s}}})$	⋯
$\mathrm{ln}K$ (m s−1)	${ \mathcal U }(-5,5)$	⋯
K (m s−1)	⋯	${ \mathcal U }(-\inf ,\inf )$
e	${ \mathcal B }(0.867,3.03)$ d
ω (rad)	${ \mathcal U }(-\pi ,\pi )$
$e\cos \omega$	⋯	${ \mathcal U }(-1,1)$
$e\sin \omega$	⋯	${ \mathcal U }(-1,1)$
$\sqrt{e}\cos \omega$	${ \mathcal U }(-1,1)$ ⋯
$\sqrt{e}\sin \omega$	${ \mathcal U }(-1,1)$ ⋯

Notes. Gaussian distributions are denoted by ${ \mathcal N }$ and are parameterized by mean and standard deviation values. Uniform distributions are denoted by ${ \mathcal U }$ and bounded by the specified lower and upper limits. Beta distributions are denoted by ${ \mathcal B }$ and are parameterized by the shape parameters α and β.

^aδ is the SPOC-derived dilution factor applied to the TESS light curve. ^bThis prior in the fiducial model reflects that used in the TESS light analysis. However, its resulting posterior is used as an informative prior in the subsequent RV analysis. ^cThe transit depth of TOI-1235.01 reported by the SPOC: Z = 841 ppm. ^dKipping (2013).

Download table as:  ASCIITypeset image

We execute an MCMC to sample the joint posterior probability density function (pdf) of our full set of model parameters using the PyMC3 MCMC package (Salvatier et al. 2016) within exoplanet. The MCMC is initialized with four simultaneous chains, each with 4000 tuning steps and 3000 draws in the final sample. Point estimates of the maximum a posteriori (MAP) values from the marginalized posterior pdf's of the GP hyperparameters are selected to construct the GP predictive distribution, whose mean function is treated as our detrending model of the PDCSAP light curve. This mean detrending function and the detrended light curve are both shown in Figure 2. Similarly, we recover the MAP point estimates of the transit model parameters to construct the transit model shown in the bottom panel of Figure 2. MAP values and uncertainty point estimates from the 16th and 84th percentiles for all model parameters are reported in Table 4.

Table 4.  Point Estimates of the TOI-1235 Model Parameters

Parameter	Fiducial Model Valuesa	EXOFASTv2Model Valuesb
TESS Light-curve Parameters
Baseline flux, f0	1.000024 ± 0.000010	1.000035 ± 0.000018
$\mathrm{ln}{\omega }_{0}$	1.45 ± 0.17	⋯
$\mathrm{ln}{S}_{0}{\omega }_{0}^{4}$	$-{0.16}_{-0.57}^{+0.52}$	⋯
$\mathrm{ln}{s}_{\mathrm{TESS}}^{2}$	0.064 ± 0.006	⋯
TESS limb-darkening coefficient, u1	${0.47}_{-0.24}^{+0.32}$	${0.40}_{-0.26}^{+0.34}$
TESS limb-darkening coefficient, u2	${0.20}_{-0.35}^{+0.38}$	${0.23}_{-0.38}^{+0.37}$
Dilution	⋯	${0.09}_{-0.34}^{+0.21}$
RV Parameters
$\mathrm{ln}\lambda \,{\mathrm{day}}^{-1}$	${4.75}_{-0.10}^{+0.22}$	⋯
$\mathrm{ln}{\rm{\Gamma }}$	$-{0.04}_{-1.9}^{+1.9}$	⋯
$\mathrm{ln}{P}_{\mathrm{rot}}\,{\mathrm{day}}^{-1}$	${3.82}_{-0.11}^{+0.10}$	⋯
$\mathrm{ln}{a}_{\mathrm{HARPS} \mbox{-} {\rm{N}}}\,{{\rm{m}}}^{-1}\,{{\rm{s}}}^{-1}$	${2.94}_{-0.69}^{+0.71}$	⋯
$\mathrm{ln}{a}_{\mathrm{HIRES}}\,{{\rm{m}}}^{-1}\,{{\rm{s}}}^{-1}$	${1.46}_{-0.61}^{+0.76}$	⋯
Jitter, ${s}_{\mathrm{HARPS} \mbox{-} {\rm{N}}}$ (m s−1)	${1.18}_{-0.75}^{+0.64}$	${1.37}_{-0.40}^{+0.46}$
Jitter, ${s}_{\mathrm{HIRES}}$ (m s−1)	${0.11}_{-0.09}^{+0.61}$	${2.47}_{-0.83}^{+1.10}$
Velocity offset, ${\gamma }_{\mathrm{HARPS} \mbox{-} {\rm{N}}}$ (m s−1)	$-{0.81}_{-3.03}^{+2.81}$	${1.39}_{-0.43}^{+0.45}$
Velocity offset, ${\gamma }_{\mathrm{HIRES}}$ (m s−1)	${0.69}_{-2.70}^{+2.50}$	$-{0.33}_{-0.99}^{+0.96}$
TOI-1235 b Parameters
Orbital period, P (days)	${3.444729}_{-0.000028}^{+0.000031}$	${3.444727}_{-0.000039}^{+0.000035}$
Time of midtransit, T0 (BJD –2,457,000)	${1845.51696}_{-0.00098}^{+0.00099}$	${1845.5173}_{-0.0010}^{+0.0008}$
Transit duration D (hr)	${1.84}_{-0.16}^{+0.09}$	${1.94}_{-0.04}^{+0.05}$
Transit depth, Z (ppt)	${0.645}_{-0.044}^{+0.049}$	${0.662}_{-0.038}^{+0.039}$
Scaled semimajor axis, $a/{R}_{{\rm{s}}}$	${13.20}_{-0.40}^{+0.41}$	${13.15}_{-0.32}^{+0.34}$
Planet-to-star radius ratio, ${r}_{{\rm{p}}}/{R}_{{\rm{s}}}$	0.0254 ± 0.0009	0.0257 ± 0.0007
Impact parameter, b	${0.45}_{-0.19}^{+0.21}$	${0.33}_{-0.19}^{+0.15}$
Inclination, i (deg)	${88.1}_{-0.9}^{+0.8}$	${88.6}_{-0.6}^{+0.8}$
$e\cos \omega$	⋯	${0.00}_{-0.03}^{+0.03}$
$e\sin \omega$	⋯	${0.00}_{-0.06}^{+0.04}$
$\sqrt{e}\cos \omega$	${0.07}_{-0.15}^{+0.13}$	⋯
$\sqrt{e}\sin \omega$	$-{0.02}_{-0.23}^{+0.23}$	⋯
Eccentricity, e	<0.15c	<0.16c
Planet radius, rp (${R}_{\oplus }$)	${1.738}_{-0.076}^{+0.087}$	${1.763}_{-0.066}^{+0.071}$
Log RV semiamplitude, $\mathrm{ln}K\,{{\rm{m}}}^{-1}\,{{\rm{s}}}^{-1}$	${1.41}_{-0.13}^{+0.10}$	${1.46}_{-0.13}^{+0.11}$
RV semiamplitude, K (m s−1)	${4.11}_{-0.50}^{+0.43}$	${4.32}_{-0.51}^{+0.50}$
Planet mass, mp (${M}_{\oplus }$)	${6.91}_{-0.85}^{+0.75}$	${7.53}_{-0.89}^{+0.88}$
Bulk density, ${\rho }_{{\rm{p}}}$ (g cm−3)	${7.4}_{-1.3}^{+1.5}$	${7.5}_{-1.2}^{+1.4}$
Surface gravity, gp (m s−2)	${22.6}_{-3.4}^{+3.5}$	${23.7}_{-3.2}^{+3.3}$
Escape velocity, ${v}_{\mathrm{esc}}$ (km s−1)	${22.4}_{-1.5}^{+1.3}$	${23.1}_{-1.1}^{+1.2}$
Semimajor axis, a (au)	${0.03845}_{-0.00040}^{+0.00037}$	${0.03846}_{-0.00032}^{+0.00033}$
Insolation, F (${F}_{\oplus })$	${53.6}_{-4.7}^{+5.3}$	${53.6}_{-4.3}^{+4.2}$
Equilibrium temperature, ${T}_{\mathrm{eq}}$ (K)
Bond albedo = 0.0	754 ± 18	754 ± 18
Bond albedo = 0.3	689 ± 16	689 ± 16
Keplerian Parameters of the 22-day RV Signald
Period (days)	⋯	${21.99}_{-0.32}^{+0.47}$
Reference epoch (analogous to T0) (BJD –2,457,000)	⋯	${1835.34}_{-0.87}^{+0.89}$
Log RV semiamplitude, $\mathrm{ln}K\,{{\rm{m}}}^{-1}\,{{\rm{s}}}^{-1}$	⋯	${1.50}_{-0.14}^{+0.15}$
RV semiamplitude, K (m s−1)	⋯	${4.50}_{-0.57}^{+0.62}$
$e\cos \omega$	⋯	${0.02}_{-0.09}^{+0.11}$
$e\sin \omega$	⋯	${0.09}_{-0.10}^{+0.18}$

Notes.

^aOur fiducial model features sequential modeling of the TESS light curve followed by the RV analysis conditioned on the results of the transit analysis. ^bOur alternative analysis is a global model of the TESS and ground-based light curves, along with the RVs using the EXOFASTv2 software. ^c95% upper limit. ^dThe 22-day RV signal is modeled as an eccentric Keplerian in our EXOFASTv2 model, although we emphasize that here we do not attribute this signal to a second planet.

Download table as:  ASCIITypeset images: 1 2

We continue our fiducial analysis by jointly modeling the HARPS-N and HIRES RV time series. Here we are able to exploit the strong priors on P and T₀ derived from our analysis of the TESS light curve (Section 4.1).

The raw HARPS-N and HIRES RVs are shown in the top row of Figure 5, along with their Bayesian generalized Lomb–Scargle periodogram (BGLS; Mortier et al. 2015). The periodicity induced by TOI-1235 b is distinctly visible at 3.44 days. A preliminary RV analysis indicated that, following the removal of an optimized Keplerian solution for TOI-1235 b, the BGLS revealed a strong periodic signal at 22 days, which is seen at moderately low significance in the BGLS of the raw RVs in Figure 5. This periodicity is close to the first harmonic of the stellar rotation period at ${P}_{\mathrm{rot}}$/2 ≈ 22.3 days. As such, we interpret this signal as likely being produced by active regions on the rotating stellar surface. We note that this feature at ${P}_{\mathrm{rot}}$/2 is similar to the first harmonic of ${P}_{\mathrm{rot}}$ observed on the Sun that has been shown to have either a comparable amount or at times more power than at ${P}_{\mathrm{rot}}$ (Mortier & Collier Cameron 2017; Milbourne et al. 2019). However, we note that simulated RV time series with injected quasi-periodic magnetic activity signals have been shown to produce spurious, and sometimes long-lived, periodogram signals that can masquerade as rotation signatures (Nava et al. 2020). But given that the 22-day signal is nearly identical to the first harmonic of the measured rotation period, we proceed with treating the 22-day signal as stellar activity and opt to simultaneously fit the HARPS-N and HIRES RVs with model components for TOI-1235 b, in the form of a Keplerian orbit, plus a quasi-periodic GP regression model of stellar activity whose covariance kernel as a function of time t takes the form

$$\begin{eqnarray}&&{k}_{{ijs}}={a}_{s}^{2}\left[-\displaystyle \frac{{\left({t}_{i}-{t}_{j}\right)}^{2}}{2{\lambda }^{2}}-{{\rm{\Gamma }}}^{2}{\sin }^{2}\left(\displaystyle \frac{2\pi | {t}_{i}-{t}_{j}| }{{P}_{\mathrm{rot}}}\right)\right].\end{eqnarray} \tag{ 2 }$$

The quasi-periodic kernel is parameterized by four hyperparameters: the covariance amplitude a_s, where s is the index over the two spectrographs; the exponential timescale λ; the coherence Γ; and the periodic timescale, which we initialize to ${P}_{\mathrm{rot}}$/2 because of its apparent periodicity in the BGLS of the raw RVs. Because the temporally correlated signal that we are modeling with a GP likely originates from active regions on the rotating stellar surface, and given the fact that activity signals are inherently chromatic, we consider separate GP activity models for each spectrograph. We also maintain that the covariance hyperparameters {λ, Γ, P_rot} are identical within each spectrograph's GP activity model. We include an additive scalar jitter ${s}_{\mathrm{RV},{\rm{s}}}$ for each spectrograph to account for any excess noise in the activity model and fit for each spectrograph's unique zero-point offset γ_s.

Zoom In Zoom Out Reset image size

Our full RV model therefore consists of 14 free parameters: $\{\mathrm{ln}{a}_{\mathrm{HARPS} \mbox{-} {\rm{N}}}$, $\mathrm{ln}{a}_{\mathrm{HIRES}}$, lnλ, lnΓ, $\mathrm{ln}{P}_{\mathrm{rot}},$ $\mathrm{ln}{s}_{\mathrm{RV},\mathrm{HARPS} \mbox{-} {\rm{N}}},\mathrm{ln}{s}_{\mathrm{RV},\mathrm{HIRES}}$, ${\gamma }_{\mathrm{HARPS} \mbox{-} {\rm{N}}},{\gamma }_{\mathrm{HIRES}},P,{T}_{0},$ $\mathrm{ln}K,h=\sqrt{e}\cos \omega$, $k\,=\sqrt{e}\sin \omega \}$, where K is the RV semiamplitude of TOI-1235 b. The adopted model parameter priors are included in Table 3. We fit the RV data with our full model using the affine-invariant ensemble MCMC sampler emcee (Foreman-Mackey et al. 2013), throughout which we use the george package (Ambikasaran et al. 2014) to evaluate the marginalized likelihood of the GP activity models. MAP point estimates of the model parameters are derived from their respective marginalized posterior pdf's and are reported in Table 4.

The second row in Figure 5 depicts the activity component in our RV model after the MAP Keplerian solution for TOI-1235 b is subtracted from the raw RVs. The residual periodicity close to ${P}_{\mathrm{rot}}$/2 = 22 days becomes clearly visible in the BGLS of the RV activity signal. In our GP activity model, we measure an exponential timescale of λ = 115 ± 20 days, indicating that active regions are relatively stable over a few rotation cycles. According to detailed investigations of periodogram signals in simulated RV time series, the persistence of the maximum RV activity peak at ${P}_{\mathrm{rot}}$/2 is consistent with active region lifetimes on TOI-1235 exceeding ${P}_{\mathrm{rot}}$ (Nava et al. 2020).

In the third row of Figure 5, the BGLS of the TOI-1235 b signal is clearly dominated by the 3.44-day periodicity as expected. We measure an RV semiamplitude of $K={4.11}_{-0.50}^{+0.43}$ m s⁻¹, which is detected at 8.2σ and is clearly visible in the phase-folded RVs in Figure 5. The RV residuals, after removing each spectrograph's mean GP activity model and the MAP Keplerian solution, show no signs of any probable periodicities and have rms values of 1.90 and 1.65 m s⁻¹ for the HARPS-N and HIRES RVs, respectively. We note that these rms values exceed the typical RV measurement uncertainties of 1.2 m s⁻¹ and may be indicative of an incomplete RV model. We reserve an exploration of this prospect until Section 5.3.

To assess the robustness of the parameters derived from our fiducial modeling strategy (Sections 4.1 and 4.2), here we consider an alternative global model using the EXOFASTv2 exoplanet transit plus RV fitting package (Eastman et al. 2013, 2019).

Here we highlight a few notable differences between our fiducial analysis and the global model using EXOFASTv2. In our fiducial model of the TESS PDCSAP light curve, we simultaneously fit the data with a GP detrending model plus a transit model such that the uncertainties in the recovered planetary parameters are marginalized over our uncertainties in the detrending model. Conversely, EXOFASTv2 takes as input a pre-detrended light curve to which the transit model is fit. We construct the detrended light curve to supply to EXOFASTv2 using the mean function of the predictive GP distribution shown in Figure 2. With this method, the uncertainties in the planetary parameters of interest are not marginalized over uncertainties in the detrending model and may consequently be underestimated. Similarly, the RV model in our fiducial analysis considers temporally correlated RV activity signals and models them as a quasi-periodic GP. Conversely, modeling of the prominent 22-day signal in the RVs with EXOFAST requires one to assume a deterministic functional form for the signal in order to construct a more complete RV model. For this purpose, we model the 22-day signal as an eccentric Keplerian within EXOFASTv2. We adopt broad uniform priors on the signal's P and T₀ and adopt identical priors on its semiamplitude, e cos ω, $e\sin \omega$ as are used for TOI-1235 b (Table 3).

The EXOFASTv2 model has the important distinction of evaluating a global model that jointly considers the TESS photometry along with the HARPS-N and HIRES RVs. By virtue of this, the common planet parameters between these data sets (i.e., P, T₀, e, ω) will be self-consistent. In particular, the eccentricity of TOI-1235 b will be jointly constrained by the transit duration, the RV solution, and the stellar density, which is constrained by our priors on the stellar mass and radius (Table 3). The EXOFASTv2 software also explicitly fits for any excess photometric dilution, therefore providing an improved accuracy on the transit depth and hence on the recovered planetary radius. Within EXOFASTv2, the dilution is defined as the fractional flux contribution from neighboring stars (see Section 12 Eastman et al. 2019).

We report the results from our global model in Table 4 and compare the planetary parameters to those derived from our fiducial analysis. All planetary parameters are consistent between our two analysis strategies at <1σ. In particular, in our fiducial and EXOFASTv2 analyses, we measure consistent values for the observables r_p/R_s = 0.0254 ± 0.0009 and 0.0257 ± 0.0007 and $\mathrm{ln}K={1.41}_{-0.13}^{+0.10}$ and ${1.46}_{-0.13}^{+0.11}$. Given the identical stellar parameter priors in each analysis, this consistency directly translates into consistent measures of TOI-1235 b's fundamental planet parameters.

5.1.1. Orbital Separation, Mass, and Radius

Our analysis of the TESS PDCSAP light curve reveals that TOI-1235 b has an orbital period of $P={3.444729}_{-0.000028}^{+0.000031}$ days and a planetary radius of ${r}_{{\rm{p}}}={1.738}_{-0.076}^{+0.087}$ ${R}_{\oplus }$. The corresponding semimajor axis for TOI-1235 b is $a={0.03845}_{-0.00040}^{+0.00037}$ au, where it receives ${53.6}_{-4.7}^{+5.3}$ times Earth's insolation. Assuming uniform heat redistribution and a Bond albedo of zero, TOI-1235 b has an equilibrium temperature of ${T}_{\mathrm{eq}}$ = 754 ± 18 K.

From our RV analysis, we obtain an 8.1σ planetary mass measurement of ${m}_{{\rm{p}}}={6.91}_{-0.85}^{+0.75}$ ${M}_{\oplus }$. Taken together, the mass and radius of TOI-1235 b give a bulk density of ${\rho }_{{\rm{p}}}={7.4}_{-1.3}^{+1.5}$ g cm⁻³. In Figure 6 we add TOI-1235 b to the mass–radius diagram of small M-dwarf planets with ≥3σ mass measurements. Comparing TOI-1235 b's mass and radius to internal structure models of two-layer, fully differentiated planet interiors (Zeng & Sasselov 2013) reveals that the bulk composition of TOI-1235 b is consistent with an Earth-like composition of 33% iron plus 67% silicate rock by mass.

Zoom In Zoom Out Reset image size

Intriguingly, the mass and radius of TOI-1235 b are nearly identical to those of LHS 1140 b, despite LHS 1140 b having a wider 25-day orbit around a mid-M dwarf, thus making it much more temperate than TOI-1235 b (${T}_{\mathrm{eq}}$ = 230 K; Dittmann et al. 2017; Ment et al. 2019). Both planets are situated within the radius valley around low-mass stars (CM19) and have masses that appear to represent the upper limit of terrestrial planet masses in a planetary mass regime where rocky Earth-like planets are inherently rare (i.e., 5–10 ${M}_{\oplus }$; Figure 6). These planets offer unique opportunities to study nature's largest terrestrial planets, whose tectonic and outgassing processes may differ significantly from those on Earth-sized terrestrial planets (Valencia et al. 2007).

With the planetary mass measurement presented herein, TOI-1235 adds to the growing list of small planets transiting M dwarfs with precise RV masses (GJ 3470, Bonfils et al. 2012; GJ 1214, Charbonneau et al. 2009; GJ 1132, Bonfils et al. 2018; K2-3, Damasso et al. 2018; K2-18, Cloutier et al. 2019b; LHS 1140, Ment et al. 2019), which has been rapidly expanding since the launch of TESS (GJ 357, Luque et al. 2019; GJ 1252, Shporer et al. 2020; L98-59, Cloutier et al. 2019a; L168-9, Astudillo-Defru et al. 2020; LTT 3780, Cloutier et al. 2020; Nowak et al. 2020). Notably, TOI-1235 b also directly contributes to the completion of the TESS level-one science requirement of obtaining precise masses for 50 planets smaller than four Earth radii.

5.1.2. Iron and Envelope Mass Fractions

We wish to place self-consistent limits on the iron mass fraction X_Fe and envelope mass fraction X_env of TOI-1235 b. Here the iron mass fraction is defined as the ratio of the total mass of the core and mantle that is composed of iron, with the remainder in magnesium silicate. The envelope mass fraction is then defined as the fraction of the planet's total mass that is in its gaseous envelope. However, it is important to note that these values are degenerate such that we cannot derive a unique solution given only the planet's mass and radius. For example, the bulk composition of TOI-1235 b is consistent with being Earth-like, thus suggesting a small envelope mass fraction,⁵⁰ but one could also imagine a more exotic scenario that is consistent with the planet's mass and radius of a planetary core with X_Fe = 1, surrounded by an extended H/He envelope. In the simplest case, we assume that magnesium silicate and iron are the only major constituents of TOI-1235 b's bulk composition such that X_env = 0. Under this assumption, we derive X_Fe by Monte Carlo sampling the uncorrelated marginalized posterior pdf's of m_p and r_p and use the analytical rock/iron mass–radius relation from Fortney et al. (2007) to recover X_Fe. We find that TOI-1235 b has an iron mass fraction of ${X}_{\mathrm{Fe}}={20}_{-12}^{+15}$% that is <46% at 90% confidence.

To infer the distribution of envelope mass fractions that are consistent with the data, we first impose a physically motivated prior on X_Fe of ${ \mathcal N }(0.33,0.10)$. The relatively narrow width of this Gaussian prior is qualitatively supported by observations of nearby Sun-like stellar metallicities that show that the abundance ratios of Mg/Fe, Si/Fe, and Mg/Si at similar ages and metallicities vary by less than 10%. This indicates a low level of compositional diversity in the refractory building blocks of planets (Bedell et al. 2018). The width of our X_Fe Gaussian prior is chosen in an ad hoc way to approximately reflect this level of chemical diversity. The homogeneity of refractory chemical abundances among Sun-like stars, coupled with their similar condensation temperatures (Lodders 2003), suggests a narrow range in iron mass fractions among close-in terrestrial planets. This assertion is supported by the locus of terrestrial planets with r_p ≲ 1.8 ${R}_{\oplus }$ that are consistent with an Earth-like bulk composition (Figure 6). This concept of similar X_Fe values is particularly compelling for the most massive terrestrial planets (e.g., TOI-1235 b), for which a significant increase in X_Fe by collisional mantle stripping is energetically infeasible owing to the large binding energies of such planets (Marcus et al. 2010).

To proceed with deriving the distribution of TOI-1235 b envelope mass fractions assuming an Earth-like core, we extend the solid two-layer interior structure model to include an H/He envelope with a mean molecular weight equal to that of a solar-metallicity gas (μ = 2.35). Our adopted planetary model is commonly used for sub-Neptune-sized planets (e.g., Rafikov 2006; Lee & Chiang 2015; Ginzburg et al. 2016; Owen & Wu 2017; Gupta & Schlichting 2019). This model features a solid core surrounded by an H/He gaseous envelope that, depending on the planetary parameters, either is fully radiative throughout or may be convective in the deep interior up to the height of the radiative–convective boundary (RCB), above which the atmosphere becomes radiative and isothermal with temperature ${T}_{\mathrm{eq}}$. The latter scenario represents the general case, whereas the former is only invoked when the planetary parameters result in a height of the RCB that is less than the atmospheric pressure scale height at ${T}_{\mathrm{eq}}$. To first order, the height of the RCB above the planetary surface r_RCB, and hence X_env, is determined by $\{{T}_{\mathrm{eq}},{m}_{{\rm{p}}},{r}_{{\rm{p}}},{X}_{\mathrm{Fe}}\}$. Each of ${T}_{\mathrm{eq}}$, m_p, and r_p are directly constrained by our data if we assume a Bond albedo to infer ${T}_{\mathrm{eq}}$. We derive r_RCB and X_env by Monte Carlo sampling X_Fe from its prior, along with the zero-albedo ${T}_{\mathrm{eq}}$, m_p, and r_p from their respective marginalized posterior pdf's. We then rescale each ${T}_{\mathrm{eq}}$ draw by ${(1-{A}_{B})}^{1/4}$, where A_B is the Bond albedo. Super-Earth bond albedos have poor empirical constraints, so we opt to condition a broad uniform prior on A_B of ${ \mathcal U }(0.1,0.8)$ based on the solar system planets. Lastly, although we expect r_RCB to shrink over time as the H/He envelope cools and contracts, this effect on X_env is known to be a weak function of planet age (Owen & Wu 2017) such that we fix the age of TOI-1235 to 5 Gyr in our calculations.

We use a customized version of the EvapMass software (Owen & Campos Estrada 2020) to self-consistently solve for r_RCB and X_env given samples of $\{{T}_{\mathrm{eq}},{m}_{{\rm{p}}},{r}_{{\rm{p}}},{X}_{\mathrm{Fe}}\}$. We attempt to sample these parameters in 10⁵ realizations, although not all parameter combinations are physically capable of producing a self-consistent solution. In practice, our Monte Carlo sampling results in 94,131 successful planetary model realizations (i.e., 94.1% success rate). The resulting distributions of X_Fe and X_env that are consistent with our measurements of TOI-1235 b are shown in Figure 7. We find that 41.1% of successful planet model realizations have fully radiative atmospheres, with the remaining 58.9% being convective in the lower atmosphere. These models produce largely disparate results with radiative atmospheres being favored for increasingly smaller X_Fe and always having X_env ≲ 10⁻³. Conversely, atmospheres with a deep convective region are more extended, thus requiring a more compressed core (i.e., large X_Fe) and larger X_env. Overall we see the positive correlation between X_Fe and X_env because at a fixed m_p the core radius must shrink with increasing X_Fe, which requires the envelope to become extended to match the observed radius. Extending the envelope increases the limits of integration over the atmospheric density profile from the planetary surface to the top of the atmosphere, consequently increasing the envelope mass. With our models, we find that TOI-1235 b has a maximum envelope mass fraction of 2.3%. Marginalizing over all other model parameters and both atmospheric equations of state, we find that X_env must be <0.5% at 90% confidence.

Zoom In Zoom Out Reset image size

Observational studies of the occurrence rate of close-in planets around Sun-like stars have revealed a bimodality in the distribution of planetary radii known as the radius valley (e.g., Fulton et al. 2017; Fulton & Petigura 2018; Mayo et al. 2018). This dearth of planets between 1.7 and 2.0 ${R}_{\oplus }$ around Sun-like stars likely marks the transition between rocky planets and larger planets that host extended gaseous envelopes. Physical models of the emergence of the radius valley from thermally driven atmospheric mass loss (i.e., photoevaporation or core-powered mass loss), and from terrestrial planet formation in a gas-poor environment, make distinct predictions regarding the slope of the radius valley in period–radius space. The slope of the radius valley around Sun-like stars with ${T}_{\mathrm{eff}}$ > 4700 K was measured by Martinez et al. (2019) using the stellar sample from Fulton et al. (2017). The recovered slope was shown to be consistent with model predictions from thermally driven atmospheric mass loss (${r}_{{\rm{p}},\mathrm{valley}}\propto {P}^{-0.15};$ Lopez & Rice 2018). On the other hand, the slope around lower-mass dwarfs with ${T}_{\mathrm{eff}}$ < 4700 K (i.e., mid-K to mid-M dwarfs) was measured by CM19 and was shown to have a flipped sign that instead was consistent with predictions from gas-poor formation (${r}_{{\rm{p}},\mathrm{valley}}\propto {P}^{0.11};$ Lopez & Rice 2018). One interpretation of this is that the dominant mechanism for sculpting the radius valley is stellar mass dependent and that thermally driven mass loss becomes less efficient toward mid- to late M dwarfs, where a new formation pathway of terrestrial planets in a gas-poor environment emerges (CM19). The stellar mass at which this proposed transition occurs is not well resolved by occurrence rate measurements, but it may be addressed by the detailed characterization of individual planets that span the model predictions in period–radius space (e.g., TOI-1235 b).

Differences in the slopes of the radius valley around Sun-like and lower-mass stars naturally carve out a subset of the period–radius space in which the models make opposing predictions for the bulk compositions of planets. This subspace around low-mass stars cooler than 4700 K was quantified by CM19 and is highlighted in Figure 8. At periods less than 23.5 days, planets within the highlighted subspace are expected to be rocky according to models of thermally driven hydrodynamic escape. Conversely, gas-poor formation models predict that those planets should instead be nonrocky with envelope mass fractions of at least a few percent depending on their composition. TOI-1235 b falls within this region of interest and therefore provides direct constraints on the efficiency of the competing physical processes on close-in planets around early M dwarfs.

Zoom In Zoom Out Reset image size

Our transit and RV analyses revealed that TOI-1235 b is a predominantly rocky planet with an iron mass fraction of ${20}_{-12}^{+15}$% and an envelope mass fraction that is <0.5% at 90% confidence. Given its period and radius, this finding is consistent with models of thermally driven mass loss but is inconsistent with the gas-poor formation scenario. Indeed, based on the photoevaporation-driven hydrodynamic escape simulations by Lopez & Fortney (2013), the mass of TOI-1235 b places its insolation flux ($F={53.6}_{-4.7}^{+5.3}$ F_⊕) right at the threshold insolation required for the planet to lose its gaseous envelope: F_threshold = 52 ± 14 F_⊕.⁵¹ These results suggest that thermally driven mass loss continues to be an efficient process for sculpting the radius valley around early M dwarfs like TOI-1235. CM19 suggested that although thermally driven mass loss seems to be prevalent around Sun-like stars, evolution in the structure of the radius valley with stellar mass suggests that this prevalence weakens with decreasing stellar mass and that gas-poor formation may emerge as the dominant mechanism for sculpting the radius valley around early to mid-M dwarfs. Although the stellar mass at which this proposed transition occurs has yet to be resolved, the rocky nature of TOI-1235 b further suggests that the stellar mass at which this transition occurs is likely less than that of TOI-1235 (0.640 ± 0.016 ${M}_{\odot }$).

As an aside, we note that distinguishing between photoevaporation and core-powered mass loss cannot be achieved with the data presented herein. Fortunately, the distinction can be addressed at the planet population level by investigating the radius valley's dependence with time and with stellar mass (Gupta & Schlichting 2020).

Recall that after removing the TOI-1235 b signal from our RV time series, a strong residual periodicity emerges at about 22 days (second row in Figure 5). We initially interpreted this signal as being likely related to rotationally induced stellar activity because of its proximity to the first harmonic of the probable stellar rotation period inferred from ground-based photometric monitoring (${P}_{\mathrm{rot}}$ = 44.7 ± 4.5 days; Figure 3). Although the measurement of ${P}_{\mathrm{rot}}$ makes the 22-day RV signal suggestive of being related to stellar activity, here we conduct a suite of tests that instead favor a planetary origin.

The treatment of the 22-day RV signal as either a quasi-periodic GP in our fiducial model or an eccentric Keplerian in our EXOFASTv2 global model (see Sections 4.2 and 4.3) gives an activity semiamplitude of ≈5 m s⁻¹. This value appears to be at odds with reasonable predictions of the RV signal based on the star's long-term photometric variability from ground-based monitoring (Section 3.2). Using the ${FF}^{\prime}$ model to predict the activity-induced RV variations from photometric variability (Aigrain et al. 2012), we would expect the semiamplitude of the TOI-1235 RV activity signal to be at the level of 1–2 m s⁻¹ instead of the observed value of 5 m s⁻¹ under the single-planet model. However, it is important to note that photometry is not a perfect predictor of RV variations because (i) stellar activity undergoes cycles and there is no guarantee that the level of activity is constant between the epochs of photometric monitoring and the RV observations, (ii) photometry is not sensitive to all spot distributions (Aigrain et al. 2012), and (iii) bright chromospheric plages can produce RV variations with amplitudes similar to those induced by spots of the same size, but with potentially 10 times less flux variations (Dumusque et al. 2014). Therefore, the discrepancy between the observed RV activity variations and the ${FF}^{\prime}$ model predictions is merely suggestive that our RV activity models are overpredicting the amplitude of the RV activity signal, which would then require an additional RV component to model the excess signal in the RV residuals.

Rotationally induced RV signals from active regions arise from the temperature difference between the active regions and the surrounding stellar surface. As such, the active region contrast has an inherent wavelength dependence that increases toward shorter wavelengths such that RV activity signals should be larger at bluer wavelengths (Reiners et al. 2010). We elected to investigate the chromatic dependence of the 22-day RV signal by considering sets of "blue" and "red" RVs from HARPS-N. We rederived the HARPS-N RVs using the same methodology as in our fiducial analysis but focused separately on the spectral orders 0–45 (388–550 nm) and 46–68 (550–689 nm) to derive sets of blue and red RVs, respectively. Each range of orders was selected to achieve a comparable median RV measurement uncertainty in each time series of 1.98 and 2.04 m s⁻¹. We then investigated the chromatic dependence of the probability of the 3.44- and 22-day periodicities in the BGLS. While the 22-day signal was marginally more probable in the red RVs, we found that the planetary signal varied by many more orders of magnitude than the 22-day signal. This behavior is unexpected for a planetary signal that is known to be achromatic. We therefore concluded that this chromatic analysis of our data set is unreliable, and we make no claims regarding the physical origin of the 22-day signal based on its chromatic dependence.

To explicitly test the idea that an additional RV component is required to completely model the data, we considered a two-planet RV model with components for TOI-1235 b, a second Keplerian "c" at 22 days, plus quasi-periodic GP activity models for each spectrograph with an imposed prior on its periodic timescale equal to that of ${P}_{\mathrm{rot}}$: ${ \mathcal N }(44.7,4.5)$ days. We sampled the two-planet model parameter posteriors using an identical method to what was used in our fiducial analysis of the one-planet RV model (Section 4.2). We adopted narrow uniform priors on P_c of ${ \mathcal U }(17,27)$ days and on ${T}_{0,{\rm{c}}}$ of ${ \mathcal U }(1821.5,1848.5)$ BJD –2,457,000. The resulting Keplerian model parameters on the hypothetical planet "c" are reported in Table 5. We find that the hypothetical planet would have a period of ${P}_{{\rm{c}}}={21.8}_{-0.8}^{+0.9}$ days and an RV semiamplitude of ${K}_{{\rm{c}}}={4.2}_{-1.7}^{+1.2}$, which implies a minimum mass of ${m}_{{\rm{p}},{\rm{c}}}\sin i={13.0}_{-5.3}^{+3.8}$ ${M}_{\oplus }$.

Table 5.  Point Estimates of the Hypothetical Planet "c" Keplerian Model Parameters

Parameter	Model Values
Orbital period, Pc (days)	${21.8}_{-0.8}^{+0.9}$
Time of midtransit, ${T}_{0,{\rm{c}}}$ (BJD –2,457,000)	${1835.3}_{-2.1}^{+2.2}$
Log RV semiamplitude, $\mathrm{ln}{K}_{{\rm{c}}}\,{{\rm{m}}}^{-1}\,{{\rm{s}}}^{-1}$	${1.4}_{-0.5}^{+0.3}$
$\sqrt{{e}_{{\rm{c}}}}\cos {\omega }_{{\rm{c}}}$	${0.08}_{-0.4}^{+0.3}$
$\sqrt{{e}_{{\rm{c}}}}\sin {\omega }_{{\rm{c}}}$	${0.12}_{-0.47}^{+0.37}$
Derived Parameters
RV semiamplitude, Kc (m s−1)	${4.2}_{-1.7}^{+1.2}$
Minimum planet mass, ${m}_{{\rm{p}},{\rm{c}}}\sin i$ (${M}_{\oplus }$)	${13.0}_{-5.3}^{+3.8}$
Semimajor axis, ac (au)	${0.1319}_{-0.0043}^{+0.0046}$
Insolation, Fc (${F}_{\oplus })$	${4.6}_{-0.5}^{+0.6}$
Equilibrium temperature, ${T}_{\mathrm{eq},c}$ (K)
Bond albedo = 0.0	407 ± 12
Bond albedo = 0.3	373 ± 11

Note. Note that we do not conclude that the hypothetical planet "c" presented in this table is a bona fide planet.

Download table as:  ASCIITypeset image

We now have one- and two-planet RV models of the HARPS-N plus HIRES RVs that both include a GP activity component whose periodic timescales are constrained to be close to ${P}_{\mathrm{rot}}$/2 and ${P}_{\mathrm{rot}}$, respectively. Therefore, we can use our models to conduct a model comparison to assess the favorability of one model over the other. We used the marginalized posterior pdf's from each model's MCMC results to estimate their Bayesian model evidences ${ \mathcal Z }$ using the estimator from Perrakis et al. (2014). We estimate model evidences of $\mathrm{ln}{{ \mathcal Z }}_{1}=-110.0$ and $\mathrm{ln}{{ \mathcal Z }}_{2}=-91.0$, which gives a model evidence ratio of ${{ \mathcal Z }}_{2}/{{ \mathcal Z }}_{1}={10}^{8}$. This result strongly favors the two-planet model, although we caution that Bayesian model evidences are notoriously difficult to accurately calculate and their interpretation is dependent on the assumed model parameter priors (Nelson et al. 2020). Alternatively, we also compute the Bayesian information criterion (BIC) and the Akaike information criterion (AIC) to perform model comparisons that are independent of the model priors. We measure BIC₁ = 225.0 and BIC₂ = 205.2 such that the two-planet model is again strongly favored since ${\rm{\Delta }}{\mathrm{BIC}}_{12}\equiv {\mathrm{BIC}}_{1}-{\mathrm{BIC}}_{2}=19.8\gt 10$. This is further supported by AIC₁ = 202.1 and AIC₂ = 174.1, whereby the two-planet model remains strongly favored as ${\rm{\Delta }}{\mathrm{AIC}}_{12}\equiv {\mathrm{AIC}}_{1}-{\mathrm{AIC}}_{2}\,=28.0\gt 10$.

Encouraged by the prospect of a second planet orbiting TOI-1235, we used its measured orbital period ${P}_{{\rm{c}}}={21.8}_{-0.8}^{+0.9}$ days and its time of inferior conjunction ${T}_{0,{\rm{c}}}=2,458,{835.3}_{-2.1}^{+2.2}$ BJD to search for transit-like events in the TESS PDCSAP and archival MEarth-North light curves. With TESS we conducted the search for periodic transit-like signals close to P_c using the implementation of the Box Least Squares (BLS) algorithm (Kovács et al. 2002) in Cloutier (2019). We conducted a complementary BLS search on the full MEarth-North light curve following the methods outlined in Ment et al. (2019). We do not find any significant transit-like signals other than those associated with TOI-1235 b. Therefore, if the 22-day signal is truly a planet, then it is unlikely to be transiting. This result is perhaps unsurprising given that if the hypothetical planet "c" is coplanar with TOI-1235 b at 881, then "c" would not have a transiting configuration at its separation of ${a}_{{\rm{c}}}/{R}_{{\rm{s}}}={45.1}_{-1.9}^{+2.0}$.

We emphasize that while the aforementioned lines of evidence are suggestive of a second, nontransiting planet around TOI-1235, the data presented herein are not sufficient to firmly distinguish between planetary and stellar activity origins of the 22-day RV signal. Ongoing spectroscopic monitoring of TOI-1235 over many rotation cycles may help to solve this ambiguity by testing for temporal correlations of the signal's amplitude over the star's evolving magnetic activity cycle. A more secure detection of the stellar rotation period from continued photometric monitoring would also be beneficial.

Following the announcement of the TOI-1235.01 level-one planet candidate in 2019 October, multiple precision RV instrument teams began pursing its mass characterization through TFOP. This study has presented the subset of those efforts from HARPS-N and HIRES, but we acknowledge that another collaboration has also submitted a paper presenting their own RV time series and analysis (Bluhm et al. 2020). Although the submissions of these complementary studies were coordinated between the two groups, their respective data, analyses, and write-ups were intentionally conducted independently.