# TOI-1338: TESS' First Transiting Circumbinary Planet

Veselin B. Kostov, Jerome A. Orosz, Adina D. Feinstein, William F. Welsh, Wolf Cukier, Nader Haghighipour, Billy Quarles, David V. Martin, Benjamin T. Montet, Guillermo Torres, et al.

► To cite this version:

Veselin B. Kostov, Jerome A. Orosz, Adina D. Feinstein, William F. Welsh, Wolf Cukier, et al. TOI-1338: TESS' First Transiting Circumbinary Planet. The Astronomical Journal, 2020, 159, 10.3847/1538-3881/ab8a48. insu-03667098

HAL Id: insu-03667098

https://insu.hal.science/insu-03667098v1

Submitted on 4 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

TOI-1338: TESS' First Transiting Circumbinary Planet

Veselin B. Kostov¹,² ID, Jerome A. Orosz³ ID, Adina D. Feinstein⁴,<fnref n="42" /> ID, William F. Welsh³ ID, Wolf Cukier⁵,
Nader Haghighipour⁶ ID, Billy Quarles⁷ ID, David V. Martin⁴,<fnref n="43" /> ID, Benjamin T. Montet⁴,⁸ ID, Guillermo Torres⁹ ID,
Amaury H. M. J. Triaud¹⁰ ID, Thomas Barclay¹,¹¹ ID, Patricia Boyd¹ ID, Cesar Briceno¹² ID, Andrew Collier Cameron¹³,
Alexandre C. M. Correia¹⁴, Emily A. Gilbert¹,⁴,¹⁵,¹⁶ ID, Samuel Gill¹⁷ ID, Michaël Gillon¹⁸ ID, Jacob Haqq-Misra¹⁹ ID,
Coel Hellier²⁰ ID, Courtney Dressing²¹ ID, Daniel C. Fabrycky⁴ ID, Gabor Furesz²², Jon M. Jenkins²³ ID, Stephen R. Kane²⁴ ID,
Ravi Kopparapu¹ ID, Vedad Kunovac Hodžić¹⁰ ID, David W. Latham⁹ ID, Nicholas Law²⁵ ID, Alan M. Levine²², Gongjie Li⁷ ID,
Chris Lintott²⁶ ID, Jack J. Lissauer²³ ID, Andrew W. Mann²⁷ ID, Tsevi Mazeh²⁸, Rosemary Mardling²⁹, Pierre F. L. Maxted²⁰ ID,
Nora Eisner²⁶ ID, Francesco Pepe³⁰, Joshua Pepper²²,³¹ ID, Don Pollacco¹⁷ ID, Samuel N. Quinn⁹ ID, Elisa V. Quintana¹,
Jason F. Rowe³² ID, George Ricker²², Mark E. Rose²³, S. Seager²²,³³,³⁴ ID, Alexandre Santerne³⁵ ID, Damien Ségransan³⁰ ID,
Donald R. Short³ ID, Jeffrey C. Smith³⁶ ID, Matthew R. Standing¹⁰ ID, Andrei Tokovinin¹² ID, Trifon Trifonov³⁷ ID, Oliver Turner³⁰,
Joseph D. Twicken²,²³ ID, Stéphane Udry³⁰ ID, Roland Vanderspek²² ID, Joshua N. Winn³⁸ ID, Eric T. Wolf³⁹ ID, Carl Ziegler⁴⁰ ID,
Peter Ansorge⁴¹, Frank Barnet⁴¹, Joel Bergeron⁴¹, Marc Huten⁴¹, Giuseppe Pappa⁴¹, and Timo van der Straeten⁴¹

¹ NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA; veselin.b.kostov@nasa.gov

² SETI Institute, 189 Bernardo Avenue, Suite 200, Mountain View, CA 94043, USA

³ Department of Astronomy, San Diego State University, 5500 Campanile Drive, San Diego, CA 92182, USA

⁴ Department of Astronomy and Astrophysics, University of Chicago, 5640 S. Ellis Avenue, Chicago, IL 60637, USA

⁵ Scarsdale High School, 1057 Post Road, Scarsdale, NY 10583, USA

⁶ Institute for Astronomy, University of Hawaii-Manoa, Honolulu, HI 96822, USA

⁷ Center for Relativistic Astrophysics, School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA

⁸ School of Physics, University of New South Wales, Sydney NSW 2052, Australia

⁹ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

¹⁰ School of Physics & Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK

¹¹ University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250, USA

¹² Cerro Tololo Inter-American Observatory, CTIO/AURA Inc., Casilla 603, La Serena, Chile

¹³ Centre for Exoplanet Science, SUPA, School of Physics and Astronomy, University of St Andrews, St Andrews KY16 9SS, UK

¹⁴ CFisUC, Department of Physics, University of Coimbra, 3004-516 Coimbra, Portugal

¹⁵ The Adler Planetarium, 1300 South Lakeshore Drive, Chicago, IL 60605, USA

¹⁶ GSFС Sellers Exoplanet Environments Collaboration, USA

¹⁷ Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK

¹⁸ Astrobiology Research Unit, Université de Liège, 19C Allée du 6 Août, B-4000 Liège, Belgium

¹⁹ Blue Marble Space Institute of Science, Seattle, Washington, USA

²⁰ Astrophysics Group, Keele University, Staffordshire, ST5 5BG, UK

²¹ Department of Astronomy, University of California at Berkeley, Berkeley, CA 94720, USA

²² Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

²³ NASA Ames Research Center, Moffett Field, CA, 94035, USA

²⁴ Department of Earth and Planetary Sciences, University of California, Riverside, CA 92521, USA

²⁵ Department of Physics and Astronomy, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599-3255, USA

²⁶ Oxford University, UK

²⁷ Department of Physics and Astronomy, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA

²⁸ Department of Astronomy and Astrophysics, Tel Aviv University, 69978 Tel Aviv, Israel

²⁹ School of Physics & Astronomy, Monash University, Victoria, 3800, Australia

³⁰ Observatoire Astronomique de l'Universit de Genève, 51 chemin des maillettes, CH-1290 Sauverny, Switzerland

³¹ Department of Physics, Lehigh University, 16 Memorial Drive East, Bethlehem, PA 18015, USA

³² Bishops University, 2600 College Street, Sherbrooke, QC J1M 1Z7, Canada

³³ Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

³⁴ Department of Aeronautics and Astronautics, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA

³⁵ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

³⁶ SETI Institute/NASA Ames Research Center, 189 Bernardo Avenue, Suite 200, Mountain View, CA, 94043, US

³⁷ Max Planck Institut fuer Astronomie, Heidelberg, Germany

³⁸ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

³⁹ Department of Atmospheric and Oceanic Sciences, Laboratory for Atmospheric and Space Physics, University of Colorado Boulder, Boulder, CO, USA

⁴⁰ Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, Ontario M5S 3H4, Canada

⁴¹ Planet Hunters TESS

Received 2019 October 30; revised 2020 April 12; accepted 2020 April 15; published 2020 May 7

Abstract

We report the detection of the first circumbinary planet (CBP) found by Transiting Exoplanet Survey Satellite (TESS). The target, a known eclipsing binary, was observed in sectors 1 through 12 at 30 minute cadence and in sectors 4 through 12 at 2 minute cadence. It consists of two stars with masses of 1.1 $M_\odot$ and 0.3 $M_\odot$ on a slightly eccentric (0.16), 14.6 day orbit, producing prominent primary eclipses and shallow secondary eclipses. The planet has a radius of ~6.9 $R_\odot$ and was observed to make three transits across the primary star of roughly equal depths

<footnote n="42">⁴² NSF Graduate Research Fellow.</footnote>

<footnote n="43">⁴³ Fellow of the Swiss National Science Foundation.</footnote>
(~0.2%) but different durations—a common signature of transiting CBPs. Its orbit is nearly circular ($e \approx 0.09$) with an orbital period of 95.2 days. The orbital planes of the binary and the planet are aligned to within ~1°. To obtain a complete solution for the system, we combined the TESS photometry with existing ground-based radial-velocity observations in a numerical photometric-dynamical model. The system demonstrates the discovery potential of TESS for CBPs and provides further understanding of the formation and evolution of planets orbiting close binary stars.

*Unified Astronomy Thesaurus concepts: Exoplanet astronomy (486); Eclipsing binary stars (444)*

# 1. Introduction

One of the most exciting breakthroughs from the Kepler mission was the discovery of circumbinary planets (CBPs). Four years of continuous observations of several thousand eclipsing binary stars (EBs; Prša et al. 2011; Slawson et al. 2011; Kirk et al. 2016) led to the discovery of 13 transiting CBPs orbiting 11 Kepler EBs (Doyle et al. 2011; Orosz et al. 2012a, 2012b, 2019; Welsh et al. 2012, 2015; Kostov et al. 2013, 2014, 2016; Schwamb et al. 2013; Socia et al. 2020). These discoveries spanned a number of firsts—e.g., the first transiting CBP, the first CBP in the Habitable Zone (HZ), the first CBP in a quadruple star system, and the first transiting multi-planet CBP system. In addition to opening a new chapter in the studies of extrasolar planets, Kepler's CBPs have confirmed theoretical predictions that planet formation in circumbinary configurations is a robust process and suggest that many such planetary systems must exist (e.g., Pierens & Nelson 2013; Kley & Haghighipour 2015). Indeed, recent studies argue that the occurrence rate of giant, Kepler-like CBPs is comparable to that of giant planets in single-star systems (~10%; Armstrong et al. 2014; Martin & Triaud 2014; Li et al. 2016; Martin et al. 2019).

As exciting as the CBPs discovered from Kepler are, however, the present sample is small, likely hindered by observational biases, and leaves a vast gap in our understanding of this new class of worlds. This is not unlike the state of exoplanet science 20 yr ago, when only a handful of hot Jupiter exoplanets were known. Pressing questions remain regarding the formation and migration efficiency of CBPs, their orbital architectures and occurrence rates, and the formation, evolution and population characteristics of their host binary stars. Addressing these questions requires more CBP discoveries— which require continuous observations of a large number of EBs for prolonged periods of time. NASA's Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) will assist CBP discovery by observing roughly half a million EBs continuously for timespans between one month and one year for the nominal mission (Sullivan et al. 2015). This motivated us to continue our search for transiting CBPs by examining the light curves of EBs observed by TESS.

Here, we report the discovery of the first CBP from TESS—TOI-1338, a Saturn-sized planet orbiting the known eclipsing binary star EBLM J0608-59⁴⁴ approximately every 95 days. At the time of this writing, this is the longest-period confirmed planet discovered by TESS. Below, we present the details of our discovery and discuss some of the characteristics of this newly found CBP that allow us to place its discovery in a broader context.

This paper is organized as follows. In Section 2, we briefly outline the discovery of the system and describe the TESS data for the target star and the detection of the CBP transits. In Section 3, we present the complementary observations and data analysis. Section 4 outlines the photometric-dynamical analysis of the system, and Section 5 presents a discussion of the results. We draw our conclusions in Section 6.

# 2. Discovery

## 2.1. TESS Mission

The primary goal of the TESS mission is to identify transiting planets around nearby bright stars that are amenable to follow-up characterization. TESS will observe about 85% of the sky during its two-year primary mission (Ricker et al. 2015), using four cameras that provide a 24° × 96° field of view (FOV); a sector is a ~month-long observation of a single FOV. Most of the stars in the full-frame images (FFIs) will be observed at 30 minute cadence, and 200,000 pre-selected stars (spread over the whole sky) will be observed at 2 minute cadence. TESS observers 13 sectors per hemisphere, per year for at least ≈27 days; where the sectors overlap near the ecliptic poles, in the two continuous viewing zones (CVZs), TESS observes for up to ≈350 days. The CVZs are especially valuable places to search for exoplanets, as the longer baseline enables the detection of smaller and/or longer-period planets— like the CBP presented here—and also overlaps with the James Webb Space Telescope CVZ (Ricker et al. 2015). Half way through its primary mission, TESS has already discovered a number of confirmed planets and identified more than a thousand planet candidates (e.g., Huang et al. 2018; Kostov et al. 2019; Vanderspek et al. 2019, and references therein), vetted by the TESS Data Validation initiative (Twicken et al. 2018; Li et al. 2019; N. Guerrero et al. 2020, in preparation).

## 2.2. Discovery of the Host Eclipsing Binary

TOI-1338 was identified as an eclipsing binary in 2009 as a part of the Eclipsing Binary Low Mass (EBLM; Triaud et al. 2013) project, a survey constructed using the false-positives of the WASP survey for transiting hot-Jupiters (Pollacco et al. 2006; Collier Cameron et al. 2007; Triaud 2011). As the observed eclipse depth is comparable to that of an inflated hot Jupiter, at the time, it was not possible to distinguish between an eclipsing late-type M dwarf and a transiting hot Jupiter from the photometric signature alone. Follow-up observations with the CORALIE<fnref n="45" /> high-resolution spectrograph obtained near quadrature revealed a semi-amplitude of 21.6 km s⁻¹, confirming that the target is a low-mass eclipsing binary (Triaud et al. 2017a). As part of the

<footnote n="44">⁴⁴ The target also has the designations TIC 260128333, TYC 8533-00950-1, and Gaia DR2 5494443978353833088. Its R.A. and decl. are 06:08:31.97 and −59:32:28.08, respectively. It has TESS magnitude $T = 11.45 \pm 0.02$ mag and $V = 11.72 \pm 0.02$ mag.</footnote>

<footnote n="45">⁴⁵ CORALIE is a fiber-fed échelle spectrograph, mounted on the Swiss Euler 1.2 m telescope at La Silla, Chile. It has a resolving power of $R = 55,000$ and achieves long-term stability through a thermally stabilized housing and nightly calibrations with respect to a Thorium–Argon reference spectrum (Lovis & Pepe 2007).</footnote>
EBLM project to improve the M-dwarf mass-radius-temperature-luminosity relation, TESS observations at 2 minute cadence were obtained for the target under a Cycle 1 Guest Observer program for Sectors 4–12 (G011278—PI O. Turner).

## 2.3. Detection of the CBP

Based on the 2 minute cadence observations from TESS, the target was flagged as an eclipsing binary on the Planet Hunters TESS platform<sup>46</sup> by user Pappa on subject 31326051. This is a citizen science project where, in addition to primarily flagging transit-like features, volunteers may flag targets as various phenomena including eclipsing binaries, variable stars, etc. Planet Hunters has had a successful contribution to the field of CBPs through the independent discovery of Kepler-64 (also known as Planet Hunters-1; Kostov et al. 2013; Schwamb et al. 2013).

As part of our search for CBPs around eclipsing binaries from TESS, we have been performing a visual inspection of targets tagged as potential eclipsing binaries (with a hashtag “eclipsingbinary”) on the Planet Hunters TESS Talk. We note that TOI-1338 was also listed on exo.MAST<sup>47</sup> as two Threshold Crossing Events with a period of ~14.61 days.

During the examination of the light curves on Planet Hunters TESS Talk, one of us (W.C.) noticed that in one of the partial light curves of TOI-1338, there was both a prominent primary eclipse and an additional feature. The latter was of similar depth to the secondary eclipse but not at the expected time according to exo.MAST. To confirm that the feature is real and not a false positive caused by instrumental artifacts, we next extracted the 2 minute cadence light curve using the Lightkurve software package (Lightkurve Collaboration et al. 2018) and confirmed that there are indeed two genuine transit-like features in Sectors 6 and 10 that are not associated with the secondary eclipses, separated by ≈95 days, and with different durations (~0.3 and ~0.6 days, respectively). Further analysis of the 30 minute cadence data, extracted with eleanor (Feinstein et al. 2019), revealed the presence of a third transit in Sector 3, ≈93 days before the transit in Sector 6, with a duration of ~0.4 days, further strengthening the CBP interpretation. Overall, the three transits exhibit the trademark “smoking gun” signatures of transiting CBPs (Welsh & Orosz 2018) where (i) the transit durations vary depending on the orbital phase of the host EB such that transits across the primary star occurring near primary eclipses have shorter durations than transits occurring near secondary eclipses, and (ii) the transit times vary significantly from a linear ephemeris where specifically in this case the interval between the first and second transits (~93 days) is significantly different than the interval between the second and third transits (~95 days).

