Fundamental effective temperature measurements for eclipsing binary stars -- III. SPIRou near-infrared spectroscopy and CHEOPS photometry of the benchmark G0V star EBLM J0113+31

EBLM J0113+31 is moderately bright (V=10.1), metal-poor ([Fe/H]$\approx-0.3$) G0V star with a much fainter M dwarf companion on a wide, eccentric orbit (=14.3 d). We have used near-infrared spectroscopy obtained with the SPIRou spectrograph to measure the semi-amplitude of the M dwarf's spectroscopic orbit, and high-precision photometry of the eclipse and transit from the CHEOPS and TESS space missions to measure the geometry of this binary system. From the combined analysis of these data together with previously published observations we obtain the following model-independent masses and radii: $M_1 = 1.029 \pm 0.025 M_{\odot}$, $M_2 = 0.197 \pm 0.003 M_{\odot}$, $R_1 = 1.417 \pm 0.014 R_{\odot}$, $R_2 = 0.215 \pm 0.002 R_{\odot}$. Using $R_1$ and the parallax from Gaia EDR3 we find that this star's angular diameter is $\theta = 0.0745 \pm 0.0007$ mas. The apparent bolometric flux of the G0V star corrected for both extinction and the contribution from the M dwarf ($<0.2$ per cent) is ${\mathcal F}_{\oplus,0} = (2.62\pm 0.05)\times10^{-9}$ erg.cm$^{-2}$.s$^{-1}$. Hence, this G0V star has an effective temperature $T_{\rm eff,1} = 6124{\rm\,K} \pm 40{\rm \,K\,(rnd.)} \pm 10 {\rm \,K\,(sys.)}$. EBLM J0113+31 is an ideal benchmark star that can be used for"end-to-end"tests of the stellar parameters measured by large-scale spectroscopic surveys, or stellar parameters derived from asteroseismology with PLATO. The techniques developed here can be applied to many other eclipsing binaries in order to create a network of such benchmark stars.


