© ESO 2019

1. Introduction

Rotation is a ubiquitous characteristic of stars; it has many connections to convection, stellar pulsations, and magnetic fields (e.g. Tassoul 2000). It can also be related to the evolution of close-by planets through tidal effects (e.g. Privitera et al. 2016). Despite this central role that it plays in stellar physics, a full understanding of its behaviour remains elusive. Only for the Sun are measurements precise enough to provide a proper insight on rotation (e.g. Thompson et al. 2003), owing to the large number of observed solar pulsation eigenmodes (Hill et al. 1996). The spherical-harmonic degrees of the detected modes reach values as high as l ≃ 200 at low radial orders (Larson & Schou 2008). Higher values, up to ≳1500, can be reached by fitting ridges in the power spectrum (Couvidat 2013). In a spherically symmetric, non-rotating star, the eigenfrequencies of the non-radial modes are degenerate with respect to the azimuthal order. However, perturbations to the velocity field caused by the rotational flow can lift this degeneracy. The resulting frequency splitting is related to the rotational flow by a linear integral equation. It is possible to invert this relation to obtain information on the solar rotation rate (see e.g. Thompson 1991). In contrast, the degrees of the eigenmodes observed in other stars are typically in the range l = 0 − 2, with some observations detecting l = 3 modes (Bouchy & Carrier 2001; Bazot et al. 2012; Appourchaux et al. 2012; Metcalfe et al. 2012). As a consequence, inverting the eigenfrequencies becomes far more challenging.

There nevertheless exist sources of information on rotation in stars other than Sun. Spectroscopic measurements offer estimates of the surface velocity as it broadens absorption lines through the Doppler effect. In general, the velocity is known up to a sin i factor, with i the inclination of the spin axis of the star with respect to the observer’s line of sight; this cannot be estimated from spectroscopic measurements alone. Recently, however, progress has been made through the advent of asteroseismology, and in particular, through the space missions CoRoT (Baglin et al. 2009) and Kepler (Borucki et al. 2010).

Improvements have come in two ways. First, it is now possible to measure frequency splittings with precision (Gizon et al. 2013; Appourchaux et al. 2014; Nielsen et al. 2014). These relate to the average rotation rate of a star. Furthermore, the inferred value of the rotational splitting is geometrically tied to the inclination of the star (Gizon & Solanki 2003) through its effect on the observed mode amplitudes. Adequate estimates of the former therefore require estimates of the latter. Combined with spectroscopic measurements and radius determinations, splitting and inclination determinations have provided consistency checks of the asteroseismically estimated surface rotational velocity (Gizon et al. 2013). It has also been argued that discrepancies between the two measurements can be interpreted as a signature of radial differential rotation (Benomar et al. 2015; Deheuvels et al. 2015).

The other major source of information on stellar rotation rates comes from photometric light curves and is related to stellar activity. It is well known that spots and plages transiting the stellar surface induce drops and increases in the integrated flux, respectively. If the activity signal is sufficiently long and coherent, this produces a quasi-periodic modulation of the intensity that is correlated with the stellar rotation rate. Such signatures have been extensively studied for large samples of Kepler stars (Reinhold et al. 2013; Walkowicz & Basri 2013; García et al. 2014; McQuillan et al. 2014; Nielsen et al. 2013). Some caution is needed with the interpretation of these rotation rates because they depend on the latitude of the spots through the variation with latitude of the surface rotation rate.

Latitudinal differential rotation is indeed a major characteristic of the solar rotation profile. It is well established that the Sun rotates faster at the equator than at the poles by a factor of roughly 1.4. It has also been shown (Schou et al. 1998; Chaplin et al. 1999; Thompson et al. 2003) that this latitudinal differential rotation remains approximately constant as a function of radius throughout the convective zone. In the radiative zone, at least outside its inner 20%, the rotation rate becomes that of a solid body. For this reason, latitudinal differential rotation is believed to be the result of the interplay between rotation and convective flows.

Much uncertainty currently remains as to which convective scale is the main driver of this phenomenon. On the one hand, it has been speculated that a mean-field approach of turbulent convection can explain differential rotation. The basic picture consists of describing the convective flow as a stationary component plus a time-dependent turbulent term that after insertion in the Navier-Stokes equations for a rotating fluid, gives rise to Reynolds stresses. The latter is in turn related to the radial and latitudinal gradients of the rotation rates, that is, to differential rotation (Rüdiger 1974, 1980, 1982). This type of model can reproduce the observations made for the Sun in a rather satisfactory way (Kitchatinov & Rüdiger 1995; Kitchatinov & Nepomnyashchikh 2017). On the other hand, the global approach consists of solving the full set of the equations of hydrodynamics using the proper spherical coordinates and boundary conditions. This work was initiated by Gilman & Glatzmaier (1981). Given the spherical arrangement of the system, differential rotation is caused by the Coriolis force acting on large-scale convective motion (Tassoul 2000). Recent work has allowed modelling enough spatial convective scales so that these simulations can reproducethe solar observations to some extent (Guerrero et al. 2013, 2016). These two explanations need not be mutually exclusive. Finally, hydrodynamic simulations for different stellar masses and rotation rates have been carried out (see e.g. Brun et al. 2017). These show a wide range of morphologies for rotation flows in Sun-like stars. It is therefore clear that observing differential rotation is key to understanding the interplay between convection and rotation. It is our goal to provide such information for Sun-like stars.

Our refined picture of differential rotation in the Sun is due to the large number of modes observed in this star, which greatly helps the inversion process. The turning point of a given mode in the stellar interior depends on its degree. Consequently, the wide range of l values obtained for the Sun allows us to probe many layers of its internal structure. Until recently, the relatively low number of frequencies detected in other stars (typically about 20 modes) had been an obstacle to inversion. Information on latitudinal differential rotation was gathered from three main sources. First, it is possible to invert the relation described above between the spot rotation rates and differential rotation to infer the latter (Donahue et al. 1996; Reinhold & Reiners 2013; Lanza et al. 2014; Davenport et al. 2015). This is done through spot modelling (e.g. Dumusque et al. 2014), which involves an implicit physical model for the stellar spots. In the same spirit, some studies have tried to identify variations in the measurements of CaII emission lines, which are tied to active stellar regions (Bertello et al. 2012). Interesting peculiar cases of photometric spot modelling are encountered for planet-hosting stars when the companion transits in front of stellar spots. The resulting decrease in transit depth may allow a precise characterisation of differential rotation (Valio et al. 2017). Another method involves the use of Doppler maps obtained from spectroscopy or spectropolarimetry (e.g. Donati & Collier Cameron 1997; Marsden et al. 2011). The last approach available consists of analysing the Fourier transforms of spectral absorption lines, whose side lobes are sensitive to the variation in latitude of the rotation rate (Huang 1961; Gray 1977). This strategy has been used extensively to detect differential rotation in A- and F-type stars (Reiners & Schmitt 2002a,b, 2003; Ammler-von Eiff & Reiners 2012). In a recent study, however, Benomar et al. (2018) demonstrated the possibility of measuring the magnitude of differential rotation in 13 Sun-like stars using the techniques of asteroseismology on Kepler time series. We extend this work here.

We report seismic detection of latitudinal differential rotation for 16 Cyg A and B. These stars have physical characteristics close to those of the Sun; in particular, the rotation rates obtained from spectroscopy suggest roughly solar values (Takeda et al. 2005). For this reason, they have sometimes been dubbed “solar analogues”. It is therefore of prime interest to assess whether the measured latitudinal differential rotation is also close to solar-like or if a difference exists. Most of the stars with measured differential rotation reported in Benomar et al. (2018) are leaning towards the F type. A handful of them have effective temperatures closer to the Sun, but still hotter by roughly 300 K. Therefore, 16 Cyg A and B are the best candidates to study rotation under nearly solar conditions, which is supported by the occasional classification of these stars as “solar twins” (e.g. King et al. 1997 ; see also the introduction of Bazot et al. 2018 for a note on the class of solar twins).

Another important fact about 16 Cyg A and B is that they are some of the brightest stars that the Kepler mission observed. Therefore, their oscillation frequencies were estimated with extremely good precision (Metcalfe et al. 2012; Davies et al. 2015). This gives us an opportunity to address the problem of how the estimates of latitudinal differential rotation depend on the uncertainties in the stellar parameters (mass, age, initial chemical composition, and mixing-length parameter).

In Sect. 2 we describe the general procedure for data fitting, including the modelling of differential rotation effects. We present its outcome in terms of coefficient estimates for a functional expansion of the rotational splitting. In Sect. 3 we explain the inversion methods used for both stars, with an emphasis on the more difficult case of 16 Cyg B. We discuss in Sect. 4 the impact of the uncertainties on the stellar parameters on our results. We also explore the implication of the measured stellar deformation on the magnetic field of the stars.

2. Acoustic mode fitting with differential rotation

2.1. Data

The sources 16 Cyg A and B have been observed by Kepler from 13 September 2010 to 8 March 2013 (covering Kepler quarters Q7 to Q16). Their magnitudes, V = 5.95 and 6.20, respectively place them beyond the saturation limit of the on-board CCDs. For this reason, it was necessary to produce a specific pixel mask that allows measuring the flux with fewer pixels (Metcalfe et al. 2012). The raw data were treated using the procedure described in García et al. (2011). It corrects for the instrumental perturbations (outliers, jumps, drifts, etc.) and merges time series for observations spanning multiple quarters.

The resulting precision on the measured flux, and subsequently, on the derived oscillation frequencies, are the best obtained by Kepler, and to this day, for any Sun-like star in addition to the Sun itself. Overall, 54 and 56 modes have been detected for16 Cyg A and B (Davies et al. 2015), respectively, with the precision ranging from 0.03 μHz to 2.74 μHz. Interestingly, Davies et al. (2015) reported projected rotational-splitting measurements of 411 ± 13 nHz for 16 Cyg A and 274 ± 17 nHz for 16 Cyg B, that is, measurements that include a sin i factor,. However, they did not detect conclusive signatures of differential rotation. Our goal is to remodel these acoustic power spectra with a model that differs from those used in Metcalfe et al. (2012) and Davies et al. (2015) in that it takes differential rotation into account. We use the same time series as in Davies et al. (2015).

2.2. Spectrum model

The fitting of the power spectra obtained from these times series is based on principles that are commonly adopted in asteroseismic studies. The first assumption is that the power spectrum in each frequency bin we consider in the Fourier space is independent of (and therefore not correlated to) its neighbours; that is, we neglect leakage coming from the convolution of the Fourier transforms of the signal and the window function. This allows us to consider the probability density of the power in each frequency bin separately. Making the further assumption that the noise on the measurements is Gaussian, the power spectrum P_i = P(ω_i) at any frequency ω_i is exponentially distributed, with a probability density

$$\begin{array}{c}\hfill f(P_k)=\frac{1}{\mathcal{P}({\omega }_k)}exp(-\frac{P_k}{\mathcal{P}({\omega }_k)}),\end{array}$$(1)

with 𝒫(ω_k) the average value of the power spectrum at ω_k, where k is an index of the frequency bins. A common model for the average power spectrum is a sum of Lorentzian functions centred at the eigenfrequencies of the pulsation modes. This is suitable for regularly damped and re-excited modes such as those observed in the Sun (Anderson et al. 1990). The central frequencies, the widths, and the heights of the Lorentzian are free parameters of the spectrum model. Usually, multiplets resulting from a rotationally induced lifting of degeneracy are fit jointly using a relation of the form ν_n, l, m = f(ν_n, l, m; (a_j)_{1 ≤ j ≤ J}). Here ν_n, l is the “central” frequency of the multiplet, of radial order n and angular degree l, m is the azimuthal order, and (a_j)_{1 ≤ j ≤ J} is a vector of coefficients allowing us to expand the splitting (Ritzwoller & Lavely 1991) as

$$\begin{array}{c}\hfill {\displaystyle \delta {\nu }_{n,l,m}={\nu }_{n,l,m}-{\nu }_{n,l}=\sum _{j=1}^J{a}_j(n,l){\zeta }_j^{(l)}(m),}\end{array}$$(2)

with the functions ${\zeta }_j^{(l)}(m)$ forming an orthogonal basis such that $\sum _m{\zeta }_j^{(l)}(m){\zeta }_i^{(l)}(m)=0$ for i ≠ j.

In the expression of 𝒫(ω_k) we introduce the effect of differential rotation. It is typical in asteroseismology to retain only the first-order term in the above expansion. The a₁ coefficient can be interpreted as a weighted average of the rotation throughout the star (Appourchaux et al. 2014; Davies et al. 2015). In this work, however, we also consider the next term in Eq. (2), as suggested by Gizon & Solanki (2004). This leads to a frequency distribution described by

$$\begin{array}{c}\hfill {\nu }_{n,l,m}={\nu }_{n,l}+m{a}_1(n,l)+{\beta }_{n,l,m}(\nu )+{\zeta }_3^{(l)}(m)a_3(n,l).\end{array}$$(3)

In the following we only consider splittings for l = 2, therefore we use ${\zeta }_3^{(2)}(m)=(5m^3-17m)/3$. We also note that a frequency-dependent term, β_n, l, m(ν), has been added to the expression resulting from Eq. (2). It includes the perturbation to the frequencies stemming from the asphericities of the star. They may be caused by the centrifugal forces, perhaps also by a large-scale magnetic field (Gough & Thompson 1990), a tidal distortion, or a strong anisotropic stellar wind. We show in Fig. 1 a typical multiplet, chosen at n = 20, l = 2, as observed in the power spectrum of 16 Cyg A and modelled using Eq. (3). The effects of the higher-order term and departure from sphericity are smaller than the contribution of a₁ to the splitting. This is shown in the right panel of Fig. 1, where we represent the individual contributions to the non-degenerate frequencies. It is interesting to note that ma₁ and ${\zeta }_3^{(l)}(m)a_3$ are symmetric functions of the azimuthal order m, while β_n, l, m is antisymmetric in m.

Fig. 1. Left panel: spectrum of 16 Cyg A in the region of the (l, n)=(2, 20) multiplet. Black shows the observations and red is the theoretical model corresponding to the MAP estimates of the parameters θ (cf. Eq. (4)). The vertical red ticks mark the position of the non-degenerate frequencies with −l ≤ m ≤ l. Right panel: splitting diagram for the multiplet. The contributions of the terms in Eq. (3) are represented individually. The final red horizontal ticks correspond to those in the left panel.

The coefficients a₁ and a₃ now become parameters of the model to be fitted to the observed spectrum. The relative heights of the modes in a multiplet also depend on the inclination i (Gizon & Solanki 2003). We denote the parameters necessary to describe the spectrum as

$$\begin{array}{c}\hfill \mathit{\theta }=({\mathit{\theta }}_{\mathit{S}},a_1,a_3,i,{\mathit{\beta }}_{\mathit{S}},\mathcal{B}).\end{array}$$(4)

Here θ_S ∈ ℝ^3N is a vector grouping of ν_n, l, H_n, l, and Γ_n, l, which are the frequency, height, and width of the Lorentzian describing the expectation value of the line profile of mode (n, l) in the power spectrum. N oscillation frequencies have been observed. Other parameters are a₁, a₃, the stellar inclination i ∈ [0, 2π], and β_S, a vector grouping of the parameters we used in the functional form of β_n, l, m. We also used a model to describe the noise contribution to the power spectrum. Its parameters are collectively denoted as ℬ, and we refer to Benomar et al. (2009) for a discussion of the issues of background-noise fitting in seismic spectra.

These parameters were estimated in a Bayesian framework, that is, we estimated the posterior density function of the parameter vector θ, conditional on y, the data

$$\begin{array}{c}\hfill {\pi }_{\mathit{\theta }|\mathit{y}}(\mathit{\theta }|\mathit{y})\propto {\pi }_{\mathit{\theta }}(\mathit{\theta })L(\mathit{\theta }).\end{array}$$(5)

Here y = (P₁, …, P_K) is the observation vector. We always note the probability density of a variable x by π_x. Likewise, if the probability is conditional on y, for instance, the corresponding density is noted π_x|y.

The likelihood function L is the product of the individual exponential distributions for all the bins considered.

$$\begin{array}{c}\hfill {\displaystyle L(\theta )=\prod _{k=1}^K\frac{1}{\mathcal{P}({\omega }_k;\theta )}exp(-\frac{P_k}{\mathcal{P}({\omega }_k;\theta )}).}\end{array}$$(6)

We recall that this function is the probability density of the data when y is the variable and is called the likelihood when seen as a function of the parameters, θ. In that case, it is not a probability density of the parameters.

The prior π_θ adopted for the parameters can be decomposed as the product of priors on the individual parameters by assuming that they are independent of each other. The priors on the usual parameters (frequencies, mode heights, and line widths) are described in the Appendix of White et al. (2016). In addition, we chose uniform priors on a₁ and a₃. The coefficient a₁ was assumed to be positive. Furthermore, because we do not expect the stars to be fast rotators, an upper boundary was set to 5 μHz. The boundaries for the prior on a₃ were more difficult to set, and we used a test-and-trial procedure. We finally set −0.1 μHz ≤ a₃ ≤ 0.1 μHz. We chose the following functional form for β_n, l, m (see also Sect. 4):

$$\begin{array}{c}\hfill {\beta }_{n,l,m}({\nu }_{n,l,m})={\beta }_0{Q}_{l,m}{\nu }_{n,l}.\end{array}$$(7)

Here we have Q_l, m = [L² − 3m²]/[(2l − 1)(2l + 3)], with L² = l(l + 1) (Gough & Thompson 1990; Kjeldsen et al. 1998), and β_S reduces to the scalar β₀.

The corresponding posterior density was sampled using an adaptive Markov chain Monte Carlo algorithm based on the scheme described by Haario et al. (2001), with some modifications inspired by Atchadé (2006). The simulations were run for 1 000 000 iterations. We used a burn-in sequence of 100 000 samples. In Fig. 1 we represent part of the spectrum obtained for the values of the parameters corresponding to the global maximum of π_θ|y, that is, the maximum a posteriori (MAP).

The resulting joint marginal densities for a₁, a₃, i, and β₀ are displayed in Figs. 2 and 3. The parameters seem fairly uncorrelated in general. The only exception worth noting is the obvious trend that can be observed in the joint density of (a₁, i). It is well documented (Gizon et al. 2013; Nielsen et al. 2014; Benomar et al. 2015) that the inferred average rotation rate, which is the main component of the splitting δν_n, l, m, increases when the inclination decreases. This is related to the lower visibility of the modes with m ≠ 0 when i ≲ 30^°, so that only larger splittings can be distinguished at lower angles. Two effects further complicate the fit of multiplets. First, if the mode lifetimes are too short, the splitting and widths of the modes become comparable, producing the blending of the non-degenerate modes. We recall that the mode lifetime τ is related to the width of the Lorentzian at half-maximum by Γ = 2/τ (see e.g. Chaplin et al. 1999). Second, l = 0 modes may obstruct l = 2 modes, which are of smaller amplitudes, if the widths of the peaks are too large. These effects will be aggravated when the signal-to-noise ratio becomes low.

Fig. 2. Marginal densities for the spectrum parameters a1, a3, i, and β0 of 16 Cyg A. The central panels show the joint marginal densities of the paired parameters. In the side panels we plot the individual marginal densities.

In Fig. 2 all marginal densities are normal to a good approximation. The inclination of 16 Cyg A is 58.5 ° ±6.8°, which confirms the findings of Davies et al. (2015), who obtained a posterior estimate of $56{}^{°}{}_{-5{}^{°}}^{+6{}^{°}}$. As previously discussed, this is an important parameter. The visibility of the modes at such an inclination will be low for the central frequency and higher for m = ±1, as shown in Fig. 1. The estimated value for a₁ is 464 ± 43 nHz. It is in good agreement with the splitting value derived in Davies et al. (2015). After deprojecting, their result is $486.3_{-29}^{+40}$ nHz. The first main result of this study is that we obtain a probability of 86% for a₃ to be strictly positive. Specifying further, we can define a 68.3% credible interval a₃ = 11.15 ± 10.95 nHz (given the Gaussian shape of the marginal density, this can be interpreted as a 1σ interval). Finally, the asphericity parameter is β₀ = 0.1 ± 0.1, that is, it is non-zero at a 1σ level. The significance of this result is discussed in Sect. 4 in relation to the magnetic properties of the star.

The situation is different for 16 Cyg B, as shown in Fig. 3. The probability densities of two parameters, namely i and a₁, show a bimodal behaviour. The most likely explanation is a poor constraint on the inclination, which is also reported in Davies et al. (2015). According to the correlation between the inclination and the rotational splitting described above, this in turn impacts the a₁ coefficient. The resulting bimodality of many marginal joint PDFs makes interpreting the inversion results difficult. These are discussed in Sect. 3. We note here that there are some differences with the results found in Davies et al. (2015). We tried to fit the power spectrum of 16 Cyg B setting a₃ = 0 nHz. In that case, we found the same results as they did. Therefore, we interpret the differences between our study and theirs as due to the including a₃.

Fig. 3. Same as Fig. 3 for 16 Cyg B.

3. Inversion of the rotation profile

The ultimate goal of this study is to provide a map of the rotation rates of 16 Cyg A and B. The derivation of the probability densities for a₁ and a₃ allows us to use methods that initially were developed for helioseismology to do so. The entire procedure relies on expanding the splitting according to Eq. (2) and the rotation rate, expressed in rad s⁻¹, in the form

$$\begin{array}{c}\hfill \mathrm{\Omega }(r,\theta )={\displaystyle \sum _{s=0}^S{\mathrm{\Omega }}_s(r)W_s(\theta ),}\end{array}$$(8)

with r the stellar radius and θ the co-latitude, the system being considered azimuthally symmetric.

In general, the orthogonality of the functions ${\zeta }_j^{(l)}(m)$ ensures that only modes with s ≥ j contribute to a_j (Brown et al. 1989). However, a particular set of functions W_s exists such that there is a one-to-one relation between a_2j + 1 and Ω_j (Ritzwoller & Lavely 1991; Schou et al. 1994; Pijpers 1997). This is obtained from the relation between the splitting and the rotation-rate components. If the rotational velocity field can be treated as a small perturbation to the hydrostatic equilibrium (Lynden-Bell & Ostriker 1967), the antisymmetric part of the frequency splitting can be related to the rotation rate through the linear integral equation

$$\begin{array}{c}\hfill {\displaystyle \frac{{\nu }_{n,l,m}-{\nu }_{n,l,-m}}{2}={\int }_0^{\pi }{\int }_0^{R_{\star }}{K}_{n,l,m}(r,\theta )\frac{\mathrm{\Omega }(r,\theta )}{2\pi }r\mathrm{d}\theta \mathrm{d}r.}\end{array}$$(9)

Here K_n, l, m is the rotation kernel for the mode defined by n, l, and m, and expressed as (Hansen et al. 1977)

$$\begin{array}{cc}\hfill K_{n,l,m}(r,\theta )=& \frac{m}{I_{n,l}}[{\xi }_{n,l}(r)[{\xi }_{n,l}(r)-\frac{2}{L}{\eta }_{n,l}(r)]P_l^m{(x)}^2\hfill \\ \hfill & +{\left(\frac{{\eta }_{n,l}(r)}{L}\right)}^2[{\left(\frac{\mathrm{d}{P}_l^m}{\mathrm{d}x}\right)}^2(1-x^2)+2P_l^m\frac{\mathrm{d}{P}_l^m}{\mathrm{d}x}x]\hfill \\ \hfill & +\frac{m^2}{1-x^2}{P}_l^m{(x)}^2]\rho (r)rsin\theta .\hfill \end{array}$$(10)

We recall the classical notation used above: considering spherical coordinates defined by the basis (e_r, e_θ, e_ϕ), the total displacement of a fluid element from its equilibrium state is

$$\begin{array}{c}\hfill {\mathit{\delta }\mathit{r}}_{n,l,m}(\mathit{r})={\xi }_{n,l}(r)Y_l^m(\theta ,\varphi ){\mathit{e}}_{\mathit{r}}+{\eta }_{n,l}(r){\mathrm{\nabla }}_h{Y}_l^m(\theta ,\varphi ),\end{array}$$

with r the position vector and $Y_l^m(\theta ,\varphi )$ the spherical harmonics. We have used the normalised associated Legendre polynomials, $P_l^m$, and defined x = cos θ. We also introduced the mode inertia

$$\begin{array}{c}\hfill I_{n,l}={\displaystyle {\int }_0^R\rho (r)r^2[{\xi }_{n,l}{(r)}^2+{\eta }_{n,l}{(r)}^2]\mathrm{d}r.}\end{array}$$(11)

In the following, we assume that the radius-dependent coefficients Ω_s(r) are piecewise constant functions that can be written explicitly as

$$\begin{array}{c}\hfill \mathrm{\Omega }(r,\theta )=\{\begin{array}{cc}{\mathrm{\Omega }}_0\hfill & \text{if}r<R_{\text{CZ}},\hfill \\ {\mathrm{\Omega }}_0+{\mathrm{\Omega }}_1{W}_1(\theta )\hfill & \text{if}r\ge R_{\text{CZ}}.\hfill \end{array}\end{array}$$(12)

We have used W₁(θ) = 1.5(5 cos²θ − 1) for the first-order basis function (Ritzwoller & Lavely 1991; Schou et al. 1994). The kernels were computed from Eq. (10) using the output of the code for stellar structure and evolution ASTEC (Christensen-Dalsgaard 2008a) and the stellar pulsation code adipls (Christensen-Dalsgaard 2008b). The results, discussed below, are summarised in Table 1.

It may then be shown (Gizon & Solanki 2004) that

$$\begin{array}{cc}& 2\pi a_1={\mathrm{\Omega }}_0{\displaystyle {\int }_0^{\pi }{\int }_0^{R_{\star }}{K}_{2,2}(r,\theta )r\mathrm{d}\theta \mathrm{d}r,}\hfill \end{array}$$(13)

$$\begin{array}{cc}& ={\mathrm{\Omega }}_0{\mathcal{K}}_0,\hfill \end{array}$$(14)

$$\begin{array}{cc}& 2\pi a_3={\mathrm{\Omega }}_1{\displaystyle {\int }_0^{\pi }{\int }_{R_{\mathrm{CZ}}}^{R_{\star }}\frac{[K_{2,2}(r,\theta )-K_{2,1}(r,\theta )]}{5}{W}_1(\theta )r\mathrm{d}\theta \mathrm{d}r,}\hfill \end{array}$$(15)

$$\begin{array}{cc}& ={\mathrm{\Omega }}_1{\mathcal{K}}_1,\hfill \end{array}$$(16)

where R_⋆ and R_CZ are the radius of the star and of the bottom of its convective envelope, respectively. The functions K_l, m = ⟨K_n, l, m⟩_n are the rotational kernels K_n, l, m averaged, with equal weights, over the radial orders.

The derivation of Eq. (15) is straightforward using Eq. (9) and the expression for ${\zeta }_3^{(2)}(m)$. To obtain Eq. (13), we assumed that δν_n, 1, 1 ≃ δν_n, 2, 2/2. This is well justified for 16 Cyg A and B. The typical departure between the two splittings is usually below 2%, which is well below the final uncertainties obtained on Ω₀, of the order of 10% or more. This assumption allowed us to use a₁ as estimated for modes with l = 2. This could be an advantage since the splittings are easier to measure for higher degree.

The rotational kernels are thus key quantities to the forward problem, that is, computing the rotational splitting using a theoretical stellar model. They depend in particular on the density of the stratified equilibrium model, ρ(r), and on the radial and horizontal mode displacements ξ_n, l(r) and η_n, l(r). These functions can be obtained by solving the equations for stellar structure, evolution, and pulsation. In order to compute them, a prerequisite is to obtain a realistic stellar model. As described in Appendix A, the model was obtained by fitting observed stellar properties, including asteroseismic data. In this section, we consider the best stellar models we obtained from our simulations. They are defined as those that maximise the posterior density function for the stellar parameters and are given in Table A.3.

After we computed the stellar models for 16 Cyg A and B are computed, the rotation kernels were obtained using Eqs. (13) and (15). We note that these two relations are sufficient to invert the profile because we only have access to a₁ and a₃. We tested spectrum models that included a₅ (sensitive only to l = 3 modes) in the truncation of the sum of the right-hand side of Eq. (8), but did not detect any significant departure from zero for this coefficient.

3.1. 16 Cyg A

The inversion for the rotation profile can be carried out straightforwardly for 16 Cyg A. The posterior marginal densities for a₁ and a₃, π_a1|y and π_a3|y, which we plot in the side panels of Fig. 2, 𝒦₀ and 𝒦₁, are known quantities. We can therefore obtain the posterior densities for Ω₀ and Ω₁. We applied Eqs. (14) and (16) to the samples obtained from the MCMC simulations {$a_1^{(1)}$, …, $a_1^{(T)}$} and {$a_3^{(1)}$, …, $a_3^{(T)}$}, with T the number of realisations in our sample. These are distributed approximately as π_a1|y and π_a3|y. This scaling gives us two new samples $\{a_1^{(1)}/2\pi {\mathcal{K}}_0,\cdots ,a_1^{(T)}/2\pi {\mathcal{K}}_0\}\sim {\pi }_{{\mathrm{\Omega }}_0|\mathit{y}}$ and $\{a_3^{(1)}/2\pi {\mathcal{K}}_1,\cdots ,a_3^{(T)}/2\pi {\mathcal{K}}_1\}\sim {\pi }_{{\mathrm{\Omega }}_1|\mathit{y}}$ (with the symbol “∼” meaning “distributed as”), which we used to approximate the desired marginal posterior probability densities.

This simple scaling obviously preserves the general structure of the posterior probability densities of the splitting coefficients. The densities for Ω₀ and Ω₁ are thus approximately described by normal densities. We estimated, in the sense of the posterior mean (PM), Ω₀/2π = 471 ± 43 nHz and Ω₁/2π = −42.7 ± 41.5 nHz. Here, we used the posterior standard deviation as a 68.3% credible interval. As expected, the non-null 1σ-level detection of latitudinal differential rotation obtained from a₃ translates into the rotation-rate coefficient. As before, the sign of Ω₁ can be assigned a probability, and the probability for it to be negative is still 86%.

The profile corresponding to the PM of Ω₀ and Ω₁ is given in Fig. 4 (left panel). It ranges from 534 nHz at the equator to 215 nHz at the pole, that is, a ratio of 2.5, which is significantly higher than that observed for the Sun. The uncertainties on this rotation profile are also shown in Fig. 4 (right panel). Their behaviour is a good indication of the regions of the stellar surface we can probe efficiently using the current seismic data. The overall structure of the posterior standard deviation seen in Fig. 4 is due to the functional form retained for Ω. The variance of Ω(R_⋆, θ) for a given co-latitude θ results from the posterior variances of Ω₀ and Ω₁. In Eq. (12) the term depending on Ω₀ does not vary with θ, hence it implies a minimum uncertainty on the rotation rate at any latitude. In contrast, Ω₁ is modulated by W₁, and we have W₁(θ = 63.4^°)=0. At this co-latitude, the rotation rate is equal in the radiative and convective envelopes and the variance on the rotation rate is minimal since Ω₁ does not contribute.

Fig. 4. Left panel: rotation-rate profile of 16 Cyg A corresponding to the PM estimates of Ω0 and Ω1. The dashed line marks the bottom of the convective envelope in the assumed stellar model. Right panel: distribution of the surface rotation rate as a function of the stellar co-latitude. The red shade represents probability densities at each latitude. The black line marks the surface rotation model for the PM values of Ω0 and Ω1 and corresponds to the map in the right panel at r = R⋆. The black dashed lines mark the mean surface rotation model plus (right curve) or minus (left curve) the corresponding posterior standard deviation estimated at each latitude.

The immediate result is that the higher latitudes are poorly constrained. This is not surprising because even in the solar case, these regions are the most difficult to probe using seismology (see e.g. Thompson et al. 2003, and references therein). The uncertainty at co-latitude 0^° is approximately ±260 nHz, which represents an error of ∼120%. The constraint is so poor that some extreme models even show an inversion of their rotation rate between the pole and the equator. Such models, although likely nonphysical, are formally admissible when only Eq. (12) is considered. Thus, the results at high latitudes should not be overinterpreted. The main conclusion that ought to be drawen is that stronger observational constraints are needed to tighten the precision of the inversion near the pole. Such constraints could help to improve our results in two ways. First, by better constraining Ω₁, large discrepancies of the rotation rate at high latitudes might be prevented. Second, by providing constraints to more accurate models, for instance, expansions of the form of Eq. (8) at higher orders (typically including a₅). We note that using such alternative rotation rates requires kernels that are sensitive to lower co-latitudes (see e.g. Lund et al. 2014).

Our results are much more constrained near the equator and up to co-latitude ∼50^°. At the equator, the uncertainty on the mean rotation rate is ±70 nHz, amounting to 13%. It then decreases to 42 nHz at ∼65^°, that is, an 8% uncertainty. The relative statistical error reaches 50% at roughly 35^°.

3.2. 16 Cyg B

The case of 16 Cyg B is more subtle and demands greater care. As we have discussed in Sect. 2, we observe some correlations in the joint marginal densities of some parameters. At first sight, this does not concern a₃. In Fig. 3 this parameter indeed looks as if it were not correlated to a₁ or i, and its marginal density π_a3|y appears to be roughly Gaussian. A normal approximation may thus be seen as suitable for modelling this density.

An examination of Eq. (3) gives some hints as to what may cause the observed bimodality. It shows that the frequency spacings between the rotationally split m-components of the l = 2 modes are ν_n, 2, 2 − ν_{n, 2, −2} = 4(a₁ + a₃) and ν_n, 2, 1 − ν_{n, 2, −1} = 2(a₁ − a₃). This can lead to degeneracies between a₁ and a₃ that can only be lifted if the mode blending (typically defined by the ratio a₁/Γ) and the noise level are not substantial.

When we estimate the first moments of the distribution with our MCMC sample, we obtain a₃  =  13.89  ±  13.95 nHz. Likewise, when we invert a₃ using the model discussed above and Eq. (15) to obtain Ω₁, the corresponding sample gives Ω₁/2π  =   − 51.68  ±  51.94 nHz. These values clearly indicate a non-detection of latitudinal differential rotation at a 68.3% level, that is, zero is included in the 1σ credible interval when we retain the posterior mean estimator in the Gaussian approximation. This said, the probabilities for a₃ to be strictly positive, or conversely, for Ω₁ to be strictly negative, remain high, at ∼85%.

All things being equal, this conclusion remains valid only as long as the normal approximation is accurate enough to describe π_a3|y or π_Ω1|y. In the following, we consider an alternative way to model the joint density π_Ω0, Ω1|y for (Ω₀, Ω₁) derived from the joint density π_a1, a3|y. We use a semi-parametric model (see e.g. Bishop 1995) called Gaussian mixture model. It provides an analytic form that can approximate the density of a random vector x ∼ π and is defined as

$$\begin{array}{c}\hfill \pi (\mathit{x})={\displaystyle \sum _{k=1}^K{p}_k{\mathcal{N}}_k({\mathit{\mu }}_k,{\mathbf{\Sigma }}_k).}\end{array}$$(17)

Here, 𝒩_k is a multivariate normal density of mean μ_k and covariance matrix Σ_k. The coefficients p_k are such that ∑_kp_k = 1. All these quantities are parameters that need to be estimated so that they reproduce the observed density satisfactorily, in our case, as approximated by the MCMC sample. We adopted a maximum-likelihood framework to do so. A classical way to obtain estimates of μ_k, Σ_k and p_k is then the expectation-minimisation algorithm (EM, Dempster et al. 1977). It is also necessary to fix the number, K, of Gaussian components to be used. We proceeded through trial-and-error steps.

Based on these principles, we proceeded in two steps. First, setting x = (Ω₀, Ω₁) in Eq. (17), we modelled the joint density π_Ω0, Ω1|y. This can be done straightforwardly using the one-to-one mappings relating a₁ and Ω₀, on one hand, and a₃ and Ω₁, on the other. Inversion then amounts to independently scaling the components of each vector and preserves the structure of the density. We selected a three-component Gaussian mixture model, which we consider to be the best trade-off between reproducing the main features of the joint probability density and over-fitting. It also has the advantage of being a simple enough model, so that it remains easy to interpret. This approach is still very simple. There exist many subtleties to mixture model fitting, with a vast literature treating them (e.g. Frühwirth-Schnatter 2006). Our goal here was to show that an improvement of the statistical model could lead us, at the 68.3% credibility level considered in this study as the reference detection threshold, from a relatively marginal non-detection to a relatively marginal detection. It is clear that from there on, significant improvements require better data. Therefore, we did not pursue more advanced techniques for the mixture modelling.

The results of the mixture modelling are shown in Fig. 5. The left panel shows the MCMC sample in the (Ω₀, Ω₁) plane and the estimated mixture model. The components of the latter can be separated into two groups. Two of them, those with the highest mean values on the Ω₀ axis (components k = 1 and 3 in Table 1), account for the peak to the right of the distribution, with a maximum Ω₀/2π >  500 nHz. They form the bulk of the density, which can be seen from the fact that the sum of their weights is 0.67. Two components in the model were necessary to account for the slightly longer tail of the main peak at higher values of a₁. The last component, with a weight of 0.33, represents the mode to the left, peaking at Ω₀/2π <  400 nHz. The global mode of the distribution is located at (Ω₀/2π, Ω₁/2π) = (558 nHz, −56.6 nHz) and is marked by a red dot. The right panel of Fig. 5 shows the projection of π_Ω0, Ω1|y on the Ω₀ and Ω₁ axis. The mixture model accounts for the marginal densities. In particular, it reproduces the two maxima in π_Ω0|y well. In the case of the marginal density of Ω₁, we also represented the results discussed above and obtained them from a Gaussian approximation.

Fig. 5. Left panel: joint marginal posterior density for (Ω0, Ω1) in 16 Cyg B. The black dots are the MCMC sample. The continuous curves show some iso-probability levels of the corresponding three-component Gaussian mixture model. Right panel: projections on the Ω0 and Ω1 axis of the two-dimensional joint probability. The histograms show the MCMC sample and the red lines the mixture model. The dashed lines show contributions of the two main peaks in the joint density (where the main peak combines components 1 and 3, and the secondary peak is component 2; see Table 1). The dot-dashed line in the Ω1 projection shows the result using a normal distribution with the mean and variance of the MCMC sample.

We now consider all three components. A first result was obtained by comparing the marginal density for Ω₁ resulting from the projection of the mixture model to the Gaussian approximation. The former is clearly a better approximation. Importantly, when we compute a 68.3% credible interval relative to its mode, we obtain ${\mathrm{\Omega }}_1/2\pi =-54.{45}_{-50.06}^{+52.12}$ nHz. This clearly excludes Ω₁/2π = 0 nHz. We thus showed that a proper modelling of the density for Ω₁ allows us to obtain a more convincing detection of latitudinal differential rotation. After we established that there is a differential rotation signal in the frequencies of 16 Cyg B, we derived a map of the rotation rate of the star Ω(r, θ). Unlike what was done in Sect. 3.1, where we retained the posterior mean (PM) estimates of Ω₀ and Ω₁, we used the MAP estimates of these parameters here. The resulting map is shown in the left panel of Fig. 6. The surface rate varies from 604 nHz at the equator to 235 nHz at the pole, for a ratio 2.6, which is only slightly higher that the rate observed for 16 Cyg A.

Fig. 6. Left panel: rotation profile Ω(r, θ) of 16 Cyg B corresponding to the MAP estimates of Ω0 and Ω1. The dashed line marks the bottom of the convective envelope. Right panel: distribution of the surface rotation rate as a function of the stellar co-latitude. The red shade represent probability densities. They have been normalised with respect to the rotational rate at fixed latitude. The black line marks the mode of the corresponding density, obtained from a Gaussian mixture model. The long-dash line shows the associated 68.3% credible interval. The short-dash line shows the model corresponding to the MAP estimates of Ω0 and Ω1.

The right panel of Fig. 6 shows the uncertainties on the surface rotation rate. It is apparent in the band of co-latitude 40^°–75^° that the density is bimodal, reflecting the shape of π_Ω0, Ω1|y. Here, we also represent the two estimates of Ω(R_⋆, θ) that we obtained: one from the direct modelling of the density π_{Ω(R⋆, θ),|y}, the other from the MAP estimate of the parameters Ω₀ and Ω₁. The latter is an approximation to the former. It was convenient to distinguish the two because π_{Ω(R⋆, θ),|y} allowed us to compute a credible interval on the surface rotation rate, while the MAP estimates of Ω₀ and Ω₁ were used to derive the map in the left panel. The approximation is valid since the two solutions never differ by more than 25% of the total width of the credible interval (with a maximum close to the pole) and is, in general, around or below 10% at co-latitudes higher than 15^°. The advantage of the rotation profile based on the MAP estimates of the rotation parameters is that using the parameters given in Table 1, it can be cast in the close simple analytic form of Eq. (8) and compared to other studies (see Sect. 4). The credible interval displayed in Fig. 6 was not as straightforward to derive as the interval shown in Fig. 4, for which we were able to use a Gaussian approximation. In this case we modelled the density of the surface rotation rate at each latitude, π_{Ω(R⋆, θ),|y}, using a Gaussian mixture model. The equatorial rotation rate obtained using this estimate is ${565}_{-129}^{+150}$ nHz.

The modelling of π_Ω0, Ω1|y allowed us to push the analysis further. It is possible to invert rotation profiles corresponding to the two main peaks found in the density. The less-likely component gives an MAP estimate ${\mathrm{\Omega }}_1/2\pi =-43.6_{-62.5}^{+63.0}$ nHz. The corresponding a₃ coefficient is $11.7_{-16.8}^{+16.9}$ nHz, which implies a marginal non-detection at a 68.3% level if this is the solution (the probability for Ω₁ to be negative remains high, however).

The main peak corresponds to the highest values of a₁, with an MAP credible interval ${531}_{-61}^{+90}$ nHz. The corresponding estimate of the a₃ coefficient is $-15.2_{-12.2}^{+12.5}$ nHz. If this turns out to be the solution, then it would correspond to a detection better than the 68.3% level. We represent the corresponding solution in Fig. 7 to provide a sense of the results that could be achieved if the data were good enough to constrain the inclination better and, consequently, a₁ and a₃.

Fig. 7. Same as Fig. 6, but for the mode corresponding to the two components with the highest mean values in Table 1.

4. Discussion

4.1. Impact of forward modelling

Up to this point, our inversion for differential rotation has relied on stellar models obtained from the MAP estimates of the stellar parameters. However, as we show in Appendix A, we also estimated uncertainties on these parameters. A legitimate question is thus whether such uncertainties can significantly affect our measurements of latitudinal differential rotation. In particular, the inferred value of Ω₁ could, in theory, be sensitive to the errors on the age and the mixing-length parameter. This could occur through the dependence of the model provided in Eq. (12) on the location of the bottom of the convective zone. The depth of the convective-radiative transition is controlled by the mixing-length parameter, which determines the magnitude of the superadiabatic gradient in the uppermost part of the convection zone and hence the entropy and structure of the adiabatically stratified bulk of the convection zone. Moreover, on the main sequence, this depth is known to decrease with the stellar age. Therefore, we determined the accuracy of the forward modelling, from stellar parameters to Ω₁, in light of these errors.

In Appendix A we obtain approximations to the joint posterior density of the stellar parameters θ_⋆ ∼ π_θ⋆|X⋆, where X_⋆ are the observations given in Table A.1 and θ_⋆ the stellar parameters. After this, we can also obtain approximations to the densities of a function h = h(θ_⋆), that is, π_{h(θ⋆)|X⋆}(h(θ_⋆)|X_⋆). The first step is to formulate a model for the stellar parameters. This can be done using the Bayesian framework we described in Sect. 2.

Table A.3 and Figs. A.1 and A.2 give for 16 Cyg A and B the estimates of the stellar parameters in the sense of the PM and the MAP and the marginal two- and one-dimensional probability densities for the stellar parameters, respectively. All the one-dimensional marginal densities are close to Gaussian. Only in the case of X₀ is the density slightly truncated as a result of the prior we chose.

In Fig. 8 we show the distributions for Ω₁ obtained using the MAP estimates for the stellar parameters, which were used to derive the rotation profiles in Figs. 4 and 6. We also display the probability density obtained by taking into account the variability of the scaling coefficient 𝒦₁ in Eq. (15) that is induced by the uncertainty on the stellar parameters. In order to obtain this latter density, we assumed that the a₃ parameter measured from the acoustic spectrum and the scaling coefficient are statistically independent. This is justified because the effect of rotation on stellar oscillations is treated as a perturbation and we did not take into account the effect of rotation on the stellar structure. In particular, the transport of angular momentum was neglected. Rewriting Eq. (15) as Ω₁ = 2πa₃/𝒦₁, we can derive the probability density π_{Ω1|X⋆, y} for Ω₁|X_⋆, y using a standard relation from probability theory that gives the density of a ratio of random variables,

Fig. 8. Densities for the differential rotation parameter Ω1 for 16 Cyg A (left panel) and B (right panel). The black lines show the densities after marginalisation over the 𝒦1xs coefficients given by Eq. (19). The red dashed lines show the densities for our best-fit stellar model.

$$\begin{array}{cc}& {\pi }_{{\mathrm{\Omega }}_1|{\mathit{X}}_{\mathbf{\star }},\mathit{y}}({\mathrm{\Omega }}_1|{\mathit{X}}_{\mathbf{\star }},\mathit{y})={\displaystyle \int {\pi }_{a_3|{\mathit{X}}_{\mathbf{\star }},\mathit{y}}\left(\frac{{\mathrm{\Omega }}_1{\mathcal{K}}_1}{2\pi }|{\mathit{X}}_{\mathbf{\star }},\mathit{y}\right)}\hfill \\ \hfill & \times {\pi }_{{\mathcal{K}}_1|{\mathit{X}}_{\mathbf{\star }}}({\mathcal{K}}_1\mid {\mathit{X}}_{\mathbf{\star }})|{\mathcal{K}}_1|\mathrm{d}{\mathcal{K}}_1,\hfill \\ \hfill & \approx \sum _i{\pi }_{a_3|\mathit{y}}\left(\frac{{\mathrm{\Omega }}_1{\mathcal{K}}_{1,i}}{2\pi }|{\mathit{X}}_{\mathbf{\star }},\mathit{y}\right)|{\mathcal{K}}_{1,i}|.\hfill \end{array}$$(18)

The density π_𝒦K1 was obtained from the MCMC simulations described above. We note that for any value of the stellar parameters we can compute new oscillation kernels and the corresponding values for the integral 𝒦₁, thus obtaining an approximation to the density π_𝒦K1. The second line in Eq. (18) is the approximation of the preceding integral using this MCMC sample. The sum is taken over all realisations. In order to compute this term, we need to know the density for a₃|y. This was done in a similar fashion as for Ω₁. We modelled the joint density for (a₁, a₃) using a three-component mixture model, each of them being bivariate Gaussian densities 𝒩(ζ_k,Λ_k) associated with weights q_k, k = 1, …, 3. This model can be marginalised analytically over a₁

$$\begin{array}{c}\hfill {\pi }_{a_3|\mathit{y}}(a_3|\mathit{y})={\displaystyle \sum _{k=1}^3{q}_k\frac{exp[-\frac{1}{2}({\lambda }_{2,k}-\frac{{\lambda }_{12,k}^2}{{\lambda }_{1,k}}{(a_3-{\zeta }_{2,k})}^2)]}{\sqrt{2\pi {\lambda }_{1,k}}{|{\mathbf{\Lambda }}_k|}^{-1/2}},}\end{array}$$(19)

with ζ_2, k the second coefficient of ζ_k and λ_1, k, λ_2, k, and λ_12, k the coefficients of the co-variance matrix Λ_k.

In Fig. 8 we show that the resulting densities are extremely close to those obtained using the MAP estimates. At any rate, this does not change our conclusions about differential rotation for 16 Cyg A and B. We can safely assess that differential-rotation measurements depend only weakly on the exact value of the stellar parameters. This is expected to hold as long as rotation can be treated as a perturbation. It is also important that the quality of the seismic data on 16 Cyg A and B offers a very good precision on the age and mixing-length parameters; this gives us an indication of the range, in the parameter space, over which our results can be regarded as robust. It remains to be understood at which level of precision this breaks down.

Another difficulty with forward modelling is the lack of information on the terms of higher order in the expansion of Ω(r, θ). Potentially, these might counteract the effects of the leading non-constant term in cos²θ. This is difficult to assess, however, since we were unable to measure a₅. Furthermore, adopting a slightly different decomposition $\mathrm{\Omega }(r,\theta )={\displaystyle \sum {\mathrm{\Omega }}_s^{\star }}(r){\cos }^{2s}\theta$ (Ritzwoller & Lavely 1991; Schou et al. 1994), ${\mathrm{\Omega }}_1^{\star }$ and ${\mathrm{\Omega }}_2^{\star }$ might be constrained using a₃. Of course, we loose in the process the one-to-one relation between a_j and ${\mathrm{\Omega }}_s^{\star }$. The meagre information obtained from spectrum fitting makes it difficult to properly estimate these parameters. So far, the best argument in favour of the preservation of latitudinal differential rotation when higher-order terms are included is the extrapolation from the solar case. We know that the term in cos⁴θ in the above expression also decreases the rotation rate as the co-latitude decreases. A common expression for the solar surface rotation rate in the convective zone is Ω_⊙(R_⊙, θ)/2π = 454 − 55 cos²θ − 76 cos⁴θ (Gizon & Solanki 2004). At the pole, the last term of the right-hand side contributes to ∼60% of the equator-to-pole braking. The missing information on higher-order terms is thus responsible for the poor constraints at high latitudes. As discussed in Sect. 3.1, the form used to describe the rotation rate implies a minimum variance at θ = 63.4^°. This is no longer the case when we introduce a term in cos⁴θ. In the Sun, this latter dominates differential rotation at co-latitudes ≲32^°, and hence the variance of the surface rotation rate in corresponding proportions. We postulate that the missing information on higher-order terms is responsible for the poor constraining of the rotation rate at high latitudes, ≳40^°, that we see in Figs. 4, 6, and 7. A possible solution to this problem is to obtain longer time series, potentially involving several l = 3 modes, which would allow measuring the a₅ coefficients. In that case, there would be a one-to-one relation between the former and the Ω₂ coefficients. These weight the function W₂(θ) in which the cos⁴θ terms appear.

We also disregarded the effect of subsurface flows on the rotation rate. It is well established in the solar case that it increases rapidly immediately below the surface (Deubner et al. 1979). For regions located in the co-latitude range 60^° − 90^°, the angular velocity gradient remains approximately constant. The layer at r = 0.97R_⊙ rotates ∼3% faster than the surface (Corbard & Thompson 2002). At lower co-latitudes, the gradient decreases in magnitude as a function of θ and even becomes positive below ∼35^°. This may bias the results presented here. The modes we used to infer a₃ are not the most sensitive to these subsurface layers, and the values derived here for latitudinal differential rotation may be more representative of the bulk of the convective zone rather than the surface itself, or the regions immediately below. The subsurface shear layer thus remains to be properly taken into account, as discussed for instance in Lund et al. (2014), in which the change in the rotational rate is uniformly modelled using a latitude-independent gradient.

4.2. Other differential rotation measurements

To our knowledge, no previous detection of differential rotation in 16 Cyg A or B has been reported so far. It is noteworthy that these stars are included in the BCool snapshot program (Marsden et al. 2014). This survey aims at detecting the average longitudinal component of stellar magnetic fields (Semel 1989; Landi Degl’Innocenti et al. 1992). Stars with conclusive detection then undergo further modelling of their magnetic topology. As a by-product, this provides an estimate of latitudinal differential rotation. The model used to reproduce the observations is a parametrisation of the magnetic field (Hussain et al. 2001) whose output is then transformed to reproduce the Zeeman profile. The latter is deconvolved from the observed V Stokes profile (Donati et al. 1997; Asensio Ramos et al. 2016). The mapping from the magnetic field to the one-dimensional Zeeman profile involves a convolution of the theoretical spectrum by the rotation profile of the star, also described here by Eq. (12). The differential rotation parameters are obtained, alongside those describing the magnetic field, using least-square minimisation. Unfortunately, the variations in the V Stokes spectra of 16 Cyg A and B could not be attributed with confidence to the magnetic field, which so far implies that their magnetic activity is not strong enough to allow a proper characterisation of the rotational profile.

Latitudinal differential rotation is usually quantified using either ΔΩ = Ω − Ω(R_⋆, 0^°) or the so-called shear parameter ΔΩ/Ω, where we set Ω = Ω(R_⋆, 90^°), the equatorial rotation rate. In the case of 16 Cyg A and B, we measured ΔΩ/2π = 320 ± 269 nHz and ${440}_{-383}^{+363}$ nHz, respectively, and $\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }=0.{65}_{-0.50}^{+0.47}$ and $0.{76}_{-0.58}^{+0.55}$. These limits correspond, as usual, to 68.3% credible intervals. We note also that the probabilities of the shear parameter to be positive are 85% and 86% for 16 Cyg A and B.

We can gain some perspective by comparing these results to the other measurements of latitudinal differential rotation provided by Benomar et al. (2018) and to other measurements obtained using spectroscopy. The asteroseismic measurements of these quantities are of a similar order of magnitude as those found for 16 Cyg A and B. However, with the exception of KIC 10963065, all the estimated shear parameters are higher. The most extreme case is KIC 9025370, for which the shear parameter is ∼4 times higher than in 16 Cyg A. This may be the reflection of a trend of deceasing magnitude for differential rotation with age. However, several factors may be at work here. As discussed below, the effective temperature (see below) is also of importance. More precise, statistical statements on the full sample of star with asteroseismically-measured latitudinal differential rotation are beyond the scope of this paper. The correlation between differential rotation and other stellar parameters will be considered in future studies.

Some other measurements of latitudinal differential rotation have also been obtained using observational techniques other than asteroseismology. In the following we focus on results obtained using spectroscopy. Even though claims of latitudinal differential rotation detection have been made using photometry (Reinhold et al. 2013; Lanza et al. 2014), it was pointed out by Aigrain et al. (2015) that the methods considered might not be entirely reliable. Therefore we do not consider them here. The only notable exception is the planet-hosting star Kepler-17 (Valio et al. 2017), for which the particular orbital configuration of the planetary system allows a precise measurement of latitudinal differential rotation. However, we recall that the method employed in this study might not lead to many detections in the future. However, even though these investigations used Kepler data, these data were not analyzed using asteroseismic techniques.

In the case of spectroscopy, Doppler imaging and (the closely related) Zeeman Doppler imaging have been important providers of latitudinal differential rotation estimates. Vogt et al. (1987) initially stated that Doppler imaging requires fast rotators so that it is the dominant mechanism that produces spectral line broadening. The same remark applies to Zeeman Doppler imaging, which can be seen, crudely, as a transposition of Doppler imaging into V Stokes profiles (Semel 1989; Brown et al. 1991). Petit et al. (2002) have shown that the method could be applied to moderate rotators to obtain latitudinal differential rotation measurements. In both cases, the methods require stellar spots that modulate the observed spectrum.

Compared to the estimates given in Sect. 3, Zeeman Doppler imaging often leads to smaller errors on the parameters of the rotation profile. An explanation is that the spot configurations encountered on the observed stars often imply surface tracers distributed on a wide range of latitudes, including, potentially, near the pole. Such observations may therefore constrain the rotation-rate profile over the entire stellar surface, while asteroseismic inversions are controlled by the sensitivity of the observed modes to the regions of the stellar interior that form their resonant cavities. Spectropolarimetric detections are always at levels ≳1.5σ (normal densities are assumed for the uncertainties on the differential rotation parameters), and sometimes better than 40σ.

It has been suggested by Barnes et al. (2005) that a relation between the effective temperatures of stars and the magnitude of their latitudinal differential rotation exists, giving the power law $\mathrm{\Delta }\mathrm{\Omega }\propto T_{\mathrm{teff}}^{8.92\pm 0.31}$. A theoretical explanation of this trend has been advanced by Kitchatinov & Olemskoy (2012), although not all behaviour could always be accounted for (Küker et al. 2011). In Fig. 9 we show a plot similar to Fig. 2 in Barnes et al. (2005), in which we display the values of the latitudinal differential rotation for 16 Cyg A and B alongside those obtained by Benomar et al. (2018), and others from Zeeman Doppler imaging or spectrographic measurements. For most cases we used the values published in this study for the effective temperature, except for HD 197890, for which we used the value of Casagrande et al. (2011). When several measurements for ΔΩ existed, we used a weighted mean and computed the corresponding standard deviation assuming Gaussian errors. We note that by doing so, we consider that the measurements are realisations of a random process, and we might be disregarding a temporal dependence of these variations that would correlate the measured values. We also added some stars observed later using spectropolarimetry or photometry: Kepler-17 (Valio et al. 2017), 61 Cyg A (Boro Saikia et al. 2016), HN Peg (Boro Saikia et al. 2015), ξ Boo (Morgenthaler et al. 2012), HD 35296 and HD 29615 (Waite et al. 2017), EK Dra (Waite et al. 2017), HD 141943 (Marsden et al. 2011), HD10650 (Waite et al. 2011), and τ Boo (Donati et al. 2008; Fares et al. 2009). The inclusion of 16 Cyg A and B does not seem to contradict the law described above. The use of weighted means seems to modify this result more than our new measurements. Barnes et al. (2005) used all the individual measurements to fit the data.

Fig. 9. Measured latitudinal differential rotation as a function of the effective temperature for 16 Cyg A and B (in red) and stars observed using spectropolarimetry, photometric transit, and asteroseismology. The Sun is represented by the blue ⊙ symbol. The grey triangles represent the upper bound provided by Ammler-von Eiff & Reiners (2012).

Alongside these differential-rotation measurements, we also display those given by Ammler-von Eiff & Reiners (2012). They were obtained for A and F stars, which are relatively fast rotators. However, these values have to be considered with caution because they are only upper limits.

4.3. Asymmetries in the power spectrum

The last point we wish to discuss is the deformation of the star by rotation and its effect on frequencies. Because stars are spheroids in rotation, it is evident that the centrifugal force may play a major role in determining the shape of the star.

Devising a method to detect the asphericity of a star might therefore give clues about the ongoing mechanisms in the stellar interior and at the surface. Interestingly, while the centrifugal force causes the star to become oblate, this is not necessarily the case for other forces. Typically, a toroidal magnetic field counteracts the effect of the centrifugal force such that the star may become prolate (Chandrasekhar 1953; Wentzel 1961).

The term β_n, l, m defined in Eq. (7) is only due to the centrifugal force, and so we can write

$$\begin{array}{cc}& {\beta }_{n,l,m}=\frac{4\pi }{3}\frac{Q_{\mathit{lm}}}{G{\rho }_{\star }}{\nu }_{n,l}{\mathrm{\Omega }}^2(r,\theta =0)\hfill \end{array}$$(20)

$$\begin{array}{cc}& \approx \frac{4\pi }{3G{\rho }_{\odot }}\frac{\mathrm{\Delta }{\nu }_{\odot }^2}{\mathrm{\Delta }{\nu }^2}{Q}_{\mathit{lm}}\nu a_1^2={\eta }_0{Q}_{\mathit{lm}}{\nu }_{n,l}{a}_1^2,\hfill \end{array}$$(21)

where G is the gravitational constant and ρ_⋆ is the average density of the star, approximated by the measure of the frequency spacing Δν of the pressure modes, ρ_⋆ ≈ ρ_⊙ (Δν/Δν_⊙)². Here we also assumed that the a₁ coefficient is representative of the the equatorial rotation. In the light of the result from the seismic inversion, this is indeed justified.

The actual pole-to-equator distortion of the star is then (Gizon et al. 2016)

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{R}_{\mathrm{c}}}{R}=\frac{3}{8\pi }{\eta }_0{a}_1^2,\end{array}$$(22)

with ΔR_c = R_eq − R_pole. Here R_eq is the equatorial radius and R_pole the polar radius of the star.

To measure an asphericity that is not only due to the centrifugal force, the general functional form β_n, l, m = β₀ Q_lm ν_n, l could be used. This relates to an effective asphericity coefficient ΔR/R by

$$\begin{array}{c}\hfill {\frac{\mathrm{\Delta }R}{R}|}_{\mathrm{eff}}=\frac{3}{8\pi }{\beta }_0.\end{array}$$(23)

This choice allows us to straightforwardly compare the case of a pure centrifugal force with cases with additional distorting forces.

In Fig. 10 we compare the effective ashericities and those computed from a₁ for 16 Cyg A and 16 Cyg B. They were determined using Eqs. (22) and (23). In the case of 16 Cyg B, although we cannot reject the possibility that the star is oblate (24.4% probability), it is striking that the measured asphericity is only marginally consistent with the case of a pure centrifugal force. The discrepancy is even more significant when considering 16 Cyg B, as the probability of being oblate is only 3.9%. One possible interpretation is that 16 Cyg A and B have a consequent (measurable) toroidal (equatorial) field that dampens equatorial waves such that they travel in a prolate cavity.

Fig. 10. Probability densities for the asphericities of 16 Cyg A (upper panel) and 16 Cyg B (lower panel). The red line shows the density resulting from a purely centrifugal force computed from a1. The black line shows the measured effective asphericity ΔR/R|eff. The probability that the star is oblate (ΔR/R >  0) is given. The shaded areas mark the 68.3% credible intervals of the corresponding distributions.

5. Conclusions

We have reported the detection of differential rotation for the two good solar analogues 16 Cyg A and B. We followed Benomar et al. (2018), where latitudinal differential rotation was detected in a sample of stars that either rotate faster than the Sun and/or with much higher differential rotation. In this case, the inferred values for the stellar rotation rate and for differential rotation are consistent with the solar regime.

We have described a way to model the latitudinal rotational splitting of 16 Cyg A and B by taking into account the impact of differential rotation. Using a Bayesian setting, we were able to state that the a₃ coefficients in these stars have a probability ≳85% of being strictly positive. Using an expansion basis for which a one-to-one relation exists between the coefficients for the rotational splitting and rotation rate, we translated this result into probabilistic statements on Ω₁, the coefficient quantifying the amount of latitudinal differential rotation in the convective zone. Importantly, it has the same ≳85% probability of being strictly negative for both stars. This indicates that it is very likely that the azimuthal component of the flow in the convective zone undergoes an equator-to-pole braking. We also provided summary statistics for Ω₁ and associated 68.3% Bayesian credible intervals that exclude zero for both stars. These results depend only very weakly on the errors on the stellar parameters, which reinforces the robustness of our conclusions.

These results agree very well with other estimates of latitudinal differential rotation obtained either by spectropolarimetry, spectroscopy, or photometry. In particular, they seem to follow the ΔΩ − T_eff relationship suggested by Barnes et al. (2005). They are of particular significance in that they represent the first conclusive detection for solar analogues, however. So far, mostly young active stars, often still in the post-main sequence stage, have yielded convincing measurements. Together with the results of Benomar et al. (2018), this work has opened the door for the study of differential rotation in main-sequence stars using asteroseismology. More precisely, it demonstrates the feasibility of such a detection for solar analogues. Studying such objects is important since their physical states are similar to those of the Sun. Therefore, theoretical models developed for this latter are likely to still apply to these stars. A fascinating perspective would thus be to gather similar data for other solar analogues and/or solar twins, using instruments such as SONG, TESS, and PLATO, in order to be able to constrain existing theoretical models for differential rotation and even dynamo.

Acknowledgments

We thank the referee for their careful reading of the manuscript. The authors thank C. Hanson for interesting discussions. This material is based upon work supported by the NYU Abu Dhabi Institute under grant G1502. Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant DNRF106). The research was supported by the ASTERISK project (ASTERoseismic Investigations with SONG and Kepler) funded by the European Research Council (Grant agreement no.: 267864). In memory of Michael J. Thompson.

References

Appendix A: Modelling 16 Cyg A and B

A common exercise, albeit challenging, in stellar physics consists of estimating stellar parameters such as the mass, M_⋆, the age, t_⋆, the initial composition given by the initial hydrogen-mass fraction, X₀, and initial metallicity, Z₀, and the mixing-length parameter, α. In the following, we group them in a single vector θ_⋆ = (M_⋆, t_⋆, X₀, Z₀, α). The parameters are real numbers. We not only wish to estimate the value of θ_⋆, but also the uncertainty in this value due to the errors on the observational constraints.

The observational data can come from many different sources, spectroscopy, photometry, interferometry, or astrometry. In the case of 16 Cyg A and B, we used an effective temperature and surface metallicity, T_eff and [Fe/H], derived from high-precision spectroscopy measurements (Ramírez et al. 2009). The radius was obtained using interferometry (White et al. 2013). The luminosity was derived from the astrometric Hipparcos parallax (van Leeuwen 2007). The seismic data were processed from the same Kepler time series as we used in this study, published by Davies et al. (2015). We did not use the individual frequencies directly to constrain our model because this demands the introduction of heuristic surface correction to our theoretical model. Rather, we used the frequency ratios

$$\begin{array}{cc}& r_{01}(n)=\frac{{\nu }_{n-1,0}-4{\nu }_{n-1,1}+6{\nu }_{n,0}-4{\nu }_{n,1}+{\nu }_{n+1,0}}{8({\nu }_{n,1}-{\nu }_{n-1,1})},\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& r_{02}(n)=\frac{{\nu }_{n,0}-{\nu }_{n-1,2}}{{\nu }_{n,1}-{\nu }_{n-1,1}},\hfill \end{array}$$(A.2)

defined by Roxburgh & Vorontsov (2003). These are expected to be far less sensitive to the surface and thus stand out as adequate quantities for model fitting (Bazot 2013; Silva Aguirre et al. 2013). The non-seismic observational constraints are listed in Table A.1. In the following, the observations are grouped in a vector X_⋆ = (T_eff, ([Fe/H], R, L, r₀₁, r₀₂)), with r₀₁ = (r₀₁(n_01, 1),…,r₀₁(n_01, N)) and r₀₂ = (r₀₂(n_02, 1),…,r₀₂(n_02, M)) (the indices n_01, i and n_02, i represent the mode orders for which the corresponding ratio can be evaluated from the observed oscillation frequencies).

Many difficulties are present when fitting stellar data with theoretical models. In brief, the main challenges stem from the facts that the theoretically evaluated observables depend non-linearly on θ_⋆ and that the computational cost of stellar models is relatively high. The former issue implies that sophisticated methods of computational statistics may be required to solve the estimation problem. The latter problem makes it difficult to use such methods.

We are interested in obtaining approximations to the joint posterior density of the stellar parameters, θ_⋆ ∼ π_θ⋆|X⋆. We obtained an expression for this density using Bayes’ formula Eq. (5) for θ_⋆ and X_⋆. In this context, the likelihood was obtained by assuming that the observations are the sum of a deterministic and stochastic component

$$\begin{array}{c}\hfill {\mathit{X}}_{\mathbf{\star }}=\mathcal{S}({\mathit{\theta }}_{\mathbf{\star }})+\mathit{ϵ},\end{array}$$(A.3)

with 𝒮(θ_⋆) a mapping from the space of parameters to the space of observations that represents the stellar evolution code, and ϵ the realisation of a random vector.

We assumed that the uncertainties on T_eff, [Fe/H], L and R are Gaussian with respective standard deviations σ_Teff, σ_[Fe/H], σ_L, and σ_R the observational uncertainties. The components of r₀₁ and r₀₂ are correlated, therefore we treated these vector as two separate multivariate Gaussian densities 𝒩N(μ₀₁, Σ₀₁) and 𝒩N(μ₀₂, Σ₀₂). The covariance matrices were estimated numerically. Independent samples were generated for each frequency used for the evaluation of the components of r₀₁ and r₀₂, and were used to obtain samples for both random vectors. Using these samples, evaluating Σ₀₁ and Σ₀₂ is straightforward. The resulting likelihood is therefore

$$\begin{array}{cc}& \pi ({\mathit{X}}_{\mathbf{\star }}|{\mathit{\theta }}_{\mathbf{\star }})\propto exp[-\frac{1}{2}(\frac{{(T-⟨T_{\mathrm{eff}}⟩)}^2}{{\sigma }_{T_{\mathrm{eff}}}^2}+\frac{{([\mathrm{Fe}/\mathrm{H}]-⟨[\mathrm{Fe}/\mathrm{H}]⟩)}^2}{{\sigma }_{[\mathrm{Fe}/\mathrm{H}]}^2}\hfill \\ \hfill & +\frac{{(L-⟨L⟩)}^2}{{\sigma }_L^2}+\frac{{(R-⟨R⟩)}^2}{{\sigma }_R^2}+{({\mathit{r}}_{\mathbf{01}}-{\mathit{\mu }}_{\mathbf{01}})}^T{\mathbf{\Sigma }}_{\mathbf{01}}({\mathit{r}}_{\mathbf{01}}-{\mathit{\mu }}_{\mathbf{01}})\hfill \\ \hfill & +{({\mathit{r}}_{\mathbf{02}}-{\mathit{\mu }}_{\mathbf{02}})}^T{\mathbf{\Sigma }}_{\mathbf{02}}({\mathit{r}}_{\mathbf{02}}-{\mathit{\mu }}_{\mathbf{02}}))].\hfill \end{array}$$(A.4)

Here the average quantities, denoted by ⟨.⟩, are the observed quantities. We note the difference of functional form between Eq. (6) and Eq. (A.4). This stems from the difference in the underlying statistical model.

The first term in Eq. (A.3) does not have an analytic closed form. In order to express it, we must solve the equation for stellar structure and pulsations numerically. This was achieved using the Aarhus Stellar Evolution Code (ASTEC) for the former and adipls for the latter. We assumed spherical symmetry and no magnetic field. The opacities were obtained from OPAL tables (Iglesias & Rogers 1996), with low-T opacities from Ferguson et al. (2005), and the equation of state was interpolated from OPAL tables (Rogers & Nayfonov 2002). Nuclear reaction rates were taken from the NACRE collaboration (Angulo et al. 1999) and supplemented by the values given in Imbriani et al. (2005) for the ¹⁴N(p,γ)¹⁵O reaction. Convection was treated using the prescription from Böhm-Vitense (1958) for the mixing-length theory, the mean-free path of the fluid elements being proportional to the pressure scale-height.

Fig. A.1. Marginal densities for the stellar parameters M, t⋆, X0, Z0, and α of 16 Cyg A. The central panels show the joint marginal densities of the paired parameters. Individual marginal densities are plotted in the side panels.

The stellar parameters were considered to be independent. Therefore, the prior can be written π(θ_⋆) = π(M_⋆)π(t_⋆)π(X₀)π(Z₀)π(α). Priors were chosen as uniform because we do not have previous measurements on any of them. All these quantities are positive, therefore the lower bounds of these prior densities should be non-negative. To set both upper and lower bounds, we used the estimates obtained by Metcalfe et al. (2015) as first initial guesses. These were obtained using almost the exact same data. There are good indications of the range in which we expect the significant probability mass to be found. We then refined the boundaries on our priors using successive trial-and-error stages. This was done in order to avoid to sharp cuts in the domain of definition of the posterior density. This could indeed lead to numerical issues when sampling from a posterior density using an MCMC algorithm. The only notable exception to this procedure concerns the initial hydrogen-mass fraction, which cannot be higher than 0.75, which is its value after the primordial nucleosynthesis (see Bazot et al. 2012, 2016, for a discussion). The priors used in our statistical model are given in Table A.2.

Fig. A.2. Marginal densities for the stellar parameters M, t⋆, X0, Z0, and α of 16 Cyg B. The central panels show the joint marginal densities of the paired parameters. Individual marginal densities are plotted in the side panels.

The posterior density was sampled using an MCMC algorithm. The details of the algorithm can be found in Bazot et al. (2018). It was run on ten independent chains. Each chain was heated with a temperature T >  1, so that we initially sampled a posterior density of the form π(θ_⋆)π^1/T(X_⋆|θ_⋆). This procedure, known as simulated annealing according to Kirkpatrick et al. (1983), allowed us to sample the space of parameters for densities with much weaker variations than the original target. Therefore, proposals of an MCMC algorithm will tend to be accepted more often after heating. This sometimes helps avoiding that a Markov chain becomes stuck far from the real solution that is sought for, in low-probability regions. It is therefore possible to run a preliminary sequence of MCMC runs with decreasing temperatures (with the last chain having T = 1), in order to identify the regions of high probability. The parameter T was assigned a decreasing sequence 2^n_T with n_T an integer such that 0 ≤ n_T ≤ 6. The number of iteration was 1000 for n_T = 6 and 200 for 1 ≤ n_T ≤ 5. For n_T = 0, the chains were run until acceptable convergence was obtained. The chains were initialised at n_T = 6 using an overdispersed density, crudely estimated from a short test run. At subsequent stages, all chains were initialised at the MAP value of the previous one. Convergence was tested using several diagnostics: the cumulative mean, variance, and the r ratio defined by Gelman & Rubin (1992).

The two-dimensional and one-dimensional posterior densities for the stellar parameters of 16 Cyg A and B are shown in Figs. A.1 and A.2. Corresponding estimates are also given in Table A.3. We also give the initial helium-mass fraction, Y₀, which is an often-used parameter in the literature. The MAP values computed from the five-dimensional joint density are given in the first line. They are given without uncertainties. We also give the MAP values obtained from each marginalised one-dimensional posterior, as shown in the side panels of Figs. A.1 and A.2. In this case, we produced an associated credible interval. The latter is defined as the smallest interval of probability mass 0.683 that encompasses the MAP. Finally, we also give the posterior mean and the posterior variance. All three estimates agree with each other. The estimates so obtained are in fair agreement with those of Metcalfe et al. (2015), even though a detailed comparison with this work is well beyond the scope of our study. The uncertainty on the parameters is remarkably low. In particular, the age is known with a precision of roughly 200 Myr.

All Tables

All Figures

Fig. 1. Left panel: spectrum of 16 Cyg A in the region of the (l, n)=(2, 20) multiplet. Black shows the observations and red is the theoretical model corresponding to the MAP estimates of the parameters θ (cf. Eq. (4)). The vertical red ticks mark the position of the non-degenerate frequencies with −l ≤ m ≤ l. Right panel: splitting diagram for the multiplet. The contributions of the terms in Eq. (3) are represented individually. The final red horizontal ticks correspond to those in the left panel.

In the text

	Fig. 2. Marginal densities for the spectrum parameters a1, a3, i, and β0 of 16 Cyg A. The central panels show the joint marginal densities of the paired parameters. In the side panels we plot the individual marginal densities.
In the text

	Fig. 3. Same as Fig. 3 for 16 Cyg B.
In the text

Fig. 4. Left panel: rotation-rate profile of 16 Cyg A corresponding to the PM estimates of Ω0 and Ω1. The dashed line marks the bottom of the convective envelope in the assumed stellar model. Right panel: distribution of the surface rotation rate as a function of the stellar co-latitude. The red shade represents probability densities at each latitude. The black line marks the surface rotation model for the PM values of Ω0 and Ω1 and corresponds to the map in the right panel at r = R⋆. The black dashed lines mark the mean surface rotation model plus (right curve) or minus (left curve) the corresponding posterior standard deviation estimated at each latitude.

In the text

Fig. 5. Left panel: joint marginal posterior density for (Ω0, Ω1) in 16 Cyg B. The black dots are the MCMC sample. The continuous curves show some iso-probability levels of the corresponding three-component Gaussian mixture model. Right panel: projections on the Ω0 and Ω1 axis of the two-dimensional joint probability. The histograms show the MCMC sample and the red lines the mixture model. The dashed lines show contributions of the two main peaks in the joint density (where the main peak combines components 1 and 3, and the secondary peak is component 2; see Table 1). The dot-dashed line in the Ω1 projection shows the result using a normal distribution with the mean and variance of the MCMC sample.

In the text

Fig. 6. Left panel: rotation profile Ω(r, θ) of 16 Cyg B corresponding to the MAP estimates of Ω0 and Ω1. The dashed line marks the bottom of the convective envelope. Right panel: distribution of the surface rotation rate as a function of the stellar co-latitude. The red shade represent probability densities. They have been normalised with respect to the rotational rate at fixed latitude. The black line marks the mode of the corresponding density, obtained from a Gaussian mixture model. The long-dash line shows the associated 68.3% credible interval. The short-dash line shows the model corresponding to the MAP estimates of Ω0 and Ω1.

In the text

	Fig. 7. Same as Fig. 6, but for the mode corresponding to the two components with the highest mean values in Table 1.
In the text

	Fig. 8. Densities for the differential rotation parameter Ω1 for 16 Cyg A (left panel) and B (right panel). The black lines show the densities after marginalisation over the 𝒦1xs coefficients given by Eq. (19). The red dashed lines show the densities for our best-fit stellar model.
In the text

	Fig. 9. Measured latitudinal differential rotation as a function of the effective temperature for 16 Cyg A and B (in red) and stars observed using spectropolarimetry, photometric transit, and asteroseismology. The Sun is represented by the blue ⊙ symbol. The grey triangles represent the upper bound provided by Ammler-von Eiff & Reiners (2012).
In the text

Fig. 10. Probability densities for the asphericities of 16 Cyg A (upper panel) and 16 Cyg B (lower panel). The red line shows the density resulting from a purely centrifugal force computed from a1. The black line shows the measured effective asphericity ΔR/R|eff. The probability that the star is oblate (ΔR/R >  0) is given. The shaded areas mark the 68.3% credible intervals of the corresponding distributions.

In the text

	Fig. A.1. Marginal densities for the stellar parameters M, t⋆, X0, Z0, and α of 16 Cyg A. The central panels show the joint marginal densities of the paired parameters. Individual marginal densities are plotted in the side panels.
In the text

	Fig. A.2. Marginal densities for the stellar parameters M, t⋆, X0, Z0, and α of 16 Cyg B. The central panels show the joint marginal densities of the paired parameters. Individual marginal densities are plotted in the side panels.
In the text