© ESO 2019

1 Introduction

The progression of the field of exoplanets has led to more and more diverse discoveries. A particular example is planets orbniting around both stars of a tight binary, known as circumbinary planets. These planets provide insight into planet formation in different, perturbative environments (Meschiari 2012; Paardekooper et al. 2012; Rafikov 2013; Lines et al. 2014) and of planetary migration in a protoplanetary disc that is sculpted by the inner binary (Artymowicz & Lubow 1994; Pierens & Nelson 2013, 2018; Martin et al. 2013; Kley & Haghighipour 2014). Studying the abundance of circumbinary planets ellicits comparisons to that around single stars (Martin & Triaud 2014; Armstrong et al. 2014; Bonavita et al. 2016; Klagyivik et al. 2017), and the presence or even absence of planets sheds light on the formation of their host binary (Muñoz & Lai 2015; Martin et al. 2015; Hamers et al. 2016; Xu & Lai 2016). Furthermore, circumbinary planets have increased transit probabilities (Borucki & Summers 1984; Schneider 1994; Martin & Triaud 2015; Li et al. 2016; Martin 2017), making them often found in the habitable zone (Kane & Hinkel 2013; Haghighipour & Kaltenegger 2013) and potentially aiding the use of transmission and emission spectroscopy (Deming & Seager 2017).

Roughly two dozen circumbinary planets have been discovered to date, with the majority coming from only two techniques: observing the planet transit in front of one or both of its host stars (e.g. Kepler 16; Doyle et al. 2011) or inferring its existence by the measurement of eclipse timing variations (ETVs) of the binary (e.g. NN Serpentis; Qian et al. 2009; Beuermann et al. 2010). Only a handful of discoveries have come from other methods: gravitational microlensing (OGLE-2007- BLG-349; Bennett et al. 2016), direct imaging (e.g. HD 106906; Bailey et al. 2014; Lagrange et al. 2016) and pulsar timing (PSR B1620-26; Backer et al. 1993; Thorsett et al. 1993; Sigurdsson et al. 2003). Furthermore, out of the two dominant techniques only the transit discoveries are completely reliable, as the validity of the ETV planets is debated, particularly for post-common envelope binaries (Zorotovic & Schreiber 2013; Bear & Soker 2014; Hardy et al. 2016; Bours et al. 2016).

Our knowledge of circumbinary planets is largely based on a sample of transiting planets that is both small in number, 11, and impacted by observational biases. Preliminary insights into their formation and distribution have been obtained (reviewed in Welsh & Orosz 2017; Martin 2018) but a more comprehensive understanding demands not only more discoveries, but ones made with complementary observing techniques with different sensitivities.

Radial velocities (RVs) led to the first discovered exoplanet around a Sun-like star (51 Peg; Mayor & Queloz 1995) and hundreds more since. Radial velocities have also yielded dozens of circumstellar planets orbiting one component of a wider binary, for instance 16 Cygni Bb (Cochran et al. 1997), and WASP-94Ab & Bb (Neveu-VanMalle et al. 2014). There is, however, no bonafide circumbinary planet discovered by RVs. This is in spite of attempts stretching back many years, for example the TATOOINE survey of double-lined spectroscopic binaries (SB2s; Konacki et al. 2009, 2010; Hełminiak et al. 2012). There was a potential RV discovery of a circumbinary planet in the HD 202206 system by Correia et al. (2005), but astrometric data has since revealed it to be a ${17.9}_{-1.8}^{+2.9}{M}_{\text{Jup}}$ circumbinary brown dwarf (Benedict & Harrison 2017).

Despite these past difficulties, RVs still have the ability to expand our knowledge of circumbinaries. The sensitivity of RVs partially overlaps with that of transits. This means that we may use transit discoveries as a guide and motivation, but radial velocities push into new parameter spaces, with a weaker dependence on period and inclination and also having detections that are mass-dependent, rather than radius-dependent.

In 2013 we intiated a radial velocity survey named Binaries Escorted By Orbiting Planets, henceforth (BEBOP). The programme was initiated on the 1.2 m Swiss Euler Telescope using the CORALIE spectrograph. We targeted 47 known single-lined eclipsing binaries (SB1s) drawn from the EBLM programme (Triaud et al. 2017a and see Sect. 3.1), which consist of F/G/K primaries and M-dwarf secondaries. The BEBOP observations reach a precision of a few metres per second. This permits a precise characterisation of the binary orbit, such that by subtracting it from the RV signal we may then search the residuals for an orders of magnitude smaller signal of a circumbinary gas giant planet.

A fundamental design element of BEBOP is to solely observe SB1s. This means that instead of trying to solve the problemof deconvolving the two spectra in SB2s, such as in the TATOOINE survey, our approach avoids it.

This paper is structured as follows. In Sect. 2 we briefly review the present understanding of circumbinary planets, anduse this to motivate a radial velocity survey. Second, in Sect. 3 the birth and construction of the BEBOP sample is described, including its roots in the EBLM survey. We then describe the observational strategy in Sect. 4. The subsequent treatment of the data, including the reduction of the spectroscopic data, fitting of radial velocity orbits and model selection is covered in Sect. 5. In Sect. 6 we discuss the calculation of primary andsecondary masses, as well as the constraints placed on undetected orbital parameters. In Sect. 7 we present the results of the survey and the characterisation of tertiary Keplerian signals. In Sect. 8 we analyse the results and compute detection limits for each of our systems. From these detection limits, in Sect. 9 we calculate the completeness of this programme and use it to calculate abundances of circumbinary planets, circumbinary brown dwarfs and tight triple star systems. We compare these values to other surveys for gas giants around single and binary stars in Sect. 10. Finally, in Sect. 11 we briefly outline the future of the BEBOP survey, including its recent upgrade to the HARPS spectrograph, before concluding in Sect. 12.

Fig. 1 In dark blue, a histogram of the Kepler eclipsing binary catalogue for all periods longer than 1 day. Data is taken from http:// keplerebs.villanova.edu/ as of May 2017, based on a catalogue that is first outlined in Prsa et al. (2011). Red dashed lines correspond to the nine binaries known to host transiting circumbinary planets. In light blue, a histogram of the 47 BEBOP binaries.

2 Trends seen in transiting circumbinary planets

With not even a dozen confirmed transiting circumbinary planets to date, we only have a cursory knowledge of their population. To provide context for the BEBOP survey we briefly review some of the trends. For a more in depth discussion of the circumbinary planets discovered to date, the reader is directed to the reviews in Welsh & Orosz (2017) and Martin (2018).

2.1 Binary orbital periods

All circumbinary planets found as part of the Kepler mission transit eclipsing binaries. Owing to geometry, the Kepler eclipsing binary catalogue is expectedly biased towards short periods, with a median orbital period of 2.7 days. Planets have only been found around relatively long-period binaries though, with P_bin > 7 days. This is shown in Fig. 1, with a distribution of the Kepler eclipsing binary catalogue and the binary periods known to host circumbinary planets. For comparison, we overlay the distribution of the BEBOP eclipsing binary catalogue, the construction of which we discuss in Sect. 3.2.

Tighter binaries are typically accompanied by a third star (Tokovinin et al. 2006). This third star is suspected to either disrupt circumbinary planet formation (Muñoz & Lai 2015; Martin et al. 2015; Hamers et al. 2016; Xu & Lai 2016) or bias the planet sample to long-period, misaligned orbits, both of which would have been missed by the Kepler transit survey. Alternate explanations for the dearth of planets around the tightest binaries, which do not invoke the presence of a third star, have also been proposed: tidal expansion of the binary orbit, causing the planet to become unstable (Fleming et al. 2018) and UV evaporation of exoplanet atmospheres, shrinking them to an undetectable size (Sanz-Forcada et al. 2014).

2.2 Planet orbital periods

Out of the nine known transiting circumbinary systems, eight contain a planet with a period less than ten times the binary period (plotted in Fig 2, left). The scaled minimum distance between orbits is defined as r_min = [a_p(1 − e_p)]∕[a_bin(1 + e_bin)]. Seven of the nine systems contain a planet with r_min between a narrow range of 2.5 and 3.0, and hence align diagonally in a plot of planet periapse against binary apoapse in Fig. 2, right. A planet cannot orbit too close to the binary lest its orbit becomes unstable through resonant overlap (Mudryk & Wu 2006), and these seven planets all orbit within a relative distance of 50% to the stability limit (Dvorak 1986; Holman & Wiegert 1999; Chavez et al. 2015; Quarles et al. 2018)¹. The over-density of planets near the stability limit is not believed to be an observational bias (Martin & Triaud 2014; Li et al. 2016), although improved statistics are needed to draw strong conclusions.

The stability limit coincides closely with where the protoplanetary disc would have been truncated by the inner binary (Artymowicz & Lubow 1994). Planet formation is believed to be hindered this close to the binary (Meschiari 2012; Paardekooper et al. 2012; Rafikov 2013; Lines et al. 2014). Instead, the favoured explanation for a heightened frequency of planets near the stability limit is formation in the farther regions of the disc, followed by an inwards migration and then parking near the disc edge (Pierens & Nelson 2013, 2018; Kley & Haghighipour 2014).

Fig. 2 Left panel: planet and binary orbital periods, with dashed lines of constant period ratio. Right panel: planet periapse distance and binary apoapse distance, with dashed lines of constant scaled minimum distance, rmin = [ap(1 − ep)]∕[abin(1 + ebin)]. Systems with a single detected planet are shown as blue squares, whereas the three-planet Kepler-47 system is shown as grey circles.

2.3 Orbital alignment

The known transiting circumbinary planets exist on orbits that are coplanar with the binary to within ~ 4°. Li et al. (2016) concludes that this is indicative of the true underlying distribution and not an observational bias. However, the statistics are presently poor, and Martin & Triaud (2014) demonstrated that highly misaligned systems (more than just a few degrees) have a sparse, hard-to-identify transit signature, and hence could remain hidden in the Kepler data.

A nearly coplanar distribution would be indicative of either a primordially flat environment, or a re-alignment over time of the circumbinary disc (Foucart & Lai 2013) or the planet itself (Correia et al. 2016). On the other hand, misalignment may be produced by disc-warping (Facchini et al. 2013; Lodato & Facchini 2013), planet-planet scattering (Smullen et al. 2016) and tertiary star interactions (Muñoz & Lai 2015; Martin et al. 2015; Hamers et al. 2016).

2.4 Planet size and abundance

