© ESO 2022

1. Introduction

It is common in astronomy to indirectly detect a physical event or the presence of a body by searching for its signature in a data set and, more specifically, in a time series. Astronomical time series are typically corrupted by photon noise, which is uncorrelated: the noise values at two distinct times are statistically independent. In that case, as more data are acquired, the searched-for signal should emerge more clearly. However, in many cases, the data are also corrupted by correlated noise emerging from other physical events, contamination from the Earth’s atmosphere, instrumental noise, etc. In some cases, the structure of this correlated noise can mimic the signal of interest, leading to spurious detections or to a poor estimation of the model parameters.

This situation is encountered in particular when searching for exoplanets in radial velocity (RV) data. The RV of a star is the star velocity projected onto the line of sight, measured thanks to the Doppler effect. The presence of a planetary companion induces a reflex motion of the star and thus a periodic pattern in the RV time series. The latest generation of spectrographs, such as ESPRESSO (Pepe et al. 2021) or EXPRES (Blackman et al. 2020), is able to reach a RV precision of the order of 10 cm s⁻¹, and has the potential to discover Earth-like planets around Sun-like stars. However, correlated noise of stellar origin complexifies this task. The p-modes and granulation processes introduce correlated noise at different timescales (Dumusque et al. 2011, 2012). Furthermore, the random appearance of spots and faculae at the surface of the star, combined with the star’s rotation, introduces complex structure in the data, which might be difficult to disentangle from low-mass planets. At longer timescales (from hundreds to thousands of days), the stellar magnetic cycle also induces RV variations, as well as variations in activity indicators such as the flux in the calcium II H & K emission lines ($log{R}_{\text{HK}}^{\prime }$, Noyes et al. 1984) as in, for instance, (Queloz et al. 2001). The RV signals induced by these multiple physical processes of stellar origin are globally referred to as stellar activity. A common way to account for stellar activity is to model it through a Gaussian process (GP) model (Rasmussen & Williams 2006). The GP regression method allows for the modeling of complex processes by parametrizing the covariance between the measurements instead of defining a deterministic model of the physical processes. For a GP G(t) measured at times t_i and t_j, the values G(t_i) and G(t_j) are assumed to be randomly drawn from a normal distribution, with covariance C_i, j = k(t_i, t_j), where k is the chosen parametrized kernel. The GP is often assumed to be stationary, such that the kernel k only depends on the lag Δt = |t_i − t_j| between two measurements. A commonly used kernel to model stellar activity is (Aigrain et al. 2012; Haywood et al. 2014; Rajpaul et al. 2015):

$$\begin{array}{c}\hfill k(\mathrm{\Delta }t)={\sigma }^{2}exp(-\frac{\mathrm{\Delta }{t}^2}{2{\rho }^2}-\frac{{sin}^2\left(\frac{\pi \mathrm{\Delta }t}{P}\right)}{2{\eta }^2}),\end{array}$$(1)

where σ, P, ρ, and η are the hyperparameters of the GP which need to be adjusted. In the following, we refer to this kernel as the squared-exponential periodic (SEP) kernel.

Gaussian processes are known to be able to represent a wide range of signals. As such, when their hyperparameters are left free, they are prone to absorb the signal of interest (planetary signal) along with the correlated noise (stellar activity). To avoid this drawback, Rajpaul et al. (2015) proposed a framework in which the RV time series is modeled jointly with activity indicators. Building on Aigrain et al. (2012), the authors assume that the activity-induced variations of the measurements depend linearly on an underlying Gaussian process G(t) and its derivative G′(t). The evolution of the RV and indicators is modeled as:

$$\begin{array}{cc}& \mathrm{\Delta }\mathrm{RV}=V_{c}G(t)+V_r{G}^{\prime }(t),\hfill \\ & \mathrm{\Delta }\mathrm{BIS}=B_{c}G(t)+B_r{G}^{\prime }(t),\hfill \\ & \mathrm{\Delta }log{R}_{\mathrm{HK}}^{\prime }=L_{c}G(t),\hfill \end{array}$$(2)

for some constants V_c, V_r, B_c, B_r, L_c. The GP’s hyperparameters are thus constrained by the three time series, instead of the RV alone. This reduces the risk of the GP overfitting, that is, the absorbtion of planetary signals, since those signals are only present in the RV time series. This framework can be straightforwardly generalized to account for additional indicators, for the combination of several GP with different amplitudes, or even for the second order derivatives of the GP (e.g., Jones et al. 2017).

While this framework is very powerful in modeling stellar activity, it represents a challenge in terms of computational cost. Indeed, computing the likelihood (or χ²) of the model for a given set of hyperparameters requires us to solve a linear system involving the full covariance matrix of the measurements. For a time series of size n, the full covariance matrix – including RV, BIS, and $log{R}_{\text{HK}}^{\prime }$ measurements – has a size of 3n × 3n, and the cost of the solving typically scales as 𝒪((3n)³). This becomes prohibitive in terms of computer time for large data sets, especially in the context of Bayesian methods (MCMC, nested sampling, etc.), which might require billions of evaluations of the likelihood.

In the context of a GP applied to a single time series, Ambikasaran (2015) and Foreman-Mackey (2018) (see also Rybicki & Press 1995) have shown that the so-called celerite kernel,

$$\begin{array}{c}\hfill k(\mathrm{\Delta }t)={\displaystyle \sum _{s<n_{\mathrm{c}}}(a_{s}cos({\nu }_s\mathrm{\Delta }t)+b_{s}sin({\nu }_s\mathrm{\Delta }t)){\mathit{e}}^{-{\lambda }_s\mathrm{\Delta }t},}\end{array}$$(3)

where n_c is the number of components, and a_s, b_s, λ_s, and ν_s are the kernel hyperparameters, can be represented as a semiseparable matrix. As a consequence, the computational cost of evaluating the likelihood scales linearly with the number of points (𝒪(n)), allowing to apply these methods to large data sets. Delisle et al. (2020b) defined a more general class of covariance matrices with a similar linear scaling of the cost: the S+LEAF matrix, which is the sum of a semiseparable matrix and a LEAF matrix. The LEAF component, which has non-zero elements close to the diagonal, is particularly adapted to represent calibration noise (see Delisle et al. 2020b). Gordon et al. (2020) extended the celerite model to the case of two-dimensional data sets. This applies, in particular, to the case of several parallel time series (e.g., RV, BIS, and $log{R}_{\text{HK}}^{\prime }$), with measurements taken at the same times. However, Gordon et al. (2020) do not discuss the treatment of the derivatives of the GP, and thus of models similar to the Rajpaul et al. (2015) model (Eq. (2)).

In this study, we extend the celerite and S+LEAF models to account for the case of several time series, with independent calendars, modeled as a linear combination of several GP and their derivatives. This allows us to apply models similar to the model used by Rajpaul et al. (2015), but with a linear scaling of the evaluation cost of the likelihood. We call this new model S+LEAF 2, as it is a generalization of the S+LEAF model (Delisle et al. 2020b).

In Sect. 2, we recall the main properties of the celerite and S+LEAF models. We then show in Sect. 3 how to model the derivative of a celerite GP. In Sect. 4, we extend the model to the case of multiple time series. We illustrate the power of this framework by reanalyzing the HARPS RV of the nearby K2 dwarf HD 13808 in Sect. 5. Finally, we discuss our methods and results in Sect. 6. An open-source reference implementation of our algorithms as C library with python wrappers is publicly available¹.

2. The celerite and S+LEAF models for homogeneous time series

We consider a time series of measurements (t_i, y_i) (i = 1, …, n), which can be modeled by a deterministic component, a GP component, and measurement noise. In the case of radial velocities, the deterministic component encompasses the reflex motion due to companions, the systematic velocity of the system, instruments offsets, and so on. The GP might be used to model physical mechanisms that are too poorly understood or constrained to be included in the deterministic part. This is typically the case of stellar activity at different timescales (oscillations, granulation, rotation, magnetic cycles, etc.). Finally the noise component encompasses photon noise, calibration noise, and so on. The time series can thus be expressed as:

$$\begin{array}{c}\hfill y_i=m(t_i)+G(t_i)+{ϵ}_i,\end{array}$$(4)

where m is the deterministic part of the model, G is the GP, and ϵ the noise. These three components of the model might depend on a set of parameters θ. Assuming the noise to also be Gaussian (but not necessarily white), the log-likelihood function of a given set of parameters θ is

$$\begin{array}{cc}\hfill ln\mathcal{L}(\theta )& =lnp(y|\theta )\hfill \\ \hfill & =-\frac{1}{2}{(y-m_{\theta })}^{\mathrm{T}}{C}_{\theta }^{-1}(y-m_{\theta })\hfill \\ \hfill & -\frac{1}{2}lndet\left(2\pi C_{\theta }\right).\hfill \end{array}$$(5)

