A&A 427, 279-292 (2004)
DOI: 10.1051/0004-6361:20040539

Oscillations of magnetic stars

II. Axisymmetric toroidal and non-axisymmetric shear Alfvén modes in a spherical shell

D. Reese - F. Rincon - M. Rieutord

Laboratoire d'Astrophysique de Toulouse et Tarbes, Observatoire Midi-Pyrénées, 14 avenue É. Belin, 31400 Toulouse, France

Received 29 March 2004 / Accepted 16 July 2004

Abstract
We carry out numerical and mathematical investigations of shear Alfvén waves inside of a spherical shell filled with an incompressible conducting fluid, and bathed in a strong dipolar magnetic field. We focus on axisymmetric toroidal and non-axisymmetric modes, in continuation of a previous work by Rincon & Rieutord (2003, A&A, 398, 663). Analytical expressions are obtained for toroidal eigenmodes and their corresponding frequencies at low diffusivities. These oscillations behave like magnetic shear layers, in which the magnetic poles play a key role, and hence become singular when diffusivities vanish. It is also demonstrated that non-axisymmetric modes are split into two categories, namely poloidal or toroidal types, following similar asymptotic behaviours as their axisymmetric counterparts when the diffusivities become arbitrarily small.

Key words: magnetohydrodynamics (MHD) - stars: oscillations - stars: magnetic fields

1 Introduction

Numerous astrophysical systems exhibit a pulsating behaviour that can be significantly influenced by the Lorentz force when a strong magnetic field is present. This may for instance be the case in neutron stars and magnetic white dwarfs (Lou 1995). Planetary cores, which are known to sustain a strong dynamo (Stevenson 1983), are also likely to fall into this category.

One of the most exciting examples of couplings between pulsation and magnetism is given by the seismological activity of roAp stars. This class of stars, discovered by Kurtz (1978), exhibits several kG (almost) dipolar magnetic fields and luminosity variations on periods ranging from 5 to 15 min. These oscillations seem to be well approximated by a single spherical harmonic $\ell= 1$ lined up with the magnetic axis, suggesting a strong mixing between high order p-modes and Alfvénic type oscillations.

Many different models have been developed to obtain a satisfying picture of the asteroseismology of these stars. Following theoretical work by Roberts & Soward (1983); Biront et al. (1982) and Campbell & Papaloizou (1986), Dziembowski & Goode (1996) have studied acoustic star models enveloped by a layer in which magnetic effects become dominant. Using a boundary layer approximation, they came up with an outer boundary condition which was then applied for the calculation of adiabatic acoustic modes. This model has undergone a lot of refinements in order to take into account new physical processes. The latest version, proposed by Bigot & Dziembowski (2002), incorporates the centrifugal force (a non-axisymmetric effect, since the rotational and magnetic axis are often tilted in roAp stars) and suggests that the axis of the modes is not necessarily lined up with the magnetic axis. In spite of these improvements, there are still non-negligible discrepancies between the magnetically shifted eigenfrequencies computed from these models and the observed ones, showing that the precise coupling mechanism occurring in the surface layers is likely to be more complex. An important point is that a single $\ell$ value is sometimes assumed to be sufficient to describe the oscillations. This may not be the case, owing to the dipolar structure of the permanent magnetic field which induces a coupling between spherical harmonics thus producing a whole spectrum of $\ell$'s (e.g. Rincon & Rieutord 2003).

Motivated by the observation that chemical peculiarities are observed near the magnetic poles of roAp stars, Balmforth et al. (2001) have tried to determine what precise physical phenomena were occurring in the polar and equatorial regions. Since the magnetic field is almost horizontal near the equator and vertical near the poles, convection is certainly inhibited in the latter region, allowing the diffusion of different chemical elements (an excess of helium is observed on the polar caps). Oscillations triggered by a $\kappa$-mechanism may therefore preferentially be observed in this region. This approach stresses the importance of a global description of the eigenmodes in such stars.

The work presented here originates partly from the preceding remarks. It continues the study described in Rincon & Rieutord 2003 (hereafter referred to as Paper I) and should be viewed as a preliminary step to obtain global models of the more complex magneto-gravito-acoustic oscillations. We aim at understanding some of the underlying physical mechanisms potentially involved in roAp stars or similar astrophysical objects, by using an approach that rigorously treats the couplings induced by the geometry of the magnetic field. Our highly simplified model consists of a non-rotating spherical shell of incompressible magnetised fluid bathed in a dipolar magnetic field, with small magnetic and kinematic diffusivities. As a consequence of incompressibility, the modes we compute are shear Alfvén waves. In Paper I, a spherical harmonic decomposition of the linearised MHD equations was obtained and results regarding poloidal axisymmetric modes were presented (a short reminder of the classification of modes is given later on in the paper). It was shown that the least-damped modes are near the magnetic poles and exhibit internal shear layers which can potentially play a role in mode selection. In the present study, we focus on the other types of shear Alfvénic oscillations in spherical shells, namely axisymmetric toroidal and non-axisymmetric modes, and characterise numerically and mathematically their geometry, eigenspectrum, and some by-products like boundary layers.

In Sect. 2, we recall the basic physical ingredients of the model and shortly describe our numerical strategy. Section 3 is devoted to the phenomenological description and analytical study of axisymmetric toroidal modes. Section 4 covers non-axisymmetric modes and their resemblance to their axisymmetric counterparts. Finally, Sects. 5 and 6 conclude the paper. Note that most of the mathematical details involved in Sect. 3 are given in Appendices A and B to preserve the clarity of the manuscript.

2 Description of the model

2.1 Basic equations

We first give a brief description of the model that is used in calculating toroidal and non-axisymmetric modes and which was already used in Paper I. More details on the basic equations and their expansion onto the spherical harmonic base are given in Paper I.

The "star'' we work with is a spherical shell of incompressible plasma of density $\rho_o$, with a radius R and an aspect ratio $\eta$. In all the numerical examples and figures, we use $\eta = 0.35$. Within this shell is a dipolar magnetic field generated in a perfectly conducting core:

$$\vec{B} = B_o \. R^3 \left( \frac{\cos \theta}{r^3} \vec{e}_r + \frac{\sin \theta} {2r^3} \vec{e}_{\theta}\right).$$ (1)

We use the following underlined dimensionless variables:

$$\vec{r} = R\underline{\vec{r}}, \quad \vec{v} = V_{\rm A} \un... ...underline{\vec{B}}, \quad t = \frac{R}{V_{\rm A}}\underline{t},$$

where $V_{\rm A}$, the Alfvén velocity, is given by the following expression, $\mu_o$ being the magnetic permeability:

$$V_{\rm A} = \frac{B_o}{\sqrt{\rho_o \mu_o}}\cdot$$ (2)

We apply the linearised magneto-hydrodynamic equations to the kinetic $\vec{v}$ and magnetic $\vec{b}$ perturbations (we drop the underlined notation). We assume these perturbations have a time dependence of the form ${\rm e}^{\lambda t}$, where $\lambda = \tau + {\rm i} \omega$ ($\tau$ is the damping rate, $\omega$ the pulsation and i²=-1). This leads to the following set of equations:

 

$$\vec{\nabla}\cdot\vec{v} = 0,$$

$$\vec{\nabla}\cdot\vec{b} = 0,$$

$$\lambda \vec{\nabla}\times\vec{v} = \vec{\nabla}\times\left( (\ve... ...abla}\times\vec{b}) \times \vec{B} \right) + E \vec{\nabla}\times\Delta\vec{v},$$

$$\lambda \vec{b} = \vec{\nabla}\times\left(\vec{v} \times \vec{B} \right) + E_{{\rm m}}\Delta\vec{b}$$ (3)

where $\vec{B}$ denotes the permanent dipolar magnetic field. The parameters E and $E_{{\rm m}}$ are non-dimensional forms of the kinematic and magnetic diffusivities, respectively, and are given by the following expressions:

$$E = \frac{\nu}{RV_{\rm A}}, \quad E_{{\rm m}}= \frac{1}{\sigma_o \mu_o R V_{\rm A}},$$ (4)

where $\sigma_o$ is the conductivity of the fluid. In the case of roAp stars, these parameters take on the following typical values:

$$E_{{\rm m}}\sim 10^{-8}, \qquad E \sim 10^{-13}.$$

2.2 Boundary conditions

The boundary conditions used in our model are as follows. For the velocity, we use stress-free, non-penetrative, conditions to minimise the importance of boundary layers. This leads to the following conditions, valid for $r = \eta$and r = 1:

$$\frac{\partial }{\partial r} \left( \frac{v_{\theta}}{r} \rig... ... r} \left( \frac{v_{\varphi}}{r} \right) = 0, \qquad v_r = 0.$$ (5)

As for the magnetic field, different conditions apply for the inner and outer boundaries. On the interior, the perturbation to the electric field is perpendicular to the conducting core, and the perturbation to the magnetic field is tangent. This gives the following three equations:

 

$$\frac{E_{{\rm m}}}{r} \frac{\partial }{\partial r}(rb_{\theta}) = -v_{\theta} B_r,$$

$$\frac{E_{{\rm m}}}{r} \frac{\partial }{\partial r}(rb_{\varphi}) = -v_{\varphi}B_r.$$ (6)

When the other boundary conditions are taken into account, and when $\lambda$ is different from zero, these equations are not independent and only correspond to two conditions. On the outer boundary and beyond, the magnetic field is continuous and potential, since there are no currents in empty space. By using the continuity of this field, and the idea that it does not diverge toward infinity, it is possible to define two boundary conditions, which are more easily expressed in spherical harmonics (see Eq. (8)).

Equation (3), together with these boundary conditions, defines a generalised eigenvalue problem, where $\lambda$ is the eigenvalue and $(\vec{v}, \vec{b})$ is the eigenvector which can be computed numerically.

2.3 Harmonic projection

In order to solve Eq. (3), the fields $\vec{v}$ and $\vec{b}$are projected onto the harmonic base:

 

$$\vec{v} = \displaystyle \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\e... ...{\ell}\vec{R}_{\ell}^m+ v_m^{\ell}\vec{S}_{\ell}^m+ w_m^{\ell}\vec{T}_{\ell}^m,$$

$$\vec{b} = \displaystyle \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\e... ...{\ell}\vec{R}_{\ell}^m+ b_m^{\ell}\vec{S}_{\ell}^m+ c_m^{\ell}\vec{T}_{\ell}^m,$$ (7)

in which $\vec{R}_{\ell}^m$, $\vec{S}_{\ell}^m,$ and $\vec{T}_{\ell}^m$ are the normalised spherical harmonics:

$$\vec{R}_{\ell}^m= Y_{\ell}^m\vec{e}_r, \quad \vec{S}_{\ell}^m... ... \quad \vec{T}_{\ell}^m= r \vec{\nabla}\times\vec{R}_{\ell}^m.$$

The full harmonic decomposition of Eq. (3) is given in Appendix A of Paper I. The outer boundary condition on the magnetic field reads:

 

$$\frac{\partial a_m^{\ell}}{\partial r} + \frac{(\ell+2)a_m^{\ell}}{r} = 0,$$

$$c_m^{\ell} = 0.$$ (8)

2.4 Classification and symmetries

The various eigenmodes fit naturally into different categories. Firstly, as was already shown in Paper I, there is no coupling between different m's. Hence, eigenmodes will be made up of only one m. This leads to two types of modes: axisymmetric oscillations (m=0) and non-axisymmetric ones $(m \neq 0)$. Secondly, within the axisymmetric category, it is possible to distinguish between poloidal modes and toroidal ones, as the corresponding equations fully decouple for m=0 (see Appendix A of Paper I). Poloidal modes are made up of $u_0^{\ell}$, $v_0^{\ell}$, $a_0^{\ell}$, and $b_0^{\ell}$ functions which correspond to $\vec{e}_r$ and $\vec{e}_{\theta}$ components. Toroidal modes are made up of $w_0^{\ell}$ and $c_0^{\ell}$functions and are in the $\vec{e}_{\varphi}$ direction.

A certain number of symmetries are present in the physical system and lead to a few simplifications. As was pointed out in Paper I, a parity can be defined for eigenmodes. However, there was a slight confusion as to the parity of toroidal components (see Eq. (19) of Paper I) since a mode is called even when the velocity perturbation is symmetric with respect to the equator and the magnetic perturbation antisymmetric, and a mode is odd in the reverse situation. The corrected form for toroidal eigenvectors reads:

$$(\vec{v},\vec{b}) = \left(w_o^{2k+1},c_o^{2k}\right)\quad \mbox{even,}$$

$$(\vec{v},\vec{b}) = \left(w_o^{2k},c_o^{2k+1}\right) \quad \mbox{odd.}$$ (9)

Non-axisymmetric modes also have a global parity. In fact it is possible to anticipate this result without the use of spherical harmonics. If we define $\mathcal{M}$ as being the operator that gives the mirror image of a vector field with respect to the equator, and if $(\vec{v},\vec{b},\lambda)$ is a solution, then $(\mathcal{M}\vec{v},-\mathcal{M}\vec{b},\lambda)$ is also a solution. By combining these two solutions, it is possible to extract an even part and an odd part from the original solution. Eigenvectors will take on the following form:

$$(\vec{v},\vec{b})=\left(u_m^{m+2k},w_m^{m+2k+1},a_m^{m+2k+1},c_m^{m+2k}\right)\quad\mbox{even,}$$

$$(\vec{v},\vec{b})=\left(u_m^{m+2k+1},w_m^{m+2k},a_m^{m+2k},c_m^{m+2k+1}\right)\quad\mbox{odd.}$$ (10)

This symmetry enables us to work with half as many components for a given resolution, which is advantageous from the numerical point of view. Useful information about the eigenspectra can be deduced from the symmetry of $\vec{B}$ with respect to any meridian. For a given solution $(\vec{v},\vec{b},\lambda)$ it is possible to create a second solution $(\mathcal{S}\vec{v},\mathcal{S}\vec{b}, \lambda)$, where $\mathcal{S}$ is the operator that gives the mirror image of a vector field with respect to the meridian that passes through $\varphi = 0$. When applied to spherical harmonics, the azimuthal order m changes to -m. From this we conclude that the eigenspectrum of modes with azimuthal order m is the same as that of -m. This also leads to the decoupling of poloidal and toroidal components in the axisymmetric case (the poloidal modes being symmetric with respect to $\mathcal{S}$ and the toroidal modes antisymmetric). By taking into account that $\vec{B}$ is real, we deduce that the eigenspectrum for m is the conjugate of that for -m. Therefore, each eigenspectrum is symmetrical with respect to the real axis. Practically, this means that we only need to explore eigenvalues for positive azimuthal orders, and only on the upper half of the complex plane $(\omega \geq 0)$.

If rotation were taken into account, a number of these symmetries would break down. Apart from the situation where the rotation and magnetic axis are aligned, different m's become coupled and mode parity is lost. The symmetry with respect to  $\mathcal{S}$ is lost even when the rotation axis is lined up with the magnetic field.

2.5 Numerical aspects

Eigenmodes and eigenvalues are calculated numerically using two different methods. The first method, based on a QZ algorithm, gives all the eigenvalues (for the discretised problem) whereas the second method is an iterative Arnoldi-Chebyshev algorithm which only computes a selection of eigenvalues and their corresponding eigenmodes.

In our numerical calculations, we use a simplified version of Eq. (7). Since different m's are decoupled, we do not have a summation over the azimuthal order. Also, because the divergence of both perturbations is zero, the functions  $v_m^{\ell}$and  $b_m^{\ell}$ can be expressed in terms of $u_m^{\ell}$ and $a_m^{\ell}$, respectively (see Appendix A of Paper I). Furthermore, in the case of axisymmetric modes, the summation on $\ell$ will actually start at 1 instead of 0as a result of our boundary conditions. In the non-axisymmetric case, the summation on $\ell$ will start at |m| as expected. Finally, the sum on the spherical harmonics is truncated at L. This leads to the following formulas:

$$\vec{v} = \displaystyle \sum_{\ell=\ell_{min}}^{L} u_m^{\ell}\vec{R}_{\ell}^m+ v_m^{\ell}\vec{S}_{\ell}^m+ w_m^{\ell}\vec{T}_{\ell}^m,$$

$$\vec{b} = \displaystyle \sum_{\ell=\ell_{min}}^{L} a_m^{\ell}\vec{R}_{\ell}^m+ b_m^{\ell}\vec{S}_{\ell}^m+ c_m^{\ell}\vec{T}_{\ell}^m.$$ (11)

Each of the functions $u_m^{\ell}$, $w_m^{\ell}$, $a_m^{\ell}$ and $c_m^{\ell}$ is written in the form of a truncated Chebyshev series:

$$u_m^{\ell}(r) = \sum_{k=0}^{Nr} \tilde{u}_m^{\ell}(k) T_k(r).$$

Typically, we could reach a spatial resolution of L=350, Nr=930 for axisymmetric toroidal modes (when analysing boundary layers) and L=858, Nr=340 for non-axisymmetric modes.

We also define Chebyshev and harmonic spectral coefficients:

$$C_k = \frac{\displaystyle \max_{\ell} \vert\tilde{u}_m^{\ell}... ...playstyle \max_{\ell,k} \vert\tilde{u}_m^{\ell}(k)\vert}\cdot$$

They are useful for characterising certain aspects of a mode's structure, as will be seen in Sect. 3.7, as well as the convergence of the discretization.

   
3 Axisymmetric toroidal modes

3.1 Basic properties

In this section, we discuss the basic properties of axisymmetric toroidal modes. As a reminder, the fields $\vec{v}$ and $\vec{b}$ are both orientated in the direction $\vec{e}_{\varphi}$. From the point of view of spherical harmonics they are made up of $(w_0^{\ell},c_0^{\ell})$ components.