Only circumbinary planets larger than 3 R_⊕ have been found to date. Some of the larger planets have measured masses from ETVs, but for most no ETVs are detectable and hence only an upper limit may be placed. The most massive measured mass is 1.52 M_Jup (Kepler-1647; Kostov et al. 2016), but the majority are Saturn mass or smaller. The lack of Earth and super-Earth circumbinary planets is however likely to be an observational bias, owing to the unique challenges of capturing shallow planetary transits with irregular transit depths, timing and durations (Armstrong et al. 2013).

For circumbinary gas giants the studies by Martin & Triaud (2014) and Armstrong et al. (2014) provided two important results. First, it was demonstrated that the true abundance was degenerate with the mutual inclination distribution; comparing a coplanar and highly misaligned population, to produce the same number of detections the misaligned population must have a higher planetary abundance as its transit detection rate would be smaller. Second, the minimum abundance of transiting circumbinary gas giants, corresponding to a near-coplanar distribution, was found to be surprisingly similar to that around single stars (Howard et al. 2010; Mayor et al. 2011; Santerne et al. 2016). On a similar note, the imaging survey of Bonavita et al. (2016) determined that the abundance of sub-stellar companions on wide orbits did not differ significantly between single and binary stars. Overall though, these results all require verification due to the presently poor statistics, one of the primary objectives of our survey.

2.5 Binary mass ratios

In Fig. 3 we show the mass ratios of the known transiting circumbinary planets (red dashed lines). For comparison, the blue histogram shows the host mass ratios in the BEBOP sample. Planets have been found around binaries of essentially all mass ratios. Whether the abundance of circumbinary planets depends on binary mass ratio remains an open question due primarily to small number statistics², but for now it is simply re-assuring that the planets known to date are just as commonly found around small mass ratio binaries, like those probed by BEBOP. Finally, we also remind that the circumbinary brown dwarf HD 202206c orbits a binary with q = 0.08 (Correia et al. 2005; Benedict & Harrison 2017).

Fig. 3 In light blue, a histogram of the mass ratio of the 47 BEBOP binaries. The vertical red dashed lines correspond to the mass ratios of the Kepler binaries known to host transiting circumbinary planets.

3 Overview of BEBOP

The targets for the BEBOP survey were first discovered and characterised as part of the older EBLM survey for low mass eclipsing binaries. It is for this reason that they are all designated by their EBLM name. This survey has been detailed in a series of papers (Triaud et al. 2013, 2017a; Gómez Maqueo Chew et al. 2014; von Boetticher et al. 2017), so in Sect. 3.1 it is just briefly reviewed. We then in Sect. 3.2 discuss how the BEBOP survey and its target list was constructed.

3.1 A brief description of the EBLM survey for low-mass eclipsing binaries

Since 2004 the Wide Angle Search for Planets (WASP; Pollacco et al. 2006; Collier Cameron et al. 2007a,b) has been conductinga ground-based photometric survey of several million stars. Observations span both hemispheres, with sites in La Palma and the South African Astronomical Observatory. The photometric precision and observing baseline are amenable to the detection of Jupiter-sized bodies on orbits of typically less than 10 days, although some are found with periods up to 40 days. Detecting hot- and warm-Jupiter planets is the primary objective of WASP.

However, the WASP survey has also netted a large quantity of astrophysical false-positives. An ever-present challenge is the ambiguity between the transit of a giant exoplanet and the eclipse of a low-mass star. Theoretical and observational studies have demonstrated that the vast range in mass between giant planets (~ 0.1 M_Jup) and small stars (~100 M_Jup), including the brown dwarfs in between, only corresponds to a narrow range in radius of ~ 0.7−2 R_Jup, in spite of the very different physical processes taking place (Baraffe et al. 1998, 2003, 2015; Chabrier et al. 2009; Chen & Kipping 2017). When using the humble precision of WASP (compared with Hubble, Kepler, etc.), it is therefore almost impossible to distinguish between a giant exoplanet and a small star using photometry alone; spectroscopic reconnaissance is required.

It is typical for most exoplanet surveys to discard candidates that display RV amplitudes in excess of a few km s⁻¹, and to consider those as false positives. However, this is not the case for the southern (Dec < +10°) WASP candidates, which are monitored spectroscopically by the Swiss Euler Telescope in La Silla, Chile, using the CORALIE spectrograph (Queloz et al. 2001a; Wilson et al. 2008; Triaud et al. 2017b). Systems with semi-amplitudes less than 50 km s⁻¹ enter the EBLM (Eclipsing Binary Low Mass) project. The project started in 2010 as an observational probe of eclipsing binaries with low-mass, M-dwarf secondary stars. The cut of 50 km s⁻¹ is designed to concentrate our observational efforts on fully convective secondaries (< 0.35 M_⊙), for which empirical mass and radius measurements are in short supply. In other words, the EBLM project is a survey of eclipsing SB1s.

An outline of the EBLM project, and some initial results were published in Triaud et al. (2013) and Gómez Maqueo Chew et al. (2014). The spectroscopic orbits of an ensemble of 118 binaries appeared in Triaud et al. (2017a), with this sample due to double in the coming year. The most recent result of the survey is the binary EBLM J0555-57, whose secondary star comes close to the hydrogen-burning limit with a mass of 85 ± 4 M_Jup while having a radius of only ${0.84}_{-0.04}^{+0.14}$ R_Jup, comparable to Saturn (von Boetticher et al. 2017).

3.2 The birth of the BEBOP survey for circumbinary planets

A spectroscopic exoplanet survey is complementary to the work already done using transits. This can be seen in the equation for the RV semi-amplitude, defined by $$K_{\text{c}}=\frac{{\left(2\pi G\right)}^{1/3}}{\sqrt{1-e_{\text{c}}^2}}\frac{m_{\text{c}}\sin I_{\text{c}}}{{\left(m_{\text{A}}+m_{\text{B}}+m_{\text{c}}\right)}^{2/3}}\frac{1}{P_{\text{c}}^{1/3}},$$(1)

where G is Newton’s gravitational constant, P denotes the orbital period, e is the eccentricity, I the inclination compared to the plane of the sky, and m the mass. The sub-scripts stand for the primary (A), the secondary (B), and the planet (c)³. We note that the mass of the secondary star appears in the above equation because the gravitational force of the planet perturbs the barycentre of the inner binary, rather than the primary star alone. Compared to transit surveys, RVs are sensitive to a wider range of planetary orbits, with a shallower dependence on the planet’s period and inclination.

Furthermore, RVs are sensitive to mass rather than radius. The majority of the Kepler circumbinary planets do not have masses measured from ETVs, and the faintness of most Kepler stars makes them unamenable to spectroscopic follow-up. There also remains some general tension in the community between masses derived photodynamically and spectroscopically (Steffen 2016; Rajpaul et al. 2017).

Past radial-velocity surveys have not yielded any confirmed circumbinary planets. The TATOOINE survey of non-eclipsing SB2s (Konacki et al. 2009, 2010; Hełminiak et al. 2012) was the most expansive effort. One of its major successes was an improvement in the precision of radial velocity measurements of double-lined binaries by at least an order of magnitude. Nevertheless, Konacki et al. (2009) reveal a mean rms across their sample of nearly 20 m s⁻¹, which exceeds the formal uncertainties by a factor of a few, and hides most gas giants from identification. We suspect the excess noise originates from an imperfect radial-velocity extraction, as the procedure gets affected by the presence of two sets of lines. Similar effects are seen in ELODIE data (Eggenberger et al. 2004).

Konacki et al. (2010) write that maximal precision can be achieved by monitoring “single stars, or at best single-lined spectroscopic binaries where the influence of the secondary spectrum can be neglected”. However, a suitable sample of bright, short-period SB1s was not available when TATOOINE was first constructed.

The BEBOP survey started when we realised that such a binary sample did now exist: the EBLM survey. By construction, the EBLM sample is solely composed of SB1s. Indeed, thanks to their eclipsing configuration we calculate the true (not minimum) mass of the secondary and its radius. Together we can robustly estimate the level of contamination produced by the secondary (Triaud et al. 2017a). Instead of attempting to build upon the pioneering work of TATOOINE to solve the double-line binary problem, we decided to circumvent it by focusing on single-line binaries.

By avoiding the contaminating effect of a secondary set of lines, the identification of a circumbinary planet becomes equivalent to identifying a multi-planet system whose innermost object happens to have a few 100 M_Jup. We note that hot-Jupiters have dayside temperatures ranging from ~800 to 4600 K (Triaud et al. 2015; Gaudi et al. 2017), and consequently M-dwarfs and hot-Jupiters are similarly located on colour-magnitude diagrams (Triaud 2014; Triaud et al. 2014). Surveys for outer companions to hot-Jupiters are common in the literature (e.g. Knutson et al. 2014; Bryan et al. 2016; Neveu-VanMalle et al. 2016), and the BEBOP survey is conceptually similar.

In addition to being SB1s, the EBLM targets also have the following beneficial attributes:

• Consistency of the sample: the EBLM targets were all discovered and characterised using only WASP photometry and CORALIE spectroscopy, with a consistent set of procedures and sensitivities.

• Past EBLM RVs were available to be combined with new measurements taken for BEBOP, which roughly doubles our time baseline, and therefore improves our sensitivity to long-period outer companions.

• Some EBLM targets already had identified stellar activity, and hence could be avoided.

• The radial velocity amplitude of the planet (Eq. (1)) is a decreasing function of the sum of the primary andsecondary masses. Having a low-mass secondary star is therefore beneficial.

• All of our binaries eclipse, which biases the orbital orientation of any planets to maximise the RV signal.

• Another advantage of eclipsing binaries is a positive bias of the transit probability of any discovered circumbinary planet (Borucki & Summers 1984; Schneider 1994; Martin & Triaud 2014, 2015; Li et al. 2016; Martin 2017). This bias is particularly strong for small mass ratio binaries, which must have inclinations very close to 90° (further investigation in Martin, in prep.).

• The distribution of EBLM binary periods and mass ratios has significant overlap with the Kepler binaries known to host circumbinary planets (shown in Fig. 1).

• The BEBOP binaries have an average V_mag = 11, which is roughly 3.3 magnitudes brighter than the Kepler circumbinary systems.

3.3 Sample construction

The BEBOP binaries are selected from the EBLM sample, with the following protocol:

• The BEBOP binaries comply with a difference of four visual magnitudes between the primary and secondary stars, such that we avoid secondary contamination of the primary’s spectrum⁴. Almost all of the EBLM binaries naturally fulfill this criterion.

