Fundamental effective temperature measurements for eclipsing binary stars – V. The circumbinary planet system EBLM J0608 − 59

EBLM J0608 − 59 / TOI-1338 / BEBOP-1 is a 12 th -magnitude, F9 V star in an eclipsing binary with a much fainter M dwarf companion on a wide, eccentric orbit (P=14.6 d). The binary is orbited by two circumbinary planets: one transiting on a 95-day orbit and one non-transiting on a 215-day orbit. We have used high-precision photometry from the TESS mission combined with direct mass measurements for the two stars published recently to measure the following model-independent radii: 𝑅 1 = 1 . 32 ± 0 . 02 𝑅 ⊙ , 𝑅 2 = 0 . 309 ± 0 . 004 𝑅 ⊙ . Using 𝑅 1 and the parallax from Gaia EDR3 we find that this star’s angular diameter is 𝜃 = 0 . 0309 ± 0 . 0005 mas. The apparent bolometric flux of the primary star corrected for both extinction and the contribution from the M dwarf ( < 0 . 4 per cent) is F ⊕ , 0 = ( 0 . 417 ± 0 . 005 ) × 10 − 9 erg cm − 2 s − 1 . Hence, this F9 V star has an effective temperature 𝑇 eff , 1 = 6031 K ± 46 K ( rnd . ) ± 10 K ( sys . ) . EBLM J0608 − 59 is an ideal benchmark star that can be added to the sample of such systems we are establishing for “end-to-end” tests of the stellar parameters measured by large-scale spectroscopic surveys.


INTRODUCTION
EBLM J0608−59 / TOI-1338 / BEBOP-1 is an F9 V star with an Mdwarf companion star that transits the primary star once every 14.6 days.The eclipses were found in ground-based photometry by the WASP survey (Wide Angle Search for Planets, Pollacco et al. 2006).Radial velocity measurements showed that the companion is too massive to be an exoplanet and so it was flagged as an EBLM system, i.e. an eclipsing binary with a low-mass companion (Maxted et al. 2023).Triaud et al. (2017) published a spectroscopic orbit for this SB1 (single-lined) binary system along with over 100 other EBLM systems identified by the WASP project.Additional high-precision radial velocity (RV) measurements were reported as part of the BE-BOP survey (Binaries Escorted By Orbiting Planets, Martin et al. 2019).These RV measurements combined with the RV measurements reported in Triaud et al. (2017) were used to set an upper limit  ≈ 1 Jup to the mass of any circumbinary planets with periods of a few hundred days orbiting this binary system.Nevertheless, transits from a circumbinary planet with a radius ≈ 7 ⊕ were discovered in the light curve of this star (TOI-1338) from the Transiting Exoplanet Survey Satellite (TESS) mission by Kostov et al. (2020), the first such system identified from TESS photometry.Standing et al. (2023) have reported the discovery of a second circumbinary planet, BEBOP-1 c, with a minimum mass of 65 M ⊕ and an orbital period of 216 days based on additional RV measurements obtained with the HARPS and ESPRESSO spectrographs (Pepe et al. 2002(Pepe et al. , 2021)).Standing et al. found that the transiting circumbinary planet, BEBOP-1 b, has a mass < ≈ 22  ⊕ based on their analysis of all the available high-precision RV data to-date.
Although the M-dwarf companion contributes less than 0.5 per cent of the flux at optical wavelengths, Sebastian et al. (2024) were able to detect the spectrum of this faint companion star by using over 100 spectra of EBLM J0608−59 obtained with the ESPRESSO spectrograph on the 8.2-m VLT telescope.This detection combined with the spectroscopic orbit of the F9 V primary star leads to a direct mass measurement for both stars, viz  1 = 1.098 ± 0.017 ⊙ and  2 = 0.307 ± 0.003  ⊙ .
In this study, we combine the masses of the two stars measured by Sebastian et al. (2024) with an analysis of the TESS light curve to obtain precise and accurate measurements for the F9 V primary and M-type secondary star in EBLM J0608−59.We then use the parallax of the system from Gaia DR3 (Gaia Collaboration et al. 2016Collaboration et al. , 2023) ) combined with flux measurements at optical and infrared wavelengths to measure directly the effective temperature (T eff ) of the F9 V star.This direct, accurate and precise measurement of T eff for a moderately bright star ( ≈ 12) in an eclipsing binary system that has a negligible flux contribution from its M-dwarf companion makes EBLM J0608−59 an ideal benchmark star that can be used for "end-to-end" tests of stellar parameters obtained from spectroscopic surveys, i.e. a direct comparison of the T eff and log  value published by large spectroscopic surveys to accurate, independent measurements of these quantities for a benchmark star observed in the same way as other stars in the survey.

