# Radial-velocity discovery of a second planet in the TOI-1338/BEBOP-1 circumbinary system

Matthew R. Standing, Lalitha Sairam, David V. Martin, Amaury H. M. J. Triaud, Alexandre C. M. Correia, Gavin A. L. Coleman, Thomas A. Baycroft, Vedad Kunovac, Isabelle Boisse, Andrew Collier Cameron, et al.

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

To cite this version:

Matthew R. Standing, Lalitha Sairam, David V. Martin, Amaury H. M. J. Triaud, Alexandre C. M. Correia, et al.. Radial-velocity discovery of a second planet in the TOI-1338/BEBOP-1 circumbinary system. Nature Astronomy, 2023, 7, pp.702-714. 10.1038/s41550-023-01948-4. insu-04505573

HAL Id: insu-04505573

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

Submitted on 4 Oct 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Radial-velocity discovery of a second planet in the TOI-1338/BEBOP-1 circumbinary system

Received: 29 July 2022

Accepted: 16 March 2023

Published online: 12 June 2023

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

Check for updates

Matthew R. Standing Ⓓ1,2✉, Lalitha Sairam¹, David V. Martin Ⓓ3,
Amaury H. M. J. Triaud Ⓓ1, Alexandre C. M. Correia Ⓓ4,5, Gavin A. L. Coleman⁶,
Thomas A. Baycroft Ⓓ1, Vedad Kunovac Ⓓ7,8, Isabelle Boisse⁹,
Andrew Collier Cameron Ⓓ10, Georgina Dransfield¹, João P. Faria Ⓓ11,12,
Michaël Gillon¹³, Nathan C. Hara¹⁴, Coel Hellier¹⁵, Jonathan Howard¹, Ellie Lane¹,
Rosemary Mardling¹⁶, Pierre F. L. Maxted¹⁵, Nicola J. Miller¹⁵, Richard P. Nelson Ⓓ6,
Jerome A. Orosz Ⓓ17, Francesco Pepe¹⁴, Alexandre Santerne Ⓓ9,
Daniel Sebastian Ⓓ1, Stéphane Udry Ⓓ14 & William F. Welsh¹⁷

Circumbinary planets, those that orbit around both stars of a central binary star system, challenge our understanding of planet formation. With only 12 binary systems known to host circumbinary planets, identifying more of these planets, along with their physical properties, could help to discern some of the physical processes that govern planet formation. Here we analyse radial-velocity data obtained by the HARPS and ESPRESSO spectrographs and report the detection of BEBOP-1 c, a gas giant planet with a mass of 65.2 ± 11.8 Earth masses ($M_⊕$) orbiting around both stars of an eclipsing binary star system with a period of 215.5 ± 3.3 days. The system TOI-1338, hereafter referred to as BEBOP-1, which also hosts the smaller and inner transiting planet TOI-1338 b, is only the second confirmed multiplanetary circumbinary system. We do not detect TOI-1338 b with radial-velocity data alone, and we can place an upper limit on its mass of 21.8 $M_⊕$ with 99% confidence. TOI-1338 b is amenable to atmospheric characterization using JWST, so the BEBOP-1 system has the potential to act as a benchmark for circumbinary exo-atmospheric studies.

Circumbinary planets are planets that orbit both stars of a central binary. They were once confined to science fiction, but the discovery of Kepler-16 b (ref. 1) paved the way for the detection of 14 transiting planets in 12 binary systems by the Kepler² and Transiting Exoplanet Survey Satellite (TESS)³ missions. Of the 12 transiting circumbinary planet systems discovered so far, only 1 hosts multiple circumbinary planets: Kepler-47. Kepler-47 b, d and c have orbital periods of 49.5, 187.4 and 303.2 days (d) respectively⁴, with Kepler-47 c's orbit placing it within the system's habitable zone. Planet-forming disks around binaries are harsh environments for planet formation to take place in, and

circumbinary planet discoveries provide insights into the formation and migration mechanisms at play in these unique environments⁵⁻⁹.

Of the 14 confirmed transiting circumbinary planets, there are only significant mass detections (different from zero at >2σ) for six of them. These are: Kepler-34 b and Kepler-35 b (ref. 10), TIC 17290098 b (ref. 11), Kepler-16 b (refs. 1,12), Kepler-47 c and Kepler-47 d (ref. 4), the masses of which were determined from binary eclipse timing variations (ETVs) alone. For the remaining eight circumbinary planets only upper limits can be placed on their masses, because the ETVs are on the order of only seconds or minutes, and hence difficult to measure.
Their masses could therefore be much lower than expected and reveal
several inflated objects, representing ideal targets for atmospheric
transmission follow-up observations¹³. To correctly characterize these
planets, accurate masses are required.

To increase the number of known circumbinary planets, and
to provide accurate masses for systems discovered with the transit
method, we initiated a radial-velocity observing survey dedicated to
circumbinary planet detection called binaries escorted by orbiting
planets (BEBOP)¹⁴. Systems followed by BEBOP are on average four
visual magnitudes brighter than circumbinary systems identified
with Kepler¹⁴. Radial velocities are less restricted to the edge-on and
shorter orbital periods found by the transit method¹⁵. Despite initial
challenges in bringing radial-velocity precision for binaries down to
values where planets can be detected¹⁶, recent results have produced
an independent detection of Kepler-16 b with a precision of 1.5 m s⁻¹,
providing a planetary mass in agreement with ETV measurements¹².

The first circumbinary planet discovered by NASA's TESS³ mission was TOI-1338 b (ref. 17). TOI-1338 is a low-mass eclipsing binary system that consists of a 1.13 solar mass (M☉) F8-type primary star with a visual magnitude of V = 11.7 (refs. 17,18) and a 0.31 M☉ M-dwarf companion. This system was already being monitored with radial velocities with HARPS¹⁹ as part of the BEBOP project, where the target was known as EBLM J0608-59 (refs. 14,18). Radial-velocity data available at the time were unable to detect any trace of TOI-1338 b, the Saturn-sized planet announced by Kostov et al.¹⁷. To constrain the mass of this planet, we used the ESPRESSO spectrograph at the Very Large Telescope²⁰. We continued to observe with HARPS intermittently to combine both datasets more easily, and to mitigate observatory closures caused by the COVID-19 pandemic. In total, we collected 123 ESPRESSO spectra and 61 HARPS spectra over a period of 1,472 d.

Results

Analysis of the radial velocities

We analysed our radial-velocity data with the diffusive nested sampler implemented in kima version 4.0.2 (ref. 21), slightly adapted for the purpose of a circumbinary system. We included non-Keplerian effects such as gravitational redshift, light time travel<sup>22-24</sup> and tidal distortion<sup>25</sup>. We also directly fitted for the most important of the Newtonian perturbations (the apsidal precession), which we parameterized using a free parameter, $\dot{\omega}_{\text{bin}}$. Further details can be found in the Methods, but we note here that one particularly useful feature of kima is that the number of planetary signals in the data can be included as a free parameter in the analysis. The number of planets detected is determined from the ratio of evidence between models with different numbers of planets. The radial-velocity data favour a one-planet model over a zero-planet model with a ratio of probability (Bayes factor) >29,000 (where >150 is the typical threshold for detection). The same analysis provided a Bayes factor of only 1.4 in favour of a two-planet model over a one-planet model, far below the detection threshold. Therefore, only one planetary signal was formally detected in our data.

Figure 1 depicts a histogram of the posterior samples obtained
from our analysis. There are two peaks, one at the binary orbital
period of 14.6 d, and a single additional excess at 215.5 d. As this
periodicity does not correspond to the 95 d period of the transiting
planet<sup>17</sup>, and an outer orbit substantially improves the fit, we concluded
that the detected planetary signal is probably that of an additional,
outer circumbinary planet within the system. We therefore gave
the system the name BEBOP-1, making it the first entry in the
BEBOP catalogue for circumbinary planets detected with radial
velocities; the new signal was named BEBOP-1 c. The name BEBOP is
recognized as an official designation by the International Astronomical
Union, and indicates the availability of high-precision radial-
velocity measurements on the target and their ability to detect a planet.
We do not refer to TOI-1338 b as BEBOP-1 b as we are not yet able to
detect its signal using radial-velocity measurements. Figure 2 shows

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

**Fig. 1 | Histogram periodogram of posterior samples obtained from a kima run on the BEBOP-1 system.** The 14.6 d binary period and 215.5 d period of BEBOP-1 c can be clearly identified (highlighted in red and blue, respectively). The inset plots show close-ups of these two peaks, where $\delta P_{\text{bin}}$ is the period of the binary orbit, and $P_c$ is the period of BEBOP-1 c's orbit. No notable peak from the transiting planet TOI-1338 b (marked with period $P_b$) can be seen in these posterior samples.

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

**Fig. 2 | Phased Keplerian radial-velocity models of the BEBOP-1 c signal.**
Models of BEBOP-1 c after removing the binary signal are shown with ESPRESSO (blue diamonds) and HARPS (orange squares) data along with their nominal error bars. The associated residuals are displayed in the bottom panel. O, observed; C, calculated. RV data are binned by 0.1 phase units (-21.6 d) weighted by their uncertainties $\sigma_i$ and illustrated by the circular data points, with their associated errors given by $\sqrt{\frac{1}{\sum i\sigma_i^2}}$. The red Keplerian models in the top panel are based on 500 randomly drawn posterior samples from a kima run, shaded from the 50th to 99th percentiles. The shaded regions display the repeating signal.

aphase-folded plot of the radial-velocity variation caused by BEBOP-1c,
with a semi-amplitude K_c = 5.6 ± 1.0 m s⁻¹, corresponding to a
minimum mass m_c sin i_c = 0.217 ± 0.035 M_jup, where i_c is the orbital
inclination of planet c and M_jup is the mass of Jupiter. The binary's
apsidal precession was found to be marginally different from zero
(ω_bin = 66.0<sup>+65.0</sup><sub>-54.7</sub> arcsec yr⁻¹). The parameters of the binary and plane-
tary orbits obtained from our model fit to the radial velocities are
available in Table 1, while Fig. 3 provides an overview of the BEBOP-1
circumbinary planet system. We calculated the true inclusion
probability<sup>26,27</sup> of our posterior samples, which also provided a clear
>99.9% probability of the presence of a planetary signal at a period
of 215.5 d. Further details can be found in the Methods.

Verification of stellar activity

Stellar magnetic activity can produce periodic modulations in radial velocity over a wide range of timescales, sometimes mimicking a planetary companion²⁸. To verify whether this is the case for the BEBOP-1 c signal, we searched for periodic signals in five spectroscopic activity indices (see Methods for details).