INTRODUCTION
Benchmark stars have properties that have been directly and accurately measured to good precision. They play a fundamental role in stellar astrophysics because we can only ascertain the accuracy and reliability of stellar models by comparing their predictions to the observed properties of real stars. Benchmark stars can also be used to establish empirical relations between stellar properties, e.g. colour -★ E-mail: p.maxted@keele.ac.uk effective temperature ( eff ) relations (Boyajian et al. 2013;Van Belle et al. 2021;Huang et al. 2015), or equations to estimate the mass or radius of a main-sequence star from eff , log and [Fe/H] (Torres et al. 2010). Empirical relations are particularly useful in cases where stellar structure models are known to be unreliable, e.g. for low-mass stars, where stellar models tend to under-predict the radius and overpredict eff (Spada et al. 2013;Cassisi & Salaris 2019;Zhou et al. 2014;Berger et al. 2006).
Considerable effort has put into calibrating the eff scale for FGKtype dwarf stars. In recent years, this effort has been partly driven by the need for accurate eff estimates for planet host stars in order to estimate their masses and radii using stellar models (Boyajian et al. 2015;Baines et al. 2009). Much of the progress in characterising exoplanets in recent years has been due to the improved precision in measuring stellar masses and radii (Jontof-Hutter 2019).
Benchmark FGK-type stars are also essential to calibrate the level of systematic and random uncertainties in massive spectroscopic surveys such as RAVE (Steinmetz et al. 2020), the Gaia-ESO survey (Gilmore et al. 2012), LAMOST (Deng et al. 2012), GALAH (Buder et al. 2018), etc. (Blanco-Cuaresma et al. 2014aHeiter et al. 2015;Jofré et al. 2018). These surveys aim to reconstruct the formation history of the Galaxy by studying the pattern of elemental abundances in stars as a function of their mass, age and kinematics. Jofré et al. (2019) in their recent comprehensive review of "industrial scale" stellar abundance measurements suggest that it is today not possible to know the temperature of a star better than an accuracy of 50 K. This uncertainty has a direct impact on reliability of trends observed in stellar abundance patterns between different stellar populations. Errors in eff are the dominant source of uncertainty when estimating the mass, radius, composition and age of a star (Valle et al. 2018;Jofré et al. 2019).
Validation and calibration of eff estimates for FGK-type dwarf stars currently rely on angular diameter measurements for a small sample of very bright stars such as Procyon, Cet, 18 Sco, Cen A, etc. (Karovicova et al. 2022;Heiter et al. 2015;Boyajian et al. 2013). Repeated measurements of the angular diameter for the same star often show differences much larger than the quoted uncertainties, with systematic errors of 5 per cent or more being quite common. For example, the 15 repeated measurements provided in Table 9 of Karovicova et al. (2022) for 7 G-type dwarf stars require an additional "external error" of about 0.04 mas to be added to the quoted uncertainties to achieve 2 = 1 for a fit of a constant offset to these difference. Tayar et al. (2022) show that current uncertainties on measured interferometric angular diameters and bolometric fluxes set a systematic uncertainty floor of ≈ 2 per cent in eff for solar-type exoplanet host stars, i.e. ±120 K at eff = 6000 K.
Very low-mass stars (VLMSs, < ≈ 0.2 ) are attractive targets for exoplanet studies because of the possibility to detect and characterise the atmospheres of terrestrial planets in the habitable zones of these stars (Sebastian et al. 2021). There are very few well-characterised VLMSs because they are intrinsically very faint and small. For example, the recent empirical coloureff , colour -luminosity and colour -radius relations published by Boyajian et al. (2012) are based on a sample that contains only one star with a spectral type later than M4V ( ≈ 0. 2 , Mann et al. 2019). The EBLM project (Triaud et al. 2013) aims to improve our understanding of VLMSs by studying eclipsing binaries with low-mass companions that have been found by the WASP survey (Pollacco et al. 2006). These eclipsing binaries typically have a late-F-to mid-G-type primary star with an M-dwarf that contributes 1 per cent of the flux at optical wavelengths. The light curves of these EBLM systems look very similar to those of transiting hot Jupiters, which are the main targets for the WASP survey. As a result, dozens of these EBLM systems have been identified in the WASP survey. The analysis of the light curve combined with a spectroscopic orbit for the primary star and an estimate for its mass provides a direct estimate for the mass and radius of the M-dwarf companion (Von Boetticher et al. 2019;Gill et al. 2019). With very high quality photometry it is possible to detect the eclipse of the M-dwarf and, hence, its surface brightness relative to the primary star. This surface brightness ratio can then be used to infer eff for the M-dwarf given an estimate of eff for the primary star and a surface brightnesseff relation for the stars, either empirical (Graczyk et al. 2021) or based on stellar model atmospheres. The first results from an on-going programme to measure mass, radius and eff for the M-dwarf in a sample of EBLM systems using ultrahigh-precision photometry obtained as part of the guaranteed time observations (GTO) with the CHEOPS mission (Benz et al. 2021) have been published by Swayne et al. (2021). Most of the targets for this programme were selected from a sample of over 100 EBLM systems with spectroscopic orbits published by Triaud et al. (2017).
The first study of EBLM J0113+31, the target for this study, was published by Gómez Maqueo Chew et al. (2014. They used ground-based photometry of the eclipse in the J-band to infer eff ≈ 3900 K for the very low-mass companion, much higher than expected given their estimate for its mass ( 2 = 0.186 ± 0.010 ). Subsequent analysis of the TESS light curve for this binary system by Swayne et al. (2020) found no evidence for a very hot companion. Their value of eff,2 = 3208 ± 43 K is similar to that for other VLMSs. They conclude that the anomalous result from GMC+2014 was due to systematic errors in the J-band photometry, illustrating the need for very high-quality space-based photometry to make reliable measurements of eff,2 in EBLM systems.
Here we present new photometry of the transit and eclipse in EBLM J0113+31 obtained with CHEOPS, and high-resolution, phase-resolved spectroscopy obtained with the near-infrared échelle spectrograph SPIRou on the Canada-France-Hawaii telescope. We have used the SPIRou spectroscopy to directly measure the semiamplitude of the M-dwarf's spectroscopic orbit. We have used this measurement combined with the analysis of the new light curves and other published data to directly and accurately measure the mass, radius and eff of both stars in this binary system. We discuss the use of the techniques developed here to determine fundamental stellar properties for stars in EBLM systems, and conclude that it is now feasible to establish a network of well-studied EBLM systems that will be an ideal set of benchmark stars for future surveys.

CHEOPS photometry
CHEOPS is a telescope with an effective aperture of 30 cm in low Earth orbit that is designed to obtain ultrahigh precision broadband photometry of bright stars (Benz et al. 2021). To achieve this, the instrument has been purposely defocused to produce a point-spread function (PSF) with a diameter of approximately 32 . We observed two transits and one eclipse of EBLM J0113+31 with CHEOPS (Table 1). There are gaps in the observations due to occultation of the target by the Earth and passages of the satellite though the South Atlantic Anomaly.
The raw data were processed using version 13.1.0 of the CHEOPS data reduction pipeline (DRP, Hoyer et al. 2020). The DRP corrects the images for environmental and instrumental effects before performing aperture photometry of the target. The contamination of the photometric aperture during the exposure by nearby stars is estimated using simulations of the field of view based on the Gaia DR2 catalogue (Gaia Collaboration et al. 2018). The instrument response function for CHEOPS is very similar to the Gaia G band. The detector used on the CHEOPS instrument is a frame-transfer charge-coupled device (CCD), so the DRP must also account for the "smear" trails from bright stars produced during the frame transfer. Both of these effects (contamination and smear) vary from image to image because the satellite rotates continuously during its 99-minute orbit. Aperture photometry is extracted by the DRP using three different fixed aperture sizes labelled "RINF", "DEFAULT" and "RSUP" (at radii of 22.5, 25.0 and 30.0 pixels, respectively) and a further "OPTIMAL" aperture whose size is determined for each visit dependent upon the amount of contamination. The observed and processed data are made available on the Data Analysis Center for Exoplanets (DACE) web platform 1 . We downloaded our data from DACE using pycheops 2 , a module developed for the analysis of data from the CHEOPS mission (Maxted et al. 2021).
There are three stars that are 5-6 magnitudes fainter than EBLM J0113+31 within 1 of the target (Fig. 1). As a result, the OPTIMAL aperture is set to its maximum allowed value by the DRP (40 pixels = 40 ). Although this maximises the contamination of the aperture by these nearby stars, it minimises the noise due to the variations in this contamination due to changes in the fraction of the stars' PSFs inside the photometric aperture as the field of view rotates.
We used spectra extracted from the raw data provided by the observatory using data reduction system (DRS) version 0.6.131. We dealt with the data order-by-order, selecting only those orders with little contamination due to telluric features. The selected wavelength regions are listed in Table 2. The corrections for the échelle blaze function and telluric absorption provided by the DRS were applied.

TESS photometry
One transit and two secondary eclipses of EBLM J0113+31 were observed at 120 s cadence by the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) in Sector 17 of the primary mission. The TESS target pixel files were downloaded from the Mikulski Archive for Space Telescopes 3 (MAST) and processed to produce a light curve using the package L 2.0 (Lightkurve . 4 The pixels used to extract the photometry from the target pixel file are shown in Fig. 2. Instrumental noise was removed using the cotrending basis vectors (CBVs) provided by the TESS Science Processing Operations Center (SPOC) (Jenkins et al. 2016). We used 16 "Single-Scale" and 7 "Spike" CBVs to model trends present in all targets on the same CCD as EBLM J0113+31. The amplitude of each CBV was determined using only data outside the eclipses and transit. We set the L2-norm penalty to = 0.1 to achieve a balance between over-fitting the data and effectively removing instrumental trends.

Radial velocity measurements from the SPIRou data
We use synthetic spectra taken from Husser et al. (2013) 5 to produce a template for the spectrum of the G0V primary star, using linear interpolation to create a spectrum appropriate for eff = 6150 K, [Fe/H] = −0.4, log = 4.15 and [ /Fe] = 0.0. We then measure the position of the peak in the cross-correlation against this template for the observed spectra order-by-order. Low-frequency noise in the data for each order was removed prior to cross-correlation using a 5-th order high-pass Butterworth filter with a critical frequency of 16/4096 pixels −1 . We then reject measurements more than 5 km s −1 from the median before calculating the mean and standard error in the mean given in Table 3.

Pre-processing of the SPIRou data
The M dwarf contributes less than 2 per cent of the flux at 1 m so we removed the spectral features in the SPIRou data due to the G0V primary star prior to our attempt to detect the faint companion in these spectra. We use the spectroscopic orbit from GMC+2014 to shift the template spectrum for the primary star to the radial velocity corresponding to the time of mid-exposure for each SPIRou spectrum and then divide the observed spectrum by the shifted model spectrum.
The correction for telluric absorption in the observed spectra will be imperfect so we mask pixels where the telluric absorption is greater than 50 per cent. We also mask all pixels in order 47 at wavelengths > 1616 nm because there is a strong telluric absorption band at these wavelengths. The removal of spectral features from the primary star will also be imperfect so we mask pixels where absorption lines in the template spectrum are deeper than 50 per cent. We then flatten the spectrum by dividing the data by a 16th-order polynomial fit by least-squares to the unmasked data in each order.
Outliers due to cosmic ray hits on the detector and other image anomalies were then identified and removed by flagging pixels more than 4 times the inter-quartile range from the mean in 10 blocks of data per order.
The signal-to-noise is similar for each spectrum but varies quite strongly with wavelength so we use 1.25× the mean absolute devia-tion of the data across the observed spectra to assign a standard error to the pixels at each wavelength.

Detection of the M dwarf in the SPIRou spectra
The signal from the M dwarf is too weak to be detected in the individual SPIRou spectra, but it is possible to measure the semi-amplitude of M dwarf's spectroscopic orbit, 2 , by calculating the average cross-correlation function against a suitable template spectrum after shifting these CCFs to the rest frame of the M dwarf assuming a range of 2 values. The barycentric radial velocity of the M dwarf at the time of mid-exposure for each spectrum is The value of the eccentricity and the longitude of periastron 2 = 1 + are known accurately from the spectroscopic orbit of the primary star with longitude of periastron 1 , taken from GMC+2014. Similarly, the true anomaly at the time of mid-exposure, , can be accurately predicted from the values of 0 (time of mid-transit), (orbital period), and 1 , also taken from GMC+2014.
We use synthetic spectra taken from Husser et al. (2013) as a template for the spectra of the M dwarf, using linear interpolation to create a spectrum appropriate for eff = 3300 K, [Fe/H] = −0.4, log = 5.0 and [ /Fe] = 0.0. The cross-correlation function is calculated order-by-order. Low-frequency noise in the data for each order was removed prior to cross-correlation using a 5-th order highpass Butterworth filter with a critical frequency of 32/4096 pixels −1 . The data are apodized using a Gaussian filter with a standard deviation of 64 pixels applied to the data at each end of the order. The correlation coefficient for each order is calculated after shifting the template according to radial velocity computed with equation (1) includes the weights calculated from the estimated standard errors on each pixel. This is repeated for a uniform grid of 2 values. The average CCF over all orders and all exposures as a function of 2 ("stacked CCF") is shown in Fig. 3. Gómez Maqueo Chew et al.
(2014) estimate that 2 = 80.3 ± 1.5 km s −1 . There is indeed a peak in the stacked CCF near this value of 2 . To measure the position of this peak we model the stacked CCF as a Gaussian process (GP) plus a Gaussian profile. We use the package (Foreman-Mackey et al. 2017) to compute the likelihood for a GP with a kernel of the form ( ) = − and the affine-invariant Markov chain Monte Carlo sampler (Foreman-Mackey et al. 2013) to sample the posterior probability distribution for the model parameters. Based on this analysis, the peak in the stacked CCF occurs at 2 = 82.9 ± 0.7 km s −1 and has a width of 5.7 ± 0.6 km s −1 . The broad peak in the stacked CCF around 2 = 0 is due to imperfect removal of telluric features and spectral features from the G0V primary star. We compared the stacked CCF to the average CCF computed with negative values of 2 plotted against | 2 |, i.e. the mirror image of the stacked CCF. As can be seen in Fig. 3, there is no corresponding peak at 2 ≈ −83 km s −1 . This reassures us that the peak at 2 ≈ +83 km s −1 is unlikely to be due to imperfect removal of telluric features or spectral features from the G0V primary star. We used a fit to the stacked CCF done in the same way as above but excluding data around the peak at 2 = 83 km s −1 to estimate the statistical significance of this feature. Based on the GP prediction of the correlated noise in this region shown in Fig. 3, we estimate that the peak height corresponds to a detection with a significance ≈ 4-. We also verified that the height of the peak in the stacked CCF is very close to the height expected for an M-dwarf companion given the flux ratio ℓ ≈ 0.00155 inferred from the depth of the secondary eclipse in the TESS light curve. We subtracted a scaled version of the template M-dwarf spectrum from the spectra used to compute the stacked CCF based on this flux ratio in the TESS band and recomputed the stacked CCF. As can be seen in Fig. 3, the resulting stacked CCF has no peak near 2 = 83 km s −1 . Based on these three tests, we are confident that our detection of the M-dwarf is robust and that the measurement of 2 is reliable.

Initial assessment of the CHEOPS data
We used the software package (Maxted et al. 2021) to make an initial assessment of the light curve data from each of the three CHEOPS visits to the target listed in Table 1. We excluded data from the analysis where the background level in the images due to scattered light is more than 20 per cent larger than the median value during the visit. The data file provided by the data reduction pipeline (DRP, Hoyer et al. 2020) includes a quantity LC_CONTAM that is an estimate of the contamination of the photometric aperture by nearby stars. This quantity varies during the orbit because of the rotation of the spacecraft and the strongly asymmetric point spread function of the instrument. This calculation of LC_CONTAM is based on a simulation of the field of view using the mean G-band magnitudes of the target and nearby stars from Gaia DR2. Fig. 1 shows the results of this simulation for one image. This contamination estimate does not account for variability of target itself, so we added a new function to version 1.0.2 that corrects the measured flux (FLUX) by subtracting the value LC_CONTAM × 10 −0.4( − 0 ) from FLUX. The zero-point 0 is calculated from the average value of where frac is the fraction of the total flux from the target in the photometric aperture, is the mean G-band magnitude of the target, and the average is taken over data points outside of transit and eclipse.
Based on the simulated image of the field of view shown in Fig. 1, we decided to use the OPTIMAL aperture with a radius of 40 pixels for our analysis. This maximises the contamination of the photometric aperture but minimises the uncertainty in this quantity due to errors in measuring the positions of the stars in the image and, hence, the fraction of the flux from each star that is contained in the aperture. We repeated our analysis using the DEFAULT aperture with a radius of 25 pixels and found that the results are entirely consistent with those presented here. For each visit we calculate a best-fit for a transit or eclipse model to the light curve including linear decorrelation against LC_CONTAM to account for small errors in estimating the amplitude of the variations in this quantity. We then calculate the best-fit light curves including each of the other available decorrelation parameters and add them one-by-one if the Bayes factor for the parameter exceeds 1. The decorrelation parameters selected by this method are listed in Table 1. The fit to the data from a typical visit including these detrending parameters in shown in Fig. 4.

Updated transit time ephemeris
The two times of mid-transit measured from the CHEOPS data during the initial assessment of the data described above are listed in Table 4 together with the time of mid-transit from GMC+2014 and one new time of mid-transit from a least-squares fit to the TESS light curve using the transit model from . From a least-squares fit to these data we obtain the following updated ephemeris for the times of mid-transit in EBLM J0113+31:   There is no evidence for any change in orbital period greater than | / | ≈ 1 × 10 −5 from these times of mid-transit.

Combined analysis of light curve and radial velocity data
We used the light curve model ellc (Maxted 2016) to calculate synthetic light curves in the TESS and CHEOPS bands, and the spectroscopic orbit of the primary star. This model gives us more flexibility in choosing the level of numerical noise in these synthetic light curves than is possible with the qpower2 algorithm used in (Maxted & Gill 2019). For the analysis presented here we used the "default" grid size for the primary star and the "very_sparse" grid size for the companion, which gives numerical noise of only a few ppm at most orbital phases and everywhere less than 10 ppm. We also tested for the impact of the gravitational distortion of the stars by their mutual gravity on the light curve. This is less than 1.5 ppm through the transit so we assumed spherical stars for our analysis in order to speed-up the calculation.
The parameters of the binary star model are: the radii of the stars in units of the semi-major axis (fractional radii), 1 = 1 / and 2 = 2 / ; the surface brightness ratios in the TESS and CHEOPS bands, T and C , respectively; the orbital inclination, ; the time of mid-transit, 0 ; the orbital period, ; = √ sin( ) and = √ cos( ), where is the orbital eccentricity and is the longitude of periastron for the primary star; the semi-amplitude of the primary star's spectroscopic orbit, 1 ; the limb-darkening parameters assuming a power-2 limb-darkening law, ℎ 1,T and ℎ 2,T in the TESS band, and ℎ 1,C and ℎ 2,C in the CHEOPS band. The ephemeris for the time of mid-transit derived in Section 3.5 is very accurate so we fix 0 and at these values in our analysis. The curvature of the light curve between the second and third contact points is very clearly seen in the CHEOPS and TESS light curves, and is almost directly related to the parameters ℎ 1,C and ℎ 1,T , respectively, so we leave this as a free parameter in the analysis. The parameters ℎ 2,C and ℎ 2,T will have a much more subtle influence on the light curve that is almost entirely confined to the ingress and egress phases so we impose priors on these parameters based on the tabulated values of ℎ 2 in the TESS and CHEOPS bands from Maxted (2018). The width of the priors is based on the comparison of these tabulated values to the observed values of this parameter from an analysis of the Kepler light curves of transiting exoplanet systems in the same study.
Prior to the analysis of the CHEOPS data combined with the other data sets, we applied a correction for hot pixels in the photometric aperture. Quantitatively, we define hot pixels as pixels with dark current above 3 e − s −1 . Since the beginning of the mission, hot pixels have appeared regularly in the CHEOPS CCD at a rate of ∼100 new hot pixels per day. The CHEOPS Instrument Team monitors closely the number and location of hot pixels. Approximately once per week, "dark images" are acquired for that purpose (10 full frame images obtained observing a region of the sky void of stars). These images are used to produce the reference files that track the location and dark current of hot pixels. These reference file are available from the CHEOPS data archive. 6 We used hot pixel maps generated about 2 days after each visit to EBLM J0113+31 to calculate the contribution of these hot pixels to the count rate in the photometric aperture. The hot-pixel contamination is ≈ 0.6 per cent in the OPTIMAL aperture for the visit in 2020 and ≈ 1.2 per cent for the visits in 2021. The hot-pixel contamination in the DEFAULT aperture is ≈ 0.3 per cent for all visits. The hot pixel contamination is calculated for every image but the variation in this quantity is small ( < ≈ 10 per cent of its value) so we apply the correction by subtracting the mean value of the contamination during the visit from the count rate.
Our model includes the parameter ℓ 3,C that is a constant added to the synthetic CHEOPS light curve to account for contamination of the photometric aperture. We applied a correction to the light curves for contamination prior to the combined analysis so, to account for uncertainties in these corrections, we assign a Gaussian prior to ℓ 3,C with a mean value of 0 and a standard error equal to 50 per cent of the total contamination estimate. Similarly, the parameter ℓ 3,T accounts for the contamination of the photometric aperture shown in Fig. 2 used to extract the TESS light curve. We noticed that the entry TIC 400048098 in the TESS input catalogue (TIC, Stassun et al. 2019) has no counterpart in GAIA EDR3 (Gaia Collaboration et al. 2016 so we assume that this is a spurious entry and do not include it in our calculation of the contamination. The star TIC 400048094 appears near the edge of the default photometric aperture provided with the target pixel file. We added one pixel to this aperture so that there is no ambiguity over whether this star should be included in the calculation of the contamination or not. From the magnitudes listed the TIC we estimate ℓ 3,T = 0.0030. We allow this parameter to vary in the fit but assign a Gaussian prior to it, assuming an arbitrary uncertainty of 50 per cent.
The light curves produced by CHEOPS are known to have very low levels of instrumental noise after decorrelation. Similarly, the TESS light curve following correction for instrumental trends that we calculated with shows little sign of residual instrumental noise or stellar variability. We therefore adopt a white noise model 6 https://cheops-archive.astro.unige.ch/archive_browser/ for our analysis and assume that the standard deviation per point in the TESS and CHEOPS light curves -T and C , respectively -are the same for all data points from the same instrument. The logarithm of these standard errors are included as a hyperparameters in our analysis by correctly normalising the calculation of the posterior probability distribution. We only include data from the TESS light curves at orbital phases near the transit and eclipses in this analysis. For both the CHEOPS and TESS data, each section of data around a transit or eclipse is divided by a straight line fit to the data either side of the transit or eclipse prior to analysis.
We use all the radial velocities published by GMC+2014 plus the new radial velocities from Table 3 in our analysis. We see no evidence for excess noise in the radial velocities so we use their standard errors as quoted for the calculation of the posterior probability distribution.
In total, we are using nine sets of data, each of which has an uncertain zero-point that should be included in the analysis. Additionally, there are eleven basis functions that are used for the removal of instrumental noise from the CHEOPS data, each with its own coefficient that should be varied independently during the fit to the data. To avoid explicitly calculating these nuisance parameters we use the procedure described by Luger et al. (2017), in which the likelihood for any proposed set of model parameters marginalised over a set of nuisance parameters for a linear model can be calculated by modifying the covariance matrix.
We used (Foreman-Mackey et al. 2013), a implementation of an affine invariant Markov chain Monte Carlo (MCMC) ensemble sampler, to calculate the posterior probability distribution of the model parameters. The maximum-likelihood model fit to the data is shown in Fig. 5. The mean and standard error of the posterior probability distributions for each of the model parameters and various derived parameters are given in Table 5.
The parameters in Table 5 can be combined with our measured value of 2 from the analysis of the SPIRou spectra to determine the masses and radii of both stars with no additional model input. To account for the correlations between parameters, we do this using the sampled posterior probability distribution for the relevant parameters generated by together with a sample of 2 values assuming a Gaussian distribution for this parameter. The masses and radii derived are given in Table 6.

Direct measurement of the stellar effective temperature
The effective temperature for a star with Rosseland radius and total luminosity is defined by the equation = 4 2 SB T 4 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 0, = 0,1 + 0,2 = SB 4 2 1 T 4 eff,1 + 2 2 T 4 eff,2 , 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 J0113+31 provided we can accurately integrate the observed flux distributions for the two stars independently. This is possible because photometry of the combined flux from both stars is available from ultraviolet to mid-infrared wavelengths, and the flux ratio at wavelengths where the majority of the flux is emitted by the primary star has been measured from the TESS and CHEOPS light curves. 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 coloureff relations. The M-dwarf contributes less than 0.2 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 G0-type primary star. The photometry used in this analysis is given in Table 7. The NUV and FUV magnitudes are taken from GALEX data release GR7 (Bianchi et al. 2014) with a correction to the IUE flux scale based on the results from Camarota & Holberg (2014). We assume that the flux from the M-dwarf at ultraviolet wavelengths is negligible. 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). Photometry in the PanSTARRS-1 photometry system is taken from Tonry et al. (2018). Details of the zero-points and response functions used to calculate synthetic photometry from an assumed spectral energy distribution are given in Miller et al. (2020).
To estimate the reddening towards EBLM J0113+31 we use the calibration of E(B−V) versus the equivalent width of the interstellar Na I D 1 line by Munari & Zwitter (1997). To measure EW(Na I D 1 ) we used 11 spectra obtain with the FIES spectrograph on the Nordic Optical Telescope used in medium resolution mode ( = 46, 000). 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 signalto-noise spectrum of the G0V primary star. We then divided each observed spectrum by this spectrum of the G0V primary star after shifting it back to the barycentric rest frame. We then took the median of these residual spectra to obtain a high signal-to-noise spectrum of the interstellar features. The equivalent width of the Na I D 1 line measured by numerical integration is EW(Na I D 1 ) = 77.1 ± 6.0 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 the value of E(B−V) for EBLM J0113+31 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Å. Based on this analysis we obtain the estimate E(B−V) = 0.002 ± 0.012. We use this as a Gaussian prior in our analysis but exclude negative values of E(B−V).
To establish coloureff 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. Coloureff relations suitable for the primary G0V star were calculated in similar way based on stars selected from the Geneva-Copenhagen survey (Holmberg et al. 2009;Casagrande et al. 2011) with 5950 K < T eff < 6250 K, (B − V) < 0.01 and 3.5 < log < 4.5. The results are given in Table 8.
The method we have developed to measure eff for eclipsing binary stars is described fully in Miller et al. (2020). Briefly, we use (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 8, and the Gaussian prior on the reddening described above. The hyperparameters 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 zeropoint 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, i.e. 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. 7 We use linear interpolation to obtain a reference SED for the G0V star appropriate for eff,1 = 6130 K, log 1 = 4.15, [Fe/H] = −0.3 and [ /Fe] = 0.0. For the reference SED for the M dwarf companion we assume eff,1 = 3380 K, log 1 = 5.0, and the same composition. We experimented with distortion functions with 1, 2, 3, 4 coefficients per star and found the results to be very similar in all cases. The results presented here use two distortion coefficient per star because there is no improvement in the quality if the fit if we use a larger number of coefficients. The predicted apparent magnitudes including their uncertainties from errors in the zero-points for each photometric system are compared to the observed apparent magnitudes in Table 7. The posterior probability distribution for the model parameters is summarised in Table 9 and the spectral energy distribution is plotted in Fig. 6.
The random errors quoted in Table 9 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.

Abundance analysis
We have used the H-band spectrum of EBLM J0113+31 A to estimate this star's metallicity. For this abundance analysis we used the observed SPIRou spectra merged into 1-dimensional spectra provided by the observatory. We first subtracted the model spectrum for Table 7. Observed apparent magnitudes for EBLM J0113+31 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 fina column. 1.397 − 0.5812 1 ± 0.045 5.620 − 3.248 2 ± 0.23 Table 9. Results from our analysis to obtain the effective temperatures for both stars in EBLM J0113+31. N.B. eff,1 and eff,2 are subject to additional systematic uncertainty of 10 K and 7 K, respectively. the M-dwarf companion described in Section 3.3 from each of these 1-dimensional spectra, scaled such that the flux ratio in the TESS band matched the value measured from the depth of the secondary eclipse in the light curve. We then co-added the spectra in the rest frame of the primary star to produce a high signal-to-noise spectrum of the G0V star with negligible contamination from the M-dwarf. For the analysis of this spectrum we used iSpec (Blanco-Cuaresma et al. 2014b; Blanco-Cuaresma 2019) with the APOGEE line list for atomic and molecular data in the wavelength range 1500 -1700 nm (Shetrone et al. 2015). We followed Sarmento et al. (2020) in selecting Turbospectrum (Alvarez & Plez 1998;Plez 2012) for the spectrum synthesis assuming a micro-turbulent velocity mic = 1.06 km s −1 with model atmospheres from the MARCS grid (Gustafsson et al. 2008) and solar abundances from Grevesse et al. (2007). We excluded from the fit ±4 nm around the two Brackett series lines at 1681.11 nm and 1641.17 nm, and also some instrumental features that occur near the ends of the échelle orders at 1657 -1659 nm and 1622-1624 nm. We fixed the value of eff = 6124 K and log = 4.15. For the macro-turbulent velocity we used the calibration by Valenti & Fischer (2005) to estimate mac = 4.67 km s −1 . We included the rotational broadening parameter sin as a free parameter in the least-squares fit with a linear limb-darkening coefficient of 0.5 in the H-band based on the results from Claret (2018). We attempted a least-squares fit including the -element abundance as a free parameter but found that the value obtain is not accurate enough to be useful so we fixed [ /Fe] = 0 in the least-squares fit. From this least-squares fit we obtained [M/H] = −0.33 ± 0.01 and sin = 6.6 ± 0.3 km s −1 . There are several additional sources of uncertainty in this analysis, e.g. inaccurate normalisation, errors in atomic data, approximations in the stellar atmosphere models, etc., so the accuracy of our metallicity estimate will be much worse than the precision estimated from the least-squares fitting algorithm (Blanco-Cuaresma 2019; Jofré et al. 2019). Based on the results from independent analyses of APOGEE spectra by Jönsson et al. (2018), we assume an accuracy of 0.15 dex, i.e. [M/H] = −0.33 ± 0.15. The fit to the spectrum is shown in Fig. 7.
We used the co-added FIES spectra of the star to determine the stellar atmospheric parameters ( eff , log , micro-turbulence, and [Fe/H]) and chemical abundances following the methodology described in our previous works (Sousa 2014;Santos et al. 2013;Adibekyan et al. 2012Adibekyan et al. , 2015. In brief, we make use of the equivalent widths (EW) of spectral lines, as measured using the ARES v2 code 8 (Sousa et al. 2015), and we assume ionization and excitation equilibrium. The process makes use of a grid of Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973).
For the stellar spectroscopic parameters we obtained eff = 6025 ± 50 K, log = 4.10 ± 0.05, tur = 1.07 ± 0.06 km s −1 and [Fe/H] = −0.31 ± 0.04. Within the uncertainties, these values are in agreement with those presented in Table|9. In order to be consistent, and because of higher accuracy, we fixed the values of effective temperature and surface gravity to eff = 6124 ± 40 K and log  (Adibekyan et al. 2011).
Using the astrometric data from Gaia EDR3 and the radial velocity of the system (11.179 ± 0.004 km s −1 , GMC+2014) we calculated the Galactic space velocity components ( , , ) = (−17, 16, 21) km s −1 with respect to the local standard of rest (Schönrich et al. 2010). Based on these velocities, adopting the characteristics parameters of Galactic stellar populations of Reddy et al. (2006), and following Adibekyan et al. (2012) we estimated a probability of 99% that the star belongs to the Galactic thin disk, which is in agreement with our conclusion based on the composition of the star.
Based on the results from the analysis of the SPIRou and FIES spectra we adopt the value [M/H] = −0.3 ± 0.1 for the metallicity of EBLM J0113+31. The co-added SPIRou spectra corrected for the contribution from the M-dwarf and the co-added FIES spectrum are available from the supplementary online information that accompanies this article.

Astrometric noise due to binary orbital motion
The projected semi-major axis of the G0V star's orbit is 1 = 1 / = 0.11 mas, so we expect excess noise in the Gaia astrometry ≈ 0.1 mas due to the orbital motion of the primary star. Indeed, the astrometric excess noise in the Gaia EDR3 catalogue for EBLM J0113+31 is 0.163 mas. This is higher than expected for a good fit to the data for a single star with G ≈ 10, and consistent with the noise expected from the orbital motion of the G0V star. This will only lead to a systematic error in the parallax if the position angle of the binary at the times of observation are not randomly distributed around the binary star orbit. This can be checked using the parameter ipd_gof_harmonic_amplitude provided in the EDR3 catalogue (Lindegren et al. 2021). For EBLM J0113+31, this parameter takes the value 0.014, which is less than the median value of this statistic for stars with 6-parameter solutions in the magnitude range G=9-12 (0.020). Although the detection of the astrometric noise is statistically significant, it is a small contribution to the uncertainties on the parallax. The renormalised unit weight error for EBLM J0113+31 is RUWE=1.154, which is only slightly higher than the median value for stars with 6-parameter solutions in the magnitude range , and is close to the expected value ≈ 1 for "for well behaved sources".
We can therefore be confident that the orbital motion of the G0V star does not produce a systematic error in the measured Gaia parallax.

Comparison to stellar evolution models
The mass, radius and effective temperature for both stars in EBLM J0113+31 are given in Table 6, together with the derived surface gravity, mean stellar density and luminosity for both stars.
We used the software package (Maxted et al. 2015) to compare the parameters of the primary star, EBLM J0113+31 A, to a grid of stellar models computed with the stellar evolution code (Weiss & Schlattl 2008). The methods used to calculate the stellar model grid are described in Serenelli et al. (2013). uses 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]. We find a very good fit to the observed parameters of EBLM J0113+31 A for an age of 6.7 ± 0.5 Gyr, as can be seen in Fig. 8. More the 99 per cent of samples from the PPD correspond to models where EBLM J0113+311 A is a post main-sequence star that has exhausted all the hydrogen in its core. The model grid accounts for diffusion so the initial metal abundance for this star is inferred to be [Fe/H] = −0.2 ± 0.1. Isochrones for the same age and initial metal abundance from the Dartmouth stellar evolution database (Dotter et al. 2008) and the MESA Isochrones & Stellar Tracks (MIST, Choi et al. 2016) are also shown in Fig. 8. There is very good agreement between these different stellar evolution codes, as might be expected given that the properties of EBLM J0113+31 A are similar to the Sun and all three grids of stellar models are calibrated to match the observed properties of the Sun.
The same isochrones from the Dartmouth and MIST stellar model grids are compared to the properties of EBLM J0113+31 B in Fig. 9. Our grid of models does not extend to these very low masses. The agreement between the models and observations is reasonably good, which is somewhat surprising given the long-standing observation that stellar models tend to under-predict the radius and over-predict eff for low-mass stars (Spada et al. 2013;Cassisi & Salaris 2019;Zhou et al. 2014;Berger et al. 2006;Hoxie 1973;Lacy 1977). This can be seen from the mass, radius and eff measurements for six other very low-mass stars in the same figure. These six stars are members of three eclipsing binaries with orbital periods less 2 days. This complicates the interpretation of their properties in the light of the so-called "radius inflation" problem because these stars will be forced to rotate much faster than most single M-dwarf stars by tidal forces in these short-period binaries. EBLM J0113+31 B is a valuable addition to the small sample of well-characterised VLMSs because we have an independent estimate of its age and initial metal abundance based on observations of the G0V primary star to add to the accurate mass, radius and eff measurements.

EBLM systems as benchmark stars
Benchmark FGK dwarf stars with direct eff measurements based on angular diameters measured by interferometry typically have apparent magnitudes V = 1 -6 (Jofré et al. 2014). This is 5 -10 magnitudes brighter than the magnitude limits for large-scale spectroscopic surveys, so special observing modes must be employed to obtain spectra of these benchmark stars. These bright benchmark stars also tend to be single stars, so there are often no direct measurements of their mass or surface gravity. It is difficult to extend this sample because new candidates for benchmark stars will necessarily be more distant than the existing sample, i.e. they will have smaller angular diameters than the existing benchmark stars. For example, a nominal Sun-like star at distance of 10 pc will have an angular diameter = 0.465 mas, so a systematic error of only 0.04 mas, which is typical for existing Figure 7. The H-band spectrum of EBLM J0113+31 A (blue) and a synthetic spectrum fit by least-squares using iSpec (red). Residuals from the synthetic spectrum fit are shown in green offset vertically by 1.05 units. measurements (Karovicova et al. 2022), implies a systematic error of 250 K in the measured value of eff for such a star.
In contrast, EBLM J0113+31 is within the magnitude range of recent large-scale spectroscopic surveys, e.g. the TESS-HERMES survey (10 < V < 13.1, Sharma et al. 2018), LAMOST "VB mode" observations (9.0 ≤ ≤ 12. 5, Luo et al. 2015), and stars in open clusters observed as part of the Gaia-ESO survey (9 < V < 16. 5, Bragaglia et al. 2021). This makes it feasible to observe EBLM J0113+31 and other EBLM binaries in exactly the same way as other stars observed by these survey instruments as part of their routine operations. The contribution of the M-dwarf to the total flux at optical wavelengths is < ≈ 0.2 per cent for EBLM binaries, so the M-dwarf will have a completely negligible effect on the atmospheric parameters derived from the analysis of the optical spectrum. This makes it possible to make an "end-to-end" test of the accuracy of parameters derived by the combination of these survey instruments plus their data processing and analysis pipelines. Even at near-infrared wavelengths used by surveys such as APOGEE (Jönsson et al. 2018) the contribution from the M-dwarf is < ≈ 1 per cent, so the results of any analysis that includes a correction for this small contribution to the total flux will be insensitive to the details of how this correction is done.
Many EBLM binaries in the magnitude range 10 < ≈ V < ≈ 12 have been identified and have well-determined spectroscopic orbits that have been published (Triaud et al. 2017) or that are in preparation thanks to the EBLM project and BEBOP survey (Standing et al. 2022). High-quality space-based photometry is already available for many of these stars from the TESS survey and/or from our on-going CHEOPS GTO programme. Several échelle spectrographs that can provide high-resolution spectroscopy at near-infrared wavelengths are currently operational on 4 -10 m telescopes, e.g. CARMENES on the Calar Alto Observatory 3.5-m telescope (Quirrenbach et al. 2016), NIRPS on the ESO 3.6-m telescope (Grieves et al. 2021), CRIRES+ on the ESO 8.2-m VLT (Kaeufl et al. 2004), and IRD on the 8.2-m Subaru telescope (Kotani et al. 2014). We can also look forward to high-quality spectrophotometry and improved parallax measurements for these EBLM systems in future data releases from the Gaia mission. 9 In summary, the instrumentation, data and targets needed to create a network of moderately-bright FGK dwarf stars covering both hemispheres that are ideal benchmark stars for ongoing large-scale spectroscopic surveys are all now available.
Apart from their utility as benchmarks stars for large-scale spectroscopic surveys, follow-up observations of additional EBLM systems will also provide valuable data on the properties of very low-mass stars. With observations similar to those presented here we can create a sample VLMSs with precise and accurate eff , mass and radius measurements. These EBLM binaries will have independent estimates for their age and initial metallicity based on the observed properties of the primary stars in these systems. It is not feasible to obtain a direct spectrum for these very faint companion stars, but it should be possible given sufficiently high-quality data to estimate the projected rotational velocity of the star from the width of the peak in the stacked-CCF. Data of this quality will be very useful for testing and calibrating models of very low mass stars that include additional physics to account for the radius inflation problem (Mullan et al. 2018;Feiden & Chaboyer 2014).
Many of these EBLM binary systems will also be ideal benchmark stars for the upcoming PLATO mission (Rauer et al. 2014) if we can measure model-independent masses for the primary star using the techniques presented in this study. The PLATO mission will focus on bright stars (4 -11 mag) with the aim to detect and characterize planets down to Earth-size by photometric transits. Asteroseismology will be performed for these bright stars to obtain stellar parameters, Figure 8. EBLM J0113+31 A in the mass-radius and Hertzsprung-Russell diagrams compared to isochrones for an age of 6.7 ± 0.7 Gyr assuming an initial metal abundance [Fe/H] = −0.2 interpolated from a grid of stellar models. The ellipses show 1-and 2-confidence regions on the parameters of EBLM J0113+31 A. Also shown are isochrones for the same age and initial metal abundance from the Dartmouth stellar evolution database (cyan dotted line) and MIST (green dashed line).
including masses and ages. The PLATO Definition Study Report 10 ("red book") specifies that PLATO must be capable of delivering accurate stellar ages with a precision of 10 per cent. Some corrections for systematic errors in the current generation of stellar models will be needed to reach this accuracy in stellar ages (Goupil 2017). The planned observing strategy includes a step-and-stare phase that will cover about 50 per-cent of the sky. EBLM binaries can be used to perform "end-to-end" tests of the PLATO data analysis to ensure that the mass estimates delivered for these stars are accurate, and to calibrate the next generation of stellar models using direct mass, radius, and eff measurements combined with asteroseismology.

CONCLUSIONS
We have derived precise and accurate masses, radii and effective temperatures for both stars in the eclipsing binary system EBLM J0113+31. These data can be used to validate and calibrate stellar models, empirical relations for stellar properties, and to test data analysis techniques. With the techniques established here, it is feasible to create a network of moderately-bright FGK dwarf stars covering both hemispheres that are ideal benchmarks for on-going large-scale spectroscopic surveys and for the upcoming PLATO mission.

ACKNOWLEDGEMENTS
CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. The CHEOPS Consortium would like to gratefully acknowledge the support received by all the agencies, offices, universities, and industries involved. Their flexibility and willingness to explore new approaches were essential to the success of this mission.
Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) which is operated from the summit of Maunakea by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. The observations at the Canada-France-Hawaii Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site. Based on observations obtained with SPIRou, an international project led by Institut de Recherche en Astrophysique et Planétologie, Toulouse, France. PM is grateful to the observatory staff at the CHFT for their help with the planning, execution and data reduction of the SPIRou data used for this study.
Based on observations made with the Nordic Optical Telescope, owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku and the University of Oslo, representing Denmark, Finland and Norway, the University of Iceland and Stockholm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias.