© ESO, 2017

1. Introduction

The impact of the Kepler space mission (Borucki et al. 2010; Koch et al. 2010) on diverse aspects of stellar astrophysics has been enormous and revolutionary by many standards. Our understanding of stellar evolution through asteroseismology has improved dramatically, and with this improvement, new challenges have appeared.

Sun-like stars, particularly red giants, exhibit a very rich pulsation pattern (De Ridder et al. 2009; Stello et al. 2009; Hekker et al. 2009). Some of these pulsations can be associated with pressure (p) modes, which are excited stochastically by turbulent convection. These p modes propagate throughout the star with the highest sensitivity to the external convective envelope, as opposed to the internal gravity (g) modes, which propagate only throughout the radiative core and hence are beyond observational reach. The p and g propagation zones generally do not overlap, and the region between them is called the evanescent zone. Modes of mixed character, behaving like g modes in the core and p modes in the envelope, bridge the evanescent zone and have substantial amplitudes in both the core and the envelope (Beck et al. 2011). These modes are extremely useful for obtaining information about the internal rotation of the core. This is possible because rotation induces frequency splittings in the modes, which would otherwise be degenerate in a spherically symmetric star (Ledoux 1951). Rotation induces a preferential axis in the star, and if the rotation rate is much smaller than the pulsation frequencies, then the splitting δ_nlm of a mode with radial, angular, and azimuthal wavenumbers n,l,m can be computed as (see, e.g., Aerts et al. 2010,  ) ${\mathit{\delta }}_{\mathit{nlm}}\mathrm{=}\mathit{m}\hspace{0.17em}{\mathit{\beta }}_{\mathit{nl}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{K}}_{\mathit{nl}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{\Omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r,}$(1)where the kernel K_nl(r) and β_nl are functions of the vertical and horizontal material displacement eigefunctions ξ_r(r) and ξ_h(r) (see Sect. 3 for more details). Thus, rotation lifts the azimuthal wave number degeneracy of the modes.

The milestone work of Beck et al. (2012) on red giants, using data obtained by the Kepler space telescope, took advantage of the detection of rotationally split mixed modes in KIC 8366239 and concluded that the stellar core spins about ten times faster than the envelope. In the same year, Mosser et al. (2012b) presented the core rotation of a sample of 300 red giants, establishing that the cores slow down significantly during the last stages of the red giant branch. A number of other studies followed. Most notably, Deheuvels et al. (2012) determined a core/envelope rotation rate ratio of about five for a star in the lower giant branch observed by Kepler, Deheuvels et al. (2014) computed the rotation rate ratio of six subgiants and young red giants, Deheuvels et al. (2015) obtained rotations rates for seven red giants in the secondary clump, and very recently, the work by Di Mauro et al. (2016) resolved the core rotation better than previous studies. In all these cases, the rotation rate of the core does not match the expectations of current angular momentum theories. Indeed, according to our current understanding of the evolution of the angular momentum in stellar interiors, the core is expected to spin up considerably as it contracts in stars at this stage of evolution. Cantiello et al. (2014) explicitly showed the inadequacy of current models in reproducing the observed slow rotation of the core in RGB stars, even after magnetic effects were included. It is then clear that a very effective mechanism for angular momentum transport is at work. Internal gravity waves are capable of transferring considerable amounts of angular momentum (Rogers 2015; Alvan et al. 2013), which provides a suitable explanation for “anomalous” rotation rates in other types of stars, such as those reported by Kurtz et al. (2014), Saio et al. (2015), and Triana et al. (2015). However, Fuller et al. (2014) showed that this mechanism falls short of explaining the observed rotation rate ratios in stars on the red giant branch. However, mixed modes can also transport angular momentum, as demonstrated recently by Belkacem et al. (2015). According to this study, the mixed-mode wave heat flux has an appreciable effect on the mean angular momentum in the inner regions of the star.

The methods used to obtain the internal rotation rates in red giant stars have seen a number of developments well worth mentioning here. The way to proceed, after the mode detection and identification process (see, e.g., Corsaro et al. 2015) is usually to develop a seismic model of the star with oscillation frequencies that are as similar as possible to the observed frequencies. This process is computationally expensive as a large number of evolutionary tracks with different stellar parameters need to be computed. A seismic model provides oscillation kernels that allow the application of inversion techniques (see Sect. 3) to determine the approximate rotation rates in different regions of the star. This is the approach taken by Deheuvels et al. (2012, 2014), Triana et al. (2015), and Di Mauro et al. (2016).

Goupil et al. (2013) developed a powerful method that allows estimating the rotation rates of both core and envelope, without recurring to any seismic model, by considering the relative amounts of “trapping” of a set of mixed modes, which as they showed are linearly related to the rotational splittings δ (see Sect. 2). The method relies solely on observed quantities as inputs and particularly, on an estimate of the asymptotic period spacing ΔΠ₁ that pure high-order g modes would have in the asymptotic regime (Mosser et al. 2012a).

This latter approach was taken by Deheuvels et al. (2015) in their sample of seven core He-burning red giants. In that study, a seismic model was obtained for one of the stars in the sample (KIC 7581399) and was used to compute rotation rates through inversions. Then the authors adapted the method of Goupil et al. (2013), and the rotation rates thus obtained were in very good agreement with the inversions based on the seismic model. After assessing the validity of the Goupil approach in this way, no further seismic modeling was attempted for the other targets, and the corresponding rotation rates reported come solely from the use of the adapted method.

Mosser et al. (2015) provided additional insight into the relationship between the relative amount of trapping of a mode (quantified by the parameter ζ) and the observed period spacing ΔP between consecutive mixed modes. The modification of Goupil’s formula introduced by Deheuvels et al. (2015) is exactly the same expression as was found by Mosser et al. (2015) to represent the ratio ΔP/ ΔΠ₁.

The model-independent method used by Goupil et al. (2013) has been compared with inversions based on seismic models for only two targets: the core helium-burning red giant KIC 7581399 (Deheuvels et al. 2015), and the early red giant KIC 4448777 (Di Mauro et al. 2016). Although a wider comparison of the method using targets in different evolutionary stages is desirable, we offer comparisons of themodel-independent method against inversions for our 13 targets, which share similar evolutionary stageswith KIC 4448777. Our targets are a subset of the original selection of 19 low-mass, low-luminosity red giants studied by Corsaro et al. (2015) using the Bayesian inference technique Diamonds (Corsaro & De Ridder 2014) to detect and identify pulsation frequencies. Pérez Hernández et al. (2016) performed a grid-based search for models designed specifically to constrain the age, mass, and initial helium content of all 19 targets. In the present work we search for optimal seismic models for the 13 targets that exhibit rotational splittings (Sect. 5) and employ Bayesian inference and inversion techniques (Sect. 3) to obtain average rotation rates. We compare these results with those obtained bythe model-independent method as implemented by Deheuvels et al. (2015) and by Mosser et al. (2015) using their expressions for the trapping parameter ζ (Sect. 6). Additionally, we use the idea proposed recently by Klion & Quataert (2016) that provides a way to localize the differential rotation of a red giant (whether in the radiative core or in the convective envelope of the star), provided that the rotation rate of the envelope is known by other means.

Fig. 1 Rotational kernels of two dipole modes with different values of the trapping parameter ζ (from the best model of KIC 007619745). The left panel shows the central region, and the right panel shows the region near the surface.

2. Rotation rate averages using the trapping parameter ζ

The parameter ζ gives an indication of how strongly a given stellar pulsation mode is localized, or “trapped”, inside the radiative core. It is defined as the ratio of the mode inertia computed within the g-mode cavity I_g (Goupil et al. 2013) to the total mode inertia I: $\mathit{\zeta }\mathrm{\equiv }\frac{{\mathit{I}}_{\mathit{g}}}{\mathit{I}}\mathrm{=}\frac{{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{1}}}^{{\mathit{r}}_{\mathrm{2}}}\mathit{\rho }{\mathit{r}}^{\mathrm{2}}\mathrm{[}{\mathit{\xi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{+}\mathit{l}\mathrm{(}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{)}{{\mathit{\xi }}_{\mathrm{h}}^{\mathrm{2}}}^{\mathrm{]}}\hspace{0.17em}\mathrm{d}\mathit{r}}{{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{R}}_{\mathrm{\ast }}}\mathit{\rho }{\mathit{r}}^{\mathrm{2}}\mathrm{[}{\mathit{\xi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{+}\mathit{l}\mathrm{(}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{)}{{\mathit{\xi }}_{\mathrm{h}}^{\mathrm{2}}}^{\mathrm{]}}\hspace{0.17em}\mathrm{d}\mathit{r}}\mathrm{·}$(2)In the expression above, ρ is the density, r₁ and r₂ are the turning points of the g-mode cavity, l is the angular wave number of the mode, ξ_r and ξ_h are the vertical and horizontal material displacement eigenfunctions, respectively, and R_∗ is the stellar radius (Aerts et al. 2010).

Figure 1 shows the rotational kernels K_nl of two mixed modes with different ζ using a seismic model for KIC 007619745. Modes with ζ close to one are gravity-dominated mixed modes, while modes with ζ close to 0.5 correspond to pressure-dominated mixed modes.

The rotational splittings δ_i are linearly related to ζ_i for each mode i, as shown by Goupil et al. (2013), which we also refer to for further details. The coefficients of this linear relationship are related directly to the rotational averagesacross the core and the envelope through $\mathit{\delta }\mathrm{=}\left(\frac{{\mathrm{\Omega }}_{\mathit{g}}}{\mathrm{2}}\mathrm{-}{\mathrm{\Omega }}_{\mathit{p}}\right)\mathit{\zeta }\mathrm{+}{\mathrm{\Omega }}_{\mathit{p}}\mathit{,}$(3)where Ω_p represents the average rotation rate in the envelope (approximated by the p-mode cavity), and Ω_g represents the average rotation rate in the radiative core (approximated by the g-mode cavity).Following Goupil et al. (2013), they are defined as ${\mathrm{\Omega }}_{\mathit{g}}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{r}}_{\mathrm{2}}}\mathit{K}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{\Omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r}}{{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{r}}_{\mathrm{2}}}\mathit{K}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r}}$(4)for the core, and ${\mathrm{\Omega }}_{\mathit{p}}\mathrm{=}\frac{{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{2}}}^{{\mathit{R}}_{\mathrm{\ast }}}\mathit{K}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{\Omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r}}{{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{2}}}^{{\mathit{R}}_{\mathrm{\ast }}}\mathit{K}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r}}$(5)for the envelope, where r₂ is the outer turning point in the g resonant cavity and K(r) is the rotational kernel of a mixed-mode. Our tests with the rotational profiles considered in Sect. 4 show that while the core averages are essentially independent of the particular mixed-mode chosen, the envelope averages differ appreciably across modes with ζ ≳ 0.9 (i.e., gravity-dominated modes). Using kernels from mixed modes with ζ ≲ 0.85 results in envelope averages with minimal variability.

