# Improving circumbinary planet detections by fitting their binary's apsidal precession

Thomas A. Baycroft, Amaury H. M. J. Triaud, João Faria, Alexandre C. M. Correia, Matthew R. Standing

► To cite this version:

Thomas A. Baycroft, Amaury H. M. J. Triaud, João Faria, Alexandre C. M. Correia, Matthew R. Standing. Improving circumbinary planet detections by fitting their binary's apsidal precession. Monthly Notices of the Royal Astronomical Society, 2023, 521, pp.1871-1879. 10.1093/mn-ras/stad607. insu-04849173

HAL Id: insu-04849173

https://insu.hal.science/insu-04849173v1

Submitted on 19 Dec 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 Attribution 4.0 International License
# Improving circumbinary planet detections by fitting their binary's apsidal precession

Thomas A. Baycroft<sup>1,*</sup> Amaury H. M. J. Triaud<sup>1</sup>,<sup>1</sup> João Faria<sup>1,2,3</sup> Alexandre C. M. Correia<sup>1,4,5</sup> and Matthew R. Standing<sup>1,6</sup>

Accepted 2023 February 22. Received 2023 January 26; in original form 2022 November 21

## ABSTRACT

Apsidal precession in stellar binaries is the main non-Keplerian dynamical effect impacting the radial-velocities of a binary star system. Its presence can notably hide the presence of orbiting circumbinary planets because many fitting algorithms assume perfectly Keplerian motion. To first order, apsidal precession ($\dot{\omega}$) can be accounted for by adding a linear term to the usual Keplerian model. We include apsidal precession in the *kima* package, an orbital fitter designed to detect and characterize planets from radial velocity data. In this paper, we detail this and other additions to *kima* that improve fitting for stellar binaries and circumbinary planets including corrections from general relativity. We then demonstrate that fitting for $\dot{\omega}$ can improve the detection sensitivity to circumbinary exoplanets by up to an order of magnitude in some circumstances, particularly in the case of multiplanetary systems. In addition, we apply the algorithm to several real systems, producing a new measurement of apsidal precession in KOI-126 (a tight triple system), and a detection of $\dot{\omega}$ in the Kepler-16 circumbinary system. Although apsidal precession is detected for Kepler-16, it does not have a large effect on the detection limit or the planetary parameters. We also derive an expression for the precession an outer planet would induce on the inner binary and compare the value this predicts with the one we detect.

**Key words:** binaries: general – planets and satellites: dynamical evolution and stability – techniques: radial velocities – methods: data analysis.

## 1 INTRODUCTION

Exoplanets exhibit a range of configurations much vaster than is present within the Solar system. Nearly three decades of discoveries have revealed that most known exoplanets are not analogous to any Solar system planets (e.g. Winn & Fabrycky 2015). This applies to individual planets being different such as hot Jupiters (e.g. Dawson & Johnson 2018) or planets with extreme eccentricities (e.g. Angelo et al. 2022), but this can also apply to entire planetary systems having more exotic configurations, such as TRAPPIST-1 a multiplanetary resonant chain orbiting a late M-dwarf (Gillon et al. 2016, 2017). One such type of exotic planetary systems are the circumbinary exoplanets that orbit about both stars of a tight stellar binary (Schneider 1994; Doyle et al. 2011).

To date there have been only 15 fully confirmed circumbinary planets.<fnref n="1" /> All but one have transited at least one of the two stars,

and were first detected from space with *Kepler* (e.g. Doyle et al. 2011; Orosz et al. 2012; Kostov et al. 2016) and with *TESS* (Kostov et al. 2020, 2021). The pace of detections is slow, since the two most common exoplanet detection techniques (transit and radial velocity) have so far both been hamstrung, each with their own issues. Circumbinary planets will generally have longer periods than planets around single stars because they need to orbit outside of an instability region produced by the binary stars' motion (Dvorak, Froeschle & Froeschle 1989; Holman & Wiegert 1999; Doolin & Blundell 2011). Because of this extra distance, circumbinary planets are geometrically less likely to produce transits than planets orbiting single stars. However, for similarly distant planets, nodal precession makes circumbinary planets more likely to create transits (Martin & Triaud 2015), even if transits do not happen at every planetary orbit (e.g. Schneider 1994; Martin & Triaud 2014).

For the radial velocity method, interference between the spectra of both components of the binary star makes it harder to obtain precise radial velocity measurements (e.g. Konacki et al. 2009). The latter issue can be circumvented by observing single-lined binaries (Konacki et al. 2010; Martin et al. 2019). This is the observing strategy employed by the Binaries Escorted By Orbiting Planets

<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>Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="3"><sup>3</sup>Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal</footnote>

<footnote n="4"><sup>4</sup>CFisUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra, Portugal</footnote>

<footnote n="5"><sup>5</sup>IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Université, 77 avenue Denfert-Rochereau, 75014 Paris, France</footnote>

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

<footnote>*E-mail: txb187@bham.ac.uk</footnote>

<footnote n="1">¹Circumbinary planets orbiting stellar remnants are claimed using eclipse timings, however, there are doubts about their existence and as such we do not consider them as fully confirmed.</footnote>
(BEBOP) survey. BEBOP has been collecting radial velocities on eclipsing single-lined binaries for over four years and has demonstrated that it can detect circumbinary planets, notably by having independently detected Kepler-16b (Triaud et al. 2022). There is also the first circumbinary planet discovered in radial velocities BEBOP-1c (Standing et al. submitted).

One significant advantage of the radial velocity method over the transit method is that radial velocities probe the full orbit, instead of just the inferior conjunction. This leads to good precision on the eccentricity *e* (sometimes down to $10^{-4}$; Triaud et al. 2017) and the argument of periastron *ω* of the binary orbit. An issue this raises is that variation in *e* and *ω* will cause problems when fitting a static Keplerian orbit, a point raised in Konacki et al. (2010); Sybilski et al. (2013). One such variation is apsidal precession of the binary: an evolution of *ω* with time, denoted by $\dot{\omega}$. This can be caused by relativistic effects, tidal effects, or – most excitingly to exoplanet hunters – planetary perturbations (Correia et al. 2013). Because the BEBOP survey has collected data over several years, the scatter caused by this precession is slowly starting to exceed the rms scatter of the residuals on some systems (Standing et al. 2022). Accounting for this effect would improve the accuracy of the fits for the binary orbit, and in turn improve our ability to detect planets, and the precision on their physical and orbital parameters. As explored in Standing et al. (2022), in most cases the radial-velocity signal of a single planet would be detected well before a non-zero $\dot{\omega}$ is significantly detected. In this paper we show that, in some multiplanetary configurations, low-amplitude planetary signals can be hidden by the precession induced by another, heavier planet.

