© ESO 2019

1. Introduction

Low-mass stars (≤0.6 M_⊙) have historically been challenging to study because of their intrinsic dimness and the low probability of finding them in eclipsing systems. Careful observations of double-lined eclipsing binaries (SB2s) can result in (almost) model-independent mass and radius estimates to a precision better than 3% in some cases (Andersen 1991; Torres et al. 2010). Interferometry can be used to achieve similar results. Interferometrically determined visibility data can be fitted with a curve appropriate for a uniformly illuminated disc to determine an angular diameter. The resultant angular diameter must be corrected for the effects of limb-darkening to obtain the true angular diameter. This correction is subtle in the infrared but more prominent in the visible and depen on stellar atmospheric parameters (see Fig. 1 of Davis et al. 2000). Numerical values to correct a uniformly-illuminated angular diameter are presented in the form of coefficients for a given limb-darkening law (e.g Hanbury Brown et al. 1974). However, the use of finite bandwidths can introduce instrumental effects (known as bandwidth smearing or chromatic aberration in radio astronomy; Bridle & Schwab 1989) or inconsistent angular diameters across the band due to spectral lines (Hanbury Brown et al. 1974).

Boyajian et al. (2012) acquire interferometric observations at the CHARA Array in the near-infrared K′ and H bands (ten Brummelaar et al. 2005) for 21 nearby, single and bright red dwarfs. They measured radii with an uncertainty below 3% and uncertainty in T_eff below 1% and robustly demonstrate that models over-predict T_eff by ∼3%, and under-predict radii by ∼5%. Rabus et al. (2019) used near-infrared long-baseline interferometry with PIONIER at the Very Large Telescope Interferometer along with data from the second Gaia data release (Gaia DR2; Gaia Collaboration 2018) to provide estimates for their linear radii, effective temperatures and luminosities. They show that Gaia underestimates M-dwarf temperatures by ∼8.2% and visually identify a discontinuity in the T_eff-radius plane. They argue this discontinuity (M_* ∼ 0.23M_⋆) arises during the transition from partly convective M-dwarfs to the fully convective domain despite residing in a less-massive regime than what is typically consider as the convective transition. Benedict et al. (2016) use white-light interferometric observations from the Hubble Space Telescope with radial velocity data from the McDonald Observatory to obtain astrometric solutions of M-dwarfs in binary systems. They achieve mass uncertainties as low as 0.4% in some cases (median error of 1.8%) but these objects are not transiting and so the radius of each component cannot be accurately inferred. Obtaining the mass for single interferometric stars often requires photometric calibrations for red-dwarfs (e.g. Henry & McCarthy 1993; Delfosse et al. 2000; Benedict et al. 2016; Mann et al. 2019). The uncertainty attributed to mass-luminoisty relations has decreased significantly in recent years to around 2−3% in redder colours (Mann et al. 2019) but makes it challenging to assess radius inflation for single stars observed with interferometry.

A review of the literature by Chaturvedi et al. (2018) finds 90 stars less massive than 0.6 M_⊙ of which mass and radius measurements are quoted to an accuracy of 10% or better. These measurements show that some low-mass stars are hotter and larger than stellar evolutionary models predict (Zhou et al. 2014; Gómez Maqueo Chew et al. 2014). The physics of low-mass stars is an interesting problem in its own right. M-dwarfs near the hydrogen burning limit are cool enough for their interiors to approach the electron Fermi temperature, resulting in an electron-degenerate gas where a classical Maxwellian description does not apply (Baraffe et al. 2003). Further complications arise from enhanced magnetic fields which may affect the evolution of M-dwarfs (Chabrier et al. 2007). Large magnetic fields are thought to be induced by tidal interaction in close binaries, enhancing rotation and the dynamo mechanism. This inhibits convection in the core and may be responsible for inflating some stellar radii above those predicted by evolutionary models (e.g. Ribas 2006; Torres et al. 2010; Kraus et al. 2011). However, studies of single M-dwarfs with interferometry (Boyajian et al. 2012) and those in double-lined eclipsing binaries (Feiden & Chaboyer 2012) are comparably inflated by around 3% making it unclear whether tidal interactions can be blamed (Spada et al. 2013).

Using spectroscopy to measure the atmospheric composition of a low-mass star is challenging. Low surface temperatures permit molecules (such as TiO, VO, CaH and H₂O) to exist which manifest as a series of broad, and mostly blended lines in the optical part of the spectrum. This problem is eased slightly in the infrared, where there are some un-blended regions. This allows equivalent-width measurements to be made of the Na 1, Ca 1 and H₂O-K2 index in the K band (Rojas-Ayala et al. 2012). The most common technique for determining the temperature of M dwarfs is comparing spectra to model atmospheres (Gaidos et al. 2014). The infrared flux method has also been extended from FGK stars to M dwarfs (Casagrande et al. 2008), however significant statistical deviations from interferometric temperatures have been noted (Mann et al. 2015).

The Wide Angle Search for Planets (WASP; Pollacco et al. 2006) is a survey for 0.8–2R_Jup objects transiting solar-like stars. Objects in this radius range can have masses which span three orders of magnitude, from Saturn-like planets to M-dwarfs. Consequently, WASP photometry has been used to identify hundreds of FGK stars with transiting M dwarf companions as a by-product of its successful exoplanet search. These systems are termed EBLMs (eclipsing binary, low-mass). We have invested considerable effort into characterising these systems, including hundreds of hours of telescope time to measure their spectroscopic orbits. A primary aim of the EBLM project is to improve our understanding of low-mass stars using accurate mass and radius measurements for transiting companions to FGK stars (Triaud et al. 2013a, 2017; Gómez Maqueo Chew et al. 2014; von Boetticher et al. 2017, 2019). So far, 13 systems have absolute parameters measured from follow-up photometry and a further 118 EBLMs presented by Triaud et al. (2017) have secondary masses. Of these, two stars (WASP-30B and J1013+01) appears to be inflated, and a third (J0113+31) is measured to be ∼600 K hotter than expected. This latter result comes from an analysis of the secondary eclipses in this system using infrared photometry. Gómez Maqueo Chew et al. (2014) were unable to explain this result with any of the mechanisms discussed in their work, including tidal dissipation. A similar result was also seen for KIC1571511 (Ofir et al. 2012) using high-precision optical photometry from the Kepler space telescope. In general, M-dwarfs measured within the EBLM project are consistent with stellar models within a few per cent (von Boetticher et al. 2019).

It is clear that we do not fully understand stars at the bottom of the main sequence and, by implication, the planets orbiting them. With interest in low mass stars increasing as a result of recent exoplanet discoveries (e.g. TRAPPIST-1 and Proxima Centauri), more effort needs to be invested into understanding what makes some low-mass stars anomalous so that we can better understand the myriad of exoplanet systems that will be discovered with TESS (Ricker et al. 2014), MEarth (Charbonneau et al. 2008), Speculoos (Gillon et al. 2013) and eventually PLATO (Rauer et al. 2016). In this paper we present high quality lightcurves and spectroscopic orbits of five EBLM systems and used these to measure the masses and radii of the stars in these systems to a precision of a few per cent in some cases. Section 2 describes the origin and reduction of data used in this work, Sect. 3 details how system parameters we extracted and we present and discuss our results in Sect. 4. We discus predictions from evolutionary models in Sects. 5 and 6 and possible sources of systematic uncertainty in Sect. 7.

2. Observations and data reduction

This section describes the data we have used to measure the physical properties of five EBLM systems (J2349−32, J2308−46, J0218−31, J1547+39, J1436−13) discovered by the WASP project (Sect. 2.3). The quality of WASP photometry is not good enough to measure masses and radii of the components to the desired precision of a few per cent so we obtained more precise follow-up photometry to provide improved size estimates of both EBLM components. A summary of observations can be found in Table 1 and the dates and time of spectroscopic observations are detailed in Table C.1.

The larger stars in these systems are far brighter than their M-dwarf companions and only the reflex motion of the primary star can be measured. These radial velocity measurements provide a constraint on the mass ratio and the analysis of the transit provides information on the relative sizes of the stars, but one additional independent constraint is needed to uniquely determine the scale of the binary system. The additional constraint we used is the mass of the primary star based on its density, effective temperature (T_eff) and metallicity ([Fe/H]) estimated from stellar models. Values of T_eff and [Fe/H] come from analysis of the coadded spectra for each star (Sect. 3.2).

2.1. Photometric colours used for SED fitting

