An update to the EVEREST K2 pipeline: Short cadence, saturated stars, and Kepler-like photometry down to Kp = 15

TL;DR

The EVEREST 2.0 pipeline significantly enhances K2 photometric data quality by improving de-trending methods, including saturated and faint stars, achieving near-Kepler precision across a wide magnitude range, and is adaptable to other surveys.

Contribution

This paper introduces a new regularization scheme for pixel level decorrelation that includes neighboring stars and handles saturated stars, improving photometric precision over previous versions.

Findings

10-20% higher photometric precision than previous EVEREST versions

Recover Kepler-like precision up to Kp=14 and Kp=15 for certain campaigns

De-trended all short cadence targets with improved accuracy

Abstract

We present an update to the EVEREST K2 pipeline that addresses various limitations in the previous version and improves the photometric precision of the de-trended light curves. We develop a fast regularization scheme for third order pixel level decorrelation (PLD) and adapt the algorithm to include the PLD vectors of neighboring stars to enhance the predictive power of the model and minimize overfitting, particularly for faint stars. We also modify PLD to work for saturated stars and improve its performance on extremely variable stars. On average, EVEREST 2.0 light curves have 10-20% higher photometric precision than those in the previous version, yielding the highest precision light curves at all Kp magnitudes of any publicly available K2 catalog. For most K2 campaigns, we recover the original Kepler precision to at least Kp = 14, and to at least Kp = 15 for campaigns 1, 5, and 6. We also de-trend all short cadence targets observed by K2, obtaining even higher photometric precision for these stars. All light curves for campaigns 0-8 are available online in the EVEREST catalog, which will be continuously updated with future campaigns. EVEREST 2.0 is open source and is coded in a general framework that can be applied to other photometric surveys, including Kepler and the upcoming TESS mission.

Figures

Equations

$$\mathbf{m}=$$

$$\sum _i\sum _j{b}_{ij}\frac{{\mathbf{p}}_i{\mathbf{p}}_j}{(\sum _n{\mathbf{p}}_n{)}^2}+$$

$$\sum _i\sum _j\sum _k{c}_{ijk}\frac{{\mathbf{p}}_i{\mathbf{p}}_j{\mathbf{p}}_k}{(\sum _n{\mathbf{p}}_n{)}^3}.$$

$$\mathbf{m}=\mathbf{X}\cdot \mathbf{w},$$

$${\mathbf{y}}^{\prime }=\mathbf{y}-\mathbf{m}.$$

$$\log {\mathcal{L}}_0=$$

$$-\frac{1}{2}\log \left|\mathbf{K}\right|-\frac{N_{dat}}{2}\log 2\pi ,$$

$$\log \mathcal{L}=\log {\mathcal{L}}_0-\frac{1}{2}{\mathbf{w}}^{\top }\cdot {\mathbf{\Lambda }}^{-1}\cdot \mathbf{w}-\frac{1}{2}\log \left|\mathbf{\Lambda }\right|,$$

$${\Lambda }_{m,n}={\lambda }_n^2{\delta }_{mn}.$$

$$\frac{\mathrm{d}\log \mathcal{L}}{\mathrm{d}\mathbf{w}}$$

$$-\left({\mathbf{X}}^{\top }\cdot {\mathbf{K}}^{-1}\cdot \mathbf{X}+{\mathbf{\Lambda }}^{-1}\right)\cdot \mathbf{w}.$$

$$\hat{\mathbf{w}}={\left({\mathbf{X}}^{\top }\cdot {\mathbf{K}}^{-1}\cdot \mathbf{X}+{\mathbf{\Lambda }}^{-1}\right)}^{-1}\cdot {\mathbf{X}}^{\top }\cdot {\mathbf{K}}^{-1}\cdot \mathbf{y}$$

$$\mathbf{m}=\mathbf{X}\cdot {\left({\mathbf{X}}^{\top }\cdot {\mathbf{K}}^{-1}\cdot \mathbf{X}+{\mathbf{\Lambda }}^{-1}\right)}^{-1}\cdot {\mathbf{X}}^{\top }\cdot {\mathbf{K}}^{-1}\cdot \mathbf{y}.$$

$$\displaystyle\mathbf{X}=\left(\begin{array}[]{ccc}\mathbf{X_{1}}&\mathbf{X_{2}}&\mathbf{X_{3}}\end{array}\right),$$

$$\displaystyle\mathbf{\Lambda}=\left(\begin{array}[]{ccc}\mathbf{\Lambda_{1}}&&\\ &\mathbf{\Lambda_{2}}&\\ &&\mathbf{\Lambda_{3}}\\ \end{array}\right)$$

$$\mathbf{m}=\mathbf{X}\cdot \mathbf{\Lambda }\cdot {\mathbf{X}}^{\top }\cdot {\left(\mathbf{X}\cdot \mathbf{\Lambda }\cdot {\mathbf{X}}^{\top }+\mathbf{K}\right)}^{-1}\cdot \mathbf{y}.$$

$$\mathbf{X}\cdot \mathbf{\Lambda }\cdot {\mathbf{X}}^{\top }$$

$$={\lambda }_1^2{\mathbf{X}}_1\cdot {\mathbf{X}}_1^{\top }+{\lambda }_2^2{\mathbf{X}}_2\cdot {\mathbf{X}}_2^{\top }+{\lambda }_3^2{\mathbf{X}}_3\cdot {\mathbf{X}}_3^{\top }$$

$$=\sum _n{\lambda }_n^2{\mathbf{X}}_{\mathbf{n}}^2,$$

$${\mathbf{X}}_{\mathbf{n}}^2\equiv {\mathbf{X}}_{\mathbf{n}}\cdot {\mathbf{X}}_{\mathbf{n}}^{\top }.$$

$$\mathbf{m}=\sum _n{\lambda }_n^2{\mathbf{X}}_{\mathbf{n}}^2\cdot {\left(\sum _n{\lambda }_n^2{\mathbf{X}}_{\mathbf{n}}^2+\mathbf{K}\right)}^{-1}\cdot \mathbf{y}.$$

$$\displaystyle\mathbf{X}=\left(\begin{array}[]{cccccc}\mathbf{X_{1}}&\mathbf{X_{1}^{\prime}}&\mathbf{X_{2}}&\mathbf{X_{2}^{\prime}}&\mathbf{X_{3}}&\mathbf{X_{3}^{\prime}}\end{array}\right),$$

$${\mathbf{K}}_{ij}=\alpha \left(1+\frac{\sqrt{3(t_i-t_j{)}^2}}{\tau }\right)e^{-\frac{\sqrt{3(t_i-t_j{)}^2}}{\tau }}+{\sigma }^2{\delta }_{ij},$$

$${\hat{\lambda }}_1$$

$${\hat{\lambda }}_2$$

$${\hat{\lambda }}_3$$

Full text

1. Introduction

The failure of the second of four reaction wheels in 2013 brought the original Kepler mission to a premature conclusion, as the spacecraft could no longer achieve the fine pointing precision required for the groundbreaking transiting exoplanet and stellar variability science that the first four years of the mission allowed. Since 2014 the spacecraft has been operating in a new mode, known as $K2$, using the solar wind and periodic thruster firings to mitigate drift (Howell et al., 2014). Despite these measures, raw $K2$ photometry is significantly poorer than that of the original Kepler mission. In order to enable continuing precision science with $K2$, numerous pipelines have been developed (e.g., Vanderburg & Johnson, 2014; Armstrong et al., 2015; Lund et al., 2015; Crossfield et al., 2015; Foreman-Mackey et al., 2015; Huang et al., 2015; Aigrain et al., 2016), many producing light curves with precision approaching that of Kepler for bright stars.

