Modelling stochastic and quasi-periodic behaviour in stellar time-series: Gaussian process regression versus power-spectrum fitting

As the hunt for an Earth-like exoplanets has intensified in recent years, so has the effort to characterise and model the stellar signals that can hide or mimic small planetary signals. Stellar variability arises from a number of sources, including granulation, supergranulation, oscillations and activity, all of which result in quasi-periodic or stochastic behaviour in photometric and/or radial velocity observations. Traditionally, the characterisation of these signals has mostly been done in the frequency domain. However, the recent development of scalable Gaussian process regression methods makes direct time-domain modelling of stochastic processes a feasible and arguably preferable alternative, obviating the need to estimate the power spectral density of the data before modelling it. In this paper, we compare the two approaches using a series of experiments on simulated data. We show that frequency domain modelling can lead to inaccurate results, especially when the time sampling is irregular. By contrast, Gaussian process regression results are often more precise, and systematically more accurate, in both the regular and irregular time sampling regimes. While this work was motivated by the analysis of radial velocity and photometry observations of main sequence stars in the context of planet searches, we note that our results may also have applications for the study of other types of astrophysical variability such as quasi-periodic oscillations in X-ray binaries and active galactic nuclei variability.


INTRODUCTION
Stochastic processes are ubiquitous in astrophysical time-series observations.Quasi-Periodic Oscillations (QPOs) arise in accreting systems, from Active Galactic Nuclei (AGN) to compact binaries, while magnetic activity, convection and oscillations in stellar interiors all give rise to stochastic variability ranging from strongly periodic to entirely aperiodic.The resulting time-series are in general too complex to describe using parametric models, but their Power Spectral Density (PSD) can often be modelled using simple analytic functions or sums thereof.It has thus become common practice, in many astrophysical sub-fields, to model time-series observations displaying stochastic behaviour by estimating their PSDs and fitting the latter in the frequency domain.
The functions used to model the PSD are often generalized forms of the Lorentzian or Cauchy distribution: where  0 is the central frequency,  is the amplitude,  is the slope and  is the width.(For a standard Cauchy or Lorentzian distribution normalised to unity,  = 2 and  =  −1 .)More complex models can be constructed by co-adding multiple components.The inferred parameters of the individual components are frequently given a physical interpretation, so it is important to understand how reliable ★ E-mail: niamh.osullivan@physics.ox.ac.uk these inferences are.Are there degeneracies between the parameters, particularly when multiple components are included, and when the number of components is not known a priori?How robust are the parameter uncertainties derived from the fits?
While it is a very powerful and widespread approach, there are a number of drawbacks to modelling stochastic processes in the Fourier domain.The dynamic range of PSDs is generally very large, spanning many orders of magnitude in both axes, and the number of independent estimates of the power at different frequencies is a strong function of the frequency itself, which makes fitting analytic functions to the PSD notoriously challenging.Furthermore, many real-world astrophysical time-series have irregular time-sampling, which precludes using the Fast Fourier Transform (FFT) to estimate the PSD.The Generalised Lomb-Scargle (GLS) periodogram is a popular method to approximate the PSD of irregularly sampled time-series data, but the results are highly sensitive to the window function of the observations Zechmeister & Kürster (2009).When the sampling is almost regular (for example for regular sampling with gaps), it is common practice to interpolate over the missing data, but that carries its own problems.
An alternative approach, which circumvents some of these issues, is to model the data directly in the time-domain, using a stochastic process model whose PSD can be written down analytically in terms of the model's parameters.Widely-used examples include (continuous or discrete) Auto-Regressive Moving Average (ARMA) models as well as certain types of Gaussian Process (GP) models (see Aigrain & Foreman-Mackey (2023) for a recent review of GP regression for  2023).The full model (solid purple line) is the sum of quasiperiodic components corresponding to the envelope of the solar oscillations (green dashed line), and to the solar rotation period and its first two harmonics (red dashed line), aliases thereof (black dashed line), plus aperiodic components corresponding to granulation (blue dashed line) and super-granulation (orange dashed line), plus white (photon) noise (cyan dashed line).astronomical time series).For a long time, the computational cost of this approach was prohibitive for time series consisting of ≳ 10 3 observations, but recent developments have overcome this limitation (Foreman-Mackey et al. 2017;Foreman-Mackey 2018).
In this study, we use simulations to quantify the limitations of both approaches and to ascertain when one might be preferable over the other.This work was initially motivated by recent attempts to quantify the amplitudes and timescales of different contributions to stellar 'noise' in Radial Velocity (RV) planet searches, particularly the study of Al Moulla et al. (2023, hereafter A23), who used data from the HARPS and HARPS-North solar telescopes to reconstruct and model the PSD of solar disk-integrated Radial Velocity (RV) variability on timescales ranging from minutes to months.Their final fit to the solar PSD is reproduced in Figure 1, illustrating the complexity of the overall model, and the significant overlap between the timescales of the different components.These results are likely to have a significant impact on both the theoretical modelling of stellar (super-)granulation and the observing strategy of future Extreme Precision RV (EPRV) surveys.It is therefore important to understand how robust the fitted parameters are to the PSD-fitting methodology and to see if GP regression outperforms the 'standard' approach.We also note that similar types of PSD fitting are commonplace in other fields, including asteroseismology, QPOs and AGN variability, where GP modelling is also starting to emerge as a viable alternative.
The remainder of this paper is structured as follows.We start by reviewing the common practice of fitting generalized Lorentzians to the PSD of solar variability -the so called 'solar background' in Section 2. In Section 3, we discuss GP kernels that can be used to fit the solar background in the time domain.In Section 4, we describe how we simulate time-series with known PSDs, and fit them them in the frequency-and in the time-domain.We then evaluate and compare the precision and accuracy of the two approaches, for different intrinsic PSDs and time-sampling properties, in Section 5. We discuss the implications of our findings and outline our future work in Section 6.