Photometric colours for each target was extracted from the following catalogues: B_T and V_T magnitudes from the Tycho-2 catalogue (Høg et al. 2000); B, V, g′, r′ and i′ magnitudes from data release 9 of the AAVSO Photometric All Sky Survey (APASS9; Henden et al. 2016; J, H and K_s magnitudes from the Two-Micron All-Sky Survey (2MASS; Skrutskie et al. 2006; i′, J and K magnitudes from the DEep Near-Infrared Southern Sky Survey (DENIS; Epchtein et al. 1997). The reddening maps by Schlafly & Finkbeiner (2011a) were used to estimate the total line-of-sight extinction in the direction of each target, E(B − V)_map. Values of E(B − V)_map were calculated using the NASA/IPAC Extragalactic Database (NED) operated by the Jet Propulsion Laboratory, California Institute of Technology¹. Not all EBLMs have photometry in all catalogues; those that do are reported in Table 1.

2.2. Gaia colours and Interferometry used to estimate distance and evolutionary status

The second Gaia data release provides mean flux counts in three bands – G, G_BP and G_RP. The G-band has a wider wavelength coverage and is optimised to determine astrometric solutions. The mean magnitudes G_BP and G_RP provide a slice through the spectral energy distribution of stars and reveal how red or blue a star is. We obtained the mean G, G_BP and G_RP magnitudes along with parallax measurements for all nine EBLM systems from Gaia DR2 using the Gaia archive² (Table 2). There is evidence of systematic offsets in parallax measurements from Gaia DR2 (e.g. Stassun & Torres 2018) which is likely correlated with on-sky positions (α & δ), G and G_BP − G_RP (Lindegren et al. 2018). We added a systematic zero-point offset of −0.082 mas to the parallax and added an additional 0.033 mas in quadrature to the quoted parallax uncertainty (Stassun & Torres 2018). We plotted the position of all EBLMs in the M_G–G_BP − G_RP plane using data from 160 randomly selected source fields (Fig. 1) filtered using Eqs. (1) and (2) from Arenou et al. (2018).

Fig. 1. MG vs GBP − GRP plane for 160 randomly selected source fields (black) filtered using Eqs. (1) and (2) from Arenou et al. (2018). The EBLMs used in the work are also plotted and coloured appropriately.

2.3. WASP photometry for initial transit parameters and ephemerides

The WASP survey (Pollacco et al. 2006) operates two survey instruments: one at the South African Astronomical Observatory (SAAO), South Africa, and another at the Observatorio del Roque de los Muchachos, La Palma. Each instrument consists of an equatorial fork mount with eight cameras with 200 mm lenses and 2k × 2k CCD detectors. Each camera coveres approximately 64 square degrees per exposure. The data are processed by a detrending algorithm which was developed from the SYSREM algorithm of Tamuz et al. (2005) and that is described by Collier Cameron et al. (2007). In July 2012, lenses on the southern installation (WASP-South) were changed to 85 mm with f/1.2 to search for brighter exoplanet hosts (Smith & WASP Consortium 2014). Data from 85 mm lenses were not used in this study.

Photometry from the WASP cameras can suffer from a large amount of scatter due to clouds, instrumental artefacts, scattered light and other non-optimal observing conditions. We cleaned the data by removing points that were not detrended in the standard WASP reduction pipeline and removed points more than 0.5 mag. from the median magnitude of each star. Additional cleaning of the light curve was done by comparing each night of data to a phase-folded light curve binned into 500 phase bins. Any measurement 3-σ or more from the mean in each bin was excluded. The entire night of data was excluded if more than a quarter of the night’s data was excluded this way or if there are fewer than ten observations. The binned light curve is then inspected by eye to further exclude bad data points.

2.4. SAAO 1 m follow-up transit photometry

The SAAO hosts an equatorial-mounted 1-m telescope built by Grubb and Parsons that is equipped with an STE4 CCD camera with 1024 × 1024 pixels. This camera was operated in 2 × 2 binning mode to reduce readout time. J2349−32 was observed on 18 October 2016 and J2308−46 on 12 October 2016 using I (exposure time of t_exp = 50 s) and R (t_exp = 40 s) Bessel filters. J1436−13 was observed on 23 April 2017 in the R (t_exp = 40 s) Bessel filter. Photometry was extracted using standard aperture photometry routines (Southworth et al. 2009) and uncertainties were estimated from photon counting statistics. A by-eye approach was used to clean the light curve and select the best comparison star in the 5′×5′ field. A slow variation in differential magnitude with time was observed corresponding to changes in the effective airmass. To correct for this, out-of-transit regions were manually selected and we used the IDL/AMOEBA³ routine to fit a polynomial which minimised the square of the magnitude residuals. The lightcurves were divided by the fitted polynomials to normalise to zero differential magnitude.

2.5. HAO follow-up transit photometry

Optical photometry for J1847+39 was obtained at the Hereford Arizona Observatory, Arizona (HAO). Three separate transits were observed with a Meade 14-inch LX200GPS telescope. The first was obtained with the clear blue-blocking filter (CBB) on 9 October 2009 with t_exp = 100 s. The second was with a g′ filter on 18 May 2011 with t_exp = 60 s. The last was with a z′ filter on 15 June 2010 with t_exp = 60 s. All observations were made in a binned mode of operation with a pixel scale of 1.5″. The atmospheric seeing and telescope tracking produced a point-spread function with a full-width half-maximum of ∼4″ (2.5 pixels). We chose apertures with radii of 4 pixels (∼6″) to encapsulate the point-spread function of J1847+39 and comparison stars. Aperture photometry was extracted using standard photometry routines with systematic trends removed and outliers rejected. Transmission information of the telescope throughput, atmosphere, filter and CCD⁴ was used to calculate the final transmission of HAO with the CBB filter (see Fig. 2). We used the four-parameter limb-darkening look-up table for the K2 pass band instead of the CBB filter due to the similarity in final transmission since we do not have access to a four-parameter look-up table for the CBB filter. In Sect. 7.3 we fit light curves using the two-parameter quadratic limb-darkening instead of the Claret law. The final response function in Fig. 2 was used along with estimates of stellar atmospheric parameters from Sect. 3.2 to calculate quadratic coefficients using LDTK (Parviainen & Aigrain 2015).

Fig. 2. Response function of the HAO+CBB instrument. The atmospheric transmission is plotted in black, the transmission of the HAO telescope in blue-solid, the CBB filter in green and CCD response in yellow. The final response of HAO-1 with the CBB filter is plotted in red-dashed along with the Kepler transmission (blue-dashed). The atmospheric transmission line originated from equations for Rayleigh, aerosol and ozone extinction vs. wavelength for Palomar Observatory (Hayes & Latham 1975). Coefficients were adjusted until they agreed with observations of extinction at HAO over a few dates (2018, priv. comm.).

The discrepancy between the K2 and HAO+CBB pass-band differ in the blue where the limb-darkening is most significant. The validity of this assumption hinges on the fact that (1) the intensity across the stellar disc is similar in both the K2 and HAO+CBB bandpass and (2) the difference in limb-darkening coefficients for each band in negligable. The first assumption was tested using LDTK to synthesise intensity profiles for J1847+39 across the stellar disc for each pass-band and calculate the discrepancy as a function of γ (the angle between a line normal to the stellar surface and the line of sight of the observer; Fig. 3). The K2 pass-band emits 2.5% less flux than what would be observed with HAO+CBB towards the limb. The second assumption was tested by calculating quadratic limb-darkening coefficients for the K2 pass-band to be u₁ = 0.496 ± 0.050, u₂ = 0.157 ± 0.050 and for the HAO+CBB pass-band to be u₁ = 0.468 ± 0.050, u₂ = 0.148 ± 0.051. These are comparable within 1-σ and so adopting the K2 pass-band for J1847+39 will have a negligible effect on the transit shape.

Fig. 3. Difference in theoretical intensity acoss the stellar disc for the HAO+CBB filter and the Kepler/K2 Filter as a function of the angle between a line normal to the stellar surface and the line of sight of the observer (γ) for J1847+39.

2.6. CTIO follow-up transit photometry

J0218−31 was observed on 14 November 2010 with the CTIO-0.9-m telescope and Tek2K CCD camera. The detector consists of a 2K × 2K array of 15 μm pixels placed at Cassegrain focus giving a 0.4″ pixel⁻¹ plate scale. Thus the entire array projects to a 13.7′ FOV. The observed signal is fed into four amplifiers causing the raw images to have a quadrant effect with the readnoise between 3.9 and 4.5 e⁻ and gain of 2.5−2.8 e⁻/ADU, depending on the amplifier. The detector has a readout time of 39 s and a 60 000 count well depth before non-linearity sets in when using 1 × 1 binning mode. J0218−31 and the surrounding field were monitored throughout the night using the Sloan griz filter set alternating and continuously between all four filters. Exposure times were chosen to maximise the flux in the target star and nearby reference stars while keeping the peak pixel value in J0218−31 below 60 000 counts (well depth). The telescope was defocused to allow for longer exposure times to build up signal in the fainter reference stars without saturating J0218−31. We adopted an exposure time of 10 s for the g′, r′, and i′–band observations and longer exposures of 15 s in the z′ filter where the detector is less sensitive. An overall light curve cadence of ≈3.3 min was achieved in each filter accounting for the exposure times, the read out time, and filter changes. The light curves were created from approximately 75 images taken in each filter during the single observing night.

A set of 11 bias calibration frames and 11 dome flat fields in all four filters were obtained at the beginning of the observing night. The images were processed in a standard way using routines written by L. Hebb in the IDL programming language. Each of the four amplifiers was processed independently. All object and calibration frames were first overscan corrected (by subtracting a line-by-line median overscan value), bias subtracted and then trimmed. Stacked bias images were created by averaging all bias frames observed each night and subtracted from all science and flat-field frames. All dome flats were averaged into a single dome flat in each filter and then applied to the trimmed and bias-corrected science images.

Source detection and aperture photometry were performed on all processed science images using the Cambridge Astronomical Survey Unit catalogue extraction software (Irwin & Lewis 2001). The software has been compared with SExtractor (Bertin & Arnouts 1996) and found to be very similar in the completeness, astrometry and photometry tests⁵. This photometry software was applied to all processed images of J0218−31. Adopting conservative parameters to define the detection threshold, the target star and dozens of fainter stars in the field were detected in each image. Aperture photometry was performed on all detected stars using a five pixel radius circular aperture, which was selected to match the typical seeing. Five bright, non-variable reference stars were selected from the many detected stars and used to perform differential photometry on the target star. In each image, the flux from all reference stars was summed into a single super comparison star that was divided by the aperture flux from J0218−31 and converted to a differential magnitude.

2.7. CORALIE spectra used for radial velocities and atmospheric parameters

CORALIE is a fibre-fed échelle spectrograph installed on the 1.2-m Leonard Euler telescope at the ESO La Silla Observatory and has a resolving power R = 50 000–60 000 (Queloz et al. 2001; Wilson et al. 2008). The spectra used in this study were all obtained with an exposure times between t_exp = 600 − 900 s. Observations of J0218−31 include spectra obtained through the primary eclipse that show the Rossiter-McLaughlin effect. The spectra for each star were processed with the CORALIE standard reduction pipeline (Baranne et al. 1996). Radial velocity measurements were obtained using standard cross-correlation techniques (using numerical masks) and checked for obvious outliers (Triaud et al. 2017). Each spectrum was corrected into the laboratory reference frame and co-added onto a common wavelength range. Maximum and median filters were applied to identify continuum regions which were fitted with spline functions (one every nm) to normalise the spectra (a standard function within ISPEC v20161118; Blanco-Cuaresma et al. 2017).

2.8. INT spectra used for radial velocities and atmospheric parameters

Spectra for J1847+39 were obtained using the intermediate dispersion spectrograph (IDS) mounted on the 2.5-m Isaac Newton Telescope (INT) at the Roque de Los Muchachos Observatory. The 235 mm camera and EEV10 CCD detector was used with the H1800V grating to obtain spectra in a small region around the Hα line with R ≈ 10 000⁶. A total of ten spectra were obtained for J1847+39 with an exposure time t_exp = 600 − 900 s. Radial velocity measurements were extracted using cross-correlation routines provided within ISPEC. A synthetic F0 spectrum was used as a template with a mask applied to the core of the Hα line. A Gaussian function was fitted to the peak in each cross-correlation function to obtain the radial velocity measurement (the peak of the Gaussian function), and uncertainty (standard deviation of the Gaussian function). Each spectrum was corrected into a laboratory reference frame and co-added onto a common wavelength range. The relatively small wavelength range does not permit the use of maximum and median filters to normalise the spectra. Suitable continuum regions were identified by-eye to normalise the spectrum using a second-order polynomial fit by least-squares.

2.9. Lucky imaging used to identify nearby companions

The lucky-imaging technique (e.g. Law et al. 2006) was used to obtain high-resolution images of J2308−46, J2349−32, J0055−00, J1652−19 and J2217−04 in July 2017, in order to search for stars contributing contaminating light, as well as potential bound companions to the eclipsing binaries. The observations were conducted using the Two Colour Instrument (TCI) on the Danish 1.54-m Telescope at La Silla Observatory. The TCI consists of two Electron Multiplying CCDs capable of imaging simultaneously in two pass bands at a frame rate of 10 Hz, with a 40″ × 40″ field of view. The “red” arm has a pass band similar to a combined i + z filter or the Cousins I filter, whilst the “visible” arm has a mean wavelength close to that of the Johnson V filter. A detailed description of the instrument can be found in Skottfelt et al. (2015) and the lucky imaging reduction pipeline is described by Harpsøe et al. (2012).

The observations and data reduction were carried out using the method outlined in Evans et al. (2018), and is briefly described here. Both targets were observed for 170 s. The raw data were reduced automatically by the instrument pipeline, which performs bias and flat frame corrections, removes cosmic rays, and determines the quality of each frame, with the end product being ten sets of stacked frames, ordered by quality. The data were run through a custom star-detection algorithm that is described in Evans et al. (2018), which is designed to detect close companion stars that may not be fully resolved.

3. Methods

This section describes the methods we have used to analyse our data in order to measure the masses and radii of each EBLM system. Our method shares some similarities to von Boetticher et al. (2019); we used the same spectroscopic analysis routine from Gill et al. (2018) which was confirmed with SED fitting, the same light curve model and Bayesian sampling routines for the orbital solution (although some fitted parameters are different). We used a modified sampler from von Boetticher et al. (2019) to measure the masses and radii of each EBLM system from common stellar models. Our approach to quantifying inflation is also different to von Boetticher et al. (2019), both in terms of the models we used and approach. For completeness, we describe our methods in the following sections.

3.1. SED fitting

Empirical colour–effective temperature relations were used used to estimate the effective temperature of the primary star in each system. These were used to complement our spectroscopic analysis and to provide a measurement of reddening. They were not used to interpolate stellar models or inform limb-darkening coefficients. We also assume that the flux contribution from the M-dwarf companion is negligible compared to the primary star (see Sect. 3.5.3).

Our model for the observed photometry has the following parameters – $g_0^{\prime }$: the apparent g′-band magnitude corrected for extinction; T_eff, the effective temperature; E(B − V), the reddening to the system; and σ_ext, the additional systematic error added in quadrature to each measurement to account for systematic errors. For each trial combination of these parameters the empirical colour – effective temperature relations of Boyajian et al. (2013) were used to predict the apparent magnitudes of the star in each of the observed bands. The transformation between the Johnson and 2MASS photometric systems is the same as Boyajian et al. (2013). The Cousins I_C band was used as an approximation to the DENIS Gunn i′ band and the 2MASS K_s as an approximation to the DENIS K band (see Fig. 4 of Bessell 2005). Table 3 of Bessell (2000) was interpolated to transform the Johnson B, V magnitudes to Tycho-2 B_T and V_T magnitudes. This assumed that the extinction in the V band is 3.1 × E(B − V). Extinction in the SDSS and 2MASS bands is calculated using A_r = 2.770 × E(B − V) from Fiorucci & Munari (2003) and extinction coefficients relative to the r′ band taken from Davenport et al. (2014).

The reddening maps by Schlafly & Finkbeiner (2011b) were used to estimate the total line-of-sight extinction in the direction of each target, E(B − V)_map. This value is used to impose the following (un-normalised) prior on Δ = E(B − V)− E(B − V)_map:

$$\begin{array}{c}\hfill P(\mathrm{\Delta })=\{\begin{array}{cc}1\hfill & \mathrm{\Delta }\le 0\hfill \\ exp(-0.5{(\mathrm{\Delta }/0.034)}^2)\hfill & \mathrm{\Delta }>0\hfill \end{array}\end{array}$$

The constant 0.034 is taken from Maxted et al. (2014) and is based on a comparison of E(B − V)_map to E(B − V) determined using Strömgren photometry for 150 A-type stars. We used EMCEE (Foreman-Mackey et al. 2013) to sample the posterior probability distribution (PPD) for our model parameters. EMCEE uses affine-invariant ensemble sampling (parallel stretch move algorithm; Goodman & Weare 2010) to split Markov chains into sub-groups and update the position of a chain using the positions of chains in the other sub-groups. The algorithms affine-invariance can cope with skewed probability distributions and generally has shorter autocorrelation times than a classic Metropolis-Hastings algorithm. The empirical colour–temperature relations we have used are valid over the approximate range T_eff = 3450 K−8600 K. Between these limits uniform priors were used on the values of T_eff. We also used uniform priors for $g_0^{\prime }$. We evolved 10 000 steps from 100 walkers as a burn-in. A further 10 000 steps were drawn and the step with the highest likelihood value is selected, with uncertainties equal to the standard deviation of each parameter in the second chain. An example posterior probability distribution (PPD) for J2308−46 is shown in Fig. 4; the PPDs for the other targets are shown in the Appendix A. The residuals of each fit to all of the EBLMs (observed magnitudes – calculated magnitudes) are shown in Fig. 5.

Fig. 4. Posterior probability distribution of EBLM J2349−32 from photometric fitting. The 1-σ contour are shown.

Fig. 5. Difference between observed and calculated magnitudes for each EBLM.

3.2. Spectroscopic analysis

In this section we describe how we measured T_eff, [Fe/H], V sin i and log g from the CORALIE spectra or INT spectra.

3.2.1. CORALIE – wavelet analysis

We re-sample between 450−650 nm with 2¹⁷ values and co-add the spectra. We calculate the wavelet coefficients W_{i = 4 − 14, k} (see Fig. 2 of Gill et al. 2018 for visual justification of our choice of wavelet coefficients) and fit the same coefficients with model spectra in a Bayesian framework. We initiated 100 walkers and generated 100 000 draws as a burn-in phase. We generated a further 100 000 draws to sample the PPD for T_eff, [Fe/H], Vsin i and log g. Gill et al. (2018) note an [Fe/H] offset of −0.18 dex which we correct for by adding 0.18 dex to the PPD for [Fe/H]. They also note a significant trend in log g with T_eff which we also correct for using their Eq. (9). The wavelet method for CORALIE spectra can determine T_eff to a precision of 85 K, [Fe/H] to a precision of 0.06 dex and V sin i to a precision of 1.35 km s⁻¹ for stars with V sin i ≥ 5 km s⁻¹. However, measurements of log g from wavelet analysis are not reliable beyond confirming dwarf-like gravity (log g ≈ 4.5 dex). Subsequently, we fit the wings of the magnesium triplets with spectral synthesis by fixing T_eff, [Fe/H] and V sin i and changing log g until an acceptable fit was found.

3.2.2. INT – synthesis

INT observations of J1847+39 are unsuitable for wavelet analysis since there only a small wavelength region around the Hα line was observed. The spectral synthesis technique was used to measure T_eff from the wings of the Hα line and [Fe/H] from a limited number of Fe lines around the Hα line. There are no gravity sensitive lines visible in the INT spectra and so we have not attempted to estimate log g from these spectra.

3.3. Ephemerides and first estimates for transit parameters

WASP photometry was used to obtain first estimates for the parameters of the transit and to determine a prior for the ephemerides; it was not used to determine the physical properties of the M-dwarf companions as transit depths can be unreliable. Using the framework of Beatty et al. (2007), we obtain first order approximations to the ratio of semi-major axis, a, and the radius of the primary star, R_⋆, using the width of the transit, Δt_tr,

$$\begin{array}{c}\hfill \frac{R_{\star }}{a}\approx \pi \frac{\mathrm{\Delta }{t}_{\mathrm{tr}}}{P},\end{array}$$(1)

and the ratio of the radii, k, can be estimated from

$$\begin{array}{c}\hfill k=\frac{R_2}{R_{\star }}\approx \sqrt{\mathrm{\Delta }m},\end{array}$$(2)

where R₂ is the radius of the M-dwarf companion and Δm is the depth of transit caused by the M-dwarf. To measure P and T₀, we used the method of Kwee & van Woerden (1956) to compute accurately the epoch of minimum of each complete eclipse in the WASP photometry. Bayesian sampling was used to minimise the correlation between the uncertainties in T₀ and P. We generated 20 000 draws from 12 walkers and selected the reference period and time of minimum from the solution with the highest log-likliehood. The uncertainty for each parameter was estimated from the standard deviation of each parameters PPD and was used as a prior in the orbital fit (Sect. 3.5). We inspected the difference between calculated models and observed epochs to concluded that there is no evidence of transit-timing variations for any of the five EBLMs.

3.4. Rotational modulation

Each system has thousands of observations from the WASP survey which have been taken over many years. Consequently, it is possible to measure variations in the light-curve caused by spot coverage or tidal interactions. We used the method outlined in Maxted et al. (2011) to search the WASP photometry for frequencies attributed to rotational modulation. Each season of photometry is treated separately and in-transit data were excluded. We inspected the periodogram and false-alarm probabilities (FAP) for each system to assess the reliability of any detected periods. The false-alarm probabilities were calculated using the method of Press & Rybicki (1989). We also phase-folded the light-curve at any detected period to check for ellipsoidal variation. An example periodogram for J0218−31 can be seen in Fig. 6. The Lomb-Scargle periodogram for each season for all targets can be found in Appendix B.

Fig. 6. Power spectrum of J0218−31 using 3 seasons of WASP photometry.

3.5. Combined radial velocity and light curve fit

We fited all follow-up photometry (from SAAO, CTIO and HAO) and radial velocity measurements simultaneously to obtain the final orbital solution for each system. We performed a χ² fit in a Bayesian framework to estimate the PPD of each parameter in the vector model. The vector model of parameters includes photometric zero-points for each ith light-curve –zp_i, R_⋆/a, k, the impact parameter –b = acos(i)/R_⋆, T₀, P, the limb-darkening temperature –T_eff, ld, the semi-amplitude of radial velocity measurements –K, the systematic radial velocity – γ and the change in systematic radial velocity with time –d(γ)/dt. First estimates of R_⋆/a, k, T₀ and P were obtained from photometry described in Sect. 3.3. The first estimate of T_eff, ld was the spectroscopic value of T_eff from Sect. 3.2. Instead of fitting the argument of the periastron (ω) and the eccentricity (e), we choose to use $f_{\mathrm{c}}=\sqrt{e}cos\omega$ and $f_{\mathrm{s}}=\sqrt{e}sin\omega$ since these have a uniform prior probability distribution and are not strongly correlated with each other. We also include a jitter term (σ_J) to account for spot activity which can introduce noise in to the radial velocity measurements (Ford 2006). We used T_eff, ld to interpolate coefficients for the Claret limb-darkening law (provided with the python package ELLC, Maxted 2016) using fixed values of [Fe/H] and log g from Sect. 3.2. The stellar metalicity and surface gravity are fixed when interpolating limb-darkening coefficients as varying them has a second-order effect relative to the effective temperature. We used a Gaussian prior for T_eff, ld using the value of T_eff from Sect. 3.2 with width equal to the uncertainty of T_eff measurements from the wavelet method, 85 K. The follow-up photometry for each system is modelled with ELLC assuming detached and spherical star-shapes. Gaussian priors for T₀ and P from Sect. 3.3 were used to constrain the ephemerides of each system.

We compare these models to data using a Bayesian framework with the likelihood function ℒ(d|m) = exp(−χ²/2), with

$$\begin{array}{cc}\hfill {\chi }^2={\displaystyle \sum _{i=1}^{N_{\mathrm{mag}}}\frac{{(m_i-m_{\mathrm{model}})}^2}{{\sigma }_{m_i}^2}}& +{\displaystyle \sum _{i=1}^{N_{\mathrm{r}\mathrm{v}}}\frac{{(r{v}_i-r{v}_{\mathrm{model}})}^2}{{\sigma }_J^2+{\sigma }_{{\mathrm{r}\mathrm{v}}_i}^2}}\hfill \\ \hfill & +\frac{{(T_{\mathrm{eff},\mathrm{ld}}-T_{\mathrm{eff}})}^2}{{\sigma }_{T_{\mathrm{eff}}}^2}·\hfill \end{array}$$(3)

Here, m_i and rv_i represent the ith measurement of magnitude and radial velocity with standard errors σ_mi and σ_rvi, respectively. We initiated 50 walkers and generated 50 000 draws, after an initial burn-in phase of 50 000 draws using EMCEE. We initially selected the model with the highest value of ℒ(d|m) from the PPD to extract the best-fitting model parameters. For J2308−46 and J1847+39 we find these values to be up to 1-σ away from the median value of each parameters PPD (Fig. 7), and so we chose to use the median value from each parameters PPD instead. The uncertainties were calculated from the largest difference between the median and the 16th and 84th percentile of the cumulative PPD for each parameter from the second chain.

Fig. 7. Orbital solution of J2349−32 with 1-σ contours plotted. We show the median value of each parameter from the Markov chains (green) along with the the solution with the highest log-likliehood (red).

3.5.1. Rossiter-McLaughlin

We obtained radial velocity measurements of J0218−31 during the primary transit which display variations in radial velocity caused by the Rossiter-McLaughlin effect. Unlike Triaud et al. (2017) who excluded such measurements when determining the spectroscopic orbit, we fit these measurements simultaneously with out-of-transit radial velocities and follow-up transit photometry. The orbital fit for this system required two more de-correlated parameters, $\sqrt{Vcosi}sin\lambda$ and $\sqrt{Vsini}cos\lambda$, where λ is sky-projected angle between the orbital and stellar rotation angular momentum vectors.

3.5.2. Star shapes

We find ellipsoidal variations in the WASP photometry of J2308−46 which required the use of Roche geometry to estimate the initial transit parameters from the WASP photometry. The follow-up photometry of J2308−46 was fitted using a spherical star shape with the assumption that there is only a small amount of out-of-eclipse photometry, which was detrended. A caveat is that the spherical volume of the star will not be the same as the volume of the triaxial ellipsoid used to approximate its shape with ELLC. We assessed the magnitude of this problem by comparing the models for J2308−46 where both stars are described by spheres to those where both stars are described using Roche models (Fig. 8). We find a maximum difference of ≈0.1 ppm which is far below the white-noise level (a few thousand ppm) and so we do not attempt to correct for this. The final orbital solution for all EBLMs assumes detached and spherical star-shapes and does not use Roche geometry.

Fig. 8. Difference between the spherical model and Roche model of J2308−46 using ELLC.

3.5.3. M-dwarf luminosity and transit depth

The orbital solution assumes that the transit depth is not modified by the luminosity of the transiting M-dwarf. The justification of this assumption requires some foresight of the results (Table 3). The flux contribution depends on the transmission profile observations were made in and the stellar atmospheric parameters of each star. The modification of the primary transit depth further depends on the ratio of the radii, k. For the primary star, these quantities are estimated from spectral analysis. For the M-dwarf companion, we used log g₂ from Table 3 which was determined from the orbital solution and assume the metalicity is the identical to the host star under the premise they both formed in the same molecular cloud. The temperature of the M-dwarf companion was estimated from the MESA stellar evolution models (Dotter 2016; Choi et al. 2016; Paxton et al. 2011, 2013, 2015) using the nominal mass, age and composition. For each band-pass in each EBLM system, we used the following procedure:

1. We interpolated high resolution PHOENIX model spectra (Husser et al. 2013) for each star.

2. The model spectrum for the host and the M-dwarf were convolved them with transmission profiles of the filters that their respective photometric transits were observed in. The ratio of the M-dwarf flux to the primary star flux (F₂/F_⋆) could then be calculated.

3. The surface brightness ratio (k² × F₂/F_⋆) was calculated along with the change in primary transit depth accounting (Δ Depth).

This process was repeated 1000 times for each transmission profile of each EBLM system by perturbing nominal measurements by their respective uncertainties in Table 3. The median value of F₂/F_⋆ and Δ Depth were adopted as nominal measurements with uncertainties equal to the standard deviation of all draws (Table 2).

EBLMs J2349−32, J2308−46 and J0218−31 have primary transit depths which are less than 100 ppm shallower when accounting for the luminosity of M-dwarf companions. J1847+39 and J1436−13 have primary transits around 400 ppm shallower due to the redder filters we obtained transit photometry with and higher values of k. The larger relative uncertainty for Δ Depth originates from a poorly constrained ratio of radii (k) due to a high impact parameter. We assume this effect is negligible and we do not correct for the light of the M-dwarf in this work.

3.6. Mass and age estimates

Breaking the degeneracy between the mass of the primary star – M_⋆, the M-dwarf companion – M₂, and the age of the system – τ, is non-trivial. One approach by Hebb et al. (2009) uses Kepler’s equation to estimate the density of the primary,

$$\begin{array}{c}\hfill {\rho }_{\star }=\frac{3\pi }{G{P}^2}{\left(\frac{a}{R_{\star }}\right)}^3-\frac{3M_2}{4\pi R_{\star }^3},\end{array}$$(4)

and then combine it with measurements of T_eff and [Fe/H] to interpolate between stellar models for M_⋆ and τ. Typically this is repeated with a better estimate of M_⋆ until the solution converges iteratively. Another approach uses empirical mass and radius calibrations (Southworth 2011; Torres 2013) to obtain M_⋆ and R_⋆. These are combined with k and Eq. (4) to obtain M₂ and R₂. Yet another approach by Triaud et al. (2013b) is to mix the two methods while fitting alongside orbital parameters. The mass function (Hilditch 2001) can be expressed in terms of radial velocity parameters,

$$\begin{array}{c}\hfill f(m)=\frac{{(M_{2}sini)}^3}{{(M_{\star }+M_2)}^2}={(1-e^2)}^{\frac{3}{2}}\frac{P{K}^3}{2\pi G},\end{array}$$(5)

where i is the inclination of the orbit and G is the gravitational constant. The middle and right part of Eq. (5) can be equated and solved numerically for M₂ assuming an initial guess of M_⋆ from empirical calibrations. Stellar models are interpolated to give a new estimate of M_⋆. The better value of M_⋆ can be used to iteratively solve Eq. (5) and generate better estimates of M_⋆ and M₂ until a solution is converged upon (Triaud et al. 2013b). A final method relies on three assumptions: (1) the circularisation timescale (τ_circ) is much shorter than τ, (2) the rotation is synchronised (τ ≫ τ_syn) and (3) that rotational and orbital inclination are the same. Under these assumptions it is possible to directly calculate the mass and radius of both components (see Eqs. (14)–(17) of Beatty et al. 2007 or Eqs. (2)–(5) of Zhou et al. 2014). We used Eqs. (2)–(5) of Zhou et al. (2014) to estimate the masses and radii for the EBLMs presented in this work (Table 3). We find a significant discrepancy with our method likely arising from the lack of synchronicity in these systems. Furthermore, it is not enough to assume tidal circularisation and synchronisation from a short orbital period alone (Fleming et al. 2019).

To estimate the mass and age of the primary star we combined the atmospheric parameters (Sect. 3.2) and the best fitting orbital solution (Sect. 3.5) and interpolate between evolutionary models computed with the GARSTEC stellar evolution code (Weiss & Schlattl 2008). We made no assumptions regarding circularised or synchronised orbits. We used a modified version of the open-source code BAGEMASS (Maxted et al. 2015a) tailored exclusively for EBLM systems (EBLMMASS). EBLMMASS uses the jump parameters of age, primary mass (M_⋆), the initial iron abundance in dex [Fe/H]_i, M₂ and the full-width half maximum of the transit w. The vector of observed parameters is given by d = (f(m), T_eff,  log L_⋆,  [Fe/H]_s, R_⋆/a, w) where log L_⋆ is the luminosity of the primary star and [Fe/H]_s is the surface metal abundance in dex and w is the transit width. The model parameters are m = (M_⋆, M₂, τ, [Fe/H]_i, w). [Fe/H]_s differs from the initial abundance ([Fe/H]_i) due to diffusion and mixing processes throughout stellar evolution. The GARSTEC evolutionary models used here are the same as the ones used in BAGEMASS. GARSTEC uses the FreeEOS⁷ equation of state (Cassisi et al. 2003) and standard mixing length theory for convection (Kippenhahn & Weigert 1990). The mixing length parameter used to calculate the default model grid is α_MLT = 1.78. With this value of α_MLT GARSTEC reproduces the observed properties of the present day Sun assuming that the composition is that given by Grevesse & Sauval (1998), the overall initial solar metallicity is Z_⊙ = 0.01826, and the initial solar helium abundance is Y_⊙ = 0.26646. These are slightly different to the value in Serenelli et al. (2013) because we have included additional mixing below the convective zone in order reduce the effect of gravitational settling and so to better match the properties of metal-poor stars. Due to the effects of microscopic diffusion, the initial solar composition corresponds to an initial iron abundance [Fe/H]_i = +0.06 dex. The stellar model grid covers the mass range 0.6 M_⊙–2.0 M_⊙ in steps of 0.02 M_⊙. The grid of initial metallicity values covers the range [Fe/H]_i = −0.75 dex to −0.05 dex in steps of 0.1 dex and the range [Fe/H]_i = −0.05 to +0.55 in 0.05 dex steps.

To obtain M₂ from f(m), M_⋆ and P, we need to know inclination from the transit light-curve. Degeneracies between i, R_⋆/a and k are such that we choose to fit the full-width half maximum of the transit,

$$\begin{array}{c}\hfill w=\frac{R_{\star }}{a}\frac{\sqrt{1-b^2}}{\pi },\end{array}$$(6)

instead of the inclination. We implement a Gaussian prior on [Fe/H]_s from spectroscopy and used uniform priors for age, M_⋆ and M₂. We ran a burn-in chain of 100 000 draws before drawing an additional 50 000 draws to sample the PPD for M_⋆, M₂ and τ. The number of post-burn-in draws matches that of the orbital fit.

We used an up-to-date constant from IAU resolution B3 (Prša et al. 2016) to calculate a from P, M_⋆ and M₂,

$$\begin{array}{c}\hfill a=4.208278\times P^{\frac{2}{3}}{(M_{\star }+M_2)}^{\frac{1}{3}}·\end{array}$$(7)

This can then be combined with R_⋆/a and k to calculate the PPD for R_⋆ and R₂. We selected the the median value of each parameters PPD as our measurements (see Fig. 9), with uncertainty equal to the largest difference between the median and the 16th and 84th percentile of the cumulative PPD for each parameter from the second chain.

Fig. 9. Posterior probability distributions from EBLMMASS for J2349−32. Over plotted are the 1-σ contours.

4. Results

4.1. EBLM J2349−32

J2349−32 was observed over three consecutive years by the WASP project. In each season we find quasi-periodic signals at periods of 4.42 d, 4.35 d and 4.42 d with amplitudes between 3–4 mmag. Each of these signals has a false-alarm probability < 10⁻⁵ and so we assumed this is detection of the rotational period of the primary star (P_rot = 4.40 ± 0.03 d). From the WASP photometry, we measured Δt_tr = 0.126 d and Δm = 0.043 mag. corresponding to R_⋆/a ≈ 0.11 and k ≈ 0.21.

The best SED fit (${\chi }_{\text{red}}^2=1.24$) corresponds to a star with spectral type F9 with a low reddening (E(B − V) ≤ 0.034 to 1-σ). This system was included in Gaia DR2 (source ID 2314099177602409856). The G magnitude was measured to be 11.413 and the parallax is 3.769 ± 0.092 mas (265.32 ± 6.47 pc). Gaia DR2 shows a single star (G = 15.219) 48″ away at a position angle of 111° (source ID 2314099173307737088). This source is not included in the sky annulus of the SAAO 1-m photometry, but falls within the WASP aperture where it will contribute around 3% of the total flux. The proper motions of this star and J2349−32 are significantly different in right ascension and declination so we concluded that they are not associated.

J2349−32 was observed with lucky imaging on 2017-07-08 where two companion stars were detected. A close companion was found at a separation of 1.402 ± 0.013″ and position angle of 308.6 ° ±0.6° (Fig. 10). We measured the companion to be 5.55 ± 0.08 mag fainter in the TCI red-arm images; the companion was not sufficiently resolved in the TCI visible-arm images to obtain any reliable measurements. A second, distant companion was detected at a separation of 25.70 ± 0.07″, position angle of 218.6 ± 0.3°. We find that it is 9.0 ± 0.3 mag fainter with the TCI red-arm images and 8.5 ± 0.3 mag fainter in the visible-arm images. This is the same source identified by Gaia DR2 (source ID 2314099173307737088). If the closest companion is blended in the CORALIE and SAAO 1-m apertures, we estimate that it only contributes 0.6% of the light and is too faint to significantly modify the transit light-curve.

Fig. 10. Lucky imaging of J2349−32 (red arm) revealing a close companion 1.3″ (blue circle) away at a position angle of 308.6 ± 0.6° (blue circle).

The eightteen CORALIE spectra were combined to produce a spectrum with S/N = 40. The analysis of this spectrum shows that the primary is a slightly metal-deficient star with a temperature consistent with the SED fit. There is a weak Li I line at 670.7 nm from which we measured a lithium abundance log A_Li + 12 = 2.4 ± 0.08. This value was estimated by synthesising a small region around this line in ISPEC using fixed atmospheric parameters from wavelet analysis (Table 3) and manually adjusting the lithium abundance to obtain the best fit by-eye.

The RVs were fitted simultaneously with a single transit in I-band from the SAAO 1-m telescope to obtain the best fitting orbital solution (${\chi }_{\text{red}}^2=0.93$; Fig. 11). The PPD for eccentricity is consistent with a circular orbit (e ≤ 0.05 to 5-σ). We find a negligible drift in systematic velocity (≤15 m s⁻¹ yr⁻¹ to 1-σ). The best-fitting limb-darkening temperature agrees with effective temperatures measured with SED fitting and wavelet analysis to better than 1-σ.

Fig. 11. Orbital fit of J2349−32. Top panel: detrended I-band light-curve from the SAAO 1-m telescope (black) with the best fitting transit model (red). Upper-middle panel: phase-folded WASP lightcurve (black). Lower-middle panel: drift-corrected radial velocity measurements (black) with the best model (red). Bottom panel: residuals from radial velocity model measurements.

EBLMMASS predicts a primary star which has a mass and radius similar to the Sun, but is approximately 350 K hotter. This is partly due to this being a metal-poor star, but also because it is approximately half the age of the Sun. The youthfulness of this star in conjunction with a convection zone which is unable to transport lithium to the core where it would be burnt may explain why lithium is detected with spectroscopy. The secondary component’s mass is consistent with that of an M-dwarf below the fully convective limit.

4.2. EBLM J2308−46

J2308−46 has WASP photometry spanning 5 years. The last season of data had less than 400 data points so was excluded. We measured a strong P/2 signal in two seasons of data. Phase folding WASP photometry at this period reveals a moderate ellipsoidal variation with an amplitude of 5 mmag (Fig. 14). We fixed parameters associated with ellipsoidal variation to produce a good out-of-transit fit to the WASP photometry (q = M₂/M_⋆ = 0.05, gravity darkening coefficient = 0.1) to measure Δt_tr = 0.109 d and Δm = 0.018 mag and estimate the transit parameters R_⋆/a ≈ 0.20 and k ≈ 0.13.

SED fitting measured the effective temperature of the primary star to be consistent with a spectral type F7 (T_eff = 6270 ± 140 K; ${\chi }_{\text{red}}^2=0.77$) with a low reddening (E(B − V)≤0.054 to 1-σ). This system is included in the Gaia DR2 catalogue (Source ID 6539811294185397120; G = 11.361) with a parallax of 2.187 ± 0.113 mas (457.24 ± 23.63 pc). There is a close companion 22.5″ from J2308−46 at a position angle of 282° (G = 15.388; source ID 6539811500344886016). This is clearly resolved in the follow-up 1-m R-band photometry from SAAO and does not contaminate the sky annulus. It does fall within the WASP annulus, contributing ≈3% of the total flux. There is another source 48″ away from J2308−46 at a position angle of 23°. This is also included in Gaia DR2 (G = 14.919; source ID 6539817204061452544) which is on the limits of the WASP aperture and would contribute less flux than the source at position angle 282°. J2308−46 was observed with lucky imaging on 7 July 2017 with only a single, faint companion being found, located 21.38 ± 0.04″ away at a position angle of 208.2 ± 0.4° (Fig. 12). We measured magnitude differences of 8.2 ± 0.3 mag in the TCI red-arm images and 8.0 ± 0.2 mag in the TCI visible-arm images. This object is included in the Gaia DR2 catalogue with Source ID 6539811289890737408 with G = 19.538. We compared the proper motion of J2308−46 with this object and conclude they are not physically associated.

Fig. 12. Lucky imaging of J2308−46 (red arm) revealing a companion 20″ away with a position angle of 208° (blue circle).

A total of twenty two CORALIE spectra were co-added to produce a spectrum with S/N = 20. Wavelet decomposition shows that the primary star is moderately-rotating (V sin i ≈ 39 km s⁻¹) and metal poor ([Fe/H] = −0.15 dex). The primary star’s effective temperature appears to be close to the temperature at which the Kraft break (Kraft 1967) becomes apparent. The Kraft break is an abrupt reduction in surface rotation for stars with effective temperatures below ∼6200 K. This is due to the presence of a efficient magnetic dynamo which transfers angular momentum from the star through stellar winds resulting in magnetic breaking.

Fitting the follow-up photometry jointly with radial velocity measurements was non-trivial as clear systematic errors remained in the SAAO 1-m light-curve after initial detrending. We obtained an orbital solution in the same framework as EBLM J2349−32 but found an unacceptable fit around contact point 2 and the continuum prior to contact point 1 in the SAAO 1-m light-curve (see top panel of Fig. 14). We attempted further detrending of the follow-up light-curve with airmass, CCD position and time but this did not successfully remove the problem. Instead, we decided to generate a red-noise model using Gaussian processes. We used the CELERITE package (Foreman-Mackey et al. 2017) and the following kernel with the default value of ϵ = 0.01 to approximate the Matérn-3/2 covariance function:

$$\begin{array}{c}\hfill k({\tau }_k)={\sigma }^2\left[(1+1/ϵ)e^{-(1-ϵ)\sqrt{3}{\tau }_k/\rho }(1-1/ϵ)e^{-(1+ϵ)\sqrt{3}{\tau }_k/\rho }\right].\end{array}$$

Here, τ_k is the time difference between two observations, ρ is a parameter that controls the time scale over which observational errors are correlated and σ controls the amplitude of such variations. The free parameters, logρ and logσ, tended to a value that over-fitted the noise in the light-curve if it remained as a free parameter in the joint fit. Instead, we adjusted these values by-eye until we find an acceptable red-noise model that accounted for the data around the second contact point (logρ = 2 and logσ = 2). The parameters were then fixed at these values to find an acceptable orbital solution (${\chi }_{\text{red}}^2=1.32$) in the same way as J2349−32.

The primary star is close the the “blue-hook” phase of its post main-sequence evolution (Fig. 13). This results in two peaks in the PPDs for M_⋆, M₂ and τ which are consistent to within 2-σ. Both solutions could be valid and so it is a requirement to fit these systems to assess the likelihood and validity of each. We fitted double-Gaussian models to the PPDs of M_⋆, M₂, τ and a which have been sorted into 100 equal bins. We used the Levenberg-Marquardt algorithm to find the optimal model vector m = (A₁, μ₁, σ₁, A₂, μ₂, σ₂) for the double Gaussian model:

$$\begin{array}{c}\hfill y=A_1{e}^{-\frac{{(x-{\mu }_1)}^2}{2{\sigma }_1^2}}+A_2{e}^{-\frac{{(x-{\mu }_2)}^2}{2{\sigma }_2^2}},\end{array}$$(8)

Fig. 13. Posterior probability distributions for the density and temperature of the primary star in J2308−46 is shown in the top panel. The zero-age main sequence is show(black-dashed) along with the best fitting isochrone (blue-solid) and the respective isochrones for ±1-σ in [Fe/H]. The lower panels show the PPD distributions for M1, M2 and τ with best-fitting double-Gaussian models in red.

Fig. 14. Orbital fit of J2308−46. Upper panel: R-band transit obtained from the SAAO 1-m telescope (black) with the best fitting transit model (green dashed). We plot the best fitting transit model generated using Gaussian processes (red). Middle panel: phase-folded WASP observations (black) and observations binned into groups of 50 (blue). We also plot the Roche model used to approximate the out-of-transit photometry used to measure transit parameters from WASP photometry (red). Lower panel: drift-corrected radial velocity measurements (black) with the best fitting model (red) and residuals from the best fitting orbital model.

where x is the position of the bin, y is the number of models in the respective bin, μ is the measurement of the model, σ is the uncertainty associated with the model, and A represents the number of models at the peak of the of the distribution. The resulting fit for J2308−46 isn’t entirely satisfactory; the fitted values of μ do not entirely match up with the peaks of the PPD for M_⋆, M₂ and τ. This is partly due to the PPDs being poorly described by a Gaussian. Other EBLMs (e.g. J0218−31) have double-peaked PPDs which are well described by Gaussian, so we decided to add additional uncertainty rather than seeking a more complex model. To account for this, we added an additional uncertainty of 2% for M_⋆, M₂ and τ which was estimated by measuring the offset between the fitted values of μ and the peaks of the respective PPDs. Moreover, the width of each PPD (σ) is underestimated upon visual inspection leading to an additional 1% uncertainty which was determined by-eye. The total additional uncertainty for each σ is 3%. We assessed each solution using the ratio of likelihoods. We find that the younger solution (τ = 3.98 ± 0.86 Gyr) is preferred over the older solution (τ = 5.81 ± 1.0 Gyr) with a factor ℒ(3.98 Gyr)/ℒ(5.81 Gyr)≈3.07. This is moderate evidence to favour the younger solution but far from conclusive so we report both solutions in Table 3.

4.3. EBLM J0218−31

J0218−31 was observed over three years by the WASP survey. We find a tentative detection of spot-induced variation across the three seasons (P_rot = 2.30 d, 2.14 d and 2.60 d). Each have a false alarm probability < 10⁻⁵ and an amplitude around 1 mmag amongst a complex periodogram of similar (but smaller) amplitudes making it unclear whether this is due to spot-induced variations (P_rot = 2.35 ± 0.20 d) or poor-quality photometry. From the WASP photometry, we estimated Δt_tr = 0.241 d and Δm = 0.03 mag, corresponding to R₁/a ≈ 0.09 and k ≈ 0.18.

We obtained a good SED fit (${\chi }_{\text{red}}^2=0.75$) with the effective surface temperature of the primary star consistent with a spectral type F9 (T_eff ≈ 6020 K). J0218−31 is included in the second data release of Gaia (G = 9.734; Source ID 4971670729566470528) with a parallax measurement of 3.762 ± 0.092 mas (265.82 ± 6.50 pc). There are 3 close and faint companions within 22″ at position angles 266°, 332° and 92°. The brightest at position angle 332° has G = 17.140 (ΔG = 7.406) which would have a negligible flux contribution to the aperture of the WASP photometry and the follow-up R-band photometry. The neighbours at positional angles 266° (G = 18.128) and 92° (G = 19.813) are fainter still. A brighter companion (G = 16.015; source ID 4971670935725243904) is located 50″ away at a position angle of 330°. This does not overlap the sky annulus of the 1-m SAAO photometry and will have a negligible flux contribution to the WASP photometry. The proper motions of these stars are not similar to J0218−31 and so we concluded that they are not physically associated.

We co-added fifty-three out-of-transit spectra to produce a spectrum with S/N = 30. Using wavelet decomposition, we estimated T_eff = 6100 ± 85 K confirming a spectral class of F9 from SED fitting. The effective temperature is 1-σ hotter than predicted by SED fitting suggesting there could be some additional reddening that is unaccounted for. The iron content is higher than the Sun ([Fe/H] = 0.15 ± 0.12 dex). There is also a strong Li I line in the spectrum from which we measured log A_Li + 12 = 3.24 ± 0.08 suggesting that the convective shell of J0218−31 may be similar to that of J2349−32.

We fitted the Rossiter-McLaughlin measurements alongside the out-of-transit radial velocity measurements with g′, r′, i & z′ band photometry to obtain the best fitting orbital solution (${\chi }_{\text{red}}^2=1.68$; Fig. 15). We initially fitted an independent value of k to the g′, r′, i & z′ follow-up photometery. The fitted value of k for each bandpass agreed with each other to 1-σ suggesting there is no wavelength-dependent transit depths which may have indicated a source of third light. However, we do find a significant drift in systematic velocity (d(γ)/dt = −69.9 ± 4.1 m s⁻¹ yr⁻¹) which suggests there may be a faint third body in the system. With the addition of R-M measurements, we were able to calculate the sky-projected angle between the rotational and orbital axes, λ = 4 ± 7°, which is consistent with the assumption that these axes are aligned. From this we also measured V sin i = 10.28 ± 2.12 km s⁻¹ which is in agreement with the value inferred from wavelet decomposition.

Fig. 15. Orbital fit for J0218−31. Top panel: transit photometry from CTIO in g′ (blue), r′ (red), i′ (cyan) and z′ (green) with best fitting models shown in black. Upper-middle panel: phase-folded WASP lightcurve. Lower-middle panel: drift-corrected radial velocity measurements from CORALIE with best fitting model plotted in red, along with residuals. Bottom panel: drift-corrected radial velocity measurements during transit (the Rossiter–McLaughlin effect; black) with the best fitting model (red). Error bars have been omitted for clarity.

Similarly to J2308−46, J0218−31 has entered the “blue-hook” part of it’s post main-sequence evolution resulting in double-peaked PPDs of τ, M_⋆ and M₂ (see Fig 16). We used the same approach for J2308−46 to fit a double-Gaussian to the PPDs for τ, M_⋆ and M₂ and found that the younger solution (2.4 ± 0.25 Gyr, M_⋆ = 1.55 ± 0.05 M_⊙, R_⋆ = 2.13 ± 0.09 R_⊙) is favoured with almost twice the likelihood ℒ(2.35 Gyr)/ℒ(3.80 Gyr)≈3.55 of the older solution. This is moderate evidence to suggest the younger solution is favoured but we report both solutions in Table 3 as a precaution.

Fig. 16. Posterior probability distributions for the density and temperature of the primary star in J0218−31 is shown in the top panel. The zero-age main sequence is show(black-dashed) along with the best fitting isochrone (blue-solid) and the respective isochrones for ±1-σ in [Fe/H]. Lower panels: PPD distributions for M1, M2 and τ with best-fitting double-Gaussian models in red.

4.4. EBLM J1847+39

J1847+39 was observed for three years with the WASP survey. For each season, we find significant spot-induced variations at periods 7.55 d, 7.14 d and 7.17 d with amplitudes of 3–4 mmag; each of these signals has a false alarm probability < 10⁻⁵ and we assumed this to be a detection of rotational spot modulation at a period of 7.29 ± 0.19 d. We find no evidence for ellipsoidal variation in the WASP lightcurve from which we estimated Δt_tr = 0.02 d and Δm = 0.10 mag corresponding to initial estimates of R_⋆/a ≈ 0.07 and k ≈ 0.32.

The best SED fit estimated the primary star to be of spectral type F9 with temperature of 6020 ± 100 K (${\chi }_{\text{red}}^2=1.37$). This system was included in Gaia DR2 (G = 11.677; source ID 2098283457595740288) with a parallax of 3.583 ± 0.086 mas (279.10 ± 6.70 pc). The field surrounding J1847+39 is relatively more crowded compared to the other targets, with over 7 targets brighter than G = 17 within 1.5′. The closest companion is 12″ away (G = 17.694; source ID 2098283457595821440) at a position angle of 196°. This neighbour would not have been included in the aperture of J1847+39 (radius ∼6″). A magnitude difference of ΔG = 6.17 results in less than 0.3% flux contribution if it was included in the apertures of the HAO photometry. There are two bright companions 1.28′ and 1.07′ away at position angles of 48° (G = 11.577) and 54° (G = 12.135) respectively. These are beyond the sky annulus of WASP but may still contribute a small amount to the total flux. The proper motions of these objects are dissimilar to J1847+39 and so we concluded that they are not associated.

Ten INT spectra were co-added to produce a spectrum with S/N = 30. We used the spectral synthesis method on the wings of the Hα line to estimate T_eff = 6200 ± 100 K (spectral type F8) which is consistent with the SED fit. We were able to fit 11 un-blended Fe I lines in the region around Hα from which we measured [Fe/H] for each line. We took the mean of value of [Fe/H] as the iron abundance measurement with the standard deviation as the uncertainty ([Fe/H] = −0.25 ± 0.21 dex). We were unable to determine log g due to the limited wavelength coverage of the H1800V grating so we assumed log g = 4.44 for the aforementioned synthesis and interpolation of limb-darkening coefficients.

Radial velocity measurements were fitted simultaneously with single transits in CBB, g′ and z′ filters to obtain the best fitting orbital solution (${\chi }_{\text{red}}^2=1.77$; Fig. 17). J1847+39 has the most eccentric orbit of the sample (e = 0.209 ± 0.014). We attempted to fit an independent value of k for photometry in each filter and found them all to agree within 1-σ suggesting there is no significant third-light contamination. However, we do measure d(γ)/dt = −71.9 ± 21.7 km s⁻¹ yr⁻¹ suggesting that there may be a faint third-body in the system. The best-fitting limb-darkening temperature, T_eff, ld, is ∼600 K hotter than spectroscopic and photometric analysis; the reason for this is unclear.

Fig. 17. Orbital fit of J1847+39. Top panel: single transits from the HAO in filters CBB (blue), g′ (cyan) and z′ (green) with best fitting models (red). Upper-middle panel: phase-folded WASP lightcurve. Lower-middle panel: drift-corrected radial velocity measurements (black) with the best fitting model (red) and residuals (lower panel).

The best fitting solution from EBLMMASS describes a star similar to the Sun in mass and size, but a fifth of it’s age (τ = 1.10 ± 1.80 Gyr). The systems eccentricity may be primordial in origin as there would have been insufficient time for tidal interaction to circularise the orbit. The M-dwarf’s mass is in the convective transition (∼0.35 M_⊙) and provides an interesting test for low-mass stellar models in a region that is highly debated.

4.5. EBLM J1436−13

J1436−13 was observed over 3 consecutive seasons with the WASP survey. We find significant variability for each season at periods of 3.99 d, 3.98 d and 4.02 d with amplitudes between 4 and 5 mmag. Each have a false alarm probability < 10⁻⁵ and we assumed this is a detection of spot modulation corresponding to a rotational period P_rot = 4.00 ± 0.02 d. We find no evidence for ellipsoidal variation in the WASP photometry from which we measured Δt_tr = 0.188 d and Δm = 0.065 mag corresponding to initial estimates of R_⋆/a ≈ 0.148 and k ≈ 0.256.

The SED fitting measured the primary star to have a spectral type of F9 (T_eff = 6080 ± 355 K; ${\chi }_{\text{red}}^2=0.87$) with little reddining (E(B − V)≤0.055 to 1-σ). This system is included in Gaia DR2 (G = 12.334; source I.D 6323183619200685824) with a parallax of 2.063 ± 0.097 mas (484.73 ± 22.79 pc). There is a faint (G = 18.994) background star included in Gaia DR2 that is 17″ away at a position angle of 302°. This is included in the WASP aperture but contributes less than 0.1% of the total flux.

Thirteen CORALIE spectra were co-added to produce a spectrum with S/N = 30. Wavelet decomposition measured a value of T_eff that is around 300 K hotter than predicted by SED fitting suggesting that there may be some unaccounted reddening. The iron content is slightly less than the Sun ([Fe/H] = −0.10 ± 0.12 dex) and the magnesium lines are relatively narrow suggesting a low surface gravity. We were unable to identify any measurable lithium lines.

The best fitting orbital solution (${\chi }_{\text{red}}^2=1.75$) describes a transit with a high impact parameter (b = 0.86 ± 0.07; Fig. 18). The limb-darkening temperature agrees better with SED fitting than wavelet decomposition, but is consistent with both to 1-σ. Radial velocity measurements suggest the system is circularised (e ≤ 0.004 to 1-σ) and there is no significant drift in systematic velocity. J1436−13 is slightly larger and more massive than the Sun. The uncertainty in R_⋆ (5%) is the largest in the sample owing to poorly constrained values of R_⋆/a and k owing to a high impact parameter. The M-dwarf companion is the most massive of the sample (M₂ = 0.49 M_⊙).

Fig. 18. Orbital fit of J1436−13. Top panel: a single transit obtained from SAAO in R filter (black) and the best fitting transit model (red). Upper-middle panel: phase-folded WASP lightcurve. Lower-middle panel: drift-corrected radial velocity measurements (black) with best fitting model (red) along with residuals (lower panel).

5. The mass-radius diagram

In Fig. 19, we plot the 5 Gyr isochrones for [Fe/H] = 0 dex (B15; Baraffe et al. 2015) and [Fe/H] = −0.5 dex (B98; Baraffe et al. 1998) and compare them to the five EBLMs measured in this work. The B15 isochrones rectify some of the flaws in the models presented by Baraffe et al. (1998) (e.g. optical colours that are too blue). Visual inspection of the radii shows that they are broadly consistent with evolutionary models. The EBLMs in this work mostly have a sub-solar metalicity and are expected to have radii between or about the B98 and B15 models.

Fig. 19. Masses and radii of M-dwarfs in five EBLM systems. The 5 Gyr isochrone for solar metallicity from Baraffe et al. (2015) is plotted (black-solid) along with the 5 Gyr isochrone for [M/H] = −0.5 (Baraffe et al. 1998; black-dashed). We plot low-mass M-dwarfs with masses and radii known to better than 10% (from Table 4 of Chaturvedi et al. 2018, and references therein). For J2308−46 and J0218−31 we plot both solutions and label accordingly. We also plot TRAPPIST-1 (Delrez et al. 2018), Proxima Centauri (Anglada-Escudé et al. 2016) and J0555−57 (von Boetticher et al. 2017).

Two EBLMs (J0218−31 and J1436−13) have high impact parameters leading to a larger uncertainty in R_⋆/a, k and ultimately R₁ and R₂. The effect is most significant in J1436−13 where the uncertainties in R₂ span across both B98 and B15 isochrones. The primary stars of two EBLMs (J2308−46, J0218−31) have evolved into the “blue-hook” part of their post main-sequence evolution, leading to two solutions of M_⋆, M₂ and τ. Although a single solution is marginally favoured for each, both are valid and we report both in Table 3 and Fig. 19 as a precaution.

6. Bayesian measurements of radius inflation

The traditional approach of interpolating between solar B98 ([Fe/H] = 0) and B15 ([Fe/H] = −0.5) isochrones of fixed age is not sufficient to assess inflation, especially for young systems below 1 Gyr which may still be contracting. A recent and well-sampled set of isochrones for low-mass stars are required to assess if the M-dwarf in each EBLM system is consistent with the isochrone for the respective measurement of [Fe/H] and τ. For this task, we used MESA stellar evolution models. The MESA models are created using the protosolar abundances recommended by Asplund et al. (2009) as the reference scale for all metallicities; this is consistent with the grid of spectra from wavelet analysis (Gill et al. 2018). MESA uses the OPAL equation of state tables from Rogers & Nayfonov (2002) along with opacity tables from Freedman et al. (2008), Yurchenko et al. (2011) and Frommhold et al. (2010). MESA also includes complex treatments for microscopic diffusion and gravitational settling (both important for low-mass stars), radiative levitation (important for high-mass stars), rotation, convective overshooting, magnetic fields and mass-loss.

The web interpolater⁸ was used to create a grid of MESA isochrones spanning the range [Fe/H] = −2 to +0.5 dex in steps of 0.5 dex and age range 0.8−9 Gyr in steps of 0.2 Gyr. Using this grid, we created a bi-linear interpolation routine (in dimensions of τ and [Fe/H]) to obtain an expected radius, R_2, exp for a given mass. To assess inflation, the following procedure was employed for each draw in each Markov chain from EBLMMASS and the orbital solution:

1. log g₂ can be calculated from the orbital solution,

   $$\begin{array}{c}\hfill g_2=\frac{2\pi }{P}\frac{\sqrt{1-e^2}{K}_{\star }}{{\left(\frac{R_2}{a}\right)}^{2}sini}·\end{array}$$(9)

2. The corresponding draw for M₂ can be combined with log g₂ to obtain the calculated value of R₂,

   $$\begin{array}{c}\hfill R_2=\sqrt{\frac{G{M}_2}{g_2}}·\end{array}$$(10)

3. The corresponding draw for τ was used with a random value for [Fe/H] to interpolate a MESA isochrone. The random value of [Fe/H] was drawn from a Gaussian distribution of mean and width corresponding to the measurement of [Fe/H] and uncertainty of [Fe/H] reported in Table 3.

4. The corresponding draw of M₂ was used to interpolate an expected radius for the M-dwarf companion, R_2, exp.

5. R_2, exp and R₂ can be combined to calculate the fractional radius residual,

   $$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{R}_2}{R_2}=\frac{R_2-R_{2,\exp }}{R_2}·\end{array}$$(11)

By repeating the above procedure for each draw in each Markov chain from EBLMASS and the orbital fit, we were able to estimate the PPD for the fractional radius residual for each EBLM (Fig. 20). Four EBLMs have narrow-peaked PPDs for ΔR₂/R₂ (top panel of Fig. 20). For these, we calculated the nominal fractional radius by binning the PPD into 100 bins and fitted a Gaussian model; we took the mean of the fitted Gaussian to be the measurement of ΔR₂/R₂ with uncertainty equal to the standard deviation. We find that a Gaussian shape is not a perfect fit to the PPDs of ΔR₂/R₂; there are asymmetric discrepancies where one side of the Gaussian model is lower than the PPD, whilst the other is too high. On average, the under-prediction on one side and over prediction on the other are of the same magnitude and we assume the widths still accurately represent the mean uncertainty of ΔR₂/R₂.

Fig. 20. Fractional radius residual PPD for J2349−32, J2308−46, J0218−31 and J1847+39 system along with Gaussian models (upper panel). Lower panel: J1436−13 which has a broader PPD and requires a double Gaussian model. The legends denote the best-fitting value of the radius residual with uncertainty.

J1436−13 has a significantly higher impact parameter which broadens the PPD for R₂ and thus, ΔR₂/R₂ (lower panel of Fig. 20). We approximate this shape with a double Gaussian, and used an identical routine used to measure the double-peaked PPDs in Sect. 4.2. The fit for J1436−13 is not perfect and we used the peak of the PPD as the measurement of ΔR₂/R₂ with uncertainty equal to the standard deviations of each fitted Gaussian added in quadrature.

J2308−46 and J1847+39 appear deflated by at least 1-σ compared to MESA evolutionary models. There is also some evidence to suggest J0218−31 is deflated although the measured radius is consistent with predictions from evolutionary models to 1-σ. Conversely, J2349−32 appears inflated by 1-σ. The PPD for ΔR₂/R₂ for J1436−13 is broadly consistent with evolutionary models although it is not possible to assess inflation for this since the radius is so poorly constrained. J2308−46 and J0218−31 have double-peaked distributions for M₂ and τ and we expected the PPDs for ΔR₂/R₂ to be shaped similar since M₂ is used to calculate R₂, and combined with τ to estimate R_2, exp. In creating the PPD for R₂ (Eq. (10)), the division of the PPD for M₂ with the PPD for g₂ diminishes the double-peaked nature observed in the PPD M₂, leading to a Cauchy-like PPD for R₂. The interpolated value of R_2, exp is dependent on τ and M₂ which are both double peaked. R_2, exp is not expected to have a double-peaked PPD as each combination of τ and M₂ was a trial step in EBLMMASS and will correspond to a similar expected radii (i.e. higher values of τ will correspond to lower values of M₂ and vice-versa). Thus the PPD for ΔR₂/R₂ is single peaked with width controlled by the uncertainty in M₂, g₂ and [Fe/H].

7. Systematic effects on determining mass, radius and age

One major issue remains with the method employed in this paper and previous publications of the EBLM project: we are attempting to test evolutionary models of low-mass stars using the models of better-understood F-dwarfs. This method is acceptable when the uncertainty propagated by stellar models for FG stars are much smaller than the uncertainties of physical properties propagated by observational uncertainties in radial velocity and transit photometry measurements. We must also consider the possibility of unresolved, faint components in the aperture which could systematically modify the transit shape along with physical properties inferred from modelling. In the following sections we explore the scale of some potential sources of additional uncertainties on measurements of physical properties for EBLM systems. For J2308−46 and J0218−31, we only considered the most probable solution in this section.

7.1. Evolution ambiguity, α_MLT and Y_He

The default model grid used in EBLMMASS uses a mixing length parameter α_MLT = 1.78 and an initial helium abundance Y = 0.26646 + 0.984 Z, both of which have been calibrated on the Sun. As noted by Maxted et al. (2015a), these assumptions are subject to some level of uncertainty. Maxted et al. (2015b) estimated the additional uncertainty in M_⋆ and τ for 28 transiting exoplanet host stars by assuming an error of 0.2 in α_MLT and 0.02 for ΔY. They find that systematic errors in M_⋆ and τ from Y and α_MLT can be comparable to the random errors in these values for typical observational uncertainties in the input parameters. We note that the sample used in Maxted et al. (2015b) consists primarily of stars less massive than the Sun, whereas the primary stars in this work are more massive F-type stars. Three grids of models are provided with EBLMMASS: 1. α_MLT = 1.78, ΔY = 0.00, 2. α_MLT = 1.5, ΔY = 0.00 and 3. α_MLT = 1.78, ΔY = 0.02; we used grid 1 in Table 3. We re-measured the mass, radius and age of both components with the grids 2 and 3 to see how the uncertainties in α_MLT and ΔY impact our results. We used grid 2 to assess an additional uncertainty of 0.28 for α_MLT and grid 3 to assess an additional uncertainty of 0.02 for ΔY. We used the same orbital solution and atmospheric parameters reported in Sect. 4 and report our results in Table 4.

We find that an additional uncertainty in α_MLT results in up to a 4% increase in the mass uncertainty. This is mirrored by the 1−5% increase in the mass uncertainty seen for an additional uncertainty of 0.02 for ΔY. These results are consistent with those found by Maxted et al. (2015b) and the largest increase in uncertainties are seen for the host stars. The measured age of each system has a larger fractional error than the measured masses and we can expect and additional uncertainty if α_MLT and ΔY are unconstrained. Assuming that both parameters are poorly constrained, we can expect an additional 3−5% uncertainty in the mass of the M-dwarf companion. This is significant and has potential to skew the interpretations of the mass-radius diagram.

A further limitation arises when the primary star has evolved into the post-main sequence blue hook (Henyey hook). J2308−46 and J0218−31 are in this region leading to two distinct solutions for M_⋆ and τ. A single solution is preferred for both these systems but there will always be some ambiguity until further mass constraints can be obtained. A solution may lie in the increased contrast between a FGK star and an M-dwarf in the infrared. It may be possible to detect molecular lines (VO, TiO and CaH) associated with an M-dwarf using high-resolution infrared spectroscopy. Simultaneously cross-correlating templates from G-/M-dwarfs with near-infrared spectra of EBLMs should produce similar results to what Bender et al. (2012) achieved in the optical for Kepler-16. The act of turning an SB1 into and SB2 would provide a direct test of the methods used in this work. This would also place a further constraint on which mass and age solution best describes the system in cases whereby the primary star is close to the post main-sequence “blue-hook”. Inspection of where these targets sit in the colour-magnitude diagram (Fig. 1) reveal that J2308−46 and J0218−31 are significantly closer to the giant branch than the other three EBLMs. An arbitrary cut using Gaia colours could be used to pre-select cooler host stars (≤6100 K) and avoid host stars near the post-main sequence blue hook.

7.2. Third light effect

Lucky imaging provides constraints on nearby contaminating objects. For J2349−32 and J2308−46, We find close companions which do not significantly contaminate follow-up photometry. For J0218−31 and J1847+39 we can put constraints on the amount of third light from the consistency between the ratio of the radii measured from transit photometry in different pass bands. For J1436−13 we have to rely on existing surveys to identify any nearby stars which may contaminate follow-up photometry. Inspection of the Gaia survey DR2 (resolution of ≤1″) finds no evidence of blends or contamination. Ground-based lucky imaging has a limited resolution to resolve companions with a sky-projected separation of ∼0.3″⁹. The orbital separation for each EBLM corresponding to a sky-projected separation of 0.3″ was calculated using parallax measurements from Gaia DR2 (Table 5). The period of such orbits were also calculated using measurements of M_⋆ from Table 3. We find that the closest EBLM (J2349−32 at a distance of 259 ± 3 pc) would require a semi-major axis of at least 389 au with orbital period spanning decades. The three-body systems identified by Triaud et al. (2017) will have orbital periods similar of the order of decades and would be difficult or impossible to resolve through lucky imaging.

The spectrum itself can provide useful insights for potential aperture contamination. The analysis of CORALIE spectra for 118 EBLM systems presented by Triaud et al. (2017) found that 17.8% of these systems show significant evident for non-zero values of d(γ)/dt (spanning d(γ)/dt = 0.07 − 4.5 km s⁻¹ yr⁻¹). J0218−31 and J1847+39 have best-fitting values of d(γ)/dt which are at the bottom of this bracket. If these drifts are evidence of a third body, they would have separations which are challenging to resolve with lucky imaging. If they could be resolved with lucky imaging, they would have periods which would require decades of radial velocity monitoring to characterise them (Table 5). The longest EBLM system has spectroscopic observations spanning 7.5 yr (see Table 1) and so determining the nature of the systematic drift in radial velocities for each system would be difficult. The low signal-to-noise spectra from CORALIE and INT eliminates unresolved blends which contribute more than 30% of the luminosity of the primary star by inspection of cross-correlation functions (350−750 nm).

Including third light as a free parameter in the orbital fit changes the shape and depth of a light-curve leading to degeneracies between R_⋆, k and b. We assessed this by re-fitting the orbital solution for all stars assuming a 10% light contamination from a third body which does not interact with the EBLM system. From this fit, we combined best fitting values of R_⋆/a, b, and k and their uncertainties with nominal values of the remaining parameters from the original fit to re-determine R_⋆ and R₂ (first two columns in Table 4). On average, we find a 3−8% increase in R₂ when third light is fixed to 10%; with the largest uncertainty for the smallest M-dwarfs. We ignore J1436−13 from this discussion due to the grazing transit nature. This is comparable to the reported radius inflation for low-mass stars typically quoted in the literature (e.g. 3 − 5%; Spada et al. 2013). However, if we were to see radius inflation in general for the M-dwarf components of EBLM systems then the third-light effect can only be a partial explanation. This is because the majority of these systems do not have detected third bodies in the system, and the third body will often contribute much less than 10% of the total flux in such triple-star systems.

7.3. Limb darkening

To determine accurate estimates for R_⋆, k and b we required an accurate prescription for limb-darkening in our light curve model. For this work, we have used the Claret 4-parameter law (Claret 2000),

$$\begin{array}{c}\hfill \frac{I_{\mu }}{I_0}=1-{\displaystyle \sum _{i=1}^4{a}_i(1-{\mu }_i^{\frac{1}{2}}),}\end{array}$$(12)

where a_i is the ith limb-darkening coefficient. The coefficient tables we used are provided by Claret & Bloemen (2011) for Kepler, Strömgren, Johnson-Cousins and Sloan pass bands based on ATLAS stellar atmosphere models assuming a micro-turbulent velocity ξ = 2 km s⁻¹. For the SDSS pass bands (u′, g′, r′, i′, z′), we used the tables from Claret (2004) provided with ELLC. We interpolate these tables for a given T_eff, [Fe/H] and log g to obtain 4-parameter limb darkening coefficients and a gravity darkening coefficient using the interpolation routine provided with ELLC. As described in Sect. 3.5, we allow the limb-darkening temperature, T_eff, ld, to vary as a free parameter with a Gaussian prior from spectroscopy, and fix log g and [Fe/H] to values from spectroscopic analysis. An alternative is to use the quadratic limb-darkening law (Kopal 1950) with only 2 parameters,

$$\begin{array}{c}\hfill \frac{I_{\mu }}{I_0}=1-{\displaystyle \sum _{i=1}^2{a}_i{(1-\mu )}^i,}\end{array}$$(13)

and allow both coefficients to vary in a fit using the de-correlated parameters a₊ = a₁ + a₂ and a₋ = a₁ − a₂ (Brown et al. 2001). Alternate de-correlation parameters which are not used in this work have been suggested by Kipping (2013): q₁ = (a₁ + a₂)² and q₂ = 0.5 × a₁(a₁ + a₂)⁻¹.

We assessed the choice of limb-darkening law on R_⋆ and R₂ by re-fitting each system using the quadratic limb darkening law (Eq. (13)). We generate coefficients a₁ and a₂ for each pass-band using the Python package LDTK (see Table 6; Parviainen & Aigrain 2015). LDTK uses uncertainties from T_eff, [Fe/H] and log g to estimate uncertainties in the calculated values of a₁ and a₂ (σ_a1 and σ_a2). We used these uncertainties to apply Gaussian priors to a₁ and a₂ and stop the sampler tending to unrealistic values. These priors have a mean value and uncertainty calculated from LDTK. Errors on a₁ and a₂ from errors on T_eff, etc. are very small and unlikely to reflect real uncertainty due to uncertainties in the models so we add a subjective value of 0.05 in quadrature to the uncertainties on each parameter to allow for this. A new combined orbit and light curve solution was found using the same number of draws used in Sect. 4. From this solution, we used R_⋆/a, k and b with their uncertainties and combine it with nominal parameters from the orbital solution in Sect. 4 to measure the radii of components in each system. This ensures that only parameters relating to the radii of the stars were changed.

We find that the additional uncertainty introduced by the choice of limb-darkening law (Table 6) is less than that introduced by third light. The primary and secondary stars see a reduction in R_⋆ and R₂ between 0.5 − 2%. Csizmadia et al. (2013) from their study of exoplanet-host stars conclude that fixing the limb-darkening coefficients to theoretical values does not allow the determination of R₂ to better than 1–10%; a reason why we fitted a₁ & a₂. Intertwined in this is the effects caused by stellar activity, spots and faculae. These are time-dependent effects which change at each transit event and can modify the limb-darkening values far from what is predicted. One conclusion from Csizmadia et al. (2013) is that a star with 0.5% spot coverage can still introduce a 1% uncertainty on k.

8. Conclusion

We present the orbital solutions for five F+M binary systems (EBLMs) discovered by the WASP survey. The host stars for these EBLMs are of spectral type F and are predominantly metal deficient with the exception of J0218−31. J2308−46 has V sin i = 39.82 ± 1.35 km s⁻¹ and appears to be near the Kraft break which separates stars with deep convective envelopes and efficient dynamos to those without.

We find variations in the WASP lightcurves of J2349−32, J0218−31, J1847+39 and J1436−13 which are similar to the best-fitting orbital period. There is a strong p/2 signal for J2308−46 corresponding to ellipsoidal variation. We fixed parameters associated with ellipsoidal variation (q = 0.05, gravity darkening coefficient = 0.1) to estimate the starting parameters of R_⋆/a and k for the orbital solution. We have showed that assuming a spherical star-shape for J2308−46 is a good approximation and does not significantly alter the fitted transit model.

Radial velocity measurements and follow-up photometry were fitted simultaneously to obtain the best-fitting orbital solution (Table 3). All EBLMs except J1847+39 (e = 0.209) have small eccentricities (e ≤ 0.03). J0218−31 and J1847+39 have d(γ)/dt = −69.9 ± 4.1 m s⁻¹ yr⁻¹ and −71.9 ± 21.7 m s⁻¹ yr⁻¹ which suggests that these systems may be influenced by a faint, unresolved companion. Significant trends in the follow-up photometry of J2308−46 required a red-noise model to be generated using a Matern-3/2 kernel with fixed values of logρ = 2 and logσ = 2. This produced a more acceptable fit around contact points 1 and 2 than would have been achievable without a red-noise model. J1847+39 is moderately eccentric (e ≈ 0.2) and is not circularised.

The best-fitting orbital solution was combined with atmospheric parameters to interpolate evolutionary models (EBLMMASS) and measure the masses, radii and age of each component in all five systems. The masses of the primary stars span 0.99–1.55 M_⊙ with radii spanning 0.96–2.13 R_⊙. The primary stars of J2308−46 and J0218−31 have evolved into a region near the post-main sequence blue hook resulting in two distinct solutions of masses, radii and age. For both EBLMs, one solution is slightly more favourable than the other but we report both solutions in Table 3 as a precaution. The M-dwarf companions of J2308−46 and J1847+39 appear deflated by 2-σ and 1.5-σ respectively. There is moderate evidence to suggest that J0218−31 is also deflated despite being consistent with MESA stellar models to within 1-σ. J2349−32 is inflated by 1-σ. J1436−13 has a high impact parameter and it is difficult assess inflation for this target.

We made various choices to measure the masses, radii and ages of both components for each EBLM. One choice was to use the Claret 4-parameter limb-darkening law instead of the plethora of other laws used in the literature. We re-fit all systems using the quadratic limb-darkening law with theoretical coefficients calculated using LDTK. The majority of stars see a reduction of R_⋆ and R₂ below 2%. Spectroscopy used in this work can only rule out unresolved companions which contribute > 30% of the total flux. We investigated the effect of 10% third light when measuring the radii of components in EBLMS. We find that R₂ increases by 3−8% reaffirming the necessity to rule out sources of third-light. The assumptions we made regarding which evolutionary models we used have significant consequences too. Additional uncertainty in ΔY can introduce an additional mass uncertainty ≈1 − 5%, while additional uncertainties for α_MLT introduce an additional mass uncertainty up to 4%. Assuming that both parameters are poorly constrained, we can expect an additional 3−5% uncertainty in the mass of the M-dwarf companion.

Acknowledgments

SG, PM and JE gratefully acknowledge financial support from the Science and Technology Facilities Council [ST/N504348/1, ST/M001040/1, ST/N504348/1 and ST/P000495/1]. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/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. Results from an initial assessment of the data for these targets conducted by the following undergraduate astrophysics students at Keele University were used as a starting point for the analysis presented here: J. Burrows, M. Hill, C. Huggins, M. Jacobs, T. Le Vien, G. Phillips, L. Smith, A. Stay, A. Welch. PL was supported by the project ANT1656. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement n° 803193/BEBOP). MG is FNRS Senior Research Associate, and acknowledges support from an ARC grant for Concerted Research Actions, financed by the Wallonia-Brussels Federation. This paper uses observations made at the South African Astronomical Observatory (SAAO). We thank the reviewer for his/her thorough review and highly appreciate the comments and suggestions which significantly contributed to improving the quality of this work.

References

Appendix A: SED fits

Fig. A.1. Posterior probability distribution of EBLM J2308−46 (top-left), J0218−31 (top-right), J1847+39 (bottom-left) and J1439−13 (bottom-right) from photometric fitting. Over plotted are the 1-σ contours.

Appendix B: Lomb-Scargle diagrams

Fig. B.1. Lomb-scargle diagrams for each season of WASP photometry for J2349−32 (top), J2308−46 (middle) and J0218−31 (bottom).

Fig. B.2. Lomb-scargle diagrams for each season of WASP photometry for J1847+39 (top) and J1436−13 (bottom).

Appendix C: Spectroscopic observations

All Tables

All Figures

	Fig. 1. MG vs GBP − GRP plane for 160 randomly selected source fields (black) filtered using Eqs. (1) and (2) from Arenou et al. (2018). The EBLMs used in the work are also plotted and coloured appropriately.
In the text

Fig. 2. Response function of the HAO+CBB instrument. The atmospheric transmission is plotted in black, the transmission of the HAO telescope in blue-solid, the CBB filter in green and CCD response in yellow. The final response of HAO-1 with the CBB filter is plotted in red-dashed along with the Kepler transmission (blue-dashed). The atmospheric transmission line originated from equations for Rayleigh, aerosol and ozone extinction vs. wavelength for Palomar Observatory (Hayes & Latham 1975). Coefficients were adjusted until they agreed with observations of extinction at HAO over a few dates (2018, priv. comm.).

In the text

	Fig. 3. Difference in theoretical intensity acoss the stellar disc for the HAO+CBB filter and the Kepler/K2 Filter as a function of the angle between a line normal to the stellar surface and the line of sight of the observer (γ) for J1847+39.
In the text

	Fig. 4. Posterior probability distribution of EBLM J2349−32 from photometric fitting. The 1-σ contour are shown.
In the text

	Fig. 5. Difference between observed and calculated magnitudes for each EBLM.
In the text

	Fig. 6. Power spectrum of J0218−31 using 3 seasons of WASP photometry.
In the text

	Fig. 7. Orbital solution of J2349−32 with 1-σ contours plotted. We show the median value of each parameter from the Markov chains (green) along with the the solution with the highest log-likliehood (red).
In the text

	Fig. 8. Difference between the spherical model and Roche model of J2308−46 using ELLC.
In the text

	Fig. 9. Posterior probability distributions from EBLMMASS for J2349−32. Over plotted are the 1-σ contours.
In the text

	Fig. 10. Lucky imaging of J2349−32 (red arm) revealing a close companion 1.3″ (blue circle) away at a position angle of 308.6 ± 0.6° (blue circle).
In the text

	Fig. 11. Orbital fit of J2349−32. Top panel: detrended I-band light-curve from the SAAO 1-m telescope (black) with the best fitting transit model (red). Upper-middle panel: phase-folded WASP lightcurve (black). Lower-middle panel: drift-corrected radial velocity measurements (black) with the best model (red). Bottom panel: residuals from radial velocity model measurements.
In the text

	Fig. 12. Lucky imaging of J2308−46 (red arm) revealing a companion 20″ away with a position angle of 208° (blue circle).
In the text

	Fig. 13. Posterior probability distributions for the density and temperature of the primary star in J2308−46 is shown in the top panel. The zero-age main sequence is show(black-dashed) along with the best fitting isochrone (blue-solid) and the respective isochrones for ±1-σ in [Fe/H]. The lower panels show the PPD distributions for M1, M2 and τ with best-fitting double-Gaussian models in red.
In the text

Fig. 14. Orbital fit of J2308−46. Upper panel: R-band transit obtained from the SAAO 1-m telescope (black) with the best fitting transit model (green dashed). We plot the best fitting transit model generated using Gaussian processes (red). Middle panel: phase-folded WASP observations (black) and observations binned into groups of 50 (blue). We also plot the Roche model used to approximate the out-of-transit photometry used to measure transit parameters from WASP photometry (red). Lower panel: drift-corrected radial velocity measurements (black) with the best fitting model (red) and residuals from the best fitting orbital model.

In the text

Fig. 15. Orbital fit for J0218−31. Top panel: transit photometry from CTIO in g′ (blue), r′ (red), i′ (cyan) and z′ (green) with best fitting models shown in black. Upper-middle panel: phase-folded WASP lightcurve. Lower-middle panel: drift-corrected radial velocity measurements from CORALIE with best fitting model plotted in red, along with residuals. Bottom panel: drift-corrected radial velocity measurements during transit (the Rossiter–McLaughlin effect; black) with the best fitting model (red). Error bars have been omitted for clarity.

In the text

	Fig. 16. Posterior probability distributions for the density and temperature of the primary star in J0218−31 is shown in the top panel. The zero-age main sequence is show(black-dashed) along with the best fitting isochrone (blue-solid) and the respective isochrones for ±1-σ in [Fe/H]. Lower panels: PPD distributions for M1, M2 and τ with best-fitting double-Gaussian models in red.
In the text

	Fig. 17. Orbital fit of J1847+39. Top panel: single transits from the HAO in filters CBB (blue), g′ (cyan) and z′ (green) with best fitting models (red). Upper-middle panel: phase-folded WASP lightcurve. Lower-middle panel: drift-corrected radial velocity measurements (black) with the best fitting model (red) and residuals (lower panel).
In the text

	Fig. 18. Orbital fit of J1436−13. Top panel: a single transit obtained from SAAO in R filter (black) and the best fitting transit model (red). Upper-middle panel: phase-folded WASP lightcurve. Lower-middle panel: drift-corrected radial velocity measurements (black) with best fitting model (red) along with residuals (lower panel).
In the text

Fig. 19. Masses and radii of M-dwarfs in five EBLM systems. The 5 Gyr isochrone for solar metallicity from Baraffe et al. (2015) is plotted (black-solid) along with the 5 Gyr isochrone for [M/H] = −0.5 (Baraffe et al. 1998; black-dashed). We plot low-mass M-dwarfs with masses and radii known to better than 10% (from Table 4 of Chaturvedi et al. 2018, and references therein). For J2308−46 and J0218−31 we plot both solutions and label accordingly. We also plot TRAPPIST-1 (Delrez et al. 2018), Proxima Centauri (Anglada-Escudé et al. 2016) and J0555−57 (von Boetticher et al. 2017).

In the text

	Fig. 20. Fractional radius residual PPD for J2349−32, J2308−46, J0218−31 and J1847+39 system along with Gaussian models (upper panel). Lower panel: J1436−13 which has a broader PPD and requires a double Gaussian model. The legends denote the best-fitting value of the radius residual with uncertainty.
In the text

	Fig. A.1. Posterior probability distribution of EBLM J2308−46 (top-left), J0218−31 (top-right), J1847+39 (bottom-left) and J1439−13 (bottom-right) from photometric fitting. Over plotted are the 1-σ contours.
In the text

	Fig. B.1. Lomb-scargle diagrams for each season of WASP photometry for J2349−32 (top), J2308−46 (middle) and J0218−31 (bottom).
In the text

	Fig. B.2. Lomb-scargle diagrams for each season of WASP photometry for J1847+39 (top) and J1436−13 (bottom).
In the text