The equations that govern toroidal modes can be simplified. It is no longer necessary to take the curl of the Navier-Stokes equation to remove the pressure term, because this term vanishes owing to axisymmetry. This leads to the following, nearly symmetric system:

 

$$\lambda b_{\varphi} = (\vec{B} \. \vec{\nabla} ) v_{\varphi} - 3\cos\theta/2r^4v_\varphi + E_{{\rm m}}\Delta' b_\varphi,$$

$$\lambda v_{\varphi} = (\vec{B} \. \vec{\nabla} ) b_{\varphi} + 3\cos\theta/2r^4b_\varphi + E \Delta' v_\varphi.$$ (12)

Note that $\vec{\nabla}\cdot\vec{v} = 0$, $\vec{\nabla}\cdot\vec{b} = 0$ and $\Delta' = \Delta-1/r^2\sin^2\theta$.

3.2 Eigenvalue spectrum

In many ways, the eigenvalue spectrum of toroidal modes is similar to that of the poloidal modes presented in Paper I. The complex eigenvalues are grouped into "horizontal'' branches, as can be seen in Fig. 1. These "eigenbranches'', indexed by n, correspond to modes with n nodes in the radial direction (see Fig. 2). Both perturbations have n nodes, but in different positions. Depending on parity, diffusion modes $(\omega = 0,~ \tau \neq 0)$ also appear but these modes are of little interest to asteroseismology. In all that follows, we will only consider oscillatory modes.

	Figure 1: Eigenvalue spectrum for odd modes and $E = E_{{\rm m}}= 10^{-4}$. The plusses (+) correspond to oscillatory modes and the diamonds ($\diamond$) to diffusion modes, which only appear for the odd parity. The letters n and qquantify the oscillatory mode spectrum.
Open with DEXTER

Figure 2: Axisymmetric toroidal modes with different vertical structures, for $E=E_{{\rm m}}=10^{-3}$. The left quadrant in each image corresponds to magnetic dissipation and the right one to magnetic energy. From left to right, and from top to bottom, the modes represented here correspond to (n,q)=(2,0), (3,0), (4,0) and (5,0). A logarithmic scale is used in this figure and most of the other figures in this paper, as it brings out more details on the modes' structure.

Open with DEXTER

A careful look at the eigenvalues along a single branch reveals a regular structure. These values are remarkably well lined up, much better than their poloidal counterparts, and they are equally spaced out. This suggests the following empirical law:

$$\lambda_{n,q} = \alpha_n + q \cdot \beta_n$$ (13)

where $\alpha_n$ and $\beta_n ~ (<$0) depend on the value of n, and q is given in Fig. 1. By looking at the slope of the eigenbranches, it is possible to determine the ratio between the real and imaginary parts of $\beta_n$. In our case, ${\rm Re} (\beta_n) \simeq {\rm Im} (\beta_n)$.

	Figure 3: Different modes from the n=3 eigenbranch. From top to bottom, qtakes on the values 1, 3 and 9.
Open with DEXTER

We can then look at the corresponding eigenmodes. When q increases, two phenomena appear. First of all, for the least-damped modes, the number of nodes in the horizontal direction seems to increase. However, these nodes are not all conspicuous. Secondly, modes that have a higher damping rate gradually depart from the magnetic poles. Figure 3 illustrates both of these phenomena. Care must be taken with highly damped modes: the spectral coefficients indicate that the solutions are convergent, but the corresponding eigenvalues are highly unstable numerically due to round-off errors. This limits the number of modes that can be analysed safely.

3.3 Mode structure and resonant field lines

The structure of the modes can be understood if one computes the resonance frequencies of different field lines for the ideal case $E=E_{{\rm m}}=0$. As each field line can oscillate at its own frequency, the resonance frequencies form a continuous spectrum called the Alfvén continuum. (For our configuration, there will be a continuum for each value of n.) These eigenvalues are purely imaginary because the ideal MHD operator which gives the squared eigenvalues $\lambda^2$ is self-adjoint.

For a first approximation of these frequencies, we can apply the same formula as was used by Alfvén (1945), which is derived from the time it takes for a wave traveling along a field line to go from one boundary to the other. This corresponds to a WKB-type approximation, in which the curvature terms have been neglected:

$$T = \int_{{\rm Field~line}} \frac{{\rm d}l}{\Vert \vec{V}_{... ...eft( n + \frac{1}{2} \right) \frac{\pi}{\omega_{{\rm f.l.}}},$$ (14)

where n is the number of nodes along the field line. A field line will then be resonant if it satisfies the relation $\omega_{{\rm f.l.}} = \omega$where $\omega$ is the frequency of the mode. Since the density is constant, the integral can be calculated analytically and is given by:

$$\int \frac{{\rm d}l}{\Vert \vec{V}_{\!\! {\rm A}} \Vert} = ... ...if} \quad \theta_{\eta} = \theta_1 = 0, \end{array} \right.$$ (15)

with $[f(\theta)]_{\theta_{\eta}}^{\theta_1} =f(\theta_1) - f(\theta_{\eta})$. $\theta_{\eta}$ and $\theta_1$ are the colatitudes of the field line on the inner and outer boundaries and satisfy the relation $\sin^2 \theta_{\eta} = \eta \sin^2 \theta_1$. This formula is the same as Eq. (26) by Kato & Watanabe (1958), except for the bounds on the integral. We see that the oscillation rate decreases for field lines further from the magnetic poles, hence one can expect eigenmodes with lower eigenfrequencies to be further from the poles.

Taking the variations of $\vec{B}$ into account, in order to obtain better results, yields the following equations, which still only apply to individual field lines:

 

$${\rm i}\omega_{{\rm f.l.}} b = \displaystyle \frac{\cos \theta}{r^3} \frac{{\rm d}v}{{\rm d}r} -\frac{3\cos \theta}{2 r^4} v,$$

$${\rm i}\omega_{{\rm f.l.}} v= \displaystyle \frac{\cos \theta}{r^3} \frac{{\rm d}b}{{\rm d}r} +\frac{3\cos \theta}{2 r^4} b,$$

$$v(\eta) = 0,$$

b(1) = 0 (16)

where $\cos \theta = \sqrt{1 - r \sin^2 \theta_1}$. The boundary conditions are derived from a simple analytical model (see Eq. (B.9)) and were used to establish the n+1/2quantization in Eq. (14). With this approach, the resonant field line coincides much better with the mode (see Fig. 4).

Figure 4: The kinetic energy and dissipation of eigenmode (n,q) = (1,8) with resonant field lines superimposed. The diffusivities are $E=E_{{\rm m}}=4\times 10^{-5}$. The dotted line corresponds to Eq. (14), which is a WKB-type approximation, and the dashed one to Eq. (16), which takes the variations of $\vec{B}$ into account. As opposed to other figures in this paper, a linear intensity scale is used instead of a logarithmic one. This gives a thinner appearance to the mode's structure.

Open with DEXTER

   
3.4 Polar eigenspectrum

We will now give a precise calculation (in the adiabatic case $(E=E_{{\rm m}}=0)$) of the resonance frequencies corresponding to the field line along the magnetic axis. This analysis is motivated by the role of resonating polar field lines in the asymptotic limit of small diffusivities. The adiabatic eigenfrequencies will be denoted by the superscript "0''.

These frequencies are solutions of Eq. (16) with $\cos \theta = 1$. By combining the two first equations and solving for v, we obtain:

$$r^2 \frac{{\rm d}^2 v}{{\rm d}r^2} - 3 r \frac{{\rm d}v}{{\rm d}r} + \left( \frac{15}{4} + \omega^2 r^8 \right) v= 0.$$ (17)

(We have dropped the notation $\omega_{{\rm f.l.}}$.) The solution to this equation is (see Abramowitz & Stegun (1972)):

$$v= Ar^2 J_{ \frac{1}{8}} \left( \frac{\omega r^4}{4} \right) + Br^2 J_{-\frac{1}{8}} \left( \frac{\omega r^4}{4} \right), \\$$ (18)

where A and B are two constants and $J_{\nu}$ is the Bessel function of the first kind of order $\nu$. The corresponding solution for b is:

$$b = {\rm i}Br^2 J_{ \frac{7}{8}} \left( \frac{\omega r^4}{4} ... ...Ar^2 J_{-\frac{7}{8}} \left( \frac{\omega r^4}{4} \right)\cdot$$ (19)

The boundary conditions determine the eigenfrequencies, via the following relation:

$$J_{ \frac{1}{8}} \left( \frac{\omega {\eta}^4}{4} \right) ... ...\right) J_{-\frac{7}{8}} \left( \frac{\omega}{4} \right) = 0.$$ (20)

We will use the notation $\lambda_n^0 = {\rm i} \omega_n^0$ to mean the nth eigenvalue calculated with this model. Asymptotically, for high frequencies, these values are well approximated by the WKB frequencies found in Paper I: $\omega_n^0 \sim \frac{2(2n+1)\pi}{1-\eta^4}$. The corresponding solutions are written b_n⁰ and v_n⁰ and are given by the formulas:

 

$$r^2 J_{ \frac{7}{8}} \left( \frac{\omega^0_n}{4} \right) J_{ \fra... ...omega^0_n}{4} \right) J_{-\frac{1}{8}} \left( \frac{\omega^0_n r^4}{4} \right),$$

