© The Authors 2025

1. Introduction

Recent observations have revealed a diverse landscape for the magnetism of low-mass stars. This diversity is shown in the properties of their large-scale magnetic fields, which can influence not only the small-scale surface magnetic phenomena (spots, flares, faculae, etc.) but also the physical processes that govern the interiors of stars and their environments. For example, magnetic fields are involved in the transport of energy and angular momentum in stellar interiors by constraining the convection processes (Proctor & Weiss 1982). It has also been shown that the magnetic field can play a role in the star-exoplanet interactions by coupling their two fields (Lanza 2009, 2013; Saur et al. 2013; Strugarek et al. 2015) and influencing the activity of the host star (Cuntz et al. 2000; Shkolnik et al. 2005). The stellar wind and flares can also strongly affect the habitability of exoplanets through different processes of evaporation (Vidal-Madjar et al. 2004; Owen 2019). Understanding and characterising this stellar magnetism is thus a crucial open question in astrophysics since the magnetic field influences many phenomena related to stars and their environments.

In the effort to gain a better understanding of the magnetism of stars, due to its proximity, studying the Sun is especially convenient, as it enables highly detailed observations. Shortly after the Zeeman effect revealed the presence of a magnetic field on the solar surface (Hale 1908), Larmor (1919) suggested that it could be explained by the dynamo effect, which is the capacity of a magnetised fluid to maintain and amplify its own magnetic field against Ohmic dissipation (Moffatt 1978; Parker 1993). There are several phenomena involved in the dynamo loop such as the differential rotation (DR), the turbulence in the convective envelope, the Coriolis force, and/or the meridional circulation (Brun & Browning 2017; Charbonneau 2020). Depending on the stellar properties, these phenomena play various roles in the star’s dynamo. A key parameter in the dynamo description is the Rossby number, which measures the effect of the convection compared to the rotation (see Sect. 4.3 for a precise definition). In particular, the dynamo effect helps explain the Sun’s 11-year sunspot cycle and the evolution of its magnetic topology from dipolar configurations at minimum of activity to quadrupolar at maximum, which lead to the inversion of its magnetic poles with a periodicity of 22 years.

The presence of large-scale magnetic fields suggests that the dynamo effect is also at play inside other stars. However, the diversity and complexity of the topologies that have been observed is probably an indicator of different internal properties. Some fields are much stronger than the solar one, while others are predominantly non-axisymmetric and toroidal compared to the Sun, which has a more poloidal and axisymmetric field (Donati & Landstreet 2009). The activity cycles are also extremely varied, with some stars showing cycles much shorter than the solar one, while others seem to have much longer cycles (Baliunas & Vaughan 1985; Jeffers et al. 2023). Saar & Brandenburg (1999) observed that stars are divided into two branches: one for stars with long cycles and slow rotation, called the ‘inactive branch’, and another for stars with faster rotation than the Sun and shorter cycles, called the ‘active branch’. Later, Böhm-Vitense (2007) linked these two regimes to different dynamo behaviours driven by different internal regions of the stars. Indeed, numerical simulations that reproduce similar activity cycles and other magnetic behaviours indicate that the internal structures and differential rotation profiles of low-mass stars are very diverse (Brun et al. 2022), reinforcing the idea that the stellar dynamo behaviour is varied. This makes it clear that the various ingredients of stellar dynamos need to be better characterised through more observational constraints.

As a consequence of stellar dynamos, a number of magnetic phenomena appear on the stellar surface. This contributes to what is called stellar magnetic activity, which encompasses all the phenomena that induce luminosity variations in the signal produced by a star. The phenomenon of interest in this paper is the activity due to the appearance of spots on stellar surfaces (hereafter referred to as the ‘activity’ for ‘spot activity’), which may be caused by the emergence of a magnetic structure at the surface. The emergence of this structure freezes convection locally and prevents heat from rising up, thus locally cooling down certain regions and making them appear darker than the rest of the star’s surface (Solanki 2003). The spots are a good tracer of the magnetic field because they are easy to observe and their properties (number, area, position, lifetime, etc.) are linked to the magnetic cycle. Spots on the Sun have been observed and studied for centuries; however, the knowledge gained from the behaviour of sunspots may not directly apply to other stars. The effectiveness of the solar analogy for more active stars is uncertain, especially considering that they frequently exhibit larger spotted areas with potentially longer lifetimes (Valio 2017; Namekata et al. 2019). Moreover, in the case of fast rotating stars the spots can be localised at the poles of the star, which is in contrast to what is observed on the Sun (Roettenbacher et al. 2016).

A method to study stellar spots involves the use of precise photometric data. Since stars rotate, the dark spots that appear on their surfaces induce photometric modulations that vary in time and amplitude depending on the physical properties of the spots. These variations are present in the light curves (LCs) from missions such as CoRoT (Baglin 2006), Kepler (Borucki et al. 2010), TESS (Ricker et al. 2014), and soon PLATO (Rauer et al. 2014), which provide high-precision photometric data for long periods of time for thousands of stars and thus allow information about the spots of stars with different physical properties to be retrieved and work on a large data sample to be done. In particular, the Kepler mission has collected around 100 000 LCs from main-sequence stars over a period of 4 years.

There are several indicators that can be used to study the activity of stars, such as the Ca II H and K lines (Wilson 1978; Noyes et al. 1984) emitted by the chromosphere or measurement of X-ray luminosity (Vaiana et al. 1981). In photometry, the S_ph index (see for ex. García et al. 2010; Mathur et al. 2014) is an indicator that allows one to follow magnetically induced modulations. More recently, Basri et al. (2022) presented a method to retrieve mean spot lifetimes by using an autocorrelation function. Another way to extract more precise information about spot properties by using LCs is to fit these data in the temporal domain with a physical model by spot modelling (see for ex. Mosser et al. 2009). Unfortunately, this method has proven to be complex to implement due to the numerous degeneracies of the problem. New methods of analysis could help answer the remaining open questions about stellar magnetism. The idea in this paper is to try to avoid some of these degeneracies encountered in the analysis of temporal time series by shifting to the Fourier domain, where it is much easier to identify global trends about spot properties. The better scale separation provided by sorting temporal frequency in the Fourier domain brings out average information about spots (lifetimes, transit time, and coverage surface). For example, the periodic transit of spots can induce rotation period harmonics that are called ‘rotation peaks’.

The impact of activity in the Fourier domain was first studied by Harvey (1985) and Harvey et al. (1993) for the solar case by a heuristic method. Since the different modulations produced by spot granulation (convective cells visible in the photosphere) and super granulation occur on different timescales and can be seen as memory noises (i.e. the autocorrelation of the time series is a decreasing exponential), they can each be fitted by a pseudo-Lorentzian with a different width (see Sect. 2). The spot activity with the longer timescale is located at the lower end of the spectrum.

Later, Aigrain et al. (2004) showed that some characteristics of this pseudo-Lorentzian could be related to properties of the activity. More specific properties have also been studied: Santos et al. (2017) studied the ratio between rotation peaks harmonics and demonstrated that it can be linked to the spot latitude and eventually to the inclination of the star.

However, it is possible to go beyond these approaches and use all the information contained in the power spectra to their full capacity. The aim of this work is to provide a more thorough method and apply it to a large dataset to build a landscape of stellar activity based on spot properties.

We present here a new model that allows for a more appropriate fit to the power spectrum of photometric data, including the rotational peaks, in order to extract global properties of the spots, such as their lifetime and their photometric impact, which is related to the spot area and the spot temperature contrast. We begin this work by presenting the new model built on an analytical description of a spot signature in an LC in Section 2, and this is followed validation of the model based on simulated LCs in Section 3. We then detail our results by applying the method to the Kepler data in Section 4 in order to take advantage of a large dataset that allows us to identify trends and patterns. We then discuss the results in Section 5 and conclude in Section 6.

2. Modelling stellar spots in the Fourier domain

This section describes the analytical model proposed in this paper. This description is inspired by the original study of Harvey (1985,; 1993) and Aigrain et al. (2004), which have previously worked in the Fourier domain to study the activity of the Sun. To describe the shape of the Fourier power spectra density (PSD), they used a pseudo-Lorentzian expression, namely,

$$\begin{array}{c}\hfill \mathrm{PSD}(\nu )=\frac{\mathcal{H}}{1+{(2\pi \tau \nu )}^{\alpha }},\end{array}$$(1)

where ν is the frequency, ℋ is the height of the Lorentzian, τ the timescale of the activity and α the exponent (see Harvey 1985, for the demonstration). Aigrain et al. (2004) have linked certain properties of this pseudo-Lorentzian to solar activity features: they showed that ℋ is related to the level of activity of the Sun, and suggested that τ could be related to the spot evolution timescale (of the order of a few days). This model has been widely used in stellar photometry to fit the contributions of magnetic activity or granulation. It has also been extended to complementary fields such as radial velocity observations or asteroseismology (see for example Michel et al. 2009; Appourchaux et al. 2012). However, this depiction overlooks certain properties of the spectrum, such as rotational peaks, suggesting that the model can be further improved, as proposed in this section.

We start by discussing the signature of a single spot case in photometry in the temporal domain and then in the Fourier domain in Sect. 2.1. Then we present the prescriptions used for the various functions describing the spot’s evolution in Sect. 2.2. We present the general case with an ensemble of N spots in Sect. 2.3, and finally summarise the main parameters of our model in Sect. 2.4. For the rest of the paper, we use the following notation for the Fourier transform:

$$\begin{array}{c}\hfill TF[f(t)](\nu )={\int }_{-\infty }^{+\infty }f(t){\mathrm{e}}^{2\mathrm{i}\pi \nu \mathrm{t}}\mathrm{d}\mathrm{t}=\widehat{f}(\nu )=\widehat{f(t)}.\end{array}$$

2.1. Single-spot case

We consider a single stellar spot transiting a star with a rotation period P_rot, characterised by an area S_spot and a temperature T_spot. This spot is assumed to always rotate with the same rotation period and to remain at the same latitude. In this model, we do not consider the impact of differential rotation or faculae. The scope of this hypothesis is discussed later in Section 5.4.

The transit of the spot on the stellar disc can be described using three functions, which are illustrated in Fig. 1:

• Ш_Prot: a normalised Dirac comb with a spacing of one rotation period, which represents the periodicity of the transit.

• F_transit(t): a function to describe the mean shape of the transit. The shape of the transit can vary depending on the star’s rotation axis relative to the line of sight, the position of the spot on the disc and limb-darkening effects. By considering that the transit time is defined by the width at mid-height, noted τ_tr, its value is close to P_rot/4, modulated by a factor depending on these configurations. We consider here an average transit shape for the sake of simplicity.

