Issue 
A&A
Volume 623, March 2019



Article Number  A125  
Number of page(s)  16  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201834594  
Published online  19 March 2019 
Latitudinal differential rotation in the solar analogues 16 Cygni A and B
^{1}
Division of Sciences, New York University Abu Dhabi, UAE
email: mb6215@nyu.edu
^{2}
Center for Space Science, NYUAD Institute, New York University Abu Dhabi, PO Box 129188, Abu Dhabi, UAE
^{3}
Stellar Astrophysics Centre, Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, 8000 Aarhus C, Denmark
^{4}
MaxPlanckInstitut für Sonnensystemforschung, Göttingen, Germany
^{5}
Institut für Astrophysik, GeorgAugustUniversität Göttingen, Göttingen, Germany
^{6}
Department of Astronomy & Astrophysics, Tata Institute of Fundamental Research, Mumbai 400005, India
^{7}
Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse UPS–OMP / CNRS, 14, avenue Édouard Belin, 31400 Toulouse, France
Received:
7
November
2018
Accepted:
7
January
2019
Context. 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.
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 solarlike functional form.
Results. We observe that the magnitude of latitudinal differential rotation has a strong chance of being solarlike 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 ± 269 nHz and 440_{−383}^{+363} nHz for 16 Cyg A and B, respectively, confirming that the rotation rates of these stars are almost solarlike. Their equatorial rotation rates are 535 ± 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.
Key words: asteroseismology / stars: oscillations / stars: rotation / stars: solartype / methods: statistical / methods: data analysis
© 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 closeby 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 sphericalharmonic 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, nonrotating star, the eigenfrequencies of the nonradial 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 quasiperiodic 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 meanfield approach of turbulent convection can explain differential rotation. The basic picture consists of describing the convective flow as a stationary component plus a timedependent turbulent term that after insertion in the NavierStokes 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 largescale 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 Sunlike 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 Sunlike 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 planethosting 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 Ftype stars (Reiners & Schmitt 2002a,b, 2003; Ammlervon Eiff & Reiners 2012). In a recent study, however, Benomar et al. (2018) demonstrated the possibility of measuring the magnitude of differential rotation in 13 Sunlike 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 solarlike 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 mixinglength 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 onboard 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 Sunlike 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 rotationalsplitting 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
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 reexcited 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
with the functions forming an orthogonal basis such that for i ≠ j.
In the expression of 𝒫(ω_{k}) we introduce the effect of differential rotation. It is typical in asteroseismology to retain only the firstorder 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
In the following we only consider splittings for l = 2, therefore we use . We also note that a frequencydependent 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 largescale 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 higherorder 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 nondegenerate frequencies. It is interesting to note that ma_{1} and 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 nondegenerate 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. 

