{'fraction_non_alphanumeric': 0.030226377952755904, 'fraction_numerical': 0.017549212598425196, 'mean_word_length': 4.4413560411311055, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Searching for Small Circumbinary Planets I. The STANLEY Automated Algorithm and No New Planets in Existing Systems

No circumbinary planets have been discovered smaller than 3 Earth radii, yet planets of this small size comprise over 75% of the discoveries around single stars. The observations do not prove the non-existence of small circumbinary planets, but rather they are much harder to find than around single stars, because their transit timing variations are much larger than the transit durations. We present STANLEY: an automated algorithm to find small circumbinary planets. It employs custom methods to detrend eclipsing binary light curves and stack shallow transits of variable duration and interval using N-body integrations. Applied to the Kepler circumbinaries, we recover all known planets, including the three planets of Kepler-47, and constrain the absence of additional planets of similar or smaller size. We also show that we could have detected<3 Earth radii planets in half of the known systems. Our work will ultimately be applied to a broad sample of eclipsing binaries to (hopefully) produce new discoveries, and derive a circumbinary size distribution which can be compared to that for single stars.


INTRODUCTION
All transiting circumbinary planets discovered to date are larger than 3R ⊕ (Martin 2018;Welsh & Orosz 2018). This is in contrast to planets around single stars, for which smaller planets are significantly more abundant than their gas giant counterparts (Fig. 1, Petigura et al. 2013;Fulton et al. 2017). On the surface this looks like a stark difference in populations, potentially providing an insight into planet formation in different environments. However, detection limitations make such a comparison is premature.
Around single stars, planetary transits occur on regular intervals. This regularity is the basis of common detection methods such as Boxed Least Squares (BLS, Kovács et al. 2002), where the light curve is phase-folded on a periodic interval to check for significance within a folded "box". Phase-folding transits becomes essential for small planets, as their shallow transit depths might be individually insignificant, but detectable as an ensemble.
Deviations from strict periodicity are called Transit Timing Variations (TTVs). For single stars these are typically on the order of seconds or minutes (Holman & Murray 2005;Agol et al. 2005), but for circumbinary planets they can be on the order of days or even weeks (Armstrong et al. 2013). The orbital motion of the inner binary is the primary component of these large TTVs, which are on order TTV ≈ (P 2 bin P p ) 1/3 /(2π). This relative motion also induces significant transit duration variations (TDVs); a given planet may have a transit of a few hours followed by a transit of over a day in length (Schneider & Chevreton 1990;Kostov et al. 2013). There are secondary, but significant, contribu-tions to the TTVs and TDVs from three-body interactions, which occur on both long secular timescales (e.g. nodal precession and apsidal advance) and shorter orbital timescales (Schneider 1994;Mardling 2013;Martin & Triaud 2014;Kostov et al. 2014).
Since the TTVs are typically much larger than the transit durations, phase-folding photometric data on a fixed period tends to wash away circumbinary transit signals rather than stacking them coherently. This invalidates automated techniques used on single stars.
As of now, the most popular and successful means of finding circumbinary planets has been to search for their transits by eye 1 . This has accounted for all of the Kepler discoveries (starting with Kepler-16, Doyle et al. 2011) and the only TESS discovery (TOI-1338/EBLM J0608-59, ). There are two main limitations with this method. First, it requires the individual transits to be discernable by eye, inhibiting the discovery of small planets. Second, it is difficult to quantify the efficiency of this method and constrain non-detections.
There have been efforts to develop automated detection algorithms specific to circumbinary planets: Jenkins et al. (1996), Ofir (2008) (CBP-BLS), Carter & Agol (2013) (QATS), Armstrong et al. (2014), Klagyivik et al. (2017) and Windemuth et al. (2019a) (QATS-EB). Some of these algorithms have demonstrated a potential sensitivity to small planets but have not yet been applied to a large sample (Jenkins et al. 1996;Windemuth et al. 2019a). On the other hand, there have been comprehensive applications to Kepler (Armstrong et al. 2014) and CoRoT (Klagyivik et al. 2017), but only with sensitivity to gas giants. There have also been advances in modelling the dynamics of circumbinary planets (Leung & Lee 2013;Mardling 2013;Georgakarakos & Eggl 2015), but they are yet to be applied in the detection of new planets. Finally, every published circumbinary planet discovery has included a photodynamical model fitted to the data, which fully encapsulates the three-body geometry and dynamics. However, such models are only to characterise and confirm the planet after the planet has been detected by eye, meaning that the final fit can be made over a narrow, computationally-tractable parameter space.
We present a new solution to the problem of automatically finding circumbinary planets. The Stanley algorithm invokes a brute-force grid search with an N-body integrator. This most closely resembles the oldest of the above works - Jenkins et al. (1996) -whereas subse-quent work developed clever approximations to the transit timing to speed up processing. The advantage of using an N-body code is that all dynamical and geometric effects on the transit signature are naturally accounted for. The trade-off is unsurprisingly a decrease in processing speed, and so we develop an optimized search grid so that the method remains tractable. Stanley also incorperates bespoke detrending methods specific to eclipsing binaries, to best illuminate shallow transits of variable duration in the presence of stellar and instrumental variability of typically greater magnitude.
In this paper we demonstrate the new Stanley algorithm and apply it to the known Kepler circumbinary planets. We both recover the true planets and to set detection limits for each system by recovering simulated planets with variable transit depths. We also search for additional planets in each system. This paper is a prelude to a second paper, in preparation, where we will apply Stanley to a broader sample of eclipsing binaries to constrain the size distribution of circumbinary planets, potentially including new discoveries.
This paper is structured as follows. We discuss the treatment and detrending of the lightcurves in Sect. 2. The search algorithm is presented in Sect. 3. We then apply this algorithm to the detection of known Kepler planets in Sect. 4 to find the known planets, search for additional planets to place detection limits. Finally in Sect. 5 we discuss applications of our work, future improvements, and provide a more in depth comparison of our method to existing circumbinary planet techniques, before concluding in Sect. 6. Lightkurve has four settings for discriminating data based on the data quality flags: "none", "normal", "hard" and "hardest". We conservatively used the "hard" setting because from inspection of some of the light curves the default "normal" setting included a noticeable number of outlier points. The "hardest" setting is not recommended according to the Lightkurve documentation.
We only downloaded long-cadence (29.4 minutes timesteps) data for all available Kepler quarters. No Figure 1. Histogram of confirmed Kepler objects under 13R⊕ for all orbital periods (blue) and Pp > 50 days (red), where the latter corresponds to the period range of discovered circumbinary planets. Black vertical lines denote the published Kepler circumbinary planet radii, with numbers denoting their Kepler number. The circumbinary planet radii are relatively large compared with planets around single stars, both with and without the 50-day period cut. Note that the typically larger radii of the > 50 day sample is largely an observational bias due to longer period planets producing less transits. Nevertheless, roughly half of the > 50 day planets around single stars are smaller than all of the circumbinary planets. Data taken from the NASA Exoplanet Database on May 4 2020, using all Kepler Objects of Interest (KOIs) with "CONFIRMED" status.
short-cadence data were used because whilst our search and detrending algorithms are not restricted to any specific observing cadence, they are not adapted to a variable cadence. The data were stitched together using default Lightkurve tools.
Like what was done for the successful by-eye searches for circumbinary planets, we avoid using the PDC (Presearch Data Conditioning) 3 and do our own de-trending.

Orbital and physical parameters
We obtain orbital parameters from the Villanova eclipsing binary catalogue (Prša et al. 2011;Kirk et al. 2016) 4 , which contains almost 3,000 entries. The data needed are the period (P bin ), the primary eclipse time (T 0,pri ), the primary and secondary scaled eclipse dura-3 Jerry Orosz and Bill Welsh,private communication. 4 http://keplerebs.villanova.edu/ tions 5 (w pri and w sec ) and the scaled separation between primary and secondary eclipses: For a circular orbit s = 0.5. From this we calculate e bin cos ω bin from the phase offset of eclipses, and e bin sin ω bin from the relative eclipse widths, e bin sin ω bin = w sec − w pri w sec + w pri .
Typically e bin cos ω bin is better constrained than e bin sin ω bin . We then calculate e bin and ω bin from these two relations. Note that Eqs. 2 and 3 are approximations that assume the binary is inclined at exactly 90 • .
We take the masses and radii for the primary and secondary stars (m A , m B , R A and R B from the catalog produced by Windemuth et al. (2019a), which contains data for over 700 of the Kepler eclipsing binaries. The Windemuth catalog is based on the Kepler photometry and GAIA stellar characterisation. They make calibrations with the binaries with masses calculated in the traditional way as a double-lined spectroscopic binary (e.g. Matson et al. 2017) and show an accuracy to ∼ 10−20%.
Note that we use the Windemuth catalog values even when running searches on the known Kepler circumbinary planets. The discovery papers for those planets undoubtedly contain more precise stellar values, as they are based on radial velocities, eclipse timing variations and circumbinary planet transits, none of which are accounted for by Windemuth et al. (2019a). However, we use the Windemuth catalogue in this paper in anticipation of our second paper, where we blindly search for new planets, most of which around binaries for which the Windemuth et al. (2019a) parameters are the best known. We will demonstrate that even with less precise stellar parameters we recover the known planets with high significance.

Removing stellar eclipses
The fraction of time the binary spends in eclipse (primary or secondary) is ∼ (R A + R B )/a bin for circular orbits and I bin = 90 • . This is significant, particularly for short-period binaries. Ideally, one could model the stellar eclipses, subtract this model, and be left with a light curve without large chunks removed. If done properly, blended planetary transits could be discovered. In practice though, the eclipses are orders of magnitude deeper than the planetary transits we are seeking. Any slight imperfections in the binary fit might be tiny with respect to the eclipse but could be comparable to planetary transits, and in our experience attempts to model the eclipses resulted in artifact "transits". We hence avoided this, and simply removed the in-eclipse data.
In each light curve we isolate the observing cadences in eclipse according to the Villanova catalog values of the primary and secondary eclipse times (T 0,pri and T 0,sec , calculated from s in Eq. 1) and durations (w pri and w sec ). We then remove these cadences. The Villanova catalog values are fitted allowing for eccentric binaries and hence potentially substantially different primary and secondary eclipse durations.

Transit duration variations -theory
We approximate the transit duration by P P P P Figure 2. Circumbinary planet transit duration on the primary star for Kepler-47b (red solid curve) and longer planet periods (blue dashed curves). This is calculated according to the simplified Eq. 4 assuming θp = π/2 at transit, b = 0 and both orbits are circular. This variation of τ with the binary orbital phase is the dominant contribution to the TDVs.
where b is the impact parameter, R A,B is the radius of the star being transited and v p and v A,B are the planet and star transverse velocities across the sky, respectively, calculated by These calculations assume perfectly edge-on binary and planet orbits with respect to the observer. We discuss this assumption in Sect. 3.1.1. Around single stars there may be small changes to τ (TDVs) in the case of multiple planets or precession due to tides or general relativity, but these are typically on the order of seconds or minutes. Around binary stars τ may vary considerably, owing to four effects. First, as seen in Eq. 5 the velocity of the star v A,B is a function of the binary phase, which could be anywhere between 0 and 2π when the planet transits. In Fig. 2 we show the variation of τ A as a function of θ A .
Second, the phase of the planet θ p can no longer be assumed to be very close to π/2 at transit. Unlike single stars, for which the planetary velocity is therefore the same at each transit according to Eq. 5, for circumbinaries there is range of possible velocities too (albeit a relatively smaller range than for the star).
The first two effects are applicable to a non-interacting pair of static Keplerians. However, in circumbinary sys-tems the gravitational potential of the binary perturbs the planetary orbit on both short, orbital timescales and on longer secular timescales. This produces a third element to the TDVs. These variations are amplified in the presence of eccentric binary and/or planet orbits.
For the fourth and final effect, we note that all of the above is applicable to planets on coplanar orbits with respect to the binary. If the planet is misaligned, even by a few 10th's of a degree, then an additional source of TDVs will arise primarily due to variations in the impact parameter, b, which may be so large that the planet actually stops transiting for long periods of time, as explored theoretically in Schneider & Chevreton 1990;Martin & Triaud 2014;Martin 2017aand observed in Kepler-47 (Orosz et al. 2012, Kepler-413 (Kostov et al. 2014) and Kepler-1661(Socia et al. 2020.
In circumbinary transit searches these large TDVs can be beneficial. In concert with the TTVs, they form a "smoking gun" signature which is very hard to mimic by false positives (estimates for the false positive rate given in . When it comes to detrending the data, however, these TDVs present a challenge.

Transit duration variations -in the context of detrending
Every light curve contains a multitude of trends, both instrumental and astrophysical, in addition to the transits we are trying to identify. The task at hand is to remove all variations except transits, such that these transits are better illuminated. This process of detrending is the application of a high-pass filter to remove the low-frequency signals caused by instrumental and stellar variations, whilst preserving the relatively sharp, short-timescale transit features. Detrending algorithms, whether it be the commonly-used Savitsky-Golay filter (Savitzky & Golay 1964), simpler mean and median filters, or the cosine (Mazeh & Faigler 2010) and Tukey's biweight (Mosteller & Tukey 1977) filters which we will be using, all act according to a "window length". In brief, the algorithms are designed to remove signals on a timescale longer than some multiple of the window length, and preserve those on a shorter timescale.
For a planet of a given orbital period transiting a single star its constant duration can be calculated, and then this would dictate the detrending window length such that transits are preserved. For a circumbinary planet, even at a fixed orbital period, the durations vary by almost an order of magnitude (e.g. Fig. 2).
A detrending window adapted for the short circumbinary transits risks removing the longer duration ones, which are the transits which would carry the most signal-to-noise weight when it comes to detecting the planet. On the other hand, if the detrending window is widened to account for transits of all durations, then we risk a light curve full of remaining, unwanted variations.
Our method in Stanley is to take each time in the light curve and calculate how long a hypothetical planet transit would be if it occurred then. This allows us to vary the detrending window accordingly.
The binary and stellar parameters are known from the Villanova and Windemuth catalogues. To predict the parameters of the as of now undetected planet we use P p = 6.1P bin , which roughly corresponds to most of the known Kepler circumbinary planets (Martin 2018) and is optimized for the detection of the shortest possible period planets 6 , which are also going to be the most detectable since they will have the highest number of transits. As demonstrated in Fig. 2, τ is not a sharp function of P p , but we accept that some of the longer transits of longer-period planets may be negatively impacted by the detrending.
The planet's starting phase is unknown, so we arbitrarily start it at θ p = 0. We set both the binary and planet inclinations to strictly be edge-on (I bin = I p = π/2). All of the known planets have a mutual inclination (∆I) less than 4 • so the assumption of coplanarity is reasonable. However, even small misalignments will create transit impact parameters which are both nonzero and variable, and so assuming ∆I = 0 guarantees that we overestimate the transit duration by at least a small amount (see Sect. 3.1.1). The binary eccentricity and apsidal alignment use the values from the Villanova catalogue, whereas the planet is set to circular 7 We calculate τ (t) with the N-body package Rebound (Rein & Liu 2012) 8 , using its IAS15 integrator (Rein & Spiegel 2015) to measure the star and planet velocities at each time step, which are converted to τ (t) using Eq. 4.

Cosine detrending
All of our light curve detrending is done using WOTAN 9 , which is a simple to use python package containing about two dozen detrending algorithms, and a companion paper comparing the methods' effi-6 Planets closer than ∼ 4P bin develop an unstable orbit on short timescales (Dvorak 1984;Holman & Wiegert 1999;Mardling & Aarseth 2001;Quarles et al. 2018). 7 We only assume a circular planet for the sake of detrending, and not in the planet detection aspect of Stanley where we search over a grid of planet eccentricities. The eccentricity of a planet will affect its transit duration, but this is a smaller effect than that of the binary phase. We also note that the eccentricities of the known circumbinary planets are all below < 0.15. 8 https://github.com/hannorein/rebound cacy . We use two of the Wotan algorithms in this paper. The first is an iterative cosine filter, which was originally developed by Mazeh & Faigler (2010). A series of sines and cosines is iteratively fitted and subtracted from the data until a convergence is found in the fit.
Our first detrending step uses sinusoids because a lot of the stellar variation seen in close binary light curves (ellipsoidal variation, Doppler beaming/boosting and reflection Faigler & Mazeh 2011), as well as typical star spot-induced rotation modulations, is inherently sinusoidal.  note that the iterative cosine method is similar to the CoFiAM (Cosine Filtering with Autocorrelation Minimization) algorithm developed by Kipping et al. (2013) and then reproduced and published by Rodenbeck et al. (2018). CoFiAM was developed for the search of exomoons, which like our target small circumbinary planets exhibit shallow transits with a complex transit timing signature (see also Martin 2017b). We favour the Mazeh & Faigler (2010) iterative cosine algorithm though since  demonstrated it was the most robust (as a function of exoplanet transits found) of the filters based on fitting splines, polynomials and sinusoids.
Our method is as follows: 1. Using the variable transit duration τ (t) calculated in Sect. 2.5, calculate the maximum transit duration, τ max .
2. Detrend the data using the cosine filter with a window length of 3 × τ max . All other parameters in the detrending are set to the Wotan defaults.
4. Compute the periodogram power corresponding to a False Alarm Probability (FAP) of 1%. This requires roughly white noise to be valid, which is the case here after the detrending.

5.
Check if there is any periodogram power above a 1% FAP at a period longer than τ max . If so, return to Step 2, but use now a window length of 2.5 × τ max . Repeat the subsequent steps, subtracting from the the window length multiplier by 0.5 each time until the periodogram loses all power above a 1% FAP for periods longer than τ max . Note that whenever we need to decrease the window length we return to the original light curve (with eclipses removed), rather than cumulatively detrending the light curve with ever-decreasing window lengths. This is such that we are always detrending an inherently sinusoidal signal.
6. If the periodogram still has power above a 1% FAP after the window length had decreased from 3 to 1 × τ max , return to Step 1 but instead of calculating τ max , calculate τ 75% . This value is defined as the transit duration for which 75% of all possible binary phases at transit should produce a shorter transit.
7. If the periodogram still has power above a 1% FAP by the time the window length reaches 1 × τ 75% then we stop the cosine detrending, lest we risk filtering out a significant number of transits.
In Fig. 3 we see the cosine filter being applied to Kepler-47. Most of the transits are preserved, but here we show here a set of three which include one long transit that gets removed by the filter. This illustrates the main challenge with detrending: filtering away nuissance signals whilst preserving the planet transits. This challenge is exacerbated by the variable and often long duration circumbinary transits. We will see later though in Fig. 11 that of the over two dozen transits of Kepler-47b almost all of them are preserved, and in Sect. 4.1 that we strongly detect the signal of the planet.

Variable-window biweight detrending
In our second stage of detrending we apply a Tukey's biweight filter (Mosteller & Tukey 1977) to the cosinedetrended light curve. This was found to be the most effective sliding filter by , and significantly better than the much more commonly used Savitzky-Golay filter (Savitzky & Golay 1964). By using a sliding filter we should be more sensitive to trends that are not characterised by simple sinusoids or polynomials, and hence might have survived the cosine detrending. Of course though, with this we must remain vigilant to avoid removing transits.
We create 48 copies of the cosine-detrended light curve and to each one we independently apply a Tukey's biweight filter. The window length used decreases linearly from the 3τ min for the 1st light curve to 0.75 × 3τ max for the 48th, where τ min,max are the bounds of the expected transit duration distribution calculated in Sect. 2.5. The factor of 3 comes from ; the biweight filter should almost completely preserve a signal with a timespan 1/3rd of the window length. The factor of 0.75 was chosen by us to make the detrending slightly more aggressive, even if some of the longest duration transits might be partially affected. All other parameters in the detrending are set to the Wotan defaults. The black curve is the cosine filter trend. The transit is highlighted by the red shading, with a variable duration (a baseline of 5 days is shown in all three plots). In all three cases the stellar activity is over an order of magnitude greater than the transit depths. The first and third transit are preserved by the cosine detrending, whereas the second transit is significantly longer and gets removed by the cosine detrending. This transit removal is an unfortunate consequence of aggressive detrending, yet since the large majority of the Kepler-47 transits are untouched we will still be able to recover the planet (see Fig. 11).
We then create a single light curve by calculating the expected transit duration at each time step (Eq 4) and then choosing the corresponding detrended light curve from the 48 above. In Fig. 4 we illustrate this procedure for the known circumbinary system Kepler-47. We plot a Lomb-Scargle periodogram of the light curve after applying different layers of the detrending. Initially there a lot of periodicity in the light curve, most prominently near the binary period at 7.4 days, and also at longer periods of ∼ 200 days, most likely corresponding to quarter by quarter variations in the Kepler spacecraft instruments. The cosine filter removes all of this long-timescale power. What remains is a small peak of stellar variability at a timescale of roughly half a day. This is then removed by the variable biweight filter, leaving us with a flat light curve plus the exoplanet transits.

Final cuts from the data
We apply a few final cuts to remove any artifacts from the detrended data, where we define an artifact as any significant deviation from a flat light curve that could not possibly be a transit. This process has a mixture of automation and human intervention. There are four things we look for.

Gaps in time
A challenge for detrending eclipsing binary data is dealing with frequent gaps in the data from where those eclipses were, in addition to the standard data gaps. The Wotan algorithm typically is very good at treating time gaps; for small gaps the detrender fits the data over the gap without distortion, and for large gaps the detrender creates independent fits to either side of the gap. However, in some exceptional circumstances the detrender creates a poor single fit across a gap in the data. The typical result in the detrended light curve is an upwards or downwards hook, on either one or both sides of the gap in time.
We detect such artefacts automatically by finding any temporal gap longer than 2 hours and calculating the mean flux for the 2 hours before and 2 hours after the gap. A difference of greater than 4 standard deviations (calculated from the light curve detrended up until this point) is considered significant . In this case we cut out all data on the preceding and proceeding 0.8 days of the gap.

Jumps in flux
In Kepler data a large gap in time may often be associated with a large jump in flux. An example event would be the re-orientation of the spacecraft between quarters, causing the target to fall onto a different CCD. Some of these undesirable parts of the lightcurves are removed by the Lightkurve quality flags, but not all. Such an event is typically well-handled by Wotan, or as a backup is identified based on the gap in time in Sect. 2.8.1.
An event not well-handed by Wotan is an offset in the flux over a short timespan of what is otherwise a smoothly varying light curve. Because there is little gap in time, Wotan attempts detrending with a single fit across the jump. The fit before the jump is distorted downwards, whereas the fit to the right of the jump is distorted upwards. Subtracting this away yields a light curve with differential "hooks" on both sides of the data Power Period (days)  , where the transit duration is primarily a function of the binary's orbital phase at transit (see Fig. 2). The minimum and maximum lines corrspond to the shortest and longest possible transits. The 50% line is a transit duration that is longer than half of all possible transit durations randomised over the binary phase. The 75% line is calculated similarly. On the left is for the light curve after only the eclipses have been removed, i.e. before any detrending. The largest peak is at 7.78 days, corresponding to the rotation period of the primary star ) and slightly longer than the binary orbital period of 7.45 days. The second peak is at just less than 200 days, most likely corresponding to a trend between quarters of Kepler data. The middle panel is the light curve after applying the cosine detrending filter (Sect. 2.6). The vertical scale has shrunk significantly because those large spikes of long-period variability have been removed. The horizontal scale is zoomed into shorter periods, where there is now the highest periodogram power. Very little stellar variability remains, except for one spike at 0.7 days which is above a 1% false alarm probability. The right panel is the light curve after a subsequent application of the variable Tukey's biweight filter, which removes this spike. What remains is a flat (white noise) light curve plus planet transits.
(unlike in Sect. 2.8.1 where there may be a hook on just one side of the data). We identify and remove these events automatically. We loop through the light curve and for each timestep we calculate the mean flux over the preceeding and proceeding 2.1 hours. A "jump" is identified at a given timestep if the following three conditions are met. First, the mean flux both before and after the time step must be individually more than two standard deviations from the mean flux of the detrended light curve. Second, the absolute difference between the mean flux before and after must be greater than four standard deviations (like in Sect. 2.8.1). Third, the mean flux before the time step must have a different sign to the mean after the time step, i.e. there must be one upwards hook and one downwards hook. This last step helps us avoid removing short transits.
As in Sect. 2.8.1 we remove the preceding and proceeding 0.8 days of data around any time step identified as a jump.
2.8.3. Times of low flux common to many different binaries By this stage in the detrending most non-transit signals have been removed. However, it was discovered that sometimes even transit-looking signals were undesirable.
Cross-talk between different Kepler pixels is an instrumental defect which can make transit-like signature appear in multiple different targets (not just binaries) at the same time.
To find such events we ran the detrending algorithm on 270 of the binaries from the Windemuth catalog, with cuts made to only include binaries between 7 and 50 day periods (similar to the known Kepler circumbinary planets, Welsh & Orosz 2018;Martin 2018). For each binary we outputted the times corresponding to the 9 points of lowest flux, where we avoided listing multiple times within a 1 day bin. We then created a histogram of the deepest points across all these binaries, using daywide bins, in Fig. 5.
We discovered 35 low flux times that were common in more than 5 of the 270 eclipsing binaries, which we determined to be more likely instrumental than random clumping of real astrophysical events via the analysis of Fig. 6. In Fig. 5 we demonstrate one event which was seen at time 1011.8 days (BJD 2,455,000) and is transitlike in appearance. It was impressively seen in 134/270 binaries.
We created a list of the 35 low flux times which is published online to warn other planet-hunters. This list is incorporated into the detrending algorithm such that for every target all data is removed at these times. It is possible that by removing ∼ 35 days of data we remove some legitimate, unique transits of circumbinary planets, however it would be difficult to be convinced of their authenticity if many other binaries had a dip at the same time.

A final semi-automatic cut
A final check is made of the light curve for any other events which are clearly spurious. The code automatically produces a list of potential events and then human input is used to decide which (if any) are to be removed. Potential artifacts are produced by the code in two ways. First, the 9 deepest flux points are listed as in Sect. 2.8.3. Second, any times are listed where there is an eight standard deviation change in flux from one data point to another. Such an instance could correspond to a deep transit (or a less deep transit with a sharp in/egress), which is why the code simply displays these sections of the light curve rather than removing them automatically. A record is made of any sections of the light curve which are manually cut. When testing Stanley and applying it to the known circumbinary systems we typically perform less than five manual cuts (and often zero). The duration of these manual cuts is typically similar to that of a transit, i.e. hours, and hence representing less than 0.1% of the total light curve.

SEARCH ALGORITHM
The main aspect of the Stanley code is to identify small transiting planets. We conduct a brute-force search over a static grid. For every set of parameters (masses, radii and orbital elements) we use an N-body integrator to calculate a mask of transit times and durations over the four year Kepler mission. We then phasefold the photometric data on these variable transit windows. Finally, we determine if the phase-folded light curve corresponds to a significant dip in the data.
Unsurprisingly, this method is slow and computationally expensive. We would encourage the application of more sophisticated algorithms such as genetic algorithms and simulated annealing. However, by indiscriminately sampling the entire grid we at least know that the best solution we find will be a global one, not local.

Defining the search grid
Both the binary and planetary orbits have six Keplerian orbital elements: period T ; eccentricity e; argument of periapse ω, inclination I, longitude of ascending node Ω and starting orbital phase, which we characterise with the starting true longitude θ 0 . We therefore have 12 orbital elements. Both of the stars and the planet have a mass and radius, bringing us to 18 different parameters to define our search grid.
Stanley is written such that a brute-force search over this 18-dimensional grid is possible. However, for a blind planet search, this is highly unrecommended. In particular, using an N-body algorithm at each point would be prohibitively computationally expensive. We therefore need to have a set value for many of these 18 parameters to reduce the search grid dimensionality, and to optimise the search grid for the remaining parameters.

Parameters we ignore
Some of these 18 parameters can be "ignored", by which we make an informed simplification of that parameter to a set value.
• R p = 0. Our algorithm does not fit transit shapes, but rather stacks together flux within a calculated transit interval. The depth of that flux will of course depend on the planetary radius, but it is not something we fit for or need to know a priori. The width of the interval, which is important for us, is dependent on R p , but this is a small factor since R p R A,B .
• m P = 0. The mass of the planet affects the dynamics of the 3-body system, which in turn affects the timing of transits and eclipses. However, only in about half of the known circumbinary systems were eclipse timing variations measurable. These eclipse timing variations were caused by gas giants and only on the order of minutes (less than one observing cadence), and hence can be safetly ignored, particularly when searching for low-mass planets.
Perturbations on the binary's orbit could also have a back-reaction on the planet's orbit, varying its transit timing, however this may be ignored for detection.
• Ω bin = 0 • . When considering transits we are not sensitive to individual values of Ω, but rather just the difference ∆Ω. It is typical to set Ω bin = 0 • , reducing the dimensionality by one.
• I bin = I p = 90 • and ∆Ω = 0 • . We assume completely flat, edge-on systems. The binary is known to be near 90 • since it eclipses, and any planets will have an inclination near 90 • when transiting. However, this flat assumption can be problematic.
A misalignment as small as one degree can affect the transit signature. The 3D orbital orientation of any misaligned planet precesses (a rotation of ∆Ω). This changes the transit impact parameter b, even bringing it above 1, i.e. the planet may . Poisson distributions are shown for comparison, which assume the events are random in time, as they would be if they are astrophysical. The bin from 0-1 shows the number of 1-day intervals with no low points, the bin from 1-2 corresponds to 1 low point, etc. Over the 1467 bins of daily data, 2430 low points were identified (9 in each of 270 light curves), for a mean rate of λ = 1.656 per bin. The Poisson distribution using this rate, normalized to 2430 low points, is dotted. It is ill-fitting, showing we can rule out that all the hypothesis that all the points are from a random-in-time distribution. Instead, if we suppose that some points are random and others indicate instrumental problems at particular times, then we can use the first two bins to estimate a mean rate of λ = 478/599 = 0.798 over 1306 low points, which is shown as the dashed Poisson distribution. On this basis, we identify time bins with 5 or more low points (e.g. 1011.8 days in Fig. 5) as more likely clumped than random-in-time, hence more likely instrumental than astrophysical.
stop transiting (Schneider 1994;Martin & Triaud 2014;Kostov et al. 2014). By creating a transit mask with a flat system we therefore assume transits on every passing. We therefore risk folding noise from a transit epoch with b > 1 on top of true transits, diluting the signal. To model and hence try to avoid this effect we would have to expand the search grid dimensionality by two (I p and ∆Ω), and use a very fine grid to precisely calculate the transit and non-transit epochs. This would be computationally prohibitive. Furthermore, for small circumbinary planets even if you could model the inclination effects it may still be impossible to detect it with significance if lots of its shallow transits are missed. With this in mind, we pursue this flat model and accept we may be less sensitive to misaligned planets. In Sect. 4.1 we show that with a flat model we can still re-discover misaligned planets which only have partial transitability (Kepler-47d, -413b and -1661b).

Parameters with a measured value
As described in Sect. 2.2, the Villanova and Windemuth catalogs combined fully characterize the two binary stars and their Keplerian orbit, i.e. m A , m B , R A , R B , T bin , e bin , ω bin and θ 0,bin . We only take the bestfitting value of each parameter. It is understood that the photometrically-derived stellar masses from Windemuth et al. (2019a) are only accurate to ∼ 20%, and these masses affect the dynamics and transit timing of the planet. However, we will demonstrate the ability to recover the known circumbinary planets despite using Windemuth et al. (2019a) stellar masses which are known to differ slightly from the true values. Stanley is also set up to search over a range of stellar and binary parameters if desired.

Parameters we search over
The remaining parameters in our blind planet search are P p , e p , ω p and θ 0,p . We take care to optimize this grid such that it is both a) as computationally-efficient as possible and b) as unbiased as possible.
We use a grid of planet periods that is neither constant in linear or log space. We instead optimize it based on the expected transit duration, similar to what Ofir (2014) did for single stars. The planet period grid spacing is described by δP p (P p ), i.e. the step size varies with P p . To calculate δT (P p ), we first calculate the relative velocity between the primary star and planet using Eq. 5. The shortest possible transit duration is corresponding to the planet and primary star moving in opposite directions with respect to the observer at transit. We set δP p = 3τ min . Our period grid is therefore derived with the precision to hit even the smallest duration transit within a factor of three. The factor of three is because the code ultimately allows slight deviations of the transit times from the N-body calculations up to three times the predicted duration. This small "sliding" effect is largely to account for imprecise measured stellar parameters and is further described in Sect. 3.3. The minimum period searched corresponds to a semimajor axis of where interior orbits would be almost guaranteed to be unstable (Dvorak 1984;Holman & Wiegert 1999;Mardling & Aarseth 2001;Quarles et al. 2018). The longest period searched is two years, since we demand at least three primary transits within the four years of Kepler data. This will inhibit detection of planets such as Kepler-1647b, which passed the binary twice during the four years of Kepler (Kostov et al. 2016). This is a gas giant though, and we imagine that finding Earthlike circumbinary planets with only two transits would be nearly impossible. In Sect. 4.1.3 we loosen the detection requirements and demonstrate we can then recover Kepler-1647b.
The spacing in the starting planet phase, θ 0,p also varies with the planet period P p , and is calculated by The θ 0,p grid becomes finer at greater periods. For the planet eccentricity and apsidal alignment traditional search methods either have uniform grids in e and ω or uniform grids in e sin ω and e cos ω (or possibly a slight variation such as √ e instead of e). We propose that both are problematic.
If δe and δω are constant then the effects of a large eccentricity are relatively under-sampled (and those of a small eccentricity or even circular orbit are oversampled). On the other hand, whilst using a square grid of δ(e sin ω) and δ(e cos ω) does avoid that problem, it means that the range of e sampled varies with ω varies by a factor of √ 2 10 We apply a novel technique of a circular grid in e cos ω and e sin ω. We have constant steps in eccentricity out to a maximum value that is independent of ω. For a given eccentricity we then define δω = δe/e in radians, which is the equivalent to stepping around the circumfrence of a circle in e cos ω vs e sin ω space, where the step size along the circumfrence is the same for every e.
For our search we use δe = 0.067 between e = 0 and e = 0.2. This leads to four steps of eccentricity and 34 different combinations of (e, ω), where higher eccentricities having more steps in ω according to our circular grid. Our range in eccentricity covers all of the known planets, where all but one has a small eccentricity less than 0.1. The grid spacing δe = 0.067 is a somewhat arbitrary choice, but was found to re-discover the known circumbinary planets. A finer grid would significantly increase computation time. On the other hand, the assumption of circular orbits was found to be ineffective, even for known planets with very small eccentricities. We believe the cause to be a non-negligible amount of apsidal advance over four years of Kepler observations causing the transit times to shift, even for the most circular orbits. Three-body interactions also prevent planets from maintaining perfectly circular orbits.
Overall, we typically use between 10 and 100 million individual transit models when searching for a Kepler circumbinary planet around any given target.

Running the N-body simulation
10 For a square grid of δ(e sin ω) and δ(e cos ω) the distance from the origin to the grid point is equal to the eccentricity. At say ω = 45 • , i.e. in the top right corner of the square grid, the maximum eccentricity value will be √ 2 times that at ω = 0 • along the horizontal axis of the grid.
For a given lightcurve we run the Rebound IAS-15 N-body integrator (Rein & Liu 2012;Rein & Spiegel 2015) on each set of parameters in the search grid to produce a set of transit times and durations. This is the slowest part of the algorithm, and is often best suited to a computing cluster, which is covered in Sect. 3.6.
Each N-body integration is run over the entire timespan of the light curve, which is typically about four years for Kepler. Transit times are calculated using an iterative procedure, where the code notes a sign change in the difference in the horizontal positions between the planet and primary star, and then steps backwards and forwards with increasingly small steps until the solution is converged upon. This is the default means of calculating TTVs in Rebound 11 . Note that we ignore light travel time effects, which should be no more than seconds in our examples ) and hence can be ignored for our purposes. We only calculate transits on the primary, and discuss this limitation in Sect. 5.4.
The transit duration is then calculated as in Eq. 4 based on the relative velocities of the primary star and planet, and the size of the primary star. Recall that we assume edge-on, coplanar systems so b = 0, and the planet radius is considered negligible and hence set to zero.
Even though a stability limit criterion was invoked in creating the range of planet periods tested (Eq. 7), such an equation is simply an approximation and does not take into account finer complexities such as the relative planet and binary phases, or structures in the stability map such as mean motion resonances. Consequently, some tested planets may still be unstable. We invoke a very simple test for stability in the N-body integration. At 50 equally spaced timesteps along the integration we calculate the planet's osculating Keplerian eccentricity and compare it to its starting value. If it varies by more than 0.1 we consider this system to be unstable and stop the integration.
We emphasise that this is a very rough instability approximation. In a four year integration we are typically not seeing instability fully unfold (i.e. the planet being ejected) but rather a change of 0.1 is seen as a sign that instability may come. A true test of instability would involve integrations for thousands or even millions of years (e.g. Quarles et al. 2018), which is not feasible with an N-body brute-force search algorithm. We recall that the purpose of this algorithm is an initial identification of transiting circumbinary planets, and any candidates found would be expected to receive further scrutiny, including of their long-term stability.

