© The Authors 2023

1. Introduction

The observations from the space-borne instruments MOST (Walker et al. 2003), CoRoT (Baglin et al. 2006), Kepler (Borucki et al. 2010), and TESS (Ricker et al. 2014) have been important to advancing our knowledge of stellar interiors. This is particularly true for Kepler, which could continuously observe tens of thousands of stars for nearly four uninterrupted years, enabling asteroseismic measurements that almost rival disc-integrated helioseismic measurements from a decade ago. Those precise measurements have, for example, enabled us to better understand stellar rotation and its impact on stellar evolution. The asteroseismology of Sun-like stars is based on the study of the pressure modes that are excited by turbulent convection in the outer layers of these stars. The acoustic modes may reach deep into the core or may be localised close to the stellar surface, giving access to the internal structure and dynamics of the stars. Rotation has a critical impact on the stellar structure and evolution, as it induces a material mixing process (Maeder & Meynet 2008). Rotation is also an essential ingredient of the dynamo effect (Thompson et al. 2003) and can lead to a distortion of the star’s shape due to the centrifugal force (Chandrasekhar 1969). For fast rotators, the rotational flattening must be taken into account. These stars show a variety of pulsation modes that are not present in slower rotators (e.g. Lignières et al. 2006).

Despite the Sun slow rotation rate, the solar asphericity can be measured by helioseismology. While the Sun is seen as oblate by acoustic modes during the quiet phases of its activity cycle, the modes sense a more complicated shape during activity maxima. During solar maxima, the frequencies of acoustic modes increase slightly due to the presence of active regions. The perturbation occurs near the surface and is stronger for modes that sense the active latitudes (below ≃30° for the Sun). Physically, the magnetic perturbations consist of several components that are not easy to disentangle, including stratification and wave speed perturbations (Libbrecht & Woodard 1990; Antia et al. 2000; Dziembowski et al. 2000). The magnetic perturbations affect the acoustic modes near the surface over only a few hundred kilometres, that is, on the order of 10⁻⁵ R_⊙, but it is significant enough to be measured.

The origin of stellar activity is not well understood, as it depends on the complex interplay between rotation, convection, and the magnetic field (Brun & Browning 2017). Stellar magnetic cycles have been observed in most cool stars (Simon et al. 2002) over a large range of the electromagnetic spectrum, and evidence of activity cycles has been observed in X-ray (e.g. Catura et al. 1975), radio waves (White 1999; White et al. 2017), and chromospheric emission lines (Vaughan 1983; Oláh et al. 2009, 2016) as well as through luminosity variation due to surface magnetic activity in the visible light (Hartmann & Rosner 1979; Silva-Valio et al. 2010; Ceillier et al. 2017). On Sun-like stars, the level of activity is often observed to be cyclic, with activity periods ranging from a few years to decades. Although there are relationships between the stellar age, the rotation period, and the level of activity (van Saders et al. 2016), the underlying mechanisms at play are not understood.

Since the advent of space-borne photometry and the observation of Sun-like pulsators by CoRoT, it became evident that helioseismic methods used for the Sun may be applied to asteroseismic observations. This has resulted in robust estimates of the average rotational splittings, which, in combination with the surface rotation rates inferred from photometric variability, indicate that main sequence stars have nearly uniform internal rotation rates (Gizon et al. 2013; Benomar et al. 2015; Nielsen et al. 2017). All the recent seismic studies of radial differential rotation agree that angular momentum transport in the radiative zone is much more efficient than considered in theory, even in the case of stars more massive than Sun-like stars, such as γ-Doradus stars (e.g. Mosser et al. 2012; Gehan et al. 2018; Ouazzani et al. 2019). For the best Kepler observations, asteroseismology has shown evidence of latitudinal differential rotation for main sequence stars (Benomar et al. 2018a) as well as of the radial differential rotation in red-giant-branch (RGB; e.g. Deheuvels et al. 2012, 2014).

This paper aims to provide a framework to study stellar activity and its latitude by analysing its effect on pulsation frequencies. The proposed method involves the use of the a-coefficient decomposition (Schou et al. 1994; Pijpers 1997, 1998) on the stellar power spectrum to conveniently separate the perturbations caused by rotation and asphericity. The method is tested on Sun-as-a-star data and on the solar analogues 16 Cyg A and B, which are the two brightest stars in the initial Kepler observation field.

Only a few successful measurements of the asphericity of stars other than the Sun have been made so far. Using ultra-precise measurements of the frequency splittings of time-harmonic (i.e. non-stochastic) low-degree p modes, Gizon et al. (2016) inferred the oblateness of the γ Doradus–δ Scuti star KIC 11145123 to be ΔR = (1.8 ± 0.6)×10⁻⁶R ≃ 3 ± 1 km, which is smaller than expected from rotational oblateness alone, suggesting the presence of magnetic activity at low latitudes. Bazot et al. (2019) measured the asphericity of the solar-like pulsators 16 Cyg A and B and found that 16 Cyg A is likely prolate, implying that this star may have low-latitude magnetic activity on its surface.

In the spirit of the study by Gizon (2002), we perform Monte Carlo simulations to demonstrate the possibility of inferring the even a-coefficients from simulated oscillation power spectra in order to constrain the latitude of activity. Unlike Gizon (2002), who included only the a₂ coefficient in the parametric model, we infer both the a₂ and a₄ coefficients.

In Sect. 2, we start by introducing the effects of rotation on pulsation frequencies and discuss the effect of the centrifugal force and of the activity on the mode cavities. In Sect. 3, we present the assumptions required for the asteroseismic measurement of stellar activity. We discuss the achievable accuracy of the seismic observables in Sect. 4, while in Sects. 5 and 6 we present the results on solar data and for 16 Cyg A and B. This is followed by a discussion and conclusion in Sect. 7.

2. The effects of rotation and magnetic activity

This section presents the effects of rotation and magnetic activity on pulsation modes. It also introduces the frequency model used for asteroseismic data analysis.

2.1. Frequency splittings

Slowly rotating stars without significant magnetic activity are approximately spherical, and it is common to describe the family of modes travelling inside them using spherical harmonics (see e.g. Unno et al. 1989). If the departure from spherical shape remains small enough, it is convenient to keep the spherical representation for the equilibrium model and account for distortions through a perturbation analysis. All pulsations can then be described using a set of integers (n, l, m), namely, the radial order, the mode degree, and the azimuthal order. Acoustic pressure modes observed in Sun-like stars can be identified using their frequencies ν_nlm.

In a non-rotating, non-active star, m-components are degenerate and cannot be resolved. When rotation or magnetic activity sets in, this degeneracy is lifted. The resulting frequency is treated as a perturbation to the degenerate frequency without rotation and activity,

$$\begin{array}{c}\hfill {\nu }_{\mathit{nlm}}={\nu }_{\mathit{nl}}+\delta {\nu }_{\mathit{nlm}},\end{array}$$(1)

with ν_nl as the equilibrium eigenfrequency without rotation and activity, and δν_nlm as the frequency splitting. The splitting may depend on multiple physical effects perceived by the modes (Libbrecht & Woodard 1990) within their cavity of propagation. It can be a term of the order 𝒪(Ω), with Ω the rotation rate estimated at the equator, Ω = Ω(r, θ = π/2). The splitting depends directly on the stellar rotation profile Ω(r, θ). Higher-order perturbations pertaining to the shape of the mode cavity can also exist. We note that ν_nl differs from ν_{nl, m = 0}, as the m = 0 components may have their frequency shifted by perturbations, such as the magnetic activity (see e.g. Fig. 1).

Fig. 1. Example of Lorentzian mode profiles showing the a-coefficients and their relationship with frequency spacings for l = 1 (a) and l = 2 modes (b). The orange spacing is Tn22 − Tn21, νnl is the m-averaged frequency of each multiplet, and the stellar inclination is 40°.

Splitting can be described using the Clebsch-Gordan a-coefficient decomposition (Ritzwoller & Lavely 1991), which corresponds to a representation of the splittings on a basis of polynomials ${\mathcal{P}}_j^{(l)}(m)$ of degree j in m,

$$\begin{array}{c}\hfill \delta {\nu }_{\mathit{nlm}}={\displaystyle \sum _{j=1}^{j_{max}}}{a}_j(n,l){\mathcal{P}}_j^{(l)}(m),\end{array}$$(2)

with the polynomials such that

$$\begin{array}{c}\hfill {\displaystyle \sum _{m=-l}^l}{\mathcal{P}}_i^{(l)}(m){\mathcal{P}}_j^{(l)}(m)=0\mathrm{when}i\ne j.\end{array}$$(3)

Here, a_j(n, l) is the a-coefficient of order j, and j_max = 2l is the maximum order to which the decomposition must be carried for a given degree. The standard set of polynomials used in this expansion are those normalised as per described by Schou et al. (1994).

This decomposition is extensively used in helioseismology and was used on Sun-like stars by Benomar et al. (2018a). Asteroseismology has been so far unable to observe modes of degree higher than l = 3 so that in the following, the discussion is restricted to l ≤ 3 and j_max = 2l = 6 due to the selection rule of the ${\mathcal{P}}_j^{(l)}(m)$. This limitation is the consequence of full-disc integrated photometric observations. An example of splitting including odd and even a_j coefficients is given in Fig. 1 for l = 1, 2.

This theoretical model leads to a natural interpretation of the observed splittings. One may decompose the splittings into their symmetrical, S_nlm, and antisymmetrical, T_nlm, parts. These components can then be respectively described as sums over the odd and even a-coefficients,

$$\begin{array}{cc}& S_{\mathit{nlm}}=\frac{{\nu }_{n,l,m}-{\nu }_{n,l,-m}}{2m}=\frac{1}{m}{\displaystyle \sum _{j=1}^{j_{max}/2}}{a}_{2j-1}(n,l){\mathcal{P}}_{2j-1}^{(l)}(m),\hfill \end{array}$$(4)

$$\begin{array}{cc}& T_{\mathit{nlm}}=\frac{{\nu }_{n,l,m}+{\nu }_{n,l,-m}}{2}-{\nu }_{n,l,0}={\displaystyle \sum _{j=1}^{j_{max}/2}}{a}_{2j}(n,l)({\mathcal{P}}_{2j}^{(l)}(m)-{\mathcal{P}}_{2j}^{(l)}(0)).\hfill \end{array}$$(5)

This arises from the parity relation ${\mathcal{P}}_j^{(l)}(-m)={(-1)}^j{\mathcal{P}}_j^{(l)}(m)$. These equations can be used to express the a-coefficients with S_nlm and T_nlm. They also provide relations between the ν_nlm and the a_j(n, l) (see Appendix A.2).

It can be seen from Appendix A.1 that the sums in Eqs. (4) and (5) respectively involve odd and even functions of m. Physically, this means that the symmetrical components of the splittings (Gough & Thompson 1990) result from large-scale perturbations sensitive to the prograde or retrograde nature or the waves, such as advection or the Coriolis force. However, the antisymmetric splittings are caused by processes that are not affected by the propagation direction of waves. This may include the centrifugal force that scales as 𝒪(Ω²) and whose effect on the oscillation frequencies varies as m². Magnetic fields or non-spherical deformations of the equilibrium structure also contribute to the antisymmetric splittings.

We further decompose the antisymmetric splittings into a term depending on centrifugal force-induced distortions and another term accounting for activity-related distortions,

$$\begin{array}{c}\hfill T_{\mathit{nlm}}=\delta {\nu }_{\mathit{nlm}}^{(\mathrm{CF})}+\delta {\nu }_{\mathit{nlm}}^{(\mathrm{AR})}.\end{array}$$(6)