TESS photometry
EBLM J0608−59 has been observed at 120-s cadence by TESS (Ricker et al. 2015) in 30 sectors over 3 years.The light curves available from the TESS Science Processing Operations Centre (SPOC) show systematic errors in some sectors that appear to be due to the transit interfering with the de-trending of the light curves.We therefore decided to produce our own light curves from the target pixel files (TPFs) using the package Lightkurve 2.0 (Lightkurve Collaboration et al. 2018). 1 The TPFs were downloaded from the Mikulski Archive for Space Telescopes2 (MAST) and light curves computed using the same target and background apertures used for the generation of the SPOC light curves.The main source of systematic error in these light curves appears to be inaccurate correction for the background flux.We therefore decided to use RegressionCorrector function in Lightkurve to perform a linear regression sector-bysector against the estimated background flux excluding transits and eclipses from the calculation of the model.This model was then evaluated at all times of observation in each sector and subtracted from the total flux in the target aperture.We then extracted the sections of the light curve within one full eclipse width of either a primary or secondary eclipse and divided each of these sections by a straight line fit by least squares to the data either side of the eclipse.These sections of the light curve were grouped into 49 subsets containing 1 or 2 primary eclipses and 1 secondary eclipse.
We performed a least-squares fit of the Nelson-Davis-Etzel (NDE) light curve model (Nelson & Davis 1972) independently to each of the 49 subsets of TESS photometry data using jktebop3 (Southworth 2010).The stars are very nearly spherical and we have not used data between the eclipses so the ellipsoidal effect was ignored in the analysis of the light curves.The reflection effect was also ignored.Limb darkening was modelled using the power-2 law recently implemented in jktebop (Southworth 2023).The values of the parameters ℎ 1 and ℎ 2 were estimated by interpolation within the relevant table from Maxted (2018).The effect of the assumed value of ℎ 1 for the primary star can be seen in the curvature of the light curve at the bottom of the primary eclipse so this parameter was allowed to vary in the analysis of the light curve.The effect of ℎ 2 on the light curve model is very subtle so this parameter was fixed at the value obtained from Maxted (2018).The assumed limb darkening of the secondary star has a negligible effect on the light curve so we used the fixed nominal values ℎ 1 = 0.7 and ℎ 2 = 0.4.
The free parameters in the least-squares fit to the light curve are: the sum of the stellar radii in units of the semi-major axis (fractional radii),  1 + 2 = ( 1 +  2 )/; the ratio of the stellar radii,  =  2 / 1 ; the ratio of the surface brightness at the centre of each stellar disc,  0 ; the orbital inclination, ; the time of mid-primary eclipse,  0 ;  sin() and  cos(), where  is the orbital eccentricity and  is the longitude of periastron for the primary star; "third light", ℓ 3 .The orbital period was fixed at the value  = 14.6085577.We also included priors on the values of  sin() = 0.137609 ± 0.000034 and  cos() = −0.072464± 0.000043 based on the spectroscopic orbit of the primary star from Sebastian et al. (2024).The value for the contamination of the photometric aperture provided in the meta data for the SPOC light curves is 1 per cent so we assume ℓ 3 = 0.01 ± 0.01 as a prior in the fit.

Eclipse ephemeris
We combined the 49 times of mid-primary transit from our analysis of the TESS photometry with 4 times of mid-primary transit determined from a least-squares fit to data from the WASP project to update our estimate of the orbital period.The data from the WASP project were fit with the same model as the TESS data but with only  0 ,  0 , ℓ 3 as free parameters, the other parameters being fixed at values determined from a fit to the TESS light curve.All the data from one observing season were fit together to determine each of the 4 times of mid-primary eclipse.All the times of mid-transit used in this analysis are provided in the supplementary material available with the on-line version of this manuscript.

