BEBOP V. Homogeneous Stellar Analysis of Potential Circumbinary Planet Hosts

Planets orbiting binary systems are relatively unexplored compared to those around single stars. Detections of circumbinary planets and planetary systems offer a first detailed view into our understanding of circumbinary planet formation and dynamical evolution. The BEBOP (Binaries Escorted by Orbiting Planets) radial velocity survey plays a special role in this adventure as it focuses on eclipsing single-lined binaries with an FGK dwarf primary and M dwarf secondary allowing for the highest-radial velocity precision using the HARPS and SOPHIE spectrographs. We obtained 4512 high-resolution spectra for the 179 targets in the BEBOP survey which we used to derive the stellar atmospheric parameters using both equivalent widths and spectral synthesis. We furthermore derive stellar masses, radii, and ages for all targets. With this work, we present the first homogeneous catalogue of precise stellar parameters for these eclipsing single-lined binaries.


INTRODUCTION
Circumbinary exoplanets orbit around both stars of a close binary system.The first confirmed systems have been discovered in the ★ E-mail: AXF859@student.bham.ac.uk last decade using space-based transit observations e.g.Kepler-16 b (Doyle et al. 2011) and TOI-1338 b (Kostov et al. 2020).Radial velocity (RV) detections of circumbinary planets have proven difficult because of stellar contamination, which has limited the RV precision that can be obtained in double-lined binaries (e.g.Konacki et al. 2010).Nevertheless, Triaud et al. (2022) confirmed the radial-velocity signal of the circumbinary planet Kepler-16 b.This was possible because the host binary is composed of a solar-type mainsequence star and a low-mass M-dwarf companion.This combination allowed us to treat the primary star spectroscopically as a single-lined star, for which achieving the necessary m s −1 precision to detect planetary signals is routinely done (e.g.Faria et al. 2022).
The BEBOP survey for circumbinary planets (Binaries Escorted by Orbiting Planets) focuses on the EBLM project (Eclipsing Binaries with Low-Mass companions).Root-mean-square (RMS) scatter in residuals after removing the RV signal caused by the binary can reach down to about 3 m s −1 (Standing et al. 2022).The survey has been designed as a blind all-sky survey and its sample was constructed from EBLM binaries detected through transit surveys (e.g.Triaud et al. 2013Triaud et al. , 2017;;von Boetticher et al. 2019;Lendl et al. 2020).Recently, BEBOP made its first discovery of a circumbinary planet solely based on RV measurements, EBLM J0608-59/TOI-1338/BEBOP-1 c (Standing et al. 2023, BEBOP IV).
Planet parameters, as measured through the transit or RV method, are inherently relative to their host star properties and such properties are often inferred from stellar evolution models (e.g.Dotter et al. 2008).Systematics introduced from different models are thus carried forward to analyses of the other planetary bodies in the system.It is therefore important to use a homogeneous and reliable set of stellar parameters for our BEBOP target stars to minimise any systematics in stellar as well as planetary parameters and ensure the statistical validity of the full survey.Several similar efforts have been done for other surveys (e.g.Sousa et al. 2011;Buchhave et al. 2014).
Throughout the BEBOP survey, we have already assembled a large archive containing thousands of high-resolution spectra of the primary FGK main-sequence stars.These spectra are first and foremost used to build up our RV time series, but they can also be used to measure stellar atmospheric parameters.Furthermore, thanks to the relative brightness of this sample (V ∼ 8 − 13 mag) precise stellar parallaxes from Gaia DR3 (Gaia Collaboration et al. 2016aCollaboration et al. , 2023) ) as well as broadband photometry from the 2MASS (Skrutskie et al. 2006) and AllWISE (Cutri et al. 2021) surveys are available for all targets.This allows us to use these in combination with our spectroscopic parameters to derive homogeneous masses, radii and ages for all targets.
In this work, we homogeneously measure stellar parameters for BEBOP's primary stars, using the same methodology overall with consistent input data.We analyse high-resolution spectra to derive the effective temperature (T eff ), surface gravity (log  ★ ), metallicity ([Fe/H]), and projected rotational velocity ( sin  ★ ) of each star.We then use these parameters to derive stellar physical parameters such as mass ( ★ ), radius ( ★ ), and stellar age.These homogeneous parameters are of interest to produce accurate secondary star masses and radii as part of the EBLM survey (Triaud et al. 2017), and accurate circumbinary masses in the context of the BEBOP survey (Martin et al. 2019).The paper is structured as followed: In Section 2 we introduce the BEBOP sample and in Section 3 the spectroscopic data.Our analysis method is described in Section 4 with results in Section 5. Finally, we conclude in Section 6.

THE BEBOP SAMPLE
A total of 179 systems are analysed in this paper with a magnitude range m v = 8.31 to m v = 12.96 (see Figure1).The sample is split into a Northern sample (93 systems observed with the SOPHIE spectrograph, see Section 3.1), and a Southern sample (110 systems observed with HARPS, see Section 3.2).A total of nine systems are common between the Northern and Southern samples, selected on purpose in order to compare the sensitivity of both instruments to circumbinary exoplanets.This particular sub-sample is used to cross calibrate the spectroscopic parameters produced by both instruments.

BEBOP-South Sample
The BEBOP-South sample was defined first.All systems identified as part of the EBLM sample (see Triaud et al. 2013Triaud et al. , 2017) ) were considered.Those EBLM secondaries were identified as transiting exoplanet false-positives during the WASP survey (Wide Angle Search for Planets; Pollacco et al. 2006;Triaud 2011) thanks to the CORALIE spectrograph (Udry et al. 2000).All had received at least 13 spectra and some many more (e.g.Martin et al. 2019).The BEBOP sample was selected to maximise radial-velocity precision, and minimise the contribution of the secondary stars.All visually identified double-lined binaries within the CORALIE spectra were removed.A Keplerian model was adjusted to all systems to obtain preliminary masses for the secondaries.All systems with binary period  bin > 4.1 days (because no circumbinary has been found where  bin < 5 days Martin 2018), where the variance in the span of the bisector slope of individual spectra (as defined in Queloz et al. 2001) is below 177 m s −1 , and where the full width at half maximum of the absorption lines is below 28 km s −1 were kept (both to maximise radial-velocity precision and improve sensitivity to circumbinary exoplanets).This resulted in a sample of 56 eclipsing binaries, all expected to have a Δmag > 4 between primary and secondary stars.
In addition, we selected 22 systems identified as likely EBLMs by the KELT survey (Kilodegree Extremely Little Telescope; Pepper et al. 2012;Collins et al. 2018) with  < +10 • ,  bin > 5 days, with spectral types later than F4, and  < 10.5 for F-types,  < 11 for Gtypes and all K-types.Typically orbital parameters for these systems were more poorly determined than for the EBLM sample.
Finally, since TESS was launched (Transiting Exoplanets Survey Satellite; Ricker et al. 2014), 32 new eclipsing systems consistent with EBLMs were added in 2021.These are typically brighter than the original EBLM sample, but usually without any prior radialvelocity information except in rare cases such as TOI-222 (Lendl et al. 2020).Their other properties, such as  bin , are consistent with the rest of the sample.

BEBOP-North Sample
The EBLM project's northern counterpart used SuperWASP to find candidates, and then SOPHIE to identify EBLM false positives (e.g.Gómez Maqueo Chew et al. 2014).However, observations and classifications were not as systematic as in the South.As such, the BEBOP-North sample was selected first from the KELT catalog (Collins et al. 2018), cross-matched with SuperWASP/SOPHIE observations for confirmation.In addition all SuperWASP/SOPHIE false positives were reviewed selecting likely EBLMs from the eclipse depth and the absence of visible secondary eclipses.As in the South, only binaries with  bin > 5 days, with spectral types later than F4 were kept.Only objects in the Northern hemisphere ( > 0 • ) were selected, resulting in a sample of 120 systems.A first reconnaissance campaign was conducted on SOPHIE in 2018 to remove double-lined binaries and confirm the binary nature of each system, reducing the sample to 93 binaries.First all systems with  > 11.5 were observed with the High Efficiency mode and all others in the High Resolution mode.After a few observations were collected, all systems with line widths > 15 km s −1 were moved to High Efficiency mode since for those, spectral resolution is not as much of an issue.
In both the Northern and Southern samples, systems were divided into a primary and a secondary sample.Typically, systems in the primary sample have more precisely measured RVs with the goal to collect of order 40 to 50 spectra and detect circumbinary planets.Systems in the secondary sample typically receive of order 10-15 measurements only.

OBSERVATIONS
The data set analysed in this paper consists of 4512 high-resolution spectra obtained with the SOPHIE and HARPS spectrographs from 2013 to 2023, the majority of which (70%) were observed with a 1800s exposure time.

SOPHIE spectroscopy
The SOPHIE échelle spectrograph (Perruchot et al. 2008;Bouchy et al. 2009) is mounted on the 193 cm reflector telescope at the Haute-Provence Observatory.SOPHIE has a wavelength range from 387.2 nm to 694.3 nm.Two observation modes are available: highresolution (HR) and high-efficiency (HE), respectively having resolutions of  = 75 000 and  = 40 000.Using the HE mode allows for a throughput increase equivalent to one magnitude.The versatility of the two modes is well demonstrated in the data set.Out of the 93 targets in the BEBOP Northern sample (taken using the SOPHIE spectrograph), 54 were observed in HE mode and 56 in HR mode (meaning 17 were observed in both modes).The median signal-tonoise ratio (SNR) achieved for individual spectra in the HE mode is 27, and the median SNR for the HR mode is 33.The sample is presented in Figure 1.When considering the combined spectra used for analysis, the median SNR for the HR mode is 88, and 87 for the HE mode.
All spectra are acquired as an échelle onto a CCD camera.The instrument is calibrated with Tungsten lamps at the start of every night to locate where each spectral order is, as well as to perform a flat field.Biases and darks are also obtained daily.In addition, Thorium-Argon lamps and a Fabry-Pérot etalon are used to establish an accurate wavelength solution.A number of reference stars are observed every night in both the HR and HE mode to track the stability of the instrument (typically of order 2 m s −1 ; Bouchy et al. 2013;Courcol et al. 2015;Hara et al. 2020).Additional Fabry-Pérot calibrations are obtained roughly every two hours throughout the night.
SOPHIE operates with two fibres.All observations obtained for BEBOP use the objAB mode where one fibre is on target and the other is on the sky so as to remove any contamination from e.g.Moonlight.Only one system was observed with the wavesimult mode where the second fibre is instead illuminated by a Fabry-Pérot etalon to produce a simultaneous calibration.This mode is usually reserved for systems where the most extreme radial-velocity precision is needed, which is not the case for the BEBOP stars.In the case of EBLM J0626+29 however, the sky fibre, which cannot be moved, coincided with another star.
Using the calibration frames, each spectral order is extracted from the CCD by the SOPHIE Data Reduction Software (DRS), producing an e2ds file.Each order is corrected for the instrumental blaze function and stitched together to create a one-dimensional spectrum, Figure 1.The median SNR of individual spectra is plotted against the apparent V band magnitude for each BEBOP sample target.The North Sample data is split into the two observation modes available with SOPHIE.The HR mode is shown as purple circles, and the HE mode as orange triangles.The South Sample, coming from HARPS, is shown as pink crosses.Note here that HARPS spectra have, on average, a shorter exposure time -this plot should not be taken as a representation of the instrumental performance.the s1d files, which is what we use for our analysis.Each individual s1d's wavelengths solution is corrected to the barycentre of the Solar system (Bouchy et al. 2009;Courcol et al. 2015).

