A&A 421, 305-322 (2004)
DOI: 10.1051/0004-6361:20040121

Eigenoscillations of the differentially rotating Sun

II. Generalization of the Laplace tidal equation

N. S. Dzhalilov^1,2,3 - J. Staude²

1 - Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation of the Russian Academy of Sciences, Troitsk City, Moscow Region, 142190 Russia
2 - Astrophysikalisches Institut Potsdam, Sonnenobservatorium Einsteinturm, 14473 Potsdam, Germany
3 - Shamakhi Astrophysical Observatory of the Azerbaidjan Academy of Sciences, Baku 370096, Azerbaidjan Rep.

Received 15 May 2002 / Accepted 13 October 2003

Abstract
The general partial differential equation governing linear adiabatic nonradial oscillations in a spherical, differentially and slowly rotating non-magnetic star is derived. This equation describes mainly low-frequency and high-degree g-modes, convective g-modes, rotational Rossby-like vorticity modes, and their mutual interaction for arbitrarily given radial and latitudinal gradients of the rotation rate. Applying to this equation the "traditional approximation'' of geophysics results in a separation into radial- and angular-dependent parts of the physical variables, each of which is described by an ordinary differential equation. The angular parts of the eigenfunctions are described by the Laplace tidal equation generalized here to take into account differential rotation. It is shown that there appears a critical latitude in the sphere where the frequencies of eigenmodes coincide with the frequencies of inertial modes. The resonant transformation of the modes into the inertial waves acts as a resonant damping mechanism of the modes. Physically this mechanism is akin to the Alfvén resonance damping mechanism for MHD waves. It applies even if the rotation is rigid. The exact solutions of the Laplace equation for low frequencies and rigid rotation are obtained. The eigenfunctions are expressed by Jacobi polynomials which are polynomials of higher order than the Legendre polynomials for spherical harmonics. In this ideal case there exists only a retrograde wave spectrum. The modes are subdivided into two branches: fast and slow modes. The long fast waves carry energy opposite to the rotation direction, while the shorter slow-mode group velocity is in the azimuthal plane along the direction of rotation. It is shown that the slow modes are concentrated around the equator, while the fast modes are concentrated around the poles. The band of latitude where the mode energy is concentrated is narrow, and the spatial location of these band depends on the wave numbers (l, m).

Key words: hydrodynamics - Sun: activity - Sun: interior - Sun: oscillations - Sun: rotation - stars: oscillations

1 Introduction

In a recent paper Dzhalilov et al. (2002; Paper I) investigated which lowest-frequency eigenoscillations can occur in the real Sun and what role they might play in redistributing angular momentum and causing solar activity. We found that such waves could only be Rossby-like vorticity modes induced by differential rotation. However, the general nonradial pulsation theory adopted from stellar rotation has some difficulties for such low-frequency oscillations. For slow rotation, when the sphericity of the star is not violated seriously, the high-frequency spherical p- and g-modes are no longer degenerate with respect to the azimuthal number m (Unno et al. 1989). Independent of the spherical modes, non-rotating toroidal flows (called "trivial'' modes with a zero frequency) become quasi-toroidal with rotation (called r-modes with a nonzero frequency; Ledoux 1951; Papaloizou & Pringle 1978; Provost et al. 1981; Smeyers et al. 1981; Wolff 1998). Although rotation removes the degeneracy of the modes, it also couples the modes with the same azimuthal order, and this makes the problem more difficult. For the high-frequency modes (Rossby number $\varepsilon_{{\mbox{\tiny R}}}=\omega/2\Omega$ $\ge$ 1, where $\omega$ and $\Omega$ are the angular frequencies of oscillations and of stellar rotation, respectively) this difficulty is resolved more or less successfully. For this case the small perturbation rotation theory is applied: the eigenfunctions are represented by power series, the angular parts of which are expressed by spherical harmonic functions Y^m_l (Unno et al. 1989). These power series are well truncated, unless $\varepsilon_{{\mbox{\tiny R}}}<1$, when the role of the Coriolis force increases.

Most pulsating stars point out the low-frequency instability (Cox 1980; Unno et al. 1989). Rotation causes strong coupling of the high-order and the convective g-modes, and the r-modes with $\varepsilon_{{\mbox{\tiny R}}}<1$ and with the same m, but different l (Lee & Saio 1986). Generally the matrix of the coupling coefficients to be determined is singular (e.g., Townsend 1997). In all papers on the eigenvalue problem of nonradially pulsating stars, there exists a "truncation problem'' for the serial eigenfunctions, the angular parts of which are represented by spherical harmonics (e.g., Lee & Saio 1997; Clement 1998).

The governing partial differential equations (PDEs) of the eigenoscillations of rotating stars are complicated from the point of view of mathematical treatment, even if the motions are adiabatic. This difficulty arises because in a spherical geometry an eigenvalue problem with a singular boundary condition has to be solved. The singularity of the governing equations on the rotation axis is not a result of the applied physical approximations (small amplitude linearization, adiabaticity of the processes, incompressibility of motion, rigid rotation approximation, etc.), but is a property of the spherical geometry, permitting special classes of motion which may become the global eigenmotion of the star. Due to these difficulties it is very important to obtain some analytical solution of the eigenvalue problem with rotation, even if rather special approximations are used.

For $\omega<\Omega$ the Coriolis force is the dominant term in the equation of motion. In this case this force causes the motion to be preferentially horizontal, $v_r<v_{\theta,\phi}$. With decreasing wave frequency, $\omega\ll\Omega$, the flows become practically horizontal, $v_r\ll v_{\theta,\phi}$. Here v_r is the radial component and  $v_{\theta,\phi}$are the surface components of the velocity. This is why in many investigations of global linear or nonlinear motions in differently rotating astrophysical objects the radial velocity component and the radial changes are ignored, i.e., v_r=0, $\frac{\partial}{\partial r}=0$. That is, 2D flows are considered, e.g. Gilman & Fox (1997), Levin & Ushomirsky (2001), Cally (2001), Charbonneu et al. (1999). Nevertheless, even such a 2D approximation describes the basic physics of the global r-modes of rotating stars, e.g., Levin & Ushomirsky (2001).

However, real situations (many stars show nonradial pulsations with rather low frequencies, e.g. Unno et al. 1989) require the inclusion of $v_r\ne 0$. It is important to study, e.g., the g-mode rotation interaction, the transport of rotation angular momentum, the mixing problems of the stellar interior, etc. The consideration of a small nonzero radial component of velocity is the next important step for improving the 2D approximation. Essentially, in geophysics the traditional approximation takes into account $v_r\ll v_{\theta,\phi}$ in the governing equations. This approximation is also often applied to global motions in astrophysical objects, for example, Bildstein et al. (1996), Lou (2000), Berthomieu et al. (1978), Lee & Saio (1997), Dziembowski & Kosovichev (1987). To avoid misunderstanding, we will give here a short excursion into this approximation; see, for example, Unno et al. (1989), Shore (1992). To show the meaning of the traditional approximation clearly let us use a local analysis procedure.

In the case of a homogeneous rotating fluid, $\Omega$ is included into the wave equation only in the term

$$\begin{eqnarray*}(\vec{\Omega}\cdot\nabla)^2=\Omega_y^2\frac{\partial^2}{\partia... ...partial y\partial z}+ \Omega_z^2\frac{\partial^2}{\partial z^2}, \end{eqnarray*}$$

where $\vec{\Omega}=\{0,\Omega_y,\Omega_z\}= \{0,\Omega\sin\theta,\Omega\cos\theta\}$, the $\vec z$ axis is along the radius (local vertical), and $\vec y$ and $\vec x$ are in the meridional and the azimuthal directions, respectively. Here we can neglect $\Omega_y$, if the condition $\vert\Omega_y\frac{\partial}{\partial y}\vert \ll \vert\Omega_z\frac{\partial}{\partial z}\vert$ is satisfied. This is the traditional approximation. For a wave motion let us use $\frac{\partial}{\partial y}=ik_y$ and $\frac{\partial}{\partial z}=ik_z$. Then the condition for the traditional approximation becomes

$$\begin{eqnarray*}\vert k_y\vert\ll k_z\cot\theta \end{eqnarray*}$$

or for the horizontal $\lambda_y$ and vertical $\lambda_z$ scales

$$\begin{eqnarray*}\frac{\sin\theta}{\cos\theta}\lambda_z \ll \lambda_y. \end{eqnarray*}$$

For a stratified star the situation is the same. In the short-wave and low-frequency limit we have a dispersion relation for the rotational-gravity waves:

$$\begin{eqnarray*}\omega^2 \approx \frac{N^2k_{\perp}^2+(2\vec{\Omega}\cdot\vec{k})^2}{k^2}, \end{eqnarray*}$$

where N is the Brunt-Väisälä frequency, $\vec{k}=k_r\vec{e_r}+\vec{k}_{\perp}$, and $k=\vert\vec{k}\vert$. Because the waves are nearly incompressible, $\vec{k}\cdot\vec{v}= k_rv_r+k_{\perp}v_{\perp}\approx 0$. Here $\vec{v}=\{v_r,v_{\theta},v_{\phi}\}= \{v_r,v_{\perp} \}$. From here

$$\begin{eqnarray*}\frac{v_r}{v_{\perp}}\approx -\frac{k_{\perp}}{k_r}\cdot \end{eqnarray*}$$

If |N²| is much larger than $\omega^2$ and $4\Omega^2$ (this is valid for most stars), we get $k_{\perp}^2/k^2\sim 0$ from the dispersion relation, which means that the horizontal component $k_{\perp}$ is much smaller than the radial one, k_r. This is the same situation as in the homogeneous case. From this we get $v_{\perp}\gg v_r$. Because the wave vector $\vec k$ is almost radial, we have $(2\vec{\Omega}\cdot\vec{k})\sim 2\Omega_rk_r$. In this way we conclude that only the radial component of  $\vec{\Omega}$ is important for determining the frequency $\omega$. This argument suggests that the neglect of  $\Omega_{\theta}$ in the governing equations should hardly alter the eigenfrequencies. This limitation, widely used in geophysical hydrodynamics (e.g. Eckart 1960), is called the "traditional approximation'' and has been used first by Laplace (1778) to study tidal waves (Lindzen & Chapman 1969). Laplace's equation (or the traditional approximation) for $\varepsilon_{{\mbox{\tiny R}}}<1$ is applicable to the stellar case too. The main advantage of this approximation is that it decomposes the initial system of equations into a set of ordinary differential equations (ODEs) (e.g., Lindzen & Chapman 1969; Berthomieu et al. 1978; Bildsten et al. 1996; Lee & Saio 1997). The angular parts of the eigenfunctions are described by Laplace's tidal equation. Solving this equation numerically by using a relaxation method, Lee & Saio (1997) first avoided the representation of the solutions by $Y_l^m(\cos\theta)$ functions for the $\Omega=$ const. case, and they had no problem with the truncation of the series.

In the present work for the non-magnetic and non-convective cases we obtain one PDE in a spherical geometry for the adiabatic pressure oscillations in the differentially rotating star ( $\Omega=\Omega(r, \theta)$) with arbitrary spatial gradients of rotation (Sect. 2). This general equation is split into $\theta$- and r-component ODEs, if the traditional approximation is applied (Sect. 3). The $\theta$-component equation is Laplace's tidal equation generalized for the differentially rotating case. In Sect. 4 we analyse this equation more qualitatively. We find the general condition for the shear instability due to differential rotation in latitude. We find that the smallest rotation gradient is responsible for the prograde (seen in the rotating frame) vorticity wave instability, while a stronger gradient causes the retrograde wave instability. For solar data (small rotation gradients) the m=1 prograde mode instability is possible (Sect. 4.4). The possible existence of such a global horizontal shear instability on the Sun has been investigated by Watson (1981) and Gilman & Fox (1997), that of shear and other dynamic instabilities and of thermal-type instabilities in stars also by Knobloch & Spruit (1982) and others. Laplace's tidal equation for low frequencies in the rigid-rotation case is investigated in detail in Sect. 5. It is shown that the eigenfunctions are defined by Jacobi polynomials which are of higher order than the Legendre polynomials.

2 Basic equations

The motion of the fluid in a self-gravitating star, neglecting magnetic field and viscosity, may be described in an inertial frame by the hydrodynamic equations. These equations in conventional definitions are written as

    

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{V})=0,$$ (1)

$$\rho ~\left(\frac{\partial}{\partial t} + \vec{V} \cdot \nabla \right) ~\vec{V}= - \nabla p - \rho \nabla \Phi,$$ (2)

$$\rho~T\left(\frac{{\partial}{}}{{\partial}{t}}+\vec{V} \cdot \nabla \right)s= \rho ~\varepsilon_{{\mbox{\tiny N}}} - \nabla \cdot \vec{F},$$ (3)

$$\nabla^2 ~\Phi= 4\pi G~\rho.$$ (4)

2.1 Equilibrium state

We suppose that the equilibrium state (zero indices of the variables) of the star is stationary and that its differential rotation is axially symmetric:

$$% \vec{V}_0(r,\theta) = \vec{\Omega}\times\vec{r} = \Omega ~r ~\sin\theta ~\vec{e}_{\phi},$$ (5)

where $\vec{\Omega}(r,\theta) = \Omega_r ~\vec{e}_r + \Omega_{\theta} ~ \vec{e}_{\theta}$, with the components $\Omega_r = \Omega\cos\theta, \ \ \Omega_{\theta} = - \Omega\sin\theta,$ and $\Omega_{\phi} \equiv 0$, is the stellar angular velocity of rotation described in spherical polar coordinates, $(r, \theta, \phi)$. Here $\vec{e}_{i}$ with $i=r, \theta,$ or $\phi$ are the unit vectors. We will not include convective motion and meridional flows into the initial steady state. In that case we may obtain, in particular from the Eq. (2) of motion, the hydrostatic equilibrium relation

 

$$% \frac{\nabla ~p_0}{\rho_0} \nabla\left(-\Phi_0 + \frac{1}{2} ~\vert\vec{\Omega}\times \vec{r}\vert^2 \right) - \Omega ~r^2 \sin^2 \theta ~ \nabla \Omega$$ =

