# EBLM XV - Revised dynamical masses for the circumbinary planet host Kepler-16 AB, using the SOPHIE spectrograph

D. Sebastian, I. Boisse, A. Santerne, A. H. M. J. Triaud, T. A. Baycroft, Y. T. Davis, M. Deleuil, S. Grouffal, G. Hébrard, N. Heidari, et al.

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

## ► To cite this version:

D. Sebastian, I. Boisse, A. Santerne, A. H. M. J. Triaud, T. A. Baycroft, et al., EBLM XV - Revised dynamical masses for the circumbinary planet host Kepler-16 AB, using the SOPHIE spectrograph. Monthly Notices of the Royal Astronomical Society, 2025, 9, p. 729-740. 10.1093/mnras/staf863. insu-05101652

HAL Id: insu-05101652

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

Submitted on 4 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# EBLM XV – revised dynamical masses for the circumbinary planet host Kepler-16 AB, using the SOPHIE spectrograph

D. Sebastian<sup>1,★</sup> I. Boisse,<sup>2</sup> A. Santerne<sup>1,2</sup> A. H.M.J. Triaud<sup>1</sup> T. A. Baycroft<sup>1</sup> Y. T. Davis<sup>1</sup>

M. Deleuil<sup>1,2</sup> S. Grouffal,<sup>2</sup> G. Hébrard,<sup>3,4</sup> N. Heidari<sup>1,3</sup> D. V. Martin<sup>1,5</sup> P. F.L. Maxted<sup>1,6</sup>

R. P. Nelson,<sup>7</sup> S. Lalitha<sup>1,8</sup> M. G. Scott<sup>1</sup>,<sup>1</sup> O. J. Scutt<sup>1</sup> and M. Standing<sup>1,9</sup>

Accepted 2025 May 22. Received 2025 May 22; in original form 2025 April 7

## ABSTRACT

Eclipsing binaries are perfect laboratories to measure precise, accurate and model-independent stellar radii and stellar masses, so long as both components are spectroscopically resolved. Resolving both components is difficult in high-contrast binaries, for instance, those composed of an FGK main-sequence star with an M-type companion. In those cases, the secondary can contribute < 1 per cent of the total flux in optical wavelengths. This makes measuring dynamical masses challenging and has typically only been attempted with large-aperture telescopes (8–10 m). The High-Resolution Cross-Correlation Spectroscopy (HRCCS) method was developed to extract weak emission and transmission spectra for exoplanet atmospheres. This method was recently adapted and applied to measure dynamical masses in high-contrast binaries. In this work, we apply the HRCCS method to optical spectra of the high-contrast binary and circumbinary planet host Kepler-16 AB, obtained with the SOPHIE spectrograph at the 1.93-m telescope at the Observatoire de Haute-Provence. The secondary, which has a contrast ratio of ~6 × 10<sup>-3</sup>, is resolved with a detection significance of 9.5σ. We derive dynamical masses with a precision of 1.5 per cent and 0.9 per cent for the primary and secondary, respectively. These are comparable, but slightly higher (within 2–7 per cent), to previous mass-measurements, which has -within the uncertainties- no implication for the mass of the known circumbinary planet. This work demonstrates that dynamical mass measurements of high-contrast binaries can be done with 2-m class telescopes. We also investigate different analysis protocols to ensure we derive robust uncertainties for dynamical masses.

**Key words:** techniques: spectroscopic – planets and satellites: atmospheres – binaries: eclipsing – binaries: spectroscopic – stars: fundamental parameters – stars: low-mass.

## 1 INTRODUCTION

Very low-mass main-sequence stars ($M < 0.35 M_⊙$) with precise absolute mass and absolute radius measurements play a crucial role to calibrate stellar evolution models at the bottom of the main sequence. This incomplete view about the physical properties of these fully convective stars have far-reaching implications. For instance exoplanets' physical parameters (as well as their atmospheric parameters) directly depend on accurate stellar parameters, which are typically based on stellar evolution models. However, stellar evolution models do not all reproduce observations well. For

instance, the 'radius inflation problem' (e.g. Casagrande, Flynn & Bessell 2008; Torres, Andersen & Giménez 2010; Spada et al. 2013; Kesseli et al. 2018) causes stellar radii from evolutionary models to be under predicted by a few percent. This discrepancy is still not fully understood.

To empirically characterize very low-mass main-sequence stars, a large sample of eclipsing binaries is needed. However few low-mass double-lined eclipsing binaries are known (e.g. DEBCAT; Southworth 2015). Observations have instead focused on F, G, and K-type primaries that host eclipsing very low-mass stellar companions. Creating a mass-radius-luminosity-metallicity relation is the goal of the EBLM project (Eclipsing Binaries–Low Mass; Triaud et al. 2013; Maxted, Triaud & Martin 2023). However, by focusing on eclipsing single-lined binaries (SB1s), the mass and radius of the

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

<footnote n="2"><sup>2</sup>Aix Marseille Univ, CNRS, CNES, Institut Origines, Marseille LAM 13388, France</footnote>

<footnote n="3"><sup>3</sup>Institut d'astrophysique de Paris, UMR 7095 CNRS université pierre et marie curie, 98 bis, boulevard Arago, F-75014, Paris, France</footnote>

<footnote n="4"><sup>4</sup>Observatoire de Haute-Provence, CNRS, Université d'Aix-Marseille, 04870, Saint-Michel-l'Observatoire, France</footnote>

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

<footnote n="6"><sup>6</sup>Astrophysics Group, Keele University, Keele ST5 5BG, UK</footnote>

<footnote n="7"><sup>7</sup>Astronomy Unit, Queen Mary University of London, Mile End Road, London E1 4NS, UK</footnote>

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

<footnote n="9"><sup>9</sup>European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Madrid, Spain</footnote>

