CETRA: a fast, sensitive exoplanet transit detection algorithm implemented for GPUs

ABSTRACT

We present the Cambridge Exoplanet Transit Recovery Algorithm (cetra), a fast and sensitive transit detection algorithm, optimized for GPUs. cetra separates the task into a search for transit signals across linear time space, followed by a phase-folding of the former to enable a periodic signal search, using a physically motivated transit model to improve detection sensitivity. It outperforms traditional methods like Box Least Squares and Transit Least Squares in both sensitivity and speed. Tests on synthetic light curves demonstrate that cetra can identify at least 20 per cent more low-SNR transits than Transit Least Squares in the same data, particularly those of long period planets. It is also shown to be up to a few orders of magnitude faster for high cadence light curves, enabling rapid large-scale searches. Through application of cetra to Transiting Exoplanet Survey Satellite short cadence data, we recover the three planets in the HD 101581 system with improved significance. In particular, the transit signal of the previously unvalidated planet TOI-6276.03 is enhanced from |${\rm SNR}=7.9$| to |${\rm SNR}=16.0$|⁠, which means it may now meet the criteria for statistical validation. cetra’s speed and sensitivity make it well-suited for current and future exoplanet surveys, particularly in the search for Earth analogues. Our implementation of this algorithm uses NVIDIA’s CUDA platform and requires an NVIDIA GPU, it is open-source and available from GitHub and PyPI.

1 INTRODUCTION

Three quarters of the current list of several thousand confirmed exoplanets were discovered through the observation of a reduction in the apparent flux of their host star as it is transited by the planet. The retired Kepler space telescope is responsible for the majority of these discoveries, through the Kepler and K2 missions (Borucki et al. 2010; Howell et al. 2014). The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) is responsible for several hundred of them, but it has generated many times more candidates that are awaiting independent confirmation. The upcoming PLAnetary Transits and Oscillations of stars (PLATO; Rauer et al. 2024) mission, with a planned 2026 launch, is designed specifically to detect and characterize terrestrial planets around Sun-like stars.

The depth of an exoplanet transit is approximately proportional to the square of the ratio of the planet and star radii. The fractional reduction in flux from a Sun-like star due to the transit of an Earth-like planet is of the order of |$10^{-4}$|⁠, or 100 parts-per-million (ppm). This level of photometric precision is essentially unachievable from the ground, but measurements made by space-based instruments highlight a further problem: the amplitude of the intrinsic variability of the star is typically many times larger than such a transit signal (McQuillan, Aigrain & Roberts 2012). Many large-scale exoplanet transit detection pipelines rely on a pre-whitening filter, i.e. they attempting to remove correlated noise (intrinsic stellar variability but often also instrumental artefacts) such that only the transit signal and uncorrelated (white) noise remain, before running a transit detector. For targets whose characteristic variability is particularly challenging for transit detection (amplitude and/or frequency-wise) it can be beneficial to attempt a simultaneous fit of the correlated noise (Garcia et al. 2024). However, this approach is relatively expensive computationally, and doubt has been cast on the need for it in most situations (Kovács, Hartman & Bakos 2016).

The standard transit detection algorithm has for many years been the Box-fitting Least-Squares algorithm (BLS) of Kovács, Zucker & Mazeh (2002). This approach phase-folds the light curve over a grid of periods and finds the least-squares box-shaped signal for each using a moving window approach. Hippke & Heller (2019) demonstrated that the sensitivity of this type of detector could be improved upon by using a transit-shaped profile, which is particularly important for low signal-to-noise ratio (SNR) transits, and they presented their Transit Least Squares (TLS) algorithm. Even modest improvements to sensitivity at low SNR are valuable since this is where prized Earth–Sun analogues typically reside. Despite the additional computational complexity over BLS, TLS is not significantly more expensive as a result of optimizations made to the period and duration search grids.

Another design choice is whether to search for periodic signals in phase-folded data, or to search for signals in linear time data and phase-fold the results. Both BLS and TLS use the former approach, which requires an expensive sorting step in phase space for each trial period. If many light curves with the same observation times are available, one might perform sorting once per period for all and share the sorted indices. This optimization is used by e.g. Aigrain & Irwin (2004). The linear-then-periodic approach offers reduced computational expense but tends to be used primarily in algorithms that incorporate simultaneous detrending (e.g. Berta et al. 2012; Foreman-Mackey et al. 2015; Garcia et al. 2024; with the notable exception of the one used in the Kepler pipeline and described by Jenkins et al. 2010). A preliminary search in linear space also enables a search in these results for single transits (a.k.a. monotransits), which are early indicators of the presence of long period planets.

The increasing volumes of data being supplied by current and future space-based instruments encourage the use of highly efficient transit search algorithms. Graphics processing units (GPUs) have potential here, since detection algorithms typically comprise highly parallel operations, operations which GPUs are ideally suited to perform.

We present a new transit search algorithm, the Cambridge Exoplanet Transit Recovery Algorithm (cetra). It performs the majority of the processing on NVIDIA GPUs using the CUDA parallel computing platform, and hence requires a compatible device. cetra incorporates the ideas behind several existing algorithms, while also including additional optimizations in order to try to improve on their sensitivity and/or compute efficiency. It is designed to be run on detrended light curves, uses an exoplanet transit shaped model, and is a two-stage linear-then-periodic algorithm. In this article, we detail the design and implementation of cetra (Section 2), and evaluate its computational performance and detection sensitivity relative to established algorithms (Section 3). We then reanalyse the HD 101581 multiplanet system in the TESS high cadence data, resulting in an increased transit signal detection significance relative to that of the existing literature (Section 4). cetra is published on github¹ and PyPI under an open-source MIT license.

2 ALGORITHM DESCRIPTION

cetra searches light curves for transit-like features using a true transit model, since this offers improved sensitivity relative to a box model, but with detection separated into linear and then periodic signal searches. The linear search component finds the maximum-likelihood model depths over a 2D grid, with trial durations on one dimension and an array of |$t_0$|s that is linear in time on the other. The periodic search component then phase-folds the output from the linear search over a grid of periods and returns the maximum joint-likelihood for each. cetra assumes that light curves it is given have already been detrended. More detail about each component, and also about some important setup stages, is provided below.

Internally, cetra works in units of normalized flux offset from baseline, i.e. baseline flux is ‘0.0’ and if flux from the star were blocked completely then its normalized flux offset would be ‘1.0’. Time units are days, and where appropriate, those relative to the light curve start time are used. These conventions mean that the use of single precision floating point numbers introduces smaller errors than otherwise. In the case of the flux for example, using a baseline of ‘1.0’ introduces a constant single precision floating point error of the order of |$0.1~{\rm ppm}$|⁠, where using ‘0.0’ incurs an error that is proportional to the flux deviation and only reaches that level when no flux is observed. For times, using this convention incurs a floating point error of only a few seconds over light curves that are several years long. These are acceptable for the purposes of transit detection, and the speed difference between single and double precision operations can be significant on some GPU models. However, these conventions are broadly transparent to the user.

2.1 Setup

2.1.1 Light-curve resampling

The first setup stage is a resampling of the input light curve into a regular, ordered time grid. For example, PLATO Level 1 light curve points are averages between camera groups (Rauer et al. 2024), which leads to an irregular cadence. Even a strictly regular on-board cadence will yield an irregular cadence when converted to Barycentric dynamical time. The main advantage of the resampling is that the resampled light-curve array element corresponding to any point in time can be determined with ease. For example, in the context of an exoplanet transit, the first and last resampled light-curve element indices of the window covering the transit are simply |$\lceil {}t/{\rm cadence}\rceil {}$| and |$\lfloor {}(t + {\rm duration})/{\rm cadence}\rfloor {}$|⁠, respectively, where t is the start time of the transit relative to the start time of the resampled light curve, duration is the duration of the transit, and cadence is the cadence of the resampled light curve. Resampling in this way avoids any need for repeated, expensive sorting of the light curve, and makes the algorithm robust against gaps and irregular cadence in the input data.

The frequency of points on the resampled grid is a choice for the user. Observed fluxes are simply assigned to the nearest output point. Where multiple input observations contribute to single output points they are combined using an inverse variance weighted mean and their uncertainties are propagated formally. Resampled grid points that have no associated observational data are assigned null fluxes with infinite error. If no resample cadence is supplied by the user then cetra will use the median cadence of the input data by default.

It is generally preferable to resample to a higher frequency than the input observations, since this limits morphological distortions of the light curve. Kipping (2010) demonstrates the potential impact of distortions in the context of binning. In terms of computational expense, the frequency of the resampled light curve only impacts the linear component of the transit search, and we demonstrate in Section 3.2 that a relatively high frequency is not particularly costly.

An additional detail is that the resampled light curve is also padded with null data at the beginning and end, with data length equal to half the maximum trial duration (see Section 2.1.3). This is a straightforward way to ensure that the linear search can detect transits that straddle the beginning and end of the data.