Mass, radius and luminosity
The mass, radius of the two stars with their standard errors have been computed in nominal solar units (Mamajek et al. 2015) using a Monte Carlo simulation assuming independent Gaussian error on the values of  1 +  2 ,  =  2 / 1 ,  and  from the fit the TESS light curve and the values of  1 and  2 from Sebastian et al. (2024).The error on  = 14.60855755 d is assumed to be negligible.The same Monte Carlo simulation is used to calculate the luminosity of the two stars and their standard errors taking random values of the effective temperatures, T eff,1 and T eff,2 , sampled from the posterior probability distribution computed as described in the Section 2.4.

Effective temperature measurements
The effective temperature for a star with Rosseland radius  and total luminosity  is defined by the equation eff , where  SB is the Stefan-Boltzmann constant.For a binary star at distance , i.e. with parallax  = 1/, the flux corrected for extinction observed at the top of Earth's atmosphere is where  1 = 2 1  is the angular diameter of star 1, and similarly for star 2. All the quantities are known or can be measured for EBLM J0608−59 provided we can accurately integrate the observed Table 1.Observed apparent magnitudes for EBLM J0608−59 and predicted values based on our synthetic photometry.The predicted magnitudes are shown with error estimates from the uncertainty on the zero-points for each photometric system.The pivot wavelength for each band pass is shown in the column headed  pivot .The magnitudes of the primary G0V star alone corrected for the contribution to the total flux from the M-dwarf are shown in the column headed  1 .The flux ratio in each band is shown in the final column. Band flux distributions for the two stars independently.This is possible because photometry of the combined flux from both stars is available from near-ultraviolet to mid-infrared wavelengths, and the flux ratio in the TESS band is measured accurately from our analysis of the light curve.Although we have no direct measurement of the flux ratio at infrared wavelengths, we can make a reasonable estimate for the small contribution of the M-dwarf to the measured total infrared flux using empirical colour - eff relations.The M-dwarf contributes less than 0.5 per cent to the total flux so it is not necessary to make a very accurate estimate of the M-dwarf flux distribution in order to derive an accurate value of  eff for the F9 V primary star.
The photometry used in this analysis is given in Table 1.The Gaia photometry is from Gaia data release EDR3.J, H and Ks magnitudes are from the 2MASS survey (Skrutskie et al. 2006).WISE magnitudes are from the All-Sky Release Catalog (Cutri & et al. 2012) with corrections to Vega magnitudes made as recommended by Jarrett et al. (2011).Details of the zero-points and response functions used to calculate synthetic photometry for these surveys from an assumed spectral energy distribution are given in Miller et al. (2020).The u-and v-band photometry are taken from the SkyMapper survey DR4 (Onken et al. 2024).
To estimate the reddening towards EBLM J0608−59 we use the calibration of E(B−V) versus the equivalent width of the interstellar Na I D 1 line by Munari & Zwitter (1997).This is an semi-empirical relation based on an analytic relation between equivalent width and column density assuming a constant gas-to-dust ratio with the constant of proportionality calibrated using using 18 stars that show a single interstellar Na I absorption line.To measure EW(Na I D 1 ) we selected 69 of the spectra obtain with the ESPRESSO spectrograph observed through an airmass < 1.4.We first shifted these spectra into the rest frame of the primary star and then took the median value at each wavelength to obtain a high signal-to-noise spectrum of the F9 V primary star.We then divided each observed spectrum by this spectrum of the F9 V primary star after shifting it back to the barycentric rest frame.We then took the median of these residual spectra to obtain the high signal-to-noise spectrum of the interstellar features shown in Fig. 2. The equivalent width of the Na I D 1 line measured by numerical integration is EW(Na I D 1 ) = 14 ± 2 mÅ.This value is less than the values of EW(Na I D 1 ) for all the stars in the calibration sample of Munari & Zwitter (1997).To estimate the uncertainty on this value of E(B−V) we take the sample standard deviation for the 5 stars in the calibration sample with the lowest values of EW(Na I D 1 ) ≈ 250 mÅ(Fig.2).As can be seen in Fig. 2, this is quite pessimistic given the very low reddening we infer for EBLM J0608−59.Based on this analysis we obtain the estimate E(B−V) = 0.0036 ± 0.0018.We use this as a Gaussian prior in our analysis but exclude negative values of E(B−V).The co-added ESPRESSO spectra can be obtained from the supplementary online information that accompanies this article or from the corresponding author's website. 4o establish colour - eff relations suitable for dwarf stars with 3100 K < T eff < 3500 K we use a robust linear fit to the stars listed in Table 6 of Fouqué et al. (2018) within this  eff range.Photometry for these stars is taken from the TESS input catalogue.To estimate a suitable standard error for a Gaussian prior based on this fit we use 1.25× the mean absolute deviation of the residuals from the fit.Colour - eff relations suitable for the primary F9 V star were calculated in a similar way based on stars selected from the Geneva-Copenhagen survey (Holmberg et al. 2009;Casagrande et al. 2011) with 5750 K < T eff < 6250 K,  (B − V) < 0.05 and 3.5 < log  < 4.5.The results are given in Table 2.
The method we have developed to measure  eff for eclipsing binary stars is described fully in Miller et al. (2020).Briefly, we use emcee (Foreman-Mackey et al. 2013) to sample the posterior probability distribution (Θ|) ∝ (|Θ)(Θ) for the model parameters Θ with prior (Θ) given the data,  (observed apparent magnitudes and flux ratios).The model parameters are The prior (Θ) is calculated using the angular diameters  1 and  2 derived from the radii  1 and  2 and the parallax , the priors on the flux ratio at infrared wavelengths based on the colour -Teff relations in Table 2, and the Gaussian prior on the reddening described above.The hyper-parameters  ext and  ℓ account for additional uncertainties in the synthetic magnitudes and flux ratio, respectively, due to errors in zero-points, inaccurate response functions, stellar variability, etc.The parallax is taken from Gaia EDR3 with corrections to the zero-point from Flynn et al. (2022).
To calculate the synthetic photometry for a given value of  eff we used a model spectral energy distribution (SED) multiplied by a distortion function, Δ().The distortion function is a linear superposition of Legendre polynomials in log wavelength.The coefficients of the distortion function for star 1 are  1,1 ,  1,2 , . . ., and similarly for star 2. The distorted SED for each star is normalized so that the total apparent flux prior to applying reddening is  SB  2 T 4 eff /4.These distorted SEDs provide a convenient function that we can integrate to calculate synthetic photometry that has realistic stellar absorption features, and where the overall shape can be adjusted to match the observed magnitudes from ultraviolet to infrared wavelengths.This means that the effective temperatures we derive are based on the integrated stellar flux and the star's angular diameter, not SED fitting.
For this analysis we use model SEDs computed from BT-Settl model atmospheres (Allard et al. 2013) obtained from the Spanish Virtual Observatory. 5For the F9 V star appropriate we use the model with  eff,1 = 6000 K, log  1 = 4.0, [Fe/H] = 0.0 and [/Fe] = 0.0.For the reference SED for the M dwarf companion we assume  eff,1 = 3000 K, log  1 = 5.0, and the same composition.We experimented with distortion functions with 3, 4 and 5 coefficients per star and found the results to be very similar in all cases.The results presented here use three distortion coefficient per star because there is no improvement in the quality of the fit if we use a larger number of coefficients.The predicted apparent magnitudes including their un- certainties from errors in the zero-points for each photometric system are compared to the observed apparent magnitudes in Table 1.