$${\rm i}r^2 J_{-\frac{7}{8}} \left( \frac{\omega^0_n}{4} \right) J... ...a^0_n}{4} \right) J_{-\frac{7}{8}} \left( \frac{\omega^0_n r^4}{4} \right)\cdot$$ (21)

It is worth noting the striking similarity between these expressions and Eq. (30) by Kato & Watanabe (1958), which corresponds to the radial variation of geomagnetic poloidal modes (instead of toroidal ones) calculated for a dipolar background in which the $\theta$-dependence of $\Vert \vec{B} \Vert$ has been neglected.

In order to show the physical meaning of these eigenvalues, we took an eigenspectrum and prolongated the lines formed by the eigenbranches (see Fig. 1). The intersections between the lines and the imaginary axis $(\tau = 0)$ correspond to what could be called "numerical polar eigenfrequencies''. The agreement between these values and $\lambda ^0_n$ is rather good, as can be seen in Table 1. This leads to the conclusion that the polar eigenfrequencies are a good indicator of the eigenbranches' positions.

Table 1: A comparison between numerical polar eigenfrequencies (Num.), adiabatic ones ( $\omega ^0_n$) and a WKB approximation. The numerical values were obtained using a least square method on the eigenvalues corresponding to $E = E_{{\rm m}}= 10^{-4}$ and correspond to the intersection between the eigenbranches and the imaginary axis.

   
3.5 Asymptotic behaviour for small diffusivities

The next step is to look at the behaviour of eigenmodes and eigenvalues when the diffusivities E and $E_{{\rm m}}$ approach zero. There are two important reasons for this. First of all, realistic astrophysical values of E and $E_{{\rm m}}$ are very small ($\la$10^-8) and still out of reach for numerical solvers. Being able to extrapolate the behaviour of eigenmodes gives an educated guess as to what they would actually be for low diffusivities. A second reason is that it may be possible this way to refine analytically the adiabatic eigenvalues $\lambda_n^0$ defined by Eq. (20) to obtain a better approximation of the true eigenvalues.

We therefore analysed the behaviour of the mode n=1, q=0 for different values of the diffusivities. Figures 5 and 6 show how its eigenvalue, its position and its "thickness'' vary. In order to calculate the mode's position and thickness, we took a profile of the magnetic energy $\vert \vec{b} \vert^2$ along a meridional cut with a fixed radius of 0.5 and a horizontal extent going from $\theta=0$ to $\theta=\pi/2$. From this profile, we calculated a mean value $\overline {\theta }$ and a standard deviation $\sigma _{\theta }$. These correspond to estimates of the mode's position and thickness, respectively.

Figure 5: Behaviour of eigenvalue (n,q)=(1,0) for different values of E and $E_{{\rm m}}$. $\Delta \omega$ corresponds to the difference between the theoretical polar eigenfrequency, solution to Eq. (20), and the actual frequency of the mode. When $E_{{\rm m}}=10^{-3}$, the eigenvalue never attains $\lambda ^0_n$. When $E=E_{{\rm m}}$, the slope of the $\Delta \omega$ line is 1/2, meaning that $\Delta \omega$ is proportional to E1/2.

Open with DEXTER

	Figure 6: Behaviour of an eigenmode for different values of E and $E_{{\rm m}}$. $\overline {\theta }$ is an estimate of the mode's position and $\sigma _{\theta }$an estimate of the thickness. In the case $E=E_{{\rm m}}$ the position and thickness seem to be approaching zero; the slope of both lines is 1/4.
Open with DEXTER

From these graphs it is possible to deduce some of the physical phenomena that are taking place. In the case where $E=E_{{\rm m}}$, we observe that the eigenvalue approaches the n=1 polar eigenfrequency given by Eq. (20). At the same time, the mode seems to get thinner and closer to the magnetic poles. We observe from Figs. 5 and 6 that:

 

$$\lambda_{1,0} = \lambda^0_1 + E^{1/2} (\tau^1 + {\rm i} \omega^1) + E \tau^{2},$$

$$\overline{\theta} \propto E^{1/4},$$

$$\sigma_{\theta} \propto E^{1/4},$$ (22)

in which $\lambda^0_1$ is the polar eigenvalue, and $\tau^{1}$, $\omega^{1}$and $\tau^{2}$ are constants. It turns out that $\tau^1 = \omega^1$. These laws suggest that toroidal modes become singular when $E=E_{{\rm m}}$approach zero thus revealing their true nature, namely that of magnetic shear layers.

By identifying Eqs. (13) with (22), we find that:

$$\alpha_1 = \lambda_1^0 + E^{1/2} (\tau^1 + {\rm i} \omega^1) + E \tau^2.$$ (23)

A comparison between $\beta_1$ and $\tau^1+{\rm i} \omega^1$ permits the identification:

$$\beta_1 = E^{1/2} (\tau^1 + {\rm i}\omega^1).$$ (24)

Therefore, we can cast Eq. (13) into a new form, which separates the adiabatic contribution from the non-adiabatic ones:

$$\lambda_{n,q} = \lambda^0_n + E^{1/2} (q+1) (\tau^1_n + {\rm i} \omega^1_n) + E \tau^2_n. \\$$ (25)

The third term, $\tau^2_n$, is roughly proportional to the square of the frequency and can therefore become noticeable when n is sufficiently large. If the eigenvalues were reduced to the first two terms, then for arbitrary values of E, they would all fall on the lines given by:

$$\frac{{\rm Im}(\lambda)-\omega_n^0}{\omega_n^1} = \frac{{\rm Re}(\lambda)}{\tau_n^1}·\cdot$$ (26)

When $E_{{\rm m}}$ is fixed, the situation is different. The eigenvalue does not appear to converge toward its corresponding polar eigenvalue. At the same time, the mode structure does not become singular, which is consistent with the behaviour of the eigenvalue. The same experiment can be repeated, but with E fixed instead of $E_{{\rm m}}$. The same phenomenon appears. However, the limit eigenvalue for $(E_{{\rm m}}= 0,E=10^{-3})$ is not the same as for $(E_{{\rm m}}= 10^{-3},E=0)$.

3.6 Asymptotic solutions

Having separated the adiabatic and non-adiabatic contributions to the eigenvalue, we will now give a justification of Eq. (25), based on an analytical model. We assume that the diffusivities are of the form $E=K\varepsilon$ and $E_{{\rm m}}=K_{{\rm m}}\varepsilon$, with $\varepsilon$ approaching zero. We find that to zeroth order in $\varepsilon$, eigenmodes can be put in the following form (see Appendix A):

$$b^{0}(r,\hat{\nu}) = b^0_n(r)f(\hat{\nu}),$$

$$v^{0}(r,\hat{\nu}) = v^0_n(r)f(\hat{\nu}),$$ (27)

where (b⁰_n,v⁰_n) are given by Eq. (21), and $\hat{\nu} = \varepsilon^{-1/4} \sin \theta /\sqrt{r}$ remains constant along field lines. We prove in Appendix A that the function f satisfies the following differential equation (see Eq. (A.10)):

$$\lambda_n^1 \mathcal{C}_1 f = - \frac{\lambda_n^0 \hat{\nu}... ...{{\rm d}f}{{\rm d}\hat{\nu}} - \frac{f}{\hat{\nu}^2} \right],$$ (28)

where $\lambda_n^0$ and $\lambda_n^1$ are the two first orders of $\lambda$, and $\mathcal{C}_1$, $\mathcal{C}_2$ and $\mathcal{C}_3$ are constants given by the following formulas:

$$\mathcal{C}_1 = \int_{\eta}^1 r^3 \left( \vert b^0_n\vert^2 + \vert v^0_n\vert^2 \right) ~ {\rm d}r,$$

$$\mathcal{C}_2 = \displaystyle \int_{\eta}^1 r^4 \left( \vert b^0_n\vert^2 + \vert v^0_n\vert^2 \right) ~ {\rm d}r,$$

$$\mathcal{C}_3 = \displaystyle \int_{\eta}^1 K_{{\rm m}}\vert b^0_n\vert^2 + K \vert v^0_n\vert^2 ~ {\rm d}r.$$ (29)

The function f vanishes at the pole since $\vec{b}$ and $\vec{v}$ are toroidal. In order to mimic the confinement of eigenmodes near the poles, we also impose the condition that f approaches zero as $\hat{\nu}$ goes to infinity. We then look for the eigenfunctions f and eigenvalues $\lambda^{1}_n$that satisfy Eq. (28). The variable $\hat{\nu}$ is rescaled to $\rho = \hat{\nu} {\rm e}^{{\rm i}\pi /8} \sqrt[4] {\frac{\omega_n^0 \mathcal{C}_2}{2 \mathcal{C}_3}}$ where $\lambda_n^0 = {\rm i} \omega_n^0$. This leads to the following equation:

$$\frac{{\rm d}^2 f}{{\rm d}\rho^2} + \frac{1}{\rho}\frac{{\r... ...rt{\omega_n^0 \mathcal{C}_2 \mathcal{C}_3/2}} \right) f = 0.$$ (30)

	Figure 7: Comparison between different numerical profiles of the magnetic field and theoretical ones. These profiles are calculated along a meridional cut at a radius $r \simeq 0.5$, and $E=E_{{\rm m}}=4\times 10^{-5}$.
Open with DEXTER

A general solution f of Eq. (30) can be expressed via a linear combination of the following functions:

$$s_1(\rho) = \rho ~ {\rm e}^{-\rho^2/2} M \left( a,2,\rho^2 \right),$$

$$s_2(\rho) = \rho ~ {\rm e}^{-\rho^2/2} U \left( a,2,\rho^2 \right),$$

$$a = \displaystyle 1 + \frac{\lambda_n^{1} {\rm e}^{-{\rm i}\pi/4} \mathcal{C}_1} {\sqrt{8 \omega_n^0 \mathcal{C}_2 \mathcal{C}_3}},$$ (31)

where M and U are confluent hypergeometric functions, solutions to Kummer's equation (see Abramowitz & Stegun 1972, M is known as Kummer's function). When ais not a negative integer, it can be shown that s₁ diverges when $\rho$ goes to infinity and s₂ diverges when $\rho$ approaches 0. Hence, the only way to satisfy our set of boundary conditions is to assume that a is a negative integer, namely a=-q. In that case, the functions M and U both become proportional to the generalised Laguerre polynomials (see Abramowitz & Stegun 1972):

$$M(-q,2,x) = \frac{L_q^{(1)}(x)}{q+1} = \frac{x^{-1}{\rm e}^x}... ...rac{{\rm d}^q}{{\rm d}x^q} \left( x^{1+q}{\rm e}^{-x} \right).$$ (32)

In summary, we obtain the following expressions for toroidal axisymmetric eigenvalues:

 

$$\lambda_{n,q} = \lambda^{0}_n + \varepsilon^{1/2} \lambda^{1}_{n,q} + \mathcal{O}(\varepsilon),$$

$$\lambda^{0}_n = \displaystyle {\rm i}\omega^0_n \simeq {\rm i}\frac{2(2n+1)\pi} {1-\eta^4},$$

$$\lambda^{1}_{n,q} = \displaystyle -2(1+q)(1+{\rm i}) \frac{\sqrt{\omega^0_n \mathcal{C}_2 \mathcal{C}_3}}{\mathcal{C}_1}\cdot$$ (33)

The corresponding modes are given by:

 

$$b = b^0 + \mathcal{O}\left(\varepsilon^{1/2}\right),$$

$$v= v^0 + \mathcal{O}\left(\varepsilon^{1/2}\right),$$

$$b^{0}(r,\theta) = b^0_n(r) \rho ~ {\rm e}^{-\rho^2/2} L_q^{(1)}(\rho^2),$$

$$v^{0}(r,\theta) = v^0_n(r) \rho ~ {\rm e}^{-\rho^2/2} L_q^{(1)}(\rho^2),$$ (34)

in which the variable $\rho$ is given by the following expression:

$$\rho = \varepsilon^{-1/4} {\rm e}^{{\rm i}\pi/8} \left\{ \fra... ...mathcal{C}_3} \right\}^{1/4} \frac{\sin \theta}{\sqrt{r}}\cdot$$ (35)

The form of Eq. (33) corresponds perfectly to that of Eq. (25). A quantitative comparison reveals the accuracy of these formulas (see Table 2). It is particularly interesting to compare the profiles along meridional cuts with these theoretical predictions. Figure 7 shows such a comparison, for different values of q, a comparison which turns out to be excellent. Naturally, the analytical formulas are more accurate when both diffusivities take on small values.

Table 2: Comparison between numerical and theoretical first order eigenvalues ( $\varepsilon ^{1/2} \lambda ^1$). The two sets of numerical values (Cols. 2 and 3) are based on the numerical eigenspectrum given in Fig. 1 and are calculated using a least square method. These values should be very close to each other. The last column is based on Eq. (33).

   
3.7 Boundary layers

To complete our study of axisymmetric toroidal eigenmodes, we also present some results regarding the existence of boundary layers. In our case, we expected the presence of Hartmann layers, which are very thin. Typically their non-dimensional thickness is given by $B_r^{-1} \sqrt{EE_{{\rm m}}}$ in which B_r is non dimensional (e.g. Pothérat et al. 2002). This can be penalising for numerical calculations as a high resolution is needed to properly resolve them.

Two different methods were used to study boundary layers. The first approach consists in coming up with a highly simplified analytical model. In this model, developed in Appendix B for clarity reasons, the permanent magnetic field is constant and vertical, and the fluid domain is enclosed between two parallel planes. This model leads to the conclusion that there should be a Hartmann layer on the lower boundary and nothing on the upper plane. This Hartmann layer only produces a finite discontinuity of the normal gradient of both perturbations when the diffusivities vanish (see Eq. (B.3)). The model also justifies a posteriori the effective boundary conditions $v(\eta)=0$ and b(1)=0 (see Eq. (B.9)) which were used in the calculation of polar eigenvalues (see Sect. 3.4).

The second approach is numerical. By looking at the Chebyshev spectral coefficients, we can see the signature of Hartmann layers. Assuming the layer is described by an exponential variation as proposed by Pothérat et al. (2002), it is possible to come up with a corresponding theoretical prediction, which can then be compared with the actual Chebyshev spectral coefficients of the eigenmodes.

	Figure 8: Comparison between numerical $(E=E_{{\rm m}}=10^{-4})$ and theoretical spectral coefficients which confirms the presence of boundary layers.
Open with DEXTER

Figure 8 shows such a comparison. The values used for B_r are those at the magnetic poles, because of the mode's localization. In other words, the theoretical inner boundary has a thickness of $\eta^3 \sqrt{EE_{{\rm m}}}$ and the outer one a thickness of $\sqrt {EE_{{\rm m}}}$ (see Eq. (1)). The theoretical boundary layers are multiplied by a constant so as to match the numerical ones. From this figure, we clearly deduce the presence of two boundary layers, one on the inside and one on the outside, unlike what was predicted by the analytical model. This is so probably because in the simplified model, the two outer boundary conditions become redundant when E and $E_{{\rm m}}$ approach zero, whereas they do not in the spherical setup.

For larger diffusivities, the agreement between the numerical and theoretical results is not as good. This is not too surprising as the eigenmodes are in general further from the poles. At smaller diffusivities, the radial resolution rapidly becomes insufficient for good comparisons. The theoretical spectral coefficients for the inner boundary diminish very little for the first several hundred Chebyshev components.

A final comment can be made about the numerical approach. By looking at the relative amplitude of the layer on the boundary, it is possible to determine what type of discontinuity it produces when both diffusivities approach zero. In our case, the amplitudes for both the inner and outer layers, obtained by comparing the numerical and theoretical spectral coefficients, were proportional to E in the case where $E=E_{{\rm m}}$, which suggests that there will be a discontinuity on the gradient of the perturbations when both diffusivities approach zero. This is the same behaviour as the single boundary layer in the analytical model.

   
4 Non-axisymmetric modes

We now focus on non-axisymmetric modes $(m \neq 0)$. From the numerical point of view, they are twice as demanding as their axisymmetric counterparts for a given resolution: the $(u_m^{\ell},a_m^{\ell})$ components are coupled to the $(w_m^{\ell},c_m^{\ell})$components, which gives birth to eigenvectors twice as large.

4.1 Eigenvalue spectrum

As was the case for poloidal and toroidal eigenspectra, non-axisymmetric eigenvalues are located along "horizontal'' branches. However, where there used to be one branch in the poloidal or toroidal case, there are now two branches next to each other for each $m \neq 0$ (see Fig. 9). This is perhaps not surprising as non-axisymmetric modes contain poloidal and toroidal components at the same time.

	Figure 9: Several non-axisymmetric eigenvalue spectra. The diffusivities Eand $E_{{\rm m}}$ both take on the value 10-3. These branches overlap with the axisymmetric ones.
Open with DEXTER

As diffusivities become small, the two branches start to differ. When compared with poloidal or toroidal branches, it becomes immediately obvious that a first group of eigenvalues resembles the toroidal branch and a second group matches its poloidal counterpart (see Fig. 10).

Figure 10: Comparison between axisymmetric poloidal, axisymmetric toroidal and non-axisymmetric (m=1) eigenvalues. The diffusivities E and $E_{{\rm m}}$ both take on the values $4\times 10^{-5}$. The lower non-axisymmetric branch corresponds to poloidal-like modes. The upper branch contains toroidal-like modes when $\tau$ is sufficiently negative. The letters "a'', "b'', "c'' and "d'' label non-axisymmetric modes that are represented in Fig. 12.

Open with DEXTER

It is then interesting to look at the structure of the corresponding modes. Not surprisingly, modes along the lower branch have a very similar appearance to that of poloidal modes. The radial nodes are in the same positions and the horizontal nodes look the same. As for the modes on the upper branches, it appears that the least-damped ones are poloidal-like, and the most damped toroidal-like. Figure 11 shows a comparison between axisymmetric poloidal and toroidal modes and non-axisymmetric modes with similar appearances. Figure 12 shows the transition from poloidal properties to toroidal ones, when looking at successive modes on the upper branch.

	Figure 11: Comparison between axisymmetric poloidal and toroidal modes and non-axisymmetric modes with similar properties. The axisymmetric modes are above, and the "corresponding'' non-axisymmetric ones below.
Open with DEXTER

	Figure 12: Transition from poloidal to toroidal characteristics in non-axisymmetric modes. The diffusivities are $E=E_{{\rm m}}=4\times 10^{-5}$ and m is equal to 1. The letters underneath each diagram correspond to labels in Fig. 10.
Open with DEXTER

One way to characterise whether a mode is more poloidal or toroidal is to calculate the following ratio:

$$\chi =\frac{H(u_m^{\ell})+H(v_m^{\ell})+H(a_m^{\ell})+H(b_m^{\ell})}{H(w_m^{\ell})+H(c_m^{\ell})},$$ (36)

where H(x) represents the energy of component x. $\chi$ is the ratio of "poloidal'' energy (both kinetic and magnetic) to "toroidal'' energy. Explicitly, it is given by:

$$\chi =\frac{\displaystyle \sum_{\ell=\ell_{\min}}^{L} \!\! \i... ...l}\vert^2 + \vert c_m^{\ell}\vert^2 \right) r^2 {\rm d}r}\cdot$$ (37)

When this test is applied to the modes of Fig. 12, we confirm that "poloidal'' modes have strong $a_m^{\ell}$, $b_m^{\ell}$, $u_m^{\ell}$ and $v_m^{\ell}$ components whereas "toroidal'' ones have strong $w_m^{\ell}$ and $c_m^{\ell}$ components (see Table 3).

4.2 Asymptotic numerical behaviour

We can take a more detailed look at the behaviour of an entire pair of branches for different values of E and $E_{{\rm m}}$ (keeping $E=E_{{\rm m}}$), to see whether the eigenvalues have poloidal or toroidal characteristics. Non-axisymmetric "poloidal'' eigenvalues appear to obey the following empirical law:

$$\lambda = i\omega^{0} + {\rm i}\omega^1 E^{\gamma} + \tau^1 E^{\delta},$$ (38)

where $\gamma \sim 1.5$ and $\delta \sim 1$. This is different from what was observed in the case of axisymmetric toroidal modes, which show a slower convergence toward an asymptotic adiabatic value. In addition each non-axisymmetric poloidal eigenvalue seems to have its own asymptotic frequency, as may also be the case for axisymmetric poloidal modes. A plot of the behaviour of different eigenvalues with respect to diffusivities is given in Fig. 13. $\omega_i$ was estimated using a least-square method for the poloidal modes and differs for each one of them.

Non-axisymmetric "toroidal'' eigenvalues behave roughly like their axisymmetric counterparts. All the values on a single branch seem to converge toward one frequency, the corresponding polar eigenvalue (see Eq. (20)). This value is used in Fig. 13 as the asymptotic frequency of the toroidal modes. The rate at which eigenvalues converge is about the same as in the axisymmetric case, namely $\lambda_i - \lambda \propto E^{1/2}$.

It is also interesting to note that the shape of the toroidal eigenbranch remains more or less the same in the axisymmetric and non-axisymmetric cases (see Fig. 14, bottom). This indicates that the quantization of these modes in the horizontal direction is very similar to the one of the axisymmetric modes. However, the spacing between each mode is not the same (see Fig. 10).

	Figure 13: Behaviour of $\tau$ and $\Delta \omega$ ( $=\omega _{i} - \omega$) for four non-axisymmetric modes. The toroidal modes start off with a poloidal structure, and undergo a metamorphosis into a toroidal appearance. The poloidal non-axisymmetric modes approach their inviscid values much faster.
Open with DEXTER

Table 3: Energy ratios as defined in Eq. (36) for modes shown in Figs. 10 and 12.

Figure 14: Behaviour of non-axisymmetric eigenbranches for different values of the diffusivities. E(= $E_{{\rm m}})$ takes on the values 10-3, $4\times 10^{-4}$, $2\times 10^{-4}$, 10-4 and $4\times 10^{-5}$. As diffusivities diminish, both eigenbranches get closer to the imaginary axis since $\vert\tau \vert$ decreases. The first figure corresponds to "poloidal'' eigenvalues and the second one to "toroidal'' ones. The classification between the two types of values was done by hand. The dotted line represents the line given by Eq. (26) for n=1.

Open with DEXTER

	Figure 15: Behaviour of non-axisymmetric modes for different values of the diffusivities. The "poloidal'' modes seem to converge toward a non-singular structure. The "toroidal'' ones, on the other hand, do not stop shrinking. The behaviour of their mean position (not shown), however, is much more irregular.
Open with DEXTER

By looking at the behaviour of the structure of different modes, it is possible to observe different characteristics (see Fig. 15). Poloidal-like modes seem to converge toward a non-degenerate structure. In other words, their widths ( $\sigma _{\theta }$) do not approach zero, nor do their positions ( $\overline {\theta }$). Their nodes become more distinct at low diffusivities.

Toroidal-like modes, on the other hand, exhibit more complex behaviour as the toroidal character appears only below some given diffusivity. We understand this behaviour as the consequence of the E^1/4-scale of toroidal modes. It is only when $E^{1/4} \ll 1$ that the toroidal character may appear but this condition is not easy to satisfy with non-axisymmetric modes for which the lowest diffusivities reached yet are $\sim$10^-5 because of memory storage requirements.

5 Discussion

The main result of this study is that there are two basic behaviours for shear Alfvén modes, both in the axisymmetric and the non-axisymmetric cases. The first one is the poloidal behaviour, for which modes remain regular as diffusivities approach zero, except for internal shear layers located on resonant field lines. The second one is the toroidal behaviour: these modes become singular as diffusivities approach zero and each branch appears to have only one eigenfrequency in the adiabatic limit.

The toroidal behaviour matches the description of the resistive Alfvén spectrum given by Kerner et al. (1986) in the context of tokamak-like configurations. The eigenvalues lie on well-defined curves in the complex plane which join the endpoints of the ideal Alfvén continua, as is suggested by Eq. (26). The slope of the eigenbranches is 45° and the damping rates of the eigenmodes are proportional to $\varepsilon^{1/2}$, which is in complete agreement with a quadratic variation of the permanent magnetic field $\theta$-profile near the magnetic axis.

Poloidal modes on the other hand, look like true eigenmodes of the ideal MHD operator with a finite number of resonances in the Alfvén continuum (resonating field lines, see Paper I). Such features are similar to those of "quasimodes" (see Goossens et al. 2002; Poedts & Kerner 1991) which combine a regular and singular behaviour; however, unlike quasimodes, the damping rate of our modes always seem to vanish in the ideal limit (e.g. Figs. 13 or 5 of Paper I).

Presently, we understand the origin of the difference between the toroidal and poloidal behaviour as the presence or absence of a coupling between magnetic field lines by total pressure. Indeed, neglecting diffusion terms, magnetic and velocity perturbations verify:

 

$$\frac{\partial \vec{v}}{\partial t} = -\vec{\nabla}p_* +(\vec{B} \. \vec{\nabla}) \vec{b} +(\vec{b} \. \vec{\nabla}) \vec{B}$$

$$\frac{\partial \vec{b}}{\partial t} = (\vec{B} \. \vec{\nabla}) \vec{v} -(\vec{v} \. \vec{\nabla}) \vec{B}$$ (39)

where $p_* = p + \vec{B} \. \vec{b}$ is the perturbed total pressure. These equations show that the time evolution of a perturbation depends on the global form of the perturbation only through the pressure gradient term. This is clear if one considers a local perturbation, namely an Alfvén wave packet. Terms like $(\vec{B} \. \vec{\nabla}) \vec{b}$ or $(\vec{B} \. \vec{\nabla}) \vec{v}$ modulate the evolution of the perturbation according to its derivative along the field lines while terms like $(\vec{b} \. \vec{\nabla}) \vec{B}$ and $(\vec{v} \. \vec{\nabla}) \vec{B}$introduce a dependence with respect to the local value of the perturbation. Last, the term $\vec{\nabla}p_*$ makes the evolution of the perturbation sensitive to the pressure perturbation which cannot be restricted to a single field line. In pure toroidal modes (axisymmetric ones) this term is zero and therefore no coupling is possible in the ideal MHD case. In non-axisymmetric modes of toroidal type it is weak and strengthens as one shifts to modes with a strong poloidal character. Modes studied by Kerner et al. (1986) have a weak, $\mathcal{O}(E_{{\rm m}}^{1/2})$, pressure term and are thus similar to our toroidal modes.

6 Conclusions

In this paper we completed the study of shear Alfvén modes of a spherical shell started in Paper I. The global picture that emerges is that both non-axisymmetric and axisymmetric modes can be divided into two subclasses: the toroidal and poloidal ones. In the limit of small diffusivities, relevant for astrophysical objects, only poloidal modes survive since toroidal ones become singular as magnetic shear layers.

However, toroidal modes may play some important role in the excitation process, especially if rotation is present. Indeed, it is well-known (e.g. Rieutord 1991) that the Coriolis force induces a coupling between toroidal and poloidal components, even in the axisymmetric case. This coupling may not strongly influence the frequency of the modes (except for some rotational splitting) but may represent a channel by which modes lose energy through dissipative shear layers; in such a case this mechanism may decide whether a mode is or is not excited by, say, a $\kappa$-mechanism.

For astrophysical applications where the compressibility is an important aspect, pure Alfvén modes coexist with magneto-acoustic modes. The incorporation of these modes will undoubtedly lead to a much more complicated spectrum of eigenvalues as well as a whole new range of physical phenomena. In particular, we can expect new resonances, namely those coming from the interaction of eigenmodes and the slow acoustic, or cusp, continuum (see Sakurai et al. 1991). Taking into account compressibility will however not affect the splitting between pure poloidal and pure toroidal axisymmetric modes; because of the divergence-less nature of purely toroidal modes, we can expect that these modes will not be affected by compressibility and remain unchanged. All the complications will concentrate on the poloidal set of eigenmodes. In the case of non-axisymmetric modes where the poloidal and toroidal characters are mixed, we may expect that only modes with a strong poloidal shape will suffer important changes. However, the importance of these changes, which incorporate the resonances with the cusp continuum, as well as its influence on the shape of the spectrum, can hardly be anticipated and will be investigated in subsequent work.

Acknowledgements

The authors wish to thank the referee for valuable suggestions and comments, which helped to improve this article.
Some of the numerical calculations were carried out on the NEC SX5 of the "Institut du Développement et des Ressources en Informatique Scientifique'' (IDRIS), which is gratefully acknowledged.

References

• Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (New York: Dover) (In the text)
• Alfvén, H. 1945, MNRAS, 105, 3 [NASA ADS] (In the text)
• Balmforth, N. J., Cunha, M. S., Dolez, N., Gough, D. O., & Vauclair, S. 2001, MNRAS, 323, 362 [NASA ADS] [CrossRef] (In the text)
• Bigot, L., & Dziembowski, W. A. 2002, A&A, 391, 235 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Biront, D., Goossens, M., Cousens, A., & Mestel, L. 1982, MNRAS, 201, 619 [NASA ADS]
• Campbell, C. G., & Papaloizou, J. C. B. 1986, MNRAS, 220, 577 [NASA ADS] (In the text)
• Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (New York: Dover), 155 (In the text)
• Dziembowski, W. A., & Goode, P. R. 1996, ApJ, 458, 338 [NASA ADS] [CrossRef] (In the text)
• Goossens, M., Andries, J., & Aschwanden, M. J. 2002, A&A, 394, L39 [EDP Sciences] [NASA ADS] [CrossRef]
• Kato, Y., & Watanabe, T. 1958, JGR, 63, 741 [NASA ADS] (In the text)
• Kerner, W., Lerbinger, K., & Riedel, K. 1986, Phys. Fluids, 29, 2975 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Kurtz, D. W. 1978, Informational Bulletin on Variable Stars, 1436, 1 [NASA ADS]
• Lou, Y. 1995, MNRAS, 275, L11 [NASA ADS] (In the text)
• Poedts, S., & Kerner, W. 1991, Phys. Rev. Lett., 66, 2871 [NASA ADS] [CrossRef]
• Pothérat, A., Sommeria, J., & Moreau, R. 2002, Phys. Fluids, 14, 403 [NASA ADS] [CrossRef] (In the text)
• Rieutord, M. 1991, Geophys. Astrophys. Fluid Dyn., 59, 185 [NASA ADS] (In the text)
• Rincon, F., & Rieutord, M. 2003, A&A, 398, 663 [EDP Sciences] [NASA ADS] [CrossRef]
• Roberts, P. H., & Soward, A. M. 1983, MNRAS, 205, 1171 [NASA ADS]
• Sakurai, T., Goossens, M., & Hollweg, J. V. 1991, Sol. Phys., 133, 227 [NASA ADS] (In the text)
• Stevenson, D. J. 1983, Rep. Progr. Phys., 46, 555 [NASA ADS] (In the text)

 

  
Online Material

   
Appendix A: $\theta$-dependence of asymptotic axisymmetric toroidal modes

We begin with the toroidal Eqs. (12), and express them explicitly in "semi-dipolar'' coordinates $(r,\nu,\varphi)$. The coordinate $\nu$ is equal to $\sin \theta/\!\sqrt{r}$, and remains constant along field lines. Choosing $\nu$ as a coordinate instead of $\theta$simplifies the expression of terms like $(\vec{B} \cdot \nabla)\vec{b}$. We obtain the following system:

$$\lambda b = \displaystyle \sqrt{1- r \nu^2} \left[ \frac{1}{r^3} \frac{\partial v}{\partial r} - \frac{3v}{2 r^4} \right] + E_{{\rm m}}\Delta' b,$$

$$\lambda v= \displaystyle \sqrt{1- r \nu^2} \left[ \frac{1}{r^3} \frac{\partial b}{\partial r} + \frac{3b}{2 r^4} \right] + E \Delta' v,$$ (A.1)

where $\Delta'$ takes on the following expression:

$$\Delta' b = \displaystyle \frac{\partial^2 b}{\partial r^2} - \fr... ... + \frac{1}{\nu} \frac{\partial b}{\partial \nu} - \frac{b}{\nu^2} \right]\cdot$$

In order to analyse the asymptotic limit of small diffusivities, we suppose E and $E_{{\rm m}}$ are of the form $K\varepsilon$ and $K_{{\rm m}}\varepsilon$ respectively, where K and $K_{{\rm m}}$ are two constants, and $\varepsilon$ a free parameter that approaches zero. We then rescale $\nu$ to $\hat{\nu} = \varepsilon^{-1/4} \nu$, in accordance with our numerical results, and neglect terms smaller than $\mathcal{O}(\varepsilon^{1/2})$. This leads to the following system:

 

$$\lambda b = \left( 1- \frac{1}{2} \varepsilon^{1/2} r \hat{\nu}^2... ...}{2 r^4} \right] \!\! + \! \frac{\varepsilon^{1/2}K_{{\rm m}}}{r^3} \Theta [b],$$

$$\lambda v= \left( 1- \frac{1}{2} \varepsilon^{1/2} r \hat{\nu}^2 ... ...+ \frac{3b}{2 r^4} \right] \!\! + \! \frac{\varepsilon^{1/2}K}{r^3} \Theta [v],$$ (A.2)

where the differential operator $\Theta$ is given by:

$$\Theta [b] = \frac{\partial^2 b}{\partial \hat{\nu}^2} + \f... ...\partial b}{\partial \hat{\nu}} - \frac{b}{\hat{\nu}^2} \cdot$$

Following the form of Eq. (25), we develop the solutions in the following manner:

$$\lambda = \lambda^0 + \varepsilon^{1/2} \lambda^{1},$$

$$b = b^{0} + \varepsilon^{1/2}b^{1},$$

$$v= v^{0} + \varepsilon^{1/2}v^{1}.$$ (A.3)

At zeroth order, we obtain the equations that give the polar eigenspectrum (see Eq. (16) for $\cos \theta = 1$). We conclude that the zeroth order solutions to Eq. (A.2) are of the form:

$$b^{0}(r,\hat{\nu}) = b^0_n(r)f(\hat{\nu}),$$

$$v^{0}(r,\hat{\nu}) = v^0_n(r)f(\hat{\nu}),$$

$$\lambda^{0} = \lambda^0_n,$$ (A.4)

where (b⁰_n,v⁰_n) is given by Eq. (21), and $\lambda ^0_n$ is the corresponding eigenvalue. $\lambda ^0_n$ gives us the position of the different eigenbranches (as already seen in Sect. 3.4). Finding a quantization on the function f would enable us to explain the different eigenvalues on a given branch. This requires looking at the first order of Eq. (A.2) (in which we substitute results from the zeroth order):

 

$$\lambda^0_n b^1 - \frac{1}{r^3} \frac{\partial v^1}{\partial r} +... ...{\lambda^0_n r\hat{\nu}^2 b^0_n f}{2} +\frac{b^0_n K_{{\rm m}}\Theta [f]}{r^3},$$

