New methods for radial-velocity measurements of double-lined binaries, and detection of a circumbinary planet orbiting TIC 172900988

Ongoing ground-based radial-velocity observations seeking to detect circumbinary planets focus on single-lined binaries even though over nine in every ten binary systems in the solar-neighbourhood are double-lined. Double-lined binaries are on average brighter, and should in principle yield more precise radial-velocities. However, as the two stars orbit one another, they produce a time-varying blending of their weak spectral lines. This makes an accurate measure of radial velocities difficult, producing a typical scatter of 10-15m/s. This extra noise prevents the detection of most orbiting circumbinary planets. We develop two new data-driven approaches to disentangle the two stellar components of a double-lined binary, and extract accurate and precise radial-velocities. Both approaches use a Gaussian Process regression, with the first one working in the spectral domain, whereas the second works on cross-correlated spectra. We apply our new methods to TIC 172900988, a proposed circumbinary system with a double-lined binary, and detect a circumbinary planet with an orbital period of 150 days, different than previously proposed. We also measure a significant residual scatter, which we speculate is caused by stellar activity. We show that our two data-driven methods outperform the traditionally used TODCOR and TODMOR, for that particular binary system.


INTRODUCTION
Planets orbiting around both stars of a binary star system are called ptype or circumbinary planets.The discovery of Kepler-16 (Doyle et al. 2011) marked the first confirmation of the existence of these longtheorised planets (Borucki & Summers 1984;Schneider 1994).Since then, the Kepler and TESS missions have identified an additional 13 transiting circumbinary planets in 11 binary star systems (e.g.Martin 2018;Kostov et al. 2020;Socia et al. 2020;Kostov et al. 2021).While the radial velocity method is a highly effective way to detect exoplanets since Mayor & Queloz (1995), only two circumbinary planets have been detected using this method: Kepler-16 b (Triaud et al. 2022) and TOI-1338/BEBOP-1 c (Standing et al. 2023).
Ground-based radial-velocity surveys of binary stars are more efficient, and less biased than transit surveys (Martin et al. 2019) and thus can construct a more insightful picture of the properties of the circumbinary planet population.The BEBOP survey (Binaries Es-★ E-mail: lalitha.sairam@ast.cam.ac.uk corted By Orbiting Planets; Martin et al. 2019) began in 2017 as a dedicated, blind, radial-velocity survey of single-lined eclipsing binary stars.The survey reaches a precision of a few metres per second (Triaud et al. 2022).However, the BEBOP survey currently remains confined to highly unequal stellar-mass pairs where the secondary's mass is usually less than 30% of the primary star (Triaud et al. 2017;Martin et al. 2019).These single-lined binaries represents only around 8% of all binaries, however, the remaining 92% of binaries are double-line binaries (Kovaleva et al. 2016).Prior to the BEBOP survey, an extensive survey of double-line binaries was carried out by the TATOOINE survey (Konacki et al. 2009(Konacki et al. , 2010)).However, TATOOINE did not yield any circumbinary planets and reported a 10 − 15 m s −1 scatter in the data, enough to hide most exoplanets.Konacki et al. (2009) had to first deconvolve their spectra to remove an iodine-cell absorption spectrum, then they use TOD-COR (Zucker & Mazeh 1994) to obtain a guess radial-velocity for each component of the binary.Then they perform their tomographic deconvolution method to accurately measure the radial-velocity for each of the stellar components.
Radial velocities are also often extracted by cross-correlating the observed spectrum with a template spectrum, or with a line-list mask (the method we use; e.g.Baranne et al. 1996).When the template matches the observed spectrum a strong signal is recorded.The crosscorrelation method is convenient to understand what might help address radial velocity scatter.
In addition to the large cross-correlation function (CCF) signal, weaker signals on either side of the main signal are also recorded.They are caused by coincidental correspondence between the template and the spectrum, which we call wiggles.The wiggles exhibit pseudo-static behaviour over time in relation to the main CCF signal.We refer to these wiggles as pseudo-static because their characteristics may vary depending on the observation conditions.As an example, Figure 1 depicts a CCF time-series showing the strong stellar signal for HD 189733 but also its weaker, and stable wiggle signals.
In the context of double-line binary stars, the presence of two bright stars leads to two distinct, and strong CCF signals that vary in time based on the respective masses and orbital parameters of the binary system.Figure 2 illustrates this time-varying signal, as well as the accompanying wiggles.Because there are two stars, there are two sets of wiggles.Each can interact with the strong CCF signal of the other star, and blending of a wiggle from star A with the CCF of star B, can lead to an error in estimating the radial-velocity of star B (and vice-versa).This effect also happens on the spectral side, where weak lines from star A can blend with strong lines of star B. If the template used for deconvolution is incomplete, this will cause a similar issue when estimating radial-velocities.This phenomenon, which we refer to as the double-lined binary problem, is likely the issue that prevents the discovery of circumbinary planets.
With this effect in mind, it is obvious that a data-driven approach needs to be taken since no template can properly reproduce every spectral feature.Any error on the spectrum/template spectrum of one star, can affect the radial-velocity measurement of the other star.
In this paper, we present two new methods to derive precise radial velocities of double-lined binaries.We treat the wiggles as a correlated signal.Such correlated signals are typically and accurately treated in a data-driven way using Gaussian processes (Aigrain et al. 2016).This method has had several successes (e.g.Czekala et al. 2017;Rajpaul et al. 2020) within the exoplanet field.
The structure of the paper is as follow: §2 we introduce the Gaussian process (GP) framework.We then describe the two new techniques we developed for deriving radial velocities of double-lined binary stars.In §3, we choose to apply our two methods to the double-lined binary system TIC 172900988 because of the presence of a planet in the system that allows to test for the recovery of a Keplerian signal.In §4, we describe how we analyse the resulting radial velocities to model a binary's Keplerian motion, search for a circumbinary planet and infer its orbital parameters.Finally, we compare the radial velocities derived from our two new methods, and compare them to traditional and publicly available methods of measuring radial velocities (TODCOR and TODMOR;Zucker & Mazeh 1994;Mazeh & Zucker 1994;Zucker et al. 2004).We also discuss the implication of detecting a circumbinary planet in TIC 172900988 binary system.We conclude the paper with a summary of our key results in §5.

