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")
../_images/32f0757360a5433a19a5be599059273336b7371677c74a60cab2db88fefe169e.png

Detrending#

Note

The following requires the installation of the following packages:

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")
../_images/128d0fecd0335981108aed5a8eb4a041ac58bcf211c698eca9084f75ad6708c3.png

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")
../_images/4245467b25b369c35b966fa269c51fa92fd422c7f8dce301aef5ea5f58a27e89.png

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")
../_images/113d044863afed095d89d803ccf0856e2ceab272918ce216e1ad34b65bb57b33.png

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")
../_images/cd87e48d0f5a8f18c10276f0c951471da2432967acfb1730e50b66b883ae2131.png

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")
../_images/49e385192fd22e68c7eded435371de16b078a8502eed05b66a2a087ac0c07ad2.png

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()
../_images/897814acf188a180ed6a4aef08b42045daf16833ae9f9d8160ab4cea94050973.png

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)")
../_images/4144ff87682ca005202c276dea7c2fd05d6ba7548f6a03c474f4469d8e621e96.png

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.