© ESO 2021

1 Introduction

High-resolution Doppler spectroscopy has been instrumental in finding and confirming the existence of planetary companions orbiting stars. Modern instruments like HARPS (Mayor et al. 2003), CARMENES (Quirrenbach et al. 2018), and ESPRESSO (Pepe et al. 2010) attain radial velocity (RV) precisions in the m s⁻¹ domain and better. This is effectively comparable to astrophysical variations even in the case of stars with moderate to low magnetic activity levels. Therefore, a precise understanding of intrinsic stellar phenomena is mandatory in order to expand the search for low-mass, rocky exoplanets.

In the cm s⁻¹ domain, stellar effects include short-term variability induced by surface granulation and stellar oscillations (see e.g. Dumusque et al. 2011; Milbourne et al. 2019). Furthermore, magnetic phenomena produce cool spots and hot faculae co-rotating with the stellar surface and with lifetimes of the order of a fraction to several rotational periods. These phenomena introduce RV signals that can distort, mimic, and even hide the Keplerian signature of a planetary companion. The filling factor (i.e. the fraction of surface covered with spots) and the location of the spots can also change over time producing long-term magnetic cycles of the order of years to decades (Suárez Mascareño et al. 2018). Magnetic effects are found to be even larger for cooler stars (M type), which are the prime targets for several RV surveys conducted over the last decade (e.g. Pepe et al. 2004; Affer et al. 2016; Reiners et al. 2018) as they offer a shortcut to finding temperate rocky planets. An overview of the impact of various stellar effects on RV measurements is provided in Dumusque et al. (2012).

The scientific community has put a lot of effort into the treatment of star-induced RV variability. Since its first use by Haywood et al. (2014), and in particular when the first results of the RV challenge by Dumusque et al. (2017) were published, Gaussian process (GP) regression (Rajpaul et al. 2015; Roberts et al. 2013) has become one of the most commonly used tools to model and mitigate correlated effects in RV time series arising from stellar activity. The flexibility of GP algorithms makes them suitable for accounting for the effects of rotating spots on the data using a kernel function defined by a few hyper-parameters and the application of optimisation techniques (e.g. maximisation of the likelihood function). However, the association of the mathematical kernel hyper-parameters with true physical phenomena has not been clearly discussed and evaluated.

We have developed the StarSim code (Herrero et al. 2016; Rosich et al. 2020) to investigate the impact of stellar activity on time-series data using a physical approach. In its forward model functionality, StarSim assumes a stellar surface configuration, which is mapped onto a grid of elements either representing the immaculate photosphere, a cool dark spot, or a hot facula surrounding the spot. Observables (photometric and spectroscopic) at a certain epoch are predicted by integrating the contributionof all elements over the full visible stellar surface, including geometric effects, convective shifts, line asymmetries, characteristics of the spots and faculae, and of the host star, including the possibility of differential rotation.

In this study, we create synthetic RV time-series data of a rotating spotted star using different realistic, evolving spot distributions and the StarSim code. To interpret the simulated data, we analyse the auto-correlation function (ACF) and fit different commonly used GP kernels. We are then able to study the connection of the introduced physical parameters with the posterior distribution of the GP hyper-parameters. We present the mathematical background to the likelihood function and the ACF in Sect. 2. In Sect. 3, we introduce the StarSim RVs of the differentspot distributions and analyse their ACFs, and in Sect. 4, we apply GP regression using the different kernels and evaluate their performance. The results are summarised in Sect. 5.

2 Auto-correlation function and kernels

2.1 Likelihood function

The goal of a regression exercise is to include a parameterised representation of the data (model) that best fits the variability measured in the observations. In the current literature, this is done by maximising a merit function called the likelihood function $\mathcal{L}$. The likelihood function represents the probability distribution of the data (a time series v and a function of time t in our case) to occur given a specific model, which can contain a number of adjustable parameters. The likelihood function arises from the common assumption that the residuals r_i, that is, the difference between the observation v_obs(t_i) and a parameterised model v_mod(t_i), are drawn from a joint multivariate Gaussian distribution and is usually defined as $$\mathcal{L}=\frac{1}{{(2\pi )}^{N/2}}|C{|}^{-1/2}\text{exp}\left(-\frac{1}{2}{\sum }_{i=1}^N{\sum }_{j=1}^N{r}_i{r}_j{\text{C}}_{ij}^{-1}\right),$$(1) $$r_i=v_{\text{obs}}(t_i)-v_{\text{mod}}(t_i),$$(2)

where N is the number of observations, and C is the covariance matrix between data residuals. This covariance matrix can be written as $$C_{ij}=({\sigma }^2+{ϵ}_i^2){\delta }_{ij}+K(t_i-t_j),$$(3)

where ϵ_i is the uncertainty on the measurement v(t_i), σ represents uncorrelated random noise (jitter), and δ_ij the Kronecker Delta, and covariances between data points are specified by the so-called stationary kernel function K (or simply kernel hereafter); for a detailed discussion see for example Baluev (2013), Anglada-Escudé et al. (2013), and Ribas et al. (2018). For numerical reasons, it is more practical to work with the logarithm of the likelihood function $\ln \mathcal{L}$, which will also be maximal for the same parameter values as the likelihood function. The maximum value of a likelihood function that depends on a number of adjustable parameters n_p will always improve when adding more. When comparing models with different numbers of parameters, the Bayesian information criteria (BIC) is often used as a way to apply a penalty to more complex models, and can be computed as $$\text{BIC}=n_{\text{p}}\ln N+2\ln \mathcal{L}.$$(4)

Although this indicator does not account for Bayesian priors nor multi-variate local minima, it is used below to additionally evaluate the relation between different models with increasing complexity.

2.2 Relation between physical models and correlated noise

It is important to note that, in this picture, our interpretation of the data can be included in two different parts of the likelihood function. Firstly, adjustable parameters can be in the directly parameterisable model v_mod, and secondly in the covariance matrix C and its kernel K. Whether a signal is better described as part of the v_mod-model (e.g. detection of a planetary signal) or inside the covariance matrix C is often used as a diagnostic to decide if the signal is induced by a real planet or is caused by stellar activity. It has been acknowledged that there must be a certain amount of degeneracy between signals described by v_mod and those in the covariances, which is the source of all sorts of controversial claims. In this section, we re-derive a basic description of correlations between data sets and show the relationship between these correlations and the signals that generate them, which should provide us with useful information for subsequent discussions.