GAUSSIAN PROCESS-BASED RADIAL VELOCITY EXTRACTION
In this section, we provide a brief overview of the GP framework, which is a non-parametric Bayesian modelling technique that we use to infer the spectra of double-lined binary star systems.A GP is a type of stochastic process that describes the distribution of a group of random variables.It can be thought of as an extension of kernel regression to probabilistic models.Using non-parametric GPs to model the unknown wiggle function of double-lined binaries is a powerful and flexible approach (Rasmussen & Williams 2006).
A normal distribution is often represented as N (,  2 ).If a random variable  is normally distributed with mean and variance, it can be expressed as with   represents mean vector and Σ   represents the covariance matrix.The covariance matrix describes the pairwise covariance between the different elements of the input data .
The likelihood that a set of observations  is drawn from GP can be written as ln L = ln (|, , ) (2) where  and  are hyperparameters of the mean and covariance functions.In Eq. 3,  refers to the covariance matrix associated with the GP.The elements of the covariance matrix depend on the chosen covariance function and the values of the hyperparameters . represents the number of elements in the vector , which contains the observations.Evaluating this likelihood provides a posterior distribution of the hyperparameters.

Method 1: Efficient Spectral Decomposition using Gaussian Processes (SD-GP)
The observed spectra of a double-line binary star system can be modelled as a GP (e.g.Czekala et al. 2017).The spectrum of a single star can be represented as a function  () where  is the wavelength.If  > 0, the observed spectrum of a single star can be modelled as a function  () with a mean function () and a covariance kernel  (,  ′ ).If the observed spectrum has finite inputs 0 <  1 <  2 < ...... <   , then the vector [  ( 1 ),  ( 2 ), .....  (  )] has a multivariate Gaussian distribution with a mean function [( 1 ), ( 2 ), .....(  )] and a covariance matrix with  (  ,   ) as its elements, where  is the kernel function and ,  = 1, 2, ..., .Therefore, the intrinsic continuous spectrum of a star can be assumed to be a function  generated from a GP: The observed spectra of a spectroscopic binary are composed of spectral lines from both stars as a composite.Due to the orbital motion of each stars in the binary, a Doppler shift is induced to the rest frame stellar spectra which are observed in the composite spectra simultaneously.For a star moving with radial velocity v, the rest frame wavelength of the observed spectra are shifted to where  is speed of light.We describe the radial velocities of binary stars as a function of time using seven parameters: the semi-amplitude of the primary star ( A ), the binary stars mass ratio q =  B  A =  A  B , binary orbital period (), the eccentricity (), the argument of periastron (), the epoch of periastron ( 0 ) and systemic velocity ().The velocity of the primary and secondary stars as a function of time is  and respectively.
For a single epoch observation of a double-line binary star, we assume that the observed composite spectrum () is a sum of realisation  for the primary star and  for the secondary star along with , the noise process realisation.
where Σ  and Σ  are covariance matrices describing the primary and secondary star.We evaluate the covariance matrix Σ  and Σ  using the wavelengths corresponding to the primary and secondary components in their rest frame with the kernel function (Eq.11).In Figure 3 we show the realisation  and  for both components of the binary star.Similar to Rajpaul et al. (2020), we choose the Matérn kernel to model our spectra.The Matérn kernel is often used when modeling spectra because it is a flexible and versatile covariance function that can capture a wide range of smoothness and correlation properties.The Matérn covariance kernel specifies the covariance between two pixels   and   as where r   has units of km s −1 with  as the speed of light,   the characteristic length scale and   the signal standard deviation.
In practice, we analyse our spectra with the following steps and assumptions: (a) The observed composite spectra are divided into smaller wavelength subsets, which we call chunks.The radial velocity shift is computed separately for each chunk.In principle this step is not necessary, but GPs are computationally intensive and without breaking each spectrum into smaller components, the calculation become intractable.(b)We allow separate values for GP hyperparameters   and   (i.e., {  ,   }  ,{  ,   }  ) for each star of the binary system.This allows an optimal reconstruction of the spectrum taking into account the different spectral types of the stars.The mean function that is obtained by drawing samples from the GP distribution is shown in Figure 3 (bottom panel).(c) Rather than forcing a single set of hyperparameters to model a spectrum, we allow different sets of hyperparameters for different regions of a spectrum.(d)To set reasonable initial values for the radial velocities of the primary and secondary stars, we first fit simple Gaussian functions to the cross-correlation functions of each component.We then use that result to define a flat/uniform prior distribution on the radial velocities of each component.A flat prior distribution assigns equal probability density to all values within a specified range.The bounds of the flat prior distribution can be defined as the minimum and maximum values of the range, which would depend on the expected range of radial velocities for each component.We chose the bounds carefully to avoid assigning unrealistic probabilities to certain values.We then use a  2 likelihood function, to refine these estimates.(e) We simultaneously explore the posterior distribution of the radial velocities and the GP hyperparameters using an MCMC (Markov Chain Monte Carlo) sampler (the emcee python package; Goodman & Weare 2010;Foreman-Mackey et al. 2013).(f) We reiterate these steps for each chunk, and then filter bad radial velocity estimates caused by telluric absorption, stellar activity contamination or instrumental systematics.(g)The radial velocities from individual chunks are combined by computing a weighted average.The weights used for computing the weighted average of the radial velocities from individual chunks are determined by the uncertainties associated with the radial velocity estimates obtained from each chunk.(h)We normalise the observed spectra, and for the GP, we set the mean function to a constant value of 1.0 (=1).