MODELLING THE SOLAR BACKGROUND IN THE FREQUENCY DOMAIN
In this section, we review the practice of characterising stellar variability on timescales ranging from hours to months by fitting a sum of modified Lorentzian functions to the PSD of the observations.This approach was first developed by Harvey (1985), in the context of helioseismic RV observations, and has since become widespread, being routinely applied to both solar and stellar RV and photometric observations (see e.g.Goldreich et al. 1994;Aigrain et al. 2004;Al Moulla et al. 2023) 2.1 The 'solar background' The term 'solar background' refers to the broad spectrum of solar variability in RV or photometry, which acts as a nuisance signal in the detection and study of -mode oscillations (see e.g.Harvey 1985).It is modelled as the sum of up to five components, in order of decreasing frequency: the broad envelope of -mode oscillations, caused by trapped sound waves in the stellar atmosphere (Kjeldsen & Bedding 1995;Goldreich et al. 1994), granulation, caused by convective upflows (Rimmele et al. 1995;Lefebvre et al. 2008;Roudier et al. 1991), super-granulation, the origin of which is still under debate (see Rieutord & Rincon 2010), the rotational modulation of active regions, caused by stellar spots rotating in and out of our field of view (Saar & Donahue 1997;Meunier et al. 2010;Haywood et al. 2016), and long term trends caused by the star's activity cycle (Santos et al. 2010).Beyond helio-(and astero-)seismology, this background is also an important nuisance signal that needs to be understood and mitigated for exoplanet searches, as it can hide and mimic planet signals in both Doppler and photometric stellar time-series.As illustrated in Figure 1, the components of the solar background have timescales ranging from minutes (in the case of -modes), to hours (granulation), to days (super-granulation), to weeks or months (rotation) and in some cases decades (magnetic activity cycles).Carefully designed observing strategies can be used to minimize the impact of the higher-frequency components on RV planet searches, for example by selecting exposure times that average out the -modes and spacing out consecutive exposures by more than the characteristic timescale of granulation (Dumusque et al. 2011).However, the super-granulation component is particularly problematic, as it has characteristic timescales of 1-2 days (Meunier et al. 2015), which limits the frequency with which a given target can be monitored without being affected by correlated noise.Meunier (2021) provides a more detailed review of the effect of solar and stellar variability on RV planet searches.
The characterisation of the components of the solar background has been traditionally done by fitting analytic functions to estimates of the PSDs of the observations.In the rest of this section, we give a brief overview of the functional forms used in the literature for modelling the different components and discuss some of the subtleties of the fitting process.Harvey (1985) introduced the notion of modelling granulation and super-granulation signals in solar RVs as a random process with exponentially decaying auto-covariance, with variance  2 and characteristic decay timescale :

Aperiodic Components
(2) The PSD of this process is (3) This is often known as a Harvey function (hence the subscript "H").
There is a difference of a factor 2 between the expression given by Harvey (1985) for the PSD of such a process and Equation (3); we can check that the latter is correct by integrating it over the full frequency domain: satisfying Parseval's theorem.This discrepancy means that some of the literature estimates of  derived from PSD fits of this type may need to be adjusted by a factor √ 2. In any case, later studies (e.g., Harvey et al. 1993;Lefebvre et al. 2008;Michel et al. 2009;Karoff et al. 2013;Kallinger et al. 2014) have shown that the power spectrum of granulation in particular is significantly steeper, and better approximated by a PSD of the form (5) This is the expression used in more recent work, including A23.The difference between Equations (3) and ( 5) is illustrated on Figure 2. The variance of this process is and thus depends on  as well as  H . Nordlund et al. (1997) provide a physical explanation for the apparently steeper power spectrum of granulation.They argue that the Harvey (1985) model corresponds to a turbulent cascade, which has power on all timescales, whereas granulation is a convection phenomenon, for which there is a minimum characteristic spatial (and thus temporal) scale, and thus a steep decay in the power spectrum.
Note that some studies, for example Harvey et al. (1993), let the index of the power spectrum be a free parameter, usually denoted as : but this tends to result in fitted values in the range 3.5 ≤  ≤ 5.5, consistent with  = 4 within the uncertainties.For lower-frequency components such as super-granulation (in RV) or faculae (in photometry, Karoff et al. 2013), even steeper PSD indices ( ∼ 6) are sometimes reported, but it is not clear how robust those estimates are, as the high-frequency tails of these components overlap strongly with the granulation component.

Periodic Components
To model (quasi-)periodic components, for example the envelope of the -modes, Harvey et al. (1993) and later Lefebvre et al. (2008) used a modified Lorentzian function of the form: where  0 is the central frequency and  L is the power at  =  0 .To reduce the number of free parameters, many recent studies that use this type of periodic model, including A23, fix  = 0 and  = 1, so that in which case Γ is the half-width at half maximum (HWHM).The variance of this process is: Once again, it depends not only on the amplitude but also on the width (or damping time-scale).Some other recent studies, for example Kallinger et al. (2014), use a Gaussian instead of a Lorentzian function to model the envelope of the -modes: where in this case  is the standard deviation of the Gaussian (not to be confused with the standard deviation of the process).This is more in line with the standard asteroseismology practice for measuring  max , the peak frequency of the stellar -modes, which is then used in scaling relations to determine the stellar mean density (see e.g.Chaplin & Miglio 2013).Neither of the above formulations have a particularly clear physical motivation, however; their use is justified empirically rather than theoretically.

