# BEBOP VII. SOPHIE discovery of BEBOP-3b, a circumbinary giant planet on an eccentric orbit

Thomas A Baycroft, Alexandre Santerne, Amaury H M J Triaud, Neda

Heidari, Daniel Sebastian, Yasmin T Davis, Alexandre C M Correia, Lalitha

Sairam, Alix V Freckelton, Aleyna Adamson, et al.

► To cite this version:

Thomas A Baycroft, Alexandre Santerne, Amaury H M J Triaud, Neda Heidari, Daniel Sebastian, et al. BEBOP VII. SOPHIE discovery of BEBOP-3b, a circumbinary giant planet on an eccentric orbit. Monthly Notices of the Royal Astronomical Society, 2025, 541 (3), pp.2801-2814. 10.1093/mn-ras/staf1184. hal-05469759

HAL Id: hal-05469759

https://hal.science/hal-05469759v1

Submitted on 21 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# BEBOP VII. SOPHIE discovery of BEBOP-3b, a circumbinary giant planet on an eccentric orbit

Thomas A. Baycroft<sup>1,★</sup> Alexandre Santerne<sup>2</sup> Amaury H. M. J. Triaud<sup>1</sup> Neda Heidari<sup>3</sup>

Daniel Sebastian<sup>1</sup> Yasmin T. Davis<sup>1</sup> Alexandre C. M. Correia<sup>4,5</sup> Lalitha Sairam<sup>6</sup>

Alix V. Freckelton<sup>1</sup> Aleyna Adamson<sup>1</sup> Isabelle Boisse<sup>2</sup> Gavin A. L. Coleman<sup>7</sup>

Georgina Dransfield,<sup>8,9</sup> João Faria<sup>10</sup> Salomé Grouffal,<sup>2</sup> Nathan Hara,<sup>2</sup> Guillaume Hébrard,<sup>3</sup>

Vedad Kunovac,<sup>11,12</sup> David V. Martin<sup>13</sup> Pierre F. L. Maxted<sup>14</sup> Richard P. Nelson,<sup>7</sup>

Madison G. Scott<sup>1</sup> Owen J. Scutt<sup>1</sup> and Matthew R. Standing<sup>15</sup>

Accepted 2025 July 16. Received 2025 July 14; in original form 2025 April 25

## ABSTRACT

Planetary systems orbiting close binaries are valuable testing grounds for planet formation and migration models. More detections with good mass measurements are needed. We present a new planet discovered during the BEBOP survey for circumbinary exoplanets using radial velocities. We use data taken with the SOPHIE spectrograph at the Observatoire de Haute-Provence, and perform a spectroscopic analysis to obtain high precision radial velocities. This planet is the first radial velocity detection of a previously unknown circumbinary system. The planet has a mass of 0.56 M<sub>Jup</sub> and orbits its host binary in 550 d with an eccentricity of 0.25. Compared to most of the previously known circumbinary planets, BEBOP-3b has a long period (relative to the binary) and a high eccentricity. There also is a candidate outer planet with a ~ 1400 d orbital period. We test the stability of potential further candidate signals inside the orbit of BEBOP-3b, and demonstrate that there are stable orbital solutions for planets near the instability region which is where the *Kepler* circumbinary planets are located. We also use our data to obtain independent dynamical masses for the two stellar components of the eclipsing binary using high resolution cross-correlation spectroscopy, and compare those results to a more traditional approach, finding them compatible with one another.

**Key words:** techniques: radial velocities – planets and satellites: detection – binaries: eclipsing – stars: individual: BEBOP-3.

## 1 INTRODUCTION

Circumbinary exoplanetary systems are a useful test of exoplanet formation and migration models, not only for binary stars but also applicable to single stars. Far from the binary the dynamical influence of the inner binary is negligible and planet formation is expected to

be comparable to around single stars. In the inner region of the circumbinary disc, the influence of the binary causes drastically different conditions for nearby bodies compared to single stars; dynamical perturbations induced by the binary cause collisions to be destructive, impairing the growth of planetesimals into protoplanets (Meschiari 2012). This has been shown to mean that circumbinary planet formation is inhibited in a region near the binary (Moriwaki & Nakagawa 2004; Paardekooper et al. 2012; Pierens, Nelson & McNally 2021), however this region is precisely where almost all

<footnote n="1"><sup>1</sup>School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="2"><sup>2</sup>Aix Marseille Univ, CNRS, CNES, LAM, F-13284 Marseille, France</footnote>

<footnote n="3"><sup>3</sup>Institut d'astrophysique de Paris, UMR 7095 CNRS université Pierre et Marie Curie, 98 bis, boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="4"><sup>4</sup>CFisUC, Departamento de Fisica, Universidade de Coimbra, P-3004-516 Coimbra, Portugal</footnote>

<footnote n="5"><sup>5</sup>IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Université, 77 avenue Denfert-Rochereau, F-75014 Paris, France</footnote>

<footnote n="6"><sup>6</sup>Institute of Astronomy, University of Cambridge, Madingley road, Cambridge CB3 0HA, UK</footnote>

<footnote n="7"><sup>7</sup>Astronomy Unit, Department of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK</footnote>

<footnote n="8"><sup>8</sup>Department of Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="9"><sup>9</sup>Magdalen College, University of Oxford, Oxford OX1 4AU, UK</footnote>

<footnote n="10"><sup>10</sup>Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="11"><sup>11</sup>Department of Physics, University of Warwick, Coventry CV4 7AL, UK</footnote>

<footnote n="12"><sup>12</sup>Centre for Exoplanets and Habitability, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK</footnote>

<footnote n="13"><sup>13</sup>Department of Physics and Astronomy, Tufts University, 574 Boston Avenue, Medford, MA 02155, USA</footnote>

<footnote n="14"><sup>14</sup>Astrophysics Group, Keele University, Keele, Staffordshire ST5 5BG, UK</footnote>

<footnote n="15"><sup>15</sup>European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, Villanueva de la Cañada, E-28692 Madrid, Spain</footnote>

<footnote>*E-mail: t.baycroft@bham.ac.uk, thomasbaycroftastro@gmail.com</footnote>
of the discovered circumbinary planets have been found to exist. These circumbinary exoplanets are therefore likely to have migrated to their current locations, meaning they are an excellent probe of migration theories in general (Pierens & Nelson 2008; Coleman, Nelson & Triaud 2023; Coleman et al. 2024).

There are currently 16 confirmed circumbinary exoplanets orbiting main-sequence binaries. Of these, 14 were discovered using the transit method (e.g. Doyle et al. 2011), 1 using eclipse timing variations (Goldberg et al. 2023), and 1 using the radial velocity method (Standing et al. 2023). Of these, few are significantly eccentric (the most eccentric is Kepler-34 b with $e \approx 0.18$; Welsh et al. 2012)<fnref n="1" /> and there are two multiplanet systems: TOI-1338/BEBOP-1 (Kostov et al. 2020; Standing et al. 2023) and Kepler-47 (Orosz et al. 2012, 2019).

Of the transiting circumbinary systems, all but one (Kepler-1647 b; Kostov et al. 2016) have their innermost planet orbiting very close to the boundary inside which orbits are unstable due to the dynamical influence of the binary (Holman & Wiegert 1999). This so-called 'pile-up' is thought to not be solely due to the observational bias of the transit method towards short-period planets (Martin & Triaud 2014; Li, Holman & Tao 2016). The pile-up has been interpreted as a result of migration through a circumbinary disc, and parking around the current location near where the disc gets truncated (Penzlin, Kley & Nelson 2021), and in fact was predicted before being observed (Pierens & Nelson 2008).

The BEBOP (Binaries Escorted By Orbiting Planets) survey is a radial velocity survey aiming to detect circumbinary planets (Martin et al. 2019). Its first phase had a focus on single-lined binaries – systems in which there are no obvious spectral lines observed from the secondary star. For a main-sequence binary this typically requires a mass ratio $q \lesssim 0.3$ (Triaud et al. 2017). So far BEBOP has made three circumbinary planet detections, with BEBOP-1c (Standing et al. 2023) being the first new planet. The other two detections were confirmations of Kepler-16b (Triaud et al. 2022) and of TIC 172900988b (Sairam et al. 2024). This paper, the seventh in the BEBOP series, presents our second circumbinary planet discovery using the radial velocity method, and the first such planet with a significant eccentricity measurement, BEBOP-3b. This newly detected circumbinary planet orbits the eclipsing binary star BD+79 230 (TIC 289949453).

There is a disadvantage to observing single-lined binaries instead of double-lined binaries. In double-lined binaries, the two sets of radial velocities mean that a model-independent dynamical mass can be derived. In single-lined systems, however, the stellar (and thus the planetary) physical parameters are ultimately dependent on the modelled parameters for the primary star. Recently, Sebastian et al. (2024b) demonstrated the use of high resolution cross-correlation spectroscopy (HRCCS) on single-lined binaries. This allows the signal of the secondary within the spectra to be retrieved, effectively converting the single-lined binary to a double-lined binary for the purposes of absolute mass determination. Sebastian et al. (2024b) apply the method to the TOI-1338/BEBOP-1 system using spectra from the HARPS and ESPRESSO spectrographs. They obtain a model-independent and more precise mass measurement for both stars in the binary. We repeat the same analysis to obtain dynamical masses for the two components of the BEBOP-3(AB) binary.

In Section 2 we describe the data collection and reduction pipeline. The various analyses performed are described in Section 3. We present the results of the radial velocity analysis and a stellar activity analysis in Section 4. We discuss the results and investigate the orbital stability of the system in Section 5, and conclude in Section 6.