• We only keep binaries on whose primary we reach a precision of 70 m s⁻¹ or better during a 30 min observation, which is the typical exposure time used for the WASP planet survey (Triaud 2011). This is sufficient to reach the planetary domain. For instance, a hypothetical 3 M_Jup planet at P_c = 50 days around a m_A + m_B = 1.2 M_⊙ binary produces a detectable radial velocity signal with semi-amplitude K_c = 146 m s⁻¹.

• We exclude systems that display signs of stellar activity, as seen in an abnormal variation of the span of the bisector (Queloz et al. 2001b; Figueira et al. 2013). While stellar activity does not prevent the large-amplitude binary orbit to be characterised, it becomes a hindrance for detecting small-amplitude planets, sometimes mimicking their signal. The identification of stellar activity in the EBLM binaries is outlined in Triaud et al. (2017a). Some of the EBLM binaries were already known to exist inside a triple star system. Outer stellar companions are thought to truncate and shorten the lifetime of the protoplanetary disc (Kraus et al. 2012; Daemgen et al. 2013; Cheetham et al. 2015), and generally be detrimental for the formation and survival of circumbinary planets (Muñoz & Lai 2015; Martin et al. 2015; Hamers et al. 2016; Xu & Lai 2016). However, such triple systems are kept in our sample for two reasons. First, the searches for circumbinary planets around the Kepler eclipsing binaries were done so without any a priori knowledge of a tertiary companion. Indeed, one example is known of Kepler-64 (Schwamb et al. 2013; Kostov et al. 2013)⁵ which has an outer stellar companion, which is itself a binary, albeit at a large separation of ~ 1000 AU. The second reason to keep triple star systems is to avoid introducing a confirmation bias into our survey.

The BEBOP binaries are typically longer period than the EBLM sample from which they were chosen. This was not chosen to match the trend seen in the Kepler results that the tightest binaries do not host planets (Fig. 1), as this would also introduce a confirmation bias. Instead, this long-period selection is a function of the obtainable RV precision. Our binaries are all (or close to it) tidally synchronised (or pseudo-synchronised if eccentric), and hence the rotation period equals the orbital period. Consequently, the tightest binaries are also the fastest rotators, which have the worst RV precision due to broadened spectral lines. An example of this can be seen by comparing EBLM J1146-42 and EBLM J1525+03. The two targets have a similar visual magnitude (V_mag = 10.29 and 10.74) and primary mass (M_A = 1.35 M_⊙ and 1.23 M_⊙). However, the RV precision is significantly different (σ_1800s = 9 and 48 m s⁻¹), which we attribute to different orbital periods (P_bin = 10.5 and 3.82 days). In fact, this 3.82-day period for EBLM J1525+03 is the shortest in the sample, and also corresponds to one of our worst precisions.

The BEBOP sample tallies 47 binaries, which we present in this paper. In Table A.1 we list some of the fundamentalobservational and physical parameters of these binaries. The calculation of the primary and secondary masses is discussed in Sect. 5.4. The primary visual magnitudes are all taken for the NOMAD survey, except for EBLM J1934-42 which did not have available data. For this exceptional case the Baraffe et al. (2015) models were used at an age of 1 Gyr. For the secondary visual magnitudes Baraffe et al. (2015) models were used in all cases. The mid-times of primary and secondary eclipse (T_pri and T_sec), are calculated based on the CORALIE spectroscopy alone and not the WASP photometry. The different σ values are the observational precisions, and are discussed in Sect. 5.3.

4 Observational strategy

All spectroscopic observations were taken at the Swiss Euler Telescope in La Silla, Chile, using the CORALIE spectrograph. CORALIE (Queloz et al. 2001a; Marmier et al. 2013) is a thermally stabilised, fibre-fed echelle spectrograph with a resolving power of R = 55 000.

The goal was to collect 20 observations of 30 min length on each binary. The flexible observing schedule of the Swiss Telescope allows for observations to be spread out over the year. This is important for probing long-period planets, like the ones we expect to find. The Kepler-discovered circumbinary planets have periods between 49.5 and 1107 days, with a median of 184 days.

Observations were instructed to be separated by at least half the binary’s orbital period. This means that the 20 observations would span at least 10P_bin, which would be long enough to cover at least one orbital revolution of a planet in eight of the nine Kepler circumbinary systems. Typically though, the observations were spaced over a longer timespan.

Other constraints on the observations were as follows:

• Separation between the target star and the Moon by at least 70°. This conservative criterion avoids contamination of the spectrum by the gentle Sun’s light reflected off the delicate Lunar surface.

• Avoidance of primary eclipses of the binary. When the secondary M dwarf passes in front of the primary star the radial velocity signal is slightly distorted by the Rossiter–McLaughlin effect (Holt 1893; Rossiter 1924; McLaughlin 1924; Queloz et al. 2000; Triaud et al. 2013). We did not instruct observers to avoid secondary eclipses of the binary, as the faintness of the secondary stars makes these phenomena negligible.

• Generally good, clear observing conditions were required. This meant an airmass of the target better than 1.5 and a seeing better than 2.0 arcsec.

Since the BEBOP sample was constructed from the existing EBLM programme, we included all available radial velocity data except those likely affected by the Rossiter–McLaughlin effect. Measurements from the EBLM programme were also removed if deemed outliers, which is explained in Sect. 5.2.

In Fig. 4 we show the calendar of observations on the 47 eclipsing binaries. The red diamonds correspond to long observations (1700+ s). These typically correspond to the BEBOP programme since late 2013, and a couple of initial observations dating back as far as 2008. The blue squares are for shorter observations, which in most cases were taken under the guise of the EBLM programme. A small number of BEBOP targets did not receive a full quota of 20 long observations, owing to limitations in available observing time, but most exceeded this.

5 Radial velocity data treatment

The radial velocity data was treated in the same way as in Triaud et al. (2017a). We therefore only provide a summary of the methods used here, and refer the reader to that paper for a more thorough discussion.

5.1 Reduction of spectroscopic data

The CORALIE Data Reduction Software (DRS) is similar to that used with the HARPS, HARPS-North and SOPHIE instruments.A cross correlation function (CCF) is created between the observed spectrum and a numerical mask (Baranne et al. 1996). The CCF is a weighted average spectral line, which contains characteristics of individual absorption lines such as their width and asymmetries, but with a heightened signal to noise. The CCF is binned in 0.5 km s⁻¹ increments, owing to theR = 55 000 resolving power of the spectrograph. Two different spectral type masks were used: G2 and K5. These were chosen based on the spectral type of the primary star, which in our sample ranges between K2 and F0. Dumusque et al. (2012) demonstrated that the closeness of the spectral type mask largely affects only the absolute radial velocity and not the radial velocity variations. Only having two spectral type masks therefore does not hinder our analysis.

The CCF of each measurement was compared with a Thorium-Argon spectrum, which was used as a wavelength-calibration reference (Lovis & Pepe 2007). This accounted for variations in the instrument which would otherwise impose a drift on the measured radial velocity of the star. The main source of variation was the pressure of the spectrograph, which followed the ambient atmospheric pressure because CORALIE, unlike HARPS and HARPS-North, is not pressurised. In 2014 a Fabry-Pérot unit was added to provide even more precise calibrations.

5.2 Outlier removal

For each CCF, which is approximately Gaussian, we measure the span of the bisector slope. The bisector is calculated by tracing vertically the midpoint of the CCF at each value of flux intensity. The span of the bisector slope is the difference between the bisector at the top and bottom of the CCF (Queloz et al. 2001b). The bisector therefore reflects any asymmetries in the absorption lines. We remove any observations with bisector positions more than three interquartile ranges below the first quartile or above the third quartile. Such outliers may be from the wrong star being observed accidentally, or an abnormally low signal to noise observation. A visual inspection was also done of the data to outliers in the CCF’s Full Width Half Maximum (FWHM).

5.3 Calculation of radial velocity uncertainties

In Table A.1 three different σ radial velocity uncertainties are listed. The first value, σ_1800s, is mean the photon noise uncertainty for all observations of 1800 s, which was the typical observation length during the BEBOP programme. The value σ_median is the median photon noise precision for all observations, that is both the 1800 s observations taken for BEBOP and earlier, shorter observations taken in the EBLM programme.

After removing observations with significant outlier bisector values there may still remain some variation in the bisector. We consider such asymmetries to be a source of error, which is shown in Table A.1 as σ_add. This value is calculated by $${\sigma }_{\text{add}}=\sqrt{\frac{{\delta }_{\text{bis}}^2}{4}}-〈{\sigma }_{\gamma }^2〉,$$(2)

where δ_bis is the root mean square of the variation of the bisector measurements around their mean and $〈{\sigma }_{\gamma }〉$ is the mean photon noise error. If $〈{\sigma }_{\gamma }〉>{\delta }_{\text{bis}}^2/4$ then we take σ_add= 0, which was the case for 30 of the 47 binaries. Otherwise, we add σ_add in quadrature to the radial velocity measurements.

Finally, the CORALIE spectrograph was historically stable to a precision of 6 m s⁻¹ (Marmier 2014). A recent change of the optical fibres from circular to hexagonal improved the stability to 3 m s⁻¹ (Triaud et al. 2017b). To each data point we add 6 m s⁻¹ of Gaussian noise in quadrature. Choosing 6 m s⁻¹ and not 3 m s⁻¹ was considered conservative, and also reasonable since most of the observations occurred before the fibre upgrade.

Fig. 4 Time-series of 1519 radial velocity observations of the 47 eclipsing binaries in the BEBOP programme. Red diamonds are observations for 1700 s or longer. Blue squares are for shorter observations. It is seen that all binaries typically receive two long observations initially before a series of short observations as part of the EBLM programme. A series of long observations typically commenced near the end of 2013 as part of the BEBOP survey.

5.4 Orbit fitting

To fit orbits to the spectroscopic data we use the YORBIT genetic algorithm, which has been developed over the years at the University of Geneva and implemented in numerous radial velocities using CORALIE and HARPS (e.g. Bonfils et al. 2013; Mayor et al. 2011; Marmier et al. 2013). Only static Keplerian orbits are fitted, that is orbital variations induced by gravitational interactions between orbits are ignored. This is a reasonable assumption except for very tight triple star systems, and in Sect. 8.4 we briefly discuss one such example. More details on YORBIT may be found in Bouchy et al. (2016).

When YORBIT is run to search for a single Keplerian orbit it will inevitably first fit that of the inner binary, as its signal is orders of magnitude higher than any potential circumbinary orbit. This binary orbit is characterised by six parameters: the period, P, semi-amplitude, K, eccentricity, e, time of periapsis passage, T₀, mass function f(m) and the argument of periapsis, ω. Error bars are calculated for each of these parameters by running 5000 Monte Carlo simulations.

5.5 Model selection

For each target we fitted five different types of model to the spectroscopic data. These are listed below, along with the number of parameters shown in parenthesis.

• 1.

  k1: a single Keplerian (6)

• 2.

  k1d1: a single Keplerian plus a linear drift (7)

• 3.

  k1d2: a single Keplerian plus a quadratic drift (8)

• 4.

  k1d3: a single Keplerian plus a cubic drift (9)

• 5.

  k2: a pair of Keplerians (12)

Models more complex than a single Keplerian are likely indicative of a tertiary companion. This tertiary companion will have its own Keplerian orbit, but if the observational timespan only covers a small fraction of this outer period then the orbit will be sufficiently modelled by a drift. A drift could alternatively be explained by an instrumental variation, but for CORALIE the temperature stabilization and nightly pressure calibrations have historically avoided this. A third explanation would be long-term stellar activity, although the binaries were all vetted for heightened activity, as described in Triaud et al. (2017a) based on the proceedures of Queloz et al. (2001b) and Figueira et al. (2013).

When attempting to fit a two-Keplerian model the genetic algorithm was restricted to searching for periods greater than four times the inner binary period. Numerous stability studies (e.g. Dvorak 1986; Holman & Wiegert 1999; Chavez et al. 2015; Quarles et al. 2018) show that circumbinary planets would be unstable with shorter orbits. Aside from this minimum period, no further restrictions are applied to the fitting. In particular, we search for and are sensitive to binaries and planets of all eccentricities. Upon the discovery of any candidate triple systems the orbital stability is then tested more carefully.

All targets in the BEBOP programme have been observed at least 16 times, with a median count of 32, which is more than there are free parameters in any of the models.

YORBITwill always retain small eccentricities like most fitting procedures (Lucy & Sweeney 1971). Therefore, for each of the above types of model we also test a fit where eccentricity is forced to zero. This allows us to test if the eccentricity fitted by YORBIT is significant. We are therefore left with a total of ten tested models, where the number of parameters for the forced circular model is always two less than the corresponding eccentric model, as both e and ω are removed. Throughout this paper we use “(circ)” and “(ecc)” to distinguish between forced circular and freely eccentric models. Note that when testing the k2 (circ) model only the binary eccentricity is forced to zero, not that of the planet.

For all ten models YORBIT outputs a χ² statistic, which is a weighted sum of the square of the residuals. The reduced χ² statistic is calculated by normalising over the number of free parameters: $${\chi }_{\text{red}}^2=\frac{{\chi }^2}{n_{\text{obs}}-k},$$(3)

where n_obs is the number of spectroscopic observations and k is the number of model parameters. A value of ${\chi }_{\text{red}}^2=1$ is indicative of an optimal fit.

To choose the most appropriate model between the ten possibilities we follow the same procedure as in Triaud et al. (2017a). For this we calculate the Bayesian information criterion (Schwarz 1978; Kass & Raftery 1995), henceforth referred to as the BIC, according to $$\text{BIC}={\chi }^2+k\ln \left(n_{\text{obs}}\right).$$(4)

The BIC is defined such that it naturally punishes complex models (large k), and hence has an inbuilt Ockham’s razor, in selecting the most parsimonious explanation possible.

In our process of model selection we start with the simplest model, which is k1 (circ) with only four parameters, and calculate the BIC. We then calculate the BIC for other models in order of complexity. To choose the next most complex model we demand that the BIC improves (decreases) by at least six. This is believed to be strong evidence for the more complex model (Kass & Raftery 1995).

We allow the model selection to “jump” levels of complexity if the BIC improves by n × 6, where n is the number of ranks of complexity moved through. For example, if BIC = 40 for k1 (circ) with four parameters, BIC = 38 for k1d1 (circ) with five parameters and BIC = 25 for k1d2 (circ) with six parameters, then the chosen model would be k1d2 (circ), even though there was only a marginal BIC improvement between k1 (circ) and k1d1 (circ).

This process is done for k1, k1d1 and k1d2 models, in both circular and eccentric flavours. These were deemed the “base” models. The k1d3 and k2 models were deemed “complex” models. Complex models are only tested if ${\chi }_{\text{red}}^2>2$ for the best base model. This criterion was an additional means of penalising overly complex models. Since the aim of the survey is to find circumbinary planets, that is k2 models, this cautious approach minimises false discoveries. We note that sometimes complex models were tested but ultimately a base model was chosen.

6 Calculating physical parameters

6.1 Primary and secondary masses

Because the BEBOP sample only contains single-lined binaries the primary and secondary masses are not directly measured, but rather only the mass function is directly measured. Primary masses are calculated the same way as in Triaud et al. (2017a), based on photometric colour fitting methodology outlined in Maxted et al. (2014).

Knowing both the primary mass and the mass function, the secondary mass is calculated by solving numerically $$f(m)=\frac{{\left(m_{\text{B}}\sin I_{\text{bin}}\right)}^3}{{\left(m_{\text{A}}+m_{\text{B}}\right)}^2}=\frac{P_{\text{bin}}{K}_{\text{pri}}^3}{2\pi G}.$$(5)

The error in m_B is calculated as $$\frac{\delta m_{\text{B}}}{m_{\text{B}}}=\frac{1}{3}\left(\frac{\delta f(m)}{f(m)}+2\frac{\delta m_{\text{A}}}{m_{\text{A}}}+3\frac{\delta \sin I_{\text{bin}}}{\sin I_{\text{bin}}}\right),$$(6)

where δ indicates a 1σ uncertainty and δf(m) is a direct output from YORBIT and δm_A is calculatedbased on Torres et al. (2010) and Maxted et al. (2014). The uncertainty in the binary inclination is calculated as δ sin I_bin = R_A∕a_pri, where a_pri is the semi-major axis of the binary and R_A is the radius of the primary star, as calculated based on Gray (2008). When m_B is calculated in Eq. (5) we take I_bin = 90° because of the existence of the binary eclipses. By adding this small inclination uncertainty we reflect our ignorance of the impact parameter of the eclipse. Thistypically adds 20% or less to the error in m_B. Based on the primary and secondary masses the semi-major axis is calculated using Kepler’s third law.

6.2 Calculating upper limits on undetected orbital parameters

We follow the same methodology as for the EBLM survey (Triaud et al. 2017a) to constrain upper limits on orbital parameters which we do not have the precision to directly measure. For all binaries where a forced circular solution was chosen by the BIC model selection we constrain the eccentricity to within zero and an upper limit. This upper limit is calculated by adding the fitted eccentricity for that model to the 1σ uncertainty on that eccentricity. The same procedure is done for drifts in the radial velocity, for example if k1 (ecc) was the chosen model then the upper limit on the coefficient of linear drift was taken as the fitted value in k1d1 (ecc) plus its 1σ uncertainty.

7 Spectroscopic results

7.1 Chosen models

Table A.2 shows all of the information pertaining to the model selection. The BIC for the selected model is highlighted in bold font. In Table 1 we count how many binaries were fitted by each of the ten possible models. The same table appears in Triaud et al. (2017a, Table 2 in that paper) for the EBLM survey. The most noticeable difference is that here we report 6 binaries fitted with k1 (circ) and 25 with k1 (ecc). Contrastingly, for the EBLM paper 58 binaries were fitted with k1 (circ) and 39 with k1 (ecc). The heightened percentage of binaries fitted with k1 (ecc) in the BEBOP survey reflects the additional long-exposure radial velocity measurements, which heighten our sensitivity to eccentricities as small as 0.001.

7.2 Innerbinary parameters

The orbital parameters and masses for the inner binary are all listed in Table A.3, all taken from the chosen model indicated in Table A.2. The 1σ uncertainty for each parameter is given in parenthesis, corresponding to the last two digits of the measured value. For example, for EBLM J0008+02 the period is P_bin= 4.7222824(48) days, which can be otherwise written as P_bin = 4.7222824 ± 0.0000048 days. For some systems upper limits are provided for the eccentricity and coefficients of drift according to Sect. 6.2. The quantity ω_bin is undefined when a forced circular model is chosen.

7.3 Discovered or potential tertiary bodies

For five of our binaries the selected model is a pair of Keplerian orbits. Unfortunately from a planetary perspective, all of the characterised tertiary orbits are within the stellar regime. The smallest characterised tertiary mass is m_c sin I_c = 0.1207 M_⊙ for EBLM J2011-71.

In Table A.4 we provide parameters for all five characterised triple star systems. Compared with the EBLM release in Triaud et al. (2017a), there is an additional system: EBLM J1038-37. This system was included in Triaud et al. (2017a) but the BIC selection criteria characterised it as a single, eccentric binary plus a cubic drift. That prior characterisation was based on 13 observations taken over 3.89 yr, with a median precision of 131 m s⁻¹. The double Keplerian characterisation presented in this paper is based on 33 observations taken over 5.01 yr, with a median precision of 74 m s⁻¹. This is an example of the improved orbital fits provided by the BEBOP survey in comparison with the original EBLM survey.

It is interesting to note that all of the minimum masses of the tertiary stars are all significantly smaller than the primary masses. However, we caution drawing too much from this result as a more massive tertiary would have dilutedthe already shallow WASP eclipse depth of the secondary star, and hence such systems may not have been detected in the first place. None of the tertiary stars are bright enough for us to directly observe their flux.

8 Analysis

8.1 Residuals as a function of orbital phase – a test for spectral contamination

The BEBOP binary sample was constructed to avoid cross contamination between two sets of stellar spectral lines by only choosing binaries with faint secondary stars. However, if there were spectral contamination then it would only be expected to affect the observations at certain binary orbital phases. The vulnerable orbital phases correspond to the primary star’s radial velocity equalling the system’s systemic velocity, as it will be also equal to the secondary star’s radial velocity (which we do not directly measure). At this point the primary and secondary spectral lines overlap, and hence the chance of contamination is maximised. It was said in Sect. 4 that observations were taken at orbital phases that avoid the primary eclipse, such that we do not observe a Rossiter–McLaughlin effect which would skew the radial velocities. This fortuitously helps us avoid spectral contamination, as the eclipse corresponds to an overlapping of the primary and secondary radial velocity signals.

The width of the spectral lines, quantified by the full width half maximum (FWHM), demarcates the range of radial velocity values that may be contaminated. A larger FWHM means that spectral contamination may occur at a greater difference in the primary and secondary radial velocities.