The symmetric part of the splitting corresponds to a term $\delta {\nu }_{\mathit{nlm}}^{(\mathrm{rot})}$ that stems from perturbations of the order 𝒪(Ω), and the total observed splitting is,

$$\begin{array}{c}\hfill \delta {\nu }_{\mathit{nlm}}=\delta {\nu }_{\mathit{nlm}}^{(\mathrm{rot})}+\delta {\nu }_{\mathit{nlm}}^{(\mathrm{CF})}+\delta {\nu }_{\mathit{nlm}}^{(\mathrm{AR})}.\end{array}$$(7)

Using the Sun as an archetype of a Sun-like star, it is possible to measure these contributions to δν_nlm, as shown in Sects. 2.2 and 2.3.

2.2. Expressions for $\delta {\nu }_{\mathit{nlm}}^{(\mathit{rot})}$, $\delta {\nu }_{\mathit{nlm}}^{(\mathit{CF})}$, and $\delta {\nu }_{\mathit{nlm}}^{(\mathit{AR})}$

To the first order, the perturbation on the frequency due to rotation is

$$\begin{array}{c}\hfill \delta {\nu }_{\mathit{nlm}}^{(\mathrm{rot})}=\frac{m}{2\pi }{\displaystyle {\int }_0^R{\int }_0^{\pi }{K}_{\mathit{nlm}}(r,\theta )\mathrm{\Omega }(r,\theta )r\mathrm{d}r\mathrm{d}\theta ,}\end{array}$$(8)

where R is the radius of the star and the kernel K_nlm(r, θ) expresses the sensitivity of a mode to the rotation at the radial point r and co-latitude θ (Hansen et al. 1977). The rotation profile of the star is represented as Ω(r, θ). It can be shown that $\delta {\nu }_{\mathit{nlm}}^{(\mathrm{rot})}$ actually only depends on symmetrical splittings (Ritzwoller & Lavely 1991), which in turn depend only on odd coefficients. For example, for l = 3, it is expressed as,

$$\begin{array}{c}\hfill \delta {\nu }_{\mathit{nlm}}^{(\mathrm{rot})}\simeq {\mathcal{P}}_1^{(l)}(m)a_1+{\mathcal{P}}_3^{(l)}(m)a_3+{\mathcal{P}}_5^{(l)}(m)a_5.\end{array}$$(9)

Centrifugal forces typically distort a spherical rotating sphere of gas into an oblate ellipsoid elongated at the equator. Functional analysis shows that the contribution of centrifugal forces to the mode splitting is proportional to Ω²R³/𝒢M. Properly integrating higher-order terms of the perturbation expansion over the aspherical volume of the star and using asymptotic expressions for the equilibrium mode eigenfunctions (assuming n is large enough), leads to the following expression for the centrifugal-force component of the frequency splitting (Gough & Taylor 1984; Gough & Thompson 1990),

$$\begin{array}{c}\hfill \delta {\nu }_{\mathit{nlm}}^{(\mathrm{CF})}=\frac{{\mathrm{\Omega }}^2{R}^3}{\mathcal{G}M}{\nu }_{\mathit{nl}}{Q}_{\mathit{lm}},\end{array}$$(10)

with the $Q_{\mathit{lm}}\approx \frac{2}{3}\frac{l(l+1)-3m^2}{(2l-1)(2l+3)}$ factor depending on the density.

Equation (10) means that

$$\begin{array}{c}\hfill \delta {\nu }_{\mathit{nlm}}^{(\mathrm{CF})}={\mathcal{P}}_2^{(l)}(m)a_2^{(\mathrm{CF})}(n,l)\end{array}$$(11)

with

$$\begin{array}{c}\hfill a_2^{(\mathrm{CF})}(n,l)=-\frac{1}{2l+3}\frac{{\mathrm{\Omega }}^2{R}^3}{\mathcal{G}M}{\nu }_{\mathit{nl}}\end{array}$$(12)

since $-(2l+3)Q_{\mathit{lm}}={\mathcal{P}}_2^{(l)}(m)$. We note that the contribution of the centrifugal force-induced deformation to the frequency splittings can be described by a linear combination of the 𝒫_2j polynomials (Gough & Taylor 1984; Gough & Thompson 1990). Since the stars are assumed to be slowly rotating, one may only retain terms of the order 𝒪(Ω²), which correspond to the contribution of the a₂ coefficient alone.

One can further approximate Eq. (12) in order to express $\delta {\nu }_{\mathit{nlm}}^{(\mathrm{CF})}$ as a function of quantities that can be obtained directly from the modelling of the acoustic power spectra of Sun-like stars (Benomar et al. 2018a). We consider first that the stellar mean density of a Sun-like star scales to a good approximation with its large separation. The large separation is the average distance in the frequency space between two modes of identical degree and consecutive orders (Ulrich 1986). In solar units, this reads as ρ = (ρ_⊙/Δν_⊙)Δν, with the solar density ρ_⊙ = (1.4060  ±  0.0005)×10³ kg m⁻³ and Δν_⊙ = 135.20  ±  0.25 μHz (García et al. 2011a). With an accuracy estimated to a few percentage points for Sun-like stars (White et al. 2011), the use of this scaling relation is thought to be a decent approximation.

The second simplification uses the fact that the Clebsch-Gordan coefficient decomposition of the frequency splitting imposes a one-to-one relationship between the a-coefficients and the coefficients associated to the decomposition of the velocity field into poloidal and toroidal components (Ritzwoller & Lavely 1991). Helioseismology suggests that the Sun rotates with a near constant angular velocity down to at least r/R_⊙ = 0.2 (e.g. Thompson et al. 2003), which is the maximum depth at which measurements from low-degree p modes are available. Furthermore, its outer convective zone shows a differential rotation of ≃30% from the equator to the pole, which leads to an a_3, ⊙(n, l)≃4 nHz and even smaller values for higher odd a-coefficients. This is significantly smaller than a_1, ⊙(n, l)≃420 nHz. Therefore, the aforementioned one-to-one relation ensures that we can retain only the leading order in the expansion of the rotation rate and treat it as an average value, given in terms of seismic observables by Ω ≃ 2πa₁. This leads to

$$\begin{array}{c}\hfill a_2^{(\mathrm{CF})}(n,l)\simeq -\frac{{\nu }_{\mathit{nl}}}{2l+3}\frac{3\pi }{\mathcal{G}{\rho }_{\odot }}\frac{\mathrm{\Delta }{\nu }_{\odot }^2}{\mathrm{\Delta }{\nu }^2}{a}_1^2.\end{array}$$(13)

Regarding $\delta {\nu }_{\mathit{nlm}}^{(\mathrm{AR})}$, there is no unambiguous theory to describe the effect of the near-surface magnetic activity on the shape of the cavity. Due to this and following Gizon (2002), a geometrical description is preferred to a physical model. This description assumes that the corresponding wave speed perturbation separates in the latitudinal and radial coordinates. The proposed form of the perturbation in frequency is

$$\begin{array}{cc}& \delta {\nu }_{\mathit{nlm}}^{(\mathrm{AR})}={\nu }_{\mathit{nl}}{ϵ}_{\mathit{nl}}{\displaystyle {\int }_0^{2\pi }{\int }_0^{\pi }}F(\theta |\mathit{x}){|Y_l^m(\theta ,\varphi )|}^{2}sin\theta \mathrm{d}\theta \mathrm{d}\varphi \hfill \\ \hfill & ={\nu }_{\mathit{nl}}{ϵ}_{\mathit{nl}}{A}_{\mathit{lm}}(\mathit{x})\hfill \\ \hfill & ={\displaystyle \sum _{j=0}^{j_{max}/2}{\mathcal{P}}_{2j}^{(l)}(m)a_{2j}^{(\mathrm{AR})}(n,l).}\hfill \end{array}$$(14)

The term $a_{2j}^{(\mathrm{AR})}(n,l)$ in Eq. (14) refers to the combined contribution of a magnetic field and other perturbations in the stellar structure (e.g. stratification, temperature). The geometrical weight function A_lm(x) describes the effect of an active zone at the co-latitude θ on a mode of degree l and azimuthal order m. It is also the product of two contributions. The first is from the normalised spherical harmonics $Y_l^m(\theta ,\varphi )$ that decompose the magnetic activity effect over each mode. These are spherical-polar coordinates defined in the inertial frame with a polar axis pointing in the direction of the rotation axis. The second is from the weight distribution (or the shape of the active region) defined by F(θ|x). Here, x refers to the parameters that are necessary to describe the function F(θ|x).

On the Sun, large active regions persist on the surface for one to two rotation periods and are randomly distributed in longitude over well-defined latitudes. We assume that the corresponding perturbation can be approximated by a function of latitude only. The general problem of distinct active regions on the differentially rotating surface goes beyond the scope of our study (see Papini & Gizon 2019, for the case of a single long-lived active region). In this paper, we only consider perturbations that are approximately steady in the inertial frame (i.e. latitudinal bands of activity).

In theory, a third integral over the radius is necessary to describe the dependence of the magnetic activity to the stellar depth. However, p modes are weakly sensitive to the deep structure inside stars. Therefore, the radial integral is replaced by ϵ_nl, the overall activity intensity. The frequency ν_nl allows for a dimensionless ϵ_nl that can be compared between stars. Section 3.2 further develops the required assumptions in order to obtain reliable information content on the active region in asteroseismology.

2.3. Modelling the frequencies

Combining Eqs. (1), (7), (8), (13), and (14) leads to

$$\begin{array}{cc}\hfill {\nu }_{\mathit{nlm}}& ={\nu }_{\mathit{nl}}(1+\frac{3\pi }{(2l+3)\mathcal{G}{\rho }_{\odot }}\frac{\mathrm{\Delta }{\nu }_{\odot }^2}{\mathrm{\Delta }{\nu }^2}{a}_1^2)\hfill \\ \hfill & +{ϵ}_{\mathit{nl}}{A}_{\mathit{lm}}(\mathit{x})+\frac{m}{2\pi }{\displaystyle {\int }_0^R{\int }_0^{\pi }{K}_{\mathit{nlm}}(r,\theta )\mathrm{\Omega }(r,\theta )r\mathrm{d}r\mathrm{d}\theta .}\hfill \end{array}$$(15)

Using Eqs. (1) and (2), ν_nlm can also be expressed using the a-coefficients without loss of generality,

$$\begin{array}{c}\hfill {\nu }_{\mathit{nlm}}={\nu }_{\mathit{nl}}+{\displaystyle \sum _{j=1}^{2l}{\mathcal{P}}_j^{(l)}(m)a_j(n,l).}\end{array}$$(16)

It is possible to use Eq. (15) to measure the activity parameters x by directly fitting the power spectrum. However, there are multiple benefits to using a two-step approach consisting of first using Eq. (16) for the fitting of the power spectrum to evaluate the a-coefficients (method detailed in Appendix B) and then fitting the coefficients obtained during the first step using solely Eqs. (13) and (14) (detailed in Appendix C). Firstly, and in the general case, the evaluation of A_lm(x) requires the precise computation of a double integral. This is a slow process that increases the time necessary to fit the power spectrum¹. Postponing to an ulterior step, the computation of A_lm(x) drastically reduces the parameter space from a few tens of parameters to just a few², effectively making the convergence rate of the fitting algorithm faster and making it easier to explore the assumptions that have to be made in Eq. (15) to have a functional approach in real cases (see the discussions in Sect. 3). Secondly, it allows us to decouple the observables from the physics, enabling the testing of various physical assumptions without having to reperform the lengthy power spectrum analysis. Finally, it eases the evaluation of the reliability of the activity determination by enabling us to pinpoint the cause of biases (if any) on either the observables (e.g. a₂, a₄; see Sect. 4) or on the interpretation of these observables in terms of physical parameters (see Sect. 3.3). The disadvantage of the two-step approach is that it requires more statistical assumptions, such as neglecting correlations and assuming Gaussian parameters. Although one can argue that it is possible to perform a hierarchical Bayesian analysis (e.g. Hogg et al. 2010; Campante et al. 2016) to partly alleviate this issue, such an approach is generally slow and may cancel the benefits of the two-step approach.