RESULTS
The TESS light curve data and best-fit models for each subset are shown in Fig. 1.The mean and standard error of the mean across the 49 subsets for each free parameter and some useful derived quantities are given in Table 3.The results from our least-squares fit of a linear ephemeris to the 49 times of mid-primary eclipse from the TESS data and 4 times of mid-primary eclipse from the WASP data are shown in the same table.
The posterior probability distribution for the model parameters from our analysis to measure the stellar effective temperatures is summarised in Table 4 and the spectral energy distribution is plotted in Fig. 3.The random errors quoted in Table 4 do not allow for the systematic error due to the uncertainty in the absolute calibration of the CALSPEC flux scale (Bohlin et al. 2014).This additional systematic error is 10 K for the G0V primary star and 7 K for the M-dwarf companion.
The mass, radius and luminosity of the two stars derived from these results are given in Table 5.The two stars compared to other stars in eclipsing binary systems with accurate mass and radius measurements are shown in the Hertzsprung-Russell diagram in Fig. 4.

DISCUSSION
We used the software package bagemass version 1.3 (Maxted et al. 2015) to compare the parameters of the primary star, EBLM J0608−59 A, to a grid of stellar models computed with the garstec stellar evolution code (Weiss & Schlattl 2008).The methods used to calculate the stellar model grid are described in Serenelli et al. (2013).A Markov-chain Monte-Carlo method to explore the posterior probability distribution (PPD) for the mass and age of a star based on its observed  eff , luminosity, mean stellar density and surface metal abundance [Fe/H].Version 1.3 uses a Gaussian prior on the stellar mass rather than a power-law prior so that we can ensure that the mass derived is consistent with the observed value.We find a good fit to the observed properties of EBLM J0608−59 A for models with an age of 6.0 ± 0.3 Gyr assuming a mixing length parameter  MLT = 1.78 (solar value) and solar helium abundance or slightly enhanced helium abundance (+0.02 dex).The quality of the fit for a reduced mixing length ( MLT = 1.5) is significantly worse.The best-fit isochrone with an age of 5.8 Gyr assuming solar helium abundance and mixing length is shown in Fig. 4. Isochrones for the same age and initial metal abundance from the Dartmouth stellar evolution database (DSEP, Dotter et al. 2008) and the MESA Isochrones & Stellar Tracks (MIST, Choi et al. 2016) are also shown in Fig. 4. EBLM J0608−59 B matches very well the mass-radius relation for low-mass stars from the MIST and DSEP stellar model grids (Fig. 4).This is unusual for stars in this mass range, which tend to be larger than predicted by stellar models.These stars also tend to be cooler than predicted, which is also the case for EBLM J0608−59 B. These measurements are a valuable addition to the small sample of stars with independent mass and radius measurements around the mass range where there are discontinuities in the relationships between the observed properties of M-dwarfs (Jao et al. 2018;Rabus et al. 2019).These discontinuities are thought to be related to the transition between stars with a radiative core and fully-convective stars, but may also be associated with mixing of 3 He during merger of envelope and core convection zones that occurs for masses ≈ 0.34  ⊙ (MacDonald & Gizis 2018).It may then be worthwhile to improve on the precision of the T eff measurement presented here using observations of the eclipse depth at near-infrared wavelengths, i.e. around the peak of the M-dwarf's SED.Kostov et al. (2020) derived the following masses and radii for the stars in EBLM J0608−59:  1 = 1.127 +0.068 −0.069 ,  1 = 1.331 +0.024 −0.026  ⊙ ,  2 = 0.313 +0.011 −0.012  2 = 0.3089 +0.0056 −0.0060  ⊙ .The value of K 2 that we have used was not available at the time of that study so constraints on the primary star's radius and mass derived from the analysis of its spectral energy distribution and spectrum were required to obtain these mass and radius estimates.The agreement between the values from this study and our results is excellent.The better precision in the values we have derived is partly because we have a precise measurement of the stellar masses from Sebastian et al. (2024), and partly because we have been able to use light curve data from 3 cycles of observations with the TESS spacecraft whereas only 1 cycle of observations was available to Kostov et al..This gives us some reassurance that the parameters of M-dwarfs and other bodies orbiting Sun-like stars (e.g.brown dwarfs and planets) can be measured reliably using techniques such as those employed by Kostov et al..
We have investigated the impact of using the reddening estimate given in the catalogue published by Anders et al. (2022) instead of the prior on E(B−V) derived from the interstellar Na I D 1 absorption line.This reddening estimate was computed with the StarHorse code (Queiroz et al. 2018) using Gaia EDR3 spectrophotometry and broad-band photometry from various optical and infrared surveys.If we use the prior E(B − V) = 0.058 ± 0.036 based on the value of A V given by Anders et al., we obtain T eff,1 = 6214 ± 115 K, T eff,2 = 3184 ± 142 K.This is consistent with the effective temperatures given in Table 4 within the rather broad errors that result from the larger uncertainty in this reddening estimate.The value of T eff = 6272 ± 173 K given in the same catalogue is consistent with the value of T eff,1 obtained using the reddening estimate from the same source.However, the value of E(B−V) from Anders et al. corresponds to a value of EW(Na I D 1 ) = 180 +67 −128 mÅ, which is an order-of-magnitude larger than the value we observe for this star.It would certainly be worthwhile to obtain suitable data to calibrate the E(B−V) -EW(Na I D 1 ) relation shown in Fig. 2 for stars with little reddening, but this current calibration is clearly preferable to the weaker constraint on E(B−V) that can currently be obtained from spectrophotometry and broad-band photometry for nearby stars.