We note that it is highly unlikely for the three CBP transits to be a false-positive scenario due to, for example, an unresolved eclipsing binary star. First, because the three transits have different durations, there would be only two plausible scenarios. One scenario, S1, involves two unresolved eclipsing binaries where one (hereafter EB1) would produce the first and third transits as primary and secondary eclipses, and the other (hereafter EB2) produces only the second transit as either a primary or a secondary eclipse. The other scenario, S2, requires an unresolved triple star system consisting of EB1 (producing

the first and the third transits as a primary and secondary eclipse) and a long-period third star (hereafter EB3) producing the second transit as an eclipse across either the primary or the secondary star of EB1. For circular orbits, EB1 needs to have an orbital period of about 400 days (twice the time between the first and third transit) regardless of the scenario, EB2 for scenario S1 should have an orbital period greater than about 200 days (so that it would not produce a second eclipse during the TESS observations), and EB3 would need to have a dynamically stable orbit around the ~400 day EB1. Second, because the three CBP transits have a depth of ~0.2%, such background EBs cannot be fainter than $T \sim 18$ mag (i.e., $\Delta T \sim 7$ mag difference compared to TOI-1338). We note that spectroscopy observations of the target do not show any signs of a second or third EB, albeit not at a contrast level of 7 s magnitudes difference. Below, we explore this further using rough approximations to estimate the order of magnitude of the probability.

While the TESS Input Catalog indicates that there are 98 contaminating sources for TOI-1338, Gaia shows that there are only six sources within $\Delta T \sim 7$ mag<sup>48</sup> inside the entire $13 \times 13$ TESS pixel array of the target, although none of them is inside the TESS aperture of the target. Assuming that these six sources are representative of the field of view, the density of the sources within $\Delta T \sim 7$ mag of the target is then ~0.036 sources pixel<sup>-1</sup>, i.e., ~$10^{-4}$ sources/sq. arcsec. The contrast sensitivity of Gaia DR2 for $\Delta G \sim 7$ mag (and thus $\Delta T \sim 7$ mag) is ~3" (Brandeker & Cataldi 2019). Thus, the probability of having one $\Delta T \sim 7$ mag source unresolved by Gaia within 9 sq. arcsec of the target is $P_{unresolved} \sim 10^{-3}$.

Using the results of Raghavan et al. (2010), we estimate that the probability of EB1 having an orbital period of 400 days is $P_{EB1,per} \sim 20\%$; the probability of EB2 having an orbital period of 200 days is $P_{EB2,per} \sim 85\%$. Using the results of Tokovinin (2014), we estimate the probability of a triple star to be $P_{triple} \sim 10\%$. The probability that EB1, EB2, and EB3 are eclipsing is roughly $(R_1 + R_2)/a$. Assuming similar stars (because of equal depth eclipses) of solar radius and mass (because larger stars are more rare and smaller stars would be too faint to produce the required contamination) for both scenarios, the corresponding probabilities are $P_{EB1,ecl} \sim (2 \times 0.0046)/1.34 \sim 0.006$ and $P_{EB2,ecl} \sim (2 \times 0.0046)/0.84 \sim 0.011$. For scenario S2, we used Rebound's IAS15 integrator (Rein & Spiegel 2015) to test that the orbital period of EB3 would need to be greater than ~2000 days to be dynamically stable. Thus, the probability that EB3 is eclipsing EB1 is $P_{triple,ecl} \sim (2 \times 0.0046)/4.48 \sim 0.002$.

We also note that the orbital phases of EB1, EB2, and EB3 would need to be such that in ~325 days (duration of TESS observations) EB1 would produce one primary and one secondary eclipse, and EB2 and EB3 would produce a single eclipse—within a window of a few days for its duration to agree with the CBP model. This introduces additional constraints. Namely, for orbital periods of 400 days, 200 days, and 2000 days for EB1, EB2, and EB3, respectively, and assuming said window is ~2 days, the corresponding probabilities are $P_{EB1,phase} \sim 80\%$ (i.e., ~325/400), $P_{EB2,phase} \sim 1\%$ (i.e., ~2/200), and $P_{EB3,phase} \sim 0.1\%$ (i.e., ~2/2000).

<footnote n="46">46 https://www.zooniverse.org/projects/nora-dot-eisner/planet-hunters-tess</footnote>

<footnote n="47">47 https://exo.mast.stsci.edu</footnote>

<footnote n="48">48 Gaia magnitude is similar to TESS magnitude, i.e., $G = T + 0.43$ (Stassun et al. 2019).</footnote>
Putting all of this together, the combined probabilities for scenarios S1 and S2 are as follows:

$$ P_{S1} = P_{\text{unresolved,EB1}} \times P_{\text{unresolved,EB2}} \times P_{\text{EB1,per}} \times P_{\text{EB2,per}} \\ \times P_{\text{EB1,ecl}} \times P_{\text{EB2,ecl}} \times P_{\text{EB1,phase}} \times P_{\text{EB2,phase}} \quad (1) $$

$$ = 10^{-3} \times 10^{-3} \times 0.2 \times 0.85 \times 0.006 \\ \times 0.011 \times 0.8 \times 0.01 \sim 10^{-13} \quad (2) $$

$$ P_{S2} = P_{\text{unresolved}} \times P_{\text{triple}} \times P_{\text{EB1,per}} \times P_{\text{EB1,ecl}} \\ \times P_{\text{EB3,ecl}} \times P_{\text{EB1,phase}} \times P_{\text{EB3,phase}} \quad (3) $$

$$ = 10^{-3} \times 0.1 \times 0.2 \times 0.006 \times 0.002 \times 0.8 \\ \times 0.001 \sim 2 \times 10^{-13}. \quad (4) $$

To convert these numbers to a false-positive probability (FPP), we compare them to the probability that TOI-1338 b is a CBP. Specifically, from the Kepler data set, the probability that any given star has a transiting CBP is ~10/200,000 = 5 × 10⁻⁵. Thus, the CBP hypothesis is ~5 × 10⁻⁵/10⁻¹³ = 5 × 10⁻⁸ times more likely than the false-positive hypothesis. The FPP is therefore about 10⁻⁸.

These false-positive scenarios can be further argued against based on the durations and depths of the three transits. Specifically, as the depths of the first and third transits are similar, then the two stars of EB1 should have comparable sizes. As discussed above, assuming sunlike stars, the semimajor axis of EB1 would be ~1.34 au, and the corresponding duration of the primary eclipse of EB1 (for circular orbit, R₁ = R₂ = R₃, and impact parameter b = 0) would be ~0.9 days, i.e., nearly twice the duration of the observed transits, thus further ruling out this specific scenario. While an eccentric orbit of EB1 might alleviate this tension to some degree, this would require special orbital elements—in addition to the requirement imposed by the special orbital phase as discussed above. Thus, overall, we consider false-positive scenarios S1 and S2 to be highly unlikely.

## 2.4. TESS Light Curve

TESS telemeters data in two modes: postage stamps, i.e., small regions, around roughly 20,000 stars at 2 minute cadence every sector<fnref n="49" /> as well as the FFIs, which contain about a million stars each (brighter than T = 15 mag), at 30 minute cadence. TOI-1338 was observed by TESS in 30 minute cadence in Sectors 1–12, and in 2 minute cadence in Sectors 4–12.

The TESS light curve of TOI-1338 is shown in Figure 1, where the 30 minute cadence data extracted with eleanor using aperture photometry are shown for Sectors 1–12 and 2 minute cadence SAPFLUX measurements from the standard processing provided by the TESS mission are shown for Sectors 4–12. The nominal units of the observation times are days from BJD 2,457,000. The eleanor software package performs background subtraction, aperture photometry, and detrending for a given source on the FFIs. It also provides the opportunity to use a custom aperture and, if desired, can use

models of the point-spread function (PSF) to extract the light curves.

Each FFI delivered by the TESS project is barycentric corrected. However, there is only a single correction applied to each CCD, each of which covers more than 100 square degrees of the sky. As a result, the barycentric-corrected times in the raw FFI data can be discrepant by up to a minute. The eleanor software corrects for this potential offset, removing the barycentric correction and applying a more accurate value given the actual position of the target in question. We verified that the eleanor timestamps were accurate, comparing them to the midpoint of the 15 2 minute images that make up a single FFI. We find that the two are consistent at the 2σ level, sufficiently precise for the photodynamical modeling we employ in Section 4.

In addition to the prominent stellar eclipses, the light curve of TOI-1338 contains several small, transit-like events that required further scrutiny. Specifically, we noticed four events near days 1390, 1391, 1403, and 1404 (Sector 3), one event near day 1484 (Sector 6), and an event near day 1579 (Sector 10). Using the PSF fitting built into eleanor, we showed that the events near days 1391, 1484, and 1579 are astrophysical in origin—namely, the three transits of the CBP—while the remaining events are artifacts caused by pointing “jitter” (mainly before reaction wheel “momentum dumps”).

The 30 minute cadence light curve from Sector 3 was particularly difficult to extract because the target fell near an edge of the detector as shown in Figure 2 (the target was well away from the detector edges in the remaining Sectors). Additionally, the light-curve extraction is especially prone to systematic errors, as occasional spacecraft pointing jitter during Sector 3 moved some of the light from the target off of the detector, thereby producing transit-like events in the observed flux. As discussed below, the PSF of the target changes during some of the events listed above, which helps us rule out an astrophysical origin.

The PSF can be modeled in eleanor using either a two-dimensional Gaussian or a Moffat profile, and in this case, both models perform equally well. Briefly, the analysis proceeded as follows. At each cadence, we fit parameters that describe the shape of the stellar PSF, assuming all stars in a 13 × 13 pixel region share the same PSF. We then optimized the flux of each star, the shared PSF parameters, and a single background level across this region. Based on this analysis, we found that during the events near days 1390, 1403, and 1404 in Sector 3, the shape of the PSF changed, affecting how many of the star falls in the optimal aperture and, therefore, how much flux is observed. It is likely that the pointing became “looser” during the times of these events, and the PSF-fitting package interpreted the resulting images as an increase in the size of the PSF. The final pointing algorithm was implemented on board the spacecraft after Sector 4, and TESS experienced sporadic pointing errors more frequently and of higher amplitude prior to that. Owing to an unfortunate set of coincidences, the spurious transit events in Sector 3 happen to have similar depths and durations as the real CBP transits and secondary eclipses. Given the PSF change, we then built a linear model from the out-of-eclipse data that predicts the flux of the star at every cadence from the PSF parameters alone. Figure 3 shows part of the Sector 3 light curve from eleanor using aperture photometry (top curve) and the model light curve predicted from changes in the PSF (bottom curve). We

<footnote n="49">⁴⁹ For sectors 1–3, slightly less than 16,000 targets received 2 minute cadence observations. The number of 2 minute cadence targets was increased to 20,000 after predicted compressibility was demonstrated in flight and the SPOC demonstrated in a ground segment test that 20,000 targets could be handled in the pipeline (Jenkins et al. 2016).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Upper panel: Sectors 1–12 of 30 minute cadence TESS *eleanor* PSF-extracted light curve (black). Some of the eclipses were missed due to data gaps. The CBP transits are highlighted in green. Middle panel: Sectors 4–12 of 2 minute cadence SAPFLUX light curve (gray). Lower panels: 5 day sections of the light curve centered on the three CBP transits (highlighted in green). The lower middle and lower right panels also show the primary and secondary eclipses near days 1483 and 1577, respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Nominal *eleanor* apertures for Sectors 3, 6, and 10. The target was on the edge of the detector in Sector 3, as seen in the left panel. Hence, in this case, a custom aperture is needed.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Aperture photometry for Sector 3 data (upper curve) with the flux model predicted from changes in the shape of the PSF at each cadence inferred by *eleanor* PSF modeling (lower curve). We expect that spurious events caused by changes in the PSF to appear in both curves, whereas astrophysical events, which occur independently of the instrumental PSF, will only appear in the upper curve. Thus, the apparent events near days 1390, 1403, and 1404 are caused by changes in the PSF shape that occur as parts of the star fall off and on the detector owing to pointing jitter. The purple- and green-shaded regions correspond to the locations of secondary eclipses and a transit of the planet, respectively. These events are apparent only in the upper curve, which is a strong indication of their validity. Finally, we note that the deep feature near day 1395 is a primary eclipse.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** A section of the TESS *eleanor* per-pixel light curve for Sector 3 near days 1402 and 1405. The panels with black light curves represent the core pixels used for light-curve extraction. The two events highlighted in green appear as transit-like in some of the core pixels and anti-transit-like in others. This demonstrates that these two events are data artifacts (see Figure 3).</figcaption></figure>

found that the PSF-model predicts transit-like features near days 1390, 1402, and 1405. However, this model does *not* predict transit-like events at the time of a secondary eclipse (days 1387 and 1402) or near day 1391 (CBP transit). This is further demonstrated in Figure 4, showing a section of the per-pixel light curve for the two events near days 1402 and 1405—the events are anticorrelated among some of the core pixels such that they appear as transit-like in some pixels and anti-transit-like in others, indicating data artifacts. We thus confirmed the reality of the day 1391 transit-like event (and also the secondary eclipses near days 1387 and 1402) and conclude that the apparent transit-like events near days 1390, 1403, and 1404 are instrumental artifacts. A similar analysis confirmed the reality of the two CBP transits in Sectors 6 and 10. We found the long- and short-cadence light curves to have different eclipse depths in some sectors because of different levels of background subtraction. In many cases, the short-cadence data overestimated the background, causing many “sky” pixels near the star to record negative flux values and the eclipses to appear artificially deep. We re-fit a background model for the short-cadence data using the FFIs, interpolating to the times of each short-cadence exposure, leading to consistent eclipse depths between the two data sets. We also tested a variety of different apertures, finding that the choice of aperture did not make a significant difference on the ultimate photometry, and in most sectors, we used the pipeline default aperture.

We achieved the best photometric precision for TOI-1338 using the PSF-based photometry built into *eleanor*. Our final adopted light curve is a combination of 30 minute cadence data extracted from the FFIs for Sectors 1, 2, and 3, and 2 minute cadence data extracted from the target pixels in the “postage stamps” for this particular target for the remaining sectors. This final light curve was detrended and normalized in the manner described in Orosz et al. (2019). As part of this iterative detrending process, we measure the durations of all eclipse and transit events. Most of the out-of-eclipse portions of the light curve are trimmed, and we keep only the out-of-eclipse regions that are within 0.25 days or $1.25\tau$ of each event (whichever is

greater), where $\tau$ is the duration of the event. The last primary eclipse in Sector 12 had unusually large residuals (for unknown reasons), and consequently, we excluded it from further analysis.

