The atmospheric parameters of FGK stars using wavelet analysis of CORALIE spectra

The WASP project (Pollacco et al. 2006, Wide Angle Search for Planets;) monitors the night sky in search of exoplanet transit signatures around bright and nearby stars. So far, WASP has discovered over 130¹¹ 1 TEPCAT (Southworth 2011) 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 10’s of km s^-1 over an orbit, suggesting that the transit is stellar in origin as opposed to Jovian. These are typically low-mass stars ($\leq 0.5\,M_{\sun}$) 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 $\leq$ M $\leq$ 0.6 $\rm M_{\sun}$) 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 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) is 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 atmospheric studies with transmission spectroscopy (Sedaghati et al. 2016; Southworth et al. 2016). 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_{\rm 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 thi star. For EBLM systems discovered by WASP, these parameters are made with CORALIE spectra using measurements of equivalent widths 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_{\rm eff}$, [Fe/H] and $V\sin i$ 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 per Ångstrom (SNR) 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 SNR ($\propto\sqrt{N_{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 SNR 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_{\rm eff}$, $\log g$ and [Fe/H], respectively for stars with a Gaia magnitude $G_{RVS}$ = 13. Using neural networks to estimate atmospheric properties has 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\,\rm K$, $0.23\,\rm dex$ and $0.16\,\rm dex$ for $T_{\rm eff}$, $\log g$ and $\rm[Fe/H]$).

Our method determines the best-fitting atmospheric parameters ($T_{\rm 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 (Demory et al. 2009; Torres et al. 2010, e.g.). We introduce wavelet decomposition as it applies to a spectrum in Sec. 2 before reviewing our Bayesian approach to determine $T_{\rm eff}$, [Fe/H], $\log g$ and $V\sin i$ in Sec. 3. We show that our method converges and is self-consistent in Sec. 4 and test against a sample of independently analysed FGK stars in Sec. 5.

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 a spectrum, $f(\lambda)$, with variable sized functions (Stumpe et al. 2012). These functions are called daughter wavelets, $\psi_{\rm a,b}(\lambda)$, which are created from a mother wavelet, $\psi(\lambda)$, using a shift-and-scale operation,

here, $a$ is a member of the dyadic sequence,

and $b=kb_{0}$, where $k$ is an integer and $b_{0}$ is chosen to ensure the recovery of $f(\lambda)$. 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$):

The likeness of a wavelet, $\psi_{i,k}$, to a section of the spectrum is given by the wavelet coefficient $WT_{f(\lambda)}(i,k)$ from Eqn. (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, $\{WT_{f(\lambda)}(0,b)\}_{k}$, which represent the likeness of a single scale across the entire spectrum. The power of each scale, $\{WT_{f(\lambda)}(i,b)\}_{k}^{2}$ can be visuallised 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):

where

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[WT_{f(\lambda)}(i,k)\right]^{2}\geq\rm 0.01$) which eliminates information that has little contribution to a signal; this is called wavelet compression.

We do not require Eqn. (4) to determine atmospheric parameters as we perform a $\chi^{2}$ fit using a subset of coefficients from Eqn. (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^{16}$ values are required to sample a spectrum over the wavelength range 450-650nm. We decided to use $2^{17}$ 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 Eqn. (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 Eqn. (4) to the two subsets of coefficients (4–12 and 13–17) is shown in Fig. 2(a). 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. 2(b)) which better represents the boundary between noise and weak lines. We do not show the reconstruction of subset $i=0-3$ in Fig. 2(a) and 2(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\,\rm nm$). We see a slow variation of some wavelet coefficients which correspond to changes in individual spectral line geometries as the each parameter changes. One thing to note is the sensitivity of each parameter; $T_{\rm eff}$ varies the most, followed by $V\sin i$ 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 frequency identification and is widely used in signal processing and data compression (Belmon et al. 2002, e.g.). A Haar wavelet, with a more step-like structure, is more suited to identifying discontinuity and is widely used in computer vision projects (Essaouabi et al. 2009, e.g.). 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.

We use the Markov chain Monte Carlo method to determine the posterior probability distribution for $T_{\rm eff}$, $\log g$, $V\sin i$ and [Fe/H] given an observed spectrum. Our method is a global $\chi^{2}$ 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 450nm–650nm 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 (Smiljanic et al. 2014, Fig. 4;). 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 and 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\,\rm kms^{-1}$. The upper limit of $50\,kms^{-1}$ would need to be extended for hotter stars beyond the Kraft break, but it suitable for this work on late-type stars. Macroturbulence are estimated using equation (5.10) from Doyle 2015 and microturbulence was accounted for at the synthesis stage using equation (3.1) from the same source. Spectra in-between grid points are extracted by trilinear interpolation, broadened to the desired value of $V\sin i$ 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 – 650 nm with 2¹⁷ values (the same as the observed spectrum) and apply Eqn. (3) to obtain the wavelet coefficients $WT_{f(\lambda)}(4-14,k)$ for the model spectra.

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

and

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(\lambda)$ and with $\sigma$ equal to the standard deviation of the spectrum, $\sigma_{f(\lambda)}$ (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 $\sigma_{f(\lambda)}$, and $\sigma_{MC}^{2}$, effectively down-weighting regions with poor blaze corrections. The free parameter $\beta$ 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 be used stellar models or empirical relations to set priors on these atmospheric parameters but we have 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, e.g., helium-rich stars.

We sample the model parameter space using the Markov chain Monte Carlo method, implemented by the python package emcee (Foreman-Mackey et al. 2013) . emcee uses affine-invariant ensemble sampling (Goodman 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_{\rm eff}$, [Fe/H] and $\log g.$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.

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_{\rm eff}$, [Fe/H], $\log g$ and $V\sin i$ selected within the limits of our grid of models. Each spectrum was then re-sampled to have $2^{17}$ 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 for 44 WASP exoplanet hosts ($\overline{\sigma_{\rm\log g}}$ = 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 Sec. 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 $\sim\,0.25\,-\,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 our grid. With no constraints on $\log g$, there were only 2 measurements of $T_{\rm 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\sin i$ when the interpolated value of $V\sin i$ below 0.5 km s^-1; 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\sin i$ values. The internal precision associated with the wavelet method is remarkably high; by taking $1\,\sigma$ values from the cumulative probability distributions we find precisions around 15 K, 0.01 dex, 0.02 dex, and 0.15 $\rm km\,s^{-1}$ for $T_{\rm eff}$, [Fe/H], $\log g$, and $V\sin i$ 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_{\rm 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^{-1}$ for $V\sin i$. 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_{\rm 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.

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. The D15 sample consists of 24 stars analysed by measurements of equivalent width and spectral fitting of high signal-to-noise and high-resolution ($\rm R\,=\,112,000$) data from the HARPS spectrograph (Queloz et al. 2001). We used lower-quality observations from the CORALIE spectrograph to determine $T_{\rm eff}$, [Fe/H], $\log g$ and $V\sin i$ 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).

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 equivalent widths from higher-quality HARPS spectra. From this we determine a precision for the parameters derived from the wavelet method of $85\,\rm K$ for $T_{\rm eff}$, $0.06\,\rm dex$ for [Fe/H] and $1.35\,\rm kms^{-1}$ for $V\sin i$. 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 as an indicator of whether or not a star is a dwarf star ($\log g\approx 4.5$) or not. 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 is applied at a correction to the [Fe/H] values from our method. We find our method is robust for échelle spectra with a SNR 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 star currently known. This method is also being used to study other EBLM systems and as part of the on-going exoplanet discovery process with for 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 SNR from other echelle spectrographs. We have implemented the wavelet method in a python module called waveletspec which is available upon request.