In addition, measuring the apsidal precession rate adds new information to our knowledge of the system. Usually, it is not possible to measure orbital inclinations from radial velocities alone. However, the binary's precession rate due to an external perturber is dependent on the mutual inclination between the binary's orbital plane and the perturber's (e.g. Correia et al. 2013; Correia, Boué & Laskar 2016). In this work, we derive an equation to calculate the apsidal precession that a third body induces on an inner binary pair which can be used to calculate the mutual inclination, this can be found in Appendix (A). In the case of BEBOP, where all binaries are known to eclipse, an upper bound on the mutual inclination directly translates into an upper bound on the orbital inclination of the planet, meaning those radial velocity data can be used to obtain not just a minimum mass $m_p \sin i_p$, but also a maximum mass. Finally, for close binaries, measuring the precession rate also provides information on the stars' internal structure (Claret & Giménez 2010).

In this paper, we first we describe a binaries-specific radial velocity model applied in the kima package in Section 2. This new model (which we will occasionally refer to as kima-binaries when comparing it with the old model) moves beyond fitting pure Keplerian orbits (as done in Faria et al. 2018), by including an apsidal precession parameter to the fitted model. The section details those changes and describes the inclusion of other tidal and relativistic effects that are known to affect orbital solutions. In Section 3, the new model is used on both simulated and observed data. The ability to accurately recover the apsidal precession rate is demonstrated, and we show how fitting for apsidal precession can improve a survey's sensitivity to circumbinary planets by producing Bayesian detection limits. Finally, in Section 4 we present a detection of the precession rate for Kepler-16, and conclude in Section 5.

## 2 A BINARY UPDATE TO KIMA

In this section we present an update to kima, developing a binary-specific radial velocity model. This model accounts for various

factors that are generally ignored when looking at radial velocities for a single star, but recommended when seeking to detect circumbinary planet signals (Konacki et al. 2009; Sybilski et al. 2013). The new model includes tidal and relativistic effects as well as, most notably, apsidal precession of the binary's orbit. The new model is also given the capability to fit double-lined binary data.

kima is an orbital fitting algorithm which makes use of diffusive nested sampling (DNest; Brewer, Pártay & Csányi 2011) to sample the posterior distribution for the model parameters. It allows for the number of Keplerian signals being fit to vary freely which is advantageous for Bayesian model comparison. There is a so called 'known-object' mode where separate priors can be defined for certain already known signals while allowing to search for further signals freely; this model is ideal to apply to circumbinary systems. As will be discussed further, this method of sampling allows for an efficient method of calculating detection limits.

### 2.1 Adding precession to kima

#### 2.1.1 A linear approximation

As a first order approximation, we add a linear precession parameter, $\dot{\omega}$ to kima. This parameter is free during a fit and its posterior is estimated. We take the usual equation for the radial velocity of a Keplerian orbit (e.g Murray & Correia 2010):

$$V = K(\cos(f + \omega) + e \cos(\omega)) + \gamma, \quad (1)$$

with $\omega$ now being time dependent:<fnref n="2" />

$$\omega(t) = \omega_0 + \dot{\omega}(t - t_0). \quad (2)$$

$K$ is the semi-amplitude of the radial velocity signal, $f$ the true anomaly, $e$ and $\omega$ the eccentricity and argument of pericentre for the orbit, $t_0$ is some reference time, and $\gamma$ is the mean velocity of the system (which can be affected by the zero-point calibration of an instrument, but does not impact other parameters). In our case, we use the mean of the times of observation for $t_0$.

#### 2.1.2 Period correction

The period of an orbit as a single value is, in any realistic scenario, not completely well defined. Various angles associated with an orbit will vary in time, such as the argument of pericentre $\omega$ or the mean anomaly $M$. Different combinations of these variations could all be called periods. We consider two of these defined as in the equations (3) and (4). We denote these as observational period, $P_{\text{obs}}$ which is the time taken peak-to-peak in radial velocity, and the anomalistic period,<fnref n="3" /> $P_{\text{ano}}$, which is the time between consecutive pericentre passages.

$$\frac{2\pi}{P_{\text{obs}}} \approx \dot{\omega} + \dot{M}, \quad (3)$$

$$\frac{2\pi}{P_{\text{ano}}} \approx \dot{M}. \quad (4)$$

$P_{\text{obs}}$ is the period that is usually referred to by observational astronomers and can be precisely measured from time between transits or eclipses. In this work, we set priors on the binary period based on eclipses so want to use $P_{\text{obs}}$ for this. When including $\dot{\omega}$ into

<footnote n="2">²We neglect terms that are order $O(t - t_0)^2$</footnote>

<footnote n="3">³This nomenclature is often used to refer to the time between two consecutive pericentre passages in precessing systems (e.g. Rosu et al. 2020; Borkovits et al. 2021).</footnote>
radial velocity fits we want to use $P_{\text{ano}}$ as the period parameter to avoid the expected correlation between $P_{\text{obs}}$ and $\dot{\omega}$. Hence we need to be able to convert from one to the other. To do this we combine the two equations to get

$$ \frac{2\pi}{P_{\text{obs}}} = \dot{\omega} + \frac{2\pi}{P_{\text{ano}}}, \qquad (5) $$

and hence, neglecting terms of order $O(\dot{\omega}P)^2$,

$$ P_{\text{ano}} = \frac{P_{\text{obs}}}{\left(1 - \frac{\dot{\omega}P_{\text{obs}}}{2\pi}\right)} \approx P_{\text{obs}} \left(1 + \frac{\dot{\omega}P_{\text{obs}}}{2\pi}\right), \qquad (6) $$

$$ P_{\text{obs}} = \frac{P_{\text{ano}}}{\left(1 + \frac{\dot{\omega}P_{\text{ano}}}{2\pi}\right)} \approx P_{\text{ano}} \left(1 - \frac{\dot{\omega}P_{\text{ano}}}{2\pi}\right). \qquad (7) $$

Our model fits for $P_{\text{obs}}$ as a parameter (i.e. the period prior is for $P_{\text{obs}}$ as is the output posterior distribution), but the model internally converts this to $P_{\text{ano}}$.

## 2.2 Other additions to the binaries model

Here, we describe the other additions made to the binary model on top of the apsidal precession described earlier, namely, we add relativistic and tidal corrections, and give the ability to fit the radial velocities for a double-lined binary.

### 2.2.1 Relativistic and tidal corrections

We include relativistic corrections for binary orbits, the main ones being light-travel time and transverse doppler (LT, TD; Sybilski et al. 2013) and gravitational redshift (GR; Zucker & Alexander 2007):<fnref n="4" />

$$ \Delta V_{\text{LT}} = \frac{K_1^2}{c} \sin^2(f + \omega)(1 + e \cos f), \qquad (8) $$

$$ \Delta V_{\text{TD}} = \frac{K_1^2}{c \sin^2 i} \left(1 + e \cos f - \frac{1 - e^2}{2}\right), \qquad (9) $$

