\section{Introduction}
Circumbinary planets are planets that orbit around both stars of a binary star system. Long postulated \citep{Borucki:1984,Schneider:1994lr}, the first unambiguous discovery of a circumbinary planet came with Kepler-16 \citep{Doyle:2011vn}, detected by identifying three transits within the lightcurve of an eclipsing binary system monitored by NASA's {\it Kepler} mission \citep{Borucki:2011qy}. {\it Kepler} went on to detect another 13 transiting circumbinary planets, in 11 systems \citep{Martin2018,Socia:2020}, with another two systems found using {\it TESS} \citep{Kostov2020,Kostov2021}. A number of circumbinary planets systems are suspected from eclipse timing variations of binaries on the main sequence \citep[e.g.][]{Borkovits2016,Getley2017}, and stellar remnants \citep[e.g.][]{Marsh:2013lr,Han2017} but most are disputed \citep[e.g.][]{Mustill2013}, and some disproven \citep[e.g.][]{Hardy2015}. Other detections include HD~106906~b, in direct imaging \citep[][]{Bailey2014} and OGLE-2007-BLG-349L(AB)c with the microlensing method \citep{Bennett:2016}. 

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

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

Overall, radial velocities are essential to create a sample of circumbinary planets that is both greater in number and less biased than the transit sample. This will allow a deeper understanding of circumbinary planets: their occurrence rate \citep[][]{Martin:2014lr,Armstrong:2014fk}, multiplicity \citep{Sutherland:2019,Orosz:2019}, formation and evolution \citep[e.g.][]{Chachan2019,Pierens2020, Penzlin2021}, and dependence on binary properties \citep{Munoz:2015,Martin:2016,Li:2016,Martin:2019b}.

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

Thanks to the advent of exoplanet transit experiments, increasing amounts of low-mass, single-lined eclipsing binaries are being identified \citep[e.g.][]{Triaud:2013lr,Triaud:2017yu,Boetticher2019,Lendl2020,Mireles2020,Acton2020}. The BEBOP survey was constructed solely using sufficiently faint secondaries, to avoid detection with spectrographs, such that we could in principle reach a radial velocity precision comparable to that around single stars of the same brightness. In principle, this ought to provide an accuracy of order ${1~\rm m\,s}^{-1}$. Our survey is ongoing, and uses the CORALIE, SOPHIE, HARPS, and ESPRESSO spectrographs. Preliminary results were published in \citet{Martin2019} (BEBOP I). In \citet[][under review, BEBOP II]{Standing2021}, we describe our observational protocol, the methods we use to detect planets, as well as how we produce detection limits.

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





\section{Velocimetric observations on Kepler-16}

Between 2016-07-08 and 2021-06-23, we collected 143 spectra using the high-resolution, high-precision, fibre-fed SOPHIE spectrograph, mounted on the 193cm at Observatoire de Haute-Provence, in France \citep{Perruchot08}. The Journal of Observations can be found in Table~\ref{tab:data}. All observations were conducted in HE mode (High Efficiency) where some of the instrumental resolution is sacrificed from 75,000 to 40,000 in favour of a $2.5\times$ greater throughput. We chose this since whilst Kepler-16 is the brightest circumbinary system, it is relatively faint for SOPHIE, with $V \sim 12.0$. 

SOPHIE has two fibres; the first stayed on target, while the other was kept on the sky in order to remove any contribution from the Moon-reflected sunlight. Standard calibrations were made at the start of night as well as roughly every two hours throughout the night to monitor the instrument's zero point. In addition, we observed one of three standards (HD~185144, HD~9407 and HD~89269~A) in HE mode nightly, which we used to track and correct for any long-term instrumental drift following procedures established in \citet{santerne2014sophie} and \citet{Courcol15}. 

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

The cross-correlation software produces two key metrics of the shape of the CCF, its FWHM (Full Width at Half Maximum), and its bisector span  \citep[as defined in][]{Queloz:2001}. In addition, we measure the $\rm H\alpha$ stellar activity indicator following \citet{boisse2009stellar}. These are provided in Table~\ref{tab:data}.

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

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

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

\section{Modelling of the radial velocities}

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


\subsection{Outlier removal}\label{sec:outliers}