HARPS spectroscopy
The HARPS spectrograph is mounted on the ESO 3.6 m telescope at La Silla observatory in Chile (Mayor et al. 2003).HARPS has a resolving power  = 115 000 over a wavelength range 378 nm − 691 nm.Individual HARPS spectra used within this paper have a median SNR of 30.Observations and recording of spectra are done in a similar fashion to SOPHIE.Like SOPHIE, HARPS utilises two fibres, with one on the target and the other used as a calibrationeither by use of a Th-Ar reference spectrum, or the sky background.
The data reduction for HARPS follows very closely the method outlined for SOPHIE.The main difference is that we only used the HR mode for observation with HARPS (called HAM).Like for SOPHIE, HARPS can either be used in objAB or in wavesimult mode.All BEBOP observations used the objAB mode except for four systems identified by the TESS mission, because their brightness and spectral properties allowed us to reach a photon noise below HARPS's long term stability (around 1 m s −1 ; Fischer et al. 2016) and a simultaneous calibration was necessary to make the best use of the instrument.HARPS is stable enough that neither reference stars nor intra-night calibrations are needed.
The radial velocity of the target spectrum can then be extracted by both the SOPHIE and HARPS reduction pipelines as the mean velocity from a Gaussian fitting to the cross-correlation function (CCF) profile.Out of all BEBOP targets, 6 were fitted with a K type mask by the HARPS and SOPHIE data reduction pipelines, with all others being fitted with G type masks.
Figure 1 shows the HARPS spectra follow the same general trend as the SOPHIE data in terms of the SNR achieved for individual spectra.After coadding, the median spectral SNR is 126, surpassing that achieved by SOPHIE as outlined in Section 3.1.All targets observed with HARPS were also observed using CORALIE, however the former was used due to the superior SNR achieved.

SPECTROSCOPIC ANALYSIS
Two commonly used methods to retrieve atmospheric parameters via spectral analysis include the curve-of-growth equivalent widths (EW) (e.g.Sousa et al. 2011;Sousa 2014) and spectral synthesis (e.g.Valenti & Piskunov 1996;Adibekyan et al. 2012;Tsantaki et al. 2018) methods.Alternative methods are also used extensively, such as the neural-network based Payne algorithm as used by the SAPP pipeline for the analysis of the PLATO core sample (Gent et al. 2022).The EW method uses the neutral and ionised absorption lines of only one element, resulting in a quick return of parameters.The atomic species used depends on the conditions of the star.For example, for young and/or active stars, it can be beneficial to use titanium (Ti) lines, as described by Baratella et al. (2020).However, for this work we employ the more commonly-used Fe I and Fe II lines, which are abundant in the spectra of main-sequence FGK stars.
On the other hand, spectral synthesis is more computationally intensive, iterating through parameters that synthesise a spectrum until the synthesised one matches the observed spectrum.Analyses of large data-sets would benefit from the speed of the EW method, however this method is unable to constrain the projected rotational velocity or macroturbulent velocity since a line's EW is conserved under these broadening parameters (e.g.Sousa et al. 2011;Santos et al. 2013).To provide the most complete stellar information, the synthesis method must also, therefore, be implemented.
For the analysis of the BEBOP survey stars, we combine the strengths of both methods to provide a full and homogeneous analysis of the spectra of these FGK main-sequence stars.We utilise the speed of the EW method to provide excellent initial parameters for the synthesis method, which then runs much quicker.We use the iSpec framework (Blanco-Cuaresma et al. 2014;Blanco-Cuaresma 2019) for our analysis.In this section we describe the details of our analysis and the tests done to ensure reliability of our parameters.In Appendix A, we describe the public python pipeline, PAWS1 , we wrote around the iSpec framework to perform our analysis.

Data Preparation
Barycentric velocity correction, made necessary due the the Earth's motion, is performed by the SOPHIE and HARPS reduction pipelines (Bouchy et al. 2009;Lovis & Pepe 2007), resulting in only needing to handle RV corrections in order to shift the BEBOP spectra into the lab frame.To perform the radial velocity correction, an atomic line mask developed from the NARVAL solar spectrum (Aurière 2003) is used as a comparison template, representing the lab frame.A cross-correlation function (CCF) is then used to determine the radial velocity of the target star, which is then corrected for in the individual spectra to shift them into the lab frame.
Consistent continuum normalisation is a crucial step to ensuring homogeneity throughout the analysis.Continuum flux in every spectrum is allocated as 1.0, with all spectral features, and therefore analysis, relative to this.To perform the normalisation, we use the fit_continuum function implemented in iSpec.Continuum fluxes are found using a median filter of window size 0.05 nm, and maximum filter of 1.0 nm.Noise is identified by the median filter, then the maximum filter is used to block the fluxes in absorption lines.A model of the continuum is created by fitting a B-Spline of 2 degrees to the spectra, every 5 nm.The spectrum is then normalized by dividing all fluxes by this model.
Individual spectra from the same target are co-added prior to analysis to form one higher-SNR spectrum.Average flux and flux error at wavelength steps of 0.001 nm are taken with a default range of 420 − 680 nm.Testing yielded no significant differences in results for varying wavelength steps.All BEBOP spectra were treated with the same wavelength range to ensure homogeneity.Spectra from the BEBOP sample are supplied without flux errors; we estimated these by dividing the flux at each pixel by spectral SNR provided in the FITS headers.

