{'fraction_non_alphanumeric': 0.04437467553336135, 'fraction_numerical': 0.014857481891671406, 'mean_word_length': 4.324032640168466, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Quasiperiodic Automated Transit Search Algorithm

We present a new algorithm for detecting transiting extrasolar planets in time-series photometry. The Quasiperiodic Automated Transit Search (QATS) algorithm relaxes the usual assumption of strictly periodic transits by permitting a variable, but bounded, interval between successive transits. We show that this method is capable of detecting transiting planets with significant transit timing variations (TTVs) without any loss of significance --"smearing"-- as would be incurred with traditional algorithms; however, this is at the cost of an slightly-increased stochastic background. The approximate times of transit are standard products of the QATS search. Despite the increased flexibility, we show that QATS has a run-time complexity that is comparable to traditional search codes and is comparably easy to implement. QATS is applicable to data having a nearly uninterrupted, uniform cadence and is therefore well-suited to the modern class of space-based transit searches (e.g., Kepler, CoRoT). Applications of QATS include transiting planets in dynamically active multi-planet systems and transiting planets in stellar binary systems.


INTRODUCTION
This paper describes a new algorithm for detecting the transits of an extrasolar planet in a time-series of stellar flux observations. Many algorithms have been developed and specialized for this exact purpose (e.g., Kovács et al. (2002), Renner et al. (2008), Defaÿ et al. (2001), Aigrain & Favata (2002); see Tingley (2003) and Moutou et al. (2005) for a comparison of many of these methods). The variety of implementation in those algorithms arose primarily from the need to optimize transit detection (in the form of increased significance) in response to the nature of the data (e.g., Grziwa et al. (2012), Jenkins et al. (2010)), or to speed or automate the detection of transiting planets in a large collection of light curve data (e.g., Weldrake & Sackett (2005)). These specializations did not, however, address the potential degradation of detection significance as a result of the non-Keplerian planetary motion (exceptions are the algorithms proposed by Ofir (2008) and Jenkins et al. (1996) to detect transiting circumbinary planets). Specifically, all of these algorithms operate on the assumption of strict periodicity of the arrival of transits.
Strict periodicity is not expected when greater than two bodies are interacting gravitationally (e.g., in many-planet systems or for planets in binary systems). In these cases, we would expect non-linearities in the transit times, commonly referred to as transit timing variations (TTVs). The detection of TTVs can suggest the presence of an unseen perturbing companion (e.g., Ballard et al. (2011), Nesvorný et al. (2012) and yield information on the masses of the bodies involved in the gravitational exchange (Agol et al. (2005), Holman & Murray (2005)). The observation of TTVs in several transiting planets in a single system has proved invaluable in confirming their planetary nature and, more importantly, has placed useful constraints on their masses (Holman et al. (2010), Lissauer et al. (2011), Cochran et al. (2011), Carter et al. (2012). Systems of this type are most likely to be detected by space-based transit surveys (e.g., Kepler, Borucki et al. (2010) and CoRoT, Baglin et al. (2006)) given their nearly continuous observing modes. García-Melendo & López-Morales (2011) have shown that the failure of traditional transit algorithms to account for TTVs would result in reduction of detection significance (transit 'smearing'). They showed that the reduction in significance is non-negligible with 'typical' system architectures comprised of many interacting planets in a single system. In particular, for a transit-timing variation with an amplitude σ T TV that is longer than the transit duration and having a period shorter than the duration of the data, the signal-to-noise is reduced by a factor of T transit /σ T TV , where T transit is the transit duration. The transit-timing variation amplitude tends to grow in proportion to orbital period, P, while the transit duration grows as P 1/3 . As such, the reduction is most significant for long period planets that are perturbed by nearby companion.
The detection algorithm described in this paper (the Quasiperiodic Automated Transit Search or QATS) presents one possible solution to the problem raised by García-Melendo & López-Morales (2011) by allowing bounded variation in the intervals between transits when determining the detection significance, accounting for possible TTVs. We pose the detection problem in § 2 and derive the algorithmic solution (QATS) in § 3 by specializing and expanding the work by Kel'Manov & Jeon (2004) who tackled the more general problem of detection and pulse-shape discrimination. In § 3.4.1, we estimate the background significance in the presence of white noise as a function of the algorithm parameters. In § 4, we provide a worked example of the application of the QATS to the Kepler data for Kepler-36 (Carter et al. (2012)).

FORMULATION OF THE PROBLEM
The QATS algorithm is capable of detecting the transits of an exoplanet in a series of observations of the relative flux of its host star.
We require -in order to simplify our discussion and to facilitate a fast algorithm -that the N observations of the normalized and median removed flux, F(t n ) for 0 < n < N − 1, follow a uniform cadence; i.e., t n − t 0 ∝ n. Given this, we may uniquely specify the time of any measurement by referring to its sequential cadence number n. We use this index exclusively in what follows as opposed to referring to the absolute time and write F(t n ) as F n .

Quasiperiodic transit light curve model
We assume that the data are completely explained by a transit light curve model, L(t n ) ≡ L n , that is contaminated by additive, independent, identically distributed Gaussian noise of width σ, such that F n = L n + ǫ with ǫ ∼ N (0, σ). Refer to Figure 1 for a qualitative representation of the model, described next.
During transit, the light curve decreases for a short amount of time when the planet occludes the star. We make the usual approximation (e.g., Kovács et al. (2002)) that the light curve profile of a single transit event is boxed-shaped with depth δ and duration q. We assume that q is an exact integral number of cadences and that the flux deficit starts precisely on a cadence.
Suppose that M transits of fixed shape (fixed q and δ) occur in the sequence {L n }. Let the collection of indices at the start of each of the M transits be (1) We may compactly define the model light curve sequence {L n } as the sum of flux deficits from all M transits as where we have defined the "root" transit deficit profile as For a strictly periodic transiting planet, we would expect the interval, ∆ m , between successive transits to obey (for m > 0) for a constant periodP. We have distinguished this period, measured as an integral number of cadences, from the absolute period, P, which is measured in units of time and may be fractional. We wish to relax the constancy of this interval to allow for possible transit timing variations. Fully unspecified intervals need not be considered; we expect, based on physical considerations, that these a priori unknown variations (in interval) are small relative to the mean interval. The smallness of this variation depends on the specifics of the physical system to be detected. For example, the amplitude of this variability is ∼ 1% of the mean orbital period for known cases of interacting planets in multi-transiting systems (e.g., Holman et al. (2010), Lissauer et al. (2011), Carter et al. (2012). For transiting circumbinary planets, the amplitude is a much larger fraction: one would expect timing variations comparable to the orbital period of the stellar binary; this may amount to a variation ∼ 10% of the mean orbital period of the planet (e.g., Doyle et al. (2011)). In either case, we may specify a variable but nearly periodic, or 'quasiperiodic' bounding, between successive transits: (Quasiperiodic Transit) where ∆ min and ∆ max are the minimum and maximum intervals considered. Following the above, it is sometimes convenient to parameterize the difference in the maximum and minimum interval to be some small fraction, f , of the minimum interval: ∆ max − ∆ min = f ∆ min .
In addition to the quasiperiodic bounding, we require that the transits be non-overlapping and that the first transit (starting at cadence n 1 ) and the last transit (starting at cadence n M ) are fully contained within the sequence. These conditions are equivalent to Subject to the constraints above, the permitted number of transits in L n lies between a minimum (M min > 1) and maximum (M max ) where it is simple to show that 2.2. Maximum-likelihood objective function The QATS algorithm is a maximum-likelihood method. In other words, QATS determines the set of starting cadences (η M,best ), the transit duration (q best ) and transit depth (δ best ) that best agrees with the data (has the maximum likelihood), subject to the quasiperiodic condition.
We have derived an expression, more suited for numerical computation than using the full likelihood, that is also maximized for the optimal parameters. This objective function, S(η M , q), is independent of σ and the transit depth when the noise is stationary. The details of this derivation have been provided in the Appendix. The result is The most likely transit depth (also derived in the Appendix) is S(η M , q) has a simple physical interpretation when maximized: where (S/N) total is the total transit signal-to-noise ratio.

THE QATS ALGORITHM
The previous section introduced an expression, S(η M , q) (Eqn. 9), that is maximized when the likelihood (Eqn. A3) is greatest. The task remains to describe the algorithm that actually maximizes S(η M , q) with respect to the transit duration, q, and the set of starting cadences, η M for a selected ∆ min and ∆ max (or ∆ min and fraction f ).
Plausible transit durations may be restricted to a reduced range of cadences for a given period subject to simple physical expectations (in particular, duration ∝ P 1/3 ). In any case, the allowed range of durations may be searched, marginalizing with respect to η M for each duration, to find the maximum likelihood and best-fitting duration. Similarly, the possible number of transits in an observation {F n } is enumerated from M min and M max , according to Eqns. (7,8), and may be searched explicitly.

Dynamic programming solution for maximization over η M by Kel'Manov & Jeon (2004)
The only non-trivial marginalization of S(η M , q) -over η M for a fixed q and M -is found using a specialization of the algorithm by Kel'Manov & Jeon (2004). We review this algorithm in what follows.
First, we further distill the objective function: for a fixed M and q, S(η M , q) is maximized whenS(η M , q) = M m=1 q−1 j=0 −F j+nm is maximized. Note thatS(η M , q) may be written asS where D n can be interpreted as the discrete convolution of the data with a box of unit depth and duration q (a box "matchedfilter"). At a transit event, D n has a symmetric triangular profile of width q, peaking at the start of the transit in the absence of noise (see Figure 2a)). Written this way, we see that the maximum ofS(η M , q) is the maximum of a sum of M numbers, drawn from D n with the selected indices subject to the quasiperiodic condition. The additive nature of this objective function along with the quasiperiodic constraints admits a 'dynamic programming' solution for the maximization, to be described shortly. Dynamic programming typically describes a class of algorithms in which calculations made at prior steps in an iteration are retained in memory for use in the current step 1 . Second, we enumerate all possible transit start times. Viz., the first transit in the sequence F n , beginning at index n 1 , must occur before the maximum interval has passed (according to Eq. 5). We consider, in turn, the possible starting indices n 1 in this set of possible indices. We refer to this set as ω 1 . For example, with q = 2, ∆ max = 4, the maximum index the first transit can begin at is n 1,max = ∆ max − q = 2 such that ω 1 = {0, 1, 2} and n 1 ∈ ω 1 . The second transit, starting at index n 2 ∈ ω 2 , has to occur between ∆ min and ∆ max after n 1 ∈ ω 1 and must also permit the inclusion of M − 2 more transits, all subject to the quasi-periodic condition, in the N observations. The location of the third, fourth, etc. transit must follow analogously. Summarily, each transit m of M total transits is restricted to a set of possible starting indices, ω m , subject to our constraints. Kel'Manov & Jeon (2004) showed that the possible sets of transit start indices may be exactly specified given ∆ min , ∆ max , N, q and M as It is also useful to consider the subset of ω m−1 that is conditioned on the knowledge of the starting index of the next transit. In particular, if the mth transit is known to occur at index n m = n, then subject to the quasiperiodic condition, the preceding transit (the (m − 1)th transit) must have started at index n m−1 such that n − ∆ max ≤ n m−1 ≤ n − ∆ min . We refer to the intersection of the set given by this inequality and the set ω m−1 of allowed transit start times of the (m − 1)th transit as γ m−1 (n). Again, Kel'Manov & Jeon (2004) provide the exact expression of this set in closed form: With the above formalism, we may detail the algorithm by Kel' Manov & Jeon (2004). We calculate the maximum of S(η M , q) by considering successive transits (restricted to the indices specified by the ω m ) starting with the first. At each transit number m and candidate transit start index n ∈ ω m , we calculate an intermediate quantityS mn ≡ maxS(η m |n m = n; q) where η m = {n 1 , ..., n m } is the subset of the first m transits in η M which is the maximum value ofS(η m , q) for this circumstance, i.e. excluding the contribution from proceeding transits but including the contribution from preceding transits. Given this definition and Eq. 14, we may recursively construct S mn for all transits m ∈ {1, ...M}: where we have utilized the set γ m−1 (n), as defined above, in the first line. It is then simple to show that the desired maximum of the objective function for fixed q and M is given by We calculateS mn following the dynamic programming method by minimally retaining in memoryS m−1,n for all n ∈ ω m−1 so as to calculate the mth row ofS mn without repeating any calculation. In Fig. 2 we show a graphic interpretation of the explicit execution of this recurrence, as would be performed by a computer, for a specific example.

The most likely indices at the start of the M transits
The most likely set of transit start indices η M , corresponding to maxS(η M , q), may be determined onceS mn has been calculated by traversing this matrix in reverse 2 . In detail, the index n = n M that gives the maximum value ofS Mn for n ∈ ω M (equal to max ηMS (η M , q) according to Eq. 24) corresponds to the most likely starting index of the final (the Mth) transit. We may then recursively determine the most likely starting indices of the earlier transits conditioned on the knowledge of the location of the proceeding transit(s): where argmax i∈I a i is the index i max such that a imax is maximum for all i ∈ I and we have made use of the conditional set γ m−1 (n), defined in the preceding discussion, in the second line. Fig. 3 shows this recursion for the example presented in Fig. 2 The most likely set of times may then be determined using q best and M best in Eqns. 15, 22, 23 applied to Eqns. 25 and 26.

Computational considerations
The QATS algorithm may be implemented as listed in Algorithm 1, however, there are some practical considerations.
The termS mn is most simply realized in code as a N × M matrix. This matrix can be memory intensive. For example, a typical Kepler light curve, observed at the nominal ∼ 30 minute cadence and collected through observing quarter 9 (750 days), contains N ≈36,000 cadences. A transiting planet with orbital period P admits M ≈ 750 days/P transits in that light curve. With 8 bytes allocated per matrix entry, we can anticipate requiring ≈ 200 MB · days/P of memory to storeS mn . Short period orbiting planets (P < 1 day) may require disk access, slowing the search dramatically. Additionally, when searching over multiple periods (see below), it is advantageous to conserve shared memory for the purpose of parallelization. It is therefore recommended during the calculation ofS best that only a 2 × M submatrix ofS mn be retained containing the (m − 1)th row and the mth row under calculation. The determination of the most likely times still requires the full matrix, however, these values are not needed for the calculation of the objective function maximum and can be calculated in a posterior operation.
The convolved data, D n , as defined in Eqn. 15, may be modified to support 'inline' data detrending with more complicated matched-filters. For example, the filter may include a linear correction outside of transit. We also note that the D n may be pre-computed and stored in memory for a number of common durations prior to the first loop in Algorithm 1.
It is likely that some cadences will be missing from the observation sequence (for example, those lost in a Kepler light curve during data downlink between quarters). Ignoring these cadences violates the preconditions of the QATS algorithm and may lead to unexpected results. We recommend filling these missing cadences with zero values. Note in this case that transits occurring during these gaps will provide no additional significance to a detection and the most likely times of these missing transits will be meaningless.

The QATS 'spectrum'
The QATS algorithm, detailed in Algorithm 1, determines the most likely duration and times for a transiting planet with 'period' specified by the pair of interval bounds, ∆ min and FIG. 2.-Here we show a graphical representation of the step-by-step execution of the QATS algorithm for a noiseless example. a) The 'data' are shown in the top time-series, containing three transits of duration q = 2 and of quasperiodic interval. The maximum and minimum intervals (∆ min and ∆max) are shown graphically. The box-filter has been applied to the data in the middle time-series (Dn, given in Eqn. 15). The grayscale bar represents Dn as well; each cell represents a single cadence and is color coded (increasing from white to black) according to the value of Dn. b) This panel shows the 'initial' state of the yet-to-be-computed matrixSnm, described in § 3.1. Each row represents an individual transit m (of three total). Each row currently has a copy of Dn; only those cadences without hashing in their respective cells are relevant during the computation. The allowed set of transit start times ωm (Eqn. 16) is indicated by a label and a range for each row m. c) There is nothing to be done for the first row. Here, and for c) through g) we sum the value of Dn and the maximum of the significance in the first transit over the possible set of transit start times for the first transit assuming the second transit starts at the indicated location (n = 5, 6, 7, 8 or 9). These sets are shown explicitly in the first row by the braced range in red (and labelled by γ m−1 (n) as defined in Eqn. 19). The maxima of those sets have been highlighted with a red boundary. h) -l) We repeat the procedure performed on the second row for the third transit iterating over the possible starts of the third transit (n = 10, 11, 12, 13 or 14). l) At this point, theSmn have been calculated. The maximum significance is then determined as the maximum of theSmn with m = 3; that cadence has been highlighted with a blue boundary. ∆ max . In practice, one does not know the period of the transiting planet (nor the variation on that period between transits) and multiple periods need be compared for significance. Traditional search algorithms establish a grid on period (typically being uniform in orbital frequency) and perform the fixed period marginalization over duration and phase at each grid point. There is no single period in the QATS algorithm, however, a number of strategies can be suggested to perform the analogous task. We suggest a couple below.
The first strategy is a uniform grid over minimum interval, ∆ min , with a fixed difference between this minimum interval and the maximum interval (such that ∆ max − ∆ min = ∆∆). In this case, a natural grid spacing is ∆∆ > 0. In this strategy, the allowed variation in the transit interval (or period) is the same for all trial 'periods' and may not reflect the underlying physical prior.
In contrast, the second strategy allows the difference between minimum and maximum interval to be some fixed fraction, f , of the minimum interval (some fraction of the 'period') which may better reflect our theoretical expectations (Agol et al. 2005) when the sought after transiting body is being perturbed by another body. In this case, we grid ∆ min such that the ith grid point is The fraction f may be tuned for a particular search or may be searched over a small grid as well; any transit-timing variations with a fractional amplitude of f /2 will be contained within one of the ranges searched within this grid. Both strategies can have coverings that are 'complete' and non-overlapping over any interval in ∆ min -our grid resolution does not affect our ability to detect a given period for a strictly periodic transiting planet, excluding signal confusion due to a significant stochastic background signal (see § 3.4.1).
The detection spectrum (i.e., max S(η M , q) or the maximum total transit signal-to-noise as a function of ∆ min ) has characteristics similar to traditional transit search spectra that are based on a box matched-filter. In particular, the QATS spec- trum is roughly proportional to (mn) −1/2 where m/n is the rational number with smallest denominator n that lies between [∆ (i) min /P, ∆ (i) max /P] with P being the actual period of the transit. This proportionality is simply explained: suppose your trial period is ∆ trial , while the actual period is P with ∆ trial /P = m/n where m and n are mutually prime integers. Then every mth transit will be detected, so the signal will be reduced by a factor of 1/m, while the noise will be changed by a factor of n/m: the signal-to-noise decreases by a factor of 1/ √ mn. This signal-to-noise is maximized when mn is minimized, giving our dependence. The numbers n and m may be efficiently calculated with use of the 'Farey' sequences by finding the lowest order Farey sequence which has a fraction, m/n, contained within the trial period search range (in units of the correct orbital period); this guarantees that m/n is the smallest mutually prime ratio and thus has the highest signal-to-noise relative to the peak. For trial periods above the correct period, the Farey sequence search should take place between [P/∆ (i) min , P/∆ (i) max ]. Figure 4 shows a plot of 1/ √ mn without the presence of noise for an interval width of f = 0.001. In practice, we find that this pattern matches the QATS spectrum of detected planet candidates well, while other forms of light curve variability show a different shaped QATS spectra. We note that the robustness of a detection of a transiting planet (or at least a periodic box signal) may be validated by comparing the measured spectrum to this prediction for a specified candidate period. See § 4 for example QATS spectra computed from noisy data.