$$\lambda^0_n v^1 - \frac{1}{r^3} \frac{\partial b^1}{\partial r} -... ... f - \frac{\lambda^0_n r\hat{\nu}^2 v^0_n f}{2} +\frac{v^0_n K\Theta [f]}{r^3},$$ (A.5)

$$b^1(r=1,\hat{\nu}) = 0,$$

$$v^1(r=\eta,\hat{\nu}) = 0.$$

(We now use the notation $\lambda_n^1$ for $\lambda^1$.) This equation is of the form $\mathcal{L}_0 \Psi_1 = \mathcal{L}_1 \Psi_0$, where $\mathcal{L}_0$, $\mathcal{L}_1$ are respectively zeroth and first order linear operators, and $\Psi_0$, $\Psi_1$ are zeroth and first order solutions. It is solvable only if its right-hand side is orthogonal to the kernel of the adjoint of the operator $\mathcal{L}_0$, with respect to the following dot product:

$$\left< \left[ \begin{array}{l} u \\ [3mm] w \end{array} \righ... ...^*(r)x(r) + w^*(r) y(r)}{\Vert \vec{B}(r,0) \Vert} ~ {\rm d}r,$$ (A.6)

where u^* is the complex conjugate of u, and $\Vert \vec{B}(r,0) \Vert = 1/r^3$. The adjoint operator $\mathcal{L}_0^\dag$ is deduced from $\mathcal{L}_0$ by using integration by parts:

$$\mathcal{L}_0^{\dag } \left[ \begin{array}{l} u \\ [3mm] w \e... ...ac{\partial u}{\partial r}+\frac{3u}{2r^4} \end{array}\right].$$ (A.7)

The adjoint boundary conditions are chosen so as to cancel out terms like $w^*(1)b(1) - w^*(\eta)b(\eta)$, which come from the integration by parts:

$$u(r=1) = 0, \qquad w(r=\eta) = 0.$$ (A.8)

The zeroth order eigenvalue $\lambda_n^0$ is purely imaginary. As a result the operator $\mathcal{L}_0$ is antisymmetric: $\mathcal{L}_0^\dag = -\mathcal{L}_0$. The boundary conditions being the same for both operators, the kernel of $\mathcal{L}_0^\dag$ is the same as that of $\mathcal{L}_0$:

$$u = b_n^0, \qquad w = v_n^0.$$ (A.9)

By calculating the dot product between (b_n⁰,v_n⁰) and the terms of Eq. (A.5) we obtain the following equation:

$$-\lambda^1_n \mathcal{C}_1 f -\frac{\lambda_n^0 \hat{\nu}^2 \mathcal{C}_2}{2}f +\mathcal{C}_3 \Theta [f] = 0$$ (A.10)

where:

 

$$\mathcal{C}_1 \left< \left[ \begin{array}{l} u \\ [3mm] w \end{array} \right], \left[ \begin{array}{l} b^0_n \\ [3mm] v^0_n \end{array} \right] \right>$$ =

$$\int_{\eta}^1 \! r^3 \left( \vert b^0_n\vert^2 + \vert v^0_n\vert^2 \right) ~ {\rm d}r,$$ =

$$\mathcal{C}_2 \left< \left[ \begin{array}{l} u \\ [3mm] w \end{array} \right], r \left[ \begin{array}{l} b^0_n \\ [3mm] v^0_n \end{array} \right] \right>$$ =

$$\int_{\eta}^1 \! r^4 \left( \vert b^0_n\vert^2 + \vert v^0_n\vert^2 \right) ~ {\rm d}r,$$ =

$$\mathcal{C}_3 \left< \left[ \begin{array}{l} u \\ [3mm] w \end{array} \right],... ...\begin{array}{l} K_{{\rm m}}b^0_n \\ [3mm] K v^0_n \end{array} \right] \right>$$ =

$$\int_{\eta}^1 \! K_{{\rm m}}\vert b^0_n\vert^2 + K \vert v^0_n\vert^2 ~ {\rm d}r.$$ = (A.11)

Equation (A.10) gives the quantization along each eigenbranch.

   
Appendix B: Simplified analytical model

The model presented here justifies the use of effective boundary conditions throughout the paper for the adiabatic calculations. In this model, we use Cartesian coordinates (x,y,z), and the fluid domain is enclosed between two infinite horizontal planes located at z=0 and z=1. The permanent magnetic field is constant and in the vertical direction $(\vec{B}=B\vec{e}_z)$. The perturbations take on the form $\vec{b} = b(z)~\vec{e}_y$ and $\vec{v} = v(z)\vec{e}_y$, in order to mimic the toroidal direction $\vec{e}_{\varphi}$. Equation (3) becomes:

 

$$\lambda b = \frac {{\rm d}v} {{\rm d}z} + E_{{\rm m}}\frac {{\rm d}^2 b}{{\rm d}z^2},$$

$$\lambda v= \frac {{\rm d}b} {{\rm d}z} + E \frac {{\rm d}^2 v}{{\rm d}z^2}\cdot$$ (B.1)

The corresponding boundary conditions are:

$$\frac{{\rm d}v}{{\rm d}z}(0) = 0, \quad E_{{\rm m}}\frac{{\rm d}b}{{\rm d}z}(0) = -v(0),$$

$$\frac{{\rm d}v}{{\rm d}z}(1) = 0, \quad b(1)=0.$$ (B.2)

Since the two conditions on the derivative of v are stress free conditions, they lack physical meaning when both diffusivities vanish. Therefore, we can expect them to disappear in that situation. The solutions to Eq. (B.1) are:

 

$$k_1^{-1}k_2 \sin k_2 \left({\rm e}^{{\rm i}k_1z}+{\rm e}^{{\rm i}k_1(2-z)} \right)$$                                    v =

$$-{\rm i}\left( 1-{\rm e}^{2{\rm i}k_1} \right) \cos (k_2(z-1)),$$

$$k_2 \sin k_2 \left( \zeta^- k_1^{-1} + {\rm i} \xi \right) \left( {\rm e}^{{\rm i}k_1z}-{\rm e}^{{\rm i}k_1(2-z)} \right)$$ b =

$$+ \left( 1 - {\rm e}^{2{\rm i}k_1} \right) \left( \zeta^- + {\rm i} \xi k_2 \right)\sin (k_2(z-1)),$$ (B.3)

where:

$$\zeta^+ = \sqrt {1 + \left( \sqrt{E} + \sqrt{E_{{\rm m}}} \right)^2 \lambda},$$

$$\zeta^- = \sqrt {1 + \left( \sqrt{E} - \sqrt{E_{{\rm m}}} \right)^2 \lambda},$$

$$\xi = \sqrt{EE_{{\rm m}}} - E,$$

$$k_1 = \displaystyle \frac{{\rm i}(\zeta^- + \zeta^+)}{2\sqrt{EE_{{\rm m}}}} \sim \frac{{\rm i}}{\sqrt{EE_{{\rm m}}}},$$

$$k_2 = \displaystyle \frac{{\rm i}(\zeta^- - \zeta^+)}{2\sqrt{EE_{{\rm m}}}} \sim -{\rm i} \lambda.$$ (B.4)

The quantities k₁, k₂, -k₁,and -k₂ are solutions of the dispersion relation of Alfvén waves in a uniform field (Chandrasekhar 1961):

$$-(\lambda + E k^2)(\lambda + E_{{\rm m}}k^2) = k^2$$ (B.5)

where the non-dimensional Alfvén velocity is 1 and the angle between the wave vector $\vec{k}$ and the permanent magnetic field is 0.

The eigenvalue $\lambda$ must satisfy the following relation:

$$\left\{ -\zeta^- \lambda + {\rm i}[(\zeta^-)^2 + \xi \lambda]k_1 \right\} \sin k_1 \cos k_2$$ (B.6)

$$\qquad = \left\{ -\zeta^- \lambda + {\rm i}[(\zeta^-)^2 + \xi \lambda]k_2 \right\} \sin k_2 \cos k_1.$$ (B.7)

If the ratio between E and $E_{{\rm m}}$ remains constant, we obtain the following formula for $\lambda$:

$$\lambda = {\rm i} \Omega_n - \left( \frac{E+E_{{\rm m}}}{2} \right) \Omega_n^2 + \mathcal{O}(E^2),$$ (B.8)

where $\Omega_n = \left(n + \frac{1}{2} \right) \pi$.

When E and $E_{{\rm m}}$ approach zero, the solutions given in Eq. (B.3) obey the following effective boundary conditions:

$$v(0) = 0 \qquad b(1) = 0,$$ (B.9)

and their gradients admit a discontinuity on the lower boundary. It is interesting to note that the stress-free condition has disappeared on the lower boundary, as was expected. On the upper boundary, it simply became redundant with the condition b(1)=0.

When comparing the eigenvalues given in Eq. (B.8) to the ones given in Eq. (33), we notice that the term $\mathcal{O}(\varepsilon^{1/2})$ is lacking in the first formula. This is because in the simplified model, eigenmodes are invariant in the $\vec{e}_x$ direction, thus removing any quantization in that direction.