Line List and Model Atmosphere
A line list must be input during spectral analysis -this provides a subset of lines in the spectrum that will be used to determine the atmospheric parameters.To ensure homogeneity throughout the analysis, the same line list was employed for both the EWs and synthesis parts of our method.The line list created for the SPECTRUM code (Gray & Corbally 1994), built on the NIST Atomic Spectra Database (Ralchenko 2005), was chosen due to its proven success and versatility in FGK dwarf analysis (Blanco-Cuaresma et al. 2014).
Absorption lines selected for analysis in the spectra are identified using a line mask.The line masks provided with iSpec contain atomic information inherently dependent on the spectral resolution of  = 47 000 for which they have been optimised.For use for the SOPHIE HE mode spectra, minimal modification was required for this line list, with only 27 lines removed due to consistently performing poorly in chi-squared testing -the wavelengths of these lines can be found in Table B1.HE mode spectra are more at risk of suffering from 'blended' lines, in which the shallower and broader lines of the lower resolution spectrum may blend together.
To preserve the higher resolution achieved in the SOPHIE HR and HARPS spectra, we developed new line masks for this work using iSpec, which are publicly available on GitHub.For use with the HARPS spectra, we identified line masks using the HARPS-N solar spectra (Dumusque et al. 2021), and the NARVAL solar spectrum at a resolution of  = 65 000 (Aurière 2003) for the SOPHIE HR spectra.Individual abundances were then calculated for each line in these masks -where these were not within 0.05 of the accepted solar values, the masks were discarded to avoid lines that would not be suitable for atmospheric parameter determination.In the case of spectral synthesis, this constraint was extended to within 0.1 dex of solar abundances to allow for more lines to be available for synthesis.
Both analysis methods require a model atmosphere grid to be input.We chose to use the ATLAS9 set of model atmospheres (Kurucz 2005).ATLAS is inclusive of a range of 4500 to 8750 K in T eff , 0.00 to 5.00 dex in log  ★ , and -5.00 to 1.00 dex in metallicity.Computation time is saved by the models only working with plane-parallel geometry atmospheres, which assume local thermodynamic equilibrium (LTE) and neglect 3D convection.This assumption breaks down for cold giant and super-giant atmospheres -although on the whole valid for FGK stars, it should be cautioned that biases may be introduced in Fe abundance and should be considered when abundances accurate to a few percent are desired, as described by Bergemann et al. (2012).
Studies such as those by Cooke et al. (2020) often demonstrate extreme difficulty in constraining log  ★ from SOPHIE HE mode spectra, with that particular study reporting their log  ★ value with a significant error of 0.22 dex.Not only is such uncertainty seen in results from lower-resolution spectra, studies such as those by Torres et al. (2012) also reveal log  ★ determination to be highly problematic in spectra of  = 46000, 48000, 51000, and 68000.Careful consideration was paid throughout this analysis to the fact that the spectroscopic log  ★ values derived for the BEBOP targets may be less reliable than those obtained via other methods such as photometry or via independently derived stellar mass and radius.

Analysis Process
Using either the EW or synthesis methods can produce reliable results.In this section, we discuss a robust approach that allows us to adopt both methods in an homogeneous way.Atmospheric parameters are first derived via the EW method, however due to this method's inability to constrain  sin  ★ and  mac this does not obtain the full set of desired parameters.These results are thus used as initial parameters for the synthesis method, which derives the final set of atmospheric parameters.

Equivalent Widths
We first use the EW method, employing only the Fe I and Fe II lines from the line list.We fit a Gaussian profile to each spectral line separately using ARES (Sousa et al. 2015).The equivalent width is defined by taking a rectangle with a height equal to that of the continuum, and varying the width of it until the rectangular area is equal to that of the line under the continuum.The WIDTH radiative transfer code (Sbordone et al. 2004) then derives individual abundances for these Fe lines based on a set of initial atmospheric parameters -in this work, these were set to the solar values collated in Blanco-Cuaresma (2019) since the BEBOP targets were chosen to be FGK dwarfs.
Through a minimisation procedure, stellar parameters are then varied with the best fitting parameters ensuring ionisation and excitation balance.We fit for T eff , log  ★ , and [Fe/H]. mic is estimated using the estimate_vmic function included in iSpec; this relation depends on T eff , log  ★ , and [Fe/H], and was derived using the results of GES (Gaia-ESO Survey) UVES data release 1 (Jofré et al. 2014;Blanco-Cuaresma et al. 2014).

Spectral Synthesis
Unlike the EW method, spectral synthesis uses every line in our line list.The WIDTH radiative transfer code is not able to perform spectral synthesis, so this part of the pipeline calls upon the SPECTRUM code (Gray & Corbally 1994), chosen due to proven speed and reliability in the analysis of FGK dwarf stars by Blanco-Cuaresma (2019).Stellar atmospheric parameters are used with SPECTRUM to generate a synthesised spectrum.Parameters are iterated and used in conjunction with a minimisation algorithm to determine the optimal fit of the synthesized spectrum to the observed one.T eff , log  ★ , [Fe/H] and  sin  ★ are fitted for, whereas  mic and  mac are calculated using estimate_vmic and estimate_vmac, with the relation for  mac again based upon the GES UVES results as described for  mic in Section 4.3.1.
The synthesis process is sensitive to its initial conditions.We therefore use the EW method to set reasonable estimates for the first iteration.This ensures that the synthesis begins in a parameter space that reflects the observed spectrum, saving considerable time compared to if beginning from a solar input for all targets.We set the maximum number of iterations to 6; Blanco-Cuaresma et al. (2014) details how this is the optimal number as more iterations can cause metallicity dispersion to be favoured disproportionately compared to other parameters.The errors on derived parameters are calculated using the covariance matrix generated by the least-squares fitting (Blanco-Cuaresma et al. 2014).In the case of T eff , the precision errors reported by iSpec are smaller than the expected accuracy of model atmospheres.To compensate for this, the errors are inflated by adding 100 K in quadrature to better reflect the uncertainty (Tayar et al. 2022).

