© ESO 2018

1 Introduction

The WASP project (Wide Angle Search for Planets; Pollacco et al. 2006) monitors the night sky in search of exoplanet transit signatures around bright and nearby stars. So far, WASP has discovered over 130¹ exoplanet systems and continues to be the origin of many exoplanet discoveries. The SuperWASP instruments are sensitive to Jovian planets transiting solar-type stars which periodically attenuates the brightness of the host. Candidate exoplanet systems are subject to radial velocity measurements from the CORALIE échelle spectrograph (Queloz et al. 2000) to determine whether the transiting companion is indeed of planetary mass. Occasionally, these measurements yield radial velocity differences of tens of km s⁻¹ over an orbit, suggesting that the transit is stellar in origin as opposed to Jovian. These are typically low-mass stars (≤0.5 M_⊙) which have radii comparable to large planets, hence mimicking an exoplanet transit signal very well. These systems are given the flag “EBLM” (eclipsing binary, low mass) in the database used to coordinate observations to exclude them from the planet-hunting process.

M-dwarfs (0.08 ≤ M ≤ 0.6 M_⊙) are the most abundant stars in the Galaxy (Henry et al. 2006). They slowly descend to the main sequence along a near-vertical Hyashi track where they can remain on the main sequence for billions of years. Those which are small enough to be fully convective restrict the build-up of helium ash in their cores, resulting in main sequence lifetimes greater than the age of the Universe (Baraffe et al. 1998).

Numerous photometric measurements around the infrared region are required to sample the peak of the spectral energy distribution(SED) of M-dwarfs. Studies of the secondary eclipses of M-dwarfs around larger FGK stars in the infrared have found temperatures much hotter than evolutionary models predict (Triaud et al. 2013; Torres 2013). Empirical measurements of radius share a similar issue, in that they are typically underestimated by at least 3% (Gómez Maqueo Chew et al. 2014; Demory et al. 2009; Spada et al. 2013; Fernandez et al. 2009; Nefs et al. 2013). The favoured explanation is that enhanced magnetic activity in M-dwarf stars inhibits the convective flow of energy to the surface, resulting in an enlarged radius (Morales et al. 2010; Fernandez et al. 2009; Zhou et al. 2014; Kraus et al. 2011). Many measurements of mass and radius are done from eclipsing binary systems, and some have suggested that tidal interactions (via rapid rotation) are the root of the problem (Kervella et al. 2016). However, there is evidence for radius inflation in single stars (Spada et al. 2013) suggesting a real problem with evolutionary models.

M-dwarfs are becoming prime targets for transiting exoplanet surveys. With all else constant, the transit of an exoplanet in the habitable zone is more likely to be observed around a star of lower mass to produce deeper eclipses that enable planetary atmosphericstudies with transmission spectroscopy (Sedaghati et al. 2016; Southworth et al. 2017). Exciting new discoveries such as the TRAPPIST-1 system (Gillon et al. 2017) and Proxima Centauri b (Anglada-Escudé et al. 2016) have attracted attention from the scientific community and public alike, bringing more focus to low-mass stars in exoplanet surveys.

We are presented with an opportunity to use EBLMs discovered by the WASP survey to measure the properties of M-dwarfs. Measurements of T_eff and [Fe/H] for the FGK star can be measured accurately from a well understood spectrum and can be combined with radial-velocity measurements and transit photometry to obtain the mass, radius, and temperature of both components in EBLM systems. Measurements of [Fe/H] for the M-dwarf can be adopted from the FGK companion with the assumption that they both formed from the same parent molecular cloud (Triaud et al. 2011).

The EBLM project (Triaud et al. 2013; Gómez Maqueo Chew et al. 2014; von Boetticher et al. 2017) is an ongoing effort to characterise transiting M-dwarfs in EBLM systems. Accurate estimates for temperature and composition are needed to estimate limb-darkening coefficients and make mass estimates for the primary star. For EBLM systems discovered by WASP, these parameters are made with CORALIE spectra using measurements of equivalent widths (EWs) and by fitting individual spectral lines (Gillon et al. 2009; Doyle et al. 2013; Doyle 2015). The individual analysis of each spectrum is acceptable for a small sample size, but we require a reliable automated procedure to measure T_eff, [Fe/H], and V sini for the entirety of the EBLM database to keep pace with future EBLM discoveries. At the same time, such a method needs to allow for noise and systematic errors present in the CORALIE spectra. The sample of EBLM spectra presented in Gómez Maqueo Chew et al. (2014) typicaly have a signal-to-noise ratio (S/N) per Ångstrom of between 3 and 7. The on-going radial velocity campaign to study EBLMs typically yields between 10 and 40 spectra per star. Co-adding spectra can increase the S/N ($\propto \sqrt{N_{\text{obs}}}$) to over 40 in some parts of the spectrum, but the regions of the spectrum near the ends of each echelle order suffer from both large photon noise and systematic errors due to inaccurate order-merging.