Method 2: Cross-Correlation Functions modelled using Gaussian process (CCF-GP)
Method 1 involves dividing the observed spectrum covering a large wavelength range into smaller chunks and applying a GP regression on each chunk.This process involves computing the kernel matrix, inverting it, and multiplying it by the training set data.The computational complexity of GP regression is O ( 3 ), where  is the number of data points in each chunk.Therefore, the overall complexity of method 1 would be O (  3 ), where  is the number of chunks in the spectrum.This can be computationally very expensive.Hence we also develop an alternative method.This alternative approach is similar to the previous one but instead of modelling the entire spectrum chunk by chunk, we instead model the cross-correlated spectra, which are a typical output of instruments such as HARPS, SOPHIE and ESPRESSO (Baranne et al. 1996;Pepe et al. 2002;Perruchot et al. 2008).We assume that the crosscorrelation functions (CCFs) are samples of GP.The mean function (()) is constructed as the sum of two Gaussian functions.The covariance kernel function is used to model the correlated wiggle signal found within the cross-correlation function.
We employ a Gaussian fit jointly with a GP model for each of the components of binary.The baseline mean-function is a Gaussian function for each of the component of the binary where A 1 , B 1 , C 1 , A 2 , B 2 , C 2 are free hyper-parameters.A 1 and A 2 correspond to the amplitude of the Gaussian, which represents the contrast of primary and secondary components.B 1 and B 2 represent the radial velocities of primary and secondary stars, while C 1 and C 2 correspond to the standard deviation of the Gaussian, which represents the full width half maximum (FWHM).
We create a custom model class that inherits from celerite.modeling.Model, which is used to define the mean function of the GP.We use a Matérn covariance kernel (Equation 11), implemented in celerite (Foreman-Mackey et al. 2017) to model the correlated wiggle signal.We set bounds on the input parameters for the model and the hyperparameters of the Matérn kernel, and then creates two Matérn kernels, one for each Gaussian component, which are combined into a single kernel.
We optimise the GP model using the L-BFGS-B method (Byrd et al. 1995), which allows us to impose bounds on the parameters while minimising the negative log-likelihood of the model.We use MCMC sampling to explore the posterior distribution of the hyperparameters.For this we also use the emcee python package (Goodman & Weare 2010;Foreman-Mackey et al. 2013) with 50 walkers and a burn-in of 100 iterations.We set broad uniform priors for each hyperparameter, and run the final MCMC with 5000 iterations to converge on a solution.We take the median of the posterior distribution as the optimum solution for each hyperparameter.We then compute the 16th, 50th (median), and 84th percentiles of the posterior distribution.The uncertainties of the hyperparameters are taken as the difference between the 84 th and 50 th percentile (upper bound) and 50 th and 16 th percentile (lower bound).
The most intensive part of this method is the optimisation of GP hyperparameters which is performed by L-BFGS-B method.However, for this method the main computational bottleneck could likely be the MCMC sampling step, which has a complexity that scales with the size of the CCF.

APPLICATION TO TIC 172900988
TIC 172900988 is an eclipsing double-lined binary system consisting of two stars with spectral types F9 and G0.The orbital period of the system is approximately ∼19.7 days.Kostov et al. (2021) reported the first discovery of a circumbinary planet via what is sometimes called the "1-2 punch technique", where multiple transits occur during one conjunction event, the planet transits once over the primary and once over the secondary.Kostov et al. (2021) used a photo-dynamical analysis but did not find a single solution.The planetary radius is constrained at R p = 11.25 ± 0.44 R ⊕ .The planetary masses are proposed within a range of 823 <  p < 981 m ⊕ , and the orbital period within 188 < P p < 204 d.
We have collected radial-velocities with SOPHIE on 10 doubleline binaries (including six from Konacki et al. 2009Konacki et al. , 2010) ) in order to test our methods to different spectral types, orbital solutions, relative velocities, etc, and observe their limitations.Given the presence of a planet within the TIC 172900988 system, it serves as a good first testbed to demonstrate our ability to extract radial-velocities without removing a Keplerian signal.By focusing first on a known circumbi-nary planet host, we can evaluate the effectiveness and accuracy of our approaches.We plan to follow this paper with another paper analysing the rest of the sample that will show how precise our new methods are.

Observation
We collected 62 epochs of high-resolution spectra between 2020 October 16 and 2023 May 05 using the SOPHIE spectrograph mounted on the 1.93m telescope at Observatoire de Haute Provence (OHP) in France (Perruchot et al. 2008).The spectra cover a wavelength range of 3872-6943 Å in 39 spectral orders, with a resolving power of / ≈ 75, 000.The exposure times ranged from 600 to 1800 s depending on the seeing conditions at OHP.They have an median signal to noise ratio of 32 at 5500 Å.These are signal to noise ratio for the composite spectra of the TIC 172900988 with an average flux fraction of 0.86.This corresponds to an signal to noise ratio of ∼17 and ∼ 14 at 5500 Å for the primary and the secondary, respectively.SOPHIE was designed to detect exoplanets with a long-term stability of 2 m s −1 .The observations were taken in objAB mode, where one fibre is used to observe the starlight and another fibre is used to observe the sky brightness to estimate the background contamination such as that produced by moonlight.The wavelength calibration was performed before the night starts using a Thorium-Argon lamp and a Fabry-Pérot, fed into both fibres.Additional Fabry-Pérot calibrations are obtained roughly every two hours within the night.The spectra are extracted using the SOPHIE automatic pipeline (Bouchy et al. 2009) and the resulting wavelength-calibrated spectra are correlated with a numerical binary mask to obtain the cross-correlation functions (Baranne et al. 1996;Pepe et al. 2002).We used a G2 mask for the correlation.