We measured the times of the eclipses and transits by fitting a simple model to the trimmed and normalized light curves for each event using the Eclipsing Light Curve (ELC) code of Orosz & Hauschildt (2000). This model has nine free parameters: the period *P*, the conjunction time $T_{conj}$, the inclination $i$, the primary radius $R_1$, the ratio of the radii $R_1/R_2$, two limb darkening parameters for the quadratic limb darkening law $q_{1,1}$, $q_{2,1}$, and two eccentricity parameters $\sqrt{e} \cos \omega$ and $\sqrt{e} \sin \omega$. The goal was to find a smooth and symmetric curve that best fit each segment, so no attempt was made to optimize more than one segment at a time. For each segment, we found the best-fitting model using the Differential Evolution Monte Carlo Markov Chain (DE-MCMC) algorithm of Ter Braak (2006) with 80 chains, which was run for 5000 generations. Using the best-fitting model, the individual uncertainties on each point were scaled in such a way to get $\chi^2_{best} = N$, where *N* is the number of points. The scale factors for the first five primary eclipses (observed in long cadence) were 0.347, 0.353, 0.626, 0.690, and 0.770, respectively. For primary eclipses observed in long cadence, the scale factors ranged from 1.042 to 1.376 with a median of 1.121. After the uncertainties on the data in each segment were scaled, the DE-MCMC code was run using 80 chains for 40,000 generations. Posterior samples were selected starting at generation 400 and sampling every 400th generation thereafter. The median of each posterior sample was adopted as the eclipse time, and the rms of the sample was taken to be the $1\sigma$ uncertainty. The measured primary and secondary cycle numbers and times are given in Table 1. The cycle numbers for the secondary eclipses are given as fractional values of the form $NN.45345$ (the mean phase of the secondary eclipses is not 0.5, owing to the orbital eccentricity).

To make the final light curve that was modeled using the photodynamical model described in Section 4 below, we simply combined the individual segments with scaled
<table><caption>**Table 1**
Observed Eclipse Times for TOI-1338</caption><thead><tr><th>Cycle</th><th>Observed Time<sup>a</sup></th><th>Model Time<sup>a</sup></th><th>O – C<sup>b</sup></th><th>Cycle</th><th>Observed Time<sup>a</sup></th><th>Model Time<sup>a</sup></th><th>O – C<sup>b</sup></th></tr><tr><td colspan="4">Primary</td><td colspan="4">Secondary</td></tr></thead><tbody><tr><td>0</td><td>...</td><td>3322.21319</td><td>...</td><td>0.45345</td><td>3328.83160 ± 0.00557</td><td>3328.83833</td><td>-9.71</td></tr><tr><td>1</td><td>3336.82137 ± 0.00027</td><td>3336.82171</td><td>-0.48</td><td>1.45345</td><td>3343.44848 ± 0.00970</td><td>3343.44687</td><td>2.32</td></tr><tr><td>2</td><td>3351.43036 ± 0.00026</td><td>3351.43028</td><td>0.12</td><td>2.45345</td><td>3358.05700 ± 0.00285</td><td>3358.05537</td><td>2.34</td></tr><tr><td>3</td><td>3366.03902 ± 0.00023</td><td>3366.03888</td><td>0.21</td><td>3.45345</td><td>3372.65843 ± 0.00410</td><td>3372.66390</td><td>-7.87</td></tr><tr><td>4</td><td>3380.64698 ± 0.00059</td><td>3380.64735</td><td>-0.52</td><td>4.45456</td><td>3387.27302 ± 0.00379</td><td>3387.27251</td><td>0.73</td></tr><tr><td>5</td><td>3395.25589 ± 0.00044</td><td>3395.25598</td><td>-0.12</td><td>5.45345</td><td>3401.86970 ± 0.01653</td><td>3401.88102</td><td>-16.29</td></tr><tr><td>6</td><td>...</td><td>3409.86452</td><td>...</td><td>6.45345</td><td>3416.49508 ± 0.00375</td><td>3416.48957</td><td>7.93</td></tr><tr><td>7</td><td>...</td><td>3424.47302</td><td>...</td><td>7.45345</td><td>3431.09816 ± 0.00340</td><td>3431.09811</td><td>0.06</td></tr><tr><td>8</td><td>3439.08165 ± 0.00020</td><td>3439.08157</td><td>0.13</td><td>8.45345</td><td>3445.70014 ± 0.00459</td><td>3445.70662</td><td>-9.34</td></tr><tr><td>9</td><td>3453.69010 ± 0.00019</td><td>3453.69017</td><td>-0.09</td><td>9.45345</td><td>3460.32094 ± 0.00364</td><td>3460.31513</td><td>8.36</td></tr><tr><td>10</td><td>...</td><td>3468.29871</td><td>...</td><td>10.45345</td><td>3474.92301 ± 0.00324</td><td>3474.92372</td><td>-1.02</td></tr><tr><td>11</td><td>3482.90704 ± 0.00022</td><td>3482.90721</td><td>-0.24</td><td>11.45345</td><td>3489.53315 ± 0.00300</td><td>3489.53225</td><td>1.30</td></tr><tr><td>12</td><td>3497.51571 ± 0.00021</td><td>3497.51582</td><td>-0.15</td><td>12.45345</td><td>...</td><td>3504.14078</td><td>...</td></tr><tr><td>13</td><td>3512.12410 ± 0.00022</td><td>3512.12432</td><td>-0.31</td><td>13.45345</td><td>3518.73613 ± 0.00491</td><td>3518.74934</td><td>-19.02</td></tr><tr><td>14</td><td>3526.73304 ± 0.00019</td><td>3526.73285</td><td>0.29</td><td>14.45345</td><td>...</td><td>3533.35786</td><td>...</td></tr><tr><td>15</td><td>3541.34140 ± 0.00022</td><td>3541.34142</td><td>-0.02</td><td>15.45345</td><td>3547.97016 ± 0.00444</td><td>3547.96636</td><td>5.48</td></tr><tr><td>16</td><td>...</td><td>3555.95001</td><td>...</td><td>16.45345</td><td>3562.57901 ± 0.00349</td><td>3562.57490</td><td>5.91</td></tr><tr><td>17</td><td>...</td><td>3570.55847</td><td>...</td><td>17.45345</td><td>3577.18768 ± 0.00349</td><td>3577.18350</td><td>6.02</td></tr><tr><td>18</td><td>3585.16712 ± 0.00032</td><td>3585.16712</td><td>0.02</td><td>18.45345</td><td>3591.78979 ± 0.00254</td><td>3591.79201</td><td>-3.19</td></tr><tr><td>19</td><td>...</td><td>3599.77564</td><td>...</td><td>19.45345</td><td>3606.40434 ± 0.00371</td><td>3606.40057</td><td>5.42</td></tr><tr><td>20</td><td>3614.38474 ± 0.00030</td><td>3614.38415</td><td>0.85</td><td>20.45345</td><td>3621.01281 ± 0.00399</td><td>3621.00911</td><td>5.32</td></tr><tr><td>21</td><td>3628.99361 ± 0.00038</td><td>3628.99271</td><td>1.31</td><td>21.45345</td><td>3635.61695 ± 0.00322</td><td>3635.61762</td><td>-0.97</td></tr><tr><td>22</td><td>3643.60114 ± 0.00036</td><td>3643.60131</td><td>-0.23</td><td>22.45345</td><td>3650.22268 ± 0.00518</td><td>3650.22613</td><td>-4.97</td></tr></tbody></table>

**Notes.**

uncertainties from all of the eclipse and transit events. Since each segment has had its uncertainties scaled individually, no one event should be given undue weight owing to underestimated uncertainties for the photometric measurements. The scale factors for the segments in long cadence were smaller than one, which suggests the photometric errors were overestimated. On the other hand, the scale factors for the event observed in short cadence were all larger than one, which suggests the photometric errors were underestimated.

### 3. Complementary Observations

#### 3.1. Radial-velocity Characterization

After its initial classification as an eclipsing binary, EBLM J0608-59 was observed with the CORALIE spectrograph to constrain the component masses. Between 2009 December and 2012 April 19, radial-velocity observations were performed to map out the Keplerian orbit of the binary (Triaud et al. 2017a). Exposures were typically 600 s long, yielding a median precision of 34 m s⁻¹. Owing to the small mass ratio of the binary, and indeed of the entire EBLM sample by construction, the system appears as a single-line spectroscopic binary. We expect a V magnitude difference of 9.8 between the primary and the secondary, given that $M_1 \sim 1.1 M_\odot$ and $M_2 \sim 0.3 M_\odot$, respectively (Triaud et al. 2017a). As a consequence of the large flux ratio between the primary and secondary, the spectral lines of the secondary are not noticeable in the observations, thereby allowing high-accuracy measurements of the primary’s radial velocity. Likewise, the secondary star is too faint for TESS to be able to detect transits of the CBP across it.

#### 3.1.1. The BEBOP Radial-velocity Search for CBPs

In late 2013, the Binaries Escorted By Orbiting Planets (BEBOP) program was created as a radial-velocity survey for the detection of CBPs. The BEBOP sample is exclusively constructed from eclipsing single-line spectroscopic binaries to avoid contamination effects that severely hinder planet detection in double-line spectroscopic binaries (Konacki et al. 2009). An initial target list of roughly 50 binaries was created from the larger EBLM sample. Selection criteria include the obtainable radial-velocity precision and the lack of stellar activity. EBLM J0608-59 was included in this initial BEBOP selection. Between 2014 November and 2015 August, we acquired 17 additional CORALIE measurements with longer exposures of 1800 s to increase the precision of the radial-velocity measurements (the median uncertainty was 25 m s⁻¹). This second set of observations coincided with a change to a new octagonal fiber. The new fiber provided greater long-term stability compared to the original circular one but at the cost of ~10% of the incoming flux. Such a fiber change may also induce a small velocity offset (Triaud et al. 2017b). Consequently, when modeling the data, we treat the CORALIE data sets as if they were acquired from two different instruments, with a free parameter for the offset.

The combined radial velocities from both the EBLM and BEBOP surveys are presented in Martin et al. (2019). No CBP was detected in the radial-velocity of this system within the sensitivity of the observations, which were adequately fit by a single-Keplerian orbit ($\chi^2 \sim 1$). This ruled out (at 95% confidence) the presence in the system of a CBP more massive than ~0.7 $M_{Jup}$ and with an orbital period roughly six times shorter than that of the binary.

<footnote><sup>a</sup> BJD—2,455,000.</footnote>

<footnote><sup>b</sup> Observed time minus model time in minutes.</footnote>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left panel: radial velocities for TOI-1338, phase-folded on the binary orbital period of 14.61 days, taken with the CORALIE instrument (blue and orange) and HARPS (red). Right panel: residuals to the single-Keplerian radial-velocity fit over time. The change in CORALIE residuals between blue and orange measurements coincide with a fiber change (see the text for details).</figcaption></figure>

<table><caption>**Table 2**

Photometric and Spectroscopic Parameters of the System Prior to the Photodynamical Solution</caption><thead><tr><td>Parameter</td><td>Value</td><td>Uncertainty</td><td>Unit</td><td>Source</td></tr></thead><tbody><tr><td>RV semi-amplitude, K<sub>1</sub></td><td>21.619</td><td>0.007</td><td>km s<sup>-1</sup></td><td>Martin et al. (2019)</td></tr><tr><td>Gravity of Primary, log g<sub>1</sub></td><td>4.0</td><td>0.08</td><td>cgs</td><td>Spectroscopy, this work</td></tr><tr><td>Metallicity of Primary, [Fe/H]<sub>1</sub></td><td>0.01</td><td>0.05</td><td>dex</td><td>Spectroscopy, this work</td></tr><tr><td>Projected Rotational Velocity of Primary, v sin i</td><td>3.6</td><td>0.6</td><td>km s<sup>-1</sup></td><td>Spectroscopy, this work</td></tr><tr><td>Reddening, E(B-V)</td><td>0.02</td><td>0.01</td><td>mag</td><td>Gaia + Photometry, this work</td></tr><tr><td>Effective Temperature of Primary, T<sub>eff,1</sub></td><td>6050</td><td>80</td><td>K</td><td>Spectroscopy, this work</td></tr><tr><td>Effective Temperature of Primary, T<sub>eff,1</sub></td><td>5990</td><td>110</td><td>K</td><td>Gaia + Photometry, this work</td></tr><tr><td>Radius of Primary, R<sub>1</sub></td><td>1.345</td><td>0.046</td><td>R<sub>☉</sub></td><td>Gaia+Photometry, this work</td></tr><tr><td>Age</td><td>4.4</td><td>0.2</td><td>Gyr</td><td>This work</td></tr></tbody></table>

Starting in 2018 April, the BEBOP survey was extended to the HARPS high-resolution spectrograph on the ESO 3.6 m telescope (Prog.ID 1101.C-0721, PI A. Triaud; Pepe et al. 2002), also at La Silla, Chile.<fnref n="50" /> Compared to CORALIE, HARPS benefits from a larger telescope aperture, a higher resolving power of R = 115,000, and greater radial-velocity stability by being both thermally stabilized and operated under vacuum. Seven HARPS spectra of the target have been acquired to date. They were reduced with the HARPS pipeline, which has been shown to achieve remarkable precision and accuracy (e.g., Mayor et al. 2009; López-Morales et al. 2014). The radial velocities were computed by using a binary mask corresponding to a G2 spectral-type template (Baranne et al. 1996). We achieve a median radial-velocity precision of 5.9 m⁻¹ in our measurements. A fit to the complete set of radial velocities (CORALIE and HARPS), shown in Figure 5, produces a fit of reduced χ²/v ~ 1 mostly because the number of HARPS measurements is close to the number of free parameters for the binary orbit. The model adjusts to the HARPS measurements first, because they have the greatest weights. CORALIE’s precision is not sufficient to detect the additional planetary signal. The HARPS measurements indicate that this system has very little activity, making it optimal for radial-velocity measurements. The BEBOP survey is ongoing, and more spectra from both HARPS and ESPRESSO will be

obtained for TOI-1338/EBLM J0608-59 and published in a subsequent paper.

## 3.2. Spectroscopic Characterization

The seven extracted HARPS spectra were co-added onto a common wavelength axis reaching a signal-to-noise ratio of approximately 73. The resulting spectrum was analyzed with the spectral analysis package ISPEC (Blanco-Cuaresma et al. 2014). We used the synthesis method to fit individual spectral lines of the co-added spectra. The radiative transfer code SPECTRUM (Gray & Corbally 1994) was used to generate model spectra with MARCS model atmospheres (Gustafsson et al. 2008), version 5 of the GAIA ESO survey atomic line list provided within ISPEC and solar abundances from Asplund et al. (2009). Macroturbulence is estimated using equation (5.10) from Doyle (2015), and microturbulence was accounted for at the synthesis stage using equation (3.1) from the same source. The Hα, Na I D, and Mg I b lines were used to infer the effective temperature T<sub>eff</sub> and gravity log g, while the Fe I and Fe II lines were used to determine the metallicity [Fe/H] and the projected rotational velocity v sin i. Trial synthetic model spectra were fit until an acceptable match to the data was found. Uncertainties were estimated by varying individual parameters until the model spectrum was no longer well-matched to the spectra of TOI-1338. For the primary star, we find an effective temperature of T<sub>eff,1</sub> = 6050 ± 80 K, a metallicity of [Fe/H]<sub>1</sub> = 0.01 ± 0.05, a gravity of log g<sub>1</sub> = 4.0 ± 0.08 dex, and a projected rotational velocity of v sin i<sub>*,1</sub> = 3.6 ± 0.6 km s<sup>-1</sup>. These measurements are summarized in Table 2.

