Disentangling Time Series Spectra with Gaussian Processes: Applications to Radial Velocity Analysis

TL;DR

This paper introduces a novel Gaussian process-based spectral disentangling method for radial velocity analysis, improving orbital parameter precision and enabling reconstruction of stellar spectra without relying on templates.

Contribution

The paper presents a new spectral disentangling technique using Gaussian processes that simultaneously estimates orbital parameters and stellar spectra from time-series data.

Findings

Achieved improved orbital parameter precision for LP661-13

Successfully reconstructed primary and secondary spectra

Demonstrated potential for broader stellar and exoplanet applications

Abstract

Measurements of radial velocity variations from the spectroscopic monitoring of stars and their companions are essential for a broad swath of astrophysics, providing access to the fundamental physical properties that dictate all phases of stellar evolution and facilitating the quantitative study of planetary systems. The conversion of those measurements into both constraints on the orbital architecture and individual component spectra can be a serious challenge, however, especially for extreme flux ratio systems and observations with relatively low sensitivity. Gaussian processes define sampling distributions of flexible, continuous functions that are well-motivated for modeling stellar spectra, enabling proficient search for companion lines in time-series spectra. We introduce a new technique for spectral disentangling, where the posterior distributions of the orbital parameters and intrinsic, rest-frame stellar spectra are explored simultaneously without needing to invoke cross-correlation templates. To demonstrate its potential, this technique is deployed on red-optical time-series spectra of the mid-M dwarf binary LP661-13. We report orbital parameters with improved precision compared to traditional radial velocity analysis and successfully reconstruct the primary and secondary spectra. We discuss potential applications for other stellar and exoplanet radial velocity techniques and extensions to time-variable spectra. The code used in this analysis is freely available as an open source Python package.

Equations

$$\mathbf{x}\sim \mathcal{N}\left({\mathbf{\mu }}_x,{\mathbf{\Sigma }}_{xx}\right)$$

$$p(\mathbf{x}|{\mathbf{\mu }}_x,{\mathbf{\Sigma }}_{xx})=\frac{1}{{\left[(2\pi {)}^N\det {\mathbf{\Sigma }}_{xx}\right]}^{1/2}}\exp \left(-\frac{1}{2}(\mathbf{x}-{\mathbf{\mu }}_x{)}^{\mathrm{T}}{\mathbf{\Sigma }}_{xx}^{-1}(\mathbf{x}-{\mathbf{\mu }}_x)\right)$$

$$\ln p(\mathbf{x}|{\mathbf{\mu }}_x,{\mathbf{\Sigma }}_{xx})=-\frac{1}{2}\left[(\mathbf{x}-{\mathbf{\mu }}_x{)}^{\mathrm{T}}{\mathbf{\Sigma }}_{xx}^{-1}(\mathbf{x}-{\mathbf{\mu }}_x)+\ln \det {\mathbf{\Sigma }}_{xx}+N\ln 2\pi \right]$$

$$\left[\begin{array}[]{c}{\bm{x}}\\ {\bm{y}}\end{array}\right]\;\sim\;{\cal N}\left(\left[\begin{array}[]{c}{\bm{\mu}_{x}}\\ {\bm{\mu}_{y}}\end{array}\right],\left[\begin{array}[]{cc}{\bm{\Sigma}}_{xx}&{\bm{\Sigma}}_{xy}\\ {\bm{\Sigma}}_{yx}&{\bm{\Sigma}}_{yy}\end{array}\right]\right),$$

$$\mathbf{x}|\mathbf{y}\sim \mathcal{N}\left({\mathbf{\mu }}_x+{\mathbf{\Sigma }}_{xy}{\mathbf{\Sigma }}_{yy}^{-1}(\mathbf{y}-{\mathbf{\mu }}_y),{\mathbf{\Sigma }}_{xx}-{\mathbf{\Sigma }}_{xy}{\mathbf{\Sigma }}_{yy}^{-1}{\mathbf{\Sigma }}_{yx}\right).$$

$$\displaystyle{\bm{\lambda}}=\left[\begin{array}[]{c}\lambda_{1}\\ \lambda_{2}\\ \vdots\\ \lambda_{w}\\ \end{array}\right]\qquad{\bm{d}}=\left[\begin{array}[]{c}d_{1}\\ d_{2}\\ \vdots\\ d_{w}\\ \end{array}\right]$$

$$f(\lambda )\sim \mathrm{GP}(\mu (\lambda ),k(\lambda ,{\lambda }^{\prime })).$$

$$\mathbf{d}=$$

$$\mathbf{d}\sim$$

$$r_{ij}=r({\lambda }_i,{\lambda }_j)=\frac{c}{2}\left|\frac{{\lambda }_i-{\lambda }_j}{{\lambda }_i+{\lambda }_j}\right|,$$

$$k_{ij}(r_{ij}|a,l)=a^2\exp \left(-\frac{r_{ij}^2}{2l^2}\right),$$

$${\mathbf{f}}_{\ast }\sim \mathcal{N}\left(1,{\mathbf{\Sigma }}_f+{\mathbf{\Sigma }}_N\right).$$

$$\displaystyle{\bm{\Lambda}}=\left[\begin{array}[]{c}{\bm{\lambda}}_{1}\\ {\bm{\lambda}}_{2}\\ \vdots\\ {\bm{\lambda}}_{W}\\ \end{array}\right],\qquad{\bm{D}}=\left[\begin{array}[]{c}{\bm{d}}_{1}\\ {\bm{d}}_{2}\\ \vdots\\ {\bm{d}}_{W}\\ \end{array}\right],\qquad{\bm{F}}=\left[\begin{array}[]{c}{\bm{f}_{1}}\\ {\bm{f}_{2}}\\ \vdots\\ {\bm{f}_{W}}\\ \end{array}\right]$$

$$\mathbf{D}=\mathbf{F}+\mathbf{N}.$$

$$\mathbf{\lambda }(v)=\sqrt{\frac{c+v}{c-v}}{\mathbf{\lambda }}_0$$

$$M(t)=E(t)-e\sin E(t)=\frac{2\pi t-T_0}{P}$$

$$\cos f(t)=\frac{\cos E(t)-e}{1-e\cos E(t)}.$$

$$v_A=K_A\left(\cos (\omega +f(t))+e\cos \omega \right)+\gamma$$

$$v_B=-\frac{K_A}{q}\left(\cos (\omega +f(t))+e\cos \omega \right)+\gamma .$$

$${\bm{\Lambda}}=\left[\begin{array}[]{c}{\bm{\lambda}}_{1}\\ {\bm{\lambda}}_{2}\\ \vdots\\ {\bm{\lambda}}_{W}\\ \end{array}\right],\qquad{\bm{D}}=\left[\begin{array}[]{c}{\bm{d}_{1}}\\ {\bm{d}_{2}}\\ \vdots\\ {\bm{d}_{W}}\\ \end{array}\right],\qquad{\bm{F}}=\left[\begin{array}[]{c}{\bm{f}_{1}}\\ {\bm{f}_{2}}\\ \vdots\\ {\bm{f}_{W}}\\ \end{array}\right],\qquad{\bm{G}}=\left[\begin{array}[]{c}{\bm{g}_{1}}\\ {\bm{g}_{2}}\\ \vdots\\ {\bm{g}_{W}}\\ \end{array}\right]$$

$$\mathbf{D}=\mathbf{F}+\mathbf{G}+\mathbf{N}.$$

$${\bm{\Lambda}}_{F}=\left[\begin{array}[]{c}{\bm{\lambda}}_{1,A}\\ {\bm{\lambda}}_{2,A}\\ \vdots\\ {\bm{\lambda}}_{W,A}\\ \end{array}\right],\qquad{\bm{\Lambda}}_{G}=\left[\begin{array}[]{c}{\bm{\lambda}}_{1,B}\\ {\bm{\lambda}}_{2,B}\\ \vdots\\ {\bm{\lambda}}_{W,B}\\ \end{array}\right].$$

$${\mathbf{\Sigma }}_{xx}={\mathbf{\Sigma }}_F+{\mathbf{\Sigma }}_G+{\mathbf{\Sigma }}_N,$$

