© The Authors 2025

1 Introduction

Stellar surfaces exhibit a variety of physical processes occurring on different timescales, including acoustic oscillations, granulation, super-granulation, and magnetic activity (Rincon & Rieutord 2018; Cegla 2019; Meunier 2023). These phenomena have signatures in photometric and spectroscopic data, which can be considered either as signal or noise depending on the scientific objective. In the context of Doppler imaging, the magnetic activity signature is the signal of interest. As the star rotates, the inhomogeneities in the stellar surface create spectral line shape changes. From high-resolution spectra taken at different stellar rotational phases, the goal is to infer a specific intensity map of the stellar surface (Deutsch 1958; Khokhlova 1976; Goncharskii et al. 1977, Goncharskii et al. 1982; Vogt & Penrod 1983; Vogt et al. 1987; Donati et al. 1997; Luger et al. 2021a). In contrast, the complex variability induced by stellar processes poses a significant challenge to exoplanet detection and characterization (Mayor et al. 2014; Sulis et al. 2020; Crass et al. 2021; Hara & Ford 2023). The goal of this series of articles is to provide a detailed, physics-based description of stellar variability signals seen as stochastic processes, useful both to model them as a noise and to gain insight on the star. In this study, we focus on stellar magnetic activity, hereafter referred to as stellar activity.

Stellar activity refers to the ensemble of processes occurring on and above a star’s surface, driven primarily by magnetic fields in the outer stellar layers. In particular, transient features might appear on the surface of stars: starspots (regions of reduced surface temperature caused by strong magnetic fields), or bright regions called faculae (see Fig. 1, and Schrijver 2002; Strassmeier 2009). Spots and faculae create an imbalance of flux on the star. Furthermore, they locally inhibit the convective motion of the gas on the stellar surface (Dravins et al. 1981). Overall, they create complex patterns in photometric and spectroscopic data (Saar & Donahue 1997; Desort et al. 2007; Meunier et al. 2010b; Boisse et al. 2012; Dumusque et al. 2014; Borgniet et al. 2015; Meunier et al. 2019; Meunier & Lagrange 2019a,b). Be they considered as signal or noise, activity-induced signals must be modeled precisely. This can be done by choosing a probability distribution of the data knowing the parameters of the signal of interest, or likelihood function.

In the context of Doppler imaging, the data consist in a time series of spectra taken over one or several rotation periods, possibly with contemporaneous photometry. The likelihood parameters describe the specific intensity map of the stellar surface, allowing one to recover the position of dark and bright magnetic regions. The parameters can be chosen as intensities on a grid covering the visible hemisphere Khokhlova (1976); Piskunov & Rice (1993) or coefficients of spherical harmonics (Luger et al. 2021a). Regardless of the parametrization, the stellar surface properties are usually retrieved by searching the parameters that fit the data best, or equivalently, that maximize the likelihood. Since several intensity maps can fit the data, it is common practice to include a regularization scheme (e.g., Goncharskii et al. 1977, Goncharskii et al. 1982; Donati et al. 1997; Petit et al. 2015). This is equivalent to choosing a prior, and maximizing the posterior distribution of the parameters describing the stellar surface. Alternatively, this posterior probability can be computed directly, which has the advantage of giving uncertainties on the reconstructed map Luger et al. (2021b,a).

In the context of exoplanets, the data are spectroscopic or photometric time series, possibly both. Most stars observed do not rotate sufficiently fast, nor do they exhibit sufficient levels of activity to resolve the stellar surface at a given time. Likelihood parameters related to stellar activity include stellar rotation period and the average lifetime of magnetic regions. The likelihood is almost always assumed to be Gaussian (Aigrain et al. 2012; Haywood et al. 2014; Rajpaul et al. 2015; Foreman-Mackey et al. 2017; Gilbertson et al. 2020; Jones et al. 2022; Hara & Ford 2023; Tran et al. 2023). Its exact parametric form has an impact on the inferred properties of the star and of the putative planets (e.g., Ahrer et al. 2021; Suárez Mascareño et al. 2021). In existing works, the likelihood parameters are not clearly related to physical quantities.

Our aim is to build a likelihood of photometric and/or spectroscopic time series given the average properties of stellar activity – the exoplanet case – which is rooted in a physical model. Conceptually, our solution consists of establishing an analytical link between the likelihood and prior used in Doppler imaging, and the one used in exoplanets. Our framework serves two purposes. First, one can expect that physics-based models disentangle more efficiently the signals due to stellar activity and planets. Second, a physics-based statistical model offers a new way to perform Doppler imaging. Indeed, Doppler imaging is possible only if several conditions are met (e.g., Rice & Strassmeier 2000; Kochukhov 2016): stellar regions must have a sufficient area and contrast to be resolved, the star needs to rotate sufficiently fast in order to see spots move across spectral lines, and the spots must be sufficiently stable during the observation time. However, one can envision that even if one of the conditions is not met and the regions cannot be resolved individually, it might be possible to perform “statistical Doppler imaging”; that is, to retrieve average properties of the stellar surface (stellar inclination, average latitude of spots, etc.).

The article is organized as follows. In Section 2, we present our general statistical formalism to transfer a physical model of the effect of magnetic regions to a likelihood function. We specify the model in Section 3 in two ways. First, we give a fully physics-based model for radial velocity (RV) and photometric observations. Second, we build a general model class that still retains the physical motivation but is more agnostic to the exact assumptions, and that can be applied to any combination of observables of a star (photometry, RV, spectroscopy, etc.). In Section 4, we describe the practical computation of the kernels. In Section 5, we apply our formalism to the analysis of HARPS-N and SORCE observations of the Sun. We compare the performances of the new and existing kernels, and we discuss the ability of physics-based FENRIR kernels to retrieve stellar inclinations as well as average properties of the magnetic regions. In Section 6, we discuss the limitations of our model. In particular, we highlight the presence of non-Gaussianity in RV observations of the Sun. We conclude in Section 7. We provide a public implementation of the FENRIR GPs with tutorial notebooks¹.

2 Statistical framework

2.1 Gaussian likelihoods: Existing approaches

We shall first briefly present how likelihoods are usually defined for exoplanet searches, with the Gaussian process (GP) framework. In general, the data consist of photometric observations, or RVs and spectral indicators, possibly both (see e.g., Hara & Ford 2023). For the sake of clarity, we first consider the case where there is only one type of observations. For example, if we have photometric observations, Y = (y(t_n))_n=1..N where y(t_n) is the photometric observation at time t_n. The likelihood is chosen as a Gaussian distribution, which is fully characterized by its mean f (η) and covariance matrix V(η), and has the form $$p(Y\mid \eta )=\frac{\exp \left[-\frac{1}{2}{(Y-f(\eta ))}^{T}V{(\eta )}^{-1}(Y-f(\eta ))\right]}{\sqrt{{(2\pi )}^N|V(\eta )|}}.$$(1)

$\underset{T}{|}V\left(\eta \right)|$ denotes the determinant of matrix V(η) and the subscript denotes matrix transposition. To simplify the notations but without loss of generality, we assume that η are the parameters of stellar activity. One can easily change Eq. (1) to account for planets, offsets etc. by adding another deterministic term in the mean (see Hara & Ford 2023), f (η) must be replaced by f (η) + h(θ) where h(θ) represents the expected signal of the planets, offsets, trends, depending on parameters θ. To add the observational uncertainties, denoting by Σ the diagonal matrix whose entries are ${\sigma }_i^2$, where σ_i is the error on observation i. The covariance V(η) must be replaced by V(η) + ς.

The value of V at row m and column n quantifies the similarity of the noise values at times t_m and t_n . Typically, the stellar noise is assumed to be stationary and is defined with a kernel function, $$V_{t_m,t_n}=k\left(\left|t_m-t_n\right|,\eta \right).$$(2)

For instance, the so-called quasi-periodic kernel (Aigrain et al. 2012; Haywood et al. 2014) is defined by a vector η of four parameters (η₁, η₂, η_з, η₄), $$k^{\text{QP}}\left(\left|t_m-t_n\right|,\eta \right)={\eta }_1^2\exp \left[-\frac{{\left(t_m-t_n\right)}^2}{2{\eta }_2^2}-\frac{2{\sin }^2\frac{\pi \left(t_m-t_n\right)}{{\eta }_3}}{{\eta }_4^2}\right].$$(3)

This kernel qualitatively expresses that the surface of the star remains similar to itself, though not identical, after one stellar rotation. The parameters η₂ and η₃ can be interpreted as the average lifetime of magnetic regions and the stellar rotation period, respectively. However, none of the parameters are clearly related to a physical model.

When several time series are available, it has been shown to be much better practice to analyze them jointly than separately (Rajpaul et al. 2015; Eastman et al. 2019; Hara & Ford 2023; Yu et al. 2024). We call these different observation types “channels”, and we suppose that we have q channels, not necessarily sampled at the same epochs. For instance, we might have observations through q = 2 channels: photometry and RVs. Channel i has N_i observations and is sampled at times $t_n^i$,n = 1..N_i, and we call the vector of all observations Y = ${\left(y_i\left(t_n^i\right)\right)}_{i=\mathrm{1..}q,n=\mathrm{1..}{N}_i}$. The likelihood of data Y can still be described by Eq. (1). In existing works, the mean is always assumed to be zero. The covariance matrix can be seen as a block matrix, each i, j block gives the covariance between channels i and j. It is characterized by the kernels $$V_{t_m^i,t_n^j}^{i,j}=k_{i,j}\left(\left|t_m^i-t_n^j\right|,\eta \right).$$(4)

Multichannel Gaussian models are often chosen such that each channel is a linear combination of a GP and its derivatives (Rajpaul et al. 2015; Gilbertson et al. 2020; Jones et al. 2022; Delisle et al. 2022). This comes down to choosing a certain parametric form for Eq. (4) (see Gilbertson et al. 2020, Eqs. (8) and (9)). Eq. (4) can also be used when the channels are fluxes at different wavelengths, thus encompasses transit spectroscopy (Gibson et al. 2012; Gordon et al. 2020; Fortune et al. 2024), as well as the statistical Doppler imaging problem introduced above. In summary, both in the single and multi channel case, the data can always be represented as a single vector, and its Gaussian likelihood is fully characterized by mean and covariance functions.

2.2 FENRIR kernels

At the end of our derivations in Section 4, we establish that for any data, Y, consisting of several time series, FENRIR GPs have a kernel associated with channels i and j is of the form $$k_{i,j}(\tau ;\eta )=k^W(\tau ;\eta )k_{i,j}^{\text{per}}(\tau ;\eta )+a_{i,j}(\eta )k^{\lambda }(\tau ;\eta ),$$(5)

where k^W(τ; η) is the part of the kernel corresponding to the intrinsic evolution of spot and faculae properties, and $k_{i,j}^{\text{per}}(\tau ;\eta )$ is the periodic part of the kernel, corresponding to the changing viewing geometry as the star rotates. The kernel k^λ(τ; η) accounts for the magnetic cycle, and its amplitude changes depending on the two channels i, j through a multiplicative coefficient, a_i,j(η).

Eq. (5) has a similar structure to Eq. (3): a decaying function times a periodic one, but in FENRIR models the exact parametric form of the kernel and the hyperparameters are linked to a physical model.

2.3 Finite energy nonlinear impulse response (FENRIR)

Our aim is to describe the statistical distribution of multichannel data, Y, knowing the stellar parameters η, p(Y | η). The first step of our reasoning is to build a physically motivated random process, without assuming it Gaussian. Below, we describe how to generate a sample from p(Y | η).

The first step consists of generating a collection of times of arrival of spots and faculae, (t_k)_{k= −∞..+∞}. The appearance rate and the properties of spots and faculae vary with the magnetic cycles of stars on the timescales of several years (~ 11 years on the Sun) (Martinez Pillet et al. 1993; Borgniet et al. 2015). We refer to transient magnetic regions (spots, faculae, or groups of those) as stellar “features”. Features appear at times drawn from a Poisson process, just like the times of arrival of photons counted on a detector. Poisson processes are characterized by the average number of events per time unit, which might vary with time, λ(t).

When a stellar feature appears at time t_k, its properties γ(t_k) are drawn from a distribution p(γ | η, t_k). This simply means that the properties of the feature depend both on the properties of the star, η, and the time at which it appears t_k . When a stellar feature appears at time t_k and has parameters γ(t_k) it is assumed to have an effect g_i(t − t_k, γ(t_k)) on channel i. For the observation geometry represented in Fig. 2, in Fig. 3 we show the effect of a point-wise dark spot with constant latitude, size and contrast on RVs and photometry.

We use the notation X ~ p(x) to indicate that the random variable X is drawn from the probability distribution p(x). Our model can be simply expressed with three equations. $$t_k\sim \text{Poisson}(\lambda (t)),$$(6) $$\gamma \left(t_k\right)\sim p\left(\gamma \mid t_k,\eta \right),$$(7) $$y_i(t)={\displaystyle \sum _{k=-\infty }^{+\infty }{g}_i}\left(t-t_k,\gamma \left(t_k\right)\right),$$(8)

where

• λ(t) is called the Poisson rate.

• p(γ | t, η) is called the distribution of feature parameters.

• ɡ_i(t, γ) is called the impulse response on channel i.

FENRIR models are filtered Poisson processes (e.g., Lefebvre & Bensalma 2015), where we make the additional assumption that the impulse response has always a finite energy.

2.4 Likelihood of FENRIR models

Regardless of the exact choice of parametrization of our model, given data Y we could try to determine exactly which features affect the dataset. We now need to search their number K, find their times of appearance (t_k)_{k =1,..,K}, and estimate their parameters (γ(t_k))_k=1,..,K. We could search feature parameters with maximum posterior probability, or equivalently minimum – log posterior probability. Using the notation (t_k) = (t_k)_{k =1,..,K} for simplicity, we would solve $$\underset{\left(t_k\right),\left({\gamma }_k\right),\eta }{\mathrm{argmin}}-\log p\left(Y\mid \left(t_k\right),\left({\gamma }_k\right),\eta \right)-\log p\left(\left(t_k\right),\left({\gamma }_k\right)\mid \eta \right).$$(9)

This is similar to the approach taken in Dumusque (2014), which fits a spot model generated with SOAP 2.0 simulations (Dumusque et al. 2014) onto RVs and spectral indicators.

In Eq. (9), the first term and second term in the sum have a counterpart in the Doppler imaging literature. Assuming a Gaussian likelihood, they correspond respectively to the χ² and the regularization term (maximum entropy or Tikhonov, e.g., Goncharskii et al. 1977, 1982; Petit et al. 2015). Usually, in Doppler imaging the fitted parameters do not describe the individual active region, but the local line depth on the stellar surface (Piskunov & Rice 1993; Luger et al. 2021a). Still, the goal is to retrieve the state of the stellar surface with a certain regularization scheme.

However, for slow rotating stars, often targeted to search for exoplanets (Lovis & Fischer 2010), it is impractical to fit the properties of individual magnetic regions. The approach taken here is to compute p(Y | η). We need to marginalize the distribution p(Y | (t_k), (γ_k), η) over (t_k), (γ_k). We found that the marginalized likelihood is impractical to compute analytically. However, we can approximate it by a Gaussian likelihood.

To obtain a Gaussian likelihood in the form of Eq. (1), we computed the mean and covariance of the FENRIR process presented in Section 2.3. The calculation is presented in Appendix A.1. The mean of the FENRIR process for channel i evaluated at time $t_m^i$ is given by $$f_{t_m^i}^i(\eta )={\displaystyle \iint {𝑔}_i}\left(t_m^i-t,\gamma \right)\lambda (t)p(\gamma \mid t,\eta )\text{d}t\text{\hspace{0.17em}d}\gamma ,$$(10)

and the covariance between the channels i and j evaluated at times $t_m^i\text{\hspace{0.17em}and\hspace{0.17em}}{t}_n^i$, respectively, is $$V_{t_m,{,}_n^j}^{i,j}(\eta )={\displaystyle \iint {𝑔}_i}\left(t_m^i-t,\gamma \right){𝑔}_j\left(t_n^j-t,\gamma \right)\lambda (t)p(\gamma \mid t,\eta )\text{d}t\text{\hspace{0.17em}d}\gamma .$$(11)

Eqs. (10), and (11) define unambiguously a Gaussian likelihood.

The Gaussian representation (1) is only an approximation of the likelihood of a FENRIR process. In particular, the higher- order cumulants of the FENRIR processes are usually non-zero. This aspect is further discussed in Section 6.2.

Eq. (11) is not yet identical to the FENRIR kernel in Eq. (5). To obtain this latter expression, as well as guaranteeing a O(N) cost of likelihood evaluation, we need further hypotheses on the impulse responses ɡ_j(t, γ), the statistical distribution of γ, p(γ | t, η), and the Poisson rate λ(t), which we now turn to.

3 Physical and data-driven models

FENRIR likelihoods come in two flavors: physics-based and data-driven. In the first case, we start from quantitative physical arguments to determine the impulse responses, distribution of feature parameters and Poisson rate (see Sections 3.1–3.3). In the data-driven case, we search those in a wide class of functions, not explicitly tied to physical parameters (Section 3.4).

3.1 Impulse response

Spots and faculae affect the photometry and spectroscopy in several ways. While spots are darker than their surroundings, faculae are brighter. Consequently, as they pass accross the visible hemisphere of the star, they have an a photometric effect. Secondly, because they break the flux imbalance of the approaching and receding limb of the star, they have a Doppler signature. Finally, their magnetic field inhibits the convection motion of the plasma in the stellar photosphere. The hot plasma has an upward motion, it cools down and moves back towards the center of the star. Because the plasma moving outwards is hotter, it represents a higher fraction of the total flux resulting in a global blueshift of the stellar light. In spots and faculae, the magnetic field tends to globally slow this motion, resulting in the so-called inhibition of the convective blueshift (Dravins et al. 1981; Meunier et al. 2010a; Dumusque et al. 2014).

To derive analytical expressions of the impulse responses, we made simplifying assumptions, whose validity is discussed in Section 6.1. First, we assume that the measured RV are the sum of the local RVs of the star weighted by their flux. Second, we assume that active regions are point-wise. We denote with ϕ and δ the latitude and longitude of the magnetic region on the stellar surface. We assume that the stellar rotation axis is well defined, and denote with $\overline{i}$ its angle respective to the sky plane (see Fig. 2). Our formalism can account for differential rotation. However, for the sake of simplicity, we assume a solid body rotation, $$\varphi (t)=\omega t+{\varphi }_0.$$(12)

In Appendix B, we show that with our assumptions, the photometric and inhibition of convective blueshift effects in RV, ɡ_ph and ɡ_cb , and the effect on photometry h are $${𝑔}_{ph}(t)={𝟙}_{\text{vis}}(t)W^{\rho }(t)l_{ph}^a(\mu )\mathrm{\Delta }F\omega R^{\star }\mu \cos \overline{i}\sin \varphi \cos \delta ,$$(13) $${𝑔}_{cb}(t)={𝟙}_{\text{vis}}(t)W^{\rho }(t)l_{cb}^a(\mu )\left(\text{Δ}F{V}_{cb}+\text{Δ}{V}_{cb}(F+\text{Δ}F)\right){\mu }^2,$$(14) $$h(t)={𝟙}_{\text{vis}}(t)W^{\rho }(t)l_f^a(\mu )\text{Δ}F\mu ,$$(15)

where R^⋆ is the stellar radius, 1_vis(t) equals one when the feature is in the visible hemisphere, and 0 otherwise, ΔF is the flux difference of the magnetic region with the quiet surface flux, and Δ V_cb is the difference between the mean velocity of the flow of the quiet surface and of the magnetic region. The function W(t) models the variation in the amplitude of the magnetic region as a function of time and is parametrized by a time scale of the appearance of spots ρ. By convention, we consider that W(t) attains its maximum (maximal area of the feature) at t = 0 for this impulse response. The quantity µ is the ratio between the surface of the magnetic region projected onto the sky and its intrinsic surface, $$\mu (\overline{i},\delta ,\varphi )=\sin \overline{i}\sin \delta +\cos \overline{i}\cos \delta \cos \varphi .$$(16)

We have µ = cos ψ, where ψ is the angle between the line of sight and the normal to the magnetic region. As such, it is the quantity classically used to express limb-darkening laws. To simplify the notation, we use µ instead of $\mu (\overline{i},\delta ,\varphi )$.

The central part of the Sun appears brighter than its edges, a phenomenon known as limb-darkening, also present in other stars. Furthermore, faculae and their counterpart in the chromosphere, plages, have a limb-brightening effect which counteracts the limb-darkening behavior of the flux (Frazier 1971; Unruh et al. 1999; Meunier et al. 2010a). The functions l_ph, l_cb, and l_f in Eqs. (13)–(15) are limb-darkening/brightening laws, which are functions of µ, and which we assume to be polynomials of the form: $$l^a(\mu )={\displaystyle \sum _{k=0}^d{a}_k}{\mu }^k.$$(17)

In Fig. 3, we show ɡ_ph, ɡ_cb and h as a function of time for different spot latitudes and stellar inclination, assuming solidbody rotation and constant limb darkening. In Fig. 4, we illustrate the effect on RV of different limb-darkening laws. The limbdarkening effect tends to smooth the effect of a spot, and shift the maximal effect towards the longitude of the observer. The term W^ρ(t) captures the increase and decrease in intensity as spots and faculae grow and vanish (see Fig. 5).

We assume that the effect of a given magnetic region in RV is a weighted sum of the photometric and convective blueshift inhibition effect, and that both effects obey a common limbdarkening law. Due to the degeneracies between the amplitude coefficients, we adopt the following parametrization for the effect of a spot or facula, $${𝑔}_1(t,\gamma )={𝟙}_{\text{vis}}(t)W^{\rho }(t)l^a(\mu )\left({\sigma }_{\text{rv-cb}}{\mu }^2+{\sigma }_{\text{rv-phot}}\mu \cos \overline{i}\sin \varphi \right),$$(18) $$h_1(t,\gamma )={𝟙}_{\text{vis}}(t)W^{\rho }(t)l^a(\mu ){\sigma }_{\text{phot}}\mu ,$$(19)

with γ = (δ, ϕ₀, ρ, σ_phot, σ_rv-cb, σ_rv-phot, ω, $\overline{i}$, a) and where ϕ(t) is given in Eq. (12), and σ_rv-cb, σ_rv-phot and σ_phot are free parameters controlling the amplitude of the respective contributions of the RV inhibition of convective blueshift, RV photometric effect, and flux effect in the photometric time series.

In the FENRIR model, it is supposed that stellar features appear potentially with a variable rate, but independently of each other. However, magnetic regions might appear in certain configurations. On the Sun, spots typically appear on longitudes shifted by 180° (Korhonen et al. 2001; Berdyugina & Usoskin 2003; Borgniet et al. 2015), and there is evidence that this behavior has effects in the Sun HARPS-N RV measurements (Hara et al. 2022). To take this into account, we consider the impulse response g as the effect of two groups of spots shifted by 180°. Overall, using the notations of Eqs. (18) and (19), the expression for the impulse response on RV and flux are $$𝑔(t,\gamma )={𝑔}_1(t)+{\alpha }_{180}{𝑔}_1\left(t-\frac{\pi }{\omega }\right),$$(20) $$h(t,\gamma )=h_1(t)+{\alpha }_{180}{h}_1\left(t-\frac{\pi }{\omega }\right),$$(21)

with γ = (δ, ϕ₀, ρ, σ_phot, σ_rv-cb, σ_rv-phot, α₁₈₀, ω, $\overline{i}$, a). Eqs. (20) and (21) can be used to model the effect of spots or faculae, or their combined effect. We can further extend this idea to groups of four spots or faculae, two symmetric about the equator and two shifted by 180° (see Section 5).

3.2 Distribution of spots and faculae parameters

In the FENRIR formalism of, each time a feature appears, its parameters γ must be drawn from the distribution of features (Eq. (7)). The parameters of the individual features γ are δ the latitude, ϕ₀ the longitude at which the feature reaches its maximal area, ρ the lifetime, σ_phot , σ_rv-cb , σ_rv-phot the amplitudes in photometry and in RV, α₁₈₀ the amplitude of the opposite magnetic region, ω the angular velocity (or the rotation period P), $\overline{i}$ the stellar inclination, a the limb-darkening coefficients. The parameters ρ, σ_phot, σ_rv-cb, σ_rv-phot, α₁₈₀, ω, $\overline{i}$, a are supposed to be identical for all magnetic regions and are included in the population parameters η. The only parameters drawn randomly from the distribution of feature parameters are (δ, ϕ₀).

In the present work, we assume that the longitude at which the feature reaches its maximal area, ϕ₀ , is uniformly distributed on [0,2π]. The latitude is drawn from two Gaussian distributions, symmetrical about the stellar equator, of means ±µ_δ and standard deviation σ_δ. Overall, we have population parameters η = (µ_δ, σ_δ, ρ, σ_phot, σ_rv-cb, σ_rv-phot, α₁₈₀, ω, $\overline{i}$, a).

The assumption that the parameters σ_rv-cb , σ_rv-phot, σ_phot and ρ are the same for all magnetic regions might seem simplistic. However this is not a huge loss of generality. As explained in Appendix C, σ_rv-cb , σ_rv-phot , σ_phot and ρ can be interpreted as the average properties of magnetic regions. Finally, we assume that p(γ|η, t) = p(γ|η), the distribution is independent of time. The model still accounts for the higher rate of apparition of magnetic regions at maximum activity.

3.3 Poisson rate

To specify the FENRIR model, we must choose the rate at which features appear on the star (λ(t) in Eq. (6)). This expresses the fact that the number of magnetic region varies with the magnetic cycle. In Fig. 6 we represent a RV time series simulated with a FENRIR process, such that there are ~50 visible spot groups per day, as observed on the Sun. The effect of inhibition of the convective blueshift is always positive, and as a consequence, a higher appearance rate of magnetic regions translates to a systematic RV effect.

In the present work, we assume that the function λ(t) is itself a stochastic process with a known mean function and covariance, similarly to the approach taken in Camacho et al. (2023). As we shall see in Section 4, this allows for great flexibility as well as a computationally efficient evaluation of the likelihood.

3.4 Formal assumptions for the physics-based and data-driven models

In Sections 3.1–3.3, we specified the key ingredients of the FENRIR model (Eqs. (6)–(8)). Their validity is discussed in Section 6.1. Our hypotheses can be summarized by the following formal assumptions.

1. The impulse response on channel number i, ɡ_i(t, γ), is represented as a product of a time window and a periodic part, each with their own variables. Denoting by γ = (γ_w, γ_per, ϕ₀), $${𝑔}_i(t,\gamma )=W\left(t,{\gamma }_W\right)P_i\left(t,{\gamma }_{\text{per}},{\varphi }_0\right).$$(22)

   In the computations, P_i is assumed to be in the form of a Fourier series $$P_i\left(t,{\gamma }_{\text{per}},{\varphi }_0\right)={\displaystyle \sum _{k=-d}^d{i}_k}\left({\gamma }_{\text{per}}\right){\text{e}}^{\text{i}k\left(\omega t+{\varphi }_0\right)},$$(23)

   where i_k(γ_per) are the k-th Fourier coefficient of P_i(t, γ_per). In the case where the signal is real-valued, as is the case here, we have $i_{-k}=\overline{i_k}\text{\hspace{0.17em}where\hspace{0.17em}}\overline{z}$ denotes the complex conjugate of z. For the sake of simplicity, the window function is the same in all the different channels.

   Fig. 5 gives an example of this representation. The black line represents W(t, γ_W) for a certain value of γ_W, the thin and thick colored lines represent P_i(t, γ_per, ϕ₀) and ɡ_i(t, γ), respectively, for three values of γ_per. Figs. 5a and b represent these quantities for two values of ϕ₀.

2. The phase ϕ₀ is distributed uniformly on [0,2π]. That is, the longitude at which the magnetic regions attain their maximal size is random on [0,2π].

3. The parameters of the window and of the periodic parts are supposed to be statistically independent, $$p(\gamma \mid \eta )=p\left({\gamma }_W\mid \eta \right)p\left({\gamma }_{\text{per\hspace{0.17em}}}\mid \eta \right).$$(24)

   In the context of stellar activity, it means that the evolution of the size of a magnetic region is independent of its longitude and latitude.

4. We assume that the properties of the distribution of feature parameters are independent of time, p(γ|η, t) = p(γ|η).

5. We assume that λ(t) is a stationary random process with expectation µ^λ(η), and covariance function k^λ(τ; η).

6. The variation timescale in λ(t) is much greater than 2π/ω (defined in Eq. (23)): the magnetic cycle has a much longer timescale than the stellar rotation period.

In the physics-based setting, we obtained analytical expressions for each of mathematical expressions appearing in assumptions 1–6. The data-driven kernel consists in considering that the real and imaginary parts of the Fourier coefficients in Eq. (23) are the parameters η of our model. As shown in Section 4.3, if we do so, there is no need to specify the distribution of feature parameters.

4 The physics-based and data-driven kernels

4.1 Computations

From assumptions 1–6 in Section 3.4, in Appendix C we show that the covariance between a measurement taken in channel i at t and channel j at t + τ is of the form of Eq. (5), which is reproduced below for convenience, $$k_{i,j}(\tau ;\eta )=k^W(\tau ;\eta )k_{i,j}^{\text{per}}(\tau ;\eta )+a_{i,j}(\eta )k^{\lambda }(\tau ;\eta ),$$

such that $$k^W(\tau ;\eta )={\displaystyle \iint W}\left(t,{\gamma }_W\right)W\left(t+\tau ,{\gamma }_W\right)p\left({\gamma }_W\mid \eta \right)\text{d}t\text{\hspace{0.17em}d}{\gamma }_W,$$(25) $$k_{i,j}^{\text{per}}(\tau ;\eta )={\displaystyle \iint P^i}\left(t,{\gamma }_{\text{per}}\right)P^j\left(t+\tau ,{\gamma }_{\text{per}}\right)p\left({\gamma }_{\text{per}}\mid \eta \right)\text{d}t\text{\hspace{0.17em}d}{\gamma }_{\text{per}},$$(26) $$a_{i,j}(\eta )={\mu }_i(\eta ){\mu }_j(\eta ),$$(27)

where µ_i(η) = E_γ(i₀(γ)), µ_j(η) = E_γ(j₀(γ)), the mathematical expectancy of i₀ and j₀ defined in Eq. (23). The kernel, k^λ(τ; η), is the covariance λ(t). It accounts for the magnetic cycle, and its amplitude changes depending on the two channels i, j through a multiplicative coefficient, a_{i, j} . The fact that this covariance term is additive is a simplification due to assumption 6. If the time-scale of stellar rotation and magnetic cycle were closer, the expression of the FENRIR kernel would be more complex (see Appendix C.6). The periodic part of the kernel can be expressed in a Fourier series, $$k_{i,j}^{\text{per}}(\tau ;\eta )={\displaystyle \sum _{k=-d}^d{C}_k^{i,j}}(\eta ){\text{e}}^{\text{i}k\omega \tau },$$(28)

where $$C_k^{i,j}(\eta )={\displaystyle \int i_{-k}}\left({\gamma }_{\text{per}}\right)j_k\left({\gamma }_{\text{per}}\right)p\left({\gamma }_{\text{per}}\mid \eta \right)\text{d}{\gamma }_{\text{per}},$$(29)

d, i_k and j_k being respectively the number of harmonics, and the Fourier coefficients of the impulse responses in channels i and j (see Eq. (23)).

In Appendix C, we present a method to evaluate the likelihood with a cost ∝ Nr_W(2d + 1) where N is the total number of data points, r_w a coefficient linked to the window (typically one to three), and d is the number of harmonics in the Fourier expansion of the periodic part. This computational scheme applies both to the physics-based and data-driven models, we now highlight their specificities.

4.2 The physics-based kernels

In the case of the physical model, the kernel $k_{i,j}^{\text{per}}$ is fully specified by the choices made in Section 3. However, for a specific choice of the distribution of spots latitudes and longitudes, the integral in Eq. (29) may not have an analytical expression. This can be circumvented in two ways. The first option consists in considering that the impulse response has always the same parameters, and that the average group of magnetic regions that it describes represents the average effect of magnetic regions. The second option is, for a given value of the hyperparameters η, we compute this integral numerically by discretizing the feature parameters γ_per. To improve performances, we precompute the integrals on a grid of values of hyperparameters, so that, when we need to evaluate the likelihood for a given value of η, we simply interpolate the integral values computed on the grid (the exact procedure is slightly more elaborate, see Appendix C).

In Section 3, we did not specify our choice for the evolution of the size of stellar spots and faculae, W(t). In Section 5, we use Matérn 1/2 or Matérn 3/2 kernels corresponding, respectively, to assuming that W(t) is a one-sided exponential or symmetric exponential (see Table C.1). Finally, for the kernel representing the magnetic cycle, k^λ , we used a Matérn kernel.

4.3 The data-driven MultiFourierKernel

The physics-based model is specified only for RV and photometry. Furthermore, the physical assumptions on the impulse response and distribution of feature parameters might be faulty, and we could want to be more agnostic to them. Interestingly, for any choice of spectroscopic indicator and photometric channels, all the possible FENRIR kernels with assumptions 1–6 can be written as in Eq. (5) where the periodic part of the kernel is what we call a MultiFourierKernel, $$k_{i,j}^{\text{per}}(\tau ;\eta )={\displaystyle \sum _{k=0}^d{a}_{k;i,j}}(\eta )\cos (k\omega \tau )+b_{k;i,j}(\eta )\sin (k\omega \tau ),$$(30)

where a_k(η) = (a_k;i,j(η))_{i,j =1..q} and b_k(η) = (b_k;i,j(η))_{i,j =1..q} are q × q matrices defined as $$a_k(\eta )={\alpha }_k{\alpha }_k^T+{\beta }_k{\beta }_k^T,$$(31) $$b_k(\eta )={\alpha }_k{\beta }_k^T-{\beta }_k{\alpha }_k^T,$$(32)

where η = (α_k, β_k)_k=0..d, and matrices α_k, β_k have a form $${\alpha }_k=\left(\begin{array}{cccc}{\alpha }_k^{11}& 0& \dots & 0\\ {\alpha }_k^{21}& {\alpha }_k^{22}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\alpha }_k^{q1}& {\alpha }_k^{q2}& \cdots & {\alpha }_k^{qq}\end{array}\right),$$(33) $${\beta }_k=\left(\begin{array}{ccccc}0& 0& \dots & \dots & 0\\ {\beta }_k^{21}& 0& \dots & \dots & 0\\ {\beta }_k^{31}& {\beta }_k^{32}& 0& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\beta }_k^{q1}& {\beta }_k^{q2}& {\beta }_k^{q3}& \dots & 0\end{array}\right).$$(34)

For real-valued data, they satisfy α_−k = α_k and β_−k = −β_k. In particular, for k = 0, β₀ = 0.

If we were to choose the expression of the Fourier coefficients i_k(γ) and the distribution p(γ | η) as generally as possible, we would not have a class of kernels more general than simply choosing as hyperparameters η the coefficients of matrices α_k and β_k, k = 0,…, d.

Here we have q observational channels. If we choose η = (α_k, β_k)_k=0..d, everything happens as if there were q FENRIR processes indexed by l = 1..q, statistically independent and with the same rate, which have a fixed impulse response of the form Eq. (23), and which we call modes. In the context of stellar activity, this means that everything happens as if there were at most q equivalent spots, always with exactly the same parameters, appearing independently of each other. The matrices α and β define the real and imaginary parts of these “equivalent spots”/modes up to a factor (see Eq. (35)).

For k > 1, α_k and β_k matrices respectively have q(q + 1) and q(q − 1) free coefficients, so in total the model has q² real-valued coefficients per harmonic. For k = 0 there are only q(q + 1)/2 real coefficients. There are d harmonics. In total, we have at most dq² + q(q + 1)/2 free coefficients in the model. However, we can make an assumption on the rank of the matrices α_k and β_k, which reduces the number of free coefficients. In the following sections, we refer to the number of columns of α_k and β_k that are non-zero as the number of modes, and denote it by m. Overall, the MultiFourierKernel is defined by the number of modes m and the maximum harmonic order d (see Eq. (23)).

For instance, if we assume a single mode, only the coefficients of the first column of are non zero. It comes down to assuming that it is always the same equivalent spot that appears randomly. Assuming that d = 3 harmonics are used, everything happens as if there was a single type of spot appearing randomly, of impulse response on channel j, $${𝑔}^j(t)=W(t){\displaystyle \sum _{k=-3}^3{𝑔}_k^j}{\text{e}}^{\text{i}k(\omega t+\varphi )},$$(35)

with ${ɡ}_0^j={\alpha }_0^{j1}$ and ${ɡ}_k^j={\overline{ɡ}}_{-k}^j=\frac{{\alpha }_k^{j1}+\text{i}{\beta }_k^{j1}}{\sqrt{2}}$ for k = 1,2,3.

It might seem surprising that in the data-driven version of FENRIR, the number of equivalent spots depends on the number of channels. This can be understood intuitively in the three- channel case outlined above, by considering that three types of magnetic regions might appear, each affecting only a certain channel. For example, spots of type one are only visible in RV, spots of type two in photometry, etc. To describe this situation, we need three independent templates. This case, although not physically motivated, shows that q templates might be needed to describe q channels. The fact that more than q templates do not introduce more expressivity of the kernels comes from degeneracies of the approximation of FENRIR processes by GPs. We refer the reader to Appendix C.7 for more details.

5 Analysis of HARPS-N solar RVs and SORCE Total Solar Irradiance

To illustrate the use of FENRIR GP kernels, we analyzed the HARPS-N solar RV time series² (Dumusque et al. 2021) together with the SORCE Total solar Irradiance time series³ (Kopp 2020). We jointly modeled the RVs and the photometric time series using a FENRIR kernel.

The SORCE photometry is available binned by 6h intervals or binned by day. In the following experiments, we are mainly interested in modeling the effect of active regions (modulated by the rotation period of the Sun), whose timescales are much longer than a day. We thus used the daily binned SORCE photometry and also binned the HARPS-N time series by day for consistency.

The SORCE time series covers a time span of 17 years from February 2003 to February 2020, while the HARPS-N data cover three years from July 2015 to July 2018. In Section 5.1, we detail implementation choices for the physics-based FENRIR models applied to the Sun. In Section 5.2, we test the ability of these models to perform statistical Doppler imaging, that is to correctly assess the physical parameters related to the Stellar obliquity and the active-region properties. In Section 5.3, we compare the predictive accuracy of FENRIR kernels (physics-based and data-driven) to other kernels with cross-validation.

5.1 Physics-based FENRIR GP for the Sun

The FENRIR model described in previous sections is a versatile framework. Here we provide more details on the specific implementation we used for the analysis of solar data. Since we used daily binned data, we assumed short-term effects, such as granulation and oscillations, to be mostly averaged out, and we modeled them with a jitter term (white noise) on each time series. We also introduced an offset for each time series, which in particular captures the mean of the FENRIR process.

5.1.1 Active regions and magnetic cycle

We tried different ways of modeling the population of spots and faculae appearing on the Sun’s surface by varying different hypotheses:

1. Spots and faculae as separated (hereafter “sep.”) or mixed populations (herefater “mix.”).

   We modeled the spots and faculae either as two independent populations (appearing at independent location and time) or by modeling an active region altogether by considering the average effect of spots and faculae (dis)appearing at the same location and time. Faculae tend to appear more often in the vicinity of spots (which would advocate the latter case), their lifetime is also typically longer than the spots’, and some active regions are also observed with faculae but without any spot (which would favor the former case). In Section 3, we showed the expression of the FENRIR model for spots and faculae. We simply note that if we use a sum of two independent FENRIR GP processes, one modeling spots the other faculae, the covariance of the sum is the sum of covariances. Alternatively, in “mix” model the impulse response is a weighted sum of the spot and the facula impulse responses (see details below).

2. Latitude distribution (hereafter “lat. dist.”) or opposite latitudes (hereafter “opp. lat.”).

   We always assumed the distribution of latitude of the active regions to be symmetrical with respect to the equator. In the first case, the distribution was assumed to follow a mixture of two Gaussian distributions with means ±µ_δ and standard deviations σ_δ. In the second case, we used an approximation of this law, where the active regions appear always at the two same opposite latitudes ±δ (i.e., σ_δ → 0). This approximation might be crude, but is also much more efficient to compute.

3. Matérn 1/2 or 3/2 decay kernel, with timescale ρ for the window function of the periodic part. The Matérn 1/2 kernel might be interpreted as a sudden appearing and exponential decay of spots/faculae while the Matérn 3/2 kernel would correspond to a symmetric exponential appearing and disappearing (see Appendix C.4).

4. Matérn 1/2 or 3/2 decay kernel for the magnetic cycle part.

   By taking all possible combinations of these four choices, we ended up testing 16 different physical models of active regions and magnetic cycle.

We always used a quadratic limb-darkening law (see Eq. (17)) with coefficients a₀ = 0.3, a₁ = 0.93, a₂ = −0.23 (e.g., Livingston 2002). Both spots and faculae are affected by the limb-darkening effect. However, the faculae are additionnaly affected by the limb-brightening effect (e.g., Dumusque et al. 2014), as their contrast with respect to the surrounding region increases when they reach the limb of the Sun. We model this effect in the same fashion as the limb-darkening, using a quadratic law: $$b(\mu )=b_0+b_1\mu +b_2{\mu }^2,$$(36)

with coefficients (see Meunier et al. 2010b): $$\begin{array}{cc}& b_{\text{plage},0}=1-b_{\text{plage},1}-b_{\text{plage},2},\\ & b_{\text{plage},1}=-\frac{407.7}{250.9-407.7+190.9},\\ & b_{\text{plage},2}=\frac{190.9}{250.9-407.7+190.9}.\end{array}$$(37)

With this choice, the brightening is equal to one when the faculae is at the center of the star. The contrast of a spot is on the contrary approximately constant (i.e., b_spot,0 = −1, b_spot,1 = b_spot,2 = 0).

These limb-darkening and limb-brightening laws enable one to model the populations of spots and faculae independently (sep. model). In order to model the overall effect of an active region including spots and faculae (mix. model), we now need to scale their contributions to the active region. While the flux decrease due to a spot is about 20 times stronger than the flux increase due to a plage of the same size at the center of the Sun, the surface covered by plages is also about 20 times larger (e.g., Meunier et al. 2010b). This means that for a typical active region, the scaling of the spot component should be about the same as that of the plage component. We thus adopt the following limb-brightening law for an active region containing spots and faculae: $$\begin{array}{cc}& b_0=-b_1-b_2,\\ & b_1=b_{\text{plage},1},\\ & b_2=b_{\text{plage},2}.\end{array}$$(38)

Since the limb-brightening effect models the change of contrast of the active region with respect to its surrounding, the overall effect of limb-brightening and limb-darkening is the product of the two polynomials l(µ)b(µ) (i.e. a polynomial of order 4), which is then to be combined with the projection effect, etc., as described in Eqs. (13)–(15). We illustrate in Fig. 7 the photometric effect of a spot, a plage, and the global photometric effect of an active region.

As explained in Sect. 4.2, computing the periodic part of the kernel requires us to evaluate the integrals of Eq. (29) numerically. We precomputed them on a grid of values for $\overline{i}$, µ_δ, and σ_δ, so that when we run a Monte Carlo Markov chain (MCMC) estimate of the posterior distribution of η, the integrals of Eq. (29) can be simply evaluated by interpolation. We developed the Fourier series up to the third harmonic of the rotation period (d = 3). We sampled $\overline{i}$ and µ_δ on a regular grid of 51 values in the range [0, π/2], and σ_δ on a regular grid of 51 values in the range [π/36, π/2]. For each point in the grid, we sampled 101 values of δ in the range $[\overline{i}-\pi /2,\pi /2]$ (latitudes that are visible from the observer) to integrate the Fourier coefficients over the population of spots/faculae (see Appendix C.5). We also computed the two coefficients (average photometeric effect and average RV convective-blueshift effect) of the terms due to the variations of the active region appearance rate (i.e., magnetic cycle, see Appendix C.6). Finally, we computed the likelihood for any value of $\overline{i}$,µ_δ, σ_δ by interpolation. In the simpler cases with only two opposite longitudes (opp. lat. model), integrals Eq. (29) were computed analytically, without interpolation on a precomputed grid.

We first included the effect of active regions at the opposite longitude (parameter α₁₈₀, see Eqs. (20) and (21)). We then multiplied the two magnetic cycle amplitudes by a factor η_mag.. This parameter measures the ratio between the amplitude of the magnetic cycle variations and the periodic variations due to stellar rotation. Finally, we normalize all the coefficients such that the total amplitude (considering both the periodic and magnetic cycle components) in photometry is σ_phot., the total amplitude of the RV convective blue-shift effect is σ_rv−cb, and the amplitude of the RV photometric effect is σ_rv−phot.. We note that in our model, there is no contribution of the RV photometric effect in the magnetic cycle since this effect averages out.

5.1.2 Model complexity and computing cost

All the 16 possible models presented above can be implemented as a S+LEAF covariance matrix, The exact rank r of the matrix depends on the model, ranging from r = 15 when using a single mixed population of active regions and Matérn 1/2 decay kernels to r = 90 when using two separated populations for the spots and the faculae and Matérn 3/2 kernels. The cost of the likelihood evaluations of the model scales as O(r²n) (see Foreman-Mackey et al. 2017; Delisle et al. 2022), where n is the total number of measurements (including RV and photometry).

5.2 Statistical Doppler imaging: Proof of concept

In Sect. 5.3.3, we compare the 16 FENRIR physics-based models presented above, together with various data-driven models, in their ability to predict missing data (cross-validation). Here we focus on the two best physics-based models, which both consider spots and faculae as independent populations (sep. model) and use Matérn 1/2 kernels for both the periodic and magnetic cycle components. The best model uses latitude distributions (lat. dist.) of active regions, while the second-best model uses two opposite latitudes (opp. lat.).

We used the samsam MCMC sampler⁴ (e.g., Delisle et al. 2018) to explore the parameter space. We provide the priors and posteriors of all the parameters explored in Tables 1 and 2. In Figs. 8 and 9, we show a corner plot of the stellar inclination $\overline{i}$ and the parameters of the distribution of spots and faculae in latitude. We additionally show the full corner plots, as well as convergence diagnostics in Appendix D (see Section 7). Since we observe the Sun from Earth, its inclination is not constant and oscillates through the year $(\overline{i}\in [-7,7]\text{deg})$. The distribution of active regions is typically symmetric with respect to the equator, with µ_δ ≈ 15° and σ_δ ≈ 8° (e.g., Mandal et al. 2017). The symmetry with respect to the equator implies that having positive or negative $\overline{i}$ values yields the same effects. The average absolute inclination through the year is about 4.5°. In both MCMC runs, the inclination is very well constrained to be close to zero, and compatible with the true average value. The average latitude of spots is also restricted to be close to the equator and compatible with the true value. The same is true for the average latitude of faculae. The width of the latitude distribution σ_δ mostly follows the priors for the spots, which means that it is not constrained by the data. In the case of the faculae, this width seem a bit more constrained towards low values. For both populations, the width is compatible with very small values. It is thus not surprising that the model assuming two opposite latitudes (equivalent to σ_δ → 0) is performing almost equally well and gives similar results.

In Eq. (18), we express the effect of a spot or faculae on RV as a function of its parameters, in particular the amplitude of their photometric and inhibition of convective blueshift effects. Since here we have a sep. model, we have amplitude terms both for the spots and faculae populations. Since spots are darker regions and faculae are brighter regions, one would expect σ_spots,phot. < 0, σ_{spots, rv-phot}. < 0, σ_faculae,phot. > 0,^and σ_{faculae, rv-phot}. > 0. Moreover, both kinds of region inhibit the convective blueshift, which would thus result in a net redshift: σ_{spots, rv-cb} > 0, σ_{faculae, rv-cb} > 0. Since we consider here a GP approximation of the FENRIR process, the overall sign of the effect is degenerate, but the relative signs of the amplitudes σ_spots,phot., σ_spots, _rv-phot., and σ_{spots, rv-cb} are meaningful. The same is true for the relative signs of the amplitudes for faculae. In Tables 1 and 2, we see that all these signs are correctly recovered, except for σ_{faculae, rv-phot.} which is found to be negative instead of positive. However, this value is constrained to be small compared to the inhibition of the convective blue-shift, and compatible with 0. This is expected since the contrast of faculae is low and their effect in RV is dominated by the inhibition of the convective blueshift.

The decay timescale of faculae (73 ± 8 d) is found to be about twice that of spots (40 ± 3 d). This trend is expected, since during the lifetime of an active region, spots typically disappear after one rotation, while faculae disappear after a few (typically 3-4) rotation periods (e.g., Xu et al. 2017). We note that the model tend to favor having spots appearing at opposite longitudes (α_{spots, 180} = 0.13 ± 0.05), while this is not true for faculae $\left({\alpha }_{\text{faculae,\hspace{0.17em}}180}={0.003}_{-0.002}^{+0.005}\right)$. Finally, the magnetic cycle component of the model is dominated by the contribution of faculae (η_{spots, mag.} ≪ η_{faculae, mag.}). As is noted in Appendix C.6, the amplitude of the magnetic cycle term is expected to increase with the typical lifetime of spots/faculae. Since the lifetime of faculae is expected to be longer than the lifetime of spots, we therefore expect their effect to dominate in the magnetic cycle (see e.g., Shapiro et al. 2016).

In Fig. 10, we show the kernel function corresponding to the maximum a posteriori set of parameters for the best physical model. It should be noted that the amplitude and timescale of the magnetic cycle component are degenerated (see the corner plot in Appendix D (see Section 7), which means that the long term component of this kernel (k^λ(τ; η) in Eq. (5)) is not well constrained. The short term part (k^W(τ; η)k^per(τ; η) in Eq. (5)), shown in Fig. 11 is much better constrained. Finally, Figs. 12 and 13 show the GP’s conditional distribution corresponding to the maximum a posteriori set of parameters for the best physical model. The second best model gives very similar results.

This analysis shows that it is possible to perform statistical Doppler imaging, that is, to retrieve the stellar inclination and the average properties of active regions from the statistical distribution of the signal.

5.3 Predictive performance for different kernels

We now turn our attention to the capability of FENRIR kernels to adequately model the data, and whether they perform better than existing models. As a performance metric, we chose the cross-validation score because it is efficient to find a model with a good trade-off between its simplicity and its agreement with the data. Besides the physics-based FENRIR kernels, we also test the following ones.

5.3.1 The FENRIR data-driven model

These kernels are very similar to the physics-based ones, and we still considered a Fourier decomposition up to the third harmonics of the rotation period (P_rot./3). However, in the data-driven case, instead of computing this Fourier decomposition from a set of physical parameters, the coefficients of the matrices α_k, β_k (see Eqs. (33) and (34)) are left free for each harmonic k = 0,…, 3. We performed tests varying the number of modes as defined in Section 4.3. The magnetic cycle is also modeled with free coefficients and always with two modes.

In summary, in this data-driven approach, the two components (active-regions and magnetic cycle) can both be represented as the product of a MultiFourierKernel and a decay kernel. For the magnetic cycle, we only consider the constant part of the Fourier decomposition (up to harmonics zero), while the active-regions part is developed up to the third harmonics. We used two modes for the magnetic cycle and varied the number of modes for the active-regions component.

5.3.2 Latent GP and derivatives

In order to compare the FENRIR models with more classical approaches, we designed a second set of data-driven models where the active region part is modeled using one or more latent GPs G_k (Aigrain et al. 2012; Rajpaul et al. 2015; Gilbertson et al. 2020; Jones et al. 2022; Tran et al. 2023). The active region component in the channel i (RV or photometry) is modeled as a linear combination of these GPs and their derivatives, $$y_{ispot}(t)={\displaystyle \sum }_{k=1}^{n_{GP}}{\alpha }_{i,k}{G}_k(t)+{\beta }_{i,k}{\dot{G}}_k(t).$$(39)

where α_i,k, β_i,k here are dummy variables, not the same as in Eqs. (33) and (34). We assume an ESP kernel for the GPs G_k. The ESP kernel is an approximation of the widely used squared- exponential periodic (SEP) kernel (e.g., Rajpaul et al. 2015): $$k_{SEP}(\tau ;P,\rho ,\eta )=\exp \left(-\frac{{\tau }^2}{2{\rho }^2}-\frac{{\sin }^2\left(\frac{\pi \tau }{P}\right)}{2{\eta }^2}\right).$$(40)

This approximation allows for efficient likelihood evaluations (i.e., with a linear scaling with the number of measurements, see Delisle et al. 2022). For the sake of consistency, we develop this approximation up to the third harmonics of the rotation period (P_rot./3). We let the hyper-parameters of each latent GP (P_k, ρ_k, η_k) free and independent of each other. The magnetic-cycle component is modeled as in the FENRIR data-driven case.

5.3.3 Cross-validation tests

In order to compare the ability of different GP kernels to model the Sun’s RV, we designed cross-validation tests. We divide the data into a training set y_train and a test set y_test. For each kernel, we run an MCMC on the training set and then compute the conditional probability of the test set given the model M: $$p\left(y_{\text{test\hspace{0.17em}}}\mid y_{\text{train\hspace{0.17em}}},M\right)={\displaystyle \int p}\left(y_{\text{test\hspace{0.17em}}}\mid y_{\text{train\hspace{0.17em}}},M,\eta \right)p\left(\eta \mid y_{\text{train\hspace{0.17em}}},M\right)\text{d}\eta ,$$(41)

where η is the set of the model’s hyper-parameters that are explored in the MCMC. The MCMC algorithm allows one to estimate this integral as $$p\left(y_{\text{test\hspace{0.17em}}}\mid y_{\text{train\hspace{0.17em}}},M\right)\approx \frac{1}{N_{\text{samples\hspace{0.17em}}}}{\displaystyle \sum _{k=1}^{N_{\text{samples\hspace{0.17em}}}}p}\left(y_{\text{test\hspace{0.17em}}}\mid y_{\text{train\hspace{0.17em}}},M,{\eta }_k\right);$$(42)

that is, the average of the conditional probability computed for each MCMC sample. Finally, the conditional probability for a given sample can easily be computed using its definition: $$p\left(y_{\text{test\hspace{0.17em}}}\mid y_{\text{train\hspace{0.17em}}},M,{\eta }_k\right)=\frac{p\left(y_{\text{test\hspace{0.17em}}},y_{\text{train\hspace{0.17em}}}\mid ,M,{\eta }_k\right)}{p\left(y_{\text{train\hspace{0.17em}}}\mid M,{\eta }_k\right)}=\frac{p\left(y_{\text{full\hspace{0.17em}}}\mid ,M,{\eta }_k\right)}{p\left(y_{\text{train\hspace{0.17em}}}\mid M,{\eta }_k\right)},$$(43)

which is the ratio of the likelihoods of the parameters η_k computed on the full dataset and on the training set.

This conditional probability allows one to compare the different models in their ability to predict the test set from the training set. Here we test our ability to disentangle planetary signatures from stellar activity in the RVs. We thus chose to keep all the photometric data in the training set, and split the RV timeseries between the training and the test sets. Moreover, we are mostly interested in comparing the active-regions component of the models. We thus decided to take eight chunks of 30 d in the RV timeseries for the test set.

We show in Table 3 the results of these cross-validation tests for the three models considered: the physics-based FEN- RIR model (see Sect. 5.1), the data-driven FENRIR model (see Sect. 5.3.1), and the classical case of latent GPs (see Sect. 5.3.2). For the physics-based FENRIR model, we tried the 16 variations of the model described in Sect. 5.1. For the data-driven FENRIR model, we tried to use one or two modes (for all harmonics) and a Matérn 1/2 or Matérn 3/2 kernel (both for the active-regions and for the magnetic cycle). Finally, for the latent GP case, we also tried one or two modes (i.e., one or two independent GPs), and varied the decay kernel for the magnetic cycle component (Matérn 1/2 or 3/2). The kernel of the latent GPs modeling the active regions is always an ESP kernel.

In Table 3, we provide the cross-validation score computed following Eq. (42) in log-scale for each model. This score is thus the log of the average conditional probability over the MCMC samples. Additionally, we provide the median and the 68.27 % interval of the log conditional probability over the MCMC samples. This interval illustrates the dependency of the conditional probability on the hyperparameters. Finally, we provide the RMS of the residuals on the train and test sets.

The best data-driven FENRIR models show slightly better performance than the best physics-based FENRIR. The difference in cross-validation score is of the same order of magnitude as the scatter in score among each model. In contrast, the best latent GP model shows significantly worse performances.

Among, the data-driven FENRIR models, the best performances are obtained when using two modes, which means that there are two independent component in the activity signal.

Among the physical-based FENRIR models, those modeling spots and faculae as two independent populations (sep.) clearly outperform the mixed populations models (mix.). The reality is in between those two extreme cases; that is, there is some correlation in the appearance of spots and faculae, and in particular some active regions consist of a spot surrounded by faculae, but this is not always the case. Our result suggest that considering spots and faculae as two independent populations is a better approximation. As discussed in Sect. 5.2, in the sep. case, while the latitude distributions of both types of active regions are similar, the decay timescale of spots is about half the decay timescale of faculae. This might explain why the sep. model is preferred over the mix. model, which inherently assumes that spots and faculae in the same active region appear and disappear at the same time. Regarding the latitude distribution (lat. dist. or opp. lat.), the two alternative models show very similar performances. As shown in Sect. 5.2, while the average latitude of active regions can be constrained, the width of of the distribution is not well constrained, which might explain the similar performances of both models.

The Matérn 1/2 kernel seems to be slightly preferred over the Matérn 3/2 kernel, both for the active regions component and for the magnetic cycle component, and both in the physics-based model and in the data-driven model.

5.3.4 Cross-validation tests using additional activity indicators

While in this article we derived the physics-based FENRIR model for RV and photometry only, the approach is more general and could be extended to additional activity indicators. Such an extension of our physics-based model is beyond the scope of this article. However, the data-driven FENRIR model, as well as the classical latent GPs approach, can readily be extended to any additional channels. We thus ran a larger series of crossvalidation tests, using the activity indicators extracted from the HARPS-N solar spectra at the same time as the RVs as provided by Dumusque et al. (2021).

We use the exact same splitting of the HARPS-N RVs between the training and test sets, and use different combinations of photometry and activity indicators as additional training data. Because we compute the conditional probability on the same RV test set, we can directly compare the cross-validation scores for any model and any set of activity indicators.

We tested all the possible combinations of photometry, full width half maximum (FWHM), bissector span (BIS), $\log {\text{R}}_{\text{HK}}^{\prime }$, both with the data-driven FENRIR model and the classical latent GPs model, with one or two modes in both cases and with Matérn 1/2 or 3/2 decays for the active regions and magnetic cycle parts.

The 30 best models, in terms of cross-validation scores, are shown in Table 4, while the full list is shown in Appendix E (see Section 7). We also show in Fig. 14 the GP’s conditional distribution corresponding to the maximum a posteriori set of parameters for the best model (data-driven FENRIR GP with a single mode and trained on the FWHM and BIS indicators).

Here are a few comments inspired by these results:

• The best cross-validation scores (Table 4) are significantly higher than the scores obtained using only the photometry as activity indicator (difference of about 16 in log, i.e., a factor of about 7 × 10⁶);

• The bissector span is used in the 29 best models, which indicates that this activity indicator is particularly helpful in predicting RVs;

• The FENRIR data-driven model seems to better predict RVs than the latent GPs model. However, it should be noted that the FENRIR data-driven model has a larger number of hyperparameters, since the amplitudes of all the harmonics of the rotation period are left free. While solar data allows one to constrain well these parameters, and to prevent overfitting (which would result in decreasing the cross-validation score), it might be more challenging to constrain them on smaller datasets.

• The models with a single mode (as defined in Section 4.3) perform well, even when multiple activity indicators (except for the photometry) are included in the models: it is sufficient to model the data as if there were only one type of “average active region”, which appears at random times.

• The Matérn 1/2 and 3/2 kernels seem to give similar predicting performances, while the computational cost of the Matérn 3/2 kernel is typically nine times larger.

6 Discussion

Our statistical model of stellar activity was built in three steps. First, we defined the FENRIR model in Section 2.3. Second, in Section 3, we made further assumptions on the impulse response, distribution of feature parameters, and rate of apparition of magnetic regions (Poisson rate). The validity of the assumptions in these two first steps is the object of Section 6.1. Third, we approximated FENRIR processes by a Gaussian distribution. This approximation is discussed in Section 6.2. In Appendix F, we present three additional discussions: the construction of relevant activity indicators, the potential relationship between phase shifts in RV and photometry with the faculae-to-spot ratio, and the application of our work to classical Doppler imaging.

6.1 Model assumptions: Limitations and novelty

In the FENRIR model, detailed in Section 2.3, the emergence times of magnetic regions follow a Poisson process (Eq. (6)). This means that the probability of a magnetic region appearing is independent of the presence of other active regions. However, this assumption does not fully hold, at least for the Sun, where active regions tend to emerge near existing ones and on active longitudes (e.g., Berdyugina & Usoskin 2003; Borgniet et al. 2015). Despite this, as discussed in Section 3.1, the issue is partially mitigated by considering impulse responses as groups of magnetic regions rather than individual ones. In our implementation, we account for spots shifted by 180°, which helps address this limitation.

We further specified the assumptions on the impulse response, distribution of feature parameters and Poisson rate in Section 3. Assumptions 1–6, outlined in Section 3.4, are common to the physics-based and data-driven models.

We assumed that the star has a solid body rotation (Eq. (23)). However, the angular velocity of the solar surface depends on the latitude, a phenomenon known as differential rotation (Gray 2005). With a preliminary analysis, we found that the firstorder effect of differential rotation is to reduce the time scale of the window kernel: the signal loses coherence more quickly.

Although the assumption that spot properties are independent of their position seems reasonable (Eq. (22)), it is known that their lifetime, size, and position depend on when they appear in the magnetic cycle (e.g., Song et al. 2024). First, at minimum solar activity, sunspots appear on average at a latitude 30-45° away from the equator. As the activity level rises, the average latitude of sunspots moves towards the equator. Second, the ratio of spots to faculae coverage varies with stellar activity (Nèmec et al. 2022), so that the parameters governing their effects should depend on time in a non-stationary way. This is further supported by the observation that if GP models are fitted to the solar HARPS-N RVs, the hyperparameters change significantly over time (Klein et al. 2024). While the general model can account for it (Eq. (7)), assumption 4 in Section 3.4 contradicts it. A possible workaround is to divide the data in shorter time intervals, such as months or years, and fit the model independently on each part.

In addition, we made several assumptions specific to the physics-based model (Sections 3.1–3.3). We computed the total RV as the sum of local stellar RVs weighted by their relative flux, which is simplistic. In Appendix B.3, we derived the expressions of RV, as well as ancillary indicators considering that the measured cross-correlation function (CCF) is a weighted sum of the local stellar CCFs, as is done in SOAP 2.0 (Dumusque et al. 2014). This formalism takes into account in particular the interdependence of velocity, contrast, and CCF width. Our model is also too simplistic for large spots, which cannot be assumed to be point-wise.

The lifetime of magnetic structures is correlated with their size, as documented by Meyer et al. (1974); Martinez Pillet et al. (1993). While our mathematical framework can accommodate this relationship, it is not currently incorporated into our physicsbased model implementation. Additionally, the influence of the magnetic network is not considered in our model (see Borgniet et al. 2015, for a discussion of its effect).

Our framework, despite its limitations, introduces several innovations. Traditionally, the representation of stellar activity in photometry and spectroscopy as a joint Gaussian Process (GP) relies on qualitative arguments (Aigrain et al. 2012; Gilbertson et al. 2020; Jones et al. 2022; Tran et al. 2023), resulting in the latent GP model in Eq. (39), possibly with higher order derivatives. In contrast, our approach offers a quantitative, physics-based model for understanding stellar activity in both photometry and spectra.

GP frameworks based on a quantitative physical models have been proposed. In Luger et al. (2021b,a) the goal is to retrieve an intensity map of the surface from photometric and spectroscopic measurements, respectively. In both cases, the GP represents the intensity pattern on the stellar surface and its uncertainties. The star rotates, but this pattern is considered fixed during the time span of the observations. Our framework focuses on modeling the evolution of stellar surfaces over time, and is more akin to Luger et al. (2021a). They start from a GP model of the stellar surface, and through a linear mapping, they derive a GP representation of photometric time series. Our framework has a different starting point, summarized by Eqs. (6)–(8). A comprehensive comparison of the two methods is beyond our current scope. We simply note that our framework extends to spectroscopy, accounts for the inhibition of the convective blueshift, incorporates variations due to the magnetic cycle, and enables likelihood evaluations scaling as O(N). Additionally, our approach features representations of non-Gaussianity, which we now explore.

6.2 Are stellar activity signals Gaussian?

A stochastic process X(t) is said to be a GP if for every finite collection of times t₁,…, t_n, then X(t₁),…, X(t_n) has a Gaussian multivariate distribution (Rasmussen & Williams 2005). The FENRIR GPs used in Section 5 satisfy this definition, but as mentioned in Section 2.4 they are only an approximation of the FENRIR models. We now discuss whether this is a good approximation.

There is a simple argument against stellar activity signals being strictly Gaussian. Indeed, if X(t) is a GP of vanishing mean, X(t) and –X(t) have exactly the same covariance, and thus the same probability. However, the signature of the inhibition of the convective blueshift on RV is not symmetrical, it always manifests as a net redshift (see Fig. 3 (b2)). Furthermore, as soon as there is an excess of the effect of regions brighter or darker than the quiet surface, the flux effect does not have a symmetric distribution. This type of non-Gaussianity could be called “spatial” and would manifest as a non-Gaussian histogram of RV and photometric observations.

Non-Gaussianity can also manifest in the temporal structure of the signal. If we were to generate several realizations of a GP, and compute the Fourier transform, it would appear that the phases are uniformly distributed on [0,2π]. In particular, phases at two different frequencies are statistically independent. In the case of stellar activity, we expect phases to be tied to each other. Suppose that only spots appear at a given latitude, and very rarely. This means that each time a spot appears, the difference between the phases of the process at the stellar rotation period and its harmonics is not only not uniform, but deterministic.

To study both spatial and temporal non-Gaussianity of stellar activity signals, we briefly introduce the notion of cumulant (for more details see Mendel 1991). Considering a stochastic process X(t), given n time stamps (t_i)_i=1,…,n the cumulant generating function of X(t₁),…, X(t_n) is $$K_{{\left(t_i\right)}_{i=1,\dots ,n}}^n\left(x_1,\dots ,x_n\right)=\ln \mathbb{E}\left\{{\text{e}}^{x_{1}X\left(t_1\right)+\dots +x_{n}X\left(t_n\right)}\right\}.$$(44)

In the Taylor expansion of $K_{{\left(t_i\right)}_{i=1,\dots n}}^n$ as a function of x_is, the cumu- lant of order n, κ_n, is the coefficient of their products, x₁ x₂…x_n. The cumulant κ_n is thus a function of (t_i)_i=1,…,n and is often referred to as the n-point correlation function. One of the properties of stationary GPs is that their cumulants of order three or greater are equal to zero.

A stationary process is such that the distribution of X(t + τ₁),…,X(t + τ_n) does not depend on t. When the FENRIR process is stationary, its n-point correlation function, or cumulant of order n is expressed as a function of time lags between n – 1 points and a reference one. In Appendix A.2, we show that the cumulant of order n can be analytically derived from the FENRIR model (Eqs. (6)–(8)), $$\begin{array}{l}{\kappa }_n\left({\tau }_1,\dots ,{\tau }_{n-1}\right)\\ =\lambda \int {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}𝑔}(t,\gamma )𝑔\left(t+{\tau }_1,\gamma \right)\dots 𝑔\left(t+{\tau }_{n-1},\gamma \right)p(\gamma \mid \eta )\text{d}t\text{\hspace{0.17em}d}\gamma .\end{array}$$(45)

Applying the Fourier transform to Eq. (45), we obtain the polyspectrum of our FENRIR process, $$\begin{array}{l}{S}_n\left({\omega }_1,\dots ,{\omega }_{n-1};\eta \right)\\ =\lambda {\displaystyle \iint \widehat{𝑔}}\left({\omega }_1,\gamma \right)\dots \widehat{𝑔}\left({\omega }_{n-1},\gamma \right)\widehat{𝑔}\left(-{\displaystyle \sum _{i=1}^{n-1}{\omega }_i},\gamma \right)p(\gamma \mid \eta )\text{d}\gamma ,\end{array}$$(46)

where $\widehat{g}(\omega ,\gamma )$ is the Fourier transform of g(t, γ) considered as a function of time t. Like cumulants, polyspectra of GPs vanish at order three and beyond. We now discuss “spatial” and “temporal” non-Gaussianity in Sections 6.2.1 and Sections 6.2.2.

6.2.1 Poisson rate and asymmetry

We might anticipate that increasing the rate of feature appearance enhances Gaussianity. At low Poisson rates, only one feature is present at a time. However, as the rate increases, multiple features begin to overlap, blurring phase information, and bringing the system closer to the conditions described by the central limit theorem.

As an illustration, we generated FENRIR processes using the same properties as in Section 3.3, with a sampling of one point per day over 20 years. We used two Poisson rates, λ: 1.2 and 0.2 magnetic regions appearing per day. The histograms of the simulated RVs are shown in Fig. 15. The RVs simulated with λ = 1.2 d^–1 appear closer to a Gaussian distribution than those with λ = 0.2 d^–1. To highlight this more clearly, we subtracted the mean from each simulated RV time series and normalized them by their standard deviations. Figure 16 displays the histogram of these normalized RVs. As anticipated, the distribution of RVs with λ = 1.2 d^–1 is closer to a normal distribution. As a remark, Fig. 15 shows that the variance of the simulated RVs increases with λ, which is expected from Eq. (11).

The FENRIR simulations, based on Section 3, are designed to approximate the characteristics of the Sun’s magnetic regions. These simulations can be compared with a three-year time series of RVs measured by HARPS-N (Dumusque et al. 2021). We binned these measurements every 15 hours and subtracted a 9th-order polynomial, thus removing frequencies lower than approximately 1 year, to avoid masking the contribution of stellar rotation. Fig. 17 displays the histogram of the residual RVs, which exhibits a similar behavior to the FENRIR processes: a sharp cut-off at negative RVs, heavier tails at positive RVs, and a maximum slightly shifted towards negative RVs. This behavior is independent of the degree of the fitted polynomial. However, as expected, when binning is done over shorter intervals, high- frequency noise tends to make the RV histogram slightly more Gaussian.

The λ = 1.2 d^–1 value is chosen to reproduce a Wolf number of ~50, and was calibrated to have the same standard deviation as the measured solar RVs, 1.23 m/s. Interestingly, while nothing was fine-tuned in the model to achieve this, the simulated and measured RVs have very similar empirical third-order moments: 1.02 m³/s³ and 0.91 m³/s³, respectively. Generating Gaussian noise with the same sampling as the solar 15-hour binned RVs, and using the same covariance (estimated via Fast Fourier Transform of solar RVs), yields a standard deviation of the third-order moment of 0.18 m³/s³. This indicates that the third-order moment of the measured RVs is non-zero with a 5σ significance and is accurately reproduced by the FENRIR model.

The fact that Gaussianity increases with λ can be qualitatively understood based on Eqs. (45) and (46). We consider a stationary FENRIR process, and its GP approximation. The two processes have identical first and second-order cumulants. Since the GP is stationary, from Eq. (45), its variance at each time σ² would be constant, proportional to λ. If we were to estimate empirically the cumulants of this GP, they would zero because of statistical fluctuation. The standar the cumulant of order n of a GP is proportional turn to λ^n/2, while all cumulants of order n are proportional to λ. Consequently, the ratio of the cumulants of order n ⩾ 2 of the FENRIR process and its GP approximation is proportional to λ/λ^n/2. This tends to zero when λ increases for all n ⩾ 3: the FENRIR process becomes more and more indistinguishable from its GP approximation.

6.2.2 Bispectrum

We now turn to “temporal” non-Gaussianity. Here, we do not study measured RVs, but only quantify the non-Gaussianity of FENRIR models. We use the polyspectrum of order 2, called the bispectrum. From Eq. (46), we obtain its expression for a FENRIR process, $$S_3\left({\omega }_1,{\omega }_2;\eta \right)=\lambda {\displaystyle \int s}\left({\omega }_1,{\omega }_2,\gamma \right)p(\gamma \mid \eta )\text{d}\gamma ,$$(47) $$s\left({\omega }_1,{\omega }_2,\gamma \right)=\widehat{𝑔}\left({\omega }_1,\gamma \right)\widehat{𝑔}\left({\omega }_2,\gamma \right)\widehat{𝑔}\left(-{\omega }_1-{\omega }_2,\gamma \right).$$(48)

The bispectrum is a function of the Fourier transform of the impulse response g, which can be expressed as a product of what we called a window function, modulating the amplitude of the signal as the region grows and decays in size, and a periodic term (see Eq. (23)). We first assume that the window term is constant (the size and brightness of the magnetic region are constant). The periodic part can be expressed as a truncated Fourier series $$g(t)={\displaystyle \sum _{k=-d}^d{a}_k}{\text{e}}^{\text{i}k\omega t}.$$(49)

This means that the Fourier transform of g is only non-zero at frequencies kω, k = –d, …, d, and the only non-zero products in the integrand of Eq. (47) are such that ω₁ = k₁ω and ω₂ = k₂ω₂ such that k₁, k₂, and –(k₁ + k₂) are all integers between –d and d.

Even if a cumulant is non-zero, its value estimated from noisy data could be compatible with zero. To evaluate this, we consider that the amplitudes of $\widehat{g}\left({\omega }_1\right),\widehat{g}\left({\omega }_2\right),\widehat{g}\left(-{\omega }_1-{\omega }_2\right)$ are estimated with an uncertainty. From there, we compute the mean squared error on s(ω₁, ω₂, γ) defined in Eq. (48).

In Fig. 18, we show the ratio of s and an upper bound on its mean squared error with the model of Section 3, assuming an inclination $\overline{i}=0$, a latitude δ = 0, and a ratio between the photometric RV and convective blueshift inhibition effects B = 1. We assume that $\widehat{g}(\omega )$ is known with a 15% accuracy for all ω, and the uncertainties on $\widehat{g}(\omega ),\widehat{g}\left({\omega }^{\prime }\right),\omega \ne {\omega }^{\prime }$, are independent. It appears that several bispectrum coefficients are more than three times greater than their mean squared error.

The bispectrum shown in Fig. 19 concerns a single magnetic region with a fixed latitude and a given stellar inclination. In Fig. 19, we show the real part of Eq. (48) for ω₁ = ω₂ = ω as a function of inclination $\overline{i}$ and latitude $\delta .\overline{i}=-{90}^{°}$ means that the one pole of the star is pointing to the observer, and $\overline{i}={90}^{°}$, the other pole, $\overline{i}=0^{°}$ means the stellar rotation axis is perpendicular to the line of sight. For each inclination, only regions satisfying $-\tan \overline{i}\tan \delta ⩽1$ are visible (see Appendix B). The bispectrum (Eq. (47)) requires to perform an integration. In the bottom of Fig. 19 we show the value of the integrand Eq. (48) averaged over latitude, assuming that all latitudes are equally probable. It appears that the real part of the bispectrum does not average out for most inclinations.

In the case where the window function (growth and decay of the spot area) is taken as non constant, thanks to the convolution theorem, the Fourier transform of g is $$\widehat{𝑔}(v)={\displaystyle \sum _{k=-d}^d{a}_k}\widehat{W}(v-k\omega ).$$(50)

Taking into account the window function would give a smoothed version of Fig. 18, but would not necessarily mask non- Gaussianity.

7 Conclusion

Our primary goal was to represent the signal induced by stellar activity as a stochastic process, derived from a quantitative physical model. We consider a general case where the data consist of observations of a given star in the form of one or several time series, which may vary in nature and are not necessarily sampled at the same epochs. These include radial velocities (RVs), photometry, activity indicators, or even timeseries of spectra used in transit spectroscopy or Doppler imaging. The FENRIR model, developed in this article, provides a joint likelihood representation of these observables: their probability distribution given the model parameters η, which include the sky-projected obliquity, average latitude of the spots, and rotation period.

The FENRIR model is built upon three key ingredients. (1) the impulse response: the effect of a given stellar feature (such as a spot or plage, or groups of them) as a function of its parameters (e.g., size, latitude). (2) The statistical distribution of feature parameters, given the average properties of the star, the parameters of the likelihood η. (3) The rate at which features appear, which may vary over the magnetic cycle.

The FENRIR model is generally not Gaussian, but it can be approximated by a Gaussian distribution. The resulting FENRIR GPs can be used like other GPs, here with hyperparameters η. We provide a public implementation in two flavors⁵, physicsbased and data-driven (see Section 3). In the first case, starting from a physical model, we specify the three ingredients of the model for the joint analysis of RV and photometry. The data- driven version considers the impulse response as a vector of free parameters, and can be used on any combination of observables.

FENRIR GPs are computed in the S+LEAF framework (Delisle et al. 2020, 2022). This enables fast likelihood evaluations, with a cost proportional to the number of observations.

Our goal is twofold: first, to improve the accuracy of stellar activity models, in order to improve exoplanet detection capabilities; and second, to constrain the statistical distribution of spots and faculae on slow rotating stars. To test the performance of FENRIR GPs, we analyzed data acquired on the Sun, SORCE photometry, and the first three years of HARPS-N spectroscopic observations (Dumusque et al. 2021). We tested the predictive accuracy of FENRIR models with cross-validation, on RV and photometry only, and on RV, FWHM, BIS, and photometry time series. In both cases, FENRIR models performed better than latent models with a GP and its derivatives. Second, using a separate process for spots and faculae, based solely on photometry and RVs, we retrieved the solar inclination with a precision of ~5°, as well as realistic spot and faculae distributions, close to the equator. This serves as a proof of concept, but further testing should be performed to ensure that our method gives reliable results on other stars.

The FENRIR framework is versatile. Here, we considered 16 variants of the physics-based model, but many more options are possible. It could be adapted to include differential rotation, and more subtle effects of spots and plages (see e.g., Dumusque et al. 2014). The FENRIR framework allows one to interpret the non-Gaussianity present in the solar data (see Section 6.2), and should, in the long run, allow one to go beyond the GP framework. In Appendix F, we also discussed several aspects: building relevant activity indicators, interpreting phase shifts between RV and photometry or log $R_{HK}^{\prime }$ as a signature of the ratio of spots and plages filling factors, using FENRIR in classical Doppler imaging. These avenues will be explored in future works.

Data availability

Appendices D (MCMC convergence tests and corner plots) and E (full table of cross-validation results) are available on Zenodo at https://doi.org/10.5281/zenodo.14919319.

Acknowledgements

This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. We acknowledge financial support from the Swiss National Science Foundation (SNSF), grant number: 205010. The authors warmly thank Xavier Dumusque, Vincent Bourrier, Christophe Lovis, Nadège Meunier, Lionel Garcia, Sophia Sulis and Olivier Flasseur for their helpful suggestions. We thank Suzanne Aigrain for her careful review of our article and her insightful comments.

Appendix A Covariance and cumulants

A.1 Covariance in the general nonstationary case

In this appendix, we compute the mean and autocovariance of FENRIR processes. These expressions are particular case of the general formula for cumulants, derived in Section A.2. Since the manipulation of cumulants might be unfamiliar to our readers, here we give an alternate calculation. For simplicity, we refer to the magnetic regions appearing on the star (spots and faculae, or groups of them) as features.

We consider that the features appear independently on a time interval [–L/2, L/2] with rate λ(t). Then the times of feature appearance follow the distribution $p(t)=\lambda (t)/{\displaystyle {\int }_{-L/2}^{L/2}\lambda }(t)\text{d}t$. The signal in a given observational channel is modeled as $$y^L𝑔(t)={\displaystyle \sum }_{k=1}^n\left(t-t_k,{\gamma }_k\right)$$(A.1)

and we denote by y(t) the limit in probability of y^L(t) as L tends to infinity (y^L(t) is a random variable). Here n is the number of features such that t_k ∈ [–L/2, L/2], and follows a Poisson distribution of parameter $I_{\lambda }^L:={\displaystyle {\int }_{-L/2}^{L/2}\lambda }(t)\text{d}t$. To compute the covariance of y, which by definition is $$\mathrm{Cov}\left(y\left(t_a\right),y\left(t_b\right)\right)=\mathbb{E}\left\{y\left(t_a\right)y\left(t_b\right)\right\}-\mathbb{E}\left\{y\left(t_a\right)\right\}\mathbb{E}\left\{y\left(t_b\right)\right\}.$$(A.2)

We first compute the expectancy of the value of y at a time t, $\mathbb{E}\{y(t)\}$, then $\mathbb{E}\left\{y\left(t_a\right)y\left(t_b\right)\right\}$. In both cases, we do the calculation for y^L and let L tend to infinity.

We first note that with our hypotheses, the probability of parameters {γ}_n = γ₁,…,γ_n knowing the times of arrival of features {t}_n = t₁,…, t_n can be written $$p\left({\left\{\gamma \right\}}_n\mid {\left\{t\right\}}_n\right)={\displaystyle \prod _{k=1}^{n}p}\left({\gamma }_k\mid t_k\right).$$

As a consequence, from the law of total probability, $$\begin{array}{l}\mathbb{E}\left\{y^L(t)\right\}={\displaystyle \sum _{n=1}^{+\infty }p}(n){\displaystyle \iint {\displaystyle \sum _{k=1}^n𝑔}}\left(t-t_k,{\gamma }_k\right)p\left({\left\{\gamma \right\}}_n\mid {\left\{t\right\}}_n\right)p\left({\left\{t\right\}}_n\right)\text{d}{\left\{}_{\gamma }\text{\hspace{0.17em}d}{\left\{}_t\\ ={\text{e}}^{-I_{\lambda }^L}{\displaystyle \sum _{n=1}^{+\infty }\frac{{\left(I_{\lambda }^L\right)}^n}{n!}}{\displaystyle \iint {\displaystyle \sum _{k=1}^n𝑔}}\left(t-t_k,{\gamma }_k\right)p\left({\gamma }_k\mid t_k\right)\frac{\lambda \left(t_k\right)}{I_{\lambda }^L}\text{\hspace{0.17em}d}{t}_k\text{\hspace{0.17em}d}{\gamma }_k\end{array}$$(A.3)

And we have $$\begin{array}{l}{\displaystyle \iint {\displaystyle \sum _{k=1}^n𝑔\left(t-t_k,{\gamma }_k\right)\frac{\lambda \left(t_k\right)}{I_{\lambda }^L}p\left({\gamma }_k|t_k\right)\text{d}{t}_k\text{d}{\gamma }_k=}}\\ \frac{n}{I_{\lambda }^L}{\displaystyle \iint 𝑔\left(t-t_0,\gamma \right)p\left(\gamma |t_0\right)\lambda \left(t_0\right)\text{d}{t}_0\text{d}\gamma .}\end{array}$$(A.4)

So $$\begin{array}{l}\mathbb{E}\left\{y^L(t)\right\}={\text{e}}^{-I_{\lambda }^L}{\displaystyle \sum _{n=1}^{+\infty }\frac{{\left(I_{\lambda }^L\right)}^{n-1}}{(n-1)!}}{{\displaystyle \iint }}_{-L/2}^{L/2}𝑔\left(t-t_0,\gamma \right)p\left(\gamma \mid t_0\right)\lambda \left(t_0\right)\text{d}{t}_0\text{\hspace{0.17em}d}\gamma \\ ={{\displaystyle \iint }}_{-L/2}^{L/2}𝑔\left(t-t_0,{\gamma }_k\right)\lambda \left(t_0\right)p\left(\gamma \mid t_0\right)\text{d}{t}_0\text{\hspace{0.17em}d}\gamma .\end{array}$$(A.5)

Provided this integral converges when L tends to infinity, $$\mathbb{E}\{y(t)\}={\displaystyle \iint 𝑔}\left(t-t_0,\gamma \right)p\left(\gamma \mid t_0\right)\lambda \left(t_0\right)\text{d}{t}_0\text{\hspace{0.17em}d}\gamma .$$(A.6)

We now compute the expectancy of the product y^L(t_a)y^L(t_b). Here also, provided the infinite sum and integral can be inverted, since for k ≠ j, t_k and t_j are statistically independent, $$\mathbb{E}\left\{y^L\left(t_a\right)y^L\left(t_b\right)\right\}={\text{e}}^{-I_{\lambda }^L}{\displaystyle \sum _{n=1}^{+\infty }\frac{{\left(I_{\lambda }^L\right)}^n}{n!}}\left({\displaystyle \sum _{k,j=1,k\ne j}^n{c}_{n,k,j}^L}+{\displaystyle \sum _{k=1}^n{d}_{n,k}^L}\right)$$(A.7) $$\begin{array}{l}{c}_{n,k,j}^L={{\displaystyle \iint }}_{-L/2}^{L/2}𝑔\left(t_a-t_k,{\gamma }_k\right)\frac{\lambda \left(t_k\right)}{I_{\lambda }^L}p\left({\gamma }_k\mid t_k\right)\text{d}{t}_k\text{\hspace{0.17em}d}{\gamma }_k\times \\ {{\displaystyle \iint }}_{-L/2}^{L/2}𝑔\left(t_b-t_j,{\gamma }_j\right)\frac{\lambda \left(t_j\right)}{I_{\lambda }^L}p\left({\gamma }_j\mid t_j\right)\text{d}{t}_j\text{\hspace{0.17em}d}{\gamma }_j\end{array}$$(A.8) $$d_{n,k}^L={{\displaystyle \iint }}_{-L/2}^{L/2}𝑔\left(t_a-t_k,{\gamma }_k\right)𝑔\left(t_b-t_k,{\gamma }_k\right)\frac{\lambda \left(t_k\right)}{{\left(I_{\lambda }^L\right)}^2}p\left({\gamma }_k\mid t_k\right)\text{d}{t}_k\text{\hspace{0.17em}d}{\gamma }_k$$(A.9)

Since there are n(n – 1) pairs of k, j where k ≠ j, $\mathbb{E}\left\{y^L\left(t_a\right)y^L\left(t_b\right)\right\}$ can be written as the sum of two terms C^L and D^L, $$\begin{array}{l}{C}^L\left(t_a,t_b\right)={\text{e}}^{-I_{\lambda }^L}{\displaystyle \sum _{n=2}^{+\infty }\frac{{\left(I_{\lambda }^L\right)}^n}{n!}}\frac{n(n-1)}{{\left(I_{\lambda }^L\right)}^2}\\ \begin{array}{c}\times {{{\displaystyle \iint }}^{\text{}}}_{-L/2}^{L/2}𝑔\left(t_a-t_0,\gamma \right)\lambda \left(t_0\right)p\left(\gamma \mid t_0\right)\text{d}{t}_0\text{d}\gamma \\ \times {{{\displaystyle \iint }}^{\text{}}}_{-L/2}^{L/2}𝑔\left(t_b-t_0\right)\lambda \left(t_0\right)p\left(\gamma \mid t_0\right)\text{d}{t}_0\text{d}\gamma \end{array}\end{array}$$(A.10)

and $$D^L\left(t_a,t_b\right)={{\displaystyle \iint }}_{-L/2}^{L/2}𝑔\left(t_a-t_0\right)𝑔\left(t_b-t_0\right)\lambda \left(t_0\right)p\left(\gamma \mid t_0\right)\text{d}{t}_0\text{\hspace{0.17em}d}\gamma .$$(A.11)

We note that as L tends to infinity, provided the integral over t_k converges, C^L(t_a, t_b) is the product of $\mathbb{E}\left\{y\left(t_a\right)\right\}$ and $\mathbb{E}\left\{y\left(t_b\right)\right\}$. If D^L(t_a, t_b) admits a limit as L tends to infinity, the covariance of y is $$\mathrm{Cov}\left(y\left(t_a\right),y\left(t_b\right)\right)={\displaystyle \iint 𝑔}\left(t_a-t_0,\gamma \right)𝑔\left(t_b-t_0,\gamma \right)\lambda \left(t_0\right)p\left(\gamma \mid t_0\right)\text{d}{t}_0\text{\hspace{0.17em}d}\gamma$$(A.12)

In this formula, the distribution of γ can depend on t_o.

We note that if we replace y(t_b) with z(t_b) where z models the effect of stellar activity on another ancillary indicator, $$ɀ(t)={\displaystyle \sum _{k=-\infty }^{+\infty }h}\left(t-t_k,\gamma \left(t_k\right)\right)$$(A.13)

the reasoning is unchanged, and provided all the limits when L tends to infinity converge, $$\mathrm{Cov}\left(y\left(t_a\right),z\left(t_b\right)\right)={\displaystyle \iint 𝑔}\left(t_a-t,\gamma \right)h\left(t_b-t,\gamma \right)\lambda (t)p(\gamma \mid t)\text{d}t\text{\hspace{0.17em}d}\gamma \text{\hspace{0.17em}d}\gamma .$$(A.14)

A.2 Cumulants and polyspectra

A.2.1 Purpose

In this Section, we compute the analytical formula of the cumulants and polyspectra of the FENRIR process y(t), defined as $$y(t)={\displaystyle \sum _{k=-\infty }^{+\infty }𝑔}\left(t-t_k,\gamma \left(t_k\right)\right).$$(A.15)

For a stationary FENRIR process, the autocovariance is of the form $k(\tau )={\displaystyle {\int }_{-\infty }^{+\infty }g}(t,\gamma )g(t+\tau ,\gamma )p(\gamma )\text{d}t$. Even if the γ can take only one value, the knowledge of k(τ) does not give unequivocally g: several functions can have the same autocorrelation. If y(t) is Gaussian, based on it the best we can do is to estimate k(τ). If y(t) is non-Gaussian we can characterize g more precisely, thus have not only less biases but also leverage more information. The purpose of cumulants is to reveal non-Gaussian behaviors.

The cumulant generating function of y(t₁), y(t₁)…, y(t_n) is $$K_{t_1,t_2,\dots ,t_n}^n\left({\alpha }_1,\dots ,{\alpha }_n\right)=\ln \mathbb{E}\left\{{\text{e}}^{{\alpha }_{1}y\left(t_1\right)+\dots +{\alpha }_{n}y\left(t_n\right)}\right\}.$$(A.16)

In the Taylor expansion of K as a function of α_is, the cumulant of order n, κ_n, is the coefficient of their products, α₁α₂…α_n. Considered as a function of t₁,t₂,…, t_n it is sometimes called the n point correlation function. For a GP, cumulants of order three and four are zero. Our main result is that

Theorem 1. For a FENRIR process y(t) defined as in Eqs. (6)- Eqs. (8), the cumulant of order n is $$\begin{array}{l}{\kappa }_n\left(y\left(t_1\right),y\left(t_2\right),\dots ,y\left(t_n\right)\right)=\\ {\displaystyle \iint 𝑔}\left(t-t_1,\gamma \right)g\left(t-t_2,\gamma \right)\dots 𝑔\left(t-t_n,\gamma \right)\lambda (t)p(\gamma ,t)\text{d}t\text{\hspace{0.17em}d}\gamma .\end{array}$$(A.17)

The covariance is the 2-point correlation function. When GPs are stationary, we have κ_n(t₁, t₂) = κ_n(|t₁ – t₂|) they are equivalently represented by the Fourier transform of the kernel κ(τ): their power spectrum density (PSD). Similarly, under the hypothesis of stationarity of the process, the quantity κ_n(t, t + τ₁,…, t + τ_n–1) does not depend on t. Considering the n point correlation function as a function of the n – 1 variables τ₁, …,τ_n–1, we can define the polyspectra as the n – 1 dimensional Fourier transform of κ_n(t, t + τ₁,…, t + τ_n–1).

The proof of Eq. (A.17) can be modified straightforwardly to prove a more general result, $$\begin{array}{l}{\kappa }_n\left(y_1\left(t_1\right),\dots ,y_n\left(t_n\right)\right)=\\ {\displaystyle \iint {𝑔}_1}\left(t-t_1,\gamma \right)\dots {𝑔}_n\left(t-t_n,\gamma \right)\lambda (t)p(\gamma ,t)\text{d}t\text{\hspace{0.17em}d}\gamma .\end{array}$$(A.18)

where g_i(t) are impulse responses associated with channel i. This expression generalizes Eq. (A.14). This expression has the advantage to be valid if the process is non stationary, and can explore non-Gaussian dependencies across channel. However, it is unlikely that all these aspects could be explored simultaneously in a practical way.

In the case where patterns appear with a constant Poisson rate, and the distribution followed by γ(t_k) does not depend on time, γ(t_k) ~ p(γ), the stochastic process y(t) is stationary. We can write Eq. (A.17) as a function of time differences τ_i = t_i – t₁. With the change of variable t ← t – t₁, $$\begin{array}{l}{\kappa }_n\left({\tau }_1,\dots ,{\tau }_{n-1};\eta \right)=\\ \lambda 𝑔(t,\gamma )𝑔\left(t+{\tau }_1,\gamma \right)\dots g\left(t+{\tau }_{n-1},\gamma \right)p(\gamma \mid \eta )\text{d}t\text{\hspace{0.17em}d}\gamma ,\end{array}$$(A.19)

In that case, the Fourier transform of κ_n as a function of τ₁, …,τ_n–1 is called the polyspectrum of order n. Here, provided the integrals on γ and τ₁,…, τ_n–1 can be inverted, it is equal to $$\begin{array}{l}{S}_n\left({\omega }_1,\dots ,{\omega }_{n-1};\eta \right)=\\ \lambda {\displaystyle \int \widehat{𝑔}}\left({\omega }_1,\gamma \right)\dots \widehat{𝑔}\left({\omega }_{n-1},\gamma \right)\widehat{𝑔}\left(-{\displaystyle \sum _{i=1}^{n-1}{\omega }_i},\gamma \right)p(\gamma \mid \eta )\text{d}\gamma .\end{array}$$(A.20)

The purpose of polyspectra can be understood by recalling that stationary GPs are completely characterized by their power spectral density, which is the Fourier transform of the kernel. Under mild conditions (Bloomfield 2004), the power spectral density has another interpretation: It is the mathematical expectancy of the squared modulus of the Fourier transform of the process. Stationary GPs thus have completely random phases. Suppose we compute the phase of a stationary GP between times t₁ and t₂ on the one hand and between times t₃ and t₄ on the other hand, such that t₃ – t₂ is much greater than the process correlation, the phase obtained on the first and second time intervals are statistically independent. For stationary GPs, for all n ⩾ 3, Eq. (A.20) should be equal to zero for all choice of frequencies ω₁,…, ω_n–1. As a consequence, significantly non-zero polyspectra indicate non-Gaussianity. We consider a quantity ν and all combinations of n – 1 frequencies such that ${\displaystyle \sum _{i=1}^{n-1}{\omega }_i}=v$. In Eq. (A.20), the phase of the Fourier transform evaluated at ν can be seen as a point of reference. If the frequencies ω₁, …,ω_n–1 and v are not independent, the poly spectrum is non-zero.

A.2.2 Proof

To establish Eq. (A.17), we use several properties of cumulants, outlined below.

Theorem 2. We consider n random variables Y₁, Y₂,…, Y_n (defined on the same σ algebra). (i) Cumulants are multilinear, for all i, for two real numbers a and b, $\kappa \left(Y_1,\dots ,a{Y}_i+b{Y}_i^{\prime },\dots ,Y_n\right)=a\kappa \left(Y_1,\dots ,Y_i,\dots ,X_n\right)+b\kappa \left(Y_1,\dots ,Y_i^{\prime },\dots ,Y_n\right)$.

(ii) If there are at least two indices i , j such that Y_i and Y_j are independent, κ(Y₁,…, Y_i,…, Y_j,…, Y_n) = 0.

(iii) If N is a random variable following a Poisson distribution of parameter Λ, κ(N, N…, N) = Λ regardless of the number of times N is repeated.

(iv) Law of total cumulance. Suppose that we have a random variable Y and we can define the conditional distribution of Y_is knowing X. Then $$\kappa \left(Y_1,\dots ,Y_n\right)={\displaystyle \sum _{\pi \in \mathbb{P}(n)}\kappa }\left(\kappa \left(Y_i:i\in B\mid X\right):B\in \pi \right)$$(A.21)

where $\pi \in \mathbb{P}(n)$ means that π runs through all the possible partitions of indices of indices {1,…, n}. B ∈ π means that B runs through the blocks of permutation π and κ(Y_i : i ∈ B | X) is the cumulant of Y_i : i ∈ B knowing X . To understand (iv), suppose we want to compute κ(Y₁, Y₂, Y₃). The partitions of our indices {1, 2, 3} are {{1}, {2}, {3}}, {{1}, {2, 3}}, {{2}, {1, 3}}, {{3}, {1, 2}}, {{1, 2, 3}}.

In the partition {{1}, {2}, {3}} we have three blocks, {1}, {2} and {3}. In the partition {{1}, {2, 3}} we have two blocks, {1} and {2, 3}. For n = 3, the law of total cumulance then writes $$\begin{array}{l}\kappa \left(Y_1,Y_2,Y_3\right)=\kappa \left(\kappa \left(Y_1,Y_2,Y_3\mid X\right)\right)\hfill \\ +\kappa \left(\kappa \left(Y_1,Y_2\mid X\right),\kappa \left(Y_3\mid X\right)\right)\hfill \\ +\kappa \left(\kappa \left(Y_1,Y_3\mid X\right),\kappa \left(Y_2\mid X\right)\right)\hfill \\ +\kappa \left(\kappa \left(Y_2,Y_3\mid X\right),\kappa \left(Y_1\mid X\right)\right)\hfill \\ +\kappa \left(\kappa \left(Y_1\mid X\right),\kappa \left(Y_2\mid X\right),\kappa \left(Y_3\mid X\right)\right)\hfill \end{array}$$

where the quantities κ(Yi : i ∈ B | X) are random variables because they are functions of the random variable X .

(v) With the same notations, the link between cumulants and non centred moments is the following $$\mathbb{E}\left\{Y_1{Y}_2\dots Y_3\right\}=\underset{\pi \in \mathbb{P}(n)}{{{\displaystyle \sum }}^{\text{}}}\underset{B\in \pi }{{{\displaystyle \prod }}^{\text{}}}\kappa \left(Y_i:i\in B\right).$$(A.22)

For instance $$\begin{array}{ll}\mathbb{E}\left\{Y_1{Y}_2{Y}_3\right\}\hfill & =\kappa \left(Y_1,Y_2,Y_3\right)\hfill \\ \hfill & +\kappa \left(Y_1,Y_2\right)\kappa \left(Y_3\right)\hfill \\ \hfill & +\kappa \left(Y_1,Y_3\right)\kappa \left(Y_2\right)\hfill \\ \hfill & +\kappa \left(Y_2,Y_3\right),\kappa \left(Y_1\right)\hfill \\ \hfill & +\kappa \left(Y_1\right)\kappa \left(Y_2\right)\kappa \left(Y_3\right)\hfill \end{array}$$

Formulae (A.21) and (A.22) are similar. In the first case one takes the cumulants over the different blocks of the distribution, in the second one simply takes the product. This will be useful in the proof of Eq. (A.17).

Proof of Eq. (A.17): We consider n times ordered increasingly T₁,…, T_n. We suppose that g has a finite support, meaning that for some T, for t outside [–T/2, T/2], g(t, γ) = 0. Because g has a finite support, computing the n-th order cumulant of y(T₁), …,y(T_n) for y as defined in Eq. (A.15) is equivalent to compute the cumulant of the variables Y_i $$Y_i={\displaystyle \sum _{k=1}^N𝑔}\left(T_i-t_k,\gamma \left(t_k\right)\right)={\displaystyle \sum _{k=1}^{N}G}\left(T_i-t_k\right).$$(A.23)

where features appear at t_k, following a Poisson process between t₁ – T/2 and t_n + T/2 with a variable rate λ(t), and we define the random variables G(T_i – t_k) = g(T_i – t_k, γ(t_k)) to simplify notations. N is the number of features in this interval and follows a Poisson distribution with parameter $\mathrm{\Lambda }={\displaystyle {\int }_{t_1-T/2}^{t_n+T/2}\lambda }(t)$. A feature appearing at time t_k has parameters drawn from the distribution p(γ, t).

In this expression, k(Gi : i ∈ B) are not random variables, but constants. Thanks to the multilinearity of cumulants (property (i) in Theorem (2)), we have $$\kappa \left(Y_1,\dots ,Y_n\mid N\right)={\displaystyle \sum _{i_1=1}^N\dots }{\displaystyle \sum _{i_n=1}^N\kappa }\left(G\left(T_1-t_{i_1}\right),\dots ,G\left(T_n-t_{i_n}\right)\mid N\right)$$(A.24)

Measurable functions of independent variables are independent. Since ti and tj are independent for i ≠, j, G(T₁ – t_i) and G(T₁ – t_j) are not independent if and only if i = j. Now, because cumulants are zero if two or more variables are independent (property (ii) in Theorem (2)), $$\kappa \left(G\left(T_1-t_{i_1}\right),\dots ,G\left(T_n-t_{i_n}\right)\mid N\right)\ne 0$$(A.25)

if and only if i₁ = i₂ = … i_n, and $$\kappa \left(Y_1,\dots ,Y_n\mid N\right)={\displaystyle \sum _{i=1}^N\kappa }\left(G\left(T_1-t_i\right),\dots ,G\left(T_n-t_i\right)\mid N\right)$$(A.26)

The t_i, i = 1..Ns are drawn independently from the same distribution, so for a given k, all the G(T_k – t_i), i = 1..N are drawn from the same distribution: κ(G(T₁ – t_i),…, G(Tn - t_i) | N) = κ(G(T₁ – t_j),…, G(T_n – t_j) | N) for all i, j = 1..N and $$\kappa \left(Y_1,\dots ,Y_n\mid N\right)=N\kappa \left(G\left(T_1-t_i\right),\dots ,G\left(T_n-t_i\right)\mid N\right)$$(A.27)

Thanks to the law of total cumulance (property (iv) in Theorem 2), $$\kappa \left(Y_1,\dots ,Y_n\right)={\displaystyle \sum _{\pi \in \mathbb{P}(n)}\kappa }\left(\kappa \left(Y_i:i\in B\mid N\right):B\in \pi \right)$$(A.28)

Injecting Eq. (A.27), denoting by G_k = G(T_k – t), we have $$\kappa \left(Y_1,\dots ,Y_n\right)={\displaystyle \sum _{\pi \in \mathbb{P}(n)}\kappa }\left(N\kappa \left(G_i:i\in B\right):B\in \pi \right)$$(A.29)

Here κ(G_i : i ∈ B) is a constant, not a random variable. Thanks to the multilinearity of cumulants (i), we have $$\kappa \left(N\kappa \left(G_i:i\in B\right):B\in \pi \right)=\left({\displaystyle \prod _{B\in \pi }\kappa }\left(G_i:i\in B\right)\right){\kappa }_{\left|B_{\pi }\right|}(N)$$(A.30)

where ${\kappa }_{\left|B_{\pi }\right|}N=\kappa (N,\dots ,N)$ where N is repeated |B_π| times, the number of blocks in partition π. For a Poisson distribution of parameter Λ, κ(N,…, N) = Λ (property (iii) in Theorem (2)). We now have $$\kappa \left(Y_1,\dots ,Y_n\right)=\mathrm{\Lambda }{\displaystyle \sum _{\pi \in \mathbb{P}(n)}{\displaystyle \prod _{B\in \pi }\kappa }}\left(G_i:i\in B\right)$$(A.31)

and thanks to the link between moments and cumulants (property (v) in Theorem (2)) we can write $$\kappa \left(Y_1,\dots ,Y_n\right)=\mathrm{\Lambda }\mathbb{E}\left\{G\left(T_1-t\right)\dots G\left(T_n-t\right)\right\}$$(A.32)

Finally, by definition of $\mathbb{E}\left\{G\left(T_1-t\right)\dots G\left(T_n-t\right)\right\}$, $$\begin{array}{l}\mathbb{E}\left\{G\left(T_1-t\right)\dots G\left(T_n-t\right)\right\}=\\ \frac{1}{\mathrm{\Lambda }}{\displaystyle \underset{t_1-T/2}{\overset{t_n+T/2}{\int }}𝑔}\left(T_1-t,\gamma (t)\right)\dots 𝑔\left(T_n-t,\gamma (t)\right)\lambda (t)p(\gamma \mid t)\text{d}t\text{\hspace{0.17em}d}\gamma \end{array}$$(A.33)

canceling Λ in Eq. (A.32), and recognizing that because of the finite support of g, the integral (A.33) is unchanged if we take as lower and upper bounds –∞ and +∞, and is true for any T, we can thus extend our results to functions with infinite support and we have the desired result, $$\begin{array}{l}\kappa \left(Y_1,\dots ,Y_n\right)=\\ {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}𝑔}\left(T_1-t,\gamma (t)\right)\dots 𝑔\left(T_n-t,\gamma (t)\right)\lambda (t)p(\gamma \mid t)\text{d}t\text{\hspace{0.17em}d}\gamma .\end{array}$$(A.34)

Appendix B Analytical approximation of the spot/facula RV effect

B.1 Weighted radial velocity (RV)

In this appendix we approximate the RV and photometric effect as follows. We assume that the RV effect on the data of the quiet stellar surface is $$R{V}_{quiet}(t)=\frac{1}{F_0}{\displaystyle \iint F}(x)V(x)\text{d}x\cdot e_{\text{obs\hspace{0.17em}}}$$(B.1)

where x is the position on the visible disk, V(x) is the local velocity, F(x) the local flux, F₀ the flux integrated on all the stellar surface and e_obs the unit vector pointing from the observer to the star. Where $$V(x)=V_{cb}{e}_r+\omega \wedge r$$(B.2)

where V_cb is the velocity modulus of the convective blueshift effect, e_r is pointing radially outwards, ω is the rotation axis times the rotational velocity and r is the vector joining the stellar center and the point at position x on the stellar surface.

Denoting by RV_mag(t) the RV when a magnetic region is present, and by V′(x) and F′ (x) the velocity and flux fields , $$R{V}_{mag}(t)=\frac{1}{F_0^{\prime }}{\displaystyle \iint F^{\prime }}(x)V^{\prime }(x)\text{d}x\cdot e_{obs}$$(B.3)

where $$V^{\prime }(x)=\left(V_{cb}+\mathrm{\Delta }{V}_{cb}\right)e_r+\omega \wedge r$$(B.4) $$F^{\prime }(x)=F(x)+\mathrm{\Delta }F(x).$$(B.5)

We are interested in the difference between RV_mag(t) and RV_quiet(t). Developing $F_0^{\text{'}}$ at first order in ΔF, we have $$\left(\frac{F^{\prime }(x)V^{\prime }(x)}{F_0^{\prime }}-\frac{F(x)V(x)}{F_0}\right)\cdot e_{obs}=y_{cb}+y_{ph}$$(B.6) $$y_{cb}=\left(\mathrm{\Delta }F{V}_{cb}\left(1-\frac{\mathrm{\Delta }F}{F_0}\right)+\mathrm{\Delta }{V}_{cb}(F+\mathrm{\Delta }F)\right)e_r\cdot e_{obs}$$(B.7) $$y_{ph}=\mathrm{\Delta }F\omega \wedge r\cdot e_{obs}$$(B.8)

The term y_cb and y_ph are respectively called the inhibition of convective blueshift term RV photometric terms.

Assuming that the magnetic region is small, the flux can be written as product of the flux per surface unit f (x) times the projected area of the magnetic region, W(x)P(x) where W(x) is the intrinsic area and P(x) the projection effect, times the limbdarkening effect l(x). We write $$F(x)=f(x)W(x)P(x)l(x),$$(B.9) $$\mathrm{\Delta }F(x)=\mathrm{\Delta }F(x)W(x)P(x)l(x).$$(B.10)

In the following, we establish the expression of the flux and RV as a function of the position of the magnetic region (longitude ϕ and latitude δ) on an inclined star. For small spots, we can neglect the ΔF/F₀ term in Eq. (B.7) (it is second order).

We model the surface of the star by a sphere and consider that a spot or facula is point-like. For the sake of simplicity, in the following we simply refer to a spot, but the reasoning for a facula is identical. We consider a direct frame x, y, ɀ such that x points in the direction of the observer and y, ɀ defines the sky plane.

We first assume that the rotation axis of the star is aligned with ɀ denote by ϕ, δ the spherical coordinate of the center of a spot, such that its position in the Cartesian frame is $$x=\cos \delta \cos \varphi$$(B.11) $$y=\cos \delta \sin \varphi$$(B.12) $$ɀ=\sin \delta$$(B.13)

The local frame (u, υ, w) at (x, y, ɀ) is such that $$u=\cos \delta \cos \varphi x+\cos \delta \sin \varphi y+\sin \delta ɀ$$(B.14) $$\upsilon =-\sin \varphi x+\cos \varphi y$$(B.15) $$w=-\sin \delta \cos \varphi x-\sin \delta \sin \varphi y+\cos \delta ɀ$$(B.16)

The effect of a spot on the RV depends on the inclination of the star with respect to the plane of the sky $\overline{i}$. To compute the position of the spot and the projection of the local frame centered at the spot in the reference frame, we apply a rotation of axis y and angle $\overline{i}\in [0,\pi /2]$. We now have $$x=\cos \overline{i}\cos \delta \cos \varphi +\sin \overline{i}\sin \delta$$(B.17) $$y=\cos \delta \sin \varphi$$(B.18) $$ɀ=-\sin \overline{i}\cos \delta \cos \varphi +\cos \overline{i}\sin \delta$$(B.19)

and $$\begin{array}{l}u=(\cos \overline{i}\cos \delta \cos \varphi +\sin \overline{i}\sin \delta )x\\ +\cos \delta \sin \varphi y+(-\sin \overline{i}\cos \delta \cos \varphi +\cos \overline{i}\sin \delta )ɀ\end{array}$$(B.20) $$\upsilon =-\cos \overline{i}\sin \varphi x+\cos \varphi y+\sin \overline{i}\sin \varphi y$$(B.21) $$\begin{array}{l}w=-\cos \overline{i}\sin \delta \cos \varphi x-\sin \delta \sin \varphi y\\ +(\sin \overline{i}\sin \delta \cos \varphi +\cos \overline{i}\cos \delta )ɀ\end{array}$$(B.22)

We assume that the effect of the spot is proportional to its projected area onto the sky plane times a velocity projected onto the x axis (more precisely -x axis, since the velocity is assumed to be positive in the direction observer - star) multiplied by a limb-darkening effect. In the case of the photometric effect, the velocity in question is the velocity of the stellar surface due to the stellar rotation in the rest frame, it is therefore the υ component of the velocity. The convective blueshift inhibition effect is due to the motion of the gas from the center of the star to the stellar surface. We are therefore interested in the u component of the velocity projected onto x. The spot position on the stellar surface changes its projected surface onto the sky plane. The correction factor is equal to the Jacobian of the projection from the (υ, w) onto the (y, ɀ) plane, from Eq. (B.21) and Eq. (B.22), the projected area of the spot is proportional to $$\begin{array}{l}J(\overline{i},\delta ,\varphi )=\cos \varphi (\sin \overline{i}\sin \delta \cos \varphi +\cos \overline{i}\cos \delta )\\ -(-\sin \delta \sin \varphi )\sin \overline{i}\sin \varphi \\ =\sin \overline{i}\sin \delta +\cos \overline{i}\cos \delta \cos \varphi .\end{array}$$(B.23)

The factor J actually coincides with the distance to the center of the star normalized by the stellar radius, µ. In the following we keep the notation µ, which is widely used for limb-darkening laws. The spot is visible if its position is such that x ⩾ 0. From Eq. (B.17), this translates to the condition $$\cos \varphi ⩾-\tan \overline{i}\tan \delta$$(B.24)

Overall, the velocity contribution of the photometric and inhibition of convective blueshift effects on RV (y_ph and y_cb , and its effect on photometry ɀ_ph is 0 when (B.24) is not satisfied, and when it is: $$y_{ph}(\overline{i},\delta ,\varphi )=W(t)l(\mu )\mathrm{\Delta }F\omega R^{\star }\mu \cos \overline{i}\sin \varphi \cos \delta$$(B.25) $$y_{cb}(\overline{i},\delta ,\varphi )=W(t)l(\mu )\left(\mathrm{\Delta }F{V}_{cb}+\mathrm{\Delta }{V}_{cb}(f+\mathrm{\Delta }F)\right){\mu }^2$$(B.26) $${ɀ}_{ph}(\overline{i},\delta ,\varphi )=W(t)l(\mu )\mathrm{\Delta }F\mu$$(B.27)

We can parametrize the intrinsic area of the feature with a timescale ρ. For instance, if we choose a one-sided exponential, W(t) = 1_t⩾0e^−t/ρ where 1_t⩾0 equals 0 for negative t, and one for positive t. In the remainder of this appendix, we ignore the W(t) term.

ΔF is the flux difference between the stellar surface and the spot and f^★ is the mean stellar flux, ω is the local rotational velocity of the star and R^★ is the stellar radius, Δ V_CB is the difference between the velocity of the plasma in the u direction with and without convective blueshift inhibition. The + sign comes from the fact that we are projecting u and v onto −x, and the inhibition of convective blueshift has a positive velocity in the −u direction. The photometric effect can be positive or negative depending on whether the area under consideration is brighter of darker than the quiet surface of the star. The symbol l represents the limb-darkening, which is typically a function of µ. For an equatorial spot, $\overline{i}=0$, and without limb-darkening, we simply have $$y_{ph}\propto \cos \varphi \sin \varphi =\frac{1}{2}\sin 2\varphi$$(B.28) $$y_{cb}\propto {\cos }^2\varphi =\frac{1}{2}(1+\cos 2\varphi )$$(B.29)

In the case of a non constant limb-darkening law l(µ), the expressions (B.25) and (B.26) are particularly simple if l is of the form $$l(\mu )={\displaystyle \sum _{k=0}^d{a}_k}{\mu }^k.$$(B.30)

In the general case the photometric and convective blueshift inhibition effect are summed. For the sake of simplicity we write their contribution up to a multiplicative factor and denote the relative contribution of the photometric effect compared to the convective blueshift inhibition as B. The combined effect on the RV is then $$𝑔(t)={\displaystyle \sum _{k=0}^d{a}_k}{\mu }^{k+1}(\mu +B\cos \overline{i}\sin \varphi )$$(B.31)

for ϕ satisfying (B.24) and 0 otherwise.

When the limb darkening law is taken as a constant, the general formula for the RV effect of a spot is $$\begin{array}{l}𝑔(\varphi )\propto {\sin }^2\overline{i}{\sin }^2\delta +\frac{{\cos }^2\overline{i}{\cos }^2\delta }{2}\\ +\frac{\sin 2\overline{i}\sin 2\delta }{2}\cos \varphi +\frac{{\cos }^2\overline{i}{\cos }^2\delta }{2}\cos 2\varphi \\ -B\frac{\sin 2\overline{i}\sin 2\delta }{4}\sin \varphi -B\frac{{\cos }^2\overline{i}{\cos }^2\delta }{2}\sin 2\varphi \end{array}$$(B.32)

For a constant rotation rate, constant δ, we have ϕ = ωt. We can rewrite g as $$\begin{array}{l}𝑔(t,\overline{i},\delta ,B)=c_0+c_1\cos \omega t+c_2\cos 2\omega t\\ +s_1\sin \omega t+s_2\sin 2\omega t\end{array}$$(B.33)

where c₀ , c₁ , c₂ , s₁ , s₂ are coefficients that depend on $\overline{i}$, δ and B. As a concluding remark, we note that because the limbdarkening is in power of µ, itself an affine function of cos ωt. Regardless of the order of the limb-darkening law chosen, we can always write g in the form $$𝑔(t,\overline{i},\delta ,B)=c_0+{\displaystyle \sum _{k=1}^d{c}_k}\cos (k\omega t)+s_k\sin (k\omega t)$$(B.34)

The next section provides a formal justification of this claim.

B.2 Fourier decomposition of the periodic part

We denote by I and J the impulse responses in channels i and j. For simplicity, we write them I(t, γ) = W(t, γ)i(t, γ) and J(t, γ) = W(t, γ) j(t, γ) where i(ϕ, γ) and j(ϕ, γ) are the periodic parts, which can be written in the form $$\begin{array}{l}i(\varphi ,\gamma )=1_{\text{vis.\hspace{0.17em}}}(\varphi ,\gamma )i_0(\varphi ,\gamma ),\\ j(\varphi ,\gamma )=1_{\text{vis.\hspace{0.17em}}}(\varphi ,\gamma )j_0(\varphi ,\gamma ),\end{array}$$(B.35)

where 1_vis. is the indicator function of the spots visibility (equals to one when the spot is on the observer’s side of the star, 0 otherwise), and i₀, j₀ describe the effect of the feature without visibility considerations. The Fourier expansions of the effect of the feature, including the limb-darkening effect, are provided in Appendix B. Overall, this can be very well approximated with a few harmonics of the rotation period (typically less than 5). We denote by i_0,k , j_0,k these expansions up to a given harmonics d: $$\begin{array}{l}{i}_0(\varphi ,\gamma )={\displaystyle \sum _{k=-d}^d{i}_{0,k}}(\gamma )\exp (\text{i}k\varphi ),\\ j_0(\varphi ,\gamma )={\displaystyle \sum _{k=-d}^d{j}_{0,k}}(\gamma )\exp (\text{i}k\varphi ).\end{array}$$(B.36)

As shown in Appendix B, the feature is visible when its longitude ϕ is such that $$\cos \varphi \ge -\tan i\tan \delta .$$(B.37)

We denote by ϕ_m the maximum allowed value of ϕ in the range [0, π]; that is, $${\varphi }_m=\{\begin{array}{c}0,\text{\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}the\hspace{0.17em}spot\hspace{0.17em}is\hspace{0.17em}never\hspace{0.17em}visible,\hspace{0.17em}i.e.,\hspace{0.17em}}-\tan i\tan \delta \ge 1\text{,\hspace{0.17em}}\\ \pi ,\text{\hspace{0.17em}\hspace{0.17em}if\hspace{0.17em}the\hspace{0.17em}spot\hspace{0.17em}is\hspace{0.17em}always\hspace{0.17em}visible,\hspace{0.17em}i.e.,\hspace{0.17em}}-\tan i\tan \delta \le -1\text{,\hspace{0.17em}}\\ \arccos (-\tan i\tan \delta ),\text{\hspace{0.17em}\hspace{0.17em}otherwise.\hspace{0.17em}}\end{array}$$

The Fourier coefficients of i are then given by $$\begin{array}{l}{i}_l(\gamma )=\frac{1}{2\pi }{\displaystyle {\int }_{-\pi }^{\pi }{1}_{\text{vis}.}}(\varphi ,\gamma )i_0(\varphi ,\gamma )\exp (-\text{i}l\varphi )\text{d}\varphi \\ =\frac{1}{2\pi }{\displaystyle {\int }_{-{\varphi }_{\text{m}}}^{{\varphi }_{\text{m}}}{i}_0}(\varphi ,\gamma )\exp (-\text{i}l\varphi )\text{d}\varphi \\ =\frac{1}{2\pi }{\displaystyle \sum _{k=-d}^d{i}_{0,k}}(\gamma ){\displaystyle {\int }_{-{\varphi }_{\text{m}}}^{{\varphi }_{\text{m}}}\exp }(\text{i}(k-l)\varphi )\text{d}\varphi \\ =\frac{{\varphi }_{\text{m}}}{\pi }{\displaystyle \sum _{k=-d}^d\mathrm{sinc}}\left((k-l){\varphi }_{\text{m}}\right)i_{0,k}(\gamma ).\end{array}$$(B.38)

Similarly, we have $$j_l(\gamma )=\frac{{\varphi }_{\text{m}}}{\pi }{\displaystyle \sum _{k=-d}^d\mathrm{sinc}}\left((k-l){\varphi }_{\text{m}}\right)j_{0,k}(\gamma ).$$(B.39)

Since the sinus-cardinal function rapidly decreases, the Fourier expansions of i and j can still be limited to a few harmonics. We note that configurations where the spots are visible for a very short fraction of the period (ϕ_m ≪ π) would require us to include more harmonics since the sinus-cardinal will decay more slowly with k − l. However, such configurations also correspond to spots that always remain close to the limb of the star, which typically lead to a lower amplitude in the RV and indicators, due to projection and limb-darkening effects. The periodic part of the kernel $k_{i,j}^{\text{per}}$ can thus be approximated with a few harmonics, $$k_{i,j}^{\text{per}}(\tau ,\gamma )={\displaystyle \sum _{k=-d}^d{\overline{i}}_k}(\gamma )j_k(\gamma )\exp (\text{i}k\omega \tau ).$$(B.40)

B.3 CCF-related channels

In Appendix B.1, we assumed that the measured RV is the sum of the local stellar RVs weighted by their relative flux. However, this is a simplistic assumption. In the present Section we derive the expressions of RV, as well as ancillary indicators considering that the measured cross correlation function (CCF) is a weighted sum of the local stellar CCFs.

When extracting RV measurements through the CCF technique, one models the CCF with a Gaussian function $$m_{\text{CCF}}(\upsilon ;\mathit{\eta })=a\left(1-c\exp \left(-\frac{1}{2}{\left(\frac{\upsilon -{\upsilon }_0}{\sigma }\right)}^2\right)\right)$$(B.41)

where a is the quiet surface flux, c is the contrast, υ₀ is the RV, $2\sqrt{2\ln (2)}\sigma$ is the FWHM. The parameters η = (a, c, υ₀, σ) are typically adjusted using a least-square estimator. We denote by η₀ the parameters obtained by fitting a CCF which is not affected by the activity contribution. We now consider the contribution δCCF of a small feature (spot/faculae) on the CCF. The impact of this feature on the parameters η can be estimated by linearizing the Gaussian model in the vicinity of η₀ $$m_{\text{CCF}}\left(\upsilon ;{\mathit{\eta }}_0+\delta \mathit{\eta }\right)=m_{\text{CCF}}\left(\upsilon ;{\mathit{\eta }}_0\right)+{\nabla }_{\mathit{\eta }}{m}_{\text{CCF}}\left(\upsilon ;{\mathit{\eta }}_0\right)\delta \mathit{\eta },$$(B.42)

with $$\begin{array}{l}{\nabla }_{\mathit{\eta }}{m}_{CCF}\left(\upsilon ;{\mathit{\eta }}_0\right)=\hfill \\ (1-cG(\upsilon )-aG(\upsilon )-\frac{ac\left(\upsilon -{\upsilon }_0\right)}{{\sigma }^2}G(v)\hfill \end{array}-\frac{ac{\left(\upsilon -{\upsilon }_0\right)}^2}{{\sigma }^3}G(\upsilon )).$$(B.43)

and $G(\upsilon )=\exp \left(-\frac{1}{2}{\left(\frac{\upsilon -{\upsilon }_0}{\sigma }\right)}^2\right)$. The parameters are then obtained by finding the value of δη minimizing $${\displaystyle {\int }_{{\upsilon }_{\min }}^{{\upsilon }_{\max }}{(\delta \text{CCF}(\upsilon )-\nabla m\cdot \delta \mathit{\eta })}^2}\text{\hspace{0.17em}d}\upsilon ,$$(B.44)

where the CCF is computed on a velocity interval [υ_min, υ_max]. This yields the parameters estimate $$\delta \mathit{\eta }={\mathit{\alpha }}^{-1}\mathit{\beta },$$(B.45)

with α the 4 × 4 matrix, and β the vector of size four given by $$\begin{array}{l}\mathit{\alpha }={\displaystyle {\int }_{{\upsilon }_{\min }}^{{\upsilon }_{\max }}\nabla }m{(}^{\upsilon }\nabla m(\upsilon )\text{d}v,\\ \mathit{\beta }={\displaystyle {\int }_{{\upsilon }_{\min }}^{{\upsilon }_{\max }}\nabla }m{(}^{\upsilon }\delta \text{CCF}(\upsilon )\text{d}\upsilon .\end{array}$$(B.46)

For on-ground spectroscopy, the quiet surface flux a is actually strongly affected by the Earth’s atmosphere and other observational circumstances, which makes this parameter useless for stellar characterization. Moreover, the correlation between the quiet surface flux a, and the other parameters can be shown to be negligible, such that we can actually ignore this parameter and approximate the variations of the three other parameters η_1: = (c, υ₀, σ) with $$\delta {\mathit{\eta }}_{1:}={\mathit{\alpha }}_{1:,1:}^{-1}{\mathit{\beta }}_{1:}.$$(B.47)

Indeed, for a sufficiently large interval [υ_min, υ_max], the quiet surface flux a is dominated by the information located on the tails of the CCF, while all other parameters use the information contained in the center of the CCF (see Eq. (B.43)). We can verify this by increasing the length of the interval [vmin, vmax] toward ] − ∞, +∞[. As we do so, all the components of the matrix α converge, except for α_0,0 (which corresponds to the quiet surface flux a), which is asymptotically proportional to the interval length (Δυ = υ_max − υ_min). We can then deduce the asymptotic behavior of the covariance matrix of the parameters α⁻¹, as well as the correlation matrix. We find that the correlation between a and the other parameters is asymptotically proportional to $1/\sqrt{\mathrm{\Delta }\upsilon }$, and thus vanishes for a sufficiently large interval.

In the limit of [υ_min, υ_max] →] − ∞, ∞[, we obtain $${\mathit{\alpha }}_{1:,1:}=a^2\sqrt{\pi }\sigma \left(\begin{array}{ccc}1& 0& \frac{c}{2\sigma }\\ 0& \frac{c^2}{2{\sigma }^2}& 0\\ \frac{c}{2\sigma }& 0& \frac{3c^2}{4{\sigma }^2}\end{array}\right),$$(B.48)

therefore we have $${\mathit{\alpha }}_{1:,1:}^{-1}=\frac{2\sigma }{a^2{c}^2\sqrt{\pi }}\left(\begin{array}{ccc}\frac{3c^2}{4{\sigma }^2}& 0& -\frac{c}{2\sigma }\\ 0& 1& 0\\ -\frac{c}{2\sigma }& 0& 1\end{array}\right),$$(B.49)

and finally $$\left(\begin{array}{c}\frac{\delta c}{\epsilon }\\ \frac{\delta \upsilon }{\sigma }\\ \frac{\delta \sigma }{\sigma }\end{array}\right)=-\frac{1}{\sqrt{2}ac}\left(\begin{array}{ccc}3& 0& -2\\ 0& 4& 0\\ -2& 0& 4\end{array}\right)\mathit{b},$$(B.50)

with b given by (k = 0, 1, 2) $$b_k=\frac{1}{\sigma \sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{+\infty }{\left(\frac{\upsilon -{\upsilon }_0}{\sigma }\right)}^k}G(\upsilon )\delta \text{CCF}(\upsilon )\text{d}\upsilon .$$(B.51)

We now need to specify the effect δCCF of a spot on the CCF to be able to compute b. We assume that in the absence of the spot, the contribution to the CCF of the stellar surface at the position of the spot would have been $${\text{CCF}}_l(\upsilon )=a_l\left(1-c_l\exp \left(-\frac{1}{2}{\left(\frac{\upsilon -{\upsilon }_l}{{\sigma }_l}\right)}^2\right)\right).$$(B.52)

Because of the presence of the spot, the local CCF is altered. As a first approximation, we assume that it remains approximately Gaussian but with altered parameters a_s , c_s , υ_s , σ_s : $${\text{CCF}}_s(\upsilon )=a_s\left(1-c_s\exp \left(-\frac{1}{2}{\left(\frac{\upsilon -{\upsilon }_s}{{\sigma }_s}\right)}^2\right)\right),$$(B.53)

where the flux a_s < a_l for spots and a_s > a_l for faculae. Due to the inhibition of the convective blue-shift, which results in a net red-shift, we expect υ_s > υ_l for both spots and faculae. Moreover, the Zeeman broadening effect in both kinds of active regions tends to decrease the contrast c_s < c_l, and to increase the FWHM (σ_s > σ_l). Finally, the effect of the spot on the CCF can be modeled as $$\delta \mathrm{CCF}(\upsilon )={\mathrm{CCF}}_s(\upsilon )-{\text{CCF}}_l(\upsilon ).$$(B.54)

Replacing this expression in Eq. (B.51), we find $$b=\left(a_s-a_l\right)\left(\begin{array}{l}1\hfill \\ 0\hfill \\ 1\hfill \end{array}\right)$$(B.55) $$\begin{array}{l}+a_l\frac{c_l{\sigma }_l}{\sqrt{{\sigma }^2+{\sigma }_l^2}}\exp \left(-\frac{1}{2}\frac{\delta {\upsilon }_l^2}{{\sigma }^2+{\sigma }_l^2}\right)\left(\begin{array}{c}\frac{\sigma }{{\sigma }^2+{\sigma }_l^2}\delta {\upsilon }_l\\ \frac{{\sigma }_l^2}{{\sigma }^2+{\sigma }_l^2}+\frac{{\sigma }^2}{{\left({\sigma }^2+{\sigma }_l^2\right)}^2}\delta {\upsilon }_l^2\end{array}\right)\\ -a_s\frac{c_s{\sigma }_s}{\sqrt{{\sigma }^2+{\sigma }_s^2}}\exp \left(-\frac{1}{2}\frac{\delta {\upsilon }_s^2}{{\sigma }^2+{\sigma }_s^2}\right)\left(\begin{array}{c}\frac{\sigma }{{\sigma }^2+{\sigma }_s^2}\delta {\upsilon }_s\\ \frac{{\sigma }_s^2}{{\sigma }^2+{\sigma }_s^2}+\frac{{\sigma }^2}{{\left({\sigma }^2+{\sigma }_s^2\right)}^2}\delta {\upsilon }_s^2\end{array}\right),\end{array}$$

with δυ_l = υ_l − υ₀, δυ_s = υ_s − υ₀. At leading order in δv, we find $$\begin{array}{l}{b}_0\approx a_s\left(1-\frac{c_s{\sigma }_s}{\sqrt{{\sigma }^2+{\sigma }_s^2}}\right)-a_l\left(1-\frac{c_l{\sigma }_l}{\sqrt{{\sigma }^2+{\sigma }_l^2}}\right),\\ b_1\approx \frac{a_l{c}_l{\sigma }_l\sigma }{{\left({\sigma }^2+{\sigma }_l^2\right)}^{3/2}}\delta {\upsilon }_l-\frac{a_s{c}_s{\sigma }_s\sigma }{{\left({\sigma }^2+{\sigma }_s^2\right)}^{3/2}}\delta {\upsilon }_s,\\ b_2\approx a_s\left(1-\frac{c_s{\sigma }_s^3}{{\left({\sigma }^2+{\sigma }_s^2\right)}^{3/2}}\right)-a_l\left(1-\frac{c_l{\sigma }_l^3}{{\left({\sigma }^2+{\sigma }_l^2\right)}^{3/2}}\right),\end{array}$$(B.56)

and thus (see Eq. (B.50)) $$\begin{array}{l}\sqrt{2}a\delta c\approx a_l\left(1-c_l\frac{{\sigma }_l\left(3{\sigma }^2+{\sigma }_l^2\right)}{{\left({\sigma }^2+{\sigma }_l^2\right)}^{3/2}}\right)-a_s\left(1-c_s\frac{{\sigma }_s\left(3{\sigma }^2+{\sigma }_s^2\right)}{{\left({\sigma }^2+{\sigma }_s^2\right)}^{3/2}}\right),\\ \frac{ac}{2\sqrt{2}{\sigma }^2}\delta \upsilon \approx \frac{a_s{c}_s{\sigma }_s}{{\left({\sigma }^2+{\sigma }_s^2\right)}^{3/2}}\delta {\upsilon }_s-\frac{a_l{c}_l{\sigma }_l}{{\left({\sigma }^2+{\sigma }_l^2\right)}^{3/2}}\delta {\upsilon }_l,\\ \frac{ac}{\sqrt{2}\sigma }\delta \sigma \approx a_l\left(1+c_l\frac{{\sigma }_l\left({\sigma }^2-{\sigma }_l^2\right)}{{\left({\sigma }^2+{\sigma }_l^2\right)}^{3/2}}\right)-a_s\left(1+c_s\frac{{\sigma }_s\left({\sigma }^2-{\sigma }_s^2\right)}{{\left({\sigma }^2+{\sigma }_s^2\right)}^{3/2}}\right).\end{array}$$(B.57)

We thus find that, in addition to the velocity variations induced by the photometric effect (proportional to a_s − a_l) and the convective blue-shift inhibition effect (υ_s − υ_l), changes in the equivalent width c_sσ_s − c_lσ_l might also affect the velocity. Moreover, we observe that at leading order the contrast variations δc and the FWHM variations δσ are mainly proportional (with negative signs) to the flux variations due to the spot (a_s − a_l), with some corrections coming from the shape variations (c_s − c_l, σ_s − σ_l). These indicators might thus be modeled in first approximation with the photometric component of our model if we neglect the correction terms for the shape variations.

Appendix C A computationally efficient GP framework

C.1 Overview of the mean and covariance calculation

To compute the mean function and covariance in an analytical way, be it for the physics-based or data-driven model. In Section 3.4, we made several assumptions reproduced below for convenience

1. The effect of a single feature is represented as a product of a time window W(t) and a periodic part, decomposed in Fourier series. Denoting by γ = (γ_W, γ_per, ϕ₀), $$g_i(t,\mathit{\gamma })=W\left(t,{\mathit{\gamma }}_W\right)P_i\left(t,{\mathit{\gamma }}_{\text{per}},{\varphi }_0\right)$$(C.1) $$P_i\left(t,{\mathit{\gamma }}_{\text{per\hspace{0.17em}}},{\varphi }_0\right)={\displaystyle \sum _{k=-d}^d{i}_k}\left({\mathit{\gamma }}_{\text{per\hspace{0.17em}}}\right){\text{e}}^{\text{i}k\left(\omega t+{\varphi }_0\right)}.$$(C.2)

2. The phase φ₀ is distributed uniformly on [0,2π].

3. The parameters of the window and of the periodic parts are supposed to be statistically independent, $$p(\mathit{\gamma }\mid \mathit{\eta })=p\left({\mathit{\gamma }}_W\mid \mathit{\eta }\right)p\left({\mathit{\gamma }}_{\text{per\hspace{0.17em}}}\mid \mathit{\eta }\right).$$(C.3)

4. We assume that the properties of the distribution are independent of time, p(γ|η, t) = p(γ|η).

5. We assume that λ(t) is a stationary random process with expectation µ^λ(η), and covariance function k^λ(τ; η)

6. The typical timescale of the variation in λ(t) is much greater than the rotation period (2π/ω in Eq. (C.2)).

Based on assumptions 1-6, we can compute the mean and covariance of the process. First, we introduce the following notations for the kernel associated with the window and periodic parts, $$k^W\left(\tau ;{\mathit{\gamma }}_W\right)={\displaystyle \int W}\left(t,{\mathit{\gamma }}_W\right)W\left(t+\tau ,{\mathit{\gamma }}_W\right)\text{d}t$$(C.4) $$k^W(\tau ;\mathit{\eta })={\displaystyle \int k^W}\left(\tau ;{\mathit{\gamma }}_W\right)p\left({\mathit{\gamma }}_W\mid \mathit{\eta }\right)\text{d}{\mathit{\gamma }}_W$$(C.5) $$k_{i,j}^{\text{per}}\left(\tau ;{\mathit{\gamma }}_{\text{per}}\right)={\displaystyle \int P_i}\left(t,{\mathit{\gamma }}_{\text{per}}\right)P_j\left(t+\tau ,{\mathit{\gamma }}_{\text{per}}\right)\text{d}t$$(C.6) $$k_{i,j}^{\text{per}}(\tau ;\mathit{\eta })={\displaystyle \int k_{i,j}^{\text{per}}}\left(\tau ;{\mathit{\gamma }}_{\text{per}}\right)p\left({\mathit{\gamma }}_{\text{per}}\mid \mathit{\eta }\right)\text{d}{\mathit{\gamma }}_{\text{per}}.$$(C.7)

We introduce the average window, its mean and covariance $${\mu }^W(t,\mathit{\eta })={\displaystyle \int W}\left(t,{\mathit{\gamma }}_W\right)p\left({\mathit{\gamma }}_W\mid \mathit{\eta }\right)\text{d}\mathit{\gamma },$$(C.8) $${\mu }^W(\mathit{\eta })={\displaystyle \int {\mu }^W}(t,\mathit{\eta })\text{d}t,$$(C.9) $$k^{{\mu }^W}(\tau ;\mathit{\eta })={\displaystyle \int {\mu }^W}(t,\mathit{\eta }){\mu }^W(t+\tau ,\mathit{\eta })\text{d}t.$$(C.10)

To establish the final form of the kernels, we proceed by step. We first suppose that λ is constant. Thanks to assumption 1, 2 and 4, as shown in Appendix C.2, Eq. (11) becomes $$k_{i,j}(\tau ;\mathit{\eta })=\lambda {\displaystyle \int k^W}(\tau ,\mathit{\gamma })k_{i,j}^{\text{per}}(\tau ,\mathit{\gamma })p(\mathit{\gamma }\mid \mathit{\eta })\text{d}\mathit{\gamma }$$(C.11)

And thanks to assumption 3, we can further simplify Eq. (C.11) $$k_{i,j}(\tau ;\mathit{\eta })=\lambda k^W(\tau ;\mathit{\eta })k_{i,j}^{\text{per}}(\tau ;\mathit{\eta }).$$(C.12)

In Appendix C.3, we show that the covariance matrix associated with the kernel $k_{i,j}^{\text{per}}$ can be put in a semi-separable form, in the particular case where the parameters of the features γ, are fixed. In Appendix C.4, we do the same for the window part. In Appendix C.5, we show that the covariance matrix is still semi-separable when in the case where γ are drawn from the distribution of feature parameters.

If we now take λ(t) as a stationary GP as in assumption 5, as shown in Appendix C.6, the covariance is $$\begin{array}{l}{k}_{i,j}(\tau ;\mathit{\eta })={\mu }^{\lambda }{k}^W(\tau ;\mathit{\eta })k_{i,j}^{\text{per\hspace{0.17em}}}(\tau ;\mathit{\eta })\hfill \\ +{\mu }_i(\mathit{\eta }){\mu }_j(\mathit{\eta })\left(k^{\lambda }\ast k^{{\mu }^W}\right)(\tau ;\mathit{\eta }),\hfill \end{array}$$(C.13)

where µ_i(η) = E_γ(i₀(γ)), µ_j(η) = E_γ(j₀(γ)), i₀ and j₀ are the zero order Fourier coefficient of the periodic part, as defined in Eq. (C.2). Since the timescale of the magnetic cycle is much longer than the magnetic region lifetime (assumption 6), we make the approximation $\left(k^{\lambda }*k^{{\mu }^W}\right)(\tau ;\mathit{\eta })\propto k^{\lambda }(\tau ;\mathit{\eta })$. Given that the $k^{\lambda }(\tau ;\mathit{\eta })$ kernel has already an amplitude term degenerate with the proportionality coefficient, we simply write $$k_{i,j}\left(\tau ;\eta \right)={\mu }^{\lambda }{k}^W\left(\tau ;\eta \right)k_{i,j}^{\text{per}}\left(\tau ;\eta \right)+{\mu }_i\left(\eta \right){\mu }_j\left(\eta \right)k^{\lambda }\left(\tau ;\eta \right)$$(C.14)

The first term in Eq. (C.13) is exactly the same as in the case of a constant rate λ = µ^λ (Eq. (C.12)) but we have an additional term capturing the variations of λ. Similarly, the mean function of the process is $$f_i(\eta )={\mu }^{\lambda }{\mu }^W(\eta ){\mu }_i(\eta )$$(C.15)

Finally, in Appendix C.7, we present in detail the data-driven kernels.

C.2 Separation of the window and the periodic components

To simplify the notations, in the reminder of Appendix C, we consider two channels y, z, with associated impulse responses I(t, γ), J(t, γ). We consider here the stationary case, where the rate of spot appearance λ is constant, and the spots properties γ do not depend on time p(γ|η, t) = p(γ|η). The covariance between two time series y(t) = ∑_k I(t − t_k, γ_k) and z(t) = ∑_k J(t − t_k, γ_k), is $$k^{y,z}(\tau ;\eta )=\lambda {{\displaystyle \iint }}^{\text{}}I(t,\gamma )J(t+\tau ,\gamma )\text{d}tp(\gamma \mid \eta )\text{d}\gamma ,$$(C.16)

where η is the set of hyper-parameters describing the spots population properties, I(t, γ) = W(t, γ)i(t, γ) and J(t, γ) = W(t, γ)j(t, γ). Up to now, we did not specify how the time of appearance of the magnetic regions, t_k is defined. In the following we define it as the time of maximum "activity" of the spot. This means that W reaches its maximum at t − t_k = 0. We assume that the longitude ϕ₀ of the spots at t = t_k (which is one of the parameters γ) is uniformly distributed in [0,2π] and independent of t_k and of the other spots properties. We thus take this parameter out of γ, and rewrite I, J as I(t, ϕ₀, γ) = W(t, γ)i(ϕ(t), γ) and J(t, ϕ₀, γ) = W(t, γ) j(ϕ(t), γ), with ϕ(t) = ϕ₀ + ωt, such that $$k^{y,z}(\tau ;\eta )=\lambda {{\displaystyle \iint }}^{\text{}}I\left(t,{\varphi }_0,\gamma \right)J\left(t+\tau ,{\varphi }_0,\gamma \right)\text{d}t\frac{\text{d}{\varphi }_0}{2\pi }p(\gamma \mid \eta )\text{d}\gamma .$$(C.17)

We first aim at computing $$k^{I,J}\left(\tau ,{\varphi }_0,\gamma \right)={\displaystyle {\int }_{-\infty }^{+\infty }I\left(t,{\varphi }_0,\gamma \right)}J\left(t+\tau ,{\varphi }_0,\gamma \right)\text{d}t.$$(C.18)

In the case I = J, the Fourier transform of k^1,I is the power spectral density of I (Wiener–Khinchin theorem) $${\widehat{k}}^{I,I}=|\widehat{I}{|}^2.$$(C.19)

In the more general case, $${\widehat{k}}^{I,J}=\overline{\widehat{I}}\widehat{J}$$(C.20)

Moreover, since I is the product of W and i, its Fourier transform is written as a product of convolution $$\widehat{I}\left(v,{\varphi }_0,\gamma \right)={\displaystyle \sum }_{k\in \mathbb{Z}}\widehat{W}(v-k\omega ,\gamma )i_k(\gamma )\exp \left(\text{i}k{\varphi }_0\right),$$(C.21)

where i_k is the Fourier coefficient $$i_k(\gamma )=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }i(\varphi ,\gamma )}\exp (-\text{i}k\varphi )\text{d}\varphi .$$(C.22)

The same is true for J and the Fourier transform of k^I,J is thus written as $$\begin{array}{cc}{\widehat{k}}^{l,J}\left(v,{\varphi }_0,\gamma \right)=& {\displaystyle \sum }_{k,l\in \mathbb{Z}}\overline{\widehat{W}}(v-k\omega ,\gamma )\widehat{W}(v-l\omega ,\gamma )\\ & \times {\overline{i}}_k(\gamma )j_l(\gamma )\exp \left(\text{i}(l-k){\varphi }_0\right).\end{array}$$(C.23)

We now integrate this over the longitude at maximum activity ϕ₀. This longitude only appears in the term exp(i(l − k)ϕ₀), and we have $${\displaystyle {\int }_0^{2\pi }\exp }\left(\text{i}(l-k){\varphi }_0\right)\frac{\text{d}{\varphi }_0}{2\pi }={\delta }_{k,l},$$(C.24)

thus the integral of ${\widehat{k}}^{I,J}$ simplifies to $${\displaystyle {\int }_0^{2\pi }{\widehat{k}}^{I,J}}\left(v,{\varphi }_0,\gamma \right)\frac{\text{d}{\varphi }_0}{2\pi }={\displaystyle \sum _{k\in \mathbb{Z}}|}\widehat{W}{|}^2(v-k\omega ,\gamma ){\overline{i}}_k(\gamma )j_k(\gamma ).$$(C.25)

Finally k^y,z simplifies to $$k^{y,z}(\tau ;\eta )=\lambda {\displaystyle \int k^W}(\tau ,\gamma )k_{i,j}^{\text{per}}(\tau ,\gamma )p(\gamma \mid \eta )\text{d}\gamma ,$$(C.26)

with $$\begin{array}{cc}{k}^W(\tau ,\gamma )& ={\displaystyle {\int }_{-\infty }^{+\infty }W(t,\gamma )}W(t+\tau ,\gamma )\text{d}t\\ k_{i,j}^{\text{per}}(\tau ,\gamma )& =\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }i(\varphi ,\gamma )}j(\varphi +\omega \tau ,\gamma )\text{d}\varphi \\ & ={\displaystyle \sum }_{k\in \mathbb{Z}}{\overline{i}}_k(\gamma )j_k(\gamma )\exp (\text{i}k\omega \tau ).\end{array}$$(C.27)

C.3 Efficient modeling with S+LEAF

The kernel of Eq. (B.40) can be very efficiently modeled using S+LEAF (Delisle et al. 2020, 2022) or celerite (Foreman-Mackey et al. 2017). Indeed, in the general case, the cost of likelihood evaluations of a GP is proportional to the number of measurements cubed. For a restricted class of kernel functions, the S+LEAF and celerite GP frameworks allow one to significantly reduce this cost, to a linear scaling with the number of points. In order to be able to use these frameworks, the covariance between the measurements Y = (Y₁,…, Y_n) must be semiseparable, that is it can be written: $$\mathrm{cov}\left(Y_a,Y_b\right)={\displaystyle \sum }_{s=1}^r{U}_{a,s}{V}_{b,s},\text{\hspace{0.17em}for\hspace{0.17em}\hspace{0.17em}}a>b,$$(C.28)

where r is the rank of the semiseparable representation. In particular, sums and products of exponential kernels, sinusoidal kernels, and Matérn kernels are semiseparable. As shown in Delisle et al. (2022), this reasoning remains valid when considering several time series affected by the same GP but with different coefficients. In this case, the different time series should be merged in a single heterogeneous time series Y, ordered by increasing time, and whose covariance should be semiseparable.

We consider q time series y⁽¹⁾(t),…, y^(q)(t), which we merge and sort by increasing time in a single heterogeneous time series Y. Considering only the periodic part k^i,j for a fixed set of spot parameters γ, and ignoring the window part k^W for now, the covariance matrix of the merged time series is $$\mathrm{cov}\left(Y_a,Y_b\right)={\displaystyle \sum }_{k=-d}^d\left(i_{I_a,k}\exp \left(\text{i}k\omega t_a\right)\right)\left({\overline{i}}_{I_b,k}\exp \left(-\text{i}k\omega t_b\right)\right),$$(C.29)

for a > b and where I_k is the index associating with each point in the merged time series Y the identifier of the original time series it comes from and i_s,k is the k-th Fourier coefficient of the time series s. The covariance of Eq. (C.29) is thus semiseparable with rank r = 2d + 1. However, to obtain Eq. (C.29), we neglected the window part k^W, and more importantly, we assumed fixed parameters γ for the spots.

C.4 Window function

Up to now we did not specify the shape W of the spots’ appearing and disappearing process, which models the evolution of the combined effect of area and temperature of the feature as a function of time. In Section 3.1, we mentioned that on the Sun, spots and faculae tend to appear faster than they decay. To have a S+LEAF representation of the noise, in the present work the window function is such that both its increase and decrease in intensity are exponential, potentially with different time-scales.

We provide in Table C.1 the kernel k^W obtained when assuming different window shapes W : sudden spot appearing and exponential decay, symmetric or asymmetric exponential appearing and disappearing. In these three cases, the obtained kernel is semiseparable with low rank (see Foreman-Mackey et al. 2017; Delisle et al. 2022). More generally, the Matérn kernel family can be used for efficient modeling of k^W with S+LEAF or celerite.

Finally, for fixed spot properties γ, the covariance matrix of the merged time series Y, including the window contribution, is semiseparable with rank r = r_Wr_per.. = r_w(2d + 1) since the kernel is the product of two semiseparable kernels (see Eq. (C.26)).

C.5 Distribution of spot properties

We now aim at integrating our kernel over the distribution of spot parameters γ. We assume that the spot parameters γ can be split in two sets of independent parameters, γ_W (only affecting W) and γ_per. only affecting the periodic part. By independent we mean $$p\left({\gamma }_W,{\gamma }_{\text{per.}}\mid \eta \right)=p\left({\gamma }_W\mid \eta \right)p\left({\gamma }_{\text{per}}\mid \eta \right),$$(C.30)

such that the integral of Eq. (C.26) can be written as $$k^{y,z}(\tau ;\eta )=\lambda k^W(\tau ;\eta )k^{i,j}(\tau ;\eta )$$(C.31)

with $$\begin{array}{cc}& k^W(\tau ;\eta )={{\displaystyle \int }}^{\text{}}{k}^W\left(\tau ,{\gamma }_W\right)p\left({\gamma }_W\mid \eta \right)\text{d}{\gamma }_W,\\ & k_{i,j}^{\text{per}}(\tau ;\eta )={{\displaystyle \int }}^{\text{}}{k}_{i,j}^{\text{per}}\left(\tau ,{\gamma }_{\text{per.}}\right)p\left({\gamma }_{\text{per.}}\mid \eta \right)\text{d}{\gamma }_{\text{per.}}.\end{array}$$(C.32)

where k^W(τ, γ_W) and $k_{i,j}^{\text{per}}\left(\tau ,{\gamma }_{\text{per.}}\right)$ are given in Eqs. (C.4) and (C.6), respectively.

In the following, we assume that the integral of the window kernel k^W(τ; η) can still be approximated with a Matérn kernel, such that it remains semiseparable. For the periodic part, $k_{i,j}^{\text{per}}(\tau ;\eta )$ can be decomposed in a Fourier series, and its k-th coefficient is $$C_k^{i,j}(\eta )={{\displaystyle \int }}^{\text{}}{\overline{i}}_k(\gamma )j_k(\gamma )p(\gamma \mid \eta )\text{d}\gamma .$$(C.33)

These integral might be performed analytically for some specific choices of distributions p(γ|η). However, this is beyond the scope of this article and we rather turn here to numerical integrals. We thus need to perform an integral for each harmonics k and each couple i, j of original time series (RV and indicators). For each harmonics k we obtain a Hermitian positive definite matrix $C_k={\left(C_k^{i,j}\right)}_{i=\mathrm{1..}q,j=\mathrm{1..}q}$ of size q × q of these integrated coefficients, where q is the number of original time series. We then compute a square root of the matrix C_k (e.g., by Cholesky decomposition or eigendecomposition) $$C_k=R_k{R}_k^{\ast },$$(C.34)

such that the coefficients, $C_k^{i,j}$, are given by $$C_k^{i,j}(\eta )={\displaystyle \sum }_{s=1}^q{R}_{k,i,s}{\overline{R}}_{k,j,s}.$$(C.35)

With such a decomposition, we obtain a semiseparable representation for $k_{i,j}^{\text{per}}$ with rank q(2d + 1), where q is the number of considered time series (RV and indicators) and d is the harmonics at which the Fourier series is truncated. Finally, the full covariance matrix – including the window part k^W and integrated over the distribution of spots parameters – is semiseparable with rank r = qr_W(2d + 1).

We note that the integrals over the distribution of spot parameters required to compute the matrices C_k and R_k might present a significant computational cost. For better performances, we precompute the matrices R_k on a grid of values for η, and then interpolate on this grid to get an estimate of R_k for a given η.

C.6 Variable rate λ

We now consider the case of a variable spot appearance rate λ(t). We assume that λ(t) is a stationary random process with expectation µ^λ(η), and covariance function k^λ(τ; η). Following previous sections, we obtain the following conditional expectations (given λ): $$\begin{array}{l}\mathbb{E}\left\{y\left(t_i\right)\mid \lambda (t),\eta \right\}=\underset{-\infty }{\overset{+\infty }{{\displaystyle \int }}}\lambda (t){\mathbb{E}}_{\gamma }\left\{I\left(t_i-t,\gamma \right)\right\}\text{d}t,\\ \begin{array}{cc}\mathbb{E}\left\{y\left(t_i\right)z\left(t_j\right)\mid \lambda (t),\eta \right\}& =\underset{-\infty }{\overset{+\infty }{{\displaystyle \int }}}\lambda (t){\mathbb{E}}_{\gamma }\left\{I\left(t_i-t,\gamma \right)J\left(t_j-t,\gamma \right)\right\}\text{d}t\\ & +\mathbb{E}\left\{y\left(t_i\right)\mid \lambda (t),\eta \right\}\mathbb{E}\left\{z\left(t_j\right)\mid \lambda (t),\eta \right\}.\end{array}\end{array}$$(C.36)

Moreover, the reasoning of Sects. C.2–C.5 holds and we have $$\begin{array}{ll}\hfill & {\mathbb{E}}_{\gamma }\left\{I\left(t_i-t,\gamma \right)\right\}={\mathbb{E}}_{\gamma }\left\{W\left(t_i-t,\gamma \right)\right\}{\mathbb{E}}_{\gamma }\left\{i_0(\gamma )\right\},\hfill \\ \hfill & {\mathbb{E}}_{\gamma }\left\{I\left(t_i-t,\gamma \right)J\left(t_j-t,\gamma \right)\right\}={\mathbb{E}}_{\gamma }\left\{W\mid \left(t_i-t,\gamma \right)W\left(t_j-t,\gamma \right)\right\}{k}_{i,j}^{per}\left(t_i-t_j;\eta \right)\hfill \end{array}$$(C.37) $$\begin{array}{ll}\hfill & \mathbb{E}\left\{y\left(t_i\right)z\left(t_j\right)\mid \eta \right\}={\mu }^{\lambda }(\eta )k^W\left(t_i-t_j;\eta \right)k_{i,j}^{per}\left(t_i-t_j;\eta \right)\hfill \\ \hfill & +\mathbb{E}\left\{y\left(t_i\right)\mid \eta \right\}\mathbb{E}\left\{z\left(t_j\right)\mid \eta \right\}\hfill \\ \hfill & +{\mathbb{E}}_{\gamma }\left\{i_0(\gamma )\right\}{\mathbb{E}}_{\gamma }\left\{j_0(\gamma )\right\}\hfill \\ \hfill & \times {\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }{k}^{\lambda }}}(t-s){\mathbb{E}}_{\gamma }\left\{W\left(t_i-t,\gamma \right)\right\}{\mathbb{E}}_{\gamma }\left\{W\left(t_j-s,\gamma \right)\right\}\text{d}t\text{d}s.\hfill \end{array}$$(C.38)

The last line is the convolution product of k^λ and the autocorrelation of the mean window function $k^{{\mu }^W}$. This autocorrelation $k^{{\mu }^W}$ differs from the mean autocorrelation k^W , since the mean over the spot properties is taken before computing the autocorrelation.

Finally, the kernel function in the case of a variable spot appearance rate is given by $$\begin{array}{cc}{k}^{y,z}(\tau ;\eta )=& {\mu }^{\lambda }(\eta )k^W(\tau ;\eta )k_{i,j}^{\text{per}}(\tau ;\eta )\\ & +{\mu }_i(\eta ){\mu }_j(\eta )\left(k^{\lambda }\ast k^{{\mu }^W}\right)(\tau ;\eta ),\end{array}$$(C.39)

where µ_i(η) = E_γ(i₀(γ)), µ_j(η) = E_γ(j₀(γ)). The first term in Eq. (C.39) is exactly the same as in the case of a constant rate λ = µ^λ, but we have an additional term capturing the variations of λ.

If we assume that both k^λ and $k^{{\mu }^W}$ can be approximated by Matérn 1/2 kernels, $$k^{\lambda }(\tau )={\sigma }_{\lambda }^2\exp \left(-\left|\tau \right|/{\rho }_{\lambda }\right),$$(C.40) $$k^{{\mu }^w}(\tau )={\sigma }_{{\mu }^w}^2\exp \left(-\left|\tau \right|/{\rho }_{{\mu }^w}\right),$$(C.41)

we obtain $$\begin{array}{ll}\left(k^{\lambda }\ast k^{{\mu }^w}\right)(\tau )=\hfill & 2{\sigma }_{\lambda }^2{\sigma }_{{\mu }^w}^2\frac{{\rho }_{\lambda }{\rho }_{{\mu }^W}^2}{{\rho }_{\lambda }^2-{\rho }_{{\mu }^w}^2}\times \hfill \\ \hfill & \left({\rho }_{\lambda }\exp \left(-\left|\tau \right|/{\rho }_{\lambda }\right)-{\rho }_{{\mu }^w}\exp \left(-\left|\tau \right|/{\rho }_{{\mu }^w}\right)\right).\hfill \end{array}$$(C.42)

In the limit ${\rho }_{{\mu }^W}\ll {\rho }_{\lambda }$, we have $$\left(k^{\lambda }\ast k^{{\mu }^w}\right)(\tau )\approx 2{\sigma }_{\lambda }^2{\sigma }_{{\mu }^w}^2{\rho }_{{\mu }^w}\exp \left(-\left|\tau \right|/{\rho }_{\lambda }\right)$$(C.43)

which is a Matérn 1/2 kernel with timescale ρ_λ and amplitude ${\sigma }_{\ast }={\sigma }_{\lambda }{\sigma }_{{\mu }^W}\sqrt{2{\rho }_{{\mu }^W}}$.

Similarly, if k^λ and $k^{{\mu }^W}$ can be approximated by Gaussian kernels, $$k^{\lambda }(\tau )={\sigma }_{\lambda }^2\exp \left(-{\tau }^2/\left(2{\rho }_{\lambda }^2\right)\right),$$(C.44) $$k^{{\mu }^W}(\tau )={\sigma }_{{\mu }^w}^2\exp \left(-{\tau }^2/\left(2{\rho }_{{\mu }^W}^2\right)\right),$$(C.45)

their convolution is also a Gaussian kernel: $$\begin{array}{l}\left(k^{\lambda }\ast k^{{\mu }^w}\right)(\tau )=\sqrt{2\pi }{\sigma }_{\lambda }^2{\sigma }_{{\mu }^W}^2\frac{{\rho }_{\lambda }{\rho }_{{\mu }^W}}{\sqrt{{\rho }_{\lambda }^2+{\rho }_{{\mu }^W}^2}}\exp \left(-{\tau }^2/\left(2\left({\rho }_{\lambda }^2+{\rho }_{{\mu }^w}^2\right)\right)\right)\\ \approx \sqrt{2\pi }{\sigma }_{\lambda }^2{\sigma }_{{\mu }^w}^2{\rho }_{{\mu }^w}\exp \left(-{\tau }^2/\left(2{\rho }_{\lambda }^2\right)\right).\end{array}$$(C.46)

In the limit ${\rho }_{{\mu }^W}\ll {\rho }_{\lambda }$, the kernel $\left(k^{\lambda }\ast k^{{\mu }^W}\right)$ can typically be approximated by a decay kernel with timescale ρ_λ and with an amplitude scaling as ${\sigma }_{\lambda }{\sigma }_{{\mu }^W}\sqrt{{\rho }_{{\mu }^W}}$.

C.7 Data-driven kernels

Replacing Eq. (C.2) in Eq. (11), we showed in Appendix C.5) that $$k_{i,j}^{\text{per}}(\tau ;\eta )={\displaystyle \sum }_{k=-d}^d{C}_k^{i,j}(\eta ){\text{e}}^{\text{i}k\omega \tau }$$(C.47)

where $$C_k^{i,j}(\eta )={{\displaystyle \int }}^{\text{}}{i}_{-k}\left({\gamma }_{\text{per}}\right)j_k\left({\gamma }_{\text{per}}\right)p\left({\gamma }_{\text{per}}\mid \eta \right)\text{d}{\gamma }_{\text{per}}.$$(C.48)

In the case of a real-valued signal, we show also in Appendix C.5 for a given k, $C_k={\left(C_k^{i,j}\right)}_{i=1\dots q,j=\mathrm{1..}q}$ considered as a q × q matrix has a Hermitian symmetry. As a consequence, it can be written with the Cholesky decomposition, as a product of a lower triangular matrix R_k and its complex transpose $R_k^{\ast }$, $$C_k(\eta )=R_k(\eta )R_k{(\eta )}^{\ast }.$$(C.49)

In practice, we need to compute the likelihood within a MCMC or nested sampling algorithm, so we need to evaluate the integrals (C.48) for millions of different values of η. For the physics-based stellar activity model of Section 3, to speed up these evaluations, we pre-compute the integrals as well as R_k(η) in Eq. (C.49) on a tight grid of η parameters, and interpolate linearly for the values of η in the MCMC chain. Another option consists in using the $R_k^{i,j}$ directly as our hyperparameters η, which we describe below.

To avoid manipulating complex numbers in the numerical implementation, instead of using complex-valued matrices R_k we use as variables their complex and imaginary parts. We denote them by A_k and B_k. Then, from Eq. (C.49), $$C_k=\left(A_k+i{B}_k\right){\left(A_k-i{B}_k\right)}^T.$$(C.50)

Since we focus on real-valued signals, we have $C_{-k}^{i,j}(\eta )=C_k^{i,j}{(\eta )}^{\ast }$ (see Eq. (C.47)). Writing k^per(τ; η) the matrix ${\left(k_{i,j}^{\text{per}}(\tau )\right)}_{i,j}$, we have $$k^{\text{per\hspace{0.17em}}}(\tau ;\eta )={\text{A}}_0(\eta )+{\displaystyle \sum _{k=1}^d\overline{C_k(\eta )}}{\text{e}}^{-ik\omega \tau }+C_k(\eta ){\text{e}}^{\text{i}k\omega \tau }$$(C.51)

where $\overline{C_k(\eta )}$ is the term-by-term conjugate of C_k(η) (not its complex transpose). Injecting Eq. (C.49) in Eq. (C.52) yields $$\begin{array}{l}{k}^{\text{per}}(\tau ;\eta )=A_0(\eta )\\ +{\displaystyle \sum _{k=1}^{d}2}\left(A_k{A}_k^T+B_k{B}_k^T\right)\cos (k\omega \tau )+2\left(A_k{B}_k^T-B_k{A}_k^T\right)\sin (k\omega \tau ).\end{array}$$(C.52)

We define our new hyper parameters η = (α_k, β_k)_k=0..d as $$\begin{array}{cc}{\alpha }_0=A_0,& {\beta }_0=0,\text{\hspace{0.17em}and\hspace{0.17em}for\hspace{0.17em}}k⩾1\end{array}$$(C.53) $$\begin{array}{cc}{\alpha }_k=\sqrt{2}{A}_k,& {\beta }_k=\sqrt{2}{B}_k\end{array}$$(C.54)

and we define the MultiFourierKernel function as $$k_{i,j}^{\text{per}}(\tau ;\eta )={\displaystyle \sum _{k=0}^d{a}_{k;i,j}}(\eta )\cos (k\omega \tau )+b_{k;i,j}(\eta )\sin (k\omega \tau )$$(C.55)

where a_k and b_k are q × q matrices defined as $$a_k={\alpha }_k{\alpha }_k^T+{\beta }_k{\beta }_k^T$$(C.56) $$b_k={\alpha }_k{\beta }_k^T-{\beta }_k{\alpha }_k^T$$(C.57)

Because the R_k matrices are lower-triangular complex matrices obtained from Cholesky decomposition, α and β have a form $${\alpha }_k=\left(\begin{array}{cccc}{\alpha }_k^{11}& 0& \dots & 0\\ {\alpha }_k^{21}& {\alpha }_k^{22}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\alpha }_k^{q1}& {\alpha }_k^{q2}& \dots & {\alpha }_k^{qq}\end{array}\right)$$(C.58) $${\beta }_k=\left(\begin{array}{ccccc}0& 0& \dots & \dots & 0\\ {\beta }_k^{21}& 0& \dots & \dots & 0\\ {\beta }_k^{31}& {\beta }_k^{32}& 0& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\beta }_k^{q1}& {\beta }_k^{q2}& {\beta }_k^{q3}& \dots & 0\end{array}\right),$$(C.59)

Appendix F Additional discussions

F.1 Phase shifts might relate to the ratio of spots to faculae

Assuming a constant limb darkening, as is noted in Aigrain et al. (2012), the RV effect due to the photometric effect (Eq. (13)) is proportional to the derivative of the inhibition of convective blueshift effect (Eq. (14)). We can write that RV impulse response g, as ɡ(t + τ) = ɡ_cb(t) + ɡ′_cb(t)τ and ɡ′_cb(t) = αɡ_ph(t). This is just another way to rewrite the FF′ approximation of Aigrain et al. (2012), and this is valid only if the limb-darkening effect is the same for the RV photometric effect, and the convective blueshift inhibition RV effect. This is incorrect, in particular because plages have a limb-brightening effect (Meunier et al. 2010a). Nonetheless, we lay out the reasoning with the assumption that the limb darkening is constant as a starting point for a more realistic model.

The opposite of the flux impulse response in Eq. (15), as apparent in Fig. 3 is in phase with the inhibition of the convective blueshift. As long as the RV photometric effect is smaller than the inhibition of convective blueshift by a factor 3 – 4, we can make the approximation $$\tau \approx -\frac{\mathrm{\Delta }F{R}^{\star }}{2\left(\mathrm{\Delta }F{V}_{cb}+\mathrm{\Delta }{V}_{cb}(F+\mathrm{\Delta }F)\right)}.$$(F.1)

The absolute value of the phase shift depends on several parameters. Its sign is more robustly defined by the ratio of spots to faculae through the sign of ΔF. In any case, ΔV_cb is positive because redshifts are positive, and ΔF is small compared to F. If ΔF is negative, then τ is positive.

We expect that the log R′_HK effect in RV is approximately proportional to the projected area of the magnetic region. As a result, it behaves approximately proportionally to the flux. If the magnetic region is bright, then RV should be late compared to photometry and if the region is dark, RV is in advance compared to the log R′_HK. In Hara et al. (2022), we showed that in the HARPS-N RV observations of the Sun (Dumusque et al. 2021), RV is consistently in advance of 20 ± 6°, which supports that the RV photometric effect is dominated by spots on the Sun.

F.2 Building relevant indicators

As mentioned in Section 2.4, we could try to find a spectral variability indicator such that the impulse response of a stellar feature h(t, γ) in this indicator is such that the effect of the feature on the signal of interest (photometry or RV), ɡ(t, γ), is proportional to h(t, γ) for all γ. However, in most cases indicators do not exhibit a linear relationship. It is known that RV and photometry, or RV and indicators might exhibit a so-called closed-loop relations (Bonfils et al. 2007; Forveille et al. 2009; Santerne et al. 2015; Lanza et al. 2018; Collier Cameron et al. 2019; Zhao et al. 2024). This means that when one is plotted against the other, the figure seems to close on itself. In Fig. F.1, we show the behavior of RV versus photometry when they are modeled as in Sec 3.1.

Another approach would be to build several indicators, such that the impulse response of stellar features in the signal of interest is a linear combination of the responses in the different channels. In Haywood et al. (2022), it is argued that the unsigned magnetic flux is a good proxy for the RV component due to the inhibition of the convective blueshift. Neglecting the limb-darkening effect, the radial magnetic field projected onto the line of sight is proportional to RV convective blueshift effect. Since this one dominates on the Sun, we expect it to be a good indicator. However, the limb-darkening law might differ for the velocity and the magnetic field might differ, so we do not expect the correlation to be exact.

F.3 Link with Doppler imaging

Our framework could be extended to perform more classical Doppler imaging. Luger et al. (2021a) map linearly the stellar surface brightness developed in spherical harmonics to the spectrum. In Lehmann & Donati (2022), the authors map the time series of LSD onto principal components, which can be seen as a data-driven way to retrieve the spherical harmonics. By defining observation channels as the spectrum projected on a basis corresponding to the spherical harmonics, we could obtain an estimate of the evolution of their coefficient with time.

References

---

All Tables

All Figures