Wavelet decomposition has been used previously as part of methods developed for spectral analysis. Manteiga et al. (2010) used multi-level wavelet decomposition in connectionist systems (artificial neural networks) to derive fundamental stellar parameters in the low-S/N domain (5–25) in preparation for spectra from the Gaia radial velocity spectrograph (RVS). This work was extended by Dafonte et al. (2016) by using a generative artificial neural network resulting in predicted uncertainties of 220 K, 0.32 dex and 0.20 dex for T_eff, log g and [Fe/H], respectively for stars with a Gaia magnitude G_RVS = 13. Using neural networks to estimate atmospheric properties entails well-known problems such as long training times and a strong dependence on the initial training set. Li et al. (2015) use wavelet decomposition in a regression framework to detect representative spectral features from a spectrum and estimate atmospheric parameters with better precision than those from neural networks (83 K, 0.23 dex and 0.16 dex for T_eff, log g and [Fe∕H]).

Our method determines the best-fitting atmospheric parameters (T_eff, [Fe/H], V sin i and log g) for FGK stars by comparing a selected subset of coefficients from a wavelet decomposition to those from a grid of stellar models. This reduces systematic errors in the estimated parameters due to poor continuum normalisation and low-quality regions of the spectrum. These measurements can then be combined with photometric follow-up observations to obtain the mass and radius of M-dwarfs in EBLMs to an accuracy of a few percent. These, in turn, provide calibratable points for empirical mass-radius relations of low-mass stars (e.g. Demory et al. 2009; Torres et al. 2010). We introduce wavelet decomposition as it applies to a spectrum in Sect. 2 before reviewing our Bayesian approach to determine T_eff, [Fe/H], log g and V sin i in Sect. 3. We show that our method converges and is self-consistent in Sect. 4 and test against a sample of independently analysed FGK stars in Sect. 5.

2 Wavelet analysis

Analysis of spectral components at different scales can be done using a discrete wavelet transform (DWT). A DWT tiles the wavelength-scale plane by convolving aspectrum, f(λ), with variable sized functions (Stumpe et al. 2012). These functions are called daughter wavelets, ψ_a,b (λ), which are created from a mother wavelet, ψ(λ), using a shift-and-scale operation, $${\psi }_{a,b}(\lambda )=\frac{1}{\sqrt{a}}\psi \left(\frac{\lambda -b}{a}\right),a,b\in \Re ,a\ne 0,$$(1) where a is a member of the dyadic sequence, $$a_i=2^i,i=0,1,2,3,\mathrm{...},n$$(2) and b = kb₀, where k is an integer and b₀ is chosen to ensure the recovery of f(λ). By employing a DWT, the appropriate values of b are selected to minimise overlap between wavelet convolutions. Following the notation in Chapter 8 of Olkkonen (2011), a discrete wavelet transform can be calculated for each dyadic scale (i) and displacement (k): $$W{T}_{f(\lambda )}(i,k)=\frac{1}{\sqrt{2^i}}{\displaystyle \int f}(\lambda )\overline{\psi \left(\frac{\lambda -k{2}^i{b}_0}{2^i}\right)}\text{d}\lambda =f(\lambda ),{\psi }_{i,k}(\lambda ).$$(3)

The likeness of a wavelet, ψ_i,k, to a section of the spectrum is given by the wavelet coefficient WT_f(λ)(i, k) from Eq. (3). Performing this calculation over the series of dyadic scales and displacements yields wavelet coefficients which represent different sized structures at different wavelengths. We split coefficients into bands with constant scales, ${\{W{T}_{f(\lambda )}(0,b)\}}_k$, which represent the likeness of a single scale across the entire spectrum. The power of each scale, ${\{W{T}_{f(\lambda )}(i,b)\}}_k^2$, can be visualised in a power Hövmoller (one value of i per row) in Fig. 1. Bands of coefficients which correspond to noise and low-order continuum artefacts (such as merged échelle orders) can then be excluded. A filtered spectrum may be reconstructed with an inverse DWT (IDWT): $$f(\lambda )={\displaystyle \sum _{i=-\infty }^{\infty }{2}^{\frac{-3i}{2}}}{\displaystyle \int W}{T}_{f(\lambda )}(i,k)\widehat{\psi }\left(\frac{\lambda -b}{2^i}\right)\text{d}b,$$(4)

where $$\widehat{\psi }\left(\frac{\lambda -b}{2^i}\right)=\frac{\psi \left(\frac{\lambda -b}{2^i}\right)}{{\displaystyle \sum _{i=0}^{i=n}{\left|\psi \left(\lambda -b\right)\right|}^2}}.$$(5)

The process of reconstructing a spectrum using a subset of wavelet coefficients is called wavelet filtering and is analogous with Fourier filtering. Alternatively, the subset of coefficients may be chosen to meet a threshold criteria (i.e. ${\left[W{T}_{f(\lambda )}(i,k)\right]}^2\ge 0.01$) which eliminates information that has little contribution to a signal; this is called wavelet compression.