$$\left[\begin{array}[]{c}{\bm{f}}_{\ast}\\ {\bm{g}}_{\ast}\\ {\bm{D}}\\ \end{array}\right]\;\sim\;{\cal N}\left(\left[\begin{array}[]{c}{\bm{\mu}}_{f}\\ {\bm{\mu}}_{g}\\ {\bm{\mu}}_{FG}\\ \end{array}\right],\left[\begin{array}[]{c@{\hspace{2em}}c@{\hspace{2em}}c}{\bm{\Sigma}}_{f}\hfil\hskip 20.00003pt&{\bm{0}}\hfil\hskip 20.00003pt&{\bm{\Sigma}}_{fF}\\ {\bm{0}}\hfil\hskip 20.00003pt&{\bm{\Sigma}}_{g}\hfil\hskip 20.00003pt&{\bm{\Sigma}}_{gG}\\ {\bm{\Sigma}}_{Ff}\hfil\hskip 20.00003pt&{\bm{\Sigma}}_{Gg}\hfil\hskip 20.00003pt&({\bm{\Sigma}}_{F}+{\bm{\Sigma}}_{G}+{\bm{\Sigma}}_{N})\\ \end{array}\right]\right),$$

$$\left.\left[\begin{array}[]{c}{\bm{f}}_{\ast}\\ {\bm{g}}_{\ast}\\ \end{array}\right]\right|{\bm{D}}\;\sim\;{\cal N}\left(\left[\begin{array}[]{c}{\bm{\mu}}_{f}\\ {\bm{\mu}}_{g}\\ \end{array}\right]+{\bm{\Sigma}}_{xy}{\bm{\Sigma}}_{yy}^{-1}({\bm{D}}-{\bm{\mu}}_{FG}),{\bm{\Sigma}}_{xx}-{\bm{\Sigma}}_{xy}{\bm{\Sigma}}_{yy}^{-1}{\bm{\Sigma}}_{yx}\right)$$

$${\bm{\Sigma}}_{xx}=\left[\begin{array}[]{cc}{\bm{\Sigma}}_{f}&{\bm{0}}\\ {\bm{0}}&{\bm{\Sigma}}_{g}\\ \end{array}\right],\quad\quad{\bm{\Sigma}}_{yy}={\bm{\Sigma}}_{F}+{\bm{\Sigma}}_{G}+{\bm{\Sigma}}_{N},\quad\quad{\bm{\Sigma}}_{xy}=\left[\begin{array}[]{c}{\bm{\Sigma}}_{fF}\\ {\bm{\Sigma}}_{gG}\\ \end{array}\right].$$

$$\hat{\mathbf{F}}=\frac{\mathbf{F}+{\mathbf{\mu }}_f}{{\mathbf{\mu }}_f}\hat{\mathbf{G}}=\frac{\mathbf{G}+{\mathbf{\mu }}_g}{{\mathbf{\mu }}_g}.$$

Full text

techniques: radial velocity — techniques: spectroscopic — binaries : spectroscopic — stars: fundamental parameters — celestial mechanics — stars: individual (LP661-13)

††journal: the Astrophysical Journal††facilities: FLWO:1.5m/TRES††software: Python, astropy, PSOAP

1 Introduction

Close binary pairs are the foundation of stellar astrophysics. Measurements of the orbital dynamics in these systems can constrain masses, the fundamental physical parameter in stellar evolution. That information is vital in our understanding of everything from star formation to death, with wide-reaching implications for issues ranging from molecular cloud collapse and fragmentation to exoplanets to cosmology with Type Ia supernovae (Torres et al., 2010). A common means of measuring binary orbital dynamics is through high resolution spectroscopic radial velocity monitoring. Traditionally, radial velocities for each stellar component are measured by cross-correlating an observed composite spectrum with various Doppler-shifted stellar templates (e.g., TODCOR; Zucker & Mazeh, 1994). The velocity for each component corresponds to the Doppler shift which delivers the maximum cross-correlation signal. While cross-correlation is commonly employed as a workhorse technique, there are several shortcomings to this limited statistical framework. First, in the case of low signal to noise data, there can be considerate uncertainty about how well the chosen templates match the true spectra of the stars. In its most straightforward form, cross-correlation is unable to meaningfully account for variable spectral lines or uncertain calibration parameters in a principled probabilistic framework, nuisances which can systematically bias a radial velocity measurement. Moreover, at no point does the cross-correlation framework reconstruct either of the component spectra, preventing a check of the suitability of the chosen template as well as any detailed photospheric analysis of the stars themselves.

Because the composite spectra of spectroscopic binaries are modulated by a deterministic Doppler shift, it is possible to “disentangle” the intrinsic spectra from multiple observations at different orbital phases. One of the first successful attempts was by Bagnuolo & Gies (1991): by viewing the composite spectra of the double-lined O star AO Cas as various projections of the intrinsic component spectra, they were able to apply iterative tomographic reconstruction techniques and recover the component spectra. Soon after, Simon & Sturm (1994) developed a sparse matrix formalism to decompose composite spectra into their components while also optimizing the orbital elements of the binary star. In contrast to traditional radial velocity techniques, spectral disentangling techniques can determine the intrinsic spectra and velocities simultaneously, simplifying the coherent propagation of uncertainties and often resulting in more precise constraints on the orbital parameters (Pavlovski & Hensberge, 2010). Spectral decomposition can provide precise radial velocities even at orbital phases where the lines from both stars overlap, an area of difficulty for cross-correlation techniques. Once disentangled, the component spectra can be further analyzed with conventional spectroscopic techniques to determine the fundamental stellar properties. For example, Rawls et al. (2016) disentangled the spectra of the double red giant eclipsing binary KIC 9246715 and then used the radiative transfer code MOOG (Sneden, 1973) to estimate the effective temperature, surface gravity, and metallicity for each star.

A major advance in spectral disentangling was the realization that the reconstruction could be performed quickly in the Fourier domain, thanks to the efficiency of the FFT (KOREL; Hadrava, 1995). However, the FFT introduces some potentially undesirable side effects. Because the FFT treats the spectrum as a continuous, periodic function, it is important to carefully choose “chunks” of the spectrum such that the edges are at the same continuum level, otherwise spectral lines will redshift and blueshift off the edge of the chunk and distort the reconstruction (Ilijic, 2004). This problem can be particularly acute when dealing with a late spectral type star with a high density of spectral lines. While Fourier techniques provide the ability to filter out high frequency noise (de-noising), they may also have difficulty constraining low-order “continuum”-like features of the spectra, making reconstructed spectra appear wavy (Hadrava, 1995). Lastly, the FFT requires that input spectra be interpolated to a log-$\lambda$ spaced wavelength grid (a process that introduces correlated noise) and treats the measurement errors as homoscedastic. This means that if the signal-to-noise of the spectrum is wavelength-dependent or there happens to be a cosmic-ray hit on a certain region of the detector, there is no straightforward way to adjust the weighting of specific portions of the spectrum.

The flexible probabilistic framework offered by Gaussian processes potentially provides a means to address many of the limitations of traditional orbit measurement and spectral disentangling techniques. In truth, stellar spectra are not actually physically generated from Gaussian processes but are rather the result of non-linear radiative transfer through complex stellar atmospheres with specific wavelength-dependent elemental and molecular opacities. However, for the purposes of inferring intrinsic spectra and stellar radial velocities, simple Gaussian processes offer an attractive framework to model stellar spectra in a purely data driven manner. Gaussian processes have been used successfully for other time-series applications in astronomy, for example modeling lensed quasar time delays (Hojjati et al., 2013; Tak et al., 2016), inferring stellar rotation periods (Angus et al., 2015), and modeling correlated noise in photometric observations of planet transits and eclipses (Evans et al., 2015; Montet et al., 2016).

The content of this paper is as follows: in §2, we introduce Gaussian processes and demonstrate how they may be used to model a stellar spectrum. Then, we model a mock double-lined spectroscopic binary–where spectral lines from both components are seen in the composite spectrum–and demonstrate how to simultaneously infer the orbital parameters and the intrinsic spectra of both stars. We also show how the precision of the orbital posteriors and quality of the reconstructed spectra respond to changes in the binary flux ratio and signal-to-noise of the dataset. In §3 we apply our technique to the mid-M binary system LP661-13, recently studied by Dittmann et al. (2016). We demonstrate precise inference of the orbital parameters and reconstruct the spectra of A and B, which are an excellent match to other mid-M spectral templates. In §4, we discuss potential extensions of the Gaussian process framework to variable stellar spectra, telluric line modeling, and precision radial velocity measurement for exoplanet detection, and in §5 we conclude the paper.

2 Gaussian Process Spectral Models

