  The radial velocity method is a very productive technique used to detect and
confirm extrasolar planets. The most recent spectrographs, such as ESPRESSO or
EXPRES, have the potential to detect Earth-like planets around Sun-like stars.
However, stellar activity can induce radial velocity variations that dilute or
even mimic the signature of a planet. A widely recognized method for
disentangling these signals is to model the radial velocity time series,
jointly with stellar activity indicators, using Gaussian processes and their
derivatives. However, such modeling is prohibitive in terms of computational
resources for large data sets, as the cost typically scales as the total number
of measurements cubed.
  Here, we present S+LEAF 2, a Gaussian process framework that can be used to
jointly model several time series, with a computational cost that scales
linearly with the data set size. This framework thus provides a
state-of-the-art Gaussian process model, with tractable computations even for
large data sets.
  We illustrate the power of this framework by reanalyzing the 246 HARPS radial
velocity measurements of the nearby K2 dwarf HD 138038, together with two
activity indicators. We reproduce the results of a previous analysis of these
data, but with a strongly decreased computational cost (more than two order of
magnitude). The gain would be even greater for larger data sets.