We do not require Eq. (4) to determine atmospheric parameters as we perform a χ² fit using a subset of coefficients from Eq. (3) to those from a grid of models (see Sect. 3). We also do not apply any threshold criterion. The nominal resolving power of the CORALIE spectrograph is R = 55 000, so at least 2¹⁶ values are required to sample a spectrum over the wavelength range 450–650 nm. We decided to use 2¹⁷ values for the wavelet decomposition to ensure no loss of information and to give us more choice in the number of wavelet bands used in our analysis. We use Eq. (3) to obtain wavelet coefficients which have information on scales in the range 0.003 nm–200 nm. Our wavelet method only uses a subset of i values. To select these, we constructed power Hövmoller diagrams (similar to Fig. 1) for a variety of regions between 450 nm and 650 nm, for different co-added spectra in our sample. We found that power associated with line absorption lies in the range 0.04–4 nm, with larger scales typically corresponding to systematic trends and shorter scales with noise. This corresponds to values of i = 4−12 (0.048–3.125 nm). The application of Eq. (4) to the two subsets of coefficients (4–12 and 13–17) is shown in Fig. 2a. We find that the subset range i = 4−12 is too restrictive to reproduce short-scale information (e.g., weak lines) and so we decide to extend this range to i = 4−14 (0.012–3.125 nm; Fig. 2b) which better represents the boundary between noise and weak lines. We do not show the reconstruction of subset i = 0−3 in Figs. 2a and b as using only 16 coefficients to reconstruct a spectrum leads to a large Daubechies-4 wavelet with some sub-structure.

We demonstrate the sensitivity of wavelet coefficients to atmospheric parameters in Fig. 3 for wavelet coefficients in the range i = 11−12 (0.04–0.09 nm). We see a slow variation of some wavelet coefficients which corresponds to changes in individual spectral line geometries as each parameter changes. One thing to note is the sensitivity of each parameter; T_eff varies the most, followed by V sini and [Fe/H]. Surface gravity is the least varying parameter in wavelet space and is dominated by a few lines sensitive to log g. As V sin i increases, we see positive and negative structures form and become stronger at higher V sin i values. This is likely to be a continuum effect as weaker lines are smeared out to average a lower continuum whilst stronger lines persist.

The choice of mother wavelet depends on the objective of the work. A Daubechies wavelet performs well for frequencyidentification and is widely used in signal processing and data compression (e.g. Belmon et al. 2002). A Haar wavelet, with a more step-like structure, is more suited to identifying discontinuity and is widely used in computer-vision projects (e.g. Essaouabi et al. 2009). We investigate the effects of wavelet choice on the determined atmospheric parameters in Sect. 5.2.1, but proceed with the Daubechies (k = 4) wavelet for the rest of this paper.

Fig. 1 Power Hövmoller of wavelet coefficients (lower panel) for a region around the Mg triplet for WASP-19 (upper panel). There is significant power (|WTi,k| from Eq. 3) for scales ~1 nm in the region of the Mg lines corresponding the wavelets likeness to spectral features. Horizontal red lines represent the scales 0.012 − 3.125 nm.

Fig. 2 Reconstruction of spectra using Eq. (4) for subsets of wavelet coefficients. (Left panel - top) Raw spectra for WASP-19 (black) and the flux reconstruction using wavelet coefficients from bands i = 4–12 using the raw spectrum (blue; offset − 0.6) and the best fitting model for WASP-19 (red; offset − 1.2). (Left panel -bottom) The reconstruction of the best-fitting model for WASP-19 (red) and the raw spectrum (blue; offset + 0.5) using coefficients i = 13–17. (Right panel) As in the left panel except with reconstructions using coefficients i = 4–14 (top) and coefficients i = 15–17 (bottom).

Fig. 3 Changes in wavelet coefficients in the range i = 11−12 for model spectra as a function of atmospheric parameters. The wavelet coefficient from a solar model has been subtracted to emphasise the subtle change in wavelet coefficients for each parameter. A similar result is seen for other values of i between 4 and 14. The colour-bar indicates the magnitude of the difference of coefficients. We note that these bars are not on the same scale and highlight the wavelet response to each coefficient.

3 Bayesian measurements

We use the Markov chain Monte Carlo method to determine the posterior probability distribution for T_eff, log g, V sin i and [Fe/H] given an observed spectrum. Our method is a global χ² fitting routine which compares subsets of wavelet coefficients (i = 4−14) to those from a pre-synthesised grid of spectra. Our grid was synthesised with the radiative transfer code SPECTRUM (Gray & Corbally 1994) using MARCS model atmospheres (Gustafsson et al. 2008), and version 5 of the GES (Gaia ESO survey) atomic line list provided within iSpec (Blanco-Cuaresma et al. 2016) with solar abundances from Asplund et al. (2009). We computed models spanning 450–650 nm over a temperature range of 4000 to 8000 K in steps of 250K, −1 to +1 dex in steps of 0.5 dex for [Fe/H] and 3.5 to 5 dex in steps of 0.5 for log g. We selected our range of [Fe/H] by looking at composition measurements of over 14 000 FGK stars from Gaia-ESO survey data release 3 (Fig. 4; Smiljanic et al. 2014). We find that 96% of stars with measurements of composition had [Fe/H] in the range −1 to 1 dex. This range in [Fe/H] is also much larger than the full range in [Fe/H] for our benchmark sample described in Sect. 5.

Spectra in the grid are calculated with zero instrumental, rotational, or macroturbulence broadening. These are accounted for in post-processing by convolving the grid spectra with the appropriate kernels. In this work, we allow V sin i to have values in the range 0–50 km s⁻¹. The upper limit of 50 km s⁻¹ would need to be extended for hotter stars beyond the Kraft break, but is suitable for this work on late-type stars. Macroturbulence are estimated using Eq. (5.10) from Doyle (2015) and microturbulence was accounted for at the synthesis stage using Eq. (3.1) from the same source. Spectra in-between grid points are extracted by trilinear interpolation, broadened to the desired value of V sini and macroturbulence, and then convolved with a Gaussian to account for instrumental broadening. For the self-consistency tests in Sect. 4 instrumental broadening was ignored, but for the CORALIE spectra in Sect. 5 we used an instrumental resolving power R = 55 000 (Queloz et al. 2001; Doyle 2015). We then re-sample between 450 and 650 nm with 2¹⁷ values (the same as the observed spectrum) and apply Eq. (3) to obtain the wavelet coefficients WT_f(λ)(4−14, k) for the model spectra.