Method 1 -SD-GP
We first obtain the spectra and and cross-correlate them using the SO-PHIE Data Reduction Software.To effectively apply Method 1, we work with two dimensional spectra (e2ds) at the instrument resolution instead of using 1D spectra (s1d), which operate at the pixel sampling level.We measure the radial velocities of both stellar components at the time of observation using a Gaussian fit to their cross-correlated spectra.Each SOPHIE spectrum covers from 3872 Å to 6943Å.We divide each observed spectrum into chunks of 5Å each, totalling to 615 chunks.For each chunk, we apply the SD-GP method to measure the radial velocities of both stars at each epoch.Using the calculated velocities and the parameters of the GP, we deconvolve the composite spectra into the individual spectrum of both individual stars for each epoch, by optimising the parameters of the model to fit the observed spectra.In Figure A1 (left panels), we give examples of this step of our analysis, where we show median of posterior predictive distribution of the predicted composite spectrum and the reconstructed spectra of each component of the binary.For better visualisation, we have arbitrarily included an offset to each spectrum.Note that the reconstructed spectrum matches the shape of the input composite spectrum.It is important to note that there may be chunks where the spectral lines are not present, possibly due to the continuum dominating the spectrum (Figure A1, right panels), resulting in large uncertainties in the RV measurements.After repeating this process for each chunk of spectra, we obtain the radial velocities for each star in the binary system at each epoch.We then apply outlier removal using a Student's-t distribution to remove any radial velocities that lie outside of the 95% confidence interval.The remaining radial velocities are assigned weights considering the associated uncertainties to estimate the weighted average radial velocities for the binary system (see §2.1).We then estimate the uncertainty of the combined radial velocity by propagating the uncertainties of each individual chunk through the weighted average.

Method 2 -CCF-GP
All 62 epochs of spectroscopic data from SOPHIE were crosscorrelated with a G2 mask.To determine the radial velocities of the primary and secondary star in the binary system TIC 172900988, we apply the CCF-GP to the resulting CCFs.To measure the radial velocities of the primary and secondary stars, we fit the cross-correlation function with a double Gaussian function, with the two Gaussians representing the primary and secondary stars, and two GPs to model each of wiggles caused by coincidental correspondence with the mask.In Figure 4 (left panel), we show the cross-correlation function time series along with the wiggles.We explore the posterior distribution of hyperparameters using MCMC sampling.We then calculate the 16 th , 50 th , and 84 th percentiles of the samples of the radial velocity, FWHM and contrast obtained from the MCMC simulation.These percentiles correspond to the lower uncertainty limit, the median value, and the upper uncertainty limit of the hyperparameters, respectively.The residuals following the subtraction of each component of the binary and the wiggles are shown in Figure 4 (right panel).

TODCOR and TODMOR
We also compute radial-velocities using traditional methods such as Two-Dimensional Correlation (TODCOR) and TODMOR (Two-Dimensional Modelling and Reconstruction).TODCOR is a method that uses a two-dimensional cross-correlation (2D-CCF) to measure the radial velocities of the primary and secondary stars in a binary system (Mazeh & Zucker 1994;Zucker & Mazeh 1994).A modern implementation of TODCOR for multi-order spectra is TODMOR (Zucker et al. 2004).TODMOR also uses a two-dimensional model of the observed spectra.To measure the radial velocities of the primary and secondary stars TODMOR compares each stellar component with a template spectrum matching their spectral type.
In the next section, we compare our two new methods to results produced by TODCOR and TODMOR.As such we apply TODCOR and TODMOR to the same observed spectra that we used for our own approaches.To apply TODCOR and TODMOR, first we need to correct the SOPHIE spectra for the instrumental blaze function, and detrend the pseudo-continuum.We then use the PHOENIX stellar model (Husser et al. 2013) to determine the best-matching theoretical spectra for the primary and secondary stars, and use these to optimise the 2D-CCF for each order of the spectra.We apply TODMOR to each SOPHIE order and determine the radial velocities of both components, discarding orders strongly influenced by telluric lines.
TODCOR and TODMOR use template spectra for primary and secondary stars and cross-correlate them to the observed composite spectrum to determine the radial velocities.Neither TODCOR nor TODMOR treat the wiggles.

Binary model
Each of the measured radial velocities obtained from method 1 (Sec.2.1 and 3.2), method 2 (Sec.2.2 and 3.3), TODCOR, and TODMOR (Sec.3.4) are independently modelled.We utilise kima, an open source software package for fitting radial velocities, to determine the physical parameters of the binary (Faria et al. 2018).Specifically, an updated kima package is used, which now includes simultaneous fitting of both components of double-lined binaries, correction for General Relativity effects, and that can fit for apsidal precession of the binary (Baycroft et al. 2023).For sampling, kima employs a diffusive nested sampling algorithm (DNest4, Brewer & Foreman-Mackey 2018).To account for stellar variability effects, a radial velocity jitter term is incorporated.Outliers are included in the procedure, and are handled by fitting with a student's t distribution.The system's derived parameters using radial velocities from both of the new approaches are provided in Table 1 (columns 1 and 2).
The precision reached thanks to method 1 and 2 means we have to take the circumbinary planet into account in order to properly compare them between one another and to TODCOR/TODMOR.As a nested sampler kima can fit for the number of orbiting objects in a system (in our case binary and planet), and naturally marginalises over all parameters, including the number of orbiting bodies and their possible orbits.We discuss the planet's parameters in section 4.3.