Open with DEXTER 
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
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_{1}, a_{3}, 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 backgroundnoise 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
Here y = (P_{1}, …, 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 π_{xy}.
The likelihood function L is the product of the individual exponential distributions for all the bins considered.
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_{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 μHz. The boundaries for the prior on a_{3} were more difficult to set, and we used a testandtrial procedure. We finally set −0.1 μHz ≤ a_{3} ≤ 0.1 μHz. We chose the following functional form for β_{n, l, m} (see also Sect. 4):
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 β_{S} reduces to the scalar β_{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 burnin 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_{1}, a_{3}, i, and β_{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 δν_{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 nondegenerate modes. We recall that the mode lifetime τ is related to the width of the Lorentzian at halfmaximum 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 signaltonoise ratio becomes low.
Fig. 2. Marginal densities for the spectrum parameters a_{1}, a_{3}, 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. 

Open with DEXTER 
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 . 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_{1} is 464 ± 43 nHz. It is in good agreement with the splitting value derived in Davies et al. (2015). After deprojecting, their result is 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 ± 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} = 0.1 ± 0.1, that is, it is nonzero 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_{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}.
Fig. 3. Same as Fig. 3 for 16 Cyg B. 

Open with DEXTER 
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^{−1}, in the form
with r the stellar radius and θ the colatitude, the system being considered azimuthally symmetric.
In general, the orthogonality of the functions 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 onetoone 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 rotationrate components. If the rotational velocity field can be treated as a small perturbation to the hydrostatic equilibrium (LyndenBell & Ostriker 1967), the antisymmetric part of the frequency splitting can be related to the rotation rate through the linear integral equation
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)
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
with r the position vector and the spherical harmonics. We have used the normalised associated Legendre polynomials, , and defined x = cos θ. We also introduced the mode inertia
In the following, we assume that the radiusdependent coefficients Ω_{s}(r) are piecewise constant functions that can be written explicitly as
We have used W_{1}(θ) = 1.5(5 cos^{2}θ − 1) for the firstorder 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 (ChristensenDalsgaard 2008a) and the stellar pulsation code adipls (ChristensenDalsgaard 2008b). The results, discussed below, are summarised in Table 1.
Differential rotation parameters for 16 Cyg A and B.
It may then be shown (Gizon & Solanki 2004) that
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 . 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 Ω_{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, ρ(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_{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 righthand 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}, π_{a1y} and π_{a3y}, which we plot in the side panels of Fig. 2, 𝒦_{0} and 𝒦_{1}, are known quantities. We can therefore obtain the posterior densities for Ω_{0} and Ω_{1}. We applied Eqs. (14) and (16) to the samples obtained from the MCMC simulations {, …, } and {, …, }, with T the number of realisations in our sample. These are distributed approximately as π_{a1y} and π_{a3y}. This scaling gives us two new samples and (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 Ω_{0} and Ω_{1} are thus approximately described by normal densities. We estimated, in the sense of the posterior mean (PM), Ω_{0}/2π = 471 ± 43 nHz and Ω_{1}/2π = −42.7 ± 41.5 nHz. Here, we used the posterior standard deviation as a 68.3% credible interval. As expected, the nonnull 1σlevel detection of latitudinal differential rotation obtained from a_{3} translates into the rotationrate coefficient. As before, the sign of Ω_{1} can be assigned a probability, and the probability for it to be negative is still 86%.
The profile corresponding to the PM of Ω_{0} and Ω_{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 Ω. The variance of Ω(R_{⋆}, θ) for a given colatitude θ results from the posterior variances of Ω_{0} and Ω_{1}. In Eq. (12) the term depending on Ω_{0} does not vary with θ, hence it implies a minimum uncertainty on the rotation rate at any latitude. In contrast, Ω_{1} is modulated by W_{1}, and we have W_{1}(θ = 63.4^{°})=0. At this colatitude, the rotation rate is equal in the radiative and convective envelopes and the variance on the rotation rate is minimal since Ω_{1} does not contribute.
Fig. 4. Left panel: rotationrate 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 colatitude. 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. 

Open with DEXTER 
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 colatitude 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 Ω_{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 colatitudes (see e.g. Lund et al. 2014).
Our results are much more constrained near the equator and up to colatitude ∼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_{3}. In Fig. 3 this parameter indeed looks as if it were not correlated to a_{1} or i, and its marginal density π_{a3y} 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 mcomponents of the l = 2 modes are ν_{n, 2, 2} − ν_{n, 2, −2} = 4(a_{1} + a_{3}) and ν_{n, 2, 1} − ν_{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}/Γ) 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 ± 13.95 nHz. Likewise, when we invert a_{3} using the model discussed above and Eq. (15) to obtain Ω_{1}, the corresponding sample gives Ω_{1}/2π = − 51.68 ± 51.94 nHz. These values clearly indicate a nondetection 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_{3} to be strictly positive, or conversely, for Ω_{1} 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 π_{a3y} or π_{Ω1y}. In the following, we consider an alternative way to model the joint density π_{Ω0, Ω1y} for (Ω_{0}, Ω_{1}) derived from the joint density π_{a1, a3y}. We use a semiparametric 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
Here, 𝒩_{k} is a multivariate normal density of mean μ_{k} and covariance matrix Σ_{k}. The coefficients p_{k} are such that ∑_{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 maximumlikelihood framework to do so. A classical way to obtain estimates of μ_{k}, Σ_{k} and p_{k} is then the expectationminimisation algorithm (EM, Dempster et al. 1977). It is also necessary to fix the number, K, of Gaussian components to be used. We proceeded through trialanderror steps.
Based on these principles, we proceeded in two steps. First, setting x = (Ω_{0}, Ω_{1}) in Eq. (17), we modelled the joint density π_{Ω0, Ω1y}. This can be done straightforwardly using the onetoone mappings relating a_{1} and Ω_{0}, on one hand, and a_{3} and Ω_{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 threecomponent Gaussian mixture model, which we consider to be the best tradeoff between reproducing the main features of the joint probability density and overfitting. 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ühwirthSchnatter 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 nondetection 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 (Ω_{0}, Ω_{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 Ω_{0} axis (components k = 1 and 3 in Table 1), account for the peak to the right of the distribution, with a maximum Ω_{0}/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_{1}. The last component, with a weight of 0.33, represents the mode to the left, peaking at Ω_{0}/2π < 400 nHz. The global mode of the distribution is located at (Ω_{0}/2π, Ω_{1}/2π) = (558 nHz, −56.6 nHz) and is marked by a red dot. The right panel of Fig. 5 shows the projection of π_{Ω0, Ω1y} on the Ω_{0} and Ω_{1} axis. The mixture model accounts for the marginal densities. In particular, it reproduces the two maxima in π_{Ω0y} well. In the case of the marginal density of Ω_{1}, 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 isoprobability levels of the corresponding threecomponent Gaussian mixture model. Right panel: projections on the Ω_{0} and Ω_{1} axis of the twodimensional 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 dotdashed line in the Ω_{1} projection shows the result using a normal distribution with the mean and variance of the MCMC sample. 

Open with DEXTER 
We now consider all three components. A first result was obtained by comparing the marginal density for Ω_{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 nHz. This clearly excludes Ω_{1}/2π = 0 nHz. We thus showed that a proper modelling of the density for Ω_{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 Ω(r, θ). Unlike what was done in Sect. 3.1, where we retained the posterior mean (PM) estimates of Ω_{0} and Ω_{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.
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 colatitude. 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 longdash line shows the associated 68.3% credible interval. The shortdash line shows the model corresponding to the MAP estimates of Ω_{0} and Ω_{1}. 

Open with DEXTER 
The right panel of Fig. 6 shows the uncertainties on the surface rotation rate. It is apparent in the band of colatitude 40^{°}–75^{°} that the density is bimodal, reflecting the shape of π_{Ω0, Ω1y}. 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 Ω_{0} and Ω_{1}. 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 Ω_{0} and Ω_{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 colatitudes 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 nHz.
The modelling of π_{Ω0, Ω1y} 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 lesslikely component gives an MAP estimate nHz. The corresponding a_{3} coefficient is nHz, which implies a marginal nondetection at a 68.3% level if this is the solution (the probability for Ω_{1} to be negative remains high, however).
The main peak corresponds to the highest values of a_{1}, with an MAP credible interval nHz. The corresponding estimate of the a_{3} coefficient is 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}.
Fig. 7. Same as Fig. 6, but for the mode corresponding to the two components with the highest mean values in Table 1. 

Open with DEXTER 
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 Ω_{1} could, in theory, be sensitive to the errors on the age and the mixinglength 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 convectiveradiative transition is controlled by the mixinglength 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 Ω_{1}, 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 onedimensional probability densities for the stellar parameters, respectively. All the onedimensional 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 Ω_{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 𝒦_{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 Ω_{1} = 2πa_{3}/𝒦_{1}, we can derive the probability density π_{Ω1X⋆, y} for Ω_{1}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 𝒦_{1}xs coefficients given by Eq. (19). The red dashed lines show the densities for our bestfit stellar model. 

Open with DEXTER 
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 𝒦_{1}, 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_{3}y. This was done in a similar fashion as for Ω_{1}. We modelled the joint density for (a_{1}, a_{3}) using a threecomponent 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_{1}
with ζ_{2, k} the second coefficient of ζ_{k} and λ_{1, k}, λ_{2, k}, and λ_{12, k} the coefficients of the covariance 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 differentialrotation 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 mixinglength 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 nonconstant term in cos^{2}θ. This is difficult to assess, however, since we were unable to measure a_{5}. Furthermore, adopting a slightly different decomposition (Ritzwoller & Lavely 1991; Schou et al. 1994), and might be constrained using a_{3}. Of course, we loose in the process the onetoone relation between a_{j} and . 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 higherorder terms are included is the extrapolation from the solar case. We know that the term in cos^{4}θ in the above expression also decreases the rotation rate as the colatitude decreases. A common expression for the solar surface rotation rate in the convective zone is Ω_{⊙}(R_{⊙}, θ)/2π = 454 − 55 cos^{2}θ − 76 cos^{4}θ (Gizon & Solanki 2004). At the pole, the last term of the righthand side contributes to ∼60% of the equatortopole braking. The missing information on higherorder 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^{4}θ. In the Sun, this latter dominates differential rotation at colatitudes ≲32^{°}, and hence the variance of the surface rotation rate in corresponding proportions. We postulate that the missing information on higherorder 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_{5} coefficients. In that case, there would be a onetoone relation between the former and the Ω_{2} coefficients. These weight the function W_{2}(θ) in which the cos^{4}θ 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 colatitude 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 colatitudes, 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_{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 latitudeindependent 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 byproduct, 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 onedimensional 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 leastsquare 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 socalled 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 nHz, respectively, and and . 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 asteroseismicallymeasured 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 planethosting star Kepler17 (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 rotationrate 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 . 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: Kepler17 (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 Ammlervon Eiff & Reiners (2012). 

Open with DEXTER 
Alongside these differentialrotation measurements, we also display those given by Ammlervon 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
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, ρ_{⋆} ≈ ρ_{⊙} (Δν/Δν_{⊙})^{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 poletoequator distortion of the star is then (Gizon et al. 2016)
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} = β_{0} Q_{lm} ν_{n, l} could be used. This relates to an effective asphericity coefficient ΔR/R by
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.
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 a_{1}. 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. 

Open with DEXTER 
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 ≳85% of being strictly positive. Using an expansion basis for which a onetoone relation exists between the coefficients for the rotational splitting and rotation rate, we translated this result into probabilistic statements on Ω_{1}, 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 equatortopole braking. We also provided summary statistics for Ω_{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 ΔΩ − 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 postmain 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 mainsequence 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
 Aigrain, S., Llama, J., Ceillier, T., et al. 2015, MNRAS, 450, 3211 [NASA ADS] [CrossRef] [Google Scholar]
 Ammlervon Eiff, M., & Reiners, A. 2012, A&A, 542, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anderson, E. R., Duvall, Jr., T. L., & Jefferies, S. M. 1990, ApJ, 364, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Appourchaux, T., Chaplin, W. J., García, R. A., et al. 2012, A&A, 543, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Appourchaux, T., Antia, H. M., Benomar, O., et al. 2014, A&A, 566, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asensio Ramos, A., de la Cruz Rodríguez, J., Martínez González, M. J., & Pastor Yabar, A. 2016, A&A, 590, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Atchadé, Y. F. 2006, Method. Comput. Appl. Probab., 8, 235 [CrossRef] [Google Scholar]
 Baglin, A., Auvergne, M., Barge, P., et al. 2009, in Transiting Planets, eds. F. Pont, D. Sasselov, & M. J. Holman, IAU Symp., 253, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, J. R., Collier Cameron, A., Donati, J.F., et al. 2005, MNRAS, 357, L1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Bazot, M. 2013, in EAS Pub. Ser., eds. G. Alecian, Y. Lebreton, O. Richard, & G. Vauclair, 63, 105 [CrossRef] [Google Scholar]
 Bazot, M., Campante, T. L., Chaplin, W. J., et al. 2012, A&A, 544, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bazot, M., ChristensenDalsgaard, J., Gizon, L., & Benomar, O. 2016, MNRAS, 460, 1254 [NASA ADS] [CrossRef] [Google Scholar]
 Bazot, M., Creevey, O., ChristensenDalsgaard, J., & Meléndez, J. 2018, A&A, 619, A172 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benomar, O., Appourchaux, T., & Baudin, F. 2009, A&A, 506, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benomar, O., Takata, M., Shibahashi, H., Ceillier, T., & García, R. A. 2015, MNRAS, 452, 2654 [NASA ADS] [CrossRef] [Google Scholar]
 Benomar, O., Bazot, M., Nielsen, M. B., et al. 2018, Science, 361, 1231 [NASA ADS] [CrossRef] [Google Scholar]
 Bertello, L., Pevtsov, A. A., & Pietarila, A. 2012, ApJ, 761, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Bishop, C. M. 1995, Neural Networks for Pattern Recognition (New York, NY, USA: Oxford University Press, Inc.) [Google Scholar]
 BöhmVitense, E. 1958, Z. Astrophys., 46, 108 [Google Scholar]
 Boro Saikia, S., Jeffers, S. V., Petit, P., et al. 2015, A&A, 573, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boro Saikia, S., Jeffers, S. V., Morin, J., et al. 2016, A&A, 594, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, in American Astronomical Society Meeting Abstracts 215, BAAS, 42, 215 [Google Scholar]
 Bouchy, F., & Carrier, F. 2001, A&A, 374, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown, S. F., Donati, J.F., Rees, D. E., & Semel, M. 1991, A&A, 250, 463 [NASA ADS] [Google Scholar]
 Brown, T. M., ChristensenDalsgaard, J., Dziembowski, W. A., et al. 1989, ApJ, 343, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Strugarek, A., Varela, J., et al. 2017, ApJ, 836, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chandrasekhar, S. 1953, Proc. R. Soc. London, Ser. A, 217, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Chaplin, W. J., Elsworth, Y., Isaak, G. R., Miller, B. A., & New, R. 1999, MNRAS, 308, 424 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J. 2008a, Ap&SS, 316, 13 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J. 2008b, Ap&SS, 316, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Corbard, T., & Thompson, M. J. 2002, Sol. Phys., 205, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Couvidat, S. 2013, Sol. Phys., 282, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Davenport, J. R. A., Hebb, L., & Hawley, S. L. 2015, ApJ, 806, 212 [NASA ADS] [CrossRef] [Google Scholar]
 Davies, G. R., Chaplin, W. J., Farr, W. M., et al. 2015, MNRAS, 446, 2959 [NASA ADS] [CrossRef] [Google Scholar]
 Deheuvels, S., Ballot, J., Beck, P. G., et al. 2015, A&A, 580, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dempster, A. P., Laird, N. M., & Rubin, D. B. 1977, J. R. Stat. Soc., Ser. B, 39, 1 [Google Scholar]
 Deubner, F.L., Ulrich, R. K., & Rhodes, Jr., E. J. 1979, A&A, 72, 177 [NASA ADS] [Google Scholar]
 Donahue, R. A., Saar, S. H., & Baliunas, S. L. 1996, ApJ, 466, 384 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., & Collier Cameron, A. 1997, MNRAS, 291, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Donati, J.F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Donati, J.F., Moutou, C., Farès, R., et al. 2008, MNRAS, 385, 1179 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Fares, R., Donati, J.F., Moutou, C., et al. 2009, MNRAS, 398, 1383 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [Google Scholar]
 FrühwirthSchnatter, S. 2006, Finite Mixture and Markov Switching Models, Springer Series in Statistics (New York: Springer) [Google Scholar]
 García, R. A., Hekker, S., Stello, D., et al. 2011, MNRAS, 414, L6 [NASA ADS] [CrossRef] [Google Scholar]
 García, R. A., Ceillier, T., Salabert, D., et al. 2014, A&A, 572, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gelman, A., & Rubin, D. B. 1992, Statist. Sci., 7, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A., & Glatzmaier, G. A. 1981, ApJS, 45, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Solanki, S. K. 2003, ApJ, 589, 1009 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Solanki, S. K. 2004, Sol. Phys., 220, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., Ballot, J., Michel, E., et al. 2013, Proc. Nat. Acad. Sci., 110, 13267 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., Sekii, T., Takata, M., et al. 2016, Sci. Adv., 2, e1601777 [NASA ADS] [CrossRef] [Google Scholar]
 Gough, D. O., & Thompson, M. J. 1990, MNRAS, 242, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Gray, D. F. 1977, ApJ, 211, 198 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Smolarkiewicz, P. K., de Gouveia Dal Pino, E. M., Kosovichev, A. G., & Mansour, N. N. 2016, ApJ, 828, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Haario, H., Saksman, E., & Tamminen, J. 2001, Bernoulli, 7, 223 [CrossRef] [MathSciNet] [Google Scholar]
 Hansen, C. J., Cox, J. P., & van Horn, H. M. 1977, ApJ, 217, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Hill, F., Stark, P. B., Stebbins, R. T., et al. 1996, Science, 272, 1292 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Huang, S.S. 1961, ApJ, 133, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Hussain, G. A. J., Jardine, M., & Collier Cameron, A. 2001, MNRAS, 322, 681 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Imbriani, G., Costantini, H., Formicola, A., et al. 2005, Eur. Phys. J. A, 25, 455 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 King, J. R., Deliyannis, C. P., Hiltgen, D. D., et al. 1997, AJ, 113, 1871 [NASA ADS] [CrossRef] [Google Scholar]
 Kirkpatrick, S., Gelatt, C. D., & Vecchi, M. P. 1983, Science, 220, 671 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Kitchatinov, L. L., & Nepomnyashchikh, A. A. 2017, Astron. Lett., 43, 332 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Olemskoy, S. V. 2012, MNRAS, 423, 3344 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Rüdiger, G. 1995, A&A, 299, 446 [NASA ADS] [Google Scholar]
 Kjeldsen, H., Arentoft, T., Bedding, T., et al. 1998, in Structure and Dynamics of the Interior of the Sun and Sunlike Stars, ed. S. Korzennik, ESA Spec. Pub., 418, 385 [NASA ADS] [Google Scholar]
 Küker, M., Rüdiger, G., & Kitchatinov, L. L. 2011, A&A, 530, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Landi Degl’Innocenti, E. 1992, in Magnetic field measurements, eds. F. Sanchez, M. Collados, & M. Vazquez (Cambridge University Press), 71 [Google Scholar]
 Lanza, A. F., Das Chagas, M. L., & De Medeiros, J. R. 2014, A&A, 564, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Larson, T. P., & Schou, J. 2008, J. Phys. Conf. Ser., 118, 012083 [NASA ADS] [CrossRef] [Google Scholar]
 Lund, M. N., Lundkvist, M., Silva Aguirre, V., et al. 2014, A&A, 570, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LyndenBell, D., & Ostriker, J. P. 1967, MNRAS, 136, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Marsden, S. C., Jardine, M. M., Ramírez Vélez, J. C., et al. 2011, MNRAS, 413, 1939 [NASA ADS] [CrossRef] [Google Scholar]
 Marsden, S. C., Petit, P., Jeffers, S. V., et al. 2014, MNRAS, 444, 3517 [NASA ADS] [CrossRef] [Google Scholar]
 McQuillan, A., Mazeh, T., & Aigrain, S. 2014, ApJS, 211, 24 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Metcalfe, T. S., Chaplin, W. J., Appourchaux, T., et al. 2012, ApJ, 748, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Metcalfe, T. S., Creevey, O. L., & Davies, G. R. 2015, ApJ, 811, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Morgenthaler, A., Petit, P., Saar, S., et al. 2012, A&A, 540, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nielsen, M. B., Gizon, L., Schunker, H., & Karoff, C. 2013, A&A, 557, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nielsen, M. B., Gizon, L., Schunker, H., & Schou, J. 2014, A&A, 568, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petit, P., Donati, J.F., & Collier Cameron, A. 2002, MNRAS, 334, 374 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Pijpers, F. P. 1997, A&A, 326, 1235 [NASA ADS] [Google Scholar]
 Privitera, G., Meynet, G., Eggenberger, P., et al. 2016, A&A, 591, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ramírez, I., Meléndez, J., & Asplund, M. 2009, A&A, 508, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reinhold, T., & Reiners, A. 2013, A&A, 557, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reiners, A., & Schmitt, J. H. M. M. 2002a, A&A, 393, L77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reiners, A., & Schmitt, J. H. M. M. 2002b, A&A, 384, 155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reiners, A., & Schmitt, J. H. M. M. 2003, A&A, 398, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reinhold, T., Reiners, A., & Basri, G. 2013, A&A, 560, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ritzwoller, M. H., & Lavely, E. M. 1991, ApJ, 369, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [NASA ADS] [CrossRef] [Google Scholar]
 Roxburgh, I. W., & Vorontsov, S. V. 2003, A&A, 411, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rüdiger, G. 1974, Astron. Nachr., 295, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G. 1980, Geophys. Astrophys. Fluid Dyn., 16, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G. 1982, Geophys. Astrophys. Fluid Dyn., 21, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Schou, J., ChristensenDalsgaard, J., & Thompson, M. J. 1994, ApJ, 433, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390 [NASA ADS] [CrossRef] [Google Scholar]
 Semel, M. 1989, A&A, 225, 456 [NASA ADS] [Google Scholar]
 Silva Aguirre, V., Basu, S., Brandão, I. M., et al. 2013, ApJ, 769, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Takeda, Y., Sato, B., Kambe, E., et al. 2005, PASJ, 57, 13 [NASA ADS] [Google Scholar]
 Tassoul, J.L. 2000, Stellar Rotation, Cambridge Astrophysics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Thompson, M. J. 1991, in Challenges to Theories of the Structure of ModerateMass Stars, eds. D. Gough, & J. Toomre (Berlin: Springer Verlag), Lect. Notes Phys., 388, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Thompson, M. J., ChristensenDalsgaard, J., Miesch, M. S., & Toomre, J. 2003, ARA&A, 41, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Valio, A., Estrela, R., Netto, Y., Bravo, J. P., & de Medeiros, J. R. 2017, ApJ, 835, 294 [NASA ADS] [CrossRef] [Google Scholar]
 van Leeuwen, F. 2007, A&A, 474, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vogt, S. S., Penrod, G. D., & Hatzes, A. P. 1987, ApJ, 321, 496 [NASA ADS] [CrossRef] [Google Scholar]
 Waite, I. A., Marsden, S. C., Carter, B. D., et al. 2011, MNRAS, 413, 1949 [NASA ADS] [CrossRef] [Google Scholar]
 Waite, I. A., Marsden, S. C., Carter, B. D., et al. 2017, MNRAS, 465, 2076 [NASA ADS] [CrossRef] [Google Scholar]
 Walkowicz, L. M., & Basri, G. S. 2013, MNRAS, 436, 1883 [NASA ADS] [CrossRef] [Google Scholar]
 Wentzel, D. G. 1961, ApJ, 133, 170 [NASA ADS] [CrossRef] [Google Scholar]
 White, T. R., Huber, D., Maestro, V., et al. 2013, MNRAS, 433, 1262 [NASA ADS] [CrossRef] [Google Scholar]
 White, T. R., Benomar, O., Silva Aguirre, V., et al. 2016, A&A, 601, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Modelling 16 Cyg A and B
Nonseismic observational properties for the two stars of the 16 Cyg system.
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 hydrogenmass fraction, X_{0}, and initial metallicity, Z_{0}, and the mixinglength parameter, α. In the following, we group them in a single vector θ_{⋆} = (M_{⋆}, t_{⋆}, X_{0}, Z_{0}, α). 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.
Lower and upper bounds used for the prior uniform densities for each stellar parameter.
Stellar parameters (and heliummass fraction) inferred using the ASTEC stellar evolution code for 16 Cyg A and B.
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 highprecision 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
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 nonseismic 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_{01}, r_{02})), with r_{01} = (r_{01}(n_{01, 1}),…,r_{01}(n_{01, N})) and r_{02} = (r_{02}(n_{02, 1}),…,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 nonlinearly 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
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_{01} and r_{02} are correlated, therefore we treated these vector as two separate multivariate Gaussian densities 𝒩N(μ_{01}, Σ_{01}) and 𝒩N(μ_{02}, Σ_{02}). The covariance matrices were estimated numerically. Independent samples were generated for each frequency used for the evaluation of the components of r_{01} and r_{02}, and were used to obtain samples for both random vectors. Using these samples, evaluating Σ_{01} and Σ_{02} is straightforward. The resulting likelihood is therefore
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 lowT 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 ^{14}N(p,γ)^{15}O reaction. Convection was treated using the prescription from BöhmVitense (1958) for the mixinglength theory, the meanfree path of the fluid elements being proportional to the pressure scaleheight.
Fig. A.1. Marginal densities for the stellar parameters M, t_{⋆}, X_{0}, Z_{0}, 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. 

Open with DEXTER 
The stellar parameters were considered to be independent. Therefore, the prior can be written π(θ_{⋆}) = π(M_{⋆})π(t_{⋆})π(X_{0})π(Z_{0})π(α). 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 nonnegative. 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 trialanderror 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 hydrogenmass 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_{⋆}, X_{0}, Z_{0}, 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. 

Open with DEXTER 
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 lowprobability 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^{nT} 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 twodimensional and onedimensional 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 heliummass fraction, Y_{0}, which is an oftenused parameter in the literature. The MAP values computed from the fivedimensional joint density are given in the first line. They are given without uncertainties. We also give the MAP values obtained from each marginalised onedimensional 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
Lower and upper bounds used for the prior uniform densities for each stellar parameter.
Stellar parameters (and heliummass fraction) inferred using the ASTEC stellar evolution code for 16 Cyg A and B.
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 nondegenerate 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. 

Open with DEXTER  
In the text 
Fig. 2. Marginal densities for the spectrum parameters a_{1}, a_{3}, 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. 

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

Open with DEXTER  
In the text 
Fig. 4. Left panel: rotationrate 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 colatitude. 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. 

Open with DEXTER  
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 isoprobability levels of the corresponding threecomponent Gaussian mixture model. Right panel: projections on the Ω_{0} and Ω_{1} axis of the twodimensional 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 dotdashed line in the Ω_{1} projection shows the result using a normal distribution with the mean and variance of the MCMC sample. 

Open with DEXTER  
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 colatitude. 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 longdash line shows the associated 68.3% credible interval. The shortdash line shows the model corresponding to the MAP estimates of Ω_{0} and Ω_{1}. 

Open with DEXTER  
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. 

Open with DEXTER  
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 𝒦_{1}xs coefficients given by Eq. (19). The red dashed lines show the densities for our bestfit stellar model. 

Open with DEXTER  
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 Ammlervon Eiff & Reiners (2012). 

Open with DEXTER  
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 a_{1}. 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. 

Open with DEXTER  
In the text 
Fig. A.1. Marginal densities for the stellar parameters M, t_{⋆}, X_{0}, Z_{0}, 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. 

Open with DEXTER  
In the text 
Fig. A.2. Marginal densities for the stellar parameters M, t_{⋆}, X_{0}, Z_{0}, 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. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.