Stellar activity has been shown to affect most, if not all, of these five
indices at once²⁹. In our case, only two show any signal. Furthermore,
<table><caption>**Table 1** | BEBOP-1 system orbital parameters from our analysis of the ESPRESSO and HARPS radial velocities after removing outliers</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td colspan="3"><strong>Binary parameters</strong></td></tr><tr><td>Orbital period, P<sub>bin</sub></td><td>d</td><td>14.6085579(57)</td></tr><tr><td>Time of periastron passage, T<sub>0,bin</sub></td><td>BJD</td><td>9287.20017(71)</td></tr><tr><td>Semi-amplitude, K<sub>1,bin</sub></td><td>kms<sup>-1</sup></td><td>21.61764(73)</td></tr><tr><td>Eccentricity, e<sub>bin</sub></td><td></td><td>0.155522(29)</td></tr><tr><td>Argument of periastron, ω<sub>bin</sub></td><td>rad</td><td>2.05549(30)</td></tr><tr><td>Apsidal precession, ω̇<sub>bin</sub></td><td>arcsec yr<sup>-1</sup></td><td>66.0<sup>+65.0</sup><sub>-54.7</sub></td></tr><tr><td>Mass of primary, M<sub>1</sub></td><td>M<sub>☉</sub></td><td>1.127(69)<sup>a</sup></td></tr><tr><td>Radius of primary, R<sub>1</sub></td><td>R<sub>☉</sub></td><td>1.345(46)<sup>a</sup></td></tr><tr><td>Mass of secondary, M<sub>2</sub></td><td>M<sub>☉</sub></td><td>0.313(12)</td></tr><tr><td>Semi-major axis, a<sub>bin</sub></td><td>au</td><td>0.1321(25)</td></tr><tr><td>Effective temperature of primary, T<sub>eff,1</sub></td><td>K</td><td>6050(80)<sup>a</sup></td></tr><tr><td>Visual magnitude of primary, V<sub>mag,1</sub></td><td></td><td>11.72(02)<sup>a</sup></td></tr><tr><td>Spectral type</td><td></td><td>F8<sup>b</sup></td></tr><tr><td colspan="3"><strong>Planet b parameters</strong></td></tr><tr><td>Orbital period, P<sub>b</sub></td><td>d</td><td>95.174(35)<sup>a</sup></td></tr><tr><td>Semi-amplitude, K<sub>b</sub></td><td>ms<sup>-1</sup></td><td>&lt;2.4(0.1)</td></tr><tr><td>Eccentricity, e<sub>b</sub></td><td></td><td>0.0880(43)<sup>a</sup></td></tr><tr><td>Mass of planet, m<sub>b</sub></td><td>M<sub>☉</sub></td><td>&lt;21.8(0.9)</td></tr><tr><td>Mass of planet, m<sub>b</sub></td><td>M<sub>Jup</sub></td><td>&lt;0.0685(29)</td></tr><tr><td>Semi-major axis, a<sub>b</sub></td><td>au</td><td>0.4607(88)<sup>a</sup></td></tr><tr><td>Mean planet density, ρ<sub>b</sub></td><td>gcm<sup>-3</sup></td><td>&lt; 0.36<sup>+0.02</sup><sub>-0.01</sub></td></tr><tr><td colspan="3"><strong>Planet c parameters</strong></td></tr><tr><td>Orbital period, P<sub>c</sub></td><td>d</td><td>215.5(3.3)</td></tr><tr><td>Time of periastron passage, T<sub>0,c</sub></td><td>BJD</td><td>9208.7(29.1)</td></tr><tr><td>Semi-amplitude, K<sub>c</sub></td><td>ms<sup>-1</sup></td><td>5.6(1.0)</td></tr><tr><td>Eccentricity, e<sub>c</sub></td><td></td><td>&lt;0.16</td></tr><tr><td>Argument of periastron, ω<sub>c</sub></td><td>rad</td><td>3.84(86)</td></tr><tr><td>Minimum planet mass, m<sub>c</sub> sin i<sub>c</sub></td><td>M<sub>☉</sub></td><td>65.2(11.8)</td></tr><tr><td>Minimum planet mass, m<sub>c</sub> sin i<sub>c</sub></td><td>M<sub>Jup</sub></td><td>0.205(37)</td></tr><tr><td>Semi-major axis, a<sub>c</sub></td><td>au</td><td>0.794(16)</td></tr><tr><td colspan="3"><strong>System parameters</strong></td></tr><tr><td>Systemic velocity, v</td><td>kms<sup>-1</sup></td><td>30.75335(99)</td></tr><tr><td colspan="3"><strong>Instrumental</strong></td></tr><tr><td>Jitter<sub>HAR</sub></td><td>ms<sup>-1</sup></td><td>4.2<sup>+1.6</sup><sub>-3.0</sub></td></tr><tr><td>Jitter<sub>ESP19</sub></td><td>ms<sup>-1</sup></td><td>4.2(2.0)</td></tr><tr><td>Jitter<sub>ESP21</sub></td><td>ms<sup>-1</sup></td><td>4.5(0.6)</td></tr><tr><td>Radial-velocity Offset<sub>ESP19</sub></td><td>ms<sup>-1</sup></td><td>-146.4<sup>+1.9</sup><sub>-1.8</sub></td></tr><tr><td>Radial-velocity Offset<sub>ESP21</sub></td><td>ms<sup>-1</sup></td><td>-150.8(1.2)</td></tr></tbody></table>

at shorter periods neither Ca II H+K nor Hα indices produced any significant peak. Combining the stellar radius and projected rotational velocity of BEBOP-1, a rotation period of 19 ± 3 d has been reported¹⁷. We did not find any statistically significant signal at the estimated rotation period.

<figure><img src="image_6.png" /><figcaption>**Fig. 3** | Overview of the BEBOP-1 system. The BEBOP-1 system is shown along with the extent of the system's habitable zone (HZ) calculated using the Multiple Star HZ website⁹⁵. The conservative habitable zone is shown by the dark green region, while the optimistic habitable zone is shown by the light green region. The binary stars are marked by the blue star symbols in the centre. The red shaded region denotes the instability region surrounding the binary stars as described by Holman and Wiegert⁹². BEBOP-1 c's orbit is shown by the red orbit models, based on 500 randomly drawn posterior samples from a kima run, shaded from the 50th to 99th percentiles. TOI-1338 b's orbit is shown by the yellow models, and is also based on 500 random samples drawn from the posterior in its discovery paper¹⁷.</figcaption></figure>

Our analysis of activity indices ruled out any association between the 215 d radial-velocity variation and stellar activity, supporting the hypothesis of BEBOP-1 c's planetary nature. Further details are provided in the Methods and Supplementary Figs. 1–11.

## Discussion

A visual inspection of TESS light curves showed no transit of BEBOP-1 c; however, owing to orbital circulation, transits are expected to occur in due time. Circumbinary orbits exhibit nodal precession, which changes the orientation of a circumbinary planet's orbital plane with respect to both the binary and the observer. This makes a planet change from a transiting to a non-transiting configuration¹⁵,³⁰ as has been seen in a few systems³¹,³². Using an analytic criterion³³,³⁴, we found that BEBOP-1 c is guaranteed to eventually transit—mainly because the binary is so well-aligned with our line of sight ($I_{bin} = 89.658°$) combined with the relatively large size of the primary star ($R_A = 1.299 R_☉$). While BEBOP-1 c will eventually transit, we are unable to predict when and how frequently. Its precession period is of order 119 years, during which time there will be two periods of transitability of a duration that depends on BEBOP-1 c's orbital inclination.

To verify the orbital parameters of our fit we carried out a global stability analysis of the system. This stability analysis showed that the system is stable when the two circumbinary planets are in nearly circular and coplanar orbits (Methods). More precisely, the eccentricities of both planets cannot exceed 0.1, in agreement with the best-fit solution (Table 1); $i_c$ is not constrained by the radial velocities, but according to orbital stability arguments it cannot be higher than 40° with respect to the orbital plane of the binary, corresponding to a maximum mass of $m_{max,e} = 0.28 M_{Jup}$.

The influence of the central binary means that circumbinary planets have unique formation pathways, and it is worth noting that the system discussed in this Article is only the second system known to host multiple planets. Hydrodynamical simulations have shown that the preference for observed circumbinary planets to be located close to the stability boundary can be explained by formation at large distances in circumbinary disks, followed by inward migration and stalling at the

<footnote n="10">10 uncertainties are provided as the last two significant digits, within parentheses except where asymmetric. Dates are given in BJD−2,450,000. <sup>a</sup>Data from ref. 17. <sup>b</sup>Data from ref. 18.</footnote>
<figure><img src="image_7.png" /></figure>

