Bayesian Methods for Joint Exoplanet Transit Detection and Systematic Noise Characterization

A standard model for systematic noise is a linear reduced basis model (Kov'acs et al. 2005) in which systematic noise vectors ${{\boldsymbol{l}}}_{i}$ are modeled as a linear combination of a set of K systematic noise basis vectors {${{\boldsymbol{v}}}_{k}$} weighted by coefficients c_i^k. Because the set of systematic noise signals {${{\boldsymbol{l}}}_{i}$} are reasonably expected to be correlated between light curves of sufficiently common instrumental origin (Stumpe et al. 2012), they may be expressed in terms of a common set of systematic noise basis vectors {${{\boldsymbol{v}}}_{k}$},

$\begin{eqnarray}&&{{\boldsymbol{l}}}_{i}=\sum _{k=1}^{K}{c}_{i}^{k}{{\boldsymbol{v}}}_{k}.\end{eqnarray} \tag{ 4 }$

Systematic basis noise vectors {${{\boldsymbol{v}}}_{k}$} may be estimated via dimensionality reduction techniques (Cunningham & Ghahramani 2015) applied on {${{\boldsymbol{y}}}_{i}$}. PCA is a commonly used technique for cotrending and produces a set of orthogonal basis vectors (Mazeh et al. 2007; Twicken et al. 2010; Petigura & Marcy 2012; Smith et al. 2012; Stumpe et al. 2012; Foreman-Mackey et al. 2015). A non-Bayesian least-squares estimation of the coefficients ${{\boldsymbol{c}}}_{i}$ is equivalent to maximum likelihood estimation (Kay 1999a) with an assumed Gaussian stellar noise model ${{\boldsymbol{s}}}_{i}$, where the samples are independent and identically distributed. However, least-squares estimation, by minimizing the net rms error, is prone to overfitting residual transit signatures (Smith et al. 2012; Stumpe et al. 2012). As shown by Smith et al. (2012), a Bayesian estimation of ${{\boldsymbol{c}}}_{i}$ can mitigate this issue by allowing the incorporation of a prior p(${{\boldsymbol{c}}}_{i}$) on the coefficients. The prior here is conditioned on the index i ∈ I; this index accommodates latent variables, such as position within the CCD and the stellar magnitude, that produce clustering in systematic noise properties as noted above (Smith et al. 2012; Stumpe et al. 2012).

By definition, the cotrending model for systematic noise in Equation (4) does not fully capture systematic noise that is temporally uncorrelated between light curves. Examples of such systematic noise effects include cosmic-ray events, sudden pixel sensitivity drop-off, and electronic image artifacts (Jenkins et al. 2010a; Stumpe et al. 2012). While our model cannot account for outlier effects, a number of residual unmodeled signatures may be treated as effectively Gaussian and absorbed into the term ${{\boldsymbol{s}}}_{i}$ in Equation (2).

In the current work, we assume that the stellar noise ${{\boldsymbol{s}}}_{i}$ (in each light curve i) may be reasonably modeled by a Gaussian noise model ${{\boldsymbol{s}}}_{i}\sim { \mathcal N }({\boldsymbol{0}},{\mathrm{Cov}}_{s,i})$. Gaussian processes are reviewed in the monograph by Gallager (2013, ch. 3). Gaussian models are powerful nonparametric models of processes with complex or unknown generating phenomena and may be justified by the central limit theorem (Grinstead & Snell 2012). Gaussian models have been shown to be effective for modeling stellar variability (Pereira et al. 2019) and have been used in the context of exoplanet detection (Carter & Winn 2009; Rajpaul et al. 2015).

Stellar flux density time series display variability on a wide range of timescales (Conroy et al. 2018), and time-correlated fluctuations significantly affect the detectability of transit signals (Borucki et al. 1985; Jenkins 2002; Pont et al. 2006). A Gaussian noise model can incorporate time-correlated structure via its covariance matrix. A review of common correlated noise estimators for exoplanet light curves is provided by Cubillos et al. (2017).

Our model further assumes that stellar noise is stationary within a fixed time window. In a stationary Gaussian colored noise model,⁶ correlations only depend on the separation between two points in time. Equivalently, the values on the stellar covariance matrix ${\mathrm{Cov}}_{s,i}$ are constant along the diagonals, and the matrix is Toeplitz in form. Such a model admits a power spectral representation and allows spectral analysis (Kay 1999a), which is considerably more efficient than time-domain analysis.

Stellar processes do, however, evolve in time and can exhibit nonstationarity as exemplified by our Sun (Borucki et al. 1985; Jenkins 2002). For this reason, we confine our stellar noise model to be stationary only within a single observational quarter of Kepler data (approximately 90 days). This assumption is motivated in part by the analysis of common stellar periodic variability in Kepler data on timescales of days to weeks (Basri et al. 2010).

The extension of our detectors to multiple Kepler observing quarters and associated implications for our statistical assumptions regarding signal covariance timescales are described in Appendix B.

A transit signal is typically approximated as a periodic box function (Kovacs et al. 2002). The periodic box function is parameterized as follows: α describes the fractional drop in light relative to the stellar signal, P is the orbital period of the planet, d is the transit duration, n is the time index, and t₀ is the phase (epoch) of the signal:

$\begin{eqnarray}{t}_{\alpha ,P,{t}_{0},d}[n]=\left\{\begin{array}{ll}\alpha , & \mathrm{if}\ (n-{t}_{0})\,{\mathrm{mod}}_{P}\leqslant d\\ 0, & \mathrm{otherwise}\end{array}\right..\end{eqnarray} \tag{ 5 }$

A typical technique to estimate transit signal parameters is to quantize the space of possible parameter values to create a discrete set of candidate transit signals ${\boldsymbol{T}}$ and to search through this set while performing detection tests (Jenkins et al. 2002). The set of possible transit signals ${\boldsymbol{T}}$ is not described herein by a Bayesian prior as, a priori, the population statistics of exoplanets are not fully known. As noted above, other parametric forms can be used in place of a periodic box function (Mandel & Agol 2002; Seager & Mallen-Ornelas 2003).

The signal model may be expressed in matrix form. A similar form is presented in Smith et al. (2012), which we adapt to include the presence of a transit signal. As before, time is indexed by n ranging from 1 to N and i is the light-curve index. ${\boldsymbol{V}}$ is an N × K matrix with columns formed from {${{\boldsymbol{v}}}_{k}$}:

$\begin{eqnarray}&&{{\boldsymbol{y}}}_{i}={\boldsymbol{t}}+{{\boldsymbol{Vc}}}_{i}+{{\boldsymbol{s}}}_{i}\end{eqnarray} \tag{ 6 }$

$\begin{eqnarray}&&\begin{array}{rcl}\left[\begin{array}{c}{y}_{i}[1]\\ {y}_{i}[2]\\ \vdots \\ {y}_{i}[N]\end{array}\right] & = & \left[\begin{array}{c}t[1]\\ t[2]\\ \vdots \\ t[N]\end{array}\right]+\left[\begin{array}{cccc}| & | & & | \\ {{\boldsymbol{v}}}_{1} & {{\boldsymbol{v}}}_{2} & \ldots & {{\boldsymbol{v}}}_{K}\\ | & | & & | \end{array}\right]\\ & & \times \,\left[\begin{array}{c}{c}_{i}^{1}\\ {c}_{i}^{2}\\ \vdots \\ {c}_{i}^{K}\end{array}\right]+\left[\begin{array}{c}{s}_{i}[1]\\ {s}_{i}[2]\\ \vdots \\ {s}_{i}[N]\end{array}\right].\end{array}\end{eqnarray} \tag{ 7 }$

Before proceeding further, a brief overview of the matched filter is provided as it plays an essential part in our detection framework. For more detail, we refer the reader to the monographs by Kay (1993) and Poor (2013).

The matched filter is the optimal Neyman–Pearson detector for a deterministic signal ${\boldsymbol{t}}$ in Gaussian noise ${\boldsymbol{n}}\sim { \mathcal N }({\boldsymbol{0}},{\mathrm{Cov}}_{n})$ for an observed signal ${\boldsymbol{y}}$. Using the form provided by Jenkins et al. (2002), the matched filter test statistic T(${\boldsymbol{y}}$) is given by

$\begin{eqnarray}&&T({\boldsymbol{y}})=\displaystyle \frac{{{\boldsymbol{y}}}^{T}{\mathrm{Cov}}_{n}^{-1}{\boldsymbol{t}}}{\sqrt{{{\boldsymbol{t}}}^{T}{\mathrm{Cov}}_{n}^{-1}{\boldsymbol{t}}}}\underset{{H}_{0}}{\overset{{H}_{1}}{\gtrless }}\tau .\end{eqnarray} \tag{ 8 }$

This form of the matched filter may be derived from the LRT formulation in Equation (3) under the assumption of Gaussianity. The matched filter is desirable from an implementation standpoint as the distribution of the matched filter test statistic conditioned on the null hypothesis is normal: $T({\boldsymbol{y}})| {H}_{0}\sim { \mathcal N }(0,1)$. Thus, a detection threshold τ should achieve the same false-alarm rate for different noise covariance matrices ${\mathrm{Cov}}_{n}$ and signals ${\boldsymbol{t}}$. Second, if ${\boldsymbol{n}}$ is wide-sense-stationary, the matched filter may be efficiently computed in the Fourier domain (Kay 1999b).

This detector computes p_h(${{\boldsymbol{y}}}_{i}$): h ∈ { 0, 1} by marginalizing over the systematic noise prior probability p(${{\boldsymbol{c}}}_{i}$). For a particular hypothesis, the conditional probability of the observed signal ${{\boldsymbol{y}}}_{i}$ conditioned on systematic noise coefficients ${{\boldsymbol{c}}}_{i}$ is given by

$\begin{eqnarray}&&{p}_{0}({{\boldsymbol{y}}}_{i}| {{\boldsymbol{c}}}_{i})=p({{\boldsymbol{s}}}_{i}={{\boldsymbol{y}}}_{i}-{{\boldsymbol{Vc}}}_{i}| {{\boldsymbol{c}}}_{i}),\end{eqnarray} \tag{ 9 }$

$\begin{eqnarray}&&{p}_{1}({{\boldsymbol{y}}}_{i}| {{\boldsymbol{c}}}_{i})=p({{\boldsymbol{s}}}_{i}={{\boldsymbol{y}}}_{i}-{\boldsymbol{t}}-{{\boldsymbol{Vc}}}_{i}| {{\boldsymbol{c}}}_{i}).\end{eqnarray} \tag{ 10 }$

Because the ${\boldsymbol{t}}$ and ${{\boldsymbol{y}}}_{i}$ are deterministic, the only stochastic term in these distributions is ${{\boldsymbol{s}}}_{i}$, which is modeled as approximately Gaussian. See Sections 2 and 2.1.1 for a discussion of contributing non-Gaussian terms in ${{\boldsymbol{s}}}_{i}$. The LRT may be computed by marginalizing over the known prior p(${{\boldsymbol{c}}}_{i}$) for both of these conditional probability functions:

$\begin{eqnarray}&&{{ \mathcal L }}_{i}({{\boldsymbol{y}}}_{i})=\displaystyle \frac{{p}_{1}({{\boldsymbol{y}}}_{i})}{{p}_{0}({{\boldsymbol{y}}}_{i})}=\displaystyle \frac{{\int }_{{{\boldsymbol{c}}}_{i}\in {{\boldsymbol{C}}}_{i}}{p}_{1}({{\boldsymbol{y}}}_{i}| {{\boldsymbol{c}}}_{i})p({{\boldsymbol{c}}}_{i})\ d{{\boldsymbol{c}}}_{i}}{{\int }_{{{\boldsymbol{c}}}_{i}\in {{\boldsymbol{C}}}_{i}}{p}_{0}({{\boldsymbol{y}}}_{i}| {{\boldsymbol{c}}}_{i})p({{\boldsymbol{c}}}_{i})\ d{{\boldsymbol{c}}}_{i}}\end{eqnarray} \tag{ 11 }$

$\begin{eqnarray}&&=\,\displaystyle \frac{{\int }_{{{\boldsymbol{c}}}_{i}\in {{\boldsymbol{C}}}_{i}}p({{\boldsymbol{y}}}_{i}={\boldsymbol{t}}+{{\boldsymbol{Vc}}}_{i}+{{\boldsymbol{s}}}_{i}| {{\boldsymbol{c}}}_{i})p({{\boldsymbol{c}}}_{i})\ d{{\boldsymbol{c}}}_{i}}{{\int }_{{{\boldsymbol{c}}}_{i}\in {{\boldsymbol{C}}}_{i}}p({{\boldsymbol{y}}}_{i}={{\boldsymbol{Vc}}}_{i}+{{\boldsymbol{s}}}_{i}| {{\boldsymbol{c}}}_{i})p({{\boldsymbol{c}}}_{i})\ d{{\boldsymbol{c}}}_{i}},\end{eqnarray} \tag{ 12 }$

where ${{\boldsymbol{C}}}_{i}$ is the domain of ${{\boldsymbol{c}}}_{i}$.

This method obtains two fixed systematic noise estimates ${\hat{{\boldsymbol{c}}}}_{1}$ and ${\hat{{\boldsymbol{c}}}}_{0}$, with and without the presence of a transit signal, respectively. These Bayesian estimates are obtained from the posterior distributions ${p}_{h}({{\boldsymbol{c}}}_{i}| {{\boldsymbol{y}}}_{i}):{\text{}}h\in \left\{0,1\right\}$ currently as MAP estimates. These can then be used to compute p_h(${{\boldsymbol{y}}}_{i}$): h ∈ {0, 1} and the LRT. Herein lies the main distinction from prior cotrending approaches; there, a single systematic noise estimate, roughly equivalent to ${\hat{{\boldsymbol{c}}}}_{0}$, is used to compute a detection test. Cotrending is performed before detection, and the systematic noise estimate is formed without modeling a particular transit signal. Such an approach favors the null hypothesis in a detection test as ${\hat{{\boldsymbol{c}}}}_{0}$ maximizes the null hypothesis posterior. Here instead we find the MAP estimates for systematic noise under both transit signal hypotheses and compare the likelihoods computed from these estimates in the form of the LRT. The conditional MAP estimates of the systematic noise under either hypothesis take the form

$\begin{eqnarray}\begin{array}{rcl}{\hat{{\boldsymbol{c}}}}_{{\boldsymbol{1}},{\boldsymbol{i}}}^{\mathrm{MAP}} & = & \mathop{\arg \,\max }\limits_{{{\boldsymbol{c}}}_{i}\in {{\boldsymbol{C}}}_{i}}{p}_{1}({{\boldsymbol{c}}}_{i}| {{\boldsymbol{y}}}_{i})\\ & = & \mathop{\arg \,\max }\limits_{{{\boldsymbol{c}}}_{i}\in {{\boldsymbol{C}}}_{i}}{p}_{1}({{\boldsymbol{y}}}_{i}| {{\boldsymbol{c}}}_{i})p({{\boldsymbol{c}}}_{i})\end{array}\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&=\mathop{\arg \,\max }\limits_{{{\boldsymbol{c}}}_{i}\in {{\boldsymbol{C}}}_{i}}\left(\mathrm{ln}[{p}_{1}({{\boldsymbol{y}}}_{i}| {{\boldsymbol{c}}}_{i})]+\mathrm{ln}[p({{\boldsymbol{c}}}_{i})]\right)\end{eqnarray} \tag{ 14 }$

$\begin{eqnarray}&&{\hat{{\boldsymbol{c}}}}_{{\boldsymbol{0}},{\boldsymbol{i}}}^{\mathrm{MAP}}=\mathop{\arg \,\max }\limits_{{{\boldsymbol{c}}}_{i}\in {{\boldsymbol{C}}}_{i}}\left(\mathrm{ln}[{p}_{0}({{\boldsymbol{y}}}_{i}| {{\boldsymbol{c}}}_{i})]+\mathrm{ln}[p({{\boldsymbol{c}}}_{i})]\right).\end{eqnarray} \tag{ 15 }$

We insert these estimates into the hypothesis test described in Equation (2) as

$\begin{eqnarray}&&{H}_{0}:{{\boldsymbol{y}}}_{i}={\boldsymbol{V}}{\hat{{\boldsymbol{c}}}}_{{\boldsymbol{0}},{\boldsymbol{i}}}^{\mathrm{MAP}}+{{\boldsymbol{s}}}_{i},\end{eqnarray} \tag{ 16 }$

$\begin{eqnarray}&&{H}_{1}:{{\boldsymbol{y}}}_{i}={\boldsymbol{t}}+{\boldsymbol{V}}{\hat{{\boldsymbol{c}}}}_{{\boldsymbol{1}},{\boldsymbol{i}}}^{\mathrm{MAP}}+{{\boldsymbol{s}}}_{i}.\end{eqnarray} \tag{ 17 }$

Because the systematic noise signal for either hypothesis is fixed, the only stochastic variable is the Gaussian stellar noise ${{\boldsymbol{s}}}_{i}$. The detection test may be equivalently restated as

$\begin{eqnarray}&&{H}_{0}:\hat{{{\boldsymbol{y}}}_{{\boldsymbol{i}}}}={{\boldsymbol{s}}}_{i}\end{eqnarray} \tag{ 18 }$

$\begin{eqnarray}&&{H}_{1}:\hat{{{\boldsymbol{y}}}_{{\boldsymbol{i}}}}={{\boldsymbol{k}}}_{i}+{{\boldsymbol{s}}}_{i}.\end{eqnarray} \tag{ 19 }$

This hypothesis test on $\hat{{{\boldsymbol{y}}}_{{\boldsymbol{i}}}}={{\boldsymbol{y}}}_{i}-{\boldsymbol{V}}{\hat{{\boldsymbol{c}}}}_{{\boldsymbol{0}},{\boldsymbol{i}}}^{\mathrm{MAP}}$ describes the detection of a known signal ${{\boldsymbol{k}}}_{i}={\boldsymbol{t}}-{\boldsymbol{V}}$ ${\hat{{\boldsymbol{c}}}}_{{\boldsymbol{0}},{\boldsymbol{i}}}^{\mathrm{MAP}}$ + ${\boldsymbol{V}}{\hat{{\boldsymbol{c}}}}_{{\boldsymbol{1}},{\boldsymbol{i}}}^{\mathrm{MAP}}$ in Gaussian noise ${{\boldsymbol{s}}}_{i}$. Therefore, as described above in Section 2.2.1, the optimal detector is the matched filter with test statistic ${T}_{i}(\hat{{{\boldsymbol{y}}}_{{\boldsymbol{i}}}})$:

$\begin{eqnarray}&&{T}_{i}(\hat{{{\boldsymbol{y}}}_{{\boldsymbol{i}}}})=\displaystyle \frac{{\hat{{\boldsymbol{y}}}}_{i}^{T}{\mathrm{Cov}}_{s,i}^{-1}{{\boldsymbol{k}}}_{i}}{\sqrt{{{\boldsymbol{k}}}_{i}^{T}{\mathrm{Cov}}_{s,i}^{-1}{{\boldsymbol{k}}}_{i}}}\underset{{H}_{0}}{\overset{{H}_{1}}{\gtrless }}\tau .\end{eqnarray} \tag{ 20 }$

A systematic noise prior p(${{\boldsymbol{c}}}_{i}$) in Gaussian form ${{\boldsymbol{c}}}_{i}\,\sim { \mathcal N }({{\boldsymbol{\mu }}}_{c,i},{\mathrm{Cov}}_{c,i})$ produces detectors in closed analytic form. A secondary motivation for adopting a Gaussian prior is that it has desirable nonparametric modeling properties when the underlying distribution is not fully known. Specifically, it is inherently regularizing and also allows correlations to be easily modeled. We motivate the choice of a Gaussian prior for ${{\boldsymbol{c}}}_{i}$ by examining the sample density of coefficient values $\left\{{\hat{{\boldsymbol{c}}}}_{i}\right\}$ obtained via least-squares fits of {${{\boldsymbol{v}}}_{k}$} to {${{\boldsymbol{y}}}_{i}$}. Notwithstanding the issues with frequentist estimation of systematic noise noted earlier, it is assumed here that the least-squares fits $\left\{{\hat{{\boldsymbol{c}}}}_{i}\right\}$ are approximately unbiased samples of the coefficient distribution of the population such that they may be used in the inference of parameters describing p(${{\boldsymbol{c}}}_{i}$). Figure 2 shows a subset of systematic basis noise vectors {${{\boldsymbol{v}}}_{k}$} obtained here for CCD module 8 data during Kepler observing quarter 2 using PCA with model order K = 20. This model order was chosen empirically as sufficient to be representative of the systematic noise while low enough as compared to the size of the light-curve population $| I| =5000$ so as to avoid overfitting. In Figure 3, a histogram of the set of sample least-squares coefficient values $\left\{{\hat{c}}_{i}^{1}\right\}$ is shown for the basis vectors in Figure 2; superscript 1 denotes the first element of each coefficient vector. A best-fit Gaussian coefficient prior is overlaid on the histogram. The histogram in Figure 3 has approximately Gaussian form but with a narrower central mode and a broader tail in the distribution; collectively, these broaden the overlaid Gaussian fit. The extremal outliers indicate excursions from our assumptions underlying the systematic noise model. We examine the impact of any non-Gaussianity in the systematic noise prior p(${{\boldsymbol{c}}}_{i}$) on detector performance in Section 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Detector A computes an LRT by marginalizing over the systematic noise distribution, as described in Section 2.2.2. A closed analytic form of this detector is derived here under the assumption of a Gaussian prior for the coefficients of the systematic noise basis vectors. If the coefficient prior p(${{\boldsymbol{c}}}_{i}$) is Gaussian ${{\boldsymbol{c}}}_{i}\sim N({\mu }_{c,i},{\mathrm{Cov}}_{c,i})$, so too are the marginal probabilities of the observed signal p_h(${{\boldsymbol{y}}}_{i}$): h ∈ {0, 1} (see Equation (9) and Equation (10)):

$\begin{eqnarray}&&{p}_{0}({\boldsymbol{y}})\sim { \mathcal N }({\boldsymbol{V}}{\mu }_{c,i},{\mathrm{Cov}}_{{\boldsymbol{s}},i}+{\boldsymbol{V}}{\mathrm{Cov}}_{{\boldsymbol{c}},i}{{\boldsymbol{V}}}^{{\boldsymbol{T}}}),\end{eqnarray} \tag{ 21 }$

$\begin{eqnarray}&&{p}_{1}({\boldsymbol{y}})\sim { \mathcal N }({\boldsymbol{t}}+{\boldsymbol{V}}{\mu }_{c,i},{\mathrm{Cov}}_{{\boldsymbol{s}},i}+{\boldsymbol{V}}{\mathrm{Cov}}_{{\boldsymbol{c}},i}{{\boldsymbol{V}}}^{{\boldsymbol{T}}}),\end{eqnarray} \tag{ 22 }$

where ${\mathrm{Cov}}_{{\boldsymbol{s}},i}$ is the covariance of the zero-mean stellar signal ${{\boldsymbol{s}}}_{i}$. The marginal detector formed from Equation (12) is therefore equivalent under these assumptions to the detection of a known signal ${\boldsymbol{t}}$ in the presence of Gaussian noise ${{\boldsymbol{z}}}_{i}\sim { \mathcal N }({\boldsymbol{0}},{\mathrm{Cov}}_{{\boldsymbol{s}},i}$+${\boldsymbol{V}}{\mathrm{Cov}}_{{\boldsymbol{c}},i}{{\boldsymbol{V}}}^{{\boldsymbol{T}}})$ within the data ${\hat{{\boldsymbol{y}}}}_{i}={{\boldsymbol{y}}}_{i}\,-{\boldsymbol{V}}{\mu }_{c,i}$. The matched filter is therefore optimal (Section 2.2.1) with test statistic ${T}_{i}(\hat{{{\boldsymbol{y}}}_{{\boldsymbol{i}}}})$:

$\begin{eqnarray}&&{T}_{i}(\hat{{{\boldsymbol{y}}}_{{\boldsymbol{i}}}})=\displaystyle \frac{{\hat{{\boldsymbol{y}}}}_{i}^{T}{{\boldsymbol{C}}}_{{\boldsymbol{z}},{\boldsymbol{i}}}^{-1}{\boldsymbol{t}}}{\sqrt{{{\boldsymbol{t}}}^{T}{{\boldsymbol{C}}}_{{\boldsymbol{z}},{\boldsymbol{i}}}^{-1}{\boldsymbol{t}}}}\underset{{H}_{0}}{\overset{{H}_{1}}{\gtrless }}\tau .\end{eqnarray} \tag{ 23 }$

Detector B produces MAP systematic estimates conditioned on the null and alternate hypothesis ${\hat{{\boldsymbol{c}}}}_{{\boldsymbol{h}},{\boldsymbol{i}}}^{\mathrm{MAP}}:{\text{}}h\in \left\{0,1\right\}$ as defined in Equations (13) and (15). These distinct conditional estimates are used as input to the detection test in Equation (20). Closed-form expressions for MAP/MMSE⁷ estimates ${\hat{{\boldsymbol{c}}}}_{{\boldsymbol{h}},{\boldsymbol{i}}}^{\mathrm{MAP}/\mathrm{MMSE}}:{\text{}}h\in \left\{0,1\right\}$ are provided based on a Gaussian systematic noise prior ${{\boldsymbol{c}}}_{i}\sim N({\mu }_{c,i},{\mathrm{Cov}}_{c,i})$ and a zero-mean stellar signal ${{\boldsymbol{s}}}_{i}$ with covariance ${\mathrm{Cov}}_{s,i}$. These estimates can be found either by computing the expectation with respect to the posterior distribution ${{\mathbb{E}}}_{h}({{\boldsymbol{c}}}_{i}| {\boldsymbol{y}})$ or by maximizing the log of the posterior as shown in Smith et al. (2012):

$\begin{eqnarray}\begin{array}{rcl}{\hat{{\boldsymbol{c}}}}_{{\boldsymbol{0}},{\boldsymbol{i}}}^{\mathrm{MAP}/\mathrm{MMSE}} & = & {({{\boldsymbol{V}}}^{{\boldsymbol{T}}}{\mathrm{Cov}}_{s,i}^{-1}{\boldsymbol{V}}+{\mathrm{Cov}}_{c,i}^{-1})}^{-1}\\ & & \times \,({{\boldsymbol{V}}}^{{\boldsymbol{T}}}{\mathrm{Cov}}_{s,i}^{-1}{{\boldsymbol{y}}}_{i}+{\mathrm{Cov}}_{c,i}^{-1}{\mu }_{c,i}),\end{array}\end{eqnarray} \tag{ 24 }$

$\begin{eqnarray}\begin{array}{rcl}{\hat{{\boldsymbol{c}}}}_{{\boldsymbol{1}},{\boldsymbol{i}}}^{\mathrm{MAP}/\mathrm{MMSE}} & = & {({{\boldsymbol{V}}}^{{\boldsymbol{T}}}{\mathrm{Cov}}_{s,i}^{-1}{\boldsymbol{V}}+{\mathrm{Cov}}_{c,i}^{-1})}^{-1}\\ & & \times \,({{\boldsymbol{V}}}^{{\boldsymbol{T}}}{\mathrm{Cov}}_{s,i}^{-1}({{\boldsymbol{y}}}_{i}-{\boldsymbol{t}})+{\mathrm{Cov}}_{c,i}^{-1}{\mu }_{c,i}).\end{array}\end{eqnarray} \tag{ 25 }$

In their implementation, the joint detectors A and B described above require several parameters characterizing the systematic and stellar noise to be estimated or held fixed. This is not trivial as there are no clean, well-separated data. The parameters that need to be estimated include those of the Gaussian prior ${{\boldsymbol{c}}}_{i}\sim N({\mu }_{c,i},{\mathrm{Cov}}_{c,i})$ for the coefficients ${{\boldsymbol{c}}}_{i}$ of the systematic noise basis vectors ${{\boldsymbol{v}}}_{k}$: k ∈ K, for fixed model order K. In addition, the parameters of the statistical distribution of the stellar signal ${{\boldsymbol{s}}}_{i}$ need to be estimated.

In the current work, we estimate the systematic noise basis vectors {${{\boldsymbol{v}}}_{k}$} using PCA as implemented in the Python module sklearn.⁸ To suppress the inclusion of stellar noise or dominating outlier light curves, the basis is constructed from 90% of the total light curves—those which have the lowest variance in absolute value. The parameter values of the Gaussian coefficient prior ${{\boldsymbol{c}}}_{i}\sim N({\mu }_{c,i},{\mathrm{Cov}}_{c,i})$ are estimated directly from the set of coefficient estimates $\left\{{\hat{{\boldsymbol{c}}}}_{i}\right\}$ obtained from least-squares fits of the systematic basis vectors {${{\boldsymbol{v}}}_{k}$} to the raw light curves {${{\boldsymbol{y}}}_{i}$}. Here we assume that for a population I the same coefficient covariance may be used, denoted by ${\mathrm{Cov}}_{c,I}$. An example sample covariance is shown in Figure 4. It can be seen that in this example the coefficient values are correlated (by the presence of nonzero off-diagonal elements). The PCA method finds an orthogonal set of basis vectors {${{\boldsymbol{v}}}_{k}$}, but there is no reason to expect independent systematic noise signals themselves to be orthogonal. Hence, the orthogonalization procedure may distribute a systematic noise signal across multiple basis vectors. This will likely produce correlations between coefficients.

Zoom In Zoom Out Reset image size

While the coefficient covariance ${\mathrm{Cov}}_{c,I}$ is estimated from the global population of least-square fits, the coefficient mean μ_{c, i} is taken to be the least-square coefficient vector ${\hat{{\boldsymbol{c}}}}_{i}$ obtained for the particular light curve under consideration.

As noted above, we assume that the stellar noise ${{\boldsymbol{s}}}_{i}$ is wide-sense-stationary (WSS) and zero-mean. This implies the stellar covariance matrix is Toeplitz and fully parameterizes the stellar noise. We therefore only need to estimate the stellar spectrum (N free parameters) as opposed to a full covariance matrix (N² free parameters). The spectral estimation technique is not prescribed by the form of the detectors; here we use a smoothed periodogram (Kay 1999b) on the least-squares cotrended light curves.

We evaluate our detection performance using single-transit injection tests (Gilliland et al. 2000; Weldrake et al. 2005; Burke et al. 2006; Burke & Catanzarite 2017) on raw Kepler SAP light curves (Jenkins et al. 2010b) selected to exclude known exoplanet detections. We use the long-cadence Kepler data in this analysis for which there is a 29.4 minute integration time (Jenkins et al. 2010a). A subset of 20,000 of such Kepler light curves were selected, comprising 5000 light curves from each of the following Kepler CCD module and observing quarter pairs [M6:Q10, M8:Q2, M14:Q9, M18:Q4].⁹ These module and quarter pairs were selected randomly over time and CCD module position. Within each pair, the light curves were sorted by angular separation from the module reference point used by the MAST Kepler data archive, and the first 5000 were selected from the sorted list. We exclude any light curves associated with exoplanets defined as confirmed by the NASA Exoplanet Archives.¹⁰ A single synthetic transit is injected once per light curve; transit signals were simulated using the Python transit¹¹ library developed by D. Foreman-Mackey. These synthetic transit signals include limb darkening (Mandel & Agol 2002; Kipping 2013a) and a complete description of the transiting Keplerian orbital elements. The injected signals are drawn from a distribution of exoplanet population parameters given in Table 1. This population parameter distribution is informed by that used by Foreman-Mackey et al. (2015) and Kipping (2013a, 2013b), but it is not identical. We adopted zero orbital eccentricity in the current work, among other changes.

As a reference detector, we adopt the standard heuristic of sequential cotrending followed by detection (Stumpe et al. 2012). We term this the standard model in what follows and provide our own implementation of this detector in the current work. In the standard model, the cotrending is performed assuming the light curve contains no transit signal: ${{\boldsymbol{y}}}_{i}$ = ${{\boldsymbol{s}}}_{i}+{{\boldsymbol{Vc}}}_{i}$, analogous to hypothesis H₀ (Equation (16)). An MAP/MMSE estimator for the systematic noise is constructed equivalent to ${\hat{{\boldsymbol{c}}}}_{0}^{\mathrm{MAP}/\mathrm{MMSE}}$ (Equation (15)) and applying the same Bayesian priors for systematic noise and stellar noise as used by detectors A and B. Detection is then performed on the cotrended light curves ${\hat{{\boldsymbol{y}}}}_{i}={{\boldsymbol{y}}}_{i}-{\boldsymbol{V}}{\hat{{\boldsymbol{c}}}}_{{\boldsymbol{i}},{\boldsymbol{0}}}^{\mathrm{MAP}/\mathrm{MMSE}}$ using the matched filter in Equation (8) with ${\mathrm{Cov}}_{n}={\mathrm{Cov}}_{s,i}$.

As discussed in Section 2.1.3, transit detection requires testing every candidate transit signal ${\boldsymbol{t}}$ ∈ ${\boldsymbol{T}}$ to find that which maximizes the test statistic T(${\boldsymbol{y}}$). The transit space ${\boldsymbol{T}}$ is typically populated by periodic box functions over a range of candidate orbital periods P, transit durations d, and epoch times t₀ in the functional form described by Equation (5). Transit depth α is omitted here as our detector forms generally do not make use of this parameter. The dimensionality of the transit parameter search space is therefore intrinsically large and the transit detection problem computationally expensive. This computational cost can be reduced sharply by constraining the range of epochs {t₀} for a candidate transit with a certain period and duration using the method of phase correlation (Averbuch & Keller 2002). We adopt this method in our numerical studies due to the significant reduction in computational cost. The phase-correlation method and its application to epoch estimation is described in Appendix A.

The phase-correlation method, however, requires the use of cotrended light curves for sufficient accuracy in the estimated epochs {t₀}. This raises the concern that the transit signal may not be detected optimally due to the use of the cotrended, as opposed to raw, light curves. To verify that this approach does not decrease detection efficiency, we compared detection results with and without phase correlation using single-transit injection tests over 5000 light curves from the broader injection test data described in Section 2.4.1 selected here from [M8:Q2]. We define detection efficiency in this context as the rate of correct detection of the known injected signals as described in Section 3.1. In each case, the standard reference detector defined in Section 2.4.2 was used; as described above, this detector comprises sequential cotrending and detection steps. In the test without phase correlation, a three-dimensional transit signal parameter search space was used as defined in Table 2; the standard detector operated on the raw light curves over this gridded search space. In the phase-correlation test, a transit epoch t₀ was estimated from each least-squares cotrended light curve using the phase-correlation method. The standard detector was then applied to the raw light curves holding {t₀} fixed to the phase-correlation estimate but searching over a residual two-dimensional search space in period and duration as tabulated in Table 2.

Table 2.  Correlation Verification Transit Search Space

Transit Parameter	Range	Step Size	Physical Units
	($\bigtriangleup {t}_{\mathrm{LC}}$)	($\bigtriangleup {t}_{\mathrm{LC}}$)
Period	PLC ∈ [20, 2125]	1	P ∈ [1, 43.4] days
Duration	dLC ∈ [3, 11]	2	d ∈ [1.4, 5.4] hr
Epoch	[0, PLC]	${d}_{{\rm{LC}}}/2$

Note. The long-cadence sample integration time is $\bigtriangleup {t}_{\mathrm{LC}}=29.4\,\mathrm{minutes}$.

Download table as:  ASCIITypeset image

The test data here comprise actual raw light curves from which confirmed exoplanets have been excluded; however, the data cannot be shown provably to exclude any hitherto undetected transit signals. As such, we define a quasi-false-alarm rate as the rate of incorrect detection with respect to the injected signal set. The detection tests using the phase-correlation method show an improvement in detection efficiency of 14% and a reduction in the quasi-false-alarm rate of 16% (at a detection threshold τ = 8.4) over the detection tests for which a direct search was performed. A comparison of detection efficiency broken down by transit parameter values is shown in Figure 5 and demonstrates no marked reduction in detection for weak signals. The phase-correlation method described in Appendix A by definition has a maximum accuracy Δt₀ in epoch t₀ of one long-cadence sample ${\rm{\Delta }}{t}_{{\rm{LC}}}$ (a single "pixel"). Nonadditive noise will reduce the accuracy of the method. By considering the minimum required correlation between a measured transit and a parameterized transit model (t₀, d, P), Jenkins et al. (2010c) provide an analysis motivating a default search spacing in the epoch of ${\rm{\Delta }}{t}_{0}=d/10$ for the Transiting Planet Search (TPS) module in the Kepler science pipeline. Our direct search here used a step size ${\rm{\Delta }}{t}_{0}=d/2$ (see Table 2), which is suboptimal relative to the maximum epoch accuracy of the phase-correlation method, thereby possibly explaining the improved detection efficiency of the latter method here. We stress here that these tests only demonstrate that the phase-correlation optimization is suitable for the transit search space considered here; generalization to broader applicability is left to future work.

Zoom In Zoom Out Reset image size

The ratio of computational cost between the test without phase correlation and that using phase correlation was ∼10²: 1. As discussed further in Section 4, the freed computational resources allow more refined searches in the remaining transit parameters and can be argued to improve overall accuracy in that sense. As a result of the positive outcome of this verification test and the significant associated reduction in computational cost, we used the phase-correlation method described in Appendix A to estimate transit epochs {t₀} in our full injection tests described in Section 2.4.1. For the full injection tests, the transit signal parameter search space is informed broadly by Jenkins et al. (2002). However, given our use of the phase-correlation method and the associated freed computational resources, we search (in units of long-cadence samples $\bigtriangleup {t}_{\mathrm{LC}}=29.4\ \mathrm{minutes}$) over the period range [20, 2125] with a 0.25 step size and within the following set of transit durations (in long-cadence samples): [2, 3, 4, 5, 6, 7, 9, 10, 12]. Detector B additionally requires a search over transit-depth parameter α. This arises from the estimation of transit-dependent systematics as in Equation (25) for ${\hat{{\boldsymbol{c}}}}_{{\boldsymbol{1}},{\boldsymbol{i}}}^{\mathrm{MAP}/\mathrm{MMSE}}$. In contrast, the matched filter depends purely on the shape of the transit signal (defined by t₀, d, P) and not the signal strength α. One can see this property by considering a scaled signal α${\boldsymbol{t}}$ in the matched filter function, Equation (8); the scaling α immediately cancels. For the rest of the detectors, the only step that is dependent on a transit signal is a matched filter step, ergo they do not depend on a transit-depth parameter. For detector B, we search over four equally spaced transit depths α ∈ {0.2, 0.5, 0.8, 1.1}, scaled by the maximum range of the least-squares cotrended light curve under consideration. Our choice of transit-depth sampling is exploratory but proved practical. We note, however, that it sets a limit on the detectability of signals with transit depths below 20% of the cotrended light curve. Future work will consider optimized sampling schemes for transit depth including estimated noise levels.

The computational complexity of the detectors and their key constituent operations acting on a single light curve are summarized in Table 3. The computational complexities are expressed in terms of the length of the light curve N and the size of the transit signal search space $| {\boldsymbol{T}}|$ as defined in the introduction of Section 2. These complexities are specific to the case of a Gaussian systematics prior described in Section 2.3.1.

Table 3.  Computational Complexity of Operations on a Single Light Curve of Length N

Operation	Complexity
Time-domain matched filter	O(N2)
Fourier-domain matched filter	$O(N\mathrm{log}N)$
MAP/MMSE systematics estimation	O(N2)
Standard detector (Gaussian prior)	$O(| {\boldsymbol{T}}| N\mathrm{log}N)$
Detector A (Gaussian prior)	$O(| {\boldsymbol{T}}| {N}^{2})$
Detector B (Gaussian prior)	$O(| {{\boldsymbol{T}}}^{* }| {N}^{2})$

Note. Where $| {\boldsymbol{T}}|$ is the search space size. For detector B, the search space size $| {{\boldsymbol{T}}}^{* }|$ is generally larger as one must additionally search over candidate transit depths $| \alpha |$.

Download table as:  ASCIITypeset image

A time-domain matched filter defined in the form of Equation (8) has a computational complexity dominated by the product of an [N × N] matrix and a length N vector; therefore, its complexity is O(N²). A matched filter implemented in the Fourier domain involves computing fast Fourier transforms (FFT) and inner products of length N vectors (Kay 1999b). Between these operations, the FFT is more computationally intensive, hence the Fourier-domain matched filter is $O(N\mathrm{log}N)$ (Bracewell & Bracewell 2000).

The computational complexity of MAP/MMSE systematics estimation is defined by the form of Equation (25). For a particular light curve, we can reuse many of the computed terms during multiple transit tests. Considering transit-dependent terms, the dominant computational term is the product of the [N × N] matrix ${\mathrm{Cov}}_{s,i}^{-1}$ and the length N vector ${\boldsymbol{t}}$; consequentially this computation is O(N²). Once per light curve, a matrix inversion of the N × N matrix ${\mathrm{Cov}}_{s,i}$ is performed, and this operation is O(N³). However, it is not leading order because generally $| {\boldsymbol{T}}| {N}^{2}\gg {N}^{3}$. As such, MAP/MMSE systematics estimation has computational complexity O(N²) per light curve per transit.

The detector complexities are determined by the form of the matched filter used and scaled by the size of the transit search space. The standard detector (Section 2.4.2) is the most computationally efficient detection strategy as per transit the only computation performed is a Fourier-domain matched filter. The standard detector searches a transit space of size $| {\boldsymbol{T}}|$ and therefore has a computational complexity $O(| {\boldsymbol{T}}| N\mathrm{log}N)$.

Detector A (Section 2.3.2) searches for a transit signal contained in non-WSS Gaussian noise; therefore, a time-domain matched filter must be used for each transit test. The search space is of size $| {\boldsymbol{T}}|$ and the net computational complexity is $O(| {\boldsymbol{T}}| {N}^{2})$.

Detector B (Section 2.3.3) uses a larger search space than the other detection strategies as it includes transit depth α; this search space was described in Section 2.4.3 and is of size $| {T}^{* }|$. In addition to this increased search space, this detector must compute an MAP/MMSE systematics estimate once per transit, which is then used as input into a Fourier-domain matched filter. The net computational complexity is therefore $O(| {{\boldsymbol{T}}}^{* }| {N}^{2})$.

Approximate elapsed wall-clock run times are summarized in Table 4. Transit detection tests were parallelized with one light curve assigned to each CPU core. All runs were performed on the Blue Waters petascale system at UIUC/NCSA (Bode et al. 2013). This is a Cray XE/XK system with a peak performance of 13.34 PF.¹²

We conducted an initial feasibility test using detectors A and B and the standard detector on a subset of Kepler light curves that did not exclude known exoplanets. These tests were designed to demonstrate initial detection feasibility only on real exoplanet transit signatures. Detection tests were performed over the transit search space identical to that used in injection tests. The transit signal parameter search space is described in Section 2.4.3. A total of 2000 light curves were analyzed in this test, 1000 light curves were selected from CCD module 2 over observing quarters [Q2, Q10, Q14], and an additional 1000 light curves were selected from CCD module 12 over observing quarters [Q3, Q7, Q15]. These CCD modules and observing quarters were chosen randomly over time and across the CCD module. As was performed for the injection tests, the light curves for each module were first sorted by angular separation from the module reference point used by the MAST Kepler data archive before the first 1000 were selected. No explicit selection for CCD module output was applied: module 2 data included outputs 3 and 4, while module 12 data included only output 4. For our detection tests, we did not use stitched quarters but instead performed separate detection tests on each of the quarters and computed an averaged test statistic (per transit over time).