C is the total covariance matrix of the time series which can be split as

$$\begin{array}{c}\hfill C=K+\mathrm{\Sigma },\end{array}$$(6)

where K is the covariance of the GP G and Σ is the covariance of the noise

$$\begin{array}{cc}& K_{i,j}=\text{cov}(G(t_i),G(t_j))=k(t_i,t_j),\hfill \\ & {\mathrm{\Sigma }}_{i,j}=\text{cov}({ϵ}_i,{ϵ}_j),\hfill \end{array}$$(7)

with k the GP kernel function.

2.1. The celerite model

The celerite model proposed by Foreman-Mackey et al. (2017) allows for very efficient computations of this model and, in particular, the likelihood and its gradient with respect to θ, in the case of white noise (diagonal covariance matrix Σ = diag(σ²)) and assuming the kernel function k to follow Eq. (3). In this case, the covariance matrix C is semiseparable (see Foreman-Mackey et al. 2017):

$$\begin{array}{c}\hfill C=\text{diag}(A+{\sigma }^2)+\text{tril}\left(U{V}^{\mathrm{T}}\right)+\text{triu}\left(V{U}^{\mathrm{T}}\right),\end{array}$$(8)

where diag(A) is the diagonal matrix built from the vector A of size n, U and V are n × r matrices (with r = 2n_c the rank of U and V) defined as:

$$\begin{array}{cc}& A_i={\displaystyle \sum _{s<n_{\mathrm{c}}}{a}_s,}\hfill \\ & U_{i,s}={\mathit{e}}^{-{\lambda }_s{t}_i}(a_{s}cos({\nu }_s{t}_i)+b_{s}sin({\nu }_s{t}_i)),\hfill \\ & U_{i,n_{\mathrm{c}}+s}={\mathit{e}}^{-{\lambda }_s{t}_i}(a_{s}sin({\nu }_s{t}_i)-b_{s}cos({\nu }_s{t}_i)),\hfill \\ & V_{i,s}={\mathit{e}}^{{\lambda }_s{t}_i}cos({\nu }_s{t}_i),\hfill \\ & V_{i,n_{\mathrm{c}}+s}={\mathit{e}}^{{\lambda }_s{t}_i}sin({\nu }_s{t}_i).\hfill \end{array}$$(9)

This low-rank representation of the covariance matrix allows us to use very efficient dedicated algorithms (Cholesky decomposition, solving, dot product, determinant, see Foreman-Mackey et al. 2017; Foreman-Mackey 2018), which are particularly useful in computing the likelihood. The memory footprint of the celerite model scales as 𝒪(nr), and the computational cost scales as 𝒪(nr²). For comparison, the memory footprint of a naive representation of the same covariance matrix scales as 𝒪(n²) and the computational costs typically scales as 𝒪(n³).

The celerite algorithm is actually not restricted to kernels of the form of Eq. (3), but it can be applied as long as the kernel admits a semiseparable representation following Eq. (8). In particular, the Matérn 3/2 and Matérn 5/2 kernels can be represented using a semiseparable representation of rank 2 and 3 respectively (see Appendix A.1).

2.2. The S+LEAF model

The S+LEAF model (Delisle et al. 2020b) extends the celerite model in the case of correlated measurement noise. While in the celerite model, the covariance matrix Σ is assumed to be diagonal, the S+LEAF model allows us to account for a more general class of noise models, where Σ can be represented by a LEAF matrix. Such LEAF matrices are sparse matrices where non-zero values are close to the diagonal (see Delisle et al. 2020b). This encompasses banded matrices, block-diagonal matrices, staircase matrices, and so on. In the context of radial velocities, LEAF matrices are especially useful to model calibration noise (see Delisle et al. 2020b). Indeed, to obtain precise RV measurements, the instrument must be calibrated periodically, typically once per night. Thus, several measurements taken with the same instrument during the same night share the same calibration noise, which introduces a block-diagonal contribution to the covariance matrix.

As in the case of the celerite model, the S+LEAF model allows for a sparse representation of the covariance matrix and for efficient dedicated algorithms to compute the likelihood and its gradient with respect to the parameters (Delisle et al. 2020b). The memory footprint of the S+LEAF model scales as $\mathcal{O}\left((r+\overline{b})n\right)$ and its computational cost as $\mathcal{O}\left((r^2+r\overline{b}+{\overline{b}}^2)n\right)$, where $\overline{b}$ is the average band width of the LEAF component.

3. Derivative of a celerite/S+LEAF Gaussian process

We consider a GP G(t) whose kernel function k is stationary and admits a semiseparable decomposition (Eq. (8)). The covariance matrix of G is thus given by:

$$\begin{array}{c}\hfill \text{cov}(G(t_i),G(t_j))=k(t_i-t_j)=U(t_i)V(t_j)=U_i{V}_j,\end{array}$$(10)

for t_i ≥ t_j, namely, for the lower triangular part. By way of symmetry, we have

$$\begin{array}{c}\hfill \text{cov}(G(t_i),G(t_j))=k(t_j-t_i)=V_i{U}_j,\end{array}$$(11)

for t_i < t_j, namely, for the upper triangular part. Assuming G to be differentiable, the covariance between G and G′ is given by the partial derivatives of the kernel function:

$$\begin{array}{cc}& \text{cov}(G^{\prime }(t_i),G(t_j))=\frac{\partial k(t_i,t_j)}{\partial t_i}=k^{\prime }(t_i-t_j)=U^{\prime }(t_i)V(t_j)=U_i^{\prime }{V}_j,\hfill \\ & \text{cov}(G(t_i),G^{\prime }(t_j))=\frac{\partial k(t_i,t_j)}{\partial t_j}=-k^{\prime }(t_i-t_j)=U_i{V}_j^{\prime },\hfill \end{array}$$(12)

where we assume t_i ≥ t_j, and the primes denote the differentiation with respect to time. For these two equations to be valid, the semiseparable representation (U, V) must verify:

$$\begin{array}{c}\hfill U\prime V^{\mathrm{T}}=-U{V}^{\prime \mathrm{T}}.\end{array}$$(13)

This relation is actually linked with the stationarity of the kernel. Indeed, since

$$\begin{array}{c}\hfill k(t,t+\mathrm{\Delta }t)=k(\mathrm{\Delta }t)=U(t)V(t+\mathrm{\Delta }t)\end{array}$$(14)

does not depend on t, we have

$$\begin{array}{c}\hfill \frac{\partial k(t,t+\mathrm{\Delta }t)}{\partial t}=0=U\prime (t)V(t+\mathrm{\Delta }t)+U(t)V\prime (t+\mathrm{\Delta }t).\end{array}$$(15)

In addition to this stationarity condition, for the GP G to be differentiable the kernel must verify k′(0) = 0 (e.g., Rasmussen & Williams 2006), thus:

$$\begin{array}{c}\hfill U_i^{\prime }{V}_i=U_i{V}_{\prime }^i=0.\end{array}$$(16)

From Eq. (12), we deduce that the covariance matrix between G and G′ admits several equivalent antisymmetric semiseparable representations

$$\begin{array}{cc}\hfill \text{cov}(G^{\prime },G)& =\text{tril}\left(U^{\prime }{V}^{\mathrm{T}}\right)-\text{triu}\left(V{U^{\prime }}^{\mathrm{T}}\right)\hfill \\ \hfill & =-\text{tril}\left(U{V^{\prime }}^{\mathrm{T}}\right)+\text{triu}\left(V^{\prime }{U}^{\mathrm{T}}\right),\hfill \\ \hfill & =\text{tril}\left(U^{\prime }{V}^{\mathrm{T}}\right)+\text{triu}\left(V^{\prime }{U}^{\mathrm{T}}\right),\hfill \\ \hfill & \dots \hfill \end{array}$$(17)

Similar representations can be deduced for cov(G,G′) by using the relation

$$\begin{array}{c}\hfill \text{cov}(G,G^{\prime })=\text{cov}{(G^{\prime },G)}^{\mathrm{T}}=-\text{cov}(G^{\prime },G).\end{array}$$(18)

For the covariance matrix of G′ itself, we have (for t_i ≥ t_j)

$$\begin{array}{cc}\hfill \text{cov}(G^{\prime }(t_i),G^{\prime }(t_j))& =\frac{{\partial }^{2}k(t_i,t_j)}{\partial t_i\partial t_j}=-k^{″}(t_i-t_j)\hfill \\ \hfill & =U_i^{\prime }{V}_j^{\prime }=-U_i^{″}{V}_j=-U_i{V}_j^{″}.\hfill \end{array}$$(19)

Therefore, the covariance matrix of G′ also admits several symmetric semiseparable representations:

$$\begin{array}{cc}\hfill \text{cov}(G^{\prime },G^{\prime })& =\text{diag}(B)+\text{tril}\left(U^{\prime }{V^{\prime }}^{\mathrm{T}}\right)+\text{triu}\left(V^{\prime }{U^{\prime }}^{\mathrm{T}}\right)\hfill \\ \hfill & =\text{diag}(B)-\text{tril}(U″V^{\mathrm{T}})-\text{triu}\left(V{U^{″}}^{\mathrm{T}}\right)\hfill \\ \hfill & =\text{diag}(B)-\text{tril}\left(U{V″}^{\mathrm{T}}\right)-\text{triu}\left(V^{″}{U}^{\mathrm{T}}\right)\hfill \\ \hfill & \dots \hfill \end{array}$$(20)

with

$$\begin{array}{c}\hfill B_i=U_i^{\prime }{V}_i^{\prime }=-U_i^{\prime }{V}_i=-U_i{V}_{″}^i.\end{array}$$(21)

Hereinafter, we use the following representations:

$$\begin{array}{cc}\hfill \text{cov}(G^{\prime },G)& =\text{tril}\left(U^{\prime }{V}^{\mathrm{T}}\right)+\text{triu}\left(V^{\prime }{U}^{\mathrm{T}}\right),\hfill \\ \hfill \text{cov}(G,G^{\prime })& =\text{tril}\left(U{V^{\prime }}^{\mathrm{T}}\right)+\text{triu}\left(V{U^{\prime }}^{\mathrm{T}}\right),\hfill \\ \hfill \text{cov}(G^{\prime },G^{\prime })& =\text{diag}(B)+\text{tril}\left(U^{\prime }{V^{\prime }}^{\mathrm{T}}\right)+\text{triu}\left(V^{\prime }{U^{\prime }}^{\mathrm{T}}\right).\hfill \end{array}$$(22)

Applying this reasoning to the celerite kernel of Eq. (3), we find (as per Eq. (9)):

$$\begin{array}{cc}\hfill k^{\prime }(\mathrm{\Delta }t)& ={\displaystyle \sum _{s<n_{\mathrm{c}}}(a_s^{\prime }cos({\nu }_s\mathrm{\Delta }t)+b_{\prime }^{s}sin({\nu }_s\mathrm{\Delta }t)){\mathit{e}}^{-{\lambda }_s\mathrm{\Delta }t},}\hfill \\ \hfill -k^{″}(\mathrm{\Delta }t)& ={\displaystyle \sum _{s<n_{\mathrm{c}}}(a_s^{″}cos({\nu }_s\mathrm{\Delta }t)+b_s^{″}sin({\nu }_s\mathrm{\Delta }t)){\mathit{e}}^{-{\lambda }_s\mathrm{\Delta }t},}\hfill \\ \hfill U_{i,s}^{\prime }& ={\mathit{e}}^{-{\lambda }_s{t}_i}(a_s^{\prime }cos({\nu }_s{t}_i)+b_{\prime }^{s}sin({\nu }_s{t}_i)),\hfill \\ \hfill U_{i,n_{\mathrm{c}}+s}^{\prime }& ={\mathit{e}}^{-{\lambda }_s{t}_i}(a_s^{\prime }sin({\nu }_s{t}_i)-b_{\prime }^{s}cos({\nu }_s{t}_i)),\hfill \\ \hfill V_{i,s}^{\prime }& ={\mathit{e}}^{{\lambda }_s{t}_i}({\lambda }_{s}cos({\nu }_s{t}_i)-{\nu }_{s}sin({\nu }_s{t}_i)),\hfill \\ \hfill V_{i,n_{\mathrm{c}}+s}^{\prime }& ={\mathit{e}}^{{\lambda }_s{t}_i}({\lambda }_{s}sin({\nu }_s{t}_i)+{\nu }_{s}cos({\nu }_s{t}_i)),\hfill \end{array}$$(23)

with

$$\begin{array}{cc}\hfill a_s^{\prime }& ={\nu }_s{b}_s-{\lambda }_s{a}_s,\hfill \\ \hfill b_{\prime }^s& =-{\nu }_s{a}_s-{\lambda }_s{b}_s,\hfill \\ \hfill a_s^{″}& =({\nu }_s^2-{\lambda }_s^2)a_s+2{\lambda }_s{\nu }_s{b}_s,\hfill \\ \hfill b_s^{″}& =({\nu }_s^2-{\lambda }_s^2)b_s-2{\lambda }_s{\nu }_s{a}_s.\hfill \end{array}$$(24)

In this case, the differentiability condition of Eq. (16) can be rewritten as

$$\begin{array}{c}\hfill k^{\prime }(0)=0={\displaystyle \sum _{s<n_{\mathrm{c}}}{a}_s^{\prime }.}\end{array}$$(25)

Thus the initial parameters must verify

$$\begin{array}{c}\hfill {\displaystyle \sum _{s<n_{\mathrm{c}}}{\nu }_s{b}_s-{\lambda }_s{a}_s=0,}\end{array}$$(26)

for the GP to be differentiable. In particular, in the case of a kernel including a single celerite component (n_c = 1), one must verify $a_s^{\prime }={\nu }_s{b}_s-{\lambda }_s{a}_s=0$ and Eq. (23) is simplified into

$$\begin{array}{cc}\hfill k^{\prime }(\mathrm{\Delta }t)& =b_s^{\prime }sin({\nu }_s\mathrm{\Delta }t){\mathit{e}}^{-{\lambda }_s\mathrm{\Delta }t},\hfill \\ \hfill U_{i,s}^{\prime }& =b_s^{\prime }{\mathit{e}}^{-{\lambda }_s{t}_i}sin({\nu }_s{t}_i),\hfill \\ \hfill U_{i,n_{\mathrm{c}}+s}^{\prime }& =-b_s^{\prime }{\mathit{e}}^{-{\lambda }_s{t}_i}cos({\nu }_s{t}_i).\hfill \end{array}$$(27)

The SHO kernel, which is a particular type of celerite kernel that only depends on three free parameters, always verifies this differentiability condition. It thus provides a smoother GP than the general celerite kernel, which is why Foreman-Mackey et al. (2017) recommended its application.

In the case of the Matérn 3/2 and 5/2 kernels, similar semiseparable decompositions can be achieved for the derivatives (see Appendix A.2). More generally, Eq. (22) provides the semiseparable decomposition of the derivatives for any differentiable semiseparable kernel.

4. S+LEAF 2: Extending celerite/S+LEAF to heterogeneous time series

Following Rajpaul et al. (2015), we assume that the time series of the radial velocities and the different indicators follow

$$\begin{array}{c}\hfill Y_{i,j}=f_i(T_{i,j})+{\displaystyle \sum _k({\alpha }_{k,i}{G}_k(T_{i,j})+{\beta }_{k,i}{G}_{\prime }^k(T_{i,j}))+{ϵ}_{i,j},}\end{array}$$(28)

where (T_1, ., Y_1, .) is the RV time series and (T_i, ., Y_i, .) are the indicators time series (i > 1), f_i is the determinist part of the model for the time series i, G_k are independent GP, and ϵ is the measurement noise (including photon noise, calibration noise, etc.).

The times and number of measurements need not be the same for all time series, which implies that T and Y are not necessarily matrices but collections of vectors of variable length. In the case of the model presented in Eq. (2), the activity indicators (BIS and $log{R}_{\text{HK}}^{\prime }$) are typically extracted from the same spectra as the RV time series and, thus, they share the same sampling. However, activity indicators can be extracted from other instruments or techniques. For instance, Haywood et al. (2014) trained a GP on the CoRoT light curve of CoRoT 7 to then use it to model the impact of stellar activity on the HARPS radial velocity time series of the same star. While both data sets were roughly contemporary, the time series did not share the same sampling. In such a case, we could define (T_1, ., Y_1, .) as the RV time series, and (T_2, ., Y_2, .) as the photometric time series (with T₁ ≠ T₂), and model both time series jointly according to Eq. (28).

For the sake of readability, we consider in the following a single GP G, such that

$$\begin{array}{c}\hfill Y_{i,j}=f_i(T_{i,j})+{\alpha }_{i}G(T_{i,j})+{\beta }_i{G}^{\prime }(T_{i,j})+{ϵ}_{i,j},\end{array}$$(29)

but the reasoning holds in the more general case of Eq. (28). The covariance matrix corresponding to Eq. (29) is

