# BEBOP VIII. SOPHIE radial velocities reveal an eccentric, circumbinary brown dwarf

Amaury Triaud, Thomas Baycroft, Neda Heidari, Alexandre Santerne, Aleyna Adamson, Isabelle Boisse, Gavin Coleman, Alexandre Correia, Yasmin Davis, Magali Deleuil, et al.

► To cite this version:

Amaury Triaud, Thomas Baycroft, Neda Heidari, Alexandre Santerne, Aleyna Adamson, et al. BEBOP VIII. SOPHIE radial velocities reveal an eccentric, circumbinary brown dwarf. Monthly Notices of the Royal Astronomical Society, 2025, 544 (2), pp.2180-2192. 10.1093/mnras/staf1884. hal-05490811

HAL Id: hal-05490811

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

Submitted on 2 Mar 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 VIII. SOPHIE radial velocities reveal an eccentric, circumbinary brown dwarf

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

Aleyna Adamson,<sup>1</sup> Isabelle Boisse,<sup>3</sup> Gavin A.L. Coleman<sup>1,4</sup> Alexandre C.M. Correia<sup>5,6</sup> Yasmin

T. Davis<sup>1</sup>,<sup>3</sup> Magali Deleuil,<sup>3</sup> Guillaume Hébrard,<sup>2</sup> David V. Martin<sup>1,7</sup> Pierre F.L. Maxted<sup>1,8</sup> Richard

P. Nelson,<sup>4</sup> Lalitha Sairam<sup>1,9</sup> Daniel Sebastian,<sup>1,10</sup> Matthew R. Standing<sup>1,11</sup> Adam T. Stevenson<sup>1,12</sup>

Amalie Stokholm<sup>1</sup>,<sup>1</sup> Mathilde Timmermans<sup>1</sup> and Stéphane Udry<sup>12</sup>

Accepted 2025 October 25. Received 2025 October 20; in original form 2025 September 25

## ABSTRACT

Circumbinary configurations offer a test of planet formation in an altered environment, where the inner binary has perturbed a protoplanetary disc. Comparisons of the physical and orbital parameters between the circumbinary planet population and the population of exoplanets orbiting single stars will reveal how these disc perturbations affect the assembly of planets. Circumbinary exoplanets detected thus far typically have masses < 3 M<sub>jup</sub> raising the question of whether high-mass circumbinary planets are possible, and also whether population features such as the brown dwarf desert would appear in circumbinary configurations like for single star systems. Here, we report observations taken with the SOPHIE high-resolution spectrograph. These observations reveal an m<sub>b</sub> sin i<sub>b</sub> = 20.9 M<sub>jup</sub> outer companion, on an eccentric (e = 0.43), 1800 d orbit, which we call BEBOP-4 (AB) b. Using dynamical arguments we constrain the true mass m<sub>b</sub> < 26.3 M<sub>jup</sub>. The inner binary's two eclipsing stellar components have masses M<sub>A</sub> = 1.51 M<sub>☉</sub>, and M<sub>B</sub> = 0.46 M<sub>☉</sub>. Their orbital period is 72 d, and their eccentricity is 0.27. This system contains the longest period binary surveyed by the BEBOP project. BEBOP-4 b is expected to be detectable using Gaia DR4 single epoch astrometric measurements. Despite a large period ratio of ~ 25 : 1, the substantial eccentricities of both orbits mean that the outer orbit is on the edge of orbital stability, and located in between two destabilizing secular resonances. Should the outer companion survive, the BEBOP-4 system appears like a precursor to several post-common envelope binaries exhibiting eclipse timing variations where very massive circumbinary companions have been proposed.

**Key words:** techniques: radial velocities – exoplanets – binaries: eclipsing – brown dwarfs.

## 1 INTRODUCTION

Only few circumbinary exoplanets are known. 14 have been discovered thanks to transits (e.g. L. R. Doyle et al. 2011; J. A. Orosz et al. 2012; V. B. Kostov et al. 2020), whilst another three have been identified from radial-velocity observations (M. R. Standing et al. 2023; T. A. Baycroft et al. 2025a, b). In addition, nine have

been proposed using direct imaging (e.g. V. Bailey et al. 2014; T. J. Dupuy et al. 2018), five from microlensing (e.g. D. P. Bennett et al. 2016; C. Han et al. 2024), one has been confirmed from eclipse timing variations (M. Goldberg et al. 2023), with another 38 being proposed (e.g. K. Beuermann et al. 2010; T. Borkovits et al. 2015; Ö. Baştırk et al. 2023; T. A. Baycroft, A. H. M. J. Triaud & P. Kervella 2023b), and one has been detected from pulsar timing variations (S. E. Thorsett, Z. Arzoumanian & J. H. Taylor 1993; S. Sigurdsson et al. 2003). Compared to over 6000 exoplanets known to orbit single stars, circumbinary planet would appear rare at first sight. However, about

<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>Institut d'astrophysique de Paris, UMR 7095 CNRS université Pierre et Marie Curie, 98 bis, boulevard Arago, F-75014 Paris, France</footnote>

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

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

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

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

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

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

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

<footnote n="10"><sup>10</sup>Thüringer Landessternwarte Tautenburg, Sternwarte 5, D-07778 Tautenburg, Germany</footnote>

<footnote n="11"><sup>11</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 n="12"><sup>12</sup>Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

<footnote>*E-mail: a.traiaud@bham.ac.uk</footnote>
10–12 per cent of close binaries host circumbinary gas giants (D. J. Armstrong et al. 2014; D. V. Martin & A. H. M. J. Triaud 2014; T. A. Baycroft, A. H. M. J. Triaud & D. Lai 2024), a number which is similar to single stars (A. Cumming et al. 2008; M. Mayor et al. 2011; A. Santerne et al. 2016; B. J. Fulton et al. 2021).

Another rare configuration is a single Sun-like star with a brown dwarf companion (e.g. D. Grether & C. H. Lineweaver 2006; B. Ma & J. Ge 2014; J. Sahlmann et al. 2014; A. H. M. J. Triaud et al. 2017; T. W. Carmichael, D. W. Latham & A. M. Vanderburg 2019; S. Dalal et al. 2021). Dubbed the *brown dwarf desert*, observations reveal that, particularly on orbital periods < 100–1000 d (F. Kiefer et al. 2019), very few brown dwarfs exists, with the rarest being of a mass in the range 25–45 M<sub>jup</sub> (J. Sahlmann et al. 2011; A. T. Stevenson et al. 2023). This pattern has been interpreted for companions nearer the Hydrogen-burning limit (~80 M<sub>jup</sub>; S. S. Kumar 1963; I. Baraffe et al. 2015) as the tail-end of stellar formation, while near the Deuterium-burning limit (~13 M<sub>jup</sub>; I. Baraffe et al. 2003) those objects might represent the high-end of planet formation, whether from core-accretion followed by runaway gas accretion (J. B. Pollack et al. 1996; G. A. L. Coleman, R. P. Nelson & A. H. M. J. Triaud 2023), or via direct collapse (D. Stamatellos & A. P. Whitworth 2009). As a result, a 25 M<sub>jup</sub> has been adopted to separate between planetary-mass and stellar-mass brown dwarfs (J. Sahlmann et al. 2011; J. Schneider et al. 2011).

