Improving circumbinary planet detections by ﬁtting their binary’s apsidal precession

Apsidal precession in stellar binaries is the main non-Keplerian dynamical eﬀect impacting the radial-velocities of a binary star system. Its presence can notably hide the presence of orbiting circumbinary planets because many ﬁtting algorithms assume perfectly Keplerian motion. To ﬁrst order, apsidal precession ( (cid:164) 𝜔 ) can be accounted for by adding a linear term to the usual Keplerian model. We include apsidal precession in the kima package, an orbital ﬁtter designed to detect and characterise planets from radial velocity data. In this paper, we detail this and other additions to kima that improve ﬁtting for stellar binaries and circumbinary planets including corrections from general relativity. We then demonstrate that ﬁtting for (cid:164) 𝜔 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 (cid:164) 𝜔 in the Kepler-16 circumbinary system. Although apsidal precession is detected for Kepler-16, it does not have a large eﬀect 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.


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(Gillon et al. , 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 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; ★ E-mail: txb187@bham.ac.uk 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 Orosz et al. 2012;Kostov et al. 2016) and with TESS (Kostov et al. 2020(Kostov et al. , 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 eclips-ing 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  (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  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 .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  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. 2013Correia et al. , 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  p sin  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.

A BINARY UPDATE TO KIMA
In this section we present an update to kima, developing a binaryspecific 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.

A linear approximation
As a first order approximation, we add a linear precession parameter,  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): with  now being time dependent2 : is the semi-amplitude of the radial velocity signal,  the true anomaly,  and  the eccentricity and argument of pericentre for the orbit,  0 is some reference time, and  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  0 .

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  or the mean anomaly .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,  obs which is the time taken peak-to-peak in radial velocity, and the anomalistic period3 ,  ano , which is the time between consecutive pericentre passages.

2𝜋
obs 2 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  obs for this.When including  into radial velocity fits we want to use  ano as the period parameter to avoid the expected correlation between  obs and .Hence we need to be able to convert from one to the other.To do this we combine the two equations to get 2 and hence, neglecting terms of order O ( ) 2 , Our model fits for  obs as a parameter (i.e. the period prior is for  obs as is the output posterior distribution), but the model internally converts this to  ano .

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.

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)4 : where ,  , , and  are respectively the eccentricity, true anomaly, argument of pericentre, and inclination of the binary orbit relative to the plane of the sky;  1 and  2 are the semi-amplitudes of the primary and secondary, respectively, and  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: where  1 ,  2 , and  1 are the mass and radius of the primary and secondary (in solar units),  the orbital period (in days) (we use  ano ),  the true anomaly, and  0 = /2 −  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  2 ) and so does not warrant the increase in computation time.

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. 2009Konacki et al. , 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  parameter to account for differing zero-point calibrations5 .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 .
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  2 =  1 − .The planetary orbits are then fit in the same way to each dataset just as kima usually does.

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  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 • , 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  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  = 1.

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  = 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. 2009Konacki et al. , 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  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  as a function of  that envelops the lower 99% of the samples.Practically speaking, the limit is produced by creating log-linear bins along the  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  variables and jitters, as well as , , , and  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.

Analytic equation for precession
In appendix A we derive an analytic equation for the expected precession rate due to an outer perturber (Eq.(A11)), and due to rotational and relativistic effects (Eq.( 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 (Eq.( 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 Eq. (A5) should be used.

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 , and that including it can greatly improve the fit.
We perform two simulations, both with a primary star of mass  = 1 M , a secondary with  = 0.37 M ,  = 21.08 days, and  = 0.16 and a (roughly Jupiter mass) planet with  = 0.001  ,  = 134.5 days, and  = 0.01.The first simulation, SIM1, has 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 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.just these three bodies, whereas the second simulation, SIM2, has an additional planet with  pl = 0.00015 M ,  = 911.2days, and  = 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.

Improving scatter and derived parameters
First, we consider SIM1.From the values of  at each point in the rebound simulation, we obtain  = 308.4± 3.7 arcsec yr −1 for the binary.Using Eq. (A11) we get a theoretical value of  = 301.1 +0.8 −1.5 arcsec yr −1 A fit using the binaries model results in a posterior estimate of  = 304 ± 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"  at various times, calculating  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 th to the 84 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 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 RMS = 6.61 m s −1 to 3.15 m s −1 .When apsidal precession is not accounted for, if we move further from the reference time  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 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).

Detecting a hidden planet
Here, we consider SIM2 and test how many planets are formally detected.To register as an -planet detection, the Bayes Factor for the -planet solution compared to the ( − 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  +1, comparing the model with  + 1 planets to that with  planets, is the ratio of the number of posterior samples with  + 1 planets  +1 to the number with i planets Should   = 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.
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  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  2,1 > 6086 > 150.The Bayes Factors for each attempted fit are found in Table 2. 12.5 6086  3,2 0.9 0.8

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 days (Yenawine et al. 2022) which corresponds to  = 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  = 21 800 ± 600 arcsec yr −1 .Equivalently this is  = 0.56 ± 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  2  from 640.9 to 1.5.This amply justifies adding an extra parameter, , 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  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 ∼ 10 4 days 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 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 success-  fully detect an apsidal precession rate of  1 = 283 +87 −85 arcsec yr −1 , which is 3.3 from 0. Using Eqs.(A11, A17) we obtain a value of  1 = 92.4+14.3 −13.8 arcsec yr −1 , which is 2.2  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  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  obs and  ano .These values for Kepler-16 are presented in Table 3 alongside values published in Triaud et al. (2022).The value we get for  obs is in statistical agreement with the previous publication, however  ano is 3.3 above this. 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  ano and  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 ≈ 2 × 10 −6 M .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 ≈ 2 × 10 −4 M 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 .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.
We show the parameters from our fit of Kepler-16 in Table B2.

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

Figure 1 .
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 .
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.

Figure 3 .
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.

Figure 4 .
Figure 4. Kepler-16: red: histogram of the density of posterior samples for fitted value of  with the median and 1 values shown in grey.blue: histogram of the density of posterior samples for the theoretically calculated value of  with the median value shown in grey.(Note that  is not cut at zero, there are in fact posteriors below zero.)

Figure 5 .
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.

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.

Table 3 .
Triaud et al. (2022)ds for Kepler-16AB.The first two values (*) are taken fromTriaud 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)

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  2  ) are taken for the best fitting model.Excepting parameters with highly asymmetric distributions, the 1 uncertainties are shown as the last few significant digits.