$$ \Delta V_{\text{GR}} = \frac{K_1(K_1 + K_2)}{c \sin^2 i} (1 + e \cos f), \qquad (10) $$

where $e, f, \omega$, and $i$ are respectively the eccentricity, true anomaly, argument of pericentre, and inclination of the binary orbit relative to the plane of the sky; $K_1$ and $K_2$ are the semi-amplitudes of the primary and secondary, respectively, and $c$ is the speed of light.

The tidal effect is calculated as in Arras et al. (2012), assuming a circular orbit. The equation for the tidally induced radial velocity signal is as follows:

$$ v_{\text{tide}} = 1184 \frac{M_2 R_1^4}{M_1(M_1 + M_2)} P^{-3} \sin^2 i \sin[2(f - \phi_0)] \text{ m s}^{-1}, \qquad (11) $$

where $M_1, M_2$, and $R_1$ are the mass and radius of the primary and secondary (in solar units), $P$ the orbital period (in days) (we use $P_{\text{ano}}$), $f$ the true anomaly, and $\phi_0 = \pi/2 - \omega$ is the observer's reference position.

These equations are incorporated as an optional feature into the model such that when a binary model is fit, these contributions to the radial velocities can be naturally accounted for. We do not include these effects for the general planet search objects as their effects will be much smaller (by orders of magnitude since these corrections scale with $M^2$) and so does not warrant the increase in computation time.

### 2.2.2 Adding in a double-lined binary model

The vast majority of spectroscopic binaries are double-lined (e.g. Kovaleva et al. 2016), and although the detection of circumbinary planets in double-lined system is problematic (as in Konacki et al. 2009, 2010), new methods to disentangle both spectral components accurately enough to detect circumbinary planets are being developed (e.g. Lalitha et al. in prep). To prepare for the time when circumbinary planets can be searched for, and be detected in double-lined binaries, we add a feature to kima to model such a configuration. As an input, the software requires files containing radial-velocities for each component of the binary. The sets of data are fit simultaneously, each with an independent $\gamma$ parameter to account for differing zero-point calibrations.<fnref n="5" /> In addition, each set has its own jitter term, added in quadrature to the RV uncertainties to account for any additional sources of white noise. Only one extra common parameter is fit, the mass ratio $q$.

Any given solution consists of a binary orbit, some number of planetary orbits, and polynomial trends up to cubic order. The binary orbit is fit to each data set with the secondary having the semi-amplitude $K$ scaled by $q$, and its argument of periastron reversed $\omega_2 = \omega_1 - \pi$. The planetary orbits are then fit in the same way to each data set just as kima usually does.

## 2.3 Using the model

The additions in the new binary model can be used in various combinations. The tidal correction and relativistic correction can each be turned on or off, they will then apply to any known objects included. A priori on $\dot{\omega}$ will need to be given for each known object as well as for general signals.

One important thing to note is that equation (11) assumes a circular orbit. We therefore recommend not using the tidal correction for eccentric binaries. We currently assume an inclination of 90°, and therefore we only consider eclipsing systems in this paper. A future update may include the inclination as a free parameter to either attempt to constrain or at the very least marginalize over.

The use of double lined binaries is also included in the options. This requires a data set (or multiple) with five columns: date; RV of primary; uncertainty on primary RV; RV of secondary; and uncertainty on secondary RV. The primary will therefore automatically be the signal placed in the second column. The mass ratio $q$ can be larger than 1 (at which point the 'secondary' is actually the more massive star), so for example in an almost equal mass case a priori can straddle $q = 1$.

## 3 PERFORMANCES OF KIMA-BINAIREs

We now show tests and applications of the binaries model using data from simulations as well as from real systems. We first show that the model is able to recover consistent values of apsidal precession, and demonstrate the improvement in the fit that ensues. We illustrate this improvement by computing detection limits.

The standard way to perform detection limits is to inject a fine grid of simulated Keplerian signals (often assuming $e = 0$) into the data where any planetary signal has been removed, and to measure which signals are recovered by the algorithm (e.g. Konacki et al. 2009, 2010; Mayor et al. 2011; Bonfils et al. 2013; Rosenthal et al. 2021). Here,

<footnote n="4">⁴Sybilski et al. (2013) also have an equation for the gravitational redshift, but it contains errors, hence we use the equation from Zucker & Alexander (2007)</footnote>

<footnote n="5">⁵Even though one would expect the same $\gamma$ for components observed with the same instrument, this may not be the case (Southworth 2013).</footnote>
instead, we use the posterior distribution of the undetected Keplerian
signals to measure the amplitudes which can still be present in the
data, as described in Standing et al. (2022).

Briefly, if the analysis indicates there are no planets in a system (circumstellar or circumbinary), we then apply a strict prior on the number of planets by fixing $N_p = 1$. One Keplerian signal is assumed to be present and the algorithm is thus forced to return all solutions that are compatible with the data, but not formally detected. We then analyse the posterior samples and compute a limit of *K* as a function of *P* that envelops the lower 99 per cent of the samples. Practically speaking, the limit is produced by creating log-linear bins along the *P* axis. It is best to ensure there are at least at least 1000 samples in each bin. Should a system have a formally detected planet (Bayes factor exceeding 150), that planet is subtracted from the data (the maximum likelihood parameters are used for this), and the detection limit is then computed as given earlier using the residual radial velocities.

The advantage of using such a method over traditional methods
is the ability to sample over all orbital parameters as finely as the
algorithm allows (indeed also, all $\gamma$ variables and *jitters*, as well
as $e, \omega, \dot{\omega}$, and $\phi_0$ which are sometimes avoided by the traditional
methods). While a traditional insertion/recovery asks the question
*can these exact signals be recovered?* our method instead asks *what
is compatible with the data?* A planet below the detection threshold is
consistent with the data, and thus, is not formally detected, whereas a
planet above the line is inconsistent with the data and would therefore
have been detected had it been there.

**3.1 Analytic equation for precession**

In Appendix A we derive an analytic equation for the expected precession rate due to an outer perturber (equation (A11)), and due to rotational and relativistic effects (equation A17). This is done in a similar way to in Correia et al. (2013) but where that was done in the invariant plane, we do the calculation in the sky plane, which is directly applicable to observations results.

The precession due to a perturber (equation A11) is under the
assumption of both an eclipsing and transiting system, such as
Kepler-16. If applying this to a system that does not conform to
these assumptions, then equation (A5) should be used.

**3.2 Testing kima-binaries with simulated data**

We begin by testing our ability to recover the apsidal precession using
simulated data, showing that we can recover a good measurement
of the apsidal precession rate $\dot{\omega}$, and that including it can greatly
improve the fit.