To test if there is wide-spread spectral contamination in our sample we analyse the residuals of the best-fitting model to each binary as a function of a scaled radial velocity value, for all 47 binaries. This scaled radial velocity value, which we denote by δ, is calculated by $$\delta =\frac{\text{RV}-{\text{RV}}_0}{FWHM}\frac{m_{\text{A}}}{m_{\text{B}}},$$(7)

where RV is an individual radial velocity measurement on the primary star and RV₀ is the systemic radial velocity mid-point for the system. The mass ratio factor m_A ∕m_B converts the radial velocities from the measured values on the primary star to the larger but not directly measured values on the secondary star. The δ value denotes how far the potentially contaminant secondary star radial velocity signal is away from RV₀.

The results are plotted in Fig. 5. In blue we plot the residuals (O–C) of all 1519 data points for the 47 BEBOP targets on top of each other as a function of δ. The two dark blue lines are the root mean squared (RMS) values at ten different values of δ, calculated separately for positive and negative values of O–C. If our sample were affected by stellar contamination then there would be a significant increase in the residuals near δ = 0, and hence the red RMS curves would deviate significantly from zero. This is not the case.

It was noted in Sect. 3.2 in the construction of the BEBOP sample that one target – EBLM J0425-46 – has a visual magnitude difference of 3.85 between the primary and secondary stars, whichis slightly smaller than the threshold of 4. Its mass ratio is the highest of the sample: m_B ∕m_A = 0.527. This may put it at risk of spectral contamination from the relatively “bright” secondary star.

In Fig. 6 we plot the 30 radial velocity measurements over the and their residuals to the fitted eccentric k1 model. This is the model that was chosen as the most appropriate by the BIC selection method, with ${\chi }_{\text{red}}^2=0.73$. We highlight with blue boxes in this plot two very marginal outliers of the fit, showing that they correspond to the two radial velocity measurements closest to the systemic velocity (0 m s⁻¹ on this plot). Having outliers here is consistent with spectral contamination. However, each value is less than two standard deviations away from the fit, and hence not statistically significant, and the fit is overall very good to the data (${\chi }_{\text{red}}^2<1$).

Fig. 5 Radial velocity residuals (“Observed minus Calculated” or O–C for short) of the best-fitting model to all BEBOP binaries forall 1519 observations, stacked on top of each other as a function of the scaled radial velocity value δ, defined according to Eq. (7). The δ value says how far the observed radial velocity observation was from the systemic velocity, scaled by the average FWHM for that target. The two roughly horizontal dark blue lines are the root mean squared (RMS) values of the residuals, split into positive and negative values. If the RMS were maximised near δ = 0 then this would be indicative of wide-spread spectral contamination from the secondary star, but this is evidently not the case.

Fig. 6 Top panel: radial velocities of the target EBLM J0425-46 over 3.15 yr and the selected eccentric single Keplerian model. Error bars are typically 10–20 m s−1 and too small to see at this scale. Bottom panel: residuals to the fitted model, with 1σ error bars. Blue boxes highlight two marginal outliers, with a dashed line connecting the radial velocity measurement, which is near0 m s−1, and the residual. Both outliers are less than 2σ from the model.

8.2 Calculating detection limits

There are two main factors that determine the detectability of a putative planet: its minimum mass, m_c sin I_c, and its period, P_c. Not only are these the main contributors to the amplitude of the radial velocity signal, but the observational timespan needs to cover a significant portion of the planetary period. The eccentricity of a planet also increases K_c, but Endl et al. (2002) demonstrate that it has a minimum effect on detectability for e_c ≲ 0.5. Therefore, whilst our search for circumbinary planets is sensitive to any eccentricity, when quantifying the detectability we only consider circular planetary orbits.

To calculate detection limits for each binary we introduce and attempt to retrieve artificial Doppler signals. We follow a similarprocedure to that described in Konacki et al. (2009), which is based on methods that are regularly employed to calculate the occurrence rates of planets by long-term Doppler surveys on single stars (e.g. Cumming et al. 2008; Mayor et al. 2011; Bonfils et al. 2013).

We first start by defining what makes a hypothetical planet detectable. For this, we use the Generalised Lomb Scargle (GLS) periodogram, which identifies periodic signals of varying strengths within data. We define a putative planet as “detectable” when a GLS periodogram displays a signal with a strength rising above a False Alarm Probability (FAP) of 1%.

Importantly, the injection of the synthetic Keplerian sinusoid must be done to data that has already been cleansed of any existing periodic signals. To do so, the main signal which we remove is that of the binary, which is multiple orders of magnitude larger than that from any planet. Additionally, in some systems we have evidence for an outer stellar companion. These additional signals can add power to the GLS periodogram and need to be removed as well. Therefore, for calculating detection limits we use the residuals to the best-fitting models, as determined in Sect. 7. On these residuals no periodic signal with a period longer than 4P_bin (the rough stability limit; Holman & Wiegert 1999) was discovered above a 1% FAP.

On the cleansed data for each of the 47 binaries we insert and retrieve artificial circumbinary signals in the following way. We create a grid of planet periods that is uniform in log n_c, where n_c = 2π∕P_c is the orbital frequency of the planet. The minimum period tested is 4P_bin, as shorter period planets would be unstable. The maximum period tested is equal to 4ΔT, where ΔT is the observing timespan of the observations for a given binary. We note that for EBLM J2046-40 the outer triple star has a period of 5557 days, which was discovered using a timespan of 1801 days, which is roughly a third in length.

For each period we insert Keplerian sinusoids with increasing radial velocity semi-amplitudes, K_c. Following that, we attempt to retrieve the artificial signal using the periodogram. The value of K_c is directly correlated to m_c sin I_c, and this takes into account our calculated values of m_A and m_B. As soon as the injected sinusoid produces a signal above a 1% FAP we define this as the minimum detectable mass for the binary at that period. For each period this process is repeated for 20 different planetary orbital phases, equally spaced between 0° and 360°, since some phases may be better illuminated by the observations than others.

In Fig. 7 we show example detection limit curves for one of the most precise BEBOP targets: EBLM J2011-71. This system is known to contain a tertiary M dwarf star, which is demarcated on the plot well above the detection limit curve. In Appendix B of this paper detection limit curves are provided for all BEBOP targets.

8.3 Genetic algorithm detection of n-body simulated radial velocity signals

Tests are run to verify that this periodogram-based definition of detectability matches our ability to detect planets with the YORBIT genetic algorithm. For two targets, J0035-69 and J0310-31, we construct a coarse grid of circumbinary planet periods and minimum masses. The grids are chosen to straddle either side of the detection limit curves which were calculated in Sect. 8.2. All other planetary parameters are set to zero, except the inclination which is taken at 90°. The binary parameters are the measured values. The reason for choosing these two targets in particular is that they have different binary parameters (including the mass ratio, period and eccentricity) and also precisions at either extreme of our programme (5–6 m s⁻¹ for J0310-31 and 50–60 m s⁻¹ for J0035-69).

At each grid point an n-body code⁶ is run to simulate the radial velocity signal of the hypothetical circumbinary system, including both the large-amplitude binary signal and the much smaller planetary signal. The radial velocity measurements are simulated at the same epochs as the actual observationswere for each target, and are given the same uncertainty. Importantly, the n-body simulation does not assume static Keplerians, and hence any dynamical perturbations by the binary on the planet’s orbit are naturally included. Contrastingly, the periodogram analysis assumes static orbits.

The YORBIT code is then run on the simulated radial velocities to search for a two-Keplerian solution. Owing to its large amplitude, the binary signal is always found easily. A second signal will always be fitted but we only consider the detection of the simulated planet to besuccessful if the YORBIT-found orbit has a period within 10% of the n-body simulated planet period. The threshold of 10% is admittedly somewhat arbitrary, but felt to be sufficient for this simple demonstration. A more thorough study of the effects of n-body interactions on RV detectability is beyond the scope of this paper.

In Fig. 8 we show the results of these tests. There is a close connection between the YORBIT detectability and the periodogram detectability. Recall that the range in the periodogram detectability is a result of testing 20 different orbital phases, whereas for YORBIT only a single phase is tested. There are a few exceptional cases where the YORBIT detectability is not a monotonic function of m_c sin I_c. There are two explanations for this. First, there is an element of randomness in any genetic algorithm. Second, the observational errors are redrawn from a normal distribution for each n-body test, and hence will randomly impact some simulations more than others.

Based on these tests, we conclude that the periodogram means of determining detectability sufficiently replicates how we actually detect circumbinary planets using YORBIT. Non-Keplerian effects, even near the stability limit for one of our most precise targets, are seemingly a negligible hindrance on detectability.

Fig. 7 Example of the detection limits calculated for EBLM J2011-71, which is one of the BEBOP targets with the highest precision. Each line represents the smallest detectable minimum circumbinary planet mass (mc sin Ic) as a functionof the planet’s period (Pc) for a different orbital phase of the planet. There are 20 orbital phases tested, all equally spaced. The blue square near the top of the plot is the detected circumbinary object in the EBLM J2011-71, which is in fact not a circumbinary planet but rather alow-mass tertiary star.

8.4 Evidence for n-body interactions

Only in one of our targets do we see likely evidence of n-body interactions: EBLM J1146-42. As seen in Table A.2 and in the orbit plots in Appendix B, there are significant residuals to the double Keplerian orbital fit: a scatter of ~ ±200 m s⁻¹ and ${\chi }_{\text{red}}^2=77.96$. In Fig. 9a we show the periodogram of the residuals to the k2 fit, which demonstrates a lack of any significant periodicitiesin the data. Indeed, attempts were made to fit additional orbits and drift parameters but none resulted in a significant improvement to the fit.

Stellar activity cannot produce residuals of that magnitude⁷, and this target shows no signs of such activity (Triaud et al. 2017a). In Fig. 9b we show a constant bisector of the radial velocity measurements. Figure 9c plots the residuals to the k2 fit as a function of the bisector. The classic indicator of activity is a negative correlation on this plot (Queloz et al. 2001b; Figueira et al. 2013), which is not apparent here. Furthermore, our spectra show no sign of contamination. This is expected since there is a 5.14 magnitude difference between the primary and secondary star flux. The tertiary body is even less massive than the secondary unless its orbital plane is misaligned by more than 55°.

Alternatively, the large residuals may be products of n-body interactions between the inner and outer orbits. Such interactions are not accounted for in the YORBIT-determined orbits, which are assumed to be static Keplerians. It is a future study to analyse such interactions in a means similar to Correia et al. (2010). This will hopefully yield a direct measurement of additional parameters such as the mutual inclination between the inner and outer orbits.

Fig. 8 Detection limit curves in black for two targets in the BEBOP programme: J0035-69(left panel) and J0310-31 (right panel). Each curve represents the smallest mc sin Ic for a putative circumbinary planet at a given period, according to the periodogram measure of detectability described in Sect. 8.2. For each target there are 20 black curves, one for each of the tested orbital phases of the injected planet. At seven discrete periods a series of tests were run to recover n-body simulated circumbinary planets with the YORBIT genetic algorithm. At these periods a green circle indicates a successful recovery, whereas a red diamond indicates a failure.

9 Abundance of circumbinary objects

9.1 Calculating the completeness of the programme

We define the completeness of our programme as a function of planet period and minimum mass, C(P_c, m_c sin I_c), as the fraction of target binaries for which a planet with such parameters is detectable. This is calculated using all of the detection limits curves calculated based on Sect. 8.2. For every binary 20 detection limit curves were calculated, each corresponding to different planetary phases. In calculating C(P_c, m_c sin I_c) these 20 curves are treated as if they were 20 individual targets. This is the same approach as was used in earlier studies such as Cumming et al. (2008) and Mayor et al. (2011).

The completeness of our programme is shown in Fig. 10. White corresponds to 0% completeness, meaning that none of our targets are sensitive to planets of such minimum mass and period. We then use red gradient contours to denote increasing completeness, with dark red corresponding to 100%, meaning that all of our targets at all 20 planetary phases are sensitive to planets at such a period and minimum mass. Our program lacks completeness at short periods less than 50 days. This is a consequence of the stability limit restriction P_c ≳ 4P_bin, although we can alternatively say that we are complete down to the stability limit. At longer periods there is a drop in completeness due to the observational timespan, which varies between targets. For periods between roughly 50 and 3000 days the completeness contours follow a rough power law $m_{\text{c}}\sin I_{\text{c}}\propto P_{\text{c}}^{1/3}$. This is expected based on the radial velocity semi-amplitude equation (Eq. (1)).

In Fig. 11 we calculate for all 47 eclipsing binaries the smallest detectable planet mass. The solid navy line is calculated across all planet phases and periods, whereas the red dashed line is calculated across all phases but assumes a period of 2 yr. It is seen that for all 47 targets we have the ability to detect a circumbinary object less massive than 2.5 M_Jup. At smaller masses, for 30/47 of the targets we are sensitive down to 0.5 M_Jup. If we consider planets with periods up to 2 yr, then 25/47 of our targets are sensitive to 0.5 M_Jup mass planets.

The smallest detectable planet across the entire survey is 0.082 M_Jup, corresponding to EBLM J2011-07. However, we note that this system contains a tertiary stellar companion at a close period of 663 days. This third star does not hinder the detectability of interior planets, but may have inhibited any from forming in the first place (Muñoz & Lai 2015; Martin et al. 2015; Hamers et al. 2016; Xu & Lai 2016). The next smallest detectable mass is 0.110 M_Jup for EBLM J0310-31, for which no tertiary star has been found. For each target’s detectability curve in Appendix B we identify the smallest detectable planet mass and the corresponding orbital period, denoting it with a yellow star.

For more massive objects, that is circumbinary brown dwarfs and triple star systems, Fig. 10 shows that we have essentially 100% completeness, aside for very short and very long orbital periods. We impose an upper limit of 500 M_Jup on the completeness of tertiary star masses. This is because more massive stars would likely be bright enough to produce detectable spectral lines, whereas our entire sample consists of solely single-lined binaries. A smaller effect would be that brighter stars would dilute the already small eclipse depths of the M-dwarf secondary stars, potentially hindering the initial discovery of the inner binary in the EBLM programme. We elaborate upon circumbinary brown dwarfs and triple stars in Sect. 9.3.

9.2 Constraining the abundance of circumbinary planets

Since we do not have any confirmed discoveries of tertiary objects in the planetary domain, we can only place upper limits on their abundance. For this we use the same process as He et al. (2017). They conducted a survey of planets transiting brown dwarfs, and similarly had no confirmed detections with a comparable number of targets. The upper limit is calculated as $${\eta }_{\text{CBP}}=\frac{1-{(1-\kappa )}^{1/n_{\text{stars}}}}{C(P_{\text{c}},m_{\text{c}}\sin I_{\text{c}})},$$(8)

where n_stars = 47 is the number of stars in the BEBOP survey and κ is the desired confidence interval, for example κ = 0.95 for a 2σ confidence interval.

This upper limit is calculated within various parameter spaces, demarcated by the lower six white boxes in Fig. 10. The period bounds are roughly evenly separated in log space: 50, 245, 1200 and 6000 days. The first periods are chosen to illicit an easy comparison with the work of Santerne et al. (2016) for single stars, which we do in Sect. 10.3. The planet minimum mass bounds are chosen to span sub-Jupiter masses up to the deuterium burning limit, which marks the lower bound of the brown dwarf regime. The values are 0.5, 1.5, 4.5 and 13.5 M_Jup, which have roughly even log spacings.

Within each parameter space the completeness C(P_c, m_c sin I_c) is taken as the mean value within the box. In Table 2 we show the abundance constraints at 50, 1σ and 2σ confidence for each of these planet parameter spaces. In Fig. 10 to be conservative we only show the 2σ constraint.

A promising result is that the abundance constraints are only a weak function of orbital period. For example, we place a < 6.6% 2σ constraint on super-Jupiter circumbinary planets between 50 and 245 days, and only a marginally inferior constraint of < 8.9% 2σ for planets of the same mass but periods between 1200 and 6000 days. This is in contrast to the transit method, which more strictly favours short orbital periods, and consequently no abundance constraints have been made for periods greater than 300 days (Armstrong et al. 2014). It was predicted by Pierens & Nelson (2013) that the most massive circumbinary planets would be far from the stability limit, not close like the sub-Jupiter planets (see Fig. 2). Kepler-1647 (the top right upwards triangle in Fig. 10) follows this trend, and indeed would have been the easiest planet to detect in our programme (see Sect. 10.1). Whilst our results at present are unable to confirm the predictions of Pierens & Nelson (2013), it is apparent that a radial velocity survey is well-suited for such a task.

Fig. 9 Data for the EBLM J1146-42 triple star system. Panel a: periodogram of the residuals to a double Keplerian fit. The dash-dotted horizontal line at the top of the plot corresponds to a false alarm probability of 1%. The dotted line below is for 10%. There is therefore no statistically significant periodicity within the residuals. See Appendix B for more plots for this system. Panel b: bisector for each radial velocity measurement. Panel c: radial velocity residuals to a double Keplerian fit as a function of the bisector. A negative correlation would be a marker of stellar activity, but is not apparent here.

9.3 Circumbinary brown dwarfs and tight triple star systems

Figure 10 shows that whilst the BEBOP survey only has partial completeness within the circumbinary planetary mass domain, it has practically 100% completeness for circumbinary brown dwarfs and tertiary stars with moderately long periods. We therefore use this information to calculate the abundances of such objects. Since we have not detected any circumbinary brown dwarfs, we can only place an upper limit. For triple stars though we have five well-characterised systems, and hence can calculate an actual abundance. The top two boxes in Fig. 10 correspond to the abundance calculations for closely-orbiting low-mass triple star systems and circumbinary brown dwarfs. These results are also included in Table 2.

We define brown dwarfs as bodies with masses within the deuterium-burning regime: 13.5–80 M_Jup. Our BEBOP survey has almost 100% completeness for such objects on orbits between 50 and 6000 days, but with no confirmed discoveries. Using the same method as for the cirucmbinary planets, we constrain the abundance within this period range to be η_BD < 6.5% to 2σ confidence. The known circumbinary brown dwarf HD 202206 (Correia et al. 2005; Benedict & Harrison 2017) interestingly falls within this parameter space, with the planet just above the deuterium burning limit.

Our result here should be considered preliminary on account of the size of the BEBOP sample. Brown dwarfs companions to single Sun-like stars are inherently rare, at a rate of < 1%. This has ledto the coining of the phrase “brown dwarf desert” to represent the paucity of companion masses in between the planetary and stellar domains (Marcy & Butler 2000; Grether & Lineweaver 2006; Sahlmann et al. 2011; Kraus et al. 2011; Cheetham et al. 2015). If this rarity extends to brown dwarfs around binary stars, then to have not discovered one in a sample of 47 binaries is not surprising.

Beyond the brown dwarf regime we use our five characterised triple star systems to calculate the tertiary abundance between 50 and 6000 days for minimum masses between 80 and 475 M_Jup. This period range corresponds to the rough limits of detectability of our program, as tighter triples would be unstable and wider triples would not be well-characterised by our observational timespan. The mass range is chosen to be equal in log space to that for the brown dwarfs, whilst remaining less than the 500 M_Jup upper limit which we impose on the programme.

We follow the work of Mayor et al. (2011) to calculate the abundance as $${\eta }_{\text{triple}}=\frac{1}{n_{\text{stars}}}{\displaystyle \sum _{i=1}^{n_{\text{det}}}}\frac{1}{C_i(P_{\text{c}},m_{\text{c}}\sin I_{\text{c}})},$$(9)

where n_det = 5 is the number of detected triple stars, C_i(P_c, m_c sin I_c) is the completeness level in the parameter space for each of the five triples. Using Poisson statistics the 1σ uncertainty is calculated as $$\sigma =2\frac{{\eta }_{\text{triple}}}{\sqrt{n_{\text{det}}}}.$$(10)

The fraction of 50–6000 day tertiary stellar companions (between 80 and 475 M_Jup) is calculated to be 12.1 ± 5.4, within a 1σ confidence interval.

It is a future task to compare our work with complementary surveys of triple star systems. These include the work by Tokovinin et al. (2006) to directly image distant stellar companions to tight spectroscopic binaries, the studies of eclipse timing variations observed by Kepler (Borkovits et al. 2015) and the imaging survey for sub-stellar companions around binaries by Bonavita et al. (2016). A goal is to construct a mass distribution of tertiary objects, that can ultimately be compared with that of secondary objects (e.g. Grether & Lineweaver 2006; Triaud et al. 2017a).

Our results for circumbinary brown dwarfs and triple stars will be improved in the future by analysing the entire EBLM programme, for which 118 binaries have been published in Triaud et al. (2017a) and the entire sample numbers over 220. Whilst the measurements on the binaries not selected for BEBOP are too imprecise to aid the abundance calculations for circumbinary planets, almost all of these binaries permit the detection of circumbinary brown dwarfs and triple star systems.

