© ESO 2018

1. Introduction

Solar activity is characterized by an 11 year cycle of the number and area of sunspots (Schwabe 1844). Its monitoring is important in many fields, such as Earth climate (Haigh 2007 or space travel studies (Pulkkinen 2007). Sunspots are regions of high concentration of the solar magnetic field (Solanki 2003), indicating that the field is the main driver of the activity cycle. In order to understand the solar magnetic behaviour, dynamo models have been developed (Charbonneau 2010). Their aim is to explain how an initially weak magnetic field can be amplified to the values observed in the Sun. A traditional observational test to which these models ought to comply is to reproduce the butterfly diagram which describes the evolution of the latitude of the solar active regions with time (Maunder 1904). Activity has also been observed in other stars. It is important in order to understand stellar evolution, for instance through magnetic braking (Thompson et al. 2003), and exoplanet habitability (Vidotto et al. 2013).

Recovering stellar butterfly diagrams is however a difficult task that requires locating individual spots or groups of spots on the stellar surface. Spot mapping using photometric data alone is known to be hampered by degeneracies in light curve models (e.g. Walkowicz et al. 2013), so that spectroscopic or interferometric data are usually favoured in order to recover stellar brightness maps (Vogt & Penrod 1983; Roettenbacher et al. 2016). A number of active targets have been monitored through long-term Doppler mapping campaigns, although due to sparse temporal sampling these time series have rarely led to actual butterfly diagrams (e.g. Hackman et al. 2011, for II Peg). However, Doppler mapping or long baseline interferometric imaging can only detect very large stellar spots, and are therefore limited to the most active stars, quite far from the Sun in terms of magnetic properties. We propose here an original method that uses Kepler time series to measure the latitudes of the active regions of stars and construct a stellar butterfly diagram. We apply it to the sun-like star HD 173701 (KIC 8006161).

2. Method

Our approach uses Kepler photometric time series (Jenkins et al. 2010) and is divided into two stages. First we wanted to obtain information on the large-scale rotational flow in the stellar interior, in particular in the convective envelope below the surface, using the information contained in the oscillation frequencies of the global acoustic pulsation modes (p modes) of the star (Sect. 2.1). To that end, we used the recent results of Benomar et al. (2018). Second, we tried to measure the rotation rates of active regions, carried by the surface rotational flow (Sect. 2.2). Third, we constructed the butterfly diagram of HD 173701 over the duration of the Kepler mission by inverting the rotation-rate measurements thanks to the asteroseismically derived rotation profile (Sect. 2.3).

We note that these measurements are uncorrelated in Sun-like stars since the observed pulsation modes of the stars are only affected by large-scale flows. In the frequency space, these two contributions are well separated, the modes having characteristic frequencies of the order of a few thousand μHz, while the rotation rates inferred from the starspots are of the order of a few hundred nHz.

2.1. Asteroseismic rotation profile

Our starting point is the derivation of a rotational profile of a Sun-like star. It is well established that the solar radiative interior rotates as a solid body, while the rotation rate of its outer convective envelope varies with latitude (Thompson et al. 2003). For the Sun, the equator rotates faster (∼476 nHz) than the pole (∼320 nHz). This differential rotation is most probably the result of an interplay between rotation and convection in the solar envelope (Rüdiger 1974; Gilman & Glatzmaier 1981).

Rotation affects stellar pulsations. The classical framework to describe stellar oscillations is built on the perturbed equations for the static stellar structure (Aerts et al. 2010). It can be shown that superimposing a small-amplitude rotational velocity field on the hydrostatic stellar structure will lift a degeneracy of the eigenfrequencies of the p modes, an effect usually termed “rotational splitting” (Lynden-Bell & Ostriker 1967). The degenerate frequency then becomes a multiplet whose distribution depends on the solid-body component of rotation rate of the star and on the magnitude of the differential rotation in the convective zone. If one defines a theoretical model Ω = Ω(θ) for the latitudinal differential rotation, with θ being the co-latitude, it is possible to relate it to the frequency splitting. This means that one can infer a differential-rotation profile if there are precise enough measurements of the frequency splitting.

In a recent work, Benomar et al. (2018) measured contributions from latitudinal differential rotation to the above-mentioned rotational splitting for a set of Sun-like stars. It includes the solar analogue HD 173701, which has a mass M = 0.95 M_⊙, with M_⊙ the solar mass, and an age t_⋆ = 4.49 Gyr (Silva Aguirre et al. 2017). This star rotates on average 1.33 times faster than the Sun, with a measured bulk rotation rate of 566 nHz. It has been observed for approximately four years by Kepler. Detailed modelling of its acoustic power spectrum led to a significant detection of non-zero latitudinal differential rotation. In the following we explain how to extend such results with rotation-rate measurements in order to construct the butterfly diagram of the HD 173701. For the sake of completeness we review the work of Benomar et al. (2018) in more detail in Appendix A.

At this point we have to make a further assumption relative to Benomar et al. (2018). It has been shown by Kiefer et al. (2017), Santos et al. (2018) and Santos et al. (2018) that the oscillation frequencies of HD 173701 shift with the activity cycle. We therefore need to take these shifts to be constant over frequency ranges of the order of a rotational splitting.

2.2. Photometric characterization of active regions

Once obtained, a relation Ω(θ) can be used to infer the co-latitude, θ_a, of a stellar active region such as a group of spots provided its rotation rate, Ω_a = Ω(θ_a), can be determined. The method for rotation rate measurements is inspired by solar observations and uses photometry. The total solar irradiance varies during the activity cycle due to bright regions (plages) and dark regions (spots) moving across the solar disk (Fröhlich & Lean 2004). A similar behaviour is observed for other stars using photometric measurements. During its lifetime, an active region will produce a quasi-periodic signal in the photometric time series. The long-cadence Kepler time series allow us to measure such modulations for HD 173701.

To measure the period of this modulation, we model the time series using Gaussian processes. Owing to the intermittent influence of the plages and spots, and to the stellar limb darkening, this signal departs significantly from pure harmonic oscillations. In addition, the noise impacting the photometric times series has a frequency-dependent part mainly caused by surface granulation. The time series model depends on several parameters which are estimated using an MCMC algorithm. These parameters are the amplitude of the modulation induced by the active region A (in ppm), the lifetime of the active regions τ_d (in days), the rotation rate Ω_a (in nHz), and an additional noise component α (in ppm), describing residuals not captured by the other component of the Gaussian process (see Appendix B for details).

We model independently the time series for each Kepler quarter which span from Q0 to Q17. The fit obtained for Q3, representative of our results, is shown in Fig. 1; it reproduces well the time series and its frequency spectrum. Our analysis is restricted to a region in the frequency spectrum surrounding the low-frequency activity peak (Fig. 1, middle panel). In fact, we observe multiple peaks potentially corresponding to several active regions. We consider that each of these peaks corresponds to a spot or a group of spots that rotates at the same latitude. We treat them as a single active region and consider the rotation rate of its barycentre. This is motivated by the fact that, when counting sunspots, a larger weight is given to groups (Hathaway 2010). The independent modelling, quarter by quarter, allows us to measure the evolution of the rotation rate with time. We note that quarters Q0, Q1, and Q17 were left out because the corresponding time series are too short to obtain robust results.

Fig. 1. Upper panel: time series for the relative flux HD 173701 during Q3. The black dots mark the measurements and associated errors. The red line shows the inferred mean value of the posterior predictive density conditional on the inferred MAP of the parameters (Gelman 2004). Middle panel: corresponding power spectrum computed using the Lomb-Scargle periodogram (Scargle 1982). The black line represents the observed power spectrum. The red line is the power spectrum of the mean value of the conditional posterior predictive density. Lower panel: marginalized 1D and 2D PDFs for the parameters of the correlation function of the Gaussian process used to model the time series.

The Gaussian processes reproduce the temporal variation of the intensity curve very well as can be seen in the upper panel of Fig. 1. Using the marginal density of the rotation rate Ω_a, we can estimate its maximum a posteriori (MAP). Such a density is seen for Q3 in the lower panel Fig. 1, corresponding to a MAP of 704 nHz. The MAP estimates of the rotation rates for all quarters are in the range 330–985 nHz, which already indicates that the spots are migrating along the stellar surface. In general the amplitudes of the signal and the lifetimes of the active regions are poorly constrained. Since they are strongly correlated, a wide range of values for the lifetime and the amplitude can lead to good models. Consequently, their estimated values cover several orders of magnitude. This does not impact the final result since these parameters are not correlated with the others.

2.3. Inversion for the latitude of the active regions

Our goal is now to invert for the latitude of the active region at each for each quarter we selected. Since the latitudinal rotation profile has been derived using the entire Kepler time series, it is necessary to assume that the properties of differential rotation do not vary on timescales comparable to the activity cycle. In the solar case, this is verified to a very good approximation (Thompson et al. 2003).

In order to invert for the latitude of the active regions, we need to take into account all the errors that may affect our final estimate of θ_a(t). These are of two kinds. The first is the error on the measurement of the rotation rate, the second is the uncertainty in the theoretical model for the rotation rate, Ω(θ), as obtained by (Benomar et al. 2018, see also Appendix A).

A generic framework to solve an inverse problem that takes into account these two sources of error is given by the concept of conjunction of states of information (Tarantola & Valette 1982, see Appendix C). In the following we use the more compact notation Ω_a(t_i)≔ ≔_i and θ(t_i) ≔θ_i, with i = 1,. . .,N and N the number of rotation rate measurements. The posterior density for (Ω_i, θ_i), which represents the state of information (or state of knowledge, see Appendix C) we have once all the errors have been considered, can be written

$${\sigma }_i({\Omega }_i,{\theta }_i)=\frac{{\rho }_i({\Omega }_i,{\theta }_i)\Theta ({\Omega }_i,{\theta }_i)}{\mu ({\Omega }_i,{\theta }_i)}\cdot$$(1)

The theoretical density Θ(Ω_i, θ_i) describes the error on the model. It is basically the relation Ω_i = Ω_i(θ_i) with an associated error bar. The prior density ρ_i(Ω_i, θ_i) is the combination of the information obtained from the observation on Ω_i (see Fig. 1, lower panel) and the prior information on the parameter θ_i (which has been chosen uniform over the range 0°−90°). The posterior density σ_i(Ω_i, θ_i) results from the conjunction of these two pieces of information. In Fig. 2 we show these three densities for quarter Q3, respectively in the left, central, and right panel. The marginal posterior on θ_i is shown in the lower right panel, and is obtained from the integration over all possible rotation rates of the joint posterior. The meaning of the homogeneous density μ(Ω_i, θ_i) is given in Appendix C.

Fig. 2. Probability density for the latitude of an active region at median time of Q3. Left panel: theoretical density for the couple (Ωa, θa). Central panel: prior density for the couple (Ωa, θa); it is Gaussian for Ωa and uniform for θa. Upper right panel: posterior density for the parameters. Lower right panel: marginal density for θa, obtained after integration over Ωa.

3. Result

We use the marginal posteriors of θ_a(t) to derive the butterfly diagram pictured in Fig. 3 for all selected quarters. The latitude of the active region is restricted to the range [0°−90°] since we do not resolve the stellar disk and thus do not know on which hemisphere the spot is located.

Fig. 3. Asteroseismic butterfly diagram of HD 173701, giving the latitude of the active region as a function of the median time of each quarter (labelled on the upper axis). The grey shades show the posterior densities obtained for each of the selected quarters and normalized to the global maximum; greyscale from 0 (white) to 1 (black). Red dots mark the maximum a posteriori estimates of the latitude. The vertical bars give the limits of the corresponding 68.3% credible intervals. The green squares show the maximum of the density of the 1024 estimates obtained from artificial time series. Lower panel: variation in amplitude of the signal, SQ, normalized to its maximum.

The credible intervals were computed as the smallest intervals that encompass the MAP and over which the probability mass is 0.683. There is a clear variation in the data. Half the active regions are found at the equator. Among the other half, five clearly exclude the 0° latitude at a 68.3% level. The active regions seen in Q6 and Q14 exclude it at a 99.7% level. It is thus ruled out that the signal seen could be due to equatorial active regions, with the high-latitude regions being outliers driven by intrinsic stochastic variability.

The estimated co-latitude (red dots in Fig. 3) rely on the parameters of the Gaussian process, which were estimated from the time series. This means that different data sets of the same star in the same activity configuration would lead to different inferred latitudes. In order to investigate the variability of the angle estimation procedure (and hence check for the robustness of these results), we made Monte Carlo simulations. We simulated 1024 time series using the MAP estimates of the parameters of the Gaussian processes for each quarter¹. We applied the latitude-estimation process to each artificial time series, resulting in a sample of 1024 MAP estimates for each such latitude of the active region. The distribution of this MAP obtained from artificial time series can be compared to the MAP obtained directly from the actual data. The location of the maxima for these distributions are shown in Fig. 3 as green squares. The MAP locations are essentially not impacted by the noise realization, which supports the fact that the co-latitude extracted from the Kepler data for this star must be close to the values shown in Fig. 3. Slightly larger deviations are seen when the uncertainty on θ_a is large, which is to be expected.

Over the four years of measurement, the active regions remain located at latitudes below 50°. The butterfly diagram clearly shows an alternance between equatorial and high-latitude active regions.

In the lower panel of Fig. 3 we plotted S _Q, which we define as the standard deviation of the light curve evaluated over quarter Q. We used this quantity as a proxy for the activity index S _ph rather than the definition given in García et al. (2014) for two reasons. First it involves an evaluation over a period that is five times the rotation rate, which becomes hard to interpret when we take into account the variability of the rotation rate. Second, if we retain the average rotation rate of KIC 8006161, then a quarter is approximately 4.5 times larger, which is commensurate with the definition of S_ph and preserves the statistical independence of the points in our time series.

We first see that low-amplitude modulations are observed for active regions at high latitude and near the equator. Let us assume that the amplitude of the signal follows the variations of the fraction of the stellar surface covered by the active regions. Then this observation in qualitative agreement with what is already known for the Sun (Hathaway 2010, Fig. 28) where small spots are observed at all latitudes (within the band in which they are confined). The two highest-amplitude modulations are seen for the regions with the highest latitudes, during Q6 and Q14. If our assumption that the amplitude relates to the surface covered by the active regions, then this departs from what is seen in the solar case. There the largest spots are seen at mid-latitudes.

Extracting a periodicity for the activity cycle from these measurements remains a difficult task. The periodogram of the signal does not show any significant peaks at low frequencies. Time series with longer spans would be needed to obtain such an estimate. It is interesting to note that the two highest-latitude spots also correspond to the longest estimated lifetimes (in the sense of the MAP), with τ_d ≲ 1 yr. This is of course an overestimation, since the signal lasts less than this characteristic time (see Appendix B). It should be noted that the lifetime of the active region and its rotation rate are uncorrelated, we can therefore safely assume that these biases do not impact the estimates of the latitude for these epochs (as can be checked in Fig. 1 for Q3). Such values nevertheless indicate that active regions are producing a coherent signal over a large fraction of Q6 and Q14. This could correspond to one or several long-lived large spots. Again this seems to depart from what is known in the Sun, for which spot lifetimes are much smaller.

Other activity measurements have been obtained for HD 173701. Karo et al. (2018) have obtained Ca II H and K line measurements that can track the activity cycle. Kiefer et al. (2017), Santos et al. (2018), and Salabert et al. (2018) have all measured consistent values for the activity-induced frequency shifts over the last eleven quarters of the of the Kepler mission. The Salabert et al. measurements correlate well with S_ph. Further work and larger data sets are needed to understand precisely how the butterfly diagram correlates with these other indicators.

4. Conclusion

In this letter we presented the first stellar butterfly diagram derived obtained by combining information inferred from asteroseismic and photometric analyses. Provided an approximate stellar model is known, the only data required to perform the inversion is a photometric times series (collected from Kepler in this study). This approach identifies a powerful link between asteroseismology and other branches of stellar physics studying stellar magnetism, for instance Zeeman-Doppler Imaging (Semel 1989). Stellar butterfly diagrams, out of reach of Zeeman-Doppler mapping when it is applied to solar analogues, nicely complement the long-term monitoring of large-scale magnetic geometries still accessible for low activity stars. Both proxies of stellar cycles offer complementary views to understand the underlying dynamo processes. The technique itself requires only moderate computational time and can be envisaged as a systematic processing for surveys of starspots in the perspective of the forthcoming space mission PLATO.

Acknowledgments

This material is based upon work supported by the NYU Abu Dhabi Institute under grant G1502. Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant DNRF106). In memory of Albert Tarantola.

Appendix A: Inference of the differential rotation profile

In this section we describe briefly the principles of asteroseismic inversion that allow us to estimate the parameters of our model for the stellar rotation rate Ω(θ). In generic form it can be written as an expansion over a basis of functions W_s(θ) depending on even powers of cos(θ) (Brown et al. 1989)

$$\Omega (r,\theta )={\displaystyle \sum _{s=0}^{s_{\text{max}}}{\Omega }_s}(r)W_s(\theta ).$$(A.1)

The functions Ω_s(r) are often chosen as piecewise continuous functions in order to account for the change in rotational regime between the radiative interior and the convective envelope. In the case of HD 173701, the Ω_s(r) are chosen constants, and these are the parameters we estimate.

The rotational splitting can also be expressed as a basis expansion of the form (Brown et al. 1989)

$$\delta {\nu }_{n,l,m}={\displaystyle \sum _{j=0}^{j_{\text{max}}}{a}_{2j+1}}(n,l){\zeta }_j^{(l)}(m),$$(A.2)

where the ${\zeta }_j^{(l)}(m)$ form an orthogonal basis obeying ${\displaystyle \sum _m{\zeta }_i^{(l)}}(m){\zeta }_j^{(l)}(m)=0$ if i ≠ j. Finally, the a_j and Ω_s are related through the integral equation

$$\delta {\nu }_{n,l,m}={\displaystyle \underset{0}{\overset{\pi }{\int }}{\displaystyle \underset{0}{\overset{R_{\star }}{\int }}{K}_{n,l,m}}}(r,\theta )\Omega (r,\theta )r\text{d}\theta \text{d}r,$$(A.3)

where K_n,l,m(r, θ) is a kernel that depends on the equilibrium stellar structure and the eigenfunction of the corresponding p modes (Hansen et al. 1977). To obtain these kernels we numerically solved the equations for the stellar structure and pulsations, using, respectively, ASTEC (Christensen-Dalsgaard 2008a) and adipls (Christensen-Dalsgaard 2008b).

Fig. A.1. Top panel: power spectrum of the short-cadence Kepler data in the region of the eigenfrequency multiplet centred around the frequency νn=24,l=2,m=0. The black and red lines show respectively the data and the best-fit model of the power spectrum. The vertical red ticks mark the frequencies of the eigenmodes of the multiplet, for −2 ≦ m ≦ 2 . Bottom panel: splitting diagram for the multiplet. The first two terms of the splitting sequence correspond to the rotational effects, while β(ν) is an additional term that describes aspherical contributions to the eigenfrequencies. The red horizontal ticks correspond to those seen in the bottom panel.

A judicious choice of the basis functions ${\zeta }_j^{(l)}$ will ensure that there is a one-to-one relation between the coefficients a_2s+1 of the splitting expansion and the Ω_s of the rotation rate expansion (Ritzwoller & Lavely 1991; Schou et al. 1994)). The resulting orthogonality condition on the W_s allows us to derive their functional form (Pijpers 1997), which is $W_s(\theta )=P_{2s+1}^1(\cos \theta )/\sin \theta$, with $P_l^m(\cos \theta )$ being the associated Legendre polynomial of degree l and order m. Using this expression and setting s_max = 1 in Eq. (A.1), we obtain the formula used for the rotation rate in this study:

$$\Omega (\theta )={\Omega }_0-1.5{\Omega }_1(5{\cos }^2\theta -1)$$(A.4)

We note that the condition on s_max is justified by the above-mentioned one-to-one relation and the fact that the current seismic data only allow us to observe a₁ and a₃ (tests to detect a₅ splitting coefficients in the data considered in this study were inconclusive). The goal is to obtain estimates of Ω₀ and Ω₁ in order to derive the theoretical density of the couple (Ω_i, θ_i).

In practice, the coefficients a₁ and a₃ were estimated using a parametric model for the power spectrum (Gizon & Solanki 2003), namely a sum of Lorentzian components and some noise. The locations of the Lorentzians are the mode eigenfrequencies. Their distribution is given by the relation

$${\nu }_{n,l,m}={\nu }_{n,l,0}+\delta {\nu }_{n,l,m}+\beta (\nu ).$$(A.5)

Here δν_n,l,m describes the rotational effects and is given by Eq. (A.2). The additional term β(ν) has been introduced to account for the effects of departures from strict sphericity of the star (centrifugal force, magnetic fields, tidal distortion, etc.) and is a linear function of the central frequency of the multiplet (Gough & Thompson 1990). To estimate a₃, the mode of degree l = 2 was used.

The existence of a one-to-one relation between (a₁, a₃) and (Ω₀, Ω₁) was then used to estimate the latter. The measured rotational splitting in the observed acoustic spectrum was modelled (Benomar et al. 2018) using Eq. (A.2) with j_max = s_max = 1. The splitting coefficients were considered as free parameters to be estimated. The inferred values are a₁ = 563 ± 69 nHz and a₃ = 28.61 ± 12.41 nHz. An example of the modelled power spectrum is given in Fig. A.1, alongside the measured splittings. A model for the rotation rate of the convective zone of the form (A.4) was used together with kernels computed from stellar models (Christensen-Dalsgaard 2008a,b; Lund et al. 2014) to invert this model. Estimates for the rotation coefficients Ω₀ = 566 ± 69 nHz and Ω₁ = −104 ± 45 nHz were inferred.

Finally, we expand slightly on the assumption made concerning the influence of the activity-induced frequency shifts on the estimated value for a₃. The assumption we have to make is that the entire multiplet caused by rotational splitting is shifted as a block and thus the frequency splitting remains the same. This allows the activity-induced frequency shifts to depend on the frequency itself, but this variation has to remain negligible at scales of the order of ∼3 μHz, which is characteristic of a rotational splitting in KIC 8006161. If this holds, then we may expect a₃ to be relatively constant with time. We note that this is related to but subtly different from the assumption of non-varying differential rotation magnitude with time made at the beginning of Sect. 2.3, which says that the magnitude of differential rotation does not change over time. The new assumption says that the measurement of this constant quantity a₃ is not biased by time-varying effects affecting the frequencies.

Appendix B: Measuring rotation from active regions

One element needed to infer the latitude of an active region is a measurement of its rotation rate. It is critical to obtain an estimate of the associated uncertainties in order to compute the conjunction of information states described in Sect. 2 (see also Appendix C).

The idea retained here is to fit the low-frequency components of the time series using a Gaussian process. This method has been applied recently to Kepler data with some success (Angus et al. 2018). The critical point that allowed the use of Gaussian processes to model long time series was the development of methods for the fast inversion of covariance matrices² (Foreman-Mackey et al. 2017).

The idea of fitting a Gaussian process to a time series with N points is to consider it as the realization of a random vector of dimension N. In principle one would then have to estimate the mean and the covariance matrix of the parent distribution, i.e. N(N + 3)/2 parameters. However, the assumption is made that the process is stationary and thus that any given term of the covariance matrix can be determined from a correlation function, k, that only depends on the time difference, k(t_i, t_j) = k(τ_ij), with τ_ij = |t_j − t_i|. Adopting models developed by Foreman-Mackey et al. (2017), we chose a function of the form

$$k({\tau }_{ij})=\frac{A}{2}{e}^{-{\tau }_{ij}/{\tau }_{\text{d}}}[\cos (2\pi {\Omega }_{\text{a}}{\tau }_{ij})+1]+{\sigma }^2{\delta }_{ij}\cdot$$(B.1)

The cosine term on the right-hand side models the rotation rate. The second term including σ² is sometimes dubbed “jitter” and may capture potential model errors or compensate for underestimated observational errors (Angus et al. 2018).

The parameters of the covariance model are λ ≔ (A,τ_d, Ω_a, σ). We used a Bayesian statistical model to describe each time series. The posterior density of λ was estimated using an MCMC algorithm (Foreman-Mackey et al. 2013). An example of the resulting sampling is shown in Fig. 1. The likelihood is given by the density of the Gaussian process. The components of λ were assumed independent and the prior could be written as a product of uniform univariate densities. After several tests we found that we could efficiently sample the posterior density if we fit A and τ_d in logarithmic space and Ω_a and α in linear space. We chose uniform priors for Ω_a and α and log-uniform priors for A and τ_d. The adopted boundaries on A and τ_d are respectively [−∞, +∞] and [−10, 12]. For Ω_a and α they are [150, 1600], [150, 4500]. We note that the prior on Ω_a was set so that our results do not get contaminated by very low- or very high-frequency signals that are not properly reproduced by our model (hence relating to filtering practices that can be found in other studies; Angus et al. 2018). These boundaries were set after several tests and trials.

The prior on the lifetime also demands some caution. It should be noted that its upper boundary is far greater than the duration of any time series we are using in this study. And indeed, some estimates of τ_d from the MCMC simulation happen to be greater than a quarter duration. This can be viewed in several ways. First, it should be remembered that we are fitting a stochastic process to the time series. Obtaining long lifetimes out of the estimation process just means that high values of τ_d are compatible with the degree of correlation between the different timescales probed by the time series and that this degree always remains high. From another perspective, the stationarity assumption under which we are working breaks down if we try to only model one active region over long series. Consequently, when we encounter large estimates for the lifetimes of the active regions, the precise value of the estimate should be seen as meaningless, or at least severely biased. What is important in these cases is that τ_d is larger than the duration of the time series, indicating an active region that is stable over the duration of the entire subsample.

Appendix C: Conjunction of information states

The method of inversion is based on the framework described by Tarantola & Valette (1982). In this section we just summarize the main points of the approach and refer to the original paper for further details. The basic postulate is that any state of knowledge on the values of a set of parameters can be described using a density function³,⁴. The approach of the conjunction of states of information consists in expressing the posterior information on a couple (d, m), with d and m being the data and model-parameter vectors, using density functions

$$\sigma (d,m)=\frac{\rho (d,m)\Theta (d,m)}{\mu (d,m)},$$(C.1)

where the densities α, ρ, Θ, and μ represent states of information. They are not necessarily probability densities, which would require that they be normalizable. The density α represents the conjunction of the two states of information described by ρ and Θ. The density μ(d,m) is sometimes called a homogeneous probability density and represents the state of null information (Tarantola 2004). This means that μ is the density that is the least informative⁵ on the values of the couple (d, m). It is present in Eq. (C.1) so that the conjunction of any state of information by the null state results in no loss of information. The functions ρ and Θ are, respectively, the prior probability density and the probability density of the theoretical model on (d, m). The function ρ contains the information on the system (parameter and data) before the inference process. This can be either the observational data or the prior information on the parameters of the model. The density Θ represents the uncertainties on the theoretical model.

We now switch back to the notations of Sect. 2 using the relations d = Ω_i and m = θ_i. The data and the parameters are independent, hence we can write μ(Ω_i, θ_i) = μ_Ωi(Ω_i)μ_θi(θ_i). In this work we choose μ_Ωi and μ_θi as uniform densities (Tarantola 2004). Using the same argument of independence, we can write the prior probability density as a product of probability densities ρ(Ω,θ) = ρ_Ωi(Ω_i)ρ_θi(θ_i), where ρ_Ωi is obtained from the Gaussian-process modelling and ρ_θi is uniform on the latitude interval [0, 90°] since we have no prior information on the latitude of the active region. The range chosen is the consequence of our lack of spatial resolution between the two stellar hemispheres (see Sect. 3). Finally we write Θ(Ω_i,θ_i) = η_Ωi(Ω_i|θ_i)μ_θi(θ_i) with η_Ωi(.|θ_i) the probability density on Ω_i conditional on Θ = θ_i, i.e. the probability density of the theoretical rotation rate taken at Θ_i (see the central panel of Fig. 2 for an illustration). These choices ensure that α is normalizable and can be treated as a probability density.

The final step, once σ(Ω_i,θ_i) has been estimated, is to obtain the posterior probability density on Θ_i only instead on (Ω_i, θ_i). This is achieved by marginalizing the rotation rate, i.e. performing the integral

$${\sigma }_{{\theta }_i}({\theta }_i)={\displaystyle \int \sigma }({\Omega }_i,{\theta }_i)d{\Omega }_i,$$(C.2)

over the relevant range of rotation rates. The density σ_θi is called the marginal probability density for θ_i. The σ_θi are the objects represented in Fig. 3 as a function of time t and θ_i and is what we call the butterfly diagram.

References

All Figures

Fig. 1. Upper panel: time series for the relative flux HD 173701 during Q3. The black dots mark the measurements and associated errors. The red line shows the inferred mean value of the posterior predictive density conditional on the inferred MAP of the parameters (Gelman 2004). Middle panel: corresponding power spectrum computed using the Lomb-Scargle periodogram (Scargle 1982). The black line represents the observed power spectrum. The red line is the power spectrum of the mean value of the conditional posterior predictive density. Lower panel: marginalized 1D and 2D PDFs for the parameters of the correlation function of the Gaussian process used to model the time series.

In the text

	Fig. 2. Probability density for the latitude of an active region at median time of Q3. Left panel: theoretical density for the couple (Ωa, θa). Central panel: prior density for the couple (Ωa, θa); it is Gaussian for Ωa and uniform for θa. Upper right panel: posterior density for the parameters. Lower right panel: marginal density for θa, obtained after integration over Ωa.
In the text

Fig. 3. Asteroseismic butterfly diagram of HD 173701, giving the latitude of the active region as a function of the median time of each quarter (labelled on the upper axis). The grey shades show the posterior densities obtained for each of the selected quarters and normalized to the global maximum; greyscale from 0 (white) to 1 (black). Red dots mark the maximum a posteriori estimates of the latitude. The vertical bars give the limits of the corresponding 68.3% credible intervals. The green squares show the maximum of the density of the 1024 estimates obtained from artificial time series. Lower panel: variation in amplitude of the signal, SQ, normalized to its maximum.

In the text

Fig. A.1. Top panel: power spectrum of the short-cadence Kepler data in the region of the eigenfrequency multiplet centred around the frequency νn=24,l=2,m=0. The black and red lines show respectively the data and the best-fit model of the power spectrum. The vertical red ticks mark the frequencies of the eigenmodes of the multiplet, for −2 ≦ m ≦ 2 . Bottom panel: splitting diagram for the multiplet. The first two terms of the splitting sequence correspond to the rotational effects, while β(ν) is an additional term that describes aspherical contributions to the eigenfrequencies. The red horizontal ticks correspond to those seen in the bottom panel.

In the text