2.1.2 Transit model

cetra includes densely sampled transit light curves produced using the batman package² (Kreidberg 2015) for three different values of impact parameter. The provided models use the following parameters:

• Planet radius |$R_p=0.03R_*$|⁠,

• impact parameter |$b=[0.32, 0.93, 0.99]$|⁠,

• and quadratic limb darkening parameters |$u_1,u_2=0.4804,0.1867$|⁠.

These are also the default model parameters used by TLS (in the |$b=0.32$| case). However, users might wish to supply their own models instead, since the default parameters are not ideal for all targets. For example, if the limb darkening parameters of the target are known, even approximately e.g. from atmospheric modelling, then a model that uses these might be supplied instead and may improve detection sensitivity. Additionally, since changes to the impact parameter can make a significant difference to the shape of the model, it may be worth running cetra across a sequence of values. An optimal set of parameters is likely to depend on survey design and sensitivity requirements, and is outside the scope of this paper.

The chosen/supplied model is ultimately resampled to a user-defined number of samples using simple linear interpolation. This is done to ensure that the model will fit into shared memory on the GPU alongside some additional arrays. The linear search kernel³ makes frequent comparisons of the light curve to the model, and hence the much lower latency of shared memory makes its use worthwhile.

Within the linear search kernel the nearest model point is selected for a given light-curve point, transit reference time, and duration being considered. This nearest neighbour selection introduces a small error in flux, which is greatest where the flux gradient is largest. The maximum and mean error introduced versus the number of samples chosen is shown in Fig. 1 for the default transit model. The maximum error occurs during ingress and egress, but the mean error is around an order of magnitude smaller. Since the volume of shared memory available is limited, and memory copies should be be reduced for efficiency where practical, it is advisable to use as few model samples as necessary to achieve an acceptable level of error. For reference, the linear search requires 256 bytes of shared memory space in addition to the transit model. The default number of samples is 1024, for which the maximum error is |$\approx {}1$| per cent and the mean error is |$\approx {}0.1$| per cent.

It is worth noting that while cetra was written with exoplanet transit searches in mind, any model can be supplied to it (e.g. exocomets, exoplanets with rings, etc.). In principle it might be used for any 1D template matching purpose, unrelated to astronomical light curves or even to astronomy, though some modification might be necessary. No assumption is made about whether the deviation of the model is positive or negative, it does not enforce a minimum depth, so e.g. stellar flares and other such flux-positive transient phenomena might also be searched for instead. For alternative use cases consideration must be made as to whether the default duration (and period) grids are appropriate (see Section 2.1.3).

2.1.3 Duration, t|$_0$| and period grids

The linear transit search is performed over grids of transit reference times, |$\lbrace t_{i}\rbrace _{i=1}^{I}$| with I times, which are the mid-point of the transit; and durations, |$\lbrace D_{j}\rbrace _{j=1}^{J}$| with J durations. Subsequent periodic signal searches are performed over a grid of periods, |$\lbrace P_{k}\rbrace _{k=1}^{K}$| with K periods.

By default, cetra uses the same TLS specification for the duration and period grids, although the user is free to supply their own of either. For a given period, durations that fall outside those expected for typical Keplerian orbits (see fig. 5 of Hippke & Heller 2019 in particular) are not evaluated, though the user can override this behaviour. A major optimization of TLS is the use of duration and period grids that are astrophysically motivated and optimally sampled (Ofir 2014; Hippke & Heller 2019). We saw no reason to deviate from these grid definitions.

The cetra reference time grid uses a single step size for all durations. The default is 1 per cent of the minimum duration, although the user is free to specify an alternative fraction of the minimum duration. TLS uses a transit reference time grid that is a function of the trial duration, which means that there are fewer trial reference times for longer durations. We might adopt this optimization for cetra. However, the fixed time-step simplifies the algorithm and its output, and means that the results of the linear search could be directly probed for transit timing variations (TTVs; adopting the approaches of e.g. Carter & Agol 2013; Leleu et al. 2021) without a reduction in sensitivity for long-duration transits. The density of the reference time grid impacts the compute time of the linear search only, and as is demonstrated in Section 3.2, increasing the run time of the linear search is not a major burden.

2.2 Linear search

