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

Lalitha Sairam, Amaury Triaud, Thomas Baycroft, Jerome Orosz, Isabelle Boisse, Neda Heidari, Daniel Sebastian, Georgina Dransfield, David Martin, Alexandre Santerne, et al.

► To cite this version:

Lalitha Sairam, Amaury Triaud, Thomas Baycroft, Jerome Orosz, Isabelle Boisse, et al. New methods for radial-velocity measurements of double-lined binaries, and detection of a circumbinary planet orbiting TIC 172900988. Monthly Notices of the Royal Astronomical Society, 2024, 527 (2), pp.2261-2278. 10.1093/mnras/stad3136. hal-04505685

HAL Id: hal-04505685

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

Submitted on 16 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# New methods for radial-velocity measurements of double-lined binaries, and detection of a circumbinary planet orbiting TIC 172900988

Lalitha Sairam<sup>D, 1,2*</sup> Amaury H. M. J. Triaud<sup>D, 1*</sup> Thomas A. Baycroft<sup>D, 1*</sup> Jerome Orosz,<sup>3</sup>

Isabelle Boisse,<sup>4</sup> Neda Heidari,<sup>5</sup> Daniel Sebastian<sup>D, 1</sup>,<sup>1</sup> Georgina Dransfield<sup>D, 1</sup>,<sup>1</sup> David V. Martin<sup>D, 6</sup>

Alexandre Santerne<sup>D, 4</sup> and Matthew R. Standing<sup>D, 7</sup>

Accepted 2023 October 10. Received 2023 September 7; in original form 2023 May 31

## ABSTRACT

Ongoing ground-based radial-velocity observations seeking to detect circumbinary planets focus on single-lined binaries even though over 9 in every 10 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–15 m s⁻<fnref n="1" />. 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 d, 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.

**Key words:** techniques: radial velocities – planets and satellites: detection – planets and satellites: gaseous planets – planets and satellites: individual: TIC172900988 – binaries: eclipsing – binaries: spectroscopic.

## 1 INTRODUCTION

Planets orbiting around both stars of a binary star system are called *p*-type or *circumbinary* planets. The discovery of Kepler-16 (Doyle et al. 2011) marked the first confirmation of the existence of these long-theorized 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 Escorted 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 per cent of the primary star (Triaud et al. 2017; Martin et al. 2019). These single-lined binaries represents only around 8 per cent of all binaries; however, the remaining 92 per cent 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, 2010). However, TATOOINE did not yield any circumbinary planets and reported a 10–15 m s⁻<fnref n="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 TODCOR (Zucker & Mazeh 1994) to obtain a guess radial velocity for each component of the

<footnote n="1"><sup>1</sup>School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

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

<footnote n="3"><sup>3</sup>Department of Astronomy, San Diego State University, 5500 Campanile Drive, San Diego, CA 92182, USA</footnote>

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

<footnote n="5"><sup>5</sup>Institut d'astrophysique de Paris, UMR 7095 CNRS Université Pierre et Marie Curie, 98 bis, Boulevard Arago, F-75014 Paris, France</footnote>

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

<footnote n="7"><sup>7</sup>School of Physical Sciences, The Open University, Milton Keynes, MK7 6AA, UK</footnote>

<footnote>*E-mail: lalitha.sairam@ast.cam.ac.uk (LS); a.triaud@bham.ac.uk (AHMJT); txb187@student.bham.ac.uk (TAB)</footnote>
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 cross-correlation 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, Fig. 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. Fig. 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 (GP; Aigrain, Parviainen & Pope 2016). This method has had several successes (e.g. Czekala et al. 2017; Rajpaul, Aigrain & Buchhave 2020) within the exoplanet field.

The structure of the paper is as follow: Section 2 we introduce the GP framework. We then describe the two new techniques we developed for deriving radial velocities of double-lined binary stars. In Section 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 Section 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; Mazeh & Zucker 1994; Zucker & Mazeh 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 Section 5.

## 2 GP-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(\mu, \sigma^2)$. If a random variable $x$ is normally distributed with mean and variance, it can be expressed as

$$x \sim N(\mu_x, \Sigma_{xx}), \quad (1)$$

with $\mu_x$ represents mean vector and $\Sigma_{xx}$ represents the covariance matrix. The covariance matrix describes the pairwise covariance between the different elements of the input data $x$.

The likelihood that a set of observations $y$ is drawn from GP can be written as

$$\ln \mathcal{L} = \ln p(y|x, \phi, \theta), \quad (2)$$

$$= -\frac{1}{2}(y - \mu)^{\mathrm{T}} K^{-1}(y - \mu) - \frac{1}{2} \ln |K| - \frac{N}{2} \ln(2\pi), \quad (3)$$

where $\phi$ and $\theta$ are hyperparameters of the mean and covariance functions. In equation (3), $K$ 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 $\theta$. $N$ represents the number of elements in the vector $y$, which contains the observations. Evaluating this likelihood provides a posterior distribution of the hyperparameters.

### 2.1 Method 1: efficient spectral decomposition using GP (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 $f(\lambda)$, where $\lambda$ is the wavelength. If $\lambda > 0$, the observed spectrum of a single star can be modelled as a function $f(\lambda)$ with a mean function $\mu(\lambda)$ and a covariance kernel $k(\lambda, \lambda')$. If the observed spectrum has finite inputs $0 < \lambda_1 < \lambda_2 < \dots < \lambda_w$, then the vector $[f(\lambda_1), f(\lambda_2), \dots, f(\lambda_w)]$ has a multivariate Gaussian distribution with a mean function $[\mu(\lambda_1), \mu(\lambda_2), \dots, \mu(\lambda_w)]$ and a covariance matrix with $k(\lambda_i, \lambda_j)$ as its elements, where $k$ is the kernel function and $i, j = 1, 2, \dots, w$. Therefore, the intrinsic continuous spectrum of a star can be assumed to be a function $f$ generated from a GP:

$$f(\lambda) \sim \mathrm{GP}(\mu(\lambda), k(\lambda, \lambda')). \quad (4)$$

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 that 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

$$\lambda(v) = \left(\frac{c+v}{c-v}\right)^{\frac{1}{2}} \lambda_0, \quad (5)$$

where $c$ 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
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Left panel: A time-series of CCF 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.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic diagram showing the double-line binary problem and our solution to the problem using method 2 (CCF-GP).</figcaption></figure>

$(K_A)$, the binary stars mass ratio ($q = \frac{M_B}{M_A} = \frac{K_A}{K_B}$), binary orbital period ($P$), the eccentricity ($e$), the argument of periastron ($\omega$), the epoch of periastron ($T_0$), and systemic velocity ($\gamma$). The velocity of

the primary and secondary stars as a function of time is

$$v_A = K_A[\cos(\omega + f(t)) + e \cos \omega] + \gamma \quad (6)$$

and

$$v_B = -\frac{K_A}{q}[\cos(\omega + f(t)) + e \cos \omega] + \gamma, \quad (7)$$

respectively.

For a single-epoch observation of a double-line binary star, we assume that the observed composite spectrum ($s$) is a sum of realization $f$ for the primary star and $g$ for the secondary star along with $N$, the noise process realization:

$$s = f + g + N, \quad (8)$$

$$\sim \mathcal{N}(\mu_f, \Sigma_f) + \mathcal{N}(\mu_g, \Sigma_g) + \mathcal{N}(0, \Sigma_N), \quad (9)$$

$$\sim \mathcal{N}(\mu_f + \mu_g, \Sigma_f + \Sigma_g + \Sigma_N), \quad (10)$$

where $\Sigma_f$ and $\Sigma_g$ are covariance matrices describing the primary and secondary star. We evaluate the covariance matrix $\Sigma_f$ and $\Sigma_g$ using the wavelengths corresponding to the primary and secondary components in their rest frame with the kernel function (equation 11). In Fig. 3, we show the realization $f$ and $g$ for both components of the binary star.

Similar to Rajpaul, Aigrain & Buchhave (2020), we choose the Matérn kernel to model our spectra. The Matérn kernel is often used when modelling 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 $\lambda_i$ and $\lambda_j$ as

$$k_{ij}(r_{ij}|\theta) = \sigma_f^2 \left(1 + \frac{\sqrt{3}r}{\sigma_l}\right) \exp\left(-\frac{\sqrt{3}r}{\sigma_l}\right), \quad (11)$$

where $r_{ij}$ has units of km s⁻¹

$$r_{ij} = r(\lambda_i, \lambda_j) = \frac{c}{2} \left| \frac{\lambda_i - \lambda_j}{\lambda_i + \lambda_j} \right|, \quad (12)$$

with $c$ as the speed of light, $\sigma_l$ is the characteristic length-scale, and $\sigma_f$ is the signal standard deviation.
<figure><img src="image_4.png" /><figcaption>Figure 3. Depicted in the top two panels are the multivariate realization vector from prior distribution draw for GP hyperparameters of the primary and the secondary stars. The realizations 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.</figcaption></figure>

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 σf and σl (i.e. {σf, σl}f, {σf, σl}l) 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 Fig. 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 CCFs 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 χ² likelihood function to refine these estimates.

(e) We simultaneously explore the posterior distribution of the radial velocities and the GP hyperparameters using a Markov chain Monte Carlo (MCMC) 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 normalize the observed spectra, and for the GP, we set the mean function to a constant value of 1.0 (μ = 1).

## 2.2 Method 2: CCFs modelled using GP (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(N³), where N is the number of data points in each chunk. Therefore, the overall complexity of method 1 would be O(MN³), where M 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 CCFs are samples of GP. The mean function (μ(x)) is constructed as the sum of two Gaussian functions. The covariance kernel function is used to model the correlated wiggle signal found within the CCF.

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

$$ \mu(x) = 1 - A_1 \exp \left( - \frac{(x - B_1)^2}{2C_1^2} \right) - A_2 \exp \left( - \frac{(x - B_2)^2}{2C_2^2} \right), \quad (13) $$

where A₁, B₁, C₁, A₂, B₂, C₂ are free hyper-parameters. A₁ and A₂ correspond to the amplitude of the Gaussian, which represents the contrast of primary and secondary components. B₁ and B₂ represent the radial velocities of primary and secondary stars, while C₁ and C₂ correspond to the standard deviation of the Gaussian, which represents the full width at 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 optimize the GP model using the L-BFGS-B method (Byrd et al. 1995), which allows us to impose bounds on the parameters
while minimizing 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 84th and 50th percentile (upper bound) and 50th and 16th percentile (lower bound).

The most intensive part of this method is the optimization 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.

## 3 APPLICATION TO TIC 172900988

TIC 172 900 988 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 d. 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 photodynamical analysis but did not find a single solution. The planetary radius is constrained at $R_p = 11.25 \pm 0.44 R_⊕$. The planetary masses are proposed within a range of $823 < m_p < 981 R_⊕$, and the orbital period within $188 < P_p < 204$ d.

We have collected radial-velocities with SOPHIE on 10 double-line binaries (including six from Konacki et al. 2009, 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 172 900 988 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 circumbinary 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.

### 3.1 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.93 m 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 $\delta\lambda/\lambda \approx 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 (SNR) of 32 at 5500 Å. These are SNR for the composite spectra of the TIC 172 900 988 with an average flux fraction of 0.86. This corresponds to an SNR 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 \text{ 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 2 h 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 CCFs (Baranne et al. 1996; Pepe et al. 2002). We used a G2 mask for the correlation.

### 3.2 Method 1 – SD-GP

We first obtain the spectra and cross-correlate them using the SOPHIE 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 optimizing the parameters of the model to fit the observed spectra. In Fig. 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 visualization, 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 (Fig. A1, right panels), resulting in large uncertainties in the radial velocity (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 per cent 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 Section 2.1). We then estimate the uncertainty of the combined radial velocity by propagating the uncertainties of each individual chunk through the weighted average.

### 3.3 Method 2 – CCF-GP

All 62 epochs of spectroscopic data from SOPHIE were cross-correlated 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 CCF 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 Fig. 4 (left panel), we show the CCF time-series along with the wiggles. We explore the posterior distribution of hyperparameters using MCMC sampling. We then calculate the 16th, 50th, and 84th 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 Fig. 4 (right panel).
<figure><img src="image_5.png" /><figcaption>**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).</figcaption></figure>

## 3.4 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 function (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 multiorder 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 optimize 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.

# 4 RESULTS

## 4.1 Binary model

Each of the measured radial velocities obtained from method 1 (Sections 2.1 and 3.2), method 2 (Sections 2.2 and 3.3), TODCOR, and TODMOR (Section 3.4) are independently modelled. We utilize 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 marginalizes over all parameters, including the number of orbiting bodies and their possible orbits. We discuss the planet's parameters in Section 4.3.

## 4.2 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⁻¹ for the primary star and 50.9 m s⁻¹ for the secondary star, while the RMS scatter of method 2 (CCF-GP) is 48.8 m s⁻¹ for the primary star, and 72.2 m s⁻¹ for the secondary star.

It is worth mentioning the RMS values achieved by our new approaches are larger compared to the current state of the art, a scatter of 10–15 m s⁻¹ 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 and 7.5 m s⁻¹ for method 1 and 4.7 and 13 m s⁻¹ for method 2. In Fig. 5, we plot radial velocities measured using method 1 against radial velocities measured using method 2. We find the
<table><caption>**Table 1.** Best-fitting 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\sigma$ uncertainties.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="4">ADOPTED</th></tr><tr><th>Method 1<br>SD-GP</th><th>Method 2<br>CCF-GP</th><th>Method 1 (with 1–2 punch)<br>SD-GP</th><th>Method 2 (with 1–2 punch)<br>CCF-GP</th></tr></thead><tbody><tr><td colspan="5"><strong>Binary parameters</strong></td></tr><tr><td>P<sub>B</sub> (d)</td><td>19.657878<sup>+0.000029</sup><sub>-0.000034</sub></td><td>19.657861 ± 0.000040</td><td>19.657874<sup>+0.000030</sup><sub>-0.000037</sub></td><td>19.657869<sup>+0.000048</sup><sub>-0.000051</sub></td></tr><tr><td>e<sub>B</sub></td><td>0.448234<sup>+0.000090</sup><sub>-0.000071</sub></td><td>0.44823 ± 0.00010</td><td>0.448230<sup>+0.000081</sup><sub>-0.000083</sub></td><td>0.44819 ± 0.00013</td></tr><tr><td>ω<sub>B</sub> (rad)</td><td>1.23147<sup>+0.00019</sup><sub>-0.00027</sub></td><td>1.23136<sup>+0.00035</sup><sub>-0.00038</sub></td><td>1.23142<sup>+0.00022</sup><sub>-0.00024</sub></td><td>1.23112<sup>+0.00049</sup><sub>-0.00036</sub></td></tr><tr><td>K<sub>B</sub> (km s<sup>-1</sup>)</td><td>58.5494<sup>+0.0079</sup><sub>-0.0054</sub></td><td>58.5279<sup>+0.0076</sup><sub>-0.0074</sub></td><td>58.580 ± 0.011</td><td>58.5428<sup>+0.0104</sup><sub>-0.0092</sub></td></tr><tr><td>q<sub>B</sub></td><td>0.97196<sup>+0.00016</sup><sub>-0.00022</sub></td><td>0.97106<sup>+0.00025</sup><sub>-0.00023</sub></td><td>0.97211 ± 0.00022</td><td>0.97121<sup>+0.00029</sup><sub>-0.00027</sub></td></tr><tr><td>ω̇<sub>B</sub>(arcsec yr<sup>-1</sup>)</td><td>215<sup>+56</sup><sub>-59</sub></td><td>279<sup>+72</sup><sub>-70</sub></td><td>229<sup>+62</sup><sub>-71</sub></td><td>250<sup>+120</sup><sub>-110</sub></td></tr><tr><td>T<sub>0, B</sub> (BJD)</td><td>2 459 566.00250<sup>+0.00044</sup><sub>-0.00061</sub></td><td>2 459 566.00210<sup>+0.00079</sup><sub>-0.00088</sub></td><td>2 459 566.00267<sup>+0.00053</sup><sub>-0.00062</sub></td><td>2 459 507.02830<sup>+0.00113</sup><sub>-0.00094</sub></td></tr><tr><td colspan="5"><strong>Planet parameters</strong></td></tr><tr><td>P<sub>pl</sub> (d)</td><td>151.2 ± 1.8</td><td>149.3 ± 2.2</td><td>151.4 ± 1.7</td><td>149.9<sup>+2.4</sup><sub>-2.2</sub></td></tr><tr><td>e<sub>pl</sub></td><td>&lt;0.11</td><td>&lt;0.11</td><td>0.1243<sup>+0.0183</sup><sub>-0.0090</sub></td><td>0.1251<sup>+0.0124</sup><sub>-0.0083</sub></td></tr><tr><td>ω<sub>pl</sub> (rad)</td><td>6.0 ± 1.8</td><td>4.9 ± 2.0</td><td>5.00<sup>+0.31</sup><sub>-0.49</sub></td><td>4.76 ± 0.41</td></tr><tr><td>K<sub>pl</sub> (m s<sup>-1</sup>)</td><td>40.1<sup>+5.1</sup><sub>-5.3</sub></td><td>44.0 ± 6.8</td><td>39.1<sup>+4.5</sup><sub>-4.3</sub></td><td>40.5<sup>+6.3</sup><sub>-6.6</sub></td></tr><tr><td>T<sub>0, pl</sub> (BJD)</td><td>2 459 442<sup>+40</sup><sub>-45</sub></td><td>2 459 425<sup>+45</sup><sub>-49</sub></td><td>2 459 419.7<sup>+7.8</sup><sub>-12.0</sub></td><td>2 459 421 ± 10</td></tr><tr><td colspan="5"><strong>Derived parameters</strong></td></tr><tr><td>M<sub>A</sub> (M<sub>☉</sub>)</td><td>1.23681<sup>+0.00037</sup><sub>-0.00039</sub></td><td>1.23765<sup>+0.00051</sup><sub>-0.00046</sub></td><td>1.23670 ± 0.00036</td><td>1.23831<sup>+0.00060</sup><sub>-0.00063</sub></td></tr><tr><td>M<sub>B</sub> (M<sub>☉</sub>)</td><td>1.20207<sup>+0.00033</sup><sub>-0.00026</sub></td><td>1.20184<sup>+0.00035</sup><sub>-0.00031</sub></td><td>1.20221<sup>+0.00032</sup><sub>-0.00028</sub></td><td>1.20264<sup>+0.00048</sup><sub>-0.00046</sub></td></tr><tr><td>m<sub>pl</sub> (M<sub>Jup</sub>)</td><td>1.90 ± 0.25</td><td>2.07 ± 0.32</td><td>1.84 ± 0.21</td><td>2.09 ± 0.27</td></tr><tr><td>a<sub>bin</sub> (au)</td><td>0.191879<sup>+0.000016</sup><sub>-0.000018</sub></td><td>0.191894 ± 0.000021</td><td>0.191879 ± 0.000016</td><td>0.191931<sup>+0.000026</sup><sub>-0.000028</sub></td></tr><tr><td>a<sub>pl</sub> (au)</td><td>0.7474<sup>+0.0057</sup><sub>-0.0061</sub></td><td>0.7410<sup>+0.0069</sup><sub>-0.0075</sub></td><td>0.7479<sup>+0.0051</sup><sub>-0.0056</sub></td><td>0.7433<sup>+0.0079</sup><sub>-0.0072</sub></td></tr><tr><td colspan="5"><strong>Other parameters</strong></td></tr><tr><td>V<sub>sys, A</sub> (km s<sup>-1</sup>)</td><td>25.9865<sup>+0.0045</sup><sub>-0.0047</sub></td><td>26.0252 ± 0.0057</td><td>26.0144 ± 0.0045</td><td>26.0573 ± 0.0067</td></tr><tr><td>V<sub>sys, B</sub> (km s<sup>-1</sup>)</td><td>26.0456<sup>+0.0062</sup><sub>-0.0060</sub></td><td>26.0969<sup>+0.0087</sup><sub>-0.0089</sub></td><td>26.0636<sup>+0.0059</sup><sub>-0.0061</sub></td><td>26.116 ± 0.010</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Top panel: Radial velocities extracted using method 1 versus the radial velocities from method 2, along with $1\sigma$ 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.</figcaption></figure>

mean absolute difference between the two approaches to be 26.9 and 29.2 m s⁻<fnref n="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.<fnref n="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. Further analysis may be necessary to fully understand and quantify the sources of the observed differences.

In Fig. 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 Fig. 6 using violin plots. Each violin plot represents the distribution of velocities clustered around mean O–C values in m s⁻<fnref n="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 (Fig. 6), producing an improvement of a factor of ~4

<footnote n="1">¹We tested method 1 and 2 on bright double-lined binaries from Konacki et al. (2009, 2010), observed with SOPHIE, and achieved accuracies of order 2–4 m s⁻¹, which will be the object of a follow-up paper.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** 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 (CCF-GP), are displayed using violin plots. Note that the vertical scale is different on every panel.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**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 per cent, 1 per cent and 0.1 per cent false alarm probabilities. The vertical dotted lines indicate the highly significant peak around 150 d and its harmonics.</figcaption></figure>

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.

## 4.3 The circumbinary planet within the TIC 172 900 988 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 generalized 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 Fig. 7, we display the resulting periodogram for method 1 SD-GP (top panel) and method 2 CCF-GP (bottom panel). Using 10 000

bootstrap randomization of the input data, we compute the false alarm probability (FAP) levels of 10, 1, 0.1, and 0.01 per cent. 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 $P_{pl} \sim 151$ d with a FAP = 0.005 per cent. After subtracting the signal $P_{pl}$, the periodogram has no significant peak (Fig. 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 Fig. 8, we show the phased radial velocity data with the best-fitting 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 $\Omega$ and $i$, but $i$ 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 $\Omega$ and $i$. 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 172 900 988 yields $BF = 2300 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 172 900 988 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\sigma$ 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 ($P_B$) we find is inconsistent with any
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The residual RVs phase folded with corresponding best-fitting Keplerian circumbinary model (Table A1). The top and bottom panels correspond to our method 1 (SD-GP) and method 2 (CCF-GP), respectively.</figcaption></figure>

of the six solutions proposed by Kostov et al. (2021), but it does lie within the range that these solutions cover.<fnref n="2" /> The argument of periastron ($\omega_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 $\omega_B$ is consistent with Kostov et al. (2021). The value we get for the apsidal precession rate ($\dot{\omega}_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. Fig. 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.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Orbital configuration of TIC 172900988 showing the orbits of the binary and the planet. The green orbits are a random sample of 1000 posteriors from *kima* fitting the radial velocities from SD-GP. The blue orbits are the six suggested solutions from Kostov et al. (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).</figcaption></figure>

For the planet, we find $P_{pl} \approx 150$ d and a mass $m_{pl} \approx 2$ M$_{jup}$ ($\approx 600$ M$_{\oplus}$). 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 semimajor 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.

Fig. 9 shows the orbital configuration of TIC 172900988, displaying the orbits of the binary and the planet. A subsample 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 Fig. 11 correspond therefore

<footnote n="2">²These six solutions consist of osculating elements and the binary periods are internally inconsistent with each other.</footnote>
<figure><img src="image_11.png" /><figcaption>**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 show 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.</figcaption></figure>

to the orbits shown in green. The six solutions from Kostov et al. (2021) are also shown for comparison.<fnref n="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. 2009, 2010; Martin et al. 2019) that allow to compute a detection limit efficiently over a large parameter space, while marginalizing over all orbital elements. The detection limits for TIC 172 900 988 are shown in Fig. 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 1 yr alias.

Finally, we run the same analysis on the TODCOR and TODMOR-produced 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 Fig. 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 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 two signals instead of the usual one signal.

## 4.4 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 172 900 988 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):

$$P_{\text{pl}} = \frac{2\pi G M_{\text{bin}}}{v^3} \left( \frac{e \sin \omega + \sin(\phi + \omega)}{\sqrt{1 - e^2}} \right)^3, \quad (14)$$

where $M_{\text{bin}}$ is the total mass of the binary, $e$ and $\omega$ the eccentricity and argument of periastron of the planet, $\phi$ is the true anomaly at the point in the orbit that we are measuring and $v$ the average velocity in the plane of the sky with which the planet moved between both transit mid-points. Since we know the planet is at conjunction, we can use $\sin(\phi + \omega) \approx 1.4$.

We alter our version of KIMA and add an extra feature, to include the ‘1–2 punch’ information as part of 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:

$$\log(\mathcal{L}) = \log(\mathcal{L}_{\text{RV}}) - \frac{1}{2} \log(2\pi\sigma_{P_{12}}^2) - \frac{(P_{12} - P_{\text{pl}})^2}{2\sigma_{P_{12}}^2}, \quad (15)$$

where $\log(\mathcal{L}_{\text{RV}})$ is the log-likelihood from the radial velocity data, $P_{\text{pl}}$ is the period for the planet proposed as part of the sampling process, $P_{12}$ is the orbital period calculated using equation (14) using all other parameters proposed in the sample (e.g. $M_{\text{bin}}$, $e$, etc.). Finally, $\sigma_{P_{12}}^2$ is the variance of $P_{12}$, which is derived from the uncertainty in the transit times propagated through equation (14).

We run the analysis on TIC 172 900 988 again, with the extra input of two transits with mid-times at 2 458 883.390879 ± 0.006188 and 2 458 888.309427 ± 0.003904.<fnref n="5" /> 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,

<footnote n="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.</footnote>

<footnote n="4">⁴We note a difference from Kostov et al. (2020) where the value of this is quoted as −1, Kostov et al. (2016) use a value of 1, this may be down to the choice of reference frame.</footnote>

<footnote n="5">⁵The transit mid-times were determined using the Eclipsing Light Curve (ELC) code. The segment of the data containing a single event was isolated, and a model was fitted to the transit (or eclipse) profile. Then we execute the DE-MCMC code, and the median along with $1\sigma$ uncertainties of the posterior sample was considered as the best-fitting transit time.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** 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 radial velocity data using SD-GP, and blue from the CCF-GP.</figcaption></figure>

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 Fig. 11 for both the CCF-GP method and the SD-GP method. We note that the crescent-shaped correlations involving the eccentricity are expected. We also note that while we report the median and $1\sigma$ in Table 1, some of the parameters have non-Gaussian distributions [in particular the eccentricity ($e_{pl}$) and argument of periastron ($\omega_{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 $e_{pl} > 0.1$ and an argument of periastron $4.35 \le \omega_{pl} \le 5.31$ rad.

Fig. A6 shows the orbital configuration of TIC 172900988. This is the counterpart to Fig. 9 generated from posterior samples obtained

from the KIMA analysis which included the 1–2 punch and using the SD-GP radial velocities.

## 5 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 \text{ 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⁻¹ 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⁻¹; 2009, 2010). We speculate that this increased scatter is most likely of stellar origin. Both stellar components are fairly high-mass stars. We note that TODCOR and TODMOR, both recognized 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 *K*<sub>pl</sub> ~ 42 m s⁻¹. The fact our approaches both manage to overcome some of that scatter emphasizes the limitations of existing techniques when dealing with systems characterized 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 optimizations to these methods may lead to even more accurate and precise radial velocity measurements, particularly with the spectral decomposition.

Firstly, we recognize 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 recognize 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 of 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, v sin i and metallicity, an important parameter to relate planet presence to planet formation (e.g. Santos, Israelian & Mayor 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 radial-velocity follow-up of circumbinary planet candidates identified with the transit method.

## ACKNOWLEDGEMENTS

We thank Tsevi Mazeh for suggesting we compare our methods to TODCOR and TODMOR, during the 2022 EAS conference in Valencia where preliminary results were presented.

This paper is based on observations collected the Observatoire de Haute Provence (OHP). Here, we would also like to acknowledge the staff, particularly the night assistants, at the OHP for their dedication and hard work, especially during the COVID pandemic.

This research received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no 803193/BEBOP) and from the Leverhulme Trust (research project grant no RPG-2018-418). The French group acknowledges financial support from the French Programme National de Planétologie (PNP, INSU). The acquisition of data was made possible by a series of allocations through the French PNP. The computations described in this paper were performed using the University of Birmingham’s BlueBEAR high performance computing service (available at http://www.birmingham.ac.uk/bear).

## DATA AVAILABILITY

The radial velocity data are included in the appendix of this paper. The reduced spectra are available at the SOPHIE archive: http://atlas.obs-hp.fr/sophie/ (Table A2).

## REFERENCES

* Adibekyan V. Z. et al., 2013, A&A, 560, A51

* Aigrain S., Parviainen H., Pope B. J. S., 2016, MNRAS, 459, 2408

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

* Baycroft T. A., Triaud A. H. M. J., Faria J., Correia A. C. M., Standing M. R., 2023, MNRAS, 521, 1871

* Borucki W. J., Summers A. L., 1984, Icarus, 58, 121

* Bouchy F. et al., 2009, A&A, 505, 853

* Brewer B. J., Foreman-Mackey D., 2018, J. Stat. Softw., 86, 1

* Byrd R. H., Lu P., Nocedal J., Zhu C., 1995, SIAM J. Sci. Comput., 16, 1190

* Czekala I., Mandel K. S., Andrews S. M., Dittmann J. A., Ghosh S. K., Montet B. T., Newton E. R., 2017, ApJ, 840, 49

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

* Faria J. P., Santos N. C., Figueira P., Brewer B. J., 2018, J. Open Source Softw., 3, 487

* Foreman-Mackey D., 2016, J. Open Source Softw., 1, 24

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Foreman-Mackey D., Agol E., Ambikasaran S., Angus R., 2017, AJ, 154, 220

* Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5, 65

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

* Husser T. O., Wende-von Berg S., Dreizler S., Homeier D., Reiners A., Barman T., Hauschildt P. H., 2013, A&A, 553, A6

* Kass R. E., Raftery A. E., 1995, J. Am. Statist. Assoc., 90, 773

* Konacki M., Muterspaugh M. W., Kulkarni S. R., Helminiak K. G., 2009, ApJ, 704, 513

* Konacki M., Muterspaugh M. W., Kulkarni S. R., Helminiak K. G., 2010, ApJ, 719, 1293

* Kostov V. B. et al., 2016, ApJ, 827, 86
Kostov V. B. et al., 2020, AJ, 159, 253

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

Kovaleva D., Malkov O., Yungelson L., Chulkov D., 2016, Baltic Astron., 25, 419

Martin D. V., 2018, Populations of Planets in Multiple Star Systems. Springer International Publishing, Cham

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

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

Mayor M., Queloz D., 1995, Nature, 378, 355

Mazeh T., Zucker S., 1994, Ap&SS, 212, 349

Pepe F., Mayor M., Galland F., Naef D., Queloz D., Santos N. C., Udry S., Burnet M., 2002, A&A, 388, 632

Perruchot S. et al., 2008, in McLean I. S., Casali M. M., eds, Proc. SPIE Conf. Ser. Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II. SPIE, Bellingham, p. 70140J

Rajpaul V. M., Aigrain S., Buchhave L. A., 2020, MNRAS, 492, 3960

Rasmussen C. E., Williams C. K. I., 2006, Gaussian Processes for Machine Learning. MIT Press, Cambridge, Massachusetts

Santos N. C., Israelian G., Mayor M., 2004, A&A, 415, 1153

Schneider J., 1994, Planet. Space Sci., 42, 539

Socia Q. J. et al., 2020, AJ, 159, 94

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

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

Triaud A. H. M. J. et al., 2017, A&A, 608, A129

Triaud A. H. M. J. et al., 2022, MNRAS, 511, 3561

Zechmeister M., Kürster M., 2009, A&A, 496, 577

Zucker S., Mazeh T., 1994, ApJ, 420, 806

Zucker S., Mazeh T., Santos N. C., Udry S., Mayor M., 2004, A&A, 426, 695

APPENDIX: SOME EXTRA MATERIAL
<figure><img src="image_13.png" /><figcaption>**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 realization drawn from posterior predictive distribution are shown 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 radial velocities from such chunks are weighted.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**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 2000 random MCMC steps.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure A3.** The residual periodgram after removing the signal $P_b$ 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 $P_b$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure A5.** The density of posterior samples from kima run on TIC172900988 with $N_p$ fixed to 1. The orange and blue lines indicate the detection limits for method 1 and method 2.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure A4.** Left panel method 1 (SD-GP) and right panel method 2 (CCF-GP): Grey the scatter plot showing the planets consistent with the posterior distribution on the apsidal precession rate, the GR precession rate is accounted for both this and the third-body dynamical precession rate are calculated using the equations in Baycroft et al. (2023). Blue the locations of the solution presented in this work. Red the 6 suggested solutions from Kostov et al. (2021). Dashed line the HW-stability limit. Blue line the detection limit with planet b removed.</figcaption></figure>
<table><caption>**Table A1.** 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 not excluded here so these parameters are from the full posterior (parameters from the posteriors excluding unstable samples can be found in Table 1) The parameters are determined along with their corresponding $1\sigma$ uncertainties.</caption><thead><tr><th>Parameters</th><th>Method 1<br>SD-GP</th><th>Method 2<br>CCF-GP</th><th>Method 1 (with 1–2 punch)<br>SD-GP</th><th>Method 2 (with 1–2 punch)<br>CCF-GP</th></tr></thead><tbody><tr><td><strong>Binary parameters</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>P<sub>B</sub> (d)</td><td>19.657878<sup>+0.000030</sup><sub>-0.000035</sub></td><td>19.657861 ± 0.000039</td><td>19.657872<sup>+0.000033</sup><sub>-0.000038</sub></td><td>19.657865<sup>+0.000037</sup><sub>-0.000040</sub></td></tr><tr><td>e<sub>B</sub></td><td>0.448234<sup>+0.000090</sup><sub>-0.000071</sub></td><td>0.44823 ± 0.00010</td><td>0.448232 ± 0.000083</td><td>0.44821 ± 0.00011</td></tr><tr><td>ω<sub>B</sub> (rad)</td><td>1.23146<sup>+0.00020</sup><sub>-0.00027</sub></td><td>1.23136<sup>+0.00034</sup><sub>-0.00038</sub></td><td>1.23142<sup>+0.00022</sup><sub>-0.00026</sub></td><td>1.23132<sup>+0.00042</sup><sub>-0.00038</sub></td></tr><tr><td>K<sub>B</sub> (km s<sup>-1</sup>)</td><td>58.5494<sup>+0.0079</sup><sub>-0.0054</sub></td><td>58.5278 ± 0.0076</td><td>58.5557<sup>+0.0067</sup><sub>-0.0075</sub></td><td>58.5320 ± 0.0085</td></tr><tr><td>q<sub>B</sub></td><td>0.97196<sup>+0.00016</sup><sub>-0.00022</sub></td><td>0.97106 ± 0.00025</td><td>0.97211<sup>±0.00022</sup></td><td>0.97118 ± 0.00026</td></tr><tr><td>ω&#775;<sub>B</sub>(arcsec yr<sup>-1</sup>)</td><td>215 ± 59</td><td>279<sup>+73</sup><sub>-69</sub></td><td>230<sup>+62</sup><sub>-68</sub></td><td>279<sup>+87</sup><sub>-77</sub></td></tr><tr><td>T<sub>0,B</sub> (BJD)</td><td>2 459 566.00248<sup>+0.00044</sup><sub>-0.00060</sub></td><td>2 459 566.00209<sup>+0.00080</sup><sub>-0.00087</sub></td><td>2 459 566.00267<sup>+0.00053</sup><sub>-0.00067</sub></td><td>2 459 566.00236<sup>+0.00089</sup><sub>-0.00086</sub></td></tr><tr><td><strong>Planet parameters</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>P<sub>pl</sub> (d)</td><td>151.3 ± 1.8</td><td>149.2<sup>+2.2</sup><sub>-2.5</sub></td><td>151.7<sup>+1.4</sup><sub>-1.7</sub></td><td>148.6<sup>+2.1</sup><sub>-2.4</sub></td></tr><tr><td>e<sub>pl</sub></td><td>&lt;0.21</td><td>&lt;0.19</td><td>0.153<sup>+0.095</sup><sub>-0.034</sub></td><td>0.141<sup>+0.059</sup><sub>-0.019</sub></td></tr><tr><td>ω<sub>pl</sub> (rad)</td><td>6.0 ± 1.5</td><td>4.9<sup>+1.8</sup><sub>-2.0</sub></td><td>5.37<sup>+0.34</sup><sub>-0.69</sub></td><td>4.57<sup>+0.75</sup><sub>-0.61</sub></td></tr><tr><td>K<sub>pl</sub> (m s<sup>-1</sup>)</td><td>40.0 ± 5.2</td><td>43.8 ± 6.9</td><td>39.1 ± 4.4</td><td>44.6 ± 5.9</td></tr><tr><td>T<sub>0,pl</sub> (BJD)</td><td>2 459 442 ± 36</td><td>2 459 425<sup>+43</sup><sub>-45</sub> ± 36</td><td>2 459 428.6<sup>+8.2</sup><sub>-16.4</sub></td><td>2 459 417<sup>+17</sup><sub>-14</sub></td></tr><tr><td><strong>Derived<br>parameters</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>M<sub>A</sub> (M<sub>☉</sub>)</td><td>1.23681<sup>+0.00037</sup><sub>-0.00039</sub></td><td>1.23766<sup>+0.00051</sup><sub>-0.00047</sub></td><td>1.23670<sup>±0.00036</sup></td><td>1.23763<sup>+0.00056</sup><sub>-0.00052</sub></td></tr><tr><td>M<sub>B</sub> (M<sub>☉</sub>)</td><td>1.20207<sup>+0.00032</sup><sub>-0.00027</sub></td><td>1.20184<sup>+0.00035</sup><sub>-0.00031</sub></td><td>1.20221<sup>+0.00032</sup><sub>-0.00029</sub></td><td>1.20196<sup>+0.00046</sup><sub>-0.00036</sub></td></tr><tr><td>m<sub>pl</sub> (M<sub>Jup</sub>)</td><td>1.88 ± 0.25</td><td>2.05 ± 0.33</td><td>1.82 ± 0.21</td><td>2.08 ± 0.28</td></tr><tr><td>a<sub>bin</sub> (au)</td><td>0.191879<sup>+0.000016</sup><sub>-0.000018</sub></td><td>0.191894<sup>+0.000021</sup><sub>-0.000019</sub></td><td>0.191879 ± 0.000017</td><td>0.191896<sup>+0.000026</sup><sub>-0.000020</sub></td></tr><tr><td>a<sub>pl</sub> (au)</td><td>0.7478<sup>+0.0056</sup><sub>-0.0060</sub></td><td>0.7408<sup>+0.0073</sup><sub>-0.0082</sub></td><td>0.7488<sup>+0.0046</sup><sub>-0.0057</sub></td><td>0.7388<sup>+0.0070</sup><sub>-0.0080</sub></td></tr><tr><td><strong>Other parameters</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>V<sub>sys,A</sub> (km s<sup>-1</sup>)</td><td>25.9867 ± 0.0047</td><td>26.0254<sup>+0.0056</sup><sub>-0.0058</sub></td><td>26.0146 ± 0.0045</td><td>26.0531 ± 0.0057</td></tr><tr><td>V<sub>sys,B</sub> (km s<sup>-1</sup>)</td><td>26.0459 ± 0.0060</td><td>26.0970<sup>+0.0088</sup><sub>-0.0090</sub></td><td>26.0639<sup>+0.0058</sup><sub>-0.0060</sub></td><td>26.1144<sup>+0.0087</sup><sub>-0.0089</sub></td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Figure A6.** Orbital configuration of TIC 172900988 showing the orbits of the binary and the planet. The green orbits are a random sample of 1000 posteriors from kima fitting the radial velocities from SD-GP. The blue orbits are the 6 suggested solutions from Kostov et al. (2021), the dashed grey line is the stability limit as calculated by Holman & Wiegert (1999). The radial velocity data are fit along with the 1–2 punch transit times.</figcaption></figure>
<table><caption>**Table A2.** Radial velocities for TIC 172900988 from method 1 (SD-GP).</caption><thead><tr><th>Time</th><th>Primary velocity (km s<sup>-1</sup>)</th><th>Secondary velocity (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>2459147.672</td><td>52.79199671 ± 0.006297356</td><td>-1.341327759 ± 0.006306834</td></tr><tr><td>2459171.636</td><td>91.97097921 ± 0.005992599</td><td>-41.7090222 ± 0.006286077</td></tr><tr><td>2459175.691</td><td>-23.52844488 ± 0.00644886</td><td>77.15182984 ± 0.006474733</td></tr><tr><td>2459181.596</td><td>7.347251678 ± 0.006376945</td><td>45.27023063 ± 0.006680279</td></tr><tr><td>2459190.691</td><td>92.70303819 ± 0.006187038</td><td>-42.5376363 ± 0.006332743</td></tr><tr><td>2459239.365</td><td>-1.129883138 ± 0.006176608</td><td>53.91700289 ± 0.006657609</td></tr><tr><td>2459241.492</td><td>14.08887308 ± 0.005837241</td><td>38.29504511 ± 0.00627873</td></tr><tr><td>2459244.379</td><td>37.5613718 ± 0.006226358</td><td>14.18090136 ± 0.005747394</td></tr><tr><td>2459266.368</td><td>60.76693256 ± 0.005976897</td><td>-9.64611818 ± 0.006047371</td></tr><tr><td>2459267.47</td><td>73.56029354 ± 0.006017367</td><td>-22.86220394 ± 0.006961391</td></tr><tr><td>2459273.489</td><td>-21.54275388 ± 0.005940911</td><td>75.02796022 ± 0.005957771</td></tr><tr><td>2459275.417</td><td>-20.68023758 ± 0.006389657</td><td>74.17527381 ± 0.006500415</td></tr><tr><td>2459277.469</td><td>-9.111063764 ± 0.006309079</td><td>62.28799085 ± 0.006159118</td></tr><tr><td>2459279.412</td><td>3.98339319 ± 0.006324967</td><td>48.82902088 ± 0.006526375</td></tr><tr><td>2459298.342</td><td>-1.020137615 ± 0.006330255</td><td>54.01401655 ± 0.006227837</td></tr><tr><td>2459303.354</td><td>37.55683484 ± 0.0074648</td><td>14.27880206 ± 0.006027881</td></tr><tr><td>2459306.305</td><td>67.86532276 ± 0.006444953</td><td>-16.94025943 ± 0.006238535</td></tr><tr><td>2459349.348</td><td>73.62773529 ± 0.006983177</td><td>-22.86640098 ± 0.006296323</td></tr><tr><td>2459500.669</td><td>44.32693159 ± 0.006540338</td><td>7.125655004 ± 0.006049167</td></tr><tr><td>2459503.676</td><td>77.25865033 ± 0.006945358</td><td>-26.80055497 ± 0.006602366</td></tr><tr><td>2459505.66</td><td>92.93966618 ± 0.005800538</td><td>-42.7710604 ± 0.006673939</td></tr><tr><td>2459520.677</td><td>47.73082394 ± 0.006196444</td><td>3.716136893 ± 0.006353854</td></tr><tr><td>2459538.668</td><td>32.37888203 ± 0.006456517</td><td>19.38534042 ± 0.00619455</td></tr><tr><td>2459540.703</td><td>51.37790001 ± 0.005821974</td><td>-0.173462113 ± 0.006350227</td></tr><tr><td>2459547.591</td><td>-3.063614773 ± 0.006373649</td><td>55.9661014 ± 0.006367542</td></tr><tr><td>2459549.642</td><td>-23.7319162 ± 0.006405175</td><td>77.21294787 ± 0.006348954</td></tr><tr><td>2459551.67</td><td>-15.38293878 ± 0.00663949</td><td>68.53497782 ± 0.005787423</td></tr><tr><td>2459554.647</td><td>4.143067667 ± 0.006584258</td><td>48.46498983 ± 0.006984549</td></tr><tr><td>2459561.567</td><td>64.5010298 ± 0.006332341</td><td>-13.53132504 ± 0.006124583</td></tr><tr><td>2459564.581</td><td>93.16248364 ± 0.006361619</td><td>-42.99898808 ± 0.005963061</td></tr><tr><td>2459593.454</td><td>0.683072938 ± 0.006046486</td><td>52.26118897 ± 0.006787273</td></tr><tr><td>2459594.52</td><td>8.129553129 ± 0.007271555</td><td>44.54464934 ± 0.006443498</td></tr></tbody></table>

<table><caption>**Table A2** – continued</caption><thead><tr><th>Time</th><th>Primary velocity (km s<sup>-1</sup>)</th><th>Secondary velocity (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>2459595.551</td><td>15.80628424 ± 0.006564392</td><td>36.64180076 ± 0.00590492</td></tr><tr><td>2459601.498</td><td>71.67462121 ± 0.006331687</td><td>-20.77559905 ± 0.006184636</td></tr><tr><td>2459603.604</td><td>93.11739363 ± 0.005906311</td><td>-42.88178153 ± 0.005967396</td></tr><tr><td>2459606.465</td><td>0.153346735 ± 0.00681785</td><td>52.83101539 ± 0.006158358</td></tr><tr><td>2459620.521</td><td>64.26088258 ± 0.00619908</td><td>-13.29834682 ± 0.006548925</td></tr><tr><td>2459622.461</td><td>87.01039359 ± 0.00627504</td><td>-36.66314116 ± 0.006265194</td></tr><tr><td>2459633.431</td><td>5.236038169 ± 0.006660511</td><td>47.48907668 ± 0.006490337</td></tr><tr><td>2459639.522</td><td>56.92154995 ± 0.006093159</td><td>-5.735897826 ± 0.006651901</td></tr><tr><td>2459641.539</td><td>80.41509092 ± 0.006584094</td><td>-29.84720911 ± 0.006356079</td></tr><tr><td>2459644.369</td><td>67.1501422 ± 0.006394344</td><td>-16.30315308 ± 0.006183854</td></tr><tr><td>2459648.487</td><td>-22.40190367 ± 0.006323709</td><td>75.85843987 ± 0.005795921</td></tr><tr><td>2459689.378</td><td>-14.79066973 ± 0.006230134</td><td>68.00513033 ± 0.00637691</td></tr><tr><td>2459700.34</td><td>78.32983042 ± 0.006419625</td><td>-27.75129827 ± 0.006155616</td></tr><tr><td>2459703.361</td><td>66.33099669 ± 0.006216234</td><td>-15.43700494 ± 0.005620311</td></tr><tr><td>2459866.605</td><td>-12.90221977 ± 0.005845783</td><td>66.08328391 ± 0.006040523</td></tr><tr><td>2459867.575</td><td>-6.810615935 ± 0.005947117</td><td>59.73342189 ± 0.005773165</td></tr><tr><td>2459883.659</td><td>-23.8159665 ± 0.00648577</td><td>77.31777816 ± 0.006012319</td></tr><tr><td>2459883.678</td><td>-23.774671 ± 0.006567066</td><td>77.31375696 ± 0.006051184</td></tr><tr><td>2459885.631</td><td>-16.60526503 ± 0.006354589</td><td>69.98110091 ± 0.005834438</td></tr><tr><td>2459887.693</td><td>-3.678472876 ± 0.006754954</td><td>56.62591044 ± 0.00647946</td></tr><tr><td>2459903.653</td><td>-23.40576746 ± 0.006284945</td><td>77.00199673 ± 0.005950399</td></tr><tr><td>2459921.561</td><td>-14.17031488 ± 0.006465102</td><td>67.36293503 ± 0.006204055</td></tr><tr><td>2459955.49</td><td>73.37103301 ± 0.006235832</td><td>-22.72623035 ± 0.006626476</td></tr><tr><td>2459956.447</td><td>84.73111962 ± 0.006153014</td><td>-34.38022166 ± 0.006084119</td></tr><tr><td>2459956.553</td><td>85.91070219 ± 0.006025253</td><td>-35.6867799 ± 0.006509933</td></tr><tr><td>2459957.585</td><td>93.21325513 ± 0.00645166</td><td>-43.09054178 ± 0.006232159</td></tr><tr><td>2460040.371</td><td>-22.81230331 ± 0.006385354</td><td>76.3388599 ± 0.006613679</td></tr><tr><td>2460041.4</td><td>-23.15597144 ± 0.00728653</td><td>76.62366086 ± 0.006686246</td></tr><tr><td>2460070.347</td><td>40.62749601 ± 0.006123952</td><td>11.0440643 ± 0.006351236</td></tr></tbody></table>

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