$$\begin{array}{cc}\hfill \text{cov}(Y_{i,j},Y_{l,m})=& {\alpha }_i{\alpha }_l\text{cov}(G(T_{i,j}),G(T_{l,m}))\hfill \\ \hfill & +{\beta }_i{\alpha }_l\text{cov}(G^{\prime }(T_{i,j}),G(T_{l,m}))\hfill \\ \hfill & +{\alpha }_i{\beta }_l\text{cov}(G(T_{i,j}),G^{\prime }(T_{l,m}))\hfill \\ \hfill & +{\beta }_i{\beta }_l\text{cov}(G^{\prime }(T_{i,j}),G^{\prime }(T_{l,m}))\hfill \\ \hfill & +\text{cov}({ϵ}_{i,j},{ϵ}_{l,k}).\hfill \end{array}$$(30)

The full covariance matrix is a n × n matrix, where n is the total number of measurements (radial velocities and indicators). The cost of evaluating the corresponding likelihood – which requires to solve a linear system involving the covariance matrix and computing its determinant – typically scales as 𝒪(n³). Therefore, a direct evaluation of the covariance matrix becomes rapidly prohibitive in terms of memory footprint and computing time.

4.1. Semiseparable representation of the model

In order to construct a semiseparable representation of the covariance matrix of Eq. (30) we merge all time vectors T_i into a single vector t, and all data vectors Y_i into a single vector y. We additionally sort the measurements by increasing time; thus, the measurements of the different time series are completely mixed in the merged time series (t, y). For a measurement (t_k, y_k) we denote by ℐ_k (index of the original time series the measurement belongs to) and 𝒥_k (index of the measurement in this original time series) the corresponding couple of indices of this measurement in T, Y. The model of Eq. (29) can be rewritten as

$$\begin{array}{c}\hfill y_k=f_{{\mathcal{I}}_k}(t_k)+{\mathit{\alpha }}_{k}G(t_k)+{\mathit{\beta }}_k{G}^{\prime }(t_k)+{ϵ}_{{\mathcal{I}}_k,{\mathcal{J}}_k},\end{array}$$(31)

with α_k = α_ℐk and β_k = β_ℐk. The corresponding covariance matrix is (as per Eq. (30)):

$$\begin{array}{cc}\hfill C=\text{cov}(y,y)=& \left(\mathit{\alpha }{\mathit{\alpha }}^{\mathrm{T}}\right)\ast \text{cov}(G(t),G(t))\hfill \\ \hfill & +\left(\mathit{\beta }{\mathit{\alpha }}^{\mathrm{T}}\right)\ast \text{cov}(G^{\prime }(t),G(t))\hfill \\ \hfill & +\left(\mathit{\alpha }{\mathit{\beta }}^{\mathrm{T}}\right)\ast \text{cov}(G(t),G^{\prime }(t))\hfill \\ \hfill & +\left(\mathit{\beta }{\mathit{\beta }}^{\mathrm{T}}\right)\ast \text{cov}(G^{\prime }(t),G^{\prime }(t))\hfill \\ \hfill & +\mathrm{\Sigma },\hfill \end{array}$$(32)

where Σ is the covariance matrix of the measurement noise

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{k,l}=\text{cov}({ϵ}_{{\mathcal{I}}_k,{\mathcal{J}}_k},{ϵ}_{{\mathcal{I}}_l,{\mathcal{J}}_l}),\end{array}$$(33)

and X​ * ​Y is the Hadamard (or element-wise) product

$$\begin{array}{c}\hfill {(X\ast Y)}_{i,j}=X_{i,j}{Y}_{i,j}.\end{array}$$(34)

We now assume that the GP G can be modeled by a semiseparable kernel (see Eq. (8)). We can thus use the merged time vector t to compute the semiseparable representation of the covariance matrix of G(t), G′(t) according to Eqs. (8), (9), (22), and (23) (see also Appendix A), and we obtain

$$\begin{array}{cc}\hfill C=& \left(\mathit{\alpha }{\mathit{\alpha }}^{\mathrm{T}}\right)\ast (\text{diag}(A)+\text{tril}\left(U{V}^{\mathrm{T}}\right)+\text{triu}\left(V{U}^{\mathrm{T}}\right))\hfill \\ \hfill & +\left(\mathit{\beta }{\mathit{\alpha }}^{\mathrm{T}}\right)\ast (\text{tril}\left(U^{\prime }{V}^{\mathrm{T}}\right)+\text{triu}\left(V^{\prime }{U}^{\mathrm{T}}\right))\hfill \\ \hfill & +\left(\mathit{\alpha }{\mathit{\beta }}^{\mathrm{T}}\right)\ast (\text{tril}\left(U{V^{\prime }}^{\mathrm{T}}\right)+\text{triu}\left(V{U^{\prime }}^{\mathrm{T}}\right))\hfill \\ \hfill & +\left(\mathit{\beta }{\mathit{\beta }}^{\mathrm{T}}\right)\ast (\text{diag}(B)+\text{tril}\left(U^{\prime }{V^{\prime }}^{\mathrm{T}}\right)+\text{triu}\left(V^{\prime }{U^{\prime }}^{\mathrm{T}}\right))\hfill \\ \hfill & +\mathrm{\Sigma }.\hfill \end{array}$$(35)

The Hadamard product (αβ^T)​ * ​M can also be rewritten as

$$\begin{array}{c}\hfill \left(\mathit{\alpha }{\mathit{\beta }}^{\mathrm{T}}\right)\ast M=\text{diag}(\mathit{\alpha })M\text{diag}(\mathit{\beta }),\end{array}$$(36)

thus

$$\begin{array}{cc}\hfill \left(\mathit{\alpha }{\mathit{\alpha }}^{\mathrm{T}}\right)\ast \text{diag}(A)& =\text{diag}({\mathit{\alpha }}^2\ast A),\hfill \\ \hfill \left(\mathit{\alpha }{\mathit{\beta }}^{\mathrm{T}}\right)\ast \text{tril}\left(U{V}^{\mathrm{T}}\right)& =\text{tril}\left((\mathit{\alpha }\ast U){(\mathit{\beta }\ast V)}^{\mathrm{T}}\right),\hfill \\ \hfill \left(\mathit{\alpha }{\mathit{\beta }}^{\mathrm{T}}\right)\ast \text{triu}\left(U{V}^{\mathrm{T}}\right)& =\text{triu}\left((\mathit{\alpha }\ast U){(\mathit{\beta }\ast V)}^{\mathrm{T}}\right),\hfill \end{array}$$(37)

where α² is the element-wise square of α (α² = α * α) and α * U is the element-wise product of each column of U by α

$$\begin{array}{c}\hfill {(\mathit{\alpha }\ast U)}_{k,s}={\alpha }_k{U}_{k,s}.\end{array}$$(38)

Using these relations, Eq. (35) can be rewritten as

$$\begin{array}{cc}\hfill C=& \text{diag}({\mathit{\alpha }}^2\ast A)\hfill \\ \hfill & +\text{tril}\left((\mathit{\alpha }\ast U){(\mathit{\alpha }\ast V)}^{\mathrm{T}}\right)+\text{triu}\left((\mathit{\alpha }\ast V){(\mathit{\alpha }\ast U)}^{\mathrm{T}}\right)\hfill \\ \hfill & +\text{tril}\left((\mathit{\beta }\ast U^{\prime }){(\mathit{\alpha }\ast V)}^{\mathrm{T}}\right)+\text{triu}\left((\mathit{\beta }\ast V^{\prime }){(\mathit{\alpha }\ast U)}^{\mathrm{T}}\right)\hfill \\ \hfill & +\text{tril}\left((\mathit{\alpha }\ast U){(\mathit{\beta }\ast V^{\prime })}^{\mathrm{T}}\right)+\text{triu}\left((\mathit{\alpha }\ast V){(\mathit{\beta }\ast U^{\prime })}^{\mathrm{T}}\right)\hfill \\ \hfill & +\text{diag}({\mathit{\beta }}^2\ast B)\hfill \\ \hfill & +\text{tril}\left((\mathit{\beta }\ast U^{\prime }){(\mathit{\beta }\ast V^{\prime })}^{\mathrm{T}}\right)+\text{triu}\left((\mathit{\beta }\ast V^{\prime }){(\mathit{\beta }\ast U^{\prime })}^{\mathrm{T}}\right)\hfill \\ \hfill & +\mathrm{\Sigma }.\hfill \end{array}$$(39)

Finally, the latter expression can be factorized to obtain the semiseparable representation