PSD estimation and fitting
For regularly sampled observations, PSD estimation is straightforward using the Fast Fourier Transform (FFT).For irregularly sampled time-series, such as ground-based RV observations, the Lomb-Scargle (LS) periodogram (Lomb 1976;Scargle 1982) and later derivatives thereof are the PSD estimation method of choice in the astronomical community.It is worth noting that care is needed when selecting the normalisation of the LS periodogram so that it gives results that are directly comparable to the FFT-estimated PSDs.See VanderPlas (2018) for a detailed discussion of the LS periodogram and its interpretation.
The fitting process itself is often challenging, owing to the wide dynamic range covered by the PSD in both frequency and power density, the difficulty in assigning relative weights to the individual PSD samples, and the overlapping nature of the background components, which leads to significant degeneracies between their parameters.Here we describe the procedure used by A23, which is fairly representative of previous studies.They estimate the PSD using the Generalised Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009), which is then binned to a regular grid in log frequency space, to give the low-frequency components adequate weight in the fit.As previously mentioned, the slope of the individual background components is fixed to improve convergence and reduce degeneracies between the components.Any white noise in the original time-series translates to a constant lower envelope in the PSD, which is incorporated as a constant term in the fit.The fitting itself is done by minimizing the sum of squared residuals using the Levenberg-Marquart algorithm, which also provides estimates of uncertainties on the fitted parameters using the curvature of the metric around the optimum.
The model of A23 (illustrated in Figure 1), includes aperiodic terms (Equation 5) to represent granulation and super-granulation, and Lorentzian periodic terms (Equation 9) to represent the envelope of the -modes and signals associated with the rotational modulation of active regions.For the latter, they include terms at the fundamental and first two harmonics of the solar rotation period, as well as 1-day aliases of the fundamental and first harmonics.The need to explicitly  3) and the steeper-sloped version of Equation (5), while the blue and pink curves respecitvely show the periodic Lorentzian PSD of Equation ( 9), and the high- SHO term PSD of Equation (13).
include terms to model the aliases is another drawback of modelling the data in the frequency domain.

MODELING THE SOLAR BACKGROUND IN THE TIME DOMAIN USING GPS
GP regression is a natural alternative to PSD fitting when modelling stochastic or quasi-periodic processes in time-series data.Indeed, stochastic process models provided the original motivation for the choice of function used to fit the solar background PSD as described in Sections 2.2 and 2.3.It is thus natural to ask whether GP regression could be used to model the same signals directly in the time domain, side-stepping some of the implementation challenges alluded to in Section 2.4.In particular, the celerite package (Foreman-Mackey et al. 2017, hereafter F17;Foreman-Mackey 2018) provides an efficient implementation of GP regression for certain classes of covariance functions that are particularly appropriate to model stellar signals, and scale well to large datasets.

The celerite SHO term
A generic celerite model is constructed by adding or multiplying together basic building blocks (known as 'terms') with covariance where  and  are complex numbers.These individual terms can then be combined additively or multiplicatively to construct more complex models.
In particular, the celerite package provides a family of kernels that approximate the behaviour of a harmonic oscillator and are particularly appropriate to model a wide range of stellar signals, from aperiodic to strongly periodic.This is implemented in the form of the 'Simple Harmonic Oscillator' (SHO) term, controlled by 3 parameters: an amplitude  0 , a characteristic (undamped) angular frequency  0 , and a quality factor .The PSD of the SHO term1 is: Around  ≈  0 , this PSD approximates a standard Lorentzian or Cauchy distribution, but it falls more steeply (with index 4) at high frequencies.
As of version 2 of the celerite package, an alternative parametrisation is also available, where the user can specify the period , damping timescale  and standard deviation  of the process instead of the angular frequency, amplitude and quality factor: This alternative parametrisation matches the variable names in Sections 2.2 and 2.3, making the relationship between the two classes of models more readily apparent.

Aperiodic components
Clearly, it is straight-forward to reproduce the PSD of Equation 3using a celerite term with real-valued  and .This is also the limiting behaviour of the SHO term when  = 1/2 and  0 → 0.
To reproduce the steeper slope at high frequencies used by more recent studies, however, a natural choice is the SHO term with  = 1/ √ 2. In this case, the PSD of the SHO term becomes: This is exactly equivalent to Equation 5 if  0 ≡ 1/2 and 2  0 ≡  H .

Periodic terms
The PSD of a generic celerite term with real-valued  but complexvalued  is (adapting Equation 11 of F17): If we set  ≡ 2  L Γ,  ≡ 2Γ and  ≡ 2 0 , we obtain which is the sum of two of the Lorentzian functions in Equation ( 9), with characteristic frequencies ± 0 .This is a fairly close match to the periodic components used by A23, but not an exact one, because of the negative frequency term.
On the other hand, we know that -modes are stochastically driven oscillations.Therefore, it makes sense to model them as such, using an SHO term with large , with PSD given by Equation ( 13).If we set  0 =  L Γ2 / 2 0 and  =  0 /2Γ, the resulting process will have the same variance and coherence time as the one defined by the Lorentzian PSD of Equation ( 9).As Figure 2 shows, the PSDs of these two processes look very similar in the vicinity of the peak at the characteristic frequency, but diverge away from it, particularly at high frequencies, where the SHO term gives rise to a steeper power-law slope, with index 4 rather than 2.

The celerite rotation term
The PSD of solar and stellar variability often displays peaks at the first harmonic of the rotation period as well as the fundamental, which need to be modelled explicitly to obtain a good fit.For example, A23 included periodic components at the first two harmonics of the solar rotation period (see Section 2.4) in their final fit.Foreman-Mackey et al. ( 2021) recommends modelling stellar rotation signals using the sum of two under-damped SHO terms: one at the fundamental period, and one at its first harmonic.This is implemented in version 2 of the celerite package as the 'rotation term'.
The rotation term is controlled by 5 parameters: the standard deviation  of the process, the period , two parameters controlling the quality factors of the fundamental and first harmonics (these are specified in such a way that the fundamental has quality factor  1 > 1/2 and the harmonic has  2 >  1 ), and one parameter corresponding to the amplitude of the first harmonic term relative to the fundamental.

