key: cord-0264580-cqlweafg
authors: Triaud, Amaury H.M.J.; Standing, Matthew R.; Heidari, Neda; Martin, David V.; Boisse, Isabelle; Santerne, Alexandre; Correia, Alexandre C.M.; Acuna, Lorana; Battley, Matthew; Bonfils, Xavier; Carmona, Andr'es; Cameron, Andrew Collier; Cort'es-Zuleta, P'ia; Dransfield, Georgina; Dalal, Shweta; Deleuil, Magali; Delfosse, Xavier; Faria, Joao; Forveille, Thierry; Hara, Nathan C.; H'ebrard, Guillaume; Hoyer, Sergio; Kiefer, Flavien; Kunovac, Vedad; Maxted, Pierre F. L.; Martioli, Eder; Miller, Nikki; Nelson, Richard P.; Poveda, Mathilde; Rein, Hanno; Sairam, Lalitha; Udry, St'ephane; Willett, Emma
title: BEBOP III. Observations and an independent mass measurement of Kepler-16 (AB) b -- the first circumbinary planet detected with radial velocities
date: 2021-12-13
journal: nan
DOI: 10.1093/mnras/stab3712
sha: a8bda9e7dd4756c29e92f2ae0460fab1866ac494
doc_id: 264580
cord_uid: cqlweafg

The radial velocity method is amongst the most robust and most established means of detecting exoplanets. Yet, it has so far failed to detect circumbinary planets despite their relatively high occurrence rates. Here, we report velocimetric measurements of Kepler-16A, obtained with the SOPHIE spectrograph, at the Observatoire de Haute-Provence's 193cm telescope, collected during the BEBOP survey for circumbinary planets. Our measurements mark the first radial velocity detection of a circumbinary planet, independently determining the mass of Kepler-16~(AB)~b to be $0.313 pm 0.039,{rm M}_{rm Jup}$, a value in agreement with eclipse timing variations. Our observations demonstrate the capability to achieve photon-noise precision and accuracy on single-lined binaries, with our final precision reaching $rm 1.5~m,s^{-1}$ on the binary and planetary signals. Our analysis paves the way for more circumbinary planet detections using radial velocities which will increase the relatively small sample of currently known systems to statistically relevant numbers, using a method that also provides weaker detection biases. Our data also contain a long-term radial velocity signal, which we associate with the magnetic cycle of the primary star.

Circumbinary planets are planets that orbit around both stars of a binary star system. Long postulated (Borucki & Summers 1984; Schneider 1994) , the first unambiguous discovery of a circumbinary planet came with Kepler-16 (Doyle et al. 2011) , detected by identifying three transits within the lightcurve of an eclipsing binary system monitored by NASA's Kepler mission (Borucki et al. 2011) . Kepler went on to detect another 13 transiting circumbinary planets, in 11 systems (Martin 2018; Socia et al. 2020) , with another two systems found using TESS (Kostov et al. 2020 (Kostov et al. , 2021 . A number of circumbinary planets systems are suspected from eclipse timing variations of binaries on the main sequence (e.g. Borkovits et al. 2016; Getley et al. 2017) , and stellar remnants (e.g. Marsh et al. 2013; Han et al. 2017) but most are disputed (e.g. Mustill et al. 2013) , and some disproven (e.g. Hardy et al. 2015) . Other detections include HD 106906 b, in direct imaging (Bailey et al. 2014 ) and OGLE-2007-BLG-349L(AB)c with the microlensing method (Bennett et al. 2016) .

Despite successes with almost every observational methods, no circumbinary planet signal has been detected using radial velocities yet. In addition, radial velocities have detected many planets with masses compatible with currently known circumbinary planets. This is remarkable since radial velocities are one of the earliest, most established and efficient method of exoplanet detection. The system closest to a circumbinary configuration identified thus far is HD 202206 (Correia et al. 2005, see Sect. 4) .