<footnote>*E-mail: D.Sebastian.1@bham.ac.uk</footnote>
M-dwarf companion are dependent on accurate estimates for the primary stars mass and radius. Despite highly reliable parameters, this means the EBLM systems are not included in catalogs such as DEBCAT (Southworth 2015), and thus not included to compare stellar models and observations.

EBLM binaries do have an important advantage. Their high-contrast between both stars in optical wavelengths (on average < 0.6 per cent; Triaud et al. 2017) mean (i) that the primary stars can be precisely characterized using standard spectroscopic analysis and stellar evolution models similar to single F, G, and K-stars (Freckelton et al. 2024), and (ii) that accurate and precise radial-velocities can be obtained to enable the detection of circumbinary exoplanets (e.g. Martin et al. 2019; Standing et al. 2022). Detecting a circumbinary planet typically requires dozens of high signal-to-noise spectra, enough that the secondary can be spectroscopically resolved via a careful analysis.

Sebastian, Triaud & Brogi (2024a, hereafter S24) recently used the HRCCS method (High-Resolution Cross-Correlation Spectroscopy; Snellen et al. 2010; Birkby 2018; Sebastian et al. 2024b) to recover the weak signal of the faint secondary star in the TOI-1338/BEBOP-1 circumbinary system. In this case, optical high-resolution spectra from ESPRESSO at the 8.2-m VLT (Pepe et al. 2021) were used, and precise dynamical masses (< 2 per cent) of the binary were obtained. Whilst regularly applied to detect and characterize exoplanet atmospheres (e.g. Brogi et al. 2012; de Kok et al. 2013), this was the first application of the HRCCS method to measure the absolute masses of a binary system. This method is highly effective at removing the primary's signal from the data (S24). This way the single-lined binary (SB1) was turned into a double-lined (SB2) eclipsing binary. Once combined with new, precise TESS data (Maxted et al. 2024) this enabled a full characterization of this binary system. It also validated the parameters obtained via the traditional SB1 approach normally used by the EBLM project (Davis et al. 2024).

In this paper, we turn our attention to Kepler-16AB, a *P* ~ 41 d EBLM binary that hosts the first unambiguously detected circumbinary exoplanet, Kepler-16 b (*P* ~ 229 d; Doyle et al. 2011). The planet was identified from transits with light curves obtained by the *Kepler* mission (Borucki et al. 2010). Kepler-16 also became the first circumbinary planet with a mass measured with radial velocities (Triaud et al. 2022). Thanks to the relatively late-type primary (K6), the system has already been measured as a double lined binary by Bender et al. (2012), who used six spectra, taken with the HRS spectrograph mounted on the 9.2-m Hobby-Eberly Telescope (Tull 1998). They applied a standard cross-correlation tool TODCOR (Zucker & Mazeh 1994) to resolve the weak secondary star. The HRS spectra covered the spectral range from 6100 - 8870 Å and the reported accuracies on the dynamical masses were 2.5 per cent on the primary and 1.5 per cent on the secondary, respectively.

In this paper, we (i) bench-mark the HRCCS method, using archival data of the well characterized binary system Kepler-16 AB, which were obtained by Triaud et al. (2022) and (ii) show that we can measure accurately and precisely dynamical masses of both stellar components using optical high-resolution spectra from a 2-m-class telescope.

The data selection and HRCCS analysis are described in Sections 3 and 4. We introduce an improved approach to robustly estimate uncertainties using the HRCCS method in Section 5 and extract phase- and wavelength dependent information of the faint secondary in Section 7.

## 2 SOPHIE DATA COLLECTION

Between 2016 August and 2023 November, 174 spectra have been obtained using the stabilised, fibre-fed high-resolution spectrograph SOPHIE at the 1.93-m telescope of the Observatoire de Haute-Provence in France (Perruchot et al. 2008). The SOPHIE data cover a slightly bluer spectral range of 3900–6900 Å, compared to the HRS spectral range (Bender et al. 2012). Hence, the expected contrast between both components (~ 0.6 %) is slightly larger.

The data were obtained in the high efficiency mode (R = 40 000), with a typical exposure time of 1800 s, and an average signal-to-noise-ratio (SNR) of 35 at 550 nm. The main fibre was placed on the star and the second fibre on the sky to allow optimal subtraction of the sky background. A radial velocity standard star has been observed every night and any possible wavelength drift during each night is monitored using a Fabry-Pérot etalon. All spectra were reduced using the SOPHIE data reduction pipeline (DRS Bouchy et al. 2009). The extracted radial velocities from the early part of this data-set have been used in Triaud et al. (2022, hereafter T22) to fit the primary's orbit and to detect the circumbinary planet Kepler-16 b. The data-set we analyse in this paper covers an additional 2.3 yr (27 new spectra).

## 3 DETRENDING AND CCF ANALYSIS

We use one dimensional 's1d' spectra, provided from the SOPHIE DRS. In these data products, all échelle orders of the spectrograph are already combined and the spectra are resampled to have a common resolution of 0.01 Å. We chose to use the s1d spectra since the combination of orders in the DRS has been optimized to reduce noise at the edges of the detector. For the data preparation, we follow closely the methods used in S24. We split the s1d spectra in 60 spectral chunks of 5119 points each and then normalize each spectral chunk -which are aligned in the primary's rest frame- using a polynomial of order 6. For this, we linearly interpolate the spectra using the orbital parameters from Triaud et al. (2022). We use this array of aligned spectra also to exclude columns with standard deviation more than 1.9-σ higher, compared to the average standard deviation in the normalized chunk. This value has been found experimentally for the data to exclude prominent outliers effectively without affecting the stellar absorption spectrum.