We perform two simulations, both with a primary star of mass $M = 1$ M$_\odot$, a secondary with $M = 0.37$M$_\odot$, $P = 21.08$ d, and $e = 0.16$ and a (roughly Jupiter mass) planet with $M = 0.001$ M$_\odot$, $P = 134.5$ d, and $e = 0.01$. The first simulation, SIM1, has just these 3 bodies, whereas the second simulation, SIM2, has an additional planet with $M_{pl} = 0.00015$M$_\odot$, $P = 911.2$ d, and $e = 0$, corresponding to about 3 times the mass of Neptune. We chose these parameters to emulate a typical circumbinary planet: the binary is similar to Kepler-16 in mass-ratio and eccentricity, with a shorter period to increase the amount of precession that will have happened across the time that we 'observe'. Planet 1 was placed between the 6:1 and 7:1 mean-motion resonances with the binary and Planet 2 at a similar period ratio again. Three different masses of Planet 2 were tried and we report here the one that had the right mass to be missed without using precession but detected when including it. The decimal places for the periods

were chosen randomly to try and avoid integer numbers of days and
potential accidental resonances.

Simulations are made using the rebound package (Rein & Liu 2012), the integrations used the IAS15 integrator (Rein & Spiegel 2015). Radial velocity simulations are taken as the velocity along the line of sight within the simulation. The simulation uses the same observational cadence as for Kepler-16 (Triaud et al. 2022), thus producing a simulated data set including all Newtonian perturbations. Both simulated data sets are given a Gaussian white noise.

**3.2.1 Improving scatter and derived parameters**

First, we consider SIM1. From the values of $\omega$ at each point in the rebound simulation, we obtain $\dot{\omega} = 308.4 \pm 3.7$ arcsec yr$^{-1}$ for the binary. Using equation (A11) we get a theoretical value of $\dot{\omega} = 301.1_{-1.5}^{+0.8}$ arcsec yr$^{-1}$. A fit using the binaries model results in a *posteriori* estimate of $\dot{\omega} = 304 \pm 16$ arcsec/yr, which is in agreement with both the simulated and theoretical values. This run is done with the apsidal precession of the binary fit for, but without the relativistic or tidal corrections. The uncertainty in the rebound value comes from 'sampling' $\omega$ at various times, calculating $\dot{\omega}$ from these and then taking its mean and variance. The uncertainty on the theoretical value is propagated in a Monte Carlo way from the posterior uncertainty on the binary and planetary parameters. The *kima-binaries* value's uncertainty is defined from 16 to 84 per cent of the posterior distribution.

Table 1 lists the parameters of the binary and planet taken from rebound and fit with precession (*kima-binaries*) and without precession (*kima*). The $1\sigma$ uncertainty in each measurement is shown in brackets as the last few significant figures, to make comparison easier these are all scaled so that each value on a row is shown to the same number of decimal places. The parameters for rebound are read out as the osculating parameters at the times of each data point and then the mean and standard deviation of the values are calculated.

We note that in many cases the fitted values are inconsistent
with the rebound values, and give a word of warning for using
the Keplerian parameters from an *N*-body fitter such as this. When
taking the Keplerian orbital parameters of a body from a rebound
simulation at a given time, these are taken from the osculating
Keplerian orbit which may not be representative of the average orbit.
Consider the planet's orbital period, rebound effectively gives us
an anomalous period as defined in Section 2.1.2. Because of the
perturbed motion this is not the time it will take to actually complete
one orbit and we get an observed period a few days shorter. This
effect cannot be fit as an apsidal precession of the planet as the orbit
is not detectably eccentric.

So a Keplerian (or quasi-Keplerian) fit does not reproduce the
osculating Keplerian parameters from a *N*-body simulation, but it
does (to a reasonable accuracy) reproduce the mass. We see in Table
1 that the mass of the binary is accurate to 3 decimal places (which
is more than the precision we usually get on the mass of the primary
star anyway) and the mass of the planet is accurately characterized
(more so when apsidal precession is take into account).

We can also compare the precision of the two fits, in the sense of
how tight a *posteriori* distribution we get for each parameter. In most
cases we can see an improvement in precision by about a factor of
two.

The reduction in residual scatter can be seen in Fig. 1 where the
rms improves from rms = 6.61m s⁻¹ to 3.15m s⁻¹. When apsidal
precession is not accounted for, if we move further from the reference
<table><caption>**Table 1.** For SIM1 we show the parameters from rebound (Rein & Liu 2012) (note these are keplerian parameters taken from a Newtonian simulation) alongside the fitted parameters both with precession (kima-binaries) and without (kima). For the planet we state the upper bound on the eccentricity and omit the angle parameters as these are not resolved (close to circular orbit). The $1\sigma$ uncertainties are shown as the last few significant digits, all of which are on the same scale as the smallest uncertainty to allow for easy comparison. Goodness-of-fit parameters are also shown to compare the two fits.</caption><thead><tr><th></th><th>rebound</th><th>kima</th><th>kima-binaries</th><th>Units</th></tr></thead><tbody><tr><td>P<sub>B</sub></td><td>21.0805330(8492)</td><td>21.0810474(52)</td><td>21.0810395(21)</td><td>days</td></tr><tr><td>M<sub>B</sub></td><td>0.37</td><td>0.3700886(153)</td><td>0.3700981(57)</td><td>M<sub>☉</sub></td></tr><tr><td>K<sub>B</sub></td><td>23415.30(45)</td><td>23419.67(80)</td><td>23420.11(31)</td><td>ms<sup>-1</sup></td></tr><tr><td>e<sub>B</sub></td><td>0.160106(53)</td><td>0.160114(34)</td><td>0.160096(15)</td><td></td></tr><tr><td>ω<sub>B</sub></td><td>4.50018(236)</td><td>4.50061(24)</td><td>4.50016(20)</td><td>rad</td></tr><tr><td>φ<sub>0,B</sub></td><td>6.27400(23628)</td><td>6.28256(24)</td><td>6.28288(26)</td><td>rad</td></tr><tr><td>ω̇<sub>B</sub></td><td>308.4(3.7)</td><td>0</td><td>304.0(15.8)</td><td>arcsec yr<sup>-1</sup></td></tr><tr><td>P<sub>pl</sub></td><td>135.084(1.147)</td><td>131.373(120)</td><td>131.392(57)</td><td>days</td></tr><tr><td>M<sub>pl</sub></td><td>1.0476</td><td>1.076(27)</td><td>1.046(14)</td><td>M<sub>J</sub></td></tr><tr><td>K<sub>pl</sub></td><td>33.62(09)</td><td>34.81(90)</td><td>33.82(43)</td><td>ms<sup>-1</sup></td></tr><tr><td>e<sub>pl</sub></td><td>0.0161(79)</td><td>&lt;0.050</td><td>&lt;0.033</td><td></td></tr><tr><td>rms</td><td></td><td>6.61</td><td>3.15</td><td>ms<sup>-1</sup></td></tr><tr><td>Jitter</td><td></td><td>6.27</td><td>2.04</td><td>ms<sup>-1</sup></td></tr><tr><td>χ<sup>2</sup><sub>v</sub></td><td></td><td>12.35</td><td>3.15</td><td></td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** A comparison of radial velocity residuals, after removing the binary's orbital solution, for SIM1 where apsidal precession is included in the kima-binaries model (in red) and not included in kima (in blue).</figcaption></figure>

