© ESO, 2012

1. Introduction

Stars rotate, and this rotation has important consequences on their evolution. On the one hand, centrifugal acceleration reduces local gravity, mimicking a lower mass. On the other hand, rotation induces meridional circulation (Eddington 1925; Mestel 1953) and shear and baroclinic instabilities (Mathis et al. 2004), which contribute to the mixing of chemical elements.

The problem of transport of angular momentum inside stars has not yet been fully understood. Several mechanisms seem to be active, the most commonly invoked being diffusion by turbulent viscosity, transport by meridional circulation (Zahn 1992; Maeder & Zahn 1998), torques due to magnetic fields (Maeder & Meynet 2004; Mathis & Zahn 2005; Strugarek et al. 2011), and transport by gravity waves (Talon & Charbonnel 2005). For stars with significant convective envelopes (M_⋆ ≲ 1.4 M_⊙), the problem is complicated further by the magnetic braking by stellar winds (Kawaler 1988). Another problem concerns the initial rotation state, which can have an influence on the rotation profile during the pre-main sequence (PMS) and early main sequence (MS). The initial rotation state depends on the presence and lifetime of a circumstellar disk during the first stages of the PMS; magnetic fields seem to effectively lock the stellar surface to the disk (Shu et al. 1994). Maeder’s (2009) textbook contains an extensive description of these processes.

Several stellar evolutionary codes already include transport of angular momentum and associated chemical element mixing (Chaboyer et al. 1995; Talon et al. 1997; Palacios et al. 2003; Eggenberger et al. 2008). They have been developed first with a focus on studying consequences on the evolution of massive stars or evolved stars (see, e.g., Maeder & Meynet 2000; Sills & Pinsonneault 2000). The Sun and solar-like stars have been studied as above (e.g., Chaboyer et al. 1995; Zahn et al. 1997; Talon 1997), and one important application was the lithium surface depletion (e.g., Talon & Charbonnel 1998; Charbonnel & Talon 1999).

One way to study the internal transport and evolution of angular momentum is to obtain seismic information on the internal rotation profile of stars in the low and intermediate mass range at different stages of evolution. The SoHO satellite for the Sun (Turck-Chièze et al. 2010) and the ultra high precision photometry (UHP) asteroseismic space missions CoRoT (Baglin 2003) and Kepler (Borucki et al. 2010) offer such an opportunity. Seismic information has being obtained for PMS and MS stars (e.g., Benomar et al. 2010; Escobar et al. 2012), and now we also have seismic observational constraints on the internal rotation of red-giant stars (e.g., Beck et al. 2012; Mosser et al. 2012a; Deheuvels et al. 2012).

Our goal is to use seismic diagnostics to test the description of transport of angular momentum processes in 1D stellar models of low-mass stars, from the PMS to the red-giant branch (RGB). Rotational splittings (differences between nonaxisymmetric oscillation modes) are such diagnostics. When the rotation is slow, the relation between rotational splittings and rotation rate is linear and therefore relatively easy to interpret.

The information that rotational splittings provide on the internal rotation, however, depends on the physical nature of the stochastically excited modes, which in turn depend on the structure of the star, hence its age. In addition, some rotation gradient is expected to develop during the evolution; the core rotation rate is expected to become higher than the surface, as the core contracts and the envelope expands. If the rotation rate of some layers becomes too high, the associated distortion of the star and/or the Coriolis acceleration invalidates the linear approximation that provides the usual (linear) rotational splittings. In such cases, nonperturbative methods must be used (such as those developed in Lignières et al. 2006; Reese et al. 2006; Ouazzani et al. 2012). It is therefore necessary to determine if linear splitting are valid for an evolutionary stage in order to interpret correctly the observations.

We have thus started a series of papers devoted to establish seismically validated processes of transport of angular momentum that play an essential role in shaping the rotation profile of stars. The present paper is the first of this series. The second and the third will focus specifically on the case of red-giant stars. They will investigate the seismic diagnostics of rotation for slowly rotating red-giant stars for which linear rotational splittings are valid, and then rapidly rotating red-giant stars with nonperturbative methods.

In the present paper, we compute the rotation profiles and their evolution with time from the PMS to the RGB with an evolutionary code where rotationally induced mixing using the prescription of Zahn (1992) as refined by Maeder & Zahn (1998) has been implemented. We then follow the evolution of the linear rotational splittings calculated with the rotation profiles obtained from our evolutionary models and discuss their validity.

The paper is organized as follows: in Sect. 2 we describe the physical inputs that are implemented in our evolutionary code, with special emphasis on the transport of angular momentum and rotation induced mixing. We developed a version CESAM2K of the code CESAM (Morel 1997; Morel & Lebreton 2008) by implementing rotation-induced transport in radiative zones with careful attention to conservation and transport of angular momentum. The numerical scheme for the rotationally induced transport had to be modified compared to the general scheme used in CESAM2K, and some differences in the physical inputs also exist. For these reasons, the modified version will be referred to as the CESTAM code hereafter. This will avoid any possible confusion with results from other versions of CESAM2K used in the international community. The acronym CESAM stands for Code d’Evolution Stellaire Adaptatif et Modulaire, and the extra T in CESTAM stands for transport. A validated standard version of CESTAM will be made freely available.

The numerical implementation is described in Sect. 3. We then show comparisons with other codes for validation (Sect. 4). During the implementation of the code, we encountered some difficulties already present but not clearly explained in previous work. We chose to describe them carefully. In Sect. 5 we follow the evolution of theoretical linear rotational splittings for a 1.3 M_⊙ sequence of stellar models evolved from the PMS to the RGB. The ℓ = 1,2 modes are computed with the ADIPLS code (Christensen-Dalsgaard 2008) in a frequency range where the modes are expected to be stochastically excited. We discuss the behavior of the splittings with frequency that can be observed during the various phases of evolution. For red-giant models, we find rotational splittings that are much higher than observed. In Sect. 6 we consider possible reasons for this disagreement. Finally, conclusions and some perspectives are given in Sect. 7.

2. Physical input

The original version of the CESAM2K code (Morel 1997; Morel & Lebreton 2008) consists of a set of routines that calculates 1D quasi-hydrostatic stellar evolution including microscopic diffusion of chemical species. The solution of the quasi-static equilibrium is performed by a collocation method based on piecewise polynomial approximations projected on a B-spline basis. For models without diffusion, the evolution of the chemical composition is solved by stiffly stable schemes of orders up to four; in the convection zones mixing and evolution of chemicals are simultaneous. The solution of the diffusion equation employs the Galerkin finite elements scheme, and the mixing of chemicals in convective zones is then performed by a strong turbulent diffusion.

The code CESAM2K allows the choice of several options for the physics. The microscopic input physics is updated regularly. The opacity tables are presently the OPAL95 data (Iglesias & Rogers 1996) complemented at low temperatures by the Wichita opacity data (Ferguson et al. 2005). Several sets of opacity tables are provided that correspond to various mixtures of chemical elements, for instance the Grevesse & Noels (1993) or Asplund et al. (2005, 2009) solar mixtures, as well as α-element enhanced mixtures. We included Pothekhin’s updated conductive opacities (see e.g. Cassisi et al. 2007). Several options are possible for the equation of state, the most commonly used being the OPAL2005 EoS (Rogers et al. 1996; Rogers & Nayfonov 2002). Several networks of nuclear reactions (and corresponding updated nuclear reaction rates) were implemented, allowing the evolution of stars to be calculated from the PMS up to helium burning. The microscopic diffusion of chemical elements includes gravitational settling, thermal, and concentration diffusion terms, but no radiative accelerations. Two formalisms are available for microscopic diffusion transport (Burgers 1969; Michaud & Proffitt 1993). Also, two options are available for treating convection, the classical MLT theory of Böhm-Vitense (1958) or the Canuto et al. (1996) so-called full spectrum of turbulence. The atmospheric boundary condition is derived either from gray model atmospheres or from Kurucz ATLAS9 models (Kurucz 2005). Overshooting is an option.

Mass loss can be considered with different prescriptions. In the models presented hereafter we use the empirical mass loss rates of Reimers (1975) scaled with metallicity according to Schaller et al. (1992).

Recently, the CESAM2K code has been involved in the ESTA activities undertaken to prepare the interpretation of the seismic observations of CoRoT. Stellar models were calculated for a range of mass, chemical composition, and the evolutionary stages corresponding to CoRoT main targets and have been compared with the results of several other evolutionary codes showing a very good general agreement (Lebreton et al. 2008a,b; Montalbán et al. 2008).

CESAM2K is freely available for download, with the details in Morel & Lebreton (2008).

2.1. Transport of angular momentum

In convective zones, although there is differential rotation in latitude, the mean rotation rate at a given radius weakly depends on the radius. Therefore, it is often assumed in 1D stellar evolution codes that convective zones rotate as solid bodies (e.g. Talon et al. 1997; Meynet & Maeder 2000; Palacios et al. 2003). In extended convective zones, however, comparisons with 3D numerical simulations suggest instead a prescription of uniform specific angular momentum (e.g. Denissenkov & Tout 2000; Palacios et al. 2006). Both options can be used in CESTAM.

In a radiative zone, we used the formalism of Zahn (1992), refined in Maeder & Zahn (1998), to model the transport of angular momentum and chemical species by meridional circulation and shear-induced turbulence. In what follows, we sketch the model to make it clear exactly which equations we adopted, as several slightly different versions have been used (e.g. Zahn 1992; Talon et al. 1997; Maeder & Zahn 1998; Denissenkov & Tout 2000; Palacios et al. 2003; Mathis & Zahn 2004).

Turbulence is expected to be highly anisotropic owing to the stable stratification in radiative zones. Turbulence would then be much stronger in the horizontal than in the vertical direction. Thus, the hypothesis of “shellular rotation” can be used: as differential rotation is presumably weak along isobars, it is treated as a perturbation. All variables f can be split into a mean value over an isobar and a perturbation (as in Zahn 1992): (1)where P₂(cosϑ) is the second-order Legendre polynomial and p is the pressure. Higher order effects, included in the formalism developed by Mathis & Zahn (2004), are not considered here for the moment.

The velocity of meridional circulation can be written in a spherical coordinate system as (2)where r is the mean radius of the isobar and ϑ the colatitude. The vertical component U₂ is given in Appendix A, following Maeder & Zahn (1998). The equation of continuity in the anelastic approximation gives then the horizontal component V₂: (3)where ρ is the density.

The transport of angular momentum obeys an advection-diffusion equation: (4)where ν_V is the vertical component of the turbulent viscosity and d/dt represents the Lagrangian time derivative.

The relative horizontal variation of the density, , and the mean molecular weight, , obey where D_h is the horizontal component of turbulent diffusion discussed in Appendix B. The evolution of Λ, Eq. (6), depends on the competition between the advection of a mean molecular weight gradient and its destruction by the horizontal turbulent diffusion.

2.2. Evolution of the chemical composition

The vertical advection of chemicals due to the large-scale meridional circulation coupled with a strong horizontal turbulent diffusion results in a vertical diffusion process (see, e.g. Chaboyer & Zahn 1992). The equation of the chemical composition evolution can then be written as (7)where X_i is the abundance by mass of the ith nuclear species and D_V = ν_V and D_eff the vertical diffusivity and the diffusion coefficient associated with the meridional circulation; D_V and D_eff are discussed in Appendix B.

A necessary condition for shear instability is the Richardson criterion as given by Talon & Zahn (1997). Another condition is that the turbulent viscosity ν_V must be greater than the molecular viscosity, ν, as expressed by the Reynolds criterion (8)where Re_c ≃ 10 is the critical Reynolds number (Schatzman et al. 2000). When condition (8) is not satisfied we use D_V = ν_V = ν.

2.3. Initial conditions

Stars are fully convective when they start their PMS evolution on the Hayashi track. Assuming that convective zones rotate like solid bodies, the star should have uniform angular velocity at the beginning of its evolution.

Several facts complicate this simple picture. First, Palla & Stahler (1991) showed that stars that are more massive than about 2 M_⊙ are no longer fully convective when they appear on the PMS (the exact mass depends mainly on the protostellar accretion rate and deuterium abundance). Second, if during the PMS the only process slowing down stellar rotation were the magnetic braking by stellar winds mentioned above, stars would rotate much more rapidly than observed on the ZAMS. The PMS is too short for this process to slow down the star significantly. An additional process is needed during the PMS, most likely disk locking (e.g., Bouvier et al. 1997).