$$\begin{array}{c}\hfill C=\text{diag}(\mathcal{A})+\text{tril}\left(\mathcal{U}{\mathcal{V}}^{\mathrm{T}}\right)+\text{triu}\left(\mathcal{V}{\mathcal{U}}^{\mathrm{T}}\right)+\mathrm{\Sigma },\end{array}$$(40)

with

$$\begin{array}{cc}& \mathcal{A}={\mathit{\alpha }}^2\ast A+{\mathit{\beta }}^2\ast B,\hfill \\ & \mathcal{U}=\mathit{\alpha }\ast U+\mathit{\beta }\ast U^{\prime },\hfill \\ & \mathcal{V}=\mathit{\alpha }\ast V+\mathit{\beta }\ast V^{\prime }.\hfill \end{array}$$(41)

The rank of this semiseparable representation of C (number of columns in 𝒰, 𝒱) is thus the same as the rank of the underlying GP (number of columns in U, V). In the case of the more general model of Eq. (28), with several independent GP, it is straightforward to compute the semiseparable representation of the covariance matrix by vertical concatenation of the matrix 𝒰 and 𝒱 corresponding to each independent GP and the total rank is the sum of all the GP ranks.

Keeping the rank of the covariance matrix as low as possible allows us to significantly improve the performances of the method. Indeed, the cost of likelihood evaluations with a semiseparable covariance matrix of rank r scales as 𝒪(nr²). It is thus remarkable to note that introducing the derivative of the GP in the model, as well as different coefficients α and β for the different time series, does not increase the rank of the semiseparable representation of the covariance matrix. The factorization performed between Eqs. (39) and (40) is the key step that allows us to keep the same rank as the underlying GP. This factorization is achieved thanks to the specific choice of semiseparable representation introduced in Eq. (22). Indeed, instead of using (U, V, U′,V′) to represent the covariance matrix, one could use (U, V, U′,U″) or (U, V, V′,V″) (see Eqs. (17) and (20)). However, with these choices of representation, the covariance matrix would not factorize in the same way, the rank of its semiseparable representation would be twice the rank of the underlying GP, and the cost would thus quadruple.

4.2. Measurement noise and LEAF component

Using the covariance matrix representation of Eq. (40), the results described in Sect. 2 in the case of a homogeneous time series can be extended to the case of heterogeneous time series depending on several independent GP and their derivatives. The celerite algorithms can be applied in the case of purely white noise (Σ diagonal) and the S+LEAF algorithms can be applied to the more general case of close-to-diagonal correlated noise. The LEAF component of the S+LEAF model simply needs to be defined on the merged time series (t, y).

4.3. Computational cost and memory footprint

The computational cost of our model and its memory footprint are the same as the underlying celerite/S+LEAF representation. Therefore, the memory footprint scales as $\mathcal{O}((r+\overline{b})n)$ and its computational cost as $\mathcal{O}((r^2+r\overline{b}+{\overline{b}}^2)n)$, where n is the total number of measurements (including radial velocities and indicators), r is the rank of 𝒰 and 𝒱, $\overline{b}$ is the average band width of the LEAF component (considering the merged time series (t, y)). This is to be compared with a naive implementation of the same model which has a memory footprint in 𝒪(n²) and a computational cost in 𝒪(n³).

4.4. Overflows and preconditioning

As explained in Ambikasaran (2015), Foreman-Mackey et al. (2017), Delisle et al. (2020b), a naive computer implementation of the semiseparable decomposition of Eq. (9) can lead to overflows and underflows due to the exponential terms in the definition of U and V. However, Foreman-Mackey et al. (2017) proposed a simple preconditioning method to circumvent this issue in the case of the celerite model, which is also valid in the case of the S+LEAF model (Delisle et al. 2020b). Instead of using directly the matrices U and V, one might use the matrices $\stackrel{\sim }{U}$, $\stackrel{\sim }{V}$, and ϕ defined as

$$\begin{array}{cc}\hfill {\stackrel{\sim }{U}}_{i,s}& =a_{s}cos({\nu }_s{t}_i)+b_{s}sin({\nu }_s{t}_i),\hfill \\ \hfill {\stackrel{\sim }{U}}_{i,n_{\mathrm{c}}+s}& =a_{s}sin({\nu }_s{t}_i)-b_{s}cos({\nu }_s{t}_i),\hfill \\ \hfill {\stackrel{\sim }{V}}_{i,s}& =cos({\nu }_s{t}_i),\hfill \\ \hfill {\stackrel{\sim }{V}}_{i,n_{\mathrm{c}}+s}& =sin({\nu }_s{t}_i),\hfill \\ \hfill {\varphi }_{i,s}& ={\varphi }_{i,n_{\mathrm{c}}+s}={\mathit{e}}^{-{\lambda }_s(t_{i+1}-t_i)},\hfill \end{array}$$(42)

such that

$$\begin{array}{c}\hfill U_{i,s}{V}_{j,s}={\stackrel{\sim }{U}}_{i,s}{\stackrel{\sim }{V}}_{j,s}{\displaystyle \prod _{k=j}^{i-1}{\varphi }_{k,s},}\end{array}$$(43)

and all the celerite/S+LEAF algorithms can be adapted to use this representation (see Foreman-Mackey et al. 2017; Foreman-Mackey 2018; Delisle et al. 2020b). Similarly, we define the matrices $\stackrel{\sim }{U}\prime$, ${\stackrel{\sim }{V}}^{\prime }$, $\stackrel{\sim }{\mathcal{U}}$, and $\stackrel{\sim }{\mathcal{V}}$ as

$$\begin{array}{cc}\hfill {{\stackrel{\sim }{U}}^{\prime }}_{i,s}& =a_s^{\prime }cos({\nu }_s{t}_i)+b_s^{\prime }sin({\nu }_s{t}_i),\hfill \\ \hfill {{\stackrel{\sim }{U}}^{\prime }}_{i,n_{\mathrm{c}}+s}& =a_s^{\prime }sin({\nu }_s{t}_i)-b_s^{\prime }cos({\nu }_s{t}_i),\hfill \\ \hfill {{\stackrel{\sim }{V}}^{\prime }}_{i,s}& ={\lambda }_{s}cos({\nu }_s{t}_i)-{\nu }_{s}sin({\nu }_s{t}_i),\hfill \\ \hfill {{\stackrel{\sim }{V}}^{\prime }}_{i,n_{\mathrm{c}}+s}& ={\lambda }_{s}sin({\nu }_s{t}_i)+{\nu }_{s}cos({\nu }_s{t}_i),\hfill \\ \hfill \stackrel{\sim }{\mathcal{U}}& =\mathit{\alpha }\ast \stackrel{\sim }{U}+\mathit{\beta }\ast {\stackrel{\sim }{U}}^{\prime },\hfill \\ \hfill \stackrel{\sim }{\mathcal{V}}& =\mathit{\alpha }\ast \stackrel{\sim }{V}+\mathit{\beta }\ast {\stackrel{\sim }{V}}^{\prime },\hfill \end{array}$$(44)

which allow to apply the overflow-proof version of celerite/S+LEAF algorithms. The same preconditioning method can be applied in the case of the Matérn 3/2 and 5/2 kernels (see Appendix A.3).

4.5. Efficient computation of the gradient

In most applications, one needs to explore the parameter space either to maximize the likelihood using optimization algorithms or to obtain samples from the posterior distribution of parameters using Bayesian methods (MCMC, nested sampling, etc.). In both cases, many algorithms have been designed that make use of the gradient of the likelihood with respect to the parameters to improve the convergence efficiency. Following Foreman-Mackey (2018) and Delisle et al. (2020b), we deduce gradient backpropagation algorithms for all the operations used to compute the likelihood. While a detailed presentation of these backpropagation algorithms would be cumbersome, we refer the reader to Delisle et al. (2020b) Appendix B for the general idea of the method, and to the reference S+LEAF 2 implementation² for further details.

5. Application: Reanalysis of HD 13808

In this section, we apply our algorithms to reanalyze the RV time series of HD 13808. This K2V dwarf is known to harbor two planet candidates (see Mayor et al. 2011) recently published as confirmed planets by Ahrer et al. (2021). In the latter study, the authors defined several alternative models of stellar activity and performed a Bayesian model comparison letting the number of planets vary. In addition to confirming the two candidates, Ahrer et al. (2021) concluded that the best stellar activity model was a GP trained simultaneously on the RV, BIS, and $log{R}_{\text{HK}}^{\prime }$ time series, following Eq. (2), as proposed by Rajpaul et al. (2015). Modeling this GP required the authors to solve for a 738 × 738 linear system billions of times, which was very demanding in terms of computational resources. Here, we reanalyze the 246 HARPS RV measurements of HD 13808, together with the indicators time series, but modeling the GP using S+LEAF 2.

5.1. Choice of kernel

The SEP kernel (see Eq. (1)), which was used by Ahrer et al. (2021) for their analysis of HD 13808, is not semiseparable and thus cannot be modeled with S+LEAF 2. However, other quasiperiodic kernels, such as the SHO kernel proposed by Foreman-Mackey et al. (2017), admit a semiseparable representation. Here, we aim to reproduce the main characteristics of the SEP kernel but using a semiseparable kernel. The SEP kernel is the product of a squared-exponential and the exponential of a sinusoidal (see Eq. (1)):

$$\begin{array}{c}\hfill k(\mathrm{\Delta }t)={\sigma }^{2}exp(-\frac{\mathrm{\Delta }{t}^2}{2{\rho }^2})exp(-\frac{{sin}^2\left(\frac{\pi \mathrm{\Delta }t}{P}\right)}{2{\eta }^2}).\end{array}$$(45)

The second part can be expanded as a power series, assuming 2η ≳ 1

$$\begin{array}{cc}\hfill k(\mathrm{\Delta }t)& ={\sigma }^{2}exp(-\frac{\mathrm{\Delta }{t}^2}{2{\rho }^2})(1-\frac{{sin}^2\left(\frac{\pi \mathrm{\Delta }t}{P}\right)}{2{\eta }^2}+\frac{{sin}^4\left(\frac{\pi \mathrm{\Delta }t}{P}\right)}{8{\eta }^4}+\mathcal{O}\left({\eta }^{-6}\right))\hfill \\ \hfill & ={\sigma }^{2}exp(-\frac{\mathrm{\Delta }{t}^2}{2{\rho }^2})\frac{1+fcos\left(\nu \mathrm{\Delta }t\right)+\frac{f^2}{4}cos\left(2\nu \mathrm{\Delta }t\right)}{1+f+\frac{f^2}{4}}+\mathcal{O}\left(f^3\right),\hfill \end{array}$$(46)

with f = (2η)⁻² and ν = 2π/P. This kernel thus introduces some correlation at the rotation period P, but also at the harmonics (P/2, P/3, etc.), with amplitudes decaying rapidly (scaling as η⁻²ⁿ for the harmonics P/n). The squared-exponential part implies that the correlations vanish over long timescales (Δt ≫ ρ). The squared-exponential kernel is not semiseparable, but the Matérn 1/2 (simple exponential decay), 3/2, and 5/2 kernels admit a semiseparable decomposition (see Appendix A), and offer a similar decay of the correlation over long timescales. The SEP kernel could thus be roughly approximated by:

$$\begin{array}{c}\hfill k(\mathrm{\Delta }t)={\sigma }^{2}exp(-\frac{\mathrm{\Delta }t}{\rho })\frac{1+fcos\left(\nu \mathrm{\Delta }t\right)+\frac{f^2}{4}cos\left(2\nu \mathrm{\Delta }t\right)}{1+f+\frac{f^2}{4}}·\end{array}$$(47)

However, a GP following this kernel would not be differentiable (k′(0) = − σ²/ρ ≠ 0). In order to ensure differentiability, we introduce a modified kernel, which is the combination of a Matérn 3/2 kernel and two underdamped SHO terms

$$\begin{array}{c}\hfill k(\mathrm{\Delta }t)={\sigma }^2\frac{k_{3/2}(\mathrm{\Delta }t)+f{k}_{\mathrm{SHO},\mathrm{fund}.}(\mathrm{\Delta }t)+\frac{f^2}{4}{k}_{\mathrm{SHO},\mathrm{harm}.}(\mathrm{\Delta }t)}{1+f+\frac{f^2}{4}},\end{array}$$(48)

where

$$\begin{array}{cc}\hfill k_{3/2}(\mathrm{\Delta }t)& =exp(-\frac{\sqrt{3}\mathrm{\Delta }t}{\rho })(1+\frac{\sqrt{3}\mathrm{\Delta }t}{\rho }),\hfill \\ \hfill k_{\mathrm{SHO},\mathrm{fund}.}(\mathrm{\Delta }t)& =exp(-\frac{\mathrm{\Delta }t}{\rho })(cos\left(\nu \mathrm{\Delta }t\right)+\frac{1}{\nu \rho }sin\left(\nu \mathrm{\Delta }t\right)),\hfill \\ \hfill k_{\mathrm{SHO},\mathrm{harm}.}(\mathrm{\Delta }t)& =exp(-\frac{\mathrm{\Delta }t}{\rho })(cos\left(2\nu \mathrm{\Delta }t\right)+\frac{1}{2\nu \rho }sin\left(2\nu \mathrm{\Delta }t\right)).\hfill \end{array}$$(49)

The kernel of Eq. (48), which we refer to as the Matérn 3/2 exponential periodic (MEP) kernel in the following, is differentiable and presents the main characteristics of the SEP kernel, while being semiseparable. The semiseparable representation of the MEP kernel is of rank r = 6. We note that similar kernels have already been used in the literature to model stellar activity in photometric time series (e.g., David et al. 2019; Gillen et al. 2020).

It should be noted that the MEP kernel is once mean square differentiable but not twice, which means that G′ is well defined but it is not itself differentiable. Since the time series (RV and activity indicators) are modeled as combinations of G and G′, we could require G′ to be differentiable to obtain a smoother model. Such a twice mean square differentiable kernel should satisfy k⁽³⁾(0) = 0, in addition to the mandatory differentiability condition (k′(0) = 0). For the MEP kernel, $k_{\mathrm{MEP}}^{(3)}(0)$ is non-zero, however, in practice this kernel seems to produce smooth time series (see Sect. 5.3). Nevertheless, it is possible to design semiseparable kernels that are rigorously twice differentiable. For instance, the Matérn 5/2 kernel is twice differentiable and semiseparable with a rank of 3 (see Appendix A). In Appendix B we present the exponential-sine (ES) and the exponential-sine periodic (ESP) kernels. Both kernels are twice differentiable and semiseparable. The ES kernel is of rank 3 and closely resemble the SE kernel, while the ESP kernel is of rank 15 and approximates the SEP kernel very well. We find very similar results when using the ESP kernel instead of the MEP kernel, while the cost of likelihood evaluations is roughly doubled because of the higher rank of the ESP kernel (see Appendix B). In the following, we thus adopt the MEP kernel as a replacement for the SEP kernel, to reduce the computational cost and since it produces smooth time series, at least in our case. In the general case, the ESP kernel would typically generate smoother times series than the MEP, but with a doubled cost.

Following Rajpaul et al. (2015) and Ahrer et al. (2021) we use the GP to model simultaneously the RV, BIS, and $log{R}_{\text{HK}}^{\prime }$ time series of HD 13808 according to (see Eq. (2)):

$$\begin{array}{cc}& \mathrm{\Delta }\mathrm{RV}={\alpha }_{\mathrm{RV}}G(t)+{\beta }_{\mathrm{RV}}{G}^{\prime }(t),\hfill \\ \hfill & \mathrm{\Delta }\mathrm{BIS}={\alpha }_{\mathrm{BIS}}G(t)+{\beta }_{\mathrm{BIS}}{G}^{\prime }(t),\hfill \\ \hfill & \mathrm{\Delta }log{R}_{\mathrm{HK}}^{\prime }={\alpha }_{log{R}_{\mathrm{HK}}^{\prime }}G(t),\hfill \end{array}$$(50)

where the coefficients α, β and the kernel’s hyperparameters (P, ρ, η) need to be determined.

5.2. Performances

To evaluate the performances of S+LEAF 2 for the model of Eqs. (48) and (50), we generated random RV, BIS, and $log{R}_{\text{HK}}^{\prime }$ time series and record the cost of the likelihood evaluation as a function of the number of generated measurements. We compared S+LEAF 2 performances using the MEP kernel with the computational cost of evaluating the likelihood using the full covariance matrix with the SEP kernel. The two implementations were run on the same computer, using a single core. The results are shown in Fig. 1 and confirm the 𝒪(n) scaling of S+LEAF 2, while the implementation of the full covariance matrix indeed scales as 𝒪(n³) for a large value of n. In the case of HD 13808, the total number of measurements is 738 (3 × 246), and based on Fig. 1, we obtain a gain of a factor ∼130 in computing time by using S+LEAF 2 instead of the naive implementation. For larger data sets, the gain would be even greater.

Fig. 1. Cost of a likelihood evaluation as a function of the total number of measurements using S+LEAF 2 or the full covariance matrix (see Sect. 5.2).

In addition to these performances tests, we also ran numerical precision tests. We used the same GP model as for the performances tests and assessed the stability of S+LEAF 2 by computing CC⁻¹x, that is, applying the solving and dot product algorithms on a random merged time series x. In theory, we should find CC⁻¹x = x, however, due to the limited machine precision, numerical errors accumulate and the results stray slightly from x. We computed the root mean square (rms) of these errors, for the S+LEAF 2 methods and for the full covariance matrix (see Fig. 2). In both cases, numerical errors grow with the total number of measurements. However, the level and growth rate are lower when using S+LEAF 2. For a given number of measurements, the number of arithmetic operations is lower when using S+LEAF 2, which could explain the improvement in precision.

Fig. 2. Precision of the linear solving operation as a function of the total number of measurements using S+LEAF 2 or the full covariance matrix (see Sect. 5.2).

5.3. Periodogram and false alarm probability

We analyzed HD 13808 data using the periodogram and false alarm probability (FAP) approach presented in Delisle et al. (2020a). For each frequency, all the linear parameters (here the RV and indicators offsets) are re-adjusted together with the amplitudes of the cosine and sine at the considered frequency (only applied to the RV time series). The framework of Delisle et al. (2020a) only requires slight adaptations to account for the joint fit of RV and activity indicators, which are detailed in Appendix C.

We first performed a fit of a base model including the offsets γ_RV, γ_BIS, ${\gamma }_{log{R}_{\text{HK}}^{\prime }}$, jitter terms added in quadrature to the measurements errorbars, σ_RV, σ_BIS, ${\sigma }_{log{R}_{\text{HK}}^{\prime }}$, the kernels hyperparameters, P_GP, ρ_GP, η_GP, and the amplitudes, α, β, by maximizing the likelihood. Then using the fitted noise parameters, we computed a first periodogram, reajusting the offsets for each considered frequency. The resulting periodogram is plotted in Fig. 3 (top). We observe a very significant peak at 14.19 d (FAP = 6.7 × 10⁻²²). We then fitted a Keplerian orbit to this signal, and reajusted all parameters. A second periodogram (Fig. 3, middle) was then computed on the residuals, still reajusting the offsets, but keeping the planetary and noise parameters fixed. A second significant peak is visible at 53.7 d (FAP = 7.7 × 10⁻¹²). As for the first planet, we performed a global fit of all the parameters after including this planet in the model. These fitted parameters are presented in Table 1, the kernel function corresponding to the fitted hyperparameters is illustrated in Fig. 5, and the residuals of each time series superimposed with the GP prediction are shown in Figs. 6 and 7. Finally, the periodogram of the residuals after fitting both planets (Fig. 3, bottom) does not show any significant peak (FAP above 1%). Our conclusions, based on this periodogram and FAP approach, thus agree with the findings of Ahrer et al. (2021), with two confirmed planets at 14.19 d and 53.7 d.

Fig. 3. Periodograms of the raw RV time series of HD 13808 (top) as well as of the residuals after subtracting the 14.19 d (center) and the 53.7 d (bottom) planets.

Fig. 4. Same as Fig. 3 but neglecting the harmonics component of the MEP kernel (see Eq. (48)).

Fig. 5. Kernel function used to model HD 13808’s activity (MEP, see Eq. (48)). The GP hyperparameters are taken from the best fit of the two planets solution (Table 1). For comparison, the SEP kernel, which the MEP kernel is design to roughly mimic, is also plotted using the same set of hyperparameters.

Fig. 6. GP prediction (conditional distribution) from the best fit of the two planets solution (Table 1). The prediction is plotted for the GP and its derivative (top) and the full GP prediction for the RV, BIS, and $log{R}_{\text{HK}}^{\prime }$ time series superimposed with the corresponding residuals (bottom three plots).

Fig. 7. Zoom of Fig. 6 around epoch 2 453 720 BJD.

We note that the first harmonics part (k_{SHO,  harm.}) in the MEP kernel (Eq. (48)) plays a key role in the GP modeling, even if its amplitude is significantly smaller than the fundamental term (k_{SHO,  fund.}). Indeed, when neglecting the harmonics term (see Fig. 4), we find a third highly significant signal around 19 d (FAP = 4.3 × 10⁻⁵), which corresponds to half the period of the GP (P_GP ≈ 38 d, see Table 1). On the contrary, when including the harmonics part, the peak around 19 d is no more significant (see Fig. 3). This is not surprising since stellar activity is expected to introduce signals in the RV and indicators time series at the rotation period, but also at its harmonics, and, in particular, the first harmonics.

By design, a kernel which has power at the rotation harmonics will have a tendency to absorb signals at this period (here 19 d), even if those are not due to stellar activity. To further test whether the 19 d signal could be due to a planet, Hara et al. (2022a) checked the consistency of this signal over the timespan of the data set. The authors showed that the 19 d signal, while not statistically significant, exhibits a stable presence accross the data set, which might motivate further observations.

5.4. Bayesian framework and false inclusion probability

In order to more directly compare our results with the study of Ahrer et al. (2021), we performed a full Bayesian evidence analysis using the nested sampling algorithm POLYCHORD (Handley et al. 2015). Recently, POLYCHORD was used for radial velocity exoplanet detection by Ahrer et al. (2021), Rajpaul et al. (2021), and Unger et al. (2021). We aimed at reproducing a similar analysis as Ahrer et al. (2021), but using S+LEAF 2 and the GP model detailed in Sect. 5.1.

The prior distributions we used for all parameters are detailed in Table 2. We kept mostly the same priors used by Ahrer et al. (2021), with only a few changes. We set the amplitude for the RV term of the GP to be strictly positive to avoid degeneracies in the sign of the amplitudes. We also shifted the prior for the mean longitudes from [0, 2π] to [ − π, π]. Indeed, the mean longitude of planet c is very close to 0, and this shift avoids splitting the peak in the posterior, thus improving the convergence efficiency. For all runs we used a number of live points equal to 50 n_dim, that is, 50 times the number of free parameters. For the precision criterion (stopping criterion), we used 10⁻⁹.

We ran POLYCHORD for models with zero up to three planets and each model was run five times to obtain estimates of the value and uncertainty of the evidence. The evidence for each model is presented in Table 3. We see a steady and significant increase in evidence up to the model with two planets. The two and three planet models present similar evidence (Δln Z = 3.5 ± 2.5, compatible at 1.4σ). Moreover, even if the three planet model is marginally favored, no clear period for a third planet emerges from the posteriors. To illustrate this, we computed the false inclusion probability (FIP) periodogram (Hara et al. 2022b) from the posteriors of all our runs with zero to three planets (see Fig. 8). The true inclusion probability (TIP) provides the probability for the system to host at least one planet in a given (small) range of period. On the contrary, the FIP (=1–TIP) is the probability that the system does not host any planet in the given range of period. The FIP periodogram of Fig. 8 is computed with a window size of 1/(t_max − t_min) in frequency. We observe in Fig. 8 two very significant peaks (FIP less than 10⁻⁸) around 14.1 d and 53.7 d, which confirms that these planets should be included in the model. Then, two smaller peaks around 12 d and 19 d are also visible but none of them is significant (FIP higher than 70%). The posterior distribution of parameters of the two-planet runs is given in Table 1. Our results are in agreement with the periodogram and FAP approach detailed in Sect. 5.3, and with Ahrer et al. (2021) findings.

Fig. 8. FIP periodogram for HD 13808. In blue we represent the FIP (false inclusion probability) and in yellow the TIP (true inclusion probability).

6. Conclusion

In this article, we have presented S+LEAF 2, a GP framework that is able to efficiently model multiple time series simultaneously. Classical GP models have a computational cost that scales as the cube of the number of measurements, which makes them prohibitive in terms of computational resources for large data sets. The computational cost of our GP framework scales linearly with the data set size, which allows for tractable computations even for large data sets.

This work builds on previous studies that provided efficient GP models for single time series (in particular the celerite and S+LEAF models, see Rybicki & Press 1995; Ambikasaran 2015; Foreman-Mackey et al. 2017; Delisle et al. 2020b). It is also inspired by a recent generalization of celerite to the case of two-dimensional data sets (Gordon et al. 2020), but extends it by accounting for the GP derivatives. These derivatives are especially important when modeling the effect of stellar activity on RV time series (see Aigrain et al. 2012; Rajpaul et al. 2015). Our framework additionally accounts for time series that do not share the same calendar, which is useful to train a GP simultaneously on RV and photometric measurements taken with two different instruments (e.g., Haywood et al. 2014).