## 2 RADIAL VELOCITY DATA

The sample selection for the southern component of the BEBOP survey is reported in Martin et al. (2019), but the northern sample has not been described much yet. The southern sample was built from a WASP (Wide Angle Search for Planets; Pollacco et al. 2006) false positives flag known as EBLM (Eclipsing Binary–Low Mass; Triaud et al. 2013; Maxted, Triaud & Martin 2023). Typically these are Sun-like stars with a stellar companion < 0.3 $M_\odot$, producing transit-shaped primary eclipses. Whilst EBLM-flagged false positives exist for the Northern component of WASP, the record is not as thorough as it was for the south (Triaud 2011). In the south, the survey benefited from a dedicated spectroscopic follow-up, but in the north, WASP depended on facilities with competitive access and not all false positives were followed-up enough to establish what type they were. This is why the northern BEBOP sample was selected not just from WASP, but also from the KELT survey's (Kilo-degree Extremely Little Telescope; Pepper et al. 2007) false positives (Collins et al. 2018). Systems selected are above the celestial equator, have orbital periods $P_{\text{bin}} > 5$ d, and primary spectral types > F4. The restriction to relatively long period eclipsing binaries was driven by obtainable precision; tight binaries are tidally locked and stars in < 5 d binaries rotate quickly resulting in the spectral lines being so broadened that finding CBPs is very difficult. Coincidentally, the known transiting CBPs are all around similarly 'long-period' binaries, potentially as a quirk of their formation process (Martin, Mazeh & Fabrycky 2015; Muñoz & Lai 2015). Only six systems were common between WASP-North and KELT. Following a short reconnaissance programme in 2018,<fnref n="2" /> we created a top priority subsample that became the BEBOP-North sample, where SOPHIE spectra typically produced radial velocity uncertainties < 10 m s⁻<fnref n="1" /> in 1800s or under. The BEBOP-3 system was identified from the KELT false positive catalogue (Collins et al. 2018).

Overall, 141 high resolution spectra were obtained between 2018 November 13 and 2024 November 18 with the high resolution échelle spectrograph SOPHIE, mounted on the T193 at the Observatoire de Haute-Provence (Perruchot et al. 2008). Data were acquired under programme ID 18B.DISC.TRIA and 19A.PNP.SANT.

The data were acquired in ObjAB mode where fibre A is placed on the target and fibre B is placed on the sky to enable the removal of solar spectral lines reflected by the Moon. The observations were taken in high resolution (HR) mode with resolution $R = 75\,000$.<fnref n="3" /> Standard calibrations were taken at the start of the night and at approximately 2 h intervals throughout to track the instrumental zero-point error.

The SOPHIE radial velocities (RVs) are computed by the SOPHIE Data Reduction System (DRS; Bouchy et al. 2009b), which cross-correlates the spectra with a weighted numerical mask corresponding to a G2-type star<fnref n="4" /> and subsequently fits a Gaussian to the cross-correlation function (CCF) (Baranne et al. 1996; Pepe et al. 2002).

<footnote n="1">¹The HD202206 system was discovered through radial-velocity (Correia et al. 2005; Benedict & Harrison 2017) but it is uncertain between a circumbinary brown dwarf and a planet orbiting a Star-Brown dwarf binary, the eccentricity of the 'planet' is ≈ 0.22.</footnote>

<footnote n="2">²Prog.ID 18B.DISC.TRIA</footnote>

<footnote n="3">³At 5550 Å.</footnote>

<footnote n="4">⁴The BEBOP-3 primary star has a F9 spectral classification, G2 is the most suitable mask within the SOPHIE pipeline.</footnote>
<table><caption>**Table 1.** Stellar and binary parameters from three different analysis methods.</caption><thead><tr><th>Parameter</th><th>HRCCS</th><th>RV + Photometry</th><th>Spectral</th></tr></thead><tbody><tr><td><i>i</i><sub>bin</sub>(°)</td><td></td><td>88.7488 ± 0.0092</td><td></td></tr><tr><td><i>M</i><sub>pri</sub> (M<sub>☉</sub>)</td><td>1.083 ± 0.026</td><td>1.163 ± 0.043</td><td>1.131±0.051*</td></tr><tr><td><i>M</i><sub>sec</sub> (M<sub>☉</sub>)</td><td>0.2615 ± 0.0039</td><td>0.2731<sup>+0.0065</sup><sub>-0.0064</sub></td><td></td></tr><tr><td><i>R</i><sub>pri</sub> (R<sub>☉</sub>)</td><td>–</td><td>1.414 ± 0.017**</td><td>1.4107±0.0089*</td></tr><tr><td><i>R</i><sub>sec</sub> (R<sub>☉</sub>)</td><td>–</td><td>0.2811 ± 0.0033</td><td></td></tr><tr><td><i>T</i><sub>eff,pri</sub> (K)</td><td>–</td><td>–</td><td>6033 ± 100</td></tr><tr><td>log <i>g</i><sub>pri</sub> (dex)</td><td>–</td><td>4.203 ± 0.017</td><td>4.27 ± 0.09</td></tr><tr><td>[Fe/H] (dex)</td><td>–</td><td>–</td><td>-0.02 ± 0.10</td></tr><tr><td><i>v sin i</i> (km s<sup>-1</sup>)</td><td>–</td><td>–</td><td>2.00 ± 0.66</td></tr></tbody></table>

* from ISOCHRONE fitting

** from GAIA

To enhance the precision of the SOPHIE RV measurements, we applied the optimized procedures described in Heidari (2022) and Heidari et al. (2024), including: (1) corrections for charge transfer inefficiency in the CCD (Bouchy et al. 2009a); (2) correction of template (colour) effects (Heidari et al. 2024); (3) removal of moonlight contamination by utilizing the simultaneous sky spectrum recorded via the second SOPHIE fibre (Pollacco et al. 2008); (4) correction for long-term instrumental drifts using the RV master constant time-series, constructed from the so-called ‘constant’ stars monitored each night by SOPHIE (Courcol et al. 2015; Heidari 2022; Heidari et al. 2024); and (5) correction for short-term instrumental drifts by interpolating drift measurements at the exact time of each observation, based on frequently measured drifts throughout the night. Radial velocity data are collated into Table A1.

The 141 spectra obtained have a median SNR<fnref n="5" /> of 72.9 and median radial velocity precision of 2.3 m s⁻¹ Radial velocity outliers exist. These are dealt with in various ways for the different analyses. Outliers identified are flagged in Table A1.

# 3 ANALYSIS AND METHODS

## 3.1 Spectral analysis

In Freckelton et al. (2024), the 91 SOPHIE spectra available at the time were coadded and the combined spectrum analysed using the PAWS pipeline. This analysis achieved a signal-to-noise ratio of 555. Based on this signal-noise-ratio being already high enough that model uncertainty is dominant, we do not repeat the analysis from Freckelton et al. (2024) with extra spectra. Atmospheric parameters for the primary star were determined by applying the PAWS pipeline to this coadded spectrum. The pipeline implements the curve-of-growth equivalent widths (EW) and spectral synthesis methods via iSpec (Blanco-Cuaresma 2019). The advantage of combining the two methods sequentially is that the parameters $T_{eff}$ (effective temperature), logg (surface gravity), [Fe/H] (iron abundance), and $v \sin i$ (projected rotational velocity), as presented in Table 1, are all allowed to freely vary to fit the observations in at least one stage of the pipeline. Using the WIDTH radiative transfer code (Sbordone et al. 2004), the EW method was used to determine an initial set of atmospheric parameters for the primary

star. PAWS then applies the spectral synthesis method to allow a determination of $v \sin i$, using the initial atmospheric parameters as the start point for the first iteration. Synthetic spectra were generated using these parameters by the use of the SPECTRUM (Gray & Corbally 1994) radiative transfer code, with the parameters iterated until the synthetic spectrum sufficiently reproduced the observed one, which returned the final atmospheric parameters shown in Table 1. Homogeneity between the two methods was ensured by employing the SPECTRUM (Gray & Corbally 1994) line list and ATLAS (Kurucz 2005) model atmosphere set throughout the pipeline.

## 3.2 HRCCS measurement of the secondary component

By combining the large number of SOPHIE observations, we can detect the extremely small contribution of the secondary star (~ 0.1 per cent for the spectral range of SOPHIE). This in effect turns the single-lined binary into a double-lined binary, which allows us to derive model-independent dynamical masses for both stars. Here we apply the high-resolution cross-correlation spectroscopy method (HRCCS; Snellen et al. 2010), which is commonly used to detect exoplanet atmospheres. It has been adapted and fully validated to high-contrast binaries (e.g. TOI-1338/BEBOP-1 and Kepler-16 AB; Sebastian et al. 2024b, 2025). The method includes a two-step process. First the absorption lines of the primary star are effectively removed from the spectra, second the residual spectra are cross-correlated and combined in the rest-frame of the secondary star. This enhances the SNR of its absorption lines and allows its detection. We follow closely the recipe outlined in Sebastian et al. (2024b).

We exclude 8 of the 141 SOPHIE spectra from our analysis. Three of which correspond to a wrong pointing and five with very low SNR (<23.5 at 500 nm), which is 3 times lower than the average SNR (71) for the primary. SOPHIE spectra excluded from the HRCCS analysis are indicated in Table A1. The 133 remaining spectra have an average SNR of 73.9 per resolution element at 500 nm. This means for the faint secondary, the SNR is as low as 0.08 only.