For many binaries the observational baseline of a few years only permits us to see a drift in the radial velocity residuals to a single Keplerian fit. For these targets we will take new observations in the coming years to extend this baseline and better constrain the period and mass of the outer body.

Fig. 10 Completeness of the BEBOP radial velocity survey of 47 low-mass eclipsing binaries, as a function of the circumbinary minimum mass and period. Six different colour contours indicate the programme completeness between 0% (white) and 100% (dark red). The green circles near the top of the plot correspond to the five BEBOP triples, that is binaries with well-characterised tertiary stellar companions. The upright blue triangles in the bottom half of the plot represent the four Kepler transiting circumbinary planets with published masses: Kepler-16, -34, -35 and -1647. The inverted blue triangle represents the circumbinary brown dwarf HD 202206 (mc = 17.9 MJup, Pc = 1261 days) discovered using a combination of RVs and astrometry (Correia et al. 2005; Benedict & Harrison 2017). There are eight white boxes covering different parameter spaces within which we constrain the abundance of tertiary objects. The number in each box is the 2σ upper limit to the circumbinary abundance, given as a per cent. An exception is the top box which covers triple star systems. Since we have detections in this box, we derive an actual value for the abundance and its 1σ uncertainty.

Fig. 11 Histogram of the smallest detectable planet mass for the 47 eclipsing binaries in the BEBOP survey. The navy blue solid line calculates the minimum across all possible circumbinary periods and phases. The red dashed line calculates across all possible phases but has a fixed orbital period of 2 yr. We note that one target, EBLM J0500-46, is excluded in the histogram for Pp = 2 yr because its observing timespan is too short to be sensitive to planets at this period.

10 Comparisons with other surveys

10.1 Detectability of known Kepler planets

Masses can be derived for transiting circumbinary planets if they induce detectable eclipse timing variations in their host binary. This has been done for four of the discovered systems: Kepler-16, -34, -35 and -1647. These planets are shown in Fig. 10 as upright blue triangles. Three of the planets are just below our limits of detectability. Furthermore, the other circumbinary planets without measured masses are most likely smaller, and hence would have been even tougher to find.

Only a planet with Kepler-1647’s properties (1.27 M_Jup, 1108 day period) falls within the completeness of the BEBOP CORALIE programme (14.1% for these parameters). We estimate the probability that we would have found such a planet in our survey using the equation $$D=1-{0.9}^{n_{\text{stars}}C(P_{\text{c}},m_{\text{c}}\sin I_{\text{c}})},$$(11)

where we assume that around each binary there is a 10% chance of a gas giant planet existing, according to the abundance studies of Martin & Triaud (2014) and Armstrong et al. (2014). For Kepler-1647 D = 50%.

Overall, we have demonstrated the ability of our survey to detect planetary-mass circumbinary objects, typically smaller than 1 M_Jup. However, the masses of the known Kepler transiting planets are almost all sub-Saturn, unfortunately placing them slightly below the detection threshold for most of our targets.

10.2 Comparison with Armstrong et al. (2014) circumbinary abundance calculations

Using different approaches, Armstrong et al. (2014) and Martin & Triaud (2014) estimated that the abundance of circumbinary gas giants is roughly 10%. This is compatible with the upper limits we derive at the end of our BEBOP survey with CORALIE. Armstrong et al. (2014) remarked that no circumbinary planets > 1 R_Jup were detected⁸, with the authors implying a low abundance for masses greater than Jupiter’s. However, the mass-radius relation is roughly flat within a range of approximately 1–100 M_Jup (Baraffe et al. 2015). This means that Armstrong et al. (2014), based on radius, were not particularly sensitive to the frequency of circumbinary gas giants as a function of mass, something that can instead be tested thanks to radial-velocities. Our constrained abundances from Sect. 9.2 are consistent with a ~ 10% gas giant abundance, although given our limited sensitivity to Kepler-like circumbinary planets, our comparisons can only be preliminary at this time.

The limits we place on planet occurence rates have implications on the distribution of planetary orbital inclination with respect to their binary host. Kepler is mostly sensitive to coplanar configuration. Should planetary orbital planes follow a distribution with mutual inclinations with the binary of several degrees, most planets would go undetected. Each Kepler detection would therefore imply a much more abundant population than if all circumbinary systems were coplanar.

Armstrong et al. (2014) derived circumbinary abundances within various radius bins. Since we are typically sensitive to planets more massive than Jupiter, we choose to compare with their results for >10 R_⊕. We convert the lower bound to mass using a Jupiter density: 0.76 M_Jup. We take the upper bound to be 13.5 M_Jup, so we are considering all massive planetary objects. We use a period range of 50–300 days to match Armstrong et al. (2014).

Within this parameter space we calculate upper limits on the abundance of circumbinary planets as a function of σ_ΔI, the standard deviation of the mutual inclination distribution. We calculate those following Armstrong et al. (2014), where the Gaussian distribution of ΔI is convolved with a uniform distribution in cos ΔI. For each value of σ_ΔI the abundance is calculated by increasing the mass limits by a factor sin σ_ΔI. The results are shown in Fig. 12.

The radial velocity-derived abundances are significantly less affected by mutual inclinations than the transit results from Armstrong et al. (2014). Whilst the detectability of transits is a sensitive function of the misalignment (Martin & Triaud 2015; Martin 2017), the radial velocity detectability is a more shallow K_c ∝ sin ΔI (Eq. (1) and knowing I_bin ≈ 90°). The comparison in Fig. 12 shows that within this parameter space the upper limits placed by Armstrong et al. (2014) are only more constraining for a strictly coplanar system. For 5° and above the BEBOP results place tighter constraints on the presence of planets more massive than Jupiter.

Based on a preliminary comparison with the work of Armstrong et al. (2014) in Fig. 12, it appears that there does not exist a numerous population of misaligned giant circumbinary planets, which would have evaded transit detection but have been spotted by BEBOP. We estimate that the spread of the mutual inclination is less than 10°, as otherwise the high planetary abundance would have practically guaranteed a BEBOP discovery. This result is compatible with the Kepler analysis by Li et al. (2016).

For smaller planets, between 0.5 and 1 M_Jup and periods between 50 and 300 days, we can only place a rudimentary 2σ abundance constraint of <29%. In Sect. 11 we briefly discuss the recent upgrade of the BEBOP survey to the HARPS instrument on the 3.6 m ESO telescope, and the SOPHIE instrument on the 1.93 m OHP telescope, and how this will advance our constraints on circumbinary abundances.

Fig. 12 Percentage abundance of massive circumbinary planets with periods between 50 and 300 days, as a function of the underlying spread of planetary inclinations according to a Gaussian distribution convolved with an isotropic distribution of cos ΔI. In all cases the thick, shorter lines are limits at a 50% confidence interval and longer, thinner lines are for a 2σ confidence interval. At each value of ΔI there are two lines: BEBOP in red on the left and Armstrong et al. (2014) in blue on the right. The lower mass limit is 0.76 MJup, which corresponds to 10 R⊕, assuming a Jupiter density (so we can roughly compare between our RV survey and the transit survey of Armstrong et al. 2014). The upper limit is the deuterium-burning limit of 13.5 MJup. The plot on the right is the same as the left but zoomed in the vertical axis, better showing how the BEBOP abundance constraints are nearly independent on the inclination.

10.3 Comparison with Santerne et al. (2016) single star abundance calculations

The Santerne et al. (2016) SOPHIE radial velocity survey of Kepler transit candidates is one of the most comprehensive works in the literature on gas giant abundances around single stars. A comparison between the populations of planets orbiting one and two stars would shed light on the fundamental process of planet formation and evolution.

Santerne et al. (2016) calculates the gas giant abundance between 50 and 245 days to be 3.69 ± 0.84%. In this work they also re-analyse the data from Mayor et al. (2011) within these period ranges, calculating a remarkably similar gas giant abundance of 3.85 ± 0.85%. The samplein Santerne et al. (2016) has planets with masses between ~ 0.3 and 9.3 M_Jup. If we compare our work over the mass range 0.5−13.5 M_Jup, that is including all planetary mass bodies down to our rough sensitivity limit, we calculated a 2σ constrained circumbinary abundance of <9% for 50−245 day planets. This value is compatible with the Mayor et al. (2011) and Santerne et al. (2016) results, but we currently lack the precision to know whether the circumbinary abundance is truly greater or smaller.

11 Future prospects

This initial CORALIE survey had the capacity to detect Jupiter-mass planets on two-year orbital periods for roughly half of our sample of binary stars. Our preliminary understanding of typical circumbinary masses, based on the Kepler results, was insufficient to know whether high-mass circumbinary planets were particularly abundant. Our results imply that such high-mass circumbinaries planets are indeed not frequent. Furthermore, constrained circumbinary abundances are compatible with high-mass planets orbiting single stars, but our detection capabilities inhibit a statistically strong comparison. The results of this initial survey are primarily limited by the stability of the CORALIE instrument, which is pressurised but not thermalised, and the photon noise typical to a one metre class telescope.

We sought to extend the BEBOP survey to the HARPS spectrograph on the ESO 3.6 m telescope. We first conducted a short pilot programme that demonstrated that HARPS can reach radial-velocities with ~ 1−2 m s⁻¹ rms on single-line eclipsing binaries, implying a sensitivity to planets with masses as low as Neptune’s. Those results are due to be published shortly. Building on this, we proposed and were awarded a large programme on HARPS, with the first data being collected in April 2018. An extension has also been granted to the northern hemisphere, with a three-year large programme using the OHP 1.93 m telescope with the SOPHIE spectrograph. These new observations will enhance the survey described in these pages, reaching an order of magnitude deeper in mass, extending the period range, and covering a greater number of systems. Our future observations will enable proper comparisons between the properties of planets orbiting single stars, to planets orbiting binary stars. We will report our results in a series of BEBOP papers, of which this is the first installment.

Despite concentrating our current efforts on single-line binaries, we are also motivated in observing double-line binaries. This is important as it would provide a larger and brighter sample of binaries, but also a greater range in mass ratio. Measuring planet abundances as function of inner binary mass ratio would be insightful (Martin, in prep.). In Sect. 8.1 we showed no strong evidence that our secondary stars contaminate the measurement of radial-velocities. Our target most at risk of spectral contamination, EBLM J0425-46 (q = 0.527), which has a visual magnitude difference of 3.85, still produced a statistically perfect fit (${\chi }_{\text{red}}^2=0.73$). In the future, we will investigate at which flux ratios contamination becomes an issue. We will then learn how to deal with it.