time $T_0$ (near the centre of the figure), the fit worsens, giving a characteristic bow-tie shape, but if precession is accounted for, the spread in the residuals is reduced.

The detection limits both with and without precession can be seen in Fig. 2. The improvement in detection limit is slightly larger at high periods where the radial-velocity signature of apsidal precession can be confused for a long-term trend. This improvement means the data would allow the detection of another planet signal within this system, almost an order of magnitude lower in mass at orbital periods between 1000 and 2, 000d. Whilst this sounds impressive, this simulation only had a very small amount of extra white noise added to maximize the effect of apsidal precision in order to reveal its importance. In other systems we may expect more marginal improvements (see Section 4).

### 3.2.2 Detecting a hidden planet

Here, we consider SIM2 and test how many planets are formally detected. To register as an *n*-planet detection, the Bayes Factor

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Detection limits for additional Keplerian signals using SIM1. The hexbins represent the density of posterior samples in each run. The purple-dashed line and solid-blue line are the 99 per cent detection limits. The black-dashed lines show where bodies of various masses would sit on this plot.</figcaption></figure>

for the *n*-planet solution compared to the (*n* − 1)-planet solution needs to be greater than 150. The sampling in kima is trans-dimensional, meaning that solutions with different numbers of planets are all searched simultaneously. Therefore, the Bayes Factor $BF_{i+1,i}$ comparing the model with *i* + 1 planets to that with *i* planets, is the ratio of the number of posterior samples with *i* + 1 planets $N_{i+1}$ to the number with *i* planets $N_i$

$$BF_{i+1,i} = \frac{N_{i+1}}{N_i} \qquad (12)$$

where $N_i = 0$, the Bayes Factor in equation (12) becomes infinite. This can happen if the BF is larger than the number of effective posterior samples, and would be solved eventually had the sampling continued. In this case we therefore choose to report $BF_{i+1,i} = N_{i+1}$, effectively setting $N_i = 1$. More information can be found in Faria et al. (2018); Standing et al. (2022); Triaud et al. (2022).

Running kima on the data from SIM2 without including precession, the outer planet is not formally detected but visible within the posterior as an over-density. With $BF_{2,1} = 12.5 < 150$, it would be
<table><caption>**Table 2.** The Bayes Factors for SIM2 (containing two circumbinary planets) comparing models with increasing numbers of planets both for the standard version of *kima* and the new *kima-binaries* version, which includes apsidal precession.</caption><thead><tr><th></th><th>kima</th><th>kima-binaries</th></tr></thead><tbody><tr><td><i>BF</i><sub>1,0</sub></td><td>538</td><td>0</td></tr><tr><td><i>BF</i><sub>2,1</sub></td><td>12.5</td><td>6086</td></tr><tr><td><i>BF</i><sub>3,2</sub></td><td>0.9</td><td>0.8</td></tr></tbody></table>

classified as a candidate planet. We can attribute this non-detection to the apsidal precession since when we do fit for the precession, the planet is formally detected with $BF_{2,1} > 6086 > 150$. The Bayes Factors for each attempted fit are found in Table 2.

## 3.3 Testing kima-binaries on data from the KOI-126 system

In this section, we test the ability to recover a value for the precession rate consistent with a previous solution from literature, as well as show the improvement in sensitivity to planets that accounting for precession could bring in a highly precessing system.

KOI-126 is a compact triply-eclipsing hierarchical triple star system. It contains a roughly circular, low-mass, and tight binary which is in an eccentric orbit about a more massive tertiary star. The system was first reported in Carter et al. (2011). There are radial velocity data as well as photometry during multiple eclipses. As such, the apsidal precession rate of the tertiary orbit is well measured with a period of 21 850 d (Yenawine et al. 2022) which corresponds to $\dot{\omega} = 21 650$ arcsec yr$^{-1}$. We used 29 radial velocity data from Yenawine et al. (2022). Our fit with the new binaries model recovers a consistent value for the apsidal precession, with $\dot{\omega} = 21 800 \pm 600$ arcsec yr$^{-1}$. Equivalently this is $\dot{\omega} = 0.56 \pm 0.009^{\circ}$ per cycle.

We run the analysis with apsidal precession fit for, as described in Section 2, but do not include either the general relativity or the tidal corrections. A more complete analysis would be to use a Newtonian model, rather than Keplerian with added precession, as in Yenawine et al. (2022), however including the precession improves the $\chi^2$ from 640.9 to 1.5. This amply justifies adding an extra parameter, $\dot{\omega}$, to the fit and suggests that a full dynamical model is not necessary with the current precision of the data, hence illustrating the importance of $\dot{\omega}$, since even in this dynamically complex triple system, the linear apsidal precession removes the majority of the excess noise. The detection limits are shown in Fig. 3 for reference. Here too, including precession improves the detection limit by an order of magnitude in semi-amplitude and removes much of the long-period noise where, as with the simulated data, the precession may be being mildly confused for a long term trend.

We do not use the analytic equation derived in Appendix (A) as KOI-126 is in a different orbital configuration where the precessing orbit is the outer (rather than inner) one.

As an interesting note, the orbital periods shown in Fig. 3 would be for putative *circumtertiary planets*, of which none are known in Nature. We can nonetheless state there are no stellar or brown dwarf mass companions within ~ $10^4$ d of the inner tertiary.

We show the parameters from our fit of KOI-126 in Table B1

## 4 APPLICATION OF KIMA-BINARIES TO KEPLER-16

The announcement of Kepler-16b marked the first unambiguous detection of a circumbinary planet (Doyle et al. 2011), made thanks

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Detection limits for additional signals around KOI-126. The hexbins show the density of posterior samples with the red being those when precession is included in the fit, blue when it is not included. The dashed-purple line and solid-blue line show the 99 per cent confidence detection limits. The dashed lines show where bodies of various masses, and where the Deuterium and Hydrogen fusing limits would sit on this plot.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Kepler-16: red: histogram of the density of posterior samples for fitted value of $\dot{\omega}$ with the median and $1\sigma$ values shown in grey. Blue: histogram of the density of posterior samples for the theoretically calculated value of $\dot{\omega}$ with the median value shown in grey. (Note that $\dot{\omega}$ is not cut at zero, there are in fact posteriors below zero.)</figcaption></figure>