Fitting the transit mask
For every set of circumbinary parameters the N-body code outputs a sequence of transit times and durations, which we then fit to the light curve. If our search grid contained an exact corresponding set of parameters then the N-body model would produce a perfect fit. In practice though, this is unlikely for two reasons. First, our stellar parameters are only precise to ∼ 20% (Windemuth et al. 2019a). Second, computing time prohibits constructing a grid of planet parameters sufficiently fine to produce transit times precise to sub-observing cadence (< 30 minutes).
To account for any slight imperfections in the N-body model, the algorithm loops through each of the transits and creates a window of the light curve centred on the N-body transit time, where the width of this window is three times the N-body transit duration (3τ , Eq. 4). We then slide a smaller window of width τ to find the time corresponding to the lowest mean in-transit flux. We then note this as the transit time. We do not vary the transit duration.
For each transit we take the in-transit flux data and subtract from it the mean out-of-transit flux (within the 3τ window). We then add these relative transit fluxes to a "global transit", i.e. a circumbinary transit which has been phase-folded on a variable transit interval with variable transit durations. By looking at a relative flux between inside the transit and slightly either side of it, as opposed to an absolute in-transit flux, we can account for any remaining stellar or instrumental variations which were not removed by the detrending process.
The SNR of this global transit is calculated to be SNR =f wheref is the mean in-transit relative flux for the global transit, N is the number of in-transit data points and s oot is the standard deviation calculated over the entire light curve excluding the in-transit data. We illustrate this procedure in Fig. 7 for four transits of Kepler-1661 and Kepler-47. The N-body code initially predicts a transit time and duration centered in each window, and then we slide a transit mask along the window to find the minimum flux. In these eight transits the initial transit time is about half a transit duration off the true centre. These corrections are on the order of ∼ 2 − 4 hours, which is small compared to the TTVs on the order of ∼ 16 hours and ∼ 50 hours for Kepler-1661 and Kepler-47b, respectively. The initial N-body fit of the transit time and duration is demarcated by the red dashed lines. Each transit window is equal to three times this duration, and hence some windows are longer than others. The red shaded region bounded by solid lines is the final transit picked out by the algorithm, which is selected within this window to minimize the in-transit flux. For Kepler-1661 the planet did not transit for the first half of the Kepler mission due to a slight but non-negligible misalignment of ∼ 1 • . In the first panel this is the epoch where Stanley predicts a coplanar planet would have transited. For the later three epochs the planet has precessed into transitability.
These two sets of transits are chosen to demonstrate some idiosyncrasies which the algorithm sometimes has to deal with. For Kepler-1661, the planet is mutually inclined by ∼ 1 • with respect to the binary. Consequently, the planet initially does not transit until the second half of the Kepler mission, and then the transits increase in depth as the impact parameter decreases. As discussed in Sect. 3.1.1, we assume coplanar orbits since scanning over a grid of 3D orientations would be computationally prohibitive. Here we see the drawbacks of this: we fold a non-existent transit near 634 days (BJD -2,455,000), and also we over-estimate the transit durations.
For Kepler-47b Fig. 7 shows evident (but shallow) dips for this 3R ⊕ planet in the second, third and fourth panels. In the first panel, the light curve looks flat. In fact, a transit does exist near 256 days, but it was unfortunately detrended away as seen earlier in Fig. 3. This planet transits 24 times and so the transit timing model is sufficiently constraining to predict a transit here.
Despite these small challenges, Sect. 4.1 shows that both planets are detected with high significance.
We note here a caution for the sliding method. Take one extreme, where we use solely the N-body algorithm to predict transit times and do no sliding. Since these times are entirely based on a circumbinary planet model, then if they line up with actual dips in the light curve we would be rather inclined to believe a planet has been found. However, as illustrated in Fig. 7 the N-body transit times can be slightly off even for known, legitimate planets, and hence such an algorithm would be inefficient. At the other extreme, one could completely avoid using an N-body model to predict transit times and simply scan through the light curve looking for the deepest flux points 12 . This would be very efficient, and indeed guaranteed to find all of the dips in the lightcurve. However, there would be no physical model for the timing of these dips, and hence we would be less convinced we had found a planet and not just artifacts of poor detrending. Indeed the QATS algorithm follows this method to an extent, and is thought to suffer from this for circumbinary planets  and see Sect. 5.1.4).

