Internal rotation of 13 low-mass low-luminosity red giants in the Kepler field
Santiago Andrés Triana Enrico Corsaro Joris De Ridder Alfio Bonanno Fernando Pérez Hernández Rafael García
 





Triana, Corsaro, et al.
Internal rotation of 13 red giants


                Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium
                
                Instituto de Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Spain
                
                Laboratoire AIM, CEA/DRF-CNRS, Université Paris 7 Diderot, IRFU/SAp, Centre de Saclay, 91191 Gif-sur-Yvette, France
                
                Departamento de Astrofísica, Universidad de La Laguna, E-38205 La Laguna, Tenerife, Spain
                
                INAF, Osservatorio Astrofisico di Catania, Via S.Sofia 78, 95123 Catania,Italy
                
                Royal Observatory of Belgium, Ringlaan 3, Brussels, Belgium
                trianas@oma.be 
          

   



 
  
  
   The Kepler space telescope has provided time series of red giants of such unprecedented
 quality that a detailed asteroseismic analysis becomes possible. For a limited set of about
 a dozen red giants, the observed oscillation frequencies obtained by peak-bagging together
 with the most recent pulsation codes allowed us to reliably determine the core/envelope
 rotation ratio. The results so far show that the current models are unable to reproduce
 the rotation ratios, predicting higher values than what is observed and thus indicating
 that an efficient angular momentum transport mechanism should be at work. Here we provide
 an asteroseismic analysis of a sample of 13 low-luminosity low-mass red giant stars 
observed by Kepler during its first nominal mission. These targets form a subsample 
of the 19 red giants studied previously <cit.>, which not only have
 a large number of extracted oscillation frequencies, but also unambiguous mode identifications.
  
   We aim to extend the sample of red giants for which internal rotation ratios obtained
 by theoretical modeling of peak-bagged frequencies are available. We also derive the rotation ratios
 using different methods, and compare the results of these methods with each other.
  
   We built seismic models using a grid search combined with a Nelder-Mead simplex
    algorithm and obtained rotation averages employing Bayesian inference and inversion methods. We compared these
    averages with those obtained using a previously developed
model-independent method.
  
   We find that the cores of the red giants in this sample are rotating 5 to 10 times faster than their
    envelopes, which is consistent with earlier results.
   The rotation rates computed from the different methods show good agreement for some targets, while
   some discrepancies exist for others.
  
   

   

   Internal rotation of 13 low-mass low-luminosity red giants in the Kepler field
    S. A. Triana <ref>,<ref>
          E. Corsaro <ref>,<ref>,<ref>
          J. De Ridder <ref>
          A. Bonanno <ref>
          F. Pérez Hernández<ref>,<ref>
          R. A. García <ref>
          
    Received ; accepted 
===========================================================================================================================================================================================================






§ INTRODUCTION


The impact of the Kepler space mission <cit.> 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 <cit.>. 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
<cit.>. 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 <cit.>. 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 <cit.>


    δ_nlm=m β_nl∫_0^R_* K_nl(r) Ω(r) dr,


where the kernel K_nl(r) and β_nl are functions of the
vertical and horizontal material displacement eigefunctions ξ_r(r) and ξ_h(r) 
(see Section 3 for more details).
Thus, rotation lifts the azimuthal wave number degeneracy of the modes.


The milestone work of <cit.> 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, <cit.> 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,
 <cit.> determined a core/envelope rotation rate ratio of about
five for a star in the lower giant branch observed by Kepler, 
<cit.> computed the rotation rate ratio of six subgiants and young
red giants, <cit.> obtained rotations rates for seven
red giants in the secondary clump, and very recently, the work by <cit.> 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. <cit.> 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 <cit.>, which
provides a suitable explanation for `anomalous' rotation rates in other types of stars,
 such as those reported by <cit.>, and <cit.>.
However, <cit.> 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 <cit.>. 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 <cit.> 
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. <ref>) to
determine the approximate rotation rates in different regions of the
star. This is the approach taken by <cit.>,
<cit.>, and <cit.>.

<cit.> 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 Section <ref>). The method relies solely on observed quantities  as inputs and particularly, on an estimate
of the asymptotic period spacing ΔΠ_1 that pure high-order g modes would have in the asymptotic
regime <cit.>.

This latter approach
was taken by <cit.> 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 <cit.>, 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.

<cit.> 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 <cit.> is 
exactly the same expression as was found by <cit.> to represent the ratio
Δ P/ΔΠ_1.




 The model-independent method used by <cit.> has been compared with inversions based
on seismic models for only two targets: the core helium-burning red giant KIC7581399 <cit.>,
  and the early red giant KIC4448777 <cit.>. Although a wider comparison
of the method using targets in different evolutionary stages is desirable, we offer
comparisons of the  model-independent method against inversions for our 13 targets, which share similar
evolutionary stages  with KIC4448777. Our targets
are a subset of the original selection of 19 low-mass, low-luminosity 
red giants studied by <cit.> using
the Bayesian inference technique Diamonds
 <cit.> to detect and identify pulsation frequencies. 
<cit.> 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 (Section <ref>) and
employ Bayesian inference and inversion techniques 
(Section <ref>) to obtain average rotation rates. We compare these results 
with those obtained by  the model-independent method as 
implemented by <cit.> and by <cit.> using their expressions
for the trapping parameter ζ (Section <ref>).
Additionally, we use the idea proposed recently by <cit.> 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.






§ 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 <cit.> 
to the total mode inertia I:



    ζ≡I_g/I = 
    ∫_r_1^r_2ρ r^2 [ξ_r^2 +l(l+1) ξ_h^2 ] dr /∫_0^R_*ρ r^2 [ξ_r^2 +l(l+1) ξ_h^2 ] dr .






In the expression above, ρ is the density, r_1 and r_2 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 <cit.>.

Figure <ref> shows the rotational kernels K_nl of two mixed modes with
different ζ using a seismic model for KIC007619745. 
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 <cit.>, which we also refer to for further details. The coefficients of this linear
relationship are related directly to the rotational averages  across the core and the envelope through  


    δ=(Ω_g/2-Ω_p)ζ+Ω_p,


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 <cit.>, they are defined as


    Ω_g=∫_0^r_2K(r) Ω(r) d r/∫_0^r_2K(r) dr


for the core, and


    Ω_p=∫_r_2^R_* K(r) Ω(r) d r/∫_r_2^R_* K(r) dr


for the envelope, where r_2 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 Section <ref> 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.







 §.§ Estimation of the trapping parameter ζ



The expression for ζ given by Eq. (<ref>)  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, <cit.> used an asymptotic analysis method based on the work of <cit.>
 to show that ζ can
in principle be estimated using observational data alone. The expression for ζ was later refined by
<cit.> and <cit.>. In what follows, we briefly recall the method and the main formulae,
we refer to the original works for further details.




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 <cit.>


    tanθ_p=q tanθ_g,


where q is the coupling constant between the p- and  g-mode cavities, and the phases θ_p,g are defined through


    θ_g=∫_r_1^r_2k_r dr,  θ_p=∫_r_3^r_4k_r dr,


where r_3 and r_4 are the inner and outer turning points of the p-mode cavity.
According to <cit.>, asymptotic analysis yields


    θ_p=π/Δν(ν-ν_p),  θ_g=π( 1/ΔΠ_1 ν-ϵ_g ),


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
 

    ν_p=ν_n,0+(1/2-d_01)Δν.


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 <cit.>: 


    ν_n,0=[n+ϵ_p+α/2(n-n_max)^2]Δν,


where n_max≡ν_max/Δν-ϵ_p. The approximate
 expression for the trapping parameter, denoted here as ζ_as , reads


    ζ_as=
    [ 1 + ν^2 ΔΠ_1/q Δνcos^2θ_g/cos^2θ_p]^-1.


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 ΔΠ_1, α, and
d_01. In practice, the optimal parameters ϵ_p, α, and Δν are determined first
by fitting the observed radial modes to Eq. <ref>. Then, the optimal parameters q, ΔΠ_1, ϵ_g,
and d_01 that best reproduce the observed l=1 mode frequencies can be found 
by a downhill simplex method. This requires solving Eq. <ref> for the mode frequencies ν at each search step
with a particular (q, ΔΠ_1, ϵ_g, d_01) 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. <ref>.

<cit.> obtained a result in their search of an expression
for the mixed-mode relative period spacings 
Δ P/ΔΠ_1  that exactly matched the expression
for ζ_as (Eq. <ref>) found by <cit.>. Thus, we have


    ζ_as=Δ P/ΔΠ_1.


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, Equation <ref> 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 ΔΠ_1. Some care must be taken because the period spacing
between two consecutive mixed dipole modes Δ P(n,n+1)=ν^-1_n+1-ν^-1_n is defined properly at 
ν=2/(ν^-1_n+1+ν^-1_n)≡ν_Δ 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. <ref>.


In Fig. <ref> 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. <ref>, or by the simpler
expression given by Eq. <ref>.

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. <ref>. 
By repeating these steps many times, we can obtain the 
distributions associated with Ω_g and 
Ω_p  and their associated errors.



 §.§ Bayesian inference


If a seismic model providing the oscillation eigenfunctions ξ_r,h is available,
 the trapping parameter can be computed from Eq. <ref>, 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. <ref> as the model. To accomplish this,
 we first compute a Gaussian log-likelihood function defined as <cit.>


    Λ( Ω_g, Ω_p) = Λ_0 - 1/2∑^N_i=1[ Δ_i(Ω_g ,
     Ω_p) /ϵ_i]^2   ,


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

    Λ_0 = - ∑^N_i=1ln√(2 π) ϵ_i  ,

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 <cit.>. 
      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 Figure <ref> for the star KIC007619745 (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.































§ 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 Section <ref>. Our treatment is based on <cit.>, 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. (<ref>). Explicitely, the kernels are computed from
the material displacement eigenfunctions via 


    K_nl=1/I[ ξ_r^2+l(l+1) ξ_h^2 -2ξ_rξ_h-ξ_h^2]r^2ρ,


where I is the mode inertia (see Eq. <ref>). The constant 
β_nl is given by


    β_nl=1/I∫_0^R_*[ ξ_r^2+l(l+1)ξ_h^2 -2ξ_rξ_h-ξ_h^2]r^2ρ dr.



The inverse problem consists of determining the unknown inversion
coefficients c_i(r) satisfying


    Ω̅(r)=∑_i=1^M c_i(r) δ_i/m β_i,


where Ω̅(r) 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. <ref>, 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 Ω̅(r) can be expressed in terms of the true profile Ω(r) by means of the 
averaging kernels 𝒦(r',r), which are related to the kernels
K_i(r) through 𝒦(r',r)=∑_i=1^M c_i(r') K_i(r) and fulfill


    Ω̅(r')=∫_0^R_*𝒦(r',r) Ω(r) dr.


The averaging kernels 𝒦(r',r) 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., <cit.> or <cit.>), so that
the variance of the predicted rotation rates can be estimated as


    σ^2[Ω̅(r)]=∑_i=1^M c_i^2(r) (ϵ_i/β_i)^2.


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.



 §.§ 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 Figure <ref>.
We can determine the inversion coefficients c_i associated with each
zone by finding the optimal Ω_g and Ω_p that minimize


    χ^2=∑_i=1^M (δ̅_i- δ_i/ϵ_i)^2,


where ϵ_i are the observation errors and δ̅_i are the predicted
splittings associated with the  two-zone rotation profile composed of Ω_g and Ω_p. 
 The averages are determined by enforcing ∂(χ^2)/∂Ω_g,p=0 after substituting
χ^2 using Eqs. <ref> and <ref>. 
Results from this method, with kernels from the seismic models discussed below, are presented on
Table 1 under the `Two-zone' heading.







 §.§ Subtractive optimally localized averaging

One of the differences between
the subtractive optimally localized averaging (SOLA) method <cit.>
and the method described above is that with SOLA we do not minimize χ^2 , but instead,
the method chooses the optimal linear combination of the
 inversion coefficients c_i such
that the averaging kernels 𝒦(r',r) resemble a given
target function
T(r',r) as closely as possible while keeping the variance
σ^2(Ω̅) low. Thus, we minimize


    ∫_0^R_*[ 𝒦(r',r)-T(r',r) ]^2 dr + μ∑_i=1^M c_i^2(r') ϵ_i^2

 
at each r', with the additional constraint ∫_0^R_*𝒦(r',r) dr=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:


    T(r',r)=N e^-(r'-r/s)^2,


N being a normalization factor.
In addition to the free parameter μ in Eq. <ref>, 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':


    ∑_k=1^M W_ik  c_k(r')=∫_0^R_* K_i(r)  T(r',r) dr,


where W_ik=∫_0^R_* K_i(r)  K_k(r) dr + μ δ_ik ϵ_i^2, together with the constraint ∑_k c_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. <ref>. We note that the observed splittings are
not involved in determining the c_i(r'), they determine the predicted rotation rate Ω̅(r') via Eq. <ref>. 





§ 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 KIC007619745, 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. <ref>. We also compute the `true'
rotational averages for both the core and the envelope for each case using Eqs. <ref> and <ref>.

For the first set of rotation profiles we adopted the functional form used by <cit.>:
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

    Ω(r)=
    Ω_c     r ≤ 1.5  r_H,
    Ω_m     1.5 r_H < r ≤ r_rcb,
    Ω_e(R_*/r)^α    r > r_rcb,


where 


    α=log(Ω_m/Ω_e)/log(R_*/r_rcb).


 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:
 

    Ω(r)=A e^-1/2(r/R_*/s)^2+B,


we set A=0.6 μHz, B=0.1 μHz and three different widths s: 0.01, 0.05 and 0.2. See Fig. <ref>.







We compute the rotational splittings of the six test profiles using Eq. <ref> 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, ΔΠ_1, ϵ_g, d_01) that best reproduces these model frequencies in order to
estimate ζ_as. As explained in Section <ref>, 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. <ref>.






Now we obtain the estimates of the rotation rates using each of the methods explained earlier.
The true averages are obtained following Eqs. <ref> and <ref>. 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. <ref>. In the case of the
rotation rates estimated via linear fits involving ζ_mod, as (see Fig. <ref>),
the resulting uncertainties are computed through a Monte Carlo simulation sampling
randomly and repeatedly for many times (10^5)
 from the normal distributions associated with each spliting δ. For the fits
of δ, δ_Δ P vs Δ P_ν/ΔΠ_1, Δ P/ΔΠ_1
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. <ref>. 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 𝒦(r) 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, <cit.> 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. <ref>. 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.

<cit.> were able to resolve three rotation rates in three distinct radial locations of KIC4448777, 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. <ref>). The reason for this
is not evident a priori, and although it deserves special attention, it is beyond the scope of the present study.




We can use the idea by <cit.> 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 
(Ω̅_c,Ω̅_m,Ω̅_e) 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 (Ω̅_c,Ω̅_m,Ω̅_e). If enough
information is indeed contained in the set of input splittings to constrain Ω̅_m, then this should reflect in a
sharply peaked distribution for it. Figure <ref> 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 Ω̅_c and Ω̅_e are close to the true values, but the probability density 
distribution for Ω̅_m 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 Ω̅_m
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 Ω̅_m can be recovered. 
All of the 13 targets in our sample exhibit this undesirable behavior.   







§ SEISMIC MODELING









From the original 19 young red giants studied by
<cit.>, 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 <cit.> together with
the GYRE pulsation code <cit.>. The MESA suite includes the `astero' module,
which implements a downhill simplex search method <cit.> 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 <cit.>) 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 <cit.> mixture together with OPAL opacity tables <cit.>. 
When computing mode frequencies, we peformed atmospheric corrections following
the method by <cit.> using a calibrated value of the exponent b=4.81 
as reported by <cit.>. 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 <ref>.







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 KIC007619745, 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 <ref>.   







































§ 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 <cit.> as described in Section <ref>.
Figure <ref> 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.

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. <ref>. 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 <ref> and in Figs. <ref> and <ref>.






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 Section <ref>. 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. <ref>.
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. <ref>
the slope and the intercept are not independent of each other.
The rotation rates using this technique are shown in Figs. <ref> and <ref> in the Appendix B as a cloud of small
violet 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 <ref> summarizes the rotation rates for all 13 targets obtained by all the methods described in this work.









§ DISCUSSION







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 Section <ref>. 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. <ref> and <ref>
are not exactly the same as the averages from Eq. <ref>, even if we had many kernels K(r) at our disposal to have
very well localized averaging kernels 𝒦(r',r).

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 <cit.>.

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. <ref> we show the averaging kernels 𝒦(r) for the two-zone
model (inner and outer zones) and SOLA (at r/R_*=10^-3 and r/R_*=1). According
to Eq. <ref>, 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. <ref> 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. <ref>. 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. <ref> 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. 












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.<ref>. We note that the predicted errors as computed from
Eq. <ref>  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 <cit.>. 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.

 







§ SUMMARY AND CONCLUSION




Building upon the mode-fitting and identification work of <cit.>, 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 <cit.> and improved later by 
<cit.> and <cit.>. This model-independent method aims 
to provide two rotational averages, one for the -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 <cit.> and a variation of it based on the result by <cit.>.

Before we applied these methods to our sample of red giants, we took a particular seismic model (that of KIC007619745) 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 <cit.> 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. 





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.



aa 
 








§ HARE-AND-HOUNDS TESTS


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 KIC007619745 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. <ref>. Table <ref> 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. <ref> (bottom, right panel). Here we used
a `true' rotation profile according to Eq. <ref> with
Ω_c=0.7 μHz, Ω_m=0.4 μHz, and Ω_e=0.1μHz. Compare also
with the top and middle panel in Fig. <ref>.