The subset of wavelet coefficients from the interpolated model, WT_m, are compared to those from the data, WT_d, in the following Bayesian framework: the probability of observing a spectrum for a given model is given by $\text{p}(m|d)\propto \mathcal{L}(d|m)\text{p}(m)$. The vector of model parameters is given by $m=\left(T_{\text{eff}},[\text{Fe}/\text{H}],\text{log}g,V\sin \text{i}\right)$ and we assume uniform prior probability for the model parameters within the grid range. We use the likelihood function $\mathcal{L}(d|m)=\exp (-{\chi }^2/2)$ where $${\chi }^2=\frac{{(W{T}_d-W{T}_m)}^2}{{\sigma }_{W{T}_d}^2},$$(6) and $${\sigma }_{W{T}_d}^2=\beta {\sigma }_{MC}^2.$$(7)

The term ${\sigma }_{MC}^2$ was calculated by generating 1000 spectra from the co-added spectrum with noise generated from a standard normal distribution centred around f(λ) and with σ equal to the standard deviation of the spectrum, σ_f(λ) (calculated from the standard deviation in co-added spectra). The blaze function is corrected prior to co-addition of the spectra and so deviations in blaze functions will result in uncertainties propagating through to σ_f(λ), and ${\sigma }_{MC}^2$, effectively down-weighting regions with poor blaze corrections. The free parameter β has been introduced to account for additional noise, incomplete atomic data, deviations from solar metallicity scaling, lines which form under non-local thermodynamic equilibrium, and other unaccounted errors. In principle, we could have used stellar models or empirical relations to set priors on these atmospheric parameters but decided not to do this for two reasons. Firstly, allowing the MCMC sampler to explore regions with a priori low probability gives a better indication of the reliability of our method than using a more constrained solution. Secondly, by imposing a prior from stellar models or empirical relations based on normal stars we may fail to identify interesting examples of anomolous stars in our sample, for example, helium-rich stars.

We sample the model parameter space using the Markov chain Monte Carlo method, implemented by the python packageemcee (Foreman-Mackey et al. 2013). emcee uses affine-invariant ensemble sampling (Goodman & Weare 2010, parallel stretch move algorithm) to split Markov chains into sub-groups and update the position of a chain using the positions of chains in the other subgroups. The algorithms affine-invariance can cope with skewed probability distributions and generally has shorter autocorrelation times than a classic Metropolis-Hastings algorithm.

We generate 12 Markov chains of 20 000 draws each to converge on the best atmospheric parameters. We found that the chains converged before the 5000 th draw, but as a precaution we discarded the first 10 000 draws. We take the median values of the model parameters in the remaining draws to determine the atmospheric parameters for a spectrum. An example posterior probability distribution for WASP-20 is plotted in Fig. 5. The parameter space is almost symmetric with small degeneracies between T_eff, [Fe/H] and log g. We note that the precision of the parameters determined from the standard deviation of each parameter in the Markov Chain is typically an underestimate of the true precision of these parameters because it does not account for systematic errors in the data or the models.

Fig. 4 Histogram of 14,681 [Fe/H] measurements for stars from Gaia-ESO data release 3 (Smiljanic et al. 2014). Plotted is the median value of [Fe/H] (solid blue), with 1σ from the median (dashed blue). The grid range used in enclosed by the dashed green lines.

Fig. 5 Posterior proability distributions for WASP-20.

4 Self consistency

We have assessed the ability of our method to recover atmospheric parameters from synthetic spectra in order to check that our results are self consistent. We interpolated 512 spectra with random values of T_eff, [Fe/H], log g and V sini selected within the limits of our grid of models. Each spectrum was then re-sampled to have 2¹⁷ values in the range 450–650 nm to match our choice of coefficients in Sect. 2 and the benchmark sample in Sect. 5. These spectra were split into two groups and analysed with the aforementioned method. The first group had log g as a free parameter to probe for any systematics, for the second group we imposed a prior on log g to simulate the effect of well constrained surface gravity measurement from transit photometry. The log g prior probability distribution was assumed to be Gaussian with a mean log g value equal to the value used to interpolate the spectrum and a dispersion equal to the average uncertainty of transit log g values from Mortier et al. (2013, hereafter referred to as M13) for 44 WASP exoplanet hosts ($\overline{{\sigma }_{\text{logg}}}$ = 0.02 dex). We decided not to add Gaussian noise to these spectra as noise profiles depend upon stellar parameters and instrumental conditions; this is assessed in Sect. 5.2.4. We find typical autocorrelation lengths are below 1000 steps for all parameters in the first chain and 12 chains in the second run typically produce an acceptance fraction between ~ 0.25 and 0.3.