CONCLUSION
EBLM J0608−59 adds to a small but growing sample of eclipsing binaries for which we have directly measured the effective temperature of a solar-type star with a much fainter M-dwarf companion (Maxted et al. 2022;Maxted 2023).These stars are complementary to benchmark stars for which the effective temperature estimate is based on the angular diameters () measured using interferometry.Interferometric measurements of  are technically challenging and so suffer from systematic errors that can be seen in the poor agreement between measurements of the same star made with different instruments.For solar-type dwarf stars that typically have angular diameters for < ≈ 1 mas, these systematic can be ≈ 10 per cent or more (Karovicova et al. 2022;Soubiran et al. 2024).The effective temperature measurements we have made from eclipsing binary stars have the advantage that the sources of systematic error are small and well-understood, and the surface gravity of the stars is known to very high accuracy.In addition, these stars are within the magnitude range (V∼ 10) that can be observed directly by instruments used for large-scale spectroscopic surveys in their standard observing modes.This makes it feasible to conduct "end-to-end" tests of the accuracy of stellar parameters derived by these surveys.This is often not the case for the very bright stars that are currently accessible to interferometric measurements.

Figure 1 .
Figure 1.TESS photometry of EBLM J0608−59 as a function of orbital phase (blue points) and best-fit light curve models for these data fitted as 49 individual subsets (red lines).The residuals from the best-fit models are shown offset vertically above the light curve data.