In this section, we describe a framework for modeling a stellar spectrum, and sums of stellar spectra, as Gaussian processes. First, we review common notation and theorems for multidimensional Gaussian random variables. Second, we explain how a Gaussian process can be used to model a spectrum of a single star which is stationary in time and shows no orbital motion due to the presence of a companion. Third, we introduce a Keplerian orbital model and extend the Gaussian process framework to model the spectral time series of a double-lined spectroscopic binary. Throughout this section we use archival observations of a single star to simulate observations and demonstrate the development of the framework.

We adopt the notation that

[TABLE]

signifies that the vector ${\bm{x}}$ is drawn from a multi-dimensional Gaussian distribution with mean vector ${\bm{\mu}}$ and covariance matrix ${\bm{\Sigma}}_{xx}$. The elements of ${\bm{\Sigma}}_{xx}$ are the covariances between pairwise elements within ${\bm{x}}$ and can be specified directly or via a functional prescription. By definition, the likelihood function associated with ${\bm{x}}$ is the multi-dimensional Gaussian

[TABLE]

where $N$ is the length of ${\bm{x}}$. For computational reasons, the natural logarithm of the likelihood is frequently used

[TABLE]

Following Rasmussen & Williams (2005, e.q. A.5), we let ${\bm{x}}$ and ${\bm{y}}$ be jointly Gaussian random vectors drawn from

[TABLE]

where ${\bm{\mu}}_{x}$ and ${\bm{\mu}}_{y}$ are the vector means. The sub-matrices ${\bm{\Sigma}}_{xx}$ and ${\bm{\Sigma}}_{yy}$ are the covariance matrices corresponding to the elements within ${\bm{x}}$ and ${\bm{y}}$, and ${\bm{\Sigma}}_{xy}$ encapsulates the covariances between the elements in ${\bm{x}}$ and ${\bm{y}}$; ${\bm{\Sigma}}_{yx}={\bm{\Sigma}}_{xy}^{\mathrm{T}}$. If and only if ${\bm{x}}$ and ${\bm{y}}$ are independent will ${\bm{\Sigma}}_{xy}={\bm{0}}$. The conditional distribution of ${\bm{x}}$ given ${\bm{y}}$ is also a Gaussian distribution

[TABLE]

These equations are the foundation for constructing the Gaussian process spectra model that follows.

2.1 A model for observations of a single star

For our spectroscopic application, the input vector is the sampled wavelengths of the detector $\{\lambda_{i}\}_{i=1}^{w}$ and the data vector is the observed fluxes $\{d_{i}\}_{i=1}^{w}$

[TABLE]

where each pixel is indexed by $i$ from $1$ to $w$, the number of pixels in the region of the spectrum under consideration. To illustrate the development of our framework, we use synthetic “mock” datasets generated by the following recipe. First, we create a template by stacking many archival high-resolution observations of the K5 star LkCa 14 and smoothing the result with a Gaussian kernel. Then, we sample the spectrum at ${\bm{\lambda}}$ and add a known amount of independent white noise to the dataset.

We will model the continuous, intrinsic stellar spectrum as a Gaussian process. A function $f(\lambda)$, where $\lambda>0$, is said to have a Gaussian process with mean function $\mu(t)$ and covariance kernel $k(\lambda,\lambda^{\prime})$ if for any finite collection of inputs $0<\lambda_{1}<\lambda_{2}<\ldots<\lambda_{w}$ the vector $\{f(\lambda_{1}),\,f(\lambda_{2}),\,\ldots,\,f(\lambda_{w})\}$ has a multivariate Gaussian distribution with mean vector $\{\mu(\lambda_{1}),\,\mu(\lambda_{2}),\,\ldots,\,\mu(\lambda_{w})\}$ and $w\times w$ covariance matrix with elements $k(\lambda_{i},\lambda_{j})$, where $i,j=1,2,\ldots,w$, and $k$ is a positive definite kernel function. A single function generated from a Gaussian process is called a realization of the Gaussian process. We consider the continuous, intrinsic stellar spectrum $f$ to be a function generated from the Gaussian process

[TABLE]

For the purposes of this paper, we will always work with finite-length samplings of the Gaussian process ${\bm{f}}$, either at the detector wavelengths ${\bm{\lambda}}$ or finely spaced vectors of our own design. To model a single epoch of a single star, we treat the data ${\bm{d}}$ as the sum of a realization ${\bm{f}}$ of the Gaussian process with a realization of the noise process ${\bm{N}}$,

[TABLE]

where ${\bm{\Sigma}}_{N}$ is a covariance matrix describing the noise of the dataset. For most spectroscopic observations, this will simply be a diagonal matrix with the Poisson uncertainty for each pixel, unless a particular telescope reduction pipeline provides additional information about inter-pixel covariances. It is assumed that the spectrum has been rectified to have a mean flux of 1, and so by design we choose the Gaussian process to have ${\bm{\mu}}_{f}={\bm{1}}$. In theory, we could allow the mean to vary but in practice it has little effect on the result as long as a reasonable number (or function of $\lambda$) is chosen to be $\sim 1$. The covariance matrix ${\bm{\Sigma}}_{f}$ is populated by evaluating a kernel function between pairs of input values in ${\bm{\lambda}}$. The choice of the kernel function determines the smoothness of the functions drawn from the Gaussian process. As a distance metric, we compute the velocity distance between two pairs of wavelengths $i$ and $j$ as

[TABLE]

where $c$ is the speed of light and $r_{ij}$ has units of $\mathrm{km\;s}^{-1}$. A popular kernel is the squared exponential kernel, which we choose to model our spectra. This kernel specifies the covariance between two pixels $\lambda_{i}$ and $\lambda_{j}$ as

[TABLE]

where $a$ sets the amplitude of the Gaussian process and has the same units as the rectified flux, and $l$ sets the length scale of the Gaussian process and has units of $\mathrm{km\;s}^{-1}$. In this paper, we refer to $a$ and $l$ as Gaussian processes hyperparameters.

Now we have fully specified the components of the sampling distribution for the Gaussian process $f$ (Equation 1), which can be thought of as a distribution over functions. We use the ∗ subscript to denote realizations of a vector, which is a draw from this distribution

[TABLE]

Several realizations of ${\bm{f}}_{*}$ are shown in the top panels of Figure 2.1, for two different choices of the Gaussian process hyperparameters $a$ and $l$. While the Gaussian process has a mean of 1, actual realizations of the Gaussian process scatter about the mean. As $l$ decreases, the covariance between distant elements decreases and ${\bm{f}}_{*}$ oscillates more rapidly. How does one choose the best values of $a$ and $l$? Fortunately, the Gaussian process framework provides a natural mechanism to determine the most probable values through the Gaussian process likelihood, given by Equation 2.111We will use Equations 2, 4, & 5 throughout this paper and the values of ${\bm{\Sigma}}_{xx}$, ${\bm{\Sigma}}_{yy}$, and ${\bm{\Sigma}}_{xy}$ will change depending on the context. In the presence of data (${\bm{x}}={\bm{d}}$), a choice of mean vector (${\bm{\mu}}_{x}={\bm{1}}$), a covariance matrix (${\bm{\Sigma}}_{xx}={\bm{\Sigma}}_{f}+{\bm{\Sigma}}_{N}$), and flat priors on the hyperparameters, the posterior probability distribution is directly specified by the marginal likelihood, Equation 2, which we can maximize with respect to $a$ and $l$.

With the hyperparameters optimized, we can now make realizations of the Gaussian process conditioned on the data by taking random draws from the posterior predictive distribution (Equation 5). These draws represent realizations of the inferred intrinsic stellar spectrum while the scatter in the draws serves to illustrate the uncertainty in our inference. In this application of Equation 5, ${\bm{x}}={\bm{f}}_{\ast}$; ${\bm{y}}={\bm{d}}$; ${\bm{\Sigma}}_{xx}$ is filled out by the kernel evaluated over the prediction wavelengths ${\bm{\lambda}}_{\ast}$ corresponding to ${\bm{f}}_{\ast}$; ${\bm{\Sigma}}_{yy}={\bm{\Sigma}}_{f}+{\bm{\Sigma}}_{N}$, where ${\bm{\Sigma}}_{f}$ is filled out by the kernel evaluated over the wavelengths ${\bm{\lambda}}$ corresponding to ${\bm{d}}$; and ${\bm{\Sigma}}_{xy}$ is filled out by the kernel evaluated over pairs of wavelengths corresponding to both ${\bm{\lambda}}_{\ast}$ and ${\bm{\lambda}}$. Random draws of the Gaussian process “snap” to the data (Figure 2.1, third panel), providing a very flexible mechanism for modeling the continuous stellar spectrum.