The recovery of atmospheric parameters for both groups is shown in Fig. 6 and summarised in Table 1. We find that all parameters are recovered well across the range of our grid. With no constraints on log g, there were only two measurements of T_eff that deviated from the input value by more than 150 K. A prior on log g significantly decreases the difference between measured and input atmospheric parameters and shows that our method is sensitve to log g. There is a small increase in residual scatter for measurements of V sini when the interpolated value of V sini below 0.5 km s⁻¹; this is seen in both groups and marginally improved with a prior on log g. This is expected as the resolution of the broadening kernel in combination with the edge of parameter space makes it difficult to determine low V sini values. The internal precision associated with the wavelet method is remarkably high; by taking 1 σ values from the cumulative probability distributions we find precisions around 15 K, 0.01 dex, 0.02 dex, and 0.15 km s⁻¹ for T_eff, [Fe/H], log g, and V sini respectively.More realistic uncertainties are determined in the following sections.

We also assess the sensitivity of our method by determining the atmospheric parameters of 9 spectra from a discrete set of grid points with different combinations of fixed parameters. We interpolate 96 spectra from the following grid points – 4800 K, 5800 K and 6500 K for T_eff; − 0.5, 0.0 and 0.5 dex for [Fe/H]; 3.75, 4.40 and 4.80 dex for log g; 5, 10 and 15 kms⁻¹ for V sini. In total, we determined the atmospheric parameters for each spectrum 15 times with the wavelet method using every combination of free and fixed parameters (see Fig. 7). We find that constraining one or more parameters increases the internal precision of the wavelet method significantly. A slight degeneracy exists between T_eff and log g resulting in a modest scatter when both parameters are left free. This also highlights the numerical noise introduced by starting walkers at different positions, since walkers explore parameter space by random jumps which may never reach the correct solution, despite prior knowledge that an exact solution lies somewhere within the grid.

Fig. 6 Differences between wavelet-determined atmospheric parameters and those used to synthesize spectra with all parameters free (black) and with priors on log g (red).

Fig. 7 Difference between atmospheric parameters determined by the wavelet method for 9 synthetic spectra when some parameters are fixed. We plot atmospheric parameters determined with with all parameters free (Teff, [Fe/H], log g, V sin i) in black; with one parameter fixed in blue; with two parameters fixed in cyan; three parameters fixed in orange.

5 Benchmark sample

Any spectral analysis technique must be tested against stars with high-quality measurements. For this we use stars from Doyle et al. (2013) and Doyle (2015, D15, hereafter). The D15 sample consists of 24 stars analysed by measurements of EW and spectral fitting of high-S/N and high-resolution (R = 112 000) data from the HARPS spectrograph (Queloz et al. 2001). We used lower-quality observations from the CORALIE spectrograph to determine T_eff, [Fe/H], log g and V sini of the same stars with the wavelet method. Only 22 stars in the D15 sample have CORALIE spectra available to use and we further exclude WASP-77A and the close (3′′) B-component as both component as they are un-resolved in the CORALIE fibre. This leaves a sample of 20 stars for use to calibrate our method (see Table 2).

5.1 CORALIE spectra

Each spectrum was processed with the CORALIE standard data-reduction pipeline (Baranne et al. 1996). The radial velocity shift was measured relative to a solar template² and corrected into a laboratory frame of reference. The spectra were then median-combined onto an identically sampled wavelength grid. Continuum regions were identified by applying maximum and median filters (Blanco-Cuaresma et al. 2014) and fitted with spline functions (1 every 10 nm) for normalisation. The wavelet method was then applied to each spectrum twice: once with no priors on log g and a second time with priors given by transit photometry. The priors on log g were set to those from M13 if quoted, or the relevant discovery papers otherwise (see Table 6).

5.2 Results

The results can be seen in Fig. 8 and are summarised in Tables 3 and 6. Our method determines T_eff to within 220 K of the value found by D15. Our measurements of [Fe/H] are systematically offset by approximately −0.18 dex from those of D15; this is discussed further in Sect. 5.2.1. It is difficult to constrain log g spectroscopically and our measurements often differ from those of D15 by up to 0.5 dex. Our measurements of V sini converge to 0 km s⁻¹ for seven stars in the sample due to an over estimation of v_mac or instrumental resolution. To mitigate this problem we repeated the analysis with v_mac = 0 km s⁻¹. This allowed these stars to converge on best fitting models without pushing against the edge of parameter space. These stars are marked with an asterisk in Table 6.

We find no benefit by using priors on log g. In most cases the use of log g priors increases the standard deviation in differences between atmospheric parameters from our method andpublished values. We investigated the level of agreement between spectroscopic values of log g from EW measurements (D15), log g from the wavelet method and those from transit photometry. Photometric surface gravity is typically measured to better precision than its spectroscopic counterpart, but relies on stellar models and correct limb-darkening parameters which in-turn rely on a constrained effective temperature, composition, and surface gravity. Recent work suggests a disagreement between spectroscopic and photometric log g which is correlated with T_eff (see Fig. 4 from Doyle et al. 2017). We compare the difference between spectroscopic and photometric log g measurements in Fig. 9. We find a statistically significant negative correlation (p-value ≤ 10⁻⁵) between Δlog g (log g_photometry −log g_wavelet) with T_eff from our method. The origin of this is unclear, but a similar trend is seen between spectroscopic and asteroseismic measurements (Fig. 6 of Mortier et al. 2014) which suggests to us that this is a problem to spectral analysis of late-type stars using plane-parallel non-LTE model atmospheres. For a few stars, we relaxed the log g prior to have a standard deviation of 0.2 dex (instead of 0.02 dex) and found almost no difference between these solutions and those with a uniform prior on log g.

In Fig. 10 we assess the Hα region for the model predicted from D15 and this work for the highest-quality spectrum in our sample, WASP-20, with S∕N = 150. The results from D15 were obtained using a custom line list, whereas we use version 5 of the GES atomic line list provided with iSpec to synthesise the D15 model of WASP-20 using atmospheric parameters, ν_mac and ν_mic from D15. We find both modelsagree well with the data, with the left wing fitting best and a underestimation in the right wing. The discrepancies between the two wings of the Hα line seen here are the result of the difficulty in calibrating the blaze function in this region of the spectrum. We find that the majority of Fe line depths are under-predicted with the wavelet method, with the D15 model better matching individual line profiles. This test demonstrates the need to benchmark against well studied stars and visually inspect the best models against the data.

Fig. 8 Difference between wavelet analysis and D15 (O–C) for each atmospheric parameter in the D15 sample. Each spectrum was measured twice, once with log g as a free parameter (black) and again with log g priors imposed from transit photometry (red). We exclude measurements of V sin i where macroturbulence, ξt, was set to 0 km s−1 to ensure models a best model was converged upon.

Fig. 9 Difference between spectroscopic log g and photometric log g (log gph - log gsp ) correlated with Teff,wavelet from this work (black) and from D15 (red).

Fig. 10 H-α region for WASP-20 (black) fitted with the best fitted model from D15 (red) and the best model from this work (blue). The near horizontal lines at flux 1.2 are the residuals between the D15 model (red) or wavelet model (blue) and the spectrum of WASP-20.

5.2.1 Systematic offset in [Fe/H]

There are many reasons why our method may be produce composition offsets compared with other established techniques. The interested reader should see Jofré et al. (2017) for an excellent review on how the specifics of spectroscopic analysis routines affect abundance measurements. One interesting result from Jofré et al. (2017) is the effect of continuum normalisation which increased the method-to-method scatter in abundance measurements by up to 0.3 dex (see their Fig. 5). Wavelet filtering in our method is an alternate approach to normalisation, and so an offset of around 0.18 dex is not entirely unexpected. We assess if there is a systematically lower continuum placement by adding an free parameter, C₀, which is a constant to add to the normalised flux of the model spectra before a discrete wavelet transform in the calculation of log-likelihood. We found values of C₀ converged to values between − 0.05 and 0.05 and did not affect measurements of [Fe/H] by more than 0.05 dex; T_eff remained the same for all stars within 150 K and log g changed by as much as 0.2 dex.

We also looked at components unique to the wavelet method. For instance, the mother wavelet used (Daubechies, k = 4) may not capture the true line depths when convolved with a spectrum. We again measured WASP-20 with three alternative wavelets (Daubechiesk = 20, Harr k =2 and bspline k= 103) across the range 450 – 650 nm (see Table 4). We find the choice of mother wavelet has little influence on the determined composition (and all other atmospheric parameters) for WASP-20 and we find similar results for the rest of the D15 sample. It is possible that the resolution of the finest wavelet convolution (2 pixels) is not sufficient to capture iron line depths. To assess this, we convolved a few iron lines with the Daubechies k = 4 kernel and assessed whether line depths were underdetermined. We found this not to be the case, suggesting no degradation of line depths owing to the choice in wavelets.

Finally, we consider the possibility that there may be instrumental effects at play with the CORALIE échelle spectrograph. A discrepancy in EW measurements for WASP-69 (see Fig. 3.19 in D15) suggests this instrument is prone to scattered light (Doyle 2015). This may be partly responsible for the systematic error in the iron abundance when combined with a low-quality spectrum.

The zero-point of the metallicity scale is a subject of on-going debate (e.g. Kraft & Ivans 2004). However, we can conclude that models using parameters found by D15 (as generated with line lists and atmospheres used in the above work) have better-fitting line depths for the majority of Fe lines in the D15 sample than our predicted models. For this reason, we advocate using the following correction for the [Fe/H] values measured with the wavelet method to make them consistent with the metallicity scale of D15: $${[\text{Fe}/\text{H}]}_{\text{corrected}}={[\text{Fe}/\text{H}]}_{\text{measured}}+\mathrm{0.18.}$$(8)

5.2.2 Systematic trend in log g

We also observe a negative correlation between residual log g measurements (wavelet – D15) and log g measured with the wavelet method (Fig. 8c). This trend is observed with and without Gaussian priors on log g from transitphotometry. We calculate a Pearson correlation coefficient of –0.501 for measurements with no log g prior, suggesting a significant negative correlation. We fit this trend with a first-order polynomial and found a gradient of − 0.692 and a y-intercept of 3.067. This correlation evaluates to zero at a wavelet log g value of 4.44. In principle, the following correction can be used to bring our log g measurements into line with those from D15, $$\text{log}{g}_{\text{corrected}}=\text{log}{g}_{\text{wavelet}}-3.067+0.692\times \text{log}{g}_{\text{wavelet}}.$$(9)