In Luger et al. (2016) (henceforth Paper I), we developed the first version of our $K2$ pipeline, EVEREST (EPIC Variability Extraction and Removal for Exoplanet Science Targets). We employed a variant of pixel level decorrelation (PLD) based on the method of Deming et al. (2015), a data-driven approach that uses a star’s own pixel-level light curve to remove instrumental effects. We showed that EVEREST recovered the original Kepler precision for stars brighter than $Kepler$-band magnitude $Kp\approx 13$, yielding higher average precision than any publicly available $K2$ catalog for unsaturated stars.

In this paper, we present an update to our pipeline, which we refer to as EVEREST 2.0. By combining PLD with spacecraft motion information obtained from nearby stars, this update improves the precision of $K2$ light curves at all magnitudes relative to version 1.0 and addresses certain issues with overfitting. EVEREST 2.0 also reliably de-trends saturated stars and stars observed in short cadence mode, obtaining comparable or even higher de-trending power for these targets.

The paper is organized as follows: in §2 we derive the mathematical framework of the EVEREST 2.0 model, and in §3 we describe the implementation of our pipeline in detail. We present our results in §4 and some additional remarks in §5. In §6 we briefly discuss how to use the EVEREST catalog and code, and in §7 we summarize the work.

2. The PLD Model

Here we describe the mathematical formulation of the EVEREST pixel level decorrelation (PLD) model. In PLD, products of the fractional fluxes in each pixel of the target aperture are used as regressors in a linear model:

[TABLE]

In the expression above, $\mathbf{m}$ is the instrumental model and $\mathbf{p}_{i}$ is the flux in the $i^{th}$ pixel; both are vector quantities defined at an array of times $\mathbf{t}$. Each term corresponds to a different PLD order (first, second, and third) resulting from a Taylor expansion of the instrumental signal. The $a_{i}$, $b_{ij}$, and $c_{ijk}$ are the linear weights of the model, which we seek to obtain below. For a detailed discussion of the theory behind PLD, see Deming et al. (2015) and Paper I. Below we simply discuss its mathematical implementation.

2.1. Regularized Regression (rPLD)

Given a timeseries $\mathbf{y}$ with $N_{dat}$ data points, we wish to find the linear combination of $N_{reg}$ regressors that best fits the instrumental component of $\mathbf{y}$. Expressed in vector form, our linear model is thus

[TABLE]

where $\mathbf{X}$ is the ($N_{dat}\times N_{reg}$) design matrix constructed from the set of regressors (the fractional pixel fluxes in Equation 2) and $\mathbf{w}$ is the ($N_{reg}\times 1$) vector of weights (the set $\{a_{i},b_{ij},c_{ijk}\}$). If $\mathbf{w}$ is known, the de-trended light curve is simply

[TABLE]

In Paper I, we obtained $\mathbf{w}$ by maximizing the likelihood function

[TABLE]

where $\mathbf{K}$ is the ($N_{dat}\times N_{dat}$) covariance matrix of the data and $\mathbf{y}$ is the ($N_{dat}\times 1$) simple aperture photometry (SAP) flux. Since the number of third order PLD regressors can be quite large (on the order of several thousand for a typical star, which is larger than the number of data points), the problem is ill-posed, meaning that a unique solution does not exist and maximizing $\log\mathcal{L}_{0}$ is likely to lead to overfitting. We thus constructed $\mathbf{X}$ from the (smaller) set of $N_{pc}$ principal components of the PLD regressors. We chose $N_{pc}$ by performing cross-validation, which aims to maximize the predictive power of the model while minimizing overfitting.

However, while principal component analysis (PCA) yields a set of components that captures the most variance among the PLD vectors, there is no guarantee that the principal components are the ideal regressors in the PLD problem. Dimensionality reduction techniques such as PCA inevitably lead to information loss, and so it is worthwhile to consider alternative regression methods to fully exploit the potential of PLD.

A common regression method for ill-posed problems is regularization, in which a prior is imposed on the values of the weights $\mathbf{w}$. Since overfitting occurs when $\mathbf{w}$ becomes very large, regularization recasts the problem by adding a penalty term to the likelihood that increases with increasing $|\mathbf{w}|$. While many forms of regularization exist, we focus on L2 regularization, since it has an analytic solution. Recently, Wang et al. (2016) used L2 regularization in a “causal pixel model” to de-trend light curves from the original Kepler mission. L2 regularization is mathematically equivalent to placing a Gaussian prior on each of the weights $\mathbf{w}$, so that the posterior likelihood function becomes

[TABLE]

where $\mathbf{\Lambda}$ is the ($N_{reg}\times N_{reg}$) regularization matrix, which we choose to be diagonal for simplicity and for computational efficiency:

[TABLE]

Each element $\lambda_{n}^{2}$ in $\mathbf{\Lambda}$ is the variance of the zero-mean Gaussian prior on the weight of the corresponding column of the design matrix, $\mathbf{X}_{*,n}$. In practice we find that if we choose the $\lambda_{n}$ correctly, this model has a higher predictive power than the PCA model adopted in Paper I.

Given this formulation, our task is to find the weights $\mathbf{\hat{w}}$ that maximize the posterior probability $\mathcal{L}$. Differentiating Equation (5) with respect to $\mathbf{w}$, we get

[TABLE]

By setting this expression equal to zero, we obtain the maximum a posteriori prediction for the weights,

[TABLE]

with corresponding model

[TABLE]

In what follows, we refer to the implementation of PLD with regularized regression as rPLD.

2.2. Cross-validation

Similarly to Paper I, we solve for $\mathbf{\Lambda}$ by cross-validation. For each value of $\mathbf{\Lambda}$, the model is trained on one part of the light curve (the training set) and used to de-trend the other part of the light curve (the validation set); see §3.7 for details. The value of $\mathbf{\Lambda}$ that results in the minimum scatter in the validation set is then chosen for the final de-trending step.

In principle, each of the $\lambda_{n}$ in $\mathbf{\Lambda}$ could take on a different value, but solving for each one requires minimizing an $N_{reg}$-dimensional function and is not computationally tractable. Instead, we simplify the problem by requiring that all regressors of the same order have the same regularization parameter $\lambda$. Provided we write the third order design matrix in the form

[TABLE]

where $\mathbf{X_{n}}$ is the matrix of $n^{\mathrm{th}}$ order regressors, we may express the regularization matrix as

[TABLE]

where $\mathbf{\Lambda_{n}}=\lambda_{n}^{2}\mathbf{I}$ is the $n^{\mathrm{th}}$ order regularization matrix and $\lambda_{n}^{2}$ is the variance of the prior on the $n^{\mathrm{th}}$ order regressors.

A typical $K2$ star with 30 aperture pixels has N_{reg}\sim\5,000 regressors and N_{dat}\sim\500 data points in each cross-validation light curve segment (see §3.7). Evaluating the matrix inverse in Equation (9) is thus computationally expensive, and becomes prohibitive during cross-validation, since this must be done for every set of $\lambda_{n}$’s. Fortunately, we can reduce the number of calculations with some linear algebra. First, we apply the Woodbury matrix identity (e.g., Golub & Van Loan, 1996) to Equation (9), obtaining

[TABLE]

Next, we note that

[TABLE]

where we have defined

[TABLE]

We may thus re-write our maximum a posteriori model as

[TABLE]