<footnote n="50">⁵⁰ BEBOP also surveys the northern skies, using SOPHIE, at the Observatoire de Haute-Provence, Prog.ID 19A.PNP.SANT, PI A. Santerne. The northern EBLM sample will also be included in the TESS 2 minute cadence under proposal G022253, PI D. Martin.</footnote>
As a check on the spectroscopic temperature of the primary, we first gathered brightness measurements from the literature in the Johnson, Tycho-2, Two Micron All Sky Survey, and Sloan systems, constructed 10 nonindependent color indices, and corrected each for reddening using the extinction law of Cardelli et al. (1989) with a value of $E(B - V) = 0.020 \pm 0.010$ derived from the extinction map of Schlafly & Finkbeiner (2011) and the Gaia (Gaia Collaboration et al. 2016; Gaia Collaboration 2018) distance of $405 \pm 3$ pc. These colors are unaffected by the secondary star because it is so faint. We then used color-temperature calibrations from Casagrande et al. (2010) and Huang et al. (2015) to infer a mean photometric temperature of $5990 \pm 110$ K, in good agreement with the spectroscopic value.

We estimated the radius of the primary star in three different ways. One method used the procedure outlined by Stassun et al. (2019) for the preparation of the TIC-8 catalog, involving the Gaia parallax and *G* magnitude, an extinction correction, the spectroscopic temperature, and the *G*-band bolometric correction. Another method we used was based on a fit to the spectral energy distribution (SED) performed with EXOFASTv2 (Eastman et al. 2019) and the MIST bolometric correction tables,<fnref n="51" /> using brightness measurements in the Gaia and Wide-field Infrared Survey Explorer systems in addition to those mentioned earlier. Suitable priors were placed on the distance, extinction, temperature, metallicity, and $log g_1$. A third method we used was also based on an SED fit but used instead the NextGen model atmospheres of Allard et al. (2012). The three procedures gave very similar results. We adopt the value $R_1 = 1.345 \pm 0.046 R_\odot$ in the following to use as a prior for the photometric-dynamical modeling described below.

### 3.2.1. Image Analysis

The relatively large sizes of the TESS pixels, approximately 21" on a side, leave the target susceptible to photometric contamination from nearby stars, including additional wide stellar companions. The nearest object (TIC 260128336) to TOI-1338, after proper motion correction, is separated by 53"/6. The SPOC Data Validation centroid offsets for the Sectors 1–12 multi-sector run indicated that the primary and secondary eclipses for TOI-1338 originate from the target itself (Jenkins et al. 2016); complementary analysis with the photocenter module of DAVE (Kostov et al. 2019) confirmed this. We also searched for nearby sources with SOAR speckle imaging (see Tokovinin et al. 2018, 2019, for details of the instrumentation) on 2019 March 17 UT, observing in a similar visible bandpass as TESS. We detected no nearby sources within 3" of TOI-1338. The detection sensitivity and speckle autocorrelation function of the SOAR observations are plotted in Figure 6.

## 4. Photometric-dynamical Analysis of the System

Due to the rich dynamical interactions between the two stars and a planet in a CBP system, the deviations of the planet's orbit from a strictly periodic one are much more pronounced compared to a single-star system. The secular evolution in a CBP system can occur on a timescale as short as a decade instead of thousands of years (like in the solar system). Thus, measurable changes in quantities such as the inclination of the planet's orbit can be observed with missions like Kepler and

TESS. Given the relatively rapid secular evolution, the orbits in a CBP system are not simple Keplerians. A complete description of a CBP system relies on a large number of parameters—e.g., masses, radii, and orbital parameters for the two stars and the planet(s), radiative parameters for the two stars, etc. Stellar eclipses, both primary and secondary, allow for precise measurements of times of conjunction. In addition, transits of the CBP across the primary and/or the secondary star can provide precise position measurements of both the stars and the planet at times other than the time of stellar conjunction. Thus, the dynamical complexity of the system—while computationally challenging—enables precise measurements of the system's parameters. For example, the stellar masses and radii of the two stars in the Kepler-16 CBP system have been measured to sub-percent precision (Doyle et al. 2011).

### 4.1. ELC Modeling

To obtain a complete solution for the TOI-1338 system, we carried out a photometric-dynamical analysis with the ELC code (Orosz & Hauschildt 2000), utilizing the photometry from TESS and the precise radial-velocity measurements from CORALIE and HARPS. For this task, the ELC code combines *N*-body simulations, modified to include tidal interaction and general relativistic effects, with a photometric model for the stellar eclipses and planetary transits, to reproduce a light curve of the system. These modifications to ELC to allow for modeling stellar triple and higher-order systems, and CBP systems, have been described in Welsh et al. (2015) and Orosz et al. (2019). This code has been used extensively for the analysis and confirmation of Kepler CBPs; for the sake of completeness, we outline it below.

Briefly, given instantaneous orbital parameters (e.g., the orbital period, the eccentricity, the inclination, etc.) at some reference epoch and the masses of each body, ELC solves the Newtonian equations of motion using a symplectic integrator, in this case a 12th-order Gaussian Runge-Kutta integrator (Hairer et al. 2002). When necessary, the Newtonian equations of motion can be modified to account for general relativistic precession and tidal effects (see Mardling & Lin 2002). The solution of the dynamical equations enables the positions of all of the bodies to be specified at any given time. Then, given the positions of the bodies on the plane of the sky, their radii, and their radiative properties, the observed flux is computed for any number of overlapping bodies using the algorithm discussed in Short et al. (2018). Likewise, the solution to the dynamical equations gives the radial velocity of each body at any given time. It is also possible to include other observable quantities that do not depend on time (for example, the surface gravity of the primary star) in the fitting process.

For TOI-1338, we initially had the following 25 free parameters. The binary orbit is specified by the time of conjunction $T_{\text{copj,bin}}$, the period $P_{\text{bin}}$, the eccentricity parameters $\sqrt{e_{\text{bin}} \cos \omega_{\text{bin}}}$, and $\sqrt{e_{\text{bin}} \sin \omega_{\text{bin}}}$ (where $e_{\text{bin}}$ is the eccentricity and $\omega_{\text{bin}}$ is the argument of periastron), and the inclination $i_{\text{bin}}$. We fix the initial nodal angle of the binary orbit to $\Omega_{\text{bin}} = 0$. The stellar masses are specified by the primary mass $M_1$ (in units of $M_\odot$) and the mass ratio $Q = M_2/M_1$, the stellar radii—by the primary radius $R_1$ and

<footnote n="51">51 http://waps.cfa.harvard.edu/MIST/model_grids.html</footnote>

<footnote n="52">52 Combinations and/or ratios of individual parameters are often more convenient to use.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The 5σ detection limits of the SOAR speckle imaging of TOI-1338; the inset is the speckle autocorrelation function. No nearby stars were detected within 3" of the target.</figcaption></figure>

the radius ratio $R_1/R_2$, and the stellar temperatures—by $T_{\text{eff},1}$ and the temperature ratio $T_{\text{eff},2}/T_{\text{eff},1}$. We assume a quadratic limb darkening law with the Kipping (2013) “triangular” resampling, for a total of four additional parameters. Finally, to account for tidal precession, we use two additional free parameters—the apsidal constants $k_{2,1}$ and $k_{2,2}$ for the primary and secondary, respectively. For the planet’s orbit, we fit for the sum and differences of the period and time of barycentric conjunction, i.e., $\zeta_p = (P_{\text{pl}} + T_{\text{conj,pl}})$, $\zeta_m = (P_{\text{pl}} - T_{\text{conj,pl}})$, the eccentricity parameters $\sqrt{e_{\text{pl}} \cos \omega_{\text{pl}}}$, $\sqrt{e_{\text{pl}} \sin \omega_{\text{pl}}}$, the inclination $i_{\text{pl}}$, and the nodal angle $\Omega_{\text{pl}}$. The planet’s mass is specified by $M_3$, where the units are $M_\oplus$. Finally, the radius of the planet is specified by the ratio $R_1/R_{\text{pl}}$.

The observables for the TOI-1338 system are the trimmed and normalized light curve from TESS, and the radial-velocity measurements from CORALIE and HARPS. The radial velocities span three distinct data sets where an offset velocity is found for each, namely the “early” CORALIE data before the fiber change (roughly days 200 to 1000 in units of BJD—2,455,000), the “late” CORALIE data after the fiber change (roughly days 1900 to 2250), and HARPS (roughly days 3200 to 3550). Another set of observables is the spectroscopically measured gravity for the primary star (log $g_1 = 4.0 \pm 0.08$ dex in cgs units) and the radius of the primary star ($R_1 = 1.345 \pm 0.046 R_\odot$) as determined from multicolor photometry and the distance. We use the usual $\chi^2$ statistic for the likelihood, where the $\chi^2$ contributions for the light curve, the velocity curve, the spectroscopic gravity, and the computed radius are combined. In general, model parameters are drawn uniformly within specified bounds. The use of the primary radius $R_1$ as both a fitting parameter and an observed parameter effectively gives that parameter a Gaussian prior. Finally, we also include estimates of the times of the planet transits with generous uncertainties (0.0050 days for all events) in the overall $\chi^2$. Using estimates of the transit times is very important early in the fitting process, as models with transits that are very far away from the observed times would have the same $\chi^2$ contribution from the light curve as models that are much closer to the observed times since the out-of-transit parts of the light curve are flat. Including observed transit times in the overall $\chi^2$ gives a larger penalty to the models that miss the

observed times by larger amounts compared to models that have smaller differences between the model transit times and the observed ones. As the models converge to the “true” solution, the $\chi^2$ contribution from the eclipse times becomes very small, leaving the actual fit to each individual transit profile in the light curve to determine the $\chi^2$ contribution to the overall total $\chi^2$.

The reference time for the osculating orbital parameters was chosen to be day 186 in units of BJD—2,455,000. The ending time of the numerical integration was day 3750 in the same units. For these integrations, we include the effects of both precession due to general relativity and precession due to tides. As some measure of the round-off errors that afflict nearly every numerical integrator, we measure the position and velocity of the system center of mass as a function of time. These two quantities start out at the origin, [e.g., at the coordinates (0, 0, 0)] by default, and should remain at the origin for the ideal integrator. We find that the center of mass at day 3750 is offset by $\approx 3.40 \times 10^{-13}$ au ($\approx 5.09$ cm) and that the velocity at this time is offset by $\approx 2.00 \times 10^{-16}$ au per day ($\approx 3.46 \times 10^{-10}$ meters per second). Thus, we conclude that numerical round-off error is not an issue for these integrations.

ELC has a number of optimizers available, and the two that proved to be the most useful for the analysis of TOI-1338 were the genetic algorithm of Charbonneau (1995) and the DE-MCMC algorithm of Ter Braak (2006). The spectroscopic parameters of the binary were well known from analytic models of the radial velocities. The other parameters, especially the orbital parameters for the planet, were initially not well constrained. Thus, our fitting proceeded by iteration. The genetic code was run where an initial “population” of 100 models was evolved for a few thousand generations. Next, the top 80 models were evolved using the DE-MCMC code for a few thousand more generations. The top few models, randomly chosen to be between one and five, were put back into the genetic code along with random models and evolved. After a few iterations of this process, a reasonably good model was found.

After an initial good model was found, we performed several pilot runs of the DE-MCMC code. The initial population of models consisted of the best model and mutated copies of the best model where several randomly chosen parameters were offset from their optimal values by small amounts drawn from a normal distribution. The standard deviation of the normal distribution was chosen to be between 0.001 and 0.01 times the parameter value. We found that the chains quickly spread out and achieved their final overall spreads after a burn-in period of typically 1000 generations. During these pilot runs, we checked to confirm that the prior ranges included support for the entire range with nontrivial likelihood and that the likelihood falls to extremely small values by the time model parameters reach any of these boundaries (with the possible exception of hard physical boundaries). We also found that the two apsidal parameters and the limb darkening coefficients for the secondary star were not constrained. Consequently, we fixed their values at representative values ($k_{2,1} = 0.01$, $k_{2,2} = 0.10$ for the apsidal parameters and $q_{1,2} = 0.06$, $q_{2,2} = 0.41$ for the secondary’s quadratic law limb darkening coefficients), reducing the number of free parameters to 21.

For the final step leading to the adopted parameters and their uncertainties, we used a brute-force “grid search” algorithm to find optimal models with the third body mass fixed at values
<table><caption>**Table 3**
Fitted Parameters for TOI-1338</caption><thead><tr><th>Parameter<sup>a</sup></th><th>Best</th><th>Mode</th><th>Median</th><th>+1σ</th><th>-1σ</th></tr></thead><tbody><tr><td>T<sub>conj,bin</sub></td><td>3336.8245</td><td>3336.8237</td><td>3336.8242</td><td>0.0025</td><td>0.0023</td></tr><tr><td>P<sub>bin</sub> (day)</td><td>14.608561</td><td>14.608559</td><td>14.608559</td><td>0.000013</td><td>0.000012</td></tr><tr><td>√e<sub>bin</sub> cos ω<sub>bin</sub></td><td>-0.18275</td><td>-0.18270</td><td>-0.18272</td><td>0.00040</td><td>0.00040</td></tr><tr><td>√e<sub>bin</sub> sin ω<sub>bin</sub></td><td>0.35015</td><td>0.35008</td><td>0.35020</td><td>0.00036</td><td>0.00035</td></tr><tr><td>M<sub>1</sub> (M<sub>☉</sub>)</td><td>1.038</td><td>1.149</td><td>1.127</td><td>0.068</td><td>0.069</td></tr><tr><td>Q ≡ M<sub>2</sub>/M<sub>1</sub></td><td>0.2865</td><td>0.2752</td><td>0.2774</td><td>0.0069</td><td>0.0062</td></tr><tr><td>i (deg)</td><td>89.658</td><td>89.649</td><td>89.696</td><td>0.178</td><td>0.114</td></tr><tr><td>T<sub>eff,1</sub> (K)</td><td>5990.7</td><td>6072.4</td><td>6040.8</td><td>98.1</td><td>91.5</td></tr><tr><td>T<sub>eff,2</sub>/T<sub>eff,1</sub></td><td>0.5537</td><td>0.5548</td><td>0.5516</td><td>0.0042</td><td>0.0047</td></tr><tr><td>R<sub>1</sub> (R<sub>☉</sub>)</td><td>1.299</td><td>1.338</td><td>1.331</td><td>0.024</td><td>0.026</td></tr><tr><td>R<sub>1</sub>/R<sub>2</sub></td><td>4.307</td><td>4.310</td><td>4.308</td><td>0.013</td><td>0.014</td></tr><tr><td>q<sub>1,1</sub></td><td>0.255</td><td>0.245</td><td>0.255</td><td>0.048</td><td>0.043</td></tr><tr><td>q<sub>2,1</sub></td><td>0.308</td><td>0.299</td><td>0.310</td><td>0.058</td><td>0.050</td></tr><tr><td>ζ<sub>p</sub> (day)</td><td>3436.57</td><td>3437.26</td><td>3437.32</td><td>0.83</td><td>0.78</td></tr><tr><td>ζ<sub>m</sub> (day)</td><td>-3246.29</td><td>-3246.84</td><td>-3246.97</td><td>0.71</td><td>0.77</td></tr><tr><td>√e<sub>pl</sub> cos ω<sub>pl</sub></td><td>-0.035</td><td>-0.051</td><td>-0.058</td><td>0.018</td><td>0.024</td></tr><tr><td>√e<sub>pl</sub> sin ω<sub>pl</sub></td><td>-0.303</td><td>-0.285</td><td>-0.290</td><td>0.009</td><td>0.011</td></tr><tr><td>i<sub>pl</sub> (deg)</td><td>89.22</td><td>89.29</td><td>89.37</td><td>0.35</td><td>0.00</td></tr><tr><td>Ω<sub>pl</sub> (deg)</td><td>0.87</td><td>0.83</td><td>0.91</td><td>0.35</td><td>0.35</td></tr><tr><td>M<sub>3</sub> (M<sub>☉</sub>)</td><td>30.2</td><td>33.2</td><td>33.0</td><td>20.3</td><td>19.6</td></tr><tr><td>R<sub>1</sub>/R<sub>3</sub></td><td>21.12</td><td>21.20</td><td>21.17</td><td>0.47</td><td>0.44</td></tr><tr><td>γ<sub>1</sub><sup>b</sup> (km s<sup>-1</sup>)</td><td>30.74769</td><td>30.74766</td><td>30.74778</td><td>0.00089</td><td>0.00084</td></tr><tr><td>γ<sub>2</sub><sup>c</sup> (km s<sup>-1</sup>)</td><td>30.74621</td><td>30.74618</td><td>30.74641</td><td>0.00128</td><td>0.00125</td></tr><tr><td>γ<sub>3</sub><sup>d</sup> (km s<sup>-1</sup>)</td><td>30.75866</td><td>30.75836</td><td>30.75857</td><td>0.00098</td><td>0.00098</td></tr></tbody></table>

**Notes.**

a Osculating parameters valid at BJD 2,455,186,0000.

b Relative velocity offset, “early” CORALIE data.

c Relative velocity offset, “late” CORALIE data.

d Relative velocity offset, HARPS data.

from 0.3 to 117.3 $M_{⊕}$ in steps of 1 $M_{⊕}$ as “seed” models for the final runs of the DE-MCMC code. We then ran the DE-MCMC code eight separate times, each with the same seed models, but with a different initial random number seed. All eight runs used 120 chains, and each was run for at least 8800 generations (the longest run had 38,300 generations). The posterior samples were drawn starting at generation 3000, with subsequent draws that skipped every 2000 generations until the chains ended. The individual posterior samples were combined into single samples (with $N=12,120$) for each fitting parameter and for several derived parameters of interest. Table 3 provides the parameters for the best-fitting model, the mode of the posterior sample (found using 50 bins), the median of the posterior sample, and the $+1\sigma$ and $-1\sigma$ uncertainties. Table 4 provides several derived parameters of astrophysical interest using a similar format as Table 3. In the discussion that follows, we use the posterior medians as the adopted parameter values. The model fits to the stellar eclipses are shown in Figure 7, and the model fits to the planet transits are shown in Figure 8.

There is a slight difference between the depths of the long-cadence primary eclipses and the model, such that the former are deeper. To address this issue, we reran the fit using a negative contamination factor for the long-cadence data in order to make the model deeper and match the data. The best-fit model is shown in Figure 9. The best-fit dilution factor is $-0.0336 \pm 0.0022$, indicating that there may be a variable flux offset between the long- and short-cadence data. To first order, this suggests that either the sky background in the long-cadence data was over-subtracted by 3.4% or the sky background in the

short-cadence data was under-subtracted by 3.4%. This dilution term seems to be somewhat on the high side of what one might expect based on the number of counts in the actual data. It is also hard to determine to which data set it has to be applied (long-cadence only, short-cadence only, or a combination of both). New TESS observations in Cycle 3 will help address this issue. Regardless of the specific reason and approach, the key parameters (mass, radius, etc.) do not change significantly and the results are consistent between the models with and without a dilution factor. The radial velocities are fit quite well, and Figure 10 shows the residuals for the “early” CORALIE, the “late” CORALIE, and the HARPS measurements. Finally, in Table 5, we give the initial dynamical parameters and Cartesian coordinates for the best-fitting model to full machine precision.

## 4.2. Planet Mass

As noted in Section 2.3, there is undoubtedly a transiting circumbinary object in the TOI-1338 system. Our analysis here shows that the mass of this object is well within the planetary regime, i.e., $M_3 = 33 \pm 20 M_⊕$. Here, we give a brief discussion of which feature(s) in the data allow us to make this determination. The top panel of Figure 11 shows the posterior distribution of the planet mass $M_3$ in units of $M_⊕$. The minimum and maximum values in the posterior sample are 0.01 $M_⊕$ and 117.05 $M_⊕$, respectively. The bottom panel of Figure 11 shows $\chi^2 - \chi^2_{min}$ for all computed models, where $\chi^2_{min} = 14, 587.93$. We see that $\Delta\chi^2$ is about 4 when $M_3 = 0$, and it is about 18 when $M_3 = 117 M_⊕$. As a reminder, the total
<table><caption>**Table 4**
Derived Parameters for TOI-1338</caption><thead><tr><th>Parameter<sup>a</sup></th><th>Best</th><th>Mode</th><th>Median</th><th>+1σ</th><th>-1σ</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center; font-weight: bold;">Bulk Properties</td></tr><tr><td>M<sub>1</sub> (M<sub>☉</sub>)</td><td>1.038</td><td>1.149</td><td>1.127</td><td>0.068</td><td>0.069</td></tr><tr><td>R<sub>1</sub> (R<sub>☉</sub>)</td><td>1.299</td><td>1.338</td><td>1.331</td><td>0.024</td><td>0.026</td></tr><tr><td>M<sub>2</sub> (M<sub>☉</sub>)</td><td>0.2974</td><td>0.3168</td><td>0.3128</td><td>0.0113</td><td>0.0118</td></tr><tr><td>R<sub>2</sub> (R<sub>☉</sub>)</td><td>0.3015</td><td>0.3102</td><td>0.3089</td><td>0.0056</td><td>0.0060</td></tr><tr><td>M<sub>3</sub> (M<sub>⊕</sub>)</td><td>30.2</td><td>33.2</td><td>33.0</td><td>20.3</td><td>19.6</td></tr><tr><td>R<sub>3</sub> (R<sub>⊕</sub>)</td><td>6.71</td><td>6.83</td><td>6.85</td><td>0.19</td><td>0.19</td></tr><tr><td>ρ<sub>3</sub> (g cm<sup>-3</sup>)</td><td>0.55</td><td>0.57</td><td>0.56</td><td>0.34</td><td>0.33</td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Binary Orbit</td></tr><tr><td>P<sub>bin</sub> (day)</td><td>14.608561</td><td>14.608559</td><td>14.608559</td><td>0.000013</td><td>0.000012</td></tr><tr><td>T<sub>conj,bin</sub></td><td>3336.8245</td><td>3336.8237</td><td>3336.8242</td><td>0.0025</td><td>0.0023</td></tr><tr><td>K<sub>bin</sub> (km s<sup>-1</sup>)</td><td>21.6248</td><td>21.6246</td><td>21.6247</td><td>0.0034</td><td>0.0032</td></tr><tr><td>e<sub>bin</sub></td><td>0.15601</td><td>0.15603</td><td>0.15603</td><td>0.00015</td><td>0.00015</td></tr><tr><td>ω<sub>bin</sub> (deg)</td><td>117.561</td><td>117.568</td><td>117.554</td><td>0.072</td><td>0.074</td></tr><tr><td>a<sub>bin</sub> (au)</td><td>0.1288</td><td>0.1319</td><td>0.1321</td><td>0.0024</td><td>0.0025</td></tr><tr><td>True anomaly (deg)</td><td>111.217</td><td>111.246</td><td>111.226</td><td>0.071</td><td>0.069</td></tr><tr><td>Mean anomaly (deg)</td><td>93.882</td><td>93.897</td><td>93.889</td><td>0.065</td><td>0.065</td></tr><tr><td>Mean longitude (deg)</td><td>228.779</td><td>228.783</td><td>228.779</td><td>0.020</td><td>0.020</td></tr><tr><td>i<sub>bin</sub> (deg)</td><td>89.658</td><td>89.649</td><td>89.696</td><td>0.178</td><td>0.114</td></tr><tr><td>Ω<sub>bin</sub> (deg)</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td></tr><tr><td colspan="6" style="text-align: center; font-weight: bold;">Planet Orbit</td></tr><tr><td>P<sub>pl</sub> (day)</td><td>95.141</td><td>95.175</td><td>95.174</td><td>0.031</td><td>0.035</td></tr><tr><td>T<sub>conj, pl</sub></td><td>3341.43</td><td>3342.11</td><td>3342.15</td><td>0.80</td><td>0.74</td></tr><tr><td>e<sub>pl</sub></td><td>0.0928</td><td>0.0861</td><td>0.0880</td><td>0.0043</td><td>0.0033</td></tr><tr><td>ω<sub>pl</sub> (deg)</td><td>263.3</td><td>260.3</td><td>258.6</td><td>3.7</td><td>4.8</td></tr><tr><td>a<sub>pl</sub> (au)</td><td>0.4491</td><td>0.4639</td><td>0.4607</td><td>0.0084</td><td>0.0088</td></tr><tr><td>True anomaly (deg)</td><td>136.0</td><td>141.0</td><td>141.7</td><td>5.2</td><td>4.5</td></tr><tr><td>Mean anomaly (deg)</td><td>128.3</td><td>134.1</td><td>135.2</td><td>6.1</td><td>5.4</td></tr><tr><td>Mean longitude (deg)</td><td>400.2</td><td>401.6</td><td>401.2</td><td>0.9</td><td>1.1</td></tr><tr><td>i<sub>pl</sub> (deg)</td><td>89.22</td><td>89.29</td><td>89.37</td><td>0.35</td><td>0.26</td></tr><tr><td>Ω<sub>pl</sub> (deg)</td><td>0.87</td><td>0.83</td><td>0.91</td><td>0.35</td><td>0.35</td></tr><tr><td>i<sup>b</sup> (deg)</td><td>0.98</td><td>0.92</td><td>0.99</td><td>0.31</td><td>0.28</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Osculating parameters valid at BJD 2,455,186.0000.

<sup>b</sup> Mutual inclination between orbital planes.

χ² has contributions from the TESS light curve, the three radial-velocity sets, the measured gravity, and the measured radius of the primary star. For the best-fitting model, these contributions are 14,542.39 (TESS), 13.88 (early CORALIE), 17.11 (late CORALIE), 5.44 (HARPS), 8.07 (log g), and 1.01 (R₁). For the best model with M₃ = 120 M₄, those values are 14,545.53 (TESS), 17.67 (early CORALIE), 20.31 (late CORALIE), 14.85 (HARPS), 9.72 (log g), and 0.02 (R₁). There is hardly any change in the χ² fit to the TESS light curve between the two models. The data set with the largest change in χ² is HARPS, where the χ² changed by about 9.4. When the planet mass is fixed at 150 M₄, the total χ² is 14,626.61 (38.68 larger than the overall best model), and the individual contributions are 14,553.20 (TESS), 17.96 (early CORALIE), 19.46 (late CORALIE), 25.32 (HARPS), 10.05 (log g), and 0.04 (R₁). Although the χ² for the TESS light curve got slightly worse, it seems that the HARPS radial-velocity measurements have the most sensitivity to the planet mass. The bottom panel of Figure 10 shows the residuals of the HARPS measurements for the best overall model, the best model with M₃ = 120 M₄, and the best model with M₃ = 150 M₄. The residuals of the first measurement near day 3220 and the last two measurements

near days 3560 and 3570 show the most variation with the changing planet mass. Thus, in the near term, the most effective way to better constrain the mass of the planet would be to obtain more radial-velocity measurements with a quality similar to or better than the HARPS measurements presented here.

A CBP can perturb the binary and give rise to eclipse time variations (ETVs). The interaction between the planet and the binary can cause the binary orbit to precess, which leads to changes in the phase difference between the primary and secondary eclipses. When one attempts to fit the primary and secondary eclipse times to a common ephemeris, the *O* − *C* (observed minus computed) values for the primary eclipses will have the opposite slope that the *O* − *C* values for the secondary eclipses have. Figure 12 shows the common period *O* − *C* diagram for the model primary and secondary eclipses over the whole time span of the radial-velocity and photometric observations. For our overall best-fitting model, the argument of periastron ω changes by 0.0005715 degrees per cycle. The contribution of this precession from General Relativity is 0.0001132 degrees per cycle, and the contribution from tides is 0.0000055 degrees per cycle. This precession causes a
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The folded primary and secondary eclipses are shown with the best-fitting model (solid red line). Upper left panel: primary eclipse in with residuals. Upper right panel: primary eclipse in and 30 minute cadence with residuals. Lower left panel: secondary eclipse in with residuals. Lower right panel: secondary eclipse in 30 minute cadence with residuals. There is a slight trend of the model being shallower than the long-cadence primary eclipses.</figcaption></figure>

divergence between the primary and secondary *O* − *C* curves of about 2 minutes over the roughly 10 yr time span of the data. When the planet mass is fixed at 150 $M_⊕$, the best-fitting model for that mass has a change in $\omega$ of 0.002129 degrees per cycle. This results in a divergence between the primary and secondary *O* − *C* curves of about 9 minutes. As a practical matter, we only have measurements of eclipse times over the last 1.5 yr or so, and the uncertainties are relatively large: 0.36 minutes for the primary eclipses and 5.40 minutes for the secondary eclipses. Unless the measurements of the eclipse times can be vastly improved, we would need many more years of eclipse time measurements before the time baseline is long enough to accumulate a measurable divergence in the common period *O* − *C* diagram.

### 4.3. Dynamical Evolution

The large tidal potential produced by the inner binary causes the orbital elements of the CBP to vary with time. Indeed, the

best-fit osculating orbital elements (Table 4) represent only a snapshot at the reference epoch. These variations have consequences for both the stability and observability of CBPs (see Section 4.4). Dynamical studies of CBPs indicate that a critical regime exists, such that CBPs with periods less than $P_{crit}$ are unstable, predominantly scattering onto an unbound orbit, or occasionally colliding with either star (Dvorak 1986; Holman & Wiegert 1999; Sutherland & Fabrycky 2016; Lam & Kipping 2018; Quarles et al. 2018). The process by which this instability occurs is resonant overlap (Mudryk & Wu 2006; Sutherland & Kratter 2019). The value of $P_{crit}$ is primarily a function of the dynamical mass ratio of the host stars $\mu = M_2/(M_1 + M_2)$, binary period $P_{bin}$, and eccentricity $e_{bin}$, but it also has a dependence on the mutual inclination between the two orbits, $\Delta I$.

To investigate the dynamical evolution of the TOI-1338 system, we integrated the orbit of the CBP for $\sim 10^5$ yr (corresponding to $\sim 7 \times 10^3$ orbits of the binary) using the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The transits of the CBP across the primary in Sector 3 (left panels), Sector 6 (middle panels), and Sector 10 (right panels) and the best-fitting model. The Sector 6 and 10 events were fit using 2 minute cadence data, but for clarity, we show the data binned to the 30 minute cadence.</figcaption></figure>

best-fit photometric-dynamical solution in Table 4 and the IAS15 integrator in the REBOUND integration package (Rein & Liu 2012; Rein & Spiegel 2015). Figure 13 shows the evolution of the system for 40,000 yr, using as the initial condition the best-fit photometric-dynamical solution from Table 5. As seen from the figure, both the CBP (black) and binary (red) semimajor axes are practically constant with time, indicating that the system is stable. Over the course of these integrations, the eccentricity of CBP varies in a small range from 0.0695 to 0.1763 (Figure 13(b)).

To explore whether the planet's eccentricity will continue to increase and if that will affect the stability of its orbit, we used a modified version of the *mercury6* integration package (Chambers et al. 2002) and integrated the system for $10^7$ yr. Our results showed that the extrema for the planetary eccentricity extend by an additional, but insignificant, amount of 0.0004, confirming that the orbit of the CBP is long-term stable. Figure 14 demonstrates a more global range of stability using the *mercury6* integrator for binaries, tracking the extrema of planetary eccentricity (Dvorak et al. 2004; Ramos et al. 2015). As shown here, the orbit of the CBP (green dot) lies between 6:1 and 7:1 mean-motion resonances (MMRs; downward ticks, top axis) with the binary. This is important for long-term orbital stability and is well below the eccentricity-dependent stability limit (dashed line; Quarles et al. 2018). This is an expected result and a consequence of the fact that CBPs form at large distances away from the binary and migrate to their current orbits (e.g., Pierens & Nelson 2013; Kley & Haghighipour 2015). Those that maintain stable orbits are trapped between two *N*:1 MMRs with the binary. This has indeed been observed in all Kepler CBPs, as it is critical for long-term orbital stability. Figure 14 also shows that, although the orbit of the CBP is stable, small changes in its semimajor axis or eccentricity may result in a more chaotic orbit by

situating the planet near a region of instability corresponding to *N*:1 MMRs with the binary.

As an independent test to examine the stability of the planet, we used the results of our numerical integrations in the context of the scheme developed by Quarles et al. (2018) and identified a region around the binary where the orbit of the CBP will certainly be unstable. Our analysis shows that the outer boundary of this unstable region corresponds to $P_{crit} = 64.3$ days ($a_{crit} = 0.36$ au). The observed planetary period is ≈50% longer than this critical value,<fnref n="53" /> once again confirming that the orbit of the CBP is stable. Additionally, we also used a frequency analysis (Laskar 1993) to obtain a quasi-periodic decomposition of the orbital perturbations of the CBP. We found these to be a combination of the five fundamental frequencies—i.e., the mean motions of the orbit of the binary, the CBP, the apsidal precession of the binary and CBP orbits, and the nodal precession—and fully consistent with the numerical simulations.

Our numerical simulations also show indications of both apsidal and nodal precessions in the orbit of the CBP. Figures 13(c) and (d) show the *x*-components of the planet's eccentricity (apsidal) and inclination (nodal) vectors. As seen here, many secular precession cycles of the planet occur within the span of 40,000 yr. The figures show a mode with a ~14,286 yr period and variations that occur on a much shorter timescale of decades. We use the Fast Fourier Transform routine within `scipy` to produce the periodograms shown in Figures 13(e) and (f), where the system was evolved for 100,000 yr. These periodograms show strong peaks at ~23 yr (8375 days) for the planetary apsidal precession period and 21.4 yr (7816 days) for the planetary nodal precession period. This nodal precession

<footnote n="53">53 The planetary semimajor axis is ≈30% larger than the critical semimajor axis.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The same as Figure 7 but for a model with a dilution factor for the long-cadence data as an additional parameter to address the trend seen in that figure. While this model suggests that there is a dilution of ~3.4% between the long- and short-cadence data, the best-fit parameters of the system do not change significantly.</figcaption></figure>

period differs slightly from the analytical result (8980 days) derived from formula given in Farago et al. (2009) since the stellar binary’s orbit is noncircular.

Similar to the CBP, the orbit of the binary also experiences nodal precession. This precession is predominantly due to the perturbation of the CBP, with the tidal precession and general relativistic effect being secondary. The planetary apsidal and nodal precessions occur with similar periods but in opposite directions, as expected from Lee & Peale (2006). The longer-period mode (~14,286 yr) for the planet is in phase with the binary secular precession period, where the binary causes an oscillation in the planetary argument of pericenter by ±28° per binary precession cycle (Andrade-Ines & Robutel 2018).

## 4.4. Transit and Detection Probabilities

The nodal precession of the planet’s orbit has important consequences for the long-term detectability of its transits, as these can only occur when the projected path of the planet on the plane of the sky intersects with that of the stars. Nodal precession alters the planet’s path, even to the extent that

transits disappear for long periods of time. This was predicted by Schneider (1994) and observationally confirmed by the transits of the CBP Kepler-413b (Kostov et al. 2014). In this context, Martin & Triaud (2014) discussed the “transitability” CBPs and found that transits of CBPs could also occur in non-eclipsing binaries.

In Figure 15, we show how the impact parameter of the planet, and therefore, its transitability, varies over time due to the orbital evolution of the planet. As seen from the figure, there are two windows of transitability per nodal precession cycle, each roughly ~1000 days wide as indicated by the points between the horizontal gray lines in the range -1 ≤ b ≤ 1. Figure 15(b) illustrates the motion of the planet on the sky plane for the 2000 days after the starting epoch, where the points (color-coded) are spaced by ~0.7 days and vary in size (larger, opaque is toward the observer). The horizontal ellipses represent the orbits of the primary (black) and secondary (gray) stars with respect to the center of mass at the origin. The vertical ellipses indicate the cross section that each stellar disk takes up on the sky. When the larger, opaque points overlap
<figure><img src="image_13.png" /><figcaption>**Figure 10.** The residuals for the radial velocities measured by CORALIE (top and middle panels, black points) and HARPS (bottom panel, black points) for the best-fitting photodynamical model. The red and blue points show the residuals of the HARPS measurements for models with the planet mass fixed at 120 $M_⊕$ and 150 $M_⊕$, respectively.</figcaption></figure>

<table><caption>**Table 5**
Initial Dynamical Parameters<sup>a</sup></caption><thead><tr><th>Parameter<sup>b</sup></th><th>Binary Orbit</th><th>Planet Orbit</th></tr></thead><tbody><tr><td>Period (days)</td><td>1.46085607280931704E+01</td><td>9.51407742682822573E+01</td></tr><tr><td>e cos ω</td><td>-7.21833507304852212E-02</td><td>-1.07723606253091915E-02</td></tr><tr><td>e sin ω</td><td>1.38301268444568498E-01</td><td>-9.22020142859296676E-02</td></tr><tr><td>i (rad)</td><td>1.56482063218343082E+00</td><td>1.55715854420219579E+00</td></tr><tr><td>Ω (rad)</td><td>0.0000000000000000E+00</td><td>1.52632127894949884E-02</td></tr><tr><td>T<sub>conj</sub> (days)<sup>c</sup></td><td>3.33682449014051008E+03</td><td>3.34142761690665247E+03</td></tr><tr><td>a (au)</td><td>1.28783121829547487E-01</td><td>4.49132971740733966E-01</td></tr><tr><td>e</td><td>1.56005374830665650E-01</td><td>9.28291720948987292E-02</td></tr><tr><td>ω (deg)</td><td>1.17561331987513597E+02</td><td>2.63336097695842795E+02</td></tr><tr><td>True anomaly (deg)</td><td>1.11217395295404202E+02</td><td>1.36038415640448932E+02</td></tr><tr><td>Mean anomaly (deg)</td><td>9.38816596366525289E+01</td><td>1.28272133285056185E+02</td></tr><tr><td>Mean longitude (deg)</td><td>2.28778727282917799E+02</td><td>4.00249031010939916E+02</td></tr><tr><td>i (deg)</td><td>8.96576179191039415E+01</td><td>8.92186126155212662E+01</td></tr><tr><td>Ω (deg)</td><td>0.0000000000000000E+00</td><td>8.74517674648163101E-01</td></tr></tbody></table><table><thead><tr><th>Parameter<sup>d</sup></th><th>Body 1</th><th>Body 2</th><th>Body 3</th></tr></thead><tbody><tr><td>Mass (M<sub>☉</sub>)</td><td>1.03784970719363567E+00</td><td>2.97388770751337850E-01</td><td>9.06017229632760055E-05</td></tr><tr><td>x (au)</td><td>1.95196590778876217E-02</td><td>-6.82335050435495666E-02</td><td>3.68709659800795730E-01</td></tr><tr><td>y (au)</td><td>1.32648896371446897E-04</td><td>-4.65900424342245798E-04</td><td>9.75628498103004414E-03</td></tr><tr><td>z (au)</td><td>2.22880532978689296E-02</td><td>-7.78747286204819755E-02</td><td>3.02645766911772141E-01</td></tr><tr><td>v<sub>x</sub> (au day<sup>-1</sup>)</td><td>-7.66586168788931707E-03</td><td>2.67578221143398472E-02</td><td>-1.61533175726252566E-02</td></tr><tr><td>v<sub>y</sub> (au day<sup>-1</sup>)</td><td>5.45640776560390306E-05</td><td>-1.90441050119945319E-04</td><td>6.31089229491667655E-05</td></tr><tr><td>v<sub>z</sub> (au day<sup>-1</sup>)</td><td>9.13006908476301365E-03</td><td>-3.18697187907663951E-02</td><td>2.27034214206392367E-02</td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> Reference time = 186.00000, integration step size = 0.05000 days.</footnote>

<footnote><sup>b</sup> Jacobian instantaneous (Keplerian) elements.</footnote>

<footnote><sup>c</sup> Times are relative to BJD 2,455,000.000.</footnote>

<footnote><sup>d</sup> Barycentric Cartesian coordinates.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Top panel: posterior distribution of the planet mass in units of $M_⊕$. The median value is 33.0 $M_⊕$, and the largest value is 117.05 $M_⊕$. Bottom panel: $\chi^2 - \chi^2_{min}$ vs. the planet mass in $M_⊕$ for all computed models.</figcaption></figure>

with the vertical ellipse, then transits are possible. As seen from Figure 15(a), the planet transits 29.7% of the time. This is valid both for the best-fit solution from the posterior of the photometric-dynamical analysis (see Section 4) and for the solutions taken from the overall sample from the posterior and propagated for 155,000 days (~20 nodal precession cycles), where the planet transits 29.6% ± 7.6% of the time. This is typical of the Kepler CBPs, for which the mean primary transitability across the first 10 discovered planets was ~36% (Martin 2017).

### 4.5. Future Transits

As an aid to enable further observations of the TOI-1338b transits, we present in Table 6 the predictions of the times, impact parameters, and durations of future transit events. These three quantities were computed using 9000 models from the posterior sample. The quoted values are the sample medians, and the quoted uncertainties are the sample rms. Transits will certainly occur on 2020 January 12, April 14, July 19, and October 19 since all 9000 models from the posterior had transits at these dates. Starting on 2021 April 26, not all models

from the posterior produce transits at that time, so the transits become less likely. At first, the fraction of missed transits is rather small (a few percent), but then, starting with the conjunction on 2023 May 20, the fraction of missed transits is large (~30%) and grows larger and larger thereafter. After the 2025 September 12 conjunction, the transits have a ≲30% chance of occurring, and if they do occur, their impact parameters will likely be close to -1.

The primary TESS mission ends in 2020 July, but fortunately, the mission has been extended through at least the year 2022.<fnref n="54" /> Depending on the exact pointing schedule, there is a good chance that TESS can observe transits again on 2020 October 19, 2021 January 23, 2021 April 26, and possibly July 30.

For the sake of completeness, we also calculate the extent of the HZ (Kasting et al. 1993) of the binary (Figure 16) using the Multiple Star HZ website developed by Müller & Haghighipour (2014). The inner boundary of orbital stability is shown in red, and the orbit of the CBP is shown in blue. TOI-1338 b is