Without knowing the exact cause of this trend, and given the sensitivity of our log g estimates to the continuum placement, we are reluctant to advise applying this correction and conclude that the wavelet method cannot reliably estimate log g beyond confirming a dwarf-like surface gravity. Obtaining log g from a spectrum is typically done through ionization balance (balancing the iron abundance measured from the Fe I and Fe II lines). It is also possible to measure log g by fitting the wings of gravity sensitive lines (e.g. Mg, Na) using model spectra (the synthesis method). This is essentially how the wavelet method operates (in wavelet space rather than normalised flux space). Accurate determinations of log g from the synthesis method requires detailed element-abundance measurements for gravity-sensitive Na and Mg lines. Estimating the abundances of these elements by scaling from the solar abundance values and applying some correction for α element enhancement will lead to a systematic error in log g that is difficult to quantify in individual cases. To investigate this further requires another set of comparison stars with independent log g values (preferably from binary systems where log g can be accurately measured and not planet transiting systems).

5.2.3 Precision of atmospheric parameters

The high precision of the parameters in Table 6 shows that the wavelet method can reliably converge to a well-determined set of atmospheric parameters, but to make use of these parameters we also require a reliable estimate of their true precision that accounts for additional uncertainties due to systematic errors in the data and the models. To obtain a realistic estimate of true precision of the parameters from the wavelet method, σ_wavelet, we compare theresults from our method with the correction to [Fe/H] described earlier to those from D15. The standard deviation of the residuals between the measured atmospheric parameters made by D15 and from the wavelet method, σ_{D 15−wavelet}, is a combinationof the uncertainties from methods added in quadrature: $${\sigma }_{\text{D}15-\text{wavelet}}^2={\sigma }_{\text{D}15}^2+{\sigma }_{\text{wavelet}}^2,$$(10) where σ_D15 is the quoted error on the atmospheric parameters from D15. There are two extreme cases: the first is that the uncertainty from D15 is negligible (or at least much better than what we can achieve) giving ${\sigma }_{\text{D}15-\text{wavelet}}^2\approx {\sigma }_{\text{wavelet}}^2$; and the second is that the inter-method discrepancy, ${\sigma }_{\text{D}15-\text{wavelet}}^2$, is negligible leaving uncertainties similar to those quoted by D15. In reality, the absolute uncertainty for the wavelet method is somewhere between these two extremes. We adopt a true precision of each parameter from Table 3 using a uniform prior on log g which is to assume that σ_{D 15}⟨⟨σ_wavelet. We suggest applying a correction of + 0.18 dex to [Fe/H] and not to apply a correction to log g. This means precision of 85 K for T_eff, 0.06 dex for [Fe/H] and 1.35 kms⁻¹ for V sini. The resulting value of log g is not likely to be reliable but is good enough to confirm dwarf-like gravity around log g = 4–5 dex. We note that these values are comparable to other methods (e.g. Bruntt et al. (2010)).

5.2.4 Spectrum quality

In Fig. 11 we plot the difference between atmospheric parameters obtained with the wavelet method (with no priors for log g) to those from D15 as a function of S/N. The sample falls into two categories of quality (those with S/N ≤ 90 and those with S/N ≥ 120). There is noticeably more scatter in the lower-quality group and suggests that the uncertainty of our atmospheric parameters decreases with a better-quality spectrum. The noise profile of a spectrum depends on observing conditions, properties of the star and the instrument used to make the observations. This is why adding Gaussian noise to a synthetic spectrum until the atmospheric parameters are no longer recoverable does not give a true reflection of a methods robustness to noise. Instead, we use 32 (out of 58) observations of the star with the highest S/N in the D15 sample – WASP-20. We dyadically split up these spectra and median combine them into different sets. The sets of splits used were 1 spectrum (1 set of 32 spectra), 2 spectra (2 sets of 16 spectra), 4 spectra (4 sets of 8 spectra), ..., 32 spectra (32 sets each containing just 1 spectra). We scale S/N from the coaddition of all 58 spectra: $$S/N=S/N_{58\text{\hspace{0.05em}spectra}}\times \sqrt{\frac{N_{\text{set}}}{58}}.$$(11)

Each set was measured with the aforementioned wavelet technique with no prior probability function for log g, and best fitting parameters adopted. The precision and accuracy as a function of S/N are shown in Figs. 12 and 13, respectively. We find that systematic errors dominate for a S/N below 40. A similar result is found by Smiljanic et al. (2014) who measured UVES-FLAMES spectra for FGK stars from the Gaia-ESO survey and found a systematics threshold of S∕N ≈ 50.

Fig. 11 Atmospheric parameters from the wavelet method, with no prior on log g, compared to those from D15 (xD15 − xwavelet) as a function of S/N.

Fig. 12 Precision of the wavelet method versus S/N for Teff (top left), [Fe/H] (top right), log g (bottom left), and V sini (bottom right) for WASP-20.

Fig. 13 Accuracy of the wavelet method versus S/N for Teff (top left), [Fe/H] (top right), log g (bottom left), and V sini (bottom right) for WASP-20 using the results from D15 as a zero-point.

6 Conclusion