Young stars are most often surrounded by a circumstellar disk left over after the main accretion phase is over. The magnetic coupling between the star and the disk slows the star down (see, e.g., Shu et al. 1994), and this effect is often modeled by assuming that the stellar surface corotates with the disk at a constant angular velocity (see Bouvier et al. 1997) as long as the disk exists. Once the disk has disappeared, the star surface rotation evolves freely. This scenario is implemented in CESTAM with the disk lifetime τ_disk and period P_disk as free parameters of the model.

2.4. Magnetic braking

Stars that are less massive than about 1.4 M_⊙ have significant outer convective zones. A solar-type dynamo operates there and generates a magnetic field. The coupling between the magnetic field and the plasma in the stellar wind strongly brakes the rotation of the star. Kawaler (1988, hereafter K88) proposed the following law for the loss of angular momentum : (9)where Ω_sat is a saturation angular velocity, above which magnetic field generation seems to saturate. The value of Ω_sat is often set at 8–14 Ω_⊙ (as in Bouvier et al. 1997). There are indications, however, that Ω_sat varies with stellar mass (see, e.g. Krishnamurthi et al. 1997; Andronov et al. 2003).

The parameter K_W in Eq. (9) is usually calibrated by requiring that calibrated solar models have Ω = Ω_⊙ = 2.86 × 10^-6   rad   s^-1. The precise value of K_W needed to spin down the Sun to its current period depends on the prescription for the transport of angular momentum adopted (see Sect. 4.2 below). The parameter K_W should also depend on stellar mass (see, e.g. Krishnamurthi et al. 1997).

Recently, Reiners & Mohanty (2012, hereafter RM12) have criticized the approach used in K88. Specifically, K88 supposed that the surface magnetic flux goes as some power of the angular velocity, whereas RM12 suggest instead that it is the magnetic field strength that obeys such a law. As a consequence, they obtain the following law: (10)The value of K_RM does not depend on the stellar mass, and Ω_sat ≃ 3 Ω_⊙. Both prescriptions, Eqs. (9) and (10), are implemented in CESTAM and tested below.

3. Numerical procedure

3.1. The overall problem

Stellar evolution depends on two interconnected problems: the problem of stellar structure (solving the stellar structure equations), and the problem of the evolution of the chemical composition. It has proved very difficult to solve the two problems simultaneously; different stellar evolution codes employ several techniques to overcome this difficulty. Some codes (e.g., Degl’Innocenti et al. 2008) compute the solutions of the structure equations and chemical composition in an independent way, where they use the structure calculated at the previous time step to evolve the chemical composition to the current time step, and then use this composition to calculate the structure at the current time step (or the other way around: first the structure, then the chemical composition). A second kind of code computes the chemical composition between each iteration of the structure problem, as in Scuflaire et al. (2008). And finally, Eggleton’s code (Eggleton 1971) solves the two problems simultaneously. Stancliffe (2006) has shown that differences between the three kinds of codes are only significant on the AGB.

The CESAM2K code (Morel 1997; Morel & Lebreton 2008) belongs to the second category above. It begins a time step by updating the chemical composition. With the new chemical composition, a new structure is obtained by performing one iteration of the algorithm used to solve the structure equations. This structure is used to update the chemical composition again, before a new iteration on the structure is performed. The procedure is repeated until convergence of the structure algorithm. CESTAM keeps this iterative scheme.

Rotation with angular momentum transport introduces a new problem interconnected with the previous two. The evolution of the chemical composition depends on the turbulent mixing induced by differential rotation, whereas the structure equations are changed by the inclusion of the centrifugal acceleration. In CESTAM, we inserted the resolution of the angular momentum transport equations at the beginning of the cycle, because the rotation profile is required to calculate the turbulent diffusion coefficients that are needed to update the chemical composition. Our cycle, then, is as follows: we update the rotation profile, then the chemical composition, and finally we iterate on the structure. The cycle is repeated until convergence.

In the absence of external torques, total angular momentum is conserved. The rotation profile obtained at this stage does not enforce angular momentum conservation, however, because the profile Ω(m) was computed with the stellar structure obtained in the iteration before convergence. The total angular momentum is and the function r(m) changes between iterations. We need to compute the rotation profile one more time after convergence to make sure that it is consistent with the stellar structure, so that angular momentum is numerically conserved.

3.2. The rotation profile

Equation (4) (with Eqs. (A.1), (5) and (6)), is a fourth-order differential equation in r. To solve it, we split it into four first-order differential equations, and solve the system using the well known relaxation method (Henyey et al. 1964; Press et al. 2007). Equation (6) is solved simultaneously with Eq. (4), so that we have a total of five finite difference equations. We chose this method for its simplicity in dealing with the complex relations between variables expressed in Eq. (4), taking Eqs. (A.1), (5), and (6) into account. Methods that require that the solutions are approximated by a linear combination of known functions (collocation methods, spectral methods) are therefore difficult to implement. In practice, the relaxation method we employed proved to be efficient, robust, and fairly stable.

The scheme can be fully implicit or semi-implicit in time (to ensure higher order accuracy in time). The time step can be subdivided if needed, a useful feature for future developments involving faster processes.

Four boundary conditions are needed. At the top of a radiative zone, we impose conservation of angular momentum and no differential rotation: (11)where m_t is the mass inside the top of the radiative zone. If the convective zone above it is at the surface, is the torque applied at the surface of the star, Eq. (9), otherwise . At the bottom of a radiative zone, similarly, (12)where m_b is the mass inside the bottom of the radiative zone. If the center of the star is radiative (m_b = 0), Eq. (12) is replaced by U₂ = 0, the requirement that there is no mass flow out of the center.

Intermediate convective regions (regions that are neither at the center nor at the surface) are treated as if they were special points. The equations between the beginning and the end of an intermediate convective zone (between m_i and m_f) are Equation (13) guarantees conservation of total angular momentum, Eqs. (14)–(15) impose no shear at the borders of the convective zone, and Eqs. (16)–(17) result from the absence of μ-gradients in a convective zone.

The four equations resulting from Eq. (4) are written in finite-difference form between pairs of points, say, at m = m_k − 1 and m = m_k. Equation (6) requires a special treatment, because it is not a differential equation in r (or m). We write Eq. (6) as an algebraic equation at point m = m_k − 1, because otherwise we would be solving for (the average between m_k − 1 and m_k instead of the values at m_k − 1 and m_k). Indeed, Λ only appears as in the other equations, since it has no space derivatives. A fifth boundary condition is needed: Λ = 0 at m = m_t (the top of the radiative zone) as in Eqs. (16)–(17).

4. Comparison with results from other evolutionary codes

The complicated nature of the equations describing the evolution of the rotation profile (and particularly the evolution of U₂) made it important to validate our approach. To do that we compared our results with those obtained with other implementations of rotational mixing in stellar evolution codes, namely STAREVOL (Palacios et al. 2003) for the solar case (see also Turck-Chièze et al. 2010) and the Geneva stellar evolution code (Talon et al. 1997) for the case of 3 M_⊙ and 5 M_⊙ stars.

4.1. Standard physics

In these comparisons, we used the OPAL equation of state (Rogers et al. 1996) and opacities (Iglesias & Rogers 1996), complemented at T < 10⁴ K by the Alexander & Ferguson (1994) opacities. We used the NACRE nuclear reaction rates of Angulo et al. (1999) except for the reaction, for which we used the reaction rates given in Imbriani et al. (2004). The solar compositions of Grevesse & Noels (1993, hereafter GN93) and Asplund et al. (2005, hereafter AGS05) were used. The Schwarzschild criterion was used to determine convective instability. Convective core overshoot fully mixes the chemical composition to a distance from the border of the convective core (where r_co is the radius of the core determined by the Schwarzschild criterion). The temperature gradient in the overshoot zone is ∇ = ∇_ad. The centrifugal acceleration is taken into account by adding the average centrifugal acceleration 2Ω²r/3 to gravity in the hydrostatic equilibrium equation. The atmosphere is computed in the gray approximation and integrated up to an optical depth of τ = 10^-4.

4.2. The Sun

We computed calibrated rotating solar models using the compositions of GN93 and AGS05, with and without magnetic braking according to the magnetic braking law, Eq. (9). Models were calibrated to within 10^-5 in luminosity, radius, and surface metallicity (Z = 0.0195 for GN93, Z = 0.0126 for AGS05). Initial rotational velocities (and parameter K_W for the models with magnetic braking) were chosen so that models have an equatorial velocity v_eq = 2.02 km s^-1 at the solar age (4.6 Gyr).

The temperature gradient in convective zones was computed using the convection model of Canuto et al. (1996) with l = α_CGMH_P. All models include microscopic diffusion and settling using the approximations proposed by Michaud & Proffitt (1993). The initial rotational velocities for the cases with braking were chosen in order to have v_eq = 20 km s^-1 at the ZAMS.

Table 1 shows the values of parameters resulting from the solar calibration in α_CGM, the initial helium, and metallicity abundances (Y_i and Z_i). The greatest differences between the cases with and without magnetic braking concern the initial abundances. This is because models with braking rotate much faster at the center (as shown in Turck-Chièze et al. 2010) and thus have a higher Ω-gradient, leading to a much higher turbulent diffusion. Turbulent diffusion tends to homogenize the radiative zone to a large extent, partially erasing the composition gradient created by gravitational settling.

Table 2 shows some characteristics of the calibrated solar models. The higher turbulent diffusion partially stops the settling of helium, causing a higher helium abundance in the convective zone. However, the strong Ω-gradient predicted by models with braking in the radiative zone is in direct contradiction with helioseismic results, which indicate a flat rotation profile to within r ≃ 0.25 R_⊙ of the center (see, e.g. Chaplin et al. 1999; Couvidat et al. 2003; Eff-Darwich et al. 2008). A new physical mechanism is needed to explain the discrepancy, such as transport of angular momentum by internal gravity waves (see Talon & Charbonnel 2005) and/or magnetic stresses (Maeder & Meynet 2004; Mathis & Zahn 2005).

Figure 1 shows differences between the sound speed c_s profiles of the seismic Sun of Turck-Chièze et al. (2001) using results from SoHO (GOLF-MDI) and our calibrated solar models (see also Goupil et al. 2011, for intermediate opacities AGS09). The seismic Sun of Turck-Chièze et al. (2001) agrees with Basu et al. (2009) using results from BiSON and MDI. A good agreement at the level of 0.3% is seen when using the old mixture of GN93 with only small discrepancies below the convection zone and in the central region. On the other hand, severe discrepancies occur when more recent mixtures are used, such as AGS05. These results are quite similar to what is found in the literature (see, e.g., Basu et al. 1997; Basu & Antia 2008; Turck-Chièze et al. 2011,for models with no rotation induced transport).

We compared our results for calibrated solar models including rotation-induced transport of Type I with the results obtained in Turck-Chièze et al. (2010) with the STAREVOL code. We found good agreement between the profiles of the angular velocity, U₂, and the two components of turbulent viscosity. Several improvements have been implemented in the code CESTAM compared to the CESAM2K version used in Turck-Chièze et al. (2010), but with no consequence for the solar case. We confirm that rotation induced transport of Type I does not help remove the discrepancies, as found by Turck-Chièze et al. (2010).

4.3. Higher mass, main sequence stellar models

We compared results for models with higher mass stars. These stars have thin convective or fully radiative surface layers and are therefore not expected to undergo magnetic braking. Calculations here are carried out assuming global conservation of angular momentum, which is obtained at the precision level of 10^-6.

The first comparison concerns the evolution of a 5 M_⊙ model computed with CESTAM and Geneva codes (as in Talon et al. 1997), using X_i = 0.73 and Z_i = 0.01 (and no magnetic braking). These models were computed without microscopic diffusion and settling. Figure 2 shows the profile of U₂(r) at the middle of the MS (when X_c = 0.35). Curves for the same quantities superimpose, showing excellent agreement between the results obtained with the two codes. There are small differences that we attribute to the different microphysics used in the codes. The rotation profiles, not shown, also agree during the course of the evolution.

We also compared the effects of rotational mixing on the evolutionary tracks on the Herzsprung-Russel diagram (HRD) between CESTAM models and those of Eggenberger et al. (2010) computed with the Geneva code (Eggenberger et al. 2008). Figure 3 shows evolutionary tracks on the HRD for a 3 M_⊙ model calculated without extra mixing (microscopic diffusion, convective overshoot, rotation), with rotational mixing only and with overshooting only.

We considered two cases for the models with overshooting, α_ov = 0.1 and α_ov = 0.2. The model with rotational mixing has v_eq = 150 km s^-1 at the ZAMS. As in Eggenberger et al. (2010), the evolution of the rotating model closely resembles the evolution with α_ov = 0.1. The difference between the rotating models and the nonrotating models at the ZAMS is due to the centrifugal acceleration. The decrease in the effective gravity (gravity minus centrifugal force) replicates a nonrotating star with a lower mass.

The evolutionary tracks shown in Fig. 3 are identical to those in Fig. 2 of Eggenberger et al. (2010).

4.4. Evolved low-mass models

The rotation profiles that we obtain for our red-giant models are similar to those computed by Eggenberger et al. (2010). Examples for M_⋆ = 1.3 M_⊙ are given in the section below.

5. Evolution of rotational splittings from PMS to RGB

Rotational splittings are useful as seismic diagnostics for measuring the stellar internal rotational profile. For slow rotation, a first order perturbation description provides rotational splittings that are linearly dependent on the rotation profile. They are therefore convenient tools, easy to compute from stellar models and theoretical oscillation codes for comparison with observations. The first order approximation provides the following expression for the linear rotational splittings (Christensen-Dalsgaard & Berthomieu 1991, and references therein): (18)where x = r/R_⋆ is the normalized radius, Ω is the angular rotation (rad/s) and the rotational kernel K_nℓ takes the form (19)where I_nℓ is the mode inertia (20)The quantities entering the equations above are the fluid vertical and horizontal displacement eigenfunctions, ξ_r and ξ_h respectively, and the density ρ. For solid-body rotation, the rotational splittings become: (21)with (22)where C_nℓ are the Ledoux coefficients.

For asymptotic pure p-modes, C_nℓ ~ 0 and β_nℓ ~ 1. For pure g-modes, C_nℓ ~ 1/ℓ(ℓ + 1) and β_nℓ ~ 1−1/ℓ(ℓ + 1), i.e. β_nℓ ~ 1/2 for ℓ = 1 modes.

To investigate the evolution of rotational splittings with stellar age, we computed a 1.3 M_⊙ evolutionary sequence including rotationally induced transport as described in Sect. 2 above. We used the same input physics as in the case AGS05 with braking described in Sect. 4.2, but without microscopic diffusion and settling, as we found that microscopic diffusion and settling do not affect the rotation profiles of our models. We chose M_⋆ = 1.3 M_⊙ because it is a typical seismic mass obtained for the Kepler and CoRoT red-giant stars.

Figure 4 shows an evolutionary track for a sequence of 1.3 M_⊙ stellar models in a HR diagram computed with CESTAM. In order to interpret the rotational splittings, we plot the evolution of the central and surface rotation rates with age in Fig. 5.

We then computed oscillation frequencies and rotational splittings for several models spanning the track (shown in Fig. 4) using the freely available ADIPLS adiabatic oscillation code (Christensen-Dalsgaard 2008). We computed frequencies of ℓ = 0,1 axisymmetric modes and rotational splittings for ℓ = 1 modes.

The 1.3 M_⊙ models have a convective envelope from PMS to the RGB, hence we expect stochastically excited modes all along the sequence. Thus, we chose the frequency range that spans an interval of a few radial orders below and above n_max, the radial order corresponding to the frequency at maximum power. It is estimated as n_max = ν_max/Δν, where the frequency at maximum power spectrum ν_max and the mean large separation Δν are given by the usual scalings relations with T_eff, ⊙  = 5777 K, ν_max = 3050 μHz and Δν_⊙ = 134.7 μHz for the Sun. Indeed, it has been conjectured, then shown observationally, that these relations predict well the location of the excited frequency range of stochastically excited modes (e.g. Brown et al. 1991; Kjeldsen & Bedding 1995; Kallinger et al. 2010).

5.1. Validity of linear approximation for the rotation splittings

The validity of the linear approximation is estimated by comparing the rotation rate to the oscillation frequency. In a perturbative approach, the parameter (25)is thus assumed smaller than unity. We evaluate this parameter at the center of each selected equilibrium model of the evolutionary sequence, where Ω is largest, and for a range of frequency spanning the radial order interval (n_max − 4,n_max + 4).

Sharp rotation gradients develop in the central layers and can also cause some departure from the linear approximation. In the oscillation equations (see, e.g. Eqs. (34.7)–(34.12) in Unno et al. 1989), the rotation gradient term appears with a factor (26)Thus, we consider that the linear approximation is valid if (27)where the first inequality follows from assuming ζ < 1.

We then evaluate the quantity (28)at the radius wherethe rotation gradient is greatest and for ω/2π = ν_max. Figure 6 shows both quantities, ζ and Δ, as functions of the effective temperature of the models. Both parameters remain much smaller than unity until the model reaches the base of the RGB. Therefore, the use of the linear approximation for the rotational splittings can be safely used for stochastically excited solar-type modes of low mass stars except for the fast rotating cores of red giants.

5.2. Evolution of rotational splittings of stochastically excited modes along an evolutionary track

In the following, we discuss the information that can be retrieved from the average value throughout the star  ⟨ Ω/2π ⟩  = δν_nℓ/β_nℓ weighted by the mode. This definition is given at fixed ℓ, and depends on the radial order n.

• The PMS regime: on the PMS, the central and surface rotationrates are equal as long as the model is fully convective. Therotation rate remains constant in time for the duration of disklocking. When the radiative core appears, the central rotationstarts to increase due to the contraction of the central layers. Whendisk locking stops, the surface rotation rate starts increasing as thestar contracts. Figure 7 shows the rotation profilesof the selected PMS models shown in Fig. 4. Thesurface rotation rate evolves from2.84 μHz (model #1) to 3.9 μHz (model #5) while the central rotation rate increases from 11.3 μHz to 21.2 μHz. The rotational splittings for ℓ = 1 modes are computed for models #1 to #6 according to Eq. (18) and are also shown in Fig. 7. The splittings are nearly independent on radial order and increase steadily with the surface rotation rate as the model evolves along the PMS towards the ZAMS. The value of the Ledoux constant almost vanishes, and β_nℓ ~ 1, because the modes are essentially pressure modes in the frequency range where they are expected to be stochastically excited. These modes have no amplitudes in the inner layers and cannot probe the rotation there. Thus the mean rotation rate  ⟨ Ω ⟩ /2π corresponds essentially to the rotation averaged over the outer layers. As the rotation increases inwards in these layers, we find that  ⟨ Ω ⟩ /2π = 3.8 μHz for Model #1, close to – but slightly larger than – its surface rotation rate. Model #6 is already on the MS. Its surface rotation is faster than that of younger models but its central rotation is slower. This is not seen in the corresponding splittings and the averaged rotation  ⟨ Ω ⟩ /2π remains close to the surface value. Results are similar for ℓ = 2 modes (not shown).

  Fig. 8 Top panel: rotation profiles as a function of the normalized mass m/M⋆ for 1.3 M⊙ MS models #6 (magenta), #7 (blue), #9 (red), #10 (green) and #11 (black) shown in Fig. 4. Second panel: rotational splittings for ℓ = 1 modes as a function of the normalized frequency ν/Δν for models #7, #9, #10 and #11 (same colors as above). Third panel: βnℓ for the same models. Bottom panel:  ⟨ Ω ⟩  = δνnℓ/βnℓ for the same models. The large separation Δν goes from 97.6 μHz (model #7) to 63.2 μHz (model #11).

  Open with DEXTER

• The MS regime: on the MS, the surface rotation rate decreases with age due to both braking at the surface and an increase of the stellar radius. The central rotation rate also decreases with time due to transport of angular momentum from the core to the surface. The rotation in the central regions first decreases from Model #6 to #9 and increases by roughly 60% from model #9 to #11 (Fig. 8). The excited modes are still in the frequency domain of high order p-modes and their rotational splittings reflect the surface behavior only. The splittings decrease with the surface angular velocity. At the end of the MS (model #11), the central rotation rate has already reached 15.4 μHz, while the surface rotates at a rate of 0.6 μHz: the core is rotating roughly 26 times faster than the surface. As can be seen in Fig. 8, the rotational splittings for this model yield a mean rotation  ⟨ Ω ⟩ /2π = 3.8 μHz, that is 1.22 μHz higher than the surface rotation rate but much lower than the central rotation rate.

  Fig. 9 Top panel: rotation profiles (logarithmic scale) as a function of the normalized mass m/M for 1.3 M⊙ subgiant models #12 (blue solid) to #17 (black solid) as indicated in Fig. 4. Second panel: rotational splittings for l = 1 modes as a function of the normalized frequency ν/Δν for model #12, #14, #15 and #16. Third panel: βnℓ for the same models. Bottom panel:  ⟨ Ω ⟩ /2π = δνnℓ/βnℓ for the same models. The large separation Δν goes from 56.5 μHz (model #12) to 38.0 μHz (model #16).

  Open with DEXTER

• The subgiant regime: when the model cools along the subgiant branch, its rotation keeps evolving differently in the inner and outer parts, slowing down at the surface and accelerating in the inner regions (Fig. 9). Substantial changes appear when the model leaves the MS and evolves as a subgiant. The Brunt-Väisälä frequency of red-giant stars is very high in the radiative interior, and the frequencies of g-modes enter the frequency domain of p-modes in the frequency range where stochastically excited modes are expected (Christensen-Dalsgaard 2012). During the subgiant phase, avoided crossings between p- and g-modes appear. As observationally shown (Deheuvels et al. 2010; Mosser et al. 2012b), some p-modes become mixed modes. These modes have amplitudes both in the inner regions, where they behave as gravity modes, and in the surface layers, where they behave as acoustic modes (Dziembowski & Pamyatnykh 1991). They are very interesting as they allow probing regions deep inside the star. The number of such mixed modes is small at the beginning of the subgiant branch and increases as the star evolves toward the RGB. During this phase, modes will eventually undergo avoided crossings, exchanging nature from p to g. Figure 9 shows that Model #12 is not evolved enough so that no mixed modes are present in the frequency domain shown. The splittings remain nearly frequency independent and provide a mean rotation rate close to its surface value. The variation in the splittings with frequency for the more evolved models #13 to #17 clearly shows the presence of mixed modes. With amplitudes in the dense interior, mixed modes have a much stronger inertia than the neighboring p-modes (Dupret et al. 2009). For the same reason, the rotational splittings exhibit the same behavior as the mode inertia although in a more pronounced way. They vary from one mode to the next. They are also larger than for pure p-modes when the rotation increases inward within the model. For subgiant stars, only one ℓ = 1 mixed mode exists between two consecutive radial modes. This explains the saw-like aspect of the variation in the splittings with frequency. This is confirmed by the same variation in β_nℓ with frequency: p-modes have β_nℓ close to 1 while mixed modes have β_nℓ close to 1/2. Mixed modes enter the frequency domain by the low part and one can see that β_nℓ is closer to 1 at high frequency where the g nature of the mixed modes is less pronounced. Because of the weighting by the eigenfunction, the values of the splittings yield a mean rotation  ⟨ Ω ⟩ /2π that is no longer dominated by the outer layer contribution, since the mixed modes have amplitudes in the central region where the rotation rate is larger. As a result, one obtains a rotation rate that is higher than the surface value but still much lower than the rotation rate of the central regions.

• The RGB: when the outer convective region progresses inward, its uniform rotation extends deeper as well. The surface rotation rate decreases with time essentially due to an increase in the stellar radius. The core accelerates, and in the intermediate region, where the H-burning shell lies, a sharp rotation gradient develops (see Palacios et al. 2006). The frequency spectrum of red giants is composed of g-dominated mixed modes and mixed modes with nearly equal p and g character. When the star evolves up the red-giant branch, the number of g-dominated modes increases and largely outnumbers the p-g mixed modes (Dupret et al. 2009). These modes mostly probe the core rotation (Beck et al. 2012; Deheuvels et al. 2012; Mosser et al. 2012a). For model #25 (the last shown in Fig. 4), the core rotation amounts to 251 μHz and the expected excited frequencies are about 90 μHz. This indicates a ratio ζ = 5.6 (Eq. (25)). At such rotation rates, the first-order perturbation is not relevant to computing the rotational splittings. Nonperturbative methods must then be used. However, observations show that several RGB stars have core rotation rates that are much lower than predicted by our models. (Mosser et al. 2012a). Indeed the observed values amount to a few hundred nHz, whereas here we find several dozen to a few hundred μHz. This leads to the conclusion that the central regions of our models rotate too fast and that rotationally induced transport of the kind included in our models is not efficient enough to slow down the central rotation of red-giant stars.

6. Seismic tests of transport of angular momentum: slowing down the rotating core of red-giant stars

As in the solar case, one can wonder whether another mechanism, such as internal gravity waves or magnetic fields, operates more efficiently to shape the rotation profile, either at the red-giant phase or before. However, several issues must be discussed before we can reach such a conclusion.

6.1. Uncertainties on stellar modeling and transport of angular momentum

A possible reason for this discrepancy is some inaccuracy in the current physical description of stellar models. Indeed several assumptions enter the rotationally induced transport as prescribed by Zahn (1992) and Maeder & Zahn (1998). Besides, several other uncertainties on the input physics of stellar models can also affect the resulting rotation profile to some extent. In what follows, we discuss the sensitivity of the rotation profile of red-giant models to several uncertainties.

6.1.1. Effect of convective overshoot

A 1.3 M_⊙ star has a convective core during the MS. Overshoot shifts the tracks on the HR diagram, imitating tracks with higher masses. As a result, at a given radius in the RGB, the core of the model with overshoot has had less time to spin up, since it started its contraction later in the evolution. An overshoot of 0.1 H_P reduces the core rotation rate by about 32% when R_⋆ = 3.73 R_⊙.

On the other hand, we need to reduce the mass of the model with overshoot by 0.04 M_⊙ to reproduce the RGB at the same location in the HR diagram. In that case, the model with overshoot and M_⋆ = 1.26 M_⊙ has a central rotation rate (again, when R_⋆ = 3.73 R_⊙) that is similar to the central rotation rate of the model without overshoot and M_⋆ = 1.3 M_⊙.

6.1.2. Effect of the initial rotation state

It is well known that, after the early stages on the MS, the surface rotation rate does not depend on the initial conditions; it only depends on the magnetic braking law (Eqs. (9) and (10); see, e.g. Maeder 2009 and references therein). We found that it is indeed the case for our models. The internal rotation profile down to the center is also independent of the initial conditions after the early stages of the MS.

6.1.3. Impact of the magnetic braking law

We computed an evolution with the braking law of K88 with a value of K_W twice the solar calibrated value to reduce the rotation rate. We found that the surface rotation rate is slowed down, as expected, but the central rotation rate changes only slightly. For instance, at the TAMS, the surface rotation rate changes by 20%, and by 3% at the base of the RGB, while the central rotation rate changes by 2% at the TAMS, 3% at the base of the RGB.

We calculated another evolutionary sequence with the RM12 magnetic braking law. With a coefficient K_RM calibrated so as to yield the solar rotation rate in a solar model, our 1.3 M_⊙ is much more efficiently slowed down compared to the evolution computed with K88. The central rotation rate, however, changes by at most 2%.

6.1.4. Stability of the rotation profile

The changes on the surface rotation rate caused by the mechanisms described above do not change the central rotation rate significantly because the extraction of angular momentum from the core is not efficient enough. The core and the surface are not sufficiently coupled. One potential way of coupling the core to the surface is through instabilities that might arise because of rotation gradient becoming very steep (as seen in Fig. 10). The Rayleigh criterion requires that r²Ω increases outwards. For our 1.3 M_⊙ model, the rotation profile becomes unstable according to the Rayleigh criterion in two places, at the shell source and just below the convective envelope. However, the stabilizing effect of the density stratification overcomes this instability, as expressed by the Solberg-Hoiland criterion: (29)where is the Rayleigh frequency in a rotating medium, given at the equator by (30)According to the Solberg-Hoiland criterion, the rotation profiles we obtain are stable throughout the evolution. However, the Goldreich-Schubert-Fricke instability (GSF, see Goldreich & Schubert 1967; Fricke 1968) may occur in stellar conditions. Hirschi & Maeder (2010) have studied its effects on pre-supernova models and find that if the GSF instability is always present regardless of the μ− and T − gradients. We implemented their prescription in CESTAM. Where , turbulent transport by the GSF instability operates with a viscosity coefficient given by the solution of their Eq. (20), (31)where K is the thermal diffusivity given by Eq. (A.7). The viscosity coefficient associated to the GSF instability is then given by D_GSF = 2x.

We found that in two regions from the end of the MS, just below the border of the outer CZ and at the shell source. At the shell source, the Ω-gradient is high because of the contraction of the layers left behind by the movement of the shell source, while the descending CZ lowers Ω just above its border. In our models, the GSF turbulent viscosity is on the order of D_V in the unstable regions. The turbulent transport by the GSF instability reduces the Ω-gradient in the unstable region at the shell source, reducing the central rotation rate by 3% in models at the base of the RGB. The turbulent transport by the GSF instability does not have time to change the Ω-gradient below the CZ before the convective instability sweeps over this region.

6.1.5. Uncertainties on the turbulent viscosity coefficients

To extract more angular momentum, we can either increase the meridional circulation or increase the vertical turbulent viscosity. To increase the meridional circulation, one option is to decrease the inhibiting effect of Λ (see the expression for the μ-currents, Eq. (A.5)). According to Eq. (6), an increase in D_h will lead to a decrease in Λ. The prescription for the horizontal coefficient of turbulence D_h is still open to discussion. Prescriptions derived by Zahn (1992), Maeder (2003), and Mathis et al. (2004) can differ by as much as two orders of magnitude. We adopted, as a test, a value of D_h that is 10² higher than that previously used (Eq. (B.2)).

The value of the critical Richardson number in Eq. (B.4) is usually assumed to lie between 1/4 and 1/6. Brüggen & Hillebrandt (2001) and Canuto (2002) indicate, however, that we can expect instability for Ri ≃ 1. To increase the efficiency of the vertical turbulent viscosity in transporting angular momentum from the core, we computed an evolutionary sequence using Ri_c = 1 to test the effect of the uncertainties. Figure 11 shows profiles of the rotation rate computed with different assumptions of the diffusion coefficients, at the TAMS and at the base of the RGB when R_⋆ = 3.73 R_⊙. Models have the same mass and radius, hence the same Δν. As expected, with an increased value of Ri_c from 1/6 to 1, the central rotation rate decreases by 18% at the TAMS and by 30% at the base of the RGB. With a D_h increased by two orders of magnitude (and Ri_c = 1/6), the central rotation rate decreases by 84% at the TAMS and the base of the RGB. When both changes are included (an increase of Ri_c and D_h as above), the central rotation rate is decreased by an order of magnitude, while we need a decrease of two orders of magnitude to reproduce the observed rotational splittings.

6.1.6. A slowly rotating red-giant model

We eventually computed an evolutionary sequence where we combined all effects causing a decrease in the rotation rate in the central region; an overshoot of 0.1 H_P and a vertical turbulent viscosity D_V computed with Ri_c = 1, together with D_h increased by a factor 10². This model was computed with M_⋆ = 1.2 M_⊙ because during the RGB the evolutionary track on the HRD of a 1.2 M_⊙ computed with all these effects lies close to the evolutionary track of the 1.3 M_⊙ models used as comparison. The resulting rotation profile for models at the base of the RGB with same Δν is significantly decreased in the central region with the central value Ω_c = 21.5 μHz and is displayed in Fig. 12.

Calculations of linear rotational splittings shows that the maximum splittings amounts to 5–6 μHz for model #20 at the base of the RGB (Fig. 13). This corresponds to a maximum value for δν/β_nl of 12 μHz. The rotation gradient and the value of the central rotation are now too low for nonperturbative calculations to give rise to significant corrections to the values of the linear rotational splittings. Such rotation splittings lead to a core rotation rate that is closer to the observations but still too high by almost one order of magnitude compared with what is seen in recent observations (Beck et al. 2012; Deheuvels et al. 2012; Mosser et al. 2012a).

The seismology of red-giant stars emphasizes that an additional mechanism is needed to achieve a slower core rotation of red-giant stars.

7. Conclusions and perspectives

We then computed stellar evolution models of low-mass stars, taking the transport of angular momentum by turbulent viscosity and meridional circulation in radiative zones into account. As found by Palacios et al. (2006) and Eggenberger et al. (2010), the core accelerates very rapidly in our models during the subgiant and giant phases. We showed that the linear approximation for computing rotational splittings remains nevertheless valid for stochastically excited solar-type modes, for all evolutionary stages until the RGB.

For red-giant stars, interpretations of recent Kepler observations lead to a ratio between the core and surface rotation rates of five to ten, far lower than in our models where the ratio between the core and surface rotation rates approaches 10³. Similar conclusions were reached by Eggenberger et al. (2012), which was published during the submission process of the present paper. Nonperturbative calculations lead to complex frequency spectra, with in particular nonsymmetric multiplets (see for instance without Coriolis force Lignières et al. 2006; Reese et al. 2006; Ouazzani et al. 2012, with both centrifugal and Coriolis force). Observations show that several red-giant stars do indeed have such complex spectra. Their core rotation is fast enough to have entangled rotational splittings and mixed-mode spacings. In such cases, rotation may have to be studied with nonperturbative methods. On the other hand, other red-giant stars (Beck et al. 2012) show frequency spectra where the rotation splittings are easily identified as symmetric patterns around axisymmetric modes. The values of the corresponding splittings are quite low, close to or smaller than 1 μHz, and using the linear approximation to derive the rotation splittings from stellar models is fully justified.

We have computed the evolution of a stellar model including rotationally induced transport taking uncertainties in the values of the parameters entering such a description into account. We achieved a decrease in the central rotation rate of a red-giant star by an order of magnitude at most, mostly due to an increase in D_h by two orders of magnitude. The modification of the parameters were chosen so as to enhance the transport of angular momentum so that the core rotation is slowed down as much as possible. The models in the red-giant phase were then slowed enough that the calculations of linear rotation splittings are valid. They were found to be one order of magnitude larger still than observed. This indicates that extraction of angular momentum from the core in our models is not efficient enough.

We need to decrease it by another order of magnitude to reproduce the observations. A similar situation is encountered in the case of the Sun. Several possibilities have been proposed, the main candidates being magnetic fields and internal gravity waves. We showed that such mechanisms must operate or keep operating during the whole subgiant and red-giant phase. In this context, Talon & Charbonnel (2008) show that transport of angular momentum by internal waves generated by the convective envelope can play an important role from the subgiant phase to the base of the RGB. They also find that this transport seems to have no major impact on stars ascending the RGB itself and later on. In the RGB phase, another mechanism must be called for to explain the observations. This will be investigated further in future work. On the observational side, measurements of rotational splittings for subgiant stars are crucially needed to constrain the efficiency for their central layers to slow down.

A standard, validated version of CESTAM will be available for download soon, together with grids of stellar models computed with several rotational velocities and adiabatic oscillation frequencies for selected models.

Acknowledgments

JPM acknowledges financial support through a 3-year CDD contract with the CNES. The authors also acknowledge financial support from the French National Research Agency (ANR) for the project ANR-07-BLAN-0226 SIROCO (SeIsmology, ROtation and COnvection with the CoRoT satellite). A.M., acknowledges the support of the Direccion General de Investigacion under project ESP2004-03855-C03-01. He also acknowledges a stay of two years at the Observatoire de Paris-Meudon.

References

Appendix A: meridional circulation

The vertical component of the meridional circulation is given by (A.1)where c_p,T,g,∇, and ∇_ad have their usual meanings (as defined in Kippenhahn & Weigert 1991). In Eq. (A.1) above, M ^⋆  is the reduced mass, (A.2)where m and ρ_m are the mass and mean density inside an isobar, respectively, and E_Ω and E_μ denote the so-called Ω- and μ-currents (following Palacios et al. 2003). The quantities ϕ and δ are obtained from the equation of state, and are defined as in Kippenhahn & Weigert (1991), (A.3)The Ω- and μ-currents are given by In Eqs. (A.4) and (A.5) we used the following quantities: and  are the nuclear and gravitational energy generation rates, ε_m = L/m the mean energy production inside a sphere of radius r, and ρ_m = 3m/(4πr³) the mean density. The thermal conductivity χ is given by (A.6)and the thermal diffusivity is (A.7)We also used (A.8)the temperature scale height; (A.9)the logarithmic derivatives of χ; (A.10)the logarithmic derivatives of ε; and (A.11)

Appendix B: turbulent coefficients of diffusivity

The coefficient D_eff is given by (Chaboyer & Zahn 1992) (B.1)We use the prescription of Mathis et al. (2004) for the horizontal component of turbulent diffusivity D_h: (B.2)where (B.3)We take the coefficient β = 1.5 × 10^-5 (Richard & Zahn 1999).

For the vertical component D_V, we use the prescription of Talon & Zahn (1997), (B.4)where Ri_c = 1/6 is the critical Richardson number, K the thermal diffusivity, N_T and N_μ the chemical and thermal parts of the Brunt-V*error*˝ais*error*˝ala˝ frequency, : (B.5)

All Tables

All Figures