Due to the high stability of the SOPHIE spectrograph, the telluric absorption lines in each spectrum should be at constant velocity to within a few m s⁻¹. In a first step, we identify spectra unusable for the HRCCS processing by analysing the telluric drift of all spectra. We use a telluric line list compiled from the HITRAN data base (Gordon et al. 2022) including H₂O, OH, and CO₂ transitions which were normalised using the Cerro Paranal Sky Model for vacuum wavelengths (Noll et al. 2012; Jones et al. 2013). We cross-correlate the SOPHIE orders, which are most affected by telluric absorption lines between 6841 and 6893 Å, with this telluric line list to measure and extract the telluric velocity. We find that most SOPHIE spectra are well calibrated, except for a few outliers. These are the first 12 observations taken before November 2016. The extracted RVs were corrected in T22, but since in this work, we re-calculate the CCF for each spectrum, we opt to not include these spectra in the analysis. We also excluded two observations at BJD = 2458003.32 and 2458398.33 d, which show an telluric offset of more than 25 km s⁻¹. This selection results in 160 SOPHIE spectra, which are used for the HRCCS processing in this work.

In the next step, we remove the absorption spectrum of the primary star. Following S24, we apply a singular-value-decomposition (SVD,
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Upper panel: Detrending for SOPHIE data. Black dots: Effective rank, Grey lines, RMS of the residual arrays as a function of wavelength. Lower panel: Atmospheric transmission. Wavelength areas with large RMS in the upper panel, match well with strong telluric lines, which are less correlated in the primary's rest frame.</figcaption></figure>

Kalman 1996). This is applied directly to each normalized chunk consisting of 160 spectra, which are aligned in the rest frame of the primary star. The SVD decomposes this data into components, which are sorted in descending order according to their correlation in the primary's rest frame. In this way it effectively separates the spectral features of the primary from the data. We apply the 'effective rank' (Roy & Vetterli 2007) to select the optimum number of SVD components for each chunk. In order to handle noisy chunks, we set a maximum number of removed SVD components to 32, which is a fifth of the total number of components (160). The number of SVD components, removed in each chunks is shown in Fig. 1. The median number of components is 12.5, which is less than 8 per cent of this total number (reduced rank, see S24). With an average contrast of 0.6 per-cent, the SNR of the secondary star, per resolution element, in each SOPHIE spectrum, is on the order of 0.02. According to table 4 in S24, we thus expect that this SVD detrending has removed less than 5 per cent of the secondary's spectral features, and thus, has kept its signal mainly intact.

To detect the signal of the secondary in the detrended data, we first
derive a CCF with a line mask, optimized for M-dwarfs. Second,
we make use of the K-focusing process (Sebastian et al. 2024a),
which enhances the secondary’s weak CCF by combining it in the
secondary’s rest frame. Since the rest velocity $V_{\text{rest,2}}$, as well as
the semi-amplitude $K_2$ of the secondary are a priori unknown, this
combination is explored in the $K_2 - V_{\text{rest}}$ plane. The remaining
orbital parameters which are the period ($P$), time of periastron
($T_0$,peri), eccentricity ($e$), and the argument of periastron ($\omega$) are
adopted from T22 and kept fixed in this analysis. For the CCF,
we use a line mask for an M5 star, commonly used to extract
precise RVs within the data reduction pipeline of the ESPRESSO
spectrograph at the ESO-VLT<fnref n="1" /> (Pepe et al. 2021). The rest velocity of
the primary has been measured to be $V_{\text{rest,1}} = -33.8$ km s⁻<fnref n="1" /> by T22.
We therefore derive CCFs in a grid from -60 to 0 km s⁻<fnref n="1" />, with a step
width corresponding to the SOPHIE pixel resolution (1.5 km s⁻<fnref n="1" />)<fnref n="2" />
We cover $K_2$ in a grid from 25 to 75 km s⁻<fnref n="1" />, again in a step width of

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Upper panel: Cross-correlation map and *Saltire* model for SOPHIE data. Black dotted lines show the positions of slices (a,b), displayed in the lower panel. Lower panel: Slice through CCF map at maximum significance. Red line: best-fitting *Saltire* model, red shaded area: $1\sigma$ Uncertainties from the MCMC. Error bars represent the MCMC jitter term ($\sigma_{jit}$) from the two dimensional fit.</figcaption></figure>

1.5 km s⁻<fnref n="1" />. Before we combine all CCFs we exclude data at orbital phases where the velocity of the secondary’s absorption lines differs less than 10 km s⁻<fnref n="1" /> from the primary’s absorption lines. In this way, we exclude 19 spectra, making sure any remaining artefacts from the primary in the detrended data do not interfere with the radial velocity analysis.

Fig. 2 shows the resulting cross-correlation map in the $K_2-V_{\text{rest}}$ plane. This map is fitted using the *Saltire* model (Sebastian et al. 2024a) to derive the parameters $K_2$ and $V_{\text{rest}}$. The model fits the CCF using a K-focused double Gaussian function. This function models symmetrical side-lobes to the main CCF with five parameters. These are the mean height ($h$) of the CCF outside the signal, the standard deviations ($\sigma_1$ and $\sigma_2$) of two Gaussian functions with intensities $A_1$, $A_2$, the quotient $A_2/A_1$, as well as the sum $A_1 + A_2$ of both intensities. Similar to S24, we also sample the posterior probability distribution, using the Markov chain Monte Carlo (MCMC) code EMCEE (Foreman-Mackey et al. 2013), which is implemented in the *Saltire* model. The signal of the secondary is well detected with a detection significance of 9.5 $\sigma$. Fit parameters are presented in Table 1.

4 GLOBAL TIME FITTING