3. Information content in low-order a-coefficients

The general theoretical formulations of Sect. 2 demonstrate how pulsation frequencies are expressed as a function of the rotation, the centrifugal distortion, and the activity of stars. However, observational limitations need to be accounted for in order to enable a viable, robust model that extracts as much of the information content within currently existing data as possible. This necessarily requires additional assumptions, for which the rational is detailed hereafter.

3.1. Factors contributing to the splitting accuracy and precision

As explained in Kamiaka et al. (2018), the pulsation height (H) and the noise background (N) are important factors that reduce the capabilities of asteroseismic analyses. These factors are a complex function of the global stellar characteristics (mass, radius, age) of the convective motion at the surface of Sun-like stars and also depend on the instrumental limitations. The mode height and the noise background are in fact difficult to evaluate a priori. However, the height-to-noise (HNR), defined as the ratio H/N, can be used to assess the quality of the spectrum of a pulsation mode. In fact, Kamiaka et al. (2018) have shown that the HNR at maximum of mode height, $\widehat{\mathit{HNR}}$, can be used to study biases on the stellar inclination. This is because all main sequence Sun-like stars show a similar dependency on the height and on the width as a function of the frequency (see e.g. Appourchaux et al. 2014). In Sect. 4, we propose the same approach but on a-coefficients.

The $\widehat{\mathit{HNR}}$ not only defines how many modes can be observed but also the maximum degree of the modes that is observed. Considering specifically the Kepler observations, stars exhibit an $\widehat{\mathit{HNR}}$ for l = 0 modes of up to 30, the highest HNR being observed for 16 Cyg A and 16 Cyg B. Kamiaka et al. (2018) showed that the capability of distinguishing rotationally split components is of great importance if one wants to obtain a reliable asteroseismic inference of the rotation characteristics and of the stellar inclination. In particular, the ratio a₁/Γ_νmax between the a₁ coefficient and the mode width determined at the maximum of mode amplitude Γ_νmax (see Fig. 2) determines the expected bias on the stellar inclination. As shown in Sect. 4, a₁/Γ_νmax also controls the importance of the bias for the other low-order a-coefficients. Finally, another important factor is the spectrum resolution. The higher the resolution, the more resolved are the modes. Thus, observations T_obs of several years are the most suitable in order to resolve and measure rotationally split components. Broadly speaking, observations exceeding a year and a₁/Γ > 0.4 are preferable to ensure a reliable measurement of a₁. All the limitations discussed above incite us to introduce assumptions in order to ensure robust measurements of a-coefficients (i.e. mitigate biases).

Fig. 2. Height-to-noise ratio (black) and a1/Γ ratio (red) for 16 Cyg A. We used 16 Cyg A to make our simulations.

3.2. Latitudinal profile of the activity

Noting that it is challenging to measure low-degree a-coefficients for the Sun (e.g. Toutain & Kosovichev 2001), we present a minimal set of assumptions on even a-coefficients that allow us to constrain the asphericity of stellar cavities. One of the first aspects that has to be considered is the form of F(θ|x), the function characterising the activity latitude θ (see Eq. (14)).

The butterfly diagram of the Sun (Fig. 3) is used as a reference for this latitudinal dependence. The data are from the Greenwich USAF/NOAA observatory³ and provide the daily area of the spots, counted manually over the period 1874-2016. Panel a of Fig. 3 shows the butterfly diagram with colours representing the area covered by the spots in units of percent of visible hemisphere. It focuses on the observations after 1985 and covers two full solar cycles. Vertical colour bands highlight three time intervals : 1985–1989 (purple), 1999–2002 (blue), and 2006–2009 (yellow). The first two are during a maximum of solar activity, while the last one is during a minimum of activity. Due to the gradual migration of the spots over time, the longer observation period (1985–1989; i.e. 4 yr) leads to a broader activity zone than the period 1999–2002. This indicates that the extension of the active region may not be trivial to measure in other stars because it will depend on the fraction of time the star is observed relative to the duration of its activity cycle. The activity cycles of Sun-like stars (if any) are a priori unknown, but Ca II H+K line emission and photometric studies (Oláh et al. 2009, 2016) suggest that it is of durations of roughly a few years to decades, which is the same as for the Sun. It indicates that over the course of several years, an active band as large as ≃40° may be expected.

Fig. 3. Solar activity over time. (a) Butterfly diagram for the Sun. Vertical colour bands highlight periods of maximum (1988–1992, 1999–2002) and minimum activity (2006–2009). (b) Averaged spot area for the highlighted periods. Blue and purple marks in (b) are the weighted mean for the activity latitude. (c) Filters F(θ|X) used for this study superimposed onto solar data (1986–2016). The parameters θ0 and δ are the latitude and the extension of the active region, respectively.

Panel b of Fig. 3 shows the cumulative area of spots as a function of the latitude as well as for the three considered periods. As noted earlier, the extension of the activity band is larger for the longest time frame. The area of the spots during the active solar phase are symmetrical towards the equator. This suggests that when the activity is strong, F(θ|x) is almost north-south symmetric. This may be inaccurate for the low activity phase, as shown for the period 2006–2009, but because $|Y_l^m(\theta ,\varphi ){|}^2$ is also symmetrical towards the equator, this has no incidence on the $a_j^{(\mathrm{AR})}(n,l)$ coefficients.

During the minimum, the total average area of the spots is a few times lower than during the maximum. It is also narrower, as the integral in Eq. (14) is small, reducing $a_j^{(\mathrm{AR})}(n,l)\simeq 0$. During the phase of the minimum of activity, the activity can effectively be considered as non-existent (see Sect. 5 for the analysis on solar data) so that a₂(n, l) is dominated by the centrifugal term $a_2^{(\text{CF})}(n,l)$ and the other even a-coefficients are null.

The overall latitudinal profile seemingly follows a bell shape, with sharp slopes during periods of activity. Figure 3c shows a comparison of data between 1986 and 2016 with three models for the active latitudes: a model using the gate function F(θ|x) = Π(θ₀, δ), a triangular function Λ(θ₀, δ), and a Gaussian function 𝒩(θ₀, δ). The parameters of these were adjusted manually to approximately match the solar spot active latitudes profile. The triangular and the Gaussian functions equivalently describe the data, while the gate function initially proposed by Gizon (2002) roughly fits the data. Because spots strictly appear at latitudes below 45°, the Λ function may be the most suitable for the Sun. Nevertheless, we retained these three functions and compare them further in this work.

3.3. Assumptions on the a-coefficients

The most direct method for measuring asphericity is to evaluate it directly for each mode (n,l), that is, measuring the terms a_j(n, l). As this is already difficult for low-order a-coefficients of the Sun (Chaplin et al. 2003), it seems unreasonable to expect an accurate measurement of all individual a_j(n, l) using asteroseismic data. These stars have a lower signal-to-noise ratio than the Sun and severely reduced visibility at l ≥ 3 due to the integrated photometry. After a trial and error process, jointly with power spectra simulations, we could identify a set of assumptions ensuring reliable and precise measurement of a-coefficients.

We considered a fictitious star rotating as a solid-body with a solar activity level (ϵ_nl ≃ 5 × 10⁻⁴; Gizon 2002) with a₁ = 1000 nHz. The observed oscillation frequencies of 16 Cyg A are used (Davies et al. 2015; Kamiaka et al. 2018). We derived split frequencies using Eqs. (1), (7), (13), and (14) and converted them into a-coefficients using Eqs. (A.8)–(A.19). The a-coefficients are linear functions of the frequency, which is expected, as F(θ|θ₀, δ) describes a single active region. Thus, a reasonable assumption is to consider those a-coefficients as pure first order polynomial functions of frequencies. However, tests on artificial spectra showed that it is often difficult to evaluate the slope of the a-coefficients. This is because the uncertainty on any a_j is at least of the same order as its variations within the range of observed frequencies. This suggests that current asteroseismic data lack the resolution and the signal-to-noise to reliably measure the frequency dependence on the a-coefficients.

Figure 4 shows $a_2^{(\text{AR})}$, $a_4^{(\text{AR})}$, $a_6^{(\text{AR})}$ for activity described by F = Π (black) or Λ (red) or 𝒩 (blue). These are the mean coefficients for the fictitious active star, as a function of θ₀ and δ. The figure indicates that there is a simple relationship between the a-coefficients, the co-latitude θ₀, and the extension of the activity zone δ, independently of the shape of activity. In the case of F = Π, the lines are cut near the pole and the equator due to the condition θ₀ ≥ δ/2 and θ₀ ≤ π − δ/2. In the other profiles, edge effects (truncation) have a noticeable impact near the equator. The definition of δ differs between the three profiles, which explains why a factor of a few in $a_j^{(\text{AR})}$ is noticeable between F = Π, Λ, and 𝒩 for a given δ. We note that adding the centrifugal effect reduces $a_2(n,l)=a_2^{(\text{AR})}(n,l)+a_2^{(\text{CF})}(n,l)$ because the coefficient $a_2^{(\text{CF})}(n,l)$ is always negative.

Fig. 4. Average $a_j^{(\text{AR})}$ (no centrifugal effect accounted for) with ϵnl = 5.10−4 in function of θ and δ and for the gate (Π), triangle (Λ), and Gaussian (𝒩) filter functions. Uniqueness is guaranteed for θ only if at least two a-coefficients are measured.

If an observation can constrain only a single a-coefficient (e.g. a₂), there often exists a degeneracy in θ₀, as multiple values of a_j can be obtained for a given θ₀. As the figure shows, measuring two a-coefficients alleviates this degeneracy, provided that uncertainties are small enough. In other words, the accuracy on the inference of the active region using a-coefficients is ensured only if we can simultaneously constrain two a-coefficients (e.g. a₂ and a₄). This is essential to distinguishing an activity near the pole (θ₀ ⪅ 30°) from one near mid-latitude (30 ° ⪅θ₀ ⪅ 60°) or near the equator (θ₀ ⪆ 60°).

Figure 4 indicates that the loss of information resulting from the averaging does not affect the accuracy⁴. Figure 4 also indicates that for a star with an activity intensity and an activity zone of extension commensurate with that of the Sun (ϵ_nl ≃ 5 × 10⁻⁴, δ ≃ 10° when F = Π), the uncertainty for $a_2^{(\text{AR})}$ must be approximately ⪅25 nHz, for $a_4^{(\text{AR})}\lesssim 10$ nHz, and for $a_6^{(\text{AR})}\lesssim 1.5$ nHz to be able to detect significant departures of the coefficients. This is to be compared with the spectral resolution of 7–14 nHz for 2–4 yr of observation, which is typical of the longest Kepler observations. In Fig. 4, $a_2^{(\text{AR})}$ is obtained by averaging l = 1, 2, and $a_4^{(\text{AR})}$ is obtained by averaging l = 2, while $a_6^{(\text{AR})}$ is only constrained by l = 3 modes. Averaging $a_2^{(\text{AR})}$ and $a_4^{(\text{AR})}$ over l = 1, 2, 3 increases the maximum range of the a-coefficients by a factor of approximately 1.5 for $a_2^{(\text{AR})}$ and approximately a factor of two for $a_4^{(\text{AR})}$ without changing the overall shape of the function. In Sun-like stars, the l = 3 modes have an HNR at least ten times lower than l = 0 modes because the height ratio between l = 3 and l = 0 modes is around 0.08 for the Sun (Toutain & Gouttebroze 1993; Toutain et al. 1998). These modes are therefore difficult to observe. Due to all the above, it is extremely challenging to measure a₆ with existing data. In the following, we thus focus on assessing the reliability domains of a₂ and a₄ only in the following sections.