TEST METHODOLOGY
In this section, we describe how we simulate time-series datasets with known PSDs including either periodic terms, aperiodic terms or a mixture of the two, and model them both in the frequency domain using analytic PSD functions, and in the time-domain using a GP, to compare the two approaches directly.We restrict ourselves to relatively simple processes with one or two components but also test the impact of irregular time-sampling.

Simulating the time-series
We first simulated a set of regularly spaced time stamps lasting 288 days with 100 points per day, corresponding to a 5 minute cadence.This corresponds to a Nyquest frequency of 452.38 day −1 .We then simulated observations with a known PSD by drawing samples from a celerite GP.For aperiodic components, we use an SHO term with  = 1/ √ 2, the PSD of which is given by Equation 15.For periodic terms, we use an SHO term with large , the PSD of which is given by Equation ( 13).An example of such a simulated time-series for an aperiodic component with  = 1/ √ 2,  0 = 2 (in arbitrary units) and  0 = 5 radians/day (such as might be expected for a supergranulation signal) is shown in the left-hand panel of Figure 3.
Finally, white noise with a standard deviation of  w = 0.5 was added to the simulated time-series to represent measurement uncertainties.
We later investigate the impact of changing the integration time by binning the simulated time-series, and of the baseline by truncating them.
To produce time series with irregular time-sampling, such as might be expected for ground-based RV observations, we used a representative set of time stamps from the publicly available dataset obtained by the HARPS-N solar telescope between 2015 and 2018 (Dumusque et al. 2021).The HARPS-N solar telescope observes the Sun at 5 minute cadence for several hours every day, with occasional interruptions caused by weather and technical issues.We take a 100 day subset of the HARPS-N time stamps and then simulated observations in the same way as described above for the regularly sampled case.As the HARPS-N data has a roughly 5 minute cadence, we expect the Nyquist frequency to be similar to the regular time sampling case.An example of the resulting irregularly sampled time-series is shown in the right-hand panel of Figure 3.
Note that all the variances, standard deviations and amplitudes discussed in this section are in arbitrary units, but could represent normalised flux units (such as parts-per-thousand or parts-per-million) or relative RV units (such as m/s or cm/s).

Estimating PSDs
For regularly sampled time-series, we used the numpy implementation of the Fast Fourier Transform (FFT) to estimate the Discrete Fourier Transform of the observations.The PSD is then given by the squared modulus of the DFT, divided by the frequency step 1/ where  is the duration of the simulated time array.The left panel of Figure 4 shows the FFT-estimated PSD for the simulated time-series shown in the left hand panel of Figure 3.
For irregularly sampled time-series, we used the astropy implementation of the Generalised Lomb Scargle GLS) periodogram (Zechmeister & Kürster 2009) with the normalisation parameter set to 'psd' to estimate the PSD of the observations on a regular grid of frequencies between  = 0 and  = 0.1 cycles/min (corresponding to the approximate effective Nyquist frequency of the HARPS-N solar observations), with frequency step 1/ where  was the duration of the observations.We found it was necessary to divide the output of the periodogram code by the frequency step and by the number of data points to obtain a properly normalised PSD.As the right-hand panel of Figure 4 shows, the PSDs obtained in this manner for irregularly sampled time-series can be strikingly different to the true model used to generate the observations, especially at high frequencies.Naturally, these differences will translate to differences between the true and fitted parameters, highlighting the limitations of PSD fitting for irregularly sampled observations.
Although the Lomb Scargle or GLS periodograms are widely used to evaluate the PSD of irregularly sampled time-series, they do not explicitly account for the window function of the data, which can significantly alter the PDS estimate.A set of finite-duration, irregularly sampled observations of a continuous function  () can be represented as the convolution of  () and a window function (), which is equal to one during each observation and zero elsewhere.If the duration of the observations is short compared to the time-scale of variations in  (), the window function can be approximated by a sum of delta functions located at the times of the observations, {  }.The Fourier transform of the observations is then the convolution of the Fourier transform of  (),  (), and of the Fourier transform of the window function,  ().We can gain some insight into the impact of the window function on the estimated PSD by computing the PSD of the window function, which (for instantaneous observations) is given by (see e.g.VanderPlas 2018): To illustrate this impact, we show in Figure 5 the PSD of a single aperiodic term estimated from data with both regular and irregular time-sampling, together with the PSD of the window function for both cases.This highlights how the window function dominates the irregularly sampled PSD, which differs substantially from that of the true PSD.Note that the peaks in the window function PSD (corresponding to periodicities in the observation times, e.g. 1 day) do not appear clearly in the PSD of the observations, because the PSD of the true signal, which it is convolved with, has a broad spectrum (if the true signal was strongly harmonic, this would not be the case).
Ideally one would seek to deconvolve the PSDs of the signal and the window function, but this is not possible in practice because the  The blue line shows the true PSD of the process used to simulate the data, which is the same in both panels.The grey line shows the PSD estimated using the FFT (left) and GLS (right), respectively, while the black dots show the binned version used in the PSD fit.The green lines show posterior samples from the PSD-fitting MCMC chain, and the pink lines show the same but from the GP fit.Fitted values can be found in Tables 1 and 2.
window function is zero for most values of  (VanderPlas 2018).As a result, most authors that fit PSD models for irregularly sampled data usually include a constant term to account for the extra white noise induced by the window function but do not explicitly account for its structure beyond that.This is one reason why we might expect that modelling stochastic processes in the time domain should be particularly advantageous over doing so in the frequency domain for irregularly sampled observations.