The main advantage of the K-focusing method is the possibility to
combine and detect very weak signals, such as planet atmospheres
or very faint M-dwarfs. The drawback of this method is that the
data in such CCF maps are highly correlated, making it difficult
to constrain uncertainties (Sebastian et al. 2024a). In the case of
Kepler-16 B, we can use the relatively high detection significance to
resolve the CCF trail directly in the time frame. For each detrended
spectrum, we first correct for the barycentric drift and compute a

<footnote n="1">¹ESO pipelines are publicly available on eso.org</footnote>

<footnote n="2">²SOPHIE instrument PSF</footnote>
<table><caption>**Table 1.** Priors and measured parameters from *Saltire* and global time MCMC fitting for detrended SOPHIE data. Errors of individual parameters are mean systematic uncertainties from time splitting, (Section 5.1)). Fit uncertainties are given in brackets for comparison.</caption><thead><tr><th colspan="4">Fixed parameters and priors</th></tr><tr><th>Parameters</th><th>Values</th><th>Saltire</th><th>Global fit</th></tr></thead><tbody><tr><td>P(d)</td><td>41.077772</td><td>fixed</td><td>±0.000051*</td></tr><tr><td>T<sub>0,peri</sub>(BJD)</td><td>2 457 573.0984</td><td>fixed</td><td>±0.0047*</td></tr><tr><td>e</td><td>0.15994</td><td>fixed</td><td>±0.00010*</td></tr><tr><td>ω(rad)</td><td>4.60194</td><td>fixed</td><td>±0.00051*</td></tr></tbody><thead><tr><th colspan="4">Measured parameters</th></tr><tr><th>Parameters</th><th>Saltire</th><th>Global fit</th></tr></thead><tbody><tr><td>P(d)</td><td>—</td><td>41.077777 ± 0.000035</td></tr><tr><td>T<sub>0,peri</sub>(BJD)</td><td>—</td><td>2 457 573.0995 ± 0.0032</td></tr><tr><td>e</td><td>—</td><td>0.15962 ± 0.00037</td></tr><tr><td>ω(rad)</td><td>—</td><td>4.60184 ± 0.00034</td></tr><tr><td>K<sub>2</sub>(km s<sup>-1</sup>)</td><td>47.16 ± 0.29 (0.08)</td><td>46.88 ± 0.28 (0.21)</td></tr><tr><td>V<sub>rest,2</sub>(km s<sup>-1</sup>)</td><td>-32.56 ± 0.30 (0.06)</td><td>-32.65 ± 0.22 (0.15)</td></tr><tr><td>A<sub>1</sub> + A<sub>2</sub></td><td>-(4.81 ± 0.61 (0.16) × 10<sup>-4</sup>)</td><td>-(5.42 ± 0.48 (0.36)) × 10<sup>-4</sup></td></tr><tr><td>A<sub>2</sub>/A<sub>1</sub></td><td>-0.421 ± 0.015 (0.017)</td><td>-0.363 ± 0.062 (0.057)</td></tr><tr><td>σ<sub>1</sub> (km s<sup>-1</sup>)</td><td>2.93 ± 0.23 (0.23)</td><td>2.85 ± 0.18 (0.22)</td></tr><tr><td>σ<sub>2</sub> (km s<sup>-1</sup>)</td><td>10.82 ± 2.6 (0.58)</td><td>8.9 ± 5.2 (1.2)</td></tr><tr><td>σ<sub>jit</sub></td><td>(5.1 ± 0.1) × 10<sup>-5</sup></td><td>(65.2 ± 0.4) × 10<sup>-5</sup></td></tr></tbody></table>

CCF with the same line mask for a M5 star – we used before – in a grid from –120 to 50 km s⁻¹, with a step width of 1.5 km s⁻¹. This ensures the secondary’s signal is covered in all orbital phases, taking the measured rest velocity and semi-amplitude into account. Similar to Section 3 we exclude 19 of the 160 spectra at orbital phases where the velocity of the secondary’s absorption lines differs less than 10 km s⁻¹ from the primary’s absorption lines. The left panel of Fig. 3 shows the resulting CCF functions as a function of the binary’s orbital phase. The CCF trail of Kepler-16 B is clearly resolved for a large part of the orbit. The CCF trail is relatively weak. The measured CCF contrast is smaller than the RMS noise of the CCFs (SNR<sub>trail</sub> = 0.8). Therefore, the CCF profile is hardly distinguishable from noise in individual CCFs.

Different to classical CCF fitting, which fits a Gaussian function to each individual CCF, we need to apply a global fit simultaneously to all CCF data, using a forward model of the CCF trail. All CCFs are simultaneously fitted to a double Gaussian function consisting of a main CCF signal as well as a side-lobe component, as used in the *Saltire* model. The position of the CCF signal in the CCF trail model is defined by its Keplerian movement following

$$V_{r,2} = K_2[\cos(\nu + \omega_2) + e \cos(\omega_2)], \quad (1)$$

where $K_2$ is the secondary’s semi-amplitude, and $\nu$ the true anomaly at the time of mid-exposure, which is determined from the orbital period $P$ and the time of periastron $T_{0,peri}$. Both were precisely measured by Triaud et al. (2022), who also obtained parameters for the eccentricity ($e$) and the argument of periastron of the primary ($\omega$). The argument of periastron of the secondary ($\omega_2$) is derived from the primary by $\omega_2 = \omega - \pi$. Different from the K-focusing method these Keplerian parameters are no longer kept fixed, and are now used as free parameters.

We sample the posterior probability distribution, of all free parameters using EMCEE. Given the low signal to noise of the CCF trail, we use tight uniform priors for the four Keplerian parameters $P$, $T_{0,peri}$, $e$, and $\omega_2$ using the uncertainties measured by T22. We use the fit of the *Saltire* model to obtain starting parameters, for $K_2$, $V_{rest}$, and the double Gaussian function which are listed in

Table 1. These are sampled with wide uniform priors over 6000 calls for 42 parallel chains, rejecting the burn-in samples (the first 1000 samples) of each walker and thin the remaining samples by a factor of 5. The resulting parameters and uncertainties are listed in Table 1. The middle panel of Fig. 3 shows the best-fitting model of CCF functions and the right panel the residual CCF’s. The relatively low SNR is reflected in a jitter term ($\sigma_{jit}$) which is about ten times higher compared to the K-focusing method. We find comparable parameters for the mean CCF shape as well as comparable uncertainties for the orbital parameters $K_2$ and $V_{rest}$ to the *Saltire* fitting. Before we compare the measured parameters, we analyse our results for systemic uncertainties in the following section.

## 5 UNCERTAINTIES FROM PARTIAL DATA

Different to the fitting of CCF maps in $K_p-V_{rest}$ plane from the K-focusing method, we expect the uncertainties derived from our global time fitting to be unaffected from correlated noise. We use the large data set of 141 (160-19) SOPHIE spectra to analyse the noise properties of our data set in two different methods and compare them with the uncertainties, derived from our MCMC fit. We employ two different methods to test our uncertainties, and ensure they are robust. When measuring absolute masses at the per-cent level, this exercise is an important step that ensures reliability.

### 5.1 Time splitting analysis

Similar to S24, we thus evaluate the uncertainties of the secondary fit by splitting the data set of 141 spectra in time into $n = 4, 6$, and 8 sets of 35, 23, and 18 spectra, respectively. We measure the parameters of each set (i) by fitting in the $K_2-V_{rest}$ plane using the *Saltire* code and (ii) by performing the global time fitting, introduced in Section 4. Each set is analysed individually and we estimate the systematic uncertainties as the standard error of the parameters of each individual set. Due to the low number of sets, we have to take the Bessel’s correction into account. The fractional uncertainty of deriving the standard error for such few sets ($n$) can be estimated as

<footnote>* From Triaud et al. (2022)</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Global time fitting of CCF data of Kepler-16 B. Left panel: CCF data from detrended SOPHIE spectra, folded to the orbital phase of the binary. Middle panel: Best fit double Gaussian model. Right panel: Residuals. Phases without data are presented as blank blue sections.</figcaption></figure>

<table><caption>**Table 2.** Systematic uncertainties of the M-dwarfs semi-amplitude ($K_2$) derived from time-splitting of the SOPHIE data in different number of sets. Right column: The statistical fractional uncertainty of the derived systematic uncertainty, due to the low number of sets. The mean systematic uncertainty ($\bar{X}_w$), is weighted by this fractional uncertainties.</caption><thead><tr><th>Number of sets</th><th>Saltire fit (km s<sup>-1</sup>)</th><th>Global time fit (km s<sup>-1</sup>)</th><th>Fractional uncertainty (%)</th></tr></thead><tbody><tr><td>4</td><td>0.23</td><td>0.21</td><td>41</td></tr><tr><td>6</td><td>0.27</td><td>0.31</td><td>32</td></tr><tr><td>8</td><td>0.35</td><td>0.30</td><td>27</td></tr><tr><td>X&#772;<sub>w</sub></td><td>0.29</td><td>0.28</td><td></td></tr></tbody></table>

$1/\sqrt{2n-2}$ (Topping 1972). As presented in Table 2, the fractional uncertainties are larger, when splitting the data in fewer sets. We thus have to take them into account when comparing the results for different numbers of sets.

The results for both fitting methods are presented in Table 2. We find that the systematic uncertainties are well in agreement (i) for the two different fitting methods and (ii) for a different number of sets, when taking the fractional uncertainties into account. We finally derive a mean systematic uncertainty ($\bar{X}_w$) from the different sets, which we weight by the inverse of the fractional uncertainties.

We find that this mean systematic uncertainty is only slightly larger (33 per cent), compared to the fit uncertainty from the global time

fitting (0.21 km s⁻¹, see Section 4). Taking the fractional uncertainty into account this a 1 σ agreement. In contrast, this mean systematic uncertainty it is about 4 times larger, compared to the fit uncertainty using *Saltire*. This means that, the global time fitting can be used to robustly estimate the uncertainties of the orbital parameters.

We conservatively use the mean systematic uncertainties of each parameter from the time splitting to estimate the uncertainties from the *Saltire* and global time fitting, as presented in Table 1. In brackets we add the fit uncertainties from the MCMC for comparison. Taking these uncertainties into account, we find that the semi-amplitude of the secondary is within a 1 σ agreement between the K-focusing and the global time fitting methods.

## 5.2 Bootstrapping analysis

Instead of splitting the sample as a function of time, we can also analyse randomly selected samples with a constant number of spectra. This method has been successfully applied to validate the measurement of exoplanet atmospheres, which are affected from correlated noise (e.g. Hoeijmakers et al. 2020). In this study, we select partial samples representing 50 per cent of the *m* = 141 spectra.

We draw 200 different partial samples arbitrarily, following the condition that partial samples of size *m* can only contain a certain amount of *s*<sub>max</sub> similar spectra to previous partial samples, with *s*<sub>max</sub> = 0.5 × (percentage/100) × *m* are allowed. With this, we en-
<table><caption>**Table 3.** Comparison of Masses for Kepler-16 AB.</caption><thead><tr><th colspan="2">Photometric-dynamical model</th></tr><tr><th>parameter</th><th>value</th></tr></thead><tbody><tr><td>q</td><td>0.2937 ± 0.0006<sup>†</sup></td></tr><tr><td>M<sub>1</sub>(M<sub>☉</sub>)</td><td>0.6897 ± 0.0035</td></tr><tr><td>M<sub>2</sub>(M<sub>☉</sub>)</td><td>0.20255 ± 0.00066</td></tr><tr><td>a(AU)</td><td>0.22431 ± 0.00005</td></tr><tr><th colspan="2">Dynamical masses from literature</th></tr><tr><th>parameter</th><th>value</th></tr><tr><td>K<sub>1</sub>(km s<sup>-1</sup>)</td><td>13.642 ± 0.010<sup>†</sup></td></tr><tr><td>K<sub>2</sub>(km s<sup>-1</sup>)</td><td>45.56 ± 0.47<sup>†</sup></td></tr><tr><td>q</td><td>0.2994 ± 0.0031<sup>†</sup></td></tr><tr><td>M<sub>1</sub>(M<sub>☉</sub>)</td><td>0.654 ± 0.017<sup>†</sup> (2.6 per cent)</td></tr><tr><td>M<sub>2</sub>(M<sub>☉</sub>)</td><td>0.1959 ± 0.0031<sup>†</sup> (1.6 per cent)</td></tr><tr><th colspan="2">Dynamical masses from SOPHIE data</th></tr><tr><th>parameter</th><th>value</th></tr><tr><td>K<sub>1</sub>(km s<sup>-1</sup>)</td><td>13.6787 ± 0.0015<sup>*</sup></td></tr><tr><td>K<sub>2</sub>(km s<sup>-1</sup>)</td><td>46.88 ± 0.28</td></tr><tr><td>q</td><td>0.2918 ± 0.0017</td></tr><tr><td>M<sub>1</sub>(M<sub>☉</sub>)</td><td>0.704 ± 0.011 (1.5 per cent)</td></tr><tr><td>M<sub>2</sub>(M<sub>☉</sub>)</td><td>0.2054 ± 0.0019 (0.9 per cent)</td></tr><tr><td>a(AU)</td><td>0.2257 ± 0.0010</td></tr></tbody></table>

Notes. § From Doyle et al. (2011)

†From Bender et al. (2012)

*From Triaud et al. (2022)

sure that individual drawn samples are different from each other. That means a drawn partial sample is not selected if it contain more than 17 spectra, which were drawn together in any previous partial sample. We derive the bootstrapping uncertainty for $K_p$ as the standard deviation of the fit results, for all 200 partial samples and correct for the reduced number of spectra by multiplying with $\sqrt{2}$. We finally repeat this experiment 3×, which results in an average uncertainty of 0.24 km s$^{-1}$ on $K_2$.

This uncertainty from 50 per cent bootstrapping is consistent with both the fit uncertainty from the global time fit and with the systematic fit uncertainty (Section 5.1). We therefore find that the bootstrapping is a good metric to estimate the uncertainties of noisy data (faint signals), which cannot be separated in several parts. Since we have shown in Section 5.1 that the global time fitting as well as fitting in the $K_2-V_{\text{rest}}$ plane (K-focusing) result in similar results for partial samples, the bootstrapping approach can also be applied to derive uncertainties for very noisy data, where the K-focusing is the only option.

## 6 DYNAMICAL MASSES

Combining the high precision measurement of the primary's semi-amplitude from T22 and the semi-amplitude of the secondary measured in this work, we can derive dynamical masses for both stars. Here, we use the inclination, measured from high-precision *Kepler* light-curves, which is $i = 90.3401 \pm 0.0019^{\circ}$ (Doyle et al. 2011). We follow the IAU recommended equations (table 3 in Prša et al. 2016) to derive the dynamical masses ($M_1, M_2$), the mass ratio ($q = M_2/M_1$), and semi-major ($a$) axis directly from the measured parameters as listed in Table 1.

The derived parameters are listed in Table 3 together with dynamical masses derived by Bender et al. (2012) and the photometric-dynamical model by Doyle et al. (2011). We find that the derived dynamical masses from SOPHIE data are on average $\sim 3\sigma$ larger,

compared to the literature, causing the binary system to be about 7 per cent more massive. This discrepancy could be caused by the low number of only six spectra, used by Bender et al. (2012). Comparing to the photometric-dynamical literature masses, the binary system appears to be only 2 per cent more massive, with the masses of both stars being compatible within $\sim 1.5\sigma$. This slight increase in mass would also affect the mass of the detected exoplanet. Given the reported uncertainties of the planet's radial velocity by T22, the planet's mass is still in perfect agreement to the reported mass, which is why we do not derive an updated planet mass in this work.

## 7 PHASE-CURVE AND LOW-RESOLUTION SPECTRUM

The detrended data allow us to actually extract the spectrum of the M-dwarf. This can be used to understand any phase-dependent intensity variations, or to identify spectral features.

We can extract the phase dependent CCF contrast from the CCF trail directly. The CCF contrast is a measure of the secondary contribution to the total light. Following S24, we derive the CCF trail in the rest frame of the secondary and fit each CCF in the time frame with a Gaussian function. The left panels of Fig. 4 show the extracted CCF contrast. Within the uncertainties the secondary's signal appears to be constant over the orbital phase, which is expected for a low mass star at visual wavelengths (for an exoplanet we would expect a day-night variation).

We combine all 141 spectra in the rest frame of the secondary and measure the combined CCF contrast for each of the 60 spectral chunks individually. In this way, we can extract the wavelength-resolved signal of the secondary, thus a low-resolution residual spectrum of the M-dwarf. The right panels of Fig. 4 show the residual spectrum. Since the CCF contrast is a function of the number of absorption lines and used line mask, the features do not necessarily correspond to individual spectral features, as seen in low-resolution spectra. This extracted spectrum can be compared to observed or to modelled spectra, to confirm its origin from the M-dwarf secondary.