Figure 2 .
Figure 2. Upper panel: The interstellar Na I D 1 absorption feature in the ESPRESSO spectrum of EBLM J0608−59.Vertical dashed lines shown the wavelength range used to measure EW(Na I D 1 ).Lower panel: EW(Na I D 1 ) -E(B−V) relation from Munari & Zwitter (1997) (solid line).The points are the calibrating data from Munari & Zwitter with E(B−V) < 0.2 that we use to assess the scatter in this relation at low reddening values.The red lines in the lower-left corner show the 1- range of our EW(Na I D 1 ) for EBLM J0608−59 and the inferred 1- range of E(B−V).

Figure 3 .
Figure 3.The spectral energy distribution (SED) of EBLM J0608−59.The observed fluxes are plotted with open circles and the predicted fluxes for the mean of the posterior probability distribution (PPD) integrated over the response functions shown in grey are plotted with filled symbols.The SED predicted by the mean of the PPD is plotted in dark blue and light blue shows the SEDs produced from 100 random samples from the PPD.The contribution to the total SED from the M dwarf (barely visible) is shown in orange.The TESS response function is shown in red.The W3 mid-infrared bands also used in the analysis is not shown here.

Figure 4 .
Figure 4. Upper panel: primary component of EBLM J0608−59 in the massradius plane.Middle panel: secondary component of EBLM J0608−59 in the mass -radius plane.Lower panel: both components of EBLM J0608−59 in the Hertzsprung-Russell diagram.The ellipses show 1- and 2- confidence regions on the parameters of EBLM J0608−59.All panels show isochrones for an age of 5.8 Gyr assuming [Fe/H] = 0.0 from bagemass (red solid line), DSEP (blue dotted line) and MIST (green dashed line).Cyan error bars show stars in eclipsing binary systems taken from DEBCat (Southworth 2015).

Table 3 .
Mean and standard error of the mean for free parameters and derived quantities from least-squares fits to 49 subsets of the TESS light curve of EBLM J0608−59.

Table 4 .
Results from our analysis to obtain the effective temperatures for both stars in EBLM J0608−59.