Stochastic background
The QATS algorithm permits the detection of transiting planets with large variations in transit interval without suppressing the significance of the transit; however, the exchangeable cost of this freedom is increased backgrounds that may 'swamp' a potential detection signal. The stochastic background -i.e., the conspiracy of noise that yields a smooth level of non-zero significance beneath the sharp transit detections -grows as the difference between ∆ max and ∆ min widens. It is difficult to obtain an exact closed form expres- sion for this background; however, in what follows we derive an upper bound with some simple assumptions and, further, propose formulae that closely approximate the true behavior of the background.
We assume a white-noise background, F n ∼ N (0, σ 2 ); i.e., the data are independent, identically distributed Gaussian variables with zero mean and variance σ 2 . We consider the distribution of the maximum of the total transit signal-to-noise of this background where, according to Eqns. 13, 9 and 14, max(S/N) total,BG = max Referencing Eqn. 15, we see that the D n qσ 2 −1/2 ≡ D ′ n are Gaussian distributed with zero mean and unit variance and are correlated whenever |n m1 − n m2 | < q. The requirement that the transits be non-overlapping (Eqn. 4) implies that the {D ′ nm } are independent. As such, it is simple to see that the sum Σ ηM ≡ M m=1 D ′ nm M −1/2 is distributed identically as the D ′ n ; these variables are also correlated as discussed below.
We may then regard the maximization of Σ ηM over η M as being approximately equal to the maximum, X, of some appropriate number of independent draws K from a normal distribution with zero mean and unit variance. The number K represents the effective number of independent transit configurations, yielding independent sums Σ ηM . It can be shown (pg. 374, Cramer (1946)) that the mean of the extremum variable X is asymptotically approximated for large K as When ∆ max = ∆ min (i.e., periodic transits), the number of draws K is equal to the number of independent choices for the start time of the first transit as there is no additional freedom in the choices of the remaining M − 1 transits. There are K = ∆ max − q + 1 available cadences for the start of the first transit; however, for q = 1, these choices are not independent (with correlation length ∼ q). We can account for this over-counting by dividing K by an effective correlation length l = q/α for some α ≥ 1. Then, the number of independent sums is K periodic ≈ α∆ max /q. Referring to Eqn. 31, for ∆ max large, we have that max(S/N) total,BG ≈ 2 log α∆ max /q . (Periodic Search) This analytic theory compares well with simulation (see Fig. 5) so long as α = α(q); α = {1.0, 2.2, 3.0} gives excellent agreement with simulated periodic QATS backgrounds when q = {1, 5, 14}, respectively. Obviously, this dependence also applies to singleton transits for any choice of ∆ max > ∆ min .
We may similarly consider the total number of ways to admit M transits in F n for ∆ max > ∆ min . As before, the first transit may start at ∆ max − q + 1 positions. Following the choice of the mth transit, the following (m + 1)th transit must start at a choice of ∆ max − ∆ min + 1 positions, according to the quasiperiodic condition. Applying the rule of products, and accounting for the correlation length at each transit, the total possible transit configurations is therefore However, it must be the case that max(S/N) total,BG ≤ X KTotal as the K total choices are not independent choices in this problem. This is because the contribution to the maximized sum, Σ ηM , is greater for m small than that for m large. We may attempt to account for this variable contribution by weighting the above product by a simple function of transit number m: where 1/r is a constant shaping parameter. We then expect the background signal in the QATS search to have a dependence according to max(S/N) total,BG = X K (Eqn. 31) with K = K QATS where, to simplify the computation, in the typical limit of K QATS ≫ 1, In Fig. 5, we plot the analytic prediction for the background as a function of ∆ max /q for a selection of ∆ max − ∆ min , in units of duration q, along with simulated backgrounds of the same parameters using the QATS algorithm. In this figure, we have assumed M ≈ (N/q)/(∆ max /q) with N/q = 2500 which is reasonable in a QATS search for transits of 7 hour duration in a long cadence Kepler light curve spanning 750 days (N = 35000, q = 14). We find that r is dependent on ∆ max − ∆ min with r ≈ 11.76, 7.14 and 6.25 giving the closest agreement between simulation and theory for ∆ max − ∆ min = 1q, 2q and 3q, respectively. These choices of r (and those for α = α(q) as provided above) have been used in generating the theory curves in Fig. 5.
The background level may be compared to the detection threshold of 7.1 set by the Kepler mission (Jenkins et al. (2010)), for example. The background exceeds this threshold for short periods, ∆ max /q 40, when ∆ max − ∆ min = 1q. For the Kepler-specific example, this corresponds to absolute periods of 10 days with quasiperiodic variability f 2.5% in 750 days of long cadence data.
The above formulae are useful to gauge the growth of the background QATS signal in presence of white noise, however, the actual background in a specific search is best represented via simulation if possible.
The starting indices of the most likely transit signal in the presence of background alone are best described by a one dimensional constrained random walk whose intervals are approximately uniform between the minimum and maximum intervals. This characteristic is specific to random noise and may be diagnostic of such; we would generically expect physical timing variations to have a non-stochastic form.