We have shown that our method accurately recovers the atmospheric parameters of synthetic spectra from a grid of models using subsets of wavelet coefficients in a Bayesian framework. The same method was applied to the CORALIE spectra of 20 FGK stars which have been analysed independently by measurements of EWs from higher-quality HARPS spectra. From this we determine a precision for the parameters derived from the wavelet method of 85 K for T_eff, 0.06 dex for [Fe/H] and 1.35 kms⁻¹ for V sini. Surface gravity, log g, can also be estimated using our method but it is difficult to assess the precision of this parameter in individual cases. Consequently, we recommend that log g estimates from our method are only used to decide whether or not a star is a dwarf (log g ≈ 4.5). We find an offset in our metallicity scale compared to the results of Doyle et al. (2013); Doyle (2015) in the sense that our values of [Fe/H] are lower by 0.18 dex, despite using a consistent solar abundance, and recommend that this offset be applied as a correction to the [Fe/H] values from our method. We find our method is robust for échelle spectra with a S/N above 40. Below this value the uncertainity in the measured atmospheric parameters increases to unusable levels. A further development of this method would include a more sophisticated weighting system for the wavelet coefficients beyond the Monte Carlo approach used here.

Our method has already been used to determine the atmospheric parameters of the EBLM J0555−57 (von Boetticher et al. 2017), which hosts one of the densest main sequence stars currently known. This method is also being used to study other EBLM systems and as part of the on-going exoplanet discovery process of the WASP survey. For both exoplanet systems and EBLM binaries, the contribution of the companion star to the optical flux is negligible (they are SB1 binaries) and so our method using models of single stars is appropriate, but it would not be suitable for cases where the companion is detectable in the spectrum (SB2 binaries). We have optimised our method for the application to spectra from CORALIE, but the same method should be equally applicable to spectra with moderate S/N from other echelle spectrographs. We have implemented the wavelet method in a python module called waveletspec which is available upon request.

Acknowledgements

SG acknowledges financial support from the Science and Technology Facilities Council (STFC). SG thanks A. H. M. J. Triaud and S. Udry for data acquisition and support. We thank the referee for their careful reading of the manuscript and suggestions that have helped to improve this paper.

References

All Tables

All Figures

	Fig. 1 Power Hövmoller of wavelet coefficients (lower panel) for a region around the Mg triplet for WASP-19 (upper panel). There is significant power (|WTi,k| from Eq. 3) for scales ~1 nm in the region of the Mg lines corresponding the wavelets likeness to spectral features. Horizontal red lines represent the scales 0.012 − 3.125 nm.
In the text

Fig. 2 Reconstruction of spectra using Eq. (4) for subsets of wavelet coefficients. (Left panel - top) Raw spectra for WASP-19 (black) and the flux reconstruction using wavelet coefficients from bands i = 4–12 using the raw spectrum (blue; offset − 0.6) and the best fitting model for WASP-19 (red; offset − 1.2). (Left panel -bottom) The reconstruction of the best-fitting model for WASP-19 (red) and the raw spectrum (blue; offset + 0.5) using coefficients i = 13–17. (Right panel) As in the left panel except with reconstructions using coefficients i = 4–14 (top) and coefficients i = 15–17 (bottom).

In the text

Fig. 3 Changes in wavelet coefficients in the range i = 11−12 for model spectra as a function of atmospheric parameters. The wavelet coefficient from a solar model has been subtracted to emphasise the subtle change in wavelet coefficients for each parameter. A similar result is seen for other values of i between 4 and 14. The colour-bar indicates the magnitude of the difference of coefficients. We note that these bars are not on the same scale and highlight the wavelet response to each coefficient.

In the text

	Fig. 4 Histogram of 14,681 [Fe/H] measurements for stars from Gaia-ESO data release 3 (Smiljanic et al. 2014). Plotted is the median value of [Fe/H] (solid blue), with 1σ from the median (dashed blue). The grid range used in enclosed by the dashed green lines.
In the text

	Fig. 5 Posterior proability distributions for WASP-20.
In the text

	Fig. 6 Differences between wavelet-determined atmospheric parameters and those used to synthesize spectra with all parameters free (black) and with priors on log g (red).
In the text

	Fig. 7 Difference between atmospheric parameters determined by the wavelet method for 9 synthetic spectra when some parameters are fixed. We plot atmospheric parameters determined with with all parameters free (Teff, [Fe/H], log g, V sin i) in black; with one parameter fixed in blue; with two parameters fixed in cyan; three parameters fixed in orange.
In the text

	Fig. 8 Difference between wavelet analysis and D15 (O–C) for each atmospheric parameter in the D15 sample. Each spectrum was measured twice, once with log g as a free parameter (black) and again with log g priors imposed from transit photometry (red). We exclude measurements of V sin i where macroturbulence, ξt, was set to 0 km s−1 to ensure models a best model was converged upon.
In the text

	Fig. 9 Difference between spectroscopic log g and photometric log g (log gph - log gsp ) correlated with Teff,wavelet from this work (black) and from D15 (red).
In the text

	Fig. 10 H-α region for WASP-20 (black) fitted with the best fitted model from D15 (red) and the best model from this work (blue). The near horizontal lines at flux 1.2 are the residuals between the D15 model (red) or wavelet model (blue) and the spectrum of WASP-20.
In the text

	Fig. 11 Atmospheric parameters from the wavelet method, with no prior on log g, compared to those from D15 (xD15 − xwavelet) as a function of S/N.
In the text

	Fig. 12 Precision of the wavelet method versus S/N for Teff (top left), [Fe/H] (top right), log g (bottom left), and V sini (bottom right) for WASP-20.
In the text

	Fig. 13 Accuracy of the wavelet method versus S/N for Teff (top left), [Fe/H] (top right), log g (bottom left), and V sini (bottom right) for WASP-20 using the results from D15 as a zero-point.
In the text