Whether circumbinary brown dwarfs exist remains a puzzle. In two systems detected with direct imaging, Delorme-1 (P. Delorme et al. 2013) and VHS J1256-1257 (T. J. Dupuy et al. 2023), the companion mass is of the same order of magnitude as the components of the binary (possibly also in 2M1510; P. Calissendorff et al. 2019; T. A. Baycroft et al. 2025b), hence giving the impression they are closer to a triple system than a true circumbinary configuration. Furthermore, in those systems, as a well as in all other directly detected circumbinary brown dwarfs, the brown dwarf is found at a large separation from the inner binary, far from the short-period brown dwarf desert seen in single Sun-like stars. The most evidenced short period circumbinary brown dwarf is in HD 202206, however its nature is debated. This system contains a Sun-like star with two companions, one with $m_b \sin i_b = 17$ M<sub>jup</sub>, and $m_c \sin i_c = 2.4$ M<sub>jup</sub> (A. C. M. Correia et al. 2005). Whilst J. Couetdic et al. (2010) use orbital stability arguments to estimate the true masses are $m_b < 34$ M<sub>jup</sub> and $m_c < 4.4$ M<sub>jup</sub>, and thus $i_b \sim i_c > 30^\circ$, G. F. Benedict & T. E. Harrison (2017) argue instead that astrometric measurements made with *Hubble*'s fine guidance sensor imply $i_b \sim i_c \sim 8^\circ - 10^\circ$, which would mean $m_b = 0.09$ M<sub>☉</sub> and $m_c = 18$ M<sub>jup</sub>, a circumbinary brown dwarf. Eventually Gaia DR4 astrometry should resolve this conflict, but the orbital stability arguments are hard to shake and more likely than not, J. Couetdic et al. (2010) are correct and HD 202206 is likely not hosting a circumbinary brown dwarf.

Meanwhile, BEBOP (Binaries Escorted By Orbiting Planets) is a project to survey nearby close binaries with radial velocities to detect circumbinary exoplanets. The project started in 2014, with a large sample of low-mass eclipsing binaries, and a limited radial-velocity precision. In the first publication of the project, we could place an upper limit on the occurrence rate of circumbinary brown dwarfs $\eta < 6.5$ per cent (D. V. Martin et al. 2019). In more recent work, we also notice a paucity of companions with mass > 3 M<sub>jup</sub>, found to be 5× rarer than in equivalent single star systems (T. A. Baycroft et al. 2024). In fact, within the confirmed transit, radial-velocity and eclipse-timing population, only one circumbinary planet has a mass > 3 M<sub>jup</sub>, Kepler-1660 b, with 4.9 M<sub>jup</sub> (M. Goldberg et al. 2023).

In this paper, we now report the clear detection of a circumbinary brown dwarf in a system we call BEBOP-4. This discovery is impor-

tant to better understand planet formation. Already planetary-mass brown dwarfs are difficult to assemble in single star systems and they are very rarely obtained in simulations of planetary system formation (e.g. G. A. L. Coleman & R. P. Nelson 2016; A. Emsenhuber et al. 2021). However, circumbinary planet formation is expected to be more difficult yet, at least in the close vicinity of the central binary, since the circumbinary protoplanetary disc is anticipated to be highly disturbed and more turbulent than for a single star. Vigorous turbulence in protoplanetary discs typically reduces the efficiency of pebble accretion (e.g. M. M. Rosenthal et al. 2018; A. Pierens, R. P. Nelson & C. P. McNally 2021), and perturbations to planetesimal orbits also slow the rate of planetesimal accretion (S.-J. Paardekooper et al. 2012).

In Section 2 we describe our data collection, and how radial velocities are extracted from the spectra. In Section 3, we discuss how we measure the masses for both components of the stellar binary. These data are modelled in Section 4. We interpret and discuss our results in Section 5 before concluding.

## 2 RADIAL VELOCITY DATA COLLECTION

Target selection for the Southern component of BEBOP is described in D. V. Martin et al. (2019). For the Northern survey, a brief description is provided in T. A. Baycroft et al. (2025a). The system, HD 111605 ($\alpha = 12:50:18.53, \delta = +26:09:34.29, V_{mag} = 9.6$), with spectral type G3, was identified by the KELT survey (Kilodegree Extremely Little Telescope; J. Pepper et al. 2007), and reported as a likely eclipsing binary (K. A. Collins et al. 2018, where it is known as KC08C10921), a false positive transiting planet signal. Its properties made it compatible with an EBLM (Eclipsing Binary – Low Mass; A. H. M. J. Triaud et al. 2013; P. F. L. Maxted, A. H. M. J. Triaud & D. V. Martin 2023), like all systems observed during the first phase of the BEBOP radial velocity search for circumbinary exoplanets (D. V. Martin et al. 2019). As the BEBOP-North survey catalogue was assembled, this system received the designation EBLM J1250+26 under which name it was observed. The Northern survey started with a short reconnaissance phase (2018 September 1 to 2019 February 28) using the SOPHIE spectrograph (S. Perruchot et al. 2008) at the Observatoire de Haute Provence (OHP), to assess the radial-velocity precision of each primary star, verify the orbital and physical parameters of the secondaries were in line with expectations, and to remove systems with high stellar activity (e.g. from a lack of variations in the bisector span as in D. Queloz et al. 2001). EBLM J1250+26 was selected as a top priority system for the BEBOP-North survey. Intense radial-velocity monitoring with SOPHIE, for BEBOP-North, started in earnest in 2019 March 11 and are still ongoing.

By construction, the BEBOP survey runs a blind survey, meaning all systems are in principle observed at the same cadence and under similar conditions. The collected data is analysed every six months, at the end of each observing semester to perform a quality control (remove systems with high stellar activity, triple systems etc.), and search for candidates. When candidates are identified, they are flagged, but not observed particularly more intensively so as to avoid introducing an observer's bias, however reporting candidates remains important as a means to convince telescope time allocation committees about the validity of continuing the survey for competitively obtained telescope time.

In the case of EBLM J1250+26, around March 2020 we had noticed a downward, quadratic drift in the residuals of the orbital motion of the binary. At that time we suspected a triple system, but still within planetary range. The COVID-19 pandemic started and
<figure><img src="image_2.png" /><figcaption>**Figure 1.** HRCCS detection of BEBOP-4 B. Left: K-focused cross-correlation map. Centre: SALTIRE model for the SOPHIE data, with parameters from MCMC. Right: residual CCF map after subtraction of the model. Data and model have been median-subtracted to allow for a common scale to be shared between the three panels. A secondary maximum is identified in the CCF map with a similar significance at $K_B \sim 62 \text{ km s}^{-1}$, $v_{\text{sys}} \sim -3 \text{ km s}^{-1}$. This false signal, at the incorrect systemic velocity, is actually less significant and only gets highlighted due to the presence of the main signal. In the residual map, the 'false' maximum becomes almost indistinguishable from other noise structures.</figcaption></figure>

observatories all over the world closed. The OHP stopped operations entirely between mid-March and mid-May 2020, reopening quite early, but with limited availability. As a consequence, priority was reduced on EBLM J1250+26 since drifts can be followed efficiently at a slower cadence. Four measurements were obtained, in 2020 May 29, 2020 December 2, 2021 February 25, and 2021 March 3, the period during which travel to OHP was difficult. When data were analysed, it became immediately apparent the last three of these measurements did not follow the expected drift, and that a companion was likely undergoing a periastron passage. An attempt to observe was made during twilight in early August 2020, but by then the target was too low on the horizon. Observations resumed at high cadence and with high priority in 2021 November 25. Regular fits to the orbit indicated a full orbital phase would be completed around 2024, which indeed happened. A new periastron passage is expected in September 2025 when the star will unfortunately be below the horizon at night.

Because of these unique circumstance, compared to other systems of the BEBOP survey, EBLM J1250+26 (BEBOP-4) is the system that has been the least blindly monitored.

As of 2025 March 1, we collected 84 high resolution spectra, obtained between 2019 January 7 and 2025 February 27, with the high resolution échelle spectrograph SOPHIE, mounted on the T193 at the Observatoire de Haute-Provence (S. Perruchot et al. 2008). Data were acquired under Prog.ID 18B.DISC.TRIA (PI Triaud) and 19A.PNP.SANT (PI Santerne). All data have been acquired in ObjAB mode where fibre A is placed on the target and fibre B is positioned on the sky to enable the removal of solar spectral lines reflected by the Moon. The observations were acquired using the high resolution (HR) mode with resolution $R = 75\,000$. Standard calibrations were taken at the start of the night and at approximately 2 h intervals throughout the night to track the instrumental zero-point.