The covariance between the values of v = v(t) and the same variable at a later time v′ = v(t′) is defined as $$C\left[v,v^{\prime }\right]=E\left[\left(v-E[v]\right)\cdot \left(v^{\prime }-E[v^{\prime }]\right)\right],$$(5)

where the expected value E of a variablev is defined as $$E[v]=Z^{-1}{\displaystyle \int v}(t)\text{d}t,$$(6)

and Z = ∫ dt normalises the integral for a relevant time interval. Formally, this integral should expand from − ∞ to ∞, but it can be defined over a finite interval in which case Z = T, where T is the time-span of the observations, as shown in subsequent examples. Assuming that the function v(t) is non-zero over a finite interval, the expected values satisfy E[v] = E[v′] = γ, which is a constant (average of the time-series). By explicitly computing the products before the expected values, substituting γ into Eq. (5), and defining K(t − t′) = C[v, v′] we obtain $$K(\tau )=E[v\cdot v^{\prime }]-{\gamma }^2,$$(7) $$\tau =t-t^{\prime },$$(8)

which represents the correlations induced by the physical processes described by v(t) as a function of the time-lag τ. Therefore, this is the same object as the kernel function in Eq. (5). From this expression we can immediately note that a non-zero value for the expected γ² will contribute to the kernel unless it is fitted as a free parameter in the explicit part of the model (e.g. a constant offset in v(t)). Another fact that becomes obvious when looking at this expression is that the kernel is directly connected to the physical process v(t) through the auto-correlation function of the generating signal, which can be explicitly computed using $$K(\tau )=\frac{1}{Z}{\displaystyle \int v}(t)v(t-\tau )\text{d}t,$$(9)

and can be computed for any generating function v given from models, or can be approximately computed from observations themselves (the integral becomes a discrete sum). This feature allows the development of physically motivated kernels by comparing simulated observations (signals generated by spots on rotating stars) to the kernels typically used in the literature.

Before investigating realistic cases, we illustrate the process of how to compute exact kernels in two simple but representative situations. The two generating functions discussed here – a linear trend and a sinusoidal signal – may not necessarily be used in an optimisation algorithm, as they can be more easily added in the explicit model, but we think they provide useful insights into the connection between physics and kernels, and provide clues as to the sources of degeneracy often encountered in the literature.

Case 1 – A linear trend

In this case, the generating function corresponds to v(t) = a + bt, where a is a constant offset and b is the slope. Substituting in Eq. (9), the corresponding kernel becomes $$K_{\text{lin}}(\tau )=\frac{1}{Z}\left[{\displaystyle {\int }_0^T\left(a+bt\right)}\left(a+b(t-\tau )\right)\text{d}t\right]$$ $$=a^2+abT+\frac{b^2}{3}{T}^2-\left(ab+\frac{b^{2}T}{2}\right)\tau ,$$(10)

which only contains constants and linear terms with τ. Therefore, asuitable parameterised kernel could simply be expressed as $$K_{\text{lin}}(\tau )=\alpha +\beta \tau ,$$(11)

where α and β would be related to the slope and offset of the time-series in a complicated way following Eq. (10).

Case 2 – A sinusoidal signal

The generating function has now the form v(t) = κsinwt. Again substituting in Eq. (9), the resulting kernel is $$K_{\text{sin}}(\tau )=\frac{1}{Z}{\displaystyle \int {\kappa }^2}\text{sin}\left(\omega t\right)\text{sin}\left(\omega (t-\tau )\right)\text{d}t$$(12) $$={\kappa }^2\text{cos}(\omega \tau )+O[{\kappa }^2/N],$$

where κ² (square of the semi-amplitude) and ω (angular frequency corresponding to a period P = 2π∕ω) would be the adjustable free parameters. To derive this expression we used trigonometric identities, and assumed that the integrals of the trigonometric functions over the time-domain go towards zero when several cycles N ~ T∕P are covered, as reflected in the last term in Eq. (13). This latter assumption may not hold when examining signals with periods close to the baseline of the observations, but the expression still provides an approximate idea of the functional shape of the kernel. The same calculation can be done adding an arbitrary phase ϕ as v(t) = Asin(ωt + ϕ), reaching an identical result independently of the value of ϕ. This leads to the very general conclusion that sinusoidal signals always produce cosine-like correlations. This example illustrates that, for a general signal, κ² is an approximation of the quadratic sum of all the terms in the Fourier decomposition of a signal at τ = 0. Therefore, it is useful to name this new quantity ${\kappa }_0=\sqrt{K(\tau =0)}$, which essentially represents the combined amplitude of the signals causing correlations.

Correlations affecting precision Doppler measurements are not necessarily as simple as these cases. Moreover, observed correlations are likely to be the result of the superposition of several generating functions. One can easily prove that, while physical signals are additive, their resulting correlations do not combine linearly. As an example, if a time-series contains two signals x(t) and y(t), an appropriate physically motivated kernel should consist of

Fig. 1 Auto-correlation functions of 144 RV data points of a rotating star (Prot = 25 days) with one single central spot. We show the dependence of KACF on stellar inclination angles i = 30, 45, 60, and 90 deg with the colours as indicated. The dependence on the radius of the dark spot R = 5, 10, 20, and 40 deg is illustrated by the different symbols as indicated. We show the repetitive part of the ACFs with τ < Prot.

Fig. 2 Square root of the maximum ACF value of our test-case M dwarf at τ = 0, κ0, as a functionof the maximum spot filling factor, ${\mathcal{F}}_{\text{max}}$ (top panel), and of the sine of the stellar inclination angle, sini (bottom panel). The curves are shown with colours as indicated for the different inclinations and different spot radii. The black dashed lines are the best fits on the curves.

$$K_{xy}(\tau )=\text{C}[x,x^{\prime }]+\text{C}[y,y^{\prime }]+\text{C}[x,y^{\prime }]+\text{C}[y,x^{\prime }],$$(13)

which contains the ACFs of both signals (first two terms on the right-hand-side), plus the covariances of the cross-terms (two last terms), thus resembling the effect of physical interference. The same principle generalises to more than two signals. That is, the total covariance of a generating signal always contains the auto-correlation terms first (which quadratically contribute to κ₀), and all pairs of cross terms which tend to chaotically cancel out at large time-lags τ if they are not strictly periodic.

3 Synthetic radial velocities of an M dwarf as case study

To generate test data, and as such stars are still the prime targets in the search for rocky exoplanets, we adopt an M-dwarf star with T_eff = 3750 K, log g = 4.5, and [Fe/H] = 0.0 rotating non-differentially with a rotation period of P_rot = 25 days. We assume dark spots to have a temperature 300 K cooler than the photosphere, roughly following Mallonn et al. (2018), and we consider no faculae. Spots in the StarSim code are defined by five characteristic parameters, namely longitude, colatitude, radius, time of appearance, and lifetime. We first run a suite of experiments to obtain a basic understanding of the ACFs induced by single spots (Sect. 3.1), and then we proceed with multi-spot configurations (Sect. 3.2). Commonly used kernels are fitted to the simulated data and discussed in Sect. 4.

3.1 Single-spot configurations

In a first numerical experiment we create RV time-series as induced by a single central dark spot on the stellar surface with a variety of stellar inclinations and spot radii in order to study their effect on the ACF. Radial velocity sets consist of 144 epochs over a 100-day time-span (every 1000 min), stellar inclinations of i = 30, 45, 60, and 90 deg, and spot radii of R = 5, 10, 20, and 40 deg. We assume no noise or measurement uncertainty. In Fig. 1, the repetitive parts (τ < P_rot) of the ACFs of these datasets are shown in different colours for the various values of i, and for the different R with the different symbols, as labelled. The correlation decreases as both parameters decrease. In general, K_ACF is maximal at multiples of P_rot (25 days) and close to zero at P_rot∕2. Furthermore, the curves are defined by two characteristic symmetric lobes of negative correlation. The time-lags for which a minimum of K and K = 0 are achieved and the depth of those lobes are defined by the description of convection and limb darkening of the StarSim code.

In Fig. 2 we show κ₀, linearly correlated with the RVs of our test case, as a function of the maximum projected spot coverage ${\mathcal{F}}_{\text{max}}$ (top panel) and sini (bottom panel) for each of the spot radii and inclination angles as indicated, and including theoretical values for i = 0 deg and R = 0 deg. The maximum filling factor of the spots ${\mathcal{F}}_{\text{max}}$ is calculated for the different spot radii R, inclinations i, and colatitudes θ (= 90 deg) as

Fig. 5 Graphic representation of the three kernels used in this study: the quasi-periodic (QP) kernel as implemented in the george code (top), the simple harmonic oscillator (SHO) as implemented in the celerite code (middle panel), and the quasi-periodic and cosine (QPC, bottom) kernel introduced in this work, for P = 25 days. The x-axis is the time-lag τ between two data points, and the vertical axis is the normalised value of the correlation represented by the kernel. We show the kernels in red for λ = 25 days (w = 0.31) and black forλ = 75 days. For the QP kernel, the light grey area marks 0 < w < 1, whereas the dark grey area marks the interval 0.2 < w < 0.5, which is more physically motivated. The QPC kernel has h2∕h1 = 1.

$${\mathcal{F}}_{\text{max}}=2(1-\text{cos}R)(\text{cos}\theta \text{cos}i+\text{sin}\theta \text{sin}i).$$(14)

The dependence of the RV amplitude on ${\mathcal{F}}_{\text{max}}$ is close to linear, with only the increasing projection distortion of larger spots causing some slight curvature. However, the dependence on sini is more parabolic. For an analytic representation of our test case, we could obtain a reasonable fit with a product of a second-order polynomial in ${\mathcal{F}}_{\text{max}}$ and a fourth-order polynomial in sini (black dashed curves in Fig. 2). However, we note that the details of such a fit depend on the spectral properties in a convoluted way. Our recommendation when analyzing a dataset from a star with at least one suspected spot is to perform a similar quick simulation to gauge the magnitude of κ₀ in the RV time-series. As shown below, the value of this κ₀ shall match the amplitude parameters of all of the used kernels.

3.2 Multi-spot configurations

To study situations closer to the more complex spot patterns of real stars, we calculate RV time-series of a variety of spot distribution models. The different evolving spot maps are used to create 576 synthetic RV data points with values every 3.47 days (5000 min) for a total time of 2000 days and without assuming noise or measurement uncertainties. By choosing this window function, we avoid complicated sampling effects, and have a realistic time baseline for single-instrument observations, and a sufficiently large number of data points. The inclination of the test star in this case is set to 90 deg (i.e. equator on).

The lifetime of all spots is fixed to T_spots = 100 days, because this parameter is crucial for the GP regression and we aim to calibrate it in this study. The spot radius R is selected, together with N_spot, to create a realistic σ_RV ~ 5 m s⁻¹. We also include a long-term variation, mimicking a magnetic cycle, by linearly varying the number of spots from N to N/2 to N every 1000 days in intervals of 100 days (see Table 1). In each interval i, each of the N_i spots appears at a random day (t_0,i < t_i < t_0,i + 100 day) and grows with a rate of 1 deg/day, allowing for a smooth spot evolution over the 2000 days. To explore the effects produced by different configurations, we create five different maps and RV sets for each spot distribution.

We consider four different configurations for the distribution in longitude and colatitude of dark spots, in line with the possibilities suggested by the relevant literature on M dwarfs (see e.g. Donati et al. 2008; Andersen & Korhonen 2015):

• Configuration ONE. This configuration is inspired by the observations of the fast-rotating early-M dwarf EY Dra by Vida et al. (2010), which shows a distribution of dark spots around one active longitude. To create such a spot map, we distribute a maximum of 40 small spots of 3 ± 1 deg radius, resulting in a total of 600 spots for a 2000-day baseline. We imitate the behaviour of solar spots, that is, the butterfly diagram, by distributing them from colatitudes of 95–125 deg and from 55 to 85 deg. Also, we place them closer to 90 deg when the number of spots is lower. The scheme for the spot distribution is shown in Table 1. One of the five created maps is shown in the top panels of Fig. 3 for the highest spot coverage at t = 1000 days (left panel) and the lowest coverage at 500 days (right panel). The active longitudes of the spots are normally distributed around 180 deg (σ_n = 60 deg).

• Configuration POL. This configuration is characterised by large spots at high latitudes (polar regions) of a star; second row of Fig. 3). This is observed in the K-type dwarf AB Dor by Jeffers et al. (2007) and using tomographic imaging of active M-type stars by Hébrard et al. (2016). For this polar model we introduce only 30 spots distributed in time as shown in Table 1. These have 45 deg radii, and are located at colatitudes from 0 to 20 deg, and longitudes normally distributed around 180 deg (σ_n = 90 deg).

• Configuration RAN. A third configuration assumes a random distribution of spots characteristic of very active, fast-rotating, and fully convective stars(>M 4, Barnes et al. 2011). In such cases, spots are found to be rather homogeneously distributed over the full stellar surface with high filling factors. This was shown to be the case through Zeeman Doppler imaging by Morin et al. (2008, 2010). For thisrandom distribution of spots, shown in the third row of Fig. 3, we use 600 spots with the same magnetic cycle as for ONE, but distribute each spot randomly over the surface. To reach the described σ_RV = 5 m s⁻¹, we increase the sizes of spots to 5 ± 1 deg.

• Configuration TWO. In this case, we use the solar example (Zhang et al. 2011) of two active longitudes on opposite stellar sides, which may also be the most common distribution in G-, K-, and early M-type stars (Järvinen et al. 2005; Lanza et al. 2009; Savanov 2014), and is also found for the fully convective late-M dwarf LHS 6351 (Savanov & Dmitrienko 2012). To create the spot maps, we use the same time and colatitude distribution of the 600 spots as described for ONE, but place them on two opposite longitudes (normally distributed around 90 and 270 deg, σ_{rm n} = 45 deg). Two representative maps (low and high spot coverage) are shown in the bottom panels of Fig. 3.

By sampling the simulated generating function at high cadence (1000 min ~ 0.7 d) and using numerical integration, example ACFs of all four configurations are computed, and illustrated in Fig. 4. In the left-hand-side panels, the ACFs are shown using time-lags covering the full time baseline of 2000 days. Although certain long-term variations can be observed, one cannot visually identify a clear feature associated to the injected magnetic cycle of 1000 days, and so the capability of commonly used kernels to determine features from long-term activity cycles is likely to be low. On the right-hand side, a close-in view of the ACFs at small time-lag is shown to highlight the signature of the stellar rotation. The injected spot lifetime of 100 days is seen as a general decay of the correlations at longer time-lags.

For theONE model we see recurring correlations for multiples of P_rot, with slightly decreasing power over time until the injected spot lifetime of 100 days. The spotted side of the star, having a width of 60 deg, smooths out the effect at multiples of P_rot∕2 described by the single-spot models as shown in Fig. 1. This is different from all the other curves. The POL distribution is most similar to the single-spot model, although the peaks at P_rot∕2 are not symmetric, there is diminishing power up to t = 100 days, and there isessentially no power beyond that time because of the phase shift induced by the new location of the dominant surface half. The same behaviour is seen for the peaks at P_rot. For the RAN distribution, the peaks at P_rot∕2 are higher, because the contrast with the less prominent side of the star is smaller, and they are therefore more strongly correlated and precisely located. The effect is even larger for the TWO distribution, where a strong correlation is seen at P_rot∕2. In this case, both peaks diminish in power until t = 100 days, but remain present due to the phase stability of the two opposite spotted regions.

In theory, if the two sides of the surface were identical, the kernel peaks would be at the same correlation level and the amplitude of these peaks would be the same. We try to connect this trend with the differences of filling factors of the prominent side (${\mathcal{F}}_{\text{prom}}$) and its opposite (${\mathcal{F}}_{\text{opp}}$) and calculate the minimum and maximum values in a constant 25-day window over all the spot maps using the relation of $\mathcal{F}$ and κ₀ found in Fig. 2. This delivers the mean minimum and maximum filling factors as an approximation of ${\mathcal{F}}_{\text{opp}}$ and ${\mathcal{F}}_{\text{prom}}$, respectively. The numbers are shown in Table 2. We note the anti-correlation of ${\mathcal{F}}_{\text{max}}={\mathcal{F}}_{\text{prom}}-{\mathcal{F}}_{\text{opp}}$ with the relative strength of the peak of the ACFs at P_rot∕2, as seen in Fig. 4.

Fig. 3 Spot configurations considered in this study. From top to bottom we show spot distributions with one active longitude, ONE (active longitude at 180 deg, σn = 60 deg), with axisymmetric polar spot distribution, POL, with randomly distributed spots, RAN, and with two active longitudes, TWO (90 and 270 deg, σn = 45 deg). The maps depicted illustrate the highest spot coverage at t = 1000 days (left panels)and the lowest coverage at 500 days (right panels). Note, how the spot colatitude is changing for the ONE and TWO models.

Fig. 4 Auto-correlation functions of one example StarSim RV data set for each of the four different spot distributions. Those are, from top to bottom, the ONE, POL, RAN, and TWO configurations. On the left-hand side the ACFs are shown for 0 < t <2000 days and on the right panel for 0 < t < 150 days to highlight the injected stellar rotation at Prot = 25 d, and the spot lifetime of 100 d. Vertical red and blue dashed line mark Prot and Prot∕2 intervals, respectively.

4 Gaussian process regression

We apply GP regression to our RV time-series data created for the various spot configurations as introduced in Sect. 3.2, and use the emcee code to explore the region of the parameter space with largest likelihoods using a Markov chain Monte Carlo procedure (MCMC, Foreman-Mackey et al. 2013). The likelihood optimization code needs to formally assign an uncertainty to each measurement. To deal with this, we assign identical dummy uncertainties of ϵ = 10⁻³ σ_v (σ_v being the variation in RV), and we force the jitter parameter σ (as in Eq. (3)) to be always positive in order to avoid numerical issues (divisions by zero and negative logarithm evaluations). In this sense, the jitter term can be interpreted as the unexplained variance that remains because of an imperfect fit, thus providing an additional figure of merit. In the following sections, we review the two most common kernels used in the literature, and describe the interpretation given to their free adjustable parameters. In Sect. 4.3, we further present an improved version of the QP kernel, namely the quasi-periodic and cosine (QPC) kernel.

4.1 Quasi-periodic kernel

In exoplanetliterature, the most commonly used GP kernel to model astrophysical noise is the quasi-periodic (QP) kernel, given as $$K_{\text{QP}}(\tau )=h^2\exp \left(-\frac{{\tau }^2}{2{\lambda }^2}-\frac{1}{2w^2}{\sin }^2\left(\frac{\pi }{P}\tau \right)\right),$$(15)

which introduces the four hyper-parameters h, P, λ, and w. The kernel is implemented in the george code by Ambikasaran et al. (2015). A graphic representation of this kernel with P = P_rot = 25 days, a range of values for w, and two values of λ (75 and 25 days) are shown in the upper part of Fig. 5.

The range of hyper-parameter w is shown in the figure in light grey covering the commonly used prior between 0 and 1. If we inspect Fig. 5, we can observe the influence of w. Firstly, w reflects the relative strength between the most prominent surface feature and whatever is in the opposite surface half. To match the observed features in the ACF, we examine the simulated ACFs (Fig. 4) in the different configurations. The value of w should be small and of the order of w < 0.5, otherwise the correlations would not approach zero between peaks as they do in the ACFs of the synthetic data. Furthermore, the value of w influences the width of the peaks of the kernel, which is related to how localised the prominent spots are on the dominant surface half, that is, a small localised spot leads to a narrow peak, periodically repeating at each P. If the positions of the most prominent spots change over time or there is an extended group of them, then the width of the peak widens, leading to an approximate lower limit of w >0.2. In general, we can translate the hyperparameter w to the physically realistic behaviour of the synthetic ACFs only in this somewhat arbitrary interval 0.2 < w < 0.5, as shown by the dark grey area in Fig. 5, including thereby w ~0.5, which was suggested by López-Morales et al. (2016).

We maximise the likelihood applying an MCMC procedure with 1000 steps on 1000 walkers. We set further prior limits to the parameters and hyper-parameters as shown in Table 3, where 3.4 days is the constant temporal distance between two observations, 1/4000 days is the Nyquist frequency, $\overline{v}$ is the mean RV of the respective data set, and σ_v its variation. We note in Table 3 the values leading to the largest likelihood values averaged over the five spot maps for each spot configuration. The hyper-parameters are given alongside with RV offset γ, the additional RV jitter σ, and the figures of merit, $\ln \mathcal{L}$ and BIC. In Fig. 6, we compare the best solutions of each kernel with the ACFs of all input RV data sets for the four different spot distributions (from top to bottom) ONE, POL, RAN, and TWO.

We find that the QP kernel is able to identify the introduced 25-day period for every data set. We also find that there is a clear relation of the hyper-parameter λ and the lifetime of spots with T_spots = 100 days = 1.7 ± 0.1λ. Due to the discrete nature of the problem, the number of data-point pairs N_τ with a certain time-lag τ is given with $N_{\tau }(\tau )=\frac{1-N}{T}\tau +N$ (with N as the total number of data points, T as the time baseline), which puts more statistical weight on the correlations at short-time-lags. This is why the GP regression always tries to match the most relevant features of the ACF close to τ = 0 (which tends to follow a Gaussian expτ²∕2λ² envelope), thus explaining why the correct value is found despite the poor match of the kernel and the ACF at large time-lags (see right side of Fig. 6). The hyper-parameter λ is slightly larger for the ONE and TWO configurations because the spots are always at the same rough longitudes, thus producing coherent correlations on longer time-scales. We find that w is in the range from 0.28 to 0.36 for all our experiments with different spot configurations. Moreover, a slight variation of w between those values does not result in significant changes of the likelihoods. As we search for a qualitative translation of the hyperparameters with our test case data, we recommend fixing this parameter to the average w = 0.311 ± 0.016 and removing it as a free parameter of the regression, or at least set priors to [0.2, 0.5], if our aim for the GP modelling is to remove the effect of spots in our RV data. The amplitude h (which shall match the ACF amplitude at zero time-lag κ_0,QP) is recovered quite consistently, which translates to filling factors ${\mathcal{F}}_{\text{max}}~2\%$. The fitted RV offsets are all consistent with zero.

In summary, we find that this kernel and its parameters provide a good option for adjusting correlations and physically interpreting them, but we also find that it fails to correctly recover some relevant information (i.e. the w parameter is not sensitive enough to spot configurations). Using the value of the maximum likelihood and the minimum jitter as figures of merit, we find that this kernel performs at its best on configurations ONE, TWO; and – to a slightly lesser degree – on RAN and POL configurations.

Fig. 6 Auto-correlation functions (ACF, grey and black) and GP kernels (QP blue, SHO green, QPC red) with largest fitted likelihoods of all RV data sets sorted by spot configuration (from top to bottom ONE, POL, RAN, and TWO). The SHO kernel has a less prominent correlation decay and puts more weight on correlations at larger τ in comparison to the QP and QPC kernels. We reiterate the fact that the hyperparameters of the kernels are specifically sensitive tothe correlations at shorter time-lags because of the linear decrease of the number of data-point pairs Nτ with time-lag τ, i.e. $N_{\tau }(\tau )=\frac{1-N}{T}\tau +N$.

4.2 Simple harmonic oscillator kernel

We also study the stochastically driven simple harmonic oscillator (SHO) kernel implemented in the celerite algorithm by Foreman-Mackey et al. (2017) and written as $$\begin{array}{l}{K}_{\text{SHO}}(\tau )=C_0{e}^{-\tau /P_{\text{life}}}\left[\cos \left(\frac{q}{P}\tau \right)+\frac{P}{q{P}_{\text{life}}}\sin \left(\frac{q}{P}\tau \right)\right],\\ q=2\pi \sqrt{{\left(2\pi \frac{P_{\text{life}}}{P}\right)}^2-1},\end{array}$$

with P < 2 π P_life. This kernel introduces one less hyper-parameter than the QP kernel, with C₀, P, and P_life, whereas the parameters that are fitted in the case of the celerite code are the logarithms of $$S_0=\frac{C_0}{2P_{\text{life}}}{\left(\frac{P}{\pi }\right)}^2;w_0=\frac{2\pi }{P};Q=\pi \frac{P_{\text{life}}}{P}.$$(18)

An illustration of the kernel is given in the middle part of Fig. 5 using hyper-parameter P = 25 days, and P_life = 25 (red) and 75 days (black). Here, P results in the sinusoidal variation while P_life is responsible for a slower exponential decay, and hence a stronger correlation for larger time-lags τ compared to the QP kernel. However, the most significant difference is the presence of negative values in this kernel, resulting in anti-correlations for opposite surface halves over time. C₀ is given as $K_{\text{SHO}}(\tau =0)={\kappa }_{0,\text{SHO}}^2$ in m² s⁻².

We apply the same MCMC approach as before and show the results of the GP regression using the SHO kernel in Table 3 and Fig. 6. The filling factors connected to the hyper-parameter C₀ are for all the different spot configurations slightly larger than for the QP kernel. The more consistent solutions for the five spot maps are given with the ONE spot configuration, where a P = P_rot = 25 day period is found, and TWO, where P = P_rot∕2 is visible. In those two cases, where the spots are always distributed around the same longitudes and the P_rot-periodicity in the ACFs is stable over time (black curves), the hyper-parameter P_life is rather large and unconstrained. For the POL and RAN configurations, the SHO kernel and its P hyper-parameter is not able to consistently distinguish between P_rot and P_rot∕2. However, asthe periodicity is lost after a certain time for both configurations, the hyper-parameter P_life is identified consistently with ~100 days. For the POL case there is even a solution with P_life = 8 days. P_life < P is a behaviour which we often find in our studies and in the literature of unevenly sampled RV data. When this happens, the value of P is no longer representative of the real P_rot and a kernelwith only an exponential or a Gaussian profile should be used instead to model correlations.

Large additional jitters σ of up to 3.6 m s⁻¹ are needed to fit the data and, despite having one less hyper-parameter than the QP kernel, BIC values are much larger. Given all these interpretation problems, and the fact that the fits are generally poorer (in terms of both lower likelihoods and higher jitter), we advise against using SHO kernels when studying correlations caused by stellar activity in RV time-series.

4.3 Quasi-periodic and cosine kernel

Given the shape of the ACFs found in our simulations, we propose a new kernel that we implemented to run in the george code (.yml file given in Appendix A), namely a combination of the QP kernel and a damped cosine function, dubbed a quasi-periodic and cosine (QPC) kernel, with, $$\begin{array}{ll}{K}_{\text{QPC}}(\tau )=\hfill & \exp \left(-2\frac{{\tau }^2}{{\lambda }^2}\right)\hfill \\ \hfill & \cdot \left[h_1^2\exp \left(-\frac{1}{2w_0^2}{\sin }^2\left(\pi \frac{\tau }{P}\right)\right)+h_2^2\cos \left(4\pi \frac{\tau }{P}\right)\right].\hfill \end{array}$$(19)

This kernel introduces four hyper-parameters, P, h₁, h₂, and λ. The periodic peaks of the QPC kernel are again given with P, whereas the cosine function adds a periodicity at P∕2. Here, h₁ and h₂ are the amplitudes of the QP component and of the cosine function, respectively, with ${\kappa }_{0,\text{QPC}}^2=h_1^2+h_2^2$. The kernel is shown in Fig. 5 with P = 25 days, h₁ ∕h₂ = 1, and λ = 25 (red), and 75 days (black). This kernel takes advantage of the functional dependence at small time-lags ~ exp(−τ²∕λ²) found in the QP kernel (Eq. (16)), and includes a cosine to account for the feature typically appearing at P/2. We adapted λ_QPC = 2λ_QP so that this hyper-parameter describing the exponential decay of the correlation represents the average spot lifetime more closely. Given its almost null impact in the quality of the fits, we further consider w as a fixed parameter with w = w₀ = 0.31. Nevertheless, we run our models also with free w, as a fifth GP hyper-parameter, and show the statistics in Table 3 (QPCw).

In Table 3 and Fig. 6, we observe that both hyper-parameters P and λ consistently match the introduced P_rot = 25 days and spot lifetime T_spot = 100 days. Offsets are consistent with zero. Additional jitters are small and behave similarly to the values found for the QP kernel. Amplitudes and filling factors are instead the largest of the three kernels, reaching almost 3% spot coverage. Whereas h₁ is relativelyconsistent, we see a very different behaviour for the amplitude h₂, which is directly connected to the contribution of the cosine at τ = P_rot∕2. The relationship between the filling factors of the prominent stellar surface half (${\mathcal{F}}_{\text{prom}}$), which is calculated by both amplitudes with ${(h_1^2+h_2^2)}^{\frac{1}{2}}={\kappa }_{0,\text{QPC}}$, and the opposite half (${\mathcal{F}}_{\text{opp}}$) with h₂ strongly depends on the spot distribution model. It yields ${\mathcal{F}}_{\text{opp}}/{\mathcal{F}}_{\text{prom}}$ = 28, 41, 45, and 64% for ONE,POL,RAN, and TWO configurations, respectively, accounting for the increasing correlation of the peak at P_rot∕2. We find a decrease in BIC compared to the QP kernel with increasing h₂ and significantly better models (dBIC < −10). The BIC values are larger for an additional parameter w (QPCw case), which again indicates that these parameters shall be set to a constant value.

5 Conclusions

In this paper we begin by deriving the basic relations between a generating signal, the covariances it induces, and their relation to Gaussian process kernels. We then calculate the auto-correlation functions of two simple analytic cases to illustrate the physical meaning of the parameters in the resulting kernels. With our StarSim code, we create synthetic RV time-series data for a rotating test-case M-dwarf star with a central single spot, and with four different types of evolving spot maps (which we call configurations) including: one active longitude, two active longitudes, a random spot distribution, and a large polar spot distribution. We use simulated data to explore the relation of the ACFs with stellar inclination and spot filling factor. We then apply GP regression with the commonly used QP and SHO kernels in order to calibrate the physical meaning of the hyper-parameters often presented in the literature. Because of the different imprints of the spot configurations on the ACFs and as a result of this study, we propose the use of a new quasi-periodic and cosine (QPC) kernel. Our main results can be summarised as follows:

• Generating signals of any nature produce covariances in time-series. These can be computed from models and/or observations using the definition of the ACF. This leads to a physical interpretation of the various adjustable kernel hyper-parameters often used in the literature.

• The ACF of single-spot configurations shows its main peaks at multiples of P_rot and is close to zero at multiples of P_rot∕2. It shows further characteristic lobes of negative correlations, which vary with different characteristics of the star and spots, and which are not accounted for in the commonly used GP kernels. Its amplitude κ₀ depends linearly on the filling factor ${\mathcal{F}}_{\text{max}}$ and quadratically on the stellar inclination angle sini.

• The auto-correlation functions of multi-spot configurations explicitly show that correlations generally diminish over acertain timescale connected to the spot lifetime. The repeatability of peaks in the ACF at P_rot and P_rot∕2 is connected to the phase stability of the spot configurations, that is, it remains very strong over very long timescales if the spot longitude remains roughly the same. Peaks in the ACF at P_rot∕2 are always found in ONE,POL,RAN, and TWO spot configurations, where the contrast of these peaks to the correlation at τ = 0 is connected to the filling factor difference between the two halves of the stellar surface.

• Quasi periodic (QP) kernel GP regression accounts for the ACF peaks at P_rot, but not the ones at P_rot∕2. The decay hyper-parameter λ and the lifetime of the simulated spots are clearly mapped to each other. We argue that, despite the fact that the kernel generally does a poor job adjusting the shape of the ACF at large time-lags, it stillidentifies the correct lifetime of the spots because of the larger number of pairs of data points at shorter time-lag contributing more to the likelihood function. As expected, we find that κ₀ from the simulations is mapped one-to-one to the kernel hyper-parameter h. The parameter w behaves as a form factor and should always have values in the range 0.2 < w < 0.5; its introduction as a free parameter does not have a physical translation, nor does it improve the likelihood, and so we recommend fixing it at a reference value of w₀ = 0.31 when modelling spot-induced activity with GPs in RV data.

• Simple harmonic oscillator (SHO) kernel GP regression produces unreliable estimates of the rotation periods, typically obtaining values between P_rot and P_rot∕2, strongly depending on the spot configurations. The amplitude C₀ matches the value of the ACF ${\kappa }_0^2$ at zero time-lag. However, we find that the decay parameter P_life in the SHO varies between large values and the introduced spot lifetime, and therefore its value is not a reliable measure of the typical lifetimes of stellar activity features. Overall, the kernel does not consistently recover the introduced parameters or the shape of the ACF at small time-lags, and produces fits of substantially poorer overall quality. This is reflected in lower maximum-likelihood values and the need for higher jitter values when compared with the other kernels discussed here.

• Our proposed quasi-periodic and cosine (QPC) kernel is designed to use the best features of the QP kernel and adds an additional term to account for the P_rot∕2 peaks in the ACFs. It has the same number of parameters as the QP kernel, and regression experiments find ${\kappa }_{0,\text{QPC}}^2=h_1^2+h_2^2$. The value of h₂ is related to the strength of the ACF peak at P_rot∕2, and and can therefore be used to distinguish between different spot configurations. As for the QP kernel, the periodicity of P_rot and the spot lifetime (λ_QPC = 2λ_QP) are found to be consistent with the corresponding simulated quantities. The better performance of the QPC kernel is demonstrated by the fact that it leads to significantly higher likelihoods (and smaller jitters) for most of our test cases.

We studied the qualitative connection between the auto-correlation functions of different spot configurations and the hyperparameters of different Gaussian process kernels. For this exercise, we used evenly sampled errorless RV data of a spotted rotating edge-on M-dwarf, with no differential rotation and no faculae, and fixed parameters such as the spot lifetime or temperature differences. We did not calculate Bayesian evidence in the form of marginal likelihoods, but relied rather on the qualitative behaviour of the synthetic curves and on the best fit of a Markov chain of Gaussian processes, including the best data likelihoods, additional jitters, and the controversial BIC, averaged over five iterations of stellar spot maps for each spot distribution. If we were to apply the exercise to data with for example realistic measurement uncertainties, uneven sampling, and so on, better statistics could certainly be achieved in some cases by including additional free parameters such as for example the form factor w, but we would lose the physical translation with the evolving dark spots on the rotating stellar surface as obtained by this qualitative study.

We planto further extend our studies from the RVs discussed here to other observable quantities (line widths, chromospheric emission, photometry, etc.), and their corresponding simulated output using StarSim. In addition to this calibration of hyper-parameters, it should be possible to produce improved kernels and to extract more physical information from combined measurements. We also plan to amplify our study to quantify the imprint on the RV data of different stellar parameters such as stellar type, where we expect the same qualitative results. The more closely the analysis of stellar activity can be related to physical processes and quantities, the more able we are to use observational and analysis techniques to mitigate (and eventually “model-out”) the noise of astrophysical origin in high-precision exoplanet searches.

Acknowledgements

The authors acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018-098153-B- C33, as well as the support of the Generalitat de Catalunya/CERCA programme.

Appendix A Quasi-periodic and cosine kernel C++ code as implemented in the george algorithm

name: ExpSine2CosineKernel
stationary: false
params: [per, h1, h2, lam]
reparams:
f1: return -2.0 /lam/lam;
f2: return M_PI/per;
ww: return 0.31*0.31;
value: |
double factor1 = exp(f1 *(x1-x2)*(x1-x2));
double factor2 = h1*h1 * exp(-sin(f2*(x1-x2)) * sin(f2*(x1-x2)) /2.0/ww);
double factor3 = h2*h2 * cos(4.0*f2*(x1-x2));
return factor1*(factor2+factor3);
grad:
per: |
double factor1 = exp(f1 *(x1-x2)*(x1-x2));
double factor2 = h1*h1 * exp(-sin(f2*(x1-x2)) * sin(f2*(x1-x2)) /2.0/ww);
double N1 = sin(f2*(x1-x2));
double N2 = cos(f2*(x1-x2));
double N3 = h2*h2*sin(4.0*f2*(x1-x2));
return factor1*f2*(x1-x2)/per*(factor2/ww*N1*N2 + 4.0*N3);
h1: |
double factor1 = exp(f1 *(x1-x2)*(x1-x2));
double factor2 = h1*h1 * exp(-sin(f2*(x1-x2)) * sin(f2*(x1-x2)) /2.0/ww);
return 2.0/h1*factor1*factor2;
h2: |
double factor1 = exp(f1 *(x1-x2)*(x1-x2));
double factor3 = h2*h2 * cos(4.0*f2*(x1-x2));
return 2.0/h2*factor1*factor3;
lam: |
double factor1 = exp(f1 *(x1-x2)*(x1-x2));
double factor2 = h1*h1 * exp(-sin(f2*(x1-x2)) * sin(f2*(x1-x2)) /2.0/ww);
double factor3 = h2*h2 * cos(4.0*f2*(x1-x2));
return -2.0*f1/lam * (x1-x2)*(x1-x2)*factor1*(factor2+factor3);
x1: |
double factor1 = exp(f1 *(x1-x2)*(x1-x2));
double factor2 = h1*h1 * exp(-sin(f2*(x1-x2)) * sin(f2*(x1-x2)) /2.0/ww);
double factor3 = h2*h2 * cos(4.0*f2*(x1-x2));
double N1 = sin(f2*(x1-x2));
double N2 = cos(f2*(x1-x2));
double N3 = h2*h2*sin(4.0*f2*(x1-x2));
return -factor1* (2.0*(x1-x2)*f1*(factor2+factor3) + f2*h1*h1*x1/ww*N1*N2- 4.0*x1*N3);
x2: |
double factor1 = exp(f1 *(x1-x2)*(x1-x2));
double factor2 = h1*h1 * exp(-sin(f2*(x1-x2)) * sin(f2*(x1-x2)) /2.0/ww);
double factor3 = h2*h2 * cos(4.0*f2*(x1-x2));
double N1 = sin(f2*(x1-x2));
double N2 = cos(f2*(x1-x2));
double N3 = h2*h2*sin(4.0*f2*(x1-x2));
return factor1* (2.0*(x1-x2)*f1*(factor2+factor3) + f2*h1*h1*x1/ww*N1*N2- 4.0*x1*N3);

References

All Tables

All Figures

Fig. 1 Auto-correlation functions of 144 RV data points of a rotating star (Prot = 25 days) with one single central spot. We show the dependence of KACF on stellar inclination angles i = 30, 45, 60, and 90 deg with the colours as indicated. The dependence on the radius of the dark spot R = 5, 10, 20, and 40 deg is illustrated by the different symbols as indicated. We show the repetitive part of the ACFs with τ < Prot.

In the text

Fig. 2 Square root of the maximum ACF value of our test-case M dwarf at τ = 0, κ0, as a functionof the maximum spot filling factor, ${\mathcal{F}}_{\text{max}}$ (top panel), and of the sine of the stellar inclination angle, sini (bottom panel). The curves are shown with colours as indicated for the different inclinations and different spot radii. The black dashed lines are the best fits on the curves.

In the text

Fig. 5 Graphic representation of the three kernels used in this study: the quasi-periodic (QP) kernel as implemented in the george code (top), the simple harmonic oscillator (SHO) as implemented in the celerite code (middle panel), and the quasi-periodic and cosine (QPC, bottom) kernel introduced in this work, for P = 25 days. The x-axis is the time-lag τ between two data points, and the vertical axis is the normalised value of the correlation represented by the kernel. We show the kernels in red for λ = 25 days (w = 0.31) and black forλ = 75 days. For the QP kernel, the light grey area marks 0 < w < 1, whereas the dark grey area marks the interval 0.2 < w < 0.5, which is more physically motivated. The QPC kernel has h2∕h1 = 1.

In the text

Fig. 3 Spot configurations considered in this study. From top to bottom we show spot distributions with one active longitude, ONE (active longitude at 180 deg, σn = 60 deg), with axisymmetric polar spot distribution, POL, with randomly distributed spots, RAN, and with two active longitudes, TWO (90 and 270 deg, σn = 45 deg). The maps depicted illustrate the highest spot coverage at t = 1000 days (left panels)and the lowest coverage at 500 days (right panels). Note, how the spot colatitude is changing for the ONE and TWO models.

In the text

Fig. 4 Auto-correlation functions of one example StarSim RV data set for each of the four different spot distributions. Those are, from top to bottom, the ONE, POL, RAN, and TWO configurations. On the left-hand side the ACFs are shown for 0 < t <2000 days and on the right panel for 0 < t < 150 days to highlight the injected stellar rotation at Prot = 25 d, and the spot lifetime of 100 d. Vertical red and blue dashed line mark Prot and Prot∕2 intervals, respectively.

In the text

Fig. 6 Auto-correlation functions (ACF, grey and black) and GP kernels (QP blue, SHO green, QPC red) with largest fitted likelihoods of all RV data sets sorted by spot configuration (from top to bottom ONE, POL, RAN, and TWO). The SHO kernel has a less prominent correlation decay and puts more weight on correlations at larger τ in comparison to the QP and QPC kernels. We reiterate the fact that the hyperparameters of the kernels are specifically sensitive tothe correlations at shorter time-lags because of the linear decrease of the number of data-point pairs Nτ with time-lag τ, i.e. $N_{\tau }(\tau )=\frac{1-N}{T}\tau +N$.

In the text