TRANSIT DETECTION IN THE MEarth SURVEY OF NEARBY M DWARFS: BRIDGING THE CLEAN-FIRST, SEARCH-LATER DIVIDE

The MEarth CCDs behave slightly nonlinearly at all count levels, increasing up from a 1%–2% nonlinearity at half of the detector full well up to 3%–4% near the onset of saturation. Because we often need to use comparison stars with different magnitudes than our target star, we must account for this nonlinearity. When setting exposure times, we avoid surpassing 50% of the detector's full well, and estimate a correction for the nonlinearity using sets of daytime dome flats taken with different exposure times. With these measures in place, we see no evidence that nonlinearity limits our photometric performance.

We scale dark exposures taken at the end of each night to remove some of the CCD dark current. However, until 2011 we operated our Peltier-cooled detectors at −20°C to −15°C, and at these warm temperatures they showed significant persistence. That is, images of bright stars would persist as excess localized dark current in subsequent images, slowly decaying with a half-hour timescale below an initial 1% fraction of the original fluence. This was a significant source of systematics: stars in incoming exposures could land on the same pixels as persistent ghost stars from previous exposures, thus gaining a hidden amount of flux that depended on how recently and strongly those pixels were illuminated. As this effect depends on the entire recent two-dimensional illumination history of the detector, correcting for it would be extremely complicated. We partially mitigated the persistence by ensuring that different target stars were observed on different regions of the detector, but could not eliminate problems due to overlap with comparison stars or due to changes in cadence. Updating our camera housings in 2011, we now operate at −30°C where the amplitude of the persistence is lower. We also adopt a detector preflash before each exposure; this increases the overall dark current but suppresses localized persistent images. Between the lower temperature and this preflash step, persistence no longer has a substantial effect on MEarth photometry.

We gather flat-fields at evening and morning twilight, typically eight per telescope per twilight, with empirically set exposure times estimated using the equations of Tyson & Gal (1993). To average out large-scale gradients in the illumination, we always take adjacent pairs of flats on opposite sides of the meridian, which has the effect of rotating the whole optical system relative to the sky (thanks to our German Equatorial mounts). Because our optical system shows high levels of centrally concentrated scattered light (10%–15% of sky before 2011 and <5% after; see Table 1) that corrupts the large-scale structure in twilight flat exposures, we estimate the sensitivity in the detector plane in two steps. First, we estimate a small-scale sensitivity map that accounts for dust donuts and pixel-to-pixel variations in the detector sensitivity by filtering out large-scale structure from the combined twilight flats. Second, we derive a large-scale map from dithered photometry of dense star fields to account for the non-uniform illumination across the field of view. Additionally, our camera's leaf shutter takes a finite time to open and close, resulting in a varying exposure time across the field of view (on a 1 s exposure, the amplitude of this effect is 5%); we apply a shutter correction estimated from sets of twilight flats. Altogether, our flat-fielding procedure achieves a precision of 1% across the entire detector.

However, because we hope to perform photometry down to the level of 0.1%, this 1% knowledge of the sensitivity across the field is still imperfect and will inevitably be a source of systematics in our light curves. One unavoidable problem is that our German Equatorial mounts require the detector to flip 180° when crossing the meridian. In light curves, this causes offsets as large as 1% between opposite sides of the meridian, as stars sample different regions of the large-scale sensitivity of the camera. Notably, the step-function morphology of this systematic can mimic a transit ingress or egress. In addition to this "meridian flip" problem, we achieve a blind rms pointing accuracy of 60''–120''. To improve on this, at each pointing we take a short binned image and use its astrometric solution to nudge the telescope to the correct pointing before science exposures, with a random error typically of 1''–2''. This minimizes the impact of these pointing errors, but does not completely remove the problem of stars sampling different pixels on an imperfectly flat-fielded detector.

We perform aperture photometry on all sources in the field of view. For each exposure, we derive a differential photometric correction from point sources in the field using an iterative, weighted, clipped fit that excludes variable stars from the comparison sample (see Irwin et al. 2007 for details). We calculate a theoretical uncertainty estimate σ_the(t), in magnitudes,⁵ for each point:

$$\begin{equation} \sigma _{\rm the}(t) = \frac{2.5}{\ln 10} \times \frac{ \sqrt{N_{\gamma } + \sigma _{\rm {sky}}^2 + \sigma _{\rm scint}^2 + \sigma _{\rm comp}^2}}{N_{\gamma }}, \end{equation} \tag{ 1 }$$

where N_γ is the number of photons from the source, σ²_sky is an empirically determined sky noise estimate for the photometric aperture that includes read and dark noise, σ²_scint is the anticipated scintillation noise (Young 1967), and σ²_comp accounts for the uncertainty in the comparison star solution. In some MEarth fields with very few comparisons, the σ²_comp term can be a significant contribution to the overall uncertainty.