• E_life(t): a function to describe the intrinsic temporal evolution of the spot, with a characteristic time τ_life, which represents the spot lifetime. This function reflects the possible variations of the size and temperature of the spot during its lifetime.

By combining these three components, the impact of the spot on the LC with a unitary maximum depth, noted S₀(t), can be reproduced. The convolution of Ш_Prot and F_transit(t) replicates the spot transits. Then, by multiplying by E_life(t), we add the modulations induced by the intrinsic evolution of the spot. For the LC, this yields

$$\begin{array}{c}\hfill S_0(t)=[{\mathrm{Ш}}_{\mathrm{Prot}}\ast F_{\mathrm{transit}}(t)]\times E_{\mathrm{life}}(t).\end{array}$$(2)

Fig. 1. Illustration of the analytical description of the spot transit proposed in this paper. The top panel shows the model for the single-spot case. The model is shown in the time domain on the left, and in the Fourier domain on the right. Each time, the distinction is made between the single-transit-spot/peakless case (orange) and the multiple-transit-spot/with peaks case (blue). The bottom of the figure illustrates the general case with both multiple-transit and single-transit spots’ contributions.

The value of the maximum transit photometric depth is given by A(T_spot, S_spot) which depends on the temperature contrast between the spot and the stellar surface, and the area of the spot. Therefore, the total impact of a spot on an LC can be expressed as S(t) = A(T_spot, S_spot)⋅S₀(t). This model is illustrated in the top panel of Fig. 1.

Two types of spectrum can be distinguished:

• Stars where τ_life > P_rot: multiple transits of the spot are observed in the LC; it is then categorised as the multiple-transit-spot case.

• Stars where τ_life < P_rot: only one transit of the spot is visible in the LC; it is then the single-transit-spot case.

When shifting to the Fourier domain, we obtain the following:

$$\begin{array}{c}\hfill \widehat{S_0}(\nu )=[\frac{{\mathrm{Ш}}_{1/\mathrm{Prot}}}{{\mathrm{P}}_{\mathrm{rot}}}\times \widehat{F_{\mathrm{transit}}}]\ast \widehat{E_{\mathrm{life}}}(\nu ).\end{array}$$(3)

The convolution between the Fourier comb Ш_1/Prot and the Fourier transform of E_life generates the rotation peaks that are seen in the LC spectrum (modulations of the black line in the bottom panel of Fig. 1). These peaks are multiplied by the Fourier transform of the transit component, $\widehat{F_{\mathrm{transit}}}$, giving the final shape of the spectrum. We can then link the spot properties in the temporal domain to the spectrum properties in the Fourier domain, as summarised in Fig. 1. The information about the spot lifetime τ_life appears thus in the rotation peaks. Depending on the spot lifetime, the transit component and the rotation harmonics will not appear the same. In a multiple-transit-spot case, we have τ_life > P_rot, and the height of the rotation peaks as a function of the frequency is given by the transit shape (see left side of Fig. 1). This can explain the results from Santos et al. (2017), where they find that the peak ratio between the first and second rotational peaks depends essentially on the transit profile of a single-spot in an LC. In the single-transit-spot case, we have τ_life < P_rot, the rotation peaks widen and overlap each other. As a result, the rotation peaks disappear and the spectrum comes close to the description of Harvey et al. (1993) for the solar case (see right side of Fig. 1). The information between τ_life and τ_tr is thus degenerated.

2.2. Gaussian case for the spot evolution

As a first step, we use Gaussians to describe E_life and F_transit. This description is used later in Section 3 for the simulated LC. These simple descriptions do not include the complexity of real star spots (spot shape, spot penombra, group of spot, etc.) but instead allow for an analytical approach in order to understand how the modelled spectra evolve when spots characteristics are modified. Indeed, similarly to Mosser et al. (2009), a Gaussian function can be used to describe the evolution of the spot contrast over time:

$$\begin{array}{c}\hfill E_{\mathrm{life}}(t)=exp\left(-\frac{4ln2·t^2}{{{\tau }_{\mathrm{life}}}^2}\right),\end{array}$$(4)

where τ_life is the characteristic lifetime of the spot seen as the width at half-height.

The F_transit is described using the same Gaussian function:

$$\begin{array}{c}\hfill F_{\mathrm{transit}}(t)=exp\left(-\frac{4ln2·t^2}{{{\tau }_{\mathrm{tr}}}^2}\right).\end{array}$$(5)

In the Fourier domain, and using Equation (3), we have

$$\begin{array}{cc}\hfill \widehat{S_0}(\nu )& =[\frac{{\mathrm{Ш}}_{1/\mathrm{Prot}}}{{\mathrm{P}}_{\mathrm{rot}}}\times \frac{\pi }{ln2}\frac{{\tau }_{\mathrm{tr}}}{2}{\mathrm{e}}^{-\frac{{(\pi \nu {\tau }_{\mathrm{tr}})}^2}{4ln2}}\ast \frac{{\tau }_{\mathrm{life}}}{2}{\mathrm{e}}^{-\frac{{(\pi \nu {\tau }_{\mathrm{life}})}^2}{4ln2}}]\hfill \end{array}$$(6)

$$\begin{array}{cc}& =\frac{\pi }{4ln2}{\tau }_{\mathrm{life}}{\tau }_{\mathrm{tr}}[\frac{{\mathrm{Ш}}_{1/\mathrm{Prot}}}{{\mathrm{P}}_{\mathrm{rot}}}\times {\mathrm{e}}^{-\frac{{(\pi \nu {\tau }_{\mathrm{tr}})}^2}{4ln2}}]\ast {\mathrm{e}}^{-\frac{{(\pi \nu {\tau }_{\mathrm{life}})}^2}{4ln2}}\hfill \end{array}$$(7)

$$\begin{array}{cc}& =B·[{\mathrm{Ш}}_{1/\mathrm{Prot}}\times {\mathrm{e}}^{-\frac{{(\pi \nu {\tau }_{\mathrm{tr}})}^2}{4ln2}}]\ast {\mathrm{e}}^{-\frac{{(\pi \nu {\tau }_{\mathrm{life}})}^2}{4ln2}},\hfill \end{array}$$(8)

where $B=\frac{\pi {\tau }_{\mathrm{tr}}}{4ln2{\mathrm{P}}_{\mathrm{rot}}}{\tau }_{\mathrm{life}}$ is a constant proportional to both the lifetime, τ_life, and the ratio τ_tr/P_rot, which is correlated with the inclination of the stars.

Knowing that the observed inclination i follows a probability distribution proportional to sin i, we can assume that the majority of stars from the Kepler mission that we use in this paper have inclinations much closer to 90° than 0° (inclination axis parallel to the line of sight). The detection methods used to find the rotation periods are sensitive to the modulations generated by the spots, which are less and less visible when the star inclination is low (see Sect. 4.1). This assertion should be treated with caution, as Bowler et al. (2023) found inclinations of down to 3.6° for stars whose rotation had been measured photometrically by periodograms, a method similar to the one used by Santos et al. (2019, 2021) (see Section 4.1 for more details). However, the very small inclinations found in this paper remain a minority: 13% of these stars correspond to inclinations below 30°, and more than 50% above 60°. Assuming that the majority of stars of the sample use in this paper have a high inclination, we have ${\tau }_{\mathrm{tr}}\sim \frac{{\mathrm{P}}_{\mathrm{rot}}}{4}$ and $B\sim \frac{\pi }{16ln2}{\tau }_{\mathrm{life}}$.

Other functions can be used to describe F_transit and E_life, that lead to the same physical interpretations. For example, we use Lorentzians or pseudo-Lorentzians in the Fourier domain when applying the model to the Kepler data in Section 4, since their shape appears to be more adapted to the observed spectra.

2.3. General case

When considering an LC including N spots (with N ≫ 1), with each associated with a maximum depth A_i(T_i, S_i) and a random appearance time t_i, the LC can be described as

$$\begin{array}{c}\hfill S(t)={\displaystyle \sum _{i=1}^N}{A}_i·S_0(t-t_i).\end{array}$$(9)

In the Fourier domain, this yields

$$\begin{array}{cc}\hfill \widehat{S}(\nu )& ={\displaystyle \sum _{i=1}^N}{A}_i·\widehat{S_0(t-t_i)},\hfill \end{array}$$(10)

$$\begin{array}{cc}& ={\displaystyle \sum _{i=1}^N}{A}_i·{\mathrm{e}}^{2\mathrm{i}\pi \nu t_i}\widehat{S_0}(\nu ).\hfill \end{array}$$(11)

This yields the following PSD:

$$\begin{array}{c}\hfill |\widehat{S}{(\nu )|}^2=|A_{\mathit{tot}}{|}^2·{|\widehat{S_0}(\nu )|}^2,\end{array}$$(12)

with

$$\begin{array}{cc}\hfill |A_{\mathrm{tot}}{|}^2& =|{\displaystyle \sum _{i=1}^N}{A}_i·{\mathrm{e}}^{2i\pi \nu t_i}{|}^2\hfill \end{array}$$(13)

$$\begin{array}{cc}& ={\displaystyle \sum _{i=1}^N}|A_i{|}^2·{|{\mathrm{e}}^{2i\pi \nu t_i}|}^2+{\displaystyle \sum _{i=1}^N}{\displaystyle \sum _{j=1}^N}{A}_i{A}_j{\mathrm{e}}^{2i\pi \nu t_i}{\mathrm{e}}^{2i\pi \nu t_j}.\hfill \end{array}$$(14)

Since the spots are considered incoherent, we have the following:

$$\begin{array}{c}\hfill |A_{\mathrm{tot}}{|}^2\simeq {\displaystyle \sum _{i=1}^N}{|A_i|}^2.\end{array}$$(15)

Hence, A_tot corresponds to the sum of the photometric contrast contributions from the N spots to the LC. Using the results from Section 2.2, the height of the spectrum ℋ is proportional to |A_tot ⋅ B|², which yields

$$\begin{array}{c}\hfill \mathcal{H}\sim |A_{\mathrm{tot}}·{\tau }_{\mathrm{life}}{|}^2.\end{array}$$(16)

Therefore, the amplitude of the power spectrum gives an indication of the coverage in terms of photometric attenuation (given by A_tot and due to the spot surface and spot temperature difference with its surroundings) during the entire observation time. The dependency of ℋ on τ_life is the result of the contribution added by each spot transit over time.

When considering the presence of both single- and multiple-transit spots, their two contributions can be separated: LC(t) = S_single(t)+S_multi(t). Because of the linearity of the Fourier transform, these two contributions are also separated in the Fourier domain, but the three terms appear when considering their power spectrum.

For the same reason as before, we can consider that the two contributions are incoherent and neglect the cross term:

$$\begin{array}{c}\hfill |\widehat{\mathit{LC}}{(\nu )|}^2={\widehat{S}}_{\mathrm{single}}^2(\nu )+{\widehat{S}}_{\mathrm{multi}}^2(\nu ).\end{array}$$(17)

There are hence two contributions to the PSD, as illustrated in the bottom of Fig. 1 and Fig. 2 for a realistic general case: one due to the single-transit spots, which are dominated by the transit time (orange box on Fig. 1 and in dotted orange on Fig. 2), and another from the multiple-transit spots, which makes the rotation peaks appear and yields information on the spot mean lifetime (blue box on Fig. 1 and in dotted light green on Fig. 2). The combination of these two contributions gives the LC fit in blue, which as we can see matches the large-scale trend of the data in grey. Depending on the number of spots with a single or multiple transits, we can then identify two categories of spectra: spectra with clear rotation peaks (when spots with multiple transits dominate), hereafter the “with peaks” spectrum or without rotation peaks (when spots with a single transit dominate), hereafter the “peakless” spectrum.

Fig. 2. Left panel: Example of a power spectra density of a Kepler LC (grey) ‘with peaks’ fitted by the model detailed in Sect.2. The transit component due to the single-transit spots is dash-dotted orange line, the multiple-transit spots component in a dashed light green line and the total fit in blue. The rotation peaks are marked by purple dashed lines. Right panel: Example of a peakless power spectra density of a star where the Harvey profile is used for the fit.

2.4. Summary of the model

In summary, there are two profiles that have to be taken into account when an LC is fitted: first, the rotation peak profile ${\widehat{S}}_{\mathrm{multi}}^2(\nu )$ in light green in Fig. 2, which is due to the multiple-transit spots and is related to the spot lifetime τ_life. Second, the transit profile ${\widehat{S}}_{\mathrm{single}}^2(\nu )$ in light blue in Fig. 2, which is due to the single-transit spots and is closely related to the transit time τ_tr. In the case of a star covered only by spots that last less than one rotation period (single-transit-spot case), this profile is degenerated between the transit (τ_tr) and the lifetime (τ_life) components, and we recover a profile that can be better modelled by a Harvey profile described by only one characteristic timescale (see right panel of Fig. 2).

Two types of spectra can therefore be inferred: those from stars covered mainly by multiple-transit spots where the two profiles are needed, and peakless ones where the two profiles are degenerated. This explains why some spectra show rotation peaks that are due to multiple-transit spots, while others do not. In the remainder of this paper, we refer to ‘with peaks’ spectra (left panel of Fig. 2) and ‘peakless’ spectra (right panel of Fig. 2).

Based on this total profile, three spot characteristics can be extracted from the spectrum of an LC:

1. The transit proxy τ_tr: the inverse of the width of the very low frequency component of the spectrum. It is the characteristic time of the transit and is thus close to P_rot/4 but can vary depending on the rotation-axis inclination or limb-darkening.

2. The lifetime proxy τ_life: the inverse of the width at half-height of the rotation peaks, if they are visible. It provides information about the mean lifetime of the stellar spots during the period of observation.

3. The spot impact proxy ℋ: the height of the spectrum, the sum of ℋ_single and ℋ_multi. This proxy gives information about the total spot contrast coverage on the stellar disc during the observation time. It is linked to two spot properties: their surface and their temperature difference compared to unspotted areas. This thus gives degenerate information about whether the spots are cold or small, or hot or large. This proxy is also sensitive to the inclination of the star. However, given that we are looking at average trends, and that these two phenomena are linked to the magnetic field, this proxy allows us to make an estimate of the global spot activity of the star over the duration of the observation.

We emphasise that the two temporal proxies (τ_life and τ_tr) are only relevant for the stars showing rotational peaks, while they are degenerated for the stars not showing peaks.

3. Validation of the model with simulations

In this section, the model described in Sect. 2 (Eq. (17)) is validated against simulated LCs based on the analytical model of Dorren (1987), which describes how a single circular spot impacts an LC as a function of different stellar and spot properties. By computing the impact of several spots and summing these contributions, this model permits us to simulate LCs. We can modify several parameters such as the number of spots, their lifetime, their distribution on the stellar disc, their radii, the rotation period of the star, the possible differential rotation (DR) and the inclination of the star rotation axis. The code also uses the Gaussian function given by Eq. (4) (Mosser et al. 2009) to describe the intrinsic evolution of the spot. The parameters explored are: the influence of the rotation period (P_rot), the spot lifetimes, their number, and their radius. The spectrum of the simulated LCs is then computed and fitted using the Maximum Likelihood Estimator (hereafter MLE) described by Toutain & Appourchaux (1994).

3.1. Method

We compute several sets of LCs for various stellar parameters. Each set is composed of an ensemble of 200 simulated LCs, using the same stellar parameters but a different random set of spots on the stellar disc in terms of time and location (following a normal distribution).

For the spot lifetime, we performed three different types of simulations: one with only spots that last less than the rotation period to produce peakless spectra (τ_life = 1/6P_rot hereafter “single transit”). It corresponds to the single-transit-spot case and is close to solar observations. Another type with only multiple-transit spots (i.e. τ_life = P_rot, hereafter “multi-transits”), where the spots have more than five visible transits with a significant depth in the LC (since the lifetime and transit time were defined as the width at half height). The last type has lifetimes that randomly vary according to a uniform law between 1/6P_rot and P_rot (hereafter ‘mixed’). The simulations with only multiple-transit spots are not realistic, since there should be shorter additional photometric contributions coming from physical phenomena other than spots (for instance, faculae or flares), but they serve as an extreme test case. These simulations also covered three rotation periods of 10, 20, and 30 days. For a given set, all the spots have the same area, and the LCs last for two years, similar to most of the subsequently used Kepler LCs.

Since the contrast of the spot is degenerated between temperature and area, we choose to only change the spot surface area to study how ℋ varies in our set of simulations. Other sets have been generated, where the size and number of spots varied, as well as the inclination of the stars. For each normalised LC, a Gaussian noise with a deviation of σ = 10⁻⁵ has been added. This yields to 23 sets of simulations, each set comprising 200 simulated LCs. A summary of the parameters associated with each set is given in Table 1.

Once generated with its associated stellar parameters, each LC spectrum has been fitted with the model proposed in this paper, meaning with a combination of a function for the transit profile (which represents the dominant contribution of the single-transit spots), plus a function for the rotation peaks (which represents the dominant contribution of spots with multiple transits, see Eq. (17)).

For these simulations, each spectral component associated with the transit and the time evolution profiles (i.e. each rotation peaks) can be expressed by the following Gaussian profile:

$$\begin{array}{c}\hfill G(t)=A{\mathrm{e}}^{-{(\frac{(\nu -{\nu }_0)\pi \tau }{ln(2)})}^2},\end{array}$$(18)

since the temporal evolution profile of the spot is a Gaussian, and their transit can also be approximated by a Gaussian. The set-up of the simulations is thus close to the description made in Sect. 2.2, but here we only fit the contribution of the multiple-transit spots ${\widehat{S}}_{\mathrm{multi}}(\nu )$ (i.e. the rotation peaks) by using each time the Gaussian form of Eq. (18). We do not take into account the transit component in ${\widehat{S}}_{\mathrm{multi}}(\nu )$ since it only appears in the amplitude of the peaks, which is a free parameter in the algorithm to help it converge. We choose to fit 2 rotation peaks plus the contribution centred at the zero frequency, because in most of the cases the higher harmonics are not clearly visible above the noise. For ${\widehat{S}}_{\mathrm{single}}(\nu )$, there is only one Gaussian that describes the transit time, since there are no clear rotation peaks. Moreover, the simulated LCs have also been fitted by the pseudo-Lorentzian model proposed in Eq. (1) of Harvey et al. (1993) in order to compare our new model with this already existing one.

3.2. Results

3.2.1. Spot impact proxy

The sets of simulations used here to test Eq. (16) are listed in Table 1.

Fig. 3 shows a clear correlation for all sets between the median estimate of ℋ and the surface covered by the spots. The size of the individual spots and the temporal behaviour of spots (mixed and multi-transit cases) do not impact significantly the trend. However, simulations with only single-transit spots show a different dependency, especially when considering different rotation periods, where ℋ tends to be overestimated. This may be due to the fact that, in the case of the peakless spectrum, the model does not lift the degeneracy between the signature of the lifetime and the transit, and the approximations made to express ℋ in Eq. (16) are less valid. The impact proxy ℋ is still an increasing function of |A_tot ⋅ τ_life|² for a given inclination, but nevertheless seems to slightly deviate from a linear relation.

Fig. 3. Median height of the spectrum (ℋ) of each set of simulations as a function of |Atotτlife|2 (see Eq. (16)). The colour of the dots indicates the rotation period of the star (10 days in green, 20 days in red and 30 days in purple), while the shape gives the type of spot lifetime (diamonds for the peakless stars, hexagons for mixed spots and 5-pointed stars for stars with rotation peaks spectra). The size is proportional to the surface of the spots. The large hexagon on the top-right corresponds to set 19 with spots radii of 10 degrees. The error bars are based on the first and third quartile. The points with a shading towards white are the results for the 3 simulations where the inclination varies (sets 20, 21 and 22), white being for 0° inclination (set 22).

Finally, we have tested the impact of the inclination of the stellar rotation axis using three more sets of simulations with mixed spots (sets 20, 21 and 22) for inclinations of 60° (light red), 30° and 0° (white) compared to the line of sight. The estimations of ℋ for an inclination of 30° and 0° are overlaid. ℋ is higher compared to the case with a 90° inclination. We observe that the closer the axis is to 0°, the higher is the value of ℋ. This deviation from Eq. (16) is due to the approximation τ_tr ∼ P_rot/4 made in the expression of B in Eq. (8), which is no longer true for inclination angles near to 0° where τ_tr is closer to P_rot. The spot impact is therefore higher in these configurations, since the mean transit time is longer as the mean of the spot signal over time is higher.

3.2.2. Lifetime proxy

The sets used in this section are 0, 1, 2, 3, 5, 6, 7, 8, 9, 20, 21 and 2. The results for this proxy are represented in Fig. 4. First, we can see that the rotation period does not affect the spot lifetime measurement much (similar results between the different columns). However, for mixed and multiple-transit spots (middle and bottom rows), the distribution is as when the rotation increases.

Fig. 4. Normalised histograms of τlife. Each column corresponds to a different rotation period, and each row to a type of lifetime (single-transit, multi-transit and mixed spots). Here, only the results for the new model are shown, since it cannot be retrieved using the Harvey model. For each panel, the red lines show the expected value. The vertical black dash-dotted line represents the lifetime value for the single-transit spots and the dashed the one for the spots that last several transits (i.e. τlife = 1Prot). In the case of mixed spots, these two lines are represented, as they show the range of possible values. This space has been coloured in blue for clarity. The full line represent the mean value of the lifetimes. For single-transit and mixed-spots simulations, the vertical black dotted line represents the value of lifetime where two transits are visible in the LC. The stepped histograms in green, blue and orange correspond to the simulations where the inclination of the rotation axis varies (respectively to 0, 30 and 60 degrees).

Concerning the multiple-transit-spot cases (bottom row), the mean lifetime is as expected equal to the rotation period (see Table 1), with a deviation of approximately 10 %.

For the mixed spots cases (middle row), the distribution is spread between the maximum value and P_rot/3, This value corresponds to the case where the spots have a two visible transits (τ_life = 2 ⋅ 1/6P_rot) in the LC and where the rotation peaks begin to appear. The median value is therefore located at the middle of this interval (2/3 ⋅ P_rot), the value around wihch the distribution is located.

In the case of single-transit-spot simulations (top row), the fit has a weaker agreement with the proxy. Since in this case there are none or very weak rotational peaks, as in the case of the transit proxy, the fitting algorithm encounters more difficulties to converge properly. However, there is still an accumulation around P_rot/3 (dotted line in Fig. 4). Here also, this is because some single-transit spots can sometimes produce two transits in the LC, typically when they appear right before passing behind the limb, and then re-appear at the opposite limb and fade very quickly. This induces a periodic signal that can be seen in the spectra, even for a single-transit spot. Despite this, the results show that the lifetime values found by the model are spread over a wide range of values, which also covers the lifetimes of multiple-transit-spot cases. In the case of peakless spectra produced by single-transit spots, the model therefore fails to converge reliably to a correct lifetime value, which is expected: in this case, both the lifetimes and transit times of the spots are degenerate.

The simulations where the inclination of the stars have been varied are also presented in Fig. 4. We can see that for the minimum inclination angle (0°), the values found for τ_life are biased towards shorter lifetimes, by a relative value of 30%. For the other two inclinations, the values found remain very close to the values without inclination.

3.2.3. Transit proxy, τ_tr

In this section, only the sets with varying lifetimes and rotation periods are used, since τ_tr is only sensitive to temporal parameters (see 0, 1, 2, 3, 5, 6, 7, 8, 9). For this proxy, we compare the results from our new model proposed in this paper to those from Harvey (1985). Indeed, the width of the pseudo-Lorentzian profile has been related to the characteristic time of the magnetic activity; but according to our model, this value is closer to the transit time. Fig. 5 shows the τ_tr distributions for each set of simulations.

Fig. 5. Same as Figure 4 but for the transit proxy τtr. The histograms in blue correspond to the value found by our model, and the ones in orange to the value found by the model of Harvey. The vertical black dashed lines indicate the expected value for each set of simulations, which corresponds to Prot/4. The coloured dashed lines represent the median value of each distribution.

For single-transit-spot cases (top rows), both models underestimate τ_tr. This discrepancy can be explained by the fact that, in this case, the transit time and the lifetime are degenerate. Although the rotation peaks do not appear clearly, they still contribute to the broadening of the spectrum and their presence is included in the fit function. Since the contribution of the transit is not distinguishable from the rotation peaks, the fit converges to a biased value.

For the mixed-spots case (bottom rows), the expected value is found by the algorithm in the three different rotation cases, but with shifted values for the Harvey model for fast rotators. The Harvey model also has an accumulation of points for higher values, in both single-transit-spot and mixed cases. This is due to the fitting algorithm, which in some cases will converge to the fundamental peak, which is located at a shorter interval than the transit value, because these are not taken into account in the Harvey model.

The case with only multiple-transit spots is not shown here because the model we use cannot recover a correct transit time without the contribution of the single-transit spots that we have chosen not to add in these simulated LC. However, according to our analytical description, the transit time is still present when following an envelope passing through the maximum height of the rotational peaks. This particular case is discussed in Appendix A.

This result shows that the Harvey profile is closely related to the transit time of the spots (and hence, to the rotation period), and not to the characteristic time of the activity, or more specifically the spot lifetime. The results obtained by Aigrain et al. (2004) can be explained by the fact that, in the case of the Sun, the spot lifetime is close to the rotation period. It is a “peakless spectra” case, where the rotation peaks are less apparent and the dominant profile is the one from the transit.

The inclination of the stellar rotation axis also modifies the estimation of τ_tr. The distribution of values is significantly wider and shifted towards higher values, as expected (the transits are longer when the star is inclined). The results are not included on Fig. 5, but the median values of the distribution are shifted to 6.2, 15 and 20 days for an inclination of 60°, 30° and 0°, respectively.

3.2.4. Summary of the simulations results

Analysing the profiles of the spectra in the Fourier domain allows us to validate the interpretation of the proxies proposed in Sect. 2.4 and determine their limitation :

• The spot impact proxy ℋ shows a good correlation with the contrast induced by the surface covered by the spots. The inclination of the rotation axis can have an impact on this proxy since the spots appear longer, but we can nevertheless assume that the catalogues of observed stars only contain a small fraction of highly inclined stars, as mentioned before and according to Bowler et al. (2023).

• The lifetime proxy τ_life is correlated with the mean value of the spot lifetime, when more than two transits per spot are observed, (i.e. when ${\tau }_{\mathrm{life}}>\frac{1}{3}{\mathrm{P}}_{\mathrm{rot}}$ according to the values of simulation). Otherwise, the model has difficulty to converge to an exact value, since τ_tr and τ_life are degenerated. It is therefore only possible to evaluate this proxy if there are significant rotation peaks in the spectrum under study. This means that the estimate of the mean lifetime is biased by the rotation rate: for fast rotators, we will be biased towards shorter τ_life (if any measured) since τ_life cannot be estimated if longer than 1/3P_rot. If there is no significant rotation peaks, the Harvey model is sufficient, but the calculated characteristic time cannot be used to estimate an average spot lifetime. The inclination of the rotation axis can also underestimate τ_life when the angle of inclination is very small.

• The transit proxy τ_tr is well correlated with the mean transit time of spots. It is underestimated for LCs with only single-transit spots since there is a degeneracy between the lifetime and the transit time in this case. For stars with multiple-transit spots, its value can be recovered through the contribution of single-transit spots to the spectrum, but also in the height of the rotational peaks. However, to simplify the fitting algorithm, each peak has been adjusted individually so this component has not been included in the fit function. The inclination of the stars may bias the estimation of τ_tr if this inclination is very pronounced. However, this proxy is not a primary focus of this paper.

For the remainder of the analysis using Kepler data, it is thus necessary to differentiate between spectra with and without observable peaks. This distinction ensures that we apply the model only when it can estimate lifetimes effectively. In the case of ‘peakless’ spectra, we apply the Harvey model, which still provides an estimate of ℋ.

4. Application to the Kepler data

4.1. Datasets

The LCs used in this work are part of the 150 000 stars observed by the Kepler mission. They have been taken from the MAST (Mikulski Archive for Space Telescopes) data archive. Although the Kepler mission made observations over a 4-year period, the LCs used here last for an average of two years, due to the data reduction employed. These data have been first detrended by the Presearch Data Conditioning-Maximum A Posteriori (PDC-MAP, e.g. Jenkins et al. 2010; Smith et al. 2012; Stumpe et al. 2012). Another treatment has been applied, described by de Assis Peralta et al. (2018)¹, where the main goal was to fill the gap between the quarters using a simple stitching method, to remove the outliers and to apply a linear trend correction.

Two catalogues of data have been used in this work. The first is from McQuillan et al. (2014), which provides an estimate of the rotation period for 33 040 Kepler stars using the autocorrelation function method (ACF, method detailed by McQuillan et al. 2013). The catalogue also contains other stellar parameters such as the effective temperature T_eff, the log g extracted from the Kepler Input catalogue (KIC) or estimated from Dressing & Charbonneau (2013), and an estimation of the mass, by using the isochrone from Baraffe et al. (1998), which was computed specially for low-mass stars. The sample is composed only of main-sequence stars, where the known eclipsing binaries and stars with a planet candidate (KOI, Kepler Objects of Interest) were removed.

The second catalogue used here is the concatenation from Santos et al. (2019) and Santos et al. (2021). This sample is made of 55 232 Kepler main-sequence stars and subgiants from spectral type M to F. For this paper, we removed the subgiants (5544 stars) from this sample to focus only on the main-sequence stars. This catalogue also provides an estimation of the rotation period by using a combination of three methods: the first one is the ACF, the second one is a global wavelet power spectrum analysis and the third one is based on a composite spectrum, which is a combination of the two first ones (see Ceillier et al. 2017). The catalogue also contains an estimate of the mass, T_eff and log g from Mathur et al. (2017). Most of McQuillan’s stars can also be found in the Santos set. A total of 25 855 new stars are added from this catalogue that are absent from McQuillan’s catalogue.

It is important to point out that the methodologies used by these two studies to measure rotation periods can induce errors in their estimates up to a factor of two. Despite the meticulous efforts in these papers to mitigate this issue, techniques such as autocorrelation function (ACF) analysis or peak identification can occasionally converge to a harmonic rather than the actual period. This error can also come from the stellar configuration where two large spots are present on the star at opposite locations. Some periods used in this article may therefore be incorrectly estimated.

The difference between these two catalogues relies on the method used to determine P_rot. Since McQuillan et al. (2014) used the ACF only, one can assume that the procedure of Santos et al. (2019) and Santos et al. (2021) is more sensitive to the variability and includes fewer active stars. For these two catalogues, we also used the Gaia mass estimation by FLAME², when available. However, since we observed a certainly biased accumulation of stars with mass around 0.7 M_⊙ in the Gaia data, we chose to only use Gaia masses higher than 0.8 M_⊙ and kept the mass value from the two catalogues otherwise. In the main part of this section, we analyse the results simultaneously for both McQuillan’s and Santos’ datasets.

4.2. Algorithms for the Kepler data

Two algorithms were developed to analyse the Kepler data. The first is a sorting algorithm that classifies spectra with and without peaks by examining the first four rotation peaks using three different methods. The first method analyses the derivative of the spectrum to detect significant changes in sign. The second method checks whether the peak heights are at least three times greater than the overall standard deviation of the spectrum. The third method evaluates whether each peak height is at least three times higher than the local standard deviation immediately after the peak. Thus, 12 tests (3 methods applied to 4 peaks) are performed on each spectrum. A spectrum is considered to have peaks if all three methods agree on at least two of the first three peaks, or if 8 out of the 12 tests are positive. Conversely, a spectrum is classified as lacking peaks if none of the methods converge on any of the peaks, or if fewer than 5 tests are positive. This algorithm was tested on a sample of 1000 spectra that had previously been sorted by visual inspection. In 87% of cases, the algorithm found the same results as visual sorting.

The second algorithm developed was the fitting algorithm. The same algorithm as before (Sect. 3.1) has been used to fit the LC power spectrum, but this time a pseudo-Lorentzian (Eq. (1)) profile has been used for S_single and a simple Lorentzian (α = 2) for the rotation peaks (i.e. S_multi in the case of the stars identified as “with peaks” by the sorting algorithm). These functions indeed allow for more degrees of freedom with the α exponent, and were thus more suited for actual observations. Since the aim is to retrieve a width at half-height and compare global trends for the dataset, and not to find exact lifetime values, changing the fit function should not impact the results, as explained in Sect. 2.2. The fitted spectrum takes into account two rotation peaks because that is the most common case observed. However, in some cases, the LC spectrum can show more than two peaks that can bias the fit. Given that, it will especially impact the estimation of the transit time (see Sect. 3) that is not the main proxy of interest here. This compromise can thus be considered acceptable. In the fitting process, the rotation period is fixed to the value extracted from the catalogues. The Harvey profile (i.e. only a pseudo-Lorentzian) has been used for the stars with peakless spectra. We recall that an example of these two fits is shown in Fig. 2.

The total catalogue that we were initially using consisted of more than 50 000 stars. After running the sorting algorithm, 57.4% of the data were classified as spectra with peaks, 37.1% without and 5.5% as uncertain for the whole dataset. We did not use the spectra classified as “uncertain” in our study which did not impact our sample since they represent a very small portion of the sample. After running the fitting algorithm on the data, we removed the spurious fit values (zero values, infinite values, negatives, etc.). We also removed the stars that do not respect the predicted criterion for the appearance of rotation peaks (see Sect. 3.2.4). In fact, we observe a small group of values below τ_life/P_rot ≈ 0.35, which, after a visual verification, corresponded to peakless spectra that the sorting algorithm did not classify properly. We end up with a sample of 33 755 stars. More technical details about this selection are given in Appendix B.

4.3. Results

In order to present the results obtained for ℋ and τ_life and present a landscape of magnetic spot activity as a function of different stellar parameters, we have chosen to start by plotting them as functions of stellar rotation and mass (Fig. 6). The position of groups with and without peaks is also shown. These distributions are compared to the stellar Rossby number:

Fig. 6. Top panel: Spot-impact proxy for the stellar sample as a function of mass and rotation period. The colour scale shows the impact proxy ℋ with a logarithmic scale. Yellow is for a large spot impact (active stars) and black for low spot impact (less active stars). The black lines represent the iso-lines of the Rossby number. The ⊙ point indicates the location of the Sun in this diagram. Middle panel: Same figure with contour lines showing the locations of the stars with ‘peakless’ or ‘with peaks’ spectra, respectively in orange and blue. The grey dots represent the whole sample. Bottom panel: Same as the top panel but with the lifetime proxy τlife. Only stars with spectra with peaks are presented.

$$\begin{array}{c}\hfill \mathrm{Ro}=\frac{{\mathrm{P}}_{\mathrm{rot}}}{{\tau }_{\mathrm{c}}},\end{array}$$(19)

where τ_c is the turn-over timescale for the convection. This number helps determine which effect is dominant between the rotation and the turbulent convective flows, which are the main two processes involved in the dynamo mechanisms responsible for the generation of the magnetic activity. This number is computed using the mass-dependent relation of Wright et al. (2011), where τ_tc has been computed using the colour-conversion relation of Pecaut et al. (2012). This relation has been established only for stars with masses between 0.2 and 1.2 M_⊙ and is no more valid for higher masses in the sample. We begin by describing the results for ℋ, then τ_life, and finish with a summary of the different trends that we observe as a function of the stellar mass and the Rossby number. Regarding the estimation of the transit time and the characteristic time associated with the Harvey profile, we refer the interested reader to Appendix C as the physical information provided by these parameters is not essential to the following discussion.

4.3.1. Spot impact proxy ℋ

The top panel of Fig. 6 shows the distribution of ℋ values. We recall that we are using the estimations of our model for spectra with peaks and the model of Harvey for peakless spectra.

Different trends emerge from this fig. First of all, there is a clear distinction between stars below and above one solar mass. For stars above 1 M_⊙, the spot impact is much lower. For the less massive stars, the spot impact evolves differently depending on the Rossby number. The most active stars (in yellow) are seen in a banana shape between the iso-Rossby lines Ro =  0.7 and 1 (hereafter Region 2). We refer to this area as the high activity regime. Then, on each side of the high activity regime, there are two zones where the spot impact proxy shows intermediate values (orange). The one on the faster rotator side (left) is located at smaller Rossby number values of Ro= 0.5 (hereafter Region 1). The other one, on the longer period side (right), corresponds to Rossby numbers greater than 1 (hereafter Region 3).

4.3.2. Localisation of peakless and with peaks spectra

The middle panel of Fig. 6 presents the location of stars with peakless spectra compared to those with peaks. The “with peaks” population (blue) concerns mainly fast rotators and stars below 1 M_⊙. They are massively present around the isoline of Ro = 1 and below, as well as on the branch of the fast rotators described by McQuillan et al. (2014), who observed that the rotation period distribution is divided into two branches for low effective temperature stars. The stars with peakless spectra (orange) are either slow rotators located above the isoline Ro = 1 or more massive stars (M > 1 M_⊙) that have small P_rot.

4.3.3. Lifetime proxy τ_life

The bottom panel of Fig. 6 presents the results for the mean lifetime of stars with rotation peaks. Similar to the spot impact, this proxy is related to the stellar mass: the more massive stars have very short lifetimes for their spots, while the low-mass stars have longer ones.

There is also a dependence on the rotation period: slow rotators have longer lifetimes than fast rotators, which is not surprising, given the Sect. 3.2.4 conclusions: the model struggles to reliably estimate τ_life for values below P_rot/3. Consequently, the slower a star rotates, the more challenging it becomes to accurately determine shorter lifetimes. This relationship between rotation and lifetime is illustrated more clearly on Fig. 7. The same trends has been observed by Basri et al. (2022) for the MacQuillan stars that we have in common, but we do not retrieve the diffuse second branch that they described. Similarly to the spot impact proxy, but less distinctly, τ_life of the less massive stars is organised according to the Rossby number: longer lifetimes appear between Ro = 0.7 and Ro = 1, in the same location Region 2 described in Sect. 4.3.1 where the higher ℋ are found. In Region 1 (below Ro = 0.7), longer lifetimes can also be seen under 0.6 M_⊙. This can also be seen on Fig. 7.

Fig. 7. Estimated mean lifetime τlife as a function of the rotation period for stars with spectra with rotation peaks.

But what can we say about the “peakless” group of our sample? The middle panel of Fig. 6 shows that for the less massive stars (M < 0.6 M_⊙), they have mainly slow rotation period. For these particular stars of our sample, we do not have estimations of τ_life but we still know that these stars have spot with lifetimes shorter than or equal to their rotation periods (see Sect. 2.4). As the definition of lifetime is the width at half-height, this means that the mean lifetime of spots on these stars would be shorter than their rotation. This means that this population located at large Rossby values (i.e. above Ro ∼ 1.2) has shorter lifetimes than the majority of those observed for the group with rotation peaks in this area. This statement is not so obvious for high-mass stars as there are very few stars in the group with peaks for comparison. However, we can intuitively expect that these starspots likely have short lifetimes, given their rapid rotation.

4.4. Summary: Activity and Rossby number

As seen in Sects. 4.3.3 and 4.3.1 the main proxies (ℋ and τ_life) used to characterise spot signature in LCs depend on stellar properties such as mass, P_rot and Rossby number, which we use to define stellar categories.

First, as seen in Fig. 6, stars with M > 1 M_⊙ show smaller values of ℋ regardless of other stellar properties, whereas stars with M < 1 M_⊙ show larger values of ℋ and τ_life presenting a clearly structured behaviour. We thus focus on stars of lower mass, which are in addition more numerous than the more massive stars.

The histograms of Fig. 8 summarise the different trends of this population, providing a clearer illustration of their statistical distributions. The Rossby number being a fundamental index when dealing with magnetic activity, we use it to define categories of stars depending on its value and Fig. 6: Region 1 with Ro < 0.7, Region 3 with Ro > 1, and Region 2 between these two values. The distributions of stars for these regions is shown in panel a of Fig. 8, as well as the distribution of peakless stars (in orange) having clearly higher Ro. This finally leads to four categories of stars.

Fig. 8. Histograms for the three different regions and the peakless sample as a function of different stellar parameters. In blue: Sample of stars with rotation peaks. The different shades of blue correspond to the three regions of Fig. 6 cut according to the Rossby number. This slicing is illustrated in panel (a) as a function of the Rossby number. The orange histograms represent the peakless group. From top to bottom, the histograms are for (a) the Rossby Number, (b) spot impact proxy, (c) the mean lifetime τlife,(d) the rotation period, and finally (e) the ratio Prot/τlife.

First, panel b shows these categories regarding the spot impact ℋ. These four categories show wide distributions of spot impact ℋ but with different median values, Region 2 having the highest median while Regions 1 and 3 have lower (and almost equal) medians, and peakless stars have the lowest median spot impact.

Second, panel c present their lifetime proxy τ_life. Lifetimes (τ_life) of the stars from the 3 regions (lifetimes are not defined for peakless stars) also show wide distributions but again with different median values: Region 1 being the shortest, Regions 2 and 3 being longer (panel c of Fig. 8) Regarding the bias towards long lifetimes for slow rotators, the distribution in Region 1 is likely the closest to the underlying distribution, while the other two regions are increasingly biased. However, examining the median values reveals that stars in Region 2 tend to have slightly longer lifetimes than region 3. This suggests that Region 2 may be less affected by methodological bias, with its distribution being closer to the underlying one. Additionally, panel e shows that stars in Region 2 have a τ_life/P_rot ratio close to 2, whereas stars in Region 3 exhibit a ratio closer to 1. This latter group, which has shorter lifetimes, is much more likely to be biased towards long lifetimes.

As mentioned before, we do not have lifetime estimations for the peakless group of our sample but the rotation period gives the upper limit of their characteristic lifetime. Fig. 8d shows that these stars have much slower rotations than the estimate mean lifetime of Region 2 and 3; they should therefore have shorter or equivalent lifetimes.

Third, panel d shows the rotation period of the three different regions. As they are based on the Rossby number, we can see that each of them corresponds to different rotation intervals that increase with their Rossby value. We observe that stars with the fastest rotation exhibit the shortest lifetimes (see Fig. 7).

Finally, panel e shows the ratio τ_life/P_rot. Because the P_rot/τ_life ratio is the parameter governing the appearance of rotation peaks, it is also interesting to represent its distribution in Fig. 8e. These distributions of τ_life/P_rot show quite different distributions than τ_life alone: Region 3 with the longest P_rot has the smallest τ_life/P_rot and Region 1 with the shortest P_rot has the largest τ_life/P_rot. In fact, this means that even if fast rotators have the shortest lifetimes, the spots of these stars transit over the stellar disc several times, creating clearly visible rotation peaks. It is precisely for this reason that fast-rotating stars often have better defined rotation peaks than slower-rotating stars like the Sun. Similarly, less visible rotation peaks in slower-rotating stars are not necessarily associated with shorter spot lifetimes.

Based on Figs. 6 and 8, the observed stars and their type of activity can be described as follows:

• The activity of stars under 1 M_⊙ shows trends depending on the Rossby number. When considering the median value Ro of the distribution shown in Fig. 8, these trends can be described as:

  • Ro < 0.7 (Region 1): These stars have intermediate values of spot impact and relatively short lifetimes, however increasing with decreasing mass. These stars correspond to the fast rotators branch of McQuillan et al. (2014).

  • 0.7 < Ro < 1 (Region 2): It corresponds to the most active stars with the highest spot impact and also the longest lifetimes.

  • Ro < 1 and Ro < 1.5 (Region 3): These stars have intermediate spot impacts as Region 1 and lifetimes shorter than those of Region 2 but longer than those of Region 1.

  • Ro > 1.2: It is the case of a large fraction of stars showing peakless spectra: they have also smaller spot impacts than the 3 other groups depicted before, which suggests a fourth type of activity. If we look at their rotation period, these stars are likely to have shorter mean lifetimes than the groups of stars discussed above.

• For more massive stars (> 1 M_⊙), the behaviour differs. Spot impact and lifetimes are significantly lower compared to lower mass stars, indicating that they are considerably less active. Additionally, stars with peakless spectra display similar low activity levels. This reduced activity in massive stars is not unexpected: these stars correspond to masses where the convective zone becomes thinner, resulting in a less efficient Parker-type dynamo effect.

5. Discussion

The main results of this work are now discussed with regard to previous studies on stellar activity and recent observations of magnetic fields by Zeeman-Doppler imaging. We address the potential insights brought by this work to stellar internal magnetism, and finally discuss the physical phenomena that can bias our methodology.

5.1. Comparison with other activity proxies

Two other common activity proxies have been studied in the past, namely the X-ray luminosity and the mean level of photometric variability S_ph.

Wright et al. (2011) studied the relation between stellar activity in X-rays and rotation period in a sample of about 800 low-mass stars. Fig. 2 of Wright et al. (2011) in particular shows a scatter plot from this sample of stars as a function of the Rossby number (from Ro = 10⁻³ to 3) and the X-ray luminosity. This scatter plot shows a very clear plateau for Ro values lower than ∼0.1 and a decrease above. We have reproduced the same figure, plotting ℋ as a function of the stars’ rotation Rossby number (left) and period (right) in Fig. 9.

Fig. 9. Scatter plot of the total sample of stars (grey dots in the background) as a function of the spot impact proxy and the rotation Rossby number (left) or period (right). On the left panel, the black line represents the median value of ℋ as a function of Ro for the whole sample. The hatched area around this line represents the interval of the quartiles. The dotted lines show the position of Ro = 0.7 and 1. For the right panel, these stars are grouped into slices of increasing masses (from 0.5 M⊙ in dark to 1 M⊙ in light green). The coloured lines represent the median values of ℋ for each of these mass samples. The typical on ℋ have been represented at the bottom left of the panel. They have been estimated by computing first quartiles around the median value. The large coloured dots mark the position of the maxima of these curves.

Our sample covers a slightly narrower range of Rossby numbers between 0.5 and 4, which gives us only a small overlap with the data from Wright et al. (2011). However, we clearly do not retrieve similar trends in that range of values. Indeed, the impact proxy ℋ exhibits instead a maximum between Ro = 0.7 and 1 (see the black line of Fig. 9). The difference between these two trends could be due to the difference in the active regions probed by each indicator. In the present work, we focus on the signature of spots, which are mostly correlated with the photosphere properties, while X-ray luminosity is more sensitive to coronal features. Moreover, it is also worth mentioning that we use a larger sample of stars than Wright et al. (2011), possibly allowing better mapping of the trends with the stellar parameters.

For the right panel of Fig. 9 we choose to represent the median value of ℋ for increasing intervals of masses from 0.5 to 1 solar masses. As illustrated in Fig. 6, ℋ exhibits distinct behaviours depending on both mass and rotation. The higher masses are not represented, as they show no particular structure. When considering the slowest rotators, a series of maxima appears (large dots). These maxima amplitudes increase with mass until it reaches its largest value for the mass interval M = 0.8 M_⊙. The value of the maxima then goes down. These maxima also depend on the rotation period. For low mass intervals, the maximum corresponds to a rather slow rotation, while for more massive stars in this subsample, it corresponds to smaller rotation periods. For stars with very low mass (< 0.8 M_⊙), a second regime of high ℋ seems to appear for the faster rotating stars. However, this second activity regime remains more uncertain due to the reduced amount of data that we have for low masses in this region.

On the other hand, the left panel of Fig. 9 can also be compared to the photometric activity index S_ph that was measured in a similar stellar sample by Santos et al. (2024); the result is plotted in Fig. 1 of this latter paper as a function of the rotation period and for different spectral types, which are also representative of different mass ranges. The results of Santos et al. (2024) show a shape similar to our results, when considering the global scatter plot. This is not surprising, as the impact proxy ℋ is defined as the mean level of variability due to spots in the considered LC, and it is thus similar to S_ph index. However, if you look at the trends by mass (or by spectral type in Santos et al. 2024), you do not see a maximum of activity as a function of rotation, just a plateau. This can be attributed to our choice of representation of our results, but also because our model in the Fourier domain allows us to retrieve additional information, namely the mean spot lifetime. This new index helped us to further characterise the stellar activity compared to these previous studies by identifying different activity regimes as a function of the Rossby number, clearly demonstrating the value of our new modelling.

5.2. Link with the magnetic field topology

Thanks to recent spectropolarimetric observations (e.g. ESPaDOn, NARVAL), the three-dimensional configuration of surface magnetic fields could be reconstructed for a few stars using Zeeman-Doppler techniques (e.g. Donati & Landstreet 2009; Donati et al. 2023). As seen in Fig. 3 of Donati & Landstreet (2009), different trends emerge, based on the iso-Rossby lines : stars with Ro > 1 have a weak magnetic field with a poloidal axisymmetric topology while stars with 0.1 ≲ Ro ≲ 1 have stronger, toroidal and non axisymmetric magnetic field. This sample is composed of about 15 stars with masses comparable to those in our stellar sample, that is, between about 0.5 and 1.2 M_⊙; stars with Ro ≲ 0.1 have generally very strong poloidal fields, but correspond to stars with mass under 0.5 M_⊙, which are outside the range of our study sample.

To put these observations in perspective with our study, we choose in Fig. 6 similar axes to Fig. 3 in Donati & Landstreet (2009). Our results only cover the stars from 0.4 to 1.4 solar masses, which correspond to the top of the Fig. from Donati & Landstreet (2009) According to this figure, and focusing on masses below about 1.2 M_⊙ as in Donati & Landstreet (2009), we see that stars with Ro < 0.7 (Regions 1) have short τ_life and intermediate spot impact ℋ for weak Rossby number, then reach a regime of intense activity with a large spot impact and long lifetimes (region 2). For stars with Ro > 1 (Region 3 and peakless group), we again observe a medium and even low impact spot for the highest Rossby numbers. Lifetimes are longer than those below Ro = 0.7.

We can therefore establish links between these proxies and the topologies described by Donati & Landstreet (2009). Firstly, short lifetimes seem to appear for toroidal and asymmetrical topologies, while long ones appear for more stable and poloidal configurations. Then, it is interesting to note that the transition at Ro ≈ 1 between the two kinds of observed topologies is near the high activity regime. Indeed, Ro ≈ 1 indicates a change of physical regime and the high observed activity index could be the signature of a high variability in the stellar magnetic structure (e.g. intermittent reversals) at the transition between two stability domains.

The relation depicted between the activity properties revealed in this work for a large sample of stars and the observed magnetic fields for a handle of close stars remains preliminary. In particular, the question of the variability over time and of the magnetic cycles is certainly of prime importance to perform relevant comparisons. Regarding magnetic cycles, See et al. (2016), Fabbian et al. (2017) have produced the same plot as Donati & Landstreet (2009), but have connected the magnetic field topologies to the magnetic cycles by looking at the active and inactive branches of Saar et al. (1992), Brandenburg et al. (1998). They have also showed the existence of a delimitation between the two branches at Ro = 1 between the long and short cycles, as it was observed in Saar et al. (1992).

5.3. Insights into the evolution of stellar dynamo

As seen in the previous section, different observations show that the magnetic properties of stars depend on the Rossby number. In particular, a transition seems to operate in the spot activity, field topology (Donati & Landstreet 2009) or magnetic cycles (Saar et al. 1992; See et al. 2016) around Ro = 1. In this work, we have shown that this transition is predominantly a maximum in the impact proxy, as well as the progressive emergence of long-lived spots as Ro increases and the decrease of the τ_life/P_rot ratio. These changes of regimes as a function of the Rossby number are also supported by 3D MHD numerical simulations that have already highlighted change of regimes with the variation of the Rossby number in the differential rotation (e.g. Brun & Browning 2017), dynamo processes (e.g. Raynaud et al. 2015; Zaire et al. 2022) and cycles (e.g. Brun et al. 2022).

Although good guides, current simulations are nevertheless still far from realistic regimes, and the constraints such as those brought in this work represent a good opportunity to test their theoretical outcomes. For instance, one could be tempted to draw a simple picture for the evolution of the dynamo mechanisms according to the simple trends observed as a function of the Rossby number. Indeed, owing to the efficient magnetic wind breaking in such low-mass stars, the Rossby number is expected to increase as they evolve on the main-sequence (e.g. Noraz et al. 2024). If we take into consideration the possible link between the Rossby number and the age of the stars.

A catalogue of Kepler stars with ages has been published by Lu et al. (2021), based on a gyro-chronological study. Fig. 10 shows the evolution of these ages for each intervals of Rossby defined on Fig. 6. This figure was made for different mass ranges, since the time they spend on the main-sequence are not the same depending on their mass. For stars with masses below one solar mass, the median values shows that the Rossby intervals do indeed increase with age. This suggests that a star may transition through the various types of activity outlined in Sect. 4.4. However, it is important to note that the distributions are quite broad and show significant overlap, particularly among stars with peakless spectra and those in Regions 2 and 3. It could be explained by an incorrect estimation of the ages, especially given that age estimation through gyrochronology is not a completely robust method yet, especially for high masses, or a bad sorting of the spectra with or without peaks. For stars with masses > 1 M_⊙, the position of the median value of the peakless group became lower than Regions 2 and 3 and the different regions are even more spread out.

Fig. 10. Distribution of star ages for increasing mass intervals. Each row corresponds to a decreasing mass interval (from 1 to 1.2 M⊙ for the first row, 0.9 to 1.0 M⊙ for the second row, 0.8 to 0.9 M⊙ for the third row and 0.6 to 0.8 M⊙ for the last row). The colours of the histograms correspond to the different regions of Figure 6. For each of these histograms, a line of the same colour is drawn to represent the median of the histogram.

In this way, Fig. 9 could suggest that the activity type evolves with time at a given mass. For low masses, stars are likely to evolve in a constant phase with moderate activity, followed by an increase of activity near Ro = 1 and then a rather rapid decrease. For stars with masses above 1 M_⊙, the stars potentially go from the activity of Region 1 (moderate ℋ and short τ_life) to low ℋ and longer τ_life, without the “bump” phase.

5.4. Limitations of the model: Influence of faculae and differential rotation

The method used in this paper is based on an analytical model detailed in Sect. 2 where two physical phenomena have been neglected: the differential rotation and the appearance of faculae on the stellar surface. The goal of this section is to discuss the validity of these assumptions regarding the results presented in this paper.

5.4.1. Influence of faculae

Faculae are bright magnetic regions that appear on the stellar surface, thus being able to contribute to photometry by inducing an increase in luminosity, unlike spots. They have been observed and characterised well on the Sun, but their impact and presence on other stars are still current topics of research.

Their effective contribution is debated. Basri (2018) explains that their impact might not be readily discernible in Kepler data if their influence is supposed to be similar to that of the solar faculae. However, the analysis by Shapiro et al. (2017) shows that the amplitude of the rotation peaks is diminished by faculae in the solar case, which may influence the estimate of ℋ but not the width of the rotation peaks. However, if we go back to the description of our model, the impact of faculae should mainly be seen in the transit-dependent part of the spectrum. Indeed, Toriumi et al. (2020) explains that the impact of faculae on brightness results in an increase in photometry at the edges of the transit function. In this case, it is the characteristic transit time and possibly ℋ that are affected. The reason why the rotation peaks of the Sun are very barely visible is because it has mainly spots with lifetime close to its rotation period.

Besides, Reinhold et al. (2019) studied the distinction between spot- and faculae-dominated stars using photometric and chromospheric data. This work showed that the faculae-dominated stars have Rossby number above 1 and P_rot > 20 days. This corresponds to the area where the stars with peakless spectra are dominant. This might seem contradictory to our previous argument. However, these stars exhibit the slowest rotation rates, making their rotational peaks more challenging to detect. The visibility of such peaks depends on the number of transits a spot can produce during its lifetime, and slower rotators naturally have fewer observable transits.

5.4.2. Influence of the differential rotation

The DR corresponds to a latitudinal shear between the pole and the equator. Since all the spots are not transiting at the same latitude, and therefore may have a different rotation period (Carrington 1859; Reinhold et al. 2013), the DR can widen rotation peaks and thus induce a bias on the spot lifetime proxy.

The difference with the model presented in Sect. 2.3 is that each spot is now associated with a Dirac comb, with a spacing between the values of the rotation period at the equator and at the pole. In the Fourier domain, we therefore have a superposition of comb-like signals shifted in frequency from each other by nΔΩ, with ΔΩ the value of the DR relatively to the equator at the latitude at which the spot has transited and n the number of the harmonics considered. This shift in frequency contributes to broadening the rotational peaks. Thus, the width of the n-th rotation peak is dominated by the spot lifetime if:

$$\begin{array}{c}\hfill n\mathrm{\Delta }\mathrm{\Omega }<\frac{2\pi }{{\tau }_{\mathrm{life}}}\end{array}$$(20)

Thus, the higher harmonics can be widened by the DR. Since no more than three rotation peaks have been observed in the Kepler data and that the algorithm fits only two of them, this potential bias is not distinguishable from noise.

Are there some stars where the DR can be observed in the first harmonics? Multiple DR estimation have been provided in the past years using different methods. Applying Doppler imaging, Barnes et al. (2005) found a DR that increases with the effective temperature. Combining these results with DR estimates through the broadening of the rotation peaks in the Fourier spectrum by Reiners (2006), Collier Cameron (2007) proposed a power-law fit of the form ΔΩ = 0.053(T_eff/5130)^8.6. Later, based on stellar models, Küker & Rüdiger (2011) showed that ΔΩ depends on two power laws, depending on the range of temperature of the stars. The broadening of rotation peaks in the power spectra have also been used to determine the DR in the Kepler data by Reinhold et al. (2013), Reinhold & Gizon (2015). However, this method has been questioned by Basri & Shah (2020), whose results show that DR is only visible for very long lifetimes (∼10 − 15P_rot). In this work, the temporal evolution of the spots is modelled by a symmetric and triangular function and their lifetime is defined as the total width of this function. Taking the width at half height as in our formulation, their result translates into a ratio τ_life/P_rot between 4 and 6 in our model.

Fig. 11 represents the estimation of the spot lifetime of this paper as a function of the effective temperature. We add the ΔΩ − T_eff relations from Collier Cameron (2007) and Küker & Rüdiger (2011). In these studies, the stars that do not respect the criteria of Eq. (20) are those with long lifetimes and mostly high temperatures and represent a very small fraction of the sample.

Fig. 11. Spot lifetime estimation expressed in radians per days as a function of the effective temperature in K. The orange line is the Teff − ΔΩ relation proposed by Collier Cameron (2007), with ΔΩ in days. The dark red line is the relation from Küker & Rüdiger (2011), with the dashed lines representing the two components of this relation. The hatched regions are the ones where Eq. (20) is not verify, for instance the DR may affect the lifetime estimation, τlife.

6. Conclusions

Stellar spots and their photometric signature represent very interesting probes of stellar magnetism for large samples of stars. However, finding information about spot properties in LCs is a complex process due to the intrinsic degeneracies of the problem. Indeed, the impact of a spot on an LC depends on many of parameters, such as its position, size, temperature, the rotation axis inclination of the star, and even its differential rotation. As a result, different stellar configurations and spot distributions and properties can lead to similar impacts on LCs in the temporal domain.

In this paper, we have proposed a new analytical model of the LCs PSD in the Fourier domain that allows some of the degeneracies to be removed while accounting for the main properties of spots. This model takes into account the rotation peaks of the LCs power spectrum and other components of the spectra at low frequencies. This LC power spectrum model introduces three proxies to estimate average spot properties: τ_tr, τ_life, and ℋ. This description shows that rotation peaks in the PSD occur when stars have spots with lifetimes exceeding the rotation period. Consequently, two spectral categories emerge: those dominated by spots with multiple transits that have apparent rotation peak harmonics and those from stars with a majority of spots that create a single transit in the LC, which thus have a peakless spectrum (as is the case for the Sun). After validation by simulations, we applied the model to a large Kepler dataset consisting of the catalogues of McQuillan et al. (2014) and Santos et al. (2019, 2021).

Our data analysis yielded several results. We showed a clear distinction in the spot activity according to the stellar mass in the Kepler data. Stars above one solar mass have a much lower spot impact than low-mass stars. For stars between 1 and 1.2 M_⊙, lifetimes also increase with Rossby number. For stars below 1 M_⊙, the proxies adopt different values depending on the Rossby number. In particular, an intense activity regime was observed for stars between Ro = 0.7 and 1. Spot lifetimes also tend to increase with the star’s rotation period. The increase, however, seems to stop or even decrease a little after Ro = 1. Since the peakless spectra are mainly present for large Rossby values, the evolution of spot lifetimes is difficult to infer for these values. For ℋ, a decrease is noted beyond Ro = 1. These variations indicate that the evolution of spot activity with Rossby number is not as linear as generally suggested in the literature.

The various types of activity described in terms of Rossby number indicate a potential connection to distinct dynamo regimes. They can further be correlated with the topology map constructed by Donati & Landstreet (2009). Moreover, stars in the high activity regime have Rossby numbers close to one. Previous studies have indicated behavioural changes near this Rossby number, such as alterations in magnetic cycles and topology. This makes the observed high spot activity near this delineation particularly interesting. This new methodology itself has the potential for further development, including the possibility of studying the impact of DR or using transit time estimates, which may be related to a star’s inclination or the appearance of faculae.

We conclude that these results offer promising insights that could shed light on the inner dynamo of stars. We plan to correlate these results with simulations of stellar interiors in order to gain a deeper understanding of the three distinct regimes we observed. A comparison with the work of Basri et al. (2022) on the spot lifetimes is planned. We also plan to put our results into perspective with the LAMOST activity observations. In the longer term, these results could aid in the interpretation of PLATO data, particularly for the fitting algorithm that might be used on future LC power spectra. They would also help in better constraining the impact of stellar spots on photometry and the information they provide about stellar magnetism.

Acknowledgments

This paper includes data collected by the Kepler mission and obtained from the MAST data archive at the Space Telescope Science Institute(STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). We warmly thank Cilia Damiani to provide us the simulation code to create LC. We also thank K. Belkacem, A. S. Brun, N. Meunier, P. Boumier and L. M. R. Lock for useful discussions. We also warmly thank the anonymous referee for all his comments and thoughtfulness, as they have enabled us to significantly improve our article.

References

Appendix A: Transit time in the case of the multiple-transit spots simulations

Fig. A.1 shows the results obtained for the mean transit time on LCs simulated with multiple-transit spots only. The fitted value is underestimated and much more widely dispersed for our model than for the Harvey profile, which was expected. In this case, there are only spots with multiple transits and therefore no contribution from those with only a single transit, making it impossible to observe this proxy. More precisely, the new model tries to fit a S_single contribution while it is not present in the simulated data for only multiple-transit spots. This component converges towards values close to P_rot/8, which correspond to the third rotation peaks counting from the fundamental, since only the first two are fitted by the rotation peaks component. In fact, Fig. A.3 shows that the single-transit spot component converges on this third peak.

Fig. A.1. Same as in Fig. 5 but for the simulations with only multiple transits spots.

For the Harvey profile, the median value of the distribution is close to the expected one. It is also an expected result: because there is many rotation peaks, the fit is driven by their amplitude (see Fig. A.2), which is linked to the transit profile (see Fig. 1). There are also several of accumulation points on higher values, for the same reason that for the mixed-spots cases: in some cases, when only one or two rotational harmonics are observable, the fit with the Harvey profile underestimates the width at half maximum.

Fig. A.2. Power spectrum of a simulated LC with only multiple transit spots. The spectrum has been fitted with the Harvey model (blue).

Fig. A.3. Same as Fig. A.2 but fitted with the new model describe in this article.

Appendix B: Technical description of data cleaning

As a first step and to simplify the methodology, the Harvey profile and our model were applied to all the data. The classification of spectra into those with or without peaks was performed subsequently, as the algorithm for this distinction was developed at a later stage of the study.

Fig. B.1 shows the distribution of the relative error estimates obtained using the method of Toutain & Appourchaux (1994) for the three parameters to fit in the new model. In the case of τ_life, there are two accumulations, one around 0.5% and one around 2.5%. As observed with the simulations, this second peak associated with relatively large uncertainty corresponds to the single transit spots case where the rotation peaks are reduced or not apparent, and where the fit has more difficulty converging. To avoid spurious estimates, we choose to remove the data with a relative error higher than 10% for ℋ and τ_life, and 20% for τ_tr, which are not essential for the present study. By adding these criteria, 8.5% and 18.2% of the McQuillan’s and Santos’ sets respectively are removed. We also remove some incorrect fit results by visual inspection of dubious points accumulation in theresults (more than 10 000 stars). The final set after this cleaning consist of approximately 27 000 stars, which remains a significant sample, and reassures us in the ability of this study to extract global trends.

Fig. B.1. Histograms of the relative errors of the fitting algorithm for each proxy: ℋ (top), τtr (middle) and τlife (bottom) for the McQuillan sample.

For the Harvey profile, the same criteria have been applied to ℋ and for the characteristic time; we kept the value with a relative error under 20%. Finally, 5% have been removed from the Santos data and 2.6% from the McQuillan sample, leading to a total cleaned sample of 41 176 stars.

The reason we obtain convergence of the parameters on less data with our model is that it depends on 7 free parameters against 4 for the Harvey profile, which makes the algorithm less robust.

After sorting the spectra, we identified 15 092 stars with peakless spectra and 18 663 with rotation peaks, for a total sample of 33 755 stars.

Appendix C: Transit proxy as a sanity check

The top panel of Fig. C.1 shows the results for the transit proxy as a function of the rotation period estimated by McQuillan et al. (2014). There is a first linear distribution around the relation, τ_tr = P_rot/4, which corresponds to the prediction of the model detailed in Section 2. This distribution is also clearly shifted downwards.

The accumulation of points around the red dotted line consist of a set of data for which the algorithm has not converged well. For example, in some cases there are more than two rotation peaks or other artefacts at higher frequencies, which may bias the estimation of the transit time and lead to an underestimation. In other cases, the first harmonic has very low amplitude, and is not taken into account by the fitting algorithm, which biases the value, as predicted by the simulations. Another bias can be an incorrect estimation of the rotation period by McQuillan et al. (2014) or Santos et al. (2019) by a factor of two, which can occurs with the methods used in these papers. Since this τ_tr is not a parameter of interest in this paper, these values have been retained.

The top panel of Fig. C.1 shows the same figure but with the times estimate by the Harvey profile. The results are significantly the same, which reinforces the idea that this parameter is close to the transit time.

Fig. C.1. Top panel: Transit proxy τtr as a function of the rotation period, estimated for the whole sample of stars. The red line shows the relation τtr = Prot/4, and the dotted red line its harmonic τtr = Prot/8. Bottom panel: Harvey characteristic time as a function of the rotation period.

All Tables

All Figures

Fig. 1. Illustration of the analytical description of the spot transit proposed in this paper. The top panel shows the model for the single-spot case. The model is shown in the time domain on the left, and in the Fourier domain on the right. Each time, the distinction is made between the single-transit-spot/peakless case (orange) and the multiple-transit-spot/with peaks case (blue). The bottom of the figure illustrates the general case with both multiple-transit and single-transit spots’ contributions.

In the text

Fig. 2. Left panel: Example of a power spectra density of a Kepler LC (grey) ‘with peaks’ fitted by the model detailed in Sect.2. The transit component due to the single-transit spots is dash-dotted orange line, the multiple-transit spots component in a dashed light green line and the total fit in blue. The rotation peaks are marked by purple dashed lines. Right panel: Example of a peakless power spectra density of a star where the Harvey profile is used for the fit.

In the text

Fig. 3. Median height of the spectrum (ℋ) of each set of simulations as a function of |Atotτlife|2 (see Eq. (16)). The colour of the dots indicates the rotation period of the star (10 days in green, 20 days in red and 30 days in purple), while the shape gives the type of spot lifetime (diamonds for the peakless stars, hexagons for mixed spots and 5-pointed stars for stars with rotation peaks spectra). The size is proportional to the surface of the spots. The large hexagon on the top-right corresponds to set 19 with spots radii of 10 degrees. The error bars are based on the first and third quartile. The points with a shading towards white are the results for the 3 simulations where the inclination varies (sets 20, 21 and 22), white being for 0° inclination (set 22).

In the text

Fig. 4. Normalised histograms of τlife. Each column corresponds to a different rotation period, and each row to a type of lifetime (single-transit, multi-transit and mixed spots). Here, only the results for the new model are shown, since it cannot be retrieved using the Harvey model. For each panel, the red lines show the expected value. The vertical black dash-dotted line represents the lifetime value for the single-transit spots and the dashed the one for the spots that last several transits (i.e. τlife = 1Prot). In the case of mixed spots, these two lines are represented, as they show the range of possible values. This space has been coloured in blue for clarity. The full line represent the mean value of the lifetimes. For single-transit and mixed-spots simulations, the vertical black dotted line represents the value of lifetime where two transits are visible in the LC. The stepped histograms in green, blue and orange correspond to the simulations where the inclination of the rotation axis varies (respectively to 0, 30 and 60 degrees).

In the text

	Fig. 5. Same as Figure 4 but for the transit proxy τtr. The histograms in blue correspond to the value found by our model, and the ones in orange to the value found by the model of Harvey. The vertical black dashed lines indicate the expected value for each set of simulations, which corresponds to Prot/4. The coloured dashed lines represent the median value of each distribution.
In the text

Fig. 6. Top panel: Spot-impact proxy for the stellar sample as a function of mass and rotation period. The colour scale shows the impact proxy ℋ with a logarithmic scale. Yellow is for a large spot impact (active stars) and black for low spot impact (less active stars). The black lines represent the iso-lines of the Rossby number. The ⊙ point indicates the location of the Sun in this diagram. Middle panel: Same figure with contour lines showing the locations of the stars with ‘peakless’ or ‘with peaks’ spectra, respectively in orange and blue. The grey dots represent the whole sample. Bottom panel: Same as the top panel but with the lifetime proxy τlife. Only stars with spectra with peaks are presented.

In the text

	Fig. 7. Estimated mean lifetime τlife as a function of the rotation period for stars with spectra with rotation peaks.
In the text

Fig. 8. Histograms for the three different regions and the peakless sample as a function of different stellar parameters. In blue: Sample of stars with rotation peaks. The different shades of blue correspond to the three regions of Fig. 6 cut according to the Rossby number. This slicing is illustrated in panel (a) as a function of the Rossby number. The orange histograms represent the peakless group. From top to bottom, the histograms are for (a) the Rossby Number, (b) spot impact proxy, (c) the mean lifetime τlife,(d) the rotation period, and finally (e) the ratio Prot/τlife.

In the text

Fig. 9. Scatter plot of the total sample of stars (grey dots in the background) as a function of the spot impact proxy and the rotation Rossby number (left) or period (right). On the left panel, the black line represents the median value of ℋ as a function of Ro for the whole sample. The hatched area around this line represents the interval of the quartiles. The dotted lines show the position of Ro = 0.7 and 1. For the right panel, these stars are grouped into slices of increasing masses (from 0.5 M⊙ in dark to 1 M⊙ in light green). The coloured lines represent the median values of ℋ for each of these mass samples. The typical on ℋ have been represented at the bottom left of the panel. They have been estimated by computing first quartiles around the median value. The large coloured dots mark the position of the maxima of these curves.

In the text

Fig. 10. Distribution of star ages for increasing mass intervals. Each row corresponds to a decreasing mass interval (from 1 to 1.2 M⊙ for the first row, 0.9 to 1.0 M⊙ for the second row, 0.8 to 0.9 M⊙ for the third row and 0.6 to 0.8 M⊙ for the last row). The colours of the histograms correspond to the different regions of Figure 6. For each of these histograms, a line of the same colour is drawn to represent the median of the histogram.

In the text

Fig. 11. Spot lifetime estimation expressed in radians per days as a function of the effective temperature in K. The orange line is the Teff − ΔΩ relation proposed by Collier Cameron (2007), with ΔΩ in days. The dark red line is the relation from Küker & Rüdiger (2011), with the dashed lines representing the two components of this relation. The hatched regions are the ones where Eq. (20) is not verify, for instance the DR may affect the lifetime estimation, τlife.

In the text

	Fig. A.1. Same as in Fig. 5 but for the simulations with only multiple transits spots.
In the text

	Fig. A.2. Power spectrum of a simulated LC with only multiple transit spots. The spectrum has been fitted with the Harvey model (blue).
In the text

	Fig. A.3. Same as Fig. A.2 but fitted with the new model describe in this article.
In the text

	Fig. B.1. Histograms of the relative errors of the fitting algorithm for each proxy: ℋ (top), τtr (middle) and τlife (bottom) for the McQuillan sample.
In the text

	Fig. C.1. Top panel: Transit proxy τtr as a function of the rotation period, estimated for the whole sample of stars. The red line shows the relation τtr = Prot/4, and the dotted red line its harmonic τtr = Prot/8. Bottom panel: Harvey characteristic time as a function of the rotation period.
In the text