Although the Gaussian process achieves promising success in this single-observation application, its real benefit accrues when there are multiple observations of the same star in a time series. Throughout this paper, we assume that all spectroscopic observations of a target star are acquired by the same telescope, and that the instrument line spread function (LSF) is stable across epochs. For the applications discussed in this paper, these criterion are satisfied by most high-resolution spectrographs, such as echelle spectrographs used for radial velocity planet-searches. Because convolution is a linear operation, it is reasonable to ignore the effects of the line spread function and assume that the intrinsic stellar spectrum we are modeling has already been convolved. For very precise applications (e.g. radial velocity search for earth-mass planets), variations in the line spread function may need to be considered; we discuss this further in §4.1.

Even though spectrographs may be very stable, such that from epoch to epoch nearly the same wavelengths are sampled by the detector in the topocentric frame, this is not the case in the barycentric frame. Due to the rotation of the earth and its orbital motion around the sun, observations taken of a source on different days will likely have different relative velocities between the topocentric frame and the barycentric frame, the so-called barycentric correction. This means that at each epoch, the pixels of the detector are sampling slightly different wavelengths of the intrinsic stellar spectrum (broadened by the LSF). This barycentric correction can be computed to high accuracy (e.g., Wright & Eastman, 2014), and when all spectra in a time series are corrected to the barycentric frame we have a very dense (albeit noisy) sampling of the convolved stellar spectrum (see Figure 2.1, bottom panel). We concatenate the $W$ epochs of data and denote new variables

[TABLE]

and

[TABLE]

The aforementioned covariance matrices are generated from the same kernel as before, with the major difference that it is now applied to all of the input wavelengths in all epochs, which creates a covariance matrix structure seen in Figure 2.1. In principle, we could sort the data in increasing wavelength to create a covariance matrix with a single, central diagonal band; instead we choose to keep adjacent wavelengths in the same epoch next to each other for later purposes. Posterior predictive draws for the multi-epoch Gaussian process are shown in the bottom panel of Figure 2.1.

2.2 A model for observations of a spectroscopic binary

Now, we extend the Gaussian process model to include multiple, time series observations of a spectroscopic binary star system. Typically, labeling a system a spectroscopic binary implies that spectral lines from both stars are seen in the composite spectra, however we will later show that a simplified version of this framework can also be used to model single-lined spectroscopic binary systems (e.g., stars hosting exoplanets). The orbital motion of each star induces a Doppler shift of the rest-frame stellar spectra, which are then simultaneously observed in composite. If a star is moving with a radial velocity $v$ relative to the barycentric frame, the rest-frame wavelengths ${\bm{\lambda}}_{0}$ of the spectrum are shifted to

[TABLE]

where a positive $v$ denotes a redshift, or increase in the values of the rest frame wavelengths. The radial velocities of binary stars as a function of time can be fully described by seven parameters $\theta$ = {$q$, $K_{A}$, $e$, $\omega$, $P$, $T_{0}$, $\gamma$ }, the mass ratio $q=M_{B}/M_{A}=K_{A}/K_{B}$, the velocity semi-amplitude of the primary star, the eccentricity, argument of periastron, orbital period, epoch of periastron, and systemic velocity, respectively (Murray & Correia, 2010). The mean anomaly $M$ is given by Kepler’s equation

[TABLE]

which must be solved to find the eccentric anomaly $E$. The true anomaly $f$ is given by

[TABLE]

Then, the velocity of the primary star as a function of time is

[TABLE]

and the secondary velocity is

[TABLE]

The resulting spectroscopic binary dataset has the same dimensionality as the single star dataset, but we now assume that it is produced from the sum of the spectrum of the primary star, denoted by a realization $f$ of a Gaussian process, and the spectrum of the secondary star, denoted by a realization $g$ of an (independent) Gaussian process. The secondary spectrum is produced from a Gaussian process similar to $f$ (e.g., Equation 15) with the same form of covariance kernel (Equation 20). We denote the various components of the spectroscopic binary dataset as

[TABLE]

[TABLE]

This spectroscopic binary application can be thought of as modeling the rest-frame spectrum of the primary star (star A), and separately that of the secondary star (star B), while only having access to datasets which represent the sum of the spectra at different orbital phases, all in the presence of noise. We are assuming that the two Gaussian processes are independent from each other and that the Doppler shifts from the orbital motion will allow us to disentangle them. For demonstration purposes, we assume that the orbital parameters $\theta$ are known, and therefore we can calculate the radial velocities of the primary and secondary at all epochs. To fill out the covariance matrix, we must first create new input vectors corresponding to $f$ and $g$ that relate the observed flux at each epoch to shifted versions of the rest-frame spectra of A and B,

[TABLE]

To produce these vectors, the original sub-components of ${\bm{\Lambda}}$ are shifted by the predicted radial velocity for each star at each epoch according to Equation 36, such that ${\bm{\Lambda}}_{F}$ and ${\bm{\Lambda}}_{G}$ should correspond to the wavelengths of star A and B in their respective rest-frames. To be explicit with notation, these input vectors are actually a function of the orbital parameters, ${\bm{\Lambda}}_{F}(\theta)$ and ${\bm{\Lambda}}_{G}(\theta)$.

The likelihood function for the sum of $f$ and $g$ is still given by Equation 2, however the covariance matrix is now a composite of two covariance matrices plus the noise matrix,

[TABLE]

where ${\bm{\Sigma}}_{F}$ is evaluated using ${\bm{\Lambda}}_{F}$ and ${\bm{\Sigma}}_{G}$ is evaluated using ${\bm{\Lambda}}_{G}$ with the kernel specified in Equation 20. We show how ${\bm{\Sigma}}_{xx}$ is constructed in Figure 2.2. The off-diagonal submatrices of ${\bm{\Sigma}}_{F}$ and ${\bm{\Sigma}}_{G}$ denote the cross-epoch covariances of each Gaussian process. Now that we have completely specified the likelihood function for disentangling the component spectra of a spectroscopic binary, we relax the assumption that $\theta$ is fixed. For appropriate values of $a$ and $l$, the likelihood function will be maximized when $\theta$ is chosen such that the Keplerian orbit delivers velocity shifts for each epoch that make the cross-matrices consistently describe the rest-frame spectrum of each star.

With the addition of reasonable priors, we can also explore the full posterior distribution of the orbital parameters $\theta$ and the Gaussian process hyperparameters. If we so desire, we can allow separate values of $a$, $l$ for each spectrum (i.e., $\{a,l\}_{f}$, $\{a,l\}_{g}$) or enforce that the Gaussian process hyperparameters are the same. If the spectra of A and B are drastically different in morphology (e.g. an A-type star coupled with an M-type star, or a rapidly rotating primary coupled with a slowly rotating secondary), using separate hyperparameters will yield a more optimal spectral reconstruction. We can also relax the assumption that the Gaussian process hyperparameters are the same at all wavelengths, for example to adapt to a star where the spectral line density changes significantly across the observed wavelength range, which we discuss further in 4.

Next, we focus on the formalism necessary to reconstruct the rest-frame spectra of A and B at a vector of rest-frame wavelengths ${\bm{\lambda}}_{\ast}$. The joint probability distribution of the (independent) latent Gaussian processes and the dataset is given by

[TABLE]

where ${\bm{\Sigma}}_{f}$ and ${\bm{\Sigma}}_{g}$ are evaluated over pairs of wavelengths in ${\bm{\lambda}}_{\ast}$; ${\bm{\Sigma}}_{F}$ and ${\bm{\Sigma}}_{G}$ are evaluated over pairs of wavelengths in ${\bm{\Lambda}}_{F}$ and ${\bm{\Lambda}}_{G}$, respectively; and ${\bm{\Sigma}}_{fF}$ and ${\bm{\Sigma}}_{gG}$ are evaluated over cross-pairs of wavelengths in (${\bm{\lambda}}_{\ast}$, ${\bm{\Lambda}}_{F}$) and (${\bm{\lambda}}_{\ast}$, ${\bm{\Lambda}}_{G}$), respectively. The joint posterior predictive distribution is given by

[TABLE]

where

[TABLE]