to the *Kepler* spacecraft (Borucki et al. 2010). This system is unique in also being the only circumbinary planet independently detected with radial velocity (Triaud et al. 2022). We re-analyse these radial-velocity data with our new model, and successfully detect an apsidal precession rate of $\dot{\omega}_1 = 283_{-85}^{+87}$ arcsec yr$^{-1}$, which is $3.3\sigma$ from 0. Using equations (A11) and (A17) we obtain a value of $\dot{\omega}_1 = 92.4_{-13.8}^{+14.3}$ arcsec yr$^{-1}$, which is $2.2\sigma$ away from the observed value. This theoretical value takes into account the planetary-induced and relativistic precessions (we do not include the rotational and tidal contributions as they would be very small in comparison and parameters like the Love numbers are not very well known). The posterior distributions for both the measured and theoretical precession rates are shown in Fig. 4.

The theoretical $\dot{\omega}$ is lower than the value that we measure; more data are required to determine how significant this discrepancy is. The difference is likely too important to be accounted for entirely by mutual inclination. An alternative (or additional) explanation could be further undetected planets contributing to the precession rate.
<table><caption>**Table 3.** Various binary periods for Kepler-16AB. The first two values denoted by asterisk (*) are taken from Triaud et al. (2022) using their two different algorithms, the other two are obtained using kima-binaries 1σ uncertainties are shown in brackets as the last 2 significant figures (3 in the last case for easy comparison with the others).</caption><thead><tr><th>Type of period and algorithm used</th><th>Value</th></tr></thead><tbody><tr><td>P (yorbit)*</td><td>41.077779(54)</td></tr><tr><td>P (kima)*</td><td>41.077772(51)</td></tr><tr><td>P<sub>obs</sub> (kima-binaries)</td><td>41.077716(55)</td></tr><tr><td>P<sub>ano</sub> (kima-binaries)</td><td>41.078737(305)</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Detection limits for Kepler-16 with and without accounting for precession. The posterior sample for each of these are plotted in the hexbins with the version including precession plotted in red and in front. Since there is no improvement the density plot for the model without precession is hard to see. The dashed-purple and solid-blue lines show the 99 per cent confidence detection limits with and without including precession.</figcaption></figure>

We explore the difference between $P_{\text{obs}}$ and $P_{\text{ano}}$. These values for Kepler-16 are presented in Table 3 alongside values published in Triaud et al. (2022). The value we get for $P_{\text{obs}}$ is in statistical agreement with the previous publication, however $P_{\text{ano}}$ is 3.3σ above this. $P_{\text{ano}}$ is the time between consecutive pericentre passages, the period that should in theory be used to compute physical parameters such as semimajor axis and planet mass, in a Keplerian context. In practice the difference is negligible due to there being a small difference between $P_{\text{ano}}$ and $P_{\text{obs}}$, as well as the uncertainty in the mass of the primary often being dominant. For Kepler-16 the difference in mass using the two periods is $\approx 2 \times 10^{-6} M_{\odot}$. It would take a case with very precise mass and a very high precession rate for this difference to be significant, even for KOI-126 (B + C), the difference is $\approx 2 \times 10^{-4} M_{\odot}$ which is about a fifth of the currently measured uncertainty.

In addition, we produce a detection limit for Kepler-16, comparing the results with and without including $\dot{\omega}$. The detection limits, plotted in the same way as the previous ones, are shown in Fig. 5. In this case, as we only get a marginal detection of apsidal precession there is no real improvement in the detection limits.

We show the parameters from our fit of Kepler-16 in Table B2.

## 5 CONCLUSIONS

We have shown that fitting for the apsidal precession of a binary's orbital parameters improves the radial velocity sensitivity to circumbinary planets. Our conclusions are in line with previous work

such as Konacki et al. (2010) and Sybilski et al. (2013), but extend theirs to a fully Bayesian framework. The improvement in the detection limits can be of up to an order of magnitude in some configurations, but in most cases, improvements are expected to be marginal as in the case of Kepler-16. Accounting for precession can also give improvements in the precision of the parameters recovered from a fit, as well as the potential to uncover planets that were hidden by precession (or instead require less data to detect the same planet).

We have derived a formula for calculating the theoretical precession induced in a binary (equation A11) and have discussed the potential use of a measurement of the apsidal precession rate as a way to constrain the mutual inclination of the planetary and binary orbital planes using this formula.

The theoretical and observed values of precession for Kepler-16 are in slight tension, this may be because of undetected planets or some other unknown mechanism.

The longer the baseline of radial velocity observations, the more important it is to account for apsidal precession. As the field progresses, and the number of data from surveys like BEBOP increases, fitting the apsidal precession of the binaries will become vital. To prepare for that time, we have presented an updated version of the kima package which is more adapted to fitting radial velocities for single and double-lined binaries. The code is made public on github.

## ACKNOWLEDGEMENTS

The authors thank the anonymous reviewer as well as Darin Ragozzine for their useful comments. This research is in part funded by the European Union's Horizon 2020 research and innovation programme (grants agreements number 803193/BEBOP). AC acknowledges support from CFisUC (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FIS-AST/7002/2020), and ENGAGE SKA (POCI-01-0145-FEDER-022217), funded by COMPETE 2020 and FCT, Portugal. MRS acknowledges support from the UK Science and Technology Facilities Council (ST/T000295/1).

## 6 DATA AVAILABILITY

The radial velocity data for Kepler-16 can be found in Triaud et al. (2022). The radial velocity data for KOI-126 can be found in Yenawine et al. (2022).

The code used for the analysis in this paper can be obtained at https://github.com/j-faria/kima

## REFERENCES

* Angelo I., Naoz S., Petigura E., MacDougall M., Stephan A. P., Isaacson H., Howard A. W., 2022, AJ, 163, 227

* Arras P., Burkart J., Quataert E., Weinberg N. N., 2012, MNRAS, 422, 1761

* Bonfils X., et al., 2013, A&A, 549, 109

* Borkovits T., et al., 2021, MNRAS, 510, 1352

* Borucki W. J., et al., 2010, Science, 327, 977

* Brewer B. J., Pártay L. B., Csányi G., 2011, Stat. Comput., 21, 649

* Carter J. A., et al., 2011, Science, 331, 562

* Claret A., Giménez A., 2010, A&A, 519, 57

* Correia A. C. M., Boué G., Laskar J., Morais M. H. M., 2013, A&A, 553, 39

* Correia A. C. M., Boué G., Laskar J., 2016, Cel. Mech. Dyn. Astron., 126, 189

* Dawson R. I., Johnson J. A., 2018, ARA&A, 56, 175

* Doolin S., Blundell K. M., 2011, MNRAS, 418, 2656

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

* Dvorak R., Froeschle C., Froeschle C., 1989, A&A, 226, 335
Farago F., Laskar J., 2010, MNRAS, 401, 1189

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

Gillon M., et al., 2016, Nature, 533, 221

Gillon M., et al., 2017, Nature, 542, 456

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

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

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

Kostov V. B., et al., 2016, ApJ, 827, L86

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., Triaud A. H. M. J., 2014, A&A, 570, 91