<footnote n="54">⁵⁴ https://tess.mit.edu/news/nasa-extends-the-tess-mission/</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** A common period O–C diagram that shows the ETVs in minutes for the model primary eclipses (black lines) and the model secondary eclipse (red lines). The two inner lines are for the overall best-fitting model with $M_3 = 33 M_\odot$, and the two outer lines are for the best-fitting model with the planet mass fixed at $M_3 = 150 M_\odot$. The black point and error bar show the median uncertainty of 0.36 minutes for the measured times of the primary eclipses, and the red point and error bar show the median uncertainty of 5.40 minutes for the measured times of the secondary eclipses. TESS observed the system between days ≈3336 and ≈3634.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Evolution of the orbital elements of the CBP using the best-fit parameters from Table 5 for 40,000 yr (BJD—2,455,000). Panel (a) shows small variations in the semimajor axis of the planet (black) and binary (red). Panel (b) shows the variations of orbital eccentricities, indicating that the system is stable. The evolution of the x-components of the eccentricity, panel (c), and inclination vectors, panel (d), illustrate apsidal and nodal precession, respectively. Panels (e) and (f) show Fourier periodograms using a 100,000 yr simulation, where the peak values are the planetary apsidal and nodal precession periods for the respective vectors.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Stability map exploring the range of planetary eccentricity oscillation ($e_{max} - e_{min}$; color-coded) considering a wide range of initial values in planetary semimajor axis $a_p$ (in au) and eccentricity $e_p$. The white cells indicate unstable initial conditions on a 100,000 yr timescale, where the vertical dips denote the locations of *N*: 1 mean-motion resonances (MMRs) with the inner binary. The current planetary parameters are indicated by a green dot, where the planet lies approximately midway between the 6:1 and 7:1 MMRs. The dashed curve shows the boundary of stability from Quarles et al. (2018), where the top axis marks the planetary period (in days) along with the location of the *N*:1 MMRs.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Panel (a): evolution of the impact parameter due to nodal precession of the planetary orbit for 20,000 days (BJD—2,455,000). Panel (b): motion of the planet on the sky plane, where each point (color-coded) represents a small increment in the planetary orbit (~0.7 days). The larger, opaque points indicate when the planet lies between the binary and observer, while the smaller, translucent points indicate when the planet lies behind the binary orbit. The horizontal ellipses show the orbits of the primary (black) and secondary (gray) star about the center of mass, while the vertical ellipses illustrate the cross section for transits to occur across each stellar disk.</figcaption></figure>
<table><caption>**Table 6**
Times, Durations, and Impact Parameters of Future Planet Transits</caption><thead><tr><th>BJD-2,455,000</th><th>Year</th><th>Month</th><th>Day</th><th>UTC</th><th>Impact Parameter</th><th>Duration (hr)</th><th>Transit Fraction</th></tr></thead><tbody><tr><td>3861.1845 ± 0.0161</td><td>2020</td><td>Jan</td><td>12</td><td>16:25:40.7</td><td>0.143 ± 0.128</td><td>6.74 ± 0.16</td><td>100.0%</td></tr><tr><td>3954.1353 ± 0.0214</td><td>2020</td><td>Apr</td><td>14</td><td>15:14:50.2</td><td>0.206 ± 0.108</td><td>9.80 ± 0.23</td><td>100.0%</td></tr><tr><td>4049.6930 ± 0.0196</td><td>2020</td><td>Jul</td><td>19</td><td>04:37:53.2</td><td>-0.006 ± 0.189</td><td>7.91 ± 0.25</td><td>100.0%</td></tr><tr><td>4142.3991 ± 0.0224</td><td>2020</td><td>Oct</td><td>19</td><td>21:34:40.9</td><td>0.055 ± 0.166</td><td>7.46 ± 0.18</td><td>100.0%</td></tr><tr><td>4237.9220 ± 0.0236</td><td>2021</td><td>Jan</td><td>23</td><td>10:07:40.6</td><td>-0.131 ± 0.236</td><td>9.73 ± 0.63</td><td>99.9%</td></tr><tr><td>4330.9984 ± 0.0264</td><td>2021</td><td>Apr</td><td>26</td><td>11:57:42.7</td><td>-0.128 ± 0.232</td><td>6.31 ± 0.38</td><td>99.9%</td></tr><tr><td>4425.7194 ± 0.0319</td><td>2021</td><td>Jul</td><td>30</td><td>05:15:54.7</td><td>-0.224 ± 0.266</td><td>11.78 ± 1.20</td><td>99.6%</td></tr><tr><td>4519.6951 ± 0.0328</td><td>2021</td><td>Nov</td><td>1</td><td>04:40:55.7</td><td>-0.315 ± 0.285</td><td>5.94 ± 0.79</td><td>98.7%</td></tr><tr><td>4708.3278 ± 0.0389</td><td>2022</td><td>May</td><td>8</td><td>19:52:05.0</td><td>-0.468 ± 0.308</td><td>6.06 ± 1.27</td><td>92.9%</td></tr><tr><td>4801.0174 ± 0.0461</td><td>2022</td><td>Aug</td><td>9</td><td>12:24:59.3</td><td>-0.417 ± 0.304</td><td>8.44 ± 1.58</td><td>95.5%</td></tr><tr><td>4896.7585 ± 0.0461</td><td>2022</td><td>Nov</td><td>13</td><td>06:12:14.7</td><td>-0.566 ± 0.310</td><td>6.78 ± 1.72</td><td>82.4%</td></tr><tr><td>4989.3540 ± 0.0509</td><td>2023</td><td>Feb</td><td>13</td><td>20:29:49.0</td><td>-0.540 ± 0.313</td><td>6.04 ± 1.46</td><td>85.2%</td></tr><tr><td>5084.8165 ± 0.0693</td><td>2023</td><td>May</td><td>20</td><td>07:35:48.3</td><td>-0.624 ± 0.310</td><td>8.19 ± 2.34</td><td>72.6%</td></tr><tr><td>5177.9790 ± 0.0600</td><td>2023</td><td>Aug</td><td>21</td><td>11:29:44.2</td><td>-0.633 ± 0.317</td><td>4.93 ± 1.45</td><td>70.0%</td></tr><tr><td>5272.3566 ± 0.1296</td><td>2023</td><td>Nov</td><td>23</td><td>20:33:30.0</td><td>-0.650 ± 0.311</td><td>9.79 ± 2.97</td><td>64.5%</td></tr><tr><td>5366.6770 ± 0.0742</td><td>2024</td><td>Feb</td><td>26</td><td>04:14:52.4</td><td>-0.670 ± 0.318</td><td>4.79 ± 1.52</td><td>54.2%</td></tr><tr><td>5459.7304 ± 0.1788</td><td>2024</td><td>May</td><td>29</td><td>05:31:47.0</td><td>-0.665 ± 0.315</td><td>8.76 ± 2.76</td><td>55.7%</td></tr><tr><td>5555.2836 ± 0.0967</td><td>2024</td><td>Sep</td><td>1</td><td>18:48:23.7</td><td>-0.690 ± 0.320</td><td>5.26 ± 1.77</td><td>42.9%</td></tr><tr><td>5647.6474 ± 0.1677</td><td>2024</td><td>Dec</td><td>3</td><td>03:32:17.7</td><td>-0.687 ± 0.322</td><td>6.27 ± 2.08</td><td>44.5%</td></tr><tr><td>5743.6224 ± 0.1383</td><td>2025</td><td>Mar</td><td>9</td><td>02:56:12.7</td><td>-0.707 ± 0.320</td><td>6.35 ± 2.21</td><td>35.5%</td></tr><tr><td>5836.0445 ± 0.1543</td><td>2025</td><td>Jun</td><td>9</td><td>13:04:05.3</td><td>-0.700 ± 0.329</td><td>4.93 ± 1.67</td><td>34.2%</td></tr><tr><td>5931.4246 ± 0.2222</td><td>2025</td><td>Sep</td><td>12</td><td>22:11:25.1</td><td>-0.713 ± 0.320</td><td>8.39 ± 3.01</td><td>31.0%</td></tr><tr><td>6024.6162 ± 0.1615</td><td>2025</td><td>Dec</td><td>15</td><td>02:47:16.6</td><td>-0.694 ± 0.333</td><td>4.71 ± 1.64</td><td>27.7%</td></tr><tr><td>6118.5043 ± 0.3348</td><td>2026</td><td>Mar</td><td>19</td><td>00:06:09.6</td><td>-0.702 ± 0.325</td><td>9.86 ± 3.45</td><td>27.4%</td></tr><tr><td>6213.1325 ± 0.1991</td><td>2026</td><td>Jun</td><td>21</td><td>15:10:46.2</td><td>-0.687 ± 0.333</td><td>5.22 ± 1.81</td><td>23.9%</td></tr><tr><td>6305.7169 ± 0.3200</td><td>2026</td><td>Sep</td><td>22</td><td>05:12:19.8</td><td>-0.698 ± 0.333</td><td>7.40 ± 2.63</td><td>23.8%</td></tr><tr><td>6401.3824 ± 0.2847</td><td>2026</td><td>Dec</td><td>26</td><td>21:10:36.5</td><td>-0.685 ± 0.332</td><td>6.52 ± 2.34</td><td>22.4%</td></tr><tr><td>6493.7147 ± 0.2750</td><td>2027</td><td>Mar</td><td>29</td><td>05:09:06.7</td><td>-0.673 ± 0.340</td><td>5.53 ± 1.88</td><td>20.7%</td></tr><tr><td>6589.0498 ± 0.4693</td><td>2027</td><td>Jul</td><td>2</td><td>13:11:42.7</td><td>-0.678 ± 0.329</td><td>8.97 ± 3.27</td><td>21.5%</td></tr><tr><td>6682.0668 ± 0.2875</td><td>2027</td><td>Oct</td><td>3</td><td>13:36:11.6</td><td>-0.665 ± 0.342</td><td>5.11 ± 1.70</td><td>19.6%</td></tr><tr><td>6775.8467 ± 0.7051</td><td>2028</td><td>Jan</td><td>5</td><td>08:19:17.2</td><td>-0.675 ± 0.332</td><td>10.95 ± 3.69</td><td>20.4%</td></tr><tr><td>6870.3964 ± 0.3605</td><td>2028</td><td>Apr</td><td>8</td><td>21:30:48.3</td><td>-0.672 ± 0.345</td><td>5.65 ± 1.98</td><td>20.7%</td></tr><tr><td>6962.8655 ± 0.5975</td><td>2028</td><td>Jul</td><td>10</td><td>08:46:15.3</td><td>-0.666 ± 0.340</td><td>7.74 ± 2.61</td><td>19.6%</td></tr><tr><td>7058.3711 ± 0.5322</td><td>2028</td><td>Oct</td><td>13</td><td>20:54:25.6</td><td>-0.672 ± 0.341</td><td>7.42 ± 2.70</td><td>23.1%</td></tr><tr><td>7150.8050 ± 0.4580</td><td>2029</td><td>Jan</td><td>14</td><td>07:19:14.6</td><td>-0.661 ± 0.351</td><td>5.64 ± 1.88</td><td>21.1%</td></tr><tr><td>7245.5142 ± 0.8968</td><td>2029</td><td>Apr</td><td>19</td><td>00:20:24.3</td><td>-0.682 ± 0.341</td><td>10.53 ± 4.07</td><td>26.0%</td></tr><tr><td>7339.0599 ± 0.4556</td><td>2029</td><td>Jul</td><td>21</td><td>13:26:12.1</td><td>-0.664 ± 0.361</td><td>5.33 ± 1.87</td><td>25.6%</td></tr><tr><td>7431.8295 ± 1.0387</td><td>2029</td><td>Oct</td><td>22</td><td>07:54:30.9</td><td>-0.685 ± 0.351</td><td>10.20 ± 3.65</td><td>27.5%</td></tr><tr><td>7527.1701 ± 0.5751</td><td>2030</td><td>Jan</td><td>25</td><td>16:04:55.6</td><td>-0.696 ± 0.361</td><td>6.19 ± 2.34</td><td>36.4%</td></tr></tbody></table>

substantially interior to the HZ, receiving more than nine times the Sun–Earth insolation.

# 5. Discussion

## 5.1. Comparison with Stellar Evolution Models

We compared the best-fit stellar masses, sizes, and temperatures of the primary and secondary stars of TOI-1338 against model isochrones from the MIST series (Choi et al. 2016) for the measured metallicity of the system. The fitted masses and sizes of both stars are in excellent agreement with a 4.4 Gyr isochrone model (see Figure 17). It is interesting to note that the M-dwarf secondary of TOI-1338 does not seem to be significantly inflated compared to standard models, as may be the case for some CBP systems with similar secondary stars (e.g., Kepler-38 and Kepler-47, Orosz et al. 2012b, 2019). However, this is in line with the rest of the objects detected during the EBLM Project, which show no systematic radius inflation for fully convective low-mass stars (von Boetticher et al. 2019). The effective temperature of the secondary star is also consistent with model predictions, within the errors.

## 5.2. Stellar Rotation

Starspots create modulations in the light curves of eclipsing binaries (typically seen in the out-of-eclipse regions), which affect measurements of eclipse times and, thus, photodynamical models (see, e.g., Kepler-453, Welsh et al. 2015). The systematic errors present in the TESS light curve of TOI-1338 (see Figure 1) preclude measurement of the intrinsic stellar variability from the TESS data alone, which, in turn, precludes determination of the stellar spin period based on starspot modulations.

However, if we assume that the stellar rotation axis is approximately perpendicular to the line of sight, then our measurement of $v \sin i_*$ combined with the stellar radius estimated above imply a rotation period for the primary star of $P_{rot} = 19 \pm 3$ days. Given that the binary orbit is eccentric and that the timescale for synchronization ($\sim 2.5$ Gyr) is comparable to the estimated age ($\sim 4.4$ Gyr), we can expect the rotation period to be closer to the pseudosynchronous period ($\sim 12.7$ days). If this star is magnetically active, then we expect to see modulation of the light curve at frequencies $1/P_{rot}$
<figure><img src="image_19.png" /><figcaption>**Figure 16.** The extent of the system HZ, calculated based on the properties shown in Table 4. The conservative HZ is shown in dark green, and the optimistic extension to the HZ is shown in light green. The known planet is substantially interior to the boundaries of the system HZ and close to the stability boundary.</figcaption></figure>

and/or $2/P_{\text{rot}}$, depending on the distribution of active regions on the stellar surface at the time of observation. To search for such modulations and periodic signals in the WASP light curve of TOI-1338, we used the sine-wave fitting method described in Maxted et al. (2011). The WASP light curve contains 26,492 observations obtained with the same CCD camera and 200 mm lens over three observing seasons. We calculated the periodogram over 32,768 uniformly spaced frequencies from 0 to $1.5 \text{ cycles day}^{-1}$. The false-alarm probability (FAP) is calculated using a boot-strap Monte Carlo method also described in Maxted et al. (2011). The periodogram is shown in Figure 18. From the boot-strap Monte Carlo simulations and the lack of any significant signal in this periodogram, we can place an upper limit of approximately 1 mmag on the semi-amplitude of any signal due to rotational modulation. A similar analysis for the three seasons of WASP data separately is consistent with this conclusion.

To further this point, we also obtained historical data from the All-Sky Automated Survey for Supernovae (ASAS-SN) project (Shappee et al. 2014; Kochanek et al. 2017) in order to search for evidence of starspot modulation. Our analysis indicates that ~1600 days of ASAS-SN V-band data shows no significant photometric modulation either (Figure 19). Additionally, the HARPS spectra have a dispersion of only a few m s⁻¹, compatible with a chromospherically quiet star. The bisectors of the cross-correlation function likewise show no variability (Queloz et al. 2001). These considerations further strengthen our assumption that the stellar activity and starspot-induced bias in the measured eclipse times is negligible and that their effect on the photodynamical solution is minimal. If

additional photometric observations should reveal the primary star to be heavily spotted, then the planet mass determination may need to be revised.

### 5.3. The Planet

With a mass of $33.0 \pm 20.0 M_{\oplus}$, a radius of $6.85 \pm 0.19 R_{\oplus}$, and a bulk density of $0.56 \pm 0.34 \text{ g cm}^{-3}$ (Table 4, where the quoted values are from the sample medians), the closest solar system analog of TOI-1338 b is perhaps Saturn, where $M = 95.16 M_{\oplus}$, $R = 9.14 R_{\oplus}$, and $\rho = 0.69 \text{ g cm}^{-3}$. Among the known CBPs, TOI-1338b has bulk properties similar to those of Kepler-16b (Doyle et al. 2011) and Kepler-34b (Welsh et al. 2012).

The large radius of TOI-1338 b is consistent with the predictions of planet-formation models in circumbinary disks, a trend that has also been observed among CBPs detected by the Kepler telescope. Combined with the small orbital inclination, this indicates that TOI-1338 b formed at larger distances from its host binary and migrated to its current orbit through planet-disk interaction (e.g., Pierens & Nelson 2013; Kley & Haghighipour 2015).

### 5.4. TOI-1338 within the Context of the Kepler CBP Systems

The TOI-1338 system follows the trends established by Kepler CBPs. Namely, these are gas giant planets (radii larger than $3R_{\oplus}$) with low-eccentricity, with nearly coplanar orbits<fnref n="55" /> with periods longer than ~50 days, and that orbit around binary stars with $P_{\text{bin}} \sim 7.5-40$ days (Welsh & Orosz 2018;

<footnote n="55">⁵⁵ $e_{\text{CBP}} < 0.15$, mutual orbital inclination $\Delta i < 5^{\circ}$.</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Upper panel: mass and radius of TOI-1338 compared against MIST isochrones. The dotted lines represent isochrones for ages between 2 and 6 Gyr in 1 Gyr increments, and the solid line is the best-fit isochrone corresponding to 4.4 Gyr. The inset shows an enlargement around the location of the secondary star, whose properties are seen to be consistent with the best-fit model. Lower panel: the same as the upper panel but for mass and temperature.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Periodogram of the WASP light curve for TIC 260128333. The dashed line indicates a false-alarm probability level of 0.01.</figcaption></figure>

Martin et al. 2019, and references therein).<fnref n="56" /> The near coplanarity of the planetary orbits is also consistent with the observational results of circumbinary disks around short-period

stellar binaries (Czekala et al. 2019). TOI-1338 is similar to the Kepler-38 CBP system (Orosz et al. 2012a) in terms of both the orbital periods and orbital period ratio (see Figure 20). The orbital precession timescale of TOI-1338 b is comparable to that of the CBP Kepler-413b, where the timescale is ~22 yr for the former compared to ~11 yr for the latter (Kostov et al. 2014).

<footnote n="56">⁵⁶ Whereas most of Kepler's EBs have orbital periods shorter than ~3 days (Martin et al. 2015; Muñoz & Lai 2015; Hamers et al. 2016; Fleming et al. 2018).</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Approximately 1600 days of ASAS-SN V-band photometry shows no significant rotationally induced photometric modulation.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Panel (a): ratio of the planet and binary orbital periods for Kepler CBPs (blue squares) and TESS (purple triangle). The diagonal black dashed lines indicate constant period ratios. Panel (b): corresponding ratio of the planet periapse and the binary apoapse, which together represent the shortest separation between the planet and binary.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 21.** Radius and period of the TESS CBP (purple triangle) and Kepler CBPs (blue squares) compared with the TESS planet candidates (yellow triangles), Kepler planet candidates (green circles), and the five innermost solar system planets (red diamonds). This circumbinary discovery is the longest-period confirmed TESS planet to date and exists in a parameter space similar to the Kepler CBPs.</figcaption></figure>
Overall, TOI-1338b has a relatively long period for a transiting planet, particularly when compared with other TESS candidates. This is demonstrated in Figure 21, where it resides on the very tail of the TESS planet candidate period distribution. We note that the current lack of small CBPs is likely an observational bias, since unique challenges have inhibited their detection to date, largely as a result of the transit timing variations induced by the barycentric binary motion and the orbital dynamics (Armstrong et al. 2014; Windemuth et al. 2019). The CBP population is yet to be constrained below $4R_⊕$ (Armstrong et al. 2014), and we expect that the large quantity and brightness of the TESS stars will enable the expansion of this parameter space.

## 6. Conclusions