The matrix that we must invert in Equation (27) has dimensions ($N_{dat}\times N_{dat}$), while that in Equation (9) is ($N_{reg}\times N_{reg}$). Since $N_{reg}\sim 10N_{dat}$, casting the model in this form can greatly speed up the computation. In practice, we pre-compute the three matrices $\mathbf{X^{2}_{n}}$ at the beginning of the cross-validation step, so the only time-consuming operation in Equation (27) is the inversion.

2.3. Neighboring Stars (nPLD)

One of the downsides of PLD is that the regressors used in the linear model tend to be noisy. This is particularly a problem for faint targets, whose PLD vectors are dominated by photon noise. Their light is also distributed over fewer pixels, resulting in a smaller set of vectors on which to regress. The effect of this is evident in Figure 10 of Paper I, which shows how EVEREST 1.0 light curves for faint ($Kp\gtrsim 15$) stars are significantly noisier than those of the original Kepler mission. This decrease in de-trending power at the faint end affects most other $K2$ pipelines as well, as those usually regress on information derived (either directly or indirectly) from the motion of the stellar image across the detector, which is similarly noisy.

While the spacecraft motion (the dominant source of instrumental noise in $K2$) is imprinted at relatively low signal-to-noise ratio (SNR) on the light curves of any one star, the collective light curves of all the stars on the detector encode this information at very high SNR. Therefore, a straightforward way to improve the performance of PLD for faint targets is to include this information in the design matrix. To this end, in EVEREST 2.0 we incorporate the PLD vectors of a set of other targets located on the same CCD module as the target of interest when performing the regression. We dub this method nPLD, for neighboring PLD, and discuss its implementation in §3.4. Our third-order design matrix (Equation 11) is now

[TABLE]

where $\mathbf{X_{n}^{\prime}}$ is the design matrix constructed from the $n^{th}$ order PLD vectors of all the neighboring targets. For computational speed, we still solve for a single prior amplitude $\lambda_{n}$ for each PLD order, but in principle one could assign different priors to the neighboring vectors. We discuss the implementation of nPLD in §3.

3. Implementation

3.1. Light Curves

As in Paper I, we downloaded all stars in the $K2$ EPIC catalog with long and/or short cadence target pixel files and adopted aperture #15 from the K2SFF catalog (Vanderburg, 2014; Vanderburg & Johnson, 2014). We masked all cadences with QUALITY flags 1-9, 11-14, and 16-17, though we still compute the model prediction on them. For campaigns 0-2, we remove the background signal as described in Paper I; for more recent campaigns, the background is removed by the Kepler team.

Next, we perform iterative sigma clipping to identify and mask outliers at 5$\sigma$. During each iteration, we compute the linear (unregularized) PLD model and smooth it with a Savitsky-Golay filter (Savitzky & Golay, 1964), then identify outliers based on a median absolute deviation (MAD) cut. We implement this outlier-clipping step at the beginning of each cross-validation step (§3.7), each time computing the model with a higher (regularized) PLD order, to progressively refine the outlier mask.

3.2. GP Optimization

In order to compute the covariance matrix $\mathbf{K}$ for each target, we use a Gaussian process (GP), as we did in Paper I. GP optimization can be costly, especially when performing model selection over a range of possible kernels and optimizing many hyperparameters simultaneously. For this reason, in Paper I we cut corners and performed kernel selection based on fits to the autocorrelation function of the light curve, which we also used to fix the timescale and/or period of those kernel(s). We then ran a nonlinear minimizer to optimize the overall amplitude of the GP. In practice, this worked reasonably well, but often failed for light curves dominated by high frequency stellar variability. After much experimentation, we decided to forego the kernel selection step in favor of using a single carefully optimized Matérn-3/2 kernel with an added white noise term:

[TABLE]

where the hyperparameters $\sigma$, $\alpha$, and $\tau$ are the white noise amplitude, red noise amplitude, and red noise timescale, respectively, and $t_{i}$ and $t_{j}$ correspond to the timestamps of cadences $i$ and $j$. We initialize the hyperparameters at random values and run a nonlinear optimizer to solve for the maximum likelihood (Equation 5), keeping the PLD model parameters fixed; we repeat this process several times and retain the highest likelihood solution. As with outlier clipping, we progressively optimize the GP at each of the three cross-validation steps, so that each time we train the GP on a light curve that is increasingly dominated by stellar variability (as opposed to instrumental systematics).

In principle, the quasi-periodic kernels used in Paper I should be better suited to handling variable stars, but in practice we find that a properly optimized Matérn-3/2 kernel is flexible enough to fully capture the variability and prevent PLD overfitting. We discuss this further in §5.

3.3. Breakpoints

Because the instrumental noise properties are quite variable over the course of $K2$ campaigns, we find a significant improvement in the de-trending power of our regularized regression model when we subdivide light curves into two or three segments. This is in contrast to the PCA approach in Paper I, where we did not find it necessary to split the timeseries. For all campaigns except 4 and 7, we add a single breakpoint in the light curve near the mid-campaign point, where the spacecraft roll is at a minimum. For campaigns 4 and 7, we find it necessary to insert two breakpoints. We cross-validate and de-trend each light curve segment separately and mend them at the end. In order to mitigate flux discontinuities at the breakpoints, we train the model in each segment on an additional 100 cadences past the breakpoint to remove potential edge effects and offset the models in each segment so that they align at the breakpoint. While this method removes flux discontinuities, it can introduce discontinuities in the derivative of the flux, showing up as spurious “kinks” in the light curve. We remove these in a post-processing step (§4.7).

3.4. Neighboring Stars

In principle, the larger the number of neighboring PLD vectors we include in the nPLD design matrix, the higher the de-trending power of our model. However, adding regressors significantly increases computing time, so we would like to instead select a small set of high SNR regressors that capture most of the spacecraft motion information. Moreover, since we employ a single prior for all $n^{th}$ order regressors, adding many foreign PLD vectors effectively dilutes the contribution of the target’s own PLD vectors, which are the only ones that can correct instrumental signals arising from local pixel sensitivity variations, and in practice results in poorer quality light curves. After much experimenting, we obtain the highest average de-trending power when the number of neighboring stars is about ten. We therefore de-trend each $K2$ target with the aid of the PLD vectors of ten randomly selected bright ($11\leq Kp\leq 13$) stars on the same detector module as the target. To minimize contamination of the target by outliers in its neighbors’ fluxes, we linearly interpolate over all neighbor data with flagged QUALITY bits. Finally, for computational reasons, we neglect all cross terms of the form $\prod_{i\neq j}p_{i}p_{j}$, where $p_{i}$ is the flux in the $i^{th}$ pixel, when computing the neighbors’ PLD vectors. Cross terms typically encode information specific to the sets of pixels from which they are computed and aid in correcting features such as thruster firing discontinuities (Luger et al., 2016). Cross terms from stars other than the target in question are therefore of little help in the de-trending and can be safely neglected.

One potential pitfall of nPLD is that if the PLD assumptions break down for any of the neighboring targets, the PLD regressors may become contaminated with astrophysical information from that neighbor. This is not in general an issue, since overfitting would only occur if an astrophysical signal in the target star and in its neighbor had the same period and the same phase. However, in the (unlikely) case that PLD fails for the neighboring star and this star happens to be an eclipsing binary or a transiting exoplanet host, it is possible that its transit signals could get imprinted onto the target star’s de-trended light curve, resulting in potential false positive planet detections down the line. The two cases relevant to $K2$ in which PLD could fail in such a way are for saturated stars and stars with bright contaminant sources in their apertures (Luger et al., 2016). As we show in §3.5 below, it is straightforward to adapt PLD to work reliably for saturated stars, thereby circumventing this issue. But while EVEREST 2.0 is more robust against overfitting of crowded stars (§5), highly crowded apertures remain an issue for PLD. When de-trending with nPLD, we therefore select neighboring stars with no other known sources in their apertures that are bright enough ($\Delta Kp<5$) to contaminate the PLD vectors.