First we split the 1D SOPHIE spectra (s1d) each into 70 chunks of 4387 pixels. We then shift the spectra into the rest frame of the primary, and perform continuum normalization, sigma clipping, and singular-value decomposition (SVD) detrending (Kalman 1996). The latter will effectively remove the spectral lines of the primary star from the series of spectra. The number of excluded SVD components is determined automatically using the effective matrix rank (Roy & Vetterli 2007) for each chunk. We set the maximum number of components to 26 and find that the median is 5. This results in a reduced rank of 0.04 (the number of removed components divided by the number of spectra). According to table 4 in Sebastian et al. (2024b), we expect that this detrending has removed ≤ 3 per cent of the secondary’s spectral features. The secondary’s signal in the data is therefore kept mainly intact.

We exclude a further 11 of the 133 spectra at orbital phases, where the velocity of the secondary’s signal differs less than 10 km s⁻¹ from the primary’s absorption lines. In this way we make sure that possible residuals of the detrending do not interfere with our measurements. We also discharge the 20 bluest chunks (keeping data from 477 nm and redder) as those add significant noise to the secondary signal. We then use the k-focusing method (Sebastian, Triaud & Brogi 2024a), which enhances the cross-correlation function (CCF) of the secondary star in the secondary’s rest frame by scanning different semi-amplitudes ($K_2$). This method produces a $K_2 - V_{rest}$ plane by combining all CCFs of the detrended data with an M2-dwarf

<footnote n="5">⁵The signal-to-noise ratio quoted here is that achieved at the centre of order 26 (out of orders 0–38), corresponding to 5550 Å (near the center of the V band)</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** HRCCS 5.5-$\sigma$ detection of the secondary. Upper panel: Cross-correlation map and SALTIRE model for SOPHIE data. Dotted lines show the positions of slices (a,b). Lower panel: Slices a & b through CCF map at maximum significance. Solid line: best-fitting SALTIRE model, shaded area: $1\sigma$ Uncertainties from the MCMC. Error bars represent the MCMC jitter term ($\sigma_{jit}$) from the 2D fit.</figcaption></figure>

line mask.<fnref n="6" /> We sample the map for a velocity range between -5 and 65 km s⁻¹. This range includes the systemic velocity of the binary, measured from the primary star. We choose a step width of 1.5 km s⁻¹, which takes into account the pixel resolution of the spectrograph. The unknown semi-amplitude of the secondary ($K_2$) is sampled from 55 – 105 km s⁻¹ in steps of 1.5 km s⁻¹, making sure the expected semi-amplitude is well covered.

The resulting CCF map in the $K_2 - V_{rest}$ plane is shown in Fig. 1. The CCF signal of the secondary is detected with a 5.5-$\sigma$ significance at the expected semi-amplitude and systemic velocity. We measure this signal using the SALTIRE model (Sebastian et al. 2024a), which accurately reproduces the CCF signal by fitting a k-focused double-Gaussian function to the data. We estimate the best-fitting parameters for $K_2$ and $V_{rest}$ by sampling the data using Markov Chain Monte Carlo (MCMC) with a posterior distribution of 21 000 samples. In Table 2, we present the best parameters, given as the 50th percentile of the MCMC sampling. We then follow Sebastian et al. (2025) to estimate the systematic uncertainties and randomly draw 50 per cent of the data and measure the best-fitting parameters. This is repeated 200 times and the systematic uncertainties are derived as the standard error of the 200 measurements. Results are presented in Section 4.1 and Tables 1 and 2.

### 3.3 Joint radial velocity and photometric analysis of the binary

The goal of the previous section is to obtain precise absolute, dynamical masses. Here we independently analyse the data, using