Martin D. V., Triaud A. H. M. J., 2015, MNRAS, 449, 781

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

Mayor M., et al., 2011, preprint (arXiv:1109.2497)

Morais M. H. M., Correia A. C. M., 2012, MNRAS, 419, 3447

Murray C. D., Correia A. C. M., 2010, Exoplanets. Keplerian Orbits and Dynamics of Exoplanets, University of Arizona Press, p. 15

Murray C. D., Dermott S. F., 1999, Solar system dynamics. Cambridge Univ. Press, Cambridge

Orosz J. A., et al., 2012, Science, 337, 1511

Rein H., Liu S.-F., 2012, A&A, 537, 128

Rein H., Spiegel D. S., 2015, MNRAS, 446, 1424

Rosenthal L. J., et al., 2021, ApJS, 255, 8

Rosu S., Rauw G., Conroy K. E., Gosset E., Manfroid J., Royer P., 2020, A&A, 635, 145

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

Southworth J., 2013, A&A, 557, 119

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

Sybilski P., Konacki M., Kozłowski S. K., Helminiak K. G., 2013, MNRAS, 431, 2024

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

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

Winn J. N., Fabrycky D. C., 2015, ARA&A, 53, 409

Yenawine M. E., et al., 2022, ApJ, 924, L66

Zucker S., Alexander T., 2007, ApJ, 654, L83

## APPENDIX A: DERIVATION OF THE PLANETARY INDUCED PRECESSION RATE

In this section, we derive the equation for the apsidal precession of an inner orbit (Binary or planet) due to an outer perturber.

The dominating term for the apsidal precession is the planet–binary gravitational interactions. We take the secular quadrupole Hamiltonian after averaging over the mean anomaly of both orbits is given by (e.g. Farago & Laskar 2010; Morais & Correia 2012)

$$ \mathcal{H} = C_2 [2 - 12e_1^2 - 6(1 - e_1^2) \cos^2 i + 30e_1^2 \cos^2 \alpha], \quad (A1) $$

where

$$ C_2 = \frac{G}{16} \frac{m_0 m_1}{m_0 + m_1} \frac{m_2}{(1 - e_2^2)^{3/2}} \frac{a_1^2}{a_2^3}, \quad (A2) $$

and

$$ \cos i = \sin I_1 \sin I_2 \cos(\Omega_1 - \Omega_2) + \cos I_1 \cos I_2, \quad (A3) $$

$$ \cos \alpha = \sin I_1 \sin \omega_1 \cos I_2 - \sin I_2 \cos \omega_1 \sin(\Omega_1 - \Omega_2) - \sin I_2 \cos I_1 \sin \omega_1 \cos(\Omega_1 - \Omega_2). \quad (A4) $$

In these equations, *a*, *e*, *I*, *ω*, and *Ω* refer to the semimajor axis, eccentricity, orbital inclination, argument of pericentre, and longitude of ascending node of an orbit, with subscripts 1 and 2 referring to

the inner and outer orbits, respectively. *G* refers to the gravitational constant, while *m*₀, *m*₁ and *m*₂ are the masses of the star A, B, and planet, respectively. Cos *i* and cos *α* are direction cosines, where the angle *i* corresponds to the true mutual inclination between the two orbital planes.

Then, we can compute the precession of the pericentre of the inner orbit using the Lagrange Planetary Equations (e.g. Murray & Dermott 1999) as

$$ \frac{d\omega_1}{dt} = -\frac{(1 - e_1^2) \partial \mathcal{H}}{e_1 G_1 \partial e_1} + \frac{\cot I_1 \partial \mathcal{H}}{G_1 \partial I_1}, \quad (A5) $$

where $G_1$ is the norm of the orbital angular momentum

$$ G_1 = \frac{m_0 m_1}{m_0 + m_1} \sqrt{G(m_0 + m_1)a_1(1 - e_1^2)}, \quad (A6) $$

$$ \frac{\partial \mathcal{H}}{\partial e_1} = C_2 [-24e_1 + 12e_1 \cos^2 i + 60e_1 \cos^2 \alpha], \quad (A7) $$

$$ \frac{\partial \mathcal{H}}{\partial I_1} = C_2 \left[ -12(1 - e_1^2) \cos i \frac{\partial \cos i}{\partial I_1} + 60e_1^2 \cos \alpha \frac{\partial \cos \alpha}{\partial I_1} \right], \quad (A8) $$

and

$$ \frac{\partial \cos i}{\partial I_1} = \cos I_1 \sin I_2 \cos(\Omega_1 - \Omega_2) - \sin I_1 \cos I_2, \quad (A9) $$

$$ \frac{\partial \cos \alpha}{\partial I_1} = \cos I_1 \sin \omega_1 \cos I_2 + \sin I_2 \sin I_1 \sin \omega_1 \cos(\Omega_1 - \Omega_2). \quad (A10) $$

In the case of an eclipsing binary and a transiting planet, such as Kepler-16 we can take $I_1 \approx I_2 \approx 90^{\circ}$, which allow us to simplify expression (A5) as

$$ \begin{aligned} \frac{d\omega_1}{dt} &\approx \frac{12C_2}{G_1}(1 - e_1^2)[2 - \cos^2 i - 5 \cos^2 \alpha] \\ &\approx \frac{12C_2}{G_1}(1 - e_1^2)[1 + (1 - 5 \cos^2 \omega_1) \sin^2 i]. \end{aligned} \quad (A11) $$

Therefore, for an eclipsing and transiting system, a constraint on the precession rate can be used to measure the mutual inclination.

For close-in binaries, additional sources of apsidal precession may become relevant, such as general relativity, rotational flattening, and tidal deformation. These effects, can also be modelled using a Hamiltonian formalism as<fnref n="6" /> (e.g. Correia et al. 2013; Correia et al. 2016)

$$ \mathcal{H}' = -\frac{C_g}{(1 - e_1^2)^{1/2}} - \frac{C_{r,0} + C_{r,1}}{(1 - e_1^2)^{3/2}} - \frac{C_{t,0} + C_{t,1}}{(1 - e_1^2)^{9/2}} (1 + 3e_1^2 + \frac{3}{8}e_1^4), \quad (A12) $$

where

$$ C_g = \frac{3G^2 m_0 m_1 (m_0 + m_1)}{a_1^2 c^2}, \quad (A13) $$

corresponds to the general relativity correction (*c* is the speed-of-light),

$$ C_{r,i} = \frac{G m_0 m_1 J_{2,i} R_i^2}{2 a_1^3}, \quad (A14) $$

accounts for the rotational flattening, and

$$ C_{t,i} = k_{2,i} \frac{G m_1^2 R_i^5}{2 a_1^6}, \quad (A15) $$

<footnote n="6">⁶We assume that the spin axes of both stars are normal to the orbit.</footnote>
for the tidal contribution.

