Fitting radial-velocities

In juliet, the radial-velocity model is essentially the same as the one already introduced for the lightcurve in the Lightcurve fitting with juliet tutorial, i.e., in the absence of extra linear terms (see Incorporating linear models), is of the form (see Section 2 of the juliet paper)

\(\mathcal{M}_{i}(t) + \epsilon_i(t)\),

where \(\epsilon_i(t)\) is a noise model for instrument \(i\) (which as for the Lightcurve fitting with juliet tutorial, here we assume is white-gaussian noise — i.e., we assume \(\epsilon_i(t)\sim \mathcal{N}(0,\sqrt{\sigma(t)^2 + \sigma_{w,i}})\), where \(\sigma^2_{w,i}\) is a jitter term added to each instrument — we extend this to gaussian processes in the Incorporating Gaussian Processes tutorial), and \(\mathcal{M}_{i}(t)\) is the deterministic part of the radial-velocity model for the instrument. The form of this deterministic part of the model is given by

\(\mathcal{M}_{i}(t) = \mathcal{K}(t) + \mu_i + Q(t-t_a)^2 + A(t-t_a) + B\).

Here, \(\mathcal{K}(t)\) is a Keplerian model which models the RV perturbations on the star due to the planets orbiting around it, \(\mu_i\) is the RV of the star as measured by instrument \(i\) and the coefficients \(Q, A\) and \(B\) define an additional long-term trend useful for modelling long-period signals in the RVs that might not be well modelled by an additional Keplerian signal — \(t_a\) is just an arbitrary value substracted to the input times for numerical stability of the coefficients (by default \(t_a = 2458460\) — but this can be defined by the user). By default, no long-term trend is incorporated in the models (i.e., \(Q = A = B = 0\)).

RV fits

To showcase the capabilities juliet has for radial-velocity fitting, here we will analyze the radial-velocities of the TOI-141 system (Espinoza et al. (2019)). We already analyzed the transits of this object in the Getting started tutorial; here we use the radial-velocities (RVs) of this system as it was shown that not only the signal of the transiting planet was present in the RVs, but there is also evidence for _another_ planet in the system. We have uploaded the dataset in a juliet-friendly format [here].

Let us first try to find the RV signature of the transiting planet analyzed in the Getting started tutorial in this dataset. From that analysis, the period is \(P = 1.007917 \pm 0.000073\) days and the time-of-transit center is \(t0 = 2458325.5386 \pm 0.0011\). Let us use these as priors for a first fit to the data — let us in turn assume uniform wide priors for the systemic velocities for each instrument \(\mu_i\), jitter terms and RV semi-amplitude; let us also fix the eccentricity to zero for now:

import juliet
priors = {}

# Name of the parameters to be fit:
params = ['P_p1','t0_p1','mu_CORALIE14', \
          'mu_CORALIE07','mu_HARPS','mu_FEROS',\
          'K_p1', 'ecc_p1', 'omega_p1', 'sigma_w_CORALIE14','sigma_w_CORALIE07',\
           'sigma_w_HARPS','sigma_w_FEROS']

# Distributions:
dists = ['normal','normal','uniform', \
         'uniform','uniform','uniform',\
         'uniform','fixed', 'fixed', 'loguniform', 'loguniform',\
         'loguniform', 'loguniform']

# Hyperparameters
hyperps = [[1.007917,0.000073], [2458325.5386,0.0011], [-100,100], \
           [-100,100], [-100,100], [-100,100], \
           [0.,100.], 0., 90., [1e-3, 100.], [1e-3, 100.], \
           [1e-3, 100.], [1e-3, 100.]]

# Populate the priors dictionary:
for param, dist, hyperp in zip(params, dists, hyperps):
    priors[param] = {}
    priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp

dataset = juliet.load(priors = priors, rvfilename='rvs_toi141.dat', out_folder = 'toi141_rvs')
results = dataset.fit(n_live_points = 300)

To plot the data, one can extract the models in an analogous fashion as we did for the Lightcurve fitting with juliet tutorial: we use the results.rv.evaluate() function. As with the results.lc.evaluate() function presented in the Lightcurve fitting with juliet tutorial, the function receives an instrument name and optionally times in which one wants to evaluate the model. Because each of the RV model parts are additive, it is easy to extract, e.g., the systemic-velocity corrected keplerian signal by simply evaluating the model in an arbitrary instrument and substracting the median of the systemic-velocity for that instrument. Let us do this to plot the above defined fit to see how we did — we’ll only plot the HARPS and FEROS data, as the CORALIE data is not very constraining:

import numpy as np
import matplotlib.pyplot as plt