<table><caption>**Table 2.** Orbital parameters for the BEBOP-3 binary and general fit parameters. Parameters are taken from the KIMA analysis of RVs only unless otherwise specified. The reference time (BJD) used in the analysis is 2459766.975506.</caption><thead><tr><th>Parameter</th><th>Units</th><th>Value</th><th>Note</th></tr></thead><tbody><tr><td>P</td><td>(d)</td><td>13.2176657 ± 0.0000027</td><td></td></tr><tr><td>K</td><td>(km s<sup>-1</sup>)</td><td>19.36413<sup>+0.00091</sup><sub>-0.00102</sub></td><td></td></tr><tr><td>K<sub>2</sub></td><td>(km s<sup>-1</sup>)</td><td>80.22 ± 0.74</td><td>(a)</td></tr><tr><td>e</td><td>–</td><td>0.063255 ± 0.000054</td><td></td></tr><tr><td>ω</td><td>(rad)</td><td>4.90226<sup>+0.00097</sup><sub>-0.00090</sub></td><td>(d)</td></tr><tr><td>ω̇</td><td>(arcsec yr<sup>-1</sup>)</td><td>160 ± 100</td><td></td></tr><tr><td>λ<sub>0</sub></td><td>(rad)</td><td>3.66663 ± 0.00012</td><td>(d)</td></tr><tr><td>T<sub>peri</sub></td><td>[BJD-2450000]</td><td>9756.1098<sup>+0.0021</sup><sub>-0.0019</sub></td><td></td></tr><tr><td>M<sub>pri</sub></td><td>(M<sub>☉</sub>)</td><td>1.083 ± 0.026</td><td>(b)</td></tr><tr><td>M<sub>sec</sub></td><td>(M<sub>☉</sub>)</td><td>0.2615 ± 0.0039</td><td>(b)</td></tr><tr><td>i</td><td>(deg)</td><td>88.7488 ± 0.0092</td><td>(c)</td></tr><tr><td>a<sub>pri</sub></td><td>(AU)</td><td>0.0234853<sup>+0.0000011</sup><sub>-0.0000013</sub></td><td>(b)</td></tr><tr><td>a<sub>sec</sub></td><td>(AU)</td><td>0.09728 ± 0.00090</td><td>(b)</td></tr><tr><td>Jitter</td><td>(m s<sup>-1</sup>)</td><td>6.73<sup>+1.08</sup><sub>-0.93</sub></td><td>–</td></tr><tr><td>ν</td><td>(Student's t shape)</td><td>7.9<sup>+29.0</sup><sub>-3.6</sub></td><td>–</td></tr><tr><td>v<sub>sys</sub></td><td>(km s<sup>-1</sup>)</td><td>24.46146<sup>+0.00072</sup><sub>-0.00075</sub></td><td>–</td></tr><tr><td>v<sub>rest,2</sub></td><td>(km s<sup>-1</sup>)</td><td>25.77 ± 0.77</td><td>(a)</td></tr></tbody></table>

Notes. (a) Parameter from HRCCS.

(b) Parameter from combining KIMA results with HRCCS results and using inclination from photometric analysis.

(c) Parameter taken from joint analysis of RV and photometry.

(d) The argument of pericentre $\omega$ and the true longitude at the reference time $\lambda_0$ are those of the orbit of the primary star around the centre-of mass of the primary-secondary two-body orbit, to convert to the parameters for the secondary, $\pi$ should be subtracted.

the more standard method when dealing with single-lined eclipsing binaries (e.g. Triaud et al. 2013; Davis et al. 2024). Our goal in this section is first to obtain radii for the two stars and an inclination measurement for the binary, and secondly to be in a position to compare derived (model-dependent) stellar masses, to dynamical (model-independent) stellar masses.

Following Davis et al. (2024), we perform a simultaneous joint analysis of both the TESS (Transiting Exoplanet Survey Satellite; Ricker et al. 2015) photometric data [accessed using the Mikulski Archive for Space Telescopes (MAST) web service<fnref n="7" /> and downloaded using the LIGHTKURVE software package (Lightkurve Collaboration 2018)] and the SOPHIE spectroscopic data using the ALLES-FITTER package (Günther & Daylan 2019, 2021). ALLESFITTER uses the DYNESTY nested sampling package (Speagle 2020) and ELLC (Maxted 2016) for binary system models. All currently available sectors (up to and including Sector 74) of TESS data are included from the TESS Science Processing Operations Center (SPOC) pipeline (with 120s cadence) (Jenkins et al. 2016).

The methods outlined in the remainder of this section are based on those presented in Davis et al. (in preparation). The photometric data are chunked into individual data sets to be given to ALLESFITTER as this allows for individual hybrid spline baselines to account for any out of eclipse variability. This is more computationally efficient as well as removing unnecessary complexity that comes from employing the use of a Gaussian process. Upon visual inspection any eclipse events that were partial, showed larger fluctuations, or showed strong

<footnote n="6">⁶We use the line mask of the HARPS pipeline, which has been publicly released on https://www.eso.org</footnote>

<footnote n="7">⁷https://mast.stsci.edu</footnote>
trends in the residuals were removed<fnref n="8" /> (and the remaining data re-analysed). In this analysis we use 15 of the 19 primary eclipse events. The removed events include three partial eclipses and one eclipse which showed clear additional correlated photometric variability compared to the other primary eclipses. Similarly, we use 14 of the 19 secondary eclipses where of the remove events, three were partial, and two show additional correlated photometric variability which the baseline was unable to account for. We choose to remove the eclipse events with residual correlated variability upon visual inspection to take a cautious approach to the data modelling, mitigating for any systematic bias which could be introduced by including noisier data. All radial velocity outliers identified in Section 3.4 are removed for this analysis.

Before completing the global modelling, the spectroscopic data is first analysed in isolation. This allows for appropriate prior values to be determined for radial velocity parameters related to both the binary and planetary companions. The simultaneous global fit is then completed for the photometric binary companion, and the binary and planetary companions in the spectroscopic time-series. With the addition of photometric data in the analysis, additional parameters of the system can be determined, including the radius of the stellar companion and the inclination of the binary orbit. Furthermore, by following the parameter derivation method in Davis et al. (2024), the mass of the primary star ($M_{pri}$) can be determined alongside the mass and radius of the secondary star ($M_{sec}$, $R_{sec}$). These parameters are shown in Table 1. The phase-folded light curves are shown in Fig. A1.

## 3.4 Search for planetary signals in the radial velocities

The radial velocity time-series analysis is performed using the 'BINARIESmodel' of KIMA (Faria et al. 2018; Baycroft et al. 2023a). KIMA utilizes diffusive nested sampling (DNEST4; Brewer & Foreman-Mackey 2018) to explore the parameter space, and allows for the number of Keplerian signals to be a free parameter. The binary orbit is fit with a Keplerian, using separate priors to the planetary signals. The 'BINARIESmodel' is tailored to fitting binary orbits by extending the binary's Keplerian fit with a linear apsidal precession term and by including the radial velocity contribution from general relativity (Baycroft et al. 2023a). The parameters used for the Keplerians are *P*, *K*, *e*, *ω*, $M_0$, where *P* is the orbital period, *K* the radial velocity semi-amplitude, *e* the orbital eccentricity, and *ω* the argument of periastron. Here $M_0$ is the mean anomaly at the reference time; for the purpose of reporting the parameters this is converted into a time of pericentre passage $T_{peri}$. Priors used for the radial velocity analysis with KIMA are shown in Table A2.

Outliers are dealt with in a two-stage process. The likelihood calculation uses a Student's *t* distribution rather than a Gaussian distribution, the wide tails of which allow outliers to avoid skewing the fit while remaining accounted for. The Student's *t* distribution takes a shape parameter *v*, a small value of which leads to a heavy tail, downweighting outliers, and a large value results in a near-Gaussian likelihood. Some extreme outliers exist, likely due to incorrect pointings of the telescope. We identify and remove extreme outliers by performing a first analysis with KIMA, identifying any

individual data points where the ratio of Student's *t* and Gaussian likelihoods $L_T/L_N \ge 10^{10}$. These points are removed and the final analysis is then performed, still using a Student's *t* likelihood to account for any remaining non-extreme outliers. The radial velocity data is made available in Table A1 with the points that we identify as outliers flagged: extreme outliers are flagged with an X; non-extreme outliers are flagged with a Y. Non-extreme outliers are those which have a likelihood ratio $L_T/L_N \ge 10$ in the final analysis.

# 4 RESULTS

## 4.1 Absolute physical parameters for the binary

The radial velocity analysis is performed simultaneously for the binary and the potential planet signals searched for. The detection of a planet and criteria related to that are discussed in Section 4.2. The orbital parameters of the binary are presented in Table 2. We use the values of $K_2$ from Section 3.2 and $i_{bin}$ from Section 3.3 to derive the masses of the two stars.

Table 1 shows parameters for the stars obtained from three different analyses, which each give an independent mass measurement. There is a slight tension between the different mass measurements. The measurements from HRCCS (Section 3.2) and joint radial velocity/photometry fit (Section 3.3) have a $1.7\sigma$ agreement. This difference is not so great, and a thorough study of potential systematic differences between the methods is beyond the scope of this work. We chose to use the mass from HRCCS as it is the least model-dependent therefore we take $M_{pri} = 1.083 \pm 0.026 M_\odot$ and $M_{sec} = 0.2615 \pm 0.0039 M_\odot$.

The dynamical mass of each star is therefore obtained to a precision of 2.4 per cent and 1.5 per cent. The inclination of $88.7488 \pm 0.0092^\circ$ is derived from the joint photometric and radial velocity analysis. We report orbital and fit parameters for the stars/binary in Table 2, this contains parameters from the radial velocity fit (a few of the parameters depend on results from the HRCCS or photometric fit and these are noted).

## 4.2 Detection of the planet BEBOP-3b

The trans-dimensional nature of KIMA makes comparing models with different numbers of signals possible from a single analysis. From a single KIMA analysis we can calculate the Bayes Factor for an (n)-planet model compared to a (n-1)-planet model, $BF_{n,n-1}$. A Bayes Factor over 150 is considered 'very strong evidence' for one model over another (Trotta 2008). We use a maximum of 3 Keplerians, and find that $BF_{1,0} = 3.4 \times 10^{13}$, $BF_{2,1} = 6.0$, and $BF_{3,2} = 1.1$. This means there is overwhelming evidence for a 1-planet model over a 0-planet model, and weak-to-moderate evidence for a 2-planet model over a 1-planet model. There is no evidence supporting 3 or a higher number of planets, justifying the use of 3 as the maximum, so that computation time is not wasted. For further discussion of the use of Bayesian model comparison with KIMA as a detection metric for planet signals see e.g. Standing et al. (2022, 2023); Triaud et al. (2022); Baycroft, Triaud & Kervella (2023b).

While the model comparison favours a 1-Keplerian model, this does not necessarily mean either that the signal is planetary in nature, since some other signal could mimic a Keplerian (e.g. Baycroft et al. 2023b), or that the parameters of the planet are certain, with aliases and degeneracies being possible. Fig. 2 shows the False Inclusion Probability (FIP) periodogram (Hara et al. 2022). The FIP is calculated from the proportion of weighted samples which contain a planet within each period bin. As seen in Fig. 2 the main

<footnote n="8">⁸Partial primary eclipse events at BJD 2459022.08, 2459035.13 and 2459405.39, and partial secondary eclipse events at BJD 2459412.10, 2459755.74 and 2459768.92 are removed. Primary eclipse at BJD 2460317.43 and secondary eclipses at BJD 2459015.57 and 2459028.78 removed due to additional photometric variability.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** False inclusion probability (FIP) periodogram. The dashed line is at a 1 per cent false inclusion probability.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Phased RV plot of BEBOP-3 b, with the residuals below. The parameters for the planet shown are the maximum-likelihood solution. The binned RV data are shown to guide the eye.</figcaption></figure>

detected periodicity is the one between 500 and 600 d and a second periodicity is seen around 1400 d, just surpassing the 1 per cent threshold recommended in Hara et al. (2022). The FIP has been mathematically shown to be optimal, in the sense that it is the exoplanet detection criterion which maximizes true detections for a certain tolerance to false ones (Hara et al. 2024). The validation that these signals are not due to stellar activity is presented in Section 4.4.

A circumbinary planet with $P_{pl} = 547_{-7.6}^{+6.2}$ d and $M_{pl} = 0.558_{-0.048}^{+0.051} M_{Jup}$ is therefore detected, BEBOP-3b. The model favours a moderate eccentricity of $e = 0.247_{-0.089}^{+0.077}$. The maximum-likelihood model for the planetary orbit is shown in Fig. 3 with the individual data shown as well as the residuals below. The binned radial velocities are also shown to guide the eye. The median and $1\sigma$ uncertainties on the planet parameters are presented in Table 3. The distributions and correlations for the posterior samples with the number of planets $N_p = 1$ are shown in Fig. A2.

<table><caption>**Table 3.** Orbital parameters of BEBOP-3 b, note that these are mean orbital elements rather osculating elements. The reference time (BJD) used in the analysis is 2459766.975506.</caption><thead><tr><th>Parameter</th><th>Units</th><th>Value</th><th>Note</th></tr></thead><tbody><tr><td>P</td><td>(d)</td><td>547.0<sup>+6.2</sup><sub>-7.6</sub></td><td>—</td></tr><tr><td>K</td><td>(m s<sup>-1</sup>)</td><td>11.8 ± 1.1</td><td>—</td></tr><tr><td>e</td><td></td><td>0.247<sup>+0.077</sup><sub>-0.089</sub></td><td>—</td></tr><tr><td>ω</td><td>(rad)</td><td>0.86<sup>+0.40</sup><sub>-0.46</sub></td><td>(b)</td></tr><tr><td>λ<sub>0</sub></td><td>(rad)</td><td>2.85 ± 0.20</td><td>(b)</td></tr><tr><td>T<sub>peri</sub></td><td>[BJD-2450000]</td><td>9633<sup>+31</sup><sub>-36</sub></td><td>—</td></tr><tr><td>M<sup>(a)</sup></td><td>(MJup)</td><td>0.558<sup>+0.051</sup><sub>-0.048</sub></td><td>(a)</td></tr><tr><td>a</td><td>(AU)</td><td>1.444 ± 0.017</td><td>—</td></tr></tbody></table>

Notes. (a) True mass of the planet under the assumption it is exactly coplanar with the binary.

(b) the argument of pericentre $\omega$ and the true longitude at the reference time $\lambda_0$ are those of the orbit of the binary around the centre-of-mass of the planet-binary two-body orbit, to convert to the parameters for the planet, $\pi$ should be subtracted.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Orbital configuration of BEBOP-3 showing the orbits of each star and the planet for 600 posterior draws. The stability limit is shown as a dashed line, this was calculated given the binary orbital parameters from the maximum-likelihood solution.</figcaption></figure>

With the use of KIMA we can also report the posterior distributions of the parameters for BEBOP-3b marginalized over the candidate second signal and potential unknown signals. To do this we include the posterior samples from all of the models with different numbers of planets ($N_p = 1, N_p = 2$, and $N_p = 3$). While we prefer to report the parameters from the preferred $N_p = 1$ model, we also report these marginalized distributions in Table A3. These are consistent with the ones reported in Table 3, albeit with slightly lower median values of $K$, $e$, and $M$.

The orbital configuration of the system is shown in Fig. 4. The binary and planet orbits are shown as 600 posterior draws overlain. The region of potential transitability, discussed in Section 5.3, is highlighted in magenta for each planetary orbit shown. The direction to earth is off the right hand side of the diagram. Using the
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Detection limit for BEBOP-3. In blue, the posterior density of the KIMA run forced to fit a signal, with planet b removed. The posterior distribution from the KIMA radial velocity analysis is shown in green for planet-b ($N_p = 1$ model) and that for the candidate outer signal ($N_p = 2$ model) is shown in red, straddling the detection limit. Dashed lines show expected signals of solar-system mass planets.</figcaption></figure>

primary mass $M_{\text{pri}} = 1.083 \, \text{M}_\odot$ and the orbital parameters from the maximum-likelihood posterior sample, we calculate the stability limit using the equation (3) from Holman & Wiegert (1999). This limit is shown on the diagram as a dashed grey line. Defining the scaled semimajor axis as $a_{\text{sc}} = a_{\text{pl}}/r_{\text{stab}}$, the semimajor axis of the planet's orbit divided by the radius of the stability limit (Holman & Wiegert 1999), the value for BEBOP-3b is $a_{\text{sc}} = 4.845_{-0.081}^{+0.073}$. Of the known circumbinary planetary systems, all except for Kepler-1647b have their innermost known planet between $1 < a_{\text{sc}} < 1.8$. This showcases how BEBOP-3b is distinct from the pile-up planets (see Section 5.2).

## 4.3 Further signals

As mentioned in Section 4.2, there is weak-to-moderate evidence for a second planet. The signal corresponds to a planet with $P \approx 1400 \, \text{d}$ and $M \sin i \approx 0.2 \, M_{\text{Jup}}$. The posterior distributions and covariances for planet 1 and planet 2 within the $N_p = 2$ posteriors are shown in Figs A3 and A4.

We compute a detection limit to constrain potential other planets in this system. These are computed as described in Standing et al. (2022, 2023). The maximum-likelihood signal of the detected planet is removed, and then the residual data is analysed with the number of planets that KIMA fits for being fixed to one. The binary signal is not removed, and is then fit in this second analysis in the same way as before. This analysis then returns a posterior distribution of orbits which are consistent with the data but were not detected. The resulting posterior density is shown in Fig. 5 along with the 99 per cent confidence detection limit. Sensitivity is reached down to Neptune-mass planets at short orbital periods and to Saturn-mass planets at long orbital periods.

In the detection limit in Fig. 5, the posterior distributions for planet-b and this outer candidate signal are also displayed. The planet-b distribution is that from the 1-planet posteriors from the KIMA analysis in Section 4.2 (shown in Fig. A2). The candidate's distribution is that from the 2-planet model (shown in Fig. A4). The distribution for the outer signal is majoritarily below the detection limit because the signal is not large enough (or the data not precise/accurate enough) to be detected/refuted.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Periodograms of the activity indicators H-alpha and NaD are presented with the periodogram of the radial velocities (with the binary solution removed) plotted in the background. Dotted vertical lines denote the binary period and planetary period. Horizontal dashed lines represent FAP levels at 10 per cent, 1 per cent, and 0.1 per cent.</figcaption></figure>

Further observations are ongoing to confirm/refute this candidate signal. In Section 4.4 the highest activity signal is found at $\sim 1000 \, \text{d}$. Obtaining more data will also increase the significance and the precision of the activity peaks, and help to ascertain if the candidate signal is from activity or not.

## 4.4 Stellar activity analysis

To investigate the periodic variations in the chromospheric activity indicators, we performed a generalized Lomb–Scargle periodogram analysis of the H-alpha and NaD indices (similar to Standing et al. 2023). These indices exhibit the highest signal-to-noise ratio in our data set, making them the most reliable tracers of periodic signals. Other indicators, such as CaI, HeI, and CaII, displayed significant noise and did not reveal clear periodicities. The analysis was conducted on all the observations spanning a time baseline of 2197 d except those deemed outliers. The same outliers were removed as in the KIMA radial velocity analysis with which the planet was identified, these are marked with an X in Table A1. The periodograms are shown in Fig. 6 along with false alarm probability (FAP) levels. The FAP is calculated by determining the power levels corresponding to a given probability threshold, which quantifies the likelihood that a detected signal is due to random noise rather than a true periodic signal. In this analysis, the FAP levels are set at 10 per cent, 1 per cent, and 0.1 per cent, and their corresponding power thresholds are computed from the periodogram to assess the significance of the detected peaks.

For the H-alpha index, the periodogram identifies a peak at a frequency of $0.0179 \pm 0.0001 \, \text{d}^{-1}$, corresponding to a period of $55.57 \pm 0.35 \, \text{d}$. The maximum power of the detected signal is 0.05. However, when compared with the FAP levels, this power remains below the 10 per cent significance threshold, indicating that the detected periodicity is not statistically significant. The weighted mean of the H-alpha indices is 0.1868, with an RMS scatter of 0.006.

A similar analysis of the NaD index reveals a dominant period of $84.07 \pm 0.72 \, \text{d}$. The power of this signal is $\sim 0.1$, which also remains below the 10 per cent FAP significance threshold, suggesting that the observed periodicity is not robust. The weighted mean of the NaD indices is 0.4365, with an RMS scatter of 0.0195. For comparison, the periodogram of the radial velocities, with the binary signal removed, is shown in the background of Fig. 6.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Stability analysis of an inner planet in the BEBOP-3 system assuming $m_{\text{inner}} \sim 10 \, \text{M}_{\oplus}$ and coplanar orbits. For fixed initial conditions (Tables 2 and 3), the parameter space of the system is explored by varying the orbital period $P_{\text{inner}}$ and the eccentricity $e_{\text{inner}}$ of planet-c. The step size is 1.15 d in orbital period and 0.0025 in eccentricity. For each initial condition, the system is integrated over 5000 yr and a stability indicator is calculated which involved a frequency analysis of the mean longitude of the inner planet (see the text). Chaotic diffusion is indicated when the main frequency of the mean longitude varies in time. A stability index of -2 correspond to highly unstable orbits, while -8 correspond to orbits which are likely to be stable on a billion-years time-scales. The dashed lines indicate mean-motion resonances with the planet BEBOP-3b. The circles represent the locations of the innermost circumbinary planet in the other known systems, scaled to the parameters of the BEBOP-3 binary.</figcaption></figure>

Although both activity indices exhibit periodic variations at similar time-scales, the lack of statistical significance implies that these detections will require further observations and improved data sampling to confirm the nature of the signal. Neither H-alpha periodicity nor NaD periodicity are close to the strong signal seen around 550 d in the radial velocities. As such it is unlikely stellar photospheric processes are causing the observed radial-velocity modulation. We therefore interpret it as planetary.

# 5 DISCUSSION

## 5.1 Stability analysis

The BEBOP-3(AB)b system is hierarchical ($a_{\text{bin}}/a_{\text{b}} \ll 1$), this condition typically results in stability as it avoids mean-motion resonance overlap between the planet and binary, many numerical studies find these orbits to be stable (e.g. Holman & Wiegert 1999; Doolin & Blundell 2011). Therefore if there are no planets interior to its orbit the stability of BEBOP-3b should be assured (e.g. Correia, Boué & Laskar 2016). Motivated by the population of known circumbinary planets located near the stability limit, we investigate the stability of the system when including orbits interior to that of planet b. We simulate the BEBOP-3 system using two close-in stars in a nearly circular orbit ($a_{\text{bin}} = 0.12 \, \text{au}$, $e_{\text{bin}} = 0.06$, Table 2), together with an outer planet further away in an eccentric orbit ($a_{\text{b}} = 1.44 \, \text{au}$, $e_{\text{b}} = 0.28$;<fnref n="9" /> Table 3).

In order to get a reliable and comprehensive view of the stability of potential inner planets in the system, we performed a global

frequency analysis (Laskar 1990, 1993) in the same way as achieved for other circumbinary planetary systems (eg. Correia et al. 2005; Standing et al. 2023). The system is integrated on a regular 2D mesh of initial conditions in the vicinity of the best fit (Tables 2 and 3). We used the symplectic integrator SABAC4 (Laskar & Robutel 2001), with a step size of 0.001 yr and general relativity corrections. Each initial condition is integrated for 5000 yr, and a stability indicator, $\Delta = |1 - n_b/n_{b'}|$, is computed. Here, $n_b$ and $n_{b'}$ are the main frequency of the mean longitude of the planet over 2500 yr and 5000 yr, respectively, calculated via the frequency analysis (Laskar 1993). The results are reported in colour, where yellow represents strongly chaotic trajectories with $\Delta > 10^{-2}$, while extremely stable systems for which $\Delta < 10^{-8}$ are shown in purple/black. Orange indicates the transition between the two, with $\Delta \sim 10^{-4}$.

We test the possibility of an inner planet, between the binary stars and planet-b, with a radial-velocity signature of $1 \, \text{m s}^{-1}$. Such a planet would correspond to a mass $m_{\text{inner}} \sim 10 \, \text{M}_{\oplus}$. We assume the best-fitting solution for the already known bodies in the system (Tables 2 and 3), coplanar orbits, and vary the orbital period and the eccentricity of a putative inner planet (Fig. 7). We observed that close to the binary ($50 \lesssim P_{\text{inner}} \lesssim 200 \, \text{d}$), the system is mostly stable with small eccentricity ($e_{\text{inner}} \lesssim 0.1 - 0.2$), except close to some mean-motion resonances between the putative inner planet and the binary. Interestingly, as we move further away ($P_{\text{inner}} \gtrsim 200 \, \text{d}$), we observed the exact opposite behaviour: the system becomes globally unstable, but stability is still possible, even with moderate eccentricities, near some mean-motion resonances between the inner planet and planet-b.

In Fig. 7, the location of the known innermost planets of all the known circumbinary systems (excluding Kepler-1647b) are shown. The orbital period has been scaled by the ratio of the orbital period of the respective binaries. All of these points lie in the region of stability, and 6 of these 12 planets have a mass measurement (or non-measurement) consistent with being below the mass of Neptune, which is the approximate sensitivity limit in this region for BEBOP-

<footnote n="9">⁹The maximum-likelihood parameters from the KIMA analysis are used for the stability analysis which while within the uncertainties presented in Tables 2 and 3, do not always correspond exactly to the mean values of the posterior reported.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Eccentricity distribution of the population of currently known circumbinary exoplanets, as well as BEBOP-3b. The x-axis is the semimajor axis of the planets scaled by the Holman–Wiegert stability radius for each binary. The solid curve is an indicative stability envelope above which orbits are expected to be unstable.</figcaption></figure>

3. This suggests that an additional planet consistent with the current population would remain undetected in the data collected thus far.

## 5.2 An eccentric long-period circumbinary planet

Of the circumbinary planetary systems discovered by the transit method or radial velocities, only Kepler-1647 (Kostov et al. 2016) does not contain a planet orbiting close to the instability limit (Martin & Triaud 2014). BEBOP-3 joins Kepler-1647 as the second circumbinary planetary system in this category. Additionally, some circumbinary planets have been found further from the stability limit, detected with other methods such as microlensing (Bennett et al. 2016) and eclipse timing variations (e.g. Qian et al. 2011).

In the cases of Kepler-1647 and BEBOP-3 it is unclear whether there are additional, currently undetected, planets orbiting closer to the binary or whether these circumbinary systems followed a different formation and evolutionary pathway. The stability analysis in Section 5.1 shows that there is a significant stable region exactly where it would be expected for a planet comparable to the Kepler planets to exist. In addition, the inner planets of circumbinary systems regularly appear to have a low mass [e.g. Kepler-47b (Orosz et al. 2012), TOI-1338/BEBOP-1b (Standing et al. 2023)]. In the case of TOI-1338/BEBOP-1b, despite 184 spectra, the planet remains undetected in the RVs alone (Standing et al. 2023). Based on these considerations it could be expected that either a planet below the radial velocity detection threshold is present there, or that a planet used to be present there but has since become destabilized.

BEBOP-3b has a marginal detection of eccentricity. It is the only circumbinary planet with a non-zero measured eccentricity from the radial velocity method so far. Fig. 8 displays the location of BEBOP-3b compared to the other circumbinary planets in eccentricity and semimajor axis scaled by the radius of the stability limit (calculated from equation 3 in Holman & Wiegert 1999). The red curve shows an approximate stability envelope above which a planet's orbit would cross into the instability zone. The known circumbinary planets are mostly clustered in the bottom left corner of this plot, near the stability limit and at low eccentricities. The planets with higher eccentricities have a slightly larger semimajor axis to compensate for this, leading to a pattern of planets tracing just below the stability envelope. As

is visible in Fig. 8, BEBOP-3b is distinct from the other known circumbinary planets in being both eccentric and far from crossing into the instability region.

Despite appearing somewhat special among circumbinary planets, there are a few other considerations/caveats:

(i) Many gas giants orbiting single stars occupy eccentric orbits (e.g. Winn & Fabrycky 2015), and compared to them BEBOP-3b is not abnormal.

(ii) This is only the second RV-discovered circumbinary planet, and the second longest orbital period; only more detections can tell whether its current parameters are unusual or not.

(iii) Eccentricity measurements of planetary signals from radial velocities can sometimes be unreliable. This is made worse when there are undetected planets in the system (Hara et al. 2019), in particular when near orbital period commensurabilities (Anglada-Escudé, López-Morales & Chambers 2010). Comparing the eccentricity distributions in Figs A2 and A3 show that including the outer planet candidate slightly decreases the significance of the eccentricity measurement. More data is required to assess both the significance of the second planet and of the eccentricity measurement.

(iv) The prior used for the eccentricity is that from Kipping (2013) which is tailored to planets around single-stars and may not be suitable for circumbinary planets.

More circumbinary planet detections and better characterization of the existing planets are needed to determine whether there is a population of lone long-period circumbinary planets, and also to better constrain the eccentricity distribution.

Coleman et al. (2024) investigated the formation history of circumbinary planets in the TOI-1338/BEBOP-1 system, which has similar stellar masses and binary orbital period to BEBOP-3 (see Table A4). The parameters of BEBOP-3b presented here appear consistent with the population of planets formed in their model.

## 5.3 Transitability prospects

Circumbinary planets, if misaligned with the inner binary, will undergo nodal precession of their orbit. This means that most ($\approx$ 95 per cent) of circumbinary planets orbiting eclipsing binaries are expected to eventually come into a transitable configuration (Martin & Triaud 2015). Transitability means that the planet crosses in front of the path of the binary orbit and might therefore transit in front of one or both of the stars, however in a given conjunction the planet still might happen to not transit. More thorough discussions can be found in Martin & Triaud (2015) and Martin (2017).

From our radial velocity solution for BEBOP-3b, we cannot assess the mutual inclination between the orbits of the planet and binary, or the position of the planet in its nodal precession cycle. We therefore consider times of 'possible transitability', when the planet is near inferior conjunction. If the planet is the correct orientation to be able to transit, it would have to do so during these times.

We take 1000 posterior samples for the binary and planet parameters and calculate the times of possible transitability for each sample. We then plot a histogram of these times in Fig. 9, giving a posterior probability distribution of times of possible transitability for each conjunction of the planet. These are compared with the times that the system has been observed with TESS (TIC 289949453) (Ricker et al. 2015). Based on this, even if the planet did transit at any of these conjunctions, it is highly unlikely to be present in any TESS data. We do note that the orbital configuration of the planet, as shown in Fig. 4 is optimal for a planetary transit as the transit would occur very near pericentre where the distance between planet and star is smallest.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Distributions of potential transitability of BEBOP-3 b from the radial velocity fit. The TESS sectors in which BEBOP-3 was observed are shown, along with the range of time over which the radial velocities were taken.</figcaption></figure>

## 5.4 The BEBOP catalogue of circumbinary planets with radial velocity measurements

BEBOP-3b is now the fourth published circumbinary planet detected using radial velocities. We introduce here the BEBOP catalogue of circumbinary planets with radial velocity measurements and the rationale for the numbering. The catalogue is shown in Table A4. This catalogue is intended to collect all of the circumbinary planetary systems with parameters measured using radial velocities in one place. The catalogue does include circumbinary planets that were detected first in transit and might include planets detected by other teams using radial velocities. For each system we give a BEBOP catalogue entry number without presuming to replace the names.

Kepler-16b and TIC 172900988b were first detected in transit (Doyle et al. 2011; Kostov et al. 2021) and subsequently detected in radial velocities (Triaud et al. 2022; Sairam et al. 2024). Based on the timing of their detections relative to BEBOP-1c (Standing et al. 2023) we number them 0 and 2 in the BEBOP catalogue leaving the new planet presented here as number 3. Note that TOI-1338b which was detected in transit (Kostov et al. 2020) has not yet been detected in radial velocities<fnref n="10" /> so the inner planet is only given an upper mass-limit in the BEBOP catalogue.

## 6 CONCLUSION

In this work we present the radial velocity detection of a circumbinary planet BEBOP-3b. This planet holds a distinct position relative to the known population of circumbinary planets orbiting main-sequence binaries. BEBOP-3b orbits with a long period relative to the binary and is on a moderately eccentric orbit. A dynamical analysis reveals that there is ample space for another undetected planet near the stability limit, so BEBOP-3b may not be so outlying relative to the known population. There is an additional longer period candidate signal with insufficient evidence to claim as a detection. We extracted stellar activity indices, and none show variability on the radial-velocity modulation we observed. We consequently interpret this modulation as being produced by Doppler reflex motion of a circumbinary planet. The absolute dynamical mass of the stars are calculated using high resolution cross-correlation spectroscopy. The planet's minimum mass as derived from the radial velocities is likely close to its true mass, assuming coplanarity. Should the inclination of its orbit become known, BEBOP-3b would also have an absolute mass.

## ACKNOWLEDGEMENTS

These scientific results would not have been possible without a generous allocation of telescope time by France Programme National de Planétologie, and by the dedicated staff working at the Observatoire de Haute-Provence, particularly during the COVID pandemic. We also would like to thank the help of other observers part of the timeshare agreement between groups using the SOPHIE instruments. This research leading to these results is supported by a grant from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement n° 803193/BEBOP) and by the ERC/UKRI Frontier Research Guarantee programme (EP/Z000327/1/CandY). This work was supported by the 'Programme National de Planétologie' (PNP) of CNRS/INSU co-funded by CNES and the 'Action Thématique de Physique Stellaire' (ATPS) of CNRS/INSU PN Astro co-funded by CEA and CNES. ACMC acknowledges support from FCT – Fundação para a Ciência e a Tecnologia, I.P., Portugal, through the projects 2024.01252.BD, CFisUC (UIDB/04564/2020 and UIDP/04564/2020, with DOI identifiers 10.54499/UIDB/04564/2020 and 10.54499/UIDP/04564/2020, respectively). AVF acknowledges the support of the Institute of Physics (IOP) through the Bell Burnell Graduate Scholarship Fund. RPN acknowledges funding from the STFC grant ST/X000931/1. Computations for this research were performed using the University of Birmingham's BlueBEAR HPC service http://www.birmingham.ac.uk/bear. The stability maps were performed at the Oblivion Supercomputer at the University of Évora (https://oblivion.hpc.uevora.pt).