$$J_{2,i} = k_{2,i} \frac{\Omega_i^2 R_i^3}{3Gm_i}, \qquad (A16)$$

$k_{2,i}$ is the second Love number for potential, $\Omega_i$ is the rotation rate, and $R_i$ is the radius of the star with mass $m_i$.

Then, according to expression (A5), the correction in the apsidal precession is given by

$$\frac{d\omega_1}{dt} = \frac{C_g}{G_1(1 - e_1^2)^{1/2}} + \frac{3(C_{r,0} + C_{r,1})}{G_1(1 - e_1^2)^{3/2}} + \frac{15(C_{r,0} + C_{r,1})}{G_1(1 - e_1^2)^{9/2}} (1 + \frac{3}{2}e^2 + \frac{1}{8}e^4). \qquad (A17)$$

The first term is the relativistic contribution, the middle one is the rotational contribution, and the last term is the tidal contribution. These can be taken separately as we do in Section (4)

## APPENDIX B: TABLES OF PARAMETERS

<table><caption>**Table B1.** The fitted and derived and assumed parameters for KOI-126 for the orbit of the binary B + C around star A. * since we only fit a single orbit in the radial velocity data the mass here is the combined masses of stars B and C. The quality of fit indicators (rms, Jitter, and $\chi_v^2$) are taken for the best fitting model. The mass $M_A$ is obtained from Yenawine et al. (2022). Excepting parameters with highly asymmetric distributions, the 1$\sigma$ uncertainties are shown as the last few significant digits.</caption><thead><tr><th></th><th>kima-binaries</th><th>Units</th></tr></thead><tbody><tr><td colspan="3"><strong>Assumed parameters</strong></td></tr><tr><td>$M_A$</td><td>1.2713(47)</td><td>M<sub>☉</sub></td></tr><tr><td colspan="3"><strong>Fitted parameters</strong></td></tr><tr><td>$P_{\text{obs}}$</td><td>33.77943(33)</td><td>days</td></tr><tr><td>$P_{\text{ano}}$</td><td>33.83207(69)</td><td>days</td></tr><tr><td>$K$</td><td>21 395(25)</td><td>m s<sup>-1</sup></td></tr><tr><td>$e$</td><td>0.3113(13)</td><td></td></tr><tr><td>$\omega$</td><td>1.1794(50)</td><td>rad</td></tr><tr><td>$\phi_0$</td><td>1.0210(40)</td><td>rad</td></tr><tr><td>$\dot{\omega}$</td><td>21 800(370)</td><td>arcsec yr<sup>-1</sup></td></tr><tr><td>$\gamma$</td><td>-27852<sup>+193</sup><sub>-76</sub></td><td>m s<sup>-1</sup></td></tr><tr><td colspan="3"><strong>Derived parameters</strong></td></tr><tr><td>$M_{\text{B+C}}^*$</td><td>0.4424(11)</td><td>M<sub>☉</sub></td></tr><tr><td>$T_0$</td><td>51047.4547(27)</td><td>BJD - 2 400 000</td></tr><tr><td colspan="3"><strong>Fit indicators</strong></td></tr><tr><td>rms<sub>Tull</sub></td><td>207</td><td>m s<sup>-1</sup></td></tr><tr><td>rms<sub>Tres</sub></td><td>84.2</td><td>m s<sup>-1</sup></td></tr><tr><td>Jitter<sub>Tull</sub></td><td>0.36</td><td>m s<sup>-1</sup></td></tr><tr><td>Jitter<sub>Tres</sub></td><td>49.3</td><td>m s<sup>-1</sup></td></tr><tr><td>$\chi_v^2$</td><td>3.07</td><td></td></tr></tbody></table>

Here, we show the parameters for KOI-126 in Table B1 and for Kepler-16 in Table B2 from the fits using the binaries model. The parameters for KOI-16 are of the outer orbit of the triple, modelling the short period binary as a single massive body.

<table><caption>**Table B2.** The fitted and derived and assumed parameters for Kepler-16 subscripts B refer to the binary orbit and b to the planetary orbit. The quality of fit indicators (rms, Jitter, and $\chi_v^2$) are taken for the best fitting model. Excepting parameters with highly asymmetric distributions, the $1\sigma$ uncertainties are shown as the last few significant digits.</caption><thead><tr><th></th><th>kima-binaries</th><th>Units</th></tr></thead><tbody><tr><td colspan="3"><strong>Assumed parameters</strong></td></tr><tr><td>$M_A$</td><td>0.654(20)</td><td>M<sub>☉</sub></td></tr><tr><td colspan="3"><strong>Fitted parameters</strong></td></tr><tr><td>$P_{\text{B, obs}}$</td><td>41.077716(55)</td><td>days</td></tr><tr><td>$P_{\text{B, ano}}$</td><td>41.07874(31)</td><td>days</td></tr><tr><td>$K_B$</td><td>13678.9(1.4)</td><td>m s<sup>-1</sup></td></tr><tr><td>$e_B$</td><td>0.159925(88)</td><td></td></tr><tr><td>$\omega_B$</td><td>4.60203(79)</td><td>rad</td></tr><tr><td>$\phi_{0, B}$</td><td>1.63340(76)</td><td>rad</td></tr><tr><td>$\dot{\omega}_B$</td><td>284(86)</td><td>arcsec yr<sup>-1</sup></td></tr><tr><td>$P_B$</td><td>225.8(1.7)</td><td>days</td></tr><tr><td>$K_B$</td><td>11.7(1.6)</td><td>m s<sup>-1</sup></td></tr><tr><td>$e_B$</td><td>&lt;0.29</td><td></td></tr><tr><td>$\omega_B$</td><td>3.90(92)</td><td>rad</td></tr><tr><td>$\phi_{0, B}$</td><td>2.32(87)</td><td>rad</td></tr><tr><td>$\gamma$</td><td>-33811.5<sup>+3.1</sup><sub>-0.2</sub></td><td>m s<sup>-1</sup></td></tr><tr><td colspan="3"><strong>Derived parameters</strong></td></tr><tr><td>$M_B$</td><td>0.1965(32)</td><td>M<sub>☉</sub></td></tr><tr><td>$T_0, B$</td><td>58498.4796(52)</td><td>BJD - 2 400 000</td></tr><tr><td>$M_B$</td><td>0.308(42)</td><td>M<sub>Jup</sub></td></tr><tr><td>$T_0, B$</td><td>58 388(33)</td><td>BJD - 2 400 000</td></tr><tr><td colspan="3"><strong>Fit indicators</strong></td></tr><tr><td>rms</td><td>11.01</td><td>m s<sup>-1</sup></td></tr><tr><td>Jitter</td><td>0.91</td><td>m s<sup>-1</sup></td></tr><tr><td>$\chi_v^2$</td><td>0.92</td><td></td></tr></tbody></table>

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