3.5. Saturated Stars

As discussed in Paper I, PLD typically fails for stars with saturated pixels, resulting in overfitted light curves with artificially low scatter and suppressed astrophysical information (such as transits with significantly shallower depths). This happens because saturated pixels contain nearly no astrophysical information, as the signal overflows into adjacent pixels in the same column and is ultimately dumped into the pixels at the top and bottom of the bleed trails; these “tail” pixels ultimately contain more astrophysical information than the other pixels in the aperture. Since PLD implicitly assumes that astrophysical information is constant across the aperture, the method breaks down for these stars, and PLD vectors from pixels in the saturated columns become capable of fitting out the astrophysical information in the rest of the aperture.

An obvious workaround is to simply discard pixels in saturated columns from the set of PLD regressors. However, this does not work well in practice, since the remaining regressors often have much lower SNR than the SAP flux and thus have low de-trending power. We instead suggested in Paper I that collapsing saturated columns into single pixels — by co-adding the fluxes in each of the pixels and treating the resulting timeseries as a single PLD pixel — could reduce the effect of saturation, since charge is conserved along the bleed trail. While this ensures PLD does not overfit, it leads to the loss of some of the information about the vertical motion of the stellar point spread function (PSF) across the detector. This leads to significantly poorer de-trending, and therefore we did not employ this method in the first version of the pipeline. However, we find that including the PLD vectors of neighboring stars in the design matrix (i.e., nPLD) effectively restores the information lost when saturated columns are collapsed, leading to high quality de-trended light curves of saturated stars.

In practice, we collapse all columns containing one or more pixels whose flux comes within 10% of (or exceeds) the pixel well depth for the corresponding detector channel during more than 2.5% of the timeseries. We obtained the well depths from Table 13 of the Kepler Instrument Handbook. 111archive.stsci.edu/kepler/manuals/KSCI-19033-001.pdf To avoid aperture losses, we found it necessary to use aperture K2SFF #19 (the largest of the PSF-based apertures) for these stars. We further padded these apertures by two pixels at the top and bottom of saturated columns to ensure that none of the target flux bled out of the aperture; see §5.3.

As an example, in Figure 1 we plot the light curve of EPIC 202063160, a saturated campaign 0 eclipsing binary. The raw light curve is shown at the top and the light curve de-trended with EVEREST 1.0 is shown at the center. Because three of the columns in the aperture contain saturated pixels (right panel), EVEREST 1.0 almost completely fits out the eclipses. With column-collapsed nPLD (bottom), the eclipse is preserved and the instrumental signal is effectively removed.

3.6. Short Cadence

We treat short cadence targets in much the same way as long cadence targets, with the exception that we find it necessary to introduce more breakpoints in the light curves. This is due primarily to computational reasons (short cadence $K2$ light curves are over $10^{5}$ cadences in length; computing Equation (27) for the entire light curve is not feasible). Moreover, we find that noise on sub-30 minute timescales is only properly removed when the size of the light curve segments is kept small. In practice, we find that on the order of 30 breakpoints result in the best de-trending. This might raise concerns of overfitting, but since short cadence light curves contains 30 times more data than long cadence light curves, and we split the latter into two segments, each of the short cadence segments has about twice as many cadences as the long cadence ones.

The major downside of such a large number of segments are the discontinuities that could be introduced at each breakpoint. As before, we overcompute the model into adjacent segments and match the models at the breakpoints, but some de-trended light curves display occasional jumps in either the flux or its derivative.

A second issue with short cadence light curves concerns deep transits and eclipses. As we discussed in Paper I, PLD may attempt to fit out these features if they are not properly masked, since doing so can result in a very large (but spurious) improvement in the photometric precision. With long cadence light curves, transit masking can be done by the user by simply re-computing the model with the appropriate cadences masked, since the transits are sparse and their presence does not significantly affect the cross-validation step. Moreover, outlier clipping usually masks most deep transits anyways, so this is hardly ever a problem. However, that is not the case with short cadence light curves, where transits and eclipses span upwards of fifty contiguous cadences. Since these features are so smooth, they are not flagged as outliers. And since the transit signal is no longer sparse — as it makes up a substantial fraction of the light curve segment — it is far more likely to bias the cross-validation step. In practice, we find that this leads to substantial underfitting of short cadence light curves with deep transits. As $\lambda_{n}$ increases, PLD begins to fit out the transit and the scatter in the validation set grows, forcing the algorithm to select very low values of $\hat{\lambda}_{n}$ and resulting in de-trended light curves that still contain significant instrumental signals.

We therefore explicitly mask all deep transits and eclipses in the short cadence light curves before the cross-validation step. Since only deep transits are likely to bias the cross-validation, and since the number of short cadence light curves in each campaign is relatively small, these can easily be identified by inspection.

3.7. Cross-validation

The principal step in the de-trending process is determining the prior amplitudes $\lambda_{n}$ in Equation (27), which we do by cross-validation. Our method is analogous to that of Paper I, where we performed cross-validation to obtain the optimal number of principal components to regress on. However, here we seek to optimize a three-dimensional function $\sigma_{\mathrm{v}}(\lambda_{1},\lambda_{2},\lambda_{3})$, where $\sigma_{\mathrm{v}}$ is the scatter in the validation set; this is a far more expensive calculation to do. While we could employ a nonlinear optimization algorithm (see below), in the interest of computational speed, we perform a simplification. Since we expect the first order PLD regressors to contain most of the de-trending information, with each successive PLD order providing a small correction term to the fit, we break down the minimization problem into three separate one-dimensional problems. First, we perform cross-validation on the first order PLD model by setting $\lambda_{2}=\lambda_{3}=0$ to obtain the value of $\lambda_{1}$ that minimizes the validation scatter, $\hat{\lambda}_{1}$. We do this by computing the model for each value of $\lambda_{1}$ in a logarithmically-spaced grid with 36 points in the range $[10^{0},10^{18}]$, plus $\lambda_{1}=0$, and select the minimum (details below). We then repeat this process on the second order model by fixing $\lambda_{1}$ at this estimate and keeping $\lambda_{3}=0$. Finally, we solve for $\hat{\lambda}_{3}$ by fixing the first and second order parameters at their optimum values:

[TABLE]

It is important to note that there is no a priori reason that this method should yield the global minimum of $\sigma_{\mathrm{v}}$; in fact, it very likely does not. However, we explicitly allow for $\lambda_{n}=0$ in our grid search, and thus this approximation cannot lead to overfitting, as it will always prefer a lower-order PLD model to one with higher scatter in the validation set.

As a proof of concept, we de-trended a sample of 2,700 randomly selected campaign 6 targets with nPLD using this approximation in the cross-validation step. We then repeated the de-trending by solving for the $\hat{\lambda}_{n}$ using Powell’s method, initializing the solver at different points in the vicinity of the nPLD solution and keeping the solution with the lowest average 6 hr CDPP (combined differential photometric precision; Christiansen et al., 2012) for each target; we dub this method pPLD. In Figure 2 we plot the star-by-star CDPP difference between the two models, $\mathrm{(CDPP_{pPLD}-CDPP_{nPLD})/CDPP_{nPLD}}$. While for some stars the CDPP improves substantially with pPLD, cross-validating with Powell’s method leads to a less than one percent improvement in the CDPP on average. Given that this method is more computationally expensive, we adopt the grid search method outlined above when producing the EVEREST 2.0 catalog.