Fitting PSDs
Once the PSD corresponding to a given observation has been estimated, we proceed to fit it, using a PSD model function corresponding to the true PSD of the GP used to generate the data in the first place (namely Equation 15 for aperiodic components, and Equation 13for periodic components).The parameters of the model for each compo-nent are ln  0 and ln  0 , plus ln  for the periodic components only.Fitting for the logarithm of the parameters improves convergence and ensures that they always remain positive.Additionally, a constant term,  w , (also fit in log space) is included in the fit to account for white noise and any under-sampled signals.The constant added to the PSD model is set to  w /2, where  is the total duration of the time-series.
Following A23, we bin the PSD into 100 equally spaced bins in log frequency space before fitting, except when working with periodic terms, where we found it necessary to increase the number of bins to 200 in order to properly resolve the PSD peak.We use only bins containing a minimum of 3 samples (corresponding to frequencies above 0.28 cycles/day), and discard the PSD estimates at lower frequencies.We note that binning PSDs before fitting is normally done in part to ensure that the PSD samples used in the fit are approximately uncorrelated.This would require at least 10, rather than 3, samples per bin.However, this would translate to a  1 and 2. minimum frequency of 1 cycle/day, which is very close to the typical timescale for super-granulation in Sun-like stars, and would therefore preclude the extraction of any information about supergranulation timescales from our simulated time-series.Our decision to use bins with at least  min = 3 samples represents a compromise between maximising the frequency range and minimizing correlation between bins.All bins we assigned equal weight in the fits.Again, this follows the procedure described in A23, where no mention is made of assigning weights or uncertainties to the binned PSD samples.
We perform the fit in log-space, i.e. we minimize the sum of the squared differences between the logarithm of the binned PSD and that of the model function, and assign the same weight to all the binned samples.We also carried out some tests where the fits were performed in linear space and the samples were weighted according to the standard deviation in each bin.However, we found that standard deviation of the PSD samples in each bin systematically underestimates the uncertainty on the binned PSD samples, resulting in poor fits and unrealistically small uncertainties on the fit parameters.
We use Markov Chain Monte Carlo (MCMC) to sample the joint posterior probability distribution over the parameters of the fit, using version 3 of the emcee package (Foreman- Mackey et al. 2013Mackey et al. , 2019)).The walkers are initialised in a tight Gaussian ball (with standard deviation 0.01 dex) around a local optimum found using the minimize function in scipy's optimize module.Log-uniform priors are used for all the parameters, within the interval [−11; 11] in all cases except for the white noise standard deviation log  w , which we restrict to the interval [−2; 2].The number of walkers is set to 4 times the number of parameters for each fit, and the MCMC chains are run for 10 000 steps.
We use emcee's built-in functionality to estimate the autocorrelation length of the chains, in order to assess convergence and select an appropriate burn-in and thinning factor.Provided that the longest auto-correlation length estimate across all parameters is  max ≤ 200 steps, i.e. less than a 50 th of the chain, we consider the auto-correlation time estimates to be reliable and the chains well-converged.We then discard the first 3 max samples as part of the burn-in phase, and thin the remainder of the chains by a factor  max /4.In the few cases where the chains do not converge in 10 000 steps (predominantly when using the PSD-fitting method for multicomponent models), we arbitrarily set the burn-in time to 1000 steps and the thinning factor to 30, but note that the results should be treated with caution.
For illustrative purposes, we evaluated the PSD model for 50 random samples from the truncated, thinned MCMC chains, which are shown in green on the PSD plots in this paper (such as Figure 4).We also use the corner.pypackage (Foreman-Mackey 2016) to display 1-and 2-dimensional posteriors for each fit and compute the median estimates of each parameter and the associated 1- uncertainties.Figure 6 shows the corner plots obtained in this manner for the time-series shown in Figure 3, for which the PSDs are shown in Figure 4.

GP-modelling
We also modelled the simulated time-series directly in the time domain, using a GP model of the same form as the one used to generate them.The fit parameters, priors and details of the optimization and posterior sampling procedures were identical to those described in Section 4.3, the only difference being that the figure of merit of the fit was the GP likelihood, rather than the  2 of the PSD fit.
The results of the GP fits are shown in pink on the PSD and corner plots (such as Figures 4 and 6), while the results of the PSD fits are shown in green on the same figures, and the true model in blue.

RESULTS
We are now ready to compare the results obtained for these simulated time-series with both the PSD-fitting and the GP-modelling approach, varying the type of model used and the time-sampling.

Single-component case
We start with the simplest case, namely the single, aperiodic, (super-)granulation-like component illustrated in Figures 3, 4 and 6.The fitted parameters for all the tests described in this section are reported in Tabels 1 and 2.
Even in the regularly-sampled case, the parameter estimates derived via PSD fitting already deviate from the true values at the 2-3 level.This is not unexpected for the white noise component: the finite duration of each observation and of the total time-series result in additional white noise in the PSD coming from the window function.However, it is concerning to see a significant deviation for the other parameters, even in this idealised case.The discrepancy between injected and recovered parameters is even starker in the irregular case, as expected due to the more complex window function.The PSD fits are also quite sensitive to somewhat arbitrary modelling choices such as the number of PSD bins used.
We then explored the effect of both increasing the exposure times (i.e.binning the time-series) and reducing the duration (truncating the time-series) to test the sensitivity of the results to these factor.As shown in Figure 7, this degraded the accuracy of the results further, especially in the irregular case, where all parameters were overestimated, even when the 'knee' of the PSD appears well sampled by eye.In the regular case, the PSD fits, while generally consistent with the the true value, have larger uncertainties associated with them.
By contrast, the GP-modelling results are both more precise and consistent with the truth in all these cases, though the precision worsens (as expected) when the time-series is binned, truncated or irregularly sampled.We repeated this test multiple times with a range of  0 and  0 values and found that the results were substantially the same.
We then performed the same set of tests for periodic models, with essentially similar results as shown in Figure 8.The corner plots associated with this set of tests can be found in Figure 8.
Finally, we considered a case consisting of two periodic terms, with frequencies differing by a factor 2, as implemented in the celerite rotation term mentioned in section 3.4.The results are shown in Figure 9.While both PSD-and GP-fitting were successful in the regularly sampled case (left panel), the PSD fitting approach entirely fails to recover the true shape of the PSD in the irregular case.By contrast, the GP fit still gives tight and accurate constraints on the underlying PSD and its parameters.
In all cases, the GP-modelling substantially outperformed the PSDfitting method.Not only do the GP fits result in smaller uncertainties, but they are also substantially more accurate, and generally consistent with the ground truth, unlike the PSD-fitting results.We also note that, as the time-sampling was degraded, MCMC convergence also became increasingly slow in the PSD-fitting case.