Handling Low SNR Data
Since the SNR of spectra is known to be a critical aspect in stellar analysis, we analysed the individual results of using varying SNR spectra for the same target.Such testing was aimed primarily at producing an SNR filter value, under which spectra are not used.Sousa et al. (2008a) state that 90% of their spectra surpass an SNR of 200 for their equivalent widths analysis; suggesting this filter would mainly be in place for the first step of our method.We made use of the target EBLM J0002+47, with the primary star being a slow-rotating F-type dwarf.The target has 33 spectra available with a large range of SNR, observed using the SOPHIE HR mode.Individual spectra are not combined as part of this testing; instead the EW and synthesis methods are applied separately using solar input parameters to each spectrum.
The results of using only the EW method are shown in Figure 2a, shown as the derived T eff for each individual spectrum plotted against its SNR.Included within the plot is a comparison to literature T eff values for the target, with the TESS Input Catalogue (TIC) value (Stassun et al. 2019) represented in green with the appropriate error range of 140 K, and the Gaia DR3 value represented in purple (Gaia Collaboration et al. 2016bCollaboration et al. , 2022;;Babusiaux et al. 2022).Although lower than the TIC value, the Gaia DR3 T eff (shown as the purple horizontal line) falls within the TIC T eff errors.A very clear trend is displayed, showing a clear deviation from literature T eff decreasing as spectral SNR increases, as one would expect.An SNR of 20 marks a significant cut-off, above which the majority of results lie within the accepted range of the TIC T eff value.To ensure deviations such as those shown in this plot are kept to a minimum, individual spectra with SNR < 20 were filtered out prior to analysis.The rootmean-square deviation (RMSD) from the TIC value using all results EW-only testing is 208 K, whereas after removing spectra with SNR < 20 decreases the RMSD to 91 K -a significant improvement.An additional concern is displayed in the lowest SNR spectrum of Figure 2a having the smallest error in T eff ; this suggests that the error derived for low SNR spectra by only the EW method does not accurately reflect the uncertainty of the value.Indeed, these uncertainties are statistical only and do not reflect all systematic effects involved in the parameter determination as they are not intended to be the final product of the method.Additionally, the uncertainties determined from only the EW method are highly sensitive to the flux errors supplied (Blanco-Cuaresma et al. 2014) -if these were estimated incorrectly, the uncertainties would be affected.
Figure 2b shows the result of using only the synthesis method on the individual spectra of J0002+47.No obvious trend is observed in the T eff values derived with respect to the spectral SNR.Instead, a systematic under-estimation of between 50 and 200 K below the TIC T eff is displayed from the majority of spectra.This is likely to be influenced by the solar input parameters used by default with the synthesis method (i.e. when no input parameters are specified), with all results scattered within a much smaller region than in Figure 2a.It is interesting that the synthesis method produced an effective temperature within error-bars of both literature values using the lowest SNR spectra available, however the lower SNR is reflected by a poorer synthesis fit and thus larger uncertainties.Again investigating the RMSD of the results, the synthesis-only testing has a RMSD of 155 K from the TIC value, greatly reduced compared to the 208 K from EW-only testing.Figure 2b visibly shows much less scatter in the results than Figure 2a; if we consider the RMSD from the mean of the results, rather than from the TIC value, a value of 83 K is achieved.This demonstrates the ability of the synthesis method to derive consistent parameters regardless of spectral SNR.
T eff derived by spectral synthesis appears independent of spectral SNR, however indicates an underestimation bias.
Figure 2c demonstrates the benefits of combining both methods.Where the synthesis method displays a bias to its input parameters, and the EW method has particularly poor performance at low SNR, combining the two gives a far lower dispersion of results.Returning to the metric of RMSD from the TIC T eff , using the entire pipeline returns an RMSD of 123 K, the lowest of the three tests.Additionally, the RMSD from the mean T eff determined from the pipeline is 92 K, again showing strong consistency.In fact, all but one of the results agrees with the TIC values.Agreement is stronger with the Gaia DR3 effective temperature in Figure 2c, as is also reflected in the higher SNR region of Figure 2a.No metric is available to determine which value represents the most reliable result, so within this paper they are treated as equally likely.Given that the results in Figure 2c agree largely with both values, there is no cause for concern or necessity to prove one value as more physically correct than the other.
In cases where the SNR of the combined spectrum is extremely low (< ≈ 50), we chose to skip the EW step and purely use the synthesis method for analysis.As shown in Figure 2, the synthesis method is far less affected by low SNR spectra, whereas the EW method can produce results that differ hugely from expectations.We chose to only compare the effective temperatures as photometric metallicity may not be a reliable reference (Morrell & Naylor 2019).

Handling Fast Rotators
It is well-established that for fast-rotating stars ( sin  ★ ≳ 5 km s −1 ) the reliability of the EW method is reduced due to the blending of spectral lines (Tsantaki et al. 2014).To deal with this, we adopt the same approach as in Section 4.4, in which the EW method is skipped.The full list of targets for which the EW method was skipped is shown in Table C1.The blending of lines by high rotational velocity manifests as an increased FWHM of the CCF described in Section 4.1, which is essentially the average shape of the spectral lines.The FWHM of the CCF can therefore be used to determine an estimate for the  sin  ★ , as detailed by Rainer, M. et al. (2023).Where the FWHM indicates that the star has  sin  ★ ≳ 5 km s −1 (i.e.FWHM ≳ 20 km s −1 ), only the synthesis method is used in analysis, with solar parameters as inputs, and the initial  sin  ★ set to that estimated from the FWHM.