Comparison of radial velocities
After removing all Keplerian signals, we find that the residuals' root mean square (RMS) scatter of method 1 (SD-GP) is 39.9 m s −1 for the primary star and 50.9 m s −1 for the secondary star, while the RMS scatter of method 2 (CCF-GP) is 48.8 m s −1 for the primary star, and 72.2 m s −1 for the secondary star.
It is worth mentioning the root mean square (rms) values achieved by our new approaches are larger compared to the current state-ofthe-art, a scatter of 10 − 15 m, s −1 reported by Konacki et al. (2009).However, it is crucial to recognise that this increased scatter is likely inherent to the characteristics of the star itself.We have tested our methods on other double-lined binary systems, where we find that the scatter can reach down to photon noise.
Individual measurement uncertainties for the primary and secondary stars, measured by each method, range between 5.8 − 7.5 m s −1 for method 1 and 4.7 − 13 m s −1 for method 2. In Figure 5, we plot radial velocities measured using method 1 against radial velocities measured using method 2. We find the mean absolute difference between the two approaches to be 26.9 and 29.2 m s −1 for the primary and secondary stars, respectively.These mean differences are lower than the measured scatter, but exceed the uncertainties estimated by the GP fits, which suggests the presence of a systematic bias between them 1 .This bias could be due to various factors, such as the differences in the templates used for cross-correlation in method 2. In addition, method 1 might be more susceptible to the effects of stellar activity, which affect the accuracy of the radial velocities.
1 we tested method 1 and 2 on bright double-lined binaries from (Konacki et al. 2009(Konacki et al. , 2010)), observed with SOPHIE, and achieved accuracies of order 2 − 4 m s −1 , which will be the object of a follow-up paper.
Further analysis may be necessary to fully understand and quantify the sources of the observed differences.
In Figure A2, we show the radial velocity time series for each method, along with the binary+planet Keplerian models applied to them.We find radial velocities measured by our approaches are consistent with those measured by TODCOR and TODMOR within uncertainties.This suggests that our approaches are able to accurately measure the radial velocities of the primary and secondary stars.The distribution of residuals (Observed − Calculated; O-C) for TODCOR, TODMOR, method 1, and method 2 is presented in Figure 6 using violin plots.Each violin plot represents the distribution of velocities clustered around mean O-C values in m s −1 , and the width of each plot functions like a histogram.Our proposed methods outperforms TODCOR and TODMOR in terms of root-mean-square (RMS) scatter (Figure 6), producing an improvement of a factor ∼ 4 and ∼ 2 respectively.This indicates the effectiveness of our new approaches in measuring double-lined binary radial velocities more precisely than before.Since both our methods agree between one another, we are also confident our measurements gained in precision without compromising in accuracy.