4. Bias analysis on a₁, a₂, and a₄

In order to understand the reliability of the inference on low-degree a-coefficients, it is necessary to perform a bias analysis. This requires fitting an ensemble of emulated spectra that are representative of Sun-like stars and comparing the results with the true inputs.

4.1. Test setup

The synthetic spectra used to analyse the biases use the frequencies, heights, widths, and the noise background profile of 16 Cyg A as a template for the simulations. We wanted to specifically study the impact of the mode blending (effect of a₁/Γ_νmax), stellar inclination (i), observation duration (T_obs), and maximum height-to-noise ($\widehat{\mathit{HNR}}$) background on the accuracy of a₁, a₂, and a₄. We built grids of spectra in the case of an equatorial band of activity and of a large polar activity. The grid parameters and their ranges are provided in Table 1 and discussed further in this section. The spectra were made using a spectrum simulator code⁵ that takes the reference star (or template) and modifies its properties to match those required by the user. The template was altered in terms of HNR following a similar approach in Kamiaka et al. (2018). The main difference is that we considered a frequency-dependent noise background. We rescaled heights according to

$$\begin{array}{c}\hfill H_{n,l=0}=\frac{\widehat{\mathit{HNR}}}{{\widehat{\mathit{HNR}}}_{\mathrm{ref}}}{H}_{\mathrm{ref}}(n,l=0),\end{array}$$(17)

where $\widehat{\mathit{HNR}}$ is the maximum HNR of the synthetic star, ${\widehat{\mathit{HNR}}}_{\mathrm{ref}}$ is the maximum HNR of the reference star, and H_ref(n, l = 0) is the l = 0 heights of the reference star.

As for the mode blending factor f = a₁/Γ_νmax, we calculated it fixing a₁ and altering Γ_νmax such that

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{{\nu }_{max}}=\frac{a_1}{f}.\end{array}$$(18)

We used a fiducial value of a₁ = 1000 nHz (≃2.4 the solar rotation) in the simulation. We note that this differs from Kamiaka et al. (2018), where the splitting was modified in order to obtain the desired mode blending factor.

As shown in Table 1, three HNR cases were investigated, ranging from 10 to 30. Two observation durations were considered: 2 yr and 4 yr. This is representative of the best Kepler observations (Davies et al. 2015) and of future observations from PLATO (Rauer et al. 2014). Similar to the test cases of Gizon (2002), a₂ and a₄ were determined for an equatorial activity of extension of the same order as in the Solar case (δ ≃ 10°) and in the case of a large polar cap (δ = 45 ° ). Both situations assume an activity of the same order as the Sun. The priors were set in a manner similar to how they would be if the power spectrum was from a real star. The evaluation of the mode parameters was performed using an MCMC method, just as for actual stars (see Appendix B for further details), but on the limit spectrum (no noise). Fitting the limit spectrum allowed us to assess the systematic errors by calculating the expectation value of the probability density function. The expected uncertainty can also be known by computing the standard deviation.

4.2. Results

To appreciate the level of inaccuracy achieved when measuring a_j, j ∈ [1, 2, 4] coefficients, we calculated the bias for a_j, f = a₁/Γ_νmax, and i:

$$\begin{array}{cc}\hfill b(a_j)& =\mathrm{E}[a_j^{\mathrm{meas}}]-a_j^{\mathrm{true}},\hfill \\ \hfill b(F)& =\mathrm{E}[F^{\mathrm{meas}}]-F^{\mathrm{true}},\hfill \\ \hfill b(i)& =\mathrm{E}[i^{\mathrm{meas}}]-i^{\mathrm{true}}.\hfill \end{array}$$(19)

Here, the letter E refers to the expectation value (median) from the probability density function.

Figure 5 shows the resulting bias maps for a-coefficients as a function of the blending factor F and of the stellar inclination for a-coefficients corresponding to an equatorial activity band. Figure 6 is the same but for a large polar activity cap. The bias is represented as a projected three dimensional vector using the three quantities defined by Eq. (19) and for j ∈ [1, 2, 4]. To evaluate the importance of the bias on a_j, the plot shows b(a_j)/⟨σ_aj(i, F)⟩_i, with ⟨σ_aj(i, F)⟩_i as the average standard deviation of the probability density functions (PDFs) obtained by the MCMC sampling over i ∈ [30 ° ,90 ° ]. The value of ⟨σ_aj(i, F)⟩_i (noted as σ for convenience hereafter) and its variance along the inclination axis is also shown on the plots. A value of b(a_j)/σ greater than one indicates that the inaccuracy exceeds the typical expected uncertainty at 1σ and may lead to significantly biased results during the subsequent analyses of the a-coefficients. A negative bias (underestimation) for a_j is indicated by a cross within the circle, while a positive bias (overestimation) is indicated by a dot. The darker the colour, the smaller the norm of the normalised bias. The size of the circle symbols is also proportional to the bias but capped to three in order to avoid excessively large symbols when i < 30° (see next paragraph for explanations). The cases with an observation duration of 2 yr are in the top panels (Figs. 5a–c) of the figure, while cases with a duration of 4 yr are shown in the bottom panels (Figs. 5d–f).

Fig. 5. Bias analysis for $\widehat{\mathit{HNR}}=30$ for an equatorial activity band (θ0 = 85°, δ = 10°) with similar intensity to the Sun (ϵnl = 5 × 10−4) for Tobs = 2 years (top) and Tobs = 4 years (bottom). The colour and size of the circles indicate the modulus of the bias. The colour bar gives its scale normalised by the uncertainty, b(aj)/σ. A white cross indicates an underestimation. A white dot is for an overestimation. Below 30° of inclination (grey area and symbols), the results are not reliable.

Fig. 6. Bias analysis for $\widehat{\mathit{HNR}}=30$ for a large polar activity cap (θ0 = 22.5, δ = 45) of similar intensity to the Sun (ϵnl = 5 × 10−4) for Tobs = 2 yr (top) and Tobs = 4 yr (bottom).

We note first that the biases on inclination and on F are consistent with results from Kamiaka et al. (2018). We found regions of stellar inclination below 30° provide very unreliable results, and the median of the inclination is biased toward ≃80° when the true input is i = 90°. Meanwhile, we found biases in a₁ often exceed 50% in the grey area⁶, indicating that for an inclination below ≃30°, the median estimator of the probability density function might not be reliable. This also justifies the range for the calculation of σ. In the following paragraphs, we therefore focus the discussion of the figures for the region above 30°.

In the case of an equatorial zone of activity, we note that measurements of a₁ may have negligible inaccuracies when the blending factor F is above 0.4, but it starts to be significantly deteriorated at (and probably below) f = 0.4 because |b(a₁)|/σ exceeds 1.5 when i < 50°. Split components then overlap significantly, leading to important degeneracies between a-coefficients, mode widths, and the stellar inclination. Meanwhile, we note that the biases on a₂ and a₄ remain mild compared to the uncertainty (|b(a₂)|/σ < 1 and |b(a₄)|/σ < 1), even for f = 0.4.

In the case of a polar active region, |a₂| represents ≃20% of the a₁ coefficient. In these conditions, we note that b(a₁)/σ does not exceed the unity (see Figs. 5a,d), indicating a good accuracy. The terms |b(a₂)|/σ and |b(a₄)|/σ are also lower than 0.5, provided that the stellar inclination is between 30° and 70°. This indicates that a large a₂ coefficient is generally associated with a higher accuracy for all a-coefficients. The plot also demonstrates that a large zone in the parameter space has a moderate to small bias.

Interestingly, figures for $\widehat{\mathit{HNR}}$ of 10 and 20 (Figs. D.1–D.4) lead to similar conclusions. Therefore, even in less favourable $\widehat{\mathit{HNR}}$ conditions, the expected uncertainty on the measurement generally encompass the bias, provided that the stellar inclination exceeds ≃30° (ensuring that m ≠ 0 have significant amplitudes) and that the mode blending factor is above 0.4. This indicates that for a majority of stars observed by Kepler, the accuracy may not be a major issue when measuring a₁, a₂, and a₄ if a careful assessment of F and i is performed. We note, however, that the average uncertainty σ does usually increase when the $\widehat{\mathit{HNR}}$ and/or the observation duration decrease, which has an impact on the precision of the determination of the active region.

5. Analysis of solar data

The Sun has a well-known 11-year activity cycle that makes it an ideal subject to further test the accuracy of the method described in Sect. 4. The Sun’s activity cycle is analysed at two instants, highlighted in Fig. 3 and for which high-quality helioseismic data are available: (a) the maximum of activity between January 1999 and January 2002 and (b) the minimum of activity between January 2006 and January 2009. Data are from the Variability of Solar Irradiance and Gravity Oscillations instrument on board the SOHO spacecraft (Frohlich et al. 1997) and present very few gaps.

5.1. a-coefficients of the Sun

As for the simulated data, the analysis consists of fitting the power spectrum for cases (a) and (b) in order to measure the a₁, a₂, and a₄ coefficients using the Bayesian modelling and the MCMC method described in Appendix B. Figures of the best fits and their discussions are provided in Appendix E. As the goal is to evaluate the accuracy of the inference of the activity zone based on the a-coefficients, the presented results are for a stellar inclination fixed to 90°, instead of having it as a free parameter. This value corresponds approximately to the stellar inclination seen by the SOHO satellite. This alleviates biases on a_j coefficients that may arise due to the systematic underestimation of i when it is close to 90 degrees. However, tests we carried out with a free inclination during the maximum of activity⁷ led to a difference of −12 nHz in a₂ and to −1.3 nHz in a₄. This outcome is consistent with the expected difference from the bias map of Fig. 5. In agreement with the discussion of Sect. 4, the difference does not have a significant impact on the activity inference because the uncertainties for those parameters are larger than the observed measurement shift.

Figure 7 shows the measured probability distribution function of the relevant parameters along with their correlations during the maximum (left) and the minimum (right) of solar activity. The $a_2^{(\text{CF})}$ distribution in red represents the expected centrifugal effect on a₂ as derived from frequency shifts ${\delta }_{\mathit{nlm}}^{(\mathrm{CF})}$ of Eq. (13). The distributions are Gaussians and show a weak correlation, reflecting the quality of the data. Between the maximum and the minimum of activity, the a₂ coefficient drifted significantly, from a₂ = 80  ±  19 nHz to a₂ = 11 ± 21, falling within the 1σ confidence interval of the centrifugal term. Although a₄ may have changed, the effect is below uncertainty levels and remains close to zero. As shown in Fig. 8, other time intervals may lead to a₄ departing from zero. The figures demonstrate that a comparison of a₂ and $a_2^{(\text{CF})}$ may reveal the activity of a Sun-like star. At the maximum of activity, $a_2^{(\text{CF})}$ is inconsistent with a₂ at 4.5σ, but it is in agreement at 1σ during the minimum of activity. Interestingly, Chaplin et al. (2003) also studied the frequency asymmetry of l = 2 modes in detail using BiSON and GOLF data with a different methodology and over a period that encompasses the maximum of 1999–2002. Their measurement considers only T_n22, so a direct comparison is not straightforward. However, we note that with an 844 day-long time series starting in February 1999⁸, they detected a frequency shift T_n22 ≃ 190 nHz at a similar significance (≃4σ) as our measurement when averaging over all modes between ≃2000 − 3300 μHz (to be compared to our range of ≃2300 − 3600 μHz). Their Fig. 8 also shows that the global effect of the activity cycle between 1994 and 2000 is evident in the averaged T_n22, while the frequency dependence of T_n22 may have uncertainties too large to ascertain a frequency trend. This is in line with our own findings (see our Sect. 3.3).

Fig. 7. Probability density functions and their correlations obtained by MCMC for coefficients a1, a2, a3, and a4 for the Sun at maximum activity (1999–2002, left) and at minimum activity (2006–2009, right). The red curve is the expected $a_2^{\text{CF}}$ coefficient for a pure centrifugal distortion. The light and dark grey PDF areas show the 1σ and 2σ confidence intervals, respectively.

Fig. 8. Probability density functions and their correlations obtained by MCMC for coefficients a1, a2, a3, and a4 for the Sun between two activity cycles (2006–2011). The red curve is the expected $a_2^{\text{CF}}$ coefficient for a pure centrifugal distortion. The light and dark grey PDF areas are the 1σ and 2σ confidence intervals, respectively.

A rigorous statistical evaluation of our significance requires the joint use of a₂ and a₄. This is discussed in Sect. 5.2, along with other activity results. Finally, we note that for an inclination of ≃90 deg, it is not possible to determine a₃ because the amplitudes of the azimuthal components for l = 1, 2 are not favourable (Gizon & Solanki 2004).

5.2. Activity of the Sun

We derived the activity intensity of the Sun and its latitudinal coverage from a fit of the $a_2^{(\text{AR})}$ and $a_4^{(\text{AR})}$ coefficients. Technical details are in Appendix C. A statistical summary of the results is presented in Table 2. The statistical significance of the activity exceeds 99.99% (highly significant) at the maximum of activity but is around 50.6% (not significant) during the minimum of activity, demonstrating the possibility of detecting activity in Sun-like stars. The choice of the function F describing the active region changes the detection significance by only a few percentage points. The three models we explored fit the data equivalently, implying that the shape of the active region cannot be determined with currently available a₂ and a₄ constraints.

Figures 9 and 10 show ϵ_nl = ϵ, θ₀, and δ for the Sun at its maximum (1999–2002) and its minimum (2006–2009) of activity, respectively. There are no major differences between the three activity profiles, although we do note that the uncertainties are larger in the case of a triangular and a Gaussian activity zone. During the maximum of activity, ϵ ≃ 0 is clearly excluded (see inset in Fig. 9). With F = Π, we observed a log-normal distribution, with median 7.6 × 10⁻⁴, consistent with the value reported by Gizon (2002). However, the uncertainty is large, suggesting that only the order of magnitude can be constrained in Sun-like stars. At the minimum of activity, the log-normal distribution morphs into a 1/x law, similar to the Jeffreys prior. This is a sign of weaker statistical significance for the activity.

Fig. 9. Inferred PDFs for (a) ϵ, (b) θ0, and (c) δ during the maximum activity of the Sun (1999–2002). The inset of (a) is a zoom into the near zero ϵ values with smaller binning.

Fig. 10. Inferred PDFs for (a) ϵ, (b) θ0, and (c) δ during the minimum activity of the Sun (2006–2009). The inset of (a) is a zoom into the near zero ϵ values with smaller binning.

The co-latitude θ₀ ≃ 75 < SPSDOUBLEDOLLAR > ° is consistent with the butterfly diagram, which suggests θ₀ ≃ 80°. At the minimum of activity, and despite a non-significant detection, the probability distribution of θ₀ shows a weak indication of activity at mid and high co-latitudes. Despite the high significance of the detection during the maximum of activity, δ is poorly constrained, showing that this parameter is challenging to measure on the Sun and on other Sun-like stars.

6. Analysis of 16 Cyg A and B

Due to their brightness (magnitudes V = 5.95 and 6.20), 16 Cyg A and B have modes with the highest HNR among all of Sun-like stars observed asteroseismically so far. They constitute ideal candidates to evaluate the activity. The two stars are wide binaries; there is a confirmed planet around 16 Cyg B (Cochran et al. 1997); and they were extensively studied (e.g. Neckel 1986; King et al. 1997; Deliyannis et al. 2000; Schuler et al. 2011; Takeda 2005; Metcalfe et al. 2012, 2016; Lund et al. 2014; Verma et al. 2014; Buldgen et al. 2015, 2022; Deal et al. 2015; Roxburgh 2017; Bellinger et al. 2017; Maia et al. 2019; Bazot et al. 2019; Bazot 2020; Farnir et al. 2020; Morel et al. 2021; Nsamba et al. 2022). In this section, we use the 2.5 year observation (13 September 2010 to 8 March 2013) by the Kepler space-borne instrument to measure pulsation parameters. The data are the same as those used in Bazot et al. (2019), for which instrumental issues (outliers, jumps, trends) and the quarter stitching is performed using the procedure described in García et al. (2011b). The binary system has precisely measured angular diameters, making 16 Cyg A and B two of the few Sun-like stars with known interferometric radii. The measured radii are 1.22 ± 0.02 R_⊙ and 1.12 ± 0.02 R_⊙ for 16 Cyg A and B, respectively (White et al. 2013). The spectroscopic parameters of these stars are very close to those of the Sun. The effective temperatures for 16 Cyg A and B are, respectively, T_eff = 5825 ± 50 K and T_eff = 5750 ± 50 K, and their metallicities are [M/H] = 0.10 ± 0.09 and [M/H] = 0.05 ± 0.06 (Ramírez et al. 2009), respectively. Each with an estimated age of around 7 Gyr (e.g. Metcalfe et al. 2016; Bazot 2020), the two stars are significantly older than the Sun.

6.1. Seismic constraints for 16 Cyg A

Earlier studies of 16 Cyg A have revealed around 60 modes of pulsations, with a significance in the power spectrum⁹. The asteroseismic analysis of Davies et al. (2015) found a stellar inclination of ${56}_{-5}^{+6}$ and a rotation of $23.8_{-1.8}^{+1.5}$ (${⟨\delta {\nu }_{\mathit{nlm}}/m⟩}_{\mathit{nl}}\simeq {486}_{-29}^{+40}$ nHz). Using a refined power spectrum modelling that accounts for a₁ and a₃ and that parameterises the cavity asphericity (R_eq − R_pol)/R_eq, Bazot et al. (2019) reported values (i = 58.5 ± 6.8, a₁ = 464 ± 43 nHz) consistent with Davies et al. (2015).

Figure 11 shows the PDFs for a₁, a₂, a₃, and a₄ for 16 Cyg A, along with their correlation and with $a_2^{(\text{CF})}$ superimposed in the a₂ quadrant. Table 3 synthesises the inferred values of the coefficients. The PDFs are near Gaussian, and thus the quadratic mean of the asymmetrical uncertainties is reported in the table. The a₁ coefficient is significantly lower than past estimates but remains consistent at a 2σ confidence level. The difference may be due to the lower stellar inclination, i = 45 ± 4, which is also only consistent with an earlier determination at 2σ. The a₃ coefficient is marginally greater but with a smaller uncertainty than Bazot et al. (2019) (they reported a₃ = 11.15 ± 10.95 nHz).

Fig. 11. Probability density functions and their correlations obtained by MCMC of the power spectrum of 16 Cyg A and for coefficients a1, a2, a3, and a4 for 16 Cyg A. The red curve shows the expected $a_2^{\text{CF}}$ coefficient for a pure centrifugal distortion of the star. The light and dark grey PDF areas show the 1σ and 2σ confidence intervals, respectively.

More importantly, a₂ is positive at 1σ, which is consistent with the assertion of star prolateness from Bazot et al. (2019) and is inconsistent with the centrifugal distortion term, $a_2^{(\text{CF})}$. Finally, a₄ includes zero within 1σ.

With respect to Figs. 5 and 6, a₁/Γ_νmax = 0.45 ± 0.04 and i = 45 ± 4° place 16 Cyg A in a parameter space where the expected magnitude of the bias for a₁, a₂, and a₄ is of the order of 0.5σ, 0.5σ, and 0.8σ, respectively. This remains accurate even when accounting for the slightly higher estimates of inclination from previous publications. Translated into absolute units, this corresponds to approximately b(a₁)≃ + 10 nHz (overestimation), b(a₂)≃ − 10 nHz (underestimation), and b(a₄)≃ − 10 nHz (underestimation). These fiducial values are used in Sect. 6.3 to evaluate the effect of the potential bias on the estimates of the activity zone.

6.2. Seismic constraints for 16 Cyg B

Similar to 16 Cyg A, 16 Cyg B has around 60 observed modes. However, the reported precision for the past determination of the seismic parameters is less accurate than for 16 Cyg A, despite a similar HNR. Davies et al. (2015) and Bazot et al. (2019) both note a large degeneracy between the stellar inclination and the average rotation rate, even with a clear bimodality in the distributions obtained by Bazot et al. (2019). Their global solution of $i={36}_{-7}^{+17}$ is associated with two separate solutions for a₁, centred around ≃300 nHz and ≃550 nHz, that they use to infer the latitudinal differential rotation profile. Meanwhile, as for 16 Cyg A, the star is found to be prolate, indicating a surface activity.

Figure 12 shows the measured a-coefficients. It appears that including the a₄ coefficient removed the degeneracy issue observed by previous publications. The stellar inclination, a₁, and a₂ are precisely determined and have approximately Gaussian probability distributions, and these are reported as such in Table 3. The found rotation rate a₁ corresponds to the higher solution of rotational splitting of Davies et al. (2015) and Bazot et al. (2019; see their figures of PDFs). Meanwhile, the degeneracy observed in a₁ in earlier studies is moved to a₄, and it exhibits two solutions. As our method of activity inference assumes Gaussian distributions, the bimodality of a₄ requires separating the two apparent solutions. The mean and standard deviation of both solutions were measured using Gaussian mixture modelling (Bishop 1995; Bazot et al. 2019), enabling their separate analysis. The lowest estimate’s (a₄ = −27.9 ± 6.5 nHz) weight is 30%, and the highest estimate (a₄ = −1.0 ± 6.7 nHz) is more significant, as its weight is 70%. The inference of the activity from the two solutions is discussed in Sect. 6.4. We note that $a_3={45}_{-14}^{+13}$ nHz is significantly higher than the reported values from Bazot et al. (2019) (a₃ = 13.89 ± 13.95 nHz) and may have an impact on the rotational profile.

Fig. 12. Probability density functions and their correlations obtained by MCMC of the power spectrum of 16 Cyg B and for coefficients a1, a2, a3, and a4 for 16 Cyg B. The red curve shows the expected $a_2^{\text{CF}}$ coefficient for a pure centrifugal distortion of the star. The light and dark grey PDF fillings are the 1σ and 2σ confidence intervals, respectively.

6.3. Activity inference for 16 Cyg A

Figure 13 shows the results from the inference of the activity parameters and a statistical summary is given Table 2. Triangular and Gaussian descriptions of active latitudes give larger uncertainties than the simple gate model. However, they all suggest near-equatorial activity, with a similar detection significance level of ≃79% (without bias correction), or ≃96% with fiducial correction. As for the Sun and because the shape of the active region is a priori unknown, the average distribution of the parameters is discussed. The activity intensity has a large uncertainty, but according to the median of ϵ, it may be between the maximum and the minimum of solar activity.

Fig. 13. Inferred PDFs for 16 Cyg A for raw measures (a–c) and with bias correction (d–f). The green curve is for the average PDF with F = Π,  Λ, and 𝒩. The shaded area is its 1σ confidence interval. The insets in (a) and (d) are zooms into the near zero ϵ values with smaller binning.

The posterior probability distribution of δ does not allow us to precisely constrain the extension of the activity region. As already noted in the case of the Sun, this parameter requires stringent constraints on both the a₂ and a₄ coefficients in order gain information about the size of the active region.

