Motivation#
Depending on the noise characteristics, detrending light curves before searching for transits degrades their signal-to-noise-ratio and makes them undetectable.
Let’s show an example of that by generating a light curve containing a periodic transit signal and some correlated noise
import jax
import numpy as np
from nuance import core
import matplotlib.pyplot as plt
from tinygp import GaussianProcess, kernels
depth = 1.2e-2
error = 0.004
time = np.linspace(0, 4.0, 2000)
transit_params = dict(epoch=0.2, duration=0.05, period=0.7)
transit_model = depth * core.transit(time, **transit_params)
kernel = kernels.quasisep.Matern32(sigma=3e-2, scale=transit_params["duration"])
gp = GaussianProcess(kernel, time, diag=error**2)
flux = transit_model + gp.sample(jax.random.PRNGKey(3)) + 1.0
plt.figure(figsize=(8.5, 3))
plt.plot(time, flux, ".", c="0.8")
plt.plot(time, transit_model + 1, c="k", label="transit signal")
plt.xlabel("time (days)")
plt.legend()
_ = plt.ylabel("normalized flux")

Detrending#
We will now detrend the light curve using Tukey’s bi-weight filter, with an optimal window length 3 times that of the transit duration (see Hippke et al. (2014)).
import wotan
detrended_flux = wotan.flatten(
time, flux, method="biweight", window_length=3 * transit_params["duration"]
)
plt.figure(figsize=(8.5, 3))
plt.plot(time, detrended_flux, ".", c="0.8")
plt.xlabel("time (days)")
_ = plt.ylabel("normalized flux")

That doesn’t look great. We might be tempted to use a smaller window_length
, but we will take the risk of detrending part of the transit itself.
BLS#
Let’s try to run the Box-Least-Squares algorithm (Kovacs et al. 2022) on the detrended flux.
from astropy.timeseries import BoxLeastSquares
def bls(time, flux, error, objective=None):
model = BoxLeastSquares(time, flux, dy=error)
periodogram = model.autopower(transit_params["duration"], objective=objective)
periods = periodogram.period
power = periodogram.power
return periods, power
periods, snr = bls(time, detrended_flux, error)
plt.plot(periods, snr, "k")
plt.axvline(transit_params["period"], color="k", alpha=0.3, label="true period")
plt.legend()
plt.xlabel("period (days)")
_ = plt.ylabel("BLS power")

As might be expected, we miss the true signal. In order to detrend the light curve in the least-destructive way, let’s try different window lengths
import numpy as np
from tqdm.auto import tqdm
window_lengths = np.linspace(0.1, 4, 15)
results = {}
for window_length in tqdm(window_lengths):
detrended_flux = (
wotan.flatten(time, flux, method="biweight", window_length=window_length) + 1.0
)
periods, power = bls(time, detrended_flux, error)
results[window_length] = {
"periods": periods,
"snr": power,
}
and see which one gives us the best periodogram
best = max(
window_lengths, key=lambda window_length: results[window_length]["snr"].max()
)
plt.plot(results[best]["periods"], results[best]["snr"], "k")
plt.axvline(transit_params["period"], color="k", alpha=0.3, label="true period")
plt.legend()
plt.xlabel("period (days)")
_ = plt.ylabel("BLS power")

Not great either. In the nuance paper, we show that, depending on the noise characteristics, many detrending techniques (including using a bi-weight filter) degrade the transit signal-to-noise-ratio up to the point of not being detectable. Before we move on and apply nuance to this example, let’s try one last method.
GP detrending#
Assuming we have access to the optimal Gaussian Process (GP) kernel that perfectly describes the covariance of the data (we do have access to that here), we can use the original GP used to simulate the light curve to detrend it
import jax
def trend():
return gp.condition(flux - 1, time).gp.mean
detrended_flux = flux - jax.jit(trend)()
plt.figure(figsize=(8.5, 3))
_ = plt.plot(time, detrended_flux, ".", c="0.8")

To the eye, it looks like all transits have been removed. Let’s check the BLS periodogram
periods, power = bls(time, detrended_flux, error)
plt.plot(periods, power, "k")
plt.axvline(transit_params["period"], color="k", alpha=0.3, label="true period")
plt.legend()
plt.xlabel("period (days)")
_ = plt.ylabel("BLS power")

Ok if we were ignoring smaller periods but not great on the full range.
Warning
This is a cherry-picked example 🍒!
Detrending with a Gaussian Process can be great (see nuance paper) but only assuming you perfectly know the covariance structure of the data, which is never the case.
nuance#
Let’s now run nuance on this simulated dataset. With nuance, the transit signal is searched while simultaneously modeling correlated noise, assuming that the light curve can be modeled as a Gaussian process. The first step of nuance consists in the linear search where the likelihood of a single non-periodic transit signal is computed at different epochs and durations.
from nuance.linear_search import linear_search
epochs = time.copy()
durations = np.array([0.01, transit_params["duration"]])
ls = linear_search(time, flux, gp=gp)(epochs, durations)
Important
nuance relies on tinygp to manipulate Gaussian processes. Although gp
in the previous cell is already a tinygp.GaussianProcess
instance (from the simulated data), in practice one has to build their own.
These likelihood are then combined during the periodic search
from nuance.periodic_search import periodic_search
snr_function = jax.jit(core.snr(time, flux, gp=gp))
ps_function = periodic_search(epochs, durations, ls, snr_function)
snr, params = ps_function(periods)
plt.plot(periods, snr, "k")
plt.axvline(transit_params["period"], color="k", alpha=0.3, label="true period")
plt.xlabel("period (days)")
_ = plt.ylabel("SNR")
/opt/homebrew/Caskroom/miniforge/base/envs/nuance/lib/python3.10/site-packages/multiprocess/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
self.pid = os.fork()

Here nuance leads to the successful detection of the true transit signal. Let’s plot the found signal against the original one
best_params = params[np.argmax(snr)]
linear, found, noise = core.separate_models(time, flux, gp=gp)(*best_params)
plt.figure(figsize=(8.5, 3))
plt.plot(time, flux, ".", c="0.8", label="data")
plt.plot(time, transit_model + 1.0, "k", label="true")
plt.plot(time, found + 1.0, "C0", label="found")
plt.legend()
plt.ylabel("normalized flux")
_ = plt.xlabel("time (days)")

Note
In this example we employ a perfectly optimal GP kernel. But as we show in the paper, nuance is also robust to non-optimal kernels.
Attention
nuance is not a detrending algorithm and does not produce detrended light curves! What nuance does is to compute the likelihood of a transit being present in some correlated noise, without disentangling the two.