Testing on J2046+06
We used EBLM J2046+06 to test the output of our combined method, due to having recent and reliable literature parameters determined spectroscopically by Swayne et al. (2021).Using 22 HARPS spectra combined to an SNR ≈ 300, their analysis uses ARES+MOOG, as described by Sousa (2014) and Santos et al. (2013).Differing from our method, they employ only the EW method using the ARES code (Sousa et al. 2007(Sousa et al. , 2015) ) together with the line list described in Sousa et al. (2008b), Kurucz model atmospheres (Kurucz 1993), and the MOOG radiative transfer code (Sneden 1973).
Both SOPHIE HR and HARPS spectra are available for EBLM J2046+06, allowing additionally for testing of continuity between different instruments.
Table 1 displays the results of the pipeline testing on both SOPHIE and HARPS spectra, in addition to the parameters determined by Swayne et al. (2021) and Gaia DR3 GSP-Spec, determined using spectra from the Gaia's Radial Velocity Spectrometer (RVS).(Recio-Blanco et al. 2023).The higher SNR and resolution achieved by HARPS is demonstrated as an advantage here, being closer to the values obtained by Swayne et al. (2021) than those from using the SOPHIE HR mode spectra.With 22 available spectra having an average SNR of 43, the SOPHIE combined spectrum has SNR 131, whereas the HARPS combined spectrum has an SNR of 273, from 27 spectra with average SNR of 57.Despite this, parameters within 2 are returned in both cases when compared to Swayne et al. (2021) in all cases except  mic .Table 1 displays that atmospheric parameters derived from spectra from different instruments agree well with each other, including for lower SNR and resolution SOPHIE spectra.This strengthens the reliability of our methods.Concerning the results from Gaia DR3, our results, and those from Swayne et al. (2021), good agreement is shown.

Atmospheric parameters
We analysed the spectra for all 179 targets in the BEBOP sample with the methods described in Section 4. It took approximately 66 hours of computation time2 .If a star was analysed using spectra from multiple spectrographs, the final adopted stellar parameters are calculated as an inverse-variance weighted average from the individual results.
The results of the full sample are shown in Figure 3, in the form of a log  ★ vs. T eff diagram.SOPHIE HE mode, HR mode, and HARPS spectra are split into three separate panels on the same scale to demonstrate the dispersion of each.

Comparison with Gaia DR3
Gaia DR3 parameters were used in a comparison of our output to literature parameters due to providing results for the majority of our targets.Comparisons were performed for T eff and log  ★ , which were consistently available from both GSP-Phot (Andrae et al. 2023) and GSP-Spec (Recio-Blanco et al. 2023).Figure 4 shows our results versus the Gaia DR3 values, with GSP-Spec results limited to those with a fluxNoise quality flag of Flag 2 or lower, according to Appendix C of Recio-Blanco et al. (2023).These results were split into the North and South samples to check for potential instrumental biases.Although error bars are plotted for all points, it is important to consider that the uncertainties on the Gaia DR3 values represent precision values, whereas our T eff were inflated in quadrature.
Figure 4 shows good agreement between our results and the Gaia T eff values.Considering the North sample, the RMSD between our results and GSP-Phot is 199 K; this is smaller than the mean error on our T eff of 256 K.When comparing our results to those form GSP-Spec, the RMSD is 190 K, again being smaller than the mean error.For the South sample, the benefit of the higher resolution manifests in the RMSD between our results and GSP-Phot reducing  to 184 K, with a mean uncertainty of 199 K.We see poorer agreement when comparing our results from the South sample to those from GSP-Spec, with the RMSD increasing to 266 K in this case.
Section 4.2 discusses the difficulty in obtaining a reliable value of log  ★ from spectra, hence we did not expect to see concrete agreement when comparing our log  ★ to those from Gaia DR3.However, Figure 4 does show that we do not determine any log  ★ values that would be unphysical for our targets, in addition to showing a general positive correlation between our results and results from Gaia DR3.Furthermore, we do not believe the discrepancies to be of great concern in terms of our final results, due to the extensive testing by Mortier et al. (2014) that demonstrates the effect of a changing log  ★ to be insignificant on the derivation of other stellar parameters from spectroscopy.
As discussed by Andrae et al. (2023), Gaia GSP-Phot metallicity values suffer from a systematic underestimation.To counteract this, we used the empirical calibration relation introduced by Andrae et al. (2023) via the python package gdr3apcal3 to determine calibrated GSP-Phot metallicities. Figure 5 shows a comparison of our results with both the results from GSP-Phot and GSP-Spec.Our results have a RMSD from the GSP-Spec values of 0.20 dex, increasing to 0.26 dex when we compare to the GSP-Phot values.