In addition to being more computationally tractable, there are two major benefits to minimizing $\sigma_{v}$ in this fashion. First, since we perform cross-validation three times (once for each PLD order), we are able to progressively refine the outlier masks (§3.1) and the GP hyperparameters (§3.2) for each target in between cross-validation steps. Second, it allows for some leeway in how we determine the minimum validation scatter. In Paper I, we sought to minimize the median scatter in groups of random 13-cadence segments of the light curve (the validation set). A potential issue with this method is that the noise properties of $K2$ light curves are far from constant over the course of an observing campaign; optimizing the regression based on the median (or mean) validation scatter can still, in principle, lead to overfitting in some segments. While splitting the light curves into segments with similar noise properties (§3.3) helps with this, we also modify the cross-validation process to prevent localized overfitting. For each PLD order $n$ and for each value of $\lambda_{n}$, we split each light curve segment into three roughly equal sections. For each pair of sections, we train the model on them and compute the model prediction in the third section (the validation set). We then compute the scatter $\sigma_{v}$ as the median absolute deviation of the de-trended validation set after removing the GP prediction.

We now have three $\sigma_{v}(\lambda_{n})$ curves, one for each section. In general, the minima of these curves will occur at different values of $\lambda_{n}$, so determining the optimum value $\hat{\lambda}_{n}$ requires a compromise between overfitting and underfitting in the different segments. For each segment, we compute the minimum scatter, find the set of all $\lambda_{n}$ for which $\sigma_{v}(\lambda_{n})$ is within 5% of the minimum, and keep the largest $\lambda_{n}$. We then pick $\hat{\lambda}_{n}$ to be the smallest of these values, provided it is smaller than the largest value of $\lambda_{n}$ at the minima of the three $\sigma_{v}$ curves. This process ensures that $\hat{\lambda}_{n}$ falls between the minima of the $\sigma_{v}$ curves with the smallest and largest value of $\lambda_{n}$, and that it leads to no more than 5% overfitting in one of the segments. We illustrate this procedure in Figure 3, where we show $\sigma_{v}(\lambda_{1})$ for each of the three light curve sections for EPIC 206103150. The red arrows indicate the minimum of each of the curves, and the dashed vertical line indicates the adopted $\hat{\lambda}_{n}$ based on a compromise between slight underfitting in the first two segments and slight overfitting in the third. This results in a more conservative cross-validation process than in Paper I.

4. Results

We de-trended all campaigns 0–8 stars with nPLD to produce the EVEREST 2.0 catalog. Below we report our results, starting with injection/recovery tests and a comparison of rPLD and nPLD, followed by comparisons with other pipelines and the original Kepler light curves. We report most of our results in terms of the proxy 6 hr CDPP of the de-trended light curves, which we calculate in the same way as we did in Paper I: we smooth the light curves with a Savitsky-Golay filter, clip outliers at 5$\sigma$, and compute the median standard deviation in 13-cadence segments, normalized by $\sqrt{13}$.

4.1. Injection Tests

As in Paper I, we perform simple transit injection/recovery tests to ensure our model is not overfitting. For the same sample of 2,700 campaign 6 stars as before, we injected synthetic transits of varying depths at the raw pixel level and attempt to recover them after de-trending with nPLD. We follow the exact same procedure as in §4.1 of Paper I and plot the results in Figure 4 (compare to Figure 6 in Paper I).

Each panel displays two histograms: a blue one, showing the number of transits recovered with a certain depth after de-trending with nPLD, and a red one, corresponding to a control run in which the transits were injected into the already de-trended light curve. Each row corresponds to a different injection depth $D_{0}$ ($10^{-2}$, $10^{-3}$, and $10^{-4}$, from top to bottom), and the $x$ axis in each histogram is the recovered depth $D$ scaled to this value ($D/D_{0}$). The left column corresponds to runs in which the transits were not explicitly masked during de-trending; the right column shows runs in which they were.

As with the previous version of the pipeline, we find a ${\sim}10\%$ bias toward smaller depths for low SNR transits when the transits are not explicitly masked. This is because a small decrease in the transit depth can greatly improve the CDPP of the light curve. Because the PLD regressors are noisy, the method is capable of partially fitting out transits by exploiting linear combinations of white noise in the regressors. This overfitting does not occur for high SNR transits because these are masked during the outlier clipping step.

Conversely, when transits are explicitly masked, there is no bias in the recovered depth; the median $D/D_{0}$ is consistent with unity for all three values of the injected depth. This is the same result we obtained with EVEREST 1.0, and we conclude that our new cross-validation scheme is robust in preventing overfitting when transits are masked. As before, we urge those using EVEREST light curves containing transits or eclipses to re-compute the model with those features masked. This process is quick and straightforward — refer to the EVEREST 2.0 documentation for details. 222https://github.com/rodluger/everest

4.2. rPLD

In Figure 5 we plot a comparison of the CDPP values obtained with rPLD and EVEREST 1.0 for our sample set of 2,700 campaign 6 stars. As before, the $y$ axis corresponds to the normalized relative CDPP of each model, with negative values corresponding to lower CDPP for rPLD. Each star is plotted as a blue dot and the median relative CDPP is indicated as a black line. rPLD outperforms EVEREST 1.0 at all Kepler magnitudes by ${\sim}1-6\%$ on average. However, the scatter at any value of $Kp$ is quite large, and the two models are roughly comparable for bright stars. As we argued in §3.7, the most important feature of rPLD is its increased robustness to local overfitting (see §5).

4.3. nPLD

The greatest improvement in the CDPP comes when neighboring stars’ PLD vectors are included in the design matrix. In Figure 6 we plot the CDPP comparison between nPLD and EVEREST 1.0. nPLD outperforms regular PLD by ${\sim}10-20\%$ on average, with the largest improvement occurring for fainter stars. Faint stars have the noisiest PLD vectors and benefit the most from the inclusion of higher SNR regressors. As we showed in Paper I, regular PLD already comes close to recovering Kepler photometric precision for bright ($Kp\lesssim 13$) stars, so the improvement for these stars is naturally smaller.

4.4. Comparison to Other Pipelines

In this section we compare the EVEREST 2.0 catalog to those produced by other pipelines, beginning with the previous version of EVEREST.

Figure 7 shows the CDPP comparison between EVEREST 2.0 and EVEREST 1.0 for all campaigns 0–8 stars. As in the example shown in Figure 6, EVEREST 2.0 outperforms EVEREST 1.0 by ${\sim}20\%$ for the faintest stars and by ${\sim}10\%$ for $Kp\gtrsim 12$. For $11\lesssim Kp\lesssim 12$ the two pipelines yield comparable results, though regularized regression gives EVEREST 2.0 a slight edge. Below $Kp\approx 11$, $K2$ stars become saturated; these are plotted as red dots, and the median relative CDPP for saturated stars is indicated by the dashed red line. For these stars, EVEREST 1.0 yields much lower CDPP — for $Kp\lesssim 10$, the CDPP is over a factor of two smaller than that of EVEREST 2.0. As we discussed in Paper I, the increased performance of EVEREST 1.0 for saturated stars is spurious, since the astrophysical information content of the pixels is highly variable across the aperture, leading regular PLD to overfit. As we showed in §3.5, EVEREST 2.0 does not overfit saturated stars. In §4.9 below, we show that we approximately recover the Kepler photometric precision for these stars.

In Figure 8 we show the CDPP comparison between EVEREST 2.0 and K2SFF (Vanderburg, 2014; Vanderburg & Johnson, 2014). Our pipeline yields lower average CDPP at all magnitudes, by 40% for faint stars and 10–20% for bright stars. For saturated stars, EVEREST 2.0 outperforms K2SFF by 5–10%.