Run-time analysis
The success of the QATS algorithm is contingent on being both precise and efficient in computational time (and secondarily in memory-usage, see § 3.3). In the case of the latter, the QATS algorithm performs its most complex computational task and the only task that is distinct from fixed-period searches, that of computingS mn (described in § 3.1), in 'polynomial' time-complexity. In other words, the number of simple computational evaluations on a computer scales as some (small) power of the parameters of the problem. In particular, for fixed q and M, it can be shown (Kel'Manov & Jeon (2004)) that the computation ofS mn and the maximization of the final row (see § 3.1) has time-complexity The time to execute a search scales only linearly with the (QATS-specific) interval difference ∆ max − ∆ min . For the most typical case of ∆ max − ∆ min small (see § 2.1), the computational penalty of applying the QATS algorithm as opposed to a traditional fixed-period search is negligible.

EXAMPLE APPLICATION: KEPLER-36
In this section, we present a complete application of the QATS algorithm on Kepler light curve data for the confirmed planetary system Kepler-36 (Carter et al. (2012), planet b also detected by Ofir & Dreizler (2012) as KOI-277.02).
Kepler-36 is a planetary system with two confirmed transiting planets with orbital periods of approximately 13.84 days and 16.23 days. The proximity of these planets' orbits permits significant gravitational interactions at their closest approach which modify their orbits sufficiently to yield significant (many hour) piecewise non-linearities in both their times of transit. As may be expected (see, e.g., Ford et al. (2012)), the timing non-linearities are anti-correlated between the two planets.
The planets' radii differ by a factor of nearly 2.5 such that the transit of the smaller, interior planet b results in a loss of stellar flux that is only ≈ 17% that caused by the transit of the more massive planet c; while the transits of planet c are significant by eye, those due to b are not. While both planets may be detected with a fixed-period search (Ofir & Dreizler (2012) The broken lines show the background total transit signal-to-noise significance as determined from a Monte Carlo simulation using the QATS algorithm for a selection of transit durations q and ∆max − ∆ min . The solid lines plot theory curves, given by Eqn. 31 with K given as in Eqn. 34, for the parameters described in § 3.4.1.
the QATS algorithm produces detections at much higher significance (by factors ∼ 1.7) and provides transit times that can be seen immediately to be anti-correlated without the need for any further analysis.

Preparation of the Kepler data for Kepler-36
The available Kepler data for Kepler-36 span approximately 877 days of nearly uninterrupted observation at the 29.4 minute long cadence interval (equivalent to 43,053 cadences). We utilize the 'raw' aperture photometry light curve (SAP_AP_FLUX in the Kepler fits product, see Jenkins et al. (2010)). A simple moving average of 50 cadences in width (≈ 1 day) was divided through the data to remove astrophysical and systematic trends. The detrended light curve was then median subtracted and missing cadences were filled with zeroes. The resulting light curve is shown in the top panel of Fig. 6; the lower panel shows the light curve after removing the readily visible transits of planet c and replacing those transit cadences with zeroes. We assume the data have point-topoint scatter σ = 7.89 × 10 −5 according to Carter et al. (2012).

The detection of planet c
We first apply the QATS algorithm to the light curve shown in the top of Fig. 6, containing the transits of both planets b and c. We execute the QATS algorithm outlined in Algorithm 1 for each ∆ min,i in a grid where ∆ min,i = i for 0 ≤ i ≤ 43, 052 and for three different quasiperiodic constraints ∆ max,i − ∆ min,i = {0, 1, 2}. In other words, we find the maximum total transit signal-to-noise, max S(η M , q)/σ, for all minimum intervals (in integer cadence) plausible for detection in the available data. We also fix the transit duration to q = 14, close to the duration, measured in cadences, as reported by Carter et al. (2012) for planet c.
We plot in the top of Fig. 7 max S(η M , q)/σ as a function of ∆ min for each choice of quasiperiodic constraint. The profile of the QATS spectrum is dominated by the transits of planet c and is similar for all three choices (with profile explained analytically as described in § 3.4). However, the significance of the detection peak (around 794 cadences or ≈ 16.2 days) and the level of the stochastic background (see § 3.4.1) increase with increasing ∆ max,i − ∆ min,i . In particular, we find that the detection significance of the transits of planet c increases from ≈ 80 to ≈ 130 with the most dramatic improvement going from strictly periodic transits (∆ max,i − ∆ min,i = 0) to that with the freedom of single cadence between subsequent transits. This improvement in significance coincides with better estimates of the instants of the transits (as shown with the progression of 'riverplot' figures in the bottom panel of Fig. 7) resulting in less dilution or smearing of the transits by the featureless light curve baseline.

The detection of planet b
After removing the transits of planet c (having located their starts in the previous QATS application with ∆ max,i − ∆ min,i = 2), we perform a subsequent QATS search under the same conditions so as to identify the transits of planet b. We plot in the top of Fig. 8 the maximum signal-to-noise as a function of a (smaller range) of ∆ min for ∆ max,i − ∆ min,i = {0, 1, 2}. The peak in this spectrum corresponding to planet b (at ∆ min = 677) is the only dominant feature in this secondary spectrum outside of the sloped stochastic background. The background, in this case, is appreciable compared to the detection peak, again growing with increasing quasiperiodic interval. The detection peak remains clear (despite this background) growing from ≈ 10 for periodic transits to ≈ 17 for the most liberal constraint of ∆ max,i − ∆ min,i = 2. The significance of individual transits of planet b are close to 3, making their visual identification difficult -even when folded at the best-fitting linear ephemeris (see the bottom panels of Fig. 8). However, with the guidance of the most-likely transit instants in the final figure of the lower panel of Fig. 8, its presence becomes clear. This is made even more compelling given that the transit times are anti-correlated with those of planet c (compare the final figures in Figs. 7 and 8).
A subsequent execution of the QATS algorithm, having now removed the transits of planet b, rendered no additional detections exceeding a threshold of 7.1.

SUMMARY
This paper described the Quasiperiodic Automated Transit Search (QATS) algorithm, based upon a specialization of the algorithm suggested by Kel'Manov & Jeon (2004). This algorithm is summarized in pseudocode in Algorithm 1 with certain definitions provided in § 3.1.
This algorithm extends fixed-period box maximumlikelihood search methods (e.g., Kovács et al. (2002)) by allowing for a variable interval between transits that is bounded by user-defined minimum and maximum values. For this reason, the QATS algorithm is well-suited to the unencumbered detection of transiting planets with substantial transit timing variations resulting from N-body interactions in a multi-object planetary system (García-Melendo & López-Morales (2011)). Such transiting planets have proven to be invaluable in characterizing the low mass end of the planetary mass-radius curve (e.g., Lissauer et al. (2011), Carter et al. (2012).
The additional computational load of the QATS algorithm relative to that of fixed-period searches is negligible for most searched conditions. However, in order to produce such an efficient algorithm we require certain preconditions on the data including a uniform cadence without any missing observations. Such conditions are closely approximated for the Kepler and CoRoT missions and may be reached in any circumstance following some guidelines described in § 3.3.
The additional freedom of quasiperiodicity comes at the cost of increased backgrounds. However, we show, in § 3.4.1, that the growth of this background is sub-linear with increasing quasiperiodicity and would not likely exceed detection threshold limits in most scenarios.
The algorithm presented here is generic for any quasiperiodic box-like signal and its application may not be restricted to transiting exoplanets alone. Going further, altering the box filter to match a different pulse shape could extend the algorithms utility beyond box-shaped pulses; however, this is at the expense of the transit-oriented formalism presented in this paper. We direct the reader to the more generic work of Kel'Manov & Jeon (2004)  The likelihood that the model L n represents the data F n -which is assumed to be contaminated by additive, Gaussian white noise with characteristic width σ -is Top. The top panel shows the QATS 'spectrum,' maximum detection significance as a function of minimum interval (in cadences) ∆ min , for the data shown in the top panel of Fig. 6 for a selection of fixed ∆max − ∆ min . The highest significance peak corresponds to the period of planet c (again in cadences; each cadence is 29.4 minutes long, the period of planet c is ≈ 794 cadences or ≈ 16.2 days). The profile of the spectrum closely resembles that described analytically in § 3.4. The significance of the detection peak increases with increasing ∆max − ∆ min (as does the background). Bottom panels. These 'river plots,' normalized stellar flux as a function of transit number and time modulo the mean orbital period of planet c, correspond to the QATS spectra in the top panel. The cadences highlighted in red give the most likely instants of the transit at each transit number according the QATS algorithm. More freedom in the quasiperiodic constraint permits more accurate determination of the transit start times (and a subsequent boost to the detection significance).
We consider the natural logarithm of this likelihood and substitute the definition of L n from Eqn. 2, L(η, q, δ) ≡ log L(F n |L n ; η, q, δ) For fixed noise (fixed, finite σ), maximization of the likelihood reduces to the minimization of the final braced expression, which is independent of σ, in the above log-likelihood. Expanding this term, The first term in Eqn. A3 is constant with respect to the parameters, so the maximization of the likelihood is reduced to the maximization of S(η, q, δ). We may write the second term in S 2 (η, q, δ) as The lower and upper bounds in the sum over j may be replaced with 0 and q − 1, respectively, as u j is non-zero only within those limits (Eqn. 3) and Eqns. 4-6 imply that n m > 0 and N − 1 − n m > N − 1 − n M > q − 1; S 2 (η, q, δ) =