6.4. Activity inference for 16 Cyg B

Figures 14 and 15 show the results for the activity for 16 Cyg B in the two possible scenarios of a₄ discussed in Sect. 6.2 (see Table 2 for the statistical summary). In the case of a₄ = −27.9 ± 6.5 nHz (Fig. 14), the activity is highly significant (greater than 99.3%, with or without bias correction, whatever the activity zone model) and has a stronger intensity than in the case of the Sun. The active region is then located at θ₀ ≃ 58° (θ₀ ≃ 60°, after bias correction), that is, at comparable latitudes seen during a maximum of activity of the Sun. The bias correction has a negligible effect on the inferred activity latitude. Due to a₄ significantly departing from zero, the extension of the active region is better constrained than in 16 Cyg A or the Sun, but it remains weakly informative.

Fig. 14. Inferred PDFs of ϵ, θ0, and δ for raw results of 16 Cyg B with the lower solution a4 = −29.2 ± 5.6 nHz without bias correction (a–c) and with it (d–f). The green curve is the average PDF with F = Π,  Λ, and 𝒩. The shaded area is its 1σ confidence interval. The insets of in (a) and (d) are zooms into the near zero ϵ values with smaller binning.

Fig. 15. Same as Fig. 14 but for the higher solution a4 = −0.9 ± 6.8 nHz.

The second more likely solution (probability of 70%) corresponds to a₄ = −1.0 ± 6.7 nHz (Fig. 15). It is associated with an activity at co-latitudes above ≃40° and with lower activity. Its statistical significance is low, and the activity intensity may be of the same order or lower than 16 Cyg A. In fact, it looks similar to the solar case when approaching its minimum of activity. Accounting for the fiducial bias, the solution is more concentrated in the equatorial region and differs significantly from the lower probability a₄ solution. Here, the uncertainty on a₂ and a₄ is again too large to provide a stringent constraint on the extension of the activity zone δ.

7. Discussion and conclusion

Stellar activity is a complex phenomenon emerging from the interplay between stellar plasma, rotation, and magnetism. The activity distorts the shape of the mode cavity, perturbing the pulsation frequencies. In this work, we evaluated this perturbation observationally using an a-coefficient decomposition, as it allowed us to separate observables (the a-coefficients) from their physical interpretation (the activity modelling). The separation also enabled us to demonstrate that the measurements of the average low a-coefficient, a₂, and a₄, under some observational conditions, are sufficient to reveal the presence of a statistically significant activity of similar intensity to the Sun and to determine its latitude. The required observational conditions were analysed using a methodology similar to Kamiaka et al. (2018), that is, by constructing a grid of artificial power spectra that allows the bias for a₁, a₂, and a₄ to be determined. We found that if the HNR exceeds ten, then the mode blending factor f = a₁/Γ_νmax is greater than 0.4, the inclination is above 30°, and the observation is longer than 2 years. The inaccuracy remains mild and generally smaller than the 1σ uncertainty. The uncertainty and/or the inaccuracy may, however, become too large to reliably detect any activity beyond the above-specified conditions. In particular, and in agreement with Kamiaka et al. (2018), the stellar inclination is a decisive variable for ensuring the accuracy of the measurement. Below 30°, a₁ is often wrong by a factor two. This is likely due to the fit mistakenly identifying l = 2, m = ±1 as l = 2, m = ±2 and vice-versa. Other a-coefficients are then severely inaccurate. Therefore, such a bias analysis suggests that rotation studies of stellar ensembles require a careful star selection.

We proposed a method that uses the average a₂ and a₄ coefficients to perform a subsequent analysis of the activity, considering a geometrical model of the activity effect on the pulsation frequencies and accounting for the stellar asphericity due to the centrifugal effects. We tested this two-step approach on the Sun and showed that it is effectively able to detect the change of activity between the solar maximum of activity around 1999–2002 and the minimum of activity around 2006–2009. Although the use of averaged a-coefficients makes it difficult to evaluate the extension of the activity zone, the model successfully retrieves the mean latitude of activity during the maximum of the solar cycle.

We then applied the method to the brightest stars observed during the initial observational phase of the Kepler space instrument: 16 Cyg A and B. These stars were selected as a test bed due to the fact that they are well studied and present the highest mode signal-to-noise ratio of all the currently known main sequence stars. Davies et al. (2015) suggested that these stars have mild to no activity. However, Bazot et al. (2019), using a parametric model for describing the asphericity (R_eq − R_pol)/R_eq, found an asphericity significant at 1σ. Our current analysis, performed using the same data set, confirms this asphericity, and we found a mild (relative to the Sun) to moderate activity for both stars. It is found that 16 Cyg A has a near-equatorial band of activity during the period of observation (13 September 2010 to 8 March 2013), with a significance of the detection greater than 79.8%.

The case of 16 Cyg B is more ambiguous. A bimodality on the average splitting ⟨δν_nlm/m⟩≃a₁ and on the stellar inclination has already been reported in Davies et al. (2015) and Bazot et al. (2019). Our refined model suggests that this bimodality is in fact related to the l = 2 modes. Indeed, as we accounted for a₁, a₂, a₃, and a₄, we noted that the bimodality previously seen on a₁ in earlier studies becomes displaced to a₄. The solutions of a₄ were separated using a Gaussian process algorithm, and we found that the weight (or importance) of the highest solution, close to 0 nHz, is of 70%. The lower solution, close to −28 nHz, has a weight of 30%. The associated probability distribution for stellar inclination i = 35 ± 3 is unique (instead of being bimodal in past studies). We analysed the two solutions of a₄ independently. The lower solution (with lower weight) is associated with an overall activity that is stronger than the Sun, localised at latitudes of approximately 32° (with an uncertainty of ≃3°). The higher solution is linked to a lower activity and is weakly significant (probability greater than 67.1%). Although the uncertainty is large, our study suggests an activity closer to the equatorial region. In the Sun, and as evident in its butterfly diagram, the quiet phase is associated with a magnetic activity closer to the equator, while the transition to the active Sun is abrupt and characterised by the appearance of magnetic spots at latitudes of 30–40°. In that context, and although it is not possible to rule out the possibility of a statistical fluke, an interpretation of the bimodality is that the star was transitioning from a period of low activity to a more active period during the observation time of the Kepler instrument. To evaluate that hypothesis, we selected solar data between January 2006 and January 2011, which includes the end of a cycle and the start of a new one. The measured a-coefficients for that analysis are shown in Fig. 8. There is no visible bimodality on either a₂ or a₄. Because 16 Cyg B is evidently different than the Sun, this does not refute the hypothesis but weaken it. A firmer verification would require a follow-up observation of the star, with PLATO (Rauer et al. 2014) for example, or extensive simulations in order to attempt to reproduce the bimodality.

This work demonstrates that it is possible to determine the latitude and intensity of the activity of Sun-like stars when this activity is similar to or exceeds that of the Sun. The limited observed bias suggests that this analysis of the stellar activity could also be extended to a larger set of stars. Notably, the Kepler LEGACY stars (Lund et al. 2017) are ideal candidates, as these stars generally satisfy the criteria of reliability discussed in this conclusion and detailed in Sect. 4.

Although the a-coefficient analysis has the benefit of allowing for simplification of the bias studies and making the analysis faster, a further axis of improvement would consist of using a single analysis step, for example, by fitting the power spectrum with the model of activity directly. Contrary to the two-step analysis, such an approach does not require assumptions on the properties of probability distribution of the a-coefficients, as is currently the case. It would also make use of the full set of a-coefficients (not only their average), which may lead to smaller uncertainties. However, the two approaches would certainly need to be used in tandem, as the stability of the solution (and its accuracy) is harder to ascertain when performing a single-step fitting approach.

Acknowledgments

O.B. and L.G. acknowledge partial support from NYUAD Institute Grant G1502 “Center for Space Science”. L.G. is supported in part by the European Research Council (ERC) Synergy Grant WHOLE SUN 810218 and the Max Planck Society grant “Preparations for PLATO Science”. We thank Saskia Hekker for comments and advice.

References

Appendix A: Polynomials, a-coefficients and splittings

A.1. The first six used polynomials

The ${\mathcal{P}}_j^{(l)}(m)$ coefficients, which were introduced by Schou et al. (1994), can be obtained by normalising the $H_j^m(l)$ coefficients by Ritzwoller & Lavely (1991) so that ${\mathcal{P}}_j^{(l)}(m=l)=l$.

$$\begin{array}{cc}\hfill {\mathcal{P}}_j^{(l)}(m)& =\frac{l(2l-j)!}{(2l)!}{H}_j^m(l)\hfill \\ \hfill & ={(-1)}^{-l+m}\frac{l\sqrt{(2l-j)!(2l+j+1)!}}{(2l)!}\left(\begin{array}{ccc}j& l& l\\ 0& m& -m\end{array}\right),\hfill \end{array}$$(A.1)

where the last factor means the Wigner 3j symbol. The coefficients can also be obtained by recurrence, using Equations (2) and (3) together with the normalisation condition ${\mathcal{P}}_j^l(m)=l$ and starting with ${\mathcal{P}}_0^l(m)=l$. The first six polynomials are as follows:

$$\begin{array}{c}\hfill {\mathcal{P}}_1^{(l)}(m)=m,\end{array}$$(A.2)

$$\begin{array}{c}\hfill {\mathcal{P}}_2^{(l)}(m)=\frac{3m^2-l(l+1)}{2l-1},\end{array}$$(A.3)

$$\begin{array}{c}\hfill {\mathcal{P}}_3^{(l)}(m)=\frac{5m^3-(3l(l+1)-1)m}{(l-1)(2l-1)},\end{array}$$(A.4)

$$\begin{array}{c}\hfill {\mathcal{P}}_4^{(l)}(m)=\frac{(35m^4-5(6l(l+1)-5)m^2)+3l(l+1)(l(l+1)-2)}{2(l-1)(2l-1)(2l-3)},\end{array}$$(A.5)

$$\begin{array}{c}\hfill {\mathcal{P}}_5^{(l)}(m)=\frac{252m^5-140(2L-3)m^3+(20L(3L-10)+48)m}{8(4l^4-20l^3+35l^2-25l+6)},\end{array}$$(A.6)

$$\begin{array}{cc}\hfill {\mathcal{P}}_6^{(l)}(m)& =\frac{924m^6-420m^4(3L-7)+84m^2(5L^2-25L+14)}{N}\hfill \\ \hfill & -\frac{20L(L^2-8L+12)}{N},\hfill \end{array}$$(A.7)

where L = l(l + 1) and N = 64l⁵ − 480l⁴ + 1360l³ − 1800l² + 1096l − 240.

A.2. Relationship between a-coefficients and splittings

The relationships between a-coefficients and symmetric or antisymmetric splittings can be derived using Equations (4) and (5). The equations are used to convert model frequencies ν_nlm into a-coefficients.

A.2.1. Explicit form for low a-coefficients and l=1

The relation between a-coefficients a₁(n, l = 1),  a₂(n, l = 1) and the modes frequencies is,

$$\begin{array}{c}\hfill a_1(n,l=1)=S_{n11}=\frac{{\nu }_{n,1,1}-{\nu }_{n,1,-1}}{2}\end{array}$$(A.8)

$$\begin{array}{c}\hfill a_2(n,l=1)=T_{n11}/3=\frac{({\nu }_{n,l,-1}+{\nu }_{n,l,1})/2-{\nu }_{n,l,0}}{3}\end{array}$$(A.9)

A.2.2. Explicit form for low a-coefficients and l=2

The relation between a-coefficients a₁(n, l = 2), a₂(n, l = 2), a₃(n, l = 2), a₄(n, l = 2) and the symmetric or antisymmetric splittings is,