In Figure 9 we show the comparison to the K2SC PDC light curves (Aigrain et al., 2015, 2016). EVEREST 2.0 yields lower average CDPP at all magnitudes $K_{p}\gtrsim 9$, with a 20–25% improvement for all unsaturated stars. For $K_{p}\approx 9$, K2SC slightly overperforms EVEREST 2.0, but the scatter in the plot is quite large.

We also computed the relative CDPP for each of the campaigns individually; these are plotted in Figures 10–12. The improvement over EVEREST 1.0 is approximately the same for all campaigns despite significant differences in the noise properties and stellar populations across the nine campaigns, showcasing the robustness of the nPLD method. The same is true when compared to K2SFF and K2SC, with the exception of campaigns 0–2, for which EVEREST 2.0 outperforms K2SFF by nearly 50% (i.e., a factor of 2) at all magnitudes.

4.5. Outliers

In addition to yielding lower average CDPP than any other publicly available pipeline, EVEREST yields the largest number of usable data points for any of the $K2$ campaigns to date. These data points generally appear as outliers even in the raw data, since they are highly sensitive to inter- and intra-pixel sensitivity variations. Pipelines that regress on functions of the spacecraft motion alone therefore have trouble de-trending them, resulting in ${\sim}5-10\%$ of the data points being discarded as outliers. In contrast, since PLD uses regressors containing both inter- and intra-pixel sensitivity information, it naturally de-trends data collected during thruster firing events (see Paper I).

To show this, we calculated the number of non-outlier data points per campaign for each of the four pipelines (K2SFF, K2SC, EVEREST 1.0, and EVEREST 2.0). After removing all data points with flagged QUALITY bits 1-9, 11-14, and 16-17, we smoothed each light curve with a Savitsky-Golay filter and performed iterative sigma clipping to remove all 5$\sigma$ outliers. In Figure 13 we plot histograms of the number of remaining data points for all stars in each of the first 9 campaigns. As expected, both EVEREST 1.0 and EVEREST 2.0 have, on average, 100–300 more usable data points than the other two pipelines. This roughly corresponds to the number of thruster firings per campaign, as these happen every 6–12 hours on average.

4.6. Short Cadence

In Figure 14 we plot the relative CDPP distribution of 671 light curves that were observed in both short and long cadence. In order to compute the 6 hr CDPP of short cadence light curves, we first mask outliers and then down-bin to long cadence by taking the mean of every 30 cadences; we then compute the CDPP as usual. We achieve higher average precision in the short cadence light curves by 5-10% for unsaturated stars and by up to 25% for saturated stars. The higher information content in the short cadence light curves — particularly at sub-30 minute timescales — allows EVEREST to better de-trend those stars. We show an example in Figure 15, where we plot the light curves for EPIC 201601162 (raw short cadence at the top, de-trended short cadence in the center). In the bottom panel, we show both the down-binned short cadence de-trended light curve (black) and the long cadence de-trended light curve (red). The CDPP in the short cadence light curve is ${\sim}30\%$ smaller than that of the long cadence light curve.

As we discussed in §3.6, one pitfall of our method for de-trending short cadence targets is the large number of breakpoints (${\sim}30$) we introduce in the light curve when computing the model. Overcomputing the model into adjacent segments and aligning the models at the breakpoints works very well for high SNR light curves but can often fail for light curves that are dominated by photon noise. This is the case for very faint stars ($Kp\gtrsim 17$) observed in short cadence mode, such as EPIC 201831393, which displays visible discontinuities at many of the breakpoints. Low SNR short cadence light curves may also have segments with varying white noise amplitudes due to differences in the value of the $\mathbf{\Lambda}$ prior on the PLD weights. PLD is known to perform poorly when the white noise dominates (Deming et al., 2015); in these cases, it is often desirable to down-bin the light curve and compute the model on a higher SNR signal, then predict onto the original short cadence data. Given the relatively small number of light curves for which this is an issue, we do not do this here.

4.7. Co-trending Basis Vectors

One downside of the algorithm employed by EVEREST is that GP regression has trouble distinguishing between low frequency stellar variability and low frequency instrumental systematics. Many of the EVEREST light curves display a steady rise over the course of the campaign, with hook-like features at the beginning, end, or both. While this does not affect transits or other high-frequency astrophysical signals, it could potentially lead to biases in stellar rotation studies that rely on low frequency modulation in the light curves.

We therefore run a post-processing step on all de-trended light curves to remove these residual instrumental signals. After some experimentation, we decided to use the SysRem algorithm (Tamuz et al., 2005), which identifies and removes signals shared by many light curves. SysRem is similar to PCA but allows for weighting of the input signals and is therefore better suited to dealing with light curves of varying noise properties. Our approach is similar to the presearch data conditioning (PDC) algorithm of the Kepler pipeline, which uses co-trending basis vectors from many stars on the detector to remove common instrumental signals from light curves (Stumpe et al., 2012; Smith et al., 2012).

We apply SysRem to all of the light curves in each campaign, weighting each one by the quadrature sum of the flux measurement errors and the white noise component of its GP. Given known correlations between instrumental signals and spatial position on the detector (e.g., Petigura & Marcy, 2012; Wang et al., 2016), we attempted to apply SysRem individually on each CCD module, but found that the recovered signals were often dominated by astrophysical variability originating from the brightest star(s) in each module. We found that computing the SysRem signals for the entire detector alleviated this issue without compromising the de-trending power of the method.

We separately apply SysRem to each light curve segment, obtaining one co-trending basis vector (CBV) for each segment of each campaign. We then subtract a linear fit of this CBV from each light curve and repeat the procedure to obtain additional CBVs for each segment of each campaign. In order to prevent the CBVs from fitting out or introducing high frequency signals in the light curves, we aggressively smooth them with a third order, 1000-cadence Savitsky-Golay filter. The results are shown in Figure 16. The first CBV for each segment of each campaign is plotted in blue, and the second in red. As expected, the first CBV in each segment is dominated by a linear trend with S-like hooks on either end. The second CBVs are predominantly quadratic or cubic. The remaining CBVs (not shown) are dominated by higher order trends.

For the purposes of generating the EVEREST 2.0 catalog, we perform ordinary least squares regression to fit all de-trended light curves using only the first CBV (blue curves in the figure). We find that fitting with additional CBVs often helps to remove additional systematics — in particular the hook-like features mentioned above — but may lead to overfitting of true astrophysical variability in some light curves. We therefore include all 5 CBVs of each segment of each campaign in the EVEREST 2.0 FITS files so that targets may be further corrected by the user on a case-by-case basis (see §6).

4.8. Sample Light Curves

In Figures 17 and 18 we show two sample light curves de-trended with both K2SFF and EVEREST 2.0. Figure 17 shows EPIC 201345483, a faint $Kp=15$ planet host. The CDPP of the EVEREST light curve (bottom) is a factor of 2.4 lower than that of K2SFF (top), and the light curve has visibly fewer outliers. The folded transit is shown at the right and is similarly less noisy in the EVEREST light curve.

In contrast, Figure 18 shows a very bright ($Kp=10$) saturated planet candidate host, EPIC 201862715. With the column-collapsing scheme and the inclusion of neighboring PLD vectors, EVEREST 2.0 is able to achieve a factor of 5 lower CDPP than K2SFF, as well as considerably fewer outliers.

4.9. Comparison to Kepler