A crucial assumption of this differential photometry procedure is that atmospheric or instrumental flux losses are exactly mirrored between target and comparison stars. Our wide 715–1000 nm bandpass overlaps strong telluric absorption features due to water vapor, so as the level of precipitable water vapor (PWV) changes in the column over our telescopes, their effective wavelength response will also change. Red stars will experience a larger share of this time-variable PWV-induced extinction than stars that are blue in this wavelength range. As a typical MEarth field consists of one very red target star (median target r − J = 3.8)⁶ among much bluer comparison stars (median comparison r − J = 1.3), most MEarth M dwarfs exhibit systematic trends caused by this second-order extinction effect. This PWV problem has been noted before as a limitation for cool objects observed in the near-IR (Bailer-Jones & Lamm 2003; Blake et al. 2008). Recently, Blake & Shaw (2011) showed that GPS water vapor monitoring could be used to correct for the influence of PWV variations, improving both relative and absolute photometric accuracy of SDSS red star photometry.

While we do not have a GPS water vapor monitor, we can track the impact of PWV variations on MEarth photometry using the ensemble of observations we gather each night, observations of red stars in fields of blue comparisons. Figure 1 shows all of the M dwarf light curves gathered by MEarth over one week, after applying basic differential photometry. These light curves (of different M dwarfs observed on different telescopes) move up and down in unison, reflecting water vapor changes in the atmosphere they all share. These trends correlate strongly with ground-level humidity and ambient sky temperature, which are rough tracers of PWV in the overlying column (Figure 2). As PWV variations within a night can mimic transit signals (e.g., the first panel of Figure 1), we must account for this effect when searching for planets. Fortunately, because these trends are shared among all our targets, we can estimate a "common mode" timeseries from the data themselves and use it to correct for these trends (see Section 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The S(t) term in Equation (2) enables us to include systematic trends that show clear correlations with externally measured variables. We construct S(t) as a linear combination of N_sys relevant external templates:

$$\begin{equation} S(t) = \sum _{j=1}^{N_{\rm sys}}s_{j}E_{j}(t). \end{equation} \tag{ 3 }$$

Here, E_j(t) represent timeseries of the external variables, sampled at the times as the photometric observations, and the s_j are systematics coefficients. For MEarth, at a bare minimum, we include N_sys = 6 terms in this sum: the "common mode," the "meridian flip," and the x and the y pixel positions on either side of the meridian.

• ECM(t).  
  The common mode template is constructed from the ensemble of raw M dwarf light curves from all telescopes and accounts for photometric trends that are shared in all MEarth M dwarf photometry (due to PWV variations, see Figures 1 and 2). The effect is stronger for redder stars; for MEarth targets, the best-fit values of the coefficient s_CM correlates with stellar r−J color.
• Emerid(t).  
  To account for stars sampling different regions of the detector when observing at positive or negative hour angles with MEarth's German Equatorial mounts, we include a "meridian flip" template. This template is simply defined as 0 for observations taken in one orientation and 1 for observations in the other, thus allowing light curves on two sides of the meridian flip to have different baselines.⁸
• Ex, i(t) and Ey, i(t) for i = 0,1.  
  The pixel position templates are simply the x and y centroids of the target star on the detector, with their medians subtracted. Two sets are required, one for each side of the meridian. Correlations with these templates could arise as pointing errors allow a star to drift over uncorrected small-scale features in the sensitivity of the detector (e.g., transient dust donuts).

Additional external variables may also be used as systematics templates, such as FWHM or airmass. With MEarth, we find that these variables are correlated with the photometry for only a few fields and are usually excluded.

The V(t) term in Equation (2) describes the variability of the star throughout one night, independent of the presence of a transiting planet. Such variability includes fluctuations due to rotating spots (smoothly varying on the 0.1 to 100 day timescale of the star's rotation period) and flares (impulsively appearing, with a decay timescale typically of hours). The morphology of this variability can be quite complicated; we use a simplified model to capture its key features, writing

$$\begin{eqnarray} V(t) &=& v_{{\rm night}} + v_{\sin }\sin \left(\frac{2\pi t}{v_{P}}\right)\nonumber \\ && + v_{\cos }\cos \left(\frac{2\pi t}{v_{P}}\right)+\sum _{j=1}^{N_{\rm flares}} f_{j}(t). \end{eqnarray} \tag{ 4 }$$

The first term v_night allows each night to have its own baseline flux level. By itself, this term can capture most of the variability from stars with long rotation periods, where the flux modulation from starspots smoothly varies over timescales much longer than one night. By fitting for a different v_night for each night we can piece together the variability of the star on timescales >1 day as a series of scaled step functions. The harmonic v_sin and v_cos terms capture variability with period of v_P and become especially important for stars with shorter rotation periods. Because we fit a separate v_night for each night, there can be substantial degeneracy between the harmonic terms and the nightly offsets, especially for slowly rotating stars. We discuss this issue, as well as how we estimate v_P in Section 4.2. Although we only include one harmonic of the fundamental period v_P in these sinusoidal terms, additional harmonics could be included if the data warranted them.

The final term in Equation (4) includes contributions from N_flares hypothetical stellar flares f_j(t) that may or may not be present within the night. Flares are suppressed in MEarth's relatively red bandpass but not completely eliminated (see Tofflemire et al. 2012). The main purpose of the flare term is to identify those nights of photometry that may be corrupted due to the presence of flares. While we could in principle model a night that contained both flares and a transit, we find this to be very difficult in practice, due to the morphological complexity that flares sometimes exhibit (e.g., Kowalski et al. 2010; Schmidt et al. 2012). Rather, on each night we model simple hypothetical flares as fast rising and exponentially decaying, and perform a grid search over the start time and decay timescale. If any flares have amplitudes that are detected at >4σ, we excise that night of data from our planet search. These cuts dramatically reduce planetary false positives due to flaring activity (e.g., confusing the start of the flare with the egress of a transit), with the meager cost of ignoring 3.6% of MEarth's observations. Because both planetary transits and flares are rare in MEarth data, and their overlap even more so, the losses from this strategy are small.

The last term in Equation (2), P(t), includes the signal of a hypothetical transiting planet. We model transits as having infinitely short ingress/egress times and ignore the effects of limb darkening on the host star, so transits appear as simple boxcars. In this section, we are interested only in assessing the significance of a single transit event falling within a single night, not a periodic train of transits. With these simplifications, a lone planetary eclipse signal is completely described by a transit epoch p_E, a transit duration p_T, and a transit depth D. The signal is then simply

$$\begin{equation} P(t)= \left\lbrace \begin{array}{@{}ll}D & \mbox{if $|t - p_E| < p_T/2$} \\ \ms 0 &\mbox{otherwise}. \end{array} \right. \end{equation} \tag{ 5 }$$

This model includes only one eclipse event per night. We discuss combining these lone eclipses into periodic transit candidates in Section 4.3.

A crucial component of the model is σ(t), the photometric uncertainty of each observation. Our theoretical uncertainty estimate for a given data point σ_the(t) is a lower limit on the true uncertainty. To express this fact, we introduce a noise rescaling parameter r_{σ, w} such that

$$\begin{equation} \sigma (t) = r_{\sigma , w}\sigma _{\rm the}(t), \end{equation} \tag{ 6 }$$

where r_{σ, w} ⩾ 1. The subscript w emphasizes that this is a white noise rescaling parameter that does not account for correlations between nearby data points. If left unmodeled, such red noise could substantially bias a transit's detection significance (Pont et al. 2006); we discuss a correction for red noise in Section 4.2.6.

We describe each of the N_obs photometric observations d(t_i) within a particular night as being drawn from a Gaussian distribution centered on m(t_i) and with a variance of σ(t_i)². Assuming the observations to be independent, the likelihood can be written as

$$\begin{eqnarray*} &&P(\mathbb {D}|\mathbb {M})= \prod _{i=1}^{N_{\rm obs}}\frac{1}{\sqrt{2\pi }\sigma (t_i)}\exp {\left[-\frac{1}{2}\left(\frac{d(t_i)- m(t_i)}{\sigma (t_i)}\right)^2\right]}. \nonumber \end{eqnarray*}$$

Taking the logarithm, substituting Equation (6), and defining

$$\begin{equation} \chi ^2 = \sum _{i=1}^{N_{\rm obs}}\left[\frac{d(t_i)- m(t_i)}{\sigma _{\rm the}(t_i)}\right]^2, \end{equation} \tag{ 8 }$$

we find that the (log) likelihood simplifies to

$$\begin{equation} \ln P(\mathbb {D}|\mathbb {M}) = - N_{\rm obs} \ln r_{\sigma , w} - \frac{\chi ^2}{2r_{\sigma , w}^2} + {\rm constant}, \end{equation} \tag{ 9 }$$

where we have only explicitly included terms that depend on the parameters of the model. For fixed r_{σ, w}, maximizing Equation (9) is equivalent to minimizing the commonly used χ² figure of merit.

For any one star, a particular night of MEarth photometry may contain roughly as many light curve points as there are parameters in our model. As such, the likelihood $P(\mathbb {D}|\mathbb {M})$ from one night of data only very weakly constrains the parameter space. But of course, each night of MEarth observations is just one of many nights spanning an entire season, and we should use this season-long information when investigating a single night. To implement this holistic awareness of the context provided by a large pool of observations, we generate probability distributions for various parameters by looking at the whole season of data. We then apply them as priors $P(\mathbb {M})$ on the parameters for an individual night.

By construction, the most important parameters of our model are linear parameters. The conditional likelihood of linear parameters (a slice through $P(\mathbb {D}|\mathbb {M})$ with other parameters fixed) has a Gaussian form. For marginalization, it proves quite useful for the priors to be conjugate to this shape—that is, also to take on a Gaussian form. Referring to these linear parameters with the vector $\mathbf {c} =\lbrace D, v_{\rm night}, v_{\sin }, v_{\cos }, s_{\rm CM}, s_{\rm merid}, s_{x,i}, s_{y,i}, s_{\rm other?}\rbrace$, we parameterize the prior $P(\mathbb {M})$ as being proportional to a Gaussian distribution in c_j that is centered on an expectation value $\overline{c_{j}}$ and with a variance of $\pi _{c_{j}}^2$. Multiplying the independent priors for the N_coef coefficients and defining

$$\begin{equation} \Phi ^2 = \sum _{j=1}^{N_{\rm coef}}\left(\frac{c_j - \overline{c_j}}{\pi _{c_j}}\right)^2 \end{equation} \tag{ 10 }$$

leads to a term in the prior that looks like

$$\begin{equation} \ln P(\mathbb {M}) = -\frac{1}{2}\Phi ^2+ \cdots \end{equation} \tag{ 11 }$$

The similarity in the form of Φ² to χ² is the reason that the use of conjugate Gaussian priors is often described along the lines of "adding a prior as an extra data point in the χ² sum," because the effect is identical in the overall posterior. In this framework, the smaller values of $\pi _{c_j}$ provide tighter constraints on the parameter; we could express a flat, non-informative prior for a particular c_j by choosing a large value of $\pi _{c_j}$. We set π_D = ∞, giving a flat prior on the transit depth. Note that we allow negative transit depths (i.e., "anti-transits") to avoid skewing the null distribution of transit depths away from 0.

For most of the remaining linear parameters, we take the values of $\overline{c_{j}}$ and $\pi _{c_{j}}$ directly from the results of a simultaneous fit to the star's entire season of observations. In this prior-generating season-long fit, we fit the season-long light curve with a modified version of Equation (2) that excludes both the f_j(t) term from flares and the P(t) term from hypothetical planets. To immunize against these unmodeled flares and eclipses, we perform the fit with 4σ clipping. We prefer to explain as much of the long-term variability as possible with the harmonic terms, so we fit first including only these terms in V(t). Then, fixing the values of v_sin and v_cos, we fit again with one v_{night, j} free parameter for each night represented within the season. Thus, the values of v_{night, j} then represent the deviation of the nightly flux level from a baseline sinusoidal model.

Now, for the single night flux baseline parameter v_night, we set $\pi _{v_{\rm night}}$ equal to 1.48 × MAD (median absolute deviation) of the ensemble of v_{night, j} values from the season fit. Stars that vary unpredictably from night to night will have a broad prior for v_night, thus requiring more data within a night to determine its baseline level. Conversely, stars that remain constant from night to night or have variability that is well described by a sinusoid will have a very tight prior.

To understand the impact of $\pi _{v_{\rm night}}$, imagine the following hypothetical scenario: a night in which MEarth gathered only one observation of a star, and that observation happened to fall in the middle of a transit with a 0.01 mag depth. With what significance could we detect this transit? If $\pi _{v_{\rm night}} = 0.01$ mag, then the detection significance would be at most 1σ. But if $\pi _{v_{\rm night}} = 0.001$ mag, then the transit could in principle be detected at high significance with only the single data point, provided the photon noise limit for the observation was sufficiently precise.

We note that s_{x, i} and s_{y, i}, the coefficients for the x and y pixel position templates, would not be expected to be constant throughout a season. These terms are designed to account for flat-fielding errors, which could easily change from week to week or month to month. As such, we do not take c_j and π_j from the season-wide fit for these parameters. Rather, we fix c_j = 0 and π_j = 0.001 for all four of these parameters. This has the desired effect that an apparent 0.005 mag transit event that is associated with simultaneous 5 pixel shift away from the star's mean position on the detector would not be considered as a significant event.

The most significant nonlinear parameter is the white noise rescaling parameter r_{σ, w}. In the season-long fit, Equation (9) indicates that $P(\mathbb {M}|\mathbb {D})$ would have a shape of

$$\begin{equation} \ln P(r_{\sigma , w}|\mathbb {D}) = - N_{\rm sea} \ln r_{\sigma , w} - \frac{\chi ^2_{\rm sea}}{2r_{\sigma , w}^2}, \end{equation} \tag{ 12 }$$

where χ²_sea is the season-long χ² from the N_sea observations in the ensemble fit. This is maximized when $\overline{r} = \sqrt{\chi ^2_{\rm sea}/N_{\rm sea}}$. We want the nightly prior on r_{σ, w} to push it toward $\overline{r}$, but we also want to provide enough flexibility that nights that are substantially better or worse than typical can be identified as such. To implement this, we mimic the shape of the season-long probability distribution but artificially broaden it with an effective weighting coefficient N_eff. Propagating this loose prior

$$\begin{equation} \ln P(\mathbb {M}) = \cdots - N_{\rm eff} \ln r_{\sigma , w} - \frac{N_{\rm eff}\overline{r}^2}{2r_{\sigma , w}^2} + \cdots \end{equation} \tag{ 13 }$$

into the posterior for an individual night, the MAP value of r_{σ, w} will be

$$\begin{equation} r_{\sigma , w} = \sqrt{\frac{\chi ^2 + N_{\rm eff}\overline{r}^2}{N_{\rm obs} + N_{\rm eff}}}. \end{equation} \tag{ 14 }$$

We artificially set N_eff = 4, so on nights with fewer than four observations, the MAP value of r_{σ, w} will be weighted most toward what the rest of the season says. On nights with more than four observations, the data from the night itself will more strongly drive the MAP value.

We use a modified periodogram (Irwin et al. 2011a) as part of the season-wide fit to identify the best value of $\overline{v_P}$, the period of the harmonic terms in V(t). We fix v_P to this value in all later analysis. While this effectively places an infinitely tight prior on this parameter, the degeneracy between it and the other variability parameters, especially on the timescale of a single night of data, means that its uncertainty is usually accounted for by those terms.

The shape of $P(\mathbb {M})$ offers a big advantage to our goal of estimating the marginalized transit depth probability distribution. Accounting for all the terms in $P(\mathbb {M|D})$ (Equations (7), (9), (11), and (13)), we find that the posterior $P(\mathbb {M}|\mathbb {D})$ can indeed be maximized and marginalized analytically. For fixed p_E and p_T, the system of equations

$$\begin{equation} \frac{\partial }{\partial c_j}\ln P(\mathbb {M}|\mathbb {D}, p_E, p_T) = \frac{\partial }{\partial c_j}\left(\frac{\chi ^2}{r_{\sigma ,w}^2} +\Phi ^2\right) = 0 \end{equation} \tag{ 15 }$$

can be solved exactly for the MAP vector of values $\mathbf {c_{\rm MAP}}$ using only simple matrix operations. The procedure is directly analogous to the problem of weighted linear least-squares fitting; see Sivia & Skilling (2006, chap. 8) for details of this solution. While not strictly necessary because the priors prevent unconstrained degeneracies in the solution, we use singular value decomposition (SVD) to avoid catastrophic errors in the matrix inversions (Press 2002).

The value of r_{σ, w} sets the relative weighting between the likelihood and the prior. Thus it is important to estimate r_{σ, w} accurately. We solve for it by iterating between Equations (14) and (15); the solution typically converges to the MAP value within only a couple of iterations. We forego marginalizing over r_{σ, w}, instead fixing it to its MAP value. Solving for r_{σ, w} independently on each night is a better approximation than blindly assuming a global value.

Importantly, the matrix solution to this problem gives not only the MAP values, it also gives the covariance matrix of the parameters in the fit, which is an exact representation of the shape of $P(\mathbb {M}|\mathbb {D}, p_E, p_T)$, which is a multidimensional Gaussian. The diagonal elements of this covariance matrix give the uncertainty in each parameter marginalized over all the other linear parameters. Because we have constructed our model in such a way that the parameters that most strongly influence estimates of the transit depth D are linear, we can use this analytical solution as a robust estimator of the shape of the MarPLE. It gives us both the MAP transit depth $\overline{D}$ and the Gaussian width of the distribution σ_MarPLE.

As the priors we use to regularize our model fits are themselves derived from MEarth data, one might object that the division between the likelihood and the prior is set somewhat arbitrarily. We include data only from a single night in the likelihood and group all the information from the rest of the nights into the prior. Indeed, we could have instead organized the entire season of data into the likelihood and left the priors uninformative. The division is arbitrary, but useful.

The advantages of treating nights other than that on which a candidate transit falls as external to likelihood are two-fold. First, it is more computationally efficient: instead of recalculating the likelihood of an entire season's data when investigating individual events, we only need to calculate the likelihood over the relevant night's data points. The information provided by the entire season changes little from candidate transit to candidate transit; thus it is best to store that information as a pre-computed prior.

Second, this organization scheme allows the flexibility for individual nights to behave differently. For example, consider the pixel position E_{x, i}(t) and E_{y, i}(t) terms in the systematics model, which capture the influence of stars wandering across the detector. As the detector flat-field can change from night to night, it would be foolish to try to fit an entire season's light curve with one set of coefficients for E_{x, i}(t) and E_{y, i}(t); allowing those coefficients to vary from night to night, within a tightly constrained prior, is a more useful approach. Furthermore, dividing the weight of the likelihood and priors as we do provides a helpful degree of outlier resistance, by not forcing the model on any one night to account for strange behavior on one weird night from months before.

We calculate $\overline{D}$ and σ_MarPLE on a grid of single transit epochs p_E and durations p_T. We construct this grid for all nights with usable MEarth data. The epochs in this grid are evenly spaced by Δp_E = 10 minutes, thus subsampling the typical MEarth observational cadence. The durations are evenly spaced from 0.02 to 0.1 days, spanning the likely durations for the orbital periods to which MEarth has substantial sensitivity. We extend the grid of p_E before the first and after last observation of each night by half the maximum transit duration, thus probing partial transits.

We demonstrate this calculation graphically in Figure 4. First, we show one night of a typical MEarth light curve. To give a sense of a baseline systematics and variability model, we show it the ±1σ span of model light curves arising from a fit that contains no eclipses. Next, we show a visualization of the MarPLE, the probability distribution of hypothetical eclipse depths P(D|p_E, p_T). For any chosen value of eclipse epoch and duration, the MarPLE is Gaussian shaped; error bars in Figure 4 show its central ±1σ width over the entire grid of p_E and p_T. The width of σ_MarPLE can be seen to decrease for longer durations p_T, as more data points are included in each transit window. For epochs and durations with no in-transit points, the transit depth is unconstrained and σ_MarPLE → ∞. The transit-like dip in photometry at the end of the night does not register as significant anywhere in the MarPLE. The dip can be explained by MEarth's PWV systematic (see the first panel of Figure 1, corresponding to the same night).

Zoom In Zoom Out Reset image size

At the bottom of Figure 4, we also include a subset of the variability and systematics parameters, the "nuisance parameters" over which we marginalize. We show error bars representing the Gaussian widths of both the prior $P(\mathbb {M})$, established from the entire season of data, and the results of a fit to this one night of data, $P(\mathbb {M}|\mathbb {D})$. For v_night, the nightly out-of-transit baseline level parameter, one night's data are more influential than the relatively weak prior, so the fit is notably offset from and tighter than the prior. In contrast, for the remaining nuisance parameters, the influence of one night's data is very weak, so the fit essentially reverts to the input priors.

In this example, only 10 data points are contributing to the likelihood. As the model contains almost as many parameters, one might be concerned that we are "overfitting" the data. The bottom of Figure 4 provides an initial step to allay this concern, emphasizing that except for the transit depth, each parameter in the fit has its associated prior that provides its own independent constraint on the parameter. In a pseudo-least-squares formalism, the presence of these informative priors act as (pseudo)data points, ensuring there are always more "data" than parameters.

Furthermore and perhaps more importantly, we could indeed be in severe danger of overfitting if we were interested in the exact values of the cleaned residuals from some single estimate of a best-fit systematics and variability model, but we are safe because we care instead about the marginalized probability of only one particular parameter (the transit depth). Marginalization ignores irrelevant information, so we can include an arbitrary number of nuisance parameters in the fit (see Hogg et al. 2010 for discussion). If (and only if) the inferred transit depth at any particular p_E and p_T happens to be strongly covariant with one of these nuisance parameters, then σ_MarPLE will include a contribution from that parameter. Without the priors, degeneracies could potentially inflate σ_MarPLE to ∞, but with the informative season-long priors the nuisance parameters can only vary within the range shown in Figure 4, limiting the degree to which they can in turn contribute to σ_MarPLE. In the extreme example, if we had a single data point on a night, the cleaned residuals might easily be identically zero (i.e., "overfit"), but the MarPLE would accurately express what the night told us about the presence or absence of transits.

Estimating σ_MarPLE across the whole grid of p_E for an entire season can be performed very quickly. Each grid point requires only several SVD's of a matrix whose dimension is the sum of the number of data points within the night and the number of linear parameters being fit. For a MEarth light curve containing 1000 points and spanning 100 days, the whole grid of calculations requires several seconds on a typical desktop workstation.

The likelihood in Equation (9) assumed that adjacent light curve data points were statistically independent. If our method fails to completely correct for systematics or stellar variability, this assumption will be violated. Time-correlated noise slows the $\sqrt{N}$ improvement that would be gained by obtaining N independent Gaussian measurements. So, if we were to ignore the temporal correlations between data points, we could substantially bias our estimates of σ_MarPLE.

Specifically, correlated noise would cause us to overestimate the significance of transits that spanned multiple data points. We demonstrate this phenomenon in Figure 5, which shows the MarPLE results for two simulated light curves (generated from the real-time stamps of a typical MEarth target) with different levels of correlated noise. One light curve consists of pure white Gaussian noise. The other consists of the white light curve averaged with a smoothed version itself, scaled so both light curves have an identical rms, roughly approximating a finite red noise contribution. We then inject common mode and meridian flip trends are injected into both light curves. The toy-model simulation of time-correlated noise is very coarse but is only meant to serve an illustrative purpose.

Zoom In Zoom Out Reset image size

Because each estimate of $\overline{D}$ is drawn from a Gaussian distribution with a width σ_MarPLE, the quantity $\overline{D}/\sigma _{\rm MarPLE}$ should ideally be Gaussian distributed around 0 with a variance of 1, except when real transits are present. For the light curve with pure white noise, this is true for all transit durations in Figure 5 (see the histograms at right). For the light curve with significant correlated noise, we underestimate σ_MarPLE and the distribution of $\overline{D}/\sigma _{\rm MarPLE}$ appears broadened for some durations. The effect is most pronounced at longer durations, where more data points fall within each transit. For the shorter durations, typically only one or two light curve points fall within a transit so the red noise does not substantially affect our estimate of σ_MarPLE. This general behavior, of overestimating the significance of longer duration transits, is common among MEarth targets whose light curves show features that are poorly matched by the input model.

To account for the problem, we posit that each light curve has some additional red noise source that can be expressed as a fixed fraction of the white noise, defining r_{σ, r} as the ratio of red noise to white noise in a light curve. With this parameterization, the transit depth uncertainty associated with a transit that contains N_tra data points becomes

$$\begin{equation} \sigma _{\rm MarPLE} = \sigma _{{\rm MarPLE}, w} \times \sqrt{1 + N_{\rm tra} r_{\sigma , r}^2}, \end{equation} \tag{ 16 }$$

where σ_{MarPLE, w} is the estimate of σ_MarPLE that accounted only for white noise. To determine its optimum value, we scale r_{σ, r} until the distribution of $\overline{D}/\sigma _{\rm MarPLE}$ has a MAD of 1/1.48 (i.e., the distribution has a Gaussian width of unity). This correction is similar to the $\mathcal {V}(n)$ formalism described by Pont et al. (2006). Henceforth, when we use the term σ_MarPLE, we are referring to its red-noise-corrected value.

A more ideal solution would account for time-correlated noise directly in the likelihood (Equation (9)), but doing so would substantially decrease MISS MarPLE's computational efficiency. As such, we settle on Equation (16) as a useful ad hoc solution. Figure 6 shows the amplitude of r_{σ, r} for all stars in the MEarth survey, indicating that most stars have low red noise contributions, after accounting for our stellar variability and systematics.

Zoom In Zoom Out Reset image size

We present a couple of illustrative simulations to give a sense of how MISS MarPLE works. In each case, we use fake planets with three observed transits and periods near 10 days. While such long-periods would realistically offer this many transits only rarely, we use these hand-picked examples as a convenient way to show both what individual transits of 10 day period planets look like, and what phasing these transits into periodic candidates looks like. For simplicity's sake, we left MEarth's real-time trigger out of these simulations, showing what individual transits and phased candidates look like in low-cadence data. In reality, most of the injected transits above 3σ would have been detected by the real-time trigger, and their egresses' populated with additional high-cadence observations.

Figure 7 shows one example, a 2.5 R_⊕ radius planet with a P = 9.89 day period and b = 0.1 impact parameter injected into the raw MEarth light curve of a 0.21 R_☉ star. In this case, the 8.0σ injected S/N of the transit is well recovered by MISS MarPLE at 9.2σ, as is the inferred planet radius. Three transits fell during times of MEarth observations; they are marked in the plot of D/σ_MarPLE, the eclipse S/N. This star exhibits 0% residual red noise and the transits all fall within well-sampled nights; it is thanks to these favorable conditions that the injected and recovered S/Ns are so similar.

Zoom In Zoom Out Reset image size

For contrast, Figure 8 shows another example with a different star but broadly similar planetary parameters. Here, the recovered signal's 5.1σ significance is considerably lower than its injected 9.5σ strength. One reason for the difference is that the timescale of the intrinsic stellar variability of this star is short enough that the inferred transit depths are substantially correlated with it, thus making a larger contribution to σ_MarPLE. Additionally, our model does not completely remove all the structured features in this light curve so it exhibits a large red noise fraction (r_{σ, r} = 0.5), further suppressing the detection significance.

Zoom In Zoom Out Reset image size

We also show in Figures 7 and 8 the photometry from MEarth, before and after using the MAP values of our model parameters to subtract off systematics and stellar variability from the light curves. To emphasize that the result of MISS MarPLE is not simply one best-fit model of the systematics and variability, but rather an inferred probability distribution, we plot the swaths of light curve space that are spanned at ±1σ by this distribution of models. We note that the probability distribution $P(\mathbb {M} | p_P, p_{E0}, p_T)$ is conditional on transit period, epoch, and duration, so when we visualize the models with the light curves, we have fixed these parameters to their best values (as found in the grid search in Section 4.3). Because the transit search is entangled with the cleaning process, the models and appearance of the MAP-cleaned light curve would be different for different choices of p_P, p_E0, and p_T.

By itself, a search with BLS will give the significance of a candidate transit that is conditional on the assumption that the out-of-transit baseline flux is constant and that its noise properties are globally known. If preceded by a light curve cleaning step, the transit significance is also conditional on the assumption that the aspects of the cleaning are correct. An important question is how much the marginalized significance of candidate transits found with MISS MarPLE differs from this conditional significance. Generally, the answer to this question will depend on the time sampling of the observations; for a very well-sampled and well-behaved light curve, the BLS and MarPLE results should converge to the same answer. But for the case of the real MEarth data, with its large gaps and fickle systematics, we approach this question with simulations.

Figure 9 shows the results of a head-to-head comparison of the significance with which MISS MarPLE views phased (multiple-event) candidates with the significance that would go into a BLS calculation, based on ensemble of injected transits. A full period search was not run as part of these simulations; we calculated the detection statistics in both cases assuming the period was known. This is in line with our goal with MEarth that we wish merely to identify whether a signal of a given significance is present, not accurately determine the period of that signal from the existing data.

Zoom In Zoom Out Reset image size

For MEarth's best-behaved stars (with r_{σ, r}  <  0.25), the marginalized significance estimated by MISS MarPLE is typically 80% of that estimated by our idealized BLS. For these stars, properly accounting for all of the uncertainties in the cleaning process gets us to within 20% of the significance we could achieve in the unrealistic hypothetical that there were no uncertainties in the cleaning process. That MISS MarPLE tends to be more conservative than other methods is very important to our ultimate goal of using candidates identified by MISS MarPLE to invest limited period-finding follow-up observations. The 20% factor suppression of transit significance is comparable to the degree to which global filtering methods such as TFA suppress estimated transit depths (e.g., HATNet; see Bakos et al. 2012). However, the advantage of MISS MarPLE is more than simply knowing how much light curve cleaning suppresses transit significance on average; it is knowing what the cleaning's relative influence is on individual events and which events are more, or less, reliable. MISS MarPLE can give good events on good nights appropriately higher weight, unlike more global methods.

Figure 9 also shows that the penalty imposed by the red noise correction for those stars with r_{σ, r} > 0.25 is significant but not always debilitating. Because MEarth's cadence is so low that typically only one to two points fall within any given transit window, the influence of red noise on most transits is relatively small. However, in cases where the cadence is much higher, such as a triggered event observed in real time with MEarth, the red noise penalty could be much steeper. Also, as D_injected/σ_injected is the best we could hope to achieve for each candidate, it is an important check that very few stars show significance ratios >1.

As more nights of observations are gathered, the priors on the systematics and variability parameters associated with a particular star will tighten. As these priors tighten, the significance with which a given transit can be detected will improve. We demonstrate this phenomenon graphically in Figure 10, which shows how σ_MarPLE for single events evolves as more observations are gathered as well as the impact of this evolution on the planet detection.

Zoom In Zoom Out Reset image size

We injected transits as before but calculated the MarPLE for every individual event using only the data up to and including the event, excluding all data after the third contact. These "time-machine" simulations are an approximation to the information available to the MEarth real-time trigger system when deciding whether to gather high-cadence follow-up of a candidate transit. We show the results for stars in the best and worst halves of the MEarth sample, as judged by how their white noise rescaling factors r_{σ, w} compare to the median of the sample $\overline{r_{\sigma , w}}= 1.24$. Note that transits have a distribution of injected transit depth uncertainties (σ_injected) based on the number of points in transit and the points' relative predicted uncertainties σ_the(t).

We highlight in Figure 10 the smallest planet that could be detected at 3σ confidence in a single low-cadence event, and how this quantity evolves a function of the number of nights a star is observed before the event. Imagine that a light curve contains 99 event-less nights and one event on the 100th night; Figure 10 indicates how much the information in the event-less nights improved the reliability of the single event's detection. In each panel, we show the 25% and 75% quartiles of the distribution (spanning both multiple stars and multiple random transits). For the stars with low r_{σ, w}, initially only planets larger than 2.5–3.8 R_⊕ exhibit deep enough transits to be detectable. But as more nights of observations tighten the priors, 2.0–2.6 R_⊕ planets become detectable, approaching the injected distribution. Stars with high r_{σ, w} behave very differently, presumably because our model captures fewer of the features present in the light curves. For these stars, the minimum detectable planet sizes initially span 3.0–4.4 R_⊕ and never converge to the injected values.

We also show in Figure 10 the distribution of the ratio σ_MarPLE/σ_injected for the simulated transits. The ratio starts off well in excess of unity, but approaches it as more prior-establishing observations are gathered. The range of values it spans corresponds to transits falling at more or less opportune moments. Values of σ_MarPLE/σ_injected closer to 1 are usually associated with transits that fall in the middle of a well-behaved night. Higher values correspond to events that fall at the start of a night, events in a night with high excess scatter, or events that coincide with transit-like features in the systematics or variability models. By the end of a season, the distribution of σ_MarPLE/σ_injected for single events in Figure 10 roughly approaches that for phased candidates in Figure 9. This makes sense, as the phased S/N ratios in Figure 9 use priors established from all the nights.