Multiple components
We next considered two-component models, to test our ability to disentangle between the components and recover their parameters correctly.We ran two such tests: one with one aperiodic and one periodic component, and one with two aperiodic components.All cases were run multiple times to ensure the repeatability of results.The PSDs plots for these tests can be found in the body of the text in Figure 10, while the corner plots can be found in Appendix A. It should be noted that in all cases shown below where the MCMC did not converge is due to highly multi-modol posterior distributions.

Model Comparison
We started by considering a low-frequency, high-amplitude term combined with a higher-frequency, lower amplitude, high- periodic term.In the context of asteroseismology, the aperiodic term could represent granulation, and the periodic term the envelope of the modes.The results of the PSD fit for these tests are shown in Figure 10, and the corresponding corner plots in Appendix A. In both the regular and the irregular time sampling cases, the GP once more outperformed the PSD-fitting method.It is worth noting that, in the regular time sampling case, the PSD-fitting method struggles to   1 and 2. constrain the quality factor of the oscillations (this is partly due to the binning of the PSD prior to fitting).Meanwhile, in the irregular case, the oscillations peak is entirely absent from the PSD estimate which naturally leads to poor PSD-fitting results (and the MCMC fails to converge), while the GP fits are much less severely affected: the uncertainties are increased but the results are still consistent with the ground truth.
An important question is whether the PSD and GP fitting methods can correctly identify the number and nature of the components present in the data.To check whether this is the case we repeated the fits of our 2-component time-series and PSDs using different models, including: (i) 1 periodic component and 1 aperiodic component (true model), (ii) 1 aperiodic component only, (iii) 2 aperiodic components, (iv) 2 periodic components.In each case, we compared the Bayesian Information Criterion (BIC) of the best-fit model.The results of this model comparison are reported in Table 3.Both the single aperiodic component model (case ii) and the two aperiodic component model (case iv) where preferred when using PSD-fitting, although the single aperiodic component more so, while the true model (case i) is preferred when using the GP.
We then considered two aperiodic terms with more similar decay times, representative of super-granulation and granulation.The results of these tests are shown in Figure 11, with the corner plots found in Appendix A. Here, the "true" parameters are the same order of magnitude as the ones reported by Al Moulla et al. (2023) for granulation and super-granulation in the solar case.It can be seen that the GP fits result in timescales consistent with the ground truth for both granulation and super-granulation, but the PSD fitting approach does  1 and 2. not, especially for the granulation-like component.This implies that the timescales reported by A23 should be treated with caution since they were derived using a PSD-fitting procedure very similar to that used in the present work.We also used the BIC to test whether a 2component model was preferred over a single-component case.The results are reported in Table 3.The 1-component model is favoured in the PSD case, however the preference is not statistically significant.
We then proceeded to reduce the frequency separation between the two components, to check when each method stops being able to distinguish between them.The results of these tests are shown in Table 4, where the BIC results for models with either 2 aperiodic components or 1 aperiodic component are shown.When using a GP, the true model is preferred, and thus the two components are correctly identified, though the preference over a single component model is statistically significant only down to a frequency ratio of  1  2 ∼ 4 .By contrast, the PSD fit could only separate the two components for the largest separation.It should be noted that the preference for a 2 aperiodic model in the irregular GP case with  1 \ 2 = 3 is not statistically significant, and the result may be different depending on the initial time series sample.
We also used the BIC test for the signal component celerite rotation term to see if two component models were preferred over the signal component model in this case.The results are reported in Table 3.A single aperioidc component is preferred in the PSD case.

Parameter Accuracy
Tables 1 and 2 show the recovered parameters compared to the true values for all 5 cases studied, for both regular and irregular time sampling.Cells for which a derived parameter was more than 2 away from the true value are highlighted in gray.In every case discussed above, the GP approach outperformed the PSD-fitting approach.The GP-derived parameters are always close to the true values.However, in a few cases, the uncertainties appear under-estimated, so that the 2 credible intervals do not quite overlap with the true value.This tendency to under-estimate uncertainties is particularly pronounced   1 and 2. for the white noise parameter and the quality factor of any periodic terms, but is also sometimes observed for the characteristic frequencies.The performance of the PSD-fitting method degrades very noticeable when using irregular sampling, whereas the GPfitting remains fairly robust, except in the case with two aperiodic components, where both methods failed to recover the granulation frequency and amplitude within 2 − , though the GP values are still of the correct order of magnitude.

DISCUSSION AND CONCLUSIONS
Having compared the behaviour of PSD-fitting and GP regression on simulated datasets with a known ground truth, we now proceed to discuss the implications of our results, both for RV planet searches (the original motivation for this work) and more widely, and outline perspectives for future work.