While we will soon demonstrate the many advantages of this Gaussian process formalism for disentangling spectra, it also inherits a few of the limitations common to most disentangling techniques. First, there is the somewhat obvious limitation that one needs to observe the binary system at multiple orbital phases to obtain different morphologies of composite spectra. If the orbital period of the system is very long, or by chance all observations were obtained at or near the same orbital phase, we would not be able to separate the composite spectra into individual components. Second, because the intrinsic spectra are modeled nonparametrically, without reference to physical models of spectra, reconstruction techniques are unable to say anything about the stellar velocities in an absolute sense. This means that from spectroscopic reconstruction alone, we are unable to constrain the systemic velocity of the system, $\gamma$, and so this parameter remains fixed to zero in our analysis. Once the spectra are reconstructed, however, $\gamma$ can be easily recovered by cross-correlation with physical spectral templates. Of related consequence is that the reconstruction technique is only able to constrain the relative velocities of each star from epoch to epoch; it is not possible to measure the velocity of star B relative to star A at any epoch with the reconstruction technique alone. However this can be later measured easily by correlation with physical spectra models, or if the two stars in question are bound, inferred through their Keplerian orbits.

2.3 Application to mock data

To explore the performance of the Gaussian process spectral model, we generate mock spectroscopic binary datasets from a fiducial orbit and assess the accuracy of the recovered stellar spectra and orbital parameters. To generate the datasets, we use segments of real data of the K5 star LkCa 14 observed with the CfA/TRES spectrograph on Mt. Hopkins (obtained for a separate purpose under a program with P.I. I. Czekala). Multiple epochs of data were combined and smoothed to create a high signal-to-noise template spectrum. We chose two separate regions of this template with an average spectral line density ($5235\AA$ to $5285\AA$ for the primary and $6420\AA$ to $6480\AA$ for the secondary) to mimic spectra from two different stars, and then assigned the secondary fluxes to the same wavelength range as the primary. We assume a fiducial binary orbit with $\theta=\{K_{\mathrm{A}}=5.0\,\mathrm{km\;s}^{-1},\,q=0.2,\,e=0.2,\,\omega=10.0^{\circ},\,P=10.0\,\mathrm{days},\,T_{0}=0.0\,\mathrm{JD},\,\gamma=5.0\,\mathrm{km\;s}^{-1}\}$, sampled at 10 epochs over the course of a $\sim 25\,$day period (see Figure 2.3, left panel). To generate a composite spectrum for a single epoch, we first Doppler shift each template by its calculated radial velocity, interpolate the secondary spectrum to the same wavelength vector as the primary, multiply each component spectrum by its chosen fractional flux contribution, truncate each template to the range $5265\AA$ to $5275\AA$, sum the spectra, and then add random Gaussian noise according to a chosen signal-to-noise ratio. This process is illustrated for the first three epochs in the left panel of Figure 2.3. While this mock dataset is a small fraction ($<1/500$) of the temporal and wavelength coverage of a typical radial velocity dataset, it will serve to illustrate the salient characteristics of our technique in a compact and expedient manner.

We first assess the performance of the technique on a mock dataset constructed with moderate flux ratio components and medium-high signal-to-noise “observations” ($f_{B}/f_{A}=0.2$, SNR = 60 per resolution element; linespread function FWHM = 2.5 pixels), and then explore the robustness of the technique as these quantities are varied to more extreme values. The posterior distribution of the orbital parameters and Gaussian process hyperparameters is explored simultaneously using a simple implementation of Metropolis-Hastings Markov Chain Monte Carlo (MCMC). The inferred family of orbits is shown in the right panel of Figure 2.3, and these demonstrate good agreement with the input radial velocities. Several draws of reconstructed component spectra and the predicted composite spectrum are shown in the right panel of Figure 2.3. To explore the full distribution of inferred component spectra, for each sample we first randomly draw orbital parameters and Gaussian process hyperparameters from the posterior distribution, then draw a sample from the Gaussian process posterior predictive function (Equation 46).

A common limitation of most spectral disentangling techniques is the inability to constrain the continuum level of each star without comparison to a physical model, which our technique suffers from as well. This means that while the relative amplitudes of the variation in each spectrum are well-constrained, there is no knowledge about which star is brighter, i.e., there is a degeneracy between a very luminous companion star with shallow lines and a very faint companion with very deep spectral lines (Bagnuolo et al., 1992). As with the $\gamma$ degeneracy, the flux ratio degeneracy can be easily broken by later comparison with a flux-calibrated, physical stellar model (Pavlovski & Hensberge, 2010). In the right columns of Figure 2.3, we have vertically offset each recovered spectrum by an arbitrary constant to look similar to the left panels, which has the potential to be misleading. Therefore in the bottom panel of Figure 2.3, we show the draws of the reconstructed spectra from mean zero Gaussian processes (${\bm{\mu}}_{f}={\bm{\mu}}_{g}={\bm{0}}$) to emphasize that they contain no information about the continuum level. Rather, the input spectra have been subtracted by an arbitrary constant to match the level of the reconstructed spectra. The reconstructed spectra do an excellent job of reproducing the shape of the input spectra to a precision far surpassing the average pixel noise of the mock dataset. Interestingly, there is a somewhat appreciable scatter in the vertical offset level of the Gaussian process draws. This is because ${\bm{f}}_{\ast}$ and ${\bm{g}}_{\ast}$ are drawn from a joint posterior predictive distribution while only their sum is constrained by the data, allowing the mean levels to trade off somewhat (i.e., if a primary draw is slightly higher, then the secondary draw is slightly lower).

Next, we pursue a grid of tests to explore the performance of the Gaussian process model on mock data with a range of flux ratios $f_{B}/f_{A}=\{0.05,0.2,0.5\}$ and signal-to-noise SNR $=\{20,60,120\}$. The results of the reconstruction are shown in Figure 2.3, and a visual representation of the one-dimensional marginal posterior distributions is shown in Figure 2.3. Because of the aforementioned flux ratio ambiguity, we additionally describe our suite of mock data by the standard deviation of each spectrum ($\sigma_{A}$ and $\sigma_{B}$), which is an invariant quantity as concerns the performance of the technique.

For some binary systems (e.g., eclipsing systems), it may be possible to precisely constrain the orbital parameters via alternate means (e.g., photometric eclipses) to a higher precision than is otherwise possible with a noisey and/or extreme flux ratio spectroscopic dataset alone. Therefore in Figure 2.3 we show draws of the inferred spectra marginalized over the full orbital and Gaussian process hyperparameters (pale grey lines, background) as well as the inferred spectra extracted at the fiducial orbital values but varied over the hyperparameters (hued lines, foreground). The input spectra are shown in black, vertically offset by an arbitrary constant to match the Gaussian process draws as before. To reduce the aforementioned vertical scatter of random draws in this plot, we pin the flux level of primary spectrum to be zero at $\lambda=5268\AA$. The average noise per pixel in the dataset is shown by an error bar in the right of each panel.

Unsurprisingly, the reconstruction technique performs best for high SNR datasets with moderate (near-equal) flux ratios. The primary spectrum is always accurately inferred across our grid of tests. More importantly, the technique can infer secondary spectra for moderately high flux ratios with very noisy data ($f_{B}/f_{A}=0.2$, SNR = 20) and extreme flux ratios with moderate quality data ($f_{B}/f_{A}=0.05$, SNR = 60), albeit with substantial uncertainty. The inferred values of $a_{f}$ and $a_{g}$ approximate the values of $\sigma_{A}$ and $\sigma_{B}$, while $l_{f}$ and $l_{g}$ are a crude measure of the average width of the spectroscopic lines. For the noisiest, most extreme flux ratio test ($f_{B}/f_{A}=0.05$, SNR = 20), the secondary Gaussian process is unable to constrain the shape of the secondary other than showing that the amplitude of its variations must be below a certain level (see Figure 2.3). In this situation, there also emerges a degenerate solution with large $a_{g}$ and $l_{g}$: the data spectrum is dominated by the primary star modulo the addition of an approximately flat continuum. This degenerate solution could be avoided by the addition of a prior that favors smaller $l_{g}$ at the expected width of the secondary lines. Overall, the results from this grid of tests are extremely encouraging and demonstrate that the Gaussian process framework is able to accurately infer composite spectra to a precision far exceeding the average per-pixel noise in a spectroscopic dataset. It is important to emphasize that in these tests we have only been using a single $10\AA$ chunk of the spectrum; any analyses incorporating a wider swath of spectrum will potentially be much more sensitive.

2.4 Single-lined Spectroscopic Binaries and Exoplanet Search

In the limiting case that the flux of the secondary becomes negligible, the system reverts to a single-lined spectroscopic binary and any constraint on the mass ratio ($q$) vanishes, leaving only 5 orbital parameters. If there is a strong belief that light from the primary component dominates the composite spectrum, then one could save computational time and explore a lower-dimensional parameter space with ${\bm{\Sigma}}_{G}={\bm{0}}$. However, if there are indeed unrecognized signatures of the secondary in the composite dataset, a constrained single-lined spectroscopic model will deliver biased orbital parameters. When using the full complement of binary orbital parameters, the Gaussian process framework provides the ability to place constraints on the relative flux ratio of an unknown companion in a probabilistic setting. While other search techniques rely upon cross-correlation and subtraction of a presumed template (e.g., Gullikson et al., 2015), the Gaussian process framework provides a sensitive, model-free mechanism to detect or place limits on faint spectroscopic lines from a companion star. Of course, secondary contamination is much less an issue in the case of radial velocity measurements of exoplanetary systems around single-stars. While any constraints on the nature of the secondary spectrum might be meaningless, such a Gaussian process technique could be used to infer precision radial velocities for planet search. If a statistically significant non-zero semi-amplitude ($K_{A}$) and orbital period ($P$) were found for the primary star, this would signal the existence of an orbiting exoplanet. Lastly, the framework could be used to simply infer a high-fidelity template of the primary star, which could be used for abundance analysis or as a template for traditional cross-correlation radial velocity measurements.

In the development of this framework, we have assumed that the intrinsic spectrum of the star is constant with time, an assumption also made by planet search pipelines using the standard cross-correlation technique. However, subtle changes in the spectrum due to starspots can bias the inferred radial velocity at the level of $\Delta v\sim 10\,\mathrm{m\,s}^{-1}$ (Dumusque et al., 2011), dwarfing the operating precision of state of the art spectrographs ($\Delta v\approx 1\,\mathrm{m\,s}^{-1}$; Fischer et al., 2016). In the Gaussian process framework, the effects of starspots could be minimized by incorporating an additional time-covariance that allows for flexibility in the exact shape of the spectrum. Such an extension would be especially useful when simultaneous time-series photometric observations bracket the spectroscopic observations. Measurements of starspot modulation could be used to predict the changes in the spectral line profiles, as is currently done in advanced radial velocity analyses (e.g., Haywood et al., 2014). We further explore the consequences of time variability in Section 4.2.

3 Results: Application to the Mid-M dwarf binary LP661-13

3.1 The LP661-13  dataset and traditional orbit analysis

As part of the MEarth-South survey (Irwin et al., 2015), Dittmann et al. (2016) recently discovered the mid-M dwarf LP661-13 to be a short-period eclipsing binary. Subsequent photometric and radial velocity follow-up over two observing seasons constrained the system period to be $P=4.7043512_{-0.0000010}^{+0.0000013}\,\mathrm{days}$, the primary component mass and radius to be $M_{A}=0.30795\pm 0.00084\,M_{\odot}$ and $R_{A}=0.3226\pm 0.0033\,R_{\odot}$, the secondary component mass and radius to be $M_{B}=0.19400\pm 0.00034\,M_{\odot}$ and $R_{B}=0.2174\pm 0.0023\,R_{\odot}$, and the distance to be $24.9\pm 1.3\,\mathrm{pc}$. Given the lack of strong X-ray emission and the absence of Lithium in the spectra, LP661-13 is likely a field-age system. At these masses, both stars are fully-convective ($M<0.35M_{\odot}$; Chabrier & Baraffe, 1997). Such a system provides an interesting application for our spectroscopic disentangling technique because the orbital parameters are known to high accuracy due to the eclipsing nature of the system. For the purposes of demonstrating the disentangling framework, we ignore any eclipse constraints and utilize only the spectroscopic dataset in Dittmann et al. (2016). Later, we utilize the additional constraints from the eclipses to discuss the fundamental properties of LP661-13 A and B, in particular the stellar radii.

The spectroscopic dataset consists of 14 epochs (1 hour integration each, except for one 40-minute epoch) taken over the course of February 2014 to January 2015 with the Tillinghast Reflector Echelle Spectrograph (TRES; $R=\lambda/\Delta\lambda=44,000$) on the 1.5m Tillinghast reflector at the Fred Lawrence Whipple Observatory on Mt. Hopkins, in Arizona. The data were reduced by the standard TRES pipeline (Buchhave et al., 2010), and placed in the barycentric frame. At red wavelengths, the signal to noise of the spectra ranges from 20 to 60 per resolution element, with a median value of 43. To determine radial velocities with TODCOR (Zucker & Mazeh, 1994), Dittmann et al. (2016) also obtained high signal-to-noise observations of the following mid-M template stars from Kirkpatrick et al. (1991): Gl 273 (M3.5), Gl 699 (Barnard’s star; M4), Gl 83.1 (M4.5), and Gl 406 (M6), which we later use for comparison with the disentangled spectra of LP661-13 A and B. Dittmann et al. (2016) found the strongest signal when using the Gl 699 template over the wavelengths $7060\AA-7200\AA$, which exhibit sharp molecular features, and measured a best-fit component light-ratio of $L_{B}/L_{A}=0.434\pm 0.025$. Because one or both stars may be spotted and the fraction of eclipsed spots is unknown, Dittmann et al. (2016) did not quote a measurement from the eclipse measurements. If one relies upon the derived physical parameters for the system, the temperature scale of Mann et al. (2015), and the Allard et al. (2001) models, the predicted light ratio would be $L_{B}/L_{A}=0.406$. Dittmann et al. (2016) omitted the observation on HJD 2456743 from their orbital fit because they were unable to separate the primary and secondary peaks of the CCF, meaning that at this epoch each component had a similar radial velocity. Assuming equal measurement uncertainty for the primary and secondary radial velocities, they calculated a best-fit binary orbit with the 13 epochs of measurements and found the average residuals of the primary and secondary velocities to be $0.05\,\mathrm{km\;s}^{-1}$ and $0.20\,\mathrm{km\;s}^{-1}$, respectively.

In order to make a direct comparison with the orbital parameters we will determine with the Gaussian process framework, we re-fit the 13 radial velocities in Dittmann et al. (2016) using the orbital parameterization described in §2.2. We adopt $0.05\,\mathrm{km\;s}^{-1}$ and $0.20\,\mathrm{km\;s}^{-1}$ as the uncertainty on each radial velocity measurement for the primary and secondary, respectively, and employ a simple $\chi^{2}$ likelihood function with flat prior distributions. We explore the $7$-dimensional posterior with the emcee implementation (Foreman-Mackey et al., 2013) of the Affine Invariant Ensemble Sampler (Goodman & Weare, 2010). To assess convergence, we run multiple ensembles of 28 walkers for 10,000 iterations each, starting from different locations in parameter space and burning 5,000 iterations worth of samples. The marginalized, one-dimensional posteriors are shown in Figure 3.1 and are in excellent agreement with Dittmann et al. (2016), except for our determination of $\gamma=0.025\pm 0.011\,\mathrm{km\;s}^{-1}$. Our measurement is discrepant with the value reported by Dittmann et al. (2016) ($-0.009\pm 0.014\,\mathrm{km\;s}^{-1}$) at the $2\sigma$ level, which we suspect may be due differences in how the radial velocity uncertainties were incorporated in the fitting process. As noted in Dittmann et al. (2016), the orbit is extremely circular, which leads to the common $\omega-T_{0}$ degeneracy.

In addition to the standard spectral reduction performed by Dittmann et al. (2016), we perform additional processing of the spectra to facilitate a uniform comparison between objects. Along with the observations of LP661-13 and the 4 aforementioned M-dwarf template stars, we also download an observation of Vega taken on HJD 2456740 from the TRES archive (P.I. D. Latham) to use as a reference for telluric line contamination. As part of a typical TRES observational sequence, observations of a continuum source are taken to characterize the echelle blaze function, which has been shown to be very stable with time for most echelle orders. We divide each spectrum by its associated blaze measurement to produce approximately flat spectra. This process does not directly lead to flux-calibrated spectra, however, since the absolute brightness of the source is not known. Nor does this process necessarily lead to spectra that match the same shape as flux-calibrated spectra, since the spectrum of the continuum source may be significantly different than that of a source originating from above the atmosphere. A TRES calibration program (P.I. I. Czekala) repeatedly observed several spectrophotometric standard stars over several months and demonstrated that the TRES intra-order throughput shape is very stable with time. Since the bandpass sensitivity in each order is well-described by a 4th-order Chebyshev polynomial, this means that from night to night the linear and higher order coefficients varied less than $4\%$. The inter-order sensitivity is less stable, with nightly scatter typically $\lesssim 10\%$.