2.1. Estimation of the trapping parameter ζ

The expression for ζ given by Eq. (2)in principle requires knowing the material displacement eigenfunctions ξ_r,h(r) for each mode, which are only available after the computationally expensive process of deriving a seismic model of the star. However, Goupil et al. (2013) used an asymptotic analysis method based on the work of Shibahashi (1979) to show that ζ can in principle be estimated using observational data alone. The expression for ζ was later refined by Deheuvels et al. (2015) and Mosser et al. (2015). In what follows, we briefly recall the method and the main formulae, we refer to the original works for further details.

Fig. 2 Trapping parameter as described by the true model, ζmod (Eq. (2)) as a function of the true frequencies νmod (dashed black line and black circles). Red circles represent the trapping parameter ζas as approximated by Eqs. (11)and (8)and using νas for the mode frequencies as determined by Eq. (6). The scaled period spacing ΔP/ ΔΠ1 between any two consecutive true frequencies νmod,1 and νmod,2 is plotted at the abscissa ${\mathit{\nu }}_{\mathrm{\Delta }\mathit{P}}\mathrm{=}\mathrm{2}\mathit{/}\mathrm{(}{\mathit{\nu }}_{\mathrm{mod}\mathit{,}\mathrm{1}}^{-1}\mathrm{+}{\mathit{\nu }}_{\mathrm{mod}\mathit{,}\mathrm{2}}^{-1}\mathrm{)}$ (blue crosses). ΔPν represent the interpolated value of the two adjacent period spacings around each frequency νmod (green crosses).

The method consists of finding aproximate JWKB solutions for the material displacement eigenfunctions ξ_r,h(r) in the two separate p and g cavities of the star. Matching the solutions in the evanescent zone requires (Unno 1989) $\tan {\mathit{\theta }}_{\mathit{p}}\mathrm{=}\mathit{q}\hspace{0.17em}\tan {\mathit{\theta }}_{\mathit{g}}\mathit{,}$(6)where q is the coupling constant between the p- and g-mode cavities, and the phases θ_p,g are defined through ${\mathit{\theta }}_{\mathit{g}}\mathrm{=}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{1}}}^{{\mathit{r}}_{\mathrm{2}}}{\mathit{k}}_{\mathrm{r}}\hspace{0.17em}\mathrm{d}\mathit{r,}\hspace{0.17em}\hspace{0.17em}{\mathit{\theta }}_{\mathit{p}}\mathrm{=}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{3}}}^{{\mathit{r}}_{\mathrm{4}}}{\mathit{k}}_{\mathrm{r}}\hspace{0.17em}\mathrm{d}\mathit{r,}$(7)where r₃ and r₄ are the inner and outer turning points of the p-mode cavity. According to Mosser et al. (2012a), asymptotic analysis yields ${\mathit{\theta }}_{\mathit{p}}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{\Delta }\mathit{\nu }}\mathrm{(}\mathit{\nu }\mathrm{-}{\mathit{\nu }}_{\mathit{p}}\mathrm{)}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{\theta }}_{\mathit{g}}\mathrm{=}\mathit{\pi }\left(\frac{\mathrm{1}}{\mathrm{\Delta }{\mathrm{\Pi }}_{\mathrm{1}}\hspace{0.17em}\mathit{\nu }}\mathrm{-}{\mathit{ϵ}}_{\mathit{g}}\right)\mathit{,}$(8)where ν_p is the frequency of the theoretical l = 1 pure p modes, which are related to the radial (l = 0) modes ν_n,0 through ${\mathit{\nu }}_{\mathit{p}}\mathrm{=}{\mathit{\nu }}_{\mathit{n,}\mathrm{0}}\mathrm{+}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}{\mathit{d}}_{\mathrm{01}}\right)\mathrm{\Delta }\mathit{\nu }\mathit{.}$(9)In turn, the radial modes can be expressed as a function of the radial order n involving the parameters ϵ_p, α, and the large frequency separation Δν as follows (Mosser et al. 2013): ${\mathit{\nu }}_{\mathit{n,}\mathrm{0}}\mathrm{=}\left[\mathit{n}\mathrm{+}{\mathit{ϵ}}_{\mathit{p}}\mathrm{+}\frac{\mathit{\alpha }}{\mathrm{2}}\mathrm{(}\mathit{n}\mathrm{-}{\mathit{n}}_{\max }{\mathrm{)}}^{\mathrm{2}}\right]\mathrm{\Delta }\mathit{\nu ,}$(10)where n_max ≡ ν_max/ Δν−ϵ_p. The approximate expression for the trapping parameter, denoted here as ζ_as, reads ${\mathit{\zeta }}_{\mathrm{as}}\mathrm{=}{\left[\mathrm{1}\mathrm{+}\frac{{\mathit{\nu }}^{\mathrm{2}}\hspace{0.17em}\mathrm{\Delta }{\mathrm{\Pi }}_{\mathrm{1}}}{\mathit{q}\hspace{0.17em}\mathrm{\Delta }\mathit{\nu }}\frac{{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{g}}}{{\cos }^{\mathrm{2}}{\mathit{\theta }}_{\mathit{p}}}\right]}^{-1}\mathrm{·}$(11)A total of seven parameters are required here: the coupling constant q, the offsets ϵ_p,g, the mean large frequency separation Δν, the asymptotic period spacing ΔΠ₁, α, and d₀₁. In practice, the optimal parameters ϵ_p, α, and Δν are determined first by fitting the observed radial modes to Eq. (10). Then, the optimal parameters q, ΔΠ₁, ϵ_g, and d₀₁ that best reproduce the observed l = 1 mode frequencies can be found by a downhill simplex method. This requires solving Eq. (6)for the mode frequencies ν at each search step with a particular (q, ΔΠ₁,ϵ_g,d₀₁) combination, using a Newton method to find the roots of the equation, for example.

The observed rotational splittings δ are expected to be linearly related to ζ_as , and a linear fit of δ as a function of ζ_as leads to an estimate of the average envelope rotation and the envelope core rotation through Eq. (3).

Mosser et al. (2015) obtained a result in their search of an expression for the mixed-mode relative period spacings ΔP/ ΔΠ₁ that exactly matched the expression for ζ_as (Eq. (11)) found by Deheuvels et al. (2015). Thus, we have ${\mathit{\zeta }}_{\mathrm{as}}\mathrm{=}\frac{\mathrm{\Delta }\mathit{P}}{\mathrm{\Delta }{\mathrm{\Pi }}_{\mathrm{1}}}\mathrm{·}$(12)The above is a reflection of the fact that the rotational splittings follow the same distribution as the period spacing, because both are determined by the coupling between pressure and gravity terms. As a bonus, Eq. (12)provides a simple and direct way to estimate ζ from observations without the need of the optimal parameters mentioned above (which can be slightly time consuming computationally) with the exception of ΔΠ₁. Some care must be taken because the period spacing between two consecutive mixed dipole modes $\mathrm{\Delta }\mathit{P}\mathrm{(}\mathit{n,n}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{=}{\mathit{\nu }}_{\mathit{n}\mathrm{+}\mathrm{1}}^{-1}\mathrm{-}{\mathit{\nu }}_{\mathit{n}}^{-1}$ is defined properly at $\mathit{\nu }\mathrm{=}\mathrm{2}\mathit{/}\mathrm{(}{\mathit{\nu }}_{\mathit{n}\mathrm{+}\mathrm{1}}^{-1}\mathrm{+}{\mathit{\nu }}_{\mathit{n}}^{-1}\mathrm{)}\mathrm{\equiv }{\mathit{\nu }}_{\mathrm{\Delta }\mathit{P}}$. In order to assign a ΔP to each mode ν_n, we therefore interpolate the two adjacent period spacings ΔP(n,n + 1) and ΔP(n−1,n) linearly. Similarly, when performing linear fits of δ vs. ζ_as, it is advisable to also include the interpolated rotational splitting at each location ν_ΔP using the two correspondingly adjacent values of δ to minimize biases, see Fig. 4.

In Fig. 2 we show that the trapping parameter ζ_mod as derived from a known seismic model is indeed well approximated by either ζ_as derived using asymptotic analysis, Eq. (11), or by the simpler expression given by Eq. (12).

Assuming that the errors on the frequencies and the splittings are normally distributed, we can sample randomly from them and proceed to compute interpolated splittings and spacings as explained earlier. A linear fit to these points leads to estimates of Ω_g and Ω_p according to Eq. (3). By repeating these steps many times, we can obtain the distributions associated with Ω_g and Ω_p and their associated errors.

2.2. Bayesian inference

If a seismic model providing the oscillation eigenfunctions ξ_r,h is available, the trapping parameter can be computed from Eq. (2), which we denote now as ζ_mod. With two sets of inputs, that is, the splittings δ_i and the trapping parameters ζ_mod,i, we can set out to perform a Bayesian fit using Eq. (3)as the model. To accomplish this, we first compute a Gaussian log-likelihood function defined as (see also Corsaro et al. 2013) $\mathrm{\Lambda }\mathrm{(}{\mathrm{\Omega }}_{\mathit{g}}\mathit{,}\hspace{0.17em}{\mathrm{\Omega }}_{\mathit{p}}\mathrm{)}\mathrm{=}{\mathrm{\Lambda }}_{\mathrm{0}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\left[\frac{{\mathrm{\Delta }}_{\mathit{i}}\mathrm{(}{\mathrm{\Omega }}_{\mathit{g}}\hspace{0.17em}\mathit{,}{\mathrm{\Omega }}_{\mathit{p}}\mathrm{)}}{{\mathit{ϵ}}_{\mathit{i}}}\right]}^{\mathrm{2}}\hspace{0.17em}\mathit{,}$(13)where N is the total number of rotational splittings δ_i from the observations (one for each m-multiplet of mixed modes), Δ_i(Ω_g, Ω_p) are the residuals given as the difference between the observed and the modeled splittings, ϵ_i the corresponding uncertainty, and ${\mathrm{\Lambda }}_{\mathrm{0}}\mathrm{=}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\ln \sqrt{\mathrm{2}\mathit{\pi }}\hspace{0.17em}{\mathit{ϵ}}_{\mathit{i}}\hspace{0.17em}\mathit{,}$(14)a constant term. We multiply the likelihood distribution by the prior distributions (uniform or flat in this case, in the range Ω_g ∈ [2,8] μHz, and Ω_p ∈ [0,3] μHz), obtaining a posterior probability density distribution. By marginalizing the bidimensional posterior into two one-dimensional probability density distributions, we obtain estimates for Ω_p and Ω_g that are the medians of the two one-dimensional distributions. The corresponding error bars are the Bayesian credible intervals computed as explained in Corsaro et al. (2013). This statistical approach may provide similar results to a least-squares fit, but it is conceptually very different. One of the main differences is that it is able to incorporate any a priori knowledge on the estimated parameters that we may have, for instance, in the form of the prior distributions described above.

An example of this Bayesian fit is shown in Fig. 3 for the star KIC 007619745 (solid orange line), with 1σ error bars overlaid. The results from this method for all targets are included in Table 1 under the “Bayes” heading.

Fig. 3 Illustrative use of a Bayesian fit to compute rotation rates for KIC 007619745. The trapping parameters ζmod are obtained from the eigenmodes as computed by GYRE.

3. Inversion methods

In this section we present the methods we used to obtain rotation rate averages, which are all based on the oscillation kernels provided by the seismic models described in Sect. 5. Our treatment is based on Aerts et al. (2010), who offered an extended presentation of the methods discussed below.

The so-called forward problem states that the rotational splittings δ_nlm of an oscillation mode with radial, angular, and azimuthal wavenumbers n,l,m can be computed through Eq. (1). Explicitely, the kernels are computed from the material displacement eigenfunctions via ${\mathit{K}}_{\mathit{nl}}\mathrm{=}\frac{\mathrm{1}}{\mathit{I}}\mathrm{[}{\mathit{\xi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{+}\mathit{l}\mathrm{(}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}{\mathit{\xi }}_{\mathrm{h}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{\xi }}_{\mathrm{r}}{\mathit{\xi }}_{\mathrm{h}}\mathrm{-}{{\mathit{\xi }}_{\mathrm{h}}^{\mathrm{2}}}^{\mathrm{]}}{\mathit{r}}^{\mathrm{2}}\mathit{\rho ,}$(15)where I is the mode inertia (see Eq. (2)). The constant β_nl is given by ${\mathit{\beta }}_{\mathit{nl}}\mathrm{=}\frac{\mathrm{1}}{\mathit{I}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{R}}_{\mathrm{\ast }}}\mathrm{[}{\mathit{\xi }}_{\mathrm{r}}^{\mathrm{2}}\mathrm{+}\mathit{l}\mathrm{(}\mathit{l}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\xi }}_{\mathrm{h}}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{\xi }}_{\mathrm{r}}{\mathit{\xi }}_{\mathrm{h}}\mathrm{-}{{\mathit{\xi }}_{\mathrm{h}}^{\mathrm{2}}}^{\mathrm{]}}{\mathit{r}}^{\mathrm{2}}\mathit{\rho }\hspace{0.17em}\mathrm{d}\mathit{r}\mathit{.}$(16)The inverse problem consists of determining the unknown inversion coefficients c_i(r) satisfying $\mathrm{\Omega ̅}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{M}}{\mathit{c}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\frac{{\mathit{\delta }}_{\mathit{i}}}{\mathit{m}\hspace{0.17em}{\mathit{\beta }}_{\mathit{i}}}\mathit{,}$(17)where is the predicted internal rotation rate of the star, M is the number of observed splittings, and i denotes the collective indices (n,l,m).Clearly, the inversion coefficients c_i are not determined by Eq. (17), which just states a linear relationship between the observed splittings and the predicted rotation profile. The c_i(r) are determined by minimizing the difference between observed and predicted splittings, by minimizing of the resulting uncertainties, or by adjusting the shape of the averaging kernels, as discussed below.

The approximate rotational profile can be expressed in terms of the true profile Ω(r) by means of the averaging kernels , which are related to the kernels K_i(r) through $\mathrm{𝒦}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,r}\mathrm{\right)}\mathrm{=}{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{M}}{\mathit{c}}_{\mathit{i}}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{K}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ and fulfill $\mathrm{\Omega ̅}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{R}}_{\mathrm{\ast }}}\mathrm{𝒦}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,r}\mathrm{\right)}\hspace{0.17em}\mathrm{\Omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r}\mathit{.}$(18)The averaging kernels should be localized around r′ as much as possible, ideally resembling a delta function δ(r′,r). It is usually assumed that the observational errors ϵ_i are uncorrelated (as in, e.g., Deheuvels et al. 2012 or Di Mauro et al. 2016), so that the variance of the predicted rotation rates can be estimated as ${\mathit{\sigma }}^{\mathrm{2}}\left[\mathrm{\Omega ̅}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right]\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{M}}{\mathit{c}}_{\mathit{i}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}{\left(\frac{{\mathit{ϵ}}_{\mathit{i}}}{{\mathit{\beta }}_{\mathit{i}}}\right)}^{\mathrm{2}}\mathrm{·}$(19)The expression above accounts for the errors originating from the observations alone; it does not account for the inherent errors of the inversion process itself.

3.1. Two-zone inversion models

To obtain approximate averages of the core and envelope rotation rates, we can make use of simple two-zone models where we assume an inner zone extending from the stellar center to r/R_∗ = x_c and an outer zone extending from r/R_∗ = x_c all the way to stellar surface at r/R_∗ = 1, both zones rotating uniformly with rates Ω_g and Ω_p, respectively. Our 13 targets happen to be approximately at the same evolutionary stage, and therefore it is not surprising that our seismic models show all their evanescent zones located at approximately the same radial locations (scaled by stellar radius) We have chosen x_c to coincide with the base of the convection zone for each target, see Fig. 4. We can determine the inversion coefficients c_i associated with each zone by finding the optimal Ω_g and Ω_p that minimize ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{M}}{\left(\frac{\mathit{\delta ̅}\mathit{i}\mathrm{-}{\mathit{\delta }}_{\mathit{i}}}{{\mathit{ϵ}}_{\mathit{i}}}\right)}^{\mathrm{2}}\mathit{,}$(20)where ϵ_i are the observation errors and are the predicted splittings associated with thetwo-zone rotation profile composed of Ω_g and Ω_p.The averages are determined by enforcing ∂(χ²) /∂Ω_g,p = 0 after substituting χ² using Eqs. (1)and (20). Results from this method, with kernels from the seismic models discussed below, are presented on Table 1 under the “Two-zone” heading.

Fig. 4 Propagation diagram of the best model of KIC 007619745. The Brunt-Väisälä frequency N(r) is plotted as a solid blue curve, and the l = 1 Lamb frequency S1(r) is plotted as a dashed green curve. p-modes with frequency ν are evanescent in the region where ν < S1(r). Similarly, g-modes are evanescent wherever ν >N(r). The core-envelope boundary in our two-zone models is located at rc ≈ 0.1R∗, which usually lies in the evanescent zone common to both p- and g-modes. The peak of N(r) near r/R∗ = 0.11 is associated with a chemically inhomogeneous layer in the radiative zone, left behind by convection as the boundary moves upward.

3.2. Subtractive optimally localized averaging

One of the differences between the subtractive optimally localized averaging (SOLA) method (Pijpers & Thompson 1994) and the method described above is that with SOLA we do not minimize χ², but instead, the method chooses the optimal linear combination of the inversion coefficients c_i such that the averaging kernels resemble a given target function T(r′,r) as closely as possible while keeping the variance low. Thus, we minimize ${\mathrm{\int }}_{\mathrm{0}}^{{\mathit{R}}_{\mathrm{\ast }}}{\left[\mathrm{𝒦}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,r}\mathrm{\right)}\mathrm{-}\mathit{T}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,r}\mathrm{\right)}\right]}^{\mathrm{2}}\mathrm{d}\mathit{r}\mathrm{+}\mathit{\mu }\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{M}}{\mathit{c}}_{\mathit{i}}^{\mathrm{2}}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{ϵ}}_{\mathit{i}}^{\mathrm{2}}$(21)at each r′, with the additional constraint ${{}^{\mathrm{\int }}}_{\mathrm{0}}^{{\mathit{R}}_{\mathrm{\ast }}}\mathrm{𝒦}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r}\mathrm{=}\mathrm{1}$. The target function that we have chosen in this study is a Gaussian with unit norm, centered on r = r′ with adjustable width s: $\mathit{T}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,r}\mathrm{\right)}\mathrm{=}\mathit{N}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\left(\frac{{\mathit{r}}^{\mathrm{\prime }}\mathrm{-}\mathit{r}}{\mathit{s}}\right)}^{\mathrm{2}}}\mathit{,}$(22)N being a normalization factor. In addition to the free parameter μ in Eq. (21), we can also adjust the shape of the target function T by adjusting the width s. The problem reduces to solving the linear set of M equations (i = 1,...,M) for each radial location r′: $\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{M}}{\mathit{W}}_{\mathit{ik}}\hspace{0.17em}{\mathit{c}}_{\mathit{k}}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{K}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathit{T}\mathrm{\left(}{\mathit{r}}^{\mathrm{\prime }}\mathit{,r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r,}$(23)where ${\mathit{W}}_{\mathit{ik}}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{{\mathit{R}}_{\mathrm{\ast }}}{\mathit{K}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}{\mathit{K}}_{\mathit{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r}\mathrm{+}\mathit{\mu }\hspace{0.17em}{\mathit{\delta }}_{\mathit{ik}}\hspace{0.17em}{\mathit{ϵ}}_{\mathit{i}}^{\mathrm{2}}$, together with the constraint ∑ _kc_k(r′) = 1, which is implemented via Lagrange multipliers. Given a set of kernels K_i(r) and the two parameters (s,μ), the inversion coefficients c_i(r′) are determined by solving the set of M equations given by Eq. (23). We note that the observed splittings are not involved in determining the c_i(r′), they determine the predicted rotation rate via Eq. (17).

4. Testing the methods

Before applying the methods described above to determine rotation rates for the stars in our sample, it is desirable to have an idea of how the methods perform under controlled situations. We choose a specific seismic model (of KIC 007619745, obtained as explained in the next section) as the “true” model. Then, we consider six different rotation profiles and compute the exact rotational splittings in each case via Eq. (1). We also compute the “true” rotational averages for both the core and the envelope for each case using Eqs. (4)and (5).

For the first set of rotation profiles we adopted the functional form used by Klion & Quataert (2016): we assume that the inner region of the star rotates uniformly with a rate Ω_c from the center and up to 1.5 times r_H, the outer radius of the hydrogen burning shell. Then, from 1.5 r_H and up to the base of the convective zone (r_rcb), the star follows a uniform rotation rate Ω_m. The remainder of the star rotates according a power-law profile. Thus $\mathrm{\Omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathit{r}\mathrm{\le }\mathrm{1.5}\hspace{0.17em}{\mathit{r}}_{\mathrm{H}}\mathit{,}\\ & \mathrm{1.5}\hspace{0.17em}{\mathit{r}}_{\mathrm{H}}\mathit{<}\mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{rcb}}\mathit{,}\\ & \mathit{r}\mathit{>}{\mathit{r}}_{\mathrm{rcb}}\mathit{,}\end{array}$(24)where $\mathit{\alpha }\mathrm{=}\frac{\log \mathrm{(}{\mathrm{\Omega }}_{\mathrm{m}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{e}}\mathrm{)}}{\log \mathrm{(}{\mathit{R}}_{\mathrm{\ast }}\mathit{/}{\mathit{r}}_{\mathrm{rcb}}\mathrm{)}}\mathrm{·}$(25)The exponent α is so chosen to ensure the continuity of Ω(r) at r = r_rcb). This functional form is useful to adjust the location of the differential rotation. By setting Ω_m = Ω_e all the differential rotation in the star is localized at r_rcb, inside the radiative region, and if Ω_m = Ω_c, the differential rotation is all contained in the convective envelope. We have kept Ω_c and Ω_e fixed at 0.7 μHz and 0.1 μHz, respectively. We consider three different values for Ω_m: 0.7 μHz, 0.4 μHz, and 0.1 μHz.

The other set of rotation profiles are Gaussians of different widths s plus a constant term B: $\mathrm{\Omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{A}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{\ast }}}{\mathit{s}}\right)}^{\mathrm{2}}}\mathrm{+}\mathit{B,}$(26)we set A = 0.6 μHz, B = 0.1 μHz and three different widths s: 0.01, 0.05 and 0.2. See Fig. 5.

Fig. 5 Synthetic rotation rate profiles used for testing. Ωm refers to the rotation rate between 1.5 rH and the base convective zone, see Eq. (24). s refers to the width of a Gaussian profile as defined in Eq. (26).

Fig. 6 Rotational splittings using the best seismic model of KIC 007619745 for one of the test profiles as defined in Eq. (24)with Ωc = 0.7 μHz, Ωm = 0.4 μHz, and Ωe = 0.1 μHz. The splittings δΔP (dark blue crosses) are the interpolated splittings as explained in Sect. 2.1. Dotted lines are linear fits providing estimates of the rotation rate averages in the g and p cavities following Eq. (3).

We compute the rotational splittings of the six test profiles using Eq. (1)for all the l = 1 mode frequencies of the true model within ± 3 Δν from the frequency at maximum power ν_max. Then we compute the optimal combination of parameters (q, ΔΠ₁,ϵ_g,d₀₁) that best reproduces these model frequencies in order to estimate ζ_as. As explained in Sect. 2.1, we compute interpolated splittings δ_ΔP to correspond with each period spacing ΔP as well as interpolated period spacings ΔP_ν to correspond with each splitting δ, see Fig. 6.

Fig. 7 Comparison of the predicted rotation rate averages in the core (Ωg) and in the envelope (Ωp) using the methods described in Sects. 2 and 3 based on the best seismic model for KIC 007619745. The top row corresponds to rotation profiles following Eq. (24), and the bottom row corresponds to Gaussian rotation profiles according to Eq. (26). The black dashed lines mark the true averaged rotation rates as computed from Eqs. (4)and (5). Red represents the rotation rates as determined from linear fits of (ζmod,δ), and magenta represents the rotation rates from linear fits of (ζas,δ). Blue corresponds to linear fits of δ and δΔP vs. ΔPν/ ΔΠ1 and ΔP/ ΔΠ1. Black are the predictions from two-zone inversions, and green are the predictions from SOLA inversions. Results from Bayesian inference are not included here as they are essentially identical to the two-zone models.

Now we obtain the estimates of the rotation rates using each of the methods explained earlier. The true averages are obtained following Eqs. (4)and (5). In the case of the SOLA method, we computed the predicted rotation rates in two different locations: one at the surface of the star (r/R_∗ = 1), and the other well within the radiative core (r/R_∗ = 10^-3). These SOLA-predicted rotation rates are sensitive to the width s of the target function. To obtain “calibrated” values for s, we therefore proceed first to compute the optimal two-zone model to determine the averages Ω_g and Ω_p. Then we compute the rotational splittings associated with this two-zone model and use them as inputs for a SOLA inversion, adjusting the widths s as necessary to make the corresponding SOLA predictions at r/R_∗ = 10^-3 and r/R_∗ = 1 exactly match Ω_g and Ω_p, respectively. With the widths s determined in this way, we then proceed to compute the SOLA-predicted averages using the splittings arising from the six test profiles.

We assume that the splittings all have the same uncorrelated Gaussian error distribution whose width matches the mean 1σ error from the actual measurements. In the case of the inversions, the corresponding uncertainties on the predicted rotation rates are computed via Eq. (19). In the case of the rotation rates estimated via linear fits involving ζ_{mod, as} (see Fig. 6), the resulting uncertainties are computed through a Monte Carlo simulation sampling randomly and repeatedly for many times (10⁵) from the normal distributions associated with each spliting δ. For the fits of δ,δ_ΔP vs. ΔP_ν/ ΔΠ₁, ΔP/ ΔΠ₁ we followed the same Monte Carlo approach, except that we also sampled randomly from the distributions associated each mode frequency ν, which we set as having the same standard deviation as the actual observed errors. We note that errors on the trapping parameter or systematic errors incurred by the inversions are not considered.

The predicted rotation rates are shown in Fig. 7. Although some scatter is present, the predictions are in acceptable agreement with the true averages across all methods. The test we just performed is an ideal situation, the modes are properly identified, the rotational kernels are based on the true seismic model, and the error distributions associated with the frequencies and splittings are centered exactly on the true values. Any deviation from this ideal situation will bring additional scatter to the predicted rotation averages.

The averaging kernels from the SOLA inversions at the two radial locations mentioned above, although well resolved with respect to each other, are essentially identical to the averaging kernels of the two-zone models at the corresponding zone, which suggests already that no more than two reasonably well-defined averages can be obtained using SOLA inversions.

Recently, Klion & Quataert (2016) proposed a method aiming to determine the region in a red giant star where differential rotation is concentrated. They note that the minimum normalized splitting, min(δ/ max(δ)) can be used to distinguish between rotation profiles with differential rotation localized in the core from those with differential rotation localized mostly in the envelope. This is precisely the motivation of the functional form of the rotation profiles given by Eq. (24). For a fixed ratio Ω_c/ Ω_e in these profiles, the quantity min(δ/ max(δ)) follows a one-to-one correspondence with Ω_m that determines the location of the differential rotation either in the core or in the envelope.

Di Mauro et al. (2016) were able to resolve three rotation rates in three distinct radial locations of KIC 4448777, two of them within the radiative core, which allowed them to conclude that there is a steep gradient in the rotation there, thus localizing the differential rotation of the star inside the radiative core. The target in that study is very similar to the targets in our target selection, sharing similar evolutionary stages, therefore it may be possible in principle to resolve the rotation rate in at least two points inside the core. Unfortunately, the set of rotational kernels for all of our targets is not suitable for obtaining more than an averaged value across the core (see Fig. 11). The reason for this is not evident a priori, and although it deserves special attention, it is beyond the scope of the present study.

Fig. 8 Histogram associated with the predictions for (blue, green, and red tones, respectively). The true (test) profile is described by Eq. (24)with fixed Ωc = 0.7 μHz, Ωe = 0.1 μHz. The light colors represent the histogram obtained when Ωm = 0.1 μHz, medium colors correspond to Ωm = 0.4 μHz, and dark colors represent Ωm = 0.7 μHz. Essentially identical distributions for result, regardless of the choice of Ωm. While Ωc and Ωe are recovered properly, no good prediction for Ωm can be made.

We can use the idea by Klion & Quataert (2016) to determine whether we can indeed localize the differential rotation. We consider a rotation profile with some preestablished values for the set (Ω_c,Ω_m,Ω_e) and compute its associated splittings. We consider splittings from all l = 1 mixed modes within ± 3 Δν of ν_max. With this set of splittings as input and considering errors on them matching the actual observed errors, we proceed to find the optimum combination of predicted that minimizes the difference between the input and the predicted splittings. To minimize and obtain estimates of the three parameters, we use a downhill simplex method. Then we make use once again of a Monte Carlo simulation sampling randomly from the normal distributions associated with the input splittings to obtain distributions for . If enough information is indeed contained in the set of input splittings to constrain , then this should reflect in a sharply peaked distribution for it. Figure 8 shows the results of this experiment. We set Ω_c = 0.7 μHz, Ω_e = 0.1 μHz, and vary Ω_m ∈ { 0.1, 0.4, 0.7 } μHz. The predicted values for and are close to the true values, but the probability density distribution for is wide and practically flat, thus no reliable prediction for Ω_m is possible and the differential rotation cannot be localized properly. Essentially identical distributions of are obtained regardless of the choice of Ω_m. This is a consequence of the magnitude of errors in the splittings together with the characteristics of the rotational kernels. Only if we artificially reduce the errors by an order of magnitude or less, a reasonable value of can be recovered. All of the 13 targets in our sample exhibit this undesirable behavior.

5. Seismic modeling

From the original 19 young red giants studied by Corsaro et al. (2015), we selected only those stars that exhibit rotationally split dipole modes (triplets). In some cases, depending on the stellar inclination angle, some of the l = 1 triplets were missing their central m = 0 peaks although the split m = ± 1 components were clearly visible. In these cases we assumed a central m = 0 component in the middle of the observed m = ± 1 frequencies, and we associated with it an error equal to three times the mean frequency error of the l = 0 peaks (usually larger than the error on the m = ± 1 components). This choice is conservative given that the asymmetry present in the full triplets (i.e., those that show central peaks) is usually smaller than this error.

Assuming the presence of a central peak with this frequency uncertainty has virtually no effect on the uncertainty on the inferred rotation rates, given that in these cases the splittings are computed simply as half the distance in frequency of the m = ± 1 components, without involving the hypothetical central component.

Equipped with these sets of pulsation frequencies, we set out to find approximate seismic models for each target making use of the MESA stellar evolution suite (Paxton et al. 2011, 2013, 2015) together with the GYRE pulsation code (Townsend & Teitler 2013). The MESA suite includes the “astero” module, which implements a downhill simplex search method (Nelder & Mead 1965) to obtain the best stellar parameters given a set of pulsation and spectroscopic data. To reduce computing time during the search, we opted to include only the observed radial (l = 0) modes and the dipole (l = 1) modes. Including higher l modes in the search, however desirable, would prohibitively increase the time required to find suitable seismic models for all targets.

Our approach consisted of a combination of grid and downhill simplex searches. We set up a grid of initial metallicities [Fe/H]_ini, varying from –0.22 to 0.2 with 0.06 steps (using a reference solar metallicity Z_⊙/X_⊙ = 0.02293 (Grevesse & Sauval 1998)) and initial helium content Y_ini varying from 0.23 to 0.3 with 0.005 steps. For each pair ([Fe/H]_ini, Y_ini) we performed a downhill simplex search optimizing for initial mass and overshoot (f_ov, expressed as a fraction of the pressure scale height H_P), in addition to age. We have kept the mixing length parameter α_MLT fixed at its solar calibrated value of 1.9. We also used Eddington-gray atmospheres and adopted the Asplund et al. (2009) mixture together with OPAL opacity tables (Iglesias & Rogers 1996). When computing mode frequencies, we peformed atmospheric corrections following the method by Kjeldsen et al. (2008) using a calibrated value of the exponent b = 4.81 as reported by Ball & Gizon (2014). The exact value of b is not critical as the observed modes can still be identified one-to-one to model frequencies using slightly different values. Results are summarized in Table 1.

We performed a hare-and-hounds exercise to test our grid + downhill simplex approach. For this we extracted mode frequencies from the best model of KIC 007619745, added some noise and used them as inputs to our search algorithm. The resulting model was satisfactorily close to the original, especially regarding the rotational kernels derived from them. We present more details in Appendix A.

6. Internal rotation rates

We estimated the internal rotation rates using two-zone models, Bayesian inference, and SOLA inversions (all based on the kernels provided by the seismic models) in addition to the model-independent method of Goupil et al. (2013) as described in Sect. 2. Figure 9 illustrates the use of the trapping parameter and associated linear fits as applied to the target KIC 007619745 as an example, as well as to the hare-and-hounds exercise explained in Appendix A.

Fig. 9 Left column plots show the observed rotational splittings, normalized by their maximum value (red error bars) and the observed period spacings, normalized by the observed period spacing ΔΠ1 (blue error bars) as functions of mode frequency for the target KIC 007619745 (top), and for the hare-and-hounds exercise. Right column plots show the observed rotational splitting δ as a function of ζmod (red points and error bars), δ vs. ζas (dark green points and error bars), and both δ vs. ΔPν/ ΔΠ1 and δΔP vs. ΔP/ ΔΠ1 (blue points and cyan error bars). See Sect. 2.1.

We computed two-zone models with the inter-zone boundary in the middle of the evanescent zone for each target. For the SOLA inversions we set the trade-off parameter μ to zero since it made negligible difference on the results. In addition, for SOLA, we computed rotation rates at two radial locations, one at r/R_∗ = 10^-3, deep in the radiative cores, and the other at the surface, r/R_∗ = 1. The estimates of the rotation rates through SOLA do depend slightly on the width s of the target functions used. To select the optimal value of s, we took the two-zone rotation rates determined earlier and computed their associated splittings via Eq. (1). Then, using these synthetic splittings, we iteratively determined the optimal s values required to exactly reproduce the two-zone rotation rates. With the widths s thus determined, we then computed the rotation rates at r/R_∗ = 10^-3 and at r/R_∗ = 1 using the observed rotational splittings. The rotation rates estimated in this way are presented in Table 2 and in Figs. B.1 and B.1.

We assume that the observational errors on the splittings and on the mode frequencies are normally distributed. Since in general, a given mode frequency has two period spacings associated with it and each period spacing is associated with two mode frequencies (and therefore two rotational splittings), we opted to take interpolated values as explained in Sect. 2.1. We performed linear fits to the resulting set to estimate the average rotation rates of the core Ω_g and the envelope Ω_p according to Eq. (3). A straightforward Monte Carlo approach using the observed rotational splittings and frequencies, together with their normally distributed errors, reveals a correlation between Ω_g and Ω_p that is a reflection of the fact that in Eq. (3)the slope and the intercept are not independent of each other. The rotation rates using this technique are shown in Figs. B.1 and B.1 in Appendix B as a cloud of light blue circles, the blue box is centered on the mean, and its size corresponds to the 1σ standard deviation of the points in the cloud. Table 2 summarizes the rotation rates for all 13 targets in this study.

7. Discussion

Fig. 10 Averaging kernels for both SOLA inversions at r/R∗ = 10-3, r/R∗ = 1 and a two-zone model with xc = 0.1 (from the best model of KIC 007619745). The left panel shows the central region, and the right panel shows the surface region.

The rotation rates for most of the targets show good agreement, while a few present some scatter in the envelope averages. The rotation rates agree with each other within 2σ except for one case (see further below).

There are minor differences in the ideal case of many exactly measured splittings and exact seismic models, as described in Sect. 4. The differences in this case are attributed to the different nature of the averages as computed from each method, for example, the averages as defined by Eqs. (4)and (5)are not exactly the same as the averages from Eq. (18), even if we had many kernels K(r) at our disposal to have very well localized averaging kernels .

However, considering the scatter of the predicted averages, we can still constrain the rotation rates in our targets with an error of about 0.05 μHz except in a few cases. No more than two averages that are well resolved spatially could be obtained using inversions and the seismic models, which was our hope given the previous success reported by Di Mauro et al. (2016).

We have chosen one of the stars in the sample (KIC 007619745) as an illustrative case to explain why only two values could be obtained. In Fig. 10 we show the averaging kernels for the two-zone model (inner and outer zones) and SOLA (at r/R_∗ = 10^-3 and r/R_∗ = 1). According to Eq. (18), these are simply the weight functions involved in the average. They are localized (at least enough to resolve the core and the envelope), with the inner-zone kernel and the r/R_∗ = 10^-3 SOLA kernel concentrating in approximately the same region in the stellar core. On the envelope, the outer-zone kernel and the r/R_∗ = 1 SOLA kernel coincide roughly as well, while they oscillate rapidly around zero as the radius is varied in the core regions. These rapid oscillations around zero do not present a problem since we are only interested in the integrated averages. Indeed, as Fig. 11 shows, the cumulative integrals of the averaging kernels are fortunately insensitive to the rapid oscillations. The outer-zone and the r/R_∗ = 1 SOLA cumulative kernels start growing only at around r/R_∗ = 0.2, where the inner-zone and the r/R_∗ = 10^-3 cumulative kernels have already reached unity. These properties allow us to estimate rotation rate averages of the core separately from the envelope. The inner zone averaging kernel in the SOLA inversions is basically the same when we use other radial locations well inside the radiative core, that is to say, we obtain the same average. As we approach the hydrogen burning region near r/R_∗ = 0.01, however, there is significant contamination from the outer regions of the star, as is shown in Fig. 11. All of our targets exhibit essentially identical behavior.

We note, however, that although the two-zone and SOLA kernels coincide for most targets, there were some instances where they differ slightly. In Fig. 12 we show an example illustrating this point. The two-zone kernels exhibit some leakage that leads to appreciable differences in the rotation averages when compared with the SOLA kernels, for instance. This is the case for targets KIC 011913545 and KIC 012008916.

Fig. 11 Cumulative integrals of the averaging kernels (from the best model of KIC 007619745). The integrals corresponding to the inner two-zone model and SOLA at r/R∗ = 10-3 are almost identical, the same holds for the outer two-zone and SOLA at r/R∗ = 1. A SOLA kernel at r/R∗ = 0.01 shows already significant leakage from other regions of the star, as shown by its cummulative integral not reaching unity anywhere close to that location.

Fig. 12 Cumulative integral of averaging kernels from SOLA inversions at r/R∗ = 1 (blue curves) and the cumulative integral from the outer two-zone model kernels (red curves). The left panel shows the case of KIC 008475025, where both cumulative integrals mostly coincide, are close to zero, and lead to similar average rotation rates of the envelope without appreciable leakage from central regions. In the case of KIC 011913545 (right panel) we see a difference between the cumulative integrals, particularly, the red curve is consistently below zero in the core region, which leads to leakage from these regions affecting the envelope averages.

The larger the number of splittings observed, the better the localization properties of the kernels will be and the error on the predicted rotation rates will be lower, see Fig. 13. We note that the predicted errors as computed from Eq. (19)do not include any systematic errors arising from unmodeled physics, which means that the predicted rotation rates may be precise, but not accurate. This certainly constitutes a source for discrepancies in addition to the variability induced by poorly localized averaging kernels, as discussed above.

Accidental mode misidentification can potentially lead to an inadequate seismic model. As noted earlier, we are considering only modes with a detection probability higher than 0.99 (Corsaro et al. 2015). In principle, we therefore do not expect complications arising from this issue. Since a misidentified dipole mode might reveal itself as an outlier if its rotational splitting differs considerably from its prediction, we closely examined the rotational splittings of the star KIC 012008916 and compared them with the predicted two-zone model splittings. However, no particular mode stands out in a clear way. Still, we could prune out the modes whose splittings differ the most from the two-zone prediction. This exercise modifies the predicted rotation averages and removes the scatter bringing the predictions across the six methods in good agreement. This is still not satisfactory, however, since we have no means to a priori justify the adequacy of a two-zone rotation profile, or of any other particular profile, for that matter.

To conclude our discussion, we would like to point out that the target KIC 12008916 is coincidentally one of the few in the sample that has a very high inclination angle, that is, no central peaks were detected for it. This could be connected to the scatter in the rotation rates in some way.

Fig. 13 Error on the rotation averages using a synthetic rotation profile. Starting with the seismic model of KIC 012008916, we selected all the modes that had a trapping parameter ζ smaller than a particular cut-off. The total number of modes in the model is 45 (cutoff ζ = 1). We computed two-zone inversions using synthetic splittings from modes with trapping parameters smaller than a given ζ cutoff. More splittings used certainly improve the error on the predicted rotation averages, but this does not take into account any systematic errors that may be present.

8. Summary and conclusion

Building upon the mode-fitting and identification work of Corsaro et al. (2015), who analyzed the power spectra of 19 red giant targets recorded by the Kepler space telescope, we have estimated the average core and envelope rotation rate of 13 of these targets that showed clear rotationally split mixed modes. We employed the model-independent method developed originally by Goupil et al. (2013) and improved later by Deheuvels et al. (2015) and Mosser et al. (2015). This model-independent method aims to provide two rotational averages, one for the g-mode cavity and another for the p-mode cavity. To investigate the possibility of obtaining more detailed information about the rotation profile, we have used the traditional approach of searching for optimal seismic models with the aid of the MESA stellar evolution suite and the GYRE oscillation code.

We used a total of six different methods to compute the averages, four of them based on the optimal seismic models: SOLA inversions, two-zone inversions, Bayesian inference on two-zone models, and linear fits of the rotational splittings δ as functions of the trapping parameter ζ_mod. We also used the model-independent method as implemented by Deheuvels et al. (2015) and a variation of it based on the result by Mosser et al. (2015).

Before we applied these methods to our sample of red giants, we took a particular seismic model (that of KIC 007619745) as a “true” reference model, and in conjunction with six synthetic rotation profiles, we proceeded to compute the associated rotational splittings. Using these as inputs, we compared the predictions from all the methods against the “true” synthetic rotation profiles. We found good agreement in general, with some small differences that can be attributed to the slightly different nature of the averages produced by each method.

All the rotational kernels of our 13 targets allow the computation of two distinct rotation rates. We used the functional form of the rotation profiles used in the study by Klion & Quataert (2016) to test whether we could determine the approximate location in the star where the differential rotation takes place. The results were negative given the magnitude of the observational errors. The information contained in the splittings is insufficient to localize the region with differential rotation. This is perhaps closely related to the fact that no good averaging kernels could be found for intermediate regions of any of the targets.

The averages that we could obtain, however, are in good agreement with the average rotation rates of the g- and p-mode cavities as computed from the model-independent method for most of the targets, while a few targets present more scatter in the predicted averages across the different methods. The rotation rates agree with each other within 2σ with only one exception, that of the star KIC 012008916. We identify the poorly localized two-zone averaging kernels as a contributor to the discrepancy. However, many sources of systematic errors exist that cannot be ruled out. There are indeed many simplifying assumptions involved in the stellar evolution code that directly affect the mode kernels, which means that some seismic models may not be adequate. In this respect, the model-independent methods are more reliable since they are essentially free of such systematic errors.

At any rate, the results for the rotation of the radiative cores show better agreement across all methods than the rotation of the envelopes, which is related to the fact that for p-dominated mixed modes the trapping parameter is ζ ≈ 0.5, far from 0. The cores in this target selection are spinning about 5 to 10 times faster than their envelopes, which is consistent with previous studies and still calls for a better understanding of the angular momentum redistribution in stars in the red giant branch.

Acknowledgments

The authors express their gratitude to the referee, whose comments greatly improved the manuscript. S.T. would like to thank Ehsan Moravveji for enlightening discussions and help with the MESA code. S.T. received partial funding from the ERC Advanced Grant ROTANUT (No. 670874). E.C. acknowledges funding by the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement No. 312844 (SPACEINN). R.A.G. acknowledges the support from the CNES GOLF and PLATO grants at CEA as well as the ANR (Agence Nationale de la Recherche, France) program IDEE (n ANR-12-BS05-0008) “Interaction Des Étoiles et des Exoplanètes”. We are grateful to the Kepler team and everybody who has contributed to making this mission possible. Funding for the Kepler Mission was provided by NASA’s Science Mission Directorate. The research leading to these results has received funding from the Fund for Scientific Research of Flanders (FWO) under project O6260.

References

Appendix A: Hare-and-hounds tests

Fig. A.1 Échelle diagram comparing the noise-added true model frequencies (hollow symbols) and the best model found by a grid + simplex search approach (filled symbols). Red squares correspond to l = 0 modes, green circles to l = 1, blue triangles to l = 2, and inverted black triangles to l = 3.

The adequacy of the method we used to obtain the seismic models can be assessed by a hare-and-hounds exercise. First, we selected the best seismic model of KIC 007619745 as our “true” model, then we assumed a known internal rotation profile and derived “observations” by adding some noise to the “true” mode frequencies. The errors are assumed to be the same as in the actual observations. Then we used our grid + downhill simplex search strategy to obtain a seismic model that best reproduces these “observations”.

We have deliberately used a different mixing length parameter (α_MLT = 1.7) and a different atmospheric correction power b during the search (b = 4.9) compared to the true model (which has α_MLT = 1.9, b = 4.81), thus further guaranteeing that the best model found is not identical to the true model. The best model properly recovers both the initial metallicity [Fe/H] and the initial helium content Y. The inferences on the internal rotation averages are also in good agreement with the true profile within 2σ. However, the same cannot be said of other stellar parameters such as mass, f_ov, or radius.

An échelle diagram comparing the frequencies of the best model with the noise-added true model is shown in Fig. A.1. Table A.1 lists other stellar parameters from the two models.The results for the internal rotation rates using different methods based on this best model are displayed in Fig. B.1 (bottom, right panel). Here we used a “true” rotation profile according to Eq. (24)with Ω_c = 0.7 μHz, Ω_m = 0.4 μHz, and Ω_e = 0.1μHz. Compare also with the top and middle panel in Fig. 7.

Appendix B: Internal rotation rates

We present here the rotation rates taken from Table 2 in graphical form for easier comparison.

Fig. B.1 Rotation rates for all targets. Rotation rate of the envelopes (Ωp) on the ordinates, rotation rate of the cores (Ωg) on the abscissae. We show the results inferred from two-zone inversions (black), from SOLA inversions (green), from linear fits of δ vs. ζmod (red), from linear fits of δ vs. ζas (magenta), and from Bayesian inference (gray), all with their corresponding 1σ error margins. The cloud of light blue circles is the result of applying linear fits of δ vs. δP/ ΔΠ1 sampling randomly from the splittings δ and the observed mode frequencies ν, assuming a normal distribution for both. The blue error box indicates the 1σ standard deviation of the resulting distribution on each axis.

Fig. B.1 continued.

All Tables

All Figures

	Fig. 1 Rotational kernels of two dipole modes with different values of the trapping parameter ζ (from the best model of KIC 007619745). The left panel shows the central region, and the right panel shows the region near the surface.
In the text

Fig. 2 Trapping parameter as described by the true model, ζmod (Eq. (2)) as a function of the true frequencies νmod (dashed black line and black circles). Red circles represent the trapping parameter ζas as approximated by Eqs. (11)and (8)and using νas for the mode frequencies as determined by Eq. (6). The scaled period spacing ΔP/ ΔΠ1 between any two consecutive true frequencies νmod,1 and νmod,2 is plotted at the abscissa ${\mathit{\nu }}_{\mathrm{\Delta }\mathit{P}}\mathrm{=}\mathrm{2}\mathit{/}\mathrm{(}{\mathit{\nu }}_{\mathrm{mod}\mathit{,}\mathrm{1}}^{-1}\mathrm{+}{\mathit{\nu }}_{\mathrm{mod}\mathit{,}\mathrm{2}}^{-1}\mathrm{)}$ (blue crosses). ΔPν represent the interpolated value of the two adjacent period spacings around each frequency νmod (green crosses).

In the text

	Fig. 3 Illustrative use of a Bayesian fit to compute rotation rates for KIC 007619745. The trapping parameters ζmod are obtained from the eigenmodes as computed by GYRE.
In the text

Fig. 4 Propagation diagram of the best model of KIC 007619745. The Brunt-Väisälä frequency N(r) is plotted as a solid blue curve, and the l = 1 Lamb frequency S1(r) is plotted as a dashed green curve. p-modes with frequency ν are evanescent in the region where ν < S1(r). Similarly, g-modes are evanescent wherever ν >N(r). The core-envelope boundary in our two-zone models is located at rc ≈ 0.1R∗, which usually lies in the evanescent zone common to both p- and g-modes. The peak of N(r) near r/R∗ = 0.11 is associated with a chemically inhomogeneous layer in the radiative zone, left behind by convection as the boundary moves upward.

In the text

	Fig. 5 Synthetic rotation rate profiles used for testing. Ωm refers to the rotation rate between 1.5 rH and the base convective zone, see Eq. (24). s refers to the width of a Gaussian profile as defined in Eq. (26).
In the text

	Fig. 6 Rotational splittings using the best seismic model of KIC 007619745 for one of the test profiles as defined in Eq. (24)with Ωc = 0.7 μHz, Ωm = 0.4 μHz, and Ωe = 0.1 μHz. The splittings δΔP (dark blue crosses) are the interpolated splittings as explained in Sect. 2.1. Dotted lines are linear fits providing estimates of the rotation rate averages in the g and p cavities following Eq. (3).
In the text

Fig. 7 Comparison of the predicted rotation rate averages in the core (Ωg) and in the envelope (Ωp) using the methods described in Sects. 2 and 3 based on the best seismic model for KIC 007619745. The top row corresponds to rotation profiles following Eq. (24), and the bottom row corresponds to Gaussian rotation profiles according to Eq. (26). The black dashed lines mark the true averaged rotation rates as computed from Eqs. (4)and (5). Red represents the rotation rates as determined from linear fits of (ζmod,δ), and magenta represents the rotation rates from linear fits of (ζas,δ). Blue corresponds to linear fits of δ and δΔP vs. ΔPν/ ΔΠ1 and ΔP/ ΔΠ1. Black are the predictions from two-zone inversions, and green are the predictions from SOLA inversions. Results from Bayesian inference are not included here as they are essentially identical to the two-zone models.

In the text

Fig. 8 Histogram associated with the predictions for (blue, green, and red tones, respectively). The true (test) profile is described by Eq. (24)with fixed Ωc = 0.7 μHz, Ωe = 0.1 μHz. The light colors represent the histogram obtained when Ωm = 0.1 μHz, medium colors correspond to Ωm = 0.4 μHz, and dark colors represent Ωm = 0.7 μHz. Essentially identical distributions for result, regardless of the choice of Ωm. While Ωc and Ωe are recovered properly, no good prediction for Ωm can be made.

In the text

Fig. 9 Left column plots show the observed rotational splittings, normalized by their maximum value (red error bars) and the observed period spacings, normalized by the observed period spacing ΔΠ1 (blue error bars) as functions of mode frequency for the target KIC 007619745 (top), and for the hare-and-hounds exercise. Right column plots show the observed rotational splitting δ as a function of ζmod (red points and error bars), δ vs. ζas (dark green points and error bars), and both δ vs. ΔPν/ ΔΠ1 and δΔP vs. ΔP/ ΔΠ1 (blue points and cyan error bars). See Sect. 2.1.

In the text

	Fig. 10 Averaging kernels for both SOLA inversions at r/R∗ = 10-3, r/R∗ = 1 and a two-zone model with xc = 0.1 (from the best model of KIC 007619745). The left panel shows the central region, and the right panel shows the surface region.
In the text

Fig. 11 Cumulative integrals of the averaging kernels (from the best model of KIC 007619745). The integrals corresponding to the inner two-zone model and SOLA at r/R∗ = 10-3 are almost identical, the same holds for the outer two-zone and SOLA at r/R∗ = 1. A SOLA kernel at r/R∗ = 0.01 shows already significant leakage from other regions of the star, as shown by its cummulative integral not reaching unity anywhere close to that location.

In the text

Fig. 12 Cumulative integral of averaging kernels from SOLA inversions at r/R∗ = 1 (blue curves) and the cumulative integral from the outer two-zone model kernels (red curves). The left panel shows the case of KIC 008475025, where both cumulative integrals mostly coincide, are close to zero, and lead to similar average rotation rates of the envelope without appreciable leakage from central regions. In the case of KIC 011913545 (right panel) we see a difference between the cumulative integrals, particularly, the red curve is consistently below zero in the core region, which leads to leakage from these regions affecting the envelope averages.

In the text

Fig. 13 Error on the rotation averages using a synthetic rotation profile. Starting with the seismic model of KIC 012008916, we selected all the modes that had a trapping parameter ζ smaller than a particular cut-off. The total number of modes in the model is 45 (cutoff ζ = 1). We computed two-zone inversions using synthetic splittings from modes with trapping parameters smaller than a given ζ cutoff. More splittings used certainly improve the error on the predicted rotation averages, but this does not take into account any systematic errors that may be present.

In the text

	Fig. A.1 Échelle diagram comparing the noise-added true model frequencies (hollow symbols) and the best model found by a grid + simplex search approach (filled symbols). Red squares correspond to l = 0 modes, green circles to l = 1, blue triangles to l = 2, and inverted black triangles to l = 3.
In the text

Fig. B.1 Rotation rates for all targets. Rotation rate of the envelopes (Ωp) on the ordinates, rotation rate of the cores (Ωg) on the abscissae. We show the results inferred from two-zone inversions (black), from SOLA inversions (green), from linear fits of δ vs. ζmod (red), from linear fits of δ vs. ζas (magenta), and from Bayesian inference (gray), all with their corresponding 1σ error margins. The cloud of light blue circles is the result of applying linear fits of δ vs. δP/ ΔΠ1 sampling randomly from the splittings δ and the observed mode frequencies ν, assuming a normal distribution for both. The blue error box indicates the 1σ standard deviation of the resulting distribution on each axis.

In the text

	Fig. B.1 continued.
In the text