BEBOP colour-magnitude diagram
In Figure 6 we show t he BEBOP sample in a colour-magnitude diagram using Gaia DR3 parallaxes, G band apparent magnitudes, and BP-RP colours (Gaia Collaboration et al. 2016bCollaboration et al. , 2022;;Babusiaux et al. 2022).No reddening was included so some scatter could arise from that, however, the stars have a mean distance of 336 pc, and reddening would thus be minimal overall.We add a colour-bar representing our derived effective temperature for each target.As one would expect the sample follows the main sequence with the hottest stars represented in this figure occupying the upper left of the colourmagnitude diagram, being the bluest and brightest from the sample.The reddest and dimmest stars are shown, as expected, to have the lowest effective temperatures.
Ensuring further that the effective temperatures decrease going from bluer to redder stars can present an excellent opportunity to reveal outliers.Figure 7 separates the data displayed in Figure 6 into five bins of Gaia BP-RP colour, and uses a box-and-whisker plot to show the distribution of effective temperatures in each.Outliers from this are clearly displayed as the red points lying outside of the whiskers.Beginning with the 0.4-0.6 colour bin, the cool outlier is J0954-45.Our results for this target show a  sin  ★ of 30.01 ± 16.01 km s −1 -such a high  sin  ★ could have resulted in blended lines that interfered with the analysis.There are no Gaia DR3 parameters available for J0954-45, hence we could not do a comparison here.The hotter outlier of the first bin corresponds to J1805+09, which has an exceptionally high  sin  ★ of 49.93 ± 9.30 km s −1 , hence the target is also highly susceptible to blended lines.Our T eff for this target was calculated to be 7532 ± 111 K; we can compare this to the Gaia DR3 T eff of 7267 ± 49 K, keeping in consideration that Gaia DR3 errors are precision-only, whereas we have inflated our T eff uncertainties in quadrature to reflect inaccuracies in the atmospheric models used for analysis.The Gaia DR3 result would also place the T eff as an outlier in Figure 7.Further investigation of this target reveals the TIC radius reported to be 1.68 ± 0.07 R ⊙ -together with the pipeline-derived effective temperature for the HE mode of 7375 ± 180 K, Pecaut & Mamajek (2013) places this target in the range of late A to early F type stars.This indicates a possible unsuitability of J1805+09 as part of an FGK dataset.With an HE mode combined SNR of 65, additional observations of the target would be required to reach an SNR of 100 and reduce the uncertainty in parameters.
Moving to the 0.6-0.8bin, the single outlier displayed is J1258-58.This target has no Gaia DR3 results to compare to, and with the combined spectrum reaching an SNR of 136 we do not expect that the results are unphysical.The 0.8 -1.0 colour bin also contains a single outlier, being J1916-04.Our result for the T eff of this target is 6974 ± 206 K, which is in agreement with the Gaia DR3 T eff of 7033 ± 32 K -this gives confidence in our results for this target.The final outlier in the 1.0 -1.2 colour bin is J0525+26; given that the SNR of the combined spectrum is 42, we would require additional data to determine whether this target is an outlier due to its physical conditions or the quality of the data.

Masses, Radii and Ages
The atmospheric parameters of all BEBOP primaries allow us to derive stellar parameters such as mass and radius using stellar evolution models.We applied MIST isochrones (Dotter 2016;Choi et al. 2016) to interpolate the spectral parameters using the isochrones package (Morton 2015).This interpolator utilises a multimodal nested sampler multinest (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2019) ) which allows to sample the input spectral parameters together with photometric and distance information.We used our derived T eff and [Fe/H], separating them for each instrument and instrumental mode, together with the Gaia DR3 parallaxes (Gaia Collaboration 2022), IR colours W1, W2, & W3 from Allwise (Cutri et al. 2021), and the cross matched NIR 2MASS (Skrutskie et al. 2006) colours H, J, K s from the same catalogue.We chose not to include a value for log  ★ as an input parameter due to being unable to determine its reliability.As we use a very precise parallax and a variety of magnitudes, it is also not a crucial parameter when fitting isochrones and evolutionary tracks.We fit for stellar mass, radius, age, distance and extinction (AV).For each target, we sampled 1000 live points and extract the final parameters as median values from the posterior distributions with the errors representing the 16/84 percentile.Following the methodology outlined by Tayar et al. (2022), noise floors of 5% for masses, 4.2% for radii and 20% for ages were added to our uncertainties.
The resulting BEBOP primary mass-T eff diagram, coloured by [Fe/H], is shown in Figure 9.Each point in the plot is coloured by the metallicity of the target -this reveals an abundance of solar-andhigher metallicity targets.It is well established that the occurrence rate of giant planets is increased for metal-rich stars (Sousa et al. 2011).Although not as numerous, many metal-poor targets can also be seen in this plot.These are interesting targets for the study of planet occurrence rates; Mortier et al. (2012) demonstrate that giant planet frequencies around low metallicity stars may not be as diminished as initially expected.Future comparisons of the occurrence rates of giant planets found across the entire range of metallicity targets in the BEBOP sample with those previously studied around single stars would be highly beneficial to theories of planet formation.Having stellar mass values for primary stars in the BEBOP sample will allow for such comparisons to include good constraints on planetary masses.

Spectroscopic and Isochronal Surface Gravities
As highlighted in both Sections 4.2 and 5.1.1,the spectroscopic log  ★ we derive in this work was not expected to be reliable or physically accurate.Our log  ★ values determined using the MIST isochrones better represent the physical conditions of the star.A comparison of the these and the spectroscopic values is presented in Figure 8.
This comparison is analogous to that by Tsantaki et al. (2013), who saw an underestimated when comparison their spectroscopic log  ★ values to those estimates from parallaxes.From Figure 8, it is also apparent that the uncertainties on the MIST log  ★ s are greatly reduced compared to the spectroscopic values.Taking into account these higher precisions, and the general consensus of the literature, we would suggest that any further analysis required surface gravities is done using the MIST log  ★ values.
The full set of stellar parameters is presented in an online machinereadable table.