We applied our methods to reanalyze the RV time series of the nearby K2 dwarf HD 13808. Our results are very similar to a recent state-of-the-art study of the same system (Ahrer et al. 2021) and we confirm the two planets announced in this article. However, we have shown that using our framework allowed us to dramatically decrease (by more than two orders of magnitude) the computational cost of the GP modeling. While the data set analyzed here consists of 738 measurements (RV, BIS, and $log{R}_{\text{HK}}^{\prime }$ at 246 epochs), the gain of using S+LEAF 2 would be even greater for larger data sets. Some data sets (e.g., HARPS or HARPS-N Sun-as-a-star RV time series) that could not have been analyzed with such a GP modeling are now achievable with our GP framework.

Finally, it is worth noting that the results from the periodogram and FAP approach are consistent with the much more computer intensive Bayesian evidence calculations using nested sampling. This illustrates the power of the periodogram and FAP computation including a correlated noise model, as proposed by Delisle et al. (2020a).

Acknowledgments

We thank the referee, D. Foreman-Mackey, for his very constructive feedback that helped to improve this manuscript. We acknowledge financial support from the Swiss National Science Foundation (SNSF). This work has, in part, been carried out within the framework of the National Centre for Competence in Research PlanetS supported by SNSF.

References

Appendix A: Semiseparable representation of Matérn 3/2 and Matérn 5/2 kernels

The Matérn 3/2 and 5/2 covariances are widely used in various fields of statistics. Their kernel functions are written as

$$\begin{array}{cc}& k_{3/2}(\mathrm{\Delta }t)={\sigma }^2(1+\mathrm{\Delta }x){\mathbf{e}}^{-\mathrm{\Delta }x},\hfill \\ \hfill & k_{5/2}(\mathrm{\Delta }t)={\sigma }^2(1+\mathrm{\Delta }x+\frac{1}{3}\mathrm{\Delta }{x}^2){\mathbf{e}}^{-\mathrm{\Delta }x},\hfill \end{array}$$(A.1)

where x is the rescaled time:

$$\begin{array}{cc}\hfill x& =\sqrt{3}\frac{t}{\rho }\mathrm{for}\mathrm{the}\mathrm{Mat}\acute{\mathrm{e}}\mathrm{rn}3/2\mathrm{kernel},\hfill \\ \hfill x& =\sqrt{5}\frac{t}{\rho }\mathrm{for}\mathrm{the}\mathrm{Mat}\acute{\mathrm{e}}\mathrm{rn}5/2\mathrm{kernel},\hfill \end{array}$$(A.2)

and σ and ρ are the GP hyperparameters.

A.1. Semiseparable representation

For two times t_i > t_j, we have

$$\begin{array}{c}\hfill k_{3/2}(t_i-t_j)={\sigma }^2(\left(x_i{\mathbf{e}}^{-x_i}\right){\mathbf{e}}^{x_j}+{\mathbf{e}}^{-x_i}\left((1-x_j){\mathbf{e}}^{x_j}\right)),\end{array}$$(A.3)

which provides the semiseparable representation of rank 2

$$\begin{array}{cc}& A_i={\sigma }^2,\hfill \\ \hfill & U_{i,1}={\sigma }^2{x}_i{\mathbf{e}}^{-x_i},V_{i,1}={\mathbf{e}}^{x_i},\hfill \\ \hfill & U_{i,2}={\sigma }^2{\mathbf{e}}^{-x_i},V_{i,2}=(1-x_i){\mathbf{e}}^{x_i}\hfill \end{array}$$(A.4)

for the Matérn 3/2 kernel. Similarly, the Matérn 5/2 kernel admits the semiseparable representation of rank 3

$$\begin{array}{cc}& A_i={\sigma }^2,\hfill \\ \hfill & U_{i,1}={\sigma }^2(x_i+\frac{x_i^2}{3}){\mathbf{e}}^{-x_i},V_{i,1}={\mathbf{e}}^{x_i},\hfill \\ \hfill & U_{i,2}={\sigma }^2{\mathbf{e}}^{-x_i},V_{i,2}=(1-x_i+\frac{x_i^2}{3}){\mathbf{e}}^{x_i},\hfill \\ \hfill & U_{i,3}={\sigma }^2{x}_i{\mathbf{e}}^{-x_i},V_{i,3}=-\frac{2}{3}{x}_i{\mathbf{e}}^{x_i}.\hfill \end{array}$$(A.5)

These representations are not unique and the choice of splitting into the 2 (Matérn 3/2) or 3 (Matérn 5/2) semiseparable terms is arbitrary.

A.2. Derivative of a Matérn Gaussian process

A GP following the Matérn 3/2 or 5/2 kernel is always differentiable, independently of the hyperparameters. Following the same reasoning as in Sect. 3, we compute the derivatives of U and V, as well as the matrix B, which appear in the covariance matrix between the GP G(t) and its derivative G′(t) and in the covariance matrix of G′(t) itself (see Eq. (22)). We find

$$\begin{array}{cc}\hfill U_{i,1}^{\prime }& =\frac{\sqrt{3}{\sigma }^2}{\rho }(1-x_i){\mathbf{e}}^{-x_i},\hfill \\ \hfill U_{i,2}^{\prime }& =-\frac{\sqrt{3}{\sigma }^2}{\rho }{\mathbf{e}}^{-x_i},\hfill \\ \hfill V_{i,1}^{\prime }& =\frac{\sqrt{3}}{\rho }{\mathbf{e}}^{x_i},\hfill \\ \hfill V_{i,2}^{\prime }& =-\frac{\sqrt{3}}{\rho }{x}_i{\mathbf{e}}^{x_i},\hfill \\ \hfill B_i& =\frac{3{\sigma }^2}{{\rho }^2},\hfill \end{array}$$(A.6)

for the Matérn 3/2 kernel, and

$$\begin{array}{cc}\hfill U_{i,1}^{\prime }& =\frac{\sqrt{5}{\sigma }^2}{3\rho }(3-x_i-x_i^2){\mathbf{e}}^{-x_i},\hfill \\ \hfill U_{i,2}^{\prime }& =-\frac{\sqrt{5}{\sigma }^2}{\rho }{\mathbf{e}}^{-x_i},\hfill \\ \hfill U_{i,3}^{\prime }& =\frac{\sqrt{5}{\sigma }^2}{\rho }(1-x_i){\mathbf{e}}^{-x_i},\hfill \\ \hfill V_{i,1}^{\prime }& =\frac{\sqrt{5}}{\rho }{\mathbf{e}}^{x_i},\hfill \\ \hfill V_{i,2}^{\prime }& =\frac{\sqrt{5}}{3\rho }(x_i^2-x_i){\mathbf{e}}^{x_i},\hfill \\ \hfill V_{i,3}^{\prime }& =-\frac{2\sqrt{5}}{3\rho }(1+x_i){\mathbf{e}}^{x_i},\hfill \\ \hfill B_i& =\frac{5{\sigma }^2}{3{\rho }^2},\hfill \end{array}$$(A.7)

for the Matérn 5/2 kernel.

A.3. Overflows and preconditioning

To avoid overflows (see Sect. 4.4), these representations can be adapted to use the same preconditioning method as in the case of the classical celerite quasiperiodic terms. The preconditioning matrix ϕ is defined as

$$\begin{array}{c}\hfill {\varphi }_{i,s}={\mathbf{e}}^{-(x_{i+1}-x_i)},\end{array}$$(A.8)

and the preconditioned matrices $\stackrel{\sim }{U}$, $\stackrel{\sim }{V}$, $\stackrel{\sim }{U}\prime$, and ${\stackrel{\sim }{V}}^{\prime }$ are obtained by dropping the exponential terms from the definitions of U, V, U′, and V′. For instance, we have

$$\begin{array}{cc}& {\stackrel{\sim }{U}}_{i,1}={\sigma }^2{x}_i,{\stackrel{\sim }{V}}_{i,1}=1,\hfill \\ \hfill & {\stackrel{\sim }{U}}_{i,2}={\sigma }^2,{\stackrel{\sim }{V}}_{i,2}=1-x_i,\hfill \\ \hfill & {{\stackrel{\sim }{U}}^{\prime }}_{i,1}=\frac{\sqrt{3}{\sigma }^2}{\rho }(1-x_i),{{\stackrel{\sim }{V}}^{\prime }}_{i,1}=\frac{\sqrt{3}}{\rho },\hfill \\ \hfill & {{\stackrel{\sim }{U}}^{\prime }}_{i,2}=-\frac{\sqrt{3}{\sigma }^2}{\rho },{{\stackrel{\sim }{V}}^{\prime }}_{i,2}=-\frac{\sqrt{3}}{\rho }{x}_i,\hfill \\ \hfill & {\varphi }_{i,1}={\varphi }_{i,2}={\mathbf{e}}^{-(x_{i+1}-x_i)},\hfill \end{array}$$(A.9)