## DATA AVAILABILITY

The TESS data is available to download via the MAST portal (https://mast.stsci.edu). Spectroscopic data are available on the SOPHIE archive (http://atlas.obs-hp.fr/sophie/), and were obtained under Prog.ID 18B.DISC.TRIA and 19A.PNP.SANT. Reduced radial velocities are found in Table A1, a full version of which can be found online.

## REFERENCES

* Anglada-Escudé G., López-Morales M., Chambers J. E., 2010, ApJ, 709, 168

* Baranne A. et al., 1996, A&AS, 119, 373

* Baycroft T. A., Triaud A. H. M. J., Faria J. A., Correia A. C. M., Standing M. R., 2023a, MNRAS, 521, 1871

* Baycroft T. A., Triaud A. H. M. J., Kervella P., 2023b, MNRAS, 526, 2241

* Benedict G. F., Harrison T. E., 2017, AJ, 153, 258

* Bennett D. P. et al., 2016, AJ, 152, 125

* Blanco-Cuaresma S., 2019, MNRAS, 486, 2075

<footnote n="10">¹⁰There is a photodynamical mass measurement presented in Wang & Liu (2024) but it is unclear to what extent the measurement is influenced by the radial velocities.</footnote>
Bouchy F., Isambert J., Lovis C., Boisse I., Figueira P., Hébrard G., Pepe F., 2009a, EAS Publ. Ser., 37, 247

Bouchy F. et al., 2009b, A&A, 505, 853

Brewer B. J., Foreman-Mackey D., 2018, J. Stat. Softw., 86, 1

Coleman G. A. L., Nelson R. P., Triaud A. H. M. J., 2023, MNRAS, 522, 4352

Coleman G. A. L., Nelson R. P., Triaud A. H. M. J., Standing M. R., 2024, MNRAS, 527, 414

Collins K. A. et al., 2018, AJ, 156, 234

Correia A. C. M., Udry S., Mayor M., Laskar J., Naef D., Pepe F., Queloz D., Santos N. C., 2005, A&A, 440, 751

Correia A. C. M., Boué G., Laskar J., 2016, Celest. Mech. Dyn. Astron., 126, 189

Courcol B. et al., 2015, A&A, 581, A38

Davis Y. T. et al., 2024, MNRAS, 530, 2565

Doolin S., Blundell K. M., 2011, MNRAS, 418, 2656

Doyle L. R. et al., 2011, Science, 333, 1602

Faria J. P., Santos N. C., Figueira P., Brewer B. J., 2018, J. Open Source Softw., 3, 487

Freckelton A. V. et al., 2024, MNRAS, 531, 4085

Goldberg M., Fabrycky D., Martin D. V., Albrecht S., Deeg H. J., Nowak G., 2023, MNRAS, 525, 4628

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

Günther M. N., Daylan T., 2019, record ascl:1903.003

Günther M. N., Daylan T., 2021, ApJS, 254, 13

Hara N. C., Boué G., Laskar J., Delisle J.-B., Unger N., 2019, MNRAS, 489, 738

Hara N. C., Unger N., Delisle J.-B., Díaz R. F., Ségransan D., 2022, A&A, 663, A14

Hara N. C., de Poyferré T., Delisle J.-B., Hoffmann M., 2024, Ann. Appl. Stat., 18, 749

Heidari N., 2022, PhD Thesis, Universit'e C^note d'Azur; Shahid Beheshti University (Tehran)

Heidari N. et al., 2024, A&A, 681, A55

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

Jenkins J. M. et al., 2016, in Chiozzi G., Guzman J. C., eds, Proc. SPIE Conf. Ser. Vol. 9913, Software and Cyberinfrastructure for Astronomy IV. SPIE, Bellingham, p. 99133E

Kalman D., 1996, College Math. J., 27, 2

Kipping D. M., 2013, MNRAS, 434, L51

Kostov V. B. et al., 2016, ApJ, 827, 86

Kostov V. B. et al., 2020, AJ, 159, 253

Kostov V. B. et al., 2021, AJ, 162, 234

Kurucz R. L., 2005, Mem. Soc. Astron. Ital. Suppl., 8, 14

Laskar J., 1990, Icarus, 88, 266

Laskar J., 1993, Phys. D Nonlin. Phenomena, 67, 257

Laskar J., Robutel P., 2001, Celest. Mech. Dyn. Astron., 80, 39

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

Lightkurve Collaboration 2018, recordascl:1812.013

Martin D. V., 2017, MNRAS, 465, 3235

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

Martin D. V., Triaud A. H. M. J., 2015, MNRAS, 449, 781

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

Martin D. V. et al., 2019, A&A, 624, A68

Maxted P. F. L., 2016, A&A, 591, A111

Maxted P. F. L., Triaud A. H. M. J., Martin D. V., 2023, Universe, 9, 498

Meschiari S., 2012, ApJ, 761, L7

Moriwaki K., Nakagawa Y., 2004, ApJ, 609, 1065

Muñoz D. J., Lai D., 2015, Proc. Natl Acad. Sci., 112, 9264

Orosz J. A. et al., 2012, Science, 337, 1511

Orosz J. A. et al., 2019, AJ, 157, 174

Paardekooper S.-J., Leinhardt Z. M., Thébault P., Baruteau C., 2012, ApJ, 754, L16

Penzlin A. B. T., Kley W., Nelson R. P., 2021, A&A, 645, A68

Pepe F., Mayor M., Galland F., Naef D., Queloz D., Santos N. C., Udry S., Burnet M., 2002, A&A, 388, 632

Pepper J. et al., 2007, PASP, 119, 923

Perruchot S. et al., 2008, in McLean I. S., Casali M. M., eds, Proc. SPIE Conf. Ser. Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II. SPIE, Bellingham, p. 70140J

Pierens A., Nelson R. P., 2008, A&A, 483, 633

Pierens A., Nelson R. P., McNally C. P., 2021, MNRAS, 508, 4806

Pollacco D. L. et al., 2006, PASP, 118, 1407

Pollacco D. et al., 2008, MNRAS, 385, 1576

Qian S. B. et al., 2011, MNRAS, 414, L16

Ricker G. R. et al., 2015, J. Astron. Telesc. Instrum. Syst., 1, 014003

Roy O., Vetterli M., 2007, The effective rank: a measure of effective dimensionality, in 15th European Signal Processing Conference. IEEE, Poznan, Poland, p. 606

Sairam L. et al., 2024, MNRAS, 527, 2261

Sbordone L., Bonifacio P., Castelli F., Kurucz R. L., 2004, Mem. Soc. Astron. Ital. Suppl., 5, 93

Sebastian D., Triaud A. H. M. J., Brogi M., 2024a, MNRAS, 527, 10921

Sebastian D. et al., 2024b, MNRAS, 530, 2572

Sebastian D. et al., 2025, MNRAS, 540, 2914

Snellen I. A. G., de Kok R. J., de Mooij E. J. W., Albrecht S., 2010, Nature, 465, 1049

Speagle J. S., 2020, MNRAS, 493, 3132

Standing M. R. et al., 2022, MNRAS, 511, 3571

Standing M. R. et al., 2023, Nat. Astron., 7, 702

Triaud A. H. M. J., 2011, Ph.D. thesis, https://ui.adsabs.harvard.edu/abs/2011PhDT......22T

Triaud A. H. M. J. et al., 2013, A&A, 549, A18

Triaud A. H. M. J. et al., 2017, A&A, 608, A129

Triaud A. H. M. J. et al., 2022, MNRAS, 511, 3561

Trotta R., 2008, Contemp. Phys., 49, 71

Wang M.-T., Liu H.-G., 2024, AJ, 168, 31

Welsh W. F. et al., 2012, Nature, 481, 475

Winn J. N., Fabrycky D. C., 2015, ARA&A, 53, 409

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## APPENDIX A: APPENDIX
<figure><img src="image_11.png" /><figcaption>**Figure A1.** Left panels: Phase folded primary eclipses for the binary and the residuals for the best-fitting model. Right panels: Phase folded secondary eclipses for the binary and the residuals for the best-fitting model. In both plots, 50 random posterior samples are displayed.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure A2.** Corner plot showing posterior distributions and correlations for the 1-planet model.</figcaption></figure>

<caption>**Figure A3.** Corner plot showing posterior distributions and correlations for planet 1 in the 2-planet model.</caption>
<figure><img src="image_13.png" /><figcaption>**Figure A4.** Corner plot showing posterior distributions and correlations for planet 2 in the 2-planet model.</figcaption></figure>

<table><caption>**Table A1.** Radial velocity observations taken with SOPHIE between 2018 November 13 and 2024 November 18. Radial velocity outliers are indicated with one of two flags X: those that are vast outliers found in the preliminary run and removed, these are likely wrong pointings of the telescope, and Y: those that are not excluded from the fit but accounted for with the Student's t and have a ratio of $L_T/L_N \ge 10$ in the second analysis. Spectra excluded from the HRCCS analysis are given the flag A (wrong pointing), B (low SNR), or C (blended phases).</caption><thead><tr><th>Flag</th><th>Time<br>[BJD-2400 000]</th><th>V<sub>rad</sub><br>(km s<sup>-1</sup>)</th><th>σV<sub>rad</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="12">Y</td><td>58436.63359</td><td>43.7735</td><td>0.0039</td></tr><tr><td>58437.68801</td><td>42.9949</td><td>0.0018</td></tr><tr><td>58536.46316</td><td>5.4022</td><td>0.0024</td></tr><tr><td>58539.46867</td><td>21.9657</td><td>0.0021</td></tr><tr><td>58541.46505</td><td>40.2170</td><td>0.0018</td></tr><tr><td>58556.34359</td><td>43.7120</td><td>0.0021</td></tr><tr><td>58560.44450</td><td>16.9170</td><td>0.0029</td></tr><tr><td>58564.34291</td><td>8.8614</td><td>0.0018</td></tr><tr><td>58794.64166</td><td>42.7380</td><td>0.0026</td></tr><tr><td>58804.69818</td><td>30.9787</td><td>0.0024</td></tr><tr><td>58828.65843</td><td>8.6339</td><td>0.0025</td></tr><tr><td>58855.59869</td><td>11.8996</td><td>0.0028</td></tr><tr><td>58894.46828</td><td>7.3007</td><td>0.0020</td></tr><tr><td>58899.35178</td><td>43.9156</td><td>0.0026</td></tr><tr><td rowspan="3">C</td><td>58903.33834</td><td>22.9174</td><td>0.0020</td></tr><tr><td>58907.45466</td><td>6.4592</td><td>0.0027</td></tr><tr><td>58912.32987</td><td>43.4670</td><td>0.0129</td></tr><tr><td rowspan="3">B</td><td>58917.38073</td><td>16.4726</td><td>0.0031</td></tr><tr><td>58922.32180</td><td>17.3861</td><td>0.0028</td></tr><tr><td>59130.68250</td><td>-3.6578</td><td>0.0361</td></tr><tr><td>X,A</td><td>59155.65324</td><td>13.9372</td><td>0.0018</td></tr><tr><td>...</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<table><caption>**Table A2.** Prior distributions used in the KIMA analysis. $N$ is a Normal/Gaussian distribution, $U$ is a Uniform distribution, $LU$ is a LogUniform distribution, $MLU$ is a ModifiedLogUniform distribution (LogUniform between the two values and Uniform between 0 and the lower value), $K$ is a Kumaraswamy distribution (a computationally efficient approximation to a $\beta$ distribution).</caption><thead><tr><th>Parameter</th><th>Units</th><th>Prior distribution</th></tr></thead><tbody><tr><td>P<sub>bin</sub></td><td>(d)</td><td>N(13.2177, 0.001)</td></tr><tr><td>K<sub>bin</sub></td><td>(km s<sup>-1</sup>)</td><td>N(19.366, 0.05)</td></tr><tr><td>e<sub>bin</sub></td><td>–</td><td>N(0.063, 0.01)</td></tr><tr><td>ω<sub>bin</sub></td><td>(rad)</td><td>N(4.9, 0.1)</td></tr><tr><td>M<sub>0,bin</sub></td><td>(rad)</td><td>U(0, 2π)</td></tr><tr><td>ω̂<sub>bin</sub></td><td>(arcsec yr<sup>-1</sup>)</td><td>N(0, 1000)</td></tr><tr><td>P<sub>pl</sub></td><td>(d)</td><td>LU(50, 5490)</td></tr><tr><td>K<sub>pl</sub></td><td>(m s<sup>-1</sup>)</td><td>MLU(0.1, 1000)</td></tr><tr><td>e<sub>pl</sub></td><td>–</td><td>K(0.867, 3.03)</td></tr><tr><td>ω̂<sub>pl</sub></td><td>(rad)</td><td>U(0, 2π)</td></tr><tr><td>M<sub>0,pl</sub></td><td>(rad)</td><td>U(0, 2π)</td></tr><tr><td>V<sub>sys</sub></td><td>(km s<sup>-1</sup>)</td><td>U(23.968, 24.968)</td></tr><tr><td>Jitter</td><td>(m s<sup>-1</sup>)</td><td>MLU(0.1, 100)</td></tr><tr><td>ν</td><td>–</td><td>LU(2, 1000)</td></tr></tbody></table>

<table><caption>**Table A3.** Orbital parameters of BEBOP-3 b, marginalized over the candidate and undetected potential planets by including the posterior samples from $N_P = 1$, $N_P = 2$, and $N_P = 3$. Note that these are mean orbital elements rather osculating elements. The reference time (BJD) used in the analysis is 2459766.975506.</caption><thead><tr><th>Parameter</th><th>Units</th><th>Value</th><th>Note</th></tr></thead><tbody><tr><td>P</td><td>(d)</td><td>546.0<sup>+8.7</sup><sub>-9.2</sub></td><td>–</td></tr><tr><td>K</td><td>(m s<sup>-1</sup>)</td><td>9.8<sup>+1.8</sup><sub>-2.0</sub></td><td>–</td></tr><tr><td>e</td><td>–</td><td>0.20<sup>+0.11</sup><sub>-0.13</sub></td><td>–</td></tr><tr><td>ω</td><td>(rad)</td><td>0.90<sup>+0.59</sup><sub>-0.74</sub></td><td>(b)</td></tr><tr><td>λ<sub>0</sub></td><td>(rad)</td><td>2.79 ± 0.26</td><td>(b)</td></tr><tr><td>T<sub>peri</sub></td><td>[BJD-2450000]</td><td>59633<sup>+45</sup><sub>-62</sub></td><td>–</td></tr><tr><td>M</td><td>(MJup)</td><td>0.469<sup>+0.083</sup><sub>-0.090</sub></td><td>(a)</td></tr><tr><td>a</td><td>(AU)</td><td>1.443 ± 0.019</td><td>–</td></tr></tbody></table>

Notes. (a) True mass of the planet under the assumption it is exactly coplanar with the binary.

(b) The argument of pericentre $\omega$ and the true longitude at the reference time $\lambda_0$ are those of the orbit of the binary around the centre-of-mass of the planet- binary two body orbit, to convert to the parameters for the planet, $\pi$ should be subtracted.
<table><caption>**Table A4.** The BEBOP catalogue of circumbinary planets with radial velocity measurements. $1\sigma$ uncertainties on parameters are provided at the level of the last two significant digits. $^x$ Transit detection from Doyle et al. (2011). $^y$ Transit detection and parameters from Kostov et al. (2020). $^z$ Transit detection from Kostov et al. (2021).</caption><thead><tr><th>Catalogue entry</th><th>Other names</th><th>P<sub>bin</sub><br>(d)</th><th>M<sub>pri</sub><br>(M<sub>☉</sub>)</th><th>M<sub>sec</sub><br>(M<sub>☉</sub>)</th><th>e<sub>bin</sub></th><th>P<sub>pl</sub><br>(d)</th><th>M<sub>pl</sub><br>M<sub>sup</sub></th><th>e<sub>pl</sub></th><th>Source of RVs</th><th>TIC ID</th></tr></thead><tbody><tr><td>BEBOP-0 b</td><td>Kepler-16 b<sup>x</sup></td><td>41.077772(51)</td><td>0.654(17)</td><td>0.1964(31)</td><td>0.15994(10)</td><td>226.0(1.7)</td><td>0.313(39)</td><td>≤ 0.21</td><td>Triaud et al. (2022)</td><td>299096355</td></tr><tr><td>BEBOP-1 b</td><td>TOI-1338 b<sup>y</sup></td><td>14.6085579(57)</td><td>1.098(17)</td><td>0.307(3)</td><td>0.155522(29)</td><td>95.174(35)<sup>y</sup></td><td>≤ 0.0685(29)</td><td>0.0880(43)<sup>y</sup></td><td>Standing et al. (2023)</td><td>260128333</td></tr><tr><td>BEBOP-1 c</td><td>TOI-1338 c</td><td>14.6085579(57)</td><td>1.098(17)</td><td>0.307(3)</td><td>0.155522(29)</td><td>215.5(3.3)</td><td>0.205(37)</td><td>≤ 0.16</td><td>Standing et al. (2023)</td><td>260128333</td></tr><tr><td>BEBOP-2 b</td><td>TIC 172900988 b<sup>z</sup></td><td>19.657878(34)</td><td>1.23681(39)</td><td>1.20207(33)</td><td>0.448234(90)</td><td>151.2(1.8)</td><td>1.90(25)</td><td>≤ 0.11</td><td>Sairam et al. (2024)</td><td>172900988</td></tr><tr><td>BEBOP-3 b</td><td>—</td><td>13.2176657(27)</td><td>1.083(26)</td><td>0.3615(39)</td><td>0.063255(54)</td><td>547.0<sup>(+6.2)</sup><br>(-7.6)</td><td>0.558<sup>(+5.1)</sup><br>(-4.8)</td><td>0.247<sup>(+77)</sup><br>(-89)</td><td>This work</td><td>28994953</td></tr></tbody></table>

This paper has been typeset from a TeX/LATEX file prepared by the author.