# Plot HARPS and FEROS datasets in the same panel. For this, first select any
# of the two and substract the systematic velocity to get the Keplerian signal.
# Let's do it with FEROS. First generate times on which to evaluate the model:
min_time, max_time = np.min(dataset.times_rv['FEROS'])-30,\
                     np.max(dataset.times_rv['FEROS'])+30

model_times = np.linspace(min_time,max_time,1000)

# Now evaluate the model in those times, and substract the systemic-velocity to
# get the Keplerian signal:
keplerian = results.rv.evaluate('FEROS', t = model_times) - \
            np.median(results.posteriors['posterior_samples']['mu_FEROS'])

# Now plot the (systematic-velocity corrected) RVs:
fig = plt.figure(figsize=(12,5))
instruments = ['FEROS','HARPS']
colors = ['cornflowerblue','orangered']
for i in range(len(instruments)):
    instrument = instruments[i]
    # Evaluate the median jitter for the instrument:
    jitter = np.median(results.posteriors['posterior_samples']['sigma_w_'+instrument])
    # Evaluate the median systemic-velocity:
    mu = np.median(results.posteriors['posterior_samples']['mu_'+instrument])
    # Plot original data with original errorbars:
    plt.errorbar(dataset.times_rv[instrument]-2457000,dataset.data_rv[instrument]-mu,\
                 yerr = dataset.errors_rv[instrument],fmt='o',\
                 mec=colors[i], ecolor=colors[i], elinewidth=3, mfc = 'white', \
                 ms = 7, label=instrument, zorder=10)

    # Plot original errorbars + jitter (added in quadrature):
    plt.errorbar(dataset.times_rv[instrument]-2457000,dataset.data_rv[instrument]-mu,\
                 yerr = np.sqrt(dataset.errors_rv[instrument]**2+jitter**2),fmt='o',\
                 mec=colors[i], ecolor=colors[i], mfc = 'white', label=instrument,\
                 alpha = 0.5, zorder=5)

# Plot Keplerian model:
plt.plot(model_times-2457000, keplerian,color='black',zorder=1)
plt.ylabel('RV (m/s)')
plt.xlabel('Time (BJD - 2457000)')
plt.title('1 Planet Fit | Log-evidence: {0:.3f} $\pm$ {1:.3f}'.format(results.posteriors['lnZ'],\
       results.posteriors['lnZerr']))
plt.ylim([-20,20])
plt.xlim([1365,1435])

Interesting. We have plotted both the original data with the original errorbars, and the errorbars enlarged by the best-fit jitter term. Note how the jitter is large (specially for HARPS)? This is to explain the large variations that appear in this 1-planet-fit result. Could this be due to an additional planet? To test this hypothesis, let’s try another fit but now fitting for two planets: the 1-day transiting one, and an additional one with an unknown period from, say, 1 to 10 days. To do this, add the extra priors for this model first:

# Add second planet to the prior:
params = params + ['P_p2',   't0_p2',  'K_p2',    'ecc_p2','omega_p2']
dists = dists +   ['uniform','uniform','uniform', 'fixed', 'fixed']
hyperps = hyperps + [[1.,10.],[2458325.,2458330.],[0.,100.], 0., 90.]

# Repopulate priors dictionary:
priors = {}

for param, dist, hyperp in zip(params, dists, hyperps):
    priors[param] = {}
    priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp

And let’s perform the second juliet fit with this two-planet system:

dataset = juliet.load(priors = priors, rvfilename='rvs_toi141.dat', out_folder = 'toi141_rvs_2planets')
results2 = dataset.fit(n_live_points = 300)

Repeating the same plot as above we find:

Woah! Much better fit to the data. Note also that we have plotted the log-evidences that juliet gives for these models — and the log-evidence for the 2-planet model is much larger than the one for the 1-planet model, \(\Delta \ln Z = 114.4\) which is a huge odds ratio in favor of the two-planet model. Let’s plot the posterior distributions for the parameters of this fit using Daniel Foreman-Mackey’s corner package:

import corner

posterior_names = [r"$K_1$ (m/s)", r"$P_2$ (days)", r"$K_2$ (m/s)"]
first_time = True
for i in range(len(params)):
    if dists[i] != 'fixed' and params[i] != 'P_p1' and 't0' not in params[i] and \
    params[i][0:2] != 'mu' and params[i][0:5] != 'sigma':
        if first_time:
            posterior_data = results2.posteriors['posterior_samples'][params[i]]
            first_time = False
        else:
            posterior_data  = np.vstack((posterior_data, results2.posteriors['posterior_samples'][params[i]]))