The radial velocity method has a number of advantages over the transit method. First, it is less restrictive in term of the planet's orbital inclination thus providing a weaker bias towards short orbital periods. Furthermore, the signal can be obtained at every orbital phase, and the method is more cost-effective and easier to use over a longer term thanks to using ground-based telescopes . Additionally, radial velocities provide the planet's mass, its most fundamental parameter. While the transit method can provide a mass when eclipse timing variations are detected, most circumbinary exoplanets unfortunately remain without a robust mass determination with eclipse timing variations mostly providing upper limits (e.g. Orosz et al. 2012; Schwamb et al. 2013; Kostov et al. 2020) . Only four of the known circumbinary planets have eclipse-timing mass estimates inconsistent with 0 at > 3 . The present and future TESS and PLATO missions (Ricker et al. 2014; Rauer et al. 2014 ) are set to identify several more transiting circumbinary planet candidates (e.g. Kostov et al. 2020 Kostov et al. , 2021 . However, these are unlikely to produce many reliable mass measurements, in good part due to rather short observational timespans compared to Kepler's.

Overall, radial velocities are essential to create a sample of circumbinary planets that is both greater in number and less biased than the transit sample. This will allow a deeper understanding of circumbinary planets: their occurrence rate (Martin & Triaud 2014; Armstrong et al. 2014) , multiplicity (Sutherland & Kratter 2019; Orosz et al. 2019) , formation and evolution (e.g. Chachan et al. 2019; Pierens et al. 2020; Penzlin et al. 2021) , and dependence on binary properties (Muñoz & Lai 2015; Martin et al. 2015; Li et al. 2016; .

In 2017 we created the BEBOP survey (Binaries Escorted By Orbiting Planets; , as a blind radial velocity survey for circumbinary planets. Prior to this, the most extensive radial velocity effort had been produced by the TATOOINE survey (Konacki et al. 2009 ), but the survey unfortunately did not yield any discoveries. One issue likely affected the survey, from which BEBOP learnt a great deal: TATOOINE targeted double-lined binaries, which was logical. Double-lined binaries are brighter, both stars can have model-independent mass measurements, and one can in principle measure the Doppler displacement caused by a planet on each of the two components. However, disentangling both components from their combined spectrum accurately is a complex task, and despite photon noise uncertainties regularly reaching 2 to 4 m s −1 , the survey returned a scatter of order 15 to 20 m s −1 (Konacki et al. 2009 (Konacki et al. , 2010 . Indeed, Konacki et al. (2010) recommended single-lined binaries as a solution, but too few were known at the time. Since, radial velocities have been used to constrain the binary parameters, which helps refining the planetary parameters, but not to search for circumbinary planets themselves (e.g. Kostov et al. 2013 Kostov et al. , 2014 .

Thanks to the advent of exoplanet transit experiments, increasing amounts of low-mass, single-lined eclipsing binaries are being identified (e.g. Triaud et al. 2013 Triaud et al. , 2017 von Boetticher et al. 2019; Lendl et al. 2020; Mireles et al. 2020; Acton et al. 2020) . The BEBOP survey was constructed solely using sufficiently faint secondaries, to avoid detection with spectrographs, such that we could in principle reach a radial velocity precision comparable to that around single stars of the same brightness. In principle, this ought to provide an accuracy of order 1 m s −1 . Our survey is ongoing, and uses the CORALIE, SOPHIE, HARPS, and ESPRESSO spectrographs. Preliminary results were published in (BEBOP I). In Standing et al. (2021, under review, BEBOP II) , we describe our observational protocol, the methods we use to detect planets, as well as how we produce detection limits.

In this paper we detail a complementary project to BEBOP's blind search. Between 2016 and 2021 we monitored Kepler-16, a relatively bright (Vmag = 12) single-lined eclipsing binary system with a primary mass 1 = 0.65 M (a K dwarf), a secondary mass 2 = 0.20 M (a mid M dwarf), and an orbital period bin = 41.1 days. The system is 75 pc distant and known to host a circumbinary gas giant planet with a mass pl = 0.33 M Jup , and a period pl = 229 days. Our observations demonstrate that we can indeed recover the Doppler reflex signature of a circumbinary planet. Our results act to both validate and assist our broader search for new planets. Furthermore, we can derive a "traditional" Doppler mass measurement for the planet, to be compared with that derived from photometric eclipse and transit timings. Finally, our long baseline is sensitive to additional planets, in particular to any that would occupy an orbit misaligned to the transiting inner planet's.

Between 2016-07-08 and 2021-06-23, we collected 143 spectra using the high-resolution, high-precision, fibre-fed SOPHIE spectrograph, mounted on the 193cm at Observatoire de Haute-Provence, in France (Perruchot et al. 2008) . The Journal of Observations can be found in Table A1 . All observations were conducted in HE mode (High Efficiency) where some of the instrumental resolution is sacrificed from 75,000 to 40,000 in favour of a 2.5× greater throughput. We chose this since whilst Kepler-16 is the brightest circumbinary system, it is relatively faint for SOPHIE, with ∼ 12.0. SOPHIE has two fibres; the first stayed on target, while the other was kept on the sky in order to remove any contribution from the Moon-reflected sunlight. Standard calibrations were made at the start of night as well as roughly every two hours throughout the night to monitor the instrument's zero point. In addition, we observed one of three standards (HD 185144, HD 9407 and HD 89269 A) in HE mode nightly, which we used to track and correct for any long-term instrumental drift following procedures established in Santerne et al. (2014) and Courcol et al. (2015) .

Our radial velocities were determined by cross-correlating each spectrum with a K5V mask. These methods are described in Baranne et al. (1996) , and Courcol et al. (2015) , and have been shown to produce precisions and accuracies of a few meters per second (e.g. Bouchy et al. 2013; Hara et al. 2020) , well below what we typically obtain on this system. As in Baranne et al. (1996) , and Pollacco et al. (2008) , we correct our data from lunar contamination by first scaling the calculated CCF (cross-correlation function) on fibre A and B (to account for slightly different efficiencies between the two fibres) before subtracting the two CCFs. This is a particularly important procedure for circumbinary planet searches. Most systems observed with SOPHIE are single stars, and the scheduling software informs the observer whether sunlight reflected on the Moon would create a parasitic cross-correlation signal, with a radial velocity that varies predictably with the lunar phase. In such a situation the observation is postponed. In the case of binary observations, ours, the velocity of the primary star keeps changing by km s −1 , meaning we could not practically predict possible lunar contamination at the time of acquisition. We also correct our data from the CTI (Charge Transfer Inefficiency) effect following the procedure described in Santerne et al. (2012) .

The cross-correlation software produces two key metrics of the shape of the CCF, its FWHM (Full Width at Half Maximum), and its bisector span (as defined in Queloz et al. 2001 ). In addition, we measure the H stellar activity indicator following Boisse et al. (2009) . These are provided in Table A1 .

Two measurements are immediately excluded from our analysis. On 2017-09-06 (BJD 2,458,003.32008) and 2018-10-06 (BJD 2,458,398.33382), when the fibre was mistakenly placed onto another star. This is obvious from the FWHM we extract from these measurement, and from their radial velocity. They are appropriately flagged in Table A1 . In the end, we achieve a mean radial velocity precision of ∼ 10.6 m s −1 on the remaining 141 measurements.

Radial velocity measurements of the primary have also been obtained with TRES and Keck's HIRES (Doyle et al. 2011; Winn et al. 2011 ). These datasets were not used in this analysis for several reasons. First, the TRES data has a mean precision of 21 m s −1 , which would be insufficient to detect the planet. Second, the HIRES data, despite offering a precision of a few m s −1 , was only taken on a single night and is contaminated by the Rossiter-McLaughlin signal of the eclipse. Finally, by solely using our own SOPHIE data we may produce a near independent detection. In a similar spirit, we do not use any of Kepler's or TESS's photometric data to conduct our search and parameter estimation.

Finally, Bender et al. (2012) collected near-infrared spectra to reveal the secondary's spectral lines (i.e. observing Kepler-16 as a double-lined spectroscopic binary). Similarly we do not use these measurements in our analysis, although we do use their estimate for the primary mass.

In order to ascertain our capacity to detect circumbinary planets using radial velocities, as an independent method, we decided to use two different algorithms with different methods for measuring a detection probability. Before describing this, we will detail our procedure to remove outliers. As a reminder, we only use SOPHIE data (see above), where Kepler-16 appears as a single-lined binary and we only observe the displacement of the primary star around the system's barycentre.

We searched our 141 data for measurements 1 coinciding with a primary eclipse, and likely affected by the Rossiter-McLaughlin effect (Rossiter 1924; McLaughlin 1924; Winn et al. 2011; Triaud 2018) . Fortunately no measurements needed to be excluded for this reason.

We also realised that a number of measurements were likely taken under adverse conditions. This is apparent from an unusually low signal-to-noise ratio, but also from large values of the bisector span. Typically used as a stellar activity, or a blend indicator (Queloz et al. 2001; Santos et al. 2002) , the span of the bisector slope (bisector span, or bis span) effectively informs us that the line shape varied and therefore that the mean of the cross-correlation function is likely affected. We took the mean of the bisector span measurement and removed all measurements in excess of 3 away from the mean. Five measurements are excluded this way, all with a bisector span ±100 m s −1 . Excluded measurements are reported in the Journal of Observations in Table A1 with a flag.

For visual convenience we also exclude one measurement taken on 2018-06-02 (BJD 2,458,271.53928) on account of its very low signalto-noise and correspondingly large uncertainty, seven times greater than the semi-amplitude of Kepler-16 b. Again, this is reported in Table A1 with a flag.

Finally we remove a measurement obtained on 2017-10-30 (BJD 2,458,057.38946). This measurement is ∼ 6 away from the best fit model. It is totally unclear why this is the case since its FWHM, bisector span, and H , all appear compatible with other measurements.

After removing these seven outliers, our analysis is performed on the remaining 134 SOPHIE measurements. The exoplanet detection described just below is done twice, without, and with the seven outlying measurements. Their exclusion did not affect our conclusion but refined our parameters. We also reproduced the following fitting procedures by including the previously existing TRES and HIRES data with no discernible differences.

is a radial velocity fitting tool used for exoplanet detection, described in Ségransan et al. (2011) . It assumes Keplerian orbits (see Sect 3.5). Y first performs a Lomb-Scargle periodogram of the observations, which is then used to initiate a genetic algorithm that iterates over the orbital period , the eccentricity , the argument of periastron , the semi-amplitude , a reference time 0 , and one systemic velocity per dataset (for conventions, see Hilditch 2001) . Once the algorithm has converged on a best solution, the final parameters are estimated by using a least-square fit. The tool has been routinely used to identify small planetary companions successfully (e.g. Mayor et al. 2011; Bonfils et al. 2013) , and represents a more traditional, and possibly a more recognised way of identifying a new planetary system than the nested sampler we use subsequently (see Sect. 3.3). However Y cannot do a Bayesian model comparison. Instead it computes a False Alarm Probability (FAP) by performing a bootstrap on the data thousands of times and computing for each iteration a Lomb-Scargle periodogram.

In a first instance, just one Keplerian is adjusted to the SOPHIE 1 The first series of 14 measurements were obtained from a catalogue containing erroneous proper motions and epochs, which in turn created a 1.5 km s −1 effect in the radial velocity as the Earth's motion was over compensated. This can be corrected easily, but the error remains within the archival data. data, with Y automatically finding the most prominent signal, that of the secondary star. Following this step, we remove the secondary's signal and search the resulting residuals with a periodogram. This periodogram shows excess power around 230 day with a FAP ∼ 0.1% as well as excess power for a signal longer than the range of dates we observed for, which we later associate to a magnetic cycle. To isolate the signal of Kepler-16 b, we fit the secondary's Keplerian to the data, alongside a polynomial function, which is used to detrend that longer signal. Once Y has converged, we search the residuals again with a periodogram, which provides a FAP 0.01% (Fig. 1) , clearly detecting Kepler-16 b as an additional periodic sinusoidal signal. The FAP obtained with a cubic detrending function is one order of magnitude better than that obtained with a quadratic function, so we chose the former as our baseline detrending.

To obtain results on the system's parameters, we perform a final fit to the data assuming two Keplerians and a cubic detrending function. Results of that fit are found in Table 1 . The orbital parameters of the planet are compatible with those produced in Doyle et al. (2011) (see table 1 and section 3.5 for a discussion). Our final fit produces a reduced 2 = 1.17 ± 0.14, implying no additional complexity is needed to explain the data, and supporting our choice for a circular planetary orbit 2 . In addition, this shows that we can achieve photonnoise precision on single-lined binary to detect circumbinary planets. The model fit to the data is depicted in Fig. 2 .

Including seven outliers described in Sect. 3.1, neither the FAP nor the reduced 2 are significantly affected. (58) 0.7048(11) * adopted from Bender et al. (2012) 3.3 Analysis using the diffusive nested sampler K K is a tool developed by Faria et al. (2018) which fits a sum of Keplerian curves to radial velocity data. It samples from the posterior distribution of Keplerian model parameters using a Diffusive Nested Sampling algorithm by Brewer & Foreman-Mackey (2016) . Diffusive Nested sampling allows the sampling of multi-modal distributions, such as those typically found in exoplanetary science and radial velocity data (Brewer & Donovan 2015) , evenly and efficiently. K can treat the number of planetary signals ( p ) present in an RV data set as a free parameter in its fit. Since the tool also calculates the fully marginalised likelihood (evidence) it allows for Bayesian model comparison (Trotta 2008) between models with varying p . A measure of preference of one Bayesian model over another can be ascertained by computing the "Bayes Factor" (Kass & Raftery 1995) between the two. The Bayes factor is a ratio of probabilities between the two competing models. Once a value for the Bayes factor has been calculated, we can compare it to the so-called "Jeffreys' scale" (see Trotta (2008) for more details) to rate the strength of evidence of one model over another.

A more extensive description of our use of K in the context of the BEBOP survey can be found in Standing et al. (2021) . For the analysis of the Kepler-16 system prior distributions were chosen similarly to those used in Faria et al. (2020) with the following notable adaptions. We treat the secondary star as a known object with tight uniform priors on its orbital parameters. A log-uniform distribution was used to describe the periods of any additional signals, from 4 × bin to 1 × 10 4 days. This inner limit on period is set by the instability limit found in binary star systems (Holman & Wiegert 1999). More details, particularly on the priors we use, can be found in Standing et al. (2021) . Just like for the previous analysis using Y , K is only deployed on SOPHIE data, and excluding outliers described in Sec. 3.1.

Our K analysis of the Kepler-16 data yields a Bayes Factor BF > 10000 in favour of a three Keplerian model (secondary star, planet and cubic drift). Our posterior shows over-densities at orbital periods of ≈ 230 and ≈ 2000 days, corresponding to the signal of Kepler-16 b and the cubic drift seen in Y . We then apply the clustering algorithm HDBSCAN (McInnes et al. 2017) to isolate and extract the resulting planetary orbital parameters, which can be found in table 1.

To convert our semi-amplitudes into masses for 2 and pl , for the secondary star and planet masses (respectively), we adopt a primary star mass ( 1 ) from Bender et al. (2012) . Software written for exoplanetary usage usually assumes that pl ★ (including Y and K ), however this assumption is no longer valid when comparing 2 to 1 , and a circumbinary planet to both.

First we find 2 iteratively using the mass function, following the procedure described in Triaud et al. (2013) . Then, we estimate pl from 1,pl by using the combined mass 1 + 2 . This is because whilst we are only measuring the radial velocity signature of the primary star, the gravitational force of the planet acts on the barycentre of the binary. Had we not done this extra conversion step, we would find significant differences, with erroneous 2 = 0.165 M and pl = 0.29 M Jup for the Y results. Differences between our derived parameters (bottom part of table 1) and those from Doyle et al. (2011) are mainly explained by our adoption of the more accurate 1 mass from Bender et al. (2012) rather than using the value from Doyle et al. (2011) .

The main differences between our fitted parameters and those from Doyle et al. (2011) are caused by our parameters being akin to mean parameters (e.g. the mean orbital period), whereas Doyle et al. (2011) provide osculating parameters, which are the parameters the system had at one particular date, and which constantly evolve following three-body dynamics (Mardling 2013) . Also, the planetary signal is significantly more obvious in the Kepler transiting data than in our radial velocities. Each measurement within a planetary transit over the primary star produces an SNR = 243 (SNR = 14 when over the secondary). This is why Kepler can derive osculating elements, by solving Newton's equations of motion: the planetary motion is resolved orbit after orbit (transit after transit). Comparatively, our radial velocity observations have required multiple orbits of the planet to build up a significant detection. We can therefore only measure a mean period, and are justified in using software which can only adjust non-interacting Keplerian functions.

For example, Doyle et al. (2011) provide a highly precise value of pl = 228.776 ± 0.037 days, yet this osculating period will vary by approximately ±5 days over a timescale of just years. Our measured period of pl = 226.0 ± 1.7 days has a much higher error with this uncertainty being a combination of our radial velocities being less constraining on the period, and our assumption of a static orbit. The value we obtain with K is 1.6 compatible with the value found by Doyle et al. (2011) .

With respect to the planet mass, we can derive a value with similar precision to that of Doyle et al. (2011) because whilst the transit signature of the planet is much stronger than its radial velocity signature, the transit signature itself carries very little information about the planet's mass. The photodynamical mass derived in Doyle et al. (2011) is dictated by the eclipse timing variations, which have an amplitude of a couple of minutes and have a precision on the order of tens of seconds, producing an SNR close to our radial velocities.

Finally, to validate our assumption that we cannot measure osculating elements with our current data, we used tools described in Correia et al. (2005) and Correia et al. (2010) to perform a N-body fit to the radial velocities. This fit finds no improvements in 2 indicating that Newtonian effects indeed remain below the detectable threshold.

To assess the presence of an external companion causing the additional polynomial signal we notice in our data, we force a two-Keplerian model to fit the data (a circular orbit for the planet and a free-eccentricity orbit for the binary) and analyse the residuals. We find a signal reaching a FAP < 0.01% at orbital periods exceeding the timespan of our data ( 2, 000 days). Long-term drifts can sometimes be caused by magnetic cycles since stellar spots and faculae tend to suppress convective blue-shift, producing a net change in the apparent velocity of a star (Dravins 1985; Meunier et al. 2010; Dumusque et al. 2011) . We perform a Lomb-Scargle periodogram on the H activity indicator we extracted from the SOPHIE spectra, and find a FAP < 0.01%. We therefore interpret this long-term radial velocity drift as a magnetic cycle, with a timescale longer than the timespan of our observations. The validity of this interpretation can be tested since duration of stellar magnetic cycles scales with stellar rotation periods (for single star; Suárez Mascareño et al. 2016) . We measure the primary star's rotation from the sin 1 obtained during the Rossiter-McLaughlin effect by Winn et al. (2011) to the primary's stellar radius obtained by Bender et al. (2012) and obtain a primary rotation rot,1 = 35.68 ± 1.04 days. Following the relation of Suárez Mascareño et al. (2016) , with the observed stellar rotation we ought to expect a magnetic cycle on a timescale of 1900 to 2400 days, which is entirely compatible with the H signal and the long-term radial velocity drift.

We now use K as done in Standing et al. (2021) to compute a detection limit on the presence of additional but undetected planetary companions. We first remove the highest likelihood model with two Keplerians (the planet and the long-term drift) from the data, then force p = 1 to obtain a map of all remaining signals that are compatible with the data, but remain formally undetected (the binary is also adjusted at each step). This map is shown as a greyscale density on Fig. 3 . The 99% contour informs us that we are sensitive to companions below the mass of Kepler-16 b up to orbital periods of ∼ 3, 000 days. This complements Martin & Fabrycky (2021) 's work who placed detection limits down to Earth-radius planets but only out to periods of 500 days, re-analysing the Kepler photometry.

Overall, a picture is emerging of Kepler-16 b as a lonely planet, which has implications for the formation and migration of multiplanet systems in the presence of a potentially destabilising binary (e.g. Sutherland & Kratter 2019).

Our data clearly shows an independent detection of the circumbinary planet Kepler-16 b, the first time a circumbinary planet is detected using radial velocities, and the first time a circumbinary planet is detected using ground-based telescopes as well. We re-iterate that our model fits are solely made using radial velocities and completely ignore any Kepler or other photometric data, to emulate BEBOP's blind search. Importantly our results show we can achieve a precision close to 1 m s −1 on a planetary signal, with the 1 uncertainties on the semi-amplitudes being only just 1.5 m s −1 . This is compatible with the semi-amplitudes of super-Earths and sub-Neptunes (e.g. Mayor et al. 2011) . The closest previous detection to a circumbinary planet made from the ground was produced by Correia et al. (2005) on HD 202206, a system comprised of a Sun-like star with an inner companion with mass sin = 17.4 M Jup and orbital period = 256 days, and an outer companion with mass sin = 2.44 M Jup and orbital period = 1383 days. Benedict & Harrison (2017) claim an astrometric detection of the system that implies a nearly face-on system, with = 0.89 M and = 18 M Jup , suggesting a circumbinary brown dwarf. However a dynamical analysis produced by Couetdic et al. (2010) imply such a configuration is unstable and therefore unlikely, favouring instead a more edge-on system.

We first detect Kepler-16 b by using a classical approach to planet detection, via periodograms and false alarm probabilities (FAP), but also perform a second analysis using the diffusive nested sampler K , which allows to perform model comparison and model selection in a fully Bayesian framework. K will be the method of choice for the remainder of the BEBOP survey Standing et al. (2021) . The parameters for the planetary companion are broadly compatible with those measured at the time of detection by Doyle et al. (2011) , and which have not been revised since. With our current precision it is not possible to determine the eccentricity of the orbit, but we can place an upper limit on it.

Finally we discuss the detectability of circumbinary systems. We only take the 20 first measurements we collected and run K measuring the Bayes Factor to assess the detectability of Kepler-16 b, where individual uncertainties are similar to the semi-amplitude of the signal, 1,pl . We repeat the procedure, measuring BF for each increase of five measurements until we reach a sub-sample containing the 50 first radial velocity measurements (they roughly cover four orbital periods of the planet). We reach a BF > 150 (the formal threshold for detection) with the first 40 measurements. We repeat this procedure but instead randomly select 50 measurements within the first three years of data. We run K and measure the Bayes Factor for the first 20 epochs of this sequence of 50 random epochs. We then increase from 20 to 25 until reaching 50. We find that the BF = 150 threshold is also passed at 40 measurements. We plot our results in Fig. 4 . The results between both series of tests are broadly consistent, except when we only use 25 measurements, with the Bayes Factor growing log-linearly with increasing number of measurements. Thanks to these tests we conclude that just 40 to 45 measurements would have been to formally detect Kepler-16 b with radial velocities. acknowledges support by CFisUC projects (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal.

The radial velocity data are included in the appendix of this paper. The reduced spectra are available at the SOPHIE archive: http://atlas.obs-hp.fr/sophie/. This paper has been typeset from a T E X/L A T E X file prepared by the author.

Stellar Radial Velocities

An Introduction to Close Binary Stars

Populations of Planets in Multiple Star Systems

Proceedings of the National Academy of Science

Ground-based and Airborne Instrumentation for Astronomy II

Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. p

The Rossiter-McLaughlin Effect in Exoplanet Research

Contemporary Physics

First, we would like to thank the staff, in particular the night assistants, at the Observatoire de Haute-Provence for their dedication and hard work, particularly during the COVID pandemic. The data were obtained first as part of a DDT graciously awarded by the OHP director (Prog. EW acknowledges support from the ERC Consolidator Grant funding scheme (grant agreement n • 772293/ASTROCHRONOMETRY). Support for this work was provided by NASA through the NASA Hubble Fellowship grant number HF2-51464 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. P.C. thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation, and the Moore Foundation; her participation in the program has benefited this work. SH acknowledges CNES funding through the grant 837319. ACC acknowledges support from the Science and Technology Facilities Council (STFC) consolidated grant number ST/R000824/1. ACMC