Large scale, sloping features of M-dwarf spectra are important for spectral typing and characterization. Although we do not have observations of spectrophotometric standards contemporaneous with the LP661-13 observations, because the intra-order TRES throughput is stable we can approximate flux-calibrated LP661-13 spectra by applying the sensitivity functions derived from the calibration program, which were measured from a series of observations of the O2 spectrophotometric standard BD+28 4211 in November 2013. Because the inter-order flux calibrations are not as stable as the intra-order calibrations, we further rectify the spectra so that the median flux of each order is 1. The result of this process is spectra whose intra-order shape calibration should be comparable to a bona fide flux calibration within $4\%$. Because of the per-order rectification process, the inter-order shape calibration will not be correct. It should be noted that if same-night spectrophotometric flux observations existed, the proper shape of the entire spectrum could be recovered.

3.2 Gaussian process inference of orbital parameters

and the disentangled spectra of A and B

We apply our spectroscopic binary Gaussian process framework to the pseudo-flux-calibrated dataset of LP661-13 in a similar manner as in §2.3. For computational efficiency, we focus on two high signal to noise regions of the spectrum uncontaminated by telluric lines, $6650\AA-6770\AA$ and $7060\AA-7140\AA$. Unlike with traditional cross-correlation radial velocity analysis, composite spectra with near-equal component radial velocities still provide useful information in our Gaussian process framework and help constrain the underlying intrinsic component spectra. Therefore we use all 14 epochs of spectroscopy in our analysis, including the epoch omitted by Dittmann et al. (2016). We divide the dataset into 20 chunks of $10\AA$ each so that the evaluation of the posterior function can be parallelized across a computer cluster. The posterior distribution for each $10\AA$ chunk is evaluated on an individual core, and then gathered and multiplied together to yield a final posterior for the full spectral range. This final posterior is explored using a simple Metropolis-Hastings MCMC,222included in emcee. whose jump covariances have been tuned to the structure of the posterior for efficient exploration of the parameter space. We run multiple chains from different starting locations for 20,000 iterations each and assess their convergence by the Gelman-Rubin statistic $\hat{R}$ (Gelman et al., 2013). The one-dimensional marginalized posteriors are shown in Figure 3.1 in blue. We recover similar posteriors as the traditional cross-correlation analysis, and notably several parameters are obtained with higher precision (e.g., $q$, $K_{A}$, $e$, and $P$). With more computational resources, a larger spectral range could be included in the analysis and potentially tighten the parameter constraints even further. In Figure 3.2 we further illustrate the increased precision of the Gaussian process based technique by demonstrating that it delivers a smaller radial velocity scatter at each observational epoch compared to the traditional cross-correlation based technique.

An example of the analysis for a single chunk of LP661-13 spectrum is shown in Figure 3.2, focusing on a particularly dramatic bandhead for emphasis. The composite dataset for each epoch is shown in black, with the corresponding orbital phase labeled in the upper right. Draws of the Gaussian process are generated while varying over random draws of the orbital parameters and hyperparameters from the posterior probability distribution, showing that the distribution of reconstructed component spectra fit well within the noise. As noted in §2.3, there is the additional question of choosing the normalization level for each reconstructed spectrum. And so in Figure 3.2, we arbitrarily choose ${\bm{\mu}}_{f}=0.65$ and ${\bm{\mu}}_{g}=0.35$ to display these spectra (denoted by horizontal dashed lines). We also show the residuals computed from the mean prediction of the composite spectrum at each epoch. Because the spectra have a natural covariance with wavelength, any individual realization of residuals may show a slight covariance as well.

While the exploration of the posterior is performed in $10\AA$ chunks and using only a subset of the full spectral range, by assuming the inferred orbital parameters and Gaussian process hyperparameters are valid outside of this range, one can reconstruct a wide range of component spectra. In Figures 3.2 and 3.2 we show the reconstructed spectra over the full wavelength range $6800\AA-8040\AA$, which includes many of the notable spectral regions used for classifying M-dwarfs. We also overplot an observation of Vega, a fast rotating A0 star with a near-featureless continuum, in order to highlight regions affected by telluric lines. Any spectral features intrinsic to Vega (e.g., O I $\lambda\lambda\lambda 7771$) will be broadened by its fast rotational velocity ($v\sin i=21.9\,\mathrm{km\;s}^{-1}$; Hill et al., 2004), and contrast with the narrower telluric lines. In contrast to the carefully chosen narrow range we used for the orbital parameter inference, this expanded region includes numerous telluric absorption lines and night sky emission lines. Because these features are not consistent with the rest-frame of either star A or B, they act to distort the reconstruction process by acting as outlier flux points that cannot be reproduced by the model. In principle, it should be possible to model the effects of telluric lines as an additional Gaussian process since their location is known precisely, although their absorption depth changes from epoch to epoch with airmass and atmospheric conditions. Simultaneous modeling of the telluric features would allow us to access spectral information of the stars in regions which are currently contaminated. We note that other spectral disentangling methods have successfully modeled telluric lines (Hadrava, 2006), but we leave a detailed implementation of these concepts for future work.

Typically, night sky emission lines are subtracted by the TRES reduction pipeline, although in some instances the subtraction is sub-optimal and leaves residual flux. We denote regions affected by poor-night sky line subtraction by a $\circ$ symbol and mask these features in the dataset (Osterbrock et al., 1996). Because $\sim 5\AA$ windows surrounding these lines have been masked for all epochs, there are regions for which $f$ and $g$ are minimally constrained. In these regions, the Gaussian process mean is generally a flat line, and the draws mimic draws from the prior distribution. Were we to predict intrinsic spectra outside of the range of the detector (e.g., in the coverage gaps between echelle orders), a similar phenomenon would occur.

To properly plot the reconstructed LP661-13 A and B spectra (i.e., draws of the mean ${\bm{0}}$ Gaussian processes conditional on the composite data) in a format suitable for comparison with other rectified observations of mid-M template stars, we must first renormalize the spectra. Renormalization refers to the process whereby the spectra are brought to their correct continuum level in the frame of the data with the addition of a constant, and then rectified to a mean of 1 by dividing by that same constant

[TABLE]

For LP661-13 A and B, we assume that the mean vectors are simply flat with wavelength, but in principle their relative amplitude could be guided by beliefs about how the flux ratio changes with wavelength, (i.e., ${\bm{\mu}}_{f}(\lambda)$, and ${\bm{\mu}}_{g}(\lambda)$). An incorrect renormalization will distort the continuum level and the contrast of the spectral lines. Because the secondary spectra are intrinsically fainter, they must be scaled by a larger constant to be rectified, which has the effect of making the rectified spectrum appear noisier. We experiment varying the flux ratio until we find renormalizations of LP661-13 A and B that match the amplitude and slope of the nearby mid-M templates. We find that a value of ${\bm{\mu}}_{g}/{\bm{\mu}}_{f}=0.33$ yields a renormalization that has the spectra of LP661-13 A and B closely matching the mid-M templates over the full wavelength range. We note that flux ratios as high as 0.434 (Dittmann et al., 2016) yield acceptable results over much of the spectral range, however the contrast of the bandhead at $7090\AA$ appears too shallow in the secondary spectrum. Frémat et al. (2005) discuss an alternate renormalization strategy suitable for stars with very deep absorption lines, which must remain positive even in their deep cores, and so they set a lower limit on the continuum level of the spectrum. Unfortunately, no such lines exist in the LP661-13 spectra.

3.3 LP661-13 A and B in the mid M-dwarf context

Now, with the disentangled spectra of LP661-13 A and B known (up to a modest uncertainty in the renormalization factor), we can examine the stars in the context of the other mid-M spectral templates. Although we have assumed that the flux ratio between A and B is constant across the full wavelength range, given that both stars are evidently very close in spectral type this appears to have been a reasonable decision. Next to LP661-13 A and B, we also plot TRES observations of mid-M standard stars (Kirkpatrick et al., 1991) and highlight several spectral regions commonly used in classification, which generally include molecular features. With low to moderate resolution spectra, the flux within these bands measured relative to a nearby pseudo-continuum region serves as an index for objective spectral typing. Some of the most widely indices used are summarized by Lépine et al. (2013): CaH2, CaH3, and TiO5 (Reid et al., 1995); VO1 (Hawley et al., 2002); and TiO6 and VO2 (Lépine et al., 2003). Those most sensitive (i.e., changing most rapidly) for M3 - M6 stars are the TiO5, TiO6, and VO2 indices (Lépine et al., 2013).