$$\begin{array}{c}\hfill a_1(n,l=2)=\frac{S_{n21}+4S_{n22}}{5}\end{array}$$(A.10)

$$\begin{array}{c}\hfill a_2(n,l=2)=\frac{2T_{n22}-T_{n21}}{7}\end{array}$$(A.11)

$$\begin{array}{c}\hfill a_3(n,l=2)=\frac{S_{n22}-S_{n21}}{5}\end{array}$$(A.12)

$$\begin{array}{c}\hfill a_4(n,l=2)=\frac{T_{n22}-4T_{n21}}{70}\end{array}$$(A.13)

A.2.3. Explicit form for low a-coefficients and l=3

The relation between a-coefficients  a₁(n, l = 3),  a₂(n, l = 3),  a₃(n, l = 3),  a₄(n, l = 3),  a₅(n, l = 3),  a₆(n, l = 3) and the symmetric or antisymmetric splittings is,

$$\begin{array}{c}\hfill a_1(n,l=3)=\frac{S_{n31}+4S_{n32}+9S_{n33}}{14}\end{array}$$(A.14)

$$\begin{array}{c}\hfill a_2(n,l=3)=\frac{-15T_{n31}+25T_{n33}}{126}\end{array}$$(A.15)

$$\begin{array}{c}\hfill a_3(n,l=3)=\frac{-S_{n31}-2S_{n32}+3S_{n33}}{9}\end{array}$$(A.16)

$$\begin{array}{c}\hfill a_4(n,l=3)=\frac{T_{n31}-7T_{n32}+3T_{n33}}{77}\end{array}$$(A.17)

$$\begin{array}{c}\hfill a_5(n,l=3)=\frac{5S_{n31}-8S_{n32}+3S_{n33}}{126}\end{array}$$(A.18)

$$\begin{array}{c}\hfill a_6(n,l=3)=\frac{15T_{n31}-6T_{n32}+T_{n33}}{1386}\end{array}$$(A.19)

Appendix B: Spectrum analysis

B.1. Acoustic spectrum modelling

The analysis of the asteroseismic data is very often performed by fitting the power spectrum (e.g. Appourchaux et al. 1998, 2008). This requires devising a likelihood function and a model for the observed spectrum. As our fitting involves a Bayesian approach, the priors are also required.

B.1.1. Likelihood

The likelihood function is determined by the noise statistics of the power spectrum, which is an χ-squared with two degrees of freedom. It is also based on the assumption that the frequency bins are independent and uncorrelated. This implies that the observation duration is assumed to be much greater than the typical lifetime of the pulsation modes and that the duty cycle is sufficient to allow us to neglect any leakage induced by the window function. The observation time frame for the LEGACY sample is of the order of years, which is significantly longer than the mode lifetime. Furthermore, the duty cycle is above 95%, which was shown to be sufficient to neglect leakage (Stahn 2010). Under these conditions, the likelihood function is (Anderson et al. 1990)

$$\begin{array}{c}\hfill L_1(S({\nu }_i)|\mathit{X},M)={\displaystyle \sum _{i=0}^{N-1}}\frac{1}{S({\nu }_i)}exp(-\frac{S({\nu }_i)}{M({\nu }_i,\mathit{X})}),\end{array}$$(B.1)

where S(ν_i) is the power at the central frequency ν_i of the i^th bin. Here, N is the total number of bins, and M(ν_i, X) is the model of the power spectrum with the set of variable X for the model M.

From Bayes theorem, the posterior probability density function is defined as

$$\begin{array}{c}\hfill {\pi }_1(\mathit{X}|S(\nu ),M)=\frac{{\pi }_1(\mathit{X}|M)L_1(S(\nu )|\mathit{X},M)}{{\pi }_1(S(\nu )|M)},\end{array}$$(B.2)

where π₁(X|M) is the prior knowledge on the parameters X. We note that π₁(S(ν)|M) is a normalisation constant and is essential only when comparing the significance of models (e.g. Gregory 2005; Benomar et al. 2009). We assumed the variables to be independent from each other. We evaluated the posterior using the Tempered Adaptive MCMC code¹⁰ described in Atchadé (2006) and implemented by Benomar (2008).

B.1.2. Acoustic spectrum model

The model used for the power spectrum fitting involves a sum of asymmetrical Lorentzian superimposed to a monotonically decreasing function of frequency (pink noise). The asymmetric Lorentzian is commonly used in helioseismology (Duvall et al. 1993; Nigam & Kosovichev 1998; Georgobiani et al. 2000; Toutain et al. 1998) and in asteroseismology for Sun-like stars (Benomar et al. 2018b). The sum of Lorentzian profiles is defined as,

$$\begin{array}{c}\hfill M(\nu )={\displaystyle \sum _{n=n_0}^{n_{\mathit{max}}}}{H}_{\mathit{nlm}}\frac{{(1+B_{\mathit{nlm}}z)}^2+B_{\mathit{nlm}}^2}{1+z^2}+N(\nu ),\end{array}$$(B.3)

with z = 2(ν − ν_nlm)/Γ_nlm. Each asymmetrical Lorentzian is defined by a height H_nlm, width Γ_nlm, central frequency ν_nlm, and asymmetry B_nlm. As explained by Gizon (2006); Benomar et al. (2018b), the asymmetry coefficient depends on the mode width and frequency. The normalised asymmetry coefficient ${\chi }_{\mathit{nlm}}=2\frac{B_{\mathit{nlm}}{\nu }_{\mathit{nlm}}}{{\mathrm{\Gamma }}_{\mathit{nlm}}}$ is fitted, as it is nearly constant over the range of fitted modes.

Several prescriptions exist for describing the noise background N(ν). In this work, we assumed it to be a sum of two generalised Lorentzian, sometimes referred as Harvey-like profiles (Harvey 1985), and of a white noise:

$$\begin{array}{c}\hfill N(\nu )=\frac{A_1}{1+{({\tau }_1\nu )}^{p_1}}+\frac{A_2}{1+{({\tau }_2\nu )}^{p_2}}+N_0.\end{array}$$(B.4)

Here, N₀ is the white noise, and A_k is the maximum heights of the k^th generalised Lorentzian. The τ_k parameter is the timescale that is the inverse of the full width at half maximum of the Lorentzians, and p_k is a power exponent.

The fit is performed globally over all the statistically significant peaks visible in the power spectrum of a star. This can lead to a very large number of fitted parameters (up to a few thousands), which is practically unsuitable. A model simplification is therefore preferred (similar to e.g. Appourchaux et al. (2008), Benomar et al. (2009), Campante et al. (2011), Handberg & Campante (2011)); the main simplifications recalled hereafter:

• The m dependence on heights is controlled by its relationship with the stellar inclination (Gizon & Solanki 2003). This saves several hundreds of parameters.

• The relative height of the different degree l is constant across the fitted range and for a given l, hence $V_l^2=H_{n,l}/H_{n,l=0}=\mathit{const}$. Only the l = 0 heights H_{n, l = 0} are variables. The mode visibility $V_l^2$ replaces the H_{n, l > 0} as a variable.

• Given a degree l, Γ_n, l, m = Γ_n, l is imposed. This is justified by the fact that the width depends weakly on the frequency, as all split components are assumed to have the same width.

• Because Γ_n, l is nearly independent of the degree, it is possible to fit Γ(ν) = Γ_{n, l = 0} and to interpolate it to the frequencies of the modes with degree l > 0.

These assumptions reduce the number of variables to a few tens in the case of CoRoT, Kepler, or TESS observations. The frequencies of the modes follow Equation (14) in Section 2.3.

B.2. Priors

Priors are fundamental to the Bayesian method. This section discusses π₁(X|M, I), the prior used during the power spectrum fitting. We assumed the parameters of the vector X to be independent of each other such that the product rule is used to the define π₁(X|M, I).

Heights, widths, frequencies, asymmetry, and inclination use non-informative priors. These are either Jeffreys priors for scale parameters, noted as 𝒥(x_min, x_max), or uniform priors, noted as 𝒰(x_min, x_max). The priors on frequencies are uniform and require a visual inspection of the power spectrum in order to assign the lower and upper bound of the prior for each mode that follows the expected pattern for main sequence Sun-like stars (equally spaced p modes) and that shows an excess of power relative to the background that exceeds 80%. The excess of power is determined using a smooth spectrum for which the noise statistics is derived by Appourchaux (2003). Mode visibilities are defined by Gaussian priors (noted as 𝒢(x₀, σ)), with the mean x₀ set as the solar value (see Ballot & Barban 2011), and the standard deviation σ is 10% of the mean. Table B.1 lists the type and the prior characteristic values that are used for the parameters of the modes.

The prior on a₁ is uniform between 0 and 1500 nHz. For 16 Cyg A and B, a uniform fixed prior is set on |a₃| over the range [0, 100] nHz. However at each iteration of the optimisation process, |a₃| is not allowed to exceed 20% of a₁. Preventing an extremely large relative value of a₃/a₁ improves the fit stability and ensures a faster convergence of the algorithm. For the Sun analysis and because it is not possible to measure a₃ for i ≃ 90° (due to the lack of amplitude of l = 2, m = ±1 at that inclination), a₃ is fixed to zero.

The priors on a₂ and a₄ are also uniform. The range is defined by using the maximum range of Figure 4 (showing $a_2^{(AR)}$ and $a_4^{(AR)}$, increasing it by 50%, and adding the expected centrifugal term $a_2^{(CF)}$, assuming a₁ = 400nHz for the Sun and a₁ = 600nHz for 16 Cyg A or 16 Cyg B. The Δν reported in Table 3 is also used. For the same reason as with a₃, a₂/a₁ cannot exceed 50%, and a₄/a₁ cannot exceed 20% at each iteration step of the optimisation process.

The noise priors are obtained from a global MAP approach similar to Benomar et al. (2012). This provides the best fit values and 1σ uncertainties that are to be used as priors. The noise background of the individual model analysis is described by Equation (B.4). During the global MAP fit, the model is made of that same background model plus a Gaussian envelope to account for the power excess due to the modes (e.g. Mathur et al. 2010; Huber et al. 2011).

Appendix C: Inference of the active latitudes

The determination of the activity is performed as a second step using the marginalised posterior distribution of the average a-coefficients as observables. This section describes the method and assumptions for determining the posterior distribution and computing the significance of the detection.

C.1. Posterior distribution and likelihood

Similar to the power spectrum fitting described in Appendix B, the determination of the most likely latitudes for the activity and its significance necessitate the computation of a posterior distribution,

$$\begin{array}{c}\hfill {\pi }_2(\mathit{X}|\mathit{O},M)=\frac{{\pi }_2(\mathit{X}|M)L_2(\mathit{O}|\mathit{X},M)}{{\pi }_2(\mathit{O}|M)},\end{array}$$(C.1)

with π₂(X|O, M), π₂(X|M), and L₂(O|X, M) as the posterior distribution, the prior, and the likelihood, respectively. The denominator, π₂(O|M), is the normalisation constant used for model comparison (marginal likelihood). The posterior distributions of the a-coefficients a_2, o and a_4, o obtained by power spectrum fitting are the observables contained in O. The class of the model is identified by the variable M.

For simplicity, the observables are assumed to be distributed according to uncorrelated Gaussian functions such that O = {a_2, o, a_4, o, σ_2, o, σ_4, o}, where a_2, o and a_4, o are the mean for the distributions of a₂ and a₄, while σ_2, o and σ_4, o are the standard deviations. This leads to a log-likelihood in the form of an χ-squared,

$$\begin{array}{c}\hfill ln{L}_2(\mathit{X}|\mathit{O},M)=-\frac{{(a_{2,o}-a_{2,m}(\mathit{X}))}^2}{2{\sigma }_{2,o}^2}-\frac{{(a_{4,o}-a_{4,m}(\mathit{X}))}^2}{2{\sigma }_{4,o}^2}.\end{array}$$(C.2)

In Dziembowski et al. (2000), a-coefficients are weighted using the inverse of the mode inertia (see their Equation 4). This is due to the fact that the inertia account for most of the frequency variations of the coefficients. Accounting for the mode inertia is essentially relevant if the uncertainties are small enough to observe a trend when modes are fitted individually. However, in the case of the Sun and with either GOLF or BiSON data, it is difficult to perceive a frequency trend (Chaplin et al. 2003). Our own trials on VIRGO/SPM using a linear fit to describe the frequency dependence of an a-coefficient did not detect a slope that is significant at more than 1σ. A similar analysis on 16 Cyg A and B and on simulations showed evidence of large uncertainties when attempting to determine frequency variations for stars with HNR typical of Kepler observations. In fact, only the average of the a-coefficient over l and n is shown to be robustly determined (see also our discussion Section 3.3) so that inertia effects are neglected.

The activity model M_AR depends on the variables X = {ϵ, θ₀, δ} and a_2, m(X), whereas a_4, m(X) denotes the modelled a-coefficients. The log-likelihood is then

$$\begin{array}{cc}& ln{L}_2(a_{2,o},a_{4,o},{\sigma }_{2,o},{\sigma }_{4,o}|a_{2,m},a_{4,m},M_{\mathit{AR}})=\hfill \\ \hfill & -\frac{{(a_{2,o}-a_{2,m}(ϵ,{\theta }_0,\delta ))}^2}{2{\sigma }_{2,o}^2}\hfill \\ \hfill & -\frac{{(a_{4,o}-a_{4,m}(ϵ,{\theta }_0,\delta ))}^2}{2{\sigma }_{4,o}^2}.\hfill \end{array}$$(C.3)

In our case, the observables are the average of the fitted modes’ coefficients, $a_{2,o}=a_2^{(CF)}+a_2^{(AR)}$ and $a_{4,o}=a_4^{(AR)}=0$. However, to conveniently propagate all errors on the parameters of $a_2^{(CF)}$ (a₁, Δν, Δν_⊙, ρ_⊙, ν_nl), it is preferable to use $a_{2,o}=a_2^{(AR)}$, which is obtained by subtracting the $a_2^{(CF)}$ from the a₂ measured by the power spectrum fitting. The distribution of Δν is computed by a linear fitting of the groupwise ensemble of samples (n, l = 0) for each measured frequency ν_n0. Finally, the term $a_2^{(CF)}$ is obtained through a weighted average of all independent $a_2^{(CF)}(n,l)$ computed at the star’s posterior frequencies ν_nl. The subtraction of $a_2^{(CF)}$ is then again performed using the samples of its posterior, enabling the posterior probability distribution function of $a_2^{(AR)}$ to be constructed, from which we deduced its mean a_2, o and its standard deviation σ_2, o. To evaluate the relevance of the activity, two models are considered. First, a model (M = M_CF) without activity that accounts only for the centrifugal effects and, second, a model (M = M_AR) with activity. Their details and the choice of the priors are described in the following sections individually.

C.2. Model without activity M_CF

The model M_CF has no unknown variable requiring minimisation (see Equation (13)). For the sake of comparing this model with the model M_AR, it is however important to determine the marginal likelihood of M_CF, and it is then necessary to calculate the denominator of Equation (C.1), P(O|M_CF), which is defined as an integral,

$$\begin{array}{c}\hfill P(\mathit{O}|M_{\mathit{CF}})={\displaystyle \int }\pi (\mathit{X}|M_{\mathit{CF}})L_2(\mathit{O}|a_{2,m},a_{4,m},M_{\mathit{CF}})d\mathit{X}.\end{array}$$(C.4)

An absence of activity corresponds to the limit case where the a-coefficients of the activity are exactly zero. This corresponds to setting Dirac priors on a_2, m = δ₀ and a_4, m = δ₀ and implies that the marginal likelihood is the local value of the likelihood at X = {0, 0}. The model comparison is commonly performed in log-space such that

$$\begin{array}{c}\hfill lnP(a_{2,o},a_{4,o},{\sigma }_{2,o},{\sigma }_{4,o}|M_{\mathit{CF}})=-\frac{a_{2,o}^2}{2{\sigma }_{2,o}^2}-\frac{a_{4,o}^2}{2{\sigma }_{4,o}^2}.\end{array}$$(C.5)

C.3. Model with activity M_AR

Regarding the model with activity of Equation (C.3), priors must be set on the variables subject to optimisation. The ϵ and δ parameters are indispensable intensive parameters, and the adequate non-informative prior is then the (truncated) Jeffreys prior (Jeffreys 1961), uniform in the log-space,

$$\begin{array}{c}\hfill J(x)=\frac{ln(1+x_{\mathit{max}}/x_{\mathit{min}})}{x+x_{\mathit{min}}},\end{array}$$(C.6)

where x_min and x_max are the upper and lower bounds such that if x > x_max or x < x_min, then J(x) = 0. This guarantees that the probability density is proper (the integral over x is finite). In case of weak information content in the observables, this prior is more weighted towards a null value. We set ϵ_min = 5.10⁻⁴ and ϵ_max = 10⁻². This embraces the solar value ϵ ≃ 5.10⁻⁴ (Gizon 2002). Because δ ≃ 10° at the maximum of solar activity, we set δ_min = 10°. We also note that consistent with the discussion of Figure 3 in Section 3, δ_max = 45° is required. Finally, the location parameter θ₀ has a uniform prior in the range [0, 90]( ° ).

The marginal likelihood of the model M_AR requires evaluation of the triple integral,

$$\begin{array}{c}\hfill P(\mathit{O}|M)={\displaystyle \int }\pi (ϵ)\pi (\theta )\pi (\delta )L_2(\mathit{O}|a_{2,m},a_{4,m},M)dϵd\theta d\delta .\end{array}$$(C.7)

The MCMC process used here involves the use of parallel Metropolis-Hasting tempered chains. The chains mix with each other in order to enhance the sampling. As explained in Benomar et al. 2009 (see their Section A.3), these parallel chains can be used to approximate the Equation (C.7). This technique is used with ten parallel chains following a geometrical temperature law T_k = 1.7^k − 1, with k as the chain index such that k = 1 is the target distribution, given by equation (C.1).

Appendix D: Bias map for $\widehat{\mathit{HNR}}=10$ and $\widehat{\mathit{HNR}}=20$

This section shows the bias on the a₁, a₂, and a₄ coefficients in the case of $\widehat{\mathit{HNR}}=10$ and $\widehat{\mathit{HNR}}=20$. Figures D.1 and D.3 are for the case of an equatorial band of activity, and Figures D.2 and D.4 are for a polar cap. The inaccuracy of the fit remains noticeably smaller than the standard deviation in the majority of the parameter space and for i > 30°, a₁/Γ_νmax > 0.4. Interestingly, the relative-to-error bias on a₁ is more pronounced in the case of an equatorial band of activity when {a₁/Γ_νmax ≤ 0.4, i = [30, 60]} (Figure D.1d and D.3d). Then this bias can exceed three times the uncertainty. The accuracy of the inference of activity based on the a-coefficients is expected to weakly depend on the $\widehat{\mathit{HNR}}$ (although potentially with a large uncertainty) for most stars, except within this regime (in addition to the grey area) where the centrifugal effect will be overestimated.

Fig. D.1. Bias analysis for $\widehat{\mathit{HNR}}=10$ for an equatorial activity band (θ0 = 85,δ = 10) of similar intensity to the Sun (ϵnl = 5.10−4) for Tobs = 2 years (top) and Tobs = 4 years (bottom).

Fig. D.2. Bias analysis for $\widehat{\mathit{HNR}}=10$ for a large polar activity cap (θ0 = 22.5°,δ = 45°) of similar intensity to the Sun (ϵnl = 5.10−4) for Tobs = 2 years (top) and Tobs = 4 years (bottom).

Fig. D.3. Bias analysis for $\widehat{\mathit{HNR}}=20$ for an equatorial activity band (θ0 = 85,δ = 10) of similar intensity to the Sun (ϵnl = 5.10−4) for Tobs = 2 years (top) and Tobs = 4 years (bottom).

Fig. D.4. Bias analysis for $\widehat{\mathit{HNR}}=20$ for a large polar activity cap (θ0 = 22.5,δ = 45) of similar intensity to the Sun (ϵnl = 5.10−4) for Tobs = 2 years (top) and Tobs = 4 years (bottom).

Appendix E: Figures of the best fits

This section presents visuals of the best MCMC fits for the power spectrum of the active Sun (Figure E.1), quiet Sun (Figure E.2), 16 Cyg A (Figure E.3), and for the two a₄ solutions of 16 Cyg B (Figures E.4 and E.5). For 16 Cyg B, the solutions were separated by selecting the median for the samples only below or above a₄ = 21nHz. In the case of the Sun, the impact of the activity on the profile of l = 2 is very clear, as the shifts of the m-components introduce an asymmetry in power for the overall l = 2 mode profile during the active phase. During the quiet phase, such an asymmetry is not visible.

Fig. E.1. Solar power spectrum for observation between 1999 and 2002 after Gaussian smoothing over a kernel width s = 50nHz. Superimposed is the best MCMC fit. Top panel: Highlight of a l = 2, 0, 3, 1 mode group (left to right) and their fit (red). The power asymmetry due to non-symmetric m-components is visible in the data. Top panel, inset: Overall view of the modes. Bottom panel: Residual of the fit with two levels of Gaussian smoothing.

Fig. E.2. Same as Figure E.1 but between 2006 and 2009. Top panel: Highlight of a l = 2, 0, 3, 1 mode group (left to right) and their fit (red). The power asymmetry seen for 1999-2002 is not apparent in the data. Top panel, inset: Overall view of the power spectrum. Bottom panel: Residual of the fit with two levels of Gaussian smoothing. The residual on l = 1 is a bit high due to the visibility being fixed for all modes.

Fig. E.3. Same as Figure E.1 but for 16 Cyg A. Top panel: Highlight of a l = 2, 0, 1 mode group (left to right) and their fit (red). A very mild power asymmetry is seen in the l = 2 data. Top panel, inset: Overall view of the power spectrum. Bottom panel: Residual of the fit with two levels of Gaussian smoothing. The residual shows an excess of power due to the low HNR≃1.7 l = 3 (not fitted here).

Fig. E.4. Same as Figure E.1 but for 16 Cyg B, a4 < −21 nHz. Top panel: Highlight of a l = 2, 0, 1 mode group (left to right) and their fit (red). A mild power asymmetry is seen in the l = 2 data. Top panel, inset: Overall view of the power spectrum. Bottom panel: Residual of the fit with two levels of Gaussian smoothing. The residual shows an excess of power due to the low HNR≃1.6 l = 3 (not fitted here).

Fig. E.5. Same as Figure E.1 but for 16 Cyg B, a4 > −21 nHz. Top panel: Highlight of a l = 2, 0, 1 mode group (left to right) and their fit (red). Top panel, inset: Overall view of the power spectrum. Bottom. Residual of the fit with two level s of Gaussian smoothing. The residual show an excess of power due to the low HNR≃1.6 l = 3 (not fitted here).

Similarly, although it is less pronounced than in the Sun, in the cases of 16 Cyg A and B, a weak asymmetry is visible in l = 2 modes. For 16 Cyg B, the two solutions of a₄ provide sensibly similar mode profiles, as the difference is weakly apparent to the eye only for l = 2. The residuals do not show striking differences. This visually supports the fact that the goodness of fit is similar and explains that the two solutions coexist in statistical terms.