The circumbinary planet within the TIC 172900988 system
The circumbinary planet is naturally detected and fitted by the Nested Sampler, however, we first describe a more frequentist approach as it might be closer to methods used in the binary star literature.
We initially compute a generalised Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009) for the radial velocities measured from the primary and the secondary stars, after having removed the best fitting Keplerian motion for the binary star.In Figure 7 we display the resulting periodogram for method 1 SD-GP (top panel) and method 2 CCF-GP (bottom panel).Using 10 000 bootstrap randomisation of the input data, we compute the false alarm probability (FAP) levels of 10, 1, 0.1 and 0.01%.This calculation can be done independently for the primary and secondary radial velocities.The periodogram for radial velocities using both method 1 and method 2 show excess power at  pl ∼ 151 d with a FAP = 0.005%.After subtracting the signal  pl , the periodogram has no significant peak (Figure A3).
We perform a more thorough analysis of the data using the kima analysis package which uses diffusive nested sampling (Faria et al. 2018;Baycroft et al. 2023).kima allows for Bayesian model comparison by computing the Bayes factor between a model with a binary and one planet to one with a binary but no planet from posterior samples generated by the algorithm.Using the Jeffrey's scale (Kass & Raftery 1995), a Bayes factor (BF) over 150 is considered strong evidence in favour of the more complex model (here binary+planet).Therefore we use this value as our confident detection threshold.In Figure 8, we show the phased radial velocity data with the best-fit Keplerian model for the circumbinary planet (the binary having been removed).This is done for the data from method 1 (top panel) and method 2 (bottom panel).
The version of kima we use fits for all the orbital elements of the binary, except Ω and , but  is known from the eclipsing geometry (Kostov et al. 2021).A different systemic velocity parameter is fit for each of the two stars.Keplerian models of the planet also include all orbital parameters except Ω and .Two jitter terms are also fit by kima, one for the primary and one for the secondary.To include outliers properly in our fit, we use Student's t statistics.
kima's fit of the SOPHIE data obtained on TIC 172900988 yields BF = 2 300 000 using the SD-GP method (method 1) and BF = 16 000 using the CCF-GP method (method 2).Both approaches exceed the detection threshold and imply a confident detection of a circumbinary planet.The parameters for the planet (as well as the binary) are shown in Table 1.Since TIC 172900988 is a double-lined eclipsing binary we obtain the absolute mass of each stellar component at high precision.Since the planet's orbital plane at the time of the observations is close to perpendicular to the line of sight (Kostov et al. 2021), we measure a mass that could be considered an absolute mass as well.However, it is likely the planetary orbital plane inclination has precessed, and might be out of transitability (e.g.Martin & Triaud 2014).We use the median of posteriors from kima and their 1 confidence region to produce our fit's parameters and uncertainties.We find all binary parameters to be statistically consistent with the analysis of (Kostov et al. 2021), with a few caveats.The binary period ( B ) we find is inconsistent with any of the 6 solutions proposed by (Kostov et al. 2021), but it does lie within  TableA1).The top and bottom panels correspond to our method 1 (SD-GP) and method 2 (SD-GP), respectively.
the range that these solutions cover2 .The argument of periastron ( B ) that we measure is not consistent, at first glance.However, our measurements were taken some time after the Kostov et al. (2021) paper.If we correct for the apsidal precession of the binary orbit, our value of  B is consistent with Kostov et al. (2021).The value we get for the apsidal precession rate (  B ) is also consistent with the value quoted in Kostov et al. (2021).This precession rate exceeds that expected from General Relativity, and is attributed to the third-body perturbations produced by the planet.Figure A4 shows the area of parameter-space a third body needs to have to produce and apsidal precession rate consistent with the observations.We overplot the location of the planetary parameters presented in this work and the solutions proposed by Kostov et al. (2021).All solutions can reproduce the detected precession rate well.
For the planet, we find  pl ≈ 150 d, and a mass  pl ≈ 2 M jup (≈ 600 M ⊕ ).While these are inconsistent with any of the six solutions proposed by Kostov et al. (2021), we fit the planet with a Keplerian model and report mean parameters, where Kostov et al. (2021) fit with a dynamical model and report osculating parameters.In such a dynamically complex orbit, it is difficult to compare these parameters properly.We note an additional important caveat here: We fit a static Keplerian to the planetary orbit and obtain a mean orbital period.
Other parameters such as the semi-major axis and the mass are then calculated using Kepler's law.However, due to the proximity of this orbit to the binary it is expected that non-Keplerian effects (such as apsidal precession of the planetary orbit) are present and the orbit would not conform to Kepler's law.Hence it is possible that the planet's true orbital distance and true mass are slightly larger than stated here.
Orbital parameters between both methods (SD-GP and CCF-GP) are internally consistent and are presented in Table 1, columns 1 and 2. Any posterior samples where the proposed planet crosses into the instability region (calculated using the formula in Holman & Wiegert 1999) are excluded, as described in Standing et al. (2022).We adopt the parameters in column 1 of Table 1 as our preferred solution.We chose this solution as the SD-GP method results in a higher Bayes Factor than the CCF-GP.The method used for the parameters reported in columns 3 and 4 is described in section 4.4.Table A1 shows the parameters obtained from the posterior sample with the "unstable" solutions left in.These are therefore the parameters simply consistent with the data without dynamical stability being considered.
Figure 9 shows the orbital configuration of TIC 172900988, displaying the orbits of the binary and the planet.A sub-sample of 1000 posterior samples are drawn: if a sample crosses into the instability zone it gets shown in red, otherwise in green.The parameters in Table 1 and the distributions shown in Figure 11 correspond therefore to the orbits shown in green.The 6 solutions from Kostov et al. (2021) are also shown for comparison 3 .
The kima algorithm can generate detection limits for any further signals, following the method presented in Standing et al. (2022): first all planetary detected signals are removed from the data (but the binary's orbital signature is kept), then kima is run once more and forced to fit a planetary signal (when presumably there are none left in the data).The resulting set of posterior samples corresponds to all signals that are compatible with the data, but have no statistically detectable signals.This method is an alternative to injection-recovery tests (e.g.Konacki et al. 2009Konacki et al. , 2010;;Martin et al. 2019), that allows to compute a detection limit efficiently over a large parameter space, while marginalising over all orbital elements.The detection limits for TIC 172900988 are shown in Figure A5 and reveal that the SOPHIE data analysed using our two new methods produce very similar results and that those are sensitive to planets with masses of order Jupiter's out to periods as large as 1000 d except for orbital periods around the one-year alias.
Finally, we run the same analysis on the TODCOR and TODMORproduced radial velocities.These produced BF = 0.8 and 0.7 respectively, well below the accepted detection threshold, demonstrating that our new approaches out-performed TODCOR and TODMOR.We show a comparison of all the resulting detection limits in Figure 10.The detection limits here are a little different than in Fig. A5, since the planet is not formally detected with all approaches.To allow for a proper comparison between the detection limits generated using 3 These may not be a full representation of the solutions from Kostov et al.
(2021) since they quote osculating elements and we plot them as if they were mean elements  2021), the dashed grey line is the stability limit as calculated by Holman & Wiegert (1999).The radial velocity data alone are fit (not including the 1-2 punch transit data).
the different methods we do not remove the planetary signal when calculating the detection limits on data from the SD-GP and CCF-GP methods.To ensure that the parameter space was well-covered in these cases we then force kima to fit 2 signals instead of the usual 1 signal.

Including the 1-2 punch technique
The detection of the planet, and the posteriors of its orbital and physical parameters can be improved by combining the radial velocity data with some aspects of the transit data.Ultimately, a full photodynamical analysis would need to be performed, but this is beyond the scope of our paper.
TIC 172900988 was discovered using the "1-2 punch" method (Kostov et al. 2021).Two transits within the same conjunction give an estimate of the planet's orbital period.The distance the planet has moved can be calculated from the position in its orbit of the transited star, at the time of each transit.The time between the transits and the distance travelled then allow us to calculate an estimate of the planet's orbital period (as in Kostov et al. 2020): where  bin is the total mass of the binary,  and  the eccentricity and argument of periastron of the planet,  is the true anomaly at the point in the orbit that we are measuring and  the average velocity in the plane of the sky with which the planet moved between both Table 1.Best-fit parameters for each methods both including and not including the 1-2 punch into the fit.Solutions that cross the instability limit (Holman & Wiegert 1999) are excluded here (parameters from the full posteriors can be found in Table A1).The parameters are determined along with their corresponding 1 uncertainties.We alter our version of kima and add an extra feature, to include the "1-2 punch" information as part of the the sampling process.When a solution is proposed by the sampler, the predicted period from the "1-2 punch" is calculated and compared to the proposed period.This is then included in the likelihood calculation of the sample in the same way as an extra data point would be, assuming a Gaussian distribution.Our new log-likelihood is therefore:

ADOPTED
where log(L RV ) is the log-likelihood from the radial velocity data,  pl is the period for the planet proposed as part of the sampling process,  12 is the orbital period calculated using Eq.14 using all other parameters proposed in the sample (e.g. bin , , etc).Finally, is the variance of  12 , which is derived from the uncertainty in the transit times propagated through Eq. 14.
We run the analysis on TIC 172900988 again, with the extra input of two transits with mid-times at 2 458 883.390879 ± 0.006188 and 2 458 888.309427±0.0039045 .We fix the number of planets searched for in kima to 1.The parameters for the planet and binary obtained are shown in Table 1 (columns 3 and 4).As with the previous analysis, any posterior samples where the proposed planet crosses into the instability are excluded.
A Keplerian fit of the radial velocity data probes the average parameters of the orbit over the time baseline, notably the average orbital period.A circumbinary planet, especially one like TIC 172900988 b which is quite close to the inner binary will see its orbital parameters vary throughout its orbital path meaning that when parameters of the orbit are measured over a short time frame, they may not be representative of the average orbit.The 1-2 punch, method calculates the velocity and therefore the orbital period, over a short time frame.Using this method to constrain the average period might bring in a poorly understood bias.We therefore present the results from the combined radial velocity and 1-2 punch fit out of interest, but do not adopt these parameters as our preferred solution.
The posterior samples for the planet parameters are shown as a corner plot (Foreman-Mackey 2016) in Figure 11 for both the CCF-GP method and the SD-GP method.We note that the crescentshaped correlations involving the eccentricity are expected.We also note that while we report the median and 1- in Table 1, some of the parameters have non-Gaussian distributions (in particular the eccentricity ( pl ) and argument of periastron ( pl )).
We find these new solutions are consistent with solutions of fitting just the radial velocities or just the average orbital velocity at conjunction between the transits.The combined fit suggests the circumbinary planet's orbit must have an eccentricity  pl > 0.1 and an argument of periastron 4.35 ≤  pl ≤ 5.31 rad.
Figure A6 shows the orbital configuration of TIC 172900988.This is the counterpart to Figure 9 generated from posterior samples obtained from the kima analysis which included the 1-2 punch and using the SD-GP radial velocities.

CONCLUSIONS
In this work we focus on the development of two new data-driven approaches to accurately measure radial velocities in double-lined binary systems.Despite being brighter and more precise in principle, the time-varying blending of the two stars' spectral lines makes accurate radial velocity measurement challenging.Previous methods (Konacki et al. 2009) been shown to have a typical scatter of 10 − 15 m s −1 that prevents the detection of most orbiting circumbinary planets.
In this paper, we introduce two new methods based on GP regression inspired by Czekala et al. (2017).The first method applies the GP in the spectral domain and the second is applied on cross-correlated spectra.We compare the precision and accuracy of our radial-velocity to two widely used methods: TODCOR and TODMOR (Mazeh & Zucker 1994;Zucker et al. 2004).
To perform the comparison, we analyse 62 SOPHIE spectra of the binary TIC 172900988, a binary system which was also proposed to host a circumbinary planet (Kostov et al. 2021).We show that our two methods outperform both TODCOR and TODMOR, neither of which could recover the planet whereas both our GP approaches successfully detect a circumbinary planet.However, its parameters are found statistically different from previously published solutions (Kostov et al. 2021).TIC 172900988 b will now integrate the BEBOP catalogue for circumbinary exoplanets detected with radial-velocities as its second entry.
The RMS achieved by our new approaches are ≈ 50 and 70 m s −1 for the primary and secondary stars of TIC 172900988, respectively, are larger than our measurement uncertainties, and larger than the typical scatter reported for double-lined binaries in Konacki et al.  (10 − 15 m s −1 ; 2009, 2010).We speculate this increased scatter is most likely of stellar origin.Both stellar components are fairly highmass stars.We note that TODCOR and TODMOR, both recognised methods of radial velocity extraction, also produce a high scatter.In those cases, TODCOR and TODMOR remain unable to detect the planetary signal which has a semi-amplitude  pl ∼ 42 m s −1 .The fact our approaches both manage to overcome some of that scatter emphasises the limitations of existing techniques when dealing with systems characterised by substantial scatter.The detection of a circumbinary planet in TIC 172900988 showcases the effectiveness of our data-driven methods in uncovering planetary signals even in challenging double-lined binary systems.We highlight here that should a circumbinary planet similar to the parameters of TIC 172900988 b have been present in a quieter binary star system, traditional methods such as TODCOR, TODMOR and the tomographic disentangling method have the nominal accuracy to detect it.
Our two new methods are a step forward, but there is always room for improvement.Further refinements and optimisations to these methods may lead to even more accurate and precise radial velocity measurements, particularly with the spectral decomposition.
Firstly, we recognise that our analysis does not account for possible contamination in the radial velocities obtained from each chunk of the spectra.Specifically, we do not consider the effects of stellar activity on our results.Also the chunking can be improved to avoid areas that are poor in absorption lines, and avoid areas that include band known to be highly variable such as H.In order to improve the accuracy and precision of our measurements, we plan to develop a more sophisticated approach that can identify regions of the spectra that are affected by these factors.
Secondly, we recognise that there is still much to be learned about the astrophysical properties of the binary stars themselves.In the case of TIC 172900988, both stars are equal mass.The CCF and spectral decomposition methods might need to be adapted to non-equal mass binaries to account for the their differing spectral types.
In addition, we expect that the spectral deconvolution method will yield accurate spectra for both stars individually when all wavelength chunks are combined together.Such a spectrum could be used to constrain their properties such as their temperature,  sin  and metallicity, an important parameter to relate planet presence to planet formation (e.g.Santos et al. 2004;Adibekyan et al. 2013).
Our methods open the door to extend the search for circumbinary planets using the radial-velocity method beyond single-lined binaries.With our two new approaches, it is highly probable that the discovery of circumbinary planets will be enhanced in the future.Finally, we highlight the success of our two new methods in being the first to detect a circumbinary planet using radial-velocities in a double-lined binary.Importantly this detection is made independently of any other data.Interestingly our results produce planetary parameters different from those previously published demonstrating the need for radialvelocity follow-up of circumbinary planet candidates identified with the transit method.   .SOPHIE radial velocities as a function of time and the corresponding residuals for the method 1 (SD-GP), method 2 (CCF-GP), TODCOR and TODMOR (top to bottom).The overplotted magenta and orange curves are the medial value of the posterior distribution models for the primary and the secondary stars corresponding to 2,000 random MCMC steps.