In Figures 19 and 20 we compare the EVEREST 2.0 photometric precision to that of the original Kepler mission. Figure 19 shows the CDPP as a function of $Kp$ for each of the first 9 $K2$ campaigns and Figure 20 shows the comparison for all $K2$ stars. Because of differences in the raw photometric precision and in the stellar populations across the campaigns, the results are variable, but for all campaigns except 0, 2 and 7, we recover the original Kepler precision down to at least $Kp=14$. For campaigns 1, 5, and 6, we recover the Kepler precision down to at least $Kp=15$. This also applies to saturated stars, though in some campaigns the EVEREST 2.0 CDPP is slightly higher for these stars. For stars dimmer than $Kp=15$, the EVEREST 2.0 CDPP is within a few tens of percent — or less — than that of Kepler.

5. Additional Remarks

5.1. Variable Stars

In Paper I we discussed how EVEREST often fails to properly de-trend extremely variable stars, such as RR Lyrae variables and very short period eclipsing binaries, causing overfitting in many of these light curves. We attributed this to there being too much power in the GP model, which captured both the astrophysical and the instrumental variability, resulting in an improperly optimized PLD model. However, after considerable experimentation, we found that this behavior stemmed in large part from our cross-validation scheme. In Paper I, our cross-validation sets were 13 cadences (6.5 hours) long, and we sought to minimize the median (proxy) CDPP of all such sets. We chose this timescale because it is roughly the duration of a typical transit, and in practice it worked well to minimize overfitting of transits. However, the CDPP as defined in §4 — the normalized standard deviation in 13 cadence segments after the application of a high-pass filter — is not an adequate metric of the photometric precision for stars that are intrinsically variable on similarly short timescales. In other words, RR Lyrae and other extreme variables are dominated by astrophysical variability on the timescale at which the CDPP is computed, and therefore minimizing the CDPP is a recipe for overfitting.

In EVEREST 2.0, we modified our cross-validation scheme (§3.7) in two important ways: we increased the average size of the validation sets to ${\sim}500$ cadences, and we minimized the validation scatter after the subtraction of a properly optimized GP model. This helps to ensure that we minimize only the instrumental component of the noise. In order to assess the performance of this new scheme, we visually inspected the de-trended light curves of 100 RR Lyrae stars in campaigns 0–4, chosen as the targets with the highest probability of being RRab stars according to the K2VARCAT catalog (Armstrong et al., 2016). Among the EVEREST 1.0 light curves, 92/100 had visibly damped oscillation amplitudes or clearly overfitted stellar variability features. In contrast, only 44/100 EVEREST 2.0 RR Lyrae stars showed any signs of overfitting. For 8 of these, the overfitting occurred only with the inclusion of the second or third order PLD models. While there are still issues with how the pipeline handles extremely variable stars, the improvement over version 1.0 is substantial.

For comparison, we visually inspected the same stars in the K2SFF and K2SC catalogs. 96/100 K2SFF light curves of RRab stars showed signs of overfitting (dampened oscillation amplitudes, distorted astrophysical signal, or significant de-trending artefacts), while only 12/65 (18%) of K2SC light curves appeared to be incorrectly de-trended (35 of the 100 stars were in campaigns 0–2 and are not present in the K2SC catalog). K2SC likely overperforms the other pipelines for these stars because of its robust GP optimization scheme. Since we optimize the GP and de-trend the light curve in separate steps, the covariance matrix we use is often an improper approximation to the true covariance of the astrophysical signal. Progressive optimization of the GP (§3.2) helps to maximize the amount of astrophysical information captured by the GP model, but a better procedure would be to simultaneously fit for the GP hyperparameters and the PLD coefficients. This would ensure that the PLD model captures only instrumental signals and the GP model captures only astrophysical signals. However, such a method is computationally intractable, since the problem would no longer be linear. We therefore settle for our linear method, which works extremely well for stars that do not exhibit extreme variability and has a ${\sim}50\%$ success rate for those that do. We encourage readers interested in RR Lyrae and other extreme variables to inspect the light curves of the different pipelines on a target-by-target basis.

5.2. Crowded Apertures

In Paper I we also showed how PLD is not suited to de-trend stars in excessively crowded fields, since the algorithm will often use pixels from contaminant sources to fit out the target star’s astrophysical signals. Directly addressing this issue is far less straightforward than mitigating the effects of saturation (§3.5) or extreme astrophysical variability (§5.1), as it likely requires accurate modeling and subtraction of the PSFs of contaminant sources. However, in practice we find that our new de-trending algorithm is far more robust to overfitting when contaminant sources are present. This is in part due to the inclusion of PLD regressors from neighboring stars, which reduce the contribution from potentially contaminated pixels in the target aperture. Our new cross-validation scheme (§3.7) also helps guard against the effects of crowding, since it is better at preventing localized overfitting. Because of the large amplitude of the $K2$ spacecraft drift, the amount of contamination in the target aperture often varies considerably over the duration of the campaign, as nearby sources move in and out of the aperture. In EVEREST 1.0 light curves, this resulted in a time-variable de-trending power — and time-variable overfitting — for crowded targets. As we discussed in §3.7, our new cross-validation method strongly disfavors this behavior.

In Figure 21 we show the light curves of five eclipsing binaries whose crowded apertures resulted in severe overfitting in the EVEREST 1.0 catalog. Shown are the light curves (flux versus time) folded on the period of the binary and centered on the primary eclipse for EVEREST 1.0 (left) and EVEREST 2.0 (center). The right panel shows the Palomar Optical Sky Survey333http://stdatu.stsci.edu/cgi-bin/dss_form (POSS II) red filter image of the target postage stamp, where contaminant sources are clearly visible in each of the apertures (indicated by the red contours).

The first three targets (EPIC 202072978, EPIC 202103762, and EPIC 218803648) have relatively bright contaminant sources centered on or near the edge of the aperture. In Paper I we explained how this was the “worst case” scenario for crowding, as it results in the largest spatial variation of astrophysical information from the target across the aperture; as expected, EVEREST 1.0 almost completely fits out their eclipses. EVEREST 2.0, on the other hand, achieves high de-trending power while preserving the eclipse shapes and depths seen in the raw light curves. The fourth target (EPIC 202733088) has a bright nearby source inside its aperture, separated by less than a pixel. We discussed in Paper I how targets with co-located contaminants are typically unaffected by PLD overfitting, since the ratio of the two PSFs is roughly constant everywhere. However, in this case the contaminant is sufficiently detached to result in severe overfitting in the EVEREST 1.0 light curve. As before, EVEREST 2.0 correctly de-trends the light curve with no overfitting. The final target (EPIC 211685048) is an eclipsing binary with deep transits that are completely fit out in both versions of the pipeline. In this case, two roughly equal magnitude stars are fully contained in the aperture and separated by ${\sim}3$ pixels, causing PLD to fail despite the modifications to the algorithm.

As in Paper I, we urge caution when using the EVEREST light curves of crowded targets, but note that in most of the cases known to us in which the previous version overfitted transits and eclipses, EVEREST 2.0 succeeds in removing the instrumental noise without overfitting astrophysical information. Future updates to the pipeline will include a more careful aperture selection, which can mitigate the effects of crowding for targets like EPIC 211685048.

5.3. Aperture Losses