$$\vec{g} - \vec{g}_{\mbox{\tiny R}}=\vec{\tilde{g}},$$ = (6)

where $\vec{g}=-\nabla\Phi_0$, $\vec{g}_{\mbox{\tiny R}}=\vec\Omega\times\vec\Omega \times\vec{r}$, and $\vec{\tilde{g}}=\{\tilde{g}_r,\tilde{g}_\theta,0\}$. It follows that the effective gravity  $\vec{\hbox{\~{g}}}$ cannot be a potential field if differential rotation $\nabla\Omega\ne 0$ is present. This is important for rapidly rotating stars where the configuration is deformed by the centrifugal force as well as by differential rotation. For slowly rotating stars (like the Sun) we may assume that the initial state is only marginally disturbed by rotation, and $\vec{\tilde{g}} \approx \vec{g} \approx -g(r)\vec{e_r}$ can be applied. That is, non-sphericity is not essential for the generation of waves (Unno et al. 1989).

2.2 Equations of oscillation

Small amplitude deviations from the basic state of the star may be investigated by linearizing Eqs. (1)-(4). For Eulerian perturbations (variables with a prime) the equation of motion becomes (Unno et al. 1989)

 

$$% \vec{e}_i\frac{\partial}{\partial t^{\prime}} ~V_i + 2~ \vec{\O... ...{\prime}}{\rho_0} - \frac{1}{\rho_0} ~\nabla p^{\prime} - \nabla \Phi^{\prime}.$$ (7)

Here the operator

$$% \frac{\partial}{\partial t^{\prime}} = \frac{\partial}{\partial t} + \Omega ~\frac{\partial}{\partial \phi}$$ (8)

represents the temporal derivative referring to a local frame rotating with an angular velocity $\Omega=\Omega(r, \theta)$. For low-frequency waves Saio (1982) has shown numerically that the Cowling approximation is good enough in most cases. Thus we will neglect the perturbation of the gravitational potential in Eq. (7), $\Phi^{\prime}=0$. We are interested in very slow motions such that $v_{\rm ph} \ll c_{\rm s}$, where  $v_{\rm ph}$ is the phase velocity of the waves and $c_{\rm s}$ is the sound speed. Then the incompressible fluid motion limit, $c_{\rm s}^2\to\infty$(it is within the adiabatic approximation) may be applied, and instead of Eq. (1) we use

$$% \nabla\cdot\vec{V} = \frac{1}{r^2}\frac{{\partial}{\left(r^... ...\sin\theta}\frac{{\partial}{V_{\phi}}}{{\partial}{\phi}} = 0 .$$ (9)

We have shown in Paper 1, that nonadiabatic effects are of great importance for the dynamics of low-frequency rotation modes. However, here we shall restrict ourselves to the adiabatic case only because the mathematical treatment of wave equations in spherical geometry is rather difficult. For adiabatic waves we obtain from Eq. (3)

 

$$\frac{\partial \rho^{\prime}}{\partial t^{\prime}} + \rho_0 ~\vec... ...ec{A} = \frac{1}{c^2_{\rm s}} ~\frac{\partial p^{\prime}}{\partial t^{\prime}},$$ (10)

$$\vec{A} = \frac{\nabla\rho_0}{\rho_0}-\frac{1}{\Gamma_1} \frac{\nabla p_0}{p_0}\cdot$$ (11)

In the incompressible limit ( $c_{\rm s}^2\to\infty$) Eq. (10) reads

$$% \frac{\partial \rho^{\prime}}{\partial t^{\prime}} + \rho_0 ~\vec{V}\cdot \vec{A} =0 .$$ (12)

Thus we have a complete set of Eqs. (7), (9), (12) to describe adiabatic low-frequency non-radial oscillations in a differentially rotating star. For our axisymmetric stationary initial state we may represent all the perturbed variables $\vec V$, p', and $\rho'$ in the inertial frame as

$$% \vec{V}(r, \theta, \phi; t) \Longrightarrow \vec{V} ~(r, \theta) ~ {\rm e}^{i(m \phi-\omega_0 t)}.$$ (13)

Considering ${{\partial}{}}/{{\partial}{t'}}=-i\omega$ and ${{\partial}{}}/{{\partial}{\phi}}=im$ we find from Eq. (8) the relation between the frequencies in the inertial and rotating frames

$$% \omega = \omega_0 - m~\Omega~(r, \theta),$$ (14)

from which we get $\nabla\omega=-m\nabla\Omega$. If we separate the variable part of the rotation frequency, $\Omega(r,\theta)=\Omega_{\odot}+ \tilde\Omega(r,\theta)$, then $\omega=\omega_0 - m\Omega_{\odot} - m\tilde\Omega = \omega_{\odot} - m\tilde\Omega =\omega_{\odot} - k_{\phi}v_{\rm ph}$ ($k_{\phi}$ and $v_{\rm ph}$ are the local azimuthal wave number and the phase velocity, respectively). We will study the case $\tilde\Omega \ll \Omega_{\odot}$. From now on $\Omega\approx \Omega_{\odot}$ and $\nabla\Omega = \nabla\tilde\Omega\ne 0$ will be used. Low frequencies seen in the rotating frame mean that we are close to the resonant frequencies in the inertial frame ( $\omega_{\odot}\approx k_{\phi}v_{\rm ph}$).

Excluding $\rho'$ from Eq. (7) by using Eq. (12) and taking its r and $\theta$ components we get

  

$$V_r - \frac{2\Omega_\theta}{i\omega}V_{\phi} - \frac{\tilde{g}_r}... ...vec{V}\cdot\vec{A}) = \frac{1}{\rho_0 i \omega}~\frac{\partial p'}{\partial r},$$ (15)

$$V_{\theta} + \frac{2\Omega_r}{i\omega} V_{\phi} - \frac{\tilde{g}... ...= \frac{1}{r \rho_0 i \omega} ~\frac{\partial p^{\prime}}{\partial \theta}\cdot$$ (16)

The equations may be simplified for further analysis if we consider the limit of slow rotation: $\Delta=\Omega^2/(GM/R^3)\ll 1$. For the Sun $\Delta\approx 10^{-4.7}$. Then $g\gg g_{\mbox{\tiny R}}$ where $g= GM(r)/r^2 +O(\Delta)$, $p_0= p_0(r)+O(\Delta)$, $\rho_0=\rho_0(r)+O(\Delta)$ and

$$% \frac{\tilde{g}_\theta}{\tilde{g}_r}=-\frac{\Omega^2r\sin\t... ...{g-\Omega^2r\sin^2\theta}\approx -\Delta\sin\theta\cos\theta.$$ (17)

For slow rotation $\vec{V}\cdot\vec{A}\approx V_rA_r$ and considering that

$$% N^2 = g \left( \frac{1}{\Gamma_1} ~\frac{1}{p_0} ~\frac{{\rm d} p_0}{{\rm d} r} - \frac{1}{\rho_0} ~\frac{{\rm d} \rho_0}{{\rm d} r} \right)$$ (18)

is the squared Brunt-Väisälä frequency, Eqs. (15) and (16) become

  

$$\left(1-\frac{N^2}{\omega^2}\right)V_r - \frac{2\Omega_\theta}{i\omega}V_{\phi} = \frac{1}{\rho_0 i \omega}~\frac{\partial p'}{\partial r},$$ (19)

$$V_{\theta} + \frac{2\Omega_r}{i\omega} V_{\phi} +\Delta\sin\theta... ..._r = \frac{1}{r \rho_0 i\omega}~\frac{\partial p^{\prime}}{\partial\theta}\cdot$$ (20)

In Eq. (20) the term with $\Delta$ is due to the centrifugal force. However, this term is negligible for any low frequency $\omega\ll\Omega$, as $\Delta\ll 1$. For $\omega^2\sim N^2$ the $\Delta$ term may be omitted because $V_\theta\gg \Delta V_r$. For the other frequencies we may exclude V_r in Eq. (20) using Eq. (19). Insofar as $(N^2/\omega^2)/(1-N^2/\omega^2)\approx O(1)$ the $\Delta$ term will be negligible in comparision with the other terms of Eq. (20).

Now adding to these equations the $\phi$ component of Eqs. (7) and (9) we get our set of equations:

    

$$\left(1 - \frac{N^2}{\omega^2} \right) V_r - j \frac{2\Omega_\theta} {i\omega}V_{\phi} = \frac{1}{\rho_0 i \omega}~\frac{\partial p'}{\partial r},$$ (21)

$$V_{\theta} + \frac{2\Omega_r}{i\omega} V_{\phi}= \frac{1}{r \rho_0 i\omega} ~\frac{\partial p^{\prime}}{\partial \theta},$$ (22)

$$-i\omega V_{\phi} + 2(\Omega_r V_\theta- j \Omega_\theta V_r)$$

$$\hspace*{4mm} + \sin\theta \left(j V_r r\frac{\partial\Omega}{\pa... ...{\partial\Omega}{\partial\theta} \right) = - \frac{i m~ p'}{r\rho_0\sin\theta},$$ (23)

$$i m V_{\phi} + \frac{\sin \theta}{r} ~\frac{\partial}{\partial r}... ...r^2 V_r\right) + \frac{\partial}{\partial \theta} (\sin\theta~ V_{\theta}) = 0.$$ (24)

Here we have introduced the formal parameter j to switch to the traditional approximation: in the general case $j\equiv 1$, and to make this switch we put j=0, which is valid if

  

$$V_r \ll \cot\theta~V_\theta ,\;\;\;\; r\Omega'(r) V_r \ll \Omega'(\theta)~V_\theta,$$ (25)

$$\frac{N^2}{2\omega\Omega} V_r \gg V_\phi.$$ (26)

These conditions practically agree with those discussed in the Introduction. In most cases the second condition in Eq. (25) should be omitted. In this way the traditional approximation demands a small, but nonzero radial velocity ($V_r\ne 0$) for the low-frequency ( $\omega<2\Omega$) waves. For the solar situation $N^2/(2\Omega\omega)>10^6$, except for the very narrow regions where N²(r) is going very sharply to zero (center and bottom of the convective zone). The first condition of Eq. (25) is violated very close to the equator plane. However, this is no obstacle to calculating the eigenoscillations of the hemisphere.

For the general case the system of Eqs. (21)-(24) has been reduced in Appendix A to one PDE for the pressure perturbation:

$$% \Big[\psi_1 \check{\delta}^2_r + \psi_2 \check{\delta}^2_{\... ...k{\delta}_r \check{\delta}_{\mu} + \psi_6 \Big] \tilde{P} = 0,$$ (27)

where $p^{\prime} = \mu r \rho_0 i \omega \tilde{P}$ and $\mu=\sin\theta$ is a new independent variable. The coefficients $\psi_{1-6}(r,\theta)$ of Eq. (27) are rather complicated and we present them in Appendix A. The operators are defined as

$$% \check{\delta}_r =r\frac{{\partial}{}}{{\partial}{r}} , \ \... ...eck{\delta}^2_{\mu} =\check{\delta}_{\mu}\check{\delta}_{\mu}.$$ (28)

Equation (27) is the main equation for nonradial rotational-gravity waves in a differentially rotating star. In future papers of the present series this equation will be the basis of further studies. Here we will restrict ourselves to the traditional approximation.

3 The "traditional approximation'': $j\equiv 0$

Let

$$\varepsilon_{{\mbox{\tiny R}}} = \frac{\omega}{2 \Omega} , \ \bet... ... r} , \ \beta_{\mu} = \frac{\mu}{2\Omega} \frac{\partial \Omega}{\partial \mu},$$ (29)

$$\alpha = \varepsilon^2_{{\mbox{\tiny R}}} - \left(1 - \mu^2\right)\left(1+ \beta_{\mu}\right), \ \mu=\sin\theta,$$

where $\varepsilon_{{\mbox{\tiny R}}}$ is the Rossby number (we are interested in $\varepsilon_{{\mbox{\tiny R}}}\le 1$), and  $\beta _{\mu }$ and $\beta_r$ are the latitudinal and radial logarithmic gradients of the rotation rate, respectively. Strictly speaking, the condition j=0 in Eq. (27) is not applicable in two cases: (1) at  $\omega^2=N^2(r)$, i.e. in the turning points in the radial direction; (2) at $\alpha=0$ or $\varepsilon^2_{{\mbox{\tiny R}}}\approx\cos^2\theta$, i.e. in the latitudinal turning point. The last one is more important, because the traditional approximation filters out such important and interesting phenomena as the trapping of Rossby-like waves around the equator. In geophysics this phenomenon is investigated separately (Pedlosky 1982; Gill 1982). The applicability of the traditional approximation for rigid rotation has been checked by numerical modelling as well as by experimental verification. For minor differential rotation of the star (small $\beta_r$and  $\beta _{\mu }$) these examinations are also valid. Thus for j=0Eq. (27) becomes:

 

$$% \frac{\varepsilon_{{\mbox{\tiny R}}}}{\alpha^2 \mu^2} \left[ \a... ...l^2}{\partial r^2} + b_4 r \frac{\partial}{\partial r} + b_5 \right) p^{\prime}$$ (30)

with the parameters

$$q_3 = \varepsilon_{{\mbox{\tiny R}}} \left(1 - \mu^2\right) [\alpha (1-2 a_0) - q_1] ,$$

$$q_4 = \varepsilon_{{\mbox{\tiny R}}}\left(1-\mu^2\right) \left[\alpha\left(a_0^2-\check{\delta}_{\mu} a_0\right) +a_0 q_1\right]-q_2 ,$$

$$a_0 = 1 - \beta_{\mu} \frac{m}{\varepsilon_{{\mbox{\tiny R}}}} , ... ...tiny R}}}}\left( 1+\beta_{\mu}\frac{m}{\varepsilon_{{\mbox{\tiny R}}}}\right) ,$$

$$q_1 = 2 - \mu^2 - \varepsilon^2_{{\mbox{\tiny R}}} \frac{2 - 3\mu^2}{1 - \mu^2}$$

$$\hspace*{7mm} + \beta_{\mu} \left[ 2 \mu^2 - 1 - 2 \left(1 - \mu^... ... \beta_{\mu}\right) \left( m/\varepsilon_{{\mbox{\tiny R}}} +1 \right) \right],$$

$$q_2 = h_0 + h_1 \beta_{\mu} + h_2 \beta^2_{\mu},$$

$$h_0 = \varepsilon_{{\mbox{\tiny R}}} \left(m^2 - 1\right)\left(\varepsilon^2_{{\mbox{\tiny R}}} + \mu^2 - 1\right)$$

$$\hspace*{7mm} - \mu^2 \left[m \left(1 - \mu^2\right) + \varepsilon^2_{{\mbox{\tiny R}}} \left(m - 2\varepsilon_{{\mbox{\tiny R}}}\right)\right],$$

$$h_1 = \left(1 - \mu^2\right) \left[ - 2 m \mu^2 + \varepsilon_{{\mbox{\tiny R}}} \left(m^2 - 4 + 5 \mu^2\right) \right.$$

$$\hspace*{7mm} + \left. 4 m \varepsilon^2_{{\mbox{\tiny R}}} \right] - m \varepsilon^2_{{\mbox{\tiny R}}},$$

$$h_2 = \left(1 - \mu^2\right)^2 \left(3m - 2 \varepsilon_{{\mbox{\tiny R}}} + 2 m \beta_{\mu}\right)$$

$$\hspace*{7mm} -\left(1 - \mu^2\right) m \left(1 - 2 m \varepsilon_{{\mbox{\tiny R}}}\right),$$ b₄ = 1-2b₁ + b₂ ,

$$b_5 = b_1^2-\check{\delta}_r b_1 - b_1 b_2 + b_3,$$

$$\psi = 1 - \frac{N^2}{\omega^2} \ , \ \check{\delta}_r \psi = - 2... ...\mbox{\tiny R}}}} \beta_r - \frac{r}{\omega^2} ~ \frac{{\rm d} N^2}{{\rm d} r},$$

$$b_1 = 1 + \symbol{''1A}_{\rho} - \frac{m}{\varepsilon_{{\mbox{\ti... ...\;\; \symbol{''1A}_{\rho} = \frac{r}{\rho_0} ~\frac{{\rm d} \rho_0}{{\rm d} r},$$

$$b_2 = 3 + \symbol{''1A}_{\rho} - \beta_r \frac{m}{\varepsilon_{{\mbox{\tiny R}}}} - \frac{\check{\delta}_r \psi}{\psi},$$

$$b_3 = 2 + 3\symbol{''1A}_{\rho} - \beta_r \frac{m}{\varepsilon_{{... ...\tiny R}}}} \left( 3 + \beta_r \frac{m}{\varepsilon_{{\mbox{\tiny R}}}} \right)$$

$$\hspace*{7mm} -\left( 1 + \symbol{''1A}_{\rho} - \beta_r\frac{m}{\varepsilon_{{\mbox{\tiny R}}}} \right) \frac{\check{\delta}_r \psi}{\psi},$$

$$\check{\delta}_r b_1 = \symbol{''1A}_{\rho}-\beta_r \frac{m}{\var... ...\tiny R}}}} \left(1+\beta_r\frac{m}{\varepsilon_{{\mbox{\tiny R}}}}\right)\cdot$$

Remembering that $\tilde\Omega(r,\theta) \ll \Omega_{\odot}$, the left-hand side of Eq. (30), is a function of $\mu=\sin\theta$, while the right-hand side is a function of r only.

We assume here an extra condition for the rotation rate: $\tilde{\Omega}(r,\theta)=\Omega_1(r)+\Omega_2(\theta)$. Such a representation of the variable part of the rotation rate is no strong restriction. For the Sun, for example, $\Omega_1(r)\approx 0$ except for the solar tachocline. Depending on the problem to be solved we can simplify the calculations ignoring either $\Omega_1$ or $\Omega_2$.

In that way we may separate the variables

$$% p'(r, \mu) = \Theta(\mu)Q(r).$$ (31)

Now putting Eq. (31) into (30) we obtain the "$\theta$''- and "r''-equations:

$$% \left[\alpha\varepsilon_{{\mbox{\tiny R}}}\left(1 - \mu^2\r... ...mu^2}{\varepsilon_{{\mbox{\tiny R}}}} \right]\Theta(\mu) = 0 ,$$ (32)

$$% \left(r^2\frac{{\rm d}^2}{{\rm d}r^2}+b_4 r\frac{{\rm d}}{{\rm d}r}+b_5 -\Lambda\psi\right)Q(r)=0.$$ (33)

These two equations are connected to each other by two common spectral parameters: $\omega$, the oscillation frequency and $\Lambda$, the separation parameter. Both must be searched for as a solution of the boundary value problem. So far the logarithmic gradients of the rotation rate are arbitrary functions, $\beta_{\mu}=\beta_{\mu}(\mu)$, $\beta_r=\beta_r(r)$.

4 The generalized Laplace tidal equation

Equation (32) is the generalized Laplace equation if differential rotation is present, $\beta_{\mu} \neq 0$. For rigid rotation, $\beta_{\mu} = 0$, Eq. (32) becomes the standard Laplace equation:

 

$$% \left[\frac{1-\mu^2_{\ast}}{\varepsilon^2_{{\mbox{\tiny R}}}-\m... ...}}\right) + \frac{\Lambda}{\varepsilon^2_{{\mbox{\tiny R}}}} \right]\Theta = 0,$$ (34)

where $\mu_{\ast}=\cos\theta$. A peculiarity of this equation is the presence of three singular points: at the pole, at the equator, and between these two if $\varepsilon^2_{{\mbox{\tiny R}}}\leq 1$ (our case). Therefore it is hard to solve such an equation analytically or numerically to find the eigenvalues. Since the time of Laplace investigations in geophysics were focused on almost two-dimensional (horizontal) motions in strongly stratified fluids with $V_r\approx 0$. In this situation the r-component Eq. (33) does not appear, and one is looking for the eigenvalues $\Lambda$ in Laplace's Eq. (34) for a given $\varepsilon_{{\mbox{\tiny R}}}$ (it is expressed through the thickness of the fluid layer, e.g., in a shallow water-wave system). For the special cases and for the general case too, when the eigenfunctions are expressed through the Hough functions (on the whole these are the same infinite series of Y_l^m harmonics), see Lindzen & Chapman (1969) and references therein.

In astrophysics the r-component Eq. (33) appears, and the two equations must be solved together to find both spectral parameters. In the rigid-rotation case Lee & Saio (1997) looked for $\Lambda$ numerically in an approach similar to that in geophysics, fixing $\varepsilon_{{\mbox{\tiny R}}}$ in Eq. (34). Here we offer another approach, where we will find $\Lambda$ from the r-equation for a given $\varepsilon_{{\mbox{\tiny R}}}\leq 1$.

It is convenient to introduce into Eq. (32) the new variable $x=\mu^2=\sin^2\theta$:

$$% 4\left(1 + \beta_{\mu}\right)(1 - x)(x - a)x^2\frac{{\rm d}... ...d}x^2} - A_1 2x\frac{{\rm d}\Theta}{{\rm d}x} + A_2\Theta = 0,$$ (35)

where

  

$$A_1 = (1 - x)\left[2 - x - 3\varepsilon^2_{{\mbox{\tiny R}}} + \beta_{\mu}\left(4x - 3 - 2m\varepsilon_{{\mbox{\tiny R}}}\right)\right]$$

$$\hspace*{7mm} - 2\beta_{\mu}^2 (1 - x)^2 + \varepsilon^2_{{\mbox{\tiny R}}},$$ (36)

$$A_2 = (1 - x)\left[m^2 \left(1 - \beta_{\mu}\right) + \frac{m}{\v... ...\mbox{\tiny R}}}}\left(4x - 3 - 2\varepsilon^2_{{\mbox{\tiny R}}}\right)\right]$$

$$\hspace*{7mm} - m\varepsilon_{{\mbox{\tiny R}}}(1 - x) + m\vareps... ...da}{\varepsilon^2_{{\mbox{\tiny R}}}} x(x - a)^2\left(1 + \beta_{\mu}\right)^2,$$

$$a = 1-\frac{\varepsilon^2_{{\mbox{\tiny R}}}}{1+\beta_{\mu}}\cdot$$ (37)

This equation determines the tangential structure of the eigenfunctions, while Eq. (33) is responsible for the radial behavior. A detailed investigation of Eq. (33) is not included in the present work. A similar equation for a realistically stratified model of the Sun has been investigated by Oraevsky & Dzhalilov (1997). We recall that the radial structure of the eigenfunctions depends on the sign of $\Lambda$ (either radiative or convective modes). By the usual asymptotic methods we can find the eigenvalues $\Lambda$as a solution of Eq. (33) (which has turning points at $N^2(r=r_{1,2})=\omega^2$) fixing $\varepsilon_{{\mbox{\tiny R}}}$. Then we will have an integral dispersion relation $\Lambda^{1/2}\int_{r_1}^{r_2}\sqrt{-\Psi}{\rm d}r \sim n\pi$. Now we can already estimate an approximate value of $\Lambda$:

$$% \frac{\Lambda}{\varepsilon^2_{{\mbox{\tiny R}}}} \sim \frac... ...^2} = \frac{n^2\pi^2}{N^2_m/4\Omega^2} \sim n^2\times 10^{-6},$$ (38)

where n is the radial harmonic number, N_m is the mean value of the Brunt-Väisälä frequency in the radiative interior or in the convection zone, and  $N_m/2\Omega$ is the Prandtl number. An estimate of Eq. (38) is made for a solar model, but we think similar values of the Prandtl number are valid for most other stars too. Thus we can omit from A₂ in Eq. (36) the last term, if the radial number n is not too large. This will allow us to bring our working equation into the standard form of Heun's equation.

4.1 Fluid velocities

From Eqs. (A.2), (A.3) and (A.5) we can derive in the traditional approximation the following formulae for the components of the fluid velocity:

   

$$V_r=\frac{1}{\rho_0 i\omega\psi} \frac{{\partial}{p'}}{{\partial}{r}},$$ (39)

$$V_{\theta}= \frac{\pm 1}{r \rho_0 i\omega} \frac{\varepsilon^2_{{... ...arepsilon_{{\mbox{\tiny R}}}} - 2x\frac{{\partial}{}}{{\partial}{x}}\right) p',$$ (40)

$$V_{\phi}=\frac{1}{r \rho_0 \omega} \frac{1}{\sqrt{x}} \left[\frac... ...{{\mbox{\tiny R}}} x\frac{{\partial}{}}{{\partial}{x}} - m\right) + m\right]p'.$$ (41)

Here the different signs $\pm$ of  $V_{\theta}$ correspond to the northern and southern hemispheres, so that $\cos\theta=\pm\sqrt{1-x}$. Our further aim is to find such solutions for $p'=\Theta(x)Q(r)$ so that all the components of the velocity remain limited at the pole (x=0), at the equator (x=1), and at both turning points, where x=a and $\psi(r)=0$.

4.2 Heun's equation

Now we will impose a restriction on $2\beta_{\mu} = {{\partial}{(\ln\Omega)}}/{{\partial}{(\ln\mu)}} \approx$ const, the logarithmic latitudinal gradient of the rotation frequency. We might take the linear dependence $\beta_\mu \sim x$, but for such a profile the structure of the solutions is not changed qualitatively. Really $\beta_{\mu}=$ const. corresponds to a rotation profile

 

$$% \frac{\Omega}{\Omega_{eq}} \mu^{2\beta_{\mu}}=\left(1-\cos^2\theta\right)^{\beta_{\mu}}$$ = (42)

$$1-\beta_{\mu}\cos^2\theta+\frac{\beta_{\mu}\left(\beta_{\mu}-1\right)}{2!}\cos^4\theta - \cdots .$$ =

Here, the last expansion is not valid only at the pole $\theta=0$. Equation (42) describes rather well the solar rotation at the surface at least up to a latitude of $60^\circ$ if $\beta_\mu\approx 0.14$. The rotation gradient of the polar region is not determined from observations as exactly as that around the equator, but in any case the rotation profile given by Eq. (42) is not good enough near the pole. Thus the profile given by Eq. (42) represents our modelling of the latitudinal differential rotation which is very close to the solar surface rotation at $\beta_\mu\approx 0.14$. Helioseismological data have shown that such a rotation picture is traced downwards without strong changes up to the tachocline. In the deeper solar interior the rotation is close to that of a solid body. Here $\beta _{\mu }$ is unknown and its expected value should be much smaller, $\beta_\mu \ll 1$. The rotation law of the central region is also unknown. Thus, $\beta _{\mu }$ will be a free parameter in our task.

Let us introduce the new dependent variable

$$% \Theta = x^{\sigma}Y(x) , \ {\rm where}$$ (43)

 

$$2\sigma = \beta_{\mu}S_1 \pm\sqrt{\beta^2_{\mu}S^2_1-S_2} = 2\sigma_{1, 2},$$ (44)

$$S_1 = \frac{5+2m\varepsilon_{{\mbox{\tiny R}}} + 2\beta_{\mu}}{2\left(1+\beta_{\mu} - \varepsilon^2_{{\mbox{\tiny R}}}\right)},$$

$$S_2 = \frac{\beta_{\mu}\frac{m}{\varepsilon_{{\mbox{\tiny R}}}}\l... ...{\mbox{\tiny R}}}\right)} {1+\beta_{\mu}-\varepsilon^2_{{\mbox{\tiny R}}}}\cdot$$

Then for Y(x) we get a new equation from Eq. (35)

 

$$% x(1-x)(x-a)Y'' + \frac{1}{2}\left[4\sigma(1-x)(x-a) - \frac{A_1... ...beta_{\mu}}\right]Y' -[(x-1)\nu_0 + \varepsilon_{{\mbox{\tiny R}}}\nu_1] Y = 0,$$ (45)

with

 

$$\nu_0 = \frac{1+4\beta_{\mu}}{4\left(1+\beta_{\mu}\right)} \left(... ...\mbox{\tiny R}}}} - 2\sigma\right)+ \sigma\beta_{\mu}(S_1 - 1) - \frac{S_2}{4},$$

$$\nu_1 = \frac{m\left(m\varepsilon_{{\mbox{\tiny R}}} -1\right) + 2\sigma\varepsilon_{{\mbox{\tiny R}}}}{4\left(1+\beta_{\mu}\right)}\cdot$$ (46)

Equation (45) is Heun's equation (Heun 1889) in standard form

$$% Y''+ \left( \frac{\gamma}{x} + \frac{\delta}{x-1} + \frac{\... ...'+ \frac{\tilde{\alpha}\tilde{\beta}(x-h)}{x(x-1)(x-a)} Y = 0.$$ (47)

The Riemannian scheme for this equation is

$$\begin{eqnarray*}P\left\{ \begin{array}{ccccc} 0 & 1 & a & \infty & \\ 0 & 0 & ... ... 1-\delta & 1-\varepsilon & \tilde{\beta} & \end{array}\right\}, \end{eqnarray*}$$

where the exponents of Eq. (47) are connected by Riemann's relation

$$% \tilde{\alpha} + \tilde{\beta} -\gamma - \delta - \varepsilon + 1 = 0.$$ (48)

In the Riemann scheme the first row defines the singular points of Heun's equation, while the corresponding exponents are placed in the second and third rows. The exponents of Eq. (47) are

 

$$1-\gamma = \beta_{\mu}S_1-2\sigma , \ \ 1-\delta = \frac{1}{2} ,$$ (49)

$$1-\varepsilon = 2+\beta_{\mu}\left(1-S_1 + \frac{3/2}{1+\beta_{\mu}}\right),$$

$$2\tilde{\alpha} = S + q , \ 2\tilde{\beta} = S - q ,~~~ q = \sqrt{S^2 - 4\nu_0},$$ (50)

$$h = 1-\varepsilon_{{\mbox{\tiny R}}}\frac{\nu_1}{\nu_0} , \ \ S = 2\sigma - 2 - \beta_{\mu} + \frac{3/2}{1 + \beta_{\mu}}\cdot$$

Historically the Riemannian scheme was of great importance in the development of the theory of ordinary differential equations with singularities, see, e.g., Whittaker & Watson (1927) and Ince (1927). The essence of the Riemannian scheme is the following: if we have some ordinary second order differential equation with regular singularities (more than two), the solutions of this equation may be represented by the fundamental branches of the P-function, if the exponents of the equation satisfy the Riemannian relation. Mainly, the P-functions correspond to some power series for which the domain of convergence includes at least two singular points: a circle with one singularity at the boundary, another singularity at the centre or an elliptic limacon with two singularities as foci, and a third singularity at the circumference. These regular branches of the P-functions make it possible to find, by analytical continuation, the global solution which includes all the singularities.

In our case the physical equations are reduced to Heun's Eq. (47) for which the Riemannian relation (48) is satisfied. This means that the regular eigenfunctions in the hemisphere including the singular pole and the equator exist and we can find them. A direct numerical solution of this equation is practically impossible and we intend to use the P-branches to find the analytical eigenfunctions. However, the main difficulty of Heun's Eq. (47) is connected with the "accessory'' parameter h which is not included in the Riemannian scheme. The arbitrariness of h does not make it possible to represent the solutions by a single P-function. Therefore, hypergeometric and confluent hypergeometric equations, the differential equations of Lamé, Mathieu, Legendre, Bessel, and Weber, those of the polynomials of Jacobi, Chebyshev, Laguerre and Hermite as well as those of Bateman's k-functions are special or limiting cases of Heun's equation. Due to the parameter h the solutions of Heun's equation are usually searched for as series of P-functions. For example, Erdélyi (1942) has represented in such a series the P-branches by the hypergeometric functions. In physical tasks, e.g. the damping of MHD waves in resonant layers, the same equation arises (Dzhalilov & Zhugzhda 1990). The exact solutions of Eq. (47) will be considered in the next paper of this series. Here we will restrict ourselves to a qualitative analysis and to a limiting case.

Generally for low frequencies ( $\varepsilon_{{\mbox{\tiny R}}}<1$) three singularities of Heun's equation may be realized in a hemisphere: 0 < a < 1. Note that the singularity $x=\infty$ in the Riemann scheme does not occur in our problem. At the equator (x=1), Eq. (47) has two solutions with the exponents 0 and $1-\delta=1/2$. If we put the P-solutions with these exponents into Eqs. (39)-(41), we can see that the first solution is divergent, while the second solution provides the limited $V_{\theta}$ and $V_\phi$ at the equator. In the same manner, at the pole (x=0) one of the solutions with the exponents "0'' and ($1-\gamma$) may provide the limited $V_{\theta}$and $V_\phi$ only for some selected values of $\gamma$ (see the next subsection). In the construction of the global solution in the hemisphere, the boundary conditions at the equator and at the pole should remove the divergent solutions.

To provide regularity of the global solution at the middle singularity (x=a) the two P-solutions with the exponents 0 and ( $1-\varepsilon$) must be regular. This is because every retained regular solution at the equator or at the pole will, after its analytical continuation, be expressed as the sum of these two independent solutions at x=a. The second exponent (at x=a) $(1-\varepsilon)\to 2$ if $\varepsilon_{{\mbox{\tiny R}}}\to 0$ or if $\beta_{\mu}\to 0$. If $\beta_{\mu}\to\infty$ then $(1-\varepsilon)\to 7/2$. This implies, that the second independent solution of Eq. (47) with the exponent $(1-\varepsilon)$ is regular at the singular point x=a for all variables, as follows from Eqs. (39)-(41). For the first solution, however, with the exponent of "0" at x=a the physical variables $V_{\theta}$ and $V_\phi$ are limited only if a is complex, i.e. $\omega$ is complex. In this way we conclude that the solutions of the singular boundary value problem will be complex and that the process of wave-rotation interaction produces a mechanical instability.

In our initial adiabatic motion approximation for which the final Eq. (47) is derived, the complex $\omega$ should be related to one of the following physical mechanisms:
1) Shear instability:
The latitudinal differential rotation with a gradient  $\beta _{\mu }$ can favour the development of mechanical shear instability. This instability is akin to the classical Rayleigh instability mechanism or to the Kelvin-Helmholtz instability with gravity. The latitudinal instability cannot be prevented by gravity in the horizontal direction. The development of this instability will strongly depend on the values and the sign of  $\beta _{\mu }$ as well as on the wave frequency.
2) Resonant damping:
The existence of a critical latitude $x=\sin^2\theta=a$ may be the reason for the damping of the eigenmodes. For the solar rotation the location of this latitude in the hemisphere depends weakly on  $\beta _{\mu }$, see Eq. (37). The critical latitudes exist even if the rotation is rigid at the inertial frequency $\omega=2\Omega\cos\theta=2\Omega_r$. The appearance of the critical latitudes on the rotating spherical surface and their role in the formation of a wave guide around the equator is a well known phenomenon in geophysics (e.g., Longuet-Higgins 1965). A concentration of the wave amplitude at the singular latitude is very similar to the damping of adiabatic waves in MHD resonant layers (recall that the nature of the Coriolis force is very similar to that of the ponderomotive force). Physically the damping of MHD waves, for example, in the resonant Alfvén layer where $\omega=k_x V_{\rm A}(r)$, means the transformation of waves into Alfvén waves which are concentrated along the magnetic field in the resonant layer; here  $V_{\rm A}(r)$ is the radial dependence of the Alfvén velocity, k_x is the wave number along the magnetic field. This mechanism is an important damping mechanism of MHD waves; it has been used tentatively as a corona heating mechanism (e.g., Ionson 1978, 1984; Hollweg 1987). In our case the Alfvén waves are changed to inertial oscillations.

The transformation of rotation-gravity waves to inertial waves in the narrow range of critical latitudes acts as a resonant damping mechanism. The resonance of the local inertial oscillations (frequency $2\Omega\cos\theta$) with the eigenmodes (frequency $\omega$) occurs in those places where their frequencies are very close to each other. For low frequency modes, $\varepsilon_{{\mbox{\tiny R}}}\ll 1$, the resonant layer is located close to the equator, and for higher frequency modes, $\varepsilon_{{\mbox{\tiny R}}}\sim 1$, near the pole. It can be shown that generally wave motion with a frequency of $2\Omega_r$ must be produced to ensure conservation of potential vorticity in the rotating system. The width of the critical layer will depend on the values of the imaginary parts of the eigenfrequencies. Formally, if there exists a continuous spectrum, then the critical layers may be very wide. In reality, however, we may expect the formation of some "active'' latitude belts if only some selected discrete modes are excited. In these belts the greatest part of the mode energy will be concentrated in the inertial modes. The further interaction of these waves with rotation can change the rotation velocity. The evolution of the produced resonant inertial waves in their initial state should be investigated separately on the base of our general PD Eq. (27). For the further devolopment another dissipation effect as well as nonlinearity should be taken into consideration. We are not aware of investigations of the resonant damping of rotation waves through the critical latitudes in the physics of stellar pulsations.
3) Tunneling of waves:
It is well known that the tunnel leakage of wave energy from a cavity leads to a decrease of wave amplitudes in time, i.e. the tunnel effect can work if the wave frequency is complex. This is very effective for the p-modes in the solar atmosphere (Dzhalilov et al. 2000). Equation (35) or Eq. (47) can easily be transformed into a two-terms form: Z'' + U(x) Z'=0. Due to the singular points the behavior of the wave potential U(x) is complicated, especially if $\beta_{\mu}\ne 0$ and $\varepsilon_{{\mbox{\tiny R}}}\approx 1$ (we will omit here a detailed analysis of U(x)). However, a series of changes of the sign of U(x) shows the existence of cavity (oscillation) and tunnel (non-oscillation) zones in the area 0<x<1. Close to the equator at $x\to 1$ for all cases U>0, and for most cases near the pole ($x\to 0$) U<0. This means that the eigenmodes are formed in the cavity which is located around the equator. Near the pole we must choose only a decreasing solution. For example, in the rigid rotation case we have $\Theta\sim x^{\vert m\vert}\to 0$. However, for any $\beta_{\mu}\ne 0$ (differential rotation) and complex $\omega$ we will have a complex wave vector in the $\vec\theta$-direction. This means that we will have runnig waves toward the pole. As the amplitude of these waves is zero at the pole we have no reflected waves from the pole. Thus we may conclude that the tunnel effect can work throughout each hemisphere.

These questions are important, but a detailed discussion may be given only after Heun's equation has been solved.

4.3 Wave instability condition

Here we discuss two conditions: the regularity of the solutions at the pole and the approximate instability condition for the modes. If we put the P-solutions with the exponents 0 and  $(1-\gamma)$ into Eq. (43) we get $p'\sim\Theta\sim x^{\sigma_{1,2}}$. Then $V_{\theta,\phi}\sim x^{(2\sigma-1)/2}$ means that for the regularity of the solutions at the pole the condition ${\rm Re}(2\sigma)\ge 1$ must be obeyed.

On the other hand, an instability is possible when the eigenfrequencies are complex, implying a complex $\sigma$, as the wavenumber is complex (for tunneling) if $\sigma$ is complex. Otherwise, the solution of the eigenvalue problem gives a dispersion relation for $\omega$ which becomes complex if the parameter $\sigma$ is complex (of course, this is one possibility). This is seen, for instance, directly from the dispersion relation (68, see below) for the lower frequency limit. Here |m| should be changed to $2\sigma$. For the instability it follows from Eq. (44), that the necessary condition is S₂>0. It is clear that the axially-symmetric mode with m=0 is excluded in this case.

For lower values of the rotation gradient $\vert\beta_{\mu}\vert<1$ the necessary condition S₂>0 demands for the prograde waves ( $m\varepsilon_{{\mbox{\tiny R}}}>0$) the condition $\beta_{\mu}>0$, which is more realistic for stellar situations (equatorward spinning up at the surface with radius r). Rayleigh's necessary condition for instability (Rayleigh 1880; Watson 1981) says that the function $R_l=\partial^2(\Omega\sin^2\theta)/\partial\cos^2\theta$ (gradient of vorticity) must change its sign in the flow. Rewriting this function in our definitions we get that

$$% R_l=\frac{2\Omega}{\mu^2}\left[3\beta_{\mu}-\mu^2\left(1+4\beta_{\mu}\right)\right]$$ (51)

may change its sign if $\beta_{\mu}>0$. There are instability possibilities for large negative  $\beta _{\mu }$ which are not considered in this work.

The sufficient condition for instability is obtained from Eq. (44) and reads $\beta_{\mu}^2S_1^2 < S_2$.

The regularity condition at the pole $\beta_{\mu}S_1 \ge 1$ can be rewritten as

$$% \frac{\varepsilon_{{\mbox{\tiny R}}}}{m} \le \frac{\beta_{\... ...tiny R}}}^2-\beta_{\mu}\left(\beta_{\mu}+3/2\right)} = \chi_3.$$ (52)

This condition divides the phase space $\{\varepsilon_{{\mbox{\tiny R}}}/m, \varepsilon^2_{{\mbox{\tiny R}}}\}$into three parts, depending on the values of $\beta _{\mu }$. For $0\le\beta_{\mu}<1/2$ we have the following situation: if $\varepsilon_{{\mbox{\tiny R}}}^2 < 1-\beta_{\mu}(\beta_{\mu}+3/2)$ the condition Eq. (52) is fulfilled for prograde waves $\varepsilon_{{\mbox{\tiny R}}}/m\ge 0$ (region I); in the opposite case when $\varepsilon_{{\mbox{\tiny R}}}^2 > 1-\beta_{\mu}(\beta_{\mu}+3/2)$ the condition (52) is fulfilled for retrograde waves with $\varepsilon_{{\mbox{\tiny R}}}/m < -1/2$ (region III); for $\varepsilon_{{\mbox{\tiny R}}}^2 = 1-\beta_{\mu}(\beta_{\mu}+3/2)$ these regions are separated by the asymptote $\chi_3=\infty$.

Figure 1: The domains of validity of the solution regularity and of the instability conditions in the phase space $\{\varepsilon_{{\mbox{\tiny R}}}/m, \varepsilon^2_{{\mbox{\tiny R}}}\}$ for given values of the rotation gradient $\beta _{\mu }$. a) Shows the values $\varepsilon_{{\mbox{\tiny R}}}/m = \chi_3$, Eq. (52). In the area $\varepsilon_{{\mbox{\tiny R}}}/m \le \chi_3$ the solutions are limited at the pole. b) Shows the behavior of $\chi _1$ (dashed) and $\chi _2$ (solid). Between the solid and dashed lines with the same labels (values of  $\beta _{\mu }$) mode instability is possible.

Open with DEXTER

For $\beta_{\mu} > 1/2$ (strong gradients) the condition Eq. (52) is met only for retrograde waves in the range $-1/2 < \varepsilon_{{\mbox{\tiny R}}}/m \le 0$(region II). For $\beta_{\mu} = 1/2$ we get for any $\varepsilon_{{\mbox{\tiny R}}}$ that $\chi_3=-1/2$. This is the line between regions II and III. All three regions are shown in Fig. 1a. It is seen that the regularity condition is working for $\vert m\vert\ge 1$ and $\vert\varepsilon_{{\mbox{\tiny R}}}\vert\le 1$ if $\beta_{\mu}\ne 0$. For very small  $\beta _{\mu }$ only modes with large mare possible. The smallest $m \ga 1$ modes may appear in the limit $\varepsilon_{{\mbox{\tiny R}}}\approx 1$. These conclusions are correct only if the instability occurs.

Now let us consider the second condition, the complex frequency condition $\beta^2_{\mu}S^2_1 < S_2$. This inequality may be rewritten as

 

$$\chi_2 < \frac{\varepsilon_{{\mbox{\tiny R}}}}{m} < \chi_1,$$ (53)

$$\chi_{1, 2}= \frac{2}{\beta_{\mu}\left(5+2\beta_{\mu}\right)^2} \... ...-a_{\ast}\varepsilon_{{\mbox{\tiny R}}}^2\left(5+2\beta_{\mu}\right)^2}\right],$$

$$a_{\ast}= \left(1-\varepsilon_{{\mbox{\tiny R}}}^2\right)\left(1-\varepsilon_{{\mbox{\tiny R}}}^2-\beta_{\mu}^2\right),$$

$$b_{\ast}= \left(3+2\varepsilon_{{\mbox{\tiny R}}}^2\right)\left(1... ...right)- \varepsilon_{{\mbox{\tiny R}}}^2\beta_{\mu}\left(5+2\beta_{\mu}\right).$$

In the limiting cases we have

 

$$\chi_1 \approx \frac{12\left(1+\beta_{\mu}\right)}{\beta_{\mu}\le... ...chi_2 \approx 0 \ \ \mbox{for } \ \varepsilon_{{\mbox{\tiny R}}} \rightarrow 0,$$ (54)

$$\chi_1 \approx 0 , \ \ \chi_2 \approx -\frac{8\beta_{\mu}}{\left(... ...mu}\right)^2} \ \ \ \mbox{for } \ \varepsilon_{{\mbox{\tiny R}}} \rightarrow 1.$$ (55)

Figure 1b shows for $\vert m\vert\ge+ 1$ and $\varepsilon_{{\mbox{\tiny R}}}\le 1$ the curves $\chi _1$and $\chi _2$ versus $\varepsilon^2_{{\mbox{\tiny R}}}$ for a wide range of  $\beta _{\mu }$. A comparison of Figs. 1a and b shows that in region III with $\varepsilon_{{\mbox{\tiny R}}}/m < -1/2$ the mode instability will never appear. In region I instability is possible for prograde waves, if $\beta_{\mu}<1/2$ and $\varepsilon_{{\mbox{\tiny R}}}^2 < 1-\beta_{\mu}(\beta_{\mu}+3/2)$. With decreasing  $\beta _{\mu }$ the solid and dashed curves for the same $\beta _{\mu }$ are close to $\varepsilon_{{\mbox{\tiny R}}}\approx 1$. Instability of retrograde waves is possible for strong gradients of $\beta_{\mu} > 1/2$ only in the range $-0.2 \le \varepsilon_{{\mbox{\tiny R}}}/m<0$. Here $(\chi_2)_{\rm min}=-0.2$ is valid for $\varepsilon^2_{{\mbox{\tiny R}}}=1$ and $\beta_{\mu}=5/2$, which follows from Eq. (55).

The total condition for the existence of spatially stable but temporally unstable waves reads as follows:

$$% \chi_2 < \frac{\varepsilon_{{\mbox{\tiny R}}}}{m} < \min\left(\chi_1, \chi_3\right).$$ (56)

Figure 2: The wave instability areas (hatched) in the phase space from an overlap of Figs. 1a and b and from the condition $\vert\epsilon _R/m\vert < \sqrt {\epsilon _R^2}$ for the integer m (solid line in  a) and dashed in  b)). The labels in the areas are the $\beta _{\mu }$ values. Both for prograde and retrograde waves  $\beta _{\mu }$ is limited: for prograde waves the instability is possible if $\beta _{\mu }<0.5$, for retrograde waves the instability is possible if $\beta _{\mu }>3.5$.

Open with DEXTER

Figures 2a,b show the validity ranges of this condition for some typical values of  $\beta _{\mu }$. The hatched areas are places where mode instability is possible. These figures are obtained by overlapping Figs. 1a and b. The hatched areas of prograde waves are strongly restricted also by the condition $\epsilon_R/m < \sqrt{\epsilon_R^2}$ because the integer $m\ge 1$. For prograde waves on both sides these hatched areas become very narrow: with decreasing  $\beta _{\mu }$ the extent of the hatched area decreases and tends to the point $\{\varepsilon_{{\mbox{\tiny R}}},m\}= \{1,1\}$. We will see that this is the solar case. The hatched instability areas disappear with decreasing $\beta _{\mu }$. This means we have a lower limit $(\beta_{\mu})_{\rm min}$.

For retrograde waves it is sufficient to write the condition as $\chi_2 < \varepsilon_{{\mbox{\tiny R}}}/m < \chi_3$. In Fig. 2b $\chi _2$ is the dash-dot curve and $\chi_3$ is the solid curve. Instability is possible if $\beta _{\mu }>3.5$ for $0\ge\varepsilon_{{\mbox{\tiny R}}}/m \ge -0.2$.

In Fig. 2 only the hatched areas for given $\beta _{\mu }$ are ranges of possible solutions, if wave instability occurs. Outside these hatched areas regular solutions are impossible. The case without instability (neutral oscillations) must be investigated separately.

4.4 Solar rotation profile

Let us consider at which places we might expect mode instability in the Sun. Unfortunately, it is not clear how the core rotates. Nevertheless some rotation gradients close to the Sun's centre might exist, and we could expect mode instability there. It is known from helioseismology that the radiative interior has a very small  $\beta _{\mu }$, but the exact value is unknown. We have better information on the rotation profile of the solar envelope, including the tachocline. Helioseismology data may be described by different approximate formulae. One of these is (Charbonneau et al. 1998)

 

$$\Omega (r, \theta) = \Omega_{\rm c}+\frac{1}{2}\left[1+{\rm erf}\... ...\bar{\Delta}\right)\right] \left(\Omega_{\rm s}(\theta)- \Omega_{\rm c}\right),$$ (57)

$$\Omega_{\rm s} = \Omega_{\rm eq}+c_1\cos^2\theta+c_2\cos^4\theta , \ \ \bar{\Delta}=\frac{r-r_{\rm c}}{w},$$

$$\Omega_{\rm c}/2\pi = 432.8~{\rm nHz} , \ \ \Omega_{\rm eq}/2\pi = 460.7~{\rm nHz},$$

$$c_1 = -62.69~{\rm nHz} , \ \ c_2 = -67.13~{\rm nHz},$$

where $r_{\rm c}=0.713~R_{\odot}$ is the radius at the bottom of the convective zone, and $w=0.025~R_{\odot}$ is the tachocline thickness. We can easily check that our approximation $\tilde\Omega/\Omega_{\rm c} \ll 1$ is applicable in most cases, if Eq. (57) is represented by $\Omega =\Omega_{\rm c} + \tilde\Omega(r,\theta)$. The maximum value of $(\tilde\Omega/\Omega_{\rm c})$ is in the convection zone: $\approx$0.06 for $\theta =\pi/2$ (equator) and $\approx$-0.23 for $\theta=0$ (pole). From Eq. (57) we find the latitudinal gradient of rotation

$$% \beta_{\mu} = -\frac{1+{\rm erf}(\Delta)}{2\Omega}\sin^2\theta \left(c_1+2c_2\cos^2\theta\right),$$ (58)

which should be acceptable for estimating local values of the gradients.

Using this formula we show in Fig. 3 the $\beta_{\mu}(\theta)$dependence for different $r/R_{\odot }$. We see that in the Sun $\beta_{\mu}\le 0.15$, and the maximum is in the photosphere. Figure 2 implies that unstable retrograde waves are not present in the solar case. Instability of prograde waves in the Sun occurs in the upper right corner of Fig. 2a which is enlarged in Fig. 4. Here we see that the instability area disappears when $\beta_{\mu}< 3$ $\times$ 10^-4. This boundary is located at the bold horizontal line in Fig. 3. Thus the prograde waves become unstable in the Sun in those places where 3 $\times$ 10^-4 $\le$ $\beta _{\mu }$ $\le$ 0.15. This means that instability is possible in the area $r/R_{\odot}>0.66$ which includes the greater part of the tachocline, the convective zone, and the photosphere. With increasing r the instability zone expands from middle to high latitudes. Figures 2a and 4 show that the instability occurs at high frequencies ( $\varepsilon_{{\mbox{\tiny R}}} \la 1$) and on global scales ( $\varepsilon_{{\mbox{\tiny R}}}/m > 0.5$). Considering that $m\ne 0$ is an integer we get m=1.

However, our instability analysis is based on the general Riemann scheme of Heun's equation, which is valid only if the middle singular point x=a is far from the other edges at x=0 and x=1. Thus, the limiting cases $\varepsilon_{{\mbox{\tiny R}}} \rightarrow 0$ and $\varepsilon_{{\mbox{\tiny R}}} \rightarrow 1$ (the latter is more important for the solar case) should be considered separately. In these limiting cases the regularity condition Eq. (48) may be changed, and the curve in Fig. 4 limiting the instability areas from below may be shifted. In this case instability with higher m-modes should be expected.

	Figure 3: The local estimate of the logarithmic gradient of the solar rotation frequency  $\beta _{\mu }$ for real solar data from helioseismology depending on the co-latitude and the radial distance (the labels are values of $r/R_{\odot }$). The bold horizontal line is $\beta _{\mu }= 3$ $\times$ 10-4, above which the prograde waves become unstable (see next picture).
Open with DEXTER

	Figure 4: Enlarged part of Fig. 2a for the smallest gradients of rotation. The labels 1, 2,..., 6 correspond to $\beta _{\mu }= 3$ $\times$ 10-4, 5 $\times$ 10-4, 1 $\times$ 10-3, 2 $\times$ 10-3, 3 $\times$ 10-3, 4 $\times$ 10-3. Areas of instability exist only if $\beta _{\mu } \ge 3$ $\times$ 10-4.
Open with DEXTER

5 The low-frequency waves

5.1 General remarks

After the qualitative analysis of Heun's Eq. (47) we can start a quantitive analysis. Note that the qualitative conclusions drawn above are valid for the more general Eq. (35) with the $\Lambda$ term. Heun's equation with four singularities in the general case is solved by a series of hypergeometric Gauss functions. A similar task has been considered for the damping of MHD waves at resonance levels by Dzhalilov & Zhugzhda (1990). We will start to study Eq. (47) for some simple limiting cases. At high frequencies ( $\varepsilon^2_{{\mbox{\tiny R}}}\approx 1$, where shear instability is acting in the Sun) and at low frequencies ( $\varepsilon^2_{{\mbox{\tiny R}}} \ll 1$, where the waves are stable against shear instability in the Sun) Heun's equation is much simpler. In these cases the singular level x=a is shifted either to the pole or to the equator. For both cases the solutions are expressed by one hypergeometric function.

In the present work we consider particularly the second case. Let $\varepsilon^2_{{\mbox{\tiny R}}} \ll 1+\beta_{\mu}$. Then we have $a\approx 1$ and $h\approx 1$. Equation (47) is now the hypergeometric equation:

$$% x(1-x)Y^{\prime\prime}+\left[\gamma-\left(\tilde{\alpha}+\t... ...a}+1\right)x\right] Y^{\prime}-\tilde{\alpha}\tilde{\beta}Y=0,$$ (59)

where all parameters are defined by Eqs. (44), (46) and (50). In these definitions  $\varepsilon^2_{{\mbox{\tiny R}}}$ should be set to zero. Strictly speaking, this equation is valid for $x\le 1+O(\varepsilon^2_{\mbox{\tiny R}})$. Then the $\Theta$-part of the pressure perturbations, Eq. (43), is expressed by two Gaussian hypergeometric functions:

$$% \Theta = C_1x^{\sigma_1}Y_1(x) + C_2x^{\sigma_2}Y_2(x) .$$ (60)

$$% {\rm Here \ }Y_{1, 2}(x)=F\left(\tilde{\alpha}, \tilde{\beta}; \gamma; x\right)\vert _{\sigma=\sigma_{1, 2}},$$ (61)

and C_{1, 2} are arbitrary constants. This general solution includes wave instability for larger  $\beta _{\mu }$ too. This could be realized perhaps in other, younger stars. For the Sun we have $\beta_\mu \ll 1$. We will finish this paper by considering in detail the more popular case when rotation is uniform, $\beta_{\mu} = 0$.

5.2 Rigid rotation case

Using the conditions $\beta_{\mu} = 0$ and $\varepsilon^2_{{\mbox{\tiny R}}} \ll 1$ the parameters in the solution Eq. (60) are greatly simplified. Because $2\sigma=\pm \vert m\vert$ only a regular solution at the pole (x=0) will be left. In the standard definitions of hypergeometric functions (Abramowitz & Stegun 1984) we have

 

$$\Theta = Cx^{\vert m\vert/2}F(a, b; c; x),$$ (62)

$$a = \frac{1}{2}\left(-\frac{1}{2}+\vert m\vert+q\right), \ b=\frac{1}{2}\left(-\frac{1}{2}+\vert m\vert-q\right),$$

$$c = 1+\vert m\vert , ~~~ \ q=\sqrt{\frac{1}{4}-\frac{m}{\varepsilon_{{\mbox{\tiny R}}}}} , \ \ {\rm and} \ \ C ={\rm const}.$$

Note that c-a-b=3/2. The analytical continuation of the solution Eq. (62) to the equator, $x\to 1$, gives

 

$$% \Theta Cx^{\vert m\vert/2}\left[AF\left(a, b; -\frac{1}{2}; 1-x\right)\right. +\left.B(1-x)^{3/2} F\left(c-a, c-b; \frac{5}{2}; 1-x\right)\right].$$ (63)

Here the continuation coefficients are (Abramowitz & Stegun 1984)

  

$$A = \frac{\Gamma(c)\Gamma(3/2)}{\Gamma(c-a)\Gamma(c-b)}=F(a, b; c; 1),$$ (64)

$$B=\frac{\Gamma(c)\Gamma(-3/2)}{\Gamma(a)\Gamma(b)}\cdot$$ (65)

Equations (62) and (63) imply that pressure is bounded over the whole hemisphere. Now let us consider the velocity components. Putting Eq. (62) into Eq. (40) we get

 

$$% V_{\theta} \pm\frac{CQ\varepsilon_{{\mbox{\tiny R}}}^2}{r\rho_0 i\omega} \fr... ...frac{m}{\varepsilon_{{\mbox{\tiny R}}}}-\vert m\vert\right)F(a, b; c; x)\right. -\left. x\frac{2ab}{c} F(a+1, b+1; c+1; x)\right].$$ (66)

Taking into account that F(a,b;c;0)=1, we obtain from the regularity of  $V_{\theta}$ at x=0 that $\vert m\vert\ge 1$. Axially-symmetric modes m=0 cannot be formed. The continuation of the solution Eq. (66) to the equator, $x\to 1$, gives

 

$$V_{\theta} = \pm\frac{CQ\varepsilon_{{\mbox{\tiny R}}}^2}{r\rho_0... ...^{(\vert m\vert-1)/2} \left(\frac{A}{\sqrt{1-x}}{\cal L}_1+B{\cal L}_2\right) ,$$ (67)

$${\cal L}_1 = \left(\frac{m}{\varepsilon_{{\mbox{\tiny R}}}}-\vert m\vert\right)F\left(a, b; -\frac{1}{2}; 1-x\right)$$

$$\hspace*{7mm} -x4abF\left(a+1, b+1; \frac{1}{2}; 1-x\right) ,$$

$${\cal L}_2 = \left[\left(\frac{m}{\varepsilon_{{\mbox{\tiny R}}}}... ...m\vert\right)(1 - x) + 3x\right] F\left(c - a, c - b; \frac{5}{2}; 1 - x\right)$$

$$\hspace*{7mm} +\frac{4}{5}x(1 - x)(c - a)(c - b)F\left(c\! -\! a\!+\!1, c \!-\! b \!+\! 1; \frac{7}{2}; 1 - x\right).$$

As for x=1 the functions ${\cal L}_1$ and ${\cal L}_2$ are limited, we have the relation $A\equiv 0$ from the regularity requirement. Using a property of the gamma functions (the fact that they can be presented as an infinite product) in Eq. (64) we get the condition of quantization

$$% q^2 = \left(2l+\vert m\vert+\frac{5}{2}\right)^2, ~ ~ ~ l=0, 1, 2, ...$$ (68)

From this we get the simple dispersion relation

$$% \frac{\varepsilon_{{\mbox{\tiny R}}}}{m}=\frac{1}{m}\frac{\... ...Omega}=-\frac{1} {(2+2l+\vert m\vert)(3+2l+\vert m\vert)}\cdot$$ (69)

Since $\varepsilon_{{\mbox{\tiny R}}}/m < 0$ for any $l\ge 0$, only retrograde modes (as seen in the rotating frame) are possible.

5.3 Spectrum of retrograde modes

The new dispersion relation Eq. (69) completely differs from the dispersion relation of the almost toroidal r-modes. Their dispersion relation can be derived from Eq. (69) if we formally set $l_{\rm s} = 2(1+l)+\vert m\vert$. Then

$$% \omega=-2\Omega\frac{m}{l_{\rm s}(1+l_{\rm s})}\cdot$$ (70)

However, here the degrees $l_{\rm s}$ are functions of m. Due to the coupling of the modes in our case, the eigenfunctions can never be expressed by the associated Legendre functions  $P^m_{l_{\rm s}}$. For the r-modes the axi-symmetric modes (m=0) with $\omega=0$ are possible, while our case rejects this. Unlike the r-modes the spectrum Eq. (69) has a maximum at m² = m²₀ =2(1+l)(3+2l). If $l \gg 1$ we have $\vert\varepsilon_{{\mbox{\tiny R}}}\vert _{\rm max} \approx 1/8l$. The spectrum of retrograde waves is shown in Fig. 5. Rossby waves in geophysics (Pedlosky 1982) have a similar spectrum. Using Eq. (69) we may define the local phase and group velocities (at fixed latitude $\theta_0$ and radial distance r₀) in the azimuthal plane

  

$$v_{\rm ph} = -\frac{2\Omega r_0\sin\theta_0}{(2+2l+\vert m\vert)(3+2l+\vert m\vert)},$$ (71)

$$v_{\rm gr} = 2\Omega r_0\sin\theta_0\frac{m^2-2(1+l)(3+2l)}{(2+2l+\vert m\vert)^2(3+2l+\vert m\vert)^2}\cdot$$ (72)

We see that the group velocity changes its sign at the maximum of the spectrum, where m²=m²₀. Figure 6 shows that the l=0 mode has a maximum group velocity. Long waves carry energy opposite to the rotation direction, while a packet of short waves is carried in the same direction. The facts that the frequencies $\omega(m)$ have a maximum (two different |m| correspond to the same $\omega$) and that at the maximum of $\omega=\omega_{\rm max}$ the group velocity changes direction hint at the existence of two branches of oscillations. Solving Eq. (69) for |m| we get these branches. Let $\omega>0$ and m=-|m|. Then

 

$$\vert m\vert=\frac{1}{2\varepsilon_{{\mbox{\tiny R}}}}\left(w_1\mp\sqrt{w_2}\right)=m_{1, 2},$$ (73)

$$w_1 = 1-\varepsilon_{{\mbox{\tiny R}}}-4\varepsilon_{{\mbox{\tiny R}}}(1+l),$$

$$w_2 = (1-\varepsilon_{{\mbox{\tiny R}}})^2-8\varepsilon_{{\mbox{\tiny R}}}(1+l).$$

	Figure 5: The spectrum $\omega (l,m)$ of low-frequency retrograde modes. The numbers on the curves are the l values. For the calculation of the spectrum, Eq. (69), $\Omega =\Omega _{\rm eq}=460.7$ nHz is used.
Open with DEXTER

	Figure 6: The normalized group velocity $v_{\rm g}/2\Omega r_0\sin\theta_0$ as a function of the azimuthal numbers for given degrees l.
Open with DEXTER

Figure 7: The possible domain of the existence of eigenmodes for given Rossby numbers $\varepsilon_{{\mbox{\tiny R}}}$ (the numbers on the curves). For example, the case $\varepsilon_{{\mbox{\tiny R}}}=0.002$ (close to the 22-year modes) is emphasized. All possible values of the azimuthal numbers (m=m1 and m=m2) are on this curve for given discrete l in the range $l \le l_{\rm max}$ (see text). If $\varepsilon_{{\mbox{\tiny R}}}>0.1$, the eigenmodes disappear.

Open with DEXTER

From this we get an upper limit for l if $\varepsilon_{{\mbox{\tiny R}}}$ is given: $l\le l_{\rm max}= (1-\varepsilon_{{\mbox{\tiny R}}})^2/8\varepsilon_{{\mbox{\tiny R}}} -1$. As $l_{{\rm max}}\ge 0$ we get an approximate condition for the existence of the modes: $\varepsilon_{\mbox{\tiny R}}\la 0.1$. For such degrees of l the azimuthal numbers are also limited: $m_1\le\vert m\vert\le m_2$. For $\varepsilon_{{\mbox{\tiny R}}}\to 0$ we have $l_{\rm max}\to 1/8\varepsilon_{{\mbox{\tiny R}}}$, $m_2\to 1/\varepsilon_{{\mbox{\tiny R}}}$, and $m_1\to 0$. However, considering the regularity of the solutions Eq. (66), we must take $m_1\ge 1$. m₁=1 if $l=(\sqrt{1+4/\epsilon_R}-7)/4$. For $l=l_{\rm max}$ we have m₁=m₂. This situation is shown in Fig. 7 for different values of  $\varepsilon_{{\mbox{\tiny R}}}$. A decrease of the frequency decreases the domain of existence of the modes.

Putting Eq. (73) into Eqs. (71), (72) gives

  

$$\frac{v_{\rm ph}^{\pm}}{2\Omega r_0\sin\theta_0} = -\frac{w_1\pm \sqrt{w_2}}{4(1+l)(3+2 l)},$$ (74)

$$\frac{v_{\rm gr}^{\pm}}{2\Omega r_0\sin\theta_0} = -\frac{w_2 \pm w_1 \sqrt{w_2}}{4(1+l)(3+2 l)}\cdot$$ (75)

	Figure 8: The phase velocities of fast (solid) and slow (dashed curves) modes versus l. The numbers at the curves are $\varepsilon_{{\mbox{\tiny R}}}$ values, similar to Fig. 7. Here the velocities are normalized to ( $-2\Omega~ r_0\sin \theta_0$). Every curve is restricted to the range $l \le l_{\rm max}$.
Open with DEXTER

Here $v_{\rm ph}^+$ is the phase velocity of the fast modes and  $v_{\rm ph}^-$ that of the slow modes. For $l=l_{\rm max}$ we have $v_{\rm gr}^{\pm}=0$ and $v_{\rm ph}^+=v_{\rm ph}^-$. In Eq. (73) and in Fig. 7 the m₁ branch corresponds to the fast mode, but m₂ to the slow modes, since $m_2\ge m_1$. In Fig. 8 the normalized phase velocities (with inverse sign) for the selected values of  $\varepsilon_{{\mbox{\tiny R}}}$ in Fig. 7 are shown versus l. Both branches are retrograde modes ( $v_{\rm ph}^{\pm}<0$). Using $\Omega R_{\odot} \sim 2$ km s^-1 for the Sun, we get from Fig. 8 very slow phase velocities. The fast wave velocity (solid lines) depends more strongly on l. With increasing $\varepsilon_{{\mbox{\tiny R}}}$ both branches are accelerated.

In Fig. 9 the group velocities are presented in the same way. For fast waves the group velocity is always parallel to the phase velocity ( $v_{\rm gr}^+<0$), while for the slow waves we have the opposite behavior  $v_{\rm gr}^->0$. Slow mode packets carry off energy in the rotation direction. $\vert v_{\rm gr}^+\vert>\vert v_{\rm gr}^-\vert$ is valid everywhere. With deceasing $\varepsilon_{{\mbox{\tiny R}}}$ the range $l \le l_{\rm max}$ shifts to the right, and it is seen in Fig. 9 that  $v_{\rm gr}^-$ for such low $\varepsilon_{{\mbox{\tiny R}}}$ is almost zero.

Note that m=l modes are always fast modes.

	Figure 9: Absolute values of the group velocities of fast (solid) and slow (dashed) modes normalized to $2\Omega~ r_0\sin\theta_0$ versus l for selected  $\varepsilon_{{\mbox{\tiny R}}}$. For fast waves  $v_{\rm gr}^+<0$, for slow modes  $v_{\rm gr}^->0$. For $l=l_{\rm max}$ we have $v_{\rm gr}^+=v_{\rm gr}^-=0$.
Open with DEXTER

5.4 The eigenfunctions

Taking into account the quantization condition Eq. (68) in the solutions Eqs. (63), (67), (39), and (41) we obtain the eigenfunctions. Turning from complex velocities into the real displacements, $\vec{V}=-i\omega\vec\xi$ (remember that $\vec V$ is the velocity seen in the rotating frame), we get

    

$$p' = Q(r)\Theta(\theta)\cos\left[m\left(\phi-2\frac{\varepsilon_{{\mbox{\tiny R}}}}{m}\Omega t\right)\right],$$ (76)

$$\xi_{\theta} = \xi_{\theta}^{\star}\cos\left[m\left(\phi -2\frac{\varepsilon_{{\mbox{\tiny R}}}}{m}\Omega t\right)\right],$$ (77)

$$\xi_{\phi} = \xi_{\phi}^{\star}\sin\left[m\left(\phi -2\frac{\varepsilon_{{\mbox{\tiny R}}}}{m}\Omega t\right)\right],$$ (78)

$$\xi_{r} = \xi_{r}^{\star}\cos\left[m\left(\phi -2\frac{\varepsilon_{{\mbox{\tiny R}}}}{m}\Omega t\right)\right].$$ (79)

Here the amplitude functions are

   

$$\xi_{r}^{\star} = \frac{1}{\omega^2-N^2}\frac{\Theta}{\rho_0}\frac{{\rm d}Q(r)}{{\rm d}r},$$ (80)

$$\xi_{\theta}^{\star} = \frac{Q(r)}{4\Omega^2r\rho_0}\xi_{\theta A}(\theta),$$ (81)

$$\xi_{\phi}^{\star} = \frac{Q(r)}{4\Omega^2r\rho_0}\xi_{\phi A}(\theta) .$$ (82)

The amplitudes are expressed by Jacobi's polynomials:

   

$$\Theta = B(\sin\theta)^{\vert m\vert}\cos^3\theta ~F_1(\theta),$$ (83)

$$\xi_{\theta A} = B(\sin\theta)^{\vert m\vert-1}\cos^2\theta \left(\frac{m}{\varepsilon_{{\mbox{\tiny R}}}}F_1-\frac{4lk}{5}\sin^2\theta F_2\right),$$ (84)

$$\xi_{\phi A} = -\frac{B}{\varepsilon_{{\mbox{\tiny R}}}}\frac{2lk}{5} \sin(2\theta)\cos^2\theta(\sin\theta)^{\vert m\vert}F_2(\theta),$$ (85)

 

$$F\left(-l, k; \frac{5}{2}; \cos^2\theta\right)$$                                 F₁ = (86)

$$(-1)^l \frac{l!\Gamma(5/2)}{\Gamma(l+5/2)} P_l^{(\vert m\vert, 3/2)}(\cos 2\theta)$$ =

$$\sum_{j=0}^{l} \frac{(-l)_j(k)_j}{(5/2)_j~j!}(\cos\theta)^{2j},$$ =

 

$$F\left(1-l, 1+k; \frac{7}{2}; \cos^2\theta\right)$$                                 F₂ = (87)

$$(-1)^{l-1} \frac{(l-1)!\Gamma(7/2)}{\Gamma(l+5/2)}P_{l-1}^{(1+\vert m\vert, 5/2)}(\cos 2\theta)$$ =

$$\sum_{j=0}^{l-1} \frac{(1-l)_j(1+k)_j}{(7/2)_j~ j!}(\cos\theta)^{2j}.$$ =

$$% k=\frac{5}{2}+\vert m\vert+l , \ \ (f)_j=\frac{\Gamma(f+j)}{\Gamma(f)},$$ (88)

and $F_2\equiv 0$ if l=0. The eigenfunctions of the retrograde low-frequency modes Eqs. (76)-(79) are written for the hemisphere $0\le\theta\le\frac{\pi}{2}$. In the other hemisphere we have to change $\cos\theta$ to  $\vert\!\cos\theta\vert$ and also to change the sign of  $\xi_{\theta}$. We put C=1 because all eigenfunctions are multiplied by the solution of the r-Eq. (33). Changing the sign of m the function  $\xi_{\phi}$ keeps its sign because $\varepsilon_{{\mbox{\tiny R}}}/m < 0$. It follows from Eqs. (86) and (87) that the eigenfunctions are represented by a limited number of trigonometric functions. The number of terms is defined by l in the range $l \le l_{\rm max}$. This is qualitatively different from the representation of the rotation-mode eigenfunctions by infinite series of spherical harmonics.

Figure 10: The amplitudes of eigenfunctions versus co-latitude for given pairs (l,|m|), Eqs. (83)-(85). The first row shows the  $\Theta (\theta )$ pressure function normalized to its own maximum. From left to right all curves are related to the wave numbers first given above. The middle and last rows are similar to the first row, but now show the latitudinal ( $\xi _{\theta A}$) and azimuthal ( $\xi _{\phi A}$) amplitude functions, normalized to the maximum of ( $\xi _{\theta A}$).

Open with DEXTER

Table 1: Results for the 22-year period: frequency deviations $\delta =\omega _{22}-\omega$ for permitted quantum numbers (l, |m|).

Let us consider some particular cases.

a)

Pole ($\theta=0$): taking into account B F₁=1 we have $\Theta= p'=\xi_r=\xi_{\phi}=0$. $\xi_{\theta}=0$ if |m|>1. For |m|=1we have $\xi_{\theta A}=-2(3+2l)(2+l)$.

b)

Equator ( $\theta=\frac{\pi}{2}$): here F₁=F₂=1 and $p'=\xi_{\theta}= \xi_{\phi}=0$.

c)

l=0 case: then B=F₁=1, F₂=0,

$$\begin{eqnarray*}\Theta=(\sin\theta)^{\vert m\vert}\cos^3\theta, \end{eqnarray*}$$

$$\begin{eqnarray*}\xi_{\theta A}=-(2+\vert m\vert)(3+\vert m\vert)\cos^2\theta(\sin\theta)^{\vert m\vert-1}, \end{eqnarray*}$$

and $\xi_{\phi}=0$ (in our $\varepsilon^2_{{\mbox{\tiny R}}} \ll 1$ limit).

d)

$l \gg \vert m\vert$, |m| is not large, and $\theta\ne 0,\frac{\pi}{2}$as well.

The latter case is important because the range $l \le l_{\rm max}$ is large for small  $\varepsilon^2_{{\mbox{\tiny R}}}$. Using the asymptotic formula for Jacobi's polynomial (or the hypergeometric function) of large degree we have

$$\Theta \simeq \frac{\vert m\vert!}{\sqrt{\pi l}} \frac{\cos\lambda \cos\theta}{l^{\vert m\vert}\sqrt{\sin\theta}} ,$$ (89)

$$\xi_{\theta A} \simeq \frac{\vert m\vert!}{\sqrt{\pi l}} \frac{\c... ...}{\varepsilon_{{\mbox{\tiny R}}}} + \frac{14}{5}k\tan\theta \tan\lambda\right),$$ (90)

$$\xi_{\phi A} \simeq \frac{14}{5}\frac{\vert m\vert!k}{\varepsilon... ...ny R}}}\sqrt{\pi l}} \frac{\sin\lambda}{l^{\vert m\vert}\sqrt{\sin\theta}}\cdot$$ (91)

$$\begin{eqnarray*}\lambda=\left(2l+\frac{5}{2}+\vert m\vert\right)\theta -\frac{\pi}{2}\left(\vert m\vert+\frac{1}{2}\right)\cdot \end{eqnarray*}$$

These extremely simple asymptotic formulae for the eigenfunctions might be used in most cases. The formula for the case of large l with large |m| can also be found elsewhere, e.g. in the book of Bateman & Erdélyi (1953).

In Fig. 10 for some typical selected pairs of (l,m) the amplitude functions, Eqs. (83)-(85), are shown as functions of $\theta$. The first row is the pressure function  $\Theta (\theta )$ function normalized to its maximum. The first l=1 figure represents $\Theta$ for different |m|(|m| increases from left to right). As the eigenfunctions are multiplied by a  $\sin^{\vert m\vert}\theta$ factor, the amplitudes are strongly suppressed around the pole. Increasing |m| for a given l shifts the maxima toward the equator. l is the surface node number of the $\Theta (\theta )$ function. On the contrary, increasing l for a given |m| (the second figure of the first row with m=1 in Fig. 10) suppresses the amplitudes around the equator, and the maximum is shifted toward the pole. l=|m| is the equilibrium case. In the third figure of the first row the equilibrium latitude with maximum amplitude is defined by $\theta=20^{\circ}$ for all l=|m| modes.

The second and third rows of Fig. 10 are the latitudinal ( $\xi _{\theta A}$) and azimuthal ( $\xi _{\phi A}$) eigenfunction amplitudes, respectively, normalized to the maximum of  $\xi _{\theta A}$, see Eqs. (81) and (82). The latitudinal amplitude behaviour is similar to that of the pressure. The azimuthal amplitudes are smaller than the latitudinal amplitudes, but with a change of l a redistribution of the amplitudes will not take place. $\xi _{\phi A}$ has practically the same amplitude at all latitudes and for all l.

Figure 10 implies that we can expect an interesting behaviour of the eigenfunction amplitude, when both l and |m| are large. Suppression from two sides may give rise to a concentration of wave energy in narrow latitude bands. For example, this is the case for the 22-year solar mode.

For the 22-year modes we take $\omega_{22}=2\pi$ (1.441 nHz), for which $\varepsilon_{{\mbox{\tiny R}}}\approx 0.0016$. Then we derive from the equations after Eq. (73) the limiting values of the integer l, $11\le l\le 76$. For all lin this range we find from Eq. (73), rounding off, integer azimuthal numbers m₁ and m₂ for the fast and slow modes, respectively. Putting these integer numbers into Eq. (69) we find the deviation $\delta =\omega _{22}-\omega$ from the central frequency due to the integer azimuthal numbers. Additional slow modes are also possible in the interval l=0-10. The results are given in Table 1 in Appendix A. It is seen that the fast modes with low l have larger deviations. This table includes all possible (l,|m|) pairs which correspond to the 22-year period. For some example pairs of (l,|m|) we plot in Fig. 11 the latitude dependence of the quantity $(\xi^2_{\theta}+ \xi^2_{\phi})^{1/2}$, averaged over the wave period, which characterizes the energy density of the modes. The hemisphere is divided into two equal parts: slow modes are located around the equator (solid lines), fast modes are concentrated around the pole. Each (l,|m|) pair is located in a narrow latitude band. Note that the slow modes (the group velocity of which is in the rotation direction) with sunspot-like spatial scales are at latitudes of  $30{-}40^\circ$.

5.5 Flow patterns

The eigenfunctions Eqs. (76)-(79) allow us to discuss the flow character produced by the waves, even if the solution of the radial equation Q(r) is unknown. Excluding from these equations the time-dependent phase we can obtain the trajectory equations of the fluid elements. In the meridional  $(\vec r,\vec\theta)$ plane we have

$$% \frac{\xi_{\theta}}{\xi_{r}}=\varepsilon_{{\mbox{\tiny R}}}... ...i(r)} {r~Q'(r)} \frac{\xi_{\theta A}}{\Theta}=\tan(\alpha_r).$$ (92)

It follows from this that the motion in the meridional plane is linear. The poloidal displacement vector in the meridional plane is inclined to the radius by an angle $\alpha_r$. For $\alpha_r\approx\frac{\pi}{2}$we have meridional flows and for $\alpha_r\approx 0$ radial flows. However, since $\xi_{\theta}/\xi_r\sim 1/\sin(2\theta)$, for any r the motions around the pole and around the equator will be almost meridional ( $\vert\xi_{\theta}\vert \gg \vert\xi_r\vert$). We see from Eq. (92) also that the direction of the flow vector will be changed to $\pi/2$ if we pass along the radius from the node of the Q(r)-function to Q'(r). Note that every node of Q'(r) is located between the neighbouring two nodes of Q(r). For large radial numbers (such orders are expected) these nodes are located close to each other. Then we have a complicated motion in the meridional plane.

	Figure 11: The normalized energy density of the 22-year fast (dashed) and slow (solid) modes versus co-latitude. The numbers on the curves are (l,|m|) pairs taken from Table 1.
Open with DEXTER

At the surface of the cone over $(\vec r,\vec\phi)$ the trajectory of each fluid element is an ellipse around the equilibrium point. The cone displacement equation is

$$% \frac{\xi_{\phi}^2}{\xi_{\phi}^{\star 2}}+\frac{\xi_{r}^2}{\xi_{r}^{\star 2}}=1.$$ (93)

The ratio of the semiaxes $\xi_{\phi}^*$ and $\xi_r^*$

$$% \frac{\xi_{\phi}^{\star}}{\xi_{r}^{\star}}=-\varepsilon_{{\... ...rac{Q(r)\psi(r)}{rQ'(r)}\frac{4lk}{5}\frac{F_2}{F_1}\sin\theta$$ (94)

means that close to the axis of rotation the flow is directed along this axis. On the bigger cone through the equator, where F₁=F₂=1, the motion is mainly azimuthal since lk>1. In this case the nodes of the Q(r)-function at the radius are an exception. For large $l \gg 1$ the ratio $\xi_{\phi}^*/\xi_r^*\sim \tan(2l\theta)/\cos\theta$ implies that the character of the flows is changing at middle latitudes.

On the surface of any sphere with a radius r the motion of the fluid elements is on trajectories of the ellipse

$$% \frac{\xi_{\theta}^2}{\xi_{\theta}^{\star 2}} +\frac{\xi_{\phi}^2}{\xi_{\phi}^{\star 2}}=1.$$ (95)

Here the ratio of the semiaxes is independent of the radius:

$$% \frac{\xi_{\phi}^{\star}}{\xi_{\theta}^{\star}}=-\frac{4lk}... ...n_{{\mbox{\tiny R}}}}F_1-\frac{4lk}{5} F_2 \sin^2\theta}\cdot$$ (96)

It is seen that near the equator and near the pole the motion is mainly in the $\vec{e_{\theta}}$ direction. For large $l \gg 1$ this ratio is

$$% \frac{\xi_{\phi}^{\star}}{\xi_{\theta}^{\star}} \simeq \fra... ...n_{{\mbox{\tiny R}}}}+\frac{14}{5}k\tan\theta\tan\lambda}\cdot$$ (97)

The flow pattern considered here drifts opposite to the sense of rotation with a speed  $v_{\rm ph}$ in a frame rotating with the star.

6 Conclusions

In the present work we have derived the general PDE governing non-radial adiabatic long-period (with respect to the rotation period) linear oscillations of a slowly and differentially rotating star. This general equation includes all high-order g-modes and all possible hybrids of rotation modes as well as their mutual interaction. The main objective of the present paper was to obtain an analytic solution of this equation for the extreme low-frequency limit and to show that the solutions cannot be expressed by the associated Legendre functions. Thus, in the low-frequency limit the pulsation theory of stars meets serious difficulities. For this purpose we used some simplifying approximations. The geophysical "traditional approximation'' considerably simplifies this general equation, and we get two ODEs for the r- and $\theta$-components instead of one with arbitrary gradients of rotation  $\Omega(r,\theta)$. We have obtained a more stringent condition for the applicability of this approximation to the pulsation of stars. Only for very low frequencies is this restriction the same as that for the standard case. We also imposed some restrictions on the rotation profile  $\Omega(r,\theta)$ used. All of these restrictions have been compared with empirical solar rotation profiles and, thus, we justified their use in our modelling.

The $\theta$-equation is Laplace's equation generalized to a latitudinal differential rotation. Without solving this equation we found qualitatively the condition for the appearance of a global instability. This instability is driven by the latitudinal shear of rotation. It is not influenced by buoyancy, unlike the Kelvin-Helmholtz instability, and therefore it may easily be realized in stars, even if the horizontal gradient of rotation is small.

We have shown the possible appearance of critical latitudes at which there occurs a resonant interaction of eigenmodes with the inertial modes for a frequency $2\Omega\cos\theta=2\Omega_r$. The transformation of modes with a frequency $\omega=2\Omega_r$ in narrow latitude belts acts as a resonant damping mechanism. With decreasing mode frequency these belts are shifted toward the equator. This mechanism may play an important role in the redistribution of rotation angular momentum.

The appearance of mode instability strongly depends on the Rossby number, on the azimuthal wave numbers, and on the latitudinal rotation gradients. Very large gradients produce retrograde waves (seen in the rotating frame), while a slower rotation gradient is responsible for the prograde mode instability. The rotation gradient has a lower boundary below which an instability of the modes is impossible for any Rossby number or azimuthal number m.

We have applied the instability condition to helioseismological data. Here a global instability is possible for the m=1 mode at practically all latitudes. Radially the instability zone extends from the greater part of the tachocline up to the photosphere. The shear instability for the Sun was first obtained by Watson (1981). According to his results the instability is possible only in the photospheric layers. Later Gilman & Fox (1997) showed that such an instability is possible in the tachocline too, if strong toroidal magnetic fields are included. Our results show that the instability of the m=1 modes and other modes (m>1) are possible without magnetic fields, in contradiction to Gilman & Fox (1997; see also Charbonneau et al. 1999). This difference is probably connected with the incompleteness of the equations used by Watson (1981) and by Gilman & Fox (1997); their equations are two-dimensional only.

The exact solutions of Laplace's tidal equation for lower frequencies are expressed by Jacobi's polynomials. Just for lower frequencies the numerical calculations of stellar pulsation analysis meet great problems, when one is looking for the eigenfunctions as infinite series of Legendre functions. The eigenfunctions, defined by higher-order Jacobi polynomials, cannot be expressed by convergent series of associated Legendre functions. Every Legendre function is a particular case of a Jacobi polynomial.

It has been shown here that the retrograde (slow and fast) modes with high surface wave numbers (l,m) are energetically concentrated in narrow latitude bands. This analysis was done for the 22-year modes as an example. Such a concentration of mode energy in a narrow spatial area makes such modes vulnerable to different instability mechanisms such as the $\varepsilon$-mechanism considered in Paper 1.

All of our results have been obtained in the traditional approximation. This approximation does not work for waves propagating strictly within the equatorial plane. In our work this case is excluded. Motions which were considered in our paper can never cross the equator, and our boundary condition is $v_{\theta} \to 0$ when $\theta\to \pi/2$. In the paper we derived all the solutions which obey the traditional condition: high (l,m) - short waves and  $\vert v_r\vert\ll \vert v_{\theta,\phi}\vert$.

Appendix A: The main PDE for pressure perturbations

For the independent variable $\mu=\sin\theta$ and the definitions

$$\varepsilon_{{\mbox{\tiny R}}} = \frac{\omega}{2 \Omega} , \ \bet... ... r} , \ \beta_{\mu} = \frac{\mu}{2\Omega} \frac{\partial \Omega}{\partial \mu},$$

$$\alpha = \varepsilon^2_{{\mbox{\tiny R}}} - (1 - \mu^2)(1+ \beta_{\mu}),$$ (A.1)

$$a_1 = \frac{1 + \beta_r}{\varepsilon_{{\mbox{\tiny R}}}}, \ \ a_2... ...{{\mbox{\tiny R}}}}, \ \ a_3 = \frac{\alpha}{\varepsilon^2_{{\mbox{\tiny R}}}},$$

$$a_4 = \mu^2 \frac{1 + \beta_r}{\alpha} , \ \ a_5 = \frac{m}{\varepsilon_{{\mbox{\tiny R}}}} (1 + \beta_{\mu}) - 1;$$

$$% {\rm using} \ i V_{\phi} = j a_1 \mu V_r + a_2 V_{\theta} \... ...\theta - \frac{m}{\rho_0 i \omega} ~ \frac{p^{\prime}}{r \mu},$$ (A.2)

we get three equations from Eqs. (21)-(24) for the variables  V_r, $V_{\theta}$, and p'. From those $V_{\theta}$ is eliminated by

$$% \frac{V_{\theta}}{\cos\theta} = j a_4 \tilde{V}_r - \frac{1... ...ft( a_5 - \mu \frac{\partial}{\partial \mu} \right) \tilde{P},$$ (A.3)

$$% {\rm where} \ p^{\prime} = \mu r \rho_0 i \omega \tilde{P} , ~~ V_r = \mu \tilde{V}_r,$$ (A.4)

and we get two equations for $\tilde{V}_r$ and $\tilde{P}$:

  

$$\left( 1 - \frac{N^2}{\omega^2} - j a_4 \right) \tilde{V}_r =\lef... ... 1 - \frac{1}{a_3} \right) \mu \frac{\partial}{\partial \mu} \right] \tilde{P},$$ (A.5)

$$\left[ 2 + r \frac{\partial}{\partial r} + j \left( a_8 + a_9 \mu... ...partial \mu} \right) \right] \tilde{V}_r = \check{A}\left(\mu\right) \tilde{P},$$ (A.6)

$$\check{A} (\mu) = \frac{1}{\mu^2} \left[ a_7 + (1 - \mu^2) \mu \f... ...\left( a_5 - \mu \frac{\partial}{\partial \mu} \right) + \frac{m^2}{\mu^2}\cdot$$

Here the dimensionless coefficients are defined as

$$a_6 = \frac{1}{\alpha} \left[ m \varepsilon_{{\mbox{\tiny R}}} + ... ... \left(\frac{m}{\varepsilon_{{\mbox{\tiny R}}}} \beta_{\mu} - 1\right) \right],$$

$$a_7 = \frac{m}{\varepsilon_{{\mbox{\tiny R}}}} \left(1 + \beta_{\mu}\right) \left(1 - \mu^2\right) + 1 - 2 \mu^2,$$

$$a_8 = \frac{a^*_8}{\alpha^2} ; \ a_9 = \frac{1}{\alpha} \left(1 - \mu^2\right) \left(1 + \beta_r\right),$$ (A.7)

$$a_8^* = \left(1+\beta_r\right)\left[\alpha \left(3-4\mu^2 + m\varepsilon_{{\mbox{\tiny R}}}\right)\right.$$

$$\hspace*{7mm} \left. -\left(1-\mu^2\right)\check{\delta}_{\mu}\alpha\right]-2\alpha\left(1-\mu^2\right)\beta_r \beta_{\mu},$$

$$b_1 = 1 + \symbol{''1A}_{\rho} - \frac{m}{\varepsilon_{{\mbox{\ti... ... \symbol{''1A}_{\rho} = \frac{r}{\rho_0} ~\frac{{\rm d} \rho_0}{{\rm d} r}\cdot$$

In deriving these equations the second derivatives  $\Omega''(\theta)$ and  $\Omega''(r)$ have been omitted as very small quantities. Equation (A.5) gives one equation for p':

$$% \left[\psi_1 \check{\delta}^2_r + \psi_2 \check{\delta}^2_{... ...\delta}_r \check{\delta}_{\mu} + \psi_6 \right] \tilde{P} = 0.$$ (A.8)

This is the main singular PDE for nonradial rotation-gravity waves in a differentially rotating star. Here the operators are defined as

$$% \check{\delta}_r =r\frac{{\partial}{}}{{\partial}{r}} , \ \... ...m and} \ \check{\delta}^2_r =\check{\delta}_r\check{\delta}_r,$$ (A.9)

and the coefficients $\psi_{1-6}$ are:

$$\begin{eqnarray*}&& \psi_1 = 1-\frac{N^2}{\omega^2} - j\mu^2\frac{1+\beta_r}{\al... ...ck{\delta}_{\mu}\left(\frac{a_5}{a_3}\right)\right], {\rm where} \end{eqnarray*}$$

$$\begin{eqnarray*}&& f_1 = \psi_1\left(2+j a_8\right)-\check{\delta}_r \psi_1 - j... ... \left(\mu^2 -1\right)\left(1+\beta_{\mu}\right)/\alpha=1-1/a_3. \end{eqnarray*}$$

Taking into account $\nabla\omega=-m\nabla\Omega$ we can easily obtain all the required derivations of the parameters.

Acknowledgements

The authors gratefully acknowledge the critical and very different, but constructive comments by two anonymous referees who helped to improve the paper. Moreover, we would like to thank George Isaak and Jet Katgert for careful reviews of the paper and suggestions for improving the language. The present work has been financially supported by the German Science Foundation (DFG) under grants Nos. 436 RUS 113/560/4-1 and 436 RUS 113/689/2-1 and by the Russian Foundation of Basic Research under RFBR No. 04-02-16386.

References

• Bateman, H., & Erdélyi, A. 1953, Higher Transcendental Functions (New York: McGraw-Hill) (In the text)
• Berthomieu, G., Gonczi, G., Graff, P., Provost, J., & Rocca, A. 1978, A&A, 70, 597 [NASA ADS] (In the text)
• Bildstein, L., Ushomirsky, G., & Cutler, C. 1996, ApJ, 460, 827 [NASA ADS] [CrossRef] (In the text)
• Cally, P. S. 2001, Sol. Phys., 199, 231 [NASA ADS] [CrossRef] (In the text)
• Charbonneau, P., Tomczyk, S., Schou, J., & Thompson, M. J. 1998, ApJ, 496, 1015 [NASA ADS] [CrossRef] (In the text)
• Charbonneu, P., Dikpati, M., & Gilman, P. 1999, ApJ, 526, 523 [NASA ADS] [CrossRef] (In the text)
• Clement, M. J. 1998, ApJ, 116, 57 [NASA ADS] [CrossRef] (In the text)
• Cox, J. P. 1980, Theory of Stellar Pulsation (Princeton, USA: Princeton Univ. Press) (In the text)
• Dzhalilov, N. S., Staude, J., & Arlt, K. 2000, A&A, 361, 1127 [NASA ADS] (In the text)
• Dzhalilov, N. S., Staude, J., & Oraevsky, V. N. 2002, A&A, 384, 282 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Dzhalilov, N. S., & Zhugzhda, Y. D. 1990, AZh, 67, 561 [NASA ADS] (In the text)
• Dziembowsky, W., & Kosovichev, A. 1987, Acta Astron., 37, 341 [NASA ADS] (In the text)
• Eckart, C. 1960, Hydrodynamics of Oceans and Atmospheres (Oxford, England: Pergamon Press) (In the text)
• Erdélyi, A. 1942, Duke Math. J., 9, 48 [NASA ADS] [MathSciNet] (In the text)
• Gill, A. E. 1982, Atmosphere-Ocean Dynamics (England: Acad. Press, Univ. of Cambridge) (In the text)
• Gilman, P. A., & Fox, P. A. 1997, ApJ, 484, 439 [NASA ADS] [CrossRef] (In the text)
• Heun, K. 1889, Math. Annalen, 33, 161 [NASA ADS] (In the text)
• Hollweg, J. W. 1987, ApJ, 312, 880 [NASA ADS] [CrossRef] (In the text)
• Ince, E. L. 1927, Ordinary Differential Equations (London) (In the text)
• Ionson, J. 1978, ApJ, 226, 650 [NASA ADS] [CrossRef] (In the text)
• Ionson, J. 1984, ApJ, 276, 357 [NASA ADS] [CrossRef] (In the text)
• Knobloch, E., & Spruit, H. C. 1982, A&A, 113, 261 [NASA ADS] [MathSciNet] (In the text)
• Laplace, P. S. 1778-1779, Recherches sur Plusieurs Points du Système du Monde (extracts from Mémoires de L'Académie Royale des Sciences, 1775, 1776)
• Lee, U., & Saio, H. 1986, MNRAS, 221, 365 [NASA ADS] (In the text)
• Lee, U., & Saio, H. 1997, ApJ, 491, 839 [NASA ADS] [CrossRef] (In the text)
• Ledoux, P. 1951, ApJ, 114, 373 [NASA ADS] [CrossRef] (In the text)
• Levin, Y., & Ushomirsky, G. 2001, MNRAS, 322, 515 [NASA ADS] [CrossRef] (In the text)
• Lindzen, R. S., & Chapman, S. 1969, Space Sci. Rev., 10, 3 [NASA ADS] (In the text)
• Longuet-Higgins, M. S. 1965, Proc. R. Soc. London A, 284, 40 [NASA ADS] [MathSciNet] (In the text)
• Lou, Y. Q. 2000, ApJ, 540, 1102 [NASA ADS] [CrossRef] (In the text)
• Oraevsky, V. N., & Dzhalilov, N. S. 1997, Astron. Rep., 41, 91 [NASA ADS] (In the text)
• Papaloizou, J., & Pringle, J. 1978, MNRAS, 182, 423 [NASA ADS] (In the text)
• Pedlosky, J. 1982, Geophysical Fluid Dynamics (New York-Heidelberg-Berlin: Springer-Verlag) (In the text)
• Provost, J., Berthomieu, G., & Rocca, A. 1981, A&A, 94, 126 [NASA ADS] [MathSciNet] (In the text)
• Rayleigh, J. W. S., Lord 1880, Proc. Lond. Math. Soc., 11, 57 [NASA ADS] (In the text)
• Saio, H. 1982, ApJ, 256, 717 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Shore, S. N. 1992, An Introduction to Astrophysical hydrodynamics (Acad. Press INC), 452 (In the text)
• Smeyers, P., Craeynest, D., & Martens, L. 1981, Ap&SS, 78, 483 [NASA ADS] (In the text)
• Stix, M., & Skaley, D. 1990, A&A, 232, 234 [NASA ADS]
• Townsend, R. H. D. 1997, MNRAS, 284, 839 [NASA ADS] (In the text)
• Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars (Univ. Tokio Press), 420 (In the text)
• Watson, M. 1981, Geophys. Astrophys. Fluid Dyn., 16, 285 [NASA ADS] (In the text)
• Whittaker, E. T., & Watson, G. N. 1927, A Cource of Modern Analysis (Cambridge) (In the text)
• Wolff, C. L. 1998, ApJ, 502, 961 [NASA ADS] [CrossRef] (In the text)