**Fig. 4 | Detection limits for additional undetected planetary signals in the BEBOP-1 radial-velocity data.** Hexbin plot denoting the density of ~100,000 posterior samples obtained from four separate kima runs on BEBOP-1 radial-velocity data with $N_p$ fixed to 1. The red line shows the calculated 99% detection limit, with its associated error denoted by the lines thickness. The letters mark in orange circles mark the positions of the detected planet BEBOP-1 c and planet TOI-1338 b, which is not detected in RV data alone. The orange line shows the detection limit on samples with eccentricity $e < 0.1$ (the upper limit on TOI-1338 b's eccentricity¹⁷).

with 99% confidence. This is compatible with the $33 \pm 20 M_⊕$ from

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

**Fig. 5 | Radius versus Mass plot of all transiting circumbinary planets and planets orbiting single stars.** For the single-star planets we restricted the sample to planets with mass and radius percentage errors of less than 20%. For the circumbinary planets we show all transiting planets, even if the mass is only constrained by an upper limit (triangle symbol). Error bars show $1\sigma$. TOI-1338 b is highlighted in red, with one of the lowest densities known. The newly discovered BEBOP-1 c is not plotted as it does not transit. These data are publicly available on the NASA Exoplanet Archive⁹⁶.

edge of the inner cavities formed by the central binaries³⁵⁻³⁷. We utilized a purpose-built simulation code to study the formation of circumbinary planets as a means of understanding plausible formation scenarios for the BEBOP-1 system. Details of this code are available in the Methods. The model includes the *N*-body integrator MERCURY6, adapted to include a central binary system³⁸,³⁹, and it incorporates prescriptions for a viscous circumbinary disk that includes the effects of an eccentric, precessing central cavity and photoevaporative winds⁴⁰, pebble accretion onto planetary seeds⁴¹, planet migration⁴² and gas accretion onto growing planets⁴³,⁴⁴. Our suite of simulations produced numerous systems that were qualitatively similar to BEBOP-1 (see Supplementary Information), with TOI-1338 b and BEBOP-1 c analogues landing on stable orbits at their observed locations.

With no detection of the 95 d planet in our data, the best we could do was to calculate a detection limit for this period region, which can be seen in Fig. 4. We found that TOI-1338 b has a mass of <21.8 $M_⊕$

with 99% confidence. This is compatible with the $33 \pm 20 M_⊕$ from ETVs¹⁷; Fig. 5 shows a radius versus mass plot for transiting circumbinary planets and planets orbiting single stars. Combining this mass with a planetary radius of ~6.9 $R_⊕$ (ref. 17), TOI-1338 b has a maximum mean planetary density of <0.36 g cm⁻³. This allowed us to calculate the transmission spectroscopy metric¹³ for planet b using the upper limit on the mass (<21.8 $M_⊕$); this yielded a minimum value of >39 (Methods). Of the now 15 known circumbinary exoplanets, TOI-1338 b is the only one for which JWST transmission spectroscopy can currently be pursued. If we are to unravel the mysteries of circumbinary Tatooine-like exo-atmospheres, the BEBOP-1 system may provide a new hope.

## Methods

### Collected observations

BEBOP-1 was selected for the BEBOP programme¹⁴ from a large sample of low-mass eclipsing binaries identified by the EBLM project¹⁸ under the name EBLM J0608-59. Those binaries were detected as part of the Wide Angle Search for Planets, as candidate transiting planets that were later shown using radial velocities to be 'false positive' eclipsing binaries. To be part of the EBLM and BEBOP sample, a system needs to be a single-lined eclipsing binary. For identification and selection, the southern sample exclusively used the CORALIE spectrograph (resolution $R \approx 45,000$; mounted on the 1.2 m Euler Swiss telescope at La Silla, Chile). BEBOP-1 was first observed as part of the EBLM project in 2009.

The BEBOP sample represents a subsample of the EBLM sample, for which the binaries were selected to optimize its planet-finding capability. The main discriminator is the obtained radial-velocity precision, which is typically a function of magnitude. There is also a bias towards wider (≥5 d) binaries, as rapid rotation due to tidal locking broadens spectral lines and reduces precision. Fortunately, the Kepler results in fact show a dearth of planets transiting the tightest binaries⁴⁵. We also progressively removed binaries where a long-term radial-velocity trend revealed the presence of a third star, which is believed to have a detrimental impact on planet occurrence⁴⁵. Finally, systems were also removed on the basis of simple activity indicators such as the line bisector²⁸. Targets within the BEBOP sample received longer exposures than the rest of the EBLM sample, to improve precision. BEBOP-1 was first observed as part of the BEBOP project in 2014. Further details of the selection process can be found in ref. 14.

There are 55 radial velocities obtained with CORALIE that were not used in this Article. Their typical uncertainty is 25 m s⁻¹, much higher than the signal we identified. Some of these radial velocities were used to determine the spin-orbit angle of the binary⁴⁶.

**HARPS.** HARPS is a high-resolution, high-precision échelle spectrograph built for the detection of exoplanets⁴⁷. It has a resolution $R \approx 100,000$ and typically achieves a long-term stability under 1 m s⁻¹. It is mounted on the European Southern Observatory (ESO) 3.6 m telescope at La Silla, Chile.

The southern BEBOP sample, including BEBOP-1 (under the name J0608-59), was observed by HARPS under two ESO large observing programmes (programme IDs 1101.C-0721 and 106.212H; PIA H.M.J.T.). Between 2018 April 8 and 2022 April 18 a total of 61 spectra were obtained for BEBOP-1. Typical exposure times of 1,800 s were obtained, with a median radial-velocity precision of 5.73 m s⁻¹.

HARPS data were reduced with the HARPS Data Reduction Software (DRS) version 3.5 (which is hosted at the Observatory of Geneva). Descriptions of how the DRS works can be found in ref. 48 and in ref. 49. Spectra were correlated using a weighted numerical mask that matched the spectral type of the target producing a cross-correlation function. A Gaussian function was fitted to the cross-correlation function to find the mean radial velocity.

DRS version 3.5 is not very different from version 3.0.0 held by the European Southern Observatory (https://www.eso.org/sci/software/pipelines/harps/specifics.html), but it allowed us to recorrelate the
spectra with our own specifications. The standard DRS was built to study single, slowly rotating stars. It assumes two quantities: a common mean radial velocity for all spectra of a given system, and a correlation window of 30 km s⁻¹ on either side of that mean velocity. In the case of binaries, the velocity changes according to the binary phase. In addition, some of the BEBOP sample targets rotate fast enough that a 30 km s⁻¹ window is not sufficient. Using version 3.5 allowed us to centre the correlation window to the observed velocity at each epoch, thus ensuring the same set of absorption lines were used to produce the radial velocities, epoch after epoch. For BEBOP we used the standard 30 km s⁻¹ correlation window and a G2 mask to produce the cross-correlation window. The reduction software provided the radial velocity, its uncertainty (determined from photon noise and the line width), the line width (full-width at half-maximum (FWHM)) and the span of the bisector slope (Bis_Span) and automatically produced corrections to the barycentre of the Solar System.

All HARPS data are available at the ESO public archive (http://archive.eso.org/eso/eso_archive_main.html) by searching for EBLM J0608-59, and can be found in Supplementary Table 2.

**ESPRESSO.** Following the discovery of TOI-1338 b (ref. 17) 20 radial-velocity measurements were obtained in 2019 with the ESPRESSO spectrograph located in Paranal, Chile,²⁰ in an attempt to confirm the detection of the planet and more accurately constrain its mass (programme ID 103.2024; PI A.H.M.J.T.). Exposure times were typically 900 s in duration and yielded a median precision of 2.83 m s⁻¹. These data, combined with previously obtained HARPS spectra, yielded no detection of TOI-1338 b, but hinted at the presence of an outer companion planet with an orbital period of ~200 d. To confirm this candidate planet and to attempt to detect the inner transiting planet, a further 103 ESPRESSO spectra were obtained with exposure times of 900 s, yielding a median precision of 2.63 m s⁻¹ (programme ID 106.216B; PIM.R.S.). Of these 103 measurements, 3 were obtained during the primary eclipse of the binary and 1 during the transit of planet b. These four observations have been discarded from our analysis as they were affected by the Rossiter-McLaughlin effect, which we did not model here⁵⁰,⁵¹.

In our radial-velocity analysis, we fitted for an offset between ESPRESSO data obtained in 2019 and those obtained in 2021, post COVID closures. This 4.6 m s⁻¹ offset, although small, was introduced by calibration lamp changes during the COVID closure of telescopes in 2020⁵². Our 2019 dataset was obtained after the fibre change on the instrument in 2019. Without accounting for this offset, the planetary signal was still detected in the data. The COVID closures of both the La Silla and Paranal observatories delayed the collection of data on this system and consequently the discovery of this planet by a year. Complete phase coverage of BEBOP-1 c's orbit was finally achieved when the final data points were collected in March 2022.

ESPRESSO data were reduced using version 2.3.3 of the ESPRESSO pipeline (publicly available at ESO) with a procedure similar to that used for HARPS, but adapted to ESPRESSO. Each cross-correlation function was obtained with a G2 template spectra, with a step size of 0.5 km s⁻¹, and a 300 km s⁻¹ window centred on the radial velocity of the binary. To summarize, we obtained 123 ESPRESSO observations in total, with 20 observations in ESPRESSO 2019, and 103 in ESPRESSO 2021/2022. The total timespan of the HARPS and ESPRESSO combined is 1,472 d.

All ESPRESSO data are available at the ESO public archive by searching for J0608-59. All radial-velocity data used in our analysis can be found in Supplementary Tables 2–4.

**TESS.** BEBOP-1 was a target for the radial-velocity survey before the launch of TESS. It was then later found to contain a 95 d transiting planet (TOI-1338 b), based on four sectors of TESS data (three in short cadence, 120 s). More data have been taken in the years since. These data are mentioned here for completeness, but we did not include them in the analysis: BEBOP-1 c is an independent radial-velocity discovery.

**Outlier treatment.** Outliers were identified using two methods. The first (method 1) involved a priori finding outliers in the span of the bisector slope (Bis_Span, a measure of line shape²⁸) and FWHM, and excluding these data. The radial-velocity analysis in this Article used this method. The second method used a Student's t distribution to account for outliers as part of the sampling process (for example, ref. 53); we used this to check the results from method 1.

Preliminary outliers where the wrong star was observed or the observations occurred during an eclipse or planetary transit were removed first. The following methods to deal with outliers were then applied to the remaining data.

For a priori removal, each set of points was fitted as a mixture model with inlier and outlier populations, following the method described by Hogg et al.⁵⁴. For both FWHM and Bis_Span, the inlier model was simply a constant with a small scatter allowed in addition to the uncertainties (which were taken as double the radial-velocity uncertainty). The outlier model also fitted a constant, but with a very wide population scatter.

We used PyMC3⁵⁵ to fit a mean and scatter for both the inlier and outlier distributions, as well as a parameter f for the proportion of points that are outliers; the prior on f was a β[1.5, 9] distribution favouring a low proportion of outliers. These parameters were then used to calculate the probability of the ith point being an outlier for each posterior using the following equation:

$$P_{i,out} = \frac{fL_{i,out}}{fL_{i,out} + (1-f)L_{i,in}} \quad (1)$$

where $L_{i,out}$ and $L_{i,in}$ are the likelihoods for the ith point in the respective outlier and inlier distributions.

The power of this method lies in removing the human factor from outlier identification, and not having to resort to sigma-clipping, which does not usually consider measurement uncertainties. A total of 17 outliers were identified this way.

For the Student's t method, we ran the fits on the dataset accounting for outliers using a Student's t distribution. Outliers in terms of radial velocity could then be directly identified, rather than relying on the Bis_Span and FWHM indicators used in method 1. This method identified 6 outliers, of which 3 were shared with the 17 outliers identified in method 1; the 3 points were particularly strong outliers in the Student's t analysis.

Importantly, the resulting detection of BEBOP-1 c and non-detection of TOI-1338 b were not affected by the method of dealing with outliers, and all the parameters from the Student's t fit were consistent with those found removing the outliers a priori.

In summary, three outliers were identified in the HARPS data in terms of the FWHM, one of which was also identified using Bis_Span and is such a strong outlier that it is suspected that the wrong star was observed. One outlier was identified in the ESPRESSO 2019 data using the FWHM, and 13 outliers were identified in the ESPRESSO 2021/2022 data using the FWHM (two of which were also identified in Bis_Span). Data points identified as outliers are flagged in Supplementary Tables 2–4. We denoted outliers, excluded from the analysis, with the following flags: Wrong star W, Bis_Span B, FWHM F, During binary transit R and during planetary transit P.

## Modelling of the data

**Kima.** For the radial-velocity analysis we used the kima package²¹, which models radial-velocity data with a sum of Keplerian functions from $N_p$ orbiting planets, and estimates the posterior distributions for each of the orbital parameters.

To sample the posterior distribution kima uses a diffusive nested sampling approach⁵⁶. This provides kima with estimates on the evidence for each model, allowing for model comparison⁵⁷,⁵⁸. This model comparison can then be used to compare the ratio of evidence
between models with different numbers of Keplerian signals. It should be noted that $N_p$ is a free parameter like any other.

The radial-velocity analysis follows the procedure described in Standing et al.⁴⁹ and Triaud et al.¹² with added post-Keplerian corrections described by Baycroft et al.⁵⁰. These corrections accounted for relativistic and tidal effects on the radial-velocity measurements from the binary orbit that have been derived and tested in refs. 22–25,49. In addition, we directly fitted for apsidal precession of the binary within kima, which is the largest Newtonian effect that we expected. We found a value of $\omega_{\text{bin}} = 66.0^{+65.0}_{-54.7}$ arcsec yr⁻¹, which distinct but scarcely different from zero. Details of this implementation are provided in Baycroft et al.⁵⁰. The prior distributions used in our analysis can be found in Supplementary Table 1.

**Resulting fit.** The kima combined fit of our HARPS and ESPRESSO data favours a single-planet model in addition to the binary orbit with a Bayes factor of >29,000. The Bayes factor is the ratio of the Bayesian evidence between two competing models. In our case the Bayes factor is the number of posterior samples obtained with $N_p$ planetary signals over those with $N_p - 1$. The Bayes factor indicates the measure of support in favour of one model over another⁴⁹,⁶⁰,⁶¹.

A Bayes factor of >150 indicates very strong evidence⁴⁹,⁶² in favour of a single planet in addition to the binary orbit.

When searching for the inner transiting planet in the data we obtained a Bayes factor of 1.4 in favour of two planetary signals. This is categorized as inconclusive evidence of any further signals in the data.

Forcing a fit on a 95 d Keplerian signal yielded no clear signal corresponding to TOI-1338 b with posterior samples having a semi-amplitude consistent with 0 m s⁻¹. This fit still produced a clear detection of BEBOP-1 c.

To obtain orbital parameters for the system, we followed the same procedure detailed in Standing et al.⁴⁹. Proposed posterior samples with orbits that cross one another, or enter the instability region of the binary, were removed. The remaining samples were clustered with the HBDSCAN clustering algorithm⁶³. Clusters corresponding to the binary and BEBOP-1 c were then plotted using the Corner package version 2.0.1 (ref. 64). Orbital parameters were then determined as the 50th percentile, with 1σ uncertainties estimated from the 14th and 84th percentiles of the cluster. Corner plots for the binary and BEBOP-1 c can be seen in Supplementary Figs. 1 and 2 respectively.

**Constraints on TOI-1338 b and further planetary companions.** Using kima in the same way as Standing et al.⁴⁹ and Triaud et al.¹² we calculated a detection limit on additional undetected planetary signals in the radial-velocity data. First, we removed the highest-likelihood Keplerian model corresponding to BEBOP-1 c from the data. Following this, to calculate the detection limit, we fixed $N_p = 1$ in our fit, included the binary with tight priors on its parameters and obtained posterior samples of all remaining signals that were compatible with the data. This method is conceptually similar to that described in Tuomi et al.⁶⁵, where the sampler was forced to fit an additional signal to the data. Any signals found are compatible with the data, but yet undiscovered, as only one planetary signal was present in our original fit.

The resulting posterior samples can be seen as a greyscale density plot in Fig. 4, with a blue 99% contour corresponding to our detection limit. The blue limit demonstrates that our analysis was sensitive to additional sub-Saturn mass planets for periods out to 2,000 d, while sensitive to Neptune-mass planets near the instability limit.

We note a density of posterior samples at periods around 100 d in agreement with the presence of TOI-1338 b seen in TESS photometry¹⁷. Following the findings described in Standing et al.⁴⁹, we calculated an additional detection limit to place an upper limit on the mass of TOI-1338 b. This limit was calculated in the same manner as the blue line, although only using posterior samples with eccentricities <0.1 (the approximate upper limit on the eccentricity of TOI-1338 b¹⁷, and

a value allowed by our orbital stability analysis). A running mean was applied to the red eccentricity cut line to ensure that >1,000 samples were in each bin. The resulting red detection limit demonstrated the sensitivity of our analysis to circular planets at the orbital period of TOI-1338 b was $-2.4 \pm 0.1 \text{ m s}^{-1}$. This allowed us to place an upper limit on the mass of TOI-1338 b of $21.8 \pm 0.9 M_{\oplus}$. With a radius of $-6.9 R_{\oplus}$, we calculated the mean density of TOI-1338 b to be <0.36 g cm⁻³. The mass of planet b from Kostov et al.¹⁷ was primarily determined using the apsidal precession imprinted on the ETVs under the assumption of a single planet. Our study exclusively used radial velocities; as they are not sensitive to the planet's mass, we only placed an upper limit.

The uncertainty of the detection limit was determined by performing numerical experiments with kima. We first generated a radial-velocity timeseries, following a Gaussian distribution, with dates following a log-uniform distribution. Then we ran kima just like any other system and produced a large number of posteriors (in our case 685,000). Following that, we computed a detection limit as performed in Standing et al.⁴⁹ and called this our fiducial case. The next step was to calculate detection limits for many subsamples of the posterior and measure their fractional distanced to the fiducial detection limit. We found that the fractional error in the position of subsamples' detection limits followed a square root law. Consequently, we fitted these data with a function $y = ax^{-1/2}$, where y is the fractional uncertainty of the detection limit at a given orbital period and x is the number of posteriors within a given posterior subsamples. We found $a = 15.1 \pm 0.3$. We used this relation to represent an uncertainty in Fig. 4.

**False and true inclusion probabilities.** We computed the true inclusion probability (TIP) and false inclusion probability (FIP) for a signal being present in the data over various frequency bins²⁶. Plotting this in a periodogram gives similar information to the detection limit described above. However, where the detection limits give us information about an upper limit on the mass of a potential planetary signal, the TIP periodogram effectively gives us the probability of there being a planetary signal at that given period. It has been shown that using this metric as a detection criterion is optimal in the context of exoplanet detection²⁷. We found one peak in the TIP periodogram beyond the detected planet c at ~700 d with a TIP of 0.94. This is shown in Supplementary Fig. 3, which presents the FIP and TIP periodograms. We attributed this signal to the window function as seen in Supplementary Fig. 4, although further observations will provide greater insight into the source of the signal.

## Stellar activity

Chromospheric emission lines observed in the stellar spectrum often trace stellar magnetic activity. Here we describe the activity indices measured and used in our analysis. We used the open-source package ACTIN⁶⁶ to measure the Ca II H&K, He I, Na D doublet (a and b), Ca I and Hα activity indices.

First, we computed a generalized Lomb-Scargle periodogram⁶⁷ of each activity indicator, from individual instruments, and used the combined dataset. In Supplementary Fig. 4a,b we show the window function of each individual and combined dataset and the periodogram of the radial velocity. We computed false alarm probability levels of 10%, 1% and 0.1% using a bootstrap randomization of the dataset. We considered a signal to be significant if the false alarm probability level was <0.1%. We searched for periodic signals of chromospheric activity indicators to investigate the possibility of BEBOP-1 c's period $P_c$ being produced by magnetic activity of BEBOP-1A (the primary). Supplementary Fig. 4c-g depicts the periodogram of the Ca II, Hα, Na D, He I and Ca I indices. None of the activity indicators showed any significant signal at or near $P_c$ (-215 d). We saw a significant period at ~7.5 d in Na D (false alarm probability ≈ 0.1%) and He I (false alarm probability ≈ 1%); however, the nature of these signals was not obvious. We saw weak signals in Hα at 270 d caused by a combination of window function
and a long-period term, but the signal is not significant. In addition, the activity time series shown in Supplementary Fig. 5 does not show any obvious modulation.

We searched All-Sky Automated Survey (ASAS)⁶⁸ archival data to constrain the rotation period and to investigate the nature of the BEBOP-1 c signal. The ASAS observations span more than 10 years. The ASAS data shows a strong signal at 61 d (10% false alarm probability), which could be related to stellar rotation (-3P<sub>rot</sub>).

In addition, analysis of the radial-velocity bisector shows no correlation with the -215 d signal, yielding a Pearson correlation coefficient of <0.052.

The ASAS data and activity indicators do not show any significant signal at -215 d, therefore this signal is likely to be planetary in nature.

## System dynamics

In this section we deal with various aspects of the orbital dynamics of BEBOP-1, first looking at the stability of the system, and then estimating whether the outer planet, BEBOP-1 c, will show transits at some point into the future.

**Orbital stability.** The orbital solution given in Table 1 shows a compact two-planet circumbinary system ($a_{\text{bin}} = 0.13$ au, $a_b = 0.46$ au and $a_c = 0.79$ au). The inner planet has a maximum mass similar to Neptune ($m_b = 21.8 M_{\oplus}$) and is close to the circumbinary stability boundary¹⁷,⁶⁹, whereas the minimum mass of the outer planet is similar to Saturn ($m_c = 0.21 M_{\text{Jup}}$). As a consequence, we expected strong mutual gravitational interactions between the stars and the planets. To study the stability of the system, we performed a global frequency analysis⁷⁰,⁷¹ in the same way as conducted for other circumbinary planetary systems⁷². In our analysis, we always considered the maximum mass of the inner planet, which can be seen as a superior limit for stability. Moreover, the mass of this planet is relatively small and does not greatly impact the stability, and so smaller masses do not affect the global picture described here.

The system was integrated on a regular 2D mesh of initial conditions in the vicinity of the best fit (Table 1). We used the symplectic integrator SABAC⁴,⁷³, with a step size of 0.001 yr and general relativity corrections. Each initial condition was integrated for 5,000 yr, and a stability indicator, $\Delta = |1 - n_b/n_{\text{bin}}|$, was computed. Here, $n_b$ and $n_{\text{bin}}$ are the main frequency of the mean longitude of the inner planet over 2,500 yr and 5,000 yr, respectively, calculated via the frequency analysis⁷¹. Note that while the osculating mean motion varies over a single planetary orbit as a result of energy exchange with the binary and planet c, $n_b$ calculated over time intervals much longer than any resonant variations is constant for stable systems, while it drifts for unstable systems until one of the planets escapes. The results are reported in colour in Supplementary Figs. 6–8, where yellow represents strongly chaotic trajectories with $\Delta > 10^{-2}$, while extremely stable systems with $\Delta < 10^{-8}$ are shown in purple/black. Orange indicates the transition between the two for $\Delta \approx 10^{-4}$.

Observationally, only the inner planet's eccentricity is well constrained: $e_b = 0.088 \pm 0.004$ (ref. 17). In contrast, only an upper limit can be placed for the outer, non-transiting planet: $e_c < 0.1$ (Table 1).

Therefore, in a first experiment, we explored the stability of the system in the plane ($e_b, e_c$), assuming coplanar orbits. The results are shown in Supplementary Fig. 6. We observed that the eccentricities of both planets must be smaller than 0.1 to ensure a stable system. Moreover, as we already have $e_b \approx 0.1$ (ref. 17) for the inner planet, the only remaining possibility for the outer planet is $e_c \approx 0$. We thus adopted $e_c = 0$ in the following stability analyses.

The inner planet was not detected in the radial-velocity data, its presence can only be inferred from photometric measurements¹⁷. Therefore, in a second experiment, we explored the stability of the orbit of the inner planet by varying the orbital period and the eccentricity of this planet, assuming coplanar orbits (Supplementary Fig. 7). This

allowed us to test the compatibility of the two independent observational datasets. We observed that the inner planet lies in a small stability region that is shaped by the presence of the binary system on the left-hand side and by the presence of the massive outer planet on the right-hand side. We thus concluded that the two-planet circumbinary solution is reliable. However, we assumed $e_c = 0$ to draw the stability map obtained in Supplementary Fig. 7. As we increase the eccentricity of the outer planet, the small stable region was quickly degraded and completely disappeared for $e_c > 0.1$, in conformity with the results shown in Supplementary Fig. 6. Supplementary Fig. 7 shows that the inner planet, TOI-1338 b, is surrounded by unstable regions. If the planet migrated to its current location, it would have to pass through two resonances with the outer planet. This could indicate consequences for formation migration in the system.

The radial-velocity technique alone is unable to constrain the inclination, $I_c$, and the longitude of the node, $\Omega_c$, of the outer planet. As a result, we could only determine the minimum value of its mass, which corresponds to the coplanar system ($I_c = 90^\circ$) that we assumed in previous analyses. However, in contrast to the inner planet, the outer planet does not transit, and thus cannot lie exactly in the same orbital plane as the remaining bodies in the system. In a final experiment, we explored the stability of the system by varying the inclination and longitude of the node of the outer planet (Supplementary Fig. 8). We observed that the system can be stable within a circle centred at the coplanar solution, which corresponds to mutual inclinations smaller than 40°. A larger stability region also exists centred at ($\Omega_c = 180^\circ, I_c = 90^\circ$), but it corresponds to mutual inclinations higher than 120° (that is, to retrograde orbits, which are more unlikely from a formation point of view). As we changed $I_c$, the mass of the outer planet increased. At the boundary of stability, $I_c = 90 \pm 40^\circ$, we obtained a maximum mass of the outer planet of $m_{\text{max},c} = 0.28 M_{\text{Jup}}$. Note, however, that we could also get a mutual inclination of 40° with $I_c = 90^\circ$ and $\Omega_c = \pm 40^\circ$, for which the mass of the planet was on the low side, $m_{\text{min},c} = 0.21 M_{\text{Jup}}$. Thus, the source of instability must be mutual inclinations above the 40° threshold, and not masses higher than $0.28 M_{\text{Jup}}$. Again, we assumed $e_c = 0$ to draw the stability map obtained in Supplementary Fig. 8. As we increased the eccentricity of the outer planet, the radius of the stable circle quickly shrank and completely disappeared for $e_c > 0.1$, in conformity with the results shown in Supplementary Fig. 6.

**Possible transits of BEBOP-1 c.** Orbits of circumbinary planets exhibit nodal precession. This changes the orientation of the planet's orbit with respect to both the binary and the observer. This means that a planet changes from a transiting to non-transiting configuration¹⁵,³⁰. The precession timescale is typically on the order of decades, and hence can be observable with long baselines such as those possible with Kepler and TESS. Kepler-413 exhibited this on-off transit sequence³¹. Kepler-453 did not start transiting until half way through the Kepler mission³².

Refs. ³³,³⁴ derived an analytic criterion to determine whether a planet will ever enter transitability on the primary (A) or secondary (B) star during its precession cycle:

$$ \Delta I > \left| \frac{\pi}{2} - I_{\text{bin}} \right| - \sin^{-1} \left( \frac{a_{\text{A,B}}}{a_{\text{p}}} \sin \left| \frac{\pi}{2} - I_{\text{bin}} \right| + \frac{R_{\text{A,B}}}{a_{\text{p}}} \right), \quad (2) $$

where $\Delta I$ is the mutual inclination between the planet and binary orbit. It may seem counter-intuitive at first that a planet–binary misalignment makes transitability more likely. The reason is that the planet's sky inclination oscillates around the binary's sky inclination, with an amplitude equal to $\Delta I$. A higher $\Delta I$ therefore makes it more likely that the planet's inclination will become close to 90°, and hence show transitability. For eclipsing binaries $I_{\text{bin}}$ is near 90°, and hence the vast majority of circumbinary planets orbiting eclipsing binaries will eventually transit.

Placing the values for BEBOP-1 c into equation (2) yields $\Delta I_{\text{min}} = -0.11^\circ$ for the primary star. The negative inclination criterion
means that BEBOP-1 c is guaranteed to eventually transit regardless of its mutual inclination. This is mainly because the binary is so well aligned with our line of sight ($I_{\text{bin}} = 89.658^{\circ}$) and the primary star is relatively large ($R_{\Lambda} = 1.299 R_{\odot}$). Although we think that BEBOP-1 c is guaranteed to eventually transit, we are unable to predict when and how frequently. The precession period is:

$$P_{\text{prec}} = P_{\text{p}} \frac{16}{3} \left( \frac{a_{\text{p}}}{a_{\text{bin}}} \right)^2 \frac{1}{\cos \Delta l}, \quad (3)$$

where we assumed circular orbits$^{30,74}$. For BEBOP-1 c $P_{\text{prec}} = 119$ yr, assuming that $\Delta l$ is close to zero, as is the case for known circumbinary planets. With modern advances in medical science, there is a chance that the authors will live to see BEBOP-1 c transit.

## System formation

The coplanarity between the planet and binary orbit planes for the circumbinary planets that have been discovered so far strongly suggests that these planets were formed in circumbinary protoplanetary disks that were themselves closely aligned with the binary orbit plane. We present here the results from a suite of simulations of circumbinary planet formation that were performed to examine plausible scenarios for the origin of the TOI-1338 b and BEBOP-1 c system.

### Circumbinary disk and N-body model.
Our simulations were performed using a newly developed code designed specifically to examine the formation of circumbinary planet systems. The code employs the N-body symplectic integrator MERCURY6, adapted to include a central binary system$^{38,39}$. The code utilizes the 'close-binary' algorithm$^{39}$ that calculates the temporal evolution of the positions and velocities of each body in the simulations with respect to the centre of mass of the binary stars, subject to gravitational perturbations from both stars and other large bodies. The evolution of the circumbinary disk is calculated using a 1D viscous $\alpha$-disk model$^{75}$, solving the standard diffusion equation with additional terms that account for the possible presence of a gap-forming planet and a photoevaporative wind:

$$\frac{d\Sigma}{dt} = \frac{1}{r} \frac{d}{dr} \left[ 3r^{1/2} \frac{d}{dr} (v\Sigma r^{1/2}) - \frac{2\Lambda\Sigma r^{3/2}}{GM_*} \right] - \frac{d\Sigma_{\text{pe}}}{dt} \quad (4)$$

Here, $\frac{d\Sigma_{\text{pe}}}{dt}$ is the rate of change in surface density due the photoevaporative wind, $\Lambda$ is the torque per unit mass that operates when a planet forms and becomes massive enough to open a gap in the disk and $v$ is the disk viscosity$^{75}$, given by:

$$v = \alpha c_s^2 / \Omega, \quad (5)$$

where $c_s$ is the local isothermal sound speed, $\Omega = \sqrt{GM_* / r^3}$ is the Keplerian frequency and $\alpha$ is the viscosity parameter. The planet torque per unit mass is given by:

$$\Lambda = \text{sign}(r - r_p) q^2 \frac{GM_*}{2r} \left( \frac{r}{|\Delta_p|} \right)^4, \quad (6)$$

where $q$ is the planet/star mass ratio, $r_p$ is the planet orbital radius and $|\Delta_p| = \max(H, |r - r_p|)$; $H$ is the local disk scale height. We assumed that the disk is in thermal equilibrium, and as such we used an iterative method to solve the following equation$^{76}$:

$$Q_{\text{irr,A}} + Q_{\text{irr,B}} + Q_v + Q_{\text{cloud}} - Q_{\text{cool}} = 0, \quad (7)$$

where we balanced irradiation heating from both central stars ($Q_{\text{irr,A+B}}$), background heating from the residual molecular cloud ($Q_{\text{cloud}}$) and viscous heating ($Q_v$) with blackbody cooling ($Q_{\text{cool}}$).

Mass loss due to photoevaporative winds, resulting from high-energy photons emitted by both the central stars and nearby externally located stars, occurs on a timescale determined by the adopted flux of high-energy photons$^{77,78}$.

Hydrodynamical simulations have shown that circumbinary disks develop precessing eccentric inner cavities, with their sizes and eccentricities depending on the properties of the binary and disk$^{79,80}$. This cavity plays an important role in the evolution of planets that migrate into its vicinity. To mimic the presence of this cavity in the 1D viscous disk model, we used a variable $\alpha$ model fitted to results from 2D hydrodynamic simulations using FARGO3D$^{81}$, which showed that the mass flux through the cavity remained roughly constant within the region around the binary. For these FARGO3D simulations we adopted a computational domain that spanned 0.13–20 au, and for the numerical grid we adopted a resolution of $N_x \times N_\phi = 768 \times 512$ (radial $\times$ azimuth), with logarithmic radial spacing. We used the stellar parameters outlined in this work, and the disk mass lying within 40 au was set to equal 6% of the combined binary mass. The viscosity parameter $\alpha = 10^{-3}$ and the constant disk aspect ratio was set to 0.05. To ensure that the simulations reached equilibrium, we ran them for 30,000 binary orbits.

### Planet migration.
Planets with masses that substantially exceed a lunar mass undergo migration because of gravitational interactions with the surrounding disk. We followed ref. 82 and included type-I migration in the model via the culmination of the torque formulae that account for Lindblad and corotation torques$^{42}$, as well as those of eccentricity and inclination damping$^{83,84}$:

$$\Gamma_{\text{I,tot}} = F_{\text{L}}\Gamma_{\text{LR}} + \left\{ \Gamma_{\text{VHS}}F_{p_v}G_{p_v} + \Gamma_{\text{EHS}}F_{p_v}F_{p_v} \sqrt{G_{p_v}G_{p_v}} + \Gamma_{\text{LVCT}}(1 - K_{p_v}) \right. \\ \left. + \Gamma_{\text{LECT}}\sqrt{(1 - K_{p_v})(1 - K_{p_v})} \right\} F_e F_i \quad (8)$$

where $\Gamma_{\text{LR}}$, $\Gamma_{\text{VHS}}$, $\Gamma_{\text{EHS}}$, $\Gamma_{\text{LVCT}}$ and $\Gamma_{\text{LECT}}$ are the Lindblad torque, vorticity- and entropy-related horseshoe drag torques, and linear vorticity- and entropy-related corotation torques, respectively, as given by equations (3)–(7) in ref. 85. The functions $F_{p_v}$, $F_{p_v}$, $G_{p_v}$, $G_{p_v}$, $K_{p_v}$ and $K_{p_v}$ are related to the ratios between viscous/thermal diffusion timescales and horseshoe libration/horseshoe U-turn timescales, as given by equations (23), (30) and (31) in ref. 85. If a planet becomes massive enough to form a gap in the disk$^{86}$, its migration changes from type I to type II. When transitioning from type-I to type-II migration the model accounts for gap formation self-consistently by calculating the torque acting on the disk due to the planet, with the back-reaction driving migration$^{87}$. Planets in the vicinity of the inner cavity of the circumbinary disk will experience time-varying forces arising from its eccentricity and precession. These are not naturally captured in a 1D disk model, so we have added terms to the evolution equations to incorporate these effects. Using the results from the FARGO3D simulations described above, we calculated an azimuthally averaged eccentricity profile that, when used in conjunction with a uniformly precessing disk, can give the velocity of the gas in an eccentric disk at any azimuthal location. Following the calculated precession rates from previous works$^{80}$, we assumed that the disk precession rate was equal to 3,000 binary orbits. Using the steadily precessing 2D eccentric structure of the disk obtained from the hydrodynamic simulations, we created a 2D map of the radial and azimuthal forces that an embedded particle at any position in the disk would be subject to. This map was interpolated during the N-body simulations to find the force acting on any planet that forms.

### Planetary growth.
The growth of planetary seeds in the model occurs via pebble accretion, and once a planet exceeds an Earth mass it starts to accrete gas from the disk. Pebble accretion models have shown that as protoplanetary disks evolve, dust coagulates into pebbles and settles in the midplane, where the pebbles drift inwards$^{41}$. This creates a pebble production front that expands outwards over time. As the
pebbles drift through the disk, they encounter planetary embryos that can accrete the pebbles with an efficiency that depends on the properties of the disk and planet⁸⁸,⁸⁹. We adopted the models of refs. 41,89 to account for the production of pebbles in the disk, as well as accretion onto planetary embryos. While such models were derived for single stars, with no perturbing forces from a central binary, they should be sufficiently accurate for the circumbinary disk models we present here. This is because the system is a close binary, and as such once material is only around a few astronomical units away from the central stars the gravitational force it feels is very similar to that of a single star. Given that pebble production primarily occurs far from the binary, in the outer disk regions, we assumed that the formulation for single stars is accurate. In terms of planets accreting pebbles, the same applies with most planets accreting pebbles at distances greater than 2 au. However, some planets are able to migrate and accrete pebbles near the cavity, and here it would not be appropriate to assume that the pebbles follow circular Keplerian orbits. Results from previous work have shown that dust particles with small Stokes numbers closely follow the gas⁹⁰, and as such we assumed that the pebbles follow the same eccentric profiles used to calculate the torques acting on the planet in the vicinity of the central cavity. The relative velocities between pebbles drifting past the planet were then calculated and used in the equations for pebble accretion⁸⁹.

Gas accretion onto the planets was modelled using a recently developed empirical fit⁴⁴ to the results of detailed 1D envelope structure models⁴³ that account for local disk conditions:

$$ \left(\frac{dM_{\text{ge}}}{dt}\right)_{\text{local}} = 10^{-10.199} \left(\frac{M_{\text{\textit{⊕}}}}{\text{yr}}\right)^{f_{\text{opa}}} \left(\frac{T_{\text{local}}}{1\text{K}}\right)^{-0.963} \left(\frac{M_{\text{core}}}{M_{\text{⊕}}}\right)^{5.6549} \left(\frac{M_{\text{ge}}}{M_{\text{⊕}}}\right)^{-1.159} \\ \times \left[\exp\left(\frac{M_{\text{ge}}}{M_{\text{core}}}\right)\right]^{3.6334} \quad (9) $$

where $T_{\text{local}}$ is the local disk temperature, $f_{\text{opa}}$ is an opacity reduction factor, equal in this work to 0.01, and $M_{\text{core}}$ and $M_{\text{ge}}$ are the planet's core and gas envelope masses, respectively. Given that these models calculate orbit averaged accretion rates, we assumed that this is also valid for planets in eccentric circumbinary disks; however, further work examining the validity of this assumption should be undertaken. We note these models predict that gas accretion starts off very slowly and speeds up as the core and envelope grow in mass. The onset of gas accretion therefore does not generally result in the formation of a gas giant planet because the envelope contracts on its Kelvin-Helmholtz timescale.

**The resulting circumbinary planet systems.** The simulations were initiated by placing 42 planetary seeds throughout the disk between 2 and 20 au, with eccentricities and inclinations distributed uniformly up to $e = 0.02$ and $i = 0.5°$, respectively. The initial masses of the seeds were set to equal one-tenth of the transition mass, which is where pebble accretion becomes efficient as the planetary core begins to accrete from the entirety of its Hill sphere, as opposed to its Bondi sphere. These initial masses were between $10^{-4}$ and $10^{-3} M_{\text{⊕}}$ and are consistent with recent simulations of planetary embryo formation within protoplanetary disks based on gravitational collapse⁹¹. We initialized the gas surface density according to $\Sigma = \Sigma_0(\frac{r}{1au})^{-1}$, where $\Sigma_0$ depends on the initial disk mass. We explored disk masses between 10 and 15% of the combined stellar mass in disks that extend to 100 au, and we considered metallicities of 0.5, 1 or 2× the Solar value. We assumed a viscous $\alpha$ of $2 \times 10^{-3}$, an $f_{\text{41}}$ parameter for extreme ultraviolet photoevaporation of 100 and an external photoevaporative mass loss rate for a disk of size 100 au of $3 \times 10^{-7} M_{\text{⊕}}$ yr$^{-1}$.

Supplementary Fig. 9 shows the mass and semimajor axis evolution from an example simulation that formed a planetary system similar to BEBOP-1. Black data points show the final masses and semimajor axes of the surviving planets and the vertical black dashed line denotes

the stability boundary⁹². The red triangle and plus sign indicate the locations of TOI-1338 b and BEBOP-1c.

The formation of this system occurred as follows. As the pebble front moved outwards, the planetary embryos on the most circular orbits were able to accrete pebbles efficiently, allowing a number of them to grow to masses greater than an Earth mass. These planets began to migrate, generally inwards towards the central cavity, and continued to accrete drifting pebbles, as well as gas. The central cavity in the disk provided a positive surface density gradient in the gas, and this allowed corotation torques to balance Lindblad torques for planets arriving there, resulting in migration stalling. Hence, inward migrating planets started to congregate near the outer edge of the cavity, with some resonant chains forming. As the planets continued to grow in mass, two planets with masses of ~10 $M_{\text{⊕}}$ collided, forming a more massive core. This planet was then able to accrete gas more efficiently and entered a brief period of runaway gas accretion to become a gas giant. As the planet's mass increased, it migrated closer to the central stars, pushing planets interior to it into the cavity, and closer to the stability boundary. This resulted in the innermost planets being ejected from the system as they crossed the stability boundary and were gravitationally scattered by the central stars, although the planet destined to become the TOI-1338 b analogue remained outside of the boundary. The gas giant planet opened a gap in the disk and began to undergo slow, inwards, type-II migration while continuing to accrete gas through the gap at the viscous supply rate. After 1.8 Myr the disk reached the end of its life and was fully dispersed due to photoevaporative winds, and the giant planet stopped accreting and migrating, leaving both planets with masses and semimajor axes similar to those inferred from the observations. The system was allowed to evolve for a further 8 Myr, during which time it remained stable, yielding the final system shown in Supplementary Fig. 9.

This evolution scenario was a common outcome for a number of simulated systems, with planets growing through pebble accretion at large radii in the disk before migrating inwards towards the central cavity. Planets congregated there, leading to collisions and some ejections. Supplementary Fig. 10 shows the mass versus semimajor axis for all planets formed in the simulations, with the colour coding showing the metallicity of the system in comparison with the solar value (this determines the mass in pebbles assumed in the models). Planets shown as grey data points have been lost from the systems through collisions with more massive planets or by ejection. The black dashed line shows the stability boundary, while the black triangle and plus sign show the masses and locations of TOI-1338 b and BEBOP-1c.

Supplementary Fig. 10 shows that formation of TOI-1338 b analogues is a relatively common outcome of the simulations, with numerous super-Earth and Neptune-mass planets being formed with similar orbital periods. This is unsurprising given that these planets become trapped relatively easily around the cavity, close to the observed location of TOI-1338 b. For BEBOP-1c, this planet has fewer simulated planets in its proximity, a result of its observed mass. At this mass planets are in a runaway gas accretion regime, where they quickly grow to masses greater than that of Saturn. Indeed, in our best-fit simulation, the gas giant planet had a mass roughly twice that inferred for BEBOP-1c. In our best-fitting simulations, runaway gas accretion was initiated close to the end of the disk lifetime, which limited the amount of gas that could be accreted by the planets. This necessity for appropriate timing could indicate that BEBOP-1c formed very late in the disk lifetime.

Looking at the metallicities required to form the planets, no systems with 0.5× solar metallicity were able to form systems similar to BEBOP-1. This was because cores with sufficient mass to undergo efficient gas accretion were unable to form. On the other hand, systems with 2× solar metallicity were too effective at forming massive cores, resulting in too many gas giant planets forming and few planets with masses comparable to TOI-1338 b remaining in the system, making them poor fits for comparison. These results show that the metallicity of the system
and how efficiently dust is converted into planets are important for determining the architectures of circumbinary planet systems.

Finally, we note that the orbital configuration of BEBOP-1, with planet periods of 95 and 215 d and the more massive planet on the out-side, is consistent with the prediction of Fitzmaurice et al.⁹ in their study of migrating multiplanet circumbinary systems. They predicted that if there is a planet near the stability limit, as is the case for BEBOP-1, then any equal or more massive outer planets would be located on a period more than double that of the inner planet. Their simulations show that if an exterior planet carrying greater angular momentum reaches the 2:1 resonance with the inner planet, the two lock into resonance and the inner planet is pushed into the inner cavity, resulting in ejection. The periods in BEBOP-1 at 95 and 215 d are explained by the disk dissipat-ing (and hence stopping migration) before the outer planet reached the 2:1 resonance.

## Prospects for atmospheric follow-up with JWST

To assess the suitability of TOI-1338 b for atmospheric characteriza-tion using transmission spectroscopy, we calculated its transmis-sion spectroscopy metric (TSM) as established by Kempton et. al.¹³. A planet's TSM is proportional to the signal-to-noise ratio that could be achieved by observations using JWST, and it depends on the planet's mass, radius and equilibrium temperature and the host star's radius and J-band magnitude.

Calculation of TOI-1338 b's equilibrium temperature is non-trivial; its non-zero eccentricity and changing irradiation levels caused by the binary mean that the temperature oscillates between minimum and maximum values. Full modelling of this, as was done in Kane and Hinkel⁹³, is beyond the scope of this Article, so we estimated a mean equilibrium temperature on the basis of two extremes: when the sec-ondary star is closest to and furthest from the planet. The equilibrium temperature calculated from the irradiation of just the primary star is 501 K, and the motion of the secondary star will cause the true equilib-rium temperature to vary between 603 and 659 K.

We calculated the TSM following Kempton et al.¹³ using a mean $T_{eq}$ = 630 K and the mass upper limit of $M_b$ < 21.8 $M_⊕$ in the first instance, yielding a TSM of 39. This value is below the suggested JWST cutoff of 96 for planets in the sub-Jovian regime. However, the mass used for this calculation was the upper limit; a lower mass would imply a larger scale height and therefore a substantially larger TSM.

In Supplementary Fig. 11 we present the TSM of TOI-1338 b com-pared with other known exoplanets. Its position shows that even the lower limit of its TSM range compares favourably with other sub-Jovian sized planets in its temperature regime. We have indicated with an arrow the range of TSMs possible for a mass down to $10 M_⊕$, which would imply a mean density of 0.17 g cm⁻³, comparable to Kepler-47 c, for example⁴.

We also note that, crucially, of the now 15 known circumbinary exoplanets, TOI-1338 b is the only one for which observations of this kind can currently be pursued. Most Kepler systems are too faint (for example Kepler-47 with a $J_{mag}$ = 13.970 (ref. 94) or Kepler-34 with a $J_{mag}$ = 13.605 (ref. 94)) and others (such as Kepler-16) no longer tran-sit¹. Therefore, despite the challenges it may present, TOI-1338 b may be the only possibility to shed light on the atmospheric make-up of circumbinary planets.

## Data availability

The data that support the findings of this study are available within the Article and the Supplementary Information.

## Code availability

This work made use of the kima²¹, corner⁶⁴, ACTIN⁶⁶ Astropy⁹⁷, numpy⁹⁸, pandas⁹⁹, scipy¹⁰⁰, and matplotlib¹⁰¹. The custom formation simulation code is available upon reasonable request from co-author Gavin A.L. Coleman, email:gavin.coleman@qmul.ac.uk.

## References

1. Doyle, L. R. et al. Kepler-16: a transiting circumbinary planet. *Science* **333**, 1602–1606 (2011).

2. Borucki, W. J. et al. Kepler planet-detection mission: introduction and first results. *Science* **327**, 977–980 (2010).

3. Ricker, G. R. et al. Transiting Exoplanet Survey Satellite (TESS). *J. Astron. Telec.* **1**, 14003-1-14003-10 (2015).

4. Orosz, J. A. et al. Discovery of a third transiting planet in the Kepler-47 circumbinary system. *Astron. J.* **157**, 174 (2019).

5. Meschiari, S. Circumbinary planet formation in the Kepler-16 system. I. N-body simulations. *Astrophys. J.* **752**, 71 (2012).

6. Lines, S., Leinhardt, Z. M., Paardekooper, S., Baruteau, C. & Thebault, P. Forming circumbinary planets: N-body simulations of Kepler-34. *Astrophys. J. Lett.* **782**, L11 (2014).

7. Pierens, A., McNally, C. P. & Nelson, R. P. Hydrodynamical turbulence in eccentric circumbinary discs and its impact on the in situ formation of circumbinary planets. *Mon. Not. R. Astron. Soc.* **496**, 2849–2867 (2020).

8. Martin, D. V. & Fitzmaurice, E. Running the gauntlet - survival of small circumbinary planets migrating through destabilising resonances. *Mon. Not. R. Astron. Soc.* **512**, 602–616 (2022).

9. Fitzmaurice, E., Martin, D. V. & Fabrycky, D. C. Sculpting the circumbinary planet size distribution through resonant interactions with companion planets. *Mon. Not. R. Astron. Soc.* **512**, 5023–5036 (2022).

10. Welsh, W. F. et al. Transiting circumbinary planets Kepler-34 b and Kepler-35 b. *Nature* **481**, 475–479 (2012).

11. Kostov, V. B. et al. TIC 172900988: a transiting circumbinary planet detected in one sector of TESS data. *Astron. J.* **162**, 234 (2021).

12. Triaud, A. H. M. J. et al. BEBOP III. Observations and an independent mass measurement of Kepler-16 (AB) b – the first circumbinary planet detected with radial velocities. *Mon. Not. R. Astron. Soc.* **511**, 3561–3570 (2022).

13. Kempton, E. M. R. et al. A framework for prioritizing the TESS planetary candidates most amenable to atmospheric characterization. *Publ. Astron. Soc. Pac.* **130**, 114401 (2018).

14. Martin, D. V. et al. The BEBOP radial-velocity survey for circumbinary planets. I. Eight years of CORALIE observations of 47 single-line eclipsing binaries and abundance constraints on the masses of circumbinary planets. *Astron. Astrophys.* **624**, A68 (2019).

15. Martin, D. V. & Triaud, A. H. M. J. Planets transiting non-eclipsing binaries. *Astron. Astrophys.* **570**, A91 (2014).

16. Konacki, M., Muterspaugh, M. W., Kulkarni, S. R. & Hetminiak, K. G. The radial velocity Tatooine search for circumbinary planets: planet detection limits for a sample of double-lined binary stars—initial results from Keck I/Hires, Shane/CAT/Hamspec, and TNG/Sarg Observations. *Astrophys. J.* **704**, 513–521 (2009).

17. Kostov, V. B. et al. TOI-1338: TESS' first transiting circumbinary planet. *Astron. J.* **159**, 253 (2020).

18. Triaud, A. H. et al. The EBLM Project. IV. Spectroscopic orbits of over 100 eclipsing M dwarfs masquerading as transiting hot Jupiters. *Astron. Astrophys.* **608**, A129 (2017).

19. Pepe, F. et al. HARPS: ESO's coming planet searcher. Chasing exoplanets with the La Silla 3.6-m telescope. *Messenger* **110**, 9–14 (2002).

20. Pepe, F. et al. ESPRESSO at VLT. On-sky performance and first results. *Astron. Astrophys.* **645**, A96 (2021).

21. Faria, J. P., Santos, N. C., Figueira, P. & Brewer, B. J. kima: exoplanet detection in radial velocities. *J. Open Sour. Softw.* **3**, 487 (2018).

22. Zucker, S. & Alexander, T. Spectroscopic binary mass determination using relativity. *Astrophys. J.* **654**, L83–L86 (2007).

23. Konacki, M., Muterspaugh, M. W., Kulkarni, S. R. & Hetminiak, K. G. High-precision orbital and physical parameters of double-lined spectroscopic binary stars—HD78418, HD123999, HD160922, HD200077, and HD210027. *Astrophys. J.* **719**, 1293–1314 (2010).
24. Sybilski, P., Konacki, M., Kozłowski, S. K. & Hełminiak, K. G. Non-Keplerian effects in precision radial velocity measurements of double-line spectroscopic binary stars: numerical simulations. *Mon. Not. R. Astron. Soc.* **431**, 2024–2033 (2013).

25. Arras, P., Burkart, J., Quataert, E. & Weinberg, N. N. The radial velocity signature of tides raised in stars hosting exoplanets. *Mon. Not. R. Astron. Soc.* **422**, 1761–1766 (2012).

26. Hara, N. C., Unger, N., Delisle, J.-B., Díaz, R. & Ségransan, D. Detecting exoplanets with the false inclusion probability. Comparison with other detection criteria in the context of radial velocities. *Astron. Astrophys.* **663**, A14 (2022).

27. Hara, N. C., de Poyferré, T., Delisle, J.-B. & Hoffmann, M. A continuous multiple hypothesis testing framework for optimal exoplanet detection. Preprint at https://arXiv.org/abs/2203.04957 (2022).

28. Queloz, D. et al. No planet for HD 166435. *Astron. Astrophys.* **379**, 279–287 (2001).

29. Gomes da Silva, J. et al. Long-term magnetic activity of a sample of M-dwarf stars from the HARPS program. I. Comparison of activity indices. *Astron. Astrophys.* **534**, A30 (2011).

30. Schneider, J. On the occultations of a binary star by a circum-orbiting dark companion. *Planet. Space Sci.* **42**, 539–544 (1994).

31. Kostov, V. B. et al. Kepler-413b: a slightly misaligned, Neptune-size transiting circumbinary planet. *Astrophys. J.* **784**, 14 (2014).

32. Welsh, W. F. et al. Kepler 453 b – the 10th Kepler transiting circumbinary planet. *Astrophys. J.* **809**, 26 (2015).

33. Martin, D. V. & Triaud, A. H. Circumbinary planets – why they are so likely to transit. *Mon. Not. R. Astron. Soc.* **449**, 781–793 (2015).

34. Martin, D. V. Circumbinary planets – II. When transits come and go. *Mon. Not. R. Astron. Soc.* **465**, 3235–3253 (2017).

35. Pierens, A. & Nelson, R. P. On the migration of protoplanets embedded in circumbinary disks. *Astron. Astrophys.* **472**, 993–1001 (2007).

36. Pierens, A. & Nelson, R. P. On the evolution of multiple low mass planets embedded in a circumbinary disc. *Astron. Astrophys.* **478**, 939–949 (2008).

37. Penzlin, A. B. T., Kley, W., Audiffren, H. & Schäfer, C. M. Binary orbital evolution driven by a circumbinary disc. *Astron. Astrophys.* **660**, A101 (2022).

38. Chambers, J. E. A hybrid symplectic integrator that permits close encounters between massive bodies. *Mon. Not. R. Astron. Soc.* **304**, 793–799 (1999).

39. Chambers, J. E., Quintana, E. V., Duncan, M. J. & Lissauer, J. J. Symplectic integrator algorithms for modeling planetary accretion in binary star systems. *Astron. J.* **123**, 2884–2894 (2002).

40. Coleman, G. A. L. & Haworth, T. J. Dispersal of protoplanetary discs: how stellar properties and the local environment determine the pathway of evolution. *Mon. Not. R. Astron. Soc.* **514**, 2315–2332 (2022).

41. Lambrechts, M. & Johansen, A. Forming the cores of giant planets from the radial pebble flux in protoplanetary discs. *Astron. Astrophys.* **572**, A107 (2014).

42. Paardekooper, S.-J., Baruteau, C. & Kley, W. A torque formula for non-isothermal Type I planetary migration – II. Effects of diffusion. *Mon. Not. R. Astron. Soc.* **410**, 293–303 (2011).

43. Coleman, G. A. L., Papaloizou, J. C. B. & Nelson, R. P. In situ accretion of gaseous envelopes on to planetary cores embedded in evolving protoplanetary discs. *Mon. Not. R. Astron. Soc.* **470**, 3206–3219 (2017).

44. Poon, S. T. S., Nelson, R. P. & Coleman, G. A. L. In situ formation of hot Jupiters with companion super-Earths. *Mon. Not. R. Astron. Soc.* **505**, 2500–2516 (2021).

45. Martin, D. V., Mazeh, T. & Fabrycky, D. C. No circumbinary planets transiting the tightest Kepler binaries – a possible fingerprint of a third star. *Mon. Not. R. Astron. Soc.* **453**, 3554–3567 (2015).

46. Kunovac Hodžić, V. et al. The EBLM project – VII. Spin-orbit alignment for the circumbinary planet host EBLM J0608-59 A/TOI-1338 A. *Mon. Not. R. Astron. Soc.* **497**, 1627–1633 (2020).

47. Mayor, M. et al. Setting new standards with HARPS. *Messenger* **114**, 20–24 (2003).

48. Baranne, A. et al. ELODIE: a spectrograph for accurate radial velocity measurements. *Astron. Astrophys. Suppl. Ser.* **119**, 373–390 (1996).

49. Standing, M. R. et al. BEBOP II: sensitivity to sub-Saturn circumbinary planets using radial-velocities. *Mon. Not. R. Astron. Soc.* **511**, 3571–3583 (2022).

50. Triaud, A. H. M. J. in *Handbook of Exoplanets* (eds Deeg, H. J. & Belmonte, J. A.) 1375–1401 (Springer Cham, 2018).

51. Kunovac Hodžić, V., Triaud, A. H. M. J., Cegla, H. M. Chaplin, W. J. & Davies, G. R. Orbital misalignment of the super-Earth π Men c with the spin of its star. *Mon. Not. R. Astron. Soc.* **502**, 2893–2911 (2021).

52. Faria, J. P. et al. A candidate short-period sub-Earth orbiting Proxima Centauri. *Astron. Astrophys.* **658**, A115 (2022).

53. Agol, E. et al. Refining the transit-timing and photometric analysis of TRAPPIST-1: masses, radii, densities, dynamics, and ephemerides. *Planet. Sci.* **2**, 1 (2021).

54. Hogg, D. W., Bovy, J. & Lang, D. Data analysis recipes: fitting a model to data. Preprint at https://arXiv.org/abs/1008.4686 (2010).

55. Salvatier, J., Wiecki, T. & Fonnesbeck, C. Probabilistic programming in Python using PyMC. Preprint at https://arXiv.org/abs/1507.08050 (2015).

56. Brewer, B. J., Pártay, L. B. & Csányi, G. Diffusive nested sampling. *Stat. Comput.* **21**, 649–656 (2011).

57. Feroz, F., Balan, S. T. & Hobson, M. P. Detecting extrasolar planets from stellar radial velocities using Bayesian evidence. *Mon. Not. R. Astron. Soc.* **415**, 3462–3472 (2011).

58. Brewer, B. J. Inference for trans-dimensional Bayesian models with diffusive nested sampling. Preprint at https://arXiv.org/abs/1411.3921 (2014). 1411.3921.

59. Baycroft, T. A., Triaud, A. H. M. J., Faria, J., Correia, A. C. M. & Standing, M. R. Improving circumbinary planet detections by fitting their binary's apsidal precession. *Mon. Not. R. Astron. Soc.* **521**, 1871–1879 (2023).

60. Kass, R. E. & Raftery, A. E. Bayes factors. *J. Am. Stat. Assoc.* **90**, 773–795 (1995).

61. Trotta, R. Bayes in the sky: Bayesian inference and model selection in cosmology. *Contemp. Phys.* **49**, 71–104 (2008).

62. Jeffreys, H. *Theory of Probability* (Oxford University Press, 1961).

63. McInnes, L., Healy, J. & Astels, S. hdbscan: hierarchical density based clustering. *J. Open Sour. Softw.* **2**, 205 (2017).

64. Foreman-Mackey, D. corner.py: scatterplot matrices in python. *J. Open Sour. Softw.* **1**, 24 (2016).

65. Tuomi, M., Jones, H. R. A., Barnes, J. R., Anglada-Escudé, G. & Jenkins, J. S. Bayesian search for low-mass planets around nearby M dwarfs – estimates for occurrence rate based on global detectability statistics. *Mon. Not. R. Astron. Soc.* **441**, 1545–1569 (2014).

66. Gomes da Silva, J., Figueira, P., Santos, N. & Faria, J. ACTIN: a tool to calculate stellar activity indices. *J. Open Sour. Softw.* **3**, 667 (2018).

67. Zechmeister, M., Kürster, M. & Endl, M. The M dwarf planet search programme at the ESO VLT + UVES. A search for terrestrial planets in the habitable zone of M dwarfs. *Astron. Astrophys.* **505**, 859–871 (2009).

68. Pojmanski, G. The All Sky Automated Survey. *Acta Astron.* **47**, 467–481 (1997).

69. Quarles, B., Satyal, S., Kostov, V., Kaib, N. & Haghighipour, N. Stability limits of circumbinary planets: is there a pile-up in the Kepler CBPs? *Astrophys. J.* **856**, 150 (2018).
70. Laskar, J. The chaotic motion of the solar system : a numerical estimate of the size of the chaotic zones. *Icarus* **88**, 266–291 (1990).

71. Laskar, J. Frequency analysis for multi-dimensional systems. Global dynamics and diffusion. *Phys. D* **67**, 257–281 (1993).

72. Correia, A. C. M. et al. The CORALIE survey for southern extra-solar planets. XIII. A pair of planets around HD 202206 or a circumbinary planet? *Astron. Astrophys.* **440**, 751–758 (2005).

73. Laskar, J. & Robutel, P. High order symplectic integrators for perturbed Hamiltonian systems. *Celest. Mech. Dynam. Astron.* **80**, 39–62 (2001).

74. Farago, F. & Laskar, J. High-inclination orbits in the secular quadrupolar three-body problem. *Mon. Not. R. Astron. Soc.* **401**, 1189–1198 (2010).

75. Shakura, N. I. & Sunyaev, R. A. Black holes in binary systems. Observational appearance. *Astron. Astrophys.* **24**, 337–355 (1973).

76. D'Angelo, G. & Marzari, F. Outward migration of Jupiter and Saturn in evolved gaseous disks. *Astrophys. J.* **757**, 50 (2012).

77. Dullemond, C. P., Hollenbach, D., Kamp, I. & D'Alessio, P. in *Protostars and Planets V* (eds Reipurth, B. et al.) 555–572 (University of Arizona Press, 2007).

78. Matsuyama, I., Johnstone, D. & Hartmann, L. Viscous diffusion and photoevaporation of stellar disks. *Astrophys. J.* **582**, 893–904 (2003).

79. Pierens, A. & Nelson, R. P. Migration and gas accretion scenarios for the Kepler 16, 34, and 35 circumbinary planets. *Astron. Astrophys.* **556**, A134 (2013).

80. Thun, D., Kley, W. & Picogna, G. Circumbinary discs: numerical and physical behaviour. *Astron. Astrophys.* **604**, A102 (2017).

81. Benítez-Llambay, P. & Masset, F. S. FARGO3D: a new GPU-oriented MHD code. *Astrophys. J. Suppl. Ser.* **223**, 11 (2016).

82. Coleman, G. A. L. & Nelson, R. P. On the formation of planetary systems via oligarchic growth in thermally evolving viscous discs. *Mon. Not. R. Astron. Soc.* **445**, 479–499 (2014).

83. Daisaka, J. K., Tanaka, H. & Ida, S. Orbital evolution and accretion of protoplanets tidally interacting with a gas disk. II. Solid surface density evolution with type-I migration. *Icarus* **185**, 492–507 (2006).

84. Cresswell, P. & Nelson, R. P. Three-dimensional simulations of multiple protoplanets embedded in a protostellar disc. *Astron. Astrophys.* **482**, 677–690 (2008).

85. Paardekooper, S.-J., Baruteau, C. & Kley, W. A torque formula for non-isothermal type I planetary migration – II. Effects of diffusion. *Mon. Not. R. Astron. Soc.* **410**, 293–303 (2011).

86. Crida, A., Morbidelli, A. & Masset, F. On the width and shape of gaps in protoplanetary disks. *Icarus* **181**, 587–604 (2006).

87. Lin, D. N. C. & Papaloizou, J. On the tidal interaction between protoplanets and the protoplanetary disk. III. Orbital migration of protoplanets. *Astrophys. J.* **309**, 846 (1986).

88. Lambrechts, M. & Johansen, A. Rapid growth of gas-giant cores by pebble accretion. *Astron. Astrophys.* **544**, A32 (2012).

89. Johansen, A. & Lambrechts, M. Forming planets via pebble accretion. *Annu Rev. Earth Planet Sci.* **45**, 359–387 (2017).

90. Coleman, G. A. L., Nelson, R. P. & Triaud, A. H. M. J. Dusty circumbinary discs: inner cavity structures and stopping locations of migrating planets. *Mon. Not. R. Astron. Soc.* **513**, 2563–2580 (2022).

91. Coleman, G. A. L. From dust to planets – I. Planetesimal and embryo formation. *Mon. Not. R. Astron. Soc.* **506**, 3596–3614 (2021).

92. Holman, M. J. & Wiegert, P. A. Long-term stability of planets in binary systems. *Astron. J.* **117**, 621–628 (1999).

93. Kane, S. R. & Hinkel, N. R. On the habitable zones of circumbinary planetary systems. *Astrophys. J.* **762**, 7 (2013).

94. Cutri, R. M. et al. 2MASS All-Sky Catalog of Point Sources (Cutri+ 2003) no. II/246 (VizieR Online Data Catalog, 2003).

95. Müller, T. W. A. & Haghighipour, N. Calculating the habitable zones of multiple star systems with a new interactive web site. *Astrophys. J.* **782**, 26 (2014).

96. Akeson, R. L. et al. The NASA exoplanet archive: data and tools for exoplanet research. *Pub. Astron. Sos. Pac.* **125**, 989 (2013).

97. Astropy Collaboration et al. The Astropy Project: sustaining and growing a community-oriented open-source project and the latest major release (v5.0) of the core package. *Astrophys. J.* **935**, 167 (2020).

98. Harris, C.R. et al. Array programming with NumPy. *Nature* **585**, 357–362 (2020).

99. McKinney, W. Data structures for statistical computing in Python. In *Proc. 9th Python in Science Conference* (eds van der Walt, S. & Millman, J.) 56–61 (SciPy, 2010).

100. Pauli, V. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. *Nat. Meth.* **17**, 261–272 (2020).

101. Hunter, J. D. Matplotlib: a 2D graphics environment. *Comput. Sci. Eng.* **9**, 90–95 (2007).

## Acknowledgements

We would like to thank the ESO staff at La Silla and Paranal for their continued support throughout this work, especially through the COVID-19 pandemic, with special thanks to our ESPRESSO support astronomer M. Wittkowski. We also thank all of the observers who took part in the HARPS timeshare and were instrumental in collecting data for this project. We particularly thank X. Dumusque and F. Bouchy for their work in organizing the timeshare. This Article is based on observations collected at the European Southern Observatory under ESO programmes 103.2024, 106.216B, 1101.C-0721 and 106.212H. This research has made use of the services of the ESO Science Archive Facility. A.H.M.J.T. 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 the Leverhulme Trust (research project no. RPG-2018-418) to conduct this research. M.R.S. would like to acknowledge the support of the UK Science and Technology Facilities Council (STFC) under grant number ST/T000295/1. A.C.M.C. acknowledges support from CFISUC (grant nos. UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (grant no. PTDC/FIS-AST/7002/2020), PHOBOS (grant no. POCI-01-0145-FEDER-029932) and ENGAGE SKA (grant no. POCI-01-0145-FEDER-022217), funded by COMPETE 2020 and FCT, Portugal. The stability maps were performed at the OBLIVION Supercomputer (HPC Center – University of Évora), funded by ENGAGE SKA and by the BigData@UE project (grant no. ALT20-03-0246-FEDER-000033). Support for D.V.M. was provided by NASA through the NASA Hubble Fellowship grant no. HF2-51464 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract no. NAS5-26555. D.V.M. is a NASA Sagan Fellow. V.K. acknowledges support from NSF award no. AST2009501. A.S. received funding from the French government under the 'France 2030' investment plan managed by the French National Research Agency (reference nos. ANR-16-CONV-000X/ANR-17-EURE-00XX) and from Excellence Initiative of Aix-Marseille University-A*MIDEX (reference no. AMX-21-IET-018). This work was supported by the 'Programme National de Planétologie' (PNP) of CNRS/INSU. J.P.F. received support in the form of a work contract funded by national funds through Fundação para a Ciência e a Tecnologia (FCT) with reference no. DL57/2016/CP1364/CT0005. A.C.C. acknowledges support from STFC consolidated grant nos. ST/R000824/1 and ST/V000861/1, and UKSA grant no. ST/R003203/1. M.G. is an FNRS Senior Research Associate. R.P.N. and G.A.L.C. utilized Queen Mary's Apocrita HPC facility, supported by QMUL Research-IT (https://doi.org/10.5281/zenodo.438045). This work was performed using the DiRAC Data Intensive service at
Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grant nos. ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant no. ST/R001014/1. DiRAC is part of the National e-Infrastructure.

## Author contributions

The BEBOP project was established by A.H.M.J.T. and D.V.M., building on work by the Wide Angle Search for Planets consortium that involved A.H.M.J.T., I.B., A.C.C., M.G., C.H., P.F.L.M., F.P., A.S. and S.U. The radial-velocity data used in this manuscript were secured and obtained by M.R.S., A.H.M.J.T., L.S., D.V.M., V.K., D.S., T.A.B., P.F.L.M., N.J.M. and A.S., and their reduction involved M.R.S., A.H.M.J.T., F.P. and S.U. The observational campaign was led by M.R.S. for ESPRESSO observations and A.H.M.J.T. for HARPS. The radial-velocity analysis was led by M.R.S. with assistance from A.H.M.J.T., D.V.M. and J.P.F. Outliers in the radial-velocity data were identified by T.A.B. and M.R.S. The detection limit analysis was performed by M.R.S., J.H. and E.L. with help from J.P.F. The TIP and FIP analysis was carried out by T.A.B. and N.C.H. The analysis of the stellar activity was carried out by L.S., M.R.S. and A.H.M.J.T. An independent stability analysis on the system was carried out by A.C.M.C. and D.V.M. with input from R.M. Formation pathway simulations were carried out by G.A.L.C. and R.P.N. The TSM was calculated by G.D. The system was determined to have an inner transiting planet TOI-1338 b by D.V.M., J.A.O. and W.F.W. V.K. prepared the majority of figures in the Article. All co-authors assisted in writing and reviewing the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41550-023-01948-4.

**Correspondence and requests for materials** should be addressed to Matthew R. Standing.

**Peer review information** *Nature Astronomy* thanks Maciej Konacki and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2023

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

e-mail: matthew.standing@open.ac.uk

<footnote n="1">¹School of Physics and Astronomy, University of Birmingham, Birmingham, UK. ²School of Physical Sciences, The Open University, Milton Keynes, UK.</footnote>

<footnote n="3">³Department of Astronomy, The Ohio State University, Columbus, OH, USA. ⁴CFisUC, Departamento de Física, Universidade de Coimbra, Coimbra, Portugal. ⁵IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Université, Paris, France. ⁶Astronomy Unit, Queen Mary University of London, London, UK. ⁷Lowell Observatory, Flagstaff, AZ, USA. ⁸Department of Astronomy and Planetary Science, Northern Arizona University, Flagstaff, AZ, USA. ⁹Aix Marseille Université, CNRS, CNES, Institut Origines, LAM, Marseille, France. ¹⁰Centre for Exoplanet Science/SUPA, School of Physics and Astronomy, University of St Andrews, St Andrews, UK. ¹¹Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Porto, Portugal. ¹²Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Porto, Portugal. ¹³Astrobiology Research Unit, University of Liège, Liège, Belgium.</footnote>

<footnote n="14">¹⁴Observatoire Astronomique de l'Université de Genève, Versoix, Switzerland. ¹⁵Astrophysics Group, Keele University, Keele, UK. ¹⁶School of Physics and Astronomy, Monash University, Monash, Victoria, Australia. ¹⁷Department of Astronomy, San Diego State University, San Diego, CA, USA.</footnote>