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

Thomas A. Baycroft,1 Amaury H.M.J. Triaud,1 João Faria,2,3 Alexandre C.M.
Correia,4,5 and Matthew R. Standing6
1School of Physics and Astronomy, University of Birmingham, Edgbaston,
Birmingham B15 2TT, UK
2Instituto de Astrofíısica e Ciências do Espaço, Universidade do Porto, CAUP,
Rua das Estrelas, 4150-762 Porto, Portugal
3 Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do
Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal
4CFisUC, Departamento de Física, Universidade de Coimbra, 3004-516 Coimbra,
Portugal
5IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Université, 77 av. Denfert-
Rochereau, 75014 Paris, France
6School of Physical Sciences, The Open University, Milton Keynes, MK7 6AA, UK
E-mail<EMAIL_ADDRESS>
(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 characterise 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 multi-planetary systems. In addition, we apply the algorithm to
several real systems, producing a new measurement of aspidal 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.

###### keywords:

binaries: general – planets and satellites: dynamical evolution and stability
– techniques: radial velocities – software: data analysis

††pubyear: 2022††pagerange: Improving circumbinary planet detections by
fitting their binary’s apsidal precession–B

##  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
multi-planetary 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
planets111Circumbinary 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. 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 et al., 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 BEBOP survey.
BEBOP (Binaries Escorted By Orbiting Planets) 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 $\omega$ of the binary orbit. An issue this raises is that
variation in $e$ and $\omega$ 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 $\omega$ 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, 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_{\rm p}\sin i_{\rm 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 Sect. 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 et al., 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
later, 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,$ (1)

with $\omega$ now being time dependent222We neglect terms that are order
$\mathcal{O}(t-t_{0})^{2}$:

$\omega(t)=\omega_{0}+\dot{\omega}(t-t_{0}).$ (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 Eqs. (3) and (4). We denote these as
observational period, $P_{\rm obs}$ which is the time taken peak-to-peak in
radial velocity, and the anomalistic period333This 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), $P_{\rm ano}$, which
is the time between consecutive pericentre passages.

$\frac{2\pi}{P_{\rm obs}}\approx\dot{\omega}+\dot{M},$ (3) $\frac{2\pi}{P_{\rm
ano}}\approx\dot{M}.$ (4)

$P_{\rm 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_{\rm obs}$ for this. When including $\dot{\omega}$ into radial
velocity fits we want to use $P_{\rm ano}$ as the period parameter to avoid
the expected correlation between $P_{\rm 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_{\rm obs}}=\dot{\omega}+\frac{2\pi}{P_{\rm ano}},$ (5)

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

$\displaystyle P_{\rm ano}=\frac{P_{\rm obs}}{\left(1-\frac{\dot{\omega}P_{\rm
obs}}{2\pi}\right)}\approx P_{\rm obs}\left(1+\frac{\dot{\omega}P_{\rm
obs}}{2\pi}\right),$ (6) $\displaystyle P_{\rm obs}=\frac{P_{\rm
ano}}{\left(1+\frac{\dot{\omega}P_{\rm ano}}{2\pi}\right)}\approx P_{\rm
ano}\left(1-\frac{\dot{\omega}P_{\rm ano}}{2\pi}\right).$ (7)

Our model fits for $P_{\rm obs}$ as a parameter (i.e. the period prior is for
$P_{\rm obs}$ as is the output posterior distribution), but the model
internally converts this to $P_{\rm 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 above, 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)444Sybilski et al. (2013)
also have an equation for the gravitational redshift, but it contains errors,
hence we use the equation from Zucker & Alexander (2007):

$\displaystyle\Delta V_{\rm LT}$
$\displaystyle=\frac{K_{1}^{2}}{c}\sin^{2}(f+\omega)(1+e\cos{f}),$ (8)
$\displaystyle\Delta V_{\rm TD}$
$\displaystyle=\frac{K_{1}^{2}}{c\sin^{2}{i}}\left(1+e\cos{f}-\frac{1-e^{2}}{2}\right),$
(9) $\displaystyle\Delta V_{\rm GR}$
$\displaystyle=\frac{K_{1}(K_{1}+K_{2})}{c\sin^{2}{i}}(1+e\cos{f}),$ (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_{\rm
tide}=1184\,\frac{M_{2}R_{1}^{4}}{M_{1}(M_{1}+M_{2})}P^{-3}\sin^{2}i\sin[2(f-\phi_{0})]\>\rm{m\,s^{-1}},$
(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_{\rm{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
calibrations555Even though one would expect the same $\gamma$ for components
observed with the same instrument, this may not be the case (Southworth,
2013). 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 dataset 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 dataset 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 prior 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 Eq. (11) assumes a circular orbit. We
therefore recommend not using the tidal correction for eccentric binaries. We
currently assume an inclination of $90^{\circ}$, 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 marginalise over.

The use of double lined binaries is also included in the options. This
requires a dataset (or multiple) with 5 columns: date; RV of primary;
uncertainty on primary RV; RV of secondary; 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 prior 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, 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_{\rm 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\%$ 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 above 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 (Eq. (23)), and due to rotational and relativistic
effects (Eq. (29)). 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 (Eq. (23)) 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 Eq. (17) 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$ $\rm
M_{\odot}$, a secondary with $M=0.37~{}\rm M_{\odot}$, $P=21.08$ days, and
$e=0.16$ and a (roughly Jupiter mass) planet with $M=0.001$ $M_{\odot}$,
$P=134.5$ days, and $e=0.01$. The first simulation, SIM1, has just these three
bodies, whereas the second simulation, SIM2, has an additional planet with
$M_{\rm pl}=0.00015~{}\rm M_{\odot}$, $P=911.2$ days, 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 dataset including
all Newtonian perturbations. Both simulated datasets are given a Gaussian
white noise.

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.

| rebound | kima | kima-binaries | units
---|---|---|---|---
$P_{\rm{B}}$ | 21.0805330(8492) | 21.0810474(52) | 21.0810395(21) | days
$M_{\rm{B}}$ | 0.37 | 0.3700886(153) | 0.3700981(57) | $\mathrm{M_{\odot}}$
$K_{\rm{B}}$ | 23415.30(45) | 23419.67(80) | 23420.11(31) | $\mathrm{ms^{-1}}$
$e_{\rm{B}}$ | 0.160106(53) | 0.160114(34) | 0.160096(15) |
$\omega_{\rm{B}}$ | 4.50018(236) | 4.50061(24) | 4.50016(20) | rad
$\phi_{0,\rm{B}}$ | 6.27400(23628) | 6.28256(24) | 6.28288(26) | rad
$\dot{\omega}_{\rm{B}}$ | 308.4(3.7) | 0 | 304.0(15.8) | $\mathrm{arcsec\,yr^{-1}}$
$P_{\rm{pl}}$ | 135.084(1.147) | 131.373(120) | 131.392(57) | days
$M_{\rm{pl}}$ | 1.0476 | 1.076(27) | 1.046(14) | $\mathrm{M_{J}}$
$K_{\rm{pl}}$ | 33.62(09) | 34.81(90) | 33.82(43) | $\mathrm{ms^{-1}}$
$e_{\rm{pl}}$ | 0.0161(79) | <0.050 | <0.033 |
RMS |  | 6.61 | 3.15 | $\mathrm{ms^{-1}}$
Jitter |  | 6.27 | 2.04 | $\mathrm{ms^{-1}}$
$\chi_{\nu}^{2}$ |  | 12.35 | 3.15 |

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). 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% detection limits. The black dashed lines show
where bodies of various masses would sit on this plot.

####  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\,\rm arcsec\,yr^{-1}$
for the binary. Using Eq. (23) we get a theoretical value of
$\dot{\omega}=301.1^{+0.8}_{-1.5}\,\rm arcsec\,yr^{-1}$ A fit using the
binaries model results in a posterior 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 the $16^{\rm th}$ to the $84^{\rm
th}$ percentiles 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 datapoint 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
anomalistic period as defined in Sect. 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 characterised (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
posterior 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 Figure 1 where the Root-Mean-
Square improves from ${\rm RMS}=6.61~{}\rm m\,s^{-1}$ to $3.15~{}\rm
m\,s^{-1}$. When apsidal precession is not accounted for, if we move further
from the reference 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 Figure 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 $1,000$ and $2,000~{}\rm days$. Whilst this
sounds impressive, this simulation only had a very small amount of extra white
noise added to maximise the effect of apsidal precision in order to reveal its
importance. In other systems we may expect more marginal improvements (see
Sect. 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 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}=\dfrac{N_{i+1}}{N_{i}}$ (12)

Should $N_{i}=0$, the Bayes Factor in Eq. (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 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.

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. | kima | kima-binaries
---|---|---
$BF_{1,0}$ | 538 | 0
$BF_{2,1}$ | 12.5 | 6086
$BF_{3,2}$ | 0.9 | 0.8

###  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, 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\,\rm days$ (Yenawine et al.,
2022) which corresponds to $\dot{\omega}=21\,650\,\rm 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\,\rm arcsec\,yr^{-1}$. Equivalently this is
$\dot{\omega}=0.56\pm 0.009$ degrees per cycle.

We run the analysis with apsidal precession fit for, as described in Sect. 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}_{\nu}$ 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 Figure 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 Figure 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
$\sim 10^{4}~{}\rm days$ of the inner tertiary.

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

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% 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.

##  4 Application of kima-binaries to Kepler-16

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.)

The announcement of Kepler-16b marked the first unambiguous detection of a
circumbinary planet (Doyle et al., 2011), made thanks 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 with our new model,
and successfully detect an apsidal precession rate of
$\dot{\omega}_{1}=283^{+87}_{-85}~{}\rm arcsec\,yr^{-1}$, which is $3.3\sigma$
from 0. Using Eqs. (23, 29) we obtain a value of
$\dot{\omega}_{1}=92.4^{+14.3}_{-13.8}~{}\rm 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 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.

We explore the difference between $P_{\rm obs}$ and $P_{\rm 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_{\rm obs}$ is in statistical
agreement with the previous publication, however $P_{\rm ano}$ is $3.3\sigma$
above this. $P_{\rm 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_{\rm
ano}$ and $P_{\rm 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}\,\mathrm{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}\,\mathrm{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 Figure 5. In this case, as we
only get a marginal detection of apsidal precession there is no real
improvement in the detection limits.

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% confidence detection limits with and without including precession. Type of Period and algorithm used | Value
---|---
$P$ (yorbit)* | 41.077779(54)
$P$ (kima)* | 41.077772(51)
$P_{\rm obs}$ (kima-binaries) | 41.077716(55)
$P_{\rm ano}$ (kima-binaries) | 41.078737(305)
Table 3: Various binary periods for Kepler-16AB. The first two values (*) are
taken from Triaud et al. (2022) using their two different algorithms, the
other two are obtained using kima-binaries 1$\sigma$ uncertainties are shown
in brackets as the last 2 significant figures (3 in the last case for easy
comparison with the others)

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

##  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 (Eq. (23)) 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 n∘
803193/BEBOP). A.C. 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).

## 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 et al. (2022) Angelo I., Naoz S., Petigura E., MacDougall M., Stephan A. P., Isaacson H., Howard A. W., 2022, The Astronomical Journal, 163, 227
  * Arras et al. (2012) Arras P., Burkart J., Quataert E., Weinberg N. N., 2012, Monthly Notices of the Royal Astronomical Society, 422, 1761
  * Bonfils et al. (2013) Bonfils X., et al., 2013, Astronomy & Astrophysics, 549, A109
  * Borkovits et al. (2021) Borkovits T., et al., 2021, Monthly Notices of the Royal Astronomical Society, 510, 1352
  * Borucki et al. (2010) Borucki W. J., et al., 2010, Science, 327, 977
  * Brewer et al. (2011) Brewer B. J., Pártay L. B., Csányi G., 2011, Statistics and Computing, 21
  * Carter et al. (2011) Carter J. A., et al., 2011, Science, 331, 562
  * Claret & Giménez (2010) Claret A., Giménez A., 2010, Astronomy and Astrophysics, 519, A57
  * Correia et al. (2013) Correia A. C. M., Boué G., Laskar J., Morais M. H. M., 2013, Astronomy & Astrophysics, 553, A39
  * Correia et al. (2016) Correia A. C. M., Boué G., Laskar J., 2016, Celestial Mechanics and Dynamical Astronomy, 126, 189
  * Dawson & Johnson (2018) Dawson R. I., Johnson J. A., 2018, Annual Review of Astronomy and Astrophysics, 56, 175
  * Doolin & Blundell (2011) Doolin S., Blundell K. M., 2011, Monthly Notices of the Royal Astronomical Society, 418, 2656
  * Doyle et al. (2011) Doyle L. R., et al., 2011, Science, 333, 1602
  * Dvorak et al. (1989) Dvorak R., Froeschle C., Froeschle C., 1989, Astronomy and Astrophysics, 226, 335
  * Farago & Laskar (2010) Farago F., Laskar J., 2010, Monthly Notices of the Royal Astronomical Society, 401, 1189
  * Faria et al. (2018) Faria J. P., Santos N. C., Figueira P., Brewer B. J., 2018, The Journal of Open Source Software, 3, 487
  * Gillon et al. (2016) Gillon M., et al., 2016, Nature, 533, 221
  * Gillon et al. (2017) Gillon M., et al., 2017, Nature, 542, 456
  * Holman & Wiegert (1999) Holman M. J., Wiegert P. A., 1999, The Astronomical Journal, 117, 621
  * Konacki et al. (2009) Konacki M., Muterspaugh M. W., Kulkarni S. R., Hełminiak K. G., 2009, The Astrophysical Journal, 704, 513
  * Konacki et al. (2010) Konacki M., Muterspaugh M. W., Kulkarni S. R., Hełminiak K. G., 2010, The Astrophysical Journal, 719, 1293
  * Kostov et al. (2016) Kostov V. B., et al., 2016, The Astrophysical Journal, 827, 86
  * Kostov et al. (2020) Kostov V. B., et al., 2020, The Astronomical Journal, 159, 253
  * Kostov et al. (2021) Kostov V. B., et al., 2021, The Astronomical Journal, 162, 234
  * Kovaleva et al. (2016) Kovaleva D., Malkov O., Yungelson L., Chulkov D., 2016, Baltic Astronomy, 25, 419
  * Martin & Triaud (2014) Martin D. V., Triaud A. H. M. J., 2014, Astronomy & Astrophysics, 570, A91
  * Martin & Triaud (2015) Martin D. V., Triaud A. H. M. J., 2015, Monthly Notices of the Royal Astronomical Society, 449, 781
  * Martin et al. (2019) Martin D. V., et al., 2019, Astronomy and Astrophysics, 624, A68
  * Mayor et al. (2011) Mayor M., et al., 2011, Technical report, The HARPS search for southern extra-solar planets XXXIV. Occurrence, mass distribution and orbital properties of super-Earths and Neptune-mass planets
  * Morais & Correia (2012) Morais M. H. M., Correia A. C. M., 2012, Monthly Notices of the Royal Astronomical Society, 419, 3447
  * Murray & Correia (2010) Murray C. D., Correia A. C. M., 2010, Keplerian Orbits and Dynamics of Exoplanets
  * Murray & Dermott (1999) Murray C. D., Dermott S. F., 1999, Solar system dynamics. Cambridge University Press
  * Orosz et al. (2012) Orosz J. A., et al., 2012, Science, 337, 1511
  * Rein & Liu (2012) Rein H., Liu S.-F., 2012, Astronomy and Astrophysics, 537, A128
  * Rein & Spiegel (2015) Rein H., Spiegel D. S., 2015, Monthly Notices of the Royal Astronomical Society, 446, 1424
  * Rosenthal et al. (2021) Rosenthal L. J., et al., 2021, The Astrophysical Journal Supplement Series, 255, 8
  * Rosu et al. (2020) Rosu S., Rauw G., Conroy K. E., Gosset E., Manfroid J., Royer P., 2020, Astronomy & Astrophysics, 635, A145
  * Schneider (1994) Schneider J., 1994, Planetary and Space Science, 42, 539
  * Southworth (2013) Southworth J., 2013, Astronomy and Astrophysics, 557, A119
  * Standing et al. (2022) Standing M. R., et al., 2022, Monthly Notices of the Royal Astronomical Society, 511, 3571
  * Sybilski et al. (2013) Sybilski P., Konacki M., Kozłowski S. K., Hełminiak K. G., 2013, Monthly Notices of the Royal Astronomical Society, 431, 2024
  * Triaud et al. (2017) Triaud A. H. M. J., et al., 2017, Astronomy & Astrophysics, 608, A129
  * Triaud et al. (2022) Triaud A. H. M. J., et al., 2022, Monthly Notices of the Royal Astronomical Society, 511, 3561
  * Winn & Fabrycky (2015) Winn J. N., Fabrycky D. C., 2015, Annual Review of Astronomy and Astrophysics, 53, 409
  * Yenawine et al. (2022) Yenawine M. E., et al., 2022, The Astrophysical Journal, 924, 66
  * Zucker & Alexander (2007) Zucker S., Alexander T., 2007, The Astrophysical Journal, 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}\left[2-12e_{1}^{2}-6(1-e_{1}^{2})\cos^{2}i+30e_{1}^{2}\cos^{2}\alpha\right],$
(13)

where

$C_{2}=\frac{\mathcal{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}},$
(14)

and

$\cos i=\sin I_{1}\sin I_{2}\cos(\Omega_{1}-\Omega_{2})+\cos I_{1}\cos I_{2},$
(15) $\displaystyle\cos\alpha=$ $\displaystyle\sin I_{1}\sin\omega_{1}\cos
I_{2}-\sin I_{2}\cos\omega_{1}\sin(\Omega_{1}-\Omega_{2})$ (16)
$\displaystyle-\sin I_{2}\cos I_{1}\sin\omega_{1}\cos(\Omega_{1}-\Omega_{2}).$

In these equations, $a$, $e$, $I$, $\omega$ and $\Omega$ refer to the semi-
major 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. $\mathcal{G}$ refers to the
gravitational constant, while $m_{0}$, $m_{1}$ and $m_{2}$ are the masses of
the star A, B and planet, respectively. $\cos i$ and $\cos\alpha$ 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 percientre 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})}{e_{1}G_{1}}\frac{\partial\mathcal{H}}{\partial
e_{1}}+\frac{\cot I_{1}}{G_{1}}\frac{\partial\mathcal{H}}{\partial I_{1}},$
(17)

where $G_{1}$ is the norm of the orbital angular momentum

$G_{1}=\frac{m_{0}m_{1}}{m_{0}+m_{1}}\sqrt{\mathcal{G}(m_{0}+m_{1})a_{1}(1-e_{1}^{2})},$
(18) $\frac{\partial\mathcal{H}}{\partial
e_{1}}=C_{2}\left[-24e_{1}+12e_{1}\cos^{2}i+60e_{1}\cos^{2}\alpha\right],$
(19)

$\frac{\partial\mathcal{H}}{\partial I_{1}}=C_{2}\left[-12(1-e1^{2})\cos
i\frac{\partial\cos i}{\partial
I_{1}}+60e_{1}^{2}\cos\alpha\frac{\partial\cos\alpha}{\partial I_{1}}\right],$
(20)

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},$ (21)
$\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}).$ (22)

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 (17) as

$\displaystyle\frac{d\omega_{1}}{dt}$
$\displaystyle\approx\frac{12C_{2}}{G_{1}}(1-e_{1}^{2})\left[2-\cos^{2}i-5\cos^{2}\alpha\right]$
(23)
$\displaystyle\approx\frac{12C_{2}}{G_{1}}(1-e_{1}^{2})\left[1+(1-5\cos^{2}\omega_{1})\sin^{2}i\right].$

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 modeled using a Hamiltonian formalism
as666 We assume that the spin axes of both stars are normal to the orbit.
(e.g. Correia et al., 2013, 2016)

$\begin{split}\mathcal{H}^{\prime}=&-\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}}\left(1+3e_{1}^{2}+\tfrac{3}{8}e_{1}^{4}\right)\
,\end{split}$ (24)

where

$C_{g}=\frac{3\mathcal{G}^{2}m_{0}m_{1}(m_{0}+m_{1})}{a_{1}^{2}c^{2}}$ (25)

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

$C_{r,i}=\frac{\mathcal{G}m_{0}m_{1}J_{2,i}R_{i}^{2}}{2a_{1}^{3}}$ (26)

accounts for the rotational flattening, and

$C_{t,i}=k_{2,i}\frac{\mathcal{G}m_{1-i}^{2}R_{i}^{5}}{2a_{1}^{6}}$ (27)

for the tidal contribution.

$J_{2,i}=k_{2,i}\frac{\Omega_{i}^{2}R_{i}^{3}}{3\mathcal{G}m_{i}}\ ,$ (28)

$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 (17), the correction in the apsidal precession
is given by

$\begin{split}\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_{t,0}+C_{t,1})}{G_{1}(1-e_{1}^{2})^{9/2}}\left(1+\tfrac{3}{2}e^{2}+\tfrac{1}{8}e^{4}\right)\
.\end{split}$ (29)

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 Sect. (4)

##  Appendix B Tables of parameters

Here we show the parameters for KOI-126 in Table 4 and for Kepler-16 in Table
5 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 4: 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^{2}_{\nu}$) are taken for the best fitting model. The mass $M_{\rm 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. | kima-binaries | units
---|---|---
assumed parameters |  |
$M_{\rm A}$ | 1.2713(47) | $\mathrm{M_{\odot}}$
fitted parameters |  |
$P_{\rm obs}$ | 33.77943(33) | days
$P_{\rm ano}$ | 33.83207(69) | days
$K$ | 21395(25) | $\mathrm{m\,s^{-1}}$
$e$ | 0.3113(13) |
$\omega$ | 1.1794(50) | rad
$\phi_{0}$ | 1.0210(40) | rad
$\dot{\omega}$ | 21800(370) | $\mathrm{arcsec\,yr^{-1}}$
$\gamma$ | $-27852^{+193}_{-76}$ | $\mathrm{m\,s^{-1}}$
derived parameters |  |
$M_{\rm B+C}^{*}$ | 0.4424(11) | $\mathrm{M_{\odot}}$
$T_{0}$ | 51047.4547(27) | BJD - 2,400,000
fit indicators |  |
$\rm RMS_{Tull}$ | 207 | $\mathrm{m\,s^{-1}}$
$\rm RMS_{Tres}$ | 84.2 | $\mathrm{m\,s^{-1}}$
$\rm Jitter_{Tull}$ | 0.36 | $\mathrm{m\,s^{-1}}$
$\rm Jitter_{Tres}$ | 49.3 | $\mathrm{m\,s^{-1}}$
$\chi_{\nu}^{2}$ | 3.07 |
Table 5: 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^{2}_{\nu}$) 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. | kima-binaries | units
---|---|---
assumed parameters |  |
$M_{\rm A}$ | 0.654(20) | $\mathrm{M_{\odot}}$
fitted parameters |  |
$P_{\rm B,obs}$ | 41.077716(55) | days
$P_{\rm B,ano}$ | 41.07874(31) | days
$K_{\rm B}$ | 13678.9(1.4) | $\mathrm{m\,s^{-1}}$
$e_{\rm B}$ | 0.159925(88) |
$\omega_{\rm B}$ | 4.60203(79) | rad
$\phi_{0,{\rm B}}$ | 1.63340(76) | rad
$\dot{\omega}_{\rm B}$ | 284(86) | $\mathrm{arcsec\,yr^{-1}}$
$P_{\rm b}$ | 225.8(1.7) | days
$K_{\rm b}$ | 11.7(1.6) | $\mathrm{m\,s^{-1}}$
$e_{\rm b}$ | <0.29 |
$\omega_{\rm b}$ | 3.90(92) | rad
$\phi_{0,{\rm b}}$ | 2.32(87) | rad
$\gamma$ | $-33811.5^{+3.1}_{-0.2}$ | $\mathrm{m\,s^{-1}}$
derived parameters |  |
$M_{B}$ | 0.1965(32) | $\mathrm{M_{\odot}}$
$T_{0,{\rm B}}$ | 58498.4796(52) | BJD - 2,400,000
$M_{b}$ | 0.308(42) | $\mathrm{M_{Jup}}$
$T_{0,{\rm b}}$ | 58388(33) | BJD - 2,400,000
fit indicators |  |
$\rm RMS$ | 11.01 | $\mathrm{m\,s^{-1}}$
$\rm Jitter$ | 0.91 | $\mathrm{m\,s^{-1}}$
$\chi_{\nu}^{2}$ | 0.92 |