The SOPHIE radial velocities (RVs) are computed by the SOPHIE Data Reduction System (DRS; F. Bouchy et al. 2009b), which cross-correlates the spectra with a weighted numerical mask corresponding to a G2-type star and subsequently fits a Gaussian to the cross-correlation function (CCF) (A. Baranne et al. 1996; F. Pepe et al. 2002). To enhance the precision of the SOPHIE RV measurements, we applied the recently optimized procedures described in N. Heidari (2022) and N. Heidari et al. (2024). These include: (1) corrections for charge transfer inefficiency in the CCD (F. Bouchy et al. 2009a); (2) correction of template (colour) effects (N. Heidari et al. 2024); (3)

removal of moonlight contamination by utilizing the simultaneous sky spectrum recorded via the second SOPHIE fibre (D. 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 (B. Courcol et al. 2015; N. Heidari 2022; N. 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 84 spectra obtained have a median SNR = 46.8 (at 5550 Å) and a median radial-velocity precision of 3.4 m s⁻¹.

## 3 DETERMINING STELLAR MASSES

We could assume a mass of $M_A = 1.584 \pm 0.086 \text{ M}_\odot$ for the primary star, as was obtained from a homogeneous spectral analysis of the BEBOP sample presented in A. V. Freckelton et al. (2024). However, recently developed methods allow some of our binaries to be transformed from a single to a double-lined binary (D. Sebastian et al. 2024b, 2025), to measure absolute, dynamical masses for both components.

We use the same high resolution cross-correlation spectroscopy (HRCCS) method as in D. Sebastian et al. (2024b, 2025) and T. A. Baycroft et al. (2025a), using Principal Component Analysis (PCA) and cross-correlation in the secondary's rest frame (with the HARPS M2 mask) to extract the weak lines of the secondary star. For this system, we make an additional modification whereby any additional RV shifts caused by the circumbinary brown dwarf companion (Section 4) are first removed from the spectra before we attempt to detect the secondary star.

The CCF maps (Fig. 1) are fit with the SALTIRE model (D. Sebastian, A. H. M. J. Triaud & M. Brogi 2024a) to extract $K_B$ and $v_{\text{sys}}$. However, the SALTIRE uncertainties are on the fit parameters, and do not estimate systematic uncertainty from the presence of correlated noise. To quantify the systematic error, we use a bootstrapping analysis, where we randomly draw 50 per cent of the data and measure the best-fitting parameters (as in T. A. Baycroft et al. 2025a). This is repeated 200 times, and we impose that no more than 15 (out of 42) spectra can be common between any permutation (adapted slightly from the recommendations of D. Sebastian et al. 2025, to allow 200 random permutations to be generated). From the
200 fit results, we take the 1-σ percentiles around the median as the systematic uncertainty in the parameters. The Markov chain Monte Carlo (MCMC) fit error is ±0.207 km s⁻<fnref n="1" />, and the systematic error is ±2.02 km s⁻<fnref n="1" />, and these are added quadrature.

We detect a signal with 5.8-σ significance, and measure a secondary semi-amplitude of $K_B = 51.02 \pm 2.03$ km s⁻<fnref n="1" />, implying $M_A = 1.507 \pm 0.152$ and $M_B = 0.459 \pm 0.028$ (assuming $i_{\text{bin}} = U[88.5, 91.5]^{\circ}$, and following the IAU recommended equations; A. Prša et al. 2016).

Unfortunately no primary eclipses were observed by *Transiting Exoplanet Survey Satellite (TESS)*; G. R. Ricker et al. (2015) so we cannot perform an analysis like in Y. T. Davis et al. (2024), to measure the stellar radii accurately. This analysis will be performed once a primary eclipse has been collected by *TESS*, although no new observations of this target are scheduled through *TESS* Cycle 8.

We note that $M_A$ is measured differently by different methods (see Fig. A1). Thanks to our new measurement, we confirm the fiducial $M_A$ masses used by the BEBOP survey as produced in A. V. Freckelton et al. (2024) are fully consistent to within 1-σ for this system. Using HRCCS, we likely measure $M_A$ more accurately, though at a slightly reduced precision to that achieved through spectral analysis and isochrone fitting, owing to the relatively low detection significance of the secondary.

## 4 ANALYSIS AND RESULTS

The radial velocities were analysed with *kima* (J. P. Faria et al. 2018), using the ‘BINARIES model’ (T. A. Baycroft et al. 2023a). The usual orbital elements are used as parameters. In addition, we fit a *jitter* term to rescale the uncertainties should they need to. The model also has an option to fit for apsidal precession, by having one additional free parameter, $\dot{\omega}$, and also includes a radial-velocity correction due the effects from general relativity (T. A. Baycroft et al. 2023a, 2025b). We chose to let $\dot{\omega}$ be a free parameter. A non-zero value would provide an additional constraint on the companion’s parameters. This model, and the fitting procedure have been validated and used in several publications (M. R. Standing et al. 2023; T. A. Baycroft et al. 2023a, 2025a, b; L. Sairam et al. 2024a, b).

In our resulting fit, we find a model with a circumbinary companion is clearly preferred with a Bayes Factor BF ~ 10⁶⁴ over a model with no circumbinary companion. The best-fitting result for the companion and the binary, results in residuals with RMS scatter of 10.3 m s⁻<fnref n="1" />, this is larger than the mean uncertainty of the data set and requires a jitter of ≈ 10 m s⁻<fnref n="1" />. The companion has an orbital period $P_b = 1823.5_{-4.9}^{+5.1}$ d, and a semi-amplitude $K_b = 243.7 \pm 6.5$ m s⁻<fnref n="1" />. All fit parameters are reported in Table 1. The radial velocity curve for the companion is shown in Fig. 2. The overall highest likelihood model is plotted in Fig. A2 and the binary Keplerian in Fig. A3. The posterior on the orbital parameters of the outer companion are graphically represented in Fig. A4. The orbital configuration of the system is shown in Fig. 3, where 600 random draws from the posterior are plotted. The region shown in magenta is where the companion might pass in front of the binary and as such a transit could potentially happen if the configuration is right.

We use $M_A$ and $M_B$ from Section 3, and thus the minimum mass of the circumbinary companion is $m_b \sin i_b = 20.9 \pm 1.3$ M$_{\text{jup}}$, above the Deuterium-burning limit (S. S. Kumar 1963; I. Baraffe et al. 2003), but below the 25 M$_{\text{jup}}$ that have been taken as the separation between planetary-mass and stellar-mass brown dwarfs (J. Sahlmann et al. 2011; J. Schneider et al. 2011). Constraints on the mutual inclination between the outer brown dwarf’s orbit and the inner

<table><caption>**Table 1.** Orbital parameters for the BEBOP-4 binary, planet, and general fit parameters. Stellar masses, semimajor axes, and $K_B$ are calculated from SB2 analysis in Section 3 as well as the *kima* analysis, and all remaining parameters are taken or derived from the *kima* analysis of RVs. Note that these are mean orbital elements rather than osculating elements. Since $\omega$ evolves with time, the reference time at which the value reported is calculated is 2459 736.088894 (BJD). This reference time is also that at which $M_0$ and $\lambda_0$ are calculated. The spectral type and a variety of alternate names for the BEBOP-4 system are listed as well.</caption><thead><tr><th colspan="3" style="text-align:center;">BEBOP-4-G3 V</th></tr><tr><th colspan="3" style="text-align:center;">α = 12<sup>h</sup>50<sup>m</sup>18.52<sup>s</sup>; δ = +26°09′34.0″</th></tr><tr><th rowspan="2">Parameter</th><th>Units</th><th>Value</th></tr><tr><th colspan="2" style="text-align:center;">Orbital fit parameters</th></tr></thead><tbody><tr><td>P<sub>b</sub></td><td>(d)</td><td>71.96857 ± 0.00015</td></tr><tr><td>K<sub>A</sub></td><td>(m s<sup>-1</sup>)</td><td>15 547.3 ± 2.1</td></tr><tr><td>K<sub>B</sub></td><td>(m s<sup>-1</sup>)</td><td>51 018 ± 2031</td></tr><tr><td>e<sub>bin</sub></td><td></td><td>0.26887 ± 0.00011</td></tr><tr><td>ω<sub>bin</sub><sup>(a)</sup></td><td>(rad)</td><td>4.81629<sup>+0.00049</sup><sub>-0.00051</sub></td></tr><tr><td>ω&#x0307;<sub>bin</sub></td><td>(arcsec yr<sup>-1</sup>)</td><td>88<sup>+49</sup><sub>-46</sub></td></tr><tr><td>M<sub>0,bin</sub></td><td>(rad)</td><td>1.00142<sup>+0.00048</sup><sub>-0.00043</sub></td></tr><tr><td>P<sub>b</sub></td><td>(d)</td><td>1823.5<sup>+5.1</sup><sub>-4.9</sub></td></tr><tr><td>K<sub>b</sub></td><td>(m s<sup>-1</sup>)</td><td>243.7 ± 6.5</td></tr><tr><td>e<sub>b</sub></td><td></td><td>0.428 ± 0.016</td></tr><tr><td>ω<sub>b</sub><sup>(a)</sup></td><td>(rad)</td><td>3.134 ± 0.020</td></tr><tr><td>M<sub>0,b</sub></td><td>(rad)</td><td>2.162 ± 0.021</td></tr><tr><th colspan="3" style="text-align:center;">Assumed and derived parameters</th></tr><tr><td>i<sub>bin</sub></td><td>(rad)</td><td>[1.545, π/2]</td></tr><tr><td>M<sub>A</sub></td><td>(M<sub>☉</sub>)</td><td>1.507 ± 0.152</td></tr><tr><td>M<sub>B</sub></td><td>(M<sub>☉</sub>)</td><td>0.459 ± 0.028</td></tr><tr><td>a<sub>λ</sub></td><td>(AU)</td><td>0.099063 ± 0.000013</td></tr><tr><td>a<sub>B</sub></td><td>(AU)</td><td>0.325 ± 0.013</td></tr><tr><td>λ<sub>0,bin</sub><sup>(a)</sup></td><td>(rad)</td><td>6.34386 ± 0.00030</td></tr><tr><td>T<sub>peri,bin</sub></td><td>[BJD-2 400 000]</td><td>59 724.6183<sup>+0.0049</sup><sub>-0.0056</sub></td></tr><tr><td>i<sub>b</sub></td><td>(rad)</td><td>[0.917, π/2]</td></tr><tr><td>m<sub>b</sub> sin i<sub>b</sub></td><td>(M<sub>μup</sub>)</td><td>20.9 ± 1.3</td></tr><tr><td>a<sub>b</sub></td><td>(AU)</td><td>3.634 ± 0.011</td></tr><tr><td>λ<sub>0,b</sub><sup>(a)</sup></td><td>(rad)</td><td>5.820 ± 0.017</td></tr><tr><td>T<sub>peri,b</sub></td><td>[BJD-2 400 000]</td><td>59 108.5 ± 5.7</td></tr><tr><th colspan="3" style="text-align:center;">Other fit parameters</th></tr><tr><td>Jitter</td><td>(m s<sup>-1</sup>)</td><td>10.0 ± 1.1</td></tr><tr><td>ν</td><td>(Student's t shape)</td><td>74<sup>+354</sup><sub>-62</sub></td></tr><tr><td>ν<sub>sys</sub></td><td>(m s<sup>-1</sup>)</td><td>5277.5 ± 1.8</td></tr></tbody></table>

Note. (a) The values of the argument of pericentre $\omega$ and the true longitude at the reference time $\lambda_0$ reported are those of the orbit of the primary star around the centre-of-mass of its respective two-body orbit, this is the $\omega$ used in the radial velocity equation. To obtain the $\omega$ or $\lambda_0$ of the secondary star’s or brown dwarf’s orbit, $\pi$ should be subtracted from the values reported.

eclipsing binary can be obtained with a measurement of the binary’s apsidal precession.

The apsidal precession of the BEBOP-4(AB) binary is constrained to $\dot{\omega}_{\text{bin}} = 88_{-46}^{+49}$ arcsec yr⁻<fnref n="1" />. The contribution from general relativity.<fnref n="1" /> (using equation A17 from T. A. Baycroft et al. 2023a) is $\dot{\omega}_{\text{bin}} = 0.977 \pm 0.059$ arcsec yr⁻<fnref n="1" />, meaning that the measured apsidal rate that would be induced by a tertiary companion is $\dot{\omega}_{\text{bin}} = 87_{-46}^{+49}$ arcsec yr⁻<fnref n="1" />. Fig. A2 depicts the contribution of $\dot{\omega}_{\text{bin}}$ to

<footnote n="1">¹The apsidal precession induced by tidal deformation and rotational flattening is negligible.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Phased radial velocity Keplerian and residuals for the detected circumbinary companion. The parameters used are those from the maximum-likelihood solution.</figcaption></figure>

the fit. We use $\Delta i$ to represent the difference between the binary's and brown dwarf's inclinations, analogous to the mutual inclination under the assumption that the longitudes of ascending nodes of the two orbits are equal. Assuming the brown dwarf companion is coplanar with the binary (i.e. $\Delta i = 0^{\circ} \Rightarrow m_b = m_b \sin i_b = 20.9 \pm 1.3 \, \text{M}_{\text{jup}}$), the predicted apsidal precession (using equation A5 from T. A. Baycroft et al. 2023a) is $\dot{\omega}_{\text{bin}} = 101.2 \pm 5.4 \, \text{arcsec yr}^{-1}$, consistent with the measured value.

We can use the apsidal precession constraint to obtain an upper limit on $\Delta i$ and therefore also on the true mass of the companion. As the $\Delta i$ increases, the induced apsidal precession rate decreases and eventually becomes negative as the companion moves from a coplanar to a polar orbit (Z. Zhang & D. C. Fabrycky 2019; T. A. Baycroft et al. 2025b). The weighted distance between the measured precession rate and the predicted precession rates as a function of inclination is shown in Fig. 4, along with the corresponding masses. At<fnref n="2" /> $3\sigma$, the mutual inclination is $\Delta i < 37.5^{\circ}$, $0.917 \le i_b \le 2.225$ (rad), and $m_b \le 26.3 \, \text{M}_{\text{jup}}$. At 5σ: $\Delta i < 46.7^{\circ}$, $0.756 \le i_b \le 2.385$ (rad), and $m_b \le 30.5 \, \text{M}_{\text{jup}}$<fnref n="3" /> This constrains BEBOP-4b to being in the planetary-mass regime of brown dwarfs, as such it enters the BEBOP catalogue of circumbinary planets shown in Table A2. This constraint also indicates that the outer orbit's mutual inclination is most likely below the classical Lidov–Kozai angle (39°21; M. L. Lidov 1962; Y. Kozai 1962). We also note that the BEBOP-4 configuration is tight

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Face-on view of the orbital configuration of the BEBOP-4 system showing the orbits of each star and the brown dwarf companion for 600 posterior draws. The observer is located on the right hand side. For illustrative purposes, the stability limit is calculated using (M. J. Holman & P. A. Wiegert 1999) and the binary orbital parameters from the maximum-likelihood solution is also shown. In magenta we highlight the parts of the companion orbit where it might transit one of the stars (calculated as the orbital phases where the companion crosses in front of the area swept out by the binary orbit, assuming coplanarity).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Weighted distance between the measured apsidal precession rate (with the GR contribution removed) and the predicted apsidal precession for various mutual inclination values (i.e. inclinations of the circumbinary companion). Dashed lines show the $3\sigma$ and $5\sigma$ constraints.</figcaption></figure>

<footnote n="2">²Here we assume the inclination of the binary $i_{\text{bin}} = \pi/2$.</footnote>

<footnote n="3">³The function converting $\dot{\omega}$ to $\Delta i$ is symmetrical about $\Delta i = 90^{\circ}$. As such retrograde orbits are also possible.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Detection limit plot for the BEBOP-4 system. In blue, the posterior density of the *kima* run forced to fit a signal, with the brown-dwarf companion removed. The posterior distribution for the companion is shown as a separate density map, in green (above the 13 M<sub>jup</sub> line). Dashed lines show expected signals of solar-system mass planets. The 99 per cent confidence detection limit as well as its uncertainty is also shown.</figcaption></figure>

enough that Lidov–Kozai cycles induced by the companion on the inner binary would be unlikely anyway. Lidov–Kozai cycles happen when the ratio of angular momenta $L_{\text{out}}/L_{\text{in}} > 0.5$ (D. V. Martin & A. H. M. J. Triaud 2016) and here $L_{\text{b}}/L_{\text{bin}} = 0.16$.

We calculate detection limits for any further signals in the radial velocity data in the same way as computed in M. R. Standing et al. (2022, 2023) and T. A. Baycroft et al. (2025a). The signal of the circumbinary companion is removed from the data, using the maximum-likelihood posterior sample, and the *kima* analysis is re-run with the number of planets fixed to 1. The resulting posterior distribution corresponds to all the Keplerian signals that are consistent with the noise in the data, and a 99 per cent upper limit in semi-amplitude is calculated. The detection limit is shown in Fig. 5 along with the posterior distribution for the circumbinary companion for comparison. The uncertainty in the detection limit is $E = \frac{1.97}{\sqrt{N}}$ where E is the fractional uncertainty in the detection limit for a given period bin, and N is the number of posterior samples in that bin. This relation is based on numerical experiments with *kima* (M. R. Standing et al. 2023) and is similar to the relation used in M. R. Standing et al. (2023) except they used a different scaling constant. From this exercise we confidently exclude the presence of any other planetary companion with a mass greater than Jupiter's.

Based on the current parameters, we calculate the detectability of BEBOP-4b with *Gaia*. Using the Gaia Observation Forecast tool,<fnref n="4" /> we find that there will be up to 82 epochs of astrometric data published in *Gaia* DR4. Assuming a typical epoch precision of 40 µas (B. Holl et al. 2023), we calculate the *Gaia* SNR as in J. Sahlmann, A. H. M. J. Triaud & D. V. Martin (2015) in which a conservative *Gaia* SNR = 20 is recommended as a detection threshold. The BEBOP-4(AB) binary will have a *Gaia* SNR ≈ 100, and BEBOP-4(AB)b a *Gaia* SNR ≈ 40.<fnref n="5" /> We therefore expect that *Gaia* DR4 will be sensitive to the binary's orbit and the brown dwarf's orbit, allowing for

measurements of the true mass and the mutual inclination. Amongst all confirmed circumbinary planets discovered with the transit and/or radial velocity methods, BEBOP-4 b is currently the only one which is detectable with *Gaia*.

# 5 DISCUSSION

## 5.1 Comparisons with other planets

Compared to single star systems, BEBOP-4 b does not particularly stand out. The exoplanet encyclopaedia (exoplanet.eu; J. Schneider et al. 2011) lists 168 objects with mass > 10 M<sub>jup</sub> discovered either with radial velocities or transits, 24 of which orbit stars with masses > 1.4 M<sub>☉</sub>.

Compared to known circumbinary planets identified with transits or radial velocities, BEBOP-4 b is unique. It is the most massive by a factor > 3. In addition, BEBOP-4 both contains the longest binary period, and the largest planetary period in the known population. Thanks to its large semi-amplitude, it also has the most precise eccentricity, and is the most eccentric so far discovered.

In terms of both mass and orbital parameters, BEBOP-4 b is compatible with some of the proposed solutions for post common envelope binaries, where eclipse timing variations appear to imply a planetary companion. Some examples are: HW Vir which has claimed companions of comparable mass but much longer orbital periods than BEBOP-4b (K. Beuermann et al. 2012; T. A. Baycroft et al. 2023b); V1828 Aql (L. A. Almeida, F. Jablonski & C. V. Rodrigues 2013; M. Wolf et al. 2021) and QS Vir (S. B. Qian et al. 2010; Ö. Baştürk et al. 2023) which have claimed companions of similar orbital periods but slightly lower masses.<fnref n="6" /> As such, BEBOP-4 might be seen as a precursor to the companions claimed orbiting post-common envelope binaries. The evolution of a population of circumbinary systems from main-sequence to circum-double-white-dwarf was performed by G. Columba et al. (2023) who find that 20–30 per cent of planets survive the whole evolution. Coincidentally, BEBOP-4A is itself evolving off the main-sequence (near its terminal age main sequence), and about to ignite H-shell burning, on its way to become a sub-giant.

*Gaia*'s contribution in detecting new circumbinary exoplanets will mainly be if massive circumbinary planets > 5 M<sub>jup</sub> exist (orbiting main sequence binaries, but also in post-common envelope systems; J. Sahlmann et al. 2015). BEBOP-4 b reveals that some detections are to be expected.

Protoplanetary disc masses have been observed to scale with central mass (e.g. I. Pascucci et al. 2016; S. M. Andrews et al. 2018) and the combined stellar mass at the centre of the system is $M_A + M_B = 1.97$ M<sub>☉</sub>, uncommonly large for most planets orbiting single stars, but a fairly regular occurrence for circumbinary planets (e.g. TIC 172900988; V. B. Kostov et al. 2021; L. Sairam et al. 2024a). As such the BEBOP-4 system is a good example of why studying circumbinary planets will support a better understanding on how planet formation scales with central mass. Whilst detecting a planet orbiting a single 2 M<sub>☉</sub> star is difficult (lack of absorption lines, and fast rotation), doing so by sharing the mass between two stars makes for, in this particular case, an easier detection.

<footnote n="4">⁴https://gaia.esac.esa.int/gost/</footnote>

<footnote n="5">⁵It is not expected that all of the epochs predicted will be present in the *Gaia* data; for the SNR of the companion to get below 20 would require that only 24 of the 82 measurements are available.</footnote>

<footnote n="6">⁶We note here that many claims for circumbinary exoplanets orbiting post-common envelope binaries rely on eclipsing timing variations that may have other origins than being produced by a planetary body (e.g. M. Zorotovic & M. R. Schreiber 2013; A. Hardy et al. 2015; T. A. Baycroft et al. 2023b).</footnote>
## 5.2 Possible formation pathways

As mentioned in Section 1, the occurrence frequency of brown dwarf companions to main-sequence stars has a distinct minimum for masses between ~25 and 45 M<sub>jup</sub>, with the occurrence rate increasing both towards the hydrogen burning limit at ~80 M<sub>jup</sub> and towards the deuterium burning limit at ~13 M<sub>jup</sub> (e.g. D. Grether & C. H. Lineweaver 2006; J. Sahlmann et al. 2014; S. Dalal et al. 2021). This bimodality suggests two formation channels are at play, with higher mass brown dwarfs being the low-mass tail of the stellar population that forms through fragmentation of molecular cloud cores, and lower mass objects being the high-mass end of planet formation in protoplanetary discs. Based on its minimum mass, and the constraint on its true mass from the binary's apsidal precession, BEBOP-4b appears to be an example of the latter population that formed in a circumbinary disc.

The discussion below about possible formation scenarios assumes that during the epoch of BEBOP-4b's formation, a circumbinary disc is present around the BEBOP-4 binary system, and the binary has component masses and orbital elements that are similar to those observed today. Consideration of the earlier evolution goes beyond the scope of this discussion. A further point to note is that the central binary contains a total mass ~2 M<sub>☉</sub> and disc masses are known to scale at least linearly with the central stellar mass (e.g. I. Pascucci et al. 2016). Hence, we assume the circumbinary disc contains sufficient mass to form a planet as massive as BEBOP-4b.

Formation in a circumbinary disc could have occurred through fragmentation of the disc (e.g. D. Stamatellos & A. P. Whitworth 2009) or core accretion (J. B. Pollack et al. 1996; G. A. L. Coleman et al. 2023). In the disc fragmentation scenario, the cooling time in the disc needs to be short and on the order of the local orbital time for a gravitationally unstable disc to form one or more bound clumps. Hence, fragmentation is expected in the outer regions of a protoplanetary disc beyond ~50 au (e.g. Z. Zhu et al. 2012). Given that BEBOP-4b has a semimajor axis of $a_p = 3.63$ au, formation through disc fragmentation requires it to have migrated inwards to its currently observed location. This is generally expected, since simulations of clumps formed by disc fragmentation find that disc-driven inwards migration can occur rapidly (A. C. Boley et al. 2010; C. Baruteau, F. Meru & S.-J. Paardekooper 2011; Z. Zhu et al. 2012). In their study, Z. Zhu et al. (2012) showed that clumps that grow massive enough to open gaps in the disc dramatically slow their migration. Hence, a plausible scenario for the origin of BEBOP-4b is formation at large radius in a massive circumbinary disc followed by inwards migration to its current location. However, caution should be taken. Z. Zhu et al. (2012) considered a coarsely sampled set of initial conditions, and the objects that opened gaps and survived migration were more massive than BEBOP-4b. Furthermore, a recent study aimed at understanding the origins of the Delorme 1 circumbinary planet showed that inwards migration of the planet is often reversed through interaction with a gravitationally unstable disc (M. Teasdale & D. Stamatellos 2024). Hence, a tailored study will be required to determine whether the above scenario for BEBOP-4b's formation is actually viable.

Comparing observations of transiting planets with masses in the range $20 M_{⊕} < M_p < 20 M_{jup}$ around single stars with planetary interior models demonstrates that gas giant planets are enriched with heavy elements, and the degree of enrichment depends on the mass of the planet, consistent with expectations from the core accretion scenario for planet formation (N. Miller & J. J. Fortney 2011; D. P. Thorngren et al. 2016). Hence, it seems highly probable that at least some of the observed massive super-Jovian planets formed via core

accretion. If formation occurs sufficiently far from the central binary, where its gravitational perturbations are small enough that pebble and planetesimal accretion are not inhibited (S.-J. Paardekooper et al. 2012; A. Pierens et al. 2021), then it is plausible that circumbinary super-Jovian planets such as BEBOP-4b can also form through this mechanism, with their final orbital locations being determined by post-formation disc-driven migration (A. Pierens & R. P. Nelson 2008) and the timing of disc dispersal. The formation of such a massive planet faces significant obstacles, which may explain why these bodies are rare compared to lower mass planets. For example, gas giant planets open gaps in their protoplanetary discs and for planets that exceed Jupiter's mass, efficient tidal truncation of the disc and deep gap formation can reduce the gas accretion rate onto the planet, making it difficult to accumulate a sufficiently massive gas envelope within the disc lifetime (e.g. G. Bryden et al. 1999; R. P. Nelson, E. Lega & A. Morbidelli 2023). However, as a planet grows beyond ~3 M<sub>jup</sub>, both the disc and the planet's orbit can start to become eccentric (J. C. B. Papaloizou, R. P. Nelson & F. Masset 2001; W. Kley & G. Dirksen 2006; E. Ragusa et al. 2018), reducing the effectiveness of gap opening in quenching gas accretion onto the planet and providing the possibility of massive gas envelopes accreting on Myr time-scales. This disc-driven eccentricity growth may also help to explain the observed eccentricity of BEBOP-4b, although the inferred value of $e_b = 0.43$ is larger than the values of $e \sim 0.3$ that have been obtained in hydrodynamical simulations to date. However, larger values of $e_b$ might be achieved by varying the conditions in the disc such as the temperature and viscosity. The driving of eccentricity through disc interaction is attractive when considering the probability that BEBOP-4b is coplanar with the central binary. The most likely alternative scenario for exciting BEBOP-4b's eccentricity is planet-planet scattering, and this would be expected to increase both the eccentricity and the inclination relative to the binary orbit plane, reducing the likelihood for a transit to be detected. Finally, we note that unlike the transit-discovered circumbinary planets, BEBOP-4b is not orbiting especially close to the instability zone near the binary (see Fig. 3), and hence does not obviously fit the picture considered by A. Pierens & R. P. Nelson (2008) and A. B. T. Penzlin, W. Kley & R. P. Nelson (2021) of having migrated from further out and then being parked at the edge of the cavity formed by the binary. However, the large mass of BEBOP-4b, combined with the likelihood it experienced significant disc-driven eccentricity excitation, likely led to significant sculpting of the disc, resulting in BEBOP-4b and the central binary orbiting within a common cavity. Future simulations will be required to determine what the implications of this would have been for the evolution of BEBOP-4b.

## 5.3 Orbital stability in the system

The proximity of a circumbinary planet's orbit to the stability limit can be formulated in a variety of ways, which can lead to different interpretations. We use the parameter $a_{sc} = a_b/r_{stab}$ (T. A. Baycroft et al. 2025a), where $a_b$ is the semimajor axis of the planet's orbit and $r_{stab}$ is the radius of the instability zone/the critical semimajor axis. Using the formula for the critical semimajor axis from M. J. Holman & P. A. Wiegert (1999) which takes into account $a_{bin}$, $e_{bin}$, and $q_{bin}$ gives a value of $a_{sc} \approx 2.7$. Instead using the formula from N. Georgakarakos et al. (2024), which includes more parameters (notably the eccentricity of the planet's orbit) results in $a_{sc} \approx 1.3-1.7$ (the two values being the inner and outer stability criteria defined in N. Georgakarakos et al. 2024). Instead of $a_{sc}$ using the M. J. Holman & P. A. Wiegert (1999) formula, we can incorporate the
<figure><img src="image_7.png" /></figure>

**Figure 6.** Stability analysis of the brown dwarf in the BEBOP-4 system assuming coplanar orbits. For fixed initial conditions (Table 1), the parameter space of the system is explored by varying the orbital period $P_{\text{inner}}$ and the eccentricity $e_{\text{inner}}$ of the brown dwarf. The step size is 0.5 d in orbital period and 0.005 in eccentricity. For each initial condition, the system is integrated over $10^5$ yr and a stability indicator is calculated which involved a frequency analysis of the mean longitude of the ibrown dwarf. Chaotic diffusion is indicated when the main frequency of the mean longitude varies in time. Yellow points correspond to highly unstable orbits, while purple points correspond to orbits which are likely to be stable on a billion-years time-scales. The dashed black lines indicate mean-motion resonances with the binary stars. The dot gives the position of the best-fitting solution (Table 1).

eccentricity of the brown dwarf's orbit by calculating $r_{\text{peri,sc}} = \frac{r_{\text{peri,b}}}{r_{\text{stab}}}$, where $r_{\text{peri,b}}$ is the pericentre distance of a circumbinary planet's orbit. This yields a value of $r_{\text{peri,sc}} \approx 1.5$. BEBOP-4b can therefore be considered close or far from the stability limit depending on how its eccentricity is dealt with.

To get a more definitive view about the stability of the current system, we performed a global dynamical analysis in the same way as achieved for other circumbinary planetary systems (eg. A. C. M. Correia et al. 2005; M. R. Standing et al. 2023). We used the symplectic integrator SABAC4 (J. Laskar & P. Robutel 2001), with a step size of 0.005 yr and general relativity corrections. Each initial condition is integrated for 100 kyr, and a stability indicator, $\Delta = |1 - n'/n|$, is computed, where $n$ and $n'$ are the main frequency of the mean longitude over two consecutive time intervals of 50 kyr, calculated via frequency analysis (J. Laskar 1990, 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}$.

In Fig. 6, we vary the orbital period and the eccentricity of the brown dwarf in the vicinity of the best fit (Table 1). We observe that the system is stable up to eccentricities around 0.5, which is very close to the observed eccentricity ($e = 0.435$). Moreover, there are some high order mean motion resonances with the binary nearby, namely the 25/1, that also introduce some chaotic structures for lower eccentricities. We conclude that the best-fitting system (given by a blue dot) is stable, but indeed near its stability limit, which raises interesting questions on how this system has formed. In Fig. 7, we show how the eccentricity of the brown dwarf's orbit evolves over time. The amplitude of the variations in eccentricity exceed the statistical uncertainty on the measured value of the mean eccentricity from the kima analysis.

A recent comparison between the transit-detected circumbinary planets, and the radial velocity-detected circumbinary planets shows tension in their distribution relative to the stability limit (M. J. Holman & P. A. Wiegert 1999; E. Pilat-Lohinger, B. Funk & R. Dvorak 2003; N. Georgakarakos et al. 2024). Whilst the innermost,

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

**Figure 7.** Evolution of the eccentricity of BEBOP-4b. The system is initiated with the maximum-likelihood parameters from the kima analysis, the eccentricity of which ($e = 0.435$) is marked. The uncertainty in the measured eccentricity is shown (centred on the maximum-likelihood value) as a shaded region.

transiting planets are almost always detected near the stability limit (D. V. Martin & A. H. M. J. Triaud 2014; G. Li, M. J. Holman & M. Tao 2016), radial velocity system do not seem to show such preference (T. A. Baycroft et al. 2024). One hypothesis to explain this discrepancy is for radial velocity-detected systems to host an as-of-yet undetected, low mass (but large radius) exoplanet near the stability limit. In this context, we study here whether there are dynamically stable locations between the binary and the brown dwarf for an additional planet to subsist, which might inform further data collections.

In order to check if additional planets can be hidden in the inner system, we performed again a global dynamical analysis. We assumed the best-fitting solution for the already known bodies in the system (Table 1), and scanned all orbital periods and eccentricities of a tentative inner coplanar planet between the binary stars and the brown dwarf with a radial velocity signature of 5 m s⁻¹. We observed that all orbits are unstable, and thus conclude that no other planet can subsist in the inner system.

## 6 CONCLUDING REMARKS

Thanks to 84 spectra obtained with the SOPHIE high-resolution spectrograph, we discovered BEBOP-4b, a planetary-mass brown dwarf with $m_b \sin i_b = 20.9 \pm 1.3$ M$_{\text{Jup}}$, the most massive amongst circumbinary planets confirmed via the transit or radial velocity methods. Based on the binary's apsidal precession we also expect the true mass remains $< 26.3$ M$_{\text{Jup}}$. The circumbinary brown dwarf occupies a $1823.5_{-4.9}^{+5.1}$ d orbit, with a large eccentricity $e = 0.428 \pm 0.016$, the highest yet measured amongst circumbinary planets. The binary has masses $M_A = 1.51$ M$_{\odot}$, and $M_B = 0.46$ M$_{\odot}$, an orbital period $P_{\text{bin}} = 72$ d and $e_{\text{bin}} = 0.27$. BEBOP-4 AB is the longest period binary known to host a circumbinary planet detected with transit or radial velocities. At this stage, BEBOP-4 (AB) b is unremarkable when compared to known exoplanets orbiting single stars, but unusual within the circumbinary context. Because planet formation is expected to be affected by binarity, it is surprising that such a high mass object could assemble. As surveys for circumbinary exoplanets expand their monitoring, and thanks to Gaia, it will soon be possible to know whether such high-mass companions are a regular occurrence.
ACKNOWLEDGEMENTS

We thank the reviewer for their prompt work and thoughtful remarks. These scientific results would not have been possible without a generous allocation of telescope time by France Programme National de Planétologie, and without the dedication of staff working at the Observatoire de Haute-Provence, particularly during the COVID-19 pandemic. Four measurement were obtained during that period, which are essential for our conclusions. We would also like to thank the help of other observers part of the timeshare agreement between groups using the SOPHIE instrument. 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 no. 803193/BEBOP) and by the ERC/UKRI (United Kingdom Research and Innovation) Frontier Research Guarantee programme (EP/Z000327/1/CandY). ACMC acknowledges support from FCT–Fundação para a Ciência e a Tecnologia, I.P., Portugal, through the CFisUC projects UID/04564/2025. Computations described in this paper were performed using the University of Birmingham’s Blue-BEAR High-Performance Computing service. 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 SOPH IE archive, and were obtained under Prog.ID 18B.DISC.TRIA and 19A.PNP.SANT, and under target name J1250+26. Reduced and drift-corrected radial velocities are found in Table A1, a full version of which can be found online.

REFERENCES

Almeida L. A., Jablonski F., Rodrigues C. V., 2013, ApJ, 766, 11

Andrews S. M. et al., 2018, ApJ, 869, L41

Armstrong D. J., Osborn H. P., Brown D. J. A., Faedi F., Gómez Maqueo
Chew Y., Martin D. V., Pollacco D., Udry S., 2014, MNRAS, 444, 1873

Baştürk Ö., Esmer E. M., Demir E., Selam S. O., 2023, Bull. Soc. R. Sci. Liege, 92, 11197

Bailey V. et al., 2014, ApJ, 780, L4

Baraffe I., Chabrier G., Barman T. S., Allard F., Hauschildt P. H., 2003, A&A, 402, 701

Baraffe I., Homeier D., Allard F., Chabrier G., 2015, A&A, 577, A42

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

Baruteau C., Meru F., Paardekooper S.-J., 2011, MNRAS, 416, 1971

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

Baycroft T. A., Triaud A. H. M. J., Lai D., 2024, in Lemaitre A., Libert A.-S., eds, Proc. IAU Symp. 382, Complex Planetary Systems II: Latest Methods for an Interdisciplinary Approach. Cambridge Univ. Press, Cambridge, p. 51

Baycroft T. A. et al., 2025a, MNRAS, 541, 2801

Baycroft T. A., Sairam L., Triaud A. H. M. J., Correia A. C. M., 2025b, Sci. Adv., 11, eadu0627

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

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

Beuermann K. et al., 2010, A&A, 521, L60

Beuermann K., Dreizler S., Hessman F. V., Deller J., 2012, A&A, 543, A138

Boley A. C., Hayfield T., Mayer L., Durisen R. H., 2010, Icarus, 207, 509

Borkovits T., Rappaport S., Hajdu T., Sztakovics J., 2015, MNRAS, 448, 946

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

Bryden G., Chen X., Lin D. N. C., Nelson R. P., Papaloizou J. C. B., 1999, ApJ, 514, 344

Calissendorff P., Janson M., Asensio-Torres R., Köhler R., 2019, A&A, 627, A167

Carmichael T. W., Latham D. W., Vanderburg A. M., 2019, AJ, 158, 38

Coleman G. A. L., Nelson R. P., 2016, MNRAS, 460, 2779

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

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

Columba G., Daniëlski C., Dorozsmai A., Toonen S., Lopez Puertas M., 2023, A&A, 675, A156

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

Couetdic J., Laskar J., Correia A. C. M., Mayor M., Udry S., 2010, A&A, 519, A10

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

Cumming A., Butler R. P., Marcy G. W., Vogt S. S., Wright J. T., Fischer D. A., 2008, PASP, 120, 531

Dalal S. et al., 2021, A&A, 651, A11

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

Delorme P. et al., 2013, A&A, 553, L5

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

Dupuy T. J. et al., 2018, AJ, 156, 57

Dupuy T. J., Liu M. C., Evans E. L., Best W. M. J., Pearce L. A., Sanghi A., Phillips M. W., Bardalez Gagliuffi D. C., 2023, MNRAS, 519, 1688

Emsenhuber A., Mordasini C., Burn R., Alibert Y., Benz W., Asphaug E., 2021, A&A, 656, A69

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

Fulton B. J. et al., 2021, ApJS, 255, 14

Georgakarakos N., Eggl S., Ali-Dib M., Dobbs-Dixon I., 2024, AJ, 168, 224

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

Grether D., Lineweaver C. H., 2006, ApJ, 640, 1051

Han C. et al., 2024, A&A, 685, A16

Hardy A. et al., 2015, ApJ, 800, L24

Heidari N., 2022, PhD Thesis, Université Côte d'Azur; Shahid Beheshti University (Tehran)

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

Hidalgo S. L. et al., 2018, ApJ, 856, 125

Holl B. et al., 2023, A&A, 674, A10

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

Kiefer F. et al., 2019, A&A, 631, A125

Kley W., Dirksen G., 2006, A&A, 447, 369

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

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

Kozai Y., 1962, AJ, 67, 591

Kumar S. S., 1963, ApJ, 137, 1121

Laskar J., 1990, Icarus, 88, 266

Laskar J., 1993, Physica D Nonlin. Phenom., 67, 257

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

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

Lidov M. L., 1962, Planet. Space Sci., 9, 719

Ma B., Ge J., 2014, MNRAS, 439, 2781

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

Martin D. V., Triaud A. H. M. J., 2016, MNRAS, 455, L46

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

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

Mayor M. et al., 2011, preprint (arXiv:1109.2497)

Miller N., Fortney J. J., 2011, ApJ, 736, L29

Nelson R. P., Lega E., Morbidelli A., 2023, A&A, 670, A113

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

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

Papaloizou J. C. B., Nelson R. P., Massett F., 2001, A&A, 366, 263

Pascucci I. et al., 2016, ApJ, 831, 125
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

Pilat-Lohinger E., Funk B., Dvorak R., 2003, A&A, 400, 1085

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

Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M., Greenzweig Y., 1996, Icarus, 124, 62

Prša A. et al., 2016, AJ, 152, 41

Qian S. B., Liao W. P., Zhu L. Y., Dai Z. B., Liu L., He J. J., Zhao E. G., Li L. J., 2010, MNRAS, 401, L34

Queloz D. et al., 2001, A&A, 379, 279

Ragusa E., Rosotti G., Teyssandier J., Booth R., Clarke C. J., Lodato G., 2018, MNRAS, 474, 4460

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

Rosenthal M. M., Murray-Clay R. A., Perets H. B., Wolansky N., 2018, ApJ, 861, 74

Sahlmann J., Ségransan D., Queloz D., Udry S., 2011, in Sozzetti A., Lattanzi M. G., Boss A. P., eds, Proc. IAU Symp. 276, The Astrophysics of Planetary Systems: Formation, Structure, and Dynamical Evolution. Cambridge Univ. Press, Cambridge, p. 117

Sahlmann J., Lazorenko P. F., Ségransan D., Martín E. L., Mayor M., Queloz D., Udry S., 2014, A&A, 565, A20

Sahlmann J., Triaud A. H. M. J., Martin D. V., 2015, MNRAS, 447, 287

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

Sairam L. et al., 2024b, MNRAS, 534, 3999

Santerne A. et al., 2016, A&A, 587, A64

Schneider J., Dedieu C., Le Sidaner P., Savalle R., Zolotukhin I., 2011, A&A, 532, A79

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

Sigurdsson S., Richer H. B., Hansen B. M., Stairs I. H., Thorsett S. E., 2003, Science, 301, 193

Stamatellos D., Whitworth A. P., 2009, MNRAS, 392, 413

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

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

Stevenson A. T., Haswell C. A., Barnes J. R., Barstow J. K., 2023, MNRAS, 526, 5155

Teasdale M., Stamatellos D., 2024, MNRAS, 533, 2294

Thorngren D. P., Fortney J. J., Murray-Clay R. A., Lopez E. D., 2016, ApJ, 831, 64

Thorsett S. E., Arzoumanian Z., Taylor J. H., 1993, ApJ, 412, L33

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

Wolf M., Kučáková H., Zasche P., Hornoch K., Kára J., Merc J., Zejda M., 2021, A&A, 647, A65

Zhang Z., Fabrycky D. C., 2019, ApJ, 879, 92

Zhu Z., Hartmann L., Nelson R. P., Gammie C. F., 2012, ApJ, 746, 110

Zorotovic M., Schreiber M. R., 2013, A&A, 549, A95

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

### TableA1.txt

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: ADDITIONAL FIGURES AND TABLES

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

**Figure A1.** A comparison of the mass estimates for BEBOP-4 A. The results from Section 3 and the HRCCS analysis with SALTIRE are consistent with inference from BASTA from the observed [Fe/H], $T_{eff}$, and G magnitude using the BaSTI isochrones (S. L. Hidalgo et al. 2018) with microphysics following the recommended science case for this mass range, and with the MIST isochrone value reported by A. V. Freckelton et al. (2024) for BEBOP primary stars. The TESS input catalogue (TIC) mass is almost fully inconsistent with other estimates (with a discrepancy in the median of ~0.5 M<sub>☉</sub>), highlighting the need to derive accurate stellar parameters. Mock distributions are randomly generated from values and their $1\sigma$ uncertainties to visualize the degree of consistency between estimates.
<figure><img src="image_10.png" /><figcaption>**Figure A2.** Top panel: Radial velocity observations (dots) as a function of time, with the highest likelihood model represented (line). Middle panel: residuals after removing the highest likelihood model except for $\dot{\omega}_{\text{bin}}$ fixed to 0. The line in that panel reveals the contribution of a non-zero $\dot{\omega}_{\text{bin}}$ to the fit. Bottom plot: residuals after removing the highest likelihood model (including $\dot{\omega}_{\text{bin}}$).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure A3.** Radial velocity Keplerian for the binary star (secondary's orbit around the primary). Residuals (including companion removed) are shown below.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure A4.** Distributions and correlations of the brown dwarf's parameters. The lines show the values of the maximum-likelihood solution which is used to generate the orbital and radial velocity phase plots as well as to instantiate the N-body simulations.</figcaption></figure>

<table><caption>**Table A1.** Radial velocity observations taken with SOPHIE between 2019 January 7 and 2025 February 27. Full table can be found in online version of paper.</caption><thead><tr><th>Time<br>[BJD - 2400000]</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><th>FWHM<br>(km s<sup>-1</sup>)</th><th>Bisector span<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>58491.71259</td><td>-8.73180</td><td>0.01170</td><td>10.62970</td><td>-0.05000</td></tr><tr><td>58492.64226</td><td>-8.03090</td><td>0.00350</td><td>10.44850</td><td>0.00190</td></tr><tr><td>58537.49195</td><td>3.92270</td><td>0.00180</td><td>10.49780</td><td>0.00370</td></tr><tr><td>58538.63542</td><td>2.99850</td><td>0.00330</td><td>10.47860</td><td>-0.01230</td></tr><tr><td>58539.57288</td><td>2.24080</td><td>0.00200</td><td>10.46250</td><td>0.00260</td></tr><tr><td>58542.53512</td><td>-0.07170</td><td>0.00190</td><td>10.44340</td><td>-0.00160</td></tr><tr><td>58554.56388</td><td>-8.18770</td><td>0.00410</td><td>10.45980</td><td>0.00050</td></tr><tr><td>58559.53582</td><td>-9.71920</td><td>0.00210</td><td>10.48910</td><td>0.00480</td></tr><tr><td>...</td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>**Table A2.** 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 L. R. Doyle et al. (2011). $^y$ Transit detection and parameters from V. B. Kostov et al. (2020). $^z$ Transit detection from V. B. 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>pi</sub><br>(M<sub>☉</sub>)</th><th>M<sub>sec</sub><br>(M<sub>☉</sub>)</th><th>ε<sub>bin</sub></th><th>P<sub>pl</sub><br>(d)</th><th>M<sub>pl</sub><br>(M<sub>sup</sub>)</th><th>ε<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>A. H. M. J. Triaud et al. (2022)</td><td>299096355</td></tr><tr><td>BEBOP-1 b</td><td>TOL-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>M. R. Standing et al. (2023)</td><td>260128333</td></tr><tr><td>BEBOP-1 c</td><td>TOL-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>M. R. 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>L. Sairam et al. (2024a)</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.2615(39)</td><td>0.063255(54)</td><td>547.0<sup>(+6.2)</sup><br>(-7.6)</td><td>0.558<sup>(+51)</sup><br>(-48)</td><td>0.247<sup>(+77)</sup><br>(-89)</td><td>T. A. Baycroft et al. (2025a)</td><td>28994953</td></tr><tr><td>BEBOP-4 b</td><td>HD 111065 b</td><td>71.96857(15)</td><td>1.507(152)</td><td>0.459(28)</td><td>0.26887(11)</td><td>1823.5(5.0)</td><td>20.9(1.3)</td><td>0.428(16)</td><td>This work</td><td>356710041</td></tr></tbody></table>

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

© The Author(s) 2025.

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.