We searched our 141 data for measurements\footnote{The first series of 14 measurements were obtained from a catalogue containing erroneous proper motions and epochs, which in turn created a $1.5~\rm km\,s^{-1}$ effect in the radial velocity as the Earth's motion was over compensated. This can be corrected easily, but the error remains within the archival data.} coinciding with a primary eclipse, and likely affected by the Rossiter-McLaughlin effect \citep{Rossiter:1924qy,McLaughlin:1924uq,Winn:2011wd,Triaud2018}. Fortunately no measurements needed to be excluded for this reason.

We also realised that a number of measurements were likely taken under adverse conditions. This is apparent from an unusually low signal-to-noise ratio, but also from large values of the bisector span. Typically used as a stellar activity, or a blend indicator \citep[][]{Queloz:2001,Santos:2002fk}, the span of the bisector slope (bisector span, or bis span) effectively informs us that the line shape varied and therefore that the mean of the cross-correlation function is likely affected. We took the mean of the bisector span measurement and removed all measurements in excess of $3\sigma$ away from the mean. Five measurements are excluded this way, all with a bisector span $\gtrsim \pm 100~\rm m\,s^{-1}$. Excluded measurements are reported in the Journal of Observations in Table~\ref{tab:data} with a flag. 

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

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

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




\begin{figure}
	\includegraphics[width=\columnwidth]{periodogramsS.pdf}
    \caption{Lomb-Scargle periodogram of Kepler-16's radial velocities (top) and $\rm H\alpha$ (bottom). The radial velocities are shown after removing the binary motion, and a cubic function. The four lines are, from bottom to top, the $10\%$, $1\%$, $0.1\%$ and $0.01\%$ false alarm probabilities. There is a highly significant peak around 230 days (vertical red dotted line) that is present in the radial velocities but not in $\rm H\alpha$. The $\rm H\alpha$ measurements contain significant periodogram power at $\gtrsim 2000~\rm days$, indicating a long-term trend in the chromospheric emission from the primary star.}
    \label{fig:periodogram}
\end{figure}

\subsection{Analysis using the genetic algorithm {\sc Yorbit}}

{\sc Yorbit} is a radial velocity fitting tool used for exoplanet detection, described in \citet{Segransan:2011}. It assumes Keplerian orbits (see Sect~\ref{sec:orbital_elements}). {\sc Yorbit} first performs a Lomb-Scargle periodogram of the observations, which is then used to initiate a genetic algorithm that iterates over the orbital period $P$, the eccentricity $e$, the argument of periastron $\omega$, the semi-amplitude $K$, a reference time $T_0$, and one systemic velocity per dataset $\gamma$ \citep[for conventions, see][]{Hilditch:2001uq}. Once the algorithm has converged on a best solution, the final parameters are estimated by using a least-square fit. The tool has been routinely used to identify small planetary companions successfully \citep[e.g.][]{Mayor:2011fj,Bonfils:2013fk}, and represents a more traditional, and possibly a more recognised way of identifying a new planetary system than the nested sampler we use subsequently (see Sect.~\ref{sec:kima}). However {\sc Yorbit} cannot do a Bayesian model comparison. Instead it computes a False Alarm Probability (FAP) by performing a bootstrap on the data thousands of times and computing for each iteration a Lomb-Scargle periodogram.


In a first instance, {just one Keplerian is adjusted to the SOPHIE data, with {\sc Yorbit} automatically finding the most prominent signal, that of} the secondary star. {Following this step, we remove the secondary's signal and search the resulting residuals with a} periodogram. {This periodogram shows} excess power around 230 day with a $\rm FAP~\sim 0.1\%$ as well as excess power for a signal longer than the range of dates we observed for, which we later associate to a magnetic cycle. {To isolate the signal of Kepler-16~b,} we fit the secondary's Keplerian to the data, alongside a polynomial function, which is used to detrend {that longer signal. Once {\sc Yorbit} has converged, we search the residuals again with a periodogram, which provides} a $\rm FAP~\ll 0.01\%$ (Fig.~\ref{fig:periodogram}), clearly detecting Kepler-16~b as an additional periodic sinusoidal signal. The FAP obtained with a cubic {detrending} function is one order of magnitude better than that obtained with a quadratic function, so we chose the former as our baseline detrending. 

{To obtain results on the system's parameters, we perform a final fit to the data assuming two Keplerians and a cubic detrending function. Results of that fit are} found in Table~\ref{tab:results}. The orbital parameters of the planet are compatible with those produced in \citet{Doyle:2011vn} (see table~\ref{tab:results} and section~\ref{sec:orbital_elements} for a discussion). Our final fit produces a reduced $\chi^2_\nu = 1.17\pm0.14$, implying no additional complexity is needed to explain the data, and supporting our choice for a circular planetary orbit\footnote{Making a fit with a quadratic function we obtain $\chi^2_\nu = 1.47\pm0.15$, for one fewer parameter, justifying our choice for the cubic drift.}. In addition, this shows that we can achieve photon-noise precision on single-lined binary to detect circumbinary planets. The model fit to the data is depicted in Fig.~\ref{fig:orbit}.

Including seven outliers described in Sect.~\ref{sec:outliers}, neither the FAP nor the reduced $\chi^2_\nu$ are significantly affected. 

\begin{figure}



	\includegraphics[width=\columnwidth]{kepler16rv_newS.pdf}
    \caption{Best fit adjustment to the SOPHIE radial velocity data. {\it a:} Doppler reflex motion caused by the secondary star. {\it b:} Doppler reflex motion caused by the circumbinary planet. {\it c:} Cubic drift associated with a magnetic cycle. {\it d:} radial velocities as a function of time with a binary+planet+cubic function model. Residuals are displayed below.}
    \label{fig:orbit}
\end{figure}


\begin{table}
	\centering
	\caption{Results of our analysis of the SOPHIE radial velocities only, {after removing outliers, that show} the fit's Jacobi parameters and their derived physical parameters. They are compared to previous results with $1\sigma$ uncertainties provided in the form of the last two significant digits, within brackets. Dates are given in BJD - 2,450,000. We adopt the {\sc Kima} column as our results. }
	\label{tab:results}
	\begin{tabular}{lllll} 
		\hline
		\hline
		\multicolumn{2}{l}{Parameters \& units} & {\sc Yorbit} & {\sc Kima} & Doyle+ \citeyearpar{Doyle:2011vn}\\
		\hline
		\multicolumn{5}{l}{\it binary parameters}\\
		$P_{\rm bin}$        &day             & 41.077779(54) & 41.077772(51) & 41.079220(78)\\
		$T_{0, \rm bin}$     &BJD             & 8558.9640(44) & 7573.0984(47) & --\\
		$K_{1, \rm bin}$     &$\rm m\,s^{-1}$ & 13\,678.2(1.5)  & 13\,678.7(1.5) & --\\
		$e_{\rm bin}$        &--              & 0.15989(11)   & 0.15994(10) & 0.15944(62)\\
		$\omega_{\rm bin}$   &deg             & 263.661(40)   & 263.672(40) & 263.464(27)\\
		\\
		\multicolumn{5}{l}{\it planet parameters}\\
		$P_{\rm pl}$         &day             & 228.3(1.8)  & 226.0(1.7) & 228.776(37)\\
		$T_{0, \rm pl}$      &BJD             & 8532.5(4.4) & 7535(92) & --\\
		$K_{1, \rm pl}$      &$\rm m\,s^{-1}$  & 12.8(1.5)   & 11.8(1.5) & --\\
		$e_{\rm pl}$         &--              & 0 (fixed)   & $<$0.21   & 0.0069(15)\\
		$\omega_{\rm pl}$    &deg             & --          &   231(65)    & $318^{(+10)}_{(-22)}$\\
		\\
		\multicolumn{5}{l}{\it system parameters}\\
		$\gamma$             &$\rm km\,s^{-1}$  & -33.8137(69) & -33.8065(45) & -32.769(35) \\
		$\sigma_{\rm jitter}$ &$\rm m\,s^{-1}$  & --           & 0.070 $^{+1.104}_{-0.067}$ & -- \\
		\hline
		\multicolumn{5}{l}{\it derived parameters}\\
		$M_1$         &$\rm M_\odot$& 0.654(17)$^*$ & 0.654(17)$^*$ & 0.6897(35)\\
		$M_2$         &$\rm M_\odot$& 0.1963(31) & 0.1964(31) & 0.20255(66)\\
		$m_{\rm pl}$  &$\rm M_{Jup}$& 0.345(41) & 0.313(39) & $0.333 (16)$\\
		$a_{\rm bin}$ &$\rm AU$     & 0.2207(18) & 0.2207(17) & $0.22431 (35)$\\
		$a_{\rm pl}$ &$\rm AU$      & 0.6925(67) & 0.6880(58) & $0.7048 (11)$\\
		\hline
	\end{tabular}
	{* adopted from \citet{Bender:2012lr}}
\end{table}



\subsection{Analysis using the diffusive nested sampler {\sc Kima}}\label{sec:kima}
{\sc Kima} is a tool developed by \citet{kima} which fits a sum of Keplerian curves to radial velocity data. It samples from the posterior distribution of Keplerian model parameters using a Diffusive Nested Sampling algorithm by \citet{Brewer2016}. Diffusive Nested sampling allows the sampling of multi-modal distributions, such as those typically found in exoplanetary science and radial velocity data \citep{Brewer2015}, evenly and efficiently. 

{\sc{Kima}} can treat the number of planetary signals ($N_{\rm p}$) present in an RV data set as a free parameter in its fit. Since the tool also calculates the fully marginalised likelihood (evidence) it allows for Bayesian model comparison \citep{Trotta2008} between models with varying $N_{\rm p}$. A measure of preference of one Bayesian model over another can be ascertained by computing the "Bayes Factor" \citep{KassRaftery1995} between the two. The Bayes factor is a ratio of probabilities between the two competing models. Once a value for the Bayes factor has been calculated, we can compare it to the so-called ``Jeffreys' scale" (see \citet{Trotta2008} for more details) to rate the strength of evidence of one model over another.

A more extensive description of our use of {\sc{Kima}} in the context of the BEBOP survey can be found in \citet{Standing2021}.
For the analysis of the Kepler-16 system prior distributions were chosen similarly to those used in \citet{faria2020} with the following notable adaptions. 
We treat the secondary star as a \textit{known object} with tight uniform priors on its orbital parameters. 
A log-uniform distribution was used to describe the periods of any additional signals, from $4 \times P_{\rm bin}$ to $1\times10^4$ days. This inner limit on period is set by the instability limit found in binary star systems \citep{Holman:1999lr}. More details, particularly on the priors we use, can be found in \citet{Standing2021}. {Just like for the previous analysis using {\sc Yorbit}, {\sc Kima} is only deployed on SOPHIE data, and excluding outliers described in Sec.~\ref{sec:outliers}. }

Our {\sc{Kima}} analysis of the Kepler-16 data yields a Bayes Factor $\rm BF > 10000$ in favour of a three Keplerian model (secondary star, planet and cubic drift). Our posterior shows over-densities at orbital periods of $\approx 230$ and $\approx 2000$ days, corresponding to the signal of Kepler-16~b and the cubic drift seen in {\sc{Yorbit}}. We then apply the clustering algorithm {\sc{HDBSCAN}} \citep{McInnes2017} to isolate and extract the resulting planetary orbital parameters, which can be found in table~\ref{tab:results}.


\subsection{Note on converting fitted parameters to physical values}

To convert our semi-amplitudes into masses for $M_2$ and $m_{\rm pl}$, for the secondary star and planet masses (respectively), we adopt a primary star mass ($M_1$) from \citet{Bender:2012lr}. 
Software written for exoplanetary usage usually assumes that $m_{\rm pl} \ll M_\star$ (including {\sc Yorbit} and {\sc Kima}), however this assumption is no longer valid when comparing $M_2$ to $M_1$, and a circumbinary planet to both. 

First we find $M_2$ iteratively using the mass function, following the procedure described in \citet{Triaud:2013lr}. Then, we estimate $m_{\rm pl}$ from $K_{1,\rm pl}$ by using the combined mass $M_1 + M_2$. This is because whilst we are only measuring the radial velocity signature of the primary star, the gravitational force of the planet acts on the barycentre of the binary.
Had we not done this extra conversion step, we would find significant differences, with erroneous $M_2 = 0.165~\rm M_\odot$ and $m_{\rm pl} = 0.29~\rm M_{Jup}$ for the {\sc Yorbit} results.

Differences between our derived parameters (bottom part of table~\ref{tab:results}) and those from \citet[][]{Doyle:2011vn} are mainly explained by our adoption of the more accurate $M_1$ mass from \citet{Bender:2012lr} rather than using the value from \citet[][]{Doyle:2011vn}.

\subsection{Note on circumbinary planets' orbital elements}\label{sec:orbital_elements}

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

For example, \citet{Doyle:2011vn} provide a highly precise value of $P_{\rm pl} = 228.776 \pm 0.037$ days, yet this osculating period will vary by approximately $\pm 5$ days over a timescale of just years. Our measured period of $P_{\rm pl}=226.0\pm1.7$ days has a much higher error with this uncertainty being a combination of our radial velocities being less constraining on the period, and our assumption of a static orbit. The value we obtain with {\sc Kima} is $1.6\sigma$ compatible with the value found by \citet{Doyle:2011vn}.

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

Finally, to validate our assumption that we cannot measure osculating elements with our current data, we used tools described in \citet{Correia:2005lr} and \citet{Correia_etal_2010} to perform a N-body fit to the radial velocities. This fit finds no improvements in $\chi^2$ indicating that Newtonian effects indeed remain below the detectable threshold. 

\subsection{A magnetic cycle, and constraints on additional planets}

To assess the presence of an external companion causing the additional polynomial signal we notice in our data, we force a two-Keplerian model to fit the data (a circular orbit for the planet and a free-eccentricity orbit for the binary) and analyse the residuals. We find a signal reaching a  $\rm FAP~< 0.01\%$ at orbital periods exceeding the timespan of our data ($\gtrsim 2,000 \rm~days$). Long-term drifts can sometimes be caused by magnetic cycles since stellar spots and faculae tend to suppress convective blue-shift, producing a net change in the apparent velocity of a star \citep{Dravins:1985,Meunier:2010,Dumusque:2011}. We perform a Lomb-Scargle periodogram on the $\rm H\alpha$ activity indicator we extracted from the SOPHIE spectra, and find a $\rm FAP<0.01\%$. We therefore interpret this long-term radial velocity drift as a magnetic cycle, with a timescale longer than the timespan of our observations.

The validity of this interpretation can be tested since duration of stellar magnetic cycles scales with stellar rotation periods \citep[for single star;][]{Suarez:2016}. We measure the primary star's rotation from the $v \sin i_1$ obtained during the Rossiter-McLaughlin effect by \citet{Winn:2011wd} to the primary's stellar radius obtained by \citet{Bender:2012lr} and obtain a primary rotation $P_{\rm rot,1} = 35.68\pm 1.04~\rm days$.
Following the relation of \citet{Suarez:2016}, with the observed stellar rotation we ought to expect a magnetic cycle on a timescale of 1900 to 2400 days, which is entirely compatible with the $\rm H\alpha$ signal and the long-term radial velocity drift.

We now use {\sc Kima} as done in \citet{Standing2021} to compute a detection limit on the presence of additional but undetected planetary companions. We first remove the highest likelihood model with two Keplerians (the planet and the long-term drift) from the data, then force $N_{\rm p} = 1$ to obtain a map of all remaining signals that are compatible with the data, but remain formally undetected (the binary is also adjusted at each step). This map is shown as a greyscale density on Fig.~\ref{fig:detectionlimit}. The $99\%$ contour informs us that we are sensitive to companions below the mass of Kepler-16~b up to orbital periods of $\sim3,000$ days. 
This complements \citet{Martin2021}'s work who placed detection limits down to Earth-radius planets but only out to periods of 500 days, re-analysing the {\it Kepler} photometry. 

Overall, a picture is emerging of Kepler-16\,b as a lonely planet, which has implications for the formation and migration of multi-planet systems in the presence of a potentially destabilising binary \citep[e.g.][]{Sutherland:2019}.

\begin{figure}
	\includegraphics[width=\columnwidth]{Kepler-16_hexbin.png}
    \caption{Detection sensitivity to additional planets plotted as semi-amplitude $K_{1,\rm pl}$ as a function of $P_{\rm pl}$. The hexagonal bins depict the density of posterior samples obtained from three separate {\sc Kima} runs applied on the Kepler-16 radial velocity data data after removing two Keplerian signals shown with green dots. The faded blue lines show detection limits calculated for each of the three runs on the system. The solid blue line shows the detection limit calculated from all posterior samples combined. The solid red line is the outline of the posterior that led to the detection of Kepler-16~b and of a long-term trend associated with a magnetic cycle. The green dots represent the two signals removed from the data to compute the blue detection limit.}
    \label{fig:detectionlimit}
\end{figure}

\section{Discussion}\label{sec:discussion}

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

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

Finally we discuss the detectability of circumbinary systems.
We only take the 20 first measurements we collected and run {\sc Kima} measuring the Bayes Factor to assess the detectability of Kepler-16~b, where individual uncertainties are similar to the semi-amplitude of the signal, $K_{1,\rm pl}$. We repeat the procedure, measuring BF for each increase of five measurements until we reach a sub-sample containing the 50 first radial velocity measurements (they roughly cover four orbital periods of the planet). We reach a BF $>150$ (the formal threshold for detection) with the first 40 measurements. We repeat this procedure but instead randomly select 50 measurements within the first three years of data. We run {\sc Kima} and measure the Bayes Factor for the first 20 epochs of this sequence of 50 random epochs. We then increase from 20 to 25 until reaching 50. We find that the BF~$=150$ threshold is also passed at 40 measurements. We plot our results in Fig.~\ref{fig:test}. The results between both series of tests are broadly consistent, except when we only use 25 measurements, with the Bayes Factor growing log-linearly with increasing number of measurements. Thanks to these tests we conclude that just 40 to 45 measurements would have been to formally detect Kepler-16\,b with radial velocities.

\begin{figure}
	\includegraphics[width=\columnwidth]{Kepler-16_firstvsrandom.png}
    \caption{Bayes factors obtained for subsets of data with increasing size. The blue line and triangles indicate the Bayes factors obtained by running {\sc Kima} on the first \textit{N} data points obtained on Kepler-16. The black line and dots show the same but for randomly sampled \textit{N} data points within the first 3 years of data. Coloured regions represent various thresholds corresponding to improvements in evidence in favour of the more complex binary $+$ planet model as per the Jeffreys' scale detailed in \citet{Standing2021}.}
    \label{fig:test}
\end{figure}
\section*{Acknowledgements}

First, we would like to thank the staff, in particular the night assistants, at the Observatoire de Haute-Provence for their dedication and hard work, particularly during the COVID pandemic. 
The data were obtained first as part of a DDT graciously awarded by the OHP director (Prog.ID 16.DISC.TRIA), and then by a series of allocations through the French PNP (Prog.IDs 16B.PNP.HEB2, 17A.PNP.SANT, 17B.PNP.SAN2, 18A.PNP.SANT, 18B.PNP.SAN1, 19A.PNP.SANT). The French group acknowledges financial support from the French Programme National de Plan\'etologie (PNP, INSU).
This research received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement n$^\circ$ 803193/BEBOP) and from the Leverhulme Trust (research project grant n$^\circ$ RPG-2018-418). EW acknowledges support from the ERC Consolidator Grant funding scheme (grant agreement n$^\circ$772293/ASTROCHRONOMETRY). 
Support for this work was provided by NASA through the NASA Hubble Fellowship grant number HF2-51464 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.
P.C. thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant \#1829740, the Brinson Foundation, and the Moore Foundation; her participation in the program has benefited this work. 
SH acknowledges CNES funding through the grant 837319. ACC acknowledges support from
the Science and Technology Facilities Council (STFC) consolidated
grant number ST/R000824/1.
{ACMC acknowledges support by CFisUC projects (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal.}


\section*{Data Availability}

The radial velocity data are included in the appendix of this paper. The reduced spectra are available at the SOPHIE archive: \href{http://atlas.obs-hp.fr/sophie/}{http://atlas.obs-hp.fr/sophie/}. Data obtained under Prog.IDs 16.DISC.TRIA, 16B.PNP.HEB2, 17A.PNP.SANT, 17B.PNP.SAN2, 18A.PNP.SANT, 18B.PNP.SAN1, and 19A.PNP.SANT.





\bibliographystyle{mnras}
