© ESO, 2013

1. Introduction

Thanks to helio- and asteroseismology, we are able to extract a huge amount of information about solar and stellar structures and compositions (e.g. Turck-Chièze & Couvidat 2011; Aerts et al. 2010) and their internal differential rotation profile (Garcia et al. 2007; Beck et al. 2012; Deheuvels et al. 2012). We know that internal rotation modifies the stellar structure since it generates flows, instabilities, and chemical elements mixing, which modify stars’ evolution, for example, their lifetime and their nucleosynthetic properties (e.g. Maeder 2009, and references therein). Moreover, in order to understand the obtained solar/stellar rotation profiles and the related rotational history, it is essential to develop a complete theory incorporating the different angular momentum transport mechanisms occuring in stellar interiors (e.g. Mathis 2010, and references therein). In this context, one legitimately asks the question about the origin of the internal rotation profiles. Four main processes are responsible, in different ways, for the secular angular momentum transport in radiative interiors. First, a large-scale meridional circulation is driven by structural adjustments of stars, externally applied torques, and internal stresses (e.g. Zahn 1992; Mathis & Zahn 2004; Decressin et al. 2009). Next, rotation profiles may be subject to different hydrodynamical shear and baroclinic instabilities and turbulence (e.g. Knobloch & Spruit 1982; Talon & Zahn 1997; Maeder 2003; Mathis et al. 2004). Then, fossil magnetic fields, trapped during early phases of stellar evolution once radiation zones have been formed (Braithwaite & Spruit 2004; Duez & Mathis 2010), can transport angular momentum through large scale torques and Maxwell stresses (e.g. Gough & McIntyre 1998; Mathis & Zahn 2005; Garaud & Garaud 2008; Strugarek et al. 2011). Finally, internal gravity waves (IGW) that are excited at the convection/radiation boundaries constitute the fourth mechanism able to transport angular momentum over large distances in stellar radiation zones (e.g. Press 1981; Goldreich & Nicholson 1989; Schatzman 1993; Zahn et al. 1997; Talon & Charbonnel 2005; Mathis & de Brye 2012). We note that all these processes are not necessarily present at the same time everywhere in the Hertzsprung-Russel diagram. Moreover, they act on various characteristic timescales in stars of different masses and ages.

The object of this paper is the propagation of IGW and the way they interact with the shear (i.e. the differential rotation) of the surrounding fluid. They are common in the terrestrial atmosphere and oceans (Eckart 1961; Chapman & Lindzen 1970), which is why they are fairly well known in geophysics. We will use this advantage in the stellar case. We here draw attention to the mechanism whereby IGW exchange energy with the mean flow, independently of other dissipative processes such as thermal and viscous diffusion. Indeed, when the frequency of excited waves is of the same order as the angular velocity of the fluid (we present the accurate definition later), a phenomenon of resonance occurs, which affects the properties of both the wave and the shear of the surrounding fluid. This phenomenon is called a critical layer. Under the assumption of a perfect fluid (neither a heat conductor nor viscous), Booker & Bretherton (1967) and Lindzen & Barker (1985) have provided first results about critical layers in the geophysical case. They have shown that, depending on the value of the Richardson number of the fluid, which compares the relative strength of the shear and the stable stratification (see Eq. (3)), the waves might be either attenuated or reflected by the critical layer. This reflection may even be an over-reflection together with an over-transmission (see also Sutherland & Yewchuk (2004) for a laboratory evidence). At the same time, Koppel (1964), Hazel (1967), Baldwin & Roberts (1970), Van Duin & Kelder (1986) have completed this work by taking into account the conduction of heat and the viscosity of the fluid. Surprisingly, their conclusions about the role of critical layers are identical. In the astrophysical case, Rogava et al. (2007) have studied critical layers for gravito-alfven waves. However, all these authors have produced their studies in Cartesian coordinates, assuming that the stiffness of the domain was thin enough to neglect its curvature. In the case of stellar radiation zones where critical layers may play an important role (e.g. Barker & Ogilvie 2010; Barker 2011; Rogers et al. 2012), these equations should thus be generalized to the case of spherical coordinates to be able to treat deep spherical shells.

Therefore, after exposing our notations and assumptions (Sect. 2), we present a complete mathematical study of critical layers in the case of a perfect fluid (Sect. 3) and of a non-perfect fluid (Sect. 4) in spherical geometry. Then, we determine the related mean vertical flux of angular momentum transported by IGW (Sect. 5). Next, we implement our theoretical results in the dynamical stellar evolution code STAREVOL and apply our formalism to the evolution of a one solar-mass star (Sect. 6). Finally, we present the conclusion and perspectives of this work (Sect. 7).

2. Definition, notations, and hypotheses

The star we consider is composed of (at least) a convective and a radiative region. IGW are excited at the boundary between these two regions. They propagate in the radiative zone and are evanescent in the convective zone (e.g. Press 1981). As explained in the Introductions, the study of IGW in geophysics is usually made in Cartesian coordinates, justified by a local approach. In the stellar case, however, a global approach using the spherical coordinates (r,θ,ϕ) is necessary.

In the frame of Zahn (1992), we choose a shelular angular velocity for the studied star’s radiation zone, , considering that, because of the strong stable stratification, the shear instability decreases the horizontal gradient of the angular velocity (e.g. Talon & Zahn 1997; Maeder 2003; Mathis et al. 2004). This consequently can be considered as only dependent on radius. For the moment, we neglect the action of Coriolis and centrifugal accelerations while the Doppler shift due to differential rotation is retained. We also neglect the action of a potential magnetic field.

Now, we need to introduce some quantities to describe the properties of the fluid in the studied radiative zone. Each scalar field X is written as (1)where we have introduced its horizontal average on an isobar, , and its associated fluctuation X′. The thermodynamic variables employed are the density ρ, the pressure p, the temperature T, and the specific entropy S. Next, the stratification is described in terms of the Brunt-Vaïsälä frequency (2)where is the mean gravity in the Cowling approximation, where fluctuations of the gravific potential are neglected (see Cowling 1941), and the adiabatic exponent. The relative importance of the stable stratification restoring force and the shear destabilizing effects is quantified by the Richardson number (3)When Ri is small, the velocity shear overcomes the stabilizing buoyancy and turbulence and mixing occur (e.g. Talon & Zahn 1997). In contrast, when Ri is large, the fluid remains stable. Finally, for the case of a non-perfect fluid (Sect. 4), we introduce the viscosity ν of the fluid and the coefficient of thermal conductivity κ. These notations will be recaled at the proper moment.

IGW themselves are characterized by their relative frequency (4)where σ_w is their excitation frequency (from the base of the convective zone in low-mass stars or the top of the convective core in intermediate and high-mass stars), m corresponds to a Fourier expansion along the longitudinal direction (Eq. (10)), and is the difference between the angular velocity at the level r and at the border with the convective zone. The introduction of m leads to the definition of two classes of waves. Prograde (respectively retrograde) waves correspond to negative (respectively positive) values of m.

Finally, all these notations allow us to define properly a critical layer. It arises for a wave whose relative frequency σ(r) becomes zero. The coherence with the qualitative definition given in the introduction is respected since σ(r_c) = 0 means that the excitation frequency of the wave equals  − m times the angular velocity of the fluid. It corresponds to a corotation resonance. It is important to highlight that the position of the critical layer depends on both wave characteristics (m and σ_w) and shear properties (ΔΩ(r)). In our approximation, σ(r_c) = 0 means that is constant, so the critical levels are isobaric surfaces of the star.

3. Case of the perfect fluid

In this first approach, we treat the hydrodynamic equations assuming that the fluid is neither viscous nor a heat conductor.

3.1. Equation of propagation of IGW near a critical layer

Our aim is to calculate the Eulerian velocity field of a fluid particle. We introduce the time t and the unit vectors (ê_r,ê_θ,ê_ϕ), which are associated with the classical spherical coordinates (r, θ,ϕ). The velocity field is (5)where u is the velocity associated with the wave and the azimuthal velocity field of the differential rotation. We note that we have neglected all other large-scale velocities such as meridionnal circulation. Then, we give the linearized equations of hydrodynamics governing IGW dynamics in an inertial frame, i.e. the momentum, continuity, and energy equations: (6)where D_t = ∂_t + ΔΩ∂_ϕ. The equation of energy conservation is obtained using the linearized equation of state (7)and assuming the ideal gaz law (8)where ℛ is the gaz constant and c_P the specific heat-per-unit mass at constant pressure.

We define the Lagrangien displacement ξ as u = D_tξ and following Rieutord (1986), we expand it on the vectorial spherical harmonics basis () defined by (9)where are the spherical harmonics with l ∈ N and m ∈ ˜ − l,l¨. Thus, (10)Then, we decompose ρ′ and p′ using spherical harmonics: and we obtain a new system made up of radial equations of momentum (13)of mass conservation (14)and of energy in the adiabatic limit (15)The combination of these relations leads to the system presented by Press (1981): (16)where .

In stellar radiative regions, the transport of angular momentum is dominated by low-frequency IGW with σ ≪ N, where N is the Brunt-Vaïsälä frequency defined in Eq. (2). It allows us to apply the anelastic approximation (Press 1981), where acoustic waves are filtered out. Introducing (17)we obtain the following equation of propagation: (18)As the right-hand side of Eq. (18) is of order with H_P the characteristic pressure or density height scale, it can be neglected if (19)which is the case here. Finally, we obtain the equation of propagation of IGW in a perfect fluid: (20)with (21)In Cartesian coordinates, this equation is called the Taylor-Goldstein-Synge (TGS) equation. We observe that the value r = r_c, where σ(r_c) = 0, is a singular point for this equation. Thus, we now focus on the study of the behaviour of this equation around such a critical point. Then for r ≈ r_c, the equation of propagation becomes (22)where (23)is the value of the Richardson number at the critical level, and (24)the horizontal wavenumber at the critical layer. We saw in the introduction (Eq. (3)) that the Richardson number is relevant to distinguish between the relative importance of shear and stratification. We describe this divergence of cases from a quantitave point of view in the following part.

3.2. Mathematical resolution with the method of Frobenius

The Frobenius method offers an infinite series solution for a second-order ordinary differential equation of the form in the vicinity of the singular point z₀ = 0. For more details, a mathematical description of this method can be found in Teschl (2011).

In our case, p(r) = 0, , and the theory brings out two cases depending on the value of l, m, and Ri_c:

• Case 1: The solution is (25)where A₁ and B₁ are constants representing the wave amplitude.

• Case 2: We define the complex parameter (26)and the solution is given by (27)where A₂ and B₂ are constants.

We exclude this special value of Ri_c and choose to consider only two different cases for the rest of the paper: and . The results can then be extended to the case without recourse to the logarithmic solution (Van Duin & Kelder 1986).

We now discuss the hydrodynamical behaviour corresponding to the situations where and Applying the classical method exposed in Drazin & Reid (2004) to the generalized spherical TGS equation (Eq. (22)), we can identify that the first regime corresponds to the case where the fluid stays stable with respect to the vertical shear instability at the critical layer. In the other case, this instability and thus turbulence can develop. This clearly shows how it is necessary to go beyond the currently used formalisms for stellar evolution, where IGW and vertical shear instability are considered as uncoupled. Indeed, if a fluid becomes shear unstable, mixing occurs that modifies the local stratification and thus IGW propagation (see e.g. Brown & Sutherland 2007; Nault & Sutherland 2007). Thus, we now distinguish the case of critical layers when the fluid is stable from the one when it is unstable.

3.3. The case of stable critical layers: Ri

It is now time to understand the physical behaviour of the solutions given above. First of all, we propose to cut the physical domain into two parts: above and below the critical layer. As , η_l,m defined by (26)is a purely imaginary number. In order to clearly distinguish between real and imaginary parts, we introduce (28)The two solutions can be written as (29)where (resp. ) is available when r > r_c (resp. r < r_c). Booker & Bretherton (1967) and Ringot (1998) both explain a way to connect these solutions, considering that the term (r − r_c)^{1/2 + iα_l,m} can be compared to an upward propagative wave of the form e^ikr and respectively that (r − r_c)^{1/2 − iα_l,m} can be compared to a downward propagative wave of the form e^−ikr. We obtain the following identification: (30)In order to connect the solutions, we observe the behaviour of (r − r_c) above and below the critical level. As r − r_c decreases from positive to negative values, its complex argument changes continuously from 0 to  − π (Ringot 1998). Mathematically, we get It follows that the solutions above and below the critical layer can be written as (33)Physically, these equations can be explained this way. Starting from above the critical layer, the downward propagating wave passes through the critical layer and is attenuated by a factor equal to e ^{− α_l,mπ}. At the same time, starting from below the critical layer, the upward propagating wave is attenuated by the same factor and its amplitude becomes equal to A. We also underline that both waves take a phase difference when they cross the critical layer. In Fig. 1, we represent the attenuation rate of different waves defined by the numbers l and m passing through a critical layer. The greater the ratio , the stronger is the attenuation Att  for the same value of Ri_c. The axis scale depends on l and m because the condition of validity of this result is . We so deduce that waves of high ratio  (which do not necessarily correspond to high order) are strongly attenuated if they reach their critical layer.

Fig. 1 Attenuation rate Att =  of the wave passing through a critical layer as a function of the Richardson number. We observe that Att increases with

We have not yet discussed a latter point: the choice of the method of resolution. Most of the publications concerning IGW use another process to solve the equation of propagation (Press 1981; Zahn et al. 1997; Mathis 2009). In fact, the Wentzel-Kramers-Brillouin-Jeffreys (WKBJ) theory is particularly adapted to the resolution of this equation. However, it is not convenient in our case because it imposes a condition on the value of the Richardson number as demonstrated in appendix A. It is shown that the WKBJ approximation is available only if . Despite this restriction, we write the solution: by separating the domain into two parts, we obtain (34)As (35)it becomes: (36)It is comforting to see that both methods (Frobenius and WKBJ) give the same solution to a multiplicative constant when the Richardson number becomes high: (37)

Fig. 2 Assumed neighbourhood of an unstable critical layer for the calculation of the IGW’s reflection and transmission coefficients. We assume that the unstable region around the critical layer (in red) has a thickness given by 2δ (zone II). The surrounding regions where IGW are propagative are in beige (zones I and III).

3.4. The unstable case: Ri

In the unstable regime, the Frobenius method also gives a solution but we are not able to identify upward and downward propagating waves because of shear-induced instability and turbulence. As a consequence, we can not connect the solutions at the critical layer. In order to avoid this difficulty, we here propose to solve Eq. (22) following the method developped by Lindzen & Barker (1985). They applied it in the case of Cartesian coordinates, and our work is to generalize it to spherical coordinates. The parameter η_l,m defined in Eq. (26) is real in this case, and we introduce (38)Equation (22) becomes (39)We are seeking solutions of the form , where Φ_l,m is the solution to the Bessel equation (40)Consequently, Φ_l,m is a combination of the Bessel’s modified functions I_ηl,m(X) and I _− ηl,m(X): (41)The final solution is given by: (42)We would like to calculate the reflection and transmission coefficients of the wave passing through the unstable critical layer. We assume that the fluid has the profile described in Fig. 2. We decompose the studied region into three zones defined by the value of the quantity . In zones I and III, the Richardson number is high enough to allow us to apply the WKBJ method described in the previous section. However, around the studied critical layer, in unstable zone II, we must use the solution with the modified Bessel functions. Moreover, we assume here that the thickness 2δ of this zone is small in comparison with the characteristic length of the problem. Then, we can consider that the wavenumber k_V is constant as (43)We consider a wave coming from the overside of the critical layer (zone I). It is partly transmitted toward zone III and partly reflected backward to zone I. We thus write the solutions corresponding to the three zones: (44)where Ψ_I(r) is available if (r − r_c) > δ, Ψ_II(r) if  − δ < (r − r_c) < δ and Ψ_III(r) if (r − r_c) >  − δ. The coefficients A, B, T, and R are calculated by to the following four continuity equations:(45)They correspond to the continuity of the solution and of its first derivative, which physically means that both displacement and mechanical stresses are continuous. After some algebra, we obtain the coefficients R and T of reflection and transmission, which depend on the stiffness δ, the vertical and horizontal wavenumbers k_Vc and k_Hc (see Eqs. (43) and (24)) , and the variable η_l,m (Eq. (26)). In order to simplify the formula, we note I _± ηl,m instead of I _± ηl,m(k_Hcδ). Then, we obtain (46)with and (47)with We have now calculated the transmission and reflection coefficients of an IGW, through an unstable region around a given critical layer, where . We represent the level lines of |R| and |T| in Fig. 3 for l = 4 and m =  ± 3. They are plotted as a function of the Richardson number at the critical layer Ri_c, growing from 0 to its maximum value defined by . For instance, in Fig. 3, . The other variable is the half-thickness δ of the unstable layer (zone II), arbitrarily chosen, that points the non-local character of unstable turbulent layers. The main result is that both coefficients are greater than 1 when is small enough. Consequently, for a low Richardson number at the critical layer, the wave can be over-reflected and over-transmitted at the same time. This means that, in contrast to the first stable case, the wave takes potential energy from the unstable fluid and converts it into kinetic energy. In other words, the turbulent layer acts as an excitation region. If |R| < 1 and |T| < 1, we speak about IGW tunneling (Sutherland & Yewchuk 2004; Brown & Sutherland 2007; Nault & Sutherland 2007). Another remark concerns the dependency of |R| and |T| on m. We note that the sign of m does not matter since only its square appears in the expressions. Physically, it shows that the critical layer’s action is the same on prograde and retrograde waves. This point is of importance because we know that other dissipative processes occuring during the propagation of IGW discriminate between both types of waves. In order to visualize the action of the critical layer on different waves, Fig. 2 shows the level lines |R| = 1 for 1 ≤ l ≤ 5 and 1 ≤ m ≤ l. We previously said that it is useless to consider negative values of m since |R| depends only on m². The pair (l,m) is indicated on each line, followed by the value of . Lines in the same colour correspond to the same value of l. These lines mark out the limit to observe over-reflection for a chosen wave. We observe that the higher the value of , the stronger is the condition on Ri_c to observe an over-reflection.

Fig. 3 Level lines of reflection  |R|  and transmission  |T|  coefficients of an IGW at a critical layer as a function of the Richardson number Ric and the thickness δ of the critical layer (zone II in Fig. 2). The top panels represent the level lines for an arbitrarily chosen value of (l,m) = (4, ± 3), while the bottom panels show level lines  |R|  = 1 and  |T|  = 1 for different couples (l,m).

3.5. Choice of the method

We decided to apply different methods to solve the stable and unstable cases. However, it could be legitimate to wonder if both methods are equivalent from a mathematical point of view. In this section, we present a short comparison between the solutions obtained with the method of Frobenius and the one with Bessel functions. The modified Bessel function I_ηl,m(X) can be computed using (48)At the neighbourhood of the critical layer, X tends to 0 and the first-order expression is (49)Close to the critical layer, the global solution given in Eq. (42) is then (50)Now, we recall the expression of the solution given by the method of Frobenius (Eq. (27)), rewritten with the previous notations (51)In conclusion, for a fixed couple (l,m), and vary in the same way as function of r − r_c.

4. Case of the non-perfect fluid

Until now, we have studied the role of the critical layers assuming that the fluid was perfect. In order to make the problem more realistic, we include in the second part the viscosity ν of the fluid and the coefficient of thermal conductivity κ (e.g. Koppel 1964; Hazel 1967; Baldwin & Roberts 1970; Van Duin & Kelder 1986).

4.1. Equation of propagation of IGW near a critical layer

The linearized equations of hydrodynamics in Eq. (6) become (52)The following method for forming the propagation equation is adapted from the work of Baldwin & Roberts (1970). We assume that the mean density of the fluid takes nearly a constant value, that is to say that is small compared with the characteristic lengths of the problem. First, we project Eq. (52) onto ê_r and apply the operator ∇². Then we apply D_t − κ∇². After combination with the two other equations, we obtain (53)As in the previous section, we decompose the radial velocity on the basis of spherical harmonics (54)where σ = σ_w + mΔΩ(r). Moreover, it is easier to work with dimensionless numbers. For this reason, we introduce the notations detailed in Table 1.

Then, Eq. (53) becomes for each pair (l, m) such as l ∈ N and m ∈ [[− l,l]]: (55)where Δ_l is the scalar spherical Laplacian operator: (56)Lastly, we introduce and due to a development close to the critical layer obtain the sixth-order equation: (57)

where now . Equation (57) can be compared with the one obtained by Press (1981), who has neglected the viscosity ν. Moreover, we can see that if we take ν = κ = 0, and without forgetting that P_r also depends on ν, Eq. (57) is identical to Eq. (22) for the perfect fluid.

4.2. Mathematical resolution

The resolution of Eq. (57) requests several substitutions and is quite complex. The detailed calculation can be found in Appendix B, and we give only the main steps here. We drew our inspiration from Hazel (1967), Baldwin & Roberts (1970), Koppel (1964), Van Duin & Kelder (1986), who solved the same equation in Cartesian coordinates. However, the resolution in spherical coordinates has some differences. The aim is to rewrite the equation here under a known form: the Whittaker differential equation; then, after some algebra, we can show that Eq. (57) can be written in the following form: (58)The solutions of Eq. (58) are thus the Whittaker functions (Abramowitz & Stegun 1965): (59)where is the confluent hypergeometric function of Kummer (60)with and (k)₀ = 1, being the usual gamma function, (61)and (62)Because of this solution, we obtain the expression of the radial displacement of the wave: (63)We will clarify the function f(κ,P_r) of the thermal diffusivity coefficient later.

There is still the last step to carry out. We must find a curve along which we integrate the solution, that is to say, we determine a and b. During the calculation detailed in the Appendix B, Eq. (B.3), (64)has not been used yet. A sufficient condition to make this equality true is . Considering that t³ = |t|³e^3iθ_t (θ_t being the complex argument of t), we have (65)There are several possibilities for the choice of curves. Koppel (1964) proposed building a basis of six solutions using the curves (,,) represented in Fig. 6 (left). Therefore, the solutions of the sixth-order equation (Eq. (57)) are linear combinations of (66)and V_i(η) with i ∈  { 1,2,3 } , where V_i(η) corresponds to U_i(η) with the opposite sign for M_l,m.

Fig. 4 Left: curves (,,) defining the basis (U1,U2,U3,V1,V2,V3). Right: curves (,,) defining the basis (u1,u2,u3,v1,v2,v3).

4.3. Application to IGW

It is now time to apply these mathematical results to IGW. We recall the solution for the perfect fluid obtained in Sect. 3. For the moment, it is not necessary to distinguish between the stable regime, i.e. , and the unstable one, i.e. . The Frobenius solution is available everywhere if we take as a complex number. Remembering that , the radial Lagrangian displacement is (67)For easier reading, we have removed the indices r;l,m. Subscript P designates the solution for the perfect fluid.

Concerning the solution for non-perfect fluids, Hazel (1967) proposed using another basis for the solutions of Eq. (57). The curves (, , ) are represented in Fig. 6. The new basis (u₁,u₂,u₃,v₁,v₂,v₃) is defined by

• u₁ = U₁ + U₂ + U₃ and v₁ = V₁ + V₂ + V₃,

• and ,

• and ,

where X^∗ is the complex conjugate of X. Consequently, the solution for the non-perfect fluid (subscript NP) can be written as (68)A relation between both solutions and exists if we consider that the Reynolds number is great. In the case of the Sun (e.g. Brun & Zahn 2006), the microscopic viscosity in the radiative zone is weak. For this reason, it is appropriate to consider that the Reynolds number is much greater than 1. This assumption leads to the relation (69)Baldwin & Roberts (1970) give tables for the asymptotic behaviour of u_i and v_i (i ∈ 1,2,3). The solutions u₁, u₂, v₁ and v₂ diverge when R_e →  + ∞. They are therefore physically unacceptable and we deduce that (70)Tables 2 and 3 give the asymptotic expressions of u₃ and v₃ as a function of the sign of r − r_c.

In the stable case where , is a real number. The same findings as in the previous section can be made: after the passage through the critical layer, the wave is attenuated by a factor e ^{− α_l,mπ}. In the unstable case, where , is a complex number, so we cannot interpret the solution as upward and downward propagating waves. But the expressions given in Tables 2 and 3 remain comparable to those in the first section, and we deduce that the calculation of R and T will lead to the same result: the possibility of over-reflection and over-transmission.

4.4. Radiative and viscous dampings

4.4.1. General equations

We volontarily left aside the factor f(κ,P_r) in Eq. (66). In order to establish its expression, Zahn et al. (1997) used the equation of the propagation of IGW, taking into account heat diffusion but with a viscosity coefficient (ν) equal to zero. Here, we generalize their result for ν ≠ 0 (i.e. P_r ≠ 0) and obtain (71)where (72)We introduced the general expression for N², the Brünt-Väisälä frequency, to be able to take chemical stratification into account. Then, we have (73)with and where is the pressure heigh-scale, the temperature gradient, and the mean molecular weight (μ) gradient. Moreover, we introduced the generalized equation of state (EOS) given in Kippenhahn & Weigert (1990): (74)where and . Next, σ = σ_w + mΔΩ is the Doppler-shifted frequency of the wave relative to the fluid rotation with an excitation frequency σ_w, r_c and r_CZ are respectively the positions of the critical layer and the boundary between the studied radiative zone and the convection region where IGW are initially excited. This damping is another source of attenuation independent of the presence of a critical layer. Moreover, as shown in Sects. 3.4. and 4., we have to consider IGW reflected and transmitted by unstable critical layers in addition to those initially excited by convection. Then, we introduce the notation (75)

where r₁ and r₂ are respectively the emission point of the IGW and the position with r₁ > r₂ (in the opposite case, where r₁ < r₂, limits in the integral have to be reversed). This will enable us to describe radiative and viscous dampings in every configuration. We note that in stellar radiation zones P_r ≪ 1 (Brun & Zahn 2006) inducing that the damping is mostly radiative. We now compare it with the effects of critical layers.

4.4.2. Prograde and retrograde waves

For a same environment, prograde waves have a frequency lower than retrograde ones (e.g. Eq. (4)). We chose, for example, a couple of IGW with the same excitation frequency, σ₀, the same number l and opposite azimuthal degrees m. Thus, we compared a prograde wave of frequency σ₁(r) = σ₀ − |m| ΔΩ and a retrograde one of frequency σ₂(r) = σ₀ + |m| ΔΩ. We obtained σ₁(r) − σ₂(r) =  − 2 |m| (Ω(r) − Ω_CZ) < 0 in the presence of negative gradient of Ω and σ₁(r) − σ₂(r) > 0 if the gradient is positive. We stress that τ given in Eq. (71) varies globally as . As a consequence, assuming that a negative Ω-gradient is present near the excitation layer (see Sect. 6), the radiative damping is stronger for prograde waves than retrograde waves. Therefore, prograde IGW are absorbed by the fluid much closer to their region of excitation, while the retrograde waves are damped in a deeper region. This process is responsible for the net transport of angular momentum by IGW in stars. On the other hand, critical layers do not introduce such net bias between prograde and retrgrade IGW because their effects scale with m².

4.4.3. Dependency in l and m

The second comment concerns the variation of radiative and viscous dampings as a function of l and m. On one hand, looking at the multiplicative factor, which is in front of the integral in Eq. (72), we can roughly write that . On the other hand, the expression of the attenuation factor due to stable critical layers is proportional to , which is always greater than unity since |m| ≤ l. The comparison between radiative and viscous dampings and the one due to stable critical layers for an IGW with given will be examined in Sect. 6.2.3.

4.4.4. Location

Lastly, radiative and viscous dampings occur throughout the whole propagation of IGW. In contrast, critical layers are localised. Moreover, there is a condition for a wave to reach a critical layer: the rotation speed of the fluid must be of the same order as the wave frequency for a chance of observing σ = 0. Finally, all IGW are concerned by radiative and viscous dampings, which increase around critical layers since τ ∝ σ^-4.

In Table 4, we sum up the different cases studied in these work. Depending both on the properties of the fluid and the studied wave, the wave is submitted to different phenomena.

5. Transport of angular momentum

As stated in the introduction, our goal is to study the transport of angular momentum in stellar radiation zones and to unravel the role of critical layers. Therefore, the first step is to recall the flux of angular momentum transported by propagative IGW and by the shear-induced turbulence. To illustrate our purpose, we focus here on the case of a low-mass star, where IGW are initially excited at the border of the upper convective envelope (see Fig. 5) by turbulent convection (Garcia Lopez & Spruit 1991; Dintrans et al. 2005; Rogers & Glatzmaier 2005; Belkacem et al. 2009; Brun et al. 2011; Lecoanet & Quataert 2013) and by tides if there is a close companion (Zahn 1975; Ogilvie & Lin 2007) with an amplitude A (corresponding results for massive stars with an internal convective core can be easily deduced by reversing signs).

Fig. 5 Two configurations in a low-mass star with an external convective envelope. Left: the case of a stable critical layer (CL) where IGW are damped. Right: the case of an unstable critical layer where IGW can be over-reflected/transmitted.

5.1. Angular momentum fluxes

5.1.1. Angular momentum flux transported by propagative IGW

First, we have to calculate the angular momentum flux transported by a propagative monochromatic wave over a spherical surface. It is given by the horizontal average of the Reynolds stresses associated with the wave (e.g. Zahn et al. 1997): (76)where . Besides, Mathis (2009) shows that (77)where ℱ_E;l,m,σ is the horizontal average of the energy flux in the vertical direction expressed by Lighthill (1986) as (78)So finally, the angular momentum flux is given by (79)To calculate this angular momentum flux, we need the expressions of u_r;l,m and . The expression of u_r;l,m is immediately accessible from Eq. (33), solution of Eq. (22). It is a little more complicated for  because we need to go back in the calculation leading to Eq. (22). The expression of results from the first equation of the system given in Eq. (16), which can be reduced to (80)applying the anelastic approximation and neglecting terms of order 1/L².

The total vertical flux of angular momentum transported by propagative IGW is then given by (81)and we define the associated so-called action of angular momentum as (82)

5.1.2. Angular momentum flux transported by shear-induced turbulence

In the case of shear-unstable regions, such as region II in the regime where , IGW are unstable (e.g. Drazin & Reid 2004). After the non-linear saturation of the instability, a steady turbulent state is reached. Then, as established, for example, in Zahn (1992); Talon & Zahn (1997), and now confirmed by numerical simulations (see Prat & Lignières 2013), the vertical flux of angular momentum transported by shear-induced turbulence is given by (83)where (84)The adopted value for the critical Richardson number is Ri^c = 1/6, and ν_H is the horizontal turbulent viscosity for which we assume the prescription derived by Zahn (1992).

5.1.3. Equation of transport of angular momentum

We now refocus these results in the wider frame of the complete angular momentum transport theory. Considering the other transport mechanisms we presented in the introduction, the angular momentum transport equation that takes into account meridional flows, shear-induced turbulence, and IGW (e.g. Talon & Charbonnel 2005; Mathis 2009) becomes (85)The first term on the right-hand side corresponds to the advection of angular momentum by the meridional circulation, where U_r = U₂   P₂(cosθ) is its radial component. The Lagrangian derivative takes into account the radial contractions and dilatations of the star during its evolution, which are described by ṙê_r. According to our hypothesis, we do not take into account the transport by the Lorentz force, associated with magnetic fields in stellar radiative zones (Mathis & Zahn 2005).

The major difference with previously published equations is that IGW and shear-induced turbulence transports of angular momentum are not summed linearly since they are, as we demonstrated before, intrinsically coupled. Therefore, for stable regions, one must consider the Reynolds stresses of IGW (Eq. (82)), while for unstable regions, only the vertical turbulent flux given in Eq. (90) must be taken into account.

5.2. Stable critical layer

5.2.1. Case of the perfect fluid

This is the simplest case of a stable critical layer in a perfect fluid. We apply Eq. (79) to obtain the expressions for the transported fluxes by propagative IGW below and above the critical layer: (86)where A is the initial amplitude of the IGW at r = r_CZ, (see Eq. (28)) and (87)where are the associated Legendre polynomials. Booker & Bretherton (1967) obtained similar results in Cartesian coordinates. We see the expected attenuation of the flux by a factor e ^{− 2πα_l,m} when the wave passes through the critical layer. Moreover, ℱ_J;l,m,σ ±  depends on m (and not on m²). As a consequence, we recover the classical result that prograde waves (m < 0) and retrograde ones (m > 0) have opposite angular momentum flux (respectively a deposit and an extraction). Finally, the monochromatic action of angular momentum ℒ_J;l,m,σ(r) = 4πr²ℱ_J;l,m,σ is constant in each region because of the absence of dissipation.

5.2.2. Case of the non-perfect fluid

We saw in the previous section that the solution of the equation of propagation in the case of a non-perfect fluid is similar to the one obtained for a perfect fluid. For this reason, we are able to apply the same method for the calculation of ℱ_J;l,m,σ. We obtain (88)The difference to Eq. (86) comes from the introduction of radiative and viscous dampings (Eq. (75)).

The conclusion is that, in the stable case, the attenuation due to the passage through a critical layer is added to dampings due to dissipation. This observation leads us to simply implement the role of stable critical layers as an additional term in the damping coefficient.

5.3. Unstable critical layer

5.3.1. Region I: r_c + δ ≤ r ≤ r_CZ

Using the summary given in Fig. 5 for the unstable case and the obtained results for regions where IGW are propagative (Sect. 5.2.2), we obtain: (89)where we identify the transport induced by the incident wave, which propagates downward, and the one induced by the reflected one, which propagates upward. We see that the unstable critical layer is a second excitation source for IGW propagating in region I and that the angular momentum transport will be modified in such situation, particularly when over-reflection (|R|  > 1) occurs.

5.3.2. Region II: r_c − δ < r < r_c + δ

In this unstable region, we directly use results obtained in Sect. 5.1.2. to describe the flux of angular momentum transported by the shear-induced turbulence, i.e.: (90)

5.3.3. Region III:

Using the same methodology as for region I, we get in a straightfoward way (91)where we indentify the transport induced by the transmitted wave, which propagates downward. As in region I, we can see that the unstable critical layer constitutes a secondary excitation source for IGW propagating in region III and that the angular momentum transport will be modified, particularly when over-transmission (|T| > 1) occurs.

Since the general theoretical framework has been given, we now explore the possibility of the existence of the two different regimes (stable and unstable) along the evolution of a given star. As a first application, we chose to study the case of a solar-type star which has already been studied without critical layers by Talon & Charbonnel (2005).

6. A first application: the evolution of a solar-type star

6.1. The STAREVOL code

We used the one-dimensional hydrodynamical Lagrangian stellar evolution code STAREVOL V3.10 (see Lagarde et al. (2012) and references therein for a detailed description of the input physics). We recall the main characteristics and parameters used for the modelling that are directly relevant for the present work. We use the Schwarzschild criterion to determine the position of convective zones, and compute their temperature gradient according to the mixing length theory with a α_MLT = 1.75. The solar composition is taken from Asplund et al. (2005) with the Ne abundance from (Cunha et al. 2006). We generated the opacity table for temperatures higher than 8000 K following Iglesias & Rogers (1996) by using their website¹. The opacity table at lower temperatures follows Ferguson et al. (2005)². The mass loss rate is determined following Reimers (1975) with a parameter η_R = 0.5. The increase of mass loss due to rotation is taken into account following Maeder & Meynet (2001). However due to the small mass loss and velocity of our models this effect remains weak.

In radiative regions, we follow Mathis & Zahn (2004) formalism for the transport of angular momentum and chemicals as well as the prescription from Talon & Zahn (1997) for the vertical turbulent transport (Eq. (84)). We assume that convective regions rotate as a solid body. The treatment of IGW follows Talon & Charbonnel (2005, 2008), with the difference that the volumetric excitation by Reynolds stresses in the bulk of convective zones (e.g. Goldreich & Kumar 1990; Goldreich et al. 1994; Belkacem et al. 2009) is consistently computed at each time-step as a function of their physical properties.

We start with an initial model of 1.0 M_⊙ at solar metallicity, with an initial surface velocity of 70 km s^-1. The rotation profile is initially flat. We add magnetic braking through the following law: with a constant K = 3 × 10³⁰. This value has been calibrated to reproduce the surface velocity determined in the Hyades by Gaige (1993).

6.2. The effects of critical layers

6.2.1. Location of critical layers

We studied theoretically the impact of the critical layer for a given IGW, thus assuming that there are some radii where the relation σ(r) = 0 is satisfied. As a consequence, the first question we may answer due to the simulation concerns the existence of such critical layers and their location in the radiative zone. Figure 6 shows that critical layers do exist in the radiative core of the solar-like star studied here. In the three panels, we superimpose the rotational velocity of the star’s interior as a function of the normalized radius with the position of potential critical layers, which are indicated by colored squares that correspond to positions where σ = σ_w + mΔΩ = 0. Each panel corresponds to a given value of the excitation frequency σ_w. As expected, the positions of the critical layers only depend on the azimuthal number m of the wave, and not on the degree l. Through this plot, we confirm that critical layers exist in the solar-like star. However, we did take into account the fact that not all waves can reach these positions.

Fig. 6 Position of the critical layers for retrograde IGW for three different frequencies, superimposed with the rotational velocity profile in the star.

First of all, we only consider retrograde waves because the prograde ones are damped immediately after their excitation (see Sect. 4.4.2). Moreover, in the left-hand panel of Fig. 7, we represente in logarithmic scale the luminosity of a given IGW at the moment of its initial stochastic excitation by the turbulent convection as a function of its degrees l and m, following the spectrum adopted in Talon & Charbonnel (2005). The difference between the three plots is the value of the excitation frequency σ_w. Through this representation, we can see that the maximum of excitation lies in a domain where l and m are close and quite small. Moreover, it shows that for each excitation frequency, the amplitude of the excited wave depends on l and m and may be close to zero. As a consequence, some critical layers represented in Fig. 6 may belong to a non-excited wave or to a wave completely damped at this depth. Fortunately, results obtained show that some waves really meet their critical layer.

Fig. 7 Left: luminosity of the waves at the location of their excitation as function of l and m. Right: ratio between τCL and τtot = τCL + τrad (see Eq. (92)).

6.2.2. Interaction between waves and critical layers

Concerning the way IGW interact with the surrounding fluid, the theory predicts two possible regimes depending on the value of the Richardson number Ri_c at the critical level. In our simulation, it appears that for every detected critical layer, the relation , which correponds to the stable regime, is verified. As a consequence, we establish that the second unstable regime (with the associated possible tunneling or over-reflection and transmission) does not occur for the solar-like star of our simulation; forthcoming studies may explore other types of stars at different evolutionary stages, to see if this regime can occur. Therefore, for the solar-type star studied here, we only implement in STAREVOL the terms related to stable critical layers: each time a wave passes through a critical layer, it is damped with a coefficient  (see Eq. (88)).

6.2.3. Effect of critical layers

Since all interactions between waves and critical layers are of the same kind in the star studied here, we can concentrate on the quantitative importance of their effect on the transport of angular momentum. We know that in the stable regime (Sect. 5.2.2.), the wave passing through its critical layer is damped by the factor given in Eq. (88) that is added to radiative damping (here P_r ≪ 1 and the viscous damping is thus negligible), which has already been taken into account in previous works (e.g. Talon & Charbonnel 2005). In the right-hand panel of Fig. 7, we thus choose to represent the ratio between τ_CL, the rate of attenuation due to the passage of the wave through its stable critical layer, and the sum τ_tot of this rate and that of the radiative damping. The explicit formula is (92)In contrast to the left-hand panel, only waves are represented here which meet a critical layer. That is the reason why several white zones are seen. They correspond to the waves which have been attenuated before reaching the depth of their critical layer. In the case of low σ_w (upper panels in Fig. 7), the high-degree waves (l > 45) are simply not excited, as we can see on the right. In red zones, the role of critical layers is important in comparison with the radiative damping, while dark blue regions are those where the radiative damping dominates. Therefore, this figure shows that critical layers should be taken into account.

6.2.4. Evolution of the rotation profile

Fig. 8 Evolution of the rotation profile, where the role of the critical layers in the transport of angular momentum is taken into account. The curves are labelled according to the corresponding ages in Gyr. Parameters are indicated above where K is the braking constant. Right: comparison between the evolution of the rotation with (red) and without (black) taking into account the effect of the critical layers.

We now concentrate on the evolution of the rotation profile, when momentum deposition due to the interaction between waves and critical layers is taken into account. The other transport mechanism (IGW’s radiative damping effects, shear-induced turbulence, and meridional circulation) have been previously implemented in the code (Talon & Charbonnel 2005).

In the centre of the star, the rotation velocity is lower and increases the radiative damping of retrograde waves (see Eq. (72)). This forms an angular momentum extraction front, which propagates from the core to the surface to damp the differential rotation. Three fronts are seen in the top panel of Fig. 8 and have the same form as those already obtained in Talon & Charbonnel (2005). To isolate the action of critical layers on the evolution of the rotation profile, we superimposed in the bottom panel of Fig. 8 the surface velocity as function of the evolution time with (black line) and without (red line) critical layer effects between 2.8 × 10⁸ and 3.5 × 10⁸ years. Both curves are nearly identical, which shows that despite their local action as represented in Fig. 6, IGW’s stable critical layers do not disturb the dynamical evolution of surface velocity in the case of the studied star. This can be easily understood since the radiative damping becomes mostly efficient around the positions of critical layers because of its dependance on σ^-4. Moreover, if in Talon & Charbonnel (2005) the IGW’s action and the shear-induced turbulence have been added as uncoupled physical mechanisms, it has been demonstrated that in the upper region (r > r_c) where IGW are propagative, the coefficient ν_V is negligible, while in the inner one (r < r_c) the differential rotation has been damped, leading to the same result with a transport dominated by the meridional circulation.

This clearly indicates that unstable critical layers will lead to major modification of the transport of angular momentum in stellar interiors. A systematic exploration of different types of stars for different evolutionary stages will be undertaken in the near future. Moreover, in order to give quantitative information, it is necessary to improve the way waves are excited in this model. Future work will implement a prescription about penetrative convection processes. The best way to do this is to use a realistic numerical simulation of such a mechanism to obtain the excitation spectrum at the base of the convective zone (Rogers & Glatzmaier 2005; Brun et al. 2011; Alvan et al. 2012).

7. Conclusion

In this paper, we study in detail a new mechanism of interaction between IGW and the shear of the mean flow that occurs at corotation layers in stably stratified stellar radiation zones. Taking advantage of the work in the literature about atmospheric and oceanic fluids, we highlight the similarities with such stellar regions and propose an analytical approach adapted to the related case of deep spherical shells. The use of spherical coordinates brings differences in the equations and make their resolution more complicated, but the final results are comparable. We then demonstrate the intrinsic couplings between IGW and the shear-induced instabilities and turbulence that can thus not be added linearly, as done previously in stellar evolution literature. We highlight the existence of two regimes where the interactions between IGW and the shear at critical layers are strongly different:

• In the first case, the fluid is stable and the amplitude of IGW isoverdamped by the critical layer compared to the classical case,where only radiative and viscous dampings are taken intoaccount.

• In the second case, the fluid is unstable and turbulent and the critical layer acts as a secondary excitation region. Indeed, through over-reflection/transmission (when  |R|  > 1 and  |T|  > 1), energy is taken from the unstable shear that increases the amplitude of an incident IGW. Moreover, even in the case of simple reflection and transmission where  |R|  < 1 and  |T|  < 1, this demonstrates the existence of IGW tunneling through unstable regions as identified by Sutherland & Yewchuk (2004) in laboratory experiments and by Brown & Sutherland (2007), Nault & Sutherland (2007) in geophysics.

Therefore, these mechanisms opens up a new field of investigations concerning angular momentum tranport processes by IGW in stellar interiors.

Indeed, even if according to our first evolutionary calculations with STAREVOL, only the first stable regime exists in solar-type stars, we expect to find stars where the unstable regime and possible tunneling or over-reflection/transmission take place. Moreover, while the formalism presented in this work is general, several uncertainties remain. First, concerning the stochastic excitation of IGW by convection, the model used in our evolutionary code certainly underestimates the wave flux since it considers only the volumetric excitation in the bulk of the convective envelope, while convective penetration should also be taken into account. This will influence the measured action of critical layers since it is proportional to the initial amplitude of the IGW. Then, only retrograde waves are simulated here, considering that prograde ones are immediately damped and do not penetrate deeply in the radiation zone. This should normally not affect our results because the critical layers we detected are located in the deep radiation zone, but formal equations take both types of waves into account.

The last point to bear in mind is that no latitudinal dependence on the angular velocity is considered here. We explained the reason for this choice in the introduction, but one must not forget this approximation. Finally, since our goal is to obtain a complete and coherent picture of the transport of angular momentum in stellar radiation zones for every stellar type or evolutionary stage, it will be important to extend this work to cases of gravito-inertial waves, where the action of Coriolis and centrifugal accelerations is considered (e.g. Lee & Saio 1997; Dintrans & Rieutord 2000; Mathis 2009; Ballot et al. 2010) and magneto-gravito-inertial waves, when radiation zones are magnetized (e.g. Rudraiah & Venkatachalappa 1972; Kim & MacGregor 2003; MacGregor & Rogers 2011; Mathis & de Brye 2012).

Acknowledgments

The authors thank the referee for her/his comments that allowed them to improve the paper. This work was supported by the French Programme National de Physique Stellaire (PNPS) of CNRS/INSU, by the CNES-SOHO/GOLF grant and asteroseismology support in CEA-Saclay, by the Campus Spatial of the University Paris-Diderot, and by the TOUPIES project funded by the French National Agency for Research (ANR). S.M. and L.A. are grateful to Geneva Observatory, where part of this work has been achieved. T.D. acknowledges financial support from the Swiss National Science Foundation (FNS) and from ESF-Eurogenesis.

References

Appendix A: Validity of the WKBJ approximation

The form of the equation to solve is (A.1)The first step is to introduce the Liouville transformation (e.g. Olver 1974), i.e., (A.2)We deduce and Eq. (A.1) becomes (A.5)where

In the present case we have (A.6)The WKBJ approximation is available when N² ≫ σ². Consequently, we get

and (A.10)where Ri is the Richardson number of the fluid defined in Eq. (3).

Finally, we determine that the condition for applying the WKBJ approximation is |Φ| ≪ 1, which leads to .

Appendix B: Mathematical treatment for the non-perfect fluid case

We introduce (B.1)where a and b are the limits of a domain which will be defined later. The equation of propagation can then be written as (B.2)and (B.3)The new variable u, defined by , transforms the original equation into (B.4)Then, we introduce s = Dt³ where D ∈ C ^∗ : (B.5)Finally, leads to (B.6)To obtain the final equation, we define and we get the following Whittaker equation: (B.10)

All Tables

All Figures

	Fig. 1 Attenuation rate Att =  of the wave passing through a critical layer as a function of the Richardson number. We observe that Att increases with
In the text

	Fig. 2 Assumed neighbourhood of an unstable critical layer for the calculation of the IGW’s reflection and transmission coefficients. We assume that the unstable region around the critical layer (in red) has a thickness given by 2δ (zone II). The surrounding regions where IGW are propagative are in beige (zones I and III).
In the text

	Fig. 3 Level lines of reflection  |R|  and transmission  |T|  coefficients of an IGW at a critical layer as a function of the Richardson number Ric and the thickness δ of the critical layer (zone II in Fig. 2). The top panels represent the level lines for an arbitrarily chosen value of (l,m) = (4, ± 3), while the bottom panels show level lines  |R|  = 1 and  |T|  = 1 for different couples (l,m).
In the text

	Fig. 4 Left: curves (,,) defining the basis (U1,U2,U3,V1,V2,V3). Right: curves (,,) defining the basis (u1,u2,u3,v1,v2,v3).
In the text

	Fig. 5 Two configurations in a low-mass star with an external convective envelope. Left: the case of a stable critical layer (CL) where IGW are damped. Right: the case of an unstable critical layer where IGW can be over-reflected/transmitted.
In the text

	Fig. 6 Position of the critical layers for retrograde IGW for three different frequencies, superimposed with the rotational velocity profile in the star.
In the text

	Fig. 7 Left: luminosity of the waves at the location of their excitation as function of l and m. Right: ratio between τCL and τtot = τCL + τrad (see Eq. (92)).
In the text

Fig. 8 Evolution of the rotation profile, where the role of the critical layers in the transport of angular momentum is taken into account. The curves are labelled according to the corresponding ages in Gyr. Parameters are indicated above where K is the braking constant. Right: comparison between the evolution of the rotation with (red) and without (black) taking into account the effect of the critical layers.

In the text