We would like to directly measure the indices for the disentangled spectra of LP661-13 A and B and report an index-based spectral type for these stars; however, when we measured the indices of the template stars in common with Lépine et al. (2013), we were unable to reproduce their values. Since we were able to recover the general trends with spectral type, we suspect that this discrepancy is likely due to the higher resolution of our instrument and the fact that our spectra are not flux-calibrated following a rigorous procedure. Therefore instead of comparing indices, we directly compare the morphology of the spectra within the noted index regions. We find that LP661-13 A and B naturally fall in a sequence of mid-M standards near Gl 699 (M4) and Gl83.1 (M4.5), which agrees well with the previously determined composite spectral types of M3.5 (optical; Reid et al., 2004) and M4.27 (NIR; Terrien et al., 2012).

We emphasize that other than the procedure performed in Equation 48 for LP661-13 A and B, no post-processing (e.g., continuum or pseudo-continuum normalization) has been applied to any of the spectra presented in Figures 3.2 and 3.2. If a proper flux-calibration of the data were performed on the raw data and the flux ratio correctly specified, then the inferred shape of the reconstructed spectra would be correct. The ability of the Gaussian process technique to constrain the low order shape of inferred spectra avoids a potential downside of Fourier disentangling techniques.

The Gaussian process technique can also disentangle emission lines, and we find that both LP661-13 A and B show strong Balmer series emission—in fact, the radial velocity semi-amplitude is large enough such that the emission lines from A and B are completely separated at quadrature. The strong Balmer emission indicates that both stars are magnetically and chromospherically active (Hawley et al., 1996; West et al., 2011; Lépine et al., 2013), possibly as a result of rapid, synchronous rotation due to tidal locking (Dittmann et al., 2016; Newton et al., 2017). Because magnetic activity serves to trap flux and inflate stellar radii while decreasing effective temperature (Morales et al., 2009, 2010), it may provide an explanation for the inflated radii of LP661-13 A and B measured by Dittmann et al. (2016).

4 Discussion

4.1 Hierarchical inference of spectral chunks

In our spectroscopic binary Gaussian process framework, we made the assumption that the Gaussian process kernels were the same functional form for each spectrum, and that each $10\AA$ chunk of spectrum used the same Gaussian process hyperparameters. If the flux ratios of the spectra are approximately constant over the chosen wavelength range, and the relative amplitudes of the spectral line variations are similar, as in LP661-13, then this approach is adequate. However, if we are inferring a large region of the spectrum which might contain spectral features of very different morphologies, then a more accurate spectral extraction might be obtained by modeling the Gaussian process hierarchically, where the amplitude and length scale for an individual chunk is drawn from a general hyperparameter distribution. For example, the amplitude ratio could be guided by the ratio of fluxes expected from the hypothesized effective temperatures, with the amplitude of each individual chunk allowed to vary somewhat around this mean value. Recently, high contrast observational techniques have shown great promise to detect the signatures of exoplanet atmosphere spectra via cross-correlation with templates (e.g., Snellen et al., 2014; Birkby et al., 2017)–with increasing telescope sensitivity and an appropriately flexible Gaussian process model it may become possible to reconstruct the exoplanet spectra directly.

Although we have only dealt with a single dataset acquired by the same instrument in an identical setup, in principle, data from multiple instruments covering different wavelength ranges (e.g., optical and infrared spectra) could be inferred simultaneously if the Gaussian process hyperparameters were allowed to vary individually for each instrument, tracking changes in the line spread function. Moreover, since the convolution of a Gaussian process with a line spread function is also a Gaussian process, it may even be possible to simultaneously model data taken with different instrumental setups but covering the same wavelength range if the latent spectrum is modeled at high resolution and convolved down to the resolution of each observation.

4.2 Extensions to time-variable spectra

Thus far, we have only discussed covariances with a $\lambda$-dependence. If one or more of the components is believed to be time-variable, for example due to chromospheric activity, then we could introduce an additional $t$-dependence to the covariance kernel as well. Until now, we have essentially enforced that the Gaussian process has an infinite covariance with time: pixels with nearby wavelengths contribute equally to constraining the stellar flux, regardless of measurement epoch. If we had an additional time-variable component, then we could allow subsequent epochs to vary such that only flux measurements that were close in wavelength and close in time would constrain the inferred spectrum. This also means that we would be able to make predictions (realizations) of the stellar spectrum not only as a function of wavelength, but also as a function of time. Such a mechanism could provide radial velocity measurements for the noiseiest stars, where stellar activity has hampered analysis techniques more profitably applied to quiet stars. In particular, probing the youngest exoplanet population has been hindered by stellar activity, requiring extensive analysis to prove that a planetary signal is not simply unexplained variability (e.g., Johns-Krull et al., 2016). For more mainstream exoplanet radial velocity analysis, a time-variable, $\lambda$-specific kernel could provide a natural way to identify or downweight spectral regions which bias a traditional radial velocity signal. It may also be possible to adapt the Gaussian process framework to model the iodine absorption lines commonly imprinted on a stellar spectrum for use in precision radial velocity analysis.

While we have focused on inferring stellar photospheres, the spectrum to be inferred need not necessarily be stellar in origin. A potentially exciting application is to use a Gaussian process to infer active signs of accretion in young T Tauri stars. While the choice of squared exponential kernel has worked well for modeling the intrinsic stellar spectrum of LP661-13, for other types of stars or variable phenomena, different kernels might be more effective. For example, using the Matérn kernel for the $t$-correlation might provide a Gaussian process that more closely reflects the stochastic accretion process. Such a modeling approach could also be combined with a mean function ${\bm{\mu}}_{f}(\lambda)$ specified by flexible physics-based models (e.g., Czekala et al., 2015; Gully-Santiago et al., 2017) in order to simultaneously determine the fundamental properties of the veiled star, a traditionally difficult process (Herczeg & Hillenbrand, 2014).

5 Summary

We have demonstrated the successful application of Gaussian processes for inferring component spectra of single and double-lined spectroscopic binaries, while simultaneously exploring the posteriors of the orbital parameters and the spectra themselves. We evaluated the performance of the technique by using mock data constructed over a range of flux and signal to noise ratios, and demonstrated that the technique has desirable characteristics for a wide range of data qualities, returning unbiased posterior estimates. In contrast to Fourier-domain methods, our technique excels at reconstructing the low-order shape of the component spectra and provides a natural and well-motivated probabilistic formalism to model the data in its natural pixel space. While other $\lambda$-based techniques do exist (e.g., Simon & Sturm, 1994), our Gaussian process formalism can account for the natural $\lambda$-covariances in each spectrum, providing a natural “de-noising” of the spectra typically offered by Fourier techniques.

We applied our framework to the mid-M eclipsing binary LP661-13. Using only a limited number of observations with modest SNR, we simultaneously recovered both the A and B components and inferred the orbital parameters with a precision exceeding that of standard cross-correlation techniques. For the first time, we independently analyzed each component of LP661-13 in the context of a mid-M spectral sequence, and determined that both components are M4 systems, with A appearing slightly earlier than B. We discussed potential applications of Gaussian processes for spectral modeling, including precision exoplanet radial velocity work around active young stars. The code used in our analysis is available as a freely available open-source package.333Precision Spectroscopic Orbits A-Parametrically (PSOAP) available at https://github.com/iancze/PSOAP

IC would like to thank the following people: Guillermo Torres and Dave Latham for advice and helpful discussions about fundamental stellar properties, stellar orbits, and radial velocity techniques; Kevin Gullikson for guidance about searching time series spectra for faint companions; Hyungsuk Tak for discussions about Quasar light curve analysis with Gaussian processes; Kevin Hardegree-Ullman for references to M-dwarf spectral indices; Eric Nielsen for discussions about orbital dynamics; Greg Green and Dan Foreman-Mackey for general discussions of Gaussian processes; Charlie Conroy and the hosts of the “Fitting Stars, CMDs, and Galaxies” workshop for insightful discussion about stellar inference than influenced this work; and Allyson Bieryla, Gil Escuerdo, and Jessica Mink for their service in support of the TRES instrument. KM is supported at Harvard by NSF grants AST-1211196 and AST-156854. Work by BTM was performed under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. This material was based upon work partially supported by the National Science Foundation under Grant DMS-1127914 to the Statistical and Applied Mathematical Sciences Institute. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.