Applying quality cuts
Several cuts are made by Stanley to reduce the likelihood of false positive solutions. First, when fitting each N-body transit model we exclude any transits for which more than 25% of the data is absent (based on the expected number of data points for a given transit duration). We then also do this on a global scale for the entire transit model, and demand that at least 50% of the expected data points are in the light curve.
These cuts avoid the algorithm artificially finding low mean flux solutions that are driven by one or two very deep outlier events and then the majority of the transit mask falls in gaps, and hence not affecting the mean flux. Such situations would likely be false positives, or at best true planet transits that we would never be able to confirm with the existing data.
The second criterion is that we require there to be at least three detected transits that are roughly consistent. Specifically, we calculate the SNR for each individual transit, determine the maximum value SNR max , and demand that at two of the other transits have a SNR greater than 0.45 × SNR max . Any transits that fall in a gap are ignored for the sake of consistency.
By requiring at least three roughly consistent transits we avoid selecting a solution which is solely driven by a small amount of deep events, combined with much larger number of insignificant transits. Such inconsistent events would be likely due to a completely non-periodic phenomena such as flares or artifacts in the light curve detrending. Alternatively, they could be produced by real transiting circumbinary planets, but with periods so long that they fall out of our search grid and would be hard to confirm anyway, or at least require the differing methodology of  to confirm them. Whilst our results typically use this three transit criterion, in Sect. 4.1.3 we show that by loosening it we can recover the 1108-day Kepler-1647b.
We emphasise that the three transits only need to be roughly consistent. The 0.45 SNR factor is somewhat arbitrary but allows for variations in transit duration and depth (and hence SNR) that are common in circumbinary planets, and works to recover the known planets. Loosening how consistent the transits need to be may allow for more marginal detections, for example with grazing transits, and this criterion will be easily editable in the code for the user's preference.

The best fitting solution and the signal detection efficiency (SDE)
We construct a variant of the commonly-used Signal Detection Efficiency (SDE) 13 to pick out the best solution and inform us of how significant it is compared with other tested transit models.
First, at each tested planet period we marginalise over all tested orbital parameters (which, in our work is typically θ 0,p , e p and ω p ) to find the transit model with the maximum signal to noise ratio of each fit SNR(P p ).
There is typically a decreasing trend in SNR as a function of P p . This is because the number of in-transit points N t decreases as longer period planets transit less 14 . Following Hippke & Heller (2019) we remove this trend to avoid (to the best of our ability) a bias towards short-period planets. The Tukey's biweight filter in Wotan ) is employed with a window length of This window length is derived based on the analytic TTV calculations of Armstrong et al. (2013) and should be ∼ 4 times the width of the planet-induced peaks in the SDE curve (in Fig. 9). The SDE(P p ) is then defined at each planet period according to  as where s SNR is the standard deviation of the function SNR(P p ). This standard deviation is calculated in rolling segments (with width ξ, Eq. 10) of SNR(P p ) and then we take the median value of all of the standard deviation calculations. This means that in calculating s SNR we avoid including any spikes in SNR due to transiting planets. In Fig. 8 we illustrate both SNR(P p ) and SDE(P p ) for Kepler-34. The best-fitting period is the one with the highest SDE. The values for θ p , e p and ω p are those with the highest SNR (Eq. 9) at this period. Note that the best fitting solution works in the same way if the search grid is expanded to include any of the other orbital and/or stellar parameters.
Hippke & Heller (2019)  . Ultimately a lower SDE threshold will increase survey completeness and provide a greater number of candidates to follow-up, but with the risk of added false positives. This will be quantified more thoroughly in our second

Planet Period (days) Signal Detection
Efficiency (SDE) Signal to Noise Ratio (SNR) Figure 8. Example of the Signal to Noise Ratio (SNR) from Eq. 9 for Kepler-34 in the top panel, and the conversion to Signal Detection Efficiency (SDE) from Eq. 11 in the bottom panel. The red dashed line is the trend removed using a Tukey's biweight filter in Wotan, since shorter planet periods will systematically have higher SNRs since more data points are captured in transit.
paper when applying Stanley to a larger sample in the search for small planets.

Computation times
Stanley can be run on either a single computer or on a computing cluster. The code is set up to split the search grid over P p , which typically has the most steps, and run in parallel. The results are ultimately combined seamlessly. To recover the known planets over a period grid spanning up to two years takes typically a couple of 100 computational hours, and so for a blind search a computing cluster 15 is highly recommended. Given typically 100+ million independent N-body simulations are 15 or a good book run for each planet search, this time is not surprising. Reducing the longest tested period can significantly reduce computing time (advice we follow in Sect. 4.2) because the longest planet periods require the finest grid in θ p . Fitting circular planets would also expidite the fit, but in our experience a circular orbit often failed to fit most of the transits, no matter how fine the grid in P p and θ p,0 was, although we have not strenuously benchmarked this.
We refer the reader to Rein & Liu (2012); Rein & Spiegel (2015) for details on the Rebound IAS-15 Nbody integrator. We note here that it has an adaptive timestep which is primarily set by the shortest orbital period, i.e. of the binary. A given simulation of the 7.4day Kepler-47 binary will therefore be slower than the 41.1-day Kepler-16. However, our adaptive grid is more dense for longer period binaries, which somewhat evens out the overall computing time.

Main results
There are ten published Kepler circumbinary systems. We run the Stanley detrending and detection algorithm on all of them and catalogue our fits in Table 1. All Kepler systems yield a significant detection with an SDE above 12, except for Kepler-1647 as expected, since its 1108-day planet does not transit the mandatory three times (but see Sect. 4.1.3). In Table 1 we show our fitted parameters to each system and the detection significance.
In Fig. 9 we show the SDE for each system. A few of features stand out. For all but Kepler-1647 there is at least one prominent spike, significantly above an SDE of 8. It is noticeable that this spike is rather broad. The main cause of this breadth is that the SDE at a given period is the maximum value marginalised over all different values of θ p , e p and ω p . With these degrees of freedom, even if the planet's period is changed by a few per cent, a significantly different θ p , e p and ω p could contort the transit timing to still fit the data, and hence provide a similarly good SDE. The non-Keplerian perturbation on the planet's orbit by the binary also "smear" the best-fitting solution, as our solutions only indicate the osculating orbital elements at time t 0 , and they will change significantly over four years.
Another feature seen in the SDE plots are secondary peaks at alias periods of the true peak. Such peaks correspond to solutions where either a sub-set of the transits are found, or all of the transits are found, in addition to some noise.  Figure 9. Signal Detection Efficiency (SDE) for the re-discovery of all ten known circumbinary planet systems. Inset in each plot is a representative transit as discovered by Stanley, with red in-transit data points and blue out-of-transit points. Kepler-38 is the only case where its true period of 105 days has a significant SDE of 12.2, but it is smaller than the SDE = 14.7 at double the period. This occurs because the planet-binary period ratio is close to 5.5, and the planet's starting phase awkwardly causes its early transits to alternatively be near the primary and secondary eclipses of the binary. Its second, third and fourth transits are all blended with eclipses and hence will not be detected. Transits near secondary eclipse, even if not blended, are of maximum duration (θ A = 90 • in Fig. 2) and hence are at risk of being flattened by the detrending process. By visually comparing the transit fits at 105 and 210 days it is clear that the 105 day model does include additional, legitimate transits and hence should be favoured. This is an example where some manual intervention may be needed in what is otherwise an automated process. Overall the model predicts 12 out of 13 primary transits for Kepler-38, even though some of these transits were removing by the detrending algorithm. This demonstrates that circumbinary models are highly constraining if you are able to find them.
For each re-discovered planet in Fig. 9 we show an example transit as discovered by the algorithm, i.e. with red in-transit points and blue out-of-transit points.

Multi-planet Kepler-47
Kepler-47 is the only known transiting multi-planet circumbinary system, containing three Orosz et al. 2019). From the initial run of Stanley Kepler-47b has the most prominent SDE signal, just above 15. With a radius of 3R ⊕ it is smallest planet in both the system and indeed of all known circumbinary planets. However, it is the innermost planet with a 49 day period and so it transits many times more than the outer planets, boosting its detection signal.
For Kepler-47 we cut out these detected transits of planet 'b' from the light curve and then re-run Stanley. A clear detection of the outer-most planet 'c' 16 is then seen, also at an SDE of roughly 15. We repeat the process of removing these detected transits and rerun Stanley, at which point the middle planet 'd' is revealed. Curiously it also has an SDE of ∼ 15. Finally, after removing the transits of 'd' there are no additional detectable planets. We show all of the SDE plots in Fig. 10 and an atlas of the transits of each planet in Fig. 11. 16 Confusingly, the outermost planet is called 'c' and then the middle planet is called 'd', since the latter was not confirmed in the original paper , but only later by Orosz et al. 2019. In Fig. 11 it is noticeable that the middle planet has six detected transits which grow deeper over time due to nodal precession, and are only easily visible by eye near the end of the Kepler mission. When the first paper  was written only one "orphan" transit of this middle planet had been detected, and so it was only a candidate and not confirmed until Orosz et al. (2019). The outer planet transits are more constant in depth. Actually four transits of 'c' exist in the data but one of them is blended with an eclipse, which causes it to get removed by our algorithm. Overall, Stanley is the first algorithm with a demonstrated ability to detect multi-planet circumbinary systems.

Exceptional case of Kepler-1647
Using the standard criteria of our algorithm, i.e. three transits on the primary star, we do not detect the 1108day Kepler-1647. The highest peak of the SDE in Fig. 9 is at 7.5 with a period of 134.8 days. This corresponds to one of the obvious (0.2%) dips in the light curve, folded together with some noise.
To test how Stanley performs with looser constraints, we re-run it with a requirement of only two transits. We also change the transit detection to be from the primary star to the secondary star, as the two deep transits for Kepler-1647 are on the secondary, with only one shallow, grazing one occurring on the primary.
The resulting SDE is shown in Fig. 12. Stanley finds the true planet period of 1108 days with a large SDE spike above 30, corresponding to the two deep transits on the secondary star. Curiously, there is a slightly higher SDE peak at roughly half this orbital period. This solution contains the same two transits, plus one transit epoch occurring in a gap. There are also many other "significant" looking signals that are a result of the weakened two transit criteria.
Ultimately in this case we can say that Stanley has recovered the known planet, but with a degenerate period. The initial discovery of these transits in Armstrong et al. (2014) had the same period degeneracy, but the confirmation of the planet in Kostov et al. (2016) using a full photodynamical model and the shallow but significant primary transit confirmed the 1108 day period. Given the 554 day model does not contain additional transits, one could argue that Occam's razor favour's the longer-period model, which is why we include it in Table 1.
Ultimately, it is up to the user how they use the code and the criteria are programmed to be user-editable. A two transit criterion will allow for more long-period and marginal planet detections, but at the expense of reliability. Kepler-1647b is also larger than Jupiter and;

Inner and outer planets (b,c) removed
All three planets (b,c,d) removed Figure 10. SDE for Kepler-47 where we iteratively remove the transits corresponding to the highest SDE. The top panel shows the fit using all of the data. This process shows that we can recover all three planets with a significant SDE∼ 15. We also see spikes from aliases of each planet period, which all disappear when the transits have been removed.

Inner planet (b)
Middle planet (d) Outer planet (c) Figure 11. An atlas of transits for the inner 49 day planet b (red), the middle 187 day planet d (green) and the outer 303 day planet c (blue). The plots are colour-coded to match the SDE spikes in Fig. 10. Figure 12. SDE for Kepler-1647 with a modified detection criteria of only two transits (instead of three) on the secondary star. The true period of 1108 days corresponds to the second highest peak. The highest peak at roughly half this period is for a model containing the two real transits and a gap transit epoch. We can therefore recover the known planet, but with a degenerate period.
it is unlikely that small circumbinary planets could be detected with just two transits on typically faint Keplertargets. Finally, we note that the main aim of the Stanley code is to phase-fold transits with a variable interval, and this is only relevant for three or more transits.

Scaling planet transits
We artificially reduce the depths of known transits to test Stanley's limits of detectability. To scale the transits we first calculate the smooth, noise-free trend λ(t) of the light curve using a Tukey's biweight filter in Wotan, where the known transits have been masked to be ignored by the detrending filter. The in-transit lightcurve flux, f * (t) is then scaled at each time step according to where f is the original flux and N is Gaussian random noise with a standard deviation equal to where s oot is the standard deviation of the flux out of transit, after eclipses have been removed and the light curve has been received an initial detrending. The noise is added to counter-act the effect of the original noise being effectively decreased as the transit is scaled down. In the limiting case of R * p = 0 the noise added is equal to the noise of the entire light curve.
Note that for planets with visible transits on both primary and secondary stars they both get scaled down proportionally. For the multi-planet system Kepler-47 we are only attempting to detect the innermost planet  Figure 13. Signal Detection Efficiency (SDE) for known circumbinary planets with reduced transit depths. Each of the ten known circumbinary systems is shown with a different line. The cross is the SDE for detection of the original data (i.e. the value in Table 1). For each system the planet radius is reduced from its known value down to 0R⊕ in integer steps, which proportionally scales the transit depth according to Eq. 12. Circles indicate the SDE with scaled transits. Markers are coloured red if more than 50% (typically more than 80%) of the primary transits have been recovered by the detection algorithm, and blue otherwise. The black horizontal dashed line at SDE = 8 is a rough limit, above which a detection is considered significant. All of the circles below this line are blue, and all but one above this line are red (see text in Sect. 4.2 for outlier explanation). The SDE is therefore seen to be an accurate indicator of when the algorithm successfully detects the planet transits. Kepler-16 would have been detectable down to 1R⊕. and hence only scale down these transits, not those of the outer planets.
For each system we scale down R p in steps of 1R ⊕ and re-run the detection algorithm. Since we are only focused on the known planets, typically near the stability limit, we use a shortened period grid concentrated near the edge of the stability limit: 2.2a bin to 4.1a bin , and a grid in θ p which is fixed rather than adaptive, which will not have a large effect over a small period range, and if anything will yield a slight under-sampling and hence our results are slightly conservative. This is all for computational expediency since we are not interested in planets up to the original 2 year period limit.

Finding smaller circumbinary planets
The detrending and search algorithms are then run normally on these scaled transits. In Fig. 13 we show the signal detection efficiency (SDE) as a function of the scaled planetary radius. We consider the planet detected if > 50% of transits are detected, and this is indicated with a red circle (blue otherwise). Each system follows a similar curve, where the SDE decreases as the planet radius is decreased, until it gets below SDE ∼ 8 and the transits are no longer detected above the noise. At 0R ⊕ , i.e. the transits are completely replaced by noise, no detections are made, as should be the case.
Two small exceptions to this are seen. First, for Kepler-34 the transits of the 8R ⊕ scaled planet are anomalously not detected, whereas they are for smaller radii. This is because Kepler-34 is a roughly equal mass binary and the transits are of similar depth on both stars. For the 8R ⊕ light curve Stanley has in fact discovered transits of the planet, but curiously four of the five are on the secondary star, and hence it technically gets counted as a non-detection despite a large SDE. Accounting for both primary and secondary transits will be a future improvement (Sect. 5.4.2).
The other exception is that for Kepler-453 the 2R ⊕ and 4R ⊕ scaled transits are detected but not the 3R ⊕ . That is because this system only has three transits and the detrending algorithm happens to detrend away one of the three transits in the case of 3R ⊕ . It is forever a challenge to develop detrending techniques which remove unwanted signals yet preserve transits. In Sect. 5.4.1 we discuss a few pieces of future work which could aid this.
Overall we see that our modified SDE is a good predictor of when planets are discovered, and a limit of ∼ 8 is a reasonable threshold, similar to previous studies on single stars .
In Table 2 we catalogue the detection limits for each planetary system. We define the smallest detectable planet radius as the one corresponding to an SDE of 8, using a linear interpolation in Fig. 13. In most cases we could find planets where the transits have been reduced in depth by over 50%, and over 90% in the best cases . In four out of nine cases our algorithm pushes into the mini-neptune and super-Earth regime with R p < 2.5R ⊕ , which are the most common planets around single stars (Fig. 1). On the other hand, Kepler-38 was arguably the most difficult planet (could only detect 11% smaller planets) due to the difficulties arising from the planet's awkward phasing and near 5.5 period commensurability with the binary, as discussed in Sect. 4.1.

Detection limits as a function of period
The detection limits derived in Sect. 4.2.2 are defined for the orbital parameters of the known planet. To derive detection limits at other putative planet periods we take the empirically derived detection limit from Sect. 4.2.2 and scale it according to  Figure 14. Detection limits for known circumbinary planets as a function of planet period. We derive these curves using the emperically-derived empirical detection limit at the known planetary period (given in Table 2) and scaling it according to Eq. 14. Crosses indicate the true radius and period of each known planet, and are connected to the detection limit curves by dotted lines.
which is derived by assuming the number of data points in transit is proportional to 1/P P . This is a rough criterion, since we neglect the transit duration being a function of period. This would be simple to account for in single stars, but complicated in binaries where transit durations vary with binary phase (Sect. 2.4), and this affects both the detrending robustness and the transit durations and hence SNR of the transit signal. Overall, the detection threshold of a circumbinary system is not sharp; it is highly binary-system dependent, as was seen with Kepler-38 being surprisingly such a challenge.
In Fig. 14 we calculate Eq. 14 over a period range between 4P bin (an approximate stability limit based on Holman & Wiegert (1999)) and 4/3 years (guaranteeing at least three transits within the Kepler mission, neglecting gaps in the light curve). The resulting curves, with a shallow dependence on P p , are similar to those seen in Deeg et al. (2000). Aside for Kepler-38, -47 and -413, Stanley is sensitive to smaller planets at all stable periods up to 4/3 years.
From Fig. 13 the smallest detected hypothetical circumbinary planet is a 1R ⊕ scaled version of Kepler-16, with six out of seven transits detected and an SDE of 8.8. In Fig. 15 we show how the 1R ⊕ detection is barely above the noise limit, and compare it to the relatively easy detection of a 2R ⊕ planet. We have therefore demonstrated that not only an Earth-sized circumbinary planet, but one on a ∼ 220 day orbit in the habitable zone of its binary, is detectable with the Stanley algorithm. We temper this slightly by noting that Kepler-16 is brighter (visual magnitude of 12) and quieter than  Table 2. Detection limits for known circumbinary systems at the known planetary period based on the recovery of scaled transits. The true radius is taken from the respective discovery paper, and the detection limit is the smallest radius our Stanley algorithm can detect based on scaling the transits to a shallower depth. For Kepler-47b the detection limit is only calculated for the smallest, inner-most planet (see Fig. 16), but all three planets in this system are detected at an SDE> 15 (Fig. 10). Values are only given for the innermost planet in Kepler-47, and are not produced for the 1108 Kepler-1647, because Stanley does not find with the standard detection criteria (but see Sect. 4.1.3).
most Kepler eclipsing binaries, which is probably to be expected for what was the first circumbinary planet discovery.
The smallest detected circumbinary planet so far is Kepler-47 at 3R ⊕ . This is strongly detected with SDE = 17.2. However, in the scaled transit simulations in Fig. 13 no detection was made for 1R ⊕ or 2R ⊕ . In these simulations we scaled the depths of the inner planet but did not touch the outer planets. We re-run these simulations with the outer planets removed. We consider this reasonable since the outer transits are deeper and easily visible by eye (Fig. 11). We also take smaller step sizes of 0.1R ⊕ to see more precisely what our detection limits are. The results are shown in Fig. 16.
The smallest detectable scaled planet around Kepler-47b was 2.3R ⊕ , where 22 out of 24 transits were detected and the SDE = 8.7. The SDE vs R p curve flattens out around an SDE of 5. Curiously, for two of the solutions considered undetected by the SDE a subset of 16/24 transits were still discovered, but not in a way which distinguishes them from noise, so we do not consider these legitimate detections of the planet. We also note that Kepler-47 is significantly fainter than Kepler-16 (visual magnitude of 15.4 compared with 12), which would be the main contribution to Kepler-47's weaker detection limit.

Searching for additional planets in known systems
Out of the ten known circumbinary systems, only Kepler-47 has multiple (three) planets (Orosz et al. , 2019. Around single stars the observed planet multiplicity is a higher ∼ 20% 17 The inner stability limit requiring P p 4P bin , plus the propensity for planets to orbit binaries with periods longer than 7 days (a longer period than most Kepler eclipsing binaries) limits where additional planets may stably orbit. However, Quarles et al. (2018) show that roughly half of the known systems could host shorterperiod planet, squeezed between the known planet and the binary-induced stability limit. Kepler-1647 in particular is a prime candidate given its 1108 day known planet and 11 day binary. On the other hand, for exterior planets, all but Kepler-1647 could host an additional planet with a period less than two years, which allows us to search companion period ratios at least up to two, which covers most multi-planet systems around single stars (Steffen & Hwang 2015). Under certain conditions co-orbital circumbinary planets (i.e. a 1:1 resonance) are predicted by Penzlin et al. (2019), to which the algorithm is also sensitive.
For all ten systems we cut out all known transits by eye and then re-run the search algorithm. The results, by means of the SDE, are shown in Fig. 17. No new planets are detected.
The detection limits as a function of period in Fig. 14 show our rough sensitivity to interior and exterior additional planets, although we caution that the phasing of   Figure 15. Detection of Kepler-16 with transits scaled according to a 2R⊕ (top) and 1R⊕ (bottom) planet. The top plots shows the SDE, where the spikes correspond to the same period in both cases but for 1R⊕ the SDE of 8.8 is barely above the noise. On the left plots we show the transit fit signal to noise ratio (SNR) from Eq. 9 as a function of Pp and θp, calculated for the best-fitting ep and ωp. The plots on the right show an example transit and the SDE. The 2R⊕ planet is easily detectable. The 1R⊕ planet is detectable just above the SDE threshold of 8. For the 2R⊕ planet we annotate various features of the SNR parameter space. On the left the dark blue region corresponds to zero SNR since there are no acceptable solutions here, largely due to instability from being too close to the binary. The light streaks of high SNR correspond to fits including just one of the transits. These streaks are visible because each transit is individually significant. We use red circles to indicate where two and five transits overlap, i.e. the model transit times fit multiple transits. The five transit overlap corresponds to the peak of the SDE. For the 1R⊕ planet, we are able to discover the ensemble of five transits but the individual transits are not discernible above the noise.
transits is also important for their detectability. In each system in this figure we can rule out 4R ⊕ + interior planets. For five of the systems this interior limit is better than 3R ⊕ + and for three of the systems we were sensitive to 2R ⊕ +. For Kepler-16 we could have detected planets slightly smaller than Earth on both interior and exterior orbits out to almost 500 day periods. In every system in Fig. 14 we can rule out exterior planets larger than 6R ⊕ , and for many of the systems this limit is significantly smaller. We also rule out the commensurate planets predicted by Penzlin et al. (2019), unless they are significantly smaller than the detected planet at that period. It is possible that large companion planets do exist but on misaligned orbits that did not transit during Kepler. However, it is almost guaranteed that any such planets would eventually enter transitability (Martin & Triaud 2015;Martin 2017a), and indeed this happened with the middle planet in Kepler-47 (Orosz et al. 2019). Analysis  Fig. 13, red markers indicates more than 50% of transits being found and blue for less. The cross indicates the true 3R⊕ radius. The inset figure is the SDE for a 2.3R⊕ planet, which is the smallest planet considered detected.
of the TESS revisit of the Kepler field may reveal such new planets ).

Comparison with other past methods
Here we provide a qualitative comparison with some past methods for finding circumbinary planets. The codes for these other methods are typically not published so an in depth quantitative comparison of their efficiency is not possible at this time.

Traditional box least squares
For completeness we include traditional BLS (Kovacs & Mazeh) in this list, but the requirement of no or at least minimal transit timing variations makes it inappropriate for circumbinary planets. Without modifications, it could only be expected to find planets with coincidentally regular transit intervals, e.g. from a period commensurability with the binary or for a short period planet exhibiting many transits, some of which randomly happen to be regularly-spaced. The BLS algorithm would be aided by short binary periods, reducing the TTVs, although at this point the planet transits may actually be of longer duration than the binary orbit, adding complications. 5.1.2. Jenkins 1996 Eclipsing binaries were considered as early as Borucki & Summers (1984) as an ideal target for transit surveys, since the fact that they eclipse geometrically biases the transit probability. Schneider & Chevreton (1990);Schneider (1994) expanded on this idea, including the possibility of misaligned orbits. Jenkins et al. (1996) and its later implementation in papers such as Deeg et al. (1998), to our knowledge, created the first actual algorithm dedictated to transiting circumbinary planet detection. Despite being the oldest circumbinary-specific algorithm in this list, it is in fact conceptually the most similar to our own. Jenkins et al. (1996) uses a "Matched Filter Method", which is a crosscorrelation between the light curve and a large bank of model transit light curves. It is a brute force approach, using an N-body algorithm to produce the transit models, like in our paper. Their figures 13 and 14 are very similar to our Fig. 15. Their search grid of orbital elements is limited to the planet period and starting phase, but they do discuss the implications for neglecting other effects such as eccentricity. Unlike our Stanley algorithm, Jenkins et al. (1996) explicitly tries to match the depth of transits by including in their model a variable planet radius. The means of deducing a statistically significant detection also differ from our paper Their work was specifically targeted at 1.27 day eclipsing binary CM Draconis. It consists of two M-dwarfs,  Signal Detection Efficiency (SDE) Planet Period (days) Figure 17. Signal Detection Efficiency (SDE) for all known circumbinary systems with the transits of the known planets removed, in the search for new companions. No significant new detections are made. The occasional dips in SDE, most prominent in Kepler-35 at ∼ 80 days and Kepler-1661 at ∼ 120 days, are created by regions of parameter space where there are few or no transit models that can be fit to the data, due to stability or failing the criteria of three roughly consistent transits. That produces a large reduction in the SNR in that parameter space compared with surrounding regions, leading to a significantly negative SDE. They are not signs of "inverse" or "upwards" spikes in the data, such as might be created by flares or microlensing events.
which is a beneficial but rare configuration; small planets are more easily detectable. Jenkins et al. (1996) demonstrate a sensitivity down to 1.4R ⊕ . Their grid size consisted of about 40,000 unique transit models, which is a few orders of magnitude smaller than the grids we typically use in this paper 18 . We are not aware of any published application of this algorithm to a broader sample.

CBP-BLS -Ofir 2008
The Ofir (2008) algorithm -CBP-BLS -applies BLS to a light curve that has been treated to account for the barycentric motion of the binary. This shifts the location of potential circumbinary transits based on the binary phase. Indeed, the barycentric motion of the binary is the largest source of transit timing and duration variations, and by removing it the remaining light curve will be closer to strict periodicity. It does neglect, however, non-Keplerian perturbations which can shift transits by more than a transit duration. These unaccounted for effects will blur any detection. No application of CBP-BLS to a new planet search or the recovery of known planets has been published, although conference proceedings (Ofir 2015) state that the planets discovered at that time were all recoverable.