To attempt this, we use a simulated PHOENIX model spectra (Husser et al. 2013). Stellar parameters of the primary and the radius of the secondary were published by Doyle et al. (2011). For the secondary, we use typical parameters for main-sequence stars (Pecaut & Mamajek 2013) to estimate the stellar effective temperature. We finally use PHOENIX model spectra with $T_{\text{eff},2} = 4400$ K and $\log g_{\star} = 4.5$ for the primary and $T_{\text{eff},2} = 3100$ K and $\log g_{\star} = 5.0$ for the secondary, assuming a solar metallicity ([Fe/H] = 0.0). Following S24, we correct the model wavelengths from vacuum to air, following Morton (2000), then we match the spectral resolution to SOPHIE ($R \sim 40000$) by convolution with a Gaussian kernel and correct the continuum for the radii of both stars, as measured by Doyle et al. (2011) from precise *Kepler* light-curve modelling ($R_1 = 0.6489 R_{\odot}$ and $R_2 = 0.22623 R_{\odot}$). For the final model, we used the secondary spectrum, corrected for the continuum of the primary. This noiseless model spectrum is separated into 60 spectral chunks, normalized, and cross-correlated in the secondary rest-frame to reproduce how the observations are treated. The measured CCF contrast of this model is shown in the right lower panel of Fig. 4. While the average CCF contrast of the model is smaller, compared to the SOPHIE data, dominant features -such as the increase of CCF contrast at 6300 Å- are present in the model as well.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Upper panels, CCF functions in the secondary's rest-frame. Left: as a function of the orbital phase, Right: as a function of the Wavelength. Lower panels, measurements of individual CCF contrasts (grey dots) and binned data (black dots). Left: Phase curve of the M-dwarf secondary, Right: residual flux as a function of wavelength, showing a low resolution spectrum of the M-dwarf (Black dots) and a residual spectrum of the secondary, derived from a PHOENIX model spectrum (Blue dots).</figcaption></figure>

## 8 CONCLUSIONS

Kepler-16 AB, is an eclipsing binary hosting the first circumbinary planet, which is detected with both high-precision light-curve (Doyle et al. 2011) and radial velocity measurements (Triaud et al. 2022). Intensive follow-up in the recent years made it one of the most studied circumbinary systems known to date. The faint M-dwarf secondary makes it difficult to measure the radial velocity of the secondary component and thus model-independent dynamical masses. This has only been possible using a 9-m class telescope (Bender et al. 2012) using standard cross-correlation methods.

In this study, we use HRCCS (High-Resolution Cross-Correlation Spectroscopy), commonly used to detect and analyse exoplanet atmospheres. This method has recently been qualified to measure precise dynamical masses of high-contrast binaries (Sebastian et al. 2024b). The power of this method is the use of a singular value decomposition (SVD), which can be used to remove the atmospheric lines of the primary component completely from the data, but keeps more than 95 per cent of the secondary's atmospheric lines intact. We apply HRCCS to 160 SOPHIE spectra to detect and measure the semi-amplitude of the faint secondary star, which is the first application of this kind to data of a 2-m class telescope. We use two

different fitting methods to analyse the cross-correlation functions (CCF's) of our data. The 'K-focusing', which is fitted using the Saltire model (Sebastian et al. 2024a) and a global time fit, introduced in this work. We find that both fitting methods result in comparable results.

We furthermore analyse our measurements for systematic uncertainties, by fitting partial SOPHIE data (i) by manually splitting the data in different sets and (ii) by using a bootstrap algorithm. We find that the uncertainties, derived from the global time fitting are well in agreement to the systematic uncertainties. We therefore conclude that the global fit is superior to the 'K-focusing' method, but requires a sufficiently strong CCF contrast to actually fit the CCF in the time domain ($SNR_{trail} = 0.8$ for Kepler-16 B in SOPHIE data). We furthermore find that the bootstrapping allows to estimate accurate uncertainties for those low SNR data which can only be fitted using the 'K-focusing' method.

As listed in Table 3, we derive precise dynamical masses of Kepler-16 AB with uncertainties of 1.5 per cent, and 0.9 per cent for the primary and secondary respectively. We furthermore find that the dynamical masses, derived with HRCCS, are slightly (2–7 per cent) higher for both binary components. This is well in agreement (~1.5σ) to the photometric-dynamical model (Doyle et al. 2011). For
the dynamical masses presented in Bender et al. (2012), we only find
a ~ 3 σ agreement, which might be caused by the low number of only
six measurements, which have been used in that study. Nevertheless,
the derived masses have, within the uncertainties, no implications
to the mass of the circumbinary planet published by (Triaud et al.
2022).

We have thus successfully confirmed the application of the HRCCS method to SOPHIE data and pushed down its application for a high-contrast binaries using an optical high-precision spectrograph on a 2-m class telescope. Thanks to the relatively late spectral type of Kepler-16 A (K6), the contrast ratio allowed to measure dynamical masses for this 0.2 M⊙ M-dwarf. We expect that the HRCCS method will be applicable to measure precise secondary masses of ~ 0.35 M⊙ and above for high-contrast EBLM binaries using the SOPHIE spectrograph.

ACKNOWLEDGEMENTS

We thank the French National Programme in Planetology (PNP)
for awarding the telescope time on the 193cm telescope at Ob-
servatoire de Haute-Provence (OHP) that led to this paper. Time
was granted under proposal IDs: 2024A_PNPS9, 2023B_PNPS3,
2023A_PNPS7, 2022B_PNPS2, 2022A_PNPS7, 2018B_PNPS3,
2019A_PNP001, 2018B_PNP001, 2018A_PNP002, 2017B_PNP001
& 2017A_PNP003 (PI Santerne) as well as the DDT programme
DDT_2016A (PI Triaud). We also thank the kind staff at OHP,
particularly the night time operators.

This research is funded from the European Research Council (ERC) through the European Union's Horizon 2020 research and innovation programme (grant agreement no. 803193/BEBOP), and from the ERC/UKRI Frontier Research Guarantee programme (grant agreement EP/Z000327/1/CandY). We received funding from the French Programme National de Physique Stellaire (PNPS) and the "Programme National de Planétologie" (PNP) of CNRS/INSU co-funded by CNES. MRS acknowledges support from the European Space Agency as an ESA Research Fellow.

DATA AVAILABILITY

The reduced spectra are available at the SOPHIE archive:
http://atlas obs-hp.fr/sophie/. Data obtained under Prog.IDs
16.DISC.TRIA, 16B.PNP.HEBR, 17A.PNP.SANT, 17B.PNP.SANT,
18A.PNP.SANT, 18B.PNP.SANT, and 19A.PNP.SANT.

REFERENCES

Bender C. F. et al., 2012, ApJ, 751, L31

Birkby J. L., 2018, preprint (arXiv:1806.04617)

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

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

Brogi M., Snellen I. A. G., de Kok R. J., Albrecht S., Birkby J., de Mooij E. J. W., 2012, *Nature*, 486, 502

Casagrande L., Flynn C., Bessell M., 2008, MNRAS, 389, 585

Davis Y. T. et al., 2024, MNRAS, 530, 2565

Doyle L. R., Carter J. A., Fabrycky D. C., et al., 2011, *Science*, 333, 1602

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

Freckelton A. V. et al., 2024, MNRAS, 531, 4085

Gordon I. E. et al., 2022, *J. Quant. Spec. Radiat. Transf.*, 277, 107949

Hoeijmakers H. J. et al., 2020, A&A, 641, A123

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

Jones A., Noll S., Kausch W., Szyszka C., Kimeswenger S., 2013, A&A, 560, A91

Kalman D., 1996, Coll. Math. J., 27, 2

Kesseli A. Y., Muirhead P. S., Mann A. W., Mace G., 2018, AJ, 155, 225

de Kok R. J., Brogi M., Snellen I. A. G., Birkby J., Albrecht S., de Mooij E. J. W., 2013, A&A, 554, A82

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

Maxted P. F. L., Triaud A. H. M. J., Martin D. V., 2023, *Universe*, 9, 498

Maxted P. F. L., Miller N. J., Sebastian D., Triaud A. H. M. J., Martin D. V., Duck A., 2024, MNRAS, 531, 4577

Mayor M. et al., 2003, The Messenger, 114, 20

Morton D. C., 2000, ApJS, 130, 403

Noll S., Kausch W., Barden M., Jones A. M., Szyszka C., Kimeswenger S., Vinther J., 2012, A&A, 543, A92

Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9

Pepe F. et al., 2021, A&A, 645, A96

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

Prša A. et al., 2016, AJ, 152, 41

Roy O., Vetterli M., 2007, in 2007 15th European Signal Processing Conference. p. 606

Sebastian D., Triaud A. H. M. J., Brogi M., 2024a, MNRAS, 527, 10921

Sebastian D. et al., 2024b, MNRAS, 530, 2572

Snellen I. A. G., de Kok R. J., de Mooij E. J. W., Albrecht S., 2010, *Nature*, 465, 1049

Southworth J., 2015, in Rucinski S. M., Torres G., Zejda M., eds, ASP Conf. Ser. Vol. 496, Living Together: Planets, Host Stars and Binaries. Astron. Soc. Pac., San Francisco, p. 164

Spada F., Demarque P., Kim Y. C., Sills A., 2013, ApJ, 776, 87

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

Topping J., 1972, Some Statistical Ideas. Springer Netherlands, Dordrecht, p. 29, https://doi.org/10.1007/978-94-011-6928-8_2

Torres G., Andersen J., Giménez A., 2010, A&A Rev., 18, 67

Triaud A. H. M. J. et al., 2013, A&A, 549, A18

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

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

Tull R. G., 1998, in D'Odorico S., ed., 3355, Proc. SPIE Conf. Ser. Vol. 3355, Optical Astronomical Instrumentation. SPIE, Bellingham, p. 387

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

APPENDIX A: RADIAL VELOCITY DRIFT OF
THE PRIMARY

In the precise analysis of T22 not only the planet, but also a cubic
drift on the order of ~ 20 m s⁻¹ has been reported in the reflex motion
of the binary, which is interpreted as part of the magnetic cycle of
the primary star. We thus explore the reported drift of the primary.
In this, we cross-correlate the 160 SOHPIE spectra with a line mask
for an K5 star, which is typically used to extract precise RVs within
the data reduction pipeline of the HARPS spectrograph (Mayor et al.
2003) at the ESO-3.6-m telescope. We then use our global time fitting
to precisely measure the primary’s orbit parameters. Here, we apply
an MCMC with wide uniform priors for all orbital parameters and
add a cubic drift with two additional parameters. Fig. A1 shows the
best-fitting cubic trend as well as the 1 σ uncertainties to the systemic
velocity and cubic drift. We fit our double Gaussian function to each
spectrum to extract the individual RVs. These RVs are less precise,
compared to the data products from the SOPHIE pipeline, which is
why we do not use them for the fit. We only use them for clarity in
Fig. A1, after the binary and planet orbits as reported in T22 have
been removed. These residual RVs show clearly that the drift -also
<figure><img src="image_8.png" /><figcaption>**Figure A1.** Trend from global time fitting of the primary. Red line: Best fit of the trend, using a cubic drift. Shaded area: $1\sigma$ uncertainties from the MCMC fit. Black dots: extracted RV residuals for clarity (not used for the fit).</figcaption></figure>

reported by T22- does not continue with the new data (after BJD = 2459459.36 d). Furthermore, the drift we measure from our global time fit is within the uncertainties not significantly different from zero. This supports the hypothesis by T22 that the initial drift is an artifact from the magnetic cycle (1900 to 2400 days) which is well covered by seven years of data, used in this study. As a consequence we fix the drift to zero in all measurements, presented in this study.

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