As we have demonstrated, column-collapsed PLD works extremely well for saturated stars, provided all of the target flux in the saturated columns is used. With this in mind, we used larger apertures for these stars and extended them by two pixels at the top and bottom of all saturated columns (§3.5). However, for some extremely bright ($Kp\lesssim 8.5$) targets, bleed trails in the saturated columns can extend past the edge of the target postage stamp. Since most of the astrophysical information in saturated columns is contained in the two pixels at the top and bottom, this can lead to substantial information loss. In this case, the column collapsing procedure fails to satisfy the PLD constraint that all pixels should have the same fractional astrophysical signal strength, and overfitting can occur. This is the case for EPIC 210703831, a campaign 4 $Kp=8.1$ star, whose de-trended light curve displays several discontinuities. These occur because spacecraft drift results in aperture losses only during parts of the campaign. We therefore encourage those using the EVEREST catalog to inspect the postage stamps of extremely bright stars to ensure that aperture losses are not present.

6. Using EVEREST

All EVEREST 2.0 $K2$ light curves are available for download as FITS files from MAST.444https://archive.stsci.edu/prepds/everest As with the previous version, we urge users to interface with the catalog via the Python code, which can be installed following the instructions at https://github.com/rodluger/everest. A detailed description of the data products and how to use the code is available on the github page. Below, we provide a brief outline of the EVEREST resources available online.

6.1. FITS Files

Each FITS file contains six extensions. The primary (index 0) extension consists of a header with miscellaneous target information copied from the $K2$ target pixel file (TPF). The second (index 1) extension contains a header and a binary table. The header stores miscellaneous information about the target and the settings used in the de-trending, such as the GP hyperparameters and information on the neighboring targets used in the regression. The binary table stores arrays corresponding to the cadence number (CADN), the timestamp (TIME), the raw SAP flux (FRAW), the raw SAP flux errors (FRAW_ERR), the PLD-de-trended flux (FLUX), the five CBV regressors (CBV01 – CBV05), and the de-trended flux with the CBV correction (FCOR). This extension also includes the original $K2$ QUALITY bit array, with four additional bits that signal cadences that were masked when computing the model:

[TABLE]

When using the EVEREST light curves for science purposes, it is important to properly reject these outliers. Data flagged with bit 23 is typically affected by cosmic rays or detector anomalies, and can usually be ignored. In some cases, however, deep transits or eclipses can be mistaken for detector anomalies and are incorrectly flagged in the original TPF, so a visual inspection of the light curves is recommended. Data flagged with bit 24, on the other hand, is missing in the original TPF and can thus be safely ignored. Finally, bit 25 usually corresponds to instrumental outliers that were not properly de-trended with PLD. However, transits, flares, and other short timescale astrophysical features will also be flagged with this bit, so this data should not be blindly excluded when performing transit searches.

The four remaining FITS extensions include the PLD regressors, the target aperture, and additional information used internally by the Python code.

6.2. Data Validation Summaries

Each target in the catalog also has an associated data validation summary (DVS), a PDF document showing the raw, de-trended, and CBV-corrected light curves, as well as cross-validation diagnostic plots such as those in Figure 3 and a high resolution POSS image of the target aperture like those shown in Figure 21.

6.3. Python Code

The EVEREST code can be installed using the package managing system pip:

⬇

pip install everest-pipeline

or directly from source by following the instructions on the github page. The primary way of interfacing with the catalog is to instantiate an Everest object:

⬇

import everest

star = everest.Everest(EPIC)

where EPIC is the EPIC number of the target. These lines download the target’s FITS file and populate the star object with the de-trending information and light curve arrays (time, flux, fcor, etc.). The QUALITY flags 23, 24, 25, and 27 are used to generate the badmask, nanmask, outmask, and transitmask arrays, respectively; these arrays contain the indices of all data points whose corresponding bit is flagged.

As we discussed in §4.1, it is important that transits are properly masked when computing the PLD model; otherwise, slight overfitting may occur. Users should therefore always mask known transits and re-compute the PLD model, as follows:

⬇

star.mask_planet(time, period, duration)

star.compute()

where time is the time of first transit (in units of $\mathrm{BJD}-2454833$), period is the transit period (in days), and duration is the full transit duration (also in days). Re-computing the model takes no more than a few seconds for long cadence light curves.

Users can also easily change the number of CBVs used to correct the light curve:

⬇

star.cbv_num = n

star.compute()

where n is the desired number of CBVs (0 – 5).

The code also implements various visualization routines, which are described in the documentation linked on the github page.

7. Conclusions

We have presented EVEREST 2.0, an update to the EVEREST pipeline (Luger et al., 2016) for removing instrumental noise from $K2$ photometry. In version 1.0, we constructed a linear model from the principal components of products of the fractional pixel fluxes in the aperture of each star, a variant of a method known as pixel level decorrelation (PLD, Deming et al., 2015). Here, we regress on all PLD vectors, imposing Gaussian priors on the model weights to prevent overfitting. We additionally include the PLD vectors of bright neighboring stars to increase the signal-to-noise ratio of the regressors and enhance the predictive power of the model. We developed a fast gaussian process (GP) regression scheme to de-trend all stars in the $K2$ catalog, achieving lower combined differential photometric precision (CDPP) than in version 1.0, by ${\sim}10\%$ for bright stars and ${\sim}20\%$ for faint stars. We also adapted PLD to work for saturated stars, yielding comparable de-trending power, and stars observed in short cadence mode, yielding higher photometric precision on 6 hr timescales than their long cadence counterparts. We further find that the inclusion of neighboring PLD vectors and a more conservative cross-validation scheme enhance the pipeline’s robustness to overfitting, particularly for highly variable stars.

EVEREST 2.0 light curves have higher photometric precision than the two other publicly available catalogs, K2SFF (Vanderburg & Johnson, 2014) and K2SC (Aigrain et al., 2016), at all $Kp$ magnitudes. For faint stars, EVEREST 2.0 has ${\sim}40\%$ lower CDPP than K2SFF and ${\sim}25\%$ lower CDPP than K2SC; for bright unsaturated stars, the CDPP improvement is ${\sim}20\%$ compared to both pipelines. For saturated stars, EVEREST outperforms both pipelines, but by a smaller margin. We also find that EVEREST light curves have, on average, 100–300 fewer outliers than those of other pipelines, owing primarily to the ability of PLD to correct data collected during thruster firing events.

When compared to the original Kepler mission, EVEREST 2.0 recovers Kepler photometry on average to $Kp\approx 14.5$, and past $Kp=15$ for some campaigns. For dimmer stars, the CDPP is within a few tens of percent of that of Kepler. EVEREST light curves should thus enable continued high-precision transiting exoplanet and stellar variability science for the vast majority of $K2$ stars as if they had been observed by the original four-wheeled mission.

The EVEREST 2.0 catalog of de-trended light curves is publicly available at https://archive.stsci.edu/prepds/everest. As with the previous version of the code, EVEREST 2.0 is open source under the MIT license and available at https://github.com/rodluger/everest. The reader is encouraged to use this code to interface with the EVEREST catalog and to customize the de-trending of targets of interest, particularly for masking transits to remove biases in the depth due to overfitting. We have implemented each of the PLD models discussed above (rPLD, nPLD, pPLD) in a general framework that can be adapted to different missions, including Kepler and the upcoming TESS. For more information, refer to the documentation linked on the github page.

We would like to thank Benjamin Pope and Tsevi Mazeh for their useful comments and suggestions. R.L., R.B., and E.A. acknowledge support from NASA grant NNX14AK26G and from the NASA Astrobiology Institute’s Virtual Planetary Laboratory Lead Team, funded through the NASA Astrobiology Institute under solicitation NNH12ZDA002C and Cooperative Agreement Number NNA13AA93A. E.A. acknowledges additional support from NASA grants NNX13AF20G and NNX13AF62G. E.K. acknowledges support from an NSF Graduate Student Research Fellowship. This research used the advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington.