nuance: Efficient Detection of Planets Transiting Active Stars

The detection of planetary transits in the light curves of active stars, featuring correlated noise in the form of stellar variability, remains a challenge. Depending on the noise characteristics, we show that the traditional technique that consists of detrending a light curve before searching for transits alters their signal-to-noise ratio and hinders our capability to discover exoplanets transiting rapidly rotating active stars. We present nuance, an algorithm to search for transits in light curves while simultaneously accounting for the presence of correlated noise, such as stellar variability and instrumental signals. We assess the performance of nuance on simulated light curves as well as on the Transiting Exoplanet Survey Satellite light curves of 438 rapidly rotating M dwarfs. For each data set, we compare our method to five commonly used detrending techniques followed by a search with the Box-Least-Squares algorithm. Overall, we demonstrate that nuance is the most performant method in 93% of cases, leading to both the highest number of true positives and the lowest number of false-positive detections. Although simultaneously searching for transits while modeling correlated noise is expected to be computationally expensive, we make our algorithm tractable and available as the JAX-powered Python package nuance, allowing its use on distributed environments and GPU devices. Finally, we explore the prospects offered by the nuance formalism and its use to advance our knowledge of planetary systems around active stars, both using space-based surveys and sparse ground-based observations.


INTRODUCTION
Transiting exoplanets are keystone objects for the field of exoplanetary science, but detecting transits in light curves featuring stellar variability and instrumental signals remains a challenge (e.g.Pont et al. 2006, Howell et al. 2016 or Clarice Yaptangco et al. 2024).For this reason, known transiting exoplanets tend to be found around quieter stars, or belong to the population of close-in giants whose transit signals dominate over stellar rotational variability (Simpson et al. 2023).However, transiting exoplanets around active stars are discoveries with significant scientific value.First, as younger stars are more active (Skumanich 1972), being able to detect planets transiting active stars will favor the discovery of young planetary systems (e.g.Newton et al. 2022).Second, as stellar variability may originate from surface active regions (such as starpots), transiting exoplanets can be used to map the photosphere of active stars (e.g.Morris et al. 2017), benefiting both the study of stellar atmospheres and the concerning impact of their non-uniformity on planetary atmosphere retrievals (Rackham et al. 2018).Overall, enabling the detection of transits in light curves with high levels of correlated noises will greatly benefit the study of terrestrial exoplanets around late M-dwarfs, usually observed at lower SNR and more likely to display photometric variability (e.g.Murray et al. 2020 or Petrucci et al. 2024).
Commonly used transit-search algorithms, such as the Box-Least-Square algorithm (BLS, Kovács et al. 2002) are capable of detecting transits in light curves containing only transit signals and white noise.Using this method, the simplest way to detect transits in a light curve featuring correlated noise (either astrophysical or instrumental), is to first clean it from nuisance signals before performing the search.This strategy is widely adopted by the community, both using physically-motivated systematic models like Luger et al. (2016Luger et al. ( , 2018)), or filtering techniques (Jenkins et al. 2010, Hippke et al. 2019).However, when correlated noise starts resembling transits, this cleaning step, often called detrending, is believed to degrade their detectability (see subsection 4.3 of Hippke et al. 2019).In this case, the only alternative to search for transits is to perform a full-fledged modeling of the light curve, including both transits and correlated noise, and to compute the likelihood of the data to the transit model on a wide parameter space, an approach largely avoided due to its intractable nature.Nonetheless, Kovács et al. (2016) ask: Periodic transit and variability search with simultaneous systematics filtering: Is it worth it?.The latter study explores a handful of cases and generally discards the benefit of using a full-fledged approach.
In this paper, we identify regions of light curves morphological parameter space for which a full-fledged transit search is necessary, and we present nuance 1 , a method to search for transit signals while simultaneously modeling correlated noises in a tractable way.In section 1, we describe the effect of correlated noise on transit light curves and the effect of its detrending on transit signals detectability.In section 2, we present nuance, and the two main steps on which this method is based: the linear search and the periodic search.In section 3, we test the performance of nuance on a wide variety of cases, and compare our method to commonly used transit search algorithms.This include transits injected in synthetic datasets, but also in the TESS light curves of 438 rapidly-rotating M dwarfs.Finally, in section 4, we discuss the results and the limitations of nuance, before concluding in section 5.
1 Throughout the paper, nuance written in italics refers to the algorithm, while nuance in sans-serif refers to its implementation published as a Python package in its version 0.5.2 (see section 2.5).
A strong assumption when using the BLS algorithm to search for transits is that the searched dataset only contains transit signals and white noise, justifying the need for detrending.In this section, we explore this effect by simulating light curves containing a transit signal and correlated noise in the form of stellar variability, and study how detrending the latter affects the transit signal detectability depending on the light curve morphological characteristics.Hence, we first present how transit light curves are simulated, using a stochastic model of stellar variability, and describe the detection metric we employ to quantify transit detectability in the presence of correlated noise.

Light curve simulations
To perform this study, we want to simulate realistic light curves that contain instrumental signals, transit signals and correlated noise with characteristics that can be controlled using a set of interpretable parameters.To this end, we model light curves as realizations of a Gaussian process (GP; Rasmussen & Williams 2005, Aigrain & Foreman-Mackey 2023) with a mean containing the instrumental and transit signals, and a kernel allowing to model different forms of correlated noise controlled by its hyperparameters.
Let f be the simulated differential flux of a star sampled and arranged in the vector f associated to the vector of times t, such that f ∼ N (µ, Σ), i.e. that f is drawn from a GP of mean µ and covariance matrix Σ.For practical reasons, we model µ as a linear combination of M explanatory variables, such that where the first M − 1 columns of X are contemporaneous instrumental time series measurements, such as the position of the star on the detector, the sky background, co-trending basis vectors, or any other explanatory variables.The last column of X is a box-shaped transit signal with a fixed epoch, duration and period.This way, the transit signal is part of the mean linear model, meaning that once the design matrix X is constructed the transit is only parametrized by its depth ∆.
As we are interested in active stars whose fluxes feature correlated noise in the form of stellar variability, we choose a physically-motivated GP kernel to model the covariance matrix Σ, describing stellar variability through the covariance of a stochastically-driven damped harmonic oscillator (SHO, Foreman-Mackey et al. 2017;Foreman-Mackey 2018) parametrized by its quality factor Q, its pulsation ω and the amplitude of the kernel function σ (the full expression of the kernel function is provided in Appendix A).This choice of mean and kernel function completely defines the GP, from which light curves with different levels of correlated noise can be drawn (see an example in Figure 1).The top plot shows the simulated flux time series of an active star drawn from the GP model described in section 1.1.This dataset corresponds to an observation of 4 days with an exposure time of 2 minutes.The mean of this signal consists in a periodic transit signal of period P = 0.7 days, duration D = 1.2 hours and depth ∆ = 2% (dark gray) plus an instrumental signal (blue).Correlated noise in the form of stellar variability is simulated with an SHO kernel of hyperparameters ω = π/6D, Q = 45 and σ = ∆.Finally, white noise with a standard deviation of 0.1% is added to the diagonal of the covariance matrix.

Transit detectability
One way to quantify the detectability of a unique transit signal is to compute its signal-to-noise-ratio expressed as where ∆ is the transit depth, σ the measurements uncertainty (assuming homoscedasticity), and n the number of points within transit.Although this metric is useful to assess the strength of the transit signal given a certain photometric precision, it does not account for the presence of correlated noise.However, instrumental and other astrophysical signals will necessarily affect the detectability of transits in realistic light curves (Pont et al. 2006).For this reason, the combined differential photometric precision (CDPP; Jenkins et al. 2010) metric was developed, and used in the context of the Kepler mission to assess the level of correlated noise in light curves, affecting the detectability of transits with a given duration.The CDPP is computed by decomposing the data in the time-frequency domain using wavelets, and measures the significance of a distorted transit signal in the whitened data.For a given transit duration, the CDPP is then a measure of the noise remaining after filtering the light curve in each frequency band, taking into account the presence of non-stationary correlated noise of a given timescale (see Jenkins et al. 2010 for more details).Hence, we can estimate the significance of the transit signal accounting for the presence of correlated noise as where CDP P D is the CDPP computed for a given transit duration D and ∆ is the transit depth after detrending.For simplicity, the CDPP is computed using the method from Gilliland et al. (2011)   An hourlong transit of depth 0.5% is generated on top of white noise (standard deviation of 0.15%) as part of a 24-hours observation with an exposure time of 1 minute (top).Then, in the bottom plot, correlated noise is added, generated using a GP with an SHO kernel of hyperparameters ω = 60, Q = 0.5 and σ = ∆/4 .The SNR on the right of each light curve is computed using Equation 2. Models used to simulate these data are provided in section 1.1.
Figure 2 shows the SNR from Equation 2 computed for a unique transit observed in the absence (grey) and presence of correlated noise (red).As illustrated in this figure, the presence of correlated noise strongly decreases the transit signal SNR, which would ultimately limit its detectability.

Detrending methods and their effects
The presence of instrumental correlated noise motivated the development of systematics detrending algorithms, such as the Trend Filtering Algorithm (TFA, Kovács et al. 2005, in its primary use case), SysRem (Tamuz et al. 2005) or Pixel Level Decorrelation (PLD, Deming et al. 2015; see also Everest from Luger et al. 2016Luger et al. , 2018)).Most of these methods rely on the shared nature of instrumental signals among light curves (or neighboring pixels) such that the correction applied should not degrade the transit signal and can be modeled using contemporaneous measurements (e.g.detector's temperature, pointing error, sky background or airmass time series).But even after instrumental signals have been removed, stellar variability and other astrophysical signals remain, which gave rise to several approaches.Some of them are physically-motivated and make use of GPs (e.g.Aigrain et al. 2016), others are empirical and make use of filtering and data-driven algorithms (Jenkins et al. 2010, Hippke et al. 2019).In this section, we show how these techniques impact transits detectability, depending on the morphological characteristics of light curves.
In Figure 3, we simulate a transit signal on top of which we add photometric stellar variability with different timescales, sampled from a GP with an SHO kernel described in Appendix A. For each light curve, we reconstruct and detrend stellar variability in two ways: one using the widely-adopted Tukey's bi-weight filter, presented in Mosteller & Tukey (1977) and using the implementation from wõtan2 (Hippke et al. 2019); the other using the same GP from which the data has been sampled.We then estimate the resulting transit depth and compute the remaining transit SNR using Equation 2. Figure 3 clearly shows the effect of both detrending techniques on transits SNR, and intuitively suggests that this degradation due to detrending is strongly dependent on the correlated noise characteristics encountered.Figure 3.In each plot on the left, correlated noises simulated using a GP with an SHO kernel are added (blue line) to an original light curve containing a transit of depth ∆ = 0.8%, duration D = 0.1 day, and white noise with a standard deviation of 0.15% (top left).The hyperparameters of the kernel are fixed to Q = 10 and σ = 0.8%, with the pulsation ω increasing from top to bottom.The resulting stellar variability signals are then reconstructed using (in purple) Tukey's bi-weight filter with an optimal window size of 3×D (see Hippke et al. 2019), and (in red) a GP with the same kernel used to simulate the data.
In each case the variability is reconstructed, subtracted, and the transit SNR computed using Equation 2.
To explore the parameter space for which detrending is the most problematic, we employ the model described in section 1.1 to simulate 10 000 differential light curves with different morphological characteristics, and compute the remaining transit SNR after detrending.In order to place the stellar variability hyperparameters on a relative scale with the transit signal parameters, we reparametrize the SHO kernel with where τ is the relative timescale of the variability with respect to the transit duration and δ the relative amplitude of the variability against the transit depth, both being adimensional.Hence, for (τ, δ) = (1, 1), the expressions of ω and σ given in Equation 3 correspond to a variability signal with a period half that of the transit duration, and a correlated noise amplitude comparable to the transit depth, i.e. strongly resembling the simulated transit signal.
For each of the 10 000 light curves generated, we separately reconstruct and detrend the variability signal using the two techniques employed in Figure 3, i.e. one using an optimal bi-weight filter with a window size three times that of the transit duration (Hippke et al. 2019) and the other using a GP with an optimal kernel.We then compute the detrended transit SNR using Equation 2.  Figure 4 shows that it exists an entire region of the (τ, δ) parameter space for which the bi-weight detrending degrades transit SNR to the point of no detection (SN R < 5).
Although being more robust, the same effect is observed when detrending with an optimal GP.Hence, detrending makes transit search blind to many systems, especially when the widely adopted bi-weight filter is used.We note that, outside this problematic parameter space, both techniques perform relatively well and such a degradation of the transit SNR should not be expected.Although this study should be extended to other detrending techniques, it highlights the need for a more informed transit search algorithm able to deal with correlated noise, at least if present in the form of stellar variability.

NUANCE
When searching for exoplanets using an indirect method, such as using transits or radial velocities, the initial detection is often as important as the follow-up observations, ultimately leading to the confirmation of the planetary candidate.For that reason, the final product of most transit search algorithms consists of a periodogram: a 1D detection metric given a set of trial periods that allows to identify the presence of a potential candidate, and the period and epoch at which it should be followed up.
If we assume that a transit is defined by its period P , epoch T 0 , duration D and depth ∆, we then wish to compute the likelihood of such a transit being present in the data for a set of periods P , leading to an interpretable periodogram.As we are interested in following up transits with specific parameters, such as a well-defined epoch, the likelihood we wish to compute must not be marginalized over all parameters other than P but rather specifically computed at the maximum-likelihood parameters T0 , D and ∆, leading to the profile likelihood where f is the data.In this section, we present nuance, an algorithm to compute transit-search periodograms in a tractable way, and detect planetary transits in light curves containing correlated noise such as instrumental signals and photometric stellar variability.In section 2.1, we explain how our approach requires two separate steps in order to remain tractable, and present these steps in section 2.2 and section 2.3, leading to the transit-search periodogram described in section 2.4.Finally, we present the nuance Python package in section 2.5 and perform a control test of our implementation against BLS in section 2.6.

The approach
As in section 1.1, let assume that f is a vector of size N containing the flux of a star observed at times t such that i.e. that f is drawn from a GP of mean Xw and covariance Σ.Again, we set the first M − 1 columns of the (N × M ) design matrix X as contemporaneous measurements and the last column as a normalized box-shaped transit of period P , epoch T 0 and duration D. This way, the transit signal is part of the linear model and its depth ∆ can be solved linearly.Under this assumption, the log-likelihood of the data given the presence of a periodic transit signal of period P , epoch T 0 , duration D, and a mean linear model with coefficients w is (Rasmussen & Williams 2005) ln Given the linearity of the mean model at P , T 0 and D fixed, this likelihood is maximized for the least-square parameters As we are only interested in the transit depth, we will omit w in the reminder of this paper and simply write likelihood expressions depending on its last value ∆ = w M , such as p(f |P, T 0 , D, ∆), where all other parameters of w are taken at their maximum-likelihood values.
Hence, given a period P , computing Q(P ) boils down to a non-linear optimization involving several evaluations of the likelihood given in Equation 5.While this is tractable for few periods, it becomes highly untractable for the large set of trial periods required to properly sample a transit-search periodogram (in the order of tens of thousands).
To remain tractable, nuance employs an extension of the strategy used by Foreman-Mackey et al. ( 2015), which has significant intellectual overlap with the methods used by Aigrain & Irwin (2004) and Jenkins et al. (2010).This approach separates the transit search into two components: the linear search and the periodic search.During the linear search, the likelihood p(f |T 0 , D, ∆) of a single non-periodic transit is computed for a grid of epochs and durations, each time solving for ∆ linearly.Then, the periodic search consists of combining these likelihoods to compute the likelihood p(f |P, T 0 , D, ∆) of the data given a periodic transit signal for a range of periods P .These combined likelihoods yield p(f |P, T0 , D, ∆), a transit-search periodogram on which the periodic transit detection is based.nuance differs from Foreman-Mackey et al. ( 2015) and other existing transit search algorithms as it models the covariance of the light curve with a GP, accounting for correlated noise (especially in the form of stellar variability) while keeping the model linear and tractable.This way, nuance searches for transits while, at the same time, modeling correlated noise, avoiding the separate detrending step that degrades transit signals SNR (see section 1).

The linear search
The goal of the linear search is to compute the likelihood p(f |T 0 , D, ∆) for a grid of epochs { T i } i and durations { D j } j .For each pair (T i , D j ) the transit depth ∆ i,j is linearly solved, which leads to the set of maximum likelihoods An example of such a grid of likelihoods is shown in Figure 5.To prepare for the next step, uncertainties on the depths σ i,j are also computed using Equation 6and stored.
For each pair of indices (i, j), the likelihood ln p(f |T i , D j , ∆ i,j ) is computed using the parameters from Equation 6and the expression of Equation 5.This process yields the grid of log-likelihoods ln L (bottom plot), as well as the { ∆ i,j , σ i,j } i,j transit depths and errors inferred linearly using Equation 6.

The periodic search
We then need to combine the likelihoods computed from the linear search to obtain i.e. the probability of a periodic transit of period P , epoch T 0 , duration D and depth ∆ given the data f .For a given transit duration D, any combination of (P, T 0 ) leads to K transits, for which it is tempting to write where { T k } k are the epochs matching (T 0 , P ) and { ∆ k } k the corresponding depths, so that This is the joint likelihood of transits belonging to the same periodic signal but with varying depths { ∆ k } k .However, individual transits from a periodic signal cannot be considered independent, and should instead be found periodically and share a common transit depth ∆.To this end, it can be shown (see Appendix B) that there is a closed form expression for the joint likelihood of K individual transits with depths and errors { ∆ k , σ k } k assuming a common depth ∆, corresponding to In order to compute this joint likelihood, we must assume that the likelihoods computed during the linear search are independent.However, the individual transits combined in the periodic search are not independent.Indeed, by using a GP, we assume that points belonging to separate transits may have a non-zero covariance.
In practice, we notice that this covariance is small enough to consider each transit as independent, a reasonable assumption for most physically-realistic datasets.
While Equation 8 takes a closed form, the individual epochs matching T 0 and P are not necessarily available in the grid of epochs { T k } k .In Foreman-Mackey et al. ( 2015), a similar issue is solved by using the nearest neighbors in the epochs grid.Instead, to allow the efficient matrix computation of Equation 8, we interpolate the likelihood grid from { T i } i to a common grid of transit phases { ϕ i } i , leading to the periodic search log-likelihood ln shown for few periods in Figure 6.In the latter equation, ∆ i,j is omitted since being interpolated from the linear search using ϕ i , D j and T 0 = 0.

The transit-search periodogram
Using Equation 8, we can now compute ln P for a range of periods and phases, and build a transit search periodogram using Equation 4. This has two disadvantages: First, each likelihood p(f |T 0 , D, ∆) estimated during the linear search is computed using N measurements.Hence, combining transits in the periodic search, through ∆ k , σ k and the product of K likelihoods { L k } k (see Equation 8), artificially leads to a likelihood involving up to N ×K measurements.For this reason, one has to normalize each likelihood ln P, by keeping track of the number of points used to compute each of them, which differs from one phase to another.Second, the maximum value of the likelihood ln P is relative to a given dataset, so that a more intuitive and absolute metric must be used to relate to the transit signal detection (such as the Signal Detection Efficiency in Kovács et al. 2002).This motivates a final step to produce an interpretable transit search periodogram Q. .For each period P , the joint likelihood P(P ) is computed using Equation 8, and the value of the maximum likelihood transit SNR retained as Q(P ).
For any period P, instead of taking Q(P ) as the maximum value of ln P, we compute the maximum likelihood parameters and define Q(P ) as being equal to the SNR of the transit of period P , epoch T 0 = ϕ 0 × P , duration D and depth ∆, i.e.
where ∆ and σ are obtained using Equation 8with the last column of X containing a periodic transit signal of period P , epoch T 0 , duration D and depth 1.This process and the resulting periodogram Q are shown in Figure 7.
Hence, periodic transit of period P with the maximum SNR, i.e. maximizing Q, is adopted as the best candidate, basing the confidence in this signal through its SNR.The parameters of this transit are the period P , the epoch T 0 = ϕ 0 P and duration D (Equation 9), and the depth ∆ with error σ (given by Equation 8).where gp is a tinygp GP instance and X the design matrix of the linear model.nuance exploits the use of tinygp4 , a Python package powered by JAX5 , allowing for custom kernels to be built and highly tractable computations.We can then define a set of epochs t0s and durations Ds and run the linear search with From this search object, the best transiting candidate parameters can be computed (search.best),or the Q periodogram retrieved (search.Q snr), together with valuable information about the transit search.The Nuance object also provides methods to perform transit search on light curves from multi-planetary hosts, the advantage of nuance being that the linear search only needs to be performed once and reused for the search of several transiting candidates (cf.section 3.3).An extensive and maintained online documentation is provided at nuance.readthedocs.io.

Comparison with BLS
To start testing nuance against existing methods, a simple adimensional normalized light curve is simulated, consisting of pure white noise with a standard deviation of 5 × 10 −4 spanning 6 days with an exposure time of 2 minutes.From this signal, we produce 4000 light curves, each containing box-shaped transits with periods randomly sampled from 0.3 to 2.5 days, durations of 50 minutes, and depths randomly sampled to lead to transit SNRs ranging from 4 to 30.For each light curve, transits are searched using two different tools: nuance, using its implementation from the Python package described in the previous section; and BLS, the Box-Least-Square algorithm from Kovács et al. (2002) (using astropy's BoxLeastSquares6 implementation).
For both methods, 3000 trial periods from 0.2 to 2.6 days are searched, with a single trial duration fixed to the unique known duration of 50 minutes.A transit signal is considered detected if the absolute difference between the injected and the recovered period is less than 0.01 day.To ease the detection criteria, orbital periods recovered at half or twice the injected ones (aka aliases) are considered as being detected.For this reason, detected transit epochs are not considered (although manually vetted).Results from this injection-recovery procedure are shown in Figure 8.
These results demonstrate the qualitative match between the detection capabilities of nuance and BLS on light curves with no correlated noise, where the BLS algorithm should be optimal.Explaining the subtle differences observed between the two methods when only white noise is present is beyond the scope of this paper, and we will assume that any differences observed in the following sections are due to the different treatments of correlated noise.

PERFORMANCE
Figure 4 shows that nuance's full-fledged modeling capabilities may not always be necessary and may only be beneficial for certain noise characteristics, relative to the searched transit parameters.Here, we evaluate the performance of nuance in the relative parameter space (τ, δ) described in Equation 3 and probe when its specific treatment of correlated noise in the transit search becomes necessary.
We perform this study by comparing nuance to the approach that involves removing stellar variability from light curves before performing the search on a detrended dataset.The following detrending strategies, each followed by a search with the BLS algorithm, are compared: bi-weight+BLS: employs an optimal bi-weight filter implemented in the wõtan Python package with an optimal window size of 3 × D, i.e. three times the transit duration (e.g.Hippke et al. 2019 andDransfield et al. 2024).
GP+BLS: employs a GP conditioned on the data (e.g.Lienhard et al. 2020).The kernel of the GP and its optimization is described on a case-by-case basis.
harmonics+BLS: employs a linear harmonic detrending, where the light curve is modeled as a Fourier series including four harmonics of the stellar rotation period with coefficients found through ordinary least square.
iterative+BLS iteratively detrends the light curve with a sinusoidal signal fitted to the data, each time using the dominant period of the residuals found using a lomb-scargle periodogram (5 iterations).
Like in previous sections, we consider the planet detected if the recovered period is within 0.01 day of the true period or a direct alias (such as two times or half the true period).Again, we ignore the exact match between the injected and recovered transit epochs, although we visually vet that the found epochs are consistent with the ones injected.We consider a transit detectable if its original SNR is greater than 6.Hence we define true positives has detectable transits recovered with the correct period (or an alias) and a measured SNR greater or equal to 6, and false positives as non-detectable transits recovered with a measured SNR greater or equal to 6.
As we noticed that few methods were still affected by the remaining stellar variability after detrending, the grid of orbital periods being searched for does not contain the stellar rotation period P * and its aliases.In practice this is done by removing all orbital periods P from the search grid such that dP = P P * is less than 2% from an integer value, i.e. |dP − ⌈dP ⌉| < 0.02.

Comparisons on simulated light curves
Our first comparison dataset consists in 4000 light curves simulated using the model described in section 1.1.We simulate a common periodic transit added to all light curves, of period P = 1.1 days, epoch T 0 = 0.2 days, duration D = 0.04 days and depth ∆ = 1%.Each light curve consists in a 4 days observation with an exposure time of 2 minutes, leading to N = 2880 data points with a normal error of 0.1%.
For a given pair of (τ, δ), we simulate stellar variability using a GP with an SHO kernel of hyperparameters defined by Equation 3 computed with respect to the injected transit parameters D and ∆.The same kernel is used for the search with nuance and with the GP+BLS method, an optimal choice on equal footing with the optimal 3 × D window size of the bi-weight filter employed in the bi-weight+BLS search.The 4000 pairs of (τ, δ) are drawn from τ ∼ U(0.1, 10), δ ∼ U(0.1, 25) and Q ∼ U(10, 100) where U(a, b) denotes a uniform distribution of lower bound a and upper bound b.

RESULTS
The results of this injection-recovery procedure are shown in Figure 9 and highlight particularly well the benefit of nuance against most other methods.Except for the GP+BLS approach, nuance lead to a much higher rate of true positives for transits with relatively small depths compared to the stellar variability amplitude (i.e.δ > 2), and a duration comparable to the stellar variability period (i.e.τ < 2).On the other hand, the GP+BLS strategy seems almost as performant as nuance, recovering most of the injected transits and only lacking detection in a relatively comparable portion of the (τ, δ) parameter space.We note that these empirical statements only concerns simulated light curves with a given amount of white noise, and may vary depending on the length of the observing window or the number of transits.For this reason, quantifying for which values of (τ, δ) nuance outperforms these methods would only apply to this specific but representative example.However, we verify that our conclusions remain qualitatively valid for various simulation setups.This injection-recovery is done in a particularly optimal setup, on simulated light curves not all physically realistic and using an optimal GP kernel, hence demonstrating the performance of nuance only on a purely experimental basis.In the next section, we perform transits injection-recovery on real space-based light curves.

Comparisons on rapidly-rotating M dwarfs TESS light curves
In order to assess the performance of nuance on real datasets, we inject and recover transits into light curves from the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015).We focus this proof of concept on light curves from a list of 438 M-dwarfs found to have detectable rotation signals with periods lower than one day (Ramsay et al. 2020), which lead to a parameter space justifying the use of nuance.
For each studied target, transits are injected and recovered in the TESS 2 min cadence SPOC Simple Aperture Photometry and Pre-search Data Conditioning light curves (PDCSAP, Caldwell et al. 2020) of a single sector (the first being observed for each target) spanning on average 10 days.To our knowledge, none of these 438 targets have been searched for planetary transits before.However, we note that the presence of existing transit signals in these light curves before the injection of simulated ones is possible, but will not affect the relative comparison of one method to another.

LIGHT CURVE CLEANING AND TRANSITS INJECTION
As some of the techniques compared to nuance can be affected by gaps in the data, we only use continuous measurements from half a TESS sector.We assume that all methods (including nuance) are based on an incomplete model of the data that does not account for stellar flares.For this reason, the light curve of each target is cleaned using an iterative sigma clipping approach.For each iteration, points 3 times above the standard deviation of the full light curve (previously subtracted by its median) are identified.Then, the 30 adjacent points on each side of the found outliers are masked.This way, large flare signals are masked, using a total of 3 iterations.At each iteration, the GP kernel hyperparameters are re-optimized.As PDCSAP light curves often start with a ramp-like signal, the first 300 points (as well as the last 300 points) of each continuous observation are masked.Finally, each light curve is normalized by its median value.Such a cleaned light curve is shown in Figure 10.We note that the gaps left after sigma clipping may be problematic for some of the detrending techniques (such as bspline+BLS).However, adopting this flare cleaning step and analyzing light curves with few small gaps is a practice commonly found in the literature.
For each light curve, transits of planets with 10 different orbital periods combined with 10 planetary radii are individually considered, for a total of 100 periodic transits injection-recovery per target.Orbital periods P are sampled on a regular grid between 0.4 and 5 days, and planetary radii R p are sampled on a regular grid designed to yield a minimum transit SNR of 2 and a maximum of 30.Using Equation 2 with σ r = 0, the planetary radius leading to a transit with a desired SNR s is given by with σ equal to the mean uncertainty estimated by the SPOC pipeline, R ⋆ the radius of the star reported by Ramsay et al. (2020) and n the number of points in transit computed using a transit duration assuming a circular orbit.

STELLAR VARIABILITY KERNEL
In the GP+BLS and nuance methods we model light curves using a GP with a mixture of two SHO kernels of period P ⋆ and P ⋆ /2 where P ⋆ is the rotation period of the star.This model is representative of a wide range of stochastic variability in stellar time series9 (e.g.David et al. 2019;Gillen et al. 2020).In order to account for additional corelated noises, we complement this kernel with a short and a long-timescale exponential term, so that the full kernel can be expressed as • k 2 a SHO kernel with hyperparameters where Q 0 is the quality factor for the secondary oscillation, δQ is the difference between the quality factors of the first and the second modes, f is the fractional amplitude of the secondary mode compared to the primary and σ is the standard deviation of the process.The kernels k 3 and k 4 are expressed as with ℓ and σ the scale and standard deviation of the process.These are meant to model short and long-timescale non-periodic correlated noise.In total, the rotation kernel k has 8 hyperparameters.
The hyperparameters of this kernel are optimized on trimmed and cleaned light curves containing the injected transits, using the scipy.optimize.minimizewrapper provided by the jaxopt Python package 10 , and taking advantage of the JAX implementation of tinygp and its quasiseparable kernels 11 (Foreman-Mackey et al. 2017).
As correlated noise is expected to affect the light curve uncertainty estimates performed by SPOC, the diagonal of the full covariance matrix of the data (i.e.their uncertainty, assuming homoscedasticity) is held free, increasing the number of optimized parameters to 9. The optimization is performed using the BFGS algorithm (Fletcher 1987), minimizing the negative log-likelihood of the data as expressed in Equation 5 (without transit), i.e. accounting for a linear systematic model of the data in addition of stellar variability.For simplicity, and to adopt a uniform treatment for all target light curves, a design matrix X with a single constant column is adopted, such that the systematic model only consists in a single parameter corresponding to the mean value of the differential flux (expected to be close to 1) solved linearly.Our motivations to choose this very simplistic baseline, despite the capability of nuance to account for more complex linear models, is discussed in section 4.2.

SEARCH PARAMETERS AND TRANSIT DETECTION CRITERIA
For all techniques, we only search for transits with a duration fixed to the known duration of the injected transits.This is mainly done for computational efficiency and allows for a narrower comparison between nuance and all BLS-based techniques.
Finally, the search is done over 4000 trial periods linearly sampled from 0.3 to 6 days.A realistic transit search on a wider parameter space (e.g.multiple trial durations) is discussed in the next section.

RESULTS
An example of the transit injection-recovery and its result is shown in Appendix C for the target TIC 1019692.Figure 11 shows the global results of the injection-recovery for all targets, while Figure 12 shows the true and false positives plotted against the relative parameters τ and δ (defined in Equation 3).These figures are a synthesis of Figure 13 that shows the rate of true and false positives binned in the full parameter space (τ, δ).Compared to other techniques, we find that nuance leads to the highest number of true positives, with a successful detection of 76% of the 7008 detectable transits injected (Figure 11).While the performances of other methods strongly depend on the characteristics of the variability, nuance is the best technique in 93% of cases, leading to both the highest number of true positives and the lowest number of false positives (Figure 11).From Figure 11, we note that the number of true positives of the GP+BLS method is much worse that what could be anticipated from the results presented in Figure 9, where the GP kernel was fully optimal (given it was also used to simulate the data).In that case, our kernel might not be optimal, whether because of its form or because of the values of its hyperparameters.In any case, the fact that the same kernel performs significantly better when used with nuance shows that our method is less sensitive to the choice of kernel and its optimization compared to detrending with a GP.In order to further validate nuance on a realistic dataset, we focus this section on the multi-sector TESS light curves of TOI-540, and the search for its earth-like companion TOI-540 b (Ment et al. 2021).We download the 2 minutes cadence SPOC PDCSAP light curves of TOI-540 observed in 5 sectors (4, 5, 6, 31 and 32).Like in the previous section, we use a Lomb-Scargle periodogram and identify the 0.72 days rotation period of the star, that we use as an initial value to optimize the kernel described in section 3.2 on each sector independently.Here again, we employ an upper-sigma-clipping to mask flares out of the data.The resulting light curve for sector 4 and its mean model are shown in Figure 14.
For each sector, we perform the linear search of nuance on the cleaned light curve, using the original times as the trial epochs and 10 trial transit durations linearly sampled from 15 minutes to 1.5 hours.We then perform the periodic search on all sectors combined, using a concatenation of all linear searches.By adopting this by-sector GP modeling of the light-curve, the linear search of nuance scales linearly with the number of sectors being processed.This approach is adopted for efficiency but also to encapsulate the changing properties of stellar variability from one sector to another, often separated by year-long gaps.The periodic search is done on 20 000 trial periods ranging from 0.5 to 10 days12 .This search, using nuance, is compared to the more traditional approach that consists in detrending each sector with a bi-weight filter and then search for transits with the BLS algorithm (denoted bi-weight+BLS and described in section 3.2) on all sectors combined.Since we don't know the transit duration a-priori, we perform the detrending and BLS search using 15 filtering window sizes sampled from 30 minutes to 5 hours, and retain the search that leads to the highest transit SNR peak in the periodogram (as done e.g. in the SHERLOCK transit search pipeline described in Pozuelos et al. 2020).The results of this comparison are shown in Figure 15.
After a first periodic search, trial epochs in windows of widths 2 × D centered on the detected periodic transits are masked.In practice, this is done by masking the linear search products { ln L i,j } i,j , { ∆ i,j } i,j and { σ i,j } i,j (defined in section 2.2).
As seen in Figure 15, the SNR periodogram using bi-weight+BLS and nuance are very similar, with the known transiting exoplanet TOI-540 b detected with an orbital period P = 1.24 days.This is well expected as the relative parameter τ equal 13 for TOI-540 b13 , which lies outside the range where nuance is expected to be beneficial (see Figure 4).Nonetheless, the nuance periodogram of TOI 540 features less spurious SNR peaks, largely due to the penalty naturally occurring when single transits with different depths are periodically combined.In the second search (right panel) of Figure 15, we also notice a higher number of peaks that would lead to false positive detections of transits in the bi-weight+BLS case.The proper treatment of correlated noise in nuance, as observed in section 3.2, makes these peaks non-significant, avoiding a large number of false detections.
We note that finding a known TESS candidate that displays characteristics for which nuance is expected to be significantly beneficial proved to be very challenging during the writing of this paper, as transits with such characteristics are expected to be missed by current state-of-the-art techniques (see e.g. Figure 9).Nontheless, we verify that nuance is capable of finding a large number of already known transiting exoplanet candidates, in light curves featuring various forms of correlated noises, at least as efficiently as with commonly used techniques.We reserve the search for new transit signals to a follow-up publication.

DISCUSSION
In the previous section we demonstrated the capability of nuance to search for synthetic or known transit signals, in simulated or real datasets.Here, we discuss the caveats of this algorithm, the advantages and limitations of the nuance implementation, and future prospects for its extension.

Processing time
In Figure 16, the processing time of nuance linear and periodic search are recorded against the number of points in a simulated light curve, assuming a simple nonoptimized GP with a squared exponential kernel.These are compared to the processing time of bi-weight+BLS (cf.section 3.1) separated into the bi-weight filtering step and the BLS search. .Study of nuance (black line) processing time against bi-weight+BLS (cf.section 3.1, blue line).The gray curve shows the performance of nuance linear search when applied to chunks of 10 000 points continuous observations, instead of considering these observations all together.This study is performed on a single CPU core of an Apple M2 Max chip.While not being shown, we verify that both the BLS algorithm and nuance processing times scale linearly with the number of trial transit durations and orbital periods.
As see in Figure 16, most of the computational costs of nuance and the bi-weight+BLS method come from the linear search and the bi-weight filtering step.This is not always true and depends on the size of the trial durations and periods grids.One advantage of nuance is that the linear search can be performed separately on different continuous observations, and then combined in the periodic search.Hence, if searching for transits in separate observations with approximately similar durations, such as different TESS sectors or different ground-based nightly observations, the computational cost of nuance grows linearly with the number of observations (see grey line in Figure 16).Nonetheless, considering a variety of other detrending algorithms, nuance is expected to perform between one and two order of magnitude slower than more conventional techniques associated with BLS.
Nonetheless, searching for TOI-540 b transits with the parallelized implementation of nuance (see section 3.3) on the 12 cores of an Apple M2 Max chip, took only 5 minutes and 35 seconds, 1 minute 50 seconds for the 5 sectors linear searches (around 22 seconds per TESS sector), and 3 minutes 45 seconds for the combined periodic search.In comparison, the brute force search with bi-weight+BLS, consisting in trying 15 bi-weight windows, took a total of 1 minute 49 seconds (only 7 seconds each).
Because of its computational cost, we do not recommend using nuance in the general case, but rather when light curves contain correlated noise with specific characteristics.If employing a bi-weight filter for detrending, these characteristics correspond to the ones discussed in section 1.But as these strongly depend on the type of detrending technique employed (see Figure 9), we do not provide general guidelines as when nuance should be preferred over a specific detrending technique.To aid users in making an informed choice of algorithm, extensive benchmarks and guidelines are reserved to future developments and will be progressively shared on nuance's online documentation14 .

Systematics modeling
Throughout the paper, a single column design matrix X, corresponding to the mean differential flux (ideally unitary), was employed, hence assuming that the instrumental systematics signals were non-existent.In practice, nuance has been developed to linearly model systematic signals through more complex design matrices (as in Foreman-Mackey et al. 2015), in addition with its capability to model correlated noise while searching for transits.This feature is intentionally unexploited in the comparisons presented in section 3, as detrending light curves assuming a linear systematics model, such as PLD co-trending vectors (Deming et al. 2015), is highly incomplete if applied on data while ignoring the presence of other astrophysical signals.Comparisons involving more complex design matrices would also be sensitive to the choice of linear components, and would have unwanted repercussions on their results.
As an illustration, the NEMESIS pipeline (Feliz et al. 2021) starts processing the differential light curves by employing a linear systematics detrending using a least square fit of the data with a reduced PLD basis, before smoothing the signal from stellar variability using an approach similar to the one employed in the bi-weight+BLS approach (cf.section 3.2), hence detrending the systematics with an incomplete model that does not account for stellar variability.To account for stellar variability while fitting the linear systematics model to the data, a step further would be to use a GP, such as done in the EVEREST (Luger et al. 2018) pipeline.However this would also involve some potential degradation of the transit signals (see e.g. Figure 3), with an hardly distinguishable origin.For these reasons, and to keep our comparisons as targeted as possible, we do not compare commonly used systematics detrending approaches and decided to focus our comparisons solely on stellar variability detrending techniques (although these two aspects often overlap in the literature, e.g. in Luger et al. 2016).
Although not being demonstrated here, modeling systematics signals while searching for transits on data acquired sparsly is extremely promising for the search of transiting exoplanets, including for ground-based observations that usually suffer from daily interruptions.In this respect, we note the similarity of our linear search (cf section 2.2) to the one presented in Berta et al. 2012, that focused on the detection of single eclipses in the MEarth light curves (Irwin et al. 2009).Similarly, nuance would highly benefit the search for transiting exoplanets around M-dwarf type stars, such as the ones observed by the SPECULOOS survey (Sebastian et al. 2021) whose monitoring suffers from both increased red noise (due to atmospheric and instrumental thermal effects discussed e.g. in Berta et al. 2012 andPedersen et al. 2023) and enhanced stellar variability (Murray et al. 2020).We reserve this promising application to a future study.

The GP kernel
While not being discussed in our study, the efficiency of nuance to detect transits in correlated noise might be dependant on the design of its GP kernel.In the ensemble comparison of section 3.2, the goal was to choose a kernel and an optimization strategy suited to most of the studied light curves, leading to few outliers in the results that were indicative of a badly designed and/or optimized kernel.An alternative, recommended for more realistic blind searches, is to perform model comparison on well selected kernels, and to adapt the optimization strategy to each dataset.
When using nuance on TESS light curves for example, it must be noted that the observed light curve variability might originate from contamination due to several nearby blended stars, so that a physically-interpretable GP kernel representing a single star activity is not necessarily appropriate.On the other hand, a single squared exponential GP kernel might also be sufficient for some applications; an aspect we intend to explore in future applications.
Something to emphasize is that nuance cannot be used to produce reliable detrended light curves, as an optimal GP is often flexible enough to partially model transits (see Figure 3).In contrast, the idea behind nuance is rather to compute the likelihood of data against a model containing both transits and correlated noise, without ever trying to disentangle both signals.In practice, it means that it can be very hard to actually verify the presence of transits found by nuance visually, so that transits may be detected but remain hidden in correlated noise.This is particularly true for stars displaying very high-frequency photometric pulsations (see the example in Figure 17).

Prospects
The present implementation of nuance has the potential to be extended to be used beyond the search of periodic box-shaped transits.Here are ideas of possible use and extensions, from the most straightforward to the most ambitious: 1.In order to compare nuance to BLS-based methods, we injected and retrieved only box-shaped transits.However, similarly to the Transit-Least-Square algorithm from Hippke & Heller (2019), limb-darkened transits can serve as base model in the linear search and are expected to improve the transit search in the same way TLS provided an improvement over BLS.
2. The linear search of nuance is a single-event detection algorithm that can be used to search for single transit events, but also detect transiting exocomets and flares, by simply changing the base astrophysical model in the last colum of the design matrix X.While not being tested and benchmarked, nuance already integrates this feature.An example of the detection of known transiting exocomets in β Pictoris TESS light curves is shown in Figure 17.
3. During the periodic search, no prior about the transit duration related to the orbital period of the planet was used.This was done in order to allow the detection of grazing transiting exoplanets that would produce shorter timescale transits compared to what is expected from a circular orbit with a null impact parameter.However, adding such priors might produce less spurious periodogram peaks and be very beneficial for the automatic search of transits in large datasets.Another idea, similar to the one employed in Foreman-Mackey et al. (2015), is to leverage model comparison in order to reject transits that are better described by the GP model alone.Both ideas come at no cost given our modeling approach.
4. Finally, the formalism of nuance could be adapted and used to search for exoplanets featuring transit time variations (TTV).Indeed, this application only requires a modification of the periodic search, as maximum likelihood peaks close to linearly predicted transit epochs may be considered.This could be done either with a special nearest-neighbor algorithm or with a convolution of the computed likelihood grid with a Gaussian kernel.To maximize efficiency and interpretability, we would recommend these approaches to be explored analytically, rather than using a data-driven treatment of the linear search products.

CONCLUSION
This paper presents nuance, an algorithm designed to detect planetary transits in light curves featuring correlated noise in the form of instrumental signals and stellar variability.In this context, a conventional approach involves detrending a light curve before searching for transits using a Box-Least-Square algorithm.However, we show that this approach degrades transit signal-to-noise-ratios down to the point of not being detectable.Adopting commonly used detrending strategies, we explore the extent of this degradation on simulated light curves, and its dependance on the photometric stellar variability characteristics, showing the need for a full-fledged transit search method like nuance.
The effectiveness of nuance is tested using a synthetic dataset and further validated on real TESS light curves of 438 rapidly-rotating M dwarfs.These injection-recovery tests reveal that nuance consistently outperforms commonly used transit search techniques, especially when the timescale of stellar variability is less than twice that of the transit duration.In all cases, nuance not only leads to a higher number of true positive detections but also minimizes false positives, demonstrating its robustness and reliability.
We make nuance publicly available through the nuance open-source Python package, developped with JAX to allow its use on distributed computing environments and GPU devices.Overall, we acknowledge the limitations of nuance and its increased computational cost compared to more conventional techniques.Hence, nuance should be used as an alternative to more traditional techniques only in the presence of substantial correlated noise.As guidelines for choosing our method over other techniques are highly dependant on the type of detrending algorithm employed, we reserve this study to a future work.
Finally, we suggest future improvements and extensions of the algorithm, including its application for detecting single transits, exocomets, flares, and exoplanets featuring transit time variations (TTV), underscoring its versatility and potential for broader impact in astronomical research.
The software presented in this work is open source under the MIT License and is available at https://github.com/lgrcia/nuance,with documentation and tutorials hosted at https://nuance.readthedocs.io.
We would like to thank Julien De Wit and Prajwal Niraula for meaningful discussions at the beginning of this project, Germain Garcia for his useful insights about numerical optimization, Michaël Gillon for his overall support, the member of the Exotic group at the University of Liège, and the member of the Astronomical Data Group at the Center for Computational Astrophysics for many enriching discussions and feedback.This publication benefits from the support of the French Community of Belgium in the context of the FRIA Doctoral Grant awarded to L.J.G.Software: numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), jax (Bradbury et al. 2018), jaxopt (Blondel et al. 2021) Since each depth ∆ k is found through generalized least square, each follow a normal distribution N (∆ k , σ 2 k ), centered on ∆ k with variance σ 2 k and an amplitude L k , leading to the likelihood function As for the common transit depth ∆, it can be estimated through the joint probability of all other transit depths than ∆ k , such that Hence We can now rewrite Equation B2 as The integral in this equation is a product of gaussian integrals that can be obtained analytically, leading to the log-likelihood of the data given a periodic transit signal of period P , epoch T 0 , duration D and common depth ∆.In order to reduce the number of times Equation B3 is computed, we adopt the biased estimates so that ∆ and σ are independent of k in the last sum of Equation 8. Figure 18.Results of the transits injection-recovery on TIC 1019692 half-sector light curve.Left: cleaned light curve with computed trend overplotted in black (except for nuance where it corresponds to the mean of the GP model).Right: Results of the transit search where a black square denotes a transit signal not detected, gray a signal detected at an alias period (P/2 or 2P ), and white a signal detected with the correct period.On the right plots, secondary axes show the (τ, δ) relative parameter space.For each method, the upper right legend on the left plot indicates its ranking based on the percent of recovered transit signals (where a transit with an aliased period counts as being detected).
Figure1.The top plot shows the simulated flux time series of an active star drawn from the GP model described in section 1.1.This dataset corresponds to an observation of 4 days with an exposure time of 2 minutes.The mean of this signal consists in a periodic transit signal of period P = 0.7 days, duration D = 1.2 hours and depth ∆ = 2% (dark gray) plus an instrumental signal (blue).Correlated noise in the form of stellar variability is simulated with an SHO kernel of hyperparameters ω = π/6D, Q = 45 and σ = ∆.Finally, white noise with a standard deviation of 0.1% is added to the diagonal of the covariance matrix.

Figure 2 .
Figure2.Illustration of the effect of correlated noise on a single transit SNR.An hourlong transit of depth 0.5% is generated on top of white noise (standard deviation of 0.15%) as part of a 24-hours observation with an exposure time of 1 minute (top).Then, in the bottom plot, correlated noise is added, generated using a GP with an SHO kernel of hyperparameters ω = 60, Q = 0.5 and σ = ∆/4 .The SNR on the right of each light curve is computed using Equation2.Models used to simulate these data are provided in section 1.1.

Figure 4 .
Figure 4. Transit SNR computed for 10 000 differential light curves with varying morphological characteristics.All light curves span 2.8 days with an exposure time of 2 minutes and contain the same transit signal of duration D = 1 hour, depth of 0.5% and white noise with a standard deviation of 0.1%.Light curves at the top and right side of the central plot are shown with their corresponding τ and δ values, which again correspond to the relative timescale of stellar variability against the transit duration, and the relative amplitude of stellar variability against the transit depth.The two central plots show the transit SNR after detrending the light curves with one of two methods.On the left, light curves are detrended using Tukey's bi-weight filter with an optimal window size of 3 × D (see Hippke et al. 2019), while on the right, light curves are detrended using a GP with the same kernel used to simulate the data.

Figure 5 .
Figure 5. Principle and output of the linear search.The simulated dataset (top) corresponds to the one shown and described in Figure1.First, a set of durations and depths { T i , D j } i,j is generated.For each pair of indices (i, j), the likelihood ln p(f |T i , D j , ∆ i,j ) is computed using the parameters from Equation6and the expression of Equation5.This process yields the grid of log-likelihoods ln L (bottom plot), as well as the { ∆ i,j , σ i,j } i,j transit depths and errors inferred linearly using Equation6.

Figure 6 .
Figure 6.Periodic search likelihood P(P ) computed for different trial periods P .Notice how the maximum value of P for the alias period P = 0.35 day (left plot) is lower than for P = 2 × 0.35 = 0.7 day, a result of combining the log-likelihoods using Equation 8 instead of Equation 7, in favor of individual transits matching a common depth ∆.
Figure7.For each period P , the joint likelihood P(P ) is computed using Equation8, and the value of the maximum likelihood transit SNR retained as Q(P ).

Figure 8 .
Figure 8. Binned statistics of the injection-recovery of 4000 transit signals in a flat light curve with only white noise using BLS and nuance.The color scale indicates the recovery of transits in the corresponding (period, depth) parameter space, white for a full recovery and black for no detection.

Figure 9 .
Figure 9. Rate of true positives for the 4000 simulated light curves.The color scale represents the fraction of true positives, white if all injected transits are recovered in a given portion of the (τ, δ) parameter space, black if none are recovered.

Figure 10 .
Figure 10.Trimmed and cleaned single-sector light curve of the target TIC 1019692.The top plot shows how much of the data is truncated and sigma clipped, resulting in a quasi-continuous light curve shown in the bottom plot.On this bottom plot, the black line corresponds to the mean of the GP model (with hyperparameters optimized here on the cleaned light curve).

Figure 12 .
Figure12.Rate of true and false positives of nuance compared to other methods, as a function of the relative parameters τ and δ.The solid blue line corresponds to the maximum of true positives among all methods on the top panel (dominated by the iterative+BLS method), and the minimum of false positives on the bottom panel (waned by the GP+BLS method, all other methods being above nuance).

Figure 13 .
Figure 13.Rate of true and false positives binned over the full parameter space (τ, δ) for all methods and all targets.

Figure 14 .
Figure 14.Sector 4 light curve of TOI 540.The cleaned signal (gray points) has been masked for flares (black points), and the black line corresponds to the mean of the GP model.

Figure 15 .
Figure 15.Transit search SNR periodograms of TOI-540 using bi-weight+BLS and nuance.After a first periodic search (left panel), the epochs corresponding to the maximum-SNR transit are masked before the second search is performed (right panel).
Figure16.Study of nuance (black line) processing time against bi-weight+BLS (cf.section 3.1, blue line).The gray curve shows the performance of nuance linear search when applied to chunks of 10 000 points continuous observations, instead of considering these observations all together.This study is performed on a single CPU core of an Apple M2 Max chip.While not being shown, we verify that both the BLS algorithm and nuance processing times scale linearly with the number of trial transit durations and orbital periods.

Figure 17 .
Figure17.Demonstration of the search of transiting exocomets in the TESS light curves of the star β Pictoris (2 minutes PDCSAP data from sector 33).This star is known to display rapid δ Scuti type photometric variations with a period of about 30 minutes (Lecavelier desEtangs et al. 2022).The GP kernel and its hyperparameters are chosen and optimized as in section 3.2.Here, we simply used the linear search of nuance with a different base model, one that mimics the shape of transiting exocomets (with 20 trial durations), to compute the SNR time-series of the signal over one sector.The maximum-SNR events are displayed in the bottom of the figure and match with the ones found by Lecavelier desEtangs et al. (2022).
F.J.P acknowledges financial support from the grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033 and through projects PID2019-109522GB-C52 and PID2022-137241NB-C4. S.A. acknowledges funding from the European Research Council (ERC) under grant agreement 865624 (GPRV) and from the UK Science and Technology Facilities Council (STFC) under grants ST/S000488/1 and ST/R004846/1.This paper includes data collected by the TESS mission.Funding for the TESS mission is provided by the NASA's Science Mission Directorate and NASA Explorer Program.We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center.Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products.TESS data were obtained from the MAST data archive at the Space Telescope Science Institute (STScI).STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation.