Latitudinal differential rotation in the solar analogues 16 Cyg A and B

TL;DR

This study measures latitudinal differential rotation in the solar analogues 16 Cyg A and B using asteroseismology, revealing near-solar rotation profiles and confirming the technique's effectiveness for stellar rotation analysis.

Contribution

First conclusive measurement of latitudinal differential rotation in solar analogues using asteroseismology and Bayesian modeling.

Findings

Rotation rates are higher at the equator than at the poles.

Differential rotation magnitudes are comparable to the Sun's.

Results agree with spectropolarimetry, spectroscopy, and photometry.

Abstract

Asteroseismology has undergone a profound transformation as a scientific field following the CoRoT and Kepler space missions. The latter is now yielding the first measurements of latitudinal differential rotation obtained directly from oscillation frequencies. Differential rotation is a fundamental mechanism of the stellar dynamo effect. Our goal is to measure the amount of differential rotation in the solar analogues 16 Cyg A and B, which are the components of a binary system. These stars are the brightest observed by Kepler and have therefore been extensively observed, with exquisite precision on their oscillation frequencies. We modelled the acoustic power spectrum of 16 Cyg A and B using a model that takes into account the contribution of differential rotation to the rotational frequency splitting. The estimation was carried out in a Bayesian setting. We then inverted these results to obtain the rotation profile of both stars under the assumption of a solar-like functional form. We observe that the magnitude of latitudinal differential rotation has a strong chance of being solar-like for both stars, their rotation rates being higher at the equator than at the pole. The measured latitudinal differential rotation, defined as the difference of rotation rate between the equator and the pole, is $320\pm269$ nHz and $440^{+363}_{-383}$ nHz for 16 Cyg A and B, respectively, confirming that the rotation rates of these stars are almost solar-like. Their equatorial rotation rates are $535\pm75$ nHz and $565_{-129}^{+150}$ nHz. Our results are in good agreement with measurements obtained from spectropolarimetry, spectroscopy, and photometry. We present the first conclusive measurement of latitudinal differential rotation for solar analogues. Their rotational profiles are very close to those of the Sun. These results depend weakly on the uncertainties of the stellar parameters.

Figures

Tables

Table 1. Table 1: Differential rotation parameters for 16 Cyg A and B. For the former, the posterior mean and standard deviation are given in nHz. For the latter, we give the parameters of the three-component Gaussian mixture model used to represent the joint posterior density, as inferred using the EM algorithm. The weights are p k subscript 𝑝 𝑘 p_{k} , k = 1 , 2 , 3 𝑘 1 2 3 k=1,2,3 . The vector ( Ω 0 / 2 ​ π subscript Ω 0 2 𝜋 \Omega_{0}/2\pi , Ω 1 / 2 ​ π subscript Ω 1 2 𝜋 \Omega_{1}/2\pi ) is the mean vector of each component (in nHz). The coefficients σ i , j subscript 𝜎 𝑖 𝑗 \sigma_{i,j} , i , j = 1 , 2 formulae-sequence 𝑖 𝑗 1 2 i,j=1,2 , of the covariance matrices are given in nHz 2 . We note that σ 1 , 2 = σ 2 , 1 subscript 𝜎 1 2 subscript 𝜎 2 1 \sigma_{1,2}=\sigma_{2,1} .

16 Cyg A
${\Omega }_0/2\pi$	$471\pm 43$
${\Omega }_1/2\pi$	$-42.7\pm 41.5$
16 Cyg B
	$p_1$ = 0.40	$p_2$ = 0.33	$p_3$ = 0.27
${\Omega }_0/2\pi$	518	385	618
${\Omega }_1/2\pi$	-47.1	-43.6	-68.1
${\sigma }_{1,1}$	2562	1557	4157
${\sigma }_{2,2}$	2029	3873	1849
${\sigma }_{1,2}$	-108	60.1	-647

Table 2. Table 2: Non-seismic observational properties for the two stars of the 16 Cyg system.

Star	$T_{\text{eff}}$ (K)	[Fe/H]	$L/L_{\odot }$	$R/R_{\odot }$
A	$5825\pm 50$	$0.10\pm 0.09$	$1.56\pm 0.05$	$1.22\pm 0.02$
B	$5750\pm 50$	$0.05\pm 0.06$	$1.27\pm 0.04$	$1.12\pm 0.05$

Table 3. Table 3: Lower and upper bounds used for the prior uniform densities for each stellar parameter.

Parameter	Lower bound	Upper bound
$M$ (M⊙)	0.7	1.25
$t_{\star }$ (Gyr)	1	13
$Z_0$	0.010	0.027
$X_0$	0.525	0.750
$\alpha$	1.	3.

Table 4. Table 4: Stellar parameters (and helium-mass fraction) inferred using the ASTEC stellar evolution code for 16 Cyg A and B. For each star, the first line gives the global MAP estimate. The second line is the MAP estimate for each marginalised density. The third line gives the posterior mean estimate.

Star	$M/M_{\odot }$	$t_{\star }$ (Gyr)	$X_0$	$Z_0$	$\alpha$	$Y_0$
	$1.06$	$6.74$	$0.694$	$0.0240$	$2.12$	$0.282$
16 Cyg A	${1.07}_{-0.02}^{+0.02}$	${6.70}_{-0.17}^{+0.23}$	${0.698}_{-0.013}^{+0.013}$	${0.0237}_{-0.0018}^{+0.0019}$	${2.12}_{-0.07}^{+0.09}$	${0.279}_{-0.013}^{+0.012}$
	1.07 (0.02)	6.73 (0.19)	0.698 (0.013)	0.0238 (0.0018)	2.13 (0.08)	0.278 (0.012)
	$1.05$	$6.73$	$0.711$	$0.0234$	$2.09$	$0.266$
16 Cyg B	${1.05}_{-0.02}^{+0.02}$	${6.63}_{-0.19}^{+0.20}$	${0.706}_{-0.015}^{+0.014}$	${0.0242}_{-0.0019}^{+0.0022}$	${2.12}_{-0.07}^{+0.08}$	${0.270}_{-0.013}^{+0.013}$
	1.05 (0.02)	6.64 (0.18)	0.706 (0.014)	0.0245 (0.0020)	2.12 (0.07)	0.270 (0.013)

Equations

$$f(P_k)=\frac{1}{\mathcal{P}({\omega }_k)}\exp \left(-\frac{P_k}{\mathcal{P}({\omega }_k)}\right),$$

$$\delta {\nu }_{n,l,m}={\nu }_{n,l,m}-{\nu }_{n,l}=\sum _{j=1}^J{a}_j(n,l){\zeta }_j^{(l)}(m),$$

$${\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).$$

$$\mathbf{\theta }=({\mathbf{\theta }}_{\mathbf{S}},a_1,a_3,i,{\mathbf{\beta }}_{\mathbf{S}},\mathcal{B}).$$

$${\pi }_{\mathbf{\theta }|\mathbf{y}}(\mathbf{\theta }|\mathbf{y})\propto {\pi }_{\mathbf{\theta }}(\mathbf{\theta })L(\mathbf{\theta }).$$

$$L(\theta )=\prod _{k=1}^K\frac{1}{\mathcal{P}({\omega }_k;\theta )}\exp \left(-\frac{P_k}{\mathcal{P}({\omega }_k;\theta )}\right).$$

$${\beta }_{n,l,m}({\nu }_{n,l,m})={\beta }_0{Q}_{l,m}{\nu }_{n,l}.$$

$$\Omega (r,\theta )=\sum _{s=0}^S{\Omega }_s(r)W_s(\theta ),$$

$$\frac{{\nu }_{n,l,m}-{\nu }_{n,l,-m}}{2}={\int }_0^{\pi }{\int }_0^{R_{\star }}{K}_{n,l,m}(r,\theta )\frac{\Omega (r,\theta )}{2\pi }rd\theta dr.$$

$$\begin{gathered} K_{n,l,m}(r,\theta )=\frac{m}{I_{n,l}}[{\xi }_{n,l}(r)\left[{\xi }_{n,l}(r)-\frac{2}{L}{\eta }_{n,l}(r)\right]P_l^m(x{)}^2 \\ +{\left(\frac{{\eta }_{n,l}(r)}{L}\right)}^2\left[{\left(\frac{d{P}_l^m}{dx}\right)}^2(1-x^2)+2P_l^m\frac{d{P}_l^m}{dx}x\right] \\ +\frac{m^2}{1-x^2}{P}_l^m(x{)}^2]\rho (r)r\sin \theta . \end{gathered}$$

$${\mathbf{\delta r}}_{n,l,m}(\mathbf{r})={\xi }_{n,l}(r)Y_l^m(\theta ,\varphi ){\mathbf{e}}_{\mathbf{r}}+{\eta }_{n,l}(r){\nabla }_h{Y}_l^m(\theta ,\varphi ),$$

$$I_{n,l}={\int }_0^R\rho (r)r^2[{\xi }_{n,l}(r{)}^2+{\eta }_{n,l}(r{)}^2]dr.$$

$$\Omega (r,\theta )=\{\begin{array}{ll}{\Omega }_0& \text{if r<RCZ,}\\ {\Omega }_0+{\Omega }_1{W}_1(\theta )& \text{if r≥RCZ.}\end{array}$$

$$2\pi a_1$$

$$={\Omega }_0{\mathcal{K}}_0,$$

$$2\pi a_3$$

$$={\Omega }_1{\mathcal{K}}_1,$$

$$\pi (\mathbf{x})=\sum _{k=1}^K{p}_k{\mathcal{N}}_k({\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k).$$

$${\pi }_{{\Omega }_1|{\mathbf{X}}_{\star },\mathbf{y}}({\Omega }_1|{\mathbf{X}}_{\star },\mathbf{y})=$$

$${\pi }_{{\mathcal{K}}_1|{\mathbf{X}}_{\star }}({\mathcal{K}}_1\mid {\mathbf{X}}_{\star })|{\mathcal{K}}_1|d{\mathcal{K}}_1,$$

$$\approx$$

$${\pi }_{a_3|\mathbf{y}}(a_3|\mathbf{y})=\sum _{k=1}^3{q}_k\frac{\exp \left[-\frac{1}{2}\left({\lambda }_{2,k}-\frac{{\lambda }_{12,k}^2}{{\lambda }_{1,k}}(a_3-{\zeta }_{2,k}{)}^2\right)\right]}{\sqrt{2\pi {\lambda }_{1,k}}|{\mathbf{\Lambda }}_k{|}^{-1/2}},$$

$${\beta }_{n,l,m}$$

$$\approx \frac{4\pi }{3G{\rho }_{\odot }}\frac{\Delta {\nu }_{\odot }^2}{\Delta {\nu }^2}{Q}_{lm}\nu a_1^2={\eta }_0{Q}_{lm}{\nu }_{n,l}{a}_1^2,$$

$$\frac{\Delta R_c}{R}=\frac{3}{8\pi }{\eta }_0{a}_1^2,$$

$${\frac{\Delta R}{R}|}_{\mathrm{eff}}=\frac{3}{8\pi }{\beta }_0.$$

$$r_{01}(n)$$

$$r_{02}(n)$$

$${\mathbf{X}}_{\star }=\mathcal{S}({\mathbf{\theta }}_{\star })+\mathbf{ϵ},$$

$$\begin{split}\displaystyle\pi(\boldsymbol{X_{\star}}|\boldsymbol{\theta_{\star}})\propto\exp\Bigg{[}-\frac{1}{2}\Bigg{(}\frac{(T-\langle T_{\mathrm{eff}}\rangle)^{2}}{\sigma^{2}_{T_{\mathrm{eff}}}}+\frac{(\mathrm{[Fe/H]}-\langle\mathrm{[Fe/H]}\rangle)^{2}}{\sigma^{2}_{\mathrm{[Fe/H]}}}+&\frac{(L-\langle L\rangle)^{2}}{\sigma^{2}_{L}}+\frac{(R-\langle R\rangle)^{2}}{\sigma^{2}_{R}}\\ &+(\boldsymbol{r_{01}}-\boldsymbol{\mu_{01}})^{T}\boldsymbol{\Sigma_{01}}(\boldsymbol{r_{01}}-\boldsymbol{\mu_{01}})+(\boldsymbol{r_{02}}-\boldsymbol{\mu_{02}})^{T}\boldsymbol{\Sigma_{02}}(\boldsymbol{r_{02}}-\boldsymbol{\mu_{02}})\Bigg{)}\Bigg{]}.\end{split}$$

Full text

*Aims. *Our goal is to measure the amount of differential rotation in the solar analogues 16 Cyg A and B, which are the components of a binary system. These stars are the brightest observed by Kepler and have therefore been extensively observed, with exquisite precision on their oscillation frequencies.

*Methods. *We modelled the acoustic power spectrum of 16 Cyg A and B using a model that takes into account the contribution of differential rotation to the rotational frequency splitting. The estimation was carried out in a Bayesian setting. We then inverted these results to obtain the rotation profile of both stars under the assumption of a solar-like functional form.

*Results. *We observe that the magnitude of latitudinal differential rotation has a strong chance of being solar-like for both stars, their rotation rates being higher at the equator than at the pole. The measured latitudinal differential rotation, defined as the difference of rotation rate between the equator and the pole, is $320\pm 269$ nHz and $440^{+363}_{-383}$ nHz for 16 Cyg A and B, respectively, confirming that the rotation rates of these stars are almost solar-like. Their equatorial rotation rates are $535\pm 75$ nHz and $565_{-129}^{+150}$ nHz. Our results are in good agreement with measurements obtained from spectropolarimetry, spectroscopy, and photometry.

*Conclusions. *We present the first conclusive measurement of latitudinal differential rotation for solar analogues. Their rotational profiles are very close to those of the Sun. These results depend weakly on the uncertainties of the stellar parameters.

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\simeq 200$ at low radial orders (Larson & Schou 2008). Higher values, up to $\gtrsim 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 Section 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 Section 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 Section 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 $\mu$Hz to 2.74 $\mu$Hz. Interestingly, Davies et al. (2015) reported projected rotational-splitting measurements of $411\pm 13$ nHz for 16 Cyg A and $274\pm 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(\omega_{i})$ at any frequency $\omega_{i}$ is exponentially distributed, with a probability density

[TABLE]

with $\mathcal{P}(\omega_{k})$ the average value of the power spectrum at $\omega_{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 $\nu_{n,l,m}=f(\nu_{n,l},m;(a_{j})_{1\leq j\leq J})$. Here $\nu_{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\leq j\leq J}$ is a vector of coefficients allowing us to expand the splitting (Ritzwoller & Lavely 1991) as

[TABLE]

with the functions $\zeta^{(l)}_{j}(m)$ forming an orthogonal basis such that $\sum\limits_{m}\zeta^{(l)}_{j}(m)\zeta^{(l)}_{i}(m)=0$ for $i\neq j$.

In the expression of $\mathcal{P}(\omega_{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_{1}$ 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

[TABLE]

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, $\beta_{n,l,m}(\nu)$, 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_{1}$ 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_{1}$ and $\zeta_{3}^{(l)}(m)a_{3}$ are symmetric functions of the azimuthal order $m$, while $\beta_{n,l,m}$ is antisymmetric in $m$.

The coefficients $a_{1}$ and $a_{3}$ 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

[TABLE]

Here $\boldsymbol{\theta_{S}}\in\mathbb{R}^{3N}$ is a vector grouping of $\nu_{n,l}$, $H_{n,l}$ , and $\Gamma_{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_{1}$, $a_{3}$, the stellar inclination $i\in[0,2\pi]$, and $\boldsymbol{\beta_{S}}$, a vector grouping of the parameters we used in the functional form of $\beta_{n,l,m}$. We also used a model to describe the noise contribution to the power spectrum. Its parameters are collectively denoted as $\boldsymbol{\mathcal{B}},$ 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 $\boldsymbol{\theta}$, conditional on $\boldsymbol{y}$, the data

[TABLE]

Here $\boldsymbol{y}=(P_{1},\dots,P_{K})$ is the observation vector. We always note the probability density of a variable $x$ by $\pi_{x}$. Likewise, if the probability is conditional on $\boldsymbol{y}$, for instance, the corresponding density is noted $\pi_{x|\boldsymbol{y}}$.

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

[TABLE]

We recall that this function is the probability density of the data when $\boldsymbol{y}$ is the variable and is called the likelihood when seen as a function of the parameters, $\boldsymbol{\theta}$. In that case, it is not a probability density of the parameters.

The prior $\pi_{\theta}$ 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_{1}$ and $a_{3}$. The coefficient $a_{1}$ was assumed to be positive. Furthermore, because we do not expect the stars to be fast rotators, an upper boundary was set to 5 $\mu$Hz. The boundaries for the prior on $a_{3}$ were more difficult to set, and we used a test-and-trial procedure. We finally set $-0.1$ $\mu$Hz $\leq a_{3}\leq$ 0.1 $\mu$Hz. We chose the following functional form for $\beta_{n,l,m}$ (see also Sect. 4):

[TABLE]

Here we have $Q_{l,m}=[L^{2}-3m^{2}]/[(2l-1)(2l+3)]$, with $L^{2}=l(l+1)$ (Gough & Thompson 1990; Kjeldsen et al. 1998), and $\boldsymbol{\beta_{S}}$ reduces to the scalar $\beta_{0}$.

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 $\pi_{\theta|\boldsymbol{y}}$, that is, the maximum a posteriori (MAP).

The resulting joint marginal densities for $a_{1}$, $a_{3}$, $i,$ and $\beta_{0}$ 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_{1},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 $\delta\nu_{n,l,m}$, increases when the inclination decreases. This is related to the lower visibility of the modes with $m\neq 0$ when $i\lesssim 30\degree$, 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 $\tau$ is related to the width of the Lorentzian at half-maximum by $\Gamma=2/\tau$ (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.

In Fig. 2 all marginal densities are normal to a good approximation. The inclination of 16 Cyg A is $58.5^{\circ}\pm 6.8^{\circ}$, which confirms the findings of Davies et al. (2015), who obtained a posterior estimate of $56\degree\,^{+6\degree}_{-5\degree}$. 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=\pm 1$, as shown in Fig.1. The estimated value for $a_{1}$ is $464\pm 43$ nHz. It is in good agreement with the splitting value derived in Davies et al. (2015). After deprojecting, their result is $486.3^{+40}_{-29}$ nHz. The first main result of this study is that we obtain a probability of 86% for $a_{3}$ to be strictly positive. Specifying further, we can define a 68.3% credible interval $a_{3}=11.15\pm 10.95$ nHz (given the Gaussian shape of the marginal density, this can be interpreted as a $1\sigma$ interval). Finally, the asphericity parameter is $\beta_{0}=0.1\pm 0.1$, that is, it is non-zero at a $1\sigma$ 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_{1}$ , 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_{1}$ 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_{3}=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_{3}$.

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_{1}$ and $a_{3}$ 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

[TABLE]

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

In general, the orthogonality of the functions $\zeta^{(l)}_{j}(m)$ ensures that only modes with $s\geq 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 $\Omega_{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

[TABLE]

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)

[TABLE]

We recall the classical notation used above: considering spherical coordinates defined by the basis $(\boldsymbol{e_{r}},\boldsymbol{e_{\theta}},\boldsymbol{e_{\phi}})$, the total displacement of a fluid element from its equilibrium state is

[TABLE]

with $\boldsymbol{r}$ the position vector and $Y_{l}^{m}(\theta,\phi)$ the spherical harmonics. We have used the normalised associated Legendre polynomials, $P_{l}^{m}$, and defined $x=\cos\theta$. We also introduced the mode inertia

[TABLE]

In the following, we assume that the radius-dependent coefficients $\Omega_{s}(r)$ are piecewise constant functions that can be written explicitly as

[TABLE]

We have used $W_{1}(\theta)=1.5(5\cos^{2}\theta-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

[TABLE]

where $R_{\star}$ and $R_{\mathrm{CZ}}$ are the radius of the star and of the bottom of its convective envelope, respectively. The functions $K_{l,m}=\langle K_{n,l,m}\rangle_{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 $\delta\nu_{n,1,1}\simeq\delta\nu_{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 $\Omega_{0}$, of the order of 10% or more. This assumption allowed us to use $a_{1}$ 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, $\rho(r)$, and on the radial and horizontal mode displacements $\xi_{n,l}(r)$ and $\eta_{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 4.

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_{1}$ and $a_{3}$. We tested spectrum models that included $a_{5}$ (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_{1}$ and $a_{3}$, $\pi_{a_{1}|\boldsymbol{y}}$ and $\pi_{a_{3}|\boldsymbol{y}}$, which we plot in the side panels of Fig. 2, $\mathcal{K}_{0}$ and $\mathcal{K}_{1}$, are known quantities. We can therefore obtain the posterior densities for $\Omega_{0}$ and $\Omega_{1}$. We applied Eqs. (14) and (16) to the samples obtained from the MCMC simulations $\{a_{1}^{(1)},\dots,a_{1}^{(T)}\}$ and $\{a_{3}^{(1)},\dots,a_{3}^{(T)}\}$, with $T$ the number of realisations in our sample. These are distributed approximately as $\pi_{a_{1}|\boldsymbol{y}}$ and $\pi_{a_{3}|\boldsymbol{y}}$. This scaling gives us two new samples $\{a_{1}^{(1)}/2\pi\mathcal{K}_{0},\dots,a_{1}^{(T)}/2\pi\mathcal{K}_{0}\}\sim\pi_{\Omega_{0}|\boldsymbol{y}}$ and $\{a_{3}^{(1)}/2\pi\mathcal{K}_{1},\dots,a_{3}^{(T)}/2\pi\mathcal{K}_{1}\}\sim\pi_{\Omega_{1}|\boldsymbol{y}}$ (with the symbol ‘$\sim$’ 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 $\Omega_{0}$ and $\Omega_{1}$ are thus approximately described by normal densities. We estimated, in the sense of the posterior mean (PM), $\Omega_{0}/2\pi=471\pm 43$ nHz and $\Omega_{1}/2\pi=-42.7\pm 41.5$ nHz. Here, we used the posterior standard deviation as a 68.3% credible interval. As expected, the non-null $1\sigma$-level detection of latitudinal differential rotation obtained from $a_{3}$ translates into the rotation-rate coefficient. As before, the sign of $\Omega_{1}$ can be assigned a probability, and the probability for it to be negative is still 86%.

The profile corresponding to the PM of $\Omega_{0}$ and $\Omega_{1}$ 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 $\Omega$. The variance of $\Omega(R_{\star},\theta)$ for a given co-latitude $\theta$ results from the posterior variances of $\Omega_{0}$ and $\Omega_{1}$. In Eq. (12) the term depending on $\Omega_{0}$ does not vary with $\theta$, hence it implies a minimum uncertainty on the rotation rate at any latitude. In contrast, $\Omega_{1}$ is modulated by $W_{1}$ , and we have $W_{1}(\theta=63.4\degree)=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 $\Omega_{1}$ does not contribute.

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\degree$ is approximately $\pm 260$ nHz, which represents an error of $\sim$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 $\Omega_{1}$, 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_{5}$). 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 $\sim$ $50\degree$. At the equator, the uncertainty on the mean rotation rate is $\pm 70$ nHz, amounting to 13%. It then decreases to 42 nHz at $\sim$65$\degree$, that is, an 8% uncertainty. The relative statistical error reaches 50% at roughly 35$\degree$.

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_{3}$. In Fig. 3 this parameter indeed looks as if it were not correlated to $a_{1}$ or $i$, and its marginal density $\pi_{a_{3}|\boldsymbol{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 $\nu_{n,2,2}-\nu_{n,2,-2}=4(a_{1}+a_{3})$ and $\nu_{n,2,1}-\nu_{n,2,-1}=2(a_{1}-a_{3})$. This can lead to degeneracies between $a_{1}$ and $a_{3}$ that can only be lifted if the mode blending (typically defined by the ratio $a_{1}/\Gamma$) and the noise level are not substantial.

When we estimate the first moments of the distribution with our MCMC sample, we obtain $a_{3}=13.89\pm 13.95$ nHz. Likewise, when we invert $a_{3}$ using the model discussed above and Eq. (15) to obtain $\Omega_{1}$, the corresponding sample gives $\Omega_{1}/2\pi=-51.68\pm 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\sigma$ credible interval when we retain the posterior mean estimator in the Gaussian approximation. This said, the probabilities for $a_{3}$ to be strictly positive, or conversely, for $\Omega_{1}$ to be strictly negative, remain high, at $\sim$85%.

All things being equal, this conclusion remains valid only as long as the normal approximation is accurate enough to describe $\pi_{a_{3}|\boldsymbol{y}}$ or $\pi_{\Omega_{1}|\boldsymbol{y}}$. In the following, we consider an alternative way to model the joint density $\pi_{\Omega_{0},\Omega_{1}|\boldsymbol{y}}$ for $(\Omega_{0},\Omega_{1})$ derived from the joint density $\pi_{a_{1},a_{3}|\boldsymbol{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 $\boldsymbol{x}\sim\pi$ and is defined as

[TABLE]

Here, $\mathcal{N}_{k}$ is a multivariate normal density of mean $\boldsymbol{\mu}_{k}$ and covariance matrix $\boldsymbol{\Sigma}_{k}$. The coefficients $p_{k}$ are such that $\sum_{k}p_{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 $\boldsymbol{\mu}_{k}$, $\boldsymbol{\Sigma}_{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 $\boldsymbol{x}=(\Omega_{0},\Omega_{1})$ in Eq. (17), we modelled the joint density $\pi_{\Omega_{0},\Omega_{1}|\boldsymbol{y}}$. This can be done straightforwardly using the one-to-one mappings relating $a_{1}$ and $\Omega_{0}$, on one hand, and $a_{3}$ and $\Omega_{1}$, 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 $(\Omega_{0},\Omega_{1})$ 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 $\Omega_{0}$ axis (components $k=1$ and 3 in Table 1), account for the peak to the right of the distribution, with a maximum $\Omega_{0}/2\pi>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_{1}$. The last component, with a weight of $0.33$, represents the mode to the left, peaking at $\Omega_{0}/2\pi<400$ nHz. The global mode of the distribution is located at $(\Omega_{0}/2\pi,\Omega_{1}/2\pi)=(558\leavevmode\nobreak\ \mathrm{nHz},-56.6\leavevmode\nobreak\ \mathrm{nHz})$ and is marked by a red dot. The right panel of Fig. 5 shows the projection of $\pi_{\Omega_{0},\Omega_{1}|\boldsymbol{y}}$ on the $\Omega_{0}$ and $\Omega_{1}$ axis. The mixture model accounts for the marginal densities. In particular, it reproduces the two maxima in $\pi_{\Omega_{0}|\boldsymbol{y}}$ well. In the case of the marginal density of $\Omega_{1}$, we also represented the results discussed above and obtained them from a Gaussian approximation.

We now consider all three components. A first result was obtained by comparing the marginal density for $\Omega_{1}$ 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 $\Omega_{1}/2\pi=-54.45_{-50.06}^{+52.12}$ nHz. This clearly excludes $\Omega_{1}/2\pi=0$ nHz. We thus showed that a proper modelling of the density for $\Omega_{1}$ 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 $\Omega(r,\theta)$. Unlike what was done in Sect. 3.1, where we retained the posterior mean (PM) estimates of $\Omega_{0}$ and $\Omega_{1}$, 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.

The right panel of Fig. 6 shows the uncertainties on the surface rotation rate. It is apparent in the band of co-latitude 40$\degree$ – 75$\degree$ that the density is bimodal, reflecting the shape of $\pi_{\Omega_{0},\Omega_{1}|\boldsymbol{y}}$. Here, we also represent the two estimates of $\Omega(R_{\star},\theta)$ that we obtained: one from the direct modelling of the density $\pi_{\Omega(R_{\star},\theta),|\boldsymbol{y}}$, the other from the MAP estimate of the parameters $\Omega_{0}$ and $\Omega_{1}$. The latter is an approximation to the former. It was convenient to distinguish the two because $\pi_{\Omega(R_{\star},\theta),|\boldsymbol{y}}$ allowed us to compute a credible interval on the surface rotation rate, while the MAP estimates of $\Omega_{0}$ and $\Omega_{1}$ 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$\degree$. 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, $\pi_{\Omega(R_{\star},\theta),|\boldsymbol{y}}$, using a Gaussian mixture model. The equatorial rotation rate obtained using this estimate is $565_{-129}^{+150}$ nHz.

The modelling of $\pi_{\Omega_{0},\Omega_{1}|\boldsymbol{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 $\Omega_{1}/2\pi=-43.6_{-62.5}^{+63.0}$ nHz. The corresponding $a_{3}$ 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 $\Omega_{1}$ to be negative remains high, however).

The main peak corresponds to the highest values of $a_{1}$, with an MAP credible interval $531_{-61}^{+90}$ nHz. The corresponding estimate of the $a_{3}$ 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_{1}$ and $a_{3}$.

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 $\Omega_{1}$ 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 $\Omega_{1}$, in light of these errors.

In Appendix A we obtain approximations to the joint posterior density of the stellar parameters $\boldsymbol{\theta_{\star}}\sim\pi_{\boldsymbol{\theta_{\star}}|\boldsymbol{X_{\star}}}$, where $\boldsymbol{X_{\star}}$ are the observations given in Table 2 and $\boldsymbol{\theta_{\star}}$ the stellar parameters. After this, we can also obtain approximations to the densities of a function $h=h(\boldsymbol{\theta_{\star}})$, that is, $\pi_{h(\boldsymbol{\theta_{\star}})|\boldsymbol{X_{\star}}}(h(\boldsymbol{\theta_{\star}})|\boldsymbol{X_{\star}})$. 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 4 and Figs. 11 and 12 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_{0}$ is the density slightly truncated as a result of the prior we chose.

In Fig. 8 we show the distributions for $\Omega_{1}$ 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 $\mathcal{K}_{1}$ 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_{3}$ 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 $\Omega_{1}=2\pi a_{3}/\mathcal{K}_{1}$, we can derive the probability density $\pi_{\Omega_{1}|\boldsymbol{X_{\star}},\boldsymbol{y}}$ for $\Omega_{1}|\boldsymbol{X_{\star}},\boldsymbol{y}$ using a standard relation from probability theory that gives the density of a ratio of random variables,

[TABLE]

The density $\pi_{\mathcal{K}_{1}}$ 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 $\mathcal{K}_{1}$, thus obtaining an approximation to the density $\pi_{\mathcal{K}_{1}}$. 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_{3}|\boldsymbol{y}$. This was done in a similar fashion as for $\Omega_{1}$. We modelled the joint density for $(a_{1},a_{3})$ using a three-component mixture model, each of them being bivariate Gaussian densities $\mathcal{N}(\boldsymbol{\zeta}_{k},\boldsymbol{\Lambda}_{k})$ associated with weights $q_{k}$, $k=1,\dots,3$. This model can be marginalised analytically over $a_{1}$

[TABLE]

with $\zeta_{2,k}$ the second coefficient of $\boldsymbol{\zeta}_{k}$ and $\lambda_{1,k}$, $\lambda_{2,k}$, and $\lambda_{12,k}$ the coefficients of the co-variance matrix $\boldsymbol{\Lambda}_{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 $\Omega(r,\theta)$. Potentially, these might counteract the effects of the leading non-constant term in $\cos^{2}\theta$. This is difficult to assess, however, since we were unable to measure $a_{5}$. Furthermore, adopting a slightly different decomposition $\Omega(r,\theta)=\sum\Omega_{s}^{\star}(r)\cos^{2s}\theta$ (Ritzwoller & Lavely 1991; Schou et al. 1994), $\Omega^{\star}_{1}$ and $\Omega^{\star}_{2}$ might be constrained using $a_{3}$. Of course, we loose in the process the one-to-one relation between $a_{j}$ and $\Omega^{\star}_{s}$. 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^{4}\theta$ 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 $\Omega_{\odot}(R_{\odot},\theta)/2\pi=454-55\cos^{2}\theta-76\cos^{4}\theta$ (Gizon & Solanki 2004). At the pole, the last term of the right-hand side contributes to $\sim$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 $\theta=63.4\degree$. This is no longer the case when we introduce a term in $\cos^{4}\theta$. In the Sun, this latter dominates differential rotation at co-latitudes $\lesssim 32\degree$, 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, $\gtrsim 40\degree$, 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_{5}$ coefficients. In that case, there would be a one-to-one relation between the former and the $\Omega_{2}$ coefficients. These weight the function $W_{2}(\theta)$ in which the $\cos^{4}\theta$ 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\degree-90\degree$, the angular velocity gradient remains approximately constant. The layer at $r=0.97R_{\odot}$ rotates $\sim$3% faster than the surface (Corbard & Thompson 2002). At lower co-latitudes, the gradient decreases in magnitude as a function of $\theta$ and even becomes positive below $\sim$35$\degree$. This may bias the results presented here. The modes we used to infer $a_{3}$ 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 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 $\Delta\Omega=\Omega-\Omega(R_{\star},0\degree)$ or the so-called shear parameter $\Delta\Omega/\Omega$, where we set $\Omega=\Omega(R_{\star},90\degree)$, the equatorial rotation rate. In the case of 16 Cyg A and B, we measured $\Delta\Omega/2\pi=320\pm 269$ nHz and $440^{+363}_{-383}$ nHz, respectively, and $\Delta\Omega/\Omega=0.65^{+0.47}_{-0.50}$ and $0.76^{+0.55}_{-0.58}$. 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 $\sim$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 $\gtrsim 1.5\sigma$ (normal densities are assumed for the uncertainties on the differential rotation parameters), and sometimes better than $40\sigma$.

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 $\Delta\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 $\Delta\Omega$ 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), $\xi$ 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 $\tau$ 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.

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 $\beta_{n,l,m}$ defined in Eq. (7) is only due to the centrifugal force, and so we can write

[TABLE]

where $G$ is the gravitational constant and $\rho_{\star}$ is the average density of the star, approximated by the measure of the frequency spacing $\Delta\nu$ of the pressure modes, $\rho_{\star}\approx\rho_{\odot}\,(\Delta\nu/\Delta\nu_{\odot})^{2}$. Here we also assumed that the $a_{1}$ 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)

[TABLE]

with $\Delta R_{c}=R_{\mathrm{eq}}-R_{\mathrm{pole}}$. Here $R_{\mathrm{eq}}$ is the equatorial radius and $R_{\mathrm{pole}}$ the polar radius of the star.

To measure an asphericity that is not only due to the centrifugal force, the general functional form $\beta_{n,l,m}=\beta_{0}\,Q_{lm}\,\nu_{n,l}$ could be used. This relates to an effective asphericity coefficient $\Delta R/R$ by

[TABLE]

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_{1}$ 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.

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_{3}$ coefficients in these stars have a probability $\gtrsim 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 $\Omega_{1}$, the coefficient quantifying the amount of latitudinal differential rotation in the convective zone. Importantly, it has the same $\gtrsim 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 $\Omega_{1}$ 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 $\Delta\Omega-T_{\text{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.

Acknowledgements.

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.

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_{\star}$, the age, $t_{\star}$, the initial composition given by the initial hydrogen-mass fraction, $X_{0}$, and initial metallicity, $Z_{0}$, and the mixing-length parameter, $\alpha$. In the following, we group them in a single vector $\boldsymbol{\theta_{\star}}=(M_{\star},t_{\star},X_{0},Z_{0},\alpha)$. The parameters are real numbers. We not only wish to estimate the value of $\boldsymbol{\theta_{\star}}$, 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_{\mathrm{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

[TABLE]

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 2. In the following, the observations are grouped in a vector $\boldsymbol{X_{\star}}=(T_{\mathrm{eff}},\mathrm{(}[Fe/H],R,L,\boldsymbol{r_{01}},\boldsymbol{r_{02}}))$, with $\boldsymbol{r_{01}}=(r_{01}(n_{01,1}),\dots,r_{01}(n_{01,N}))$ and $\boldsymbol{r_{02}}=(r_{02}(n_{02,1}),\dots,r_{02}(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 $\boldsymbol{\theta_{\star}}$ 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, $\boldsymbol{\theta_{\star}}\sim\pi_{\boldsymbol{\theta_{\star}}|\boldsymbol{X_{\star}}}$. We obtained an expression for this density using Bayes’ formula (5) for $\boldsymbol{\theta_{\star}}$ and $\boldsymbol{X_{\star}}$. In this context, the likelihood was obtained by assuming that the observations are the sum of a deterministic and stochastic component

[TABLE]

with $\mathcal{S}(\boldsymbol{\theta_{\star}})$ a mapping from the space of parameters to the space of observations that represents the stellar evolution code, and $\boldsymbol{\epsilon}$ the realisation of a random vector.

We assumed that the uncertainties on $T_{\mathrm{eff}}$, [Fe/H], $L$ and $R$ are Gaussian with respective standard deviations $\sigma_{T_{\mathrm{eff}}}$, $\sigma_{\mathrm{[Fe/H]}}$, $\sigma_{L}$ , and $\sigma_{R}$ the observational uncertainties. The components of $\boldsymbol{r_{01}}$ and $\boldsymbol{r_{02}}$ are correlated, therefore we treated these vector as two separate multivariate Gaussian densities $\mathcal{N}(\boldsymbol{\mu_{01}},\boldsymbol{\Sigma_{01}})$ and $\mathcal{N}(\boldsymbol{\mu_{02}},\boldsymbol{\Sigma_{02}})$. The covariance matrices were estimated numerically. Independent samples were generated for each frequency used for the evaluation of the components of $\boldsymbol{r_{01}}$ and $\boldsymbol{r_{02}}$, and were used to obtain samples for both random vectors. Using these samples, evaluating $\boldsymbol{\Sigma}_{01}$ and $\boldsymbol{\Sigma}_{02}$ is straightforward. The resulting likelihood is therefore

[TABLE]

Here the average quantities, denoted by $\langle.\rangle$, are the observed quantities. We note the difference of functional form between Eq. (6) and Eq. (27). This stems from the difference in the underlying statistical model.

The first term in Eq. (26) 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 14N(p,$\gamma$)15O 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.

The stellar parameters were considered to be independent. Therefore, the prior can be written $\pi(\boldsymbol{\theta_{\star}})=\pi(M_{\star})\pi(t_{\star})\pi(X_{0})\pi(Z_{0})\pi(\alpha)$. 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 3.

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 $\pi(\boldsymbol{\theta_{\star}})\pi^{1/T}(\boldsymbol{X_{\star}}|\boldsymbol{\theta_{\star}})$. 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\leq n_{T}\leq 6$. The number of iteration was 1000 for $n_{T}=6$ and 200 for $1\leq n_{T}\leq 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. 11 and  12. Corresponding estimates are also given in Table 4. We also give the initial helium-mass fraction, $Y_{0}$, 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. 11 and 12. 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.