We presented the discovery of the first transiting CBP from TESS, TOI-1338. The target was observed by TESS in 12 sectors of 30 minute cadence data (Sectors 1 through 12), and nine sectors of 2 minute cadence data (Sectors 4 through 12). In addition to stellar eclipses, three transit events were observed in Sectors 3, 6, and 10. These extra transit events show the hallmark characteristics of a circumbinary object where their duration depends on the binary phase and their times have significant deviations from a simple linear ephemeris. Blending is not an issue with TOI-1338, as the nearest source is $53''$ away, and speckle imaging observations from SOAR rule out nearby sources with a magnitude difference of $\Delta I \sim 4$ down to $0''/5$ from the target. Radial-velocity measurements are available, as the host eclipsing binary has been monitored for more than 3 yr by CORALIE and HARPS as part of the EBLM project. To solve for the parameters of the system, we combined the TESS data with the radial velocities into the photometric-dynamical model ELC. Our analysis confirms that the circumbinary object is indeed a planet, with a mass of $33.0 \pm 20.0 M_⊕$, a radius of $6.85 \pm 0.19 R_⊕$, and a bulk density of $0.56 \pm 0.34 \text{ g cm}^{-3}$. The planet's orbit is within $\sim 1^\circ$ of being coplanar with the binary, has a period of 95.2 days and small eccentricity, and is safely beyond the boundary for stability. The host eclipsing binary (with $P = 14.6$ days and $e \approx 0.16$) consists of G+M stars with masses $1.1 M_⊙$ and $0.3 M_⊙$, and radii of $1.3 R_⊙$ and $0.3 R_⊙$, respectively. Based on the stellar parameters, we estimate an age of 4.4 Gyr for the system.

This manuscript includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by NASA Science Mission directorate. We acknowledge the use of public TESS data from pipeline at the TESS Science Processing Operations Center. The manuscript includes data from CORALIE, an instrument mounted on the Euler 1.2 m telescope, a project of the University of Geneva, funded by the Swiss National Science Foundation. Furthermore, our analysis includes spectra obtained with HARPS, an instrument mounted on the ESO 3.6 m telescope at La Silla. Those data were obtained under Prog.ID 1101.C-0721 (PI Triaud). They are or will become available through the ESO public archive. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. This research was supported in part through research

cyberinfrastructure resources and services provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology. This research made use of Lightkurve, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration et al. 2018). This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. We are also grateful to the observer support staff at CTIO, ESO/HARPS, and Swiss Euler Telescope/CORALIE. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement. W.F.W. and J.A.O. thank John Hood Jr. for his generous support of exoplanet research at SDSU. Support was also provided and acknowledged through NASA Habitable Worlds grant 80NSSC17K0741 and NASA XRP grant 80NSSC18K0519. This work is partly supported by NASA Habitable Worlds grant 80NSSC17K0741. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under grant No. (DGE-1746045). A.H.M.J.T. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 803193/BEBOP) and from a Leverhulme Trust Research Project grant No. RPG-2018-418. A. C. acknowledges support by CfisUC strategic project (UID/FIS/04564/2019). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. E.A.G. thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining grant No. 1829740, the Brinson Foundation, and the Moore Foundation; her participation in the program has benefited this work. E.A.G. is thankful for support from GSFC Sellers Exoplanet Environments Collaboration (SEEC), which is funded by the NASA Planetary Science Divisions Internal Scientist Funding Model.

**Facility:** TESS; SOAR 4.0 m; WASP; Swiss Euler Telescope (CORALIE); ESO 3.6 m (HARPS).

**Software:** astropy (Astropy Collaboration et al. 2013); scipy (Virtanen et al. 2019); eleanor (Feinstein et al. 2019); Lightkurve (Lightkurve Collaboration et al. 2018); ispec (Blanco-Cuaresma et al. 2014); REBOUND (Rein & Liu 2012; Rein & Spiegel 2015); Mercury6 (Chambers et al. 2002).

## ORCID iDs

Veselin B. Kostov https://orcid.org/0000-0001-9786-1031
Jerome A. Orosz https://orcid.org/0000-0001-9647-2886
Adina D. Feinstein https://orcid.org/0000-0002-9464-8101
William F. Welsh https://orcid.org/0000-0003-2381-5301
Nader Haghighipour https://orcid.org/0000-0002-5234-6375
Billy Quarles https://orcid.org/0000-0002-9644-8330
David V. Martin https://orcid.org/0000-0002-7595-6360
Benjamin T. Montet https://orcid.org/0000-0001-7516-8308
Guillermo Torres https://orcid.org/0000-0002-5286-0251
Amaury H. M. J. Triaud https://orcid.org/0000-0002-5510-8751

Thomas Barclay https://orcid.org/0000-0001-7139-2724

Patricia Boyd https://orcid.org/0000-0003-0442-4284

Cesar Briceno https://orcid.org/0000-0001-7124-4094

Emily A. Gilbert https://orcid.org/0000-0002-0388-8004

Samuel Gill https://orcid.org/0000-0002-4259-0155

Michaël Gillon https://orcid.org/0000-0003-1462-7739

Jacob Haqq-Misra https://orcid.org/0000-0003-4346-2611

Coel Hellier https://orcid.org/0000-0002-3439-1439

Courtney Dressing https://orcid.org/0000-0001-8189-0233

Daniel C. Fabrycky https://orcid.org/0000-0003-3750-0183

Jon M. Jenkins https://orcid.org/0000-0002-4715-9460

Stephen R. Kane https://orcid.org/0000-0002-7084-0529

Ravi Kopparapu https://orcid.org/0000-0002-5893-2471

Vedad Kunovac Hodžić https://orcid.org/0000-0001-9419-3736

David W. Latham https://orcid.org/0000-0001-9911-7388

Nicholas Law https://orcid.org/0000-0001-9380-6457

Gongjie Li https://orcid.org/0000-0001-8308-0808

Chris Lintott https://orcid.org/0000-0001-5578-359X

Jack J. Lissauer https://orcid.org/0000-0001-6513-1659

Andrew W. Mann https://orcid.org/0000-0003-3654-1602

Pierre F. L. Maxted https://orcid.org/0000-0003-3794-1317

Nora Eisner https://orcid.org/0000-0002-9138-9028

Joshua Pepper https://orcid.org/0000-0002-3827-8417

Don Pollacco https://orcid.org/0000-0001-9850-9697

Samuel N. Quinn https://orcid.org/0000-0002-8964-8377

Jason F. Rowe https://orcid.org/0000-0002-5904-1865

S. Seager https://orcid.org/0000-0002-6892-6948

Alexandre Santerne https://orcid.org/0000-0002-3586-1316

Damien Ségransan https://orcid.org/0000-0003-2355-8034

Donald R. Short https://orcid.org/0000-0001-5504-9512

Jeffrey C. Smith https://orcid.org/0000-0002-6148-7903

Matthew R. Standing https://orcid.org/0000-0002-7608-8905

Andrei Tokovinin https://orcid.org/0000-0002-2084-0782

Trifon Trifonov https://orcid.org/0000-0002-0236-775X

Joseph D. Twicken https://orcid.org/0000-0002-6778-7552

Stéphane Udry https://orcid.org/0000-0001-7576-6236

Roland Vanderspek https://orcid.org/0000-0001-6763-6562

Joshua N. Winn https://orcid.org/0000-0002-4265-047X

Eric T. Wolf https://orcid.org/0000-0002-7188-1648

Carl Ziegler https://orcid.org/0000-0002-0619-7639

## References

Allard, F., Homeier, D., & Freytag, B. 2012, in IAU Symp. 282, From Interacting Binaries to Exoplanets: Essential Modeling Tools, ed. M. T. Richards & I. Hubeny (Cambridge: Cambridge Univ. Press), 235

Andrade-Ines, E., & Robutel, P. 2018, CeMDA, 130, 6

Armstrong, D. J., Osborn, H. P., Brown, D. J. A., et al. 2014, MNRAS, 444, 1873

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373

Blanco-Cuaresma, S., Soubiran, C., Heiter, U., & Jofré, P. 2014, A&A, 569, A111

Brandeker, A., & Cataldi, G. 2019, A&A, 621, A86

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, in IAU Symp. 135, Interstellar Dust, ed. L. J. Allamandola & A. G. G. M. Tielens (Cambridge: Cambridge Univ. Press), 5

Casagrande, L., Ramírez, I., Meléndez, J., Bessell, M., & Asplund, M. 2010, A&A, 512, A54

Chambers, J. E., Quintana, E. V., Duncan, M. J., & Lissauer, J. J. 2002, AJ, 123, 2884

Charbonneau, P. 1995, ApJS, 101, 309

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Collier Cameron, A., Wilson, D. M., West, R. G., et al. 2007, MNRAS, 380, 1230

Czekala, I., Chiang, E., Andrews, S. M., et al. 2019, ApJ, 883, 22

Doyle, A. P. 2015, PhD thesis, Keele Univ.

Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Sci, 333, 1602

Dvorak, R. 1986, A&A, 167, 379

Dvorak, R., Pilat-Lohinger, E., Schwarz, R., & Freistetter, F. 2004, A&A, 426, L37

Eastman, J. D., Rodriguez, J. E., Agol, E., et al. 2019, arXiv:1907.09480

Farago, F., Laskar, J., & Couetdic, J. 2009, CeMDA, 104, 291

Feinstein, A. D., Montet, B. T., Foreman-Mackey, D., et al. 2019, PASP, 131, 094502

Fleming, D. P., Barnes, R., Graham, D. E., Luger, R., & Quinn, T. R. 2018, ApJ, 858, 86

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Gray, R. O., & Corbally, C. J. 1994, AJ, 107, 742

Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951

Hairer, E., Lubich, C., & Wanner, G. 2002, Springer Series in Computational Mathematics, Vol. 31, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (Berlin: Springer)

Hamers, A. S., Perets, H. B., & Portegies Zwart, S. F. 2016, MNRAS, 455, 3180

Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621

Huang, C. X., Burt, J., Vanderburg, A., et al. 2018, ApJL, 868, L39

Huang, Y., Liu, X. W., Yuan, H. B., et al. 2015, MNRAS, 454, 2863

Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, Proc. SPIE, 9913, 99133E

Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icar, 101, 108

Kipping, D. M. 2013, MNRAS, 435, 2152

Kirk, B., Conroy, K., Prša, A., et al. 2016, AJ, 151, 68

Kley, W., & Haghighipour, N. 2015, A&A, 581, A20

Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502

Konacki, M., Muterspaugh, M. W., Kulkarni, S. R., & Helminiak, K. G. 2009, ApJ, 704, 513

Kostov, V. B., McCullough, P. R., Carter, J. A., et al. 2014, ApJ, 784, 14

Kostov, V. B., McCullough, P. R., Hinse, T. C., et al. 2013, ApJ, 770, 52

Kostov, V. B., Orosz, J. A., Welsh, W. F., et al. 2016, ApJ, 827, 86

Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, AJ, 158, 32

Lam, C., & Kipping, D. 2018, MNRAS, 476, 5692

Laskar, J. 1993, CeMDA, 56, 191

Lee, M. H., & Peale, S. J. 2006, Icar, 184, 573

Li, G., Holman, M. J., & Tao, M. 2016, ApJ, 831, 96

Li, J., Tenenbaum, P., Twicken, J. D., et al. 2019, PASP, 131, 024506

Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS Time Series Analysis in Python, Astrophysics Source Code Library, ascl:1812.013

López-Morales, M., Triaud, A. H. M. J., Rodler, F., et al. 2014, ApJL, 792, L31

Lovis, C., & Pepe, F. 2007, A&A, 468, 1115

Mardling, R. A., & Lin, D. N. C. 2002, ApJ, 573, 829

Martin, D. V. 2017, MNRAS, 467, 1694

Martin, D. V., Mazeh, T., & Fabrycky, D. C. 2015, MNRAS, 453, 3554

Martin, D. V., & Triaud, A. H. M. J. 2014, A&A, 570, A91

Martin, D. V., Triaud, A. H. M. J., Udry, S., et al. 2019, A&A, 624, A68

Maxted, P. F. L., Anderson, D. R., Collier Cameron, A., et al. 2011, PASP, 123, 547

Mayor, M., Udry, S., Lovis, C., et al. 2009, A&A, 493, 639

Mudryk, L. R., & Wu, Y. 2006, ApJ, 639, 423

Müller, T. W. A., & Haghighipour, N. 2014, ApJ, 782, 26

Muñoz, D. J., & Lai, D. 2015, PNAS, 112, 9264

Orosz, J. A., & Hauschildt, P. H. 2000, A&A, 364, 265

Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012a, ApJ, 758, 87

Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012b, Sci, 337, 1511

Orosz, J. A., Welsh, W. F., Haghighipour, N., et al. 2019, AJ, 157, 174

Pepe, F., Mayor, M., Rupprecht, G., et al. 2002, Msngr, 110, 9

Pierens, A., & Nelson, R. P. 2013, A&A, 556, A134

Pollacco, D. L., Skillen, I., Collier Cameron, A., et al. 2006, PASP, 118, 1407

Prša, A., Batalha, N., Slawson, R. W., et al. 2011, AJ, 141, 83

Quarles, B., Satyal, S., Kostov, V., Kaib, N., & Haghighipour, N. 2018, ApJ, 856, 150
Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279

Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1

Ramos, X. S., Correa-Otto, J. A., & Beaugé, C. 2015, CeMDA, 123, 453

Rein, H., & Liu, S.-F. 2012, A&A, 537, A128

Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, JATIS, 1, 014003

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schneider, J. 1994, P&SS, 42, 539

Schwamb, M. E., Orosz, J. A., Carter, J. A., et al. 2013, ApJ, 768, 127

Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48

Short, D. R., Orosz, J. A., Windmiller, G., & Welsh, W. F. 2018, AJ, 156, 297

Slawson, R. W., Prša, A., Welsh, W. F., et al. 2011, AJ, 142, 160

Socia, Q. J., Welsh, W. F., Orosz, J. A., et al. 2020, AJ, 159, 94

Stassun, K. G., Oelkers, R. J., Paegert, M., et al. 2019, AJ, 158, 138

Sullivan, P. W., Winn, J. N., Berta-Thompson, Z. K., et al. 2015, ApJ, 809, 77

Sutherland, A. P., & Fabrycky, D. C. 2016, ApJ, 818, 6

Sutherland, A. P., & Kratter, K. M. 2019, MNRAS, 487, 3288

Ter Braak, C. J. F. 2006, S&C, 16, 239

Tokovinin, A., Mason, B. D., Hartkopf, W. I., Mendez, R. A., & Horch, E. P. 2018, AJ, 155, 235

Tokovinin, A., Mason, B. D., Mendez, R. A., Horch, E. P., & Briceño, C. 2019, AJ, 158, 48

Triaud, A. H. M. J. 2011, PhD thesis, Observatoire Astronomique de l' Université de Genève

Triaud, A. H. M. J., Hebb, L., Anderson, D. R., et al. 2013, A&A, 549, A18

Triaud, A. H. M. J., Martin, D. V., Ségransan, D., et al. 2017a, A&A, 608, A129

Triaud, A. H. M. J., Neveu-VanMalle, M., Lendl, M., et al. 2017b, MNRAS, 467, 1714

Twicken, J. D., Catanzarite, J. H., Clarke, B. D., et al. 2018, PASP, 130, 064502

Vanderspek, R., Huang, C. X., Vanderburg, A., et al. 2019, ApJL, 871, L24

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2019, Nat. Methods, 17, 261

von Boetticher, A., Triaud, A. H. M. J., Queloz, D., et al. 2019, A&A, 625, A150

Welsh, W. F., & Orosz, J. A. 2018, in Handbook of Exoplanets, ed. H. Deeg & J. Belmonte (Cham: Springer), 34

Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Natur, 481, 475

Welsh, W. F., Orosz, J. A., Short, D. R., et al. 2015, ApJ, 809, 26

Windemuth, D., Agol, E., Carter, J., et al. 2019, MNRAS, 490, 1313