posterior_data = posterior_data.T
figure = corner.corner(posterior_data, labels = posterior_names)

Best-fit period of this second planet is at 4.76 days — this is slightly off with the value cited in the paper (which is \(4.78503 \pm 0.0005\)), we will touch on this “mistery” in the Joint transit and radial-velocity fits tutorial. The semi-amplitudes mostly agree with the values in the paper. Judging from the errorbars, it seems there still is some unexplained variance in the data. Could it be an additional planet? Let us try fitting an extra planet — this time we will try a larger prior for the period of this third signal, going all the way from 1 to 40 days, which is about half the observing window for the FEROS and HARPS observations, which are the most constraining ones:

# Add third planet to the prior:
params3pl = params + ['P_p3',   't0_p3',  'K_p3',    'ecc_p3','omega_p3']
dists3pl = dists +   ['uniform','uniform','uniform', 'fixed', 'fixed']
hyperps3pl = hyperps + [[1.,40.],[2458325.,2458355.],[0.,100.], 0., 90.]

# Repopulate priors dictionary:
priors3pl = {}

for param, dist, hyperp in zip(params3pl, dists3pl, hyperps3pl):
    priors3pl[param] = {}
    priors3pl[param]['distribution'], priors3pl[param]['hyperparameters'] = dist, hyperp

# Run juliet:
dataset = juliet.load(priors = priors3pl, rvfilename='rvs_toi141.dat', out_folder = 'toi141_rvs_3planets')
results = dataset.fit(n_live_points = 300)

The resulting fit doesn’t look too different from the 2-planet one:

keplerian = results.rv.evaluate('FEROS', t = model_times) - \
            np.median(results.posteriors['posterior_samples']['mu_FEROS'])

# Now plot the (systematic-velocity corrected) RVs:
instruments = ['FEROS','HARPS']
colors = ['cornflowerblue','orangered']
fig = plt.figure(figsize=(12,5))
for i in range(len(instruments)):
    instrument = instruments[i]
    jitter = np.median(results.posteriors['posterior_samples']['sigma_w_'+instrument])
    mu = np.median(results.posteriors['posterior_samples']['mu_'+instrument])
    # Plot original errorbars:
    plt.errorbar(dataset.times_rv[instrument]-2457000,dataset.data_rv[instrument]-mu,\
                 yerr = dataset.errors_rv[instrument],fmt='o',\
                 mec=colors[i], ecolor=colors[i], elinewidth=3, mfc = 'white', \
                 ms = 7, label=instrument, zorder=10)
    # Plot original errorbars + jitter:
    plt.errorbar(dataset.times_rv[instrument]-2457000,dataset.data_rv[instrument]-mu,\
                 yerr = np.sqrt(dataset.errors_rv[instrument]**2+jitter**2),fmt='o',\
                 mec=colors[i], ecolor=colors[i], mfc = 'white', label=None,\
                 alpha = 0.5, zorder=5)

plt.plot(model_times-2457000, keplerian,color='black',zorder=1)
plt.ylabel('RV (m/s)')
plt.xlabel('Time (BJD - 2457000)')
plt.title('3 Planet Fit | Log-evidence: {0:.3f} $\pm$ {1:.3f}'.format(results.posteriors['lnZ'],\
           results.posteriors['lnZerr']))
plt.ylim([-20,20])
plt.xlim([1365,1435])
plt.legend()

In fact, the evidence is worse in this 3-planet fit (\(\ln Z_3 = -694\)) than in the 2-planet fit (\(\ln Z_2 = -691\)). If both models were equiprobable a-priori, these log-evidences mean that, given the data, the 2-planet model is about 20 times more likely than the 3-planet model. So it seems that if there is some extra variance in the dataset, given the data at hand, this cannot be explained by an extra, third planetary signal alone — at least not with periods between 1 and 40 days. But what if there is a longer period planet creating a trend in the data? We deal with this possibility next

Long-term trends in RV data

As mentioned above, within juliet it is possible to fit for a long-term trend in the data that is common to all the instruments, parametrized by an intercept \(B\) (rv_intercept parameter within juliet), a slope \(A\) (rv_slope parameter within juliet) and a quadratic coefficient \(Q\) (rv_quad parameter within juliet). This long-term trend is useful to constrain signals whose periods might be longer than the current time baseline, which might locally appear as long-term trends. To fit those to the data, we just need to define priors for these parameters — let us do this with the TOI-141 dataset by first trying to fit a simple linear term (i.e., let us define only the parameters rv_intercept and rv_slope). Let us give wide uniform priors for those, join those priors to the 2-planet-fit priors and perform the fit:

# Add linear trend to the prior:
paramsLT = params + ['rv_intercept',   'rv_slope']
distsLT = dists +   ['uniform','uniform']
hyperpsLT = hyperps + [[-100.,100.],[-100., 100.]]

# Repopulate priors dictionary:
priorsLT = {}

for param, dist, hyperp in zip(paramsLT, distsLT, hyperpsLT):
    priorsLT[param] = {}
    priorsLT[param]['distribution'], priorsLT[param]['hyperparameters'] = dist, hyperp

# Run juliet:
dataset = juliet.load(priors = priorsLT, rvfilename='rvs_toi141.dat', out_folder = 'toi141_rvs_lineartrend')
results = dataset.fit(n_live_points = 300)

Before plotting the results, note that when we evaluate the model using results.rv.evaluate we will get back the full model — that is, a Keplerian plus the long-term trend model in our case (plus the systemic velocity of the instrument). However, one can pass an extra flag to this function, the return_components flag, which in addition to the full model returns a dictionary that will have all the (deterministic) components of the model. Let us plot all the components of the model on top of each other using this flag:

# Return full model and the components of the model:
full_model, components = results.rv.evaluate('FEROS', t = model_times, return_components = True)
# Substract systemic RV from full model (note this is part of the components):
full_model -= components['mu']

# Now plot the (systematic-velocity corrected) RVs (same code as above):
instruments = ['FEROS','HARPS']
colors = ['cornflowerblue','orangered']
fig = plt.figure(figsize=(12,5))
for i in range(len(instruments)):
    instrument = instruments[i]
    jitter = np.median(results.posteriors['posterior_samples']['sigma_w_'+instrument])
    mu = np.median(results.posteriors['posterior_samples']['mu_'+instrument])
    # Plot original errorbars:
    plt.errorbar(dataset.times_rv[instrument]-2457000,dataset.data_rv[instrument]-mu,\
                 yerr = dataset.errors_rv[instrument],fmt='o',\
                 mec=colors[i], ecolor=colors[i], elinewidth=3, mfc = 'white', \
                 ms = 7, label=instrument, zorder=10)
    # Plot original errorbars + jitter:
    plt.errorbar(dataset.times_rv[instrument]-2457000,dataset.data_rv[instrument]-mu,\
                 yerr = np.sqrt(dataset.errors_rv[instrument]**2+jitter**2),fmt='o',\
                 mec=colors[i], ecolor=colors[i], mfc = 'white', label=None,\
                 alpha = 0.5, zorder=5)

# Plot full model:
plt.plot(model_times-2457000, full_model,color='black',zorder=1, label = 'Full model')

# Extract model components and plot them:
plt.plot(model_times-2457000, components['keplerian'],color='grey',zorder=0, alpha=0.5, label = 'Keplerian')
plt.plot(model_times-2457000, components['trend'],color='grey',zorder=0,alpha=0.5, lw = 3, label = 'Linear trend')

# Labels:
plt.ylabel('RV (m/s)')
plt.xlabel('Time (BJD - 2457000)')
plt.title('2 Planet Fit + Linear Trend | Log-evidence: {0:.3f} $\pm$ {1:.3f}'.format(results.posteriors['lnZ'],\
           results.posteriors['lnZerr']))
plt.ylim([-35,25])
plt.xlim([1365,1435])
plt.legend(ncol = 2)

As can be seen, the components dictionary extracted from the results.rv.evaluate function contains the Keplerian signal under components['keplerian'], and the trend under components['keplerian']. In addition, it also stores the Keplerians of each of the individual planets under components['p1'] and components['p2'] in our case. Note however, that the linear trend appears to not be significant in our case. So it might be that the unexplained variance could be explained by something else — in the Incorporating Gaussian Processes tutorial, we explore adding a Gaussian Process to the dataset in order to explain this.

Note

Note how in our case the components dictionary for the FEROS instrument has its systemic RV stored under components['mu'], which in general is different than taking the median of the results.posteriors['posterior_samples']['mu_FEROS'] array. This is because, as was already mentioned in the Lightcurve fitting with juliet tutorial, the results.rv.evaluate function (and the results.lc.evaluate function) evaluate the model by default on nsamples = 1000 samples of the posterior. Thus, components['mu'] is the median value of the systemic RV over the same 1000 samples as the other components, whereas results.posteriors['posterior_samples']['mu_FEROS'] contains all the samples and thus, taking the median of that array should be slightly different than components['mu']. This difference, of course, is typically much smaller than the errors, so it shouldn’t be a problem in general. One can set the all_samples flag to True in the results.rv.evaluate function to use all the samples — in this case, both should give the same results.