The linear search traverses the grid of transit reference time, |$\lbrace t_{i}\rbrace _{i=1}^{I}$|⁠, and transit duration, |$\lbrace D_{j}\rbrace _{j=1}^{J}$|⁠, element pairs. Each pair defines a region in time in which the light curve should be compared to the transit model, this time span will hereafter be referred to as the transit window. Because cetra resampled the light curve during setup, it can trivially determine the section of data that covers the transit window, see Section 2.1.1. This means that expensive whole-array traversals are not required, and only two traversals of the subset of the light-curve array corresponding to the transit window are needed. The first pass computes the maximum-likelihood transit depth, by which the transit model fluxes must be multiplied in order to optimally match the observed fluxes. The second pass computes the log-likelihood of the observed fluxes given the scaled model fluxes, and subtracts from this the log-likelihood of a constant flux model to obtain the likelihood ratio between the two.

With a transit model normalized to maximum depth |$=1$|⁠, obtaining an equivalent model with a transit depth of |$\Delta {}$| is simply a matter of multiplying each element by a scale factor of |$\Delta {}$|⁠. In a given transit window, the model scale factor (i.e. transit depth) implied by each flux element, |$f_n$|⁠, is |$\Delta {}_n = f_n/m_n$|⁠, where |$m_n$| is the appropriate model element obtained using nearest neighbour interpolation (see Section 2.1.2). The uncertainty on the model scale factor is |$\sigma {}_{\Delta {}_n} = \sigma {}_{n}/m_n$|⁠, where |$\sigma {}_{n}$| is the uncertainty on the flux. The maximum-likelihood scaling factor, |$\Delta {}_{i,j}$|⁠, over the window is the inverse variance weighted mean of all values of |$\Delta {}_n$| inside the window, and the variance on this, |$\sigma _{\Delta {}_{i,j}}^2$|⁠, is the reciprocal of the sum of their weights. These two parameters give us the maximum-likelihood transit depth within the window and its error, assuming the model is appropriate. However, determination of the likelihood of the model requires a second traversal.

Recall that cetra works in normalized flux offsets from baseline internally, and that input light curves will have been detrended. As a result, out-of-transit model fluxes are zero.

For a vector of fluxes |$\mathbf {f}$|⁠, i.e. a light curve of length N, the log-likelihood of the |$t_{i},D_{j}$| model is:

Computing the result of equation (1) explicitly would require traversal of the entire light curve. However, the contribution of all out-of-transit points does not change between a transiting planet model and a constant flux model (i.e. a model containing no transit). If instead of equation (1) we compute the likelihood ratio (⁠|$LR$|⁠) between these two models we only need to consider points within the transit window. As an added bonus, the second term inside the sum also cancels inside the transit window, so the |$\log$| operation is avoided. The likelihood ratio, given in equation (2), can therefore be computed at a much reduced cost relative to the full log-likelihood.

The output of the linear search is a set of 2D arrays of: Depth, |$\lbrace \Delta {}_{i,j}\rbrace _{i,j=1}^{I,J}$|⁠; depth variance, |$\lbrace \sigma _{\Delta _{i,j}}^2\rbrace _{i,j=1}^{I,J}$|⁠; and likelihood ratio, |$\lbrace LR_{i,j}\rbrace _{i,j=1}^{I,J}$|⁠. With |$\lbrace t_{i}\rbrace _{i=1}^{I}$| and |$\lbrace D_{j}\rbrace _{j=1}^{J}$| along the two primary array axes.

Fig. 2 shows some output from the linear search for an example light curve containing transits. The example light curve is astrophysically implausible, but suitable for demonstration purposes. The input transits have 5-d durations and 1000 ppm depth, and occur at |$t={10,30,50}$|⁠, i.e. they correspond to a signal with a 20-d period. The input model parameters are correctly recovered.

Since all resampled light-curve array elements with no contributing observations have infinite uncertainties, they have no influence on either the depth or the likelihood ratio. These may include: null points padded to the beginning and end of the resampled light curve, which enable cetra to detect transits that straddle the beginning or end of the data without any special treatment within the linear search kernel; points within gaps in the observations, e.g. between observing cycles or during daytime (if the data are ground-based); and points between observations if the resampling frequency is higher than that of the observations. In some circumstances there may be no observations within the transit window, e.g. data gaps longer than any of the trial durations, or where |$t{<}D/2$| due to the aforementioned padding. Output arrays elements are null in such circumstances, the latter case can be seen in panels (b) through (d) in Fig. 2.

If required, the log-likelihood of the whole light curve for all transit models (i.e. the result of equation 1) can be calculated from the array of likelihood ratios described above. Simply add to them the log-likelihood of the light curve for a constant flux model (i.e. equation 1 evaluated for |$\Delta {}_{i,j}=0$|⁠), which is independent of t and D and hence need only be computed once for all pairs, at negligible expense.

2.3 Periodic search

The periodic search proceeds as per Foreman-Mackey et al. (2015), albeit without allowing variation in the depths of transits contributing to a periodic signal.

In principle one might expect some low level of transit depth variation, in practice however, depth variation is likely to be negligible for the purposes of transit detection. Hence, cetra assumes that all transits contributing to a periodic signal have a common depth. If transit depths do vary, this will result in a reduction in the likelihood of the periodic signal. The likelihood function for the depth at each |$t_i,D_j$| point is Gaussian, with depth |$\Delta {}_{i,j}$|⁠, and variance |$\sigma ^2_{\Delta {}_{i,j}}$|⁠. The log-likelihood of the data for a model with some alternative depth, Z, is therefore:

And it follows that the likelihood ratio for this alternative model is:

For any sequence of non-overlapping transits, the corresponding set of depths, depth variances, and likelihood ratios can be obtained from the results of the linear search. Nearest neighbour interpolation can be used to obtain approximates of these values where a transit falls between grid elements, an approach also adopted by Foreman-Mackey et al. (2015). From these we can compute their maximum likelihood common depth and its variance, |$Z_{i,j,k}$| and |$\sigma {}^2_{Z_{i,j,k}}$| hereafter, which are their inverse variance weighted mean depth and the reciprocal of the sum of their weights.

The sum of the results of equation (4) for all transits contributing to a given periodic signal, having adopted |$Z_{i,j,k}$|⁠, gives the likelihood ratio of the data given that periodic model versus a constant flux model (⁠|$LR_{i,j,k}$|⁠). This follows logically, it is simply the log-likelihood of the data in all the transit windows for a (now multiple) transit model minus the same for a constant flux model, which itself is the likelihood ratio of the periodic signal model versus the constant flux model.

This method assumes that the linear search result contribution of each transit is independent of all others. Foreman-Mackey et al. (2015) and Garcia et al. (2024) note that this is only approximately true in their case, since they were fitting systematic trends in addition to transits. Their changing belief about the transits also changes their belief about the systematics model, thereby influencing all other transits. Since cetra requires a pre-whitened light curve, the assumption might be closer to reality, though due to the intricacies of detrending it is still unlikely to be strictly true. One point of note is that should the transit windows overlap then the results will be invalid, for this reason cetra requires that trial durations are smaller than the trial period.

For each trial period the 2D array results that are the outputs of the linear search are essentially phase-folded, but with the additional common depth stipulation. For each, a smaller set of 2D arrays of depth, depth variance, and likelihood ratios is produced. Copying each of these from the GPU back to the host machine would be a major bottleneck in the algorithm. However, the user is unlikely to be interested in anything other than the model with the maximum likelihood. Consequently, the periodic search kernel incorporates a reduction operation to identify the maximum-likelihood element in the above arrays, and they return only this and the corresponding model parameters to the host machine.

The outcome of the periodic search is an array of periods with a corresponding array of their maximum-likelihood ratios (demonstrated in the upper panel of Fig. 3), and their transit reference times, durations, depths, and depth variances (demonstrated as a signal-to-noise ratio in the lower panel of Fig. 3). Should an array of log-likelihoods be required instead of likelihood ratios, then adding to them the log-likelihood of the constant flux model will achieve this. The log-likelihood of the constant flux model is itself a constant for a given detrended light curve, and this approach is common to the results of the linear search.

3 ALGORITHM PERFORMANCE

3.1 Simulated transit recovery

To evaluate the period recovery performance of cetra versus TLS, we generated 20 000 synthetic transit light curves, scattered them using three different white noise levels, and checked to see which algorithm recovered the input periods. We did not include red noise in the light curves as the objective was to test detection algorithm performance, not detrending performance. Synthetic light curves were 173 d in length, sampled at 600 s cadence. Since real observational data regularly includes breaks in observations, particularly for long epoch baselines, we included a 7 d gap in the synthetic light curves starting at day 90. The three levels of white noise that we tested were 34, 160, and 800 ppm h|$^{-1}$|⁠.

The injected transits had periods drawn randomly from a uniform distribution between 1.0 and 86.5 d, the upper bound meaning that at least two transits would always be present in the data unless one fell in the data gap. A transit duration was then drawn randomly from a uniform distribution with upper and lower bounds appropriate for the selected period (see Section 2.1.3, above; and Hippke & Heller 2019 fig. 5). A transit depth was randomly selected from a log-uniform distribution bounded at 50 and 1 000 ppm. The injected transit model had the same impact parameter and limb darkening profile that TLS and cetra use by default.

We considered an input period to have been recovered if the output period was within 1 per cent of the input period, or 1 per cent of half or double the input period. An important detail is that both algorithms were always given exactly the same light curves, i.e. the parameters and noise properties were identical. The counts and percentages of recovered periods for each noise level and algorithm are shown in Table 1, though the exact values of these are highly test-specific. Inferences about planet recovery rates from real surveys should not be drawn, since these are dependent on e.g. real transit parameter distributions, detrending performance, etc. As a consequence of this, and that these light curves contain only white noise and transits, and that detection thresholds are a choice for the user, we did not consider the level of contamination from false positives. The key takeaway is how the algorithms perform relative to each other.

In this test cetra always outperforms TLS, but at low signal-to-noise ratios the improvement is more marked. For the test sample as a whole, in the 34 and 160 ppm h|$^{-1}$| cases it correctly recovers up to 5 per cent more periods, but in the 800 ppm h|$^{-1}$| case it recovers 21 per cent more. If we consider only the subset containing Earth–Sun analogues with transit depth of 75 to 125 ppm, and with two or three observed transits,⁴ then cetra recovers |$1.3\times {}$| and |$1.7\times {}$| as many periods as TLS for the 160 and 800 ppm h|$^{-1}$| samples, respectively. Earth–Sun analogues are highly prized, so improvements to the detection rate of such systems are valuable.

We attribute the improvement in sensitivity to the combination of a few factors:

First, Hippke & Heller (2019) note (in their section 4.2) that TLS assumes a constant cadence in phase space. While the assumption is true on average, detection of longer period signals will be hindered to a greater degree. A significant cause of non-uniformity in phase-space cadence are gaps in the observational data. Gaps are almost always present in long epoch baseline light curves, as indeed they are in the above test. cetra does not assume a constant cadence in phase (or time) space.

Secondly, TLS sweeps the light curve in phase space by the larger of one data sample or one per cent (by default) of the trial duration. This means that the true TLS |$t_0$| grid spacing may be larger than expected for long period planets. For a two-transit periodic signal in a light curve, the effective |$t_0$| grid spacing is half the observing cadence. This is potentially quite significant for long cadence (e.g. Kepler) data. cetra always uses one per cent (by default) of the minimum trial duration, it does not increase the effective grid spacing due to algorithmic constraints.

Finally, though this does not impact the above test as the flux uncertainties were uniform, TLS does not consider variable flux uncertainties when calculating the maximum-likelihood depth of the transit model in a given window. TLS uses the arithmetic mean⁵ to determine the best model depth in the transit window, where cetra uses the inverse variance weighted mean.

3.2 Compute speed

We evaluated compute performance through comparison to BLS (via the astropy timeseries module which uses a C implementation via Cython) and TLS (which is Numba-compiled python code) for a range of light-curve cadences and durations. The results are shown in Fig. 4. Tested cadences were chosen to match existing and future data from the Kepler, TESS, and PLATO missions. The 25 s cadence matches the PLATO short-cadence targets, and is comparable to the 20 s cadence of some TESS second extended mission targets. The 120 s cadence matches the majority of TESS targets. The 600 s cadence matches the PLATO long-cadence targets, and the 1800 s cadence matches that of Kepler long-cadence and K2 light curves.

The measured TLS and BLS compute times are for a single CPU core. We have typically found that performance per thread is reduced somewhat for TLS when using multiple threads internally, and the astropy implementation of BLS does not offer multithreading support. Parallelization can also be achieved through running multiple separate, single-thread instances, each processing a different light curve or with different search parameters. Many system-architecture specific factors come into play when determining the optimal parallelization method. To an approximate first order, a linear scaling with CPU core count can be applied, and consequently the compute times shown in Fig. 4 can be scaled to account for the user’s own CPU core counts.

Computing overheads (resampling, grid generation, data transfers to/from GPU, model refinement, etc) were included for the cetra, TLS, and BLS measurements where applicable. The cetra time includes both components of the search – linear and periodic. BLS was run over the CETRA (/TLS) period and duration grids, skipping periods shorter than five times the duration each time. All three algorithms were run using their default parameters.

The compute times shown in Fig. 4 were measured using a machine comprised of a pair of Intel Xeon Platinum 8358 CPUs, and a pair of NVIDIA A100 GPUs, although BLS and TLS were run single-threaded and cetra used only one GPU. On a NVIDIA RTX A5000 GPU, which is very roughly one quarter of the monetary cost of an A100 at time of writing, cetra was only around one fifth slower.

Fig. 4 demonstrates that cetra run times are largely independent of observing cadence, at least within the range of cadences that are likely to be useful for transit detection purposes. The time to complete the linear search is roughly a linear function of the number of light- curve points, which is a function of light-curve cadence and length.⁶ The periodic search is only a function of light-curve length, since this increases the highest checked period (assuming a requirement of 2 or more visible transits), and the number of transits that contribute to a signal at a given period. Importantly, the linear search is anywhere between a few times faster and a few orders of magnitude faster than the periodic search. For example, for the 25 s cadence, 30 day light curve the linear search took 0.26 s, while the periodic search took 0.94 s; and for the 1800 s cadence, 720-d light curve the linear search took 0.17 s, while the periodic search took 177.92 s. The periodic search times were essentially constant for a given light-curve epoch baseline. Hence the cadence does not significantly impact the compute time. For this reason, resampling at a high frequency (see Section 2.1.1) has a low cost.

4 HD 101 581 SYSTEM RECOVERY

No demonstration of a transit detection algorithm would be complete without a test on real data. To this end, we attempted a re-detection of the HD 101581 (GJ 435, TOI-6276, TIC 397362481) multiplanet system using cetra. Kunimoto et al. (2025) validated HD 101581 b and c, a pair of Earth-sized planets found orbiting a K5V host. They also identify a third potential Earth-sized planet, TOI-6276.03, but find the |$7.9\sigma {}$| signal too weak to statistically validate. They used the TESS 120 s cadence Presearch Data Conditioning Simple Aperture Photometry (PDCSAP; Smith et al. 2012; Stumpe et al. 2012, 2014) light curves from sectors 63 and 64, with the biweight detrending algorithm implemented in the wotan Python package (Hippke et al. 2019), and TLS. We used the 20 s cadence PDCSAP light curve, with a Notch-like (Rizzuto et al. 2017) detrending algorithm (details below), and cetra.

The following analysis is available as a Jupyter Notebook in the examples directory of the cetra GitHub repository.

4.1 Detrending

Capitalizing on the speed at which cetra is capable of checking for transit-like signals in linear light-curve data, we developed a detrending algorithm inspired by the Notch filter algorithm of Rizzuto et al. (2017). While not intended to be considered part of cetra, it is included in the code base for completeness. The algorithm traverses the light curve in linear space, searching for windows in which a transit-plus-quadratic model fits the observed flux better than a pure quadratic model based on some information criteria (Bayesian and Akaike are included). It searches over grids of |$t_0$| and duration in exactly the same manner as the linear search described in Section 2.2. With the above information, a basic model of the light curve is produced, which is then used to account for potential transits while refitting the baseline flux of the light curve. The algorithm also uses the non-transit regions of the model to measure a scaling factor to be applied to the flux uncertainties such that the standardized residuals resemble a unit Gaussian. In addition to those of the linear search (e.g. |$t_0$| and duration grids), the main tuning parameters of the algorithm are the widths of the detection and flux fitting kernel windows, the type of information criterion used, and its difference threshold.

We ran the above algorithm on the 20 s cadence HD 101581 light curves from the two TESS sectors separately. We used transit detection and flux fitting kernel window widths of 1.5 and 1.0 d, respectively, and the Bayesian Information Criterion (BIC) with a difference threshold of 10. The input light curve, potential transit locations, fitted trend, and detrended light curve are shown in Fig. 5. Detrending was performed using an NVIDIA RTX A5000 GPU and took 1 to 2 s per sector.

4.2 Transit detection

We ran the detrended HD 101581 light curve through cetra with a maximum transit duration of 0.5 d, but otherwise used the default search parameters. We ran transit detection on the light curve three times, for each run after the first we nullified the flux within and around the highest likelihood transit signal from the previous search. For illustration purposes, the array output from the first run is shown in Fig. 6, and the periodogram from the final run is shown in Fig. 7.

Details of the resultant periodic signals are given in Table 2 in the order of their detection. Fig. 8 shows the maximum-likelihood transit model over the data and binned data for each of the detected signals.

4.3 Discussion

A transit detection run on the above light curve with TLS v1.32 took of order 15 min on a single thread of the test workstation. Parallelizing over all six available CPU cores reduced the compute time to around 4 min. By contrast, a single run on one NVIDIA RTX A5000 GPU took 2.6 s, a factor 90 reduction in compute time versus TLS multithreaded. The obvious implication of this is that the user can process more targets in less time. It also enables the user to more cheaply sweep detrending and detection configurations in order to optimize searches for specific targets, and to run injection-recovery tests at moderate scales without the need for large cluster computing facilities.

Note that cetra did not detect TOI-6276.03 in the 120 s cadence TESS light curves any better than Kunimoto et al. (2025) was able to with TLS. When run on the same detrended 20 s cadence light curves as cetra, TLS found all three signals with signal detection efficiency of around 18 (18.5 for TOI-6276.01, 18.0 for TOI-6276.02, and 17.8 for TOI-6276.03). Possible reasons for the marginally improved sensitivity of cetra relative to TLS are given in Section 3.1, but the signal of TOI-6276.03 in the 120 s cadence light curve is not particularly affected by any of them. So cetra not being any more sensitive than TLS in this specific case is not surprising. It is not clear why the signal of TOI-6276.03 is more significant in the short cadence light curve than the long. Possibly the detrending performance is better with a higher frequency of observations, we certainly found it easier to tune the detrending parameters for this data than for the long cadence data, but this is speculation. It is clear that there may be value in running detection algorithms on short cadence data rather than binning to longer cadence. This renders the compute efficiency improvements of cetra over existing algorithms that we demonstrated in Section 3.2 more important.

Kunimoto et al. (2025) note that the false positive probability of TOI-6276.03 is low enough to satisfy the criteria for statistical validation when boosted due to its membership in a multiplanet system, though they did not consider it a statistically valid planet due to the low |${\rm SNR}=7.9$| in their analysis. We suspect the increase in the detection significance to |${\rm SNR}=16.0$| presented in this analysis would justify reconsideration of the statistical validity of the candidate.

5 SUMMARY

Motivated by the requirement to detect exoplanet transits in ever-increasing volumes of high cadence stellar light curves, and the apparent stalling of the trend toward increasing CPU clock speeds observed in the last few decades (Leiserson et al. 2020), we present a new transit detection algorithm that takes advantage of GPUs and the highly parallel nature of the task. cetra has been implemented with CUDA for NVIDIA GPUs, although we’re exploring the possibility of porting to other frameworks to enable processing on devices from other manufacturers. After pre-processing, it first evaluates the likelihood of a transit signal across grids of |$t_0$| and transit duration in linear, i.e. time space. It then phase-folds the results of the linear search to generate a periodogram, from which periodic signals can be identified.

In common with the well established TLS code, cetra uses a true transit (or other) model in order to maximize detection sensitivity. When tested against 20 000 synthetic light curves with varying noise levels, it was found to consistently outperform TLS, with a more marked improvement seen for low SNR planets, particularly those with long periods. The improvement is largely attributable to a more robust handling of gaps in the observational data, and overcoming of a |$t_0$| sampling frequency limitation that specifically impacts detection of long period planets.

In addition to the improved sensitivity, through benchmarking against the most commonly used existing codes, BLS and TLS, cetra was demonstrated to be dramatically faster to run. The greatest improvement in compute speed of a few orders of magnitude is seen for high cadence light curves. Should the user wish to search only for monotransits, then cetra can do this in under a second even for high cadence light curves that are several years in length. Long-period planets appear as monotransits until sufficient epoch baseline and phase coverage is achieved by the observation campaign. A detection algorithm that is sensitive to them means decisions regarding key follow up observations can be made sooner.

cetra detected the three Earth-sized planets in the HD 101581 (TOI-6276) system within the TESS high cadence light curves. Two of the planets were previously detected with TLS and validated in the lower cadence light curves by Kunimoto et al. (2025), but the third was found to be too low SNR to statistically validate. We found generally improved signal-to-noise ratios, with the signal of the previously unvalidated planet (TOI-6276.03) boosted from |${\rm SNR}=7.9$| to |${\rm SNR}=16.0$|⁠, which is potentially sufficient to reconsider the statistical validity of the planet.

The ability to rapidly search large data sets makes cetra ideally suited for current and future space-based exoplanet hunting missions (e.g. TESS, PLATO). The code is available from GitHub and PyPI under an open-source MIT license.

ACKNOWLEDGEMENTS

We are grateful to the referee, René Heller, for the careful reading and helpful comments. We acknowledge the support of the PLATO Mission Consortium, and from UKRI-STFC grants ST/X001628/1 and ST/X001571/1.

DATA AVAILABILITY

TESS light curves of HD 101581 are publicly available online. Code and analyses are available from the cetra GitHub repository at https://github.com/leigh2/cetra.

Footnotes

REFERENCES