12 Conclusion

Our BEBOP radial velocity survey for circumbinary planets has started by using the CORALIE spectrograph. Such planets possess significant intrigue, yet the statistics of a mere dozen or so confirmed cases mean that our insights to date are only preliminary. This first paper reports the results from eight years monitoring 47 single-lined eclipsing binaries. We had two primary intentions. The first was to verify whether high-mass (≳ 1 M_Jup) circumbinary planets were abundant, potentially owing to an unknown misaligned population, and to make an opportunistic discovery. Whilst we made no planetary detection, we have succeeded to place constraints on the presence of giant circumbinary planets out to orbital periods of a few years. The precision of our upper limits indicates that most massive circumbinary planets likely occupy orbits close to coplanar, with a standard deviation on the mutual inclination likely less than ~ 10°. In the process of planet-hunting, we also characterise five triple star systems. Two of these have orbits tighter than 2 AU, including one with evidence for non-Keplerian dynamical interactions.

Our second main goal was to test whether single-line binaries were amenable to circumbinary planet detection, a change in direction from the classically-observed double-line systems. We successfully demonstrated that the light from the faint secondary star can be ignored without inducing spectral contamination. This enabled statistically perfect binary orbital fits with a precision as small as 5 m s⁻¹, and a sensitivity to circumbinary planets down to 0.1 M_Jup in the best cases. Our initial programme is opening opportunities to seek smaller planets with more sensitive equipment.

Acknowledgements

The authors would like to attract attention to the help and kind attention of the ESO staff at La Silla and to the dedication of the many technicians and observers from the University of Geneva and EPFL, to upkeep the telescope, acquire the data that we present here, and bring the fondue cheese all the way from Switzerland. We would also like to acknowledgethat the Euler Swiss Telescope at La Silla is a project funded by the Swiss National Science Foundation (SNSF). Over the time required to collect and analyse the data, DVM was supported by the SNSF, the University of Chicago and the Université de Gèneve. A.H.M.J.T. has received funding from the SNSF, the University of Toronto and the University of Cambridge. P.M. gratefully acknowledges support provided by the UK Science and Technology Facilities Council through grant number ST/M001040/1. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no 803193/BEBOP). We also thank François Bouchy, Dan Fabrycky and Rosemary Mardling for reviewing an earlier version of this work that appeared in D.V.M.’s PhD thesis. Finally, we kindly thank our editor and anonymous referees for reading our paper and making comments that led to significant improvements. This publication makes use of data products from two projects, which were obtained through the Simbad and VizieR services hosted at the CDS-Strasbourg: the Two Micron All Sky Survey (2MASS), which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation (Skrutskie et al. 2006); the Naval Observatory Merged Astrometric Dataset (NOMAD), which is project of the US Naval Observatory (Monet et al. 2003).

Note added to the proof. We used the Barycentric Julian Dates in our analysis. Our results are based on the equatorial solar and jovian radii and masses taken from Allen’s Astrophysical Quantities (Cox 2000).

Appendix A Data tables of the BEBOP target list, model selection and derived orbital parameters

Appendix B Radial velocities and detection limits for the 47 BEBOP systems

Fig. B.1 Top panel: radial velocity measurements over time (red) and fitted model (black) and the residuals to the model fit (O–C). Middle panel: phase-folded velocities on the binary period, where the colour indicates the observationdate. Bottom panel: detection limits as a function of the detectable tertiary period and minimum mass, where different colours are used for 20 different tested tertiary orbital phases uniformly sampled over 360°, and the yellow star highlights the smallest detectable mass for any parameters.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

Fig. B.1 continued.

References

All Tables

All Figures

	Fig. 1 In dark blue, a histogram of the Kepler eclipsing binary catalogue for all periods longer than 1 day. Data is taken from http:// keplerebs.villanova.edu/ as of May 2017, based on a catalogue that is first outlined in Prsa et al. (2011). Red dashed lines correspond to the nine binaries known to host transiting circumbinary planets. In light blue, a histogram of the 47 BEBOP binaries.
In the text

	Fig. 2 Left panel: planet and binary orbital periods, with dashed lines of constant period ratio. Right panel: planet periapse distance and binary apoapse distance, with dashed lines of constant scaled minimum distance, rmin = [ap(1 − ep)]∕[abin(1 + ebin)]. Systems with a single detected planet are shown as blue squares, whereas the three-planet Kepler-47 system is shown as grey circles.
In the text

	Fig. 3 In light blue, a histogram of the mass ratio of the 47 BEBOP binaries. The vertical red dashed lines correspond to the mass ratios of the Kepler binaries known to host transiting circumbinary planets.
In the text

Fig. 4 Time-series of 1519 radial velocity observations of the 47 eclipsing binaries in the BEBOP programme. Red diamonds are observations for 1700 s or longer. Blue squares are for shorter observations. It is seen that all binaries typically receive two long observations initially before a series of short observations as part of the EBLM programme. A series of long observations typically commenced near the end of 2013 as part of the BEBOP survey.

In the text

Fig. 5 Radial velocity residuals (“Observed minus Calculated” or O–C for short) of the best-fitting model to all BEBOP binaries forall 1519 observations, stacked on top of each other as a function of the scaled radial velocity value δ, defined according to Eq. (7). The δ value says how far the observed radial velocity observation was from the systemic velocity, scaled by the average FWHM for that target. The two roughly horizontal dark blue lines are the root mean squared (RMS) values of the residuals, split into positive and negative values. If the RMS were maximised near δ = 0 then this would be indicative of wide-spread spectral contamination from the secondary star, but this is evidently not the case.

In the text

Fig. 6 Top panel: radial velocities of the target EBLM J0425-46 over 3.15 yr and the selected eccentric single Keplerian model. Error bars are typically 10–20 m s−1 and too small to see at this scale. Bottom panel: residuals to the fitted model, with 1σ error bars. Blue boxes highlight two marginal outliers, with a dashed line connecting the radial velocity measurement, which is near0 m s−1, and the residual. Both outliers are less than 2σ from the model.

In the text

Fig. 7 Example of the detection limits calculated for EBLM J2011-71, which is one of the BEBOP targets with the highest precision. Each line represents the smallest detectable minimum circumbinary planet mass (mc sin Ic) as a functionof the planet’s period (Pc) for a different orbital phase of the planet. There are 20 orbital phases tested, all equally spaced. The blue square near the top of the plot is the detected circumbinary object in the EBLM J2011-71, which is in fact not a circumbinary planet but rather alow-mass tertiary star.

In the text

Fig. 8 Detection limit curves in black for two targets in the BEBOP programme: J0035-69(left panel) and J0310-31 (right panel). Each curve represents the smallest mc sin Ic for a putative circumbinary planet at a given period, according to the periodogram measure of detectability described in Sect. 8.2. For each target there are 20 black curves, one for each of the tested orbital phases of the injected planet. At seven discrete periods a series of tests were run to recover n-body simulated circumbinary planets with the YORBIT genetic algorithm. At these periods a green circle indicates a successful recovery, whereas a red diamond indicates a failure.

In the text

Fig. 9 Data for the EBLM J1146-42 triple star system. Panel a: periodogram of the residuals to a double Keplerian fit. The dash-dotted horizontal line at the top of the plot corresponds to a false alarm probability of 1%. The dotted line below is for 10%. There is therefore no statistically significant periodicity within the residuals. See Appendix B for more plots for this system. Panel b: bisector for each radial velocity measurement. Panel c: radial velocity residuals to a double Keplerian fit as a function of the bisector. A negative correlation would be a marker of stellar activity, but is not apparent here.

In the text

Fig. 10 Completeness of the BEBOP radial velocity survey of 47 low-mass eclipsing binaries, as a function of the circumbinary minimum mass and period. Six different colour contours indicate the programme completeness between 0% (white) and 100% (dark red). The green circles near the top of the plot correspond to the five BEBOP triples, that is binaries with well-characterised tertiary stellar companions. The upright blue triangles in the bottom half of the plot represent the four Kepler transiting circumbinary planets with published masses: Kepler-16, -34, -35 and -1647. The inverted blue triangle represents the circumbinary brown dwarf HD 202206 (mc = 17.9 MJup, Pc = 1261 days) discovered using a combination of RVs and astrometry (Correia et al. 2005; Benedict & Harrison 2017). There are eight white boxes covering different parameter spaces within which we constrain the abundance of tertiary objects. The number in each box is the 2σ upper limit to the circumbinary abundance, given as a per cent. An exception is the top box which covers triple star systems. Since we have detections in this box, we derive an actual value for the abundance and its 1σ uncertainty.

In the text

Fig. 11 Histogram of the smallest detectable planet mass for the 47 eclipsing binaries in the BEBOP survey. The navy blue solid line calculates the minimum across all possible circumbinary periods and phases. The red dashed line calculates across all possible phases but has a fixed orbital period of 2 yr. We note that one target, EBLM J0500-46, is excluded in the histogram for Pp = 2 yr because its observing timespan is too short to be sensitive to planets at this period.

In the text

Fig. 12 Percentage abundance of massive circumbinary planets with periods between 50 and 300 days, as a function of the underlying spread of planetary inclinations according to a Gaussian distribution convolved with an isotropic distribution of cos ΔI. In all cases the thick, shorter lines are limits at a 50% confidence interval and longer, thinner lines are for a 2σ confidence interval. At each value of ΔI there are two lines: BEBOP in red on the left and Armstrong et al. (2014) in blue on the right. The lower mass limit is 0.76 MJup, which corresponds to 10 R⊕, assuming a Jupiter density (so we can roughly compare between our RV survey and the transit survey of Armstrong et al. 2014). The upper limit is the deuterium-burning limit of 13.5 MJup. The plot on the right is the same as the left but zoomed in the vertical axis, better showing how the BEBOP abundance constraints are nearly independent on the inclination.

In the text

Fig. B.1 Top panel: radial velocity measurements over time (red) and fitted model (black) and the residuals to the model fit (O–C). Middle panel: phase-folded velocities on the binary period, where the colour indicates the observationdate. Bottom panel: detection limits as a function of the detectable tertiary period and minimum mass, where different colours are used for 20 different tested tertiary orbital phases uniformly sampled over 360°, and the yellow star highlights the smallest detectable mass for any parameters.

In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text

	Fig. B.1 continued.
In the text