CONCLUSION
In this work we present the homogeneous stellar analysis of the BE-BOP sample.With effective temperatures ranging from below 5000 K to over 7000 K, the BEBOP sample stretches the limit of what can be achieved with homogeneous analysis.Our spectroscopic method uses the EW method and spectral synthesis in succession to derive all stellar atmospheric parameters, but make use of the computational speed provided by the EW method.We provide a public tool, PAWS,  .that we used to perform our analysis, but is applicable to any spectra of FGK dwarfs.
Our method has been demonstrated to produce reliable parameters for solar-type main sequence stars as was the aim during its production, with the advantage of using the EW method to generate intial parameters for spectral synthesis made clear.Although the input parameters do not have a significant impact on the final derived ones, beginning in a parameter space that is close to the physical one saves  on computation time and allows more iterations to be dedicated to the refinement of well-fitting parameters.Comparison to literature parameters for the BEBOP data set revealed strong agreement throughout this work.
Emphasis is put on a minimum SNR of 100, however in this analysis we find no targets deviate significantly from expected results purely due to low SNR despite 93 of 194 spectra not meeting this.The ability to reliably return stellar parameters for spectra of SNR < 100 allows a wide range of spectra to be analysed that may not be possible to use the equivalent widths method for.As previously discussed in Section 4.4, EW analyses require high SNR spectra.Within the BEBOP sample, only 11% of combined spectra reach an SNR of 200, indicating that a pure equivalent widths analysis of the sample would not be viable.Poor performance at low SNR indicates that the EW part of our analysis should be skipped for particularly low SNR targets, instead using only synthesis to provide stellar parameters.The use of the EW method in such a scenario is unlikely to initiate the synthesis method in a parameter space reflective of the physical situation.As a result, the time-saving aspect of the EW method would be obsolete here.
Stellar parameters for the BEBOP data set from literature are far from homogeneous.The majority of targets have stellar parameters publicly available from only Gaia DR3 and the TIC.Available parameters from these sources are limited, therefore a homogeneous picture of BEBOP targets cannot be created from current literature.In providing essential stellar parameters for every BEBOP target using the same method, with the same inputs, homogeneity is ensured for any further analysis of the BEBOP sample.gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.
This the EW parameters are fed as inputs to the synthesis method, with metallicity fixed due to the proficiency of the EW method in using the iron lines.As  mic cannot be determined via spectral synthesis, it was also necessary to fix this to the value determined as detailed in Section 4.3.1.Results of the full analysis are then saved along with the original data.The final results folder includes the files detailed in Table A1.
PAWS has been set up to offer more functionality than used for the analysis of the paper to allow users to make different choices.These changeable options include: • PAWS can also use non-solar input parameters for the equivalent widths section, depending on initial knowledge of a target's spectral type.As the majority of BEBOP targets were fitted with G type masks by instrumental data reduction, the solar inputs to the equivalent widths method were not changed during their analysis, further adding to the homogeneity of the method.However, allowing for future investigation of a wider variety of targets, K and F type input parameters have been added as an option for PAWS.
• Although recommended to keep this in place to avoid the introduction of contamination to the combined spectrum, the SNR filter can be removed when using PAWS if necessary.
• Depending on the format of the input spectra, it can be selected whether they need to be coadded, or to use a pre-coadded spectrum.Continuum normalisation will occur regardless.
• The instrumental resolution and wavelength range can be customised, to allow for not only the correct resolution to be used, but also to ensure only the desired section of the spectrum is analysed.

APPENDIX B: ADDITIONAL TABLES
Table B1 presents the lines that were removed from the standard iSpec linelist in our analysis.

APPENDIX C: SYNTHESIS-ONLY TARGETS
This paper has been typeset from a T E X/L A T E X file prepared by the author.
Table A1.Files saved by PAWS into the target folders and their descriptions.

Figure 2 .
Figure2.Across all three panels, the TIC T eff value for J0002+47 is represented by the green horizontal line, with its error bars shown as the green shaded region.The purple horizontal line represents the Gaia DR3 T eff , supplied without errors.a) T eff derived from varying SNR spectra of J0002+47 using only the equivalent widths method.b) Effective Temperature derived from varying SNR spectra of J0002+47 using only the synthesis method.c) T eff derived from varying SNR spectra of J0002+47 using both methods subsequently.

Table 1 .Figure 3 .
Figure 3.The BEBOP sample presented as a log  ★ vs T eff diagram.The data are split into the respective source of the spectra, ordered as SOPHIE HE (left), SOPHIE HR (middle), and HARPS (right).

Figure 4 .
Figure 4. Comparisons of our results (y axes) to Gaia DR3 parameters (x axes).Results from the North sample (both HE and HR mode) are the top two plots, and South sample the bottom two plots, with T eff and log  ★ comparisons on the left and right respectively.

Figure 6 .
Figure 6.Colour-Magnitude diagram of the BEBOP sample, using Gaia DR3, with the colour representing our effective temperature.The cooler stars are represented by darker colours, becoming brighter as the effective temperature increases.The expected trend of the brightest, bluest stars (top left) being the hottest, and the dimmer, redder ones (bottom right) being the coolest is seen clearly here.

Figure 7 .
Figure 7. Box-and-whisker plot of the distribution of effective temperatures derived for targets in each Gaia DR3 colour bin, as taken from Figure 6.

Figure 8 .
Figure 8. Spectroscopic log  ★ values compared with those inferred from MIST isochrones, both from this work.

Figure 9 .
Figure 9. Masses obtained from MIST isochrones plotted against T eff , coloured by the metallicity of each target.
research is supported from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement n • 803193/BEBOP), and by a Leverhulme Trust Research Project Grant (n • RPG-2018-418).SPlaSH).This work is supported by the French National Research Agency in the framework of the Investissements d'Avenir program (ANR-15-IDEX-02), through the funding of the "Origin of Life" project of the Grenoble-Alpes University.EW, acknowledges support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, G.A. n. 772293 http://www.asterochronometry.eu)IB thanks the support of the Programme National de Planétologie (PNP) of CNRS-INSU co-funded by CNES.

Table B1 .
Table of parameters and errors derived only by the equivalent widths method (i.e.those used as synthesis inputs).used_linemasks.csvTable of the wavelengths and properties of the iron lines used by the equivalent widths method.params_synth_pipeline.csvTable of the final parameters derived by the pipeline.fitted_line_params.csvTable of the wavelengths and properties of the lines used in synthesis fitted.synthesized_spectrum_pipeline.fitsFITS file of the synthesized spectrum as flux vs wavelength (nm).Lines removed from analysis due to poor chi-squared results from SOPHIE HE mode synthesis fitting.

Table C1 .
Targets for which the EW step was skipped during analysis, categorised by the reasoning for doing so.