GP regression versus PSD fitting
In this work, we have shown that the analytic functions commonly used to model stochastic variability in the light curves and RV observations of the Sun and Sun-like stars can be reproduced, either  1 and 2. exactly or approximately, using Gaussian process models.We then showed,q2 using a series of tests on simulated data, that GP regression in the time domain systematically leads to better results than estimating and modelling the PSD in the frequency domain.
While the GP method leads to more accurate and precise results when applied to regularly sampled data, the difference in performance becomes even more pronounced in the irregular sampling case, where the PSD-fitting results deviate dramatically from the truth, while the GP results are generally consistent with it within 2.The main reason for this is that GP regression does not rely on PSD estimates, which by definition assume a specific model for the signal whose PSD is being estimated.In the case of the GLS periodogram, this model consists of a single sinusoidal signal, plus white noise.When the time sampling is close to regular, individual sines and cosines with frequencies between 1/ and 1/2 are approximately mutually orthogonal, and the periodogram closely approximates the true PSD, even if the signal is more complex than a single sinusoid.However, when the sampling becomes sparser, this orthogonality breaks down, and the periodogram is dominated by the window function.This is particularly problematic for stochastic signals, which have a broad intrinsic PSD, and thus deviate strongly from the simple model implicitly assumed by the GLS periodogram.The GP regression approach avoids this problem by obviating the need for a PSD estimate altogether.
It should be noted that the computational cost of GP regression is significantly larger than that of the PSD-fitting approach.For example, for the regularly sampled, single aperiodic component case on a 2020 Mac Mini with an Apple M1 chip, the MCMC sampler completed 580 iterations/second for the PSD-fitting versus 92.5 iterations/second for the GP fitting case, a difference of a factor > 6.Nonetheless, as the evaluation of the GP likelihood scales linearly with the number of data points for celerite models, and the posterior sampling is easy to parallelize, the computational cost is an obstacle that can be overcome, even for fairly large datasets.
Therefore, we recommend that GP regression in the time domain should be used, wherever possible, in preference to PSD-fitting, to model stochastic processes in astrophysical time-series datasets.One limitation of the GP-fitting method as implemented in this work, which emerged during our tests, is that it has a tendency to underestimate the uncertainties on the parameters slightly (by a factor of order 2).At this stage, the reason for this tendency is unknown, but we recommend treating these uncertainties with a degree of caution until this behaviour is better understood.

Implications for RV planet searches
The advantages of GP regression over PSD-fitting are particularly relevant to ground-based RV observations, which are always irregularly sampled.Indeed, GP regression is already widely used to model and mitigate activity signals in RV datasets (see e.g.Haywood et al. 2014;Rajpaul et al. 2015;Barragán et al. 2022).The results presented in this paper, and in particular the finding that GPs can robustly identify multi-component signals in irregularly sampled data, indicate that GP regression may also have a role to play in the characterisation and mitigation of other stellar signals such as super-granulation.
Besides the aforementioned GP methods, numerous other approaches to model activity signals at the spectral or line-profile level have been developed in recent years (see e.g.Collier Cameron et al. 2021 ;Cretignier et al. 2022Cretignier et al. , 2023;;Liang et al. 2023).Using a combination of these approaches, it is widely expected that future Extreme Precision RV (EPRV) surveys aiming to detect Earth analogues, such as the Terra Hunting Experiment (THE, Hall et al. 2018), which will benefit from unprecedented time-sampling and high signal-to-noise ratio, will succeed in mitigating activity signals in Sun-like stars from several m/s down to the sub-m/s level.
The next major challenge is thus super-granulation (Meunier & Lagrange 2020) for which no effective correction methods have been developed yet.Indeed, Lakeland et al. (2024) recently analysed the RV variability of the quiet Sun as observed by HARPS-N (after accounting for the effects of active regions using SDO data, following the method developed by Haywood et al. 2016 andMilbourne et al. 2019).During the low-activity part of the Solar cycle, they found that the quiet photosphere contribution to RV variability is as important as, if not larger than, that from active regions.
A natural continuation of this work will thus be to use the GP methodology we have proposed here to revisit recent studies aiming to quantify the timescales and amplitudes of variability due to (super-)granulation on the Sun, including A23 and Lakeland et al. (2024) 2024) used structure functions, a non-parametric approach that, like GP regression, can be applied directly to irregularly sampled data, but doesn't explicitly use a generative model, and it would be valuable to gain a better understanding of the relationship between structure functions and GP models.One potentially interesting avenue for future progress may also be to search for and characterise any non-Gaussianity in stellar RV signals.

Asteroseismology and planetary transits
The PSD models discussed in Section 2 first arose in the context of helioseismology, and are widely used in asteroseismology, whether using RV observations or photometry.While the "background" in this context is merely a nuisance signal to the oscillations, it is conceivable that there may be some advantage to be gained by modelling this background directly in the time domain using GPs rather than by evaluating and fitting a PSD.This may be particularly helpful for later-type main sequence stars, for which the oscillation frequencies are higher and the amplitudes lower than earlier type or evolved stars (see e.g.Chaplin & Miglio 2013).As we have shown in this work, multi-component GP models can be used to detect (quasi-)periodic signals in the presence of correlated noise, even when no clear peak is visible in the PSD of the data.Such models, consisting of one or more aperiodic terms to represent (super-)granulation and a single high quality-factor term to represent the envelope of the -modes, may help push the detection threshold for -modes beyond early- spectral types, which is the current limit.
This type of model has already been used to model granulation signals and oscillations simultaneously in order to characterise planetary transits around red giants (Grunblatt et al. 2017), and may likely prove important for precise measurements of planetary radii from transits (Barros et al. 2020), e.g. in the context of the PLATO space mission (Rauer et al. 2016).

Quasi-Periodic Oscillations and accreting systems
Accreting systems on all scales, from compact binaries to Active Galactic Nuclei (AGN) display Quasi-Periodic Oscillations (QPOs) (see e.g.Reig 2011;Ingram & Motta 2019).QPOs are time-variable quasi-periodic signals observed at X-ray, UV and optical wavelengths, and are usually modelled in the frequency domain as a sum of Lorentzian functions (Belloni et al. 2002;van Straaten et al. 2002).The parameters of these models, including the frequency and amplitude of any quasi-periodic components as well as the frequency at which aperiodic components display breaks in their PSDs, are used to interpret QPOs and understand their origin and diversity.There is every reason to expect that standard QPO modelling procedures are subject to the same limitations as outlined in this work for stellar signals in RV or photometry, and that GP regression could be a more robust alternative, particularly for shorter or more sparsely sampled time-series.Indeed, GP regression has begun to be used more frequently in the context of QPOs and AGN variability (see Aigrain & Foreman-Mackey 2023, and references therein).

Figure 1 .
Figure 1.Composite model of the PSD of the solar background, from Al Moulla et al. (2023).The full model (solid purple line) is the sum of quasiperiodic components corresponding to the envelope of the solar oscillations (green dashed line), and to the solar rotation period and its first two harmonics (red dashed line), aliases thereof (black dashed line), plus aperiodic components corresponding to granulation (blue dashed line) and super-granulation (orange dashed line), plus white (photon) noise (cyan dashed line).

Figure 2 .
Figure 2. Comparison between the different types of PSD functions used in this paper.The orange and green curves respectively show the aperiodic Lorentzian PSD of Equation (3) and the steeper-sloped version of Equation (5), while the blue and pink curves respecitvely show the periodic Lorentzian PSD of Equation (9), and the high- SHO term PSD of Equation (13).

Figure 3 .
Figure 3. Left: Example of a simulated, noise-free, regularly sampled super-granulation time series created using a celerite SHO term GP with  = 1/ √ 2,  0 = 2 and  0 = 5 2  cycles/day.Right: noisy, irregularly sampled time-series generated for the same process with time-stamps from the HARPS-N solar telescope.

Figure 4 .
Figure 4. PSD estimates and fits for a single, aperiodic component with  0 = 2 and  0 = 5/2 , with regular sampling (left) and irregular sampling (right).The blue line shows the true PSD of the process used to simulate the data, which is the same in both panels.The grey line shows the PSD estimated using the FFT (left) and GLS (right), respectively, while the black dots show the binned version used in the PSD fit.The green lines show posterior samples from the PSD-fitting MCMC chain, and the pink lines show the same but from the GP fit.Fitted values can be found in Tables1 and 2.

Figure 5 .Figure 6 .
Figure 5. Impact of the window function on the PSD of a for a single, aperiodic component with  0 = 2 and  0 = 5/2 , with regular sampling (left) and irregular sampling (right).The PSD estimate of the window function is shown in green, the estimated PSD of the simulated data is shown in pink, while the true PSD is shown in blue

Figure 7 .
Figure 7. Recovered versus true parameters for a single aperiodic component with regular time sampling (top row) and irregular sampling (bottom row), as the binning is increased and the duration of the time-series is truncated.Green symbols show the results of the PSD fits, while pink symbols show those of the GP fits.The true parameter values are shown by the vertical blue lines.

Figure 8 .
Figure 8. PSD estimates and fits for a single, periodic component with  = 10,  0 = 0.1 and  0 = 20/2 , with regular sampling (left) and irregular sampling (right).The colour coding is the same as Figure 4. Fitted values can be found in Tables1 and 2.

Figure 9 .
Figure9.PSD estimates for the celerite rotation term model in the regular case (left) and irregular case (right), using the same colour convention as in Figure6.The injected parameters are as follows :  = 50,  = 0.5,  0 = 20,  = 15,  = 0.8.Fitted values can be found in Tables1 and 2.

Figure A1 .
Figure A1.Corner plots for the case of a single, periodic component with  = 10  0 = 0.1 and  0 = 20/2 , with regular sampling.PSD fits are shown in Figure 8.The GP MCMC posteriors are in pink, which the FFT/GLS posteriors are in green.

Figure A2 .
Figure A2.Corner plot for the case of a single, periodic component with  = 10  0 = 0.1 and  0 = 20/2 , with irregular sampling.PSD fits are shown in Figure 8.The GP MCMC posteriors are in pink, which the FFT/GLS posteriors are in green.

Figure A3 .
Figure A3.Corner plots for the celerite rotation term model in the regular time sampling case.The injected parameters are as follows :  = 50,  = 0.5, 0 = 20,  = 15,  = 0.8.The GP MCMC posteriors are in pink, which the FFT/GLS posteriors are in green.

Figure A4 .
Figure A4.Corner plots for the celerite rotation term model in the irregular time sampling case.The injected parameters are as follows :  = 50,  = 0.5, 0 = 20,  = 15,  = 0.8.The GP MCMC posteriors are in pink, which the FFT/GLS posteriors are in green.

Table 1 .
Injected and fitted parameters for all cases with regular time sampling.Parameters more than 2 away from the true value are highlighted in gray.

Table 2 .
Injected and fitted parameters for cases with irregular time sampling.Parameters more than 2 away from the true value are highlighted in gray.

Table 3 .
ΔBIC results for multiple component cases in the irregular time sampling regime.Positive values mean the model is preferred over the true model, negative results mean the opposite.The best fit model is highlighted in gray.

Table 4 .
BIC values for 2 aperiodic components with regular and irregular time sampling as the two components are brought closer together.The only initial parameters that are changed are  1 and  2 .Positive values mean a 1 aperidoic component model is preferred over the true model, these cases are highlighted in gray, negative results mean the opposite.
, as well as in other stars (see for example Sulis et al. 2023, who compared granulation signals in photometry and RV using simultaneous observations with CHEOPS and VLT/ESPRESSO).While A23 and Sulis et al. (2023) used PSD-fitting, Lakeland et al. (