Carter & Agol 2013
The QATS algorithm by Carter & Agol (2013) is designed to find transiting planets with large transit timing variations. It has largely been applied to multi-planet systems around single stars, where in special cases mean motion resonances can make the TTVs on the order of hours. In theory, QATS can be applied to almost arbitrarily high TTVs, such as the extreme case of circumbinary planets where the TTVs may be ∼ 5−10% as large as the transit interval (Armstrong et al. 2013).
In QATS instead of fitting a fixed transit interval (or equivalently, an orbital period) you fit a maximum and minimum transit interval. At each transit epoch the explicit transit time can vary between these bounds. The downside to this method is that there is no physical basis for the variations; the variations are simply whatever best fits the data. This hampers detection reliability, and makes the method susceptible to detrending errors, particularly for long period planets .
We do note that our transit detection algorithm does use one "QATS-esque" element. Specifically, the threetimes widened transit window within which we find the 18 Although we have the benefit of two decades worth of computing improvements.
lowest flux (Sect. 3.3). However, whilst this sliding over a widened window is purely data driven, like QATS, it is a small correction relative to the overall circumbinary TTVs, which are physically determined by the N-body code.
No large circumbinary planet search with QATS has been published, but the algorithm has made important contribution to planets around single stars (e.g. Kepler-36, see Agol & Carter 2018). Kostov et al. 2013Kostov et al. (2013 was the first to present a novel application of the BLS algorithm (Sect. 5.1.1, Kovács et al. 2002) to circumbinary planets. The light curve is split up into segments with length corresponding to roughly a few times the binary orbital period, as this was a reasonable period at which to expect transiting planets based on both the stability limit and the known discoveries at the time of Kepler-16b, -34b, -35b, -38b and -47bc. Then within each segment BLS is run to find the most significant pair of transit-like dips. Given the segment length, this would correspond to consecutive transits on the same star (most likely the primary) separated by roughly the planet's period, and not a transit on each star separated by a few days. This method allows BLS to be applicable since the irregular transit intervals of circumbinary planets are not apparent when BLS is only used to find a single interval.

5.2.
The Kostov et al. (2013) algorithm then repeats this process with the flux of the light curve flipped, such that the BLS algorithm is therefore finding upwards "antitransits". A plot is produced showing the depth of all of the transits and anti-transits in the different segments of the light curve. If there is nothing but white noise then one would expect as many of each type of event, with roughly the same range of depths. If there are significant transits then these will stand out as being of greater depth. If the number of outlier dips in the light curve is greater than some merit function 19 then the algorithm triggers a human eye search to confirm that the transit signal is real and coherent across the entire light curve and not just pairs of transits in short segments. Kostov et al. (2013) portrays this as a "semi-automated" algorithm, although realistically most transiting candidates undergo a human inspection to some degree.
This algorithm was applied to yield the first discovery of Kepler-64, which was independently discovered by Planet Hunters using solely by-eye methods (Schwamb et al. 2013). Kostov et al. (2013) also presented and independent discovery and characterisation of the then recently published Kepler-47b and c. Impressively, the algorithm was also able to detect a 1.5R ⊕ transiting super-Earth injected into the Kepler-16 light curve, with 75% of the transits recovered. Improvements were later made to the method to make it applicable to misaligned planets with gaps in the transit sequence. This led to the discovery of Kepler-413 in Kostov et al. (2014). 5.2.1. Armstrong et al. 2014 Armstrong et al. (2014) is best known as the first comprehensive calculation of the occurrence rate of circumbinary planets, using the Kepler mission. However, it also contains a novel method for detecting circumbinary planets, and indeed found three circumbinary planet candidates which were later confirmed: Kepler-453 (Welsh et al. 2015), Kepler-1647 (Kostov et al. 2016), and Kepler-1661 (Socia et al. 2020). The method is to phase-fold the light curve on a fixed period and then search for significant dips within a wide window. Some transits would coincidentally be coherently stacked, but most would only be bunched together. The window was not arbitrarily chosen, but based on the TTV model derived in the (Armstrong et al. 2013, with inspiration from Agol et al. 2005). These TTVs account for both the binary's barycentric motion and apsidal precession of the planet. Armstrong et al. (2014) is one of only two large systematic circumbinary planet searches using an automated algorithm (the other being Klagyivik et al. 2017, Sect. 5.2.2). The consistency of the method allows for reliable occurrence rate constraints. The method is not, however, particularly sensitive to small planets. Indeed, constraints could only be made on planets down to 4R ⊕ (the 3R ⊕ Kepler-47b could not be recovered), and admittedly with large error bars near this limit. Klagyivik et al. 2017 Klagyivik et al. (2017) is a similar undertaking to Armstrong et al. (2014); a new method was both developed and applied to a large data set, in this case CoRoT, and constraints were placed on the occurrence rate of circumbinary planets. Their method may be considered a hybrid between QATS (Carter & Agol 2013) and Armstrong et al. (2014). A transit model is fitted allowing for variable intervals. This variation is fitted as a free parameter, but unlike in QATS for which the TTV bounds are arbitrary, Klagyivik et al. (2017) fit TTVs within the physical limits defined by the circumbinary TTV geometry (derived in Armstrong et al. 2013). Klagyivik et al. (2017) can therefore coherently stack the transits, unlike Armstrong et al. (2014) who bunched them in within a wide window.

5.2.2.
No planets were discovered from the CoRoT data and detection limits could only be placed down to 4R ⊕ . In addition to the poorer photometry than the Kepler mission, the short baseline of CoRoT is challenging for circumbinary detection, given their tendency for relatively long period orbits (Muñoz & Lai 2015;Martin et al. 2015;Hamers et al. 2016;Fleming et al. 2018). Overall, it would be interesting to see this method applied to Kepler.

QATS EB -Windemuth et al. 2019b
The most recently published method by Windemuth et al. (2019b) combines and expands several past techniques. First, the geometric positioning of the binary and planet at each transit epoch is accounted for using two independent Keplerians, creating a "regularized" light curve. The planetary orbit is assumed to be circular. The second step is to apply QATS, which accounts for any additional TTVs caused either physically (e.g. eccentric planets and non-Keplerian variations such as precession) or due to imperfect orbital and stellar parameters in the model. Compared with vanilla QATS, the transit timing models have at least some physical basis and are not completely arbitrary.
Compared with Stanley, the physical basis for the TTVs in QATS EB (independent Keplerians) is not as robust as using an N-body algorithm, but it does make the QATS EB algorithm significantly faster to run. Windemuth et al. (2019b) quote ∼ 5 minute run times per target on a single CPU, compared with ours taking typically tens or hundreds of hours. In fact, Windemuth et al. (2019b) presents their algorithm as a "compromise between the two extremes" of vanilla QATS and a fully brute force grid search such as ours. Windemuth et al. (2019b) recover known gas giants Kepler-35 and Kepler-64, with what they define as a SNR of 26.1 and 35.6, respectively. This is essentially the same as our SDE definition, which Table 1 lists as 65.2 and 16.7 for these two planets, respectively. However, our SDE peaks, as noted in Sect. 3.5, get broadened by searching over a grid of e p and ω p , whereas Windemuth assume circular planets. A closer comparison is is to calculate our SDE for the best-fitting e p and ω p , which yields SDEs of 68.6 and 26.7 for Kepler-35 and Kepler-64, respectively. Our Kepler-35 detection is therefore more significant whereas our Kepler-64 detection is a bit less. We note that differences may not be a function of the detection algorithm but rather different detrending methods. Windemuth et al. (2019a) also demonstrate a SNR = 20 detection of an injected 1R ⊕ planet in a system with orbital parameters similar to Kepler-47b. How-ever, there are two reasons why this is not comparable to our results in Sect. 4.2.2, for which 2.3R ⊕ was the smallest detectable planet. First, rather than using the true Kepler-47 light curve and scaling the transits, a synthetic light curve is created with white noise and injected transits. Given that a large part of the challenge in detecting small circumbinary planets is filtering out both physical and instrumental variations, the assumption of white noise is unrealistic. Second, the white noise used has a standard deviation of 5 × 10 −5 , which is ∼ 8 times smaller than the actual photon noise on the Kepler-47 light curve (calculated as the standard deviation of the out of eclipse and transit flux after detrending). So whilst QATS-EB is demonstrated to have a sensitivity to Earth-sized planets, it is in a highly idealized setting. Work is presently being undertaken to apply the Windemuth et al. (2019b) algorithm to a larger sample of Kepler eclipsing binaries (Windemuth, priv. comm.).

The human eye
It would be remiss of us to neglect the one method that has actually yielded all of the transiting circumbinary planets: actually has yielded new circumbinary planet publications -the human eye. Led in particular by Jerry Orosz, Bill Welsh, Veselin Kostov and Laurence Doyle, all of the known planets were painstakingly found by visually examining thousands of light curves. There has also been success with Planet Hunters, which crowd sources many human eyes. One of their first discoveries was the circumbinary planet Kepler-64b (Schwamb et al. 2013, and independently also published by Kostov et al. 2013). The complex transit signature of circumbinary planets is well suited to this platform.
We believe that our new method has advantages, but it is yet to be seen if our method can challenge the human eye for new discoveries.

Deeper optimization of the detrending efficiency
To detrend our light curves we combine existing methods with bespoke procedures specific to circumbinary planets. Some of our methods are admittedly ad hoc, and whilst they ultimately result in an algorithm that achieves its goal of small circumbinary planet detection, our methods could no doubt be more rigorously optimized. We used the biweight and cosine filters from Wotan since they are amongst the best for single stars, but it might be worth testing different filters with different window lengths. More specifically, one could run the detection algorithm from start to finish with modified detrending and quantify planet detection efficiently, similar to the tests run in Wotan ). The multi-dimensional nature of both the eclipsing binary detrending and the planet detection would make this computationally challenging though.

Transits on both primary and secondary stars
Our algorithm as it stands only detects transits on the primary star. In all but one of the circumbinary systems the planet exhibits multiple transits on the primary star and those alone would be enough for confirmation 20 In fact, most of the known circumbinary systems have secondary transits that are either invisibly shallow or geometrically non-existent (Martin 2019). Accounting for both transits could improve the SNR for near equal mass binaries such as Kepler-34, but for most binaries it would likely be folding noise onto the signal. Given the current method of timing transits by iterating the N-body integrator back and forward, timing secondary transits would also slow down the detection algorithm. A future improvement to this code will be to at least give the user the option of folding both primary and secondary transits, preferably without a significant slow down to the algorithm. The choice of whether or not to include secondary transits may be informed by whether or not the secondary eclipse depths are comparable to that of the primary. 5.5. Future applications 5.5.1. Systematic survey for small circumbinary planets in Kepler using the Windemuth catalog Paper II in preparation is an application of this new Stanley algorithm to the over 700 eclipsing binaries with stellar masses and radii derived in the Windemuth catalog (Windemuth et al. 2019a). It is both a search for new planets and a new, tighter constraint on the circumbinary planet occurrence rate.

Application to the entire Kepler eclipsing binary catalog
The entire Villanova EB catalog totals almost 3000 targets. Whilst some of these will likely remain inappropriate for planet searches, such as heartbeat stars, ellipsoidal variables and very tight contact binaries, we hope to expand our sample beyond the initial ∼ 700 Windemuth targets. Those not in the Windemuth catalog typically do not have accurately derived masses and radii though. As demonstrated in this paper though, circumbinary transits can be found even if the stellar parameters are a little bit (∼ 10 − 20%) wrong. Furthermore, the algorithm can scan over a grid of stellar parameters too. There remains hope then that the algorithm can be applied to find circumbinary transits even when the binaries are poorly constrained.

Application to TESS
The short, typically one month observing windows of TESS are challenging for the discovery of circumbinary planets, which to date have been found on orbits of 50 days and typically much longer. Some TESS targets closer to the ecliptic pole receive longer windows, but none to the same extent as Kepler's four years. To discover planets with longer periods than the observational window we may exploit the concept of a "1-2 punch", where a planet on a single conjunction transits both primary and secondary stars, which can be exploited to better characterise the planet than one of the same period around a single star . Adapting this algorithm to discovery planets with less transits, using an instrument with more variability and a smaller aperture, is a challenging future task.

Public release
The Stanley code is in the process of being finalised for public release. This involves making the code cleaner and well-commented, adding in depth documentation and testing it on a variety of systems. Please contact the authors for early access.

CONCLUSION
We have presented Stanley: a new automated algorithm to detrend eclipsing binaries and then find circumbinary planets. By phase-folding on variable transit times and durations, which are physically modelled with an N-body code, the algorithm coherently stacks transits. This allows the detection of small planets with shallow transits which may be individually indistinguishable above the noise.
We significantly recover all known circumbinary systems, including all three planets in the Kepler-47 system and the 1108 day Kepler-1647b (but only if the detection criteria are loosened). For each system we scale down the transit depths and can find significantly smaller planets than those found by eye, with the best being Kepler-16, detectable to a sub-Earth radius. We show also search for new companion planets but find none, which means that they either do not exist or are significantly smaller than the known gas giants.
This paper is a precursor to Paper II, where we apply this algorithm to over 700 Kepler eclipsing binaries to hunt for new circumbinary planets. We will also derive improved occurrence rate constraints on the circumbinary planet population, which can then be compared with that around single stars.