Figure 1 .
Figure 1.Left panel: A time series of cross correlation function for HD 189733.The bright profile shows the stellar profile.The wiggles appear on either side of bright profile.The black profile is the median CCF of all epochs.Right panel: Residual map of time series on removal of stellar profile.Wiggles are seen as alternating static bands.Overlaid in black is the median wiggle of all epochs.

Figure 2 .
Figure 2. Schematic diagram showing the double-line binary problem and our solution to the problem using method 2 (CCF-GP)

Figure 3 .
Figure 3. Depicted in the top two panels are the multivariate realisation vector from prior distribution draw for GP hyperparameters of the primary and the secondary stars.The realisations are not constrained by the data.The bottom panel depicts the mean prediction (in green) from the posterior predictive distribution of GP drawn by conditioning on the observed data from a single epoch observation.The orange and blue in the bottom panel represents the primary and the secondary star with an arbitrary offset.The bottom panel also includes the residuals.

Figure 4 .
Figure 4. Left panel: A time-series of cross correlation for double-line binary star showing the strong signals for each components.The wiggles appear as dark bands on either sides of the strong signal.Right panel: The wiggles as modelled by method 2 (CCF-GP).

Figure 5 .
Figure5.Top panel: Radial velocities extracted using method 1 versus the radial velocities from method 2, along with 1 error bars.The red and blue colours represent the primary and the secondary stars along with their respective 1:1 identity lines.Bottom panel: The residual radial velocities after removing the binary signal for method 1 versus the residual radial velocities for method 2.

Figure 6 .
Figure6.The distribution of O-C values for the primary (red) and secondary stars (blue) in TIC172900988, as determined by TODCOR, TODMOR, method 1 (SD-GP), and method 2 (GP-CCF), are displayed using violin plots.Note that the vertical scale is different on every panel.

Figure 7 .
Figure 7. Lomb-Scargle periodogram of TIC172900988 radial velocities for method 1 (top) and method 2 (bottom).The radial velocities for the primary (red) and the secondary (blue) are plotted after removing the binary motion.The three horizontal dashed lines indicate 10%, 1% and 0.1% false alarm probabilities.The vertical dotted lines indicates the highly significant peak around 150 days and its harmonics.

Figure 8 .
Figure8.The residual RVs phase folded with corresponding best fit Keplerian circumbinary model (TableA1).The top and bottom panels correspond to our method 1 (SD-GP) and method 2 (SD-GP), respectively.

Figure 9 .
Figure 9. Orbital configuration of TIC 172900988 showing the orbits of the binary and the planet.The red orbits are a random sample of 1000 posteriors from kima fitting the radial velocities from SD-GP.The green orbits are the 6 suggested solutions from Kostov et al. (2021), the dashed grey line is the stability limit as calculated byHolman & Wiegert (1999).The radial velocity data alone are fit (not including the 1-2 punch transit data).

Figure 10 .
Figure 10.Detection sensitivity to planets plotted as semi-amplitude as a function of orbital period of planets.The density of posterior samples are depicted as grey hexagonal bins.The solid green, red, orange and blue lines shows the detection limit from posterior samples for TODCOR, TODMOR, method 1 and method 2, respectively.The diagonal lines are anticipated signals of Saturn and Jupiter mass planet.

Figure 11 .
Figure11.Corner plot showing the distributions and correlations of the planetary parameters from the simultaneous fit of the radial velocity data and the 1-2 punch transit times.The contours are the 50th and 90th percentiles.Orange shows the results using the the radial velocity data using SD-GP, and blue from the CCF-GP.

Figure A1 .
Figure A1.Seven epochs of composite spectra for two chunks between 5379-5384Å and 5437-5443 Å are shown in grey with composite model in green.The mean realisation drawn from posterior predictive distribution are show in red and blue for the primary and secondary, respectively.The residuals are shown at the bottom of each panel.The first chunk shows several lines and the second chunk shows no lines in the composite spectra.The redial velocities from such chunks are weighted.
Figure A2.SOPHIE radial velocities as a function of time and the corresponding residuals for the method 1 (SD-GP), method 2 (CCF-GP), TODCOR and TODMOR (top to bottom).The overplotted magenta and orange curves are the medial value of the posterior distribution models for the primary and the secondary stars corresponding to 2,000 random MCMC steps.

Figure A3 .
Figure A3.The residual periodgram after removing the signal P  for both method 1 (top panel) and method 2 (bottom panel).The red and blue curves represent the primary and the secondary components of the binary.The grey vertical line indicates the period of  b .