Spherical magnetic instabilities in the Earth’s core and equatorial symmetries

Physics of the Earth and Planetary Interiors 153 (2005) 83–100

Spherical magnetic instabilities in the Earth’s core and equatorial symmetries C.G. Phillipsa,∗ , D.J. Iversb b

a Mathematics Learning Centre, University of Sydney, Australia School of Mathematics & Statistics, University of Sydney, Australia

Received 22 November 2004; received in revised form 3 April 2005; accepted 4 May 2005

Abstract The Boussinesq rotating magnetohydrodynamic equations in a conducting fluid sphere, linearised about a steady basic state, are discretised using scalar and vector spherical harmonics. Vector spherical harmonic representations are used for the basic state vector fields as well as the interaction terms in the equations. The solenoidal property of the perturbation magnetic and velocity fields are imposed on the fields and equations using relations between vector spherical harmonics and toroidal–poloidal representations. The method has several distinct advantages over previous discretisation methods: it produces compact angular spectral equations which require only five core subroutines to generate all interactions and product terms, and only six core subroutines to radially discretise the equations using finite-differences; the core subroutines can be extensively tested using only simple models; and its generality permits very flexible code. We apply our approach to the magnetic instabilities of a steady axisymmetric azimuthal magnetic field. To isolate magnetic instabilities and suppress other instabilities the magnetostrophic approximation is made. The perturbation velocity and magnetic field satisfy the inviscid boundary and insulating exterior conditions, respectively. Results are presented for the simplest equatorially symmetric (quadrupole) and equatorially antisymmetric (dipole) magnetic field models which vanish on a single spherical surface. For each model the perturbation fields decouple into two subproblems, in which the perturbation magnetic field has the same or opposite symmetry to the basic state field. A strong preference is found for the equatorially symmetric perturbation magnetic field in the equatorially symmetric model. However, there is only a marginal preference for the equatorial antisymmetric perturbation magnetic field in the equatorially antisymmetric model. We also consider a class of equatorially unsymmetric basic state magnetic fields which are convex combinations of the first two models. Although the critical Elsasser number does not change appreciably, the most unstable mode changes totally from equatorially symmetric to equatorially antisymmetric. Finally, we consider the effect on the equatorially symmetric model of varying the spherical surface on which the field vanishes. The most unstable mode is always equatorially symmetric. When the dimensionless magnetic-zero radius is decreased below about 0.7 the most unstable mode switches from eastward to westward propagation. The magnetic field perturbation in the westward (eastward) case is strongly confined to the region outside (inside) the magnetic-zero radius. © 2005 Published by Elsevier B.V. Keywords: Spherical magnetic instabilities; Earth’s core; Equatorial symmetries

1. Introduction

∗ Corresponding author. Tel.: +61 2 9351 4061; fax: +61 2 9351 5797. E-mail address: [email protected] (C.G. Phillips).

0031-9201/$ – see front matter © 2005 Published by Elsevier B.V. doi:10.1016/j.pepi.2005.05.008

Present day fully-dynamical geodynamo codes produce solutions that apparently simulate aspects of the terrestrial dynamo. However, the solutions are typically expensive and of great complexity, which makes physical

84

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

interpretation and understanding difficult. Studies based on the linearised equations are simpler to implement and quicker to complete, yet they still offer valuable insight into the mechanisms active within planetary and solar dynamos and complement non-linear dynamicallyconsistent dynamo calculations. In particular, linearised magnetohydrodynamic studies can allow a more detailed examination of the interactions of a subset of dynamo mechanisms. Further, linearised computations can reach regions of parameter space inaccessible to fully dynamical dynamo codes and more easily implement important physical effects such as small viscous and thermal diffusion, anisotropic diffusion, the effects of non-spherical boundaries, non-radial gravitational fields, topographic core-mantle coupling and the anelastic approximation. The linearised equations can also help in understanding the numerical solution of the non-linear equations. Vector spherical harmonic methods are used, in conjunction with toroidal–poloidal representations of the perturbation fields, to produce compact spectral equations. This type of spectral equation, which was first derived for the magnetic induction equation by James (1974), is quite distinct from analogues of the Bullard and Gellman (1954) spectral induction equations. The vector fields of the basic state are represented by linear combinations of vector spherical harmonics, as are the perturbation vector fields in the interactions between the basic state and the perturbation state. The vector spherical harmonic coefficients of the perturbation vector fields are subsequently related to the spherical harmonic coefficients of their toroidal–poloidal potentials (see (A.5) and (A.6)) and the relevant spectral equations. The coupling integrals, which arise from the spectral interactions of products in the equations, are all evaluated in terms of 3j-, 6j- and 9j-symbols. The method has significant advantages over spectral equations of Bullard and Gellman (1954) form. The vector spherical harmonic equations are more compact and therefore simpler and less error prone to produce. The vector spectral equations and, in particular, the interaction products are coded directly into the computer program (Ivers and Phillips, 2003). The process bypasses the difficult and possibly prohibitively complicated derivations of the spectral interactions and their subsequent coding. The spectral interaction expansions are performed by the numerical code itself. Redundancy is easily eliminated with a single vector spectral expression used for all terms of the same form. For example, all cross product terms, such as the induction, Coriolis or Lorenz forces, have the same vector spectral form and use the same numerical subroutine. The resulting numerical code is extremely robust in the sense that simple

models can test most of the code used by more complicated models. Radial discretisation by finite-differences is straightforward. The vector spectral analysis and the techniques of numerically coding the vector spectral equations extends to more comprehensive linear magnetohydrodynamic stability models, which offers another significant improvement over previous spherical methods. Tensor spherical harmonics (James, 1976) may be needed in such extensions, which include anisotropic diffusion, non-radial gravity and the anelastic approximation (Phillips and Ivers, 2000; Ivers and Phillips, 2005). The method is applied to the numerical solution of the linearised magnetohydrodynamic equations of an electrically-conducting fluid sphere. We investigate the effects that discrete equatorial symmetries of the basic state have on the magnetic instabilities of the steady magnetic field models, and in particular, how these symmetries may enhance or inhibit interaction models relevant to the Earth’s core. We consider three toroidal magnetic field models, which vanish on the boundary of the conducting fluid. We also consider a fourth steady magnetic field model to examine the influence that twodimensional zeros of the basic state magnetic field have on magnetic instabilities. The preference for equatorially symmetric or antisymmetric perturbations and their eastward or westward azimuthal propagation are determined. The dependence of the instabilities on azimuthal wavenumber is also examined. The non-linear equations in a frame rotating with uniform angular velocity
, which govern the velocity v, the pressure p, the magnetic induction B and the temperature
of an electrically-conducting Boussinesq fluid, are ρ∂t v = −∇P − ρω × v − ρ2
× v + J × B − ρα

g + ρν∇ 2 v

(1.1)

∂t B = η∇ 2 B + ∇ × (v × B)

(1.2)

∂t
= −v · ∇
+ κ∇ 2
+ ∇·v = 0

∇·B = 0,

Q ρcp

(1.3) (1.4)

where P = p + (1/2)ρv2 is the modified pressure, ω = ∇ × v is the vorticity, J = ∇ × B/µ0 is the electrical current, ν, η and κ are the viscous, magnetic and thermal diffusivities, ρ is the fluid density, α
is the thermal expansivity, Q is the rate of heat production per unit volume and cp is the specific heat at constant pressure. The gravitational field may be asymmetric. The diffusivities and the specific heat are constants.

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

We consider linear instabilities of a steady axisymmetric basic state (v0 , B0 ,
0 , P0 ). Ideally, we make the crucial assumption that the fields (v0 , B0 ,
0 , P0 ) satisfy the governing Eqs. (1.1)–(1.4), at least approximately. The fields v, B,
and P can be resolved into basic state and perturbation fields, v = v0 + v , B = B0 + B ,
=
0 +
 , P = P0 + P  . Subtracting the equations governing the basic state and neglecting terms which are second order in the perturbation fields, yields the linearised rotating MHD equations, ρ∂t v = −∇P  − ρω0 × v − ρω × v0 − ρ2
× v + J0 × B + J × B0 − ρα

 g + ρν∇ 2 v (1.5) ∂t B = η∇ 2 B + ∇ × (v0 × B ) + ∇ × (v × B0 ) (1.6) ∂t
 = −v0 · ∇
 − v · ∇
0 + κ∇ 2
 + ∇·v = 0

∇·B = 0.

Q ρcp

(1.7) (1.8)

The equations must be supplemented by boundary conditions. In stability problems, Eqs. (1.5)–(1.8) and the associated boundary conditions are homogeneous, in particular Q = 0, and the perturbation fields, B , v ,
 ∝ eγt . The problem is then a generalised eigen- or critical-value problem. Primes on perturbation fields are suppressed below. Linear stability problems can be categorized into important special cases: (I) thermal convection, in which v0 = 0,
0 is the prescribed conduction temperature and B0 is omitted; (II) kinematic dynamos, in which v0 is prescribed and B0 = 0,
0 are omitted; (III) magnetoconvection, in which v0 = 0, and B0 ,
0 are prescribed; and (IV) magnetic instabilities, in which v0 = 0, B0 is prescribed and
0 is omitted. The discretisation of the linearised Eqs. (1.5)–(1.8) in angle is described in Section 2 and Appendices A and B. For both the basic and perturbation states, all vector fields are expanded in vector spherical harmonics and all scalar fields in spherical harmonics, as in Appendix A. For the perturbation solenoidal vector fields, but not the basic state fields, toroidal–poloidal representations are also used. Altogether there are five independent perturbation scalar fields: the temperature, and the toroidal and poloidal potentials of the magnetic field and the velocity. This approach produces compact hybrid angular spectral equations (Ivers and Phillips, 2005), which are less error prone to code, than full spectral expansion of product toroidal–poloidal interactions. Moreover, only

85

three subroutines are required to evaluate the three angular coupling integrals of all 10 products, in terms of 3j-, 6j- and 9j-symbols. Further, the radial discretisation, which is outlined in Appendix B, is greatly simplified, since only six distinct combinations of perturbation radial functions and their derivatives occur. The present code uses second-order finite-differences on a uniform radial grid. The discretisation produces generalised eigen- and critical-value problems with very large banded complex non-hermitian matrices. The solution of these problems using inverse iteration methods, the implicitly restarted Arnoldi method and the Newton– Raphson method is described in Section 3. We apply our approach to magnetic instabilities in the Earth’s core in Section 4. We consider a rapidly-rotating electrically-conducting fluid sphere of radius a with an electrically-insulating source-free exterior. A spherical core allows us to isolate the magnetic interactions of interest and exclude any effects of an inner core. Such a geometry is also of direct interest in the early Earth and some planets, and allows us to compare our results directly with those of Fearn and Weiglhofer (1991b). The present work extends the known solutions by exploring basic state fields with different equatorial (planar) symmetries. We consider four basic state magnetic field models of the form B0 = B(r, θ)1φ ,

(1.9)

where (r, θ, φ) are spherical polar coordinates and 1φ is the unit vector in the φ-direction. The first three magnetic field models vanish on the fluid boundary, B(a, θ) = 0; the fourth model vanishes on a sphere of radius r0 , B(r0 , θ) = 0. The basic state velocity is zero in the frame of the rotating sphere. The modified Rayleigh number is zero and the temperature is omitted. The basic state magnetic field (1.9) with v0 = 0 does not in general satisfy the governing equations, in particular the magnetic induction equation. This though is not of serious concern as we consider instability on the τs timescale, which is considerably shorter than the magnetic diffusion timescale τs over which B0 evolves. Fearn and Weiglhofer (1991a,b, 1992) considered magnetic instabilities in a sphere for a selection of basic state magnetic fields that are equatorially symmetric and equatorially antisymmetric. They found instabilities in both cases though the most unstable modes of those reported were for azimuthal wavenumber m = 2 and they progressed westward. In a subsequent study Zhang and Fearn (1994) considered magnetic instabilities in a spherical shell for several basic state magnetic fields. This later study reported more unstable modes for the m = 1 azimuthal wavenumber with solutions that

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

86

progress eastward. The substantial difference in certain features of the instability motivated the later work. The comparison of results is made more difficult by the difference of an inner core. Our first two models constitute the simplest equatorial symmetric and equatorially antisymmetric basic state magnetic fields in a sphere which are zero on the conducting boundary. These models have the property that the resulting linear stability problem decouples into two sets of field symmetries for the perturbation magnetic and perturbation velocity field for each model. The resulting linear stability problems are solved numerically. In particular, we obtain values for and critical values c of the Elsasser number, =

B02 , 2 ρηµ0

where B0 is the maximum of |B0 | in the conducting sphere V, is the rotation rate of the sphere. The values of c allow us to establish the magnitude of the basic state magnetic field necessary to produce the magnetic instabilities within planetary cores. The azimuthal drift of the instabilities is given by the angular frequency ω on the slow magnetohydrodynamic timescale, τs =

2µ0 ρa2 B02

(1.10)

where a is the natural length scale. Solutions are found for each of the decoupled symmetric sub-problems for each basic state magnetic field model. The third basic state magnetic field model is a convex (linear) combination of the first two models. The fourth model is a variation of the first equatorial symmetric model, which vanishes on a single spherical surface of specified radius. Concluding remarks are made in Section 5. 2. Spectral representation and discrete symmetries of the fields 2.1. Field representation and numerical method Discrete symmetries may be present in a physical problem and can be easily exploited by an effective mathematical formulation of the problem and its numerical implementation to reduce the problem size by decoupling it into simpler sub-problems. A poor mathematical formulation may not simply manifest or even preserve such symmetries and hence not permit independent study of the smaller decoupled sub-problems. This necessitates the treatment of a larger problem than necessary and can lead to numerical and programming difficulties. We use toroidal–poloidal fields together with scalar and vector

spherical harmonics to represent the perturbation vector fields B, v, J, ω, scalar fields
, P and the basic state fields B0 , v0 , J0 , ω0 ,
0 , P0 as described below. These enable us to exploit two important discrete symmetries: symmetry or anti-symmetry under reflection in the equatorial plane, i.e. even or odd for scalar functions and quadrupole or dipole for vector functions, and invariance under rotation about the z-axis through an angle 2π/m, where m is an integer. For details of the scalar and vector spherical harmonic and toroidal poloidal harmonic representations used herein see Appendix A. For details of the radial discretisation and boundary conditions used see Appendix B. A major advantage of the vector spherical harmonic representation is that the vector spectral harmonic equations for the momentum, and induction equation have very compact forms. After representing the axisymmetric basic state fields v0 , ω0 , B0 , J0 ,
0 , and the perturbation fields S, T, s, t,
, in terms of vector and scalar spherical harmonics, the perturbation equations can be transformed into compact hybrid spectral equations for the perturbation coefficients Snm , Tnm , snm , tnm ,
m n . These equations are five infinite sets of partial differential equations in the spherical radius r and the time, which are coupled in harmonic degree n. They are derived in detail in Ivers and Phillips (2005). The coupling integrals of the spectral interactions between scalar and vector fields can be evaluated in terms of 3j-, 6j- and 9j-coefficients. The compact forms of the interactions translate into compact programs and one subroutine can easily evaluate all the interactions of a specific type. For example, one subroutine evaluates all six cross products, which can occur in the Lorentz, vorticity and magnetic induction terms. The resulting program uses five subroutines to evaluate all the cross products, the inner products, the diffusion operators, the time derivatives and the buoyancy term. The compact nature of the code means that bench-marking against elementary models tests the code for more complicated models. Further description of the code is given in Ivers and Phillips (2003). 2.2. Reflectional symmetry in the equatorial plane We denote a quantity as symmetric in the equatorial plane if its reflection in the plane equals the un-reflected quantity. Thus, let F be the image of a vector field F under reflection in the equatorial plane and let r be the image of the point r. Then F at the point r is symmetric (under reflection) in the equatorial plane if F (r) = F(r ). If F(r) is decomposed into components parallel F (r) and perpendicular F⊥ (r) to the plane, F(r) = F (r) + F⊥ (r), then F (r ) = F (r) − F⊥ (r).

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

Hence, if F is symmetric, then F (r ) = F (r) and F⊥ (r ) = −F⊥ (r). In terms of cartesian components, Fx (x, y, −z) = Fx (x, y, z), Fy (x, y, −z) = Fy (x, y, z), Fz (x, y, −z) = −Fz (x, y, z). Similarly, if F is antisymmetric (under reflection) in the equatorial plane, then F (r ) = −F (r) and F⊥ (r ) = F⊥ (r) or in terms of cartesian components, Fx (x, y, −z) = −Fx (x, y, z), Fy (x, y, −z) = −Fy (x, y, z), Fz (x, y, −z)=Fz (x, y, z). Since an axial dipole (quadrupole) is anti-symmetric (symmetric) in the equatorial plane, anti-symmetry (symmetry) in the equatorial plane will be termed D-symmetry (Q-symmetry). In terms of spherical polar components, F is Q-symmetric in the equatorial plane, if Fr (r, π − θ, φ) = Fr (r, θ, φ), Fθ (r, π − θ, φ) = −Fθ (r, θ, φ), Fφ (r, π − θ, φ) = Fφ (r, θ, φ). F is D-symmetric if Fr (r, π − θ, φ) = −Fr (r, θ, φ), Fθ (r, π − θ, φ) = Fθ (r, θ, φ), Fφ (r, π − θ, φ) = −Fφ (r, θ, φ). A scalar function f is symmetric (even or e-symmetric) [antisymmetric (odd or o-symmetric)], if f (r, π − θ, φ) = f (r, θ, φ)[f (r, π − θ, φ) = −f (r, θ, φ)]. The magnetohydrodynamic equations linearised about the axisymmetric basic state (B0 , v0 ,
0 ) decouple into simpler sub-problems under reflection in the equatorial plane in four cases. Denote the symmetries of the basic state fields (B0 , v0 ,
0 ) and the perturbation fields (B, v,
) by Qqe, Ddo, etc., where the lowercase q, d refer to the velocity. The four decoupled sub-problems are then, (1) (B0 , v0 ,
0 ) is Qqe-symmetric and (B, v,
) is Qqe-symmetric; (2) (B0 , v0 ,
0 ) is Qqe-symmetric and (B, v,
) is Ddosymmetric; (3) (B0 , v0 ,
0 ) is Dqe-symmetric and (B, v,
) is Dqesymmetric; (4) (B0 , v0 ,
0 ) is Dqe-symmetric and (B, v,
) is Qdosymmetric. Symmetric sub-problems for the magnetic instability case (IV) given in Section 1 and considered further in Section 4 are obtained from the symmetric sub-problems listed above by omitting
0 ,
, e and o. The first and third symmetries are also symmetries of the full nonlinear equations. The second and fourth symmetries are peculiar to the linearised magnetohydrodynamic equations and are not symmetries of the full non-linear equations. The symmetry in the equatorial plane is preserved by the scalar and vector spherical harmonic representations. From the parity property of the scalar spherical harmonics, Ynm (π − θ, φ) = (−)n−m Ynm (θ, φ), the scalar field component fnm Ynm of f, f = S, T, s, t,
, is

87

symmetric (even) if n − m is even and antisymmetric (odd) if n − m is odd. Either directly or from (A.5) the m m vector poloidal components Sm n := ∇ × ∇ × Sn Yn of m m m B and sn := ∇ × ∇ × sn Yn of v are symmetric (Qsymmetric) if n − m is even and antisymmetric (Dsymmetric) if n − m is odd. Similarly, the vector toroidal m m m components Tm n := ∇ × Tn Yn of B and tn := ∇ × m m tn Yn of v are symmetric (Q-symmetric) if n − m is odd and antisymmetric (D-symmetric) if n − m is even. The scalar and vector spherical harmonic representations also exhibit symmetry under rotations about the z-axis. It takes the simplest form in the scalar case. Spectral components of order m for the scalar functions f = S, T, s, t,
are of the form fnm (t, r)Pnm (θ)eimφ . Hence, since the basic state is axisymmetric, the problem decouples into one independent subproblem for each order m ≥ 0. Negative m need not be considered since real fields can be constructed from the real and imaginary parts of complex solutions for m ≥ 0. Thus, the four symmetric sub-problems listed above correspond to the following four chains of harmonics of the perturbation scalar fields S, T, s, t,
for each order m ≥ 0, m , sm ,
m , T m , t m , S m , (1) Qqe-symmetry: Sm m m m+1 m+1 m+2 m m m m ,...; sm+2 ,
m+2 , Tm+3 , tm+3 m , S m , sm ,
m , T m , (2) Ddo-symmetry: Tmm , tm m+1 m+1 m+1 m+2 m m m m tm+2 , Sm+3 , sm+3 ,
m+3 , . . . ; m ,
m , S m , t m , T m , (3) Dqe-symmetry: Tmm , sm m m+1 m+1 m+2 m m m m sm+2 ,
m+2 , Sm+3 , tm+3 , . . . ; m , t m , T m , t m ,
m , S m , (4) Qdo-symmetry: Sm m m+1 m+1 m+1 m+2 m m m tm+2 , Tm+3 , sm+3 ,
m m+3 , . . . .

The
terms are omitted in the magnetic instability case (IV) in Sections 1 and 4. Decoupling a problem into symmetric sub-problems has another advantage apart from the reduced size of the problem. Decoupled sub-problems may possess close solutions in some sense. For example, the eigenvalues of decoupled sub-problems may be similar. Such solutions can be difficult to resolve numerically unless they are separated into distinct symmetric sub-problems. 3. Time dependence and eigen-solution methods If the basic state fields and the boundary conditions are steady, the problem admits solutions for the perturbation fields which have a discrete spectrum in time,     B Bl (r, θ, φ)  v    v (r, θ, φ)  γl t  =  l e . l

l (r, θ, φ)

88

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

The discrete spectrum is due to the boundedness of the conducting region and the vanishing of the magnetic field at infinity. Thus spatial discretisation transforms the problem into a search for solutions of the form xeγt , where x represents the perturbation fields and (γ, x) is a solution of the generalised algebraic eigenproblem, (A + RC)x = γBx,

(3.1)

where A, B, C (B here not to be confused with the magnetic field) are matrices arising from the spatial discretisation of the governing equations and boundary conditions, R is a dimensionless parameter and the eigenvector x represents the spatial discretisation of the fields. The generalised eigen-problem can differ substantially from the standard eigen-problem; in particular, even for finite matrices A, B, C, there may be none, finitely- or infinitely-many solutions, if B is a singular matrix. The code is constructed so that contributions from any term in the momentum, induction or heat equation can be written into any of the matrices A, B or C. Hence once the matrices are constructed it is convenient to search for eigen-solutions for an ensemble of problems with varying values of R. Typically the time derivative contributions will be written to B including the differenced contributions from the ∂t Dn snm term in the poloidal momentum equation, contributions that have a common nondimensionalised factor, such as the magnetic Reynolds number Rm or as below the reciprocal Elsasser number −1 are written to C the remainder to A. Eigen-solution are sought using a variety of methods. The main method is generalised inverse iteration with complex eigenvalue shift γ0 , yn = (A + RC − γ0 B)−1 Bxn−1 ,

γn = γ0 + yˆ n ,

xn = yˆ n−1 yn , where yˆ n , which may be complex, is the element of yn of largest magnitude, and (γ0 , x0 ) is an estimate of (γ, x). A variation of generalised inverse iteration is also used. Let yn and zn be two consecutive non-normalised iterates given by yn = (A + RC − γ0 B)−1 Bxn−1 , zn = (A + RC − γ0 B)−1 Byn , which are then normalised using the component of largest modulus ξn = maxi |yn,i | and ζn = maxi |zn,i |, and let the next iterate be given by xn = zn /ζn . If the two eigenvectors of B−1 (A + RC − γ0 B) with smallest moduli ξ1 and ξ2 possess relatively close moduli, that is (|ξ1 | − |ξ2 |)/|ξ2 |  1, and the other eigenvalues no longer influence the iteration process, then for b and c

such that z + by + cx 0 the eigenvalues ξ1 and ξ2 will be estimated by the roots of ξ 2 + bξ + c = 0. In practice the method of least squares is used to find the values b and c of best fit. If the two eigenvalues ξ1 and ξ2 of least modulus are such that the moduli are not relatively close, then the normalisation ξn+1 can be used as an estimate for (λ − λ0 )−1 . The method, in practice, requires little more numerical work than performing an even number of inverse iterations and has proven valuable in resolving closely clustered eigenvalue solutions. Eigenvalues of relatively close moduli occur in many physically interesting problems, for instance, non-decoupled problems where differing symmetry solutions exhibit similar time behaviour, while tracking an eigensolution in parameter space as it changes from real to complex or coalesces or splits into separate modes. The method also provides an estimate for the next closest eigenvalue. The generalised eigen-problem is also solved using the implicitly restarted Arnoldi method implemented in ARPACK (Lehoucq and Sorensen, 1998). This method is the perhaps the best current method (Latter and Ivers, 2004) for finding 10–20 eigenvalues simultaneously in large sparse eigen-problems of the type which arise herein. The final method involves solving the critical value problem (A + Rc C)x = iωc Bx for the (real) critical value Rc and the critical angular frequency ωc , where the angular frequency ω := Imγ. The problem is nonlinear in Rc and ωc , and also under-determined, since it is homogeneous in x. It is solved by using a linear normalising equation in x and Newton–Raphson iteration (Zhang, 1994). 4. Basic state magnetic field models and results The code, which constructs the coefficient matrices A, B, C, and the generalised eigenvalue codes have been tested against many models with independently determined solutions. These include kinematic dynamo models with an insulating source-free exterior, where B0 = 0,
0 = 0 and the steady flow v0 is prescribed, which includes free decay (v0 = 0). For the simple roll flow given by v0 = t1 + s2 ,

t1 =

sin(πr) √ , 3

s2 =

r sin(πr) √ 5

the codes have been checked against the results of Dudley and James (1989) with good agreement. Indeed, the present method of using vector spherical harmonic

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

spectral equations directly in numerical programs was used, in a more specialised context, to explore the simple roll models in the original investigation. Rotating thermal convection driven by constant heat production, for which the basic state is v0 = 0 and
0 = β(a2 − r 2 )/2a (B is omitted), has also been considered. For g = −ga r, the no-slip boundary condition v = 0 and a perfect thermally-conducting boundary
0 = 0 at r = 1 the results compared correctly with the previous work of Roberts (1968). Results for a rotating magnetoconvection problem have also been compared successfully to those of Fearn (1979) and Fearn and Proctor (1983). The basic state toroidal magnetic field√is given by B0 = r sin θ1φ , 0 = −i 2/3r and J 0 = 2. The basic which yields B1,1 1,0 state velocity field is absent v0 = 0, the conducting temperature given by
0 = β(a2 − r 2 )/2a and the spherically symmetric gravitational field by g = −ga r/a. The no slip velocity condition v = 0 was used at the rigid outer boundary r = 1. For a further description of these benchmark checks see Ivers and Phillips (2003). 4.1. Magnetic instabilities We apply the method outlined in Sections 2 and 3 to the magnetic instabilities of the four axisymmetric toroidal magnetic field models B0 described below. The electrically conducting fluid is at rest, v0 = 0, in a rigid rapidly-rotating sphere V of radius a. The exterior of V is electrically insulating and the magnetic diffusivity η is constant. Associated with the basic state magnetic field B0 is an electric current density J0 . The Lorenz force due to B0 and J0 may perturb the velocity field from rest with the subsequent motion v inducing a magnetic field B. The driving energy is supplied by B0 . The problem is non-dimensionalised scaling length, the magnetic fields, time and the velocity in units of the radius a of V, the maximum B0 of |B0 | in V, the slow magnetohydrodynamic time scale τs = 2µ0 ρ a2 /B02 , where is the uniform rotation rate, and a/τs . The coefficients of the nondimensionalised inertial, viscous, and buoyancy terms are the effective Rossby number Ro := (B0 /2 a)2 /µ0 ρ, the Ekman number E := ν/2 a and the modified Rayleigh number Ra := aµ0 ρgα

/B02 , where g and 
are typical values of the gravitational acceleration and the temperature. They are set to zero to isolate the magnetic terms and their effects within the equations. The temperature and heat equation are omitted. The resulting equations are given by ∂t B = −1 ∇ 2 B + ∇ × (v × B0 )

(4.1)

1z × v = −∇P + J0 × B + J × B0 ,

89

(4.2)

where := B02 /2 ηµ0 ρ, is the Elsasser number. The fields B and v are then the dimensionless perturbation magnetic and velocity fields. The magnetostrophic form of the momentum equation holds. Results of Zhang and Fearn (1994) indicate that the solutions for a spherical shell depend continuously on Ro and E at Ro = 0 and E = 0, which provides some justification for neglecting the inertial and viscous terms in the momentum equation. In the generalized eigenproblem (3.1) the parameter R is the reciprocal Elsasser number, R := −1 . The matrix B of (3.1) is singular in the magnetic instability problem, since the Rossby number is zero. However, the singular nature of this matrix does not appear to cause numerical or other difficulties for m = 0. The matrix C is also singular, since the Ekman number is zero, but the significant requirement is that the combined matrix A + RC be invertible, which certainly appears to be the case for m = 0. If m = 0 the matrix A is singular due to Taylor’s constraint. At the rigid boundary the magnetic field matches an exterior potential field and the radial component of the velocity field vanishes. Thus the spectral boundary conditions are (B.3) (without the temperature condition). In this inviscid case the lack of a condition on the toroidal velocity field at r = 1 implies that the toroidal coefficient tnm is determined at r = 1 by its interactions with the magnetic field. Moreover, the only diffusion in the system is magnetic. If the magnetic diffusion is reduced, for example by increasing the Elsasser number, the structure of the velocity becomes finer at the boundary, making it more difficult to resolve, i.e. higher truncation levels are required. The Elsasser number is defined in terms of quantities with good estimates, ρ ≈ 104 kg m−3 , η ∼ 3 m2 s−1 ,

= 7 × 10−5 s−1 , µ0 = 4π × 10−7 T√m A−1 , apart from the magnetic field B0 . Thus, B0 = 2 ηµ0 ρ 1/2 ∼ 2 × 10−3 1/2 . If ∼ 10, then B0 ∼ 7 × 10−3 T = 70 gauss, which is achievable by the dynamo process in the Earth’s core. We consider three basic state magnetic field models. As indicated in Section 2.2 the eigen-problems for each model decouple with respect to the azimuthal wavenumber m, since the basic state is axisymmetric. For each m the eigen-problem decouples further for Model 1, which is Q-symmetric, into sub-problems of symmetry type Qq and Dd and for Model 2, which is D-symmetric, into sub-problems of symmetry types Qd and Dq. We truncate the spherical harmonic expansion at N and consider J radial grid-points in the sphere.

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

90

Fig. 1. Growth rates for magnetic field Model 1 with the Qq-symmetry and, from top to bottom, m = 1,2,3.

4.2. Model 1: B0 equatorially symmetric The basic state magnetic field is √ 3 3 r(1 − r 2 ) sin θ 1φ . B0 (r, θ) = 2

(4.3)

The magnetic field (4.3), which is shown in Fig. 5 with  = 1, is Q-symmetric with respect to the equatorial plane. This particular symmetry may be important during reversals within the Earth’s core. The field vanishes at r = 1 and on the z-axis. B0 is normalised so that maxr,θ |B0 | = 1. √ The location of the maximum field strength is r = 1/ 3, where B0,φ = 1.√From (A.2) the associated toroidal potential is T = (3 3/2)r(1 − r2 ) cos θ. Hence, the only non-zero spherical harmonic coefficient of T is T10 =

3 r(1 − r 2 ), 2

√ noting Y10 (θ, φ) = 3 cos θ. From this equation and (A.5) the vector spherical harmonic representation of 0 Y0 , where Y0 (θ, φ) = B0 is found to be B0 = B1,1 1,1 1,1 √ i 3/2 sin θ1φ and √ 3 0 = −i 2T10 = −i √ r(1 − r 2 ). B1,1 2

(4.4)

0 satisfies the analytic condition (B.2) near r = B1,1 0. By (A.6) the vector spherical harmonic represen0 Y0 + tation of the basic state current is J0 = J1,0 1,0 0 0 0 J1,2 Y1,2 , where Y1,0 (θ, φ) = cos θ1r − sin θ1θ and √ 0 2Y1,2 (θ, φ) = −2 cos θ1r − sin θ1θ . The non-zero vector spherical harmonic components are √ √ 0 0 = 3(3 − 5r 2 ), J1,2 = − 6r 2 . (4.5) J1,0

According to Section 2.2 there are two symmetric sub-problems, those of Qq-symmetry and Dd-symmetry. Growth rates γ for the two cases are shown in Fig. 1 and 2 for the azimuthal wavenumbers 1,2,3 and truncation levels N = 20, J = 80. These figures give confidence in the critical Elsasser numbers for the two sub-problems shown in Tables 1 and 2. Figs. 1 and 2 also indicate that the instability is of an ideal type. Resistive instabilities, as their name suggests, rely on finite conductivity and hence will not exist in the infinite conductivity limit with corresponding → ∞. Although the limit → ∞ was not examined in detail as part of this study. Convergence behaviour with truncation levels N and J is demonstrated for the Qq-symmetric solutions in Table 1. The results suggest that the m = 1 mode is well converged even for moderate truncation levels. However, it is interesting to note that for the m = 2 mode an apparently modest increase in the radial truncation from J = 80 to J = 110

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

91

Fig. 2. Growth rates for magnetic field Model 1 with the Dd-symmetry and, from top to bottom, m = 1,2,3.

increases the critical Elsasser number from c = 38.833 to c = 85.080. The m = 1 mode with Qq-symmetry, for which c = 12.09 and ωc = −1.419, is clearly the most unstable Table 1 Critical c and ωc for the magnetic field Model 1 and the Qq-symmetry with m = 1, 2, 3 m

N

J

−1 c

1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3

20 20 20 20 40 40 20 20 20 20 20 20 40 40 40 40 60 20 40 40 40

80 200 400 800 200 400 80 100 200 300 400 800 90 110 200 400 200 80 100 200 400

0.07887 0.08256 0.08268 0.08271 0.08249 0.08266 0.02575 0.01366 0.01352 0.01443 0.01470 0.01494 0.01834 0.01175 0.01620 0.01690 0.01620 0.00816 0.00640 0.00508 0.00520

ωc

c

−1.40158 −1.41881 −1.41877 −1.41873 −1.41927 −1.41868 −1.41571 −1.11289 −0.88753 −0.90261 −0.90749 −0.91203 −1.29639 −0.99241 −0.94953 −0.97286 −0.95031 −0.68562 −0.69779 −0.49324 −0.51275

12.679 12.113 12.096 12.091 12.123 12.098 38.833 73.195 73.961 69.313 68.037 66.933 54.527 85.080 61.734 59.163 61.713 122.497 156.225 196.863 192.487

mode of Model 1. In fact it is the most unstable mode of the first three models. The cylindrical polar (s, φ, z)components of the perturbation magnetic and velocity fields for this mode are shown in Figs. 3 and 4. The magnetic field here is given by B = Re{[Bs (r, θ)1s + Bφ (r, θ)1φ + Bz (r, θ)1z ]eimφ+γt }, where Re{Bs (r, θ)eimφ+γt } = [Re{Bs (r, θ) cos(mφ+γi t)} − Im{Bs (r, θ) sin(mφ+γi t)}]eγr t , and γr and γi are the real and imaginary parts of the growth rate γ. There are similar expressions for Bφ and Bz . If γi = 0 both functions ReBs (r, θ) and ImBs (r, θ) must be plotted to represent the s-component of B, etc. The eigenvector and hence the fields have been normalised so that the maximum azimuthal root-meanTable 2 Critical c and ωc for the magnetic field Model 1 and the Dd-symmetry with m = 1, 2, 3 m

N

J

−1 c

ωc

c

1 1 2 2 3 3

20 40 20 40 20 40

80 400 80 400 80 400

0.00263 0.00261 0.00432 0.00434 0.00314 0.00317

−0.03917 −0.03866 −0.21228 −0.21221 −0.23568 −0.23761

380.179 382.593 231.433 230.497 318.002 315.624

92

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

Fig. 3. Critical magnetic field for Model 1 with the Qq-symmetry and m = 1.

Fig. 4. Critical velocity field for Model 1 with the Qq-symmetry and m = 1.

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

square (rms) of the magnetic field,  1 B2 dφ max r,θ 2π

eγr t = √ max |Bs (r, θ)|2 + |Bφ (r, θ)|2 + |Bz (r, θ)|2 2 r,θ eγr t = √ . 2 The constraint of the rotation is evident in Fig. 4 with the motion mostly confined to the equatorial boundary. The azimuthal magnetic field is also strongly confined to the same region. For the Dd-symmetric sub-problem the most unstable mode is the m = 2 mode with c = 230.5 and ωc = −0.212. The critical Elsasser number is, however, so large that this mode is unlikely to be excited. A magnetic field model with Q-symmetry, but with a more complicated sin3 θ angular dependence, was considered by Fearn and Weiglhofer (1991b), who reported significantly higher critical Elsasser numbers, namely c = 357.3 for the Qq-symmetry and c = 349.7 for the Dd-symmetry. 4.3. Model 2: B0 equatorially antisymmetric The basic state magnetic field is B0 (r, θ) = 8r2 (1 − r 2 ) sin θ cos θ 1φ .

(4.6)

The field (4.6) has D-symmetry with respect to the equatorial plane and vanishes at r = 1, on the equator and on the z-axis. The field is normalised so that maxr,θ |B0 | = 1 √ with B0,φ = 1 at r = 1/ 2. The model has been studied previously by Fearn and Proctor (1983), Fearn and Weiglhofer (1991a,b) and Zhang and Fearn (1994).√ The associated toroidal potential is T√= (8/(3 5)) r2 (1 − r 2 )Y20 (θ, φ), where Y20 (θ, φ) = 5(3 cos2 θ − 1)/2. Hence, the only non-zero spherical harmonic coefficient of T is 8 T20 = √ r 2 (1 − r 2 ). 3 5 From this equation and (A.5) the vector spherical har0 Y0 (θ, φ), monic representation of B0 is B0 = B2,2 2,2 √ 0 where Y2,2 (θ, φ) = i 15/2 sin θ cos θ 1φ and  √ 0 2 2 0 B2,2 = −i 6T2 = −i8 r (1 − r 2 ). (4.7) 15 0 satisfies the analytic condition (B.2) near Clearly, B2,2 r = 0. By (A.6) the vector spherical harmonic repre0 Y0 + sentation of the basic state current is J0 = J2,1 2,1

93

Table 3 Critical c and ωc for the magnetic field Model 2 and the Qd-symmetry with m = 1, 2, 3 m

N

J

−1 c

ωc

c

1 1 2 2 3 3

20 40 20 40 20 40

80 400 80 400 80 400

0.02679 0.02687 0.01884 0.01885 0.01012 0.01019

−0.83285 −0.83102 −0.90035 −0.89680 −0.83692 −0.83452

37.33231 37.21331 53.08860 53.04969 98.79872 98.15159

√ 0 Y0 , where 2 2Y0 (θ, φ) = 2(3 cos θ 2 − 1)1 − J2,3 r 2,3 2,1 √ 0 2 6 sin θ cos θ 1θ and 2 3Y2,3 (θ, φ) = −3(3 cos θ − 1)1r − 6 sin θ cos θ 1θ . The non-zero vector spherical harmonic components are √ √ 2 4 3 3 0 2 0 J2,1 = 8 r(5 − 7r ), J2,3 = −8 r . (4.8) 5 15 According to Section 2.2 there are two symmetric sub-problems, those of Qd-symmetry and Dq-symmetry. Corresponding critical Elsasser numbers and frequencies are shown in Tables 3 and 4. For the Qd-symmetric sub-problem of the basic state magnetic field Model 2 with m = 2, which was considered by (Fearn and Weiglhofer, 1991b, see their Table 1 with R = 0), we found critical values c = 53.05 and ωc = −0.897 for N = 40 and J = 400. This value of c is in good agreement with the result of Fearn and Weiglhofer (1991b), who obtained c = 52.57. However, the value of ωc = −0.897 is in poor agreement with Fearn and Weiglhofer (1991b), who obtained ωc = 0.783. The values of ωc differ in sign and amplitude. [Note that ω := Imγ herein, as defined at the end of Section 3, whereas Fearn and Weiglhofer (1991a) define the angular frequency to be ωFW := −Imγ.] Fearn and Weiglhofer (1991a,b) observed what they termed numerical instabilities, which they distinguished from true instabilities by their poor convergence properties with respect to truncation level. Even for lower truncation levels no such numerical instabilities were observed with the present solution method. Recall that for larger values of the effect of magnetic diffusion, Table 4 Critical c and ωc for the magnetic field Model 2 and the Dq-symmetry with m = 1, 2, 3 m

N

J

−1 c

ωc

c

1 1 2 2 3 3

20 40 20 40 20 40

80 400 80 400 80 400

0.06556 0.06577 0.03360 0.03290 0.01649 0.01592

−1.30158 −1.30104 −0.57626 −0.55863 −0.72945 −0.69963

15.25247 15.20552 29.75859 30.39112 60.63632 62.80724

94

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

Fig. 5. Contour plots of the magnetic field Model 3 with  = 0,0.5,0.75,1. The x marks maximum eastward field.

which is the only diffusion in the model, is reduced, increasing the truncation necessary to resolve the solution fields. It is possible that a truncation level, which is too low, may give spurious results. For the Qd-symmetric sub-problem of Model 2 we found that the m = 1 mode was the most unstable mode with c = 37.21 and ωc = −0.831 for N = 40 and J = 400. The critical modes propagate eastward for m = 1, 2, 3. For Dq-symmetric sub-problem of Model 2 the most unstable mode was the m = 1 mode with c = 15.21 and ωc = −1.301 for N = 40 and J = 400. This Dqsymmetric solution is more unstable than the most unstable Qd-symmetric solution and it propagates in the same eastward direction. Thus, the Qd-symmetry requires a considerably stronger magnetic field to produce instability. This result differs from those of Fearn and Weiglhofer (1991b), whose most unstable mode is for m = 2. Our result is consistent with the results reported in Zhang and Fearn (1994), who considered a spherical shell with inner-core radius 0.4. They reported c = 11.26 and ωc = −0.2439. (Their angular frequency ωZF is related to that of the present paper by ω = (r)2 ωZF / where r in the dimensionless thickness of the shell.) Zhang and Fearn (1994) also considered the case with inner core radius reduced to zero and found c = 15.37, ωc = −1.29 for the Dq-symmetric m = 1 uncoupled mode using the same code as Fearn and Weiglhofer (1991b). The present work is in good agreement with their result and offers an independent verification of that result, confirming that Fearn and Weiglhofer (1991b) missed the most unstable mode. We note that the differing fluid geometry seems to affect both the efficiency of the mechanism and the rotation rate. It appears that the absence of an insulating inner core requires a stronger magnetic field for criticality. The fact that the Dq-symmetry is more unstable than the Qd-symmetry is also consistent with Zhang and Fearn (1994). In

fact, for m = 1, 2, 3 the Dq-symmetric solution is more unstable than the corresponding Qd-symmetric solution. 4.4. Model 3: B0 equatorially non-symmetric We also examine the basic state magnetic field of mixed symmetry, √ 3 3 B0 (r, θ) = C1 r(1 − r 2 ) sin θ 2

+ 8(1 − )r 2 (1 − r 2 ) sin θ cos θ

1φ , (4.9)

which is a convex combination of Models 1 and 2. The parameter  mixes the basic state magnetic field from pure D-symmetry ( = 0) for Model 2 to pure Qsymmetry ( = 1) for Model 1. The value of C1 determines the normalisation of the basic state magnetic field. To normalise (4.9) such that maxr,θ |B0 | = 1 we use C1 ≈ 1/0.859 at  = 0.5, C1 ≈ 1/0.856 at  = 0.75 and C1 = 1 at  = 0 and 1. Contour plots of the magnetic field for  = 0, 0.5, 0.75, 1 are given in Fig. 5. The critical Elsasser number c and angular frequency ωc for each of the values  = 0, 0.5, 0.75, 1 are shown in Table 5, where the values of C1 have been adjusted so that maxr,θ |B0 | = 1, using the fact that Eqs. (4.1) and (4.2) are invariant under the scaling ˜ 0 = βB0 , B P˜ = β2 P,

˜ = βB, B

v˜ = β2 v,

˜ = β−2 ,

γ˜ = β2 γ,

ω˜ = β2 ω,

(4.10)

by a constant, possibly negative, parameter β. As  varies from 0 to 1 the critical Elsasser number increases from c = 15.25 to c ≈ 17 and then decreases back to c = 12.68. Even though the change

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100 Table 5 Critical c and ωc for the magnetic field Model 3 of mixed symmetry with m = 1 

N

J

−1 c

ωc

c

0 0.5 0.5 0.75 0.75 1

20 20 40 20 40 20

80 80 200 80 200 80

0.06556 0.05891 0.05924 0.06489 0.06711 0.07887

−1.30158 −1.2458 −1.2397 −1.2907 −1.2952 −1.40158

15.25247 16.9760 16.8800 15.4112 14.9008 12.67949

in c is small as  varies from 0 to 1, the change in the critical mode is considerable. For  = 0 the critical mode is Dq-symmetric, whereas for  = 1 the critical mode is Qq-symmetric. Hence, the two extreme solutions of Table 5 have completely opposite equatorial symmetry for the magnetic field. The interaction mechanism changes completely, from supporting an equatorial antisymmetric magnetic field at the top of Table 5 to supporting an equatorially symmetric one at the bottom. This change takes place with little appreciable change in the critical Elsasser number. 4.5. Model 4: B0 =0 on a spherical surface within the conducting fluid To investigate the influence of a zero of the basic state magnetic field we consider B0 (r, θ) = C2 r(r02 − r 2 ) sin θ 1φ ,

(4.11)

where C2 is independent of r and θ though depends on r0 . For this Model 4 the magnetic field is zero on a surface r = r0 within the conducting fluid and on the z-axis. The magnetic field (4.11) is Q-symmetric and does not vanish on the equatorial plane. The constant C2 determines the normalisation of B0 . Here, we choose to normalise B0

95

such that maxr,θ |B0 | = 1. Thus,  1     1 − r 2 , if r0 ≤ r1 ; √ 0 (4.12) C2 :=  3 3    3 , if r0 ≥ r1 ; 2r0 √ where r1 = √ 3/2. The global maximum of |B0 | occurs at r = r0 / 3, θ = π/2 for r0 ≥ r1 and on the conducting boundary at r = 1, θ = π/2 for r0 ≤ r1 . From (A.2) the toroidal potential for Model 4 is T = C2 r(r02 − r 2 ) cos θ. Hence, there is only one non-zero spherical harmonic coefficient of T given by; C2 T10 = √ r(r02 − r 2 ). 3 From this equation and (A.5) the vector spherical har0 Y0 , where monic representation of B0 is B0 = B1,1 1,1 0 (θ, φ) is given in Section 4.2 and Y1,1  2 0 B1,1 = −i C2 r(r02 − r 2 ). (4.13) 3 0 satisfies the analytic condition (B.2) near r = 0. B1,1 By (A.6), the vector spherical harmonic representation 0 Y0 + J 0 Y0 , of the basic state current is J0 = J1,0 1,0 1,2 1,2 0 0 where Y1,0 (θ, φ) and Y1,2 (θ, φ) are given in Section 4.2. The non-zero vector spherical harmonic components of J are √ 2 2 2 0 2 0 J1,0 = C2 (3r0 − 5r ), J1,2 = − C2 r 2 . 3 3 (4.14)

Contour plots of the basic state magnetic field Model 4 are given in Fig. 6 where the field vanishes on spherical surfaces of radius r0 = 0.4, 0.6, 0.8, 0.9 and the z-axis.

Fig. 6. Contour plots of the magnetic field Model 4 with r0 = 0.4, 0.6, 0.8, 0.9.

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

96

Table 6 Critical c and ωc for the magnetic field Model 4 for m = 1 and the Qq-symmetry with r0 = 0.4, 0.6, 0.8, 0.9

Table 7 Critical c and ωc for the magnetic field Model 4 for m = 1 and the Dd-symmetry with r0 = 0.4, 0.6, 0.8, 0.9

r0

N

J

−1 c

ωc

c

r0

N

J

−1 c

ωc

c

0.4 0.4 0.4 0.4 0.6 0.6 0.6 0.6 0.8 0.8 0.8 0.9 0.9

20 40 20 40 20 40 20 40 20 40 40 20 40

80 200 80 200 80 200 80 200 80 200 200 80 200

0.000119 0.000102 0.005246 0.005296 0.00183 0.00182 0.00532 0.00531 0.06529 0.06528 0.01110 0.15786 0.15904

−0.00211 −0.00151 0.78486 0.82278 −0.01482 −0.01488 0.75376 0.80311 −0.92835 −0.92544 0.54991 −2.61963 −2.62003

8399.012 9805.144 190.624 188.821 545.605 548.034 187.868 188.283 15.316 15.318 90.100 6.335 6.288

0.4 0.4 0.6 0.6 0.6 0.8 0.8 0.8 0.9 0.9 0.9

20 40 20 40 40 20 40 40 20 40 40

80 200 80 200 200 80 200 200 80 200 200

0.0033133 0.0033742 0.0053406 0.0053067 0.00018316 0.0058373 0.0060365 0.0017401 0.0070965 0.0062841 0.0034217

0.99465 0.98296 0.47695 0.47738 −0.00384 0.49458 0.50994 −0.005547 0.15594 0.16955 −0.14626

301.818 296.368 187.246 188.443 5459.603 171.312 165.658 574.667 140.915 159.132 292.249

From Section 2.2 the basic state magnetic field of Model 4 decouples the perturbation fields into two subproblems, those of Qq-symmetry and Dd-symmetry. For Model 4 we restrict our attention to the m = 1 uncoupled sub-problem. The critical Elsasser numbers and frequencies are given for the azimuthal wavenumber m = 1 at selected values of radius r0 = 0.4, 0.6, 0.8, 0.9 within the conducting sphere: for the Qq-symmetry in Table 6 and for the Dd-symmetry in Table 7. In several cases the critical Elsasser number of the most unstable mode which propagates in the opposite longitudinal direction is also given. If r0 is reduced to 0.9 or 0.8, it is clear from Tables 6 and 7 that the critical mode of Model 4 for m = 1 has Qq-symmetry and propagates eastward with c ∼ 10 and ωc ∼ −1. The critical Elsasser number of the most unstable westward Qq-mode is c ∼ 100. However, as r0 is reduced further to 0.6 or 0.4, the westward Qqmode becomes the critical mode with c ∼ 190 and ωc ∼ 1. The critical Elsasser number c of the most

unstable eastward Qq-mode increases very rapidly as r0 is decreased and the magnitude of the critical angular frequency ωc decreases rapidly toward zero. Generally, for m = 1 the modes of Dd-symmetry are more stable. However, at r0 = 0.6 the most unstable westward Ddmode has essentially the same critical Elsasser number c ≈ 188 as the Qq-mode. Fig. 7 shows meridional plots of the azimuthal-rms magnetic fields with the Qq-symmetry for m = 1 at r0 = 0.8 and r0 = 0.6. The magnetic field of the westward (eastward) propagating modes are strongly confined to the region of the conducting fluid outside (inside) the radius r0 . The corresponding velocities of the westward modes are dominated by strong azimuthal jets tightly confined to the equatorial boundary. At r0 = 0.8 the velocity of the most unstable eastward mode possesses a weaker but still dominant azimuthal component in the outer half of the fluid near the equator. However, as r0 is reduced from r0 = 0.8 to 0.7 the azimuthal flow concentrates into a strong jet near the equator.

Fig. 7. Critical azimuthal-rms magnetic fields for Model 4 with the Qq-symmetry, m = 1 and r0 = 0.6, 0.8. The two modes on the left (right) propagate westward (eastward).

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

5. Concluding remarks We have found magnetic instabilities for the simplest possible azimuthal magnetic fields (4.3) and (4.6), which are of Q- and D-symmetry, respectively, and vanish on the boundary of the conducting fluid sphere. The most unstable mode, with c = 12.098, was the m = 1 Qqsymmetric mode of the magnetic field (4.3). We have also investigated the effect that a zerosurface of the basic state magnetic field in the conducting region has on the instabilities. The most unstable mode is always equatorially symmetric. When the dimensionless magnetic-zero radius is decreased below about 0.7 the most unstable mode switches from eastward to westward propagation. The magnetic field perturbation in the westward (eastward) case is strongly confined to the region outside (inside) the magnetic-zero radius. In each of the first three models considered in Section 4, the basic state magnetic field vanishes at the boundary r = 1, whereas Model 4 does not vanish at the boundary unless r0 = 1. Future work will investigate the effects of the location of zeros in a basic state magnetic field, which also vanishes on the boundary of the conducting fluid, and the dimensionality of the zero, i.e. whether it is a surface, curve or point. The critical flows within a dimensionless radius of about 0.4 were found to be weak, which suggests that the effect of an electrically conducting inner core on the critical modes might be negligible. Further work is required, however, as there is the possibility that such an inner core might destabilise a subcritical mode and hence lower the critical Elsasser number. An insulating inner core may magnetically effect the critical modes. We note in Section 4.3, comparing our results with Zhang and Fearn (1994), that an insulating inner core apparently decreases the critical Elsasser number. Appendix A. Field representations For definitions of the scalar and vector sperical harmonics used here see Phillips and Ivers (2001) and James (1976). The scalar and vector spherical harmonics are orthonormal with respect to the inner products, 1 (F, G) : = 4π (F, G) : =

1 4π



The basic state vector fields are expanded in vector spherical harmonics,  0,m m Fn,n Y , F0 = v0 , B0 , ω0 , J0 , g0 , F0 = 1 n,n1 n,n1 ,m

(A.1) where the summations are over all non-negative integers n, n1 = n, n ± 1 and m = −n : n. In numerical applications, the degree n and order m are truncated at N so that |m| ≤ n ≤ N. The index n1 is truncated at N + 1. Toroidal–poloidal representations are used for the perturbation vector fields v, B, v = ∇ × {tr} + ∇ × ∇ × {sr}, B = ∇ × {T r} + ∇ × ∇ × {Sr}.

(A.2)

The incompressible and solenoidal conditions are satisfied identically in the toroidal–poloidal potentials t, s, T, S. The vorticity ω and current J have the derived representations, ω = ∇ × {−(∇ 2 s)r} + ∇ × ∇ × {tr}, µ0 J = ∇ × {−(∇ 2 S)r} + ∇ × ∇ × {T r}.

(A.3)

The potentials t, s, T, S and the temperature
are expanded in scalar spherical harmonics,  f = fnm Ynm f = t, s, T, S,
. (A.4) n,m

As with expansion (A.1) in numerical applications the degree n and order m are truncated at N, |m| ≤ n ≤ N. The vector spherical harmonic coefficients of the m are simply related to the coefficients magnetic field Bn,n 1 Tnm and Snm of the toroidal–poloidal magnetic potentials by   n   (n + 1) ∂nn−1 Snm , if n1 = n − 1;   2n + 1   √ m −i n(n + 1)Tnm , if n1 = n; Bn,n := 1     n + 1 n+1 m   ∂ Sn , if n1 = n + 1; n 2n + 1 n (A.5)

∗

FG d ,

97

F · G∗ d , m

respectively, where d = sin θ dθdφ; i.e. (Ynmαα , Ynββ ) = n m m n n mβ δnβα δmβα and (Ynmαα,n1α , Ynββ,n1 β ) = δnβα δn1β 1α δmα .

where ∂nn1 := ∂r + (n + 1)/r if n1 = n − 1, and ∂nn1 := ∂r − n/r if n1 = n + 1. The n1 = n ± 1 (n1 = n) components represent the poloidal (toroidal) field. Similar equations relate the velocity vector coefficients vm n,n1 to the velocity potential scalar coefficients tnm and snm . The m of the elecvector spherical harmonic coefficients Jn,n 1 trical current J are related to the coefficients Tnm and Snm

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

98

of the magnetic field potentials by   n   (n + 1) ∂nn−1 Tnm ,   2n + 1     if n = n − 1; 1 m √ µ0 Jn,n1 :=  i n(n + 1)Dn Sn , if n1 = n;         n n + 1 ∂nn+1 Tnm , if n1 = n + 1. 2n + 1 (A.6) where n := rr + 2r∂r − n(n + 1). The vector m of the vorticity spherical harmonic components ωn,n 1 are related to the velocity potentials tnm and snm by similar equations. The most convenient method to prescribe basic state magnetic and velocity fields is to specify their toroidal–poloidal potentials and then use (A.5) and (A.6) and their velocity analogues to derive the associated vector spherical harmonic components. A poor choice of the angular behaviour of the basic state fields can give poor convergence properties, e.g. an infinite number of terms are required in the expansion (A.4) for f = sin θ but only one term for f = cos θ. r2 D

r2 ∂

Appendix B. Radial discretisation and boundary conditions A major advantage of the compact hybrid spectral equations is that they contain only six types of radial derivative expression made up from coefficients of the basic state, the perturbation fields and their radial derivatives. The six types are n

n

n

f0 ∂nβ1β f, f0 Dnβ f, ∂nβ1β (f0 f ), ∂nn1 (f0 ∂nβ1β f ), ∂nn1 (f0 Dnβ f ), Dn Dn f, where f0 is a scalar or vector spherical basic state component; f is a coefficient of the perturbation fields m m m m m tnββ , snββ , Tnββ , Snββ ,
nββ ; the subscript nβ is the harmonic degree of the perturbation field; the superscript n1β = nβ ± 1; the index n is the harmonic degree of the spectral equation; and n1 = n ± 1. Some care is needed with the differentiability of the basic state and perturbation fields and their scalar and vector spherical harmonic coefficients at r = 0. There is nothing special about the origin in cartesian coordinates x, y, z, but it is a coordinate singularity in spherical polar coordinates. Thus, for any integer k ≥ 0, rn+2k Ynm (θ, φ) is an analytic function of x, y, z at the origin, since it is a homogeneous polynomial in x, y, z, and rn+2k+1 Ynm (θ, φ) has singular derivatives in x, y, z of order n + 2k + 1 at the origin, but both are analytic functions of r, θ, φ at r = 0. For our purposes differen-

tiability means with respect to cartesian coordinates. We can guarantee the differentiability of scalar and vector fields at the origin from the form of the scalar and vector spherical harmonic components. Indeed, the scalar function f = S, T, s, t,
, is an analytic function of x, y and z at the origin, if the coefficients of its spherical harmonic expansion f = n,m fnm (r)Ynm (θ, φ) are of the form, fnm (r) = rn

∞ 

fnm,k r2k ,

(B.1)

k=0

where fnm,k is independent of r. Similarly, from the m definition of the vector spherical harmonic Yn,n as a 1 linear combination of scalar spherical harmonics Yn1 of order n1 and complex cartesian unit vectors eµ (see James, vector spherical harmonic expansion 1976), the m (r)Ym (θ, φ) is an analytic function F = n,n1 ,m fn,n n,n1 1 of x, y, z at the origin, if m (r) = rn1 Fn,n 1

∞ 

m,k 2k Fn,n r , 1

(B.2)

k=0 m,k is independent of r. where Fn,n 1 The sufficient conditions (B.1) and (B.2) on the scalar and vector spherical harmonic components are important as they allow us to specify basic state fields v0 , B0 , ω0 , J0 , which are guaranteed to be analytic at the origin. It is possible to specify basic state fields with unbounded derivatives in x, y, z at the origin, which may affect numerical convergence with higher resolution in the radial direction. Such fields can give misleading results due the non-analytic behaviour of the basic state at the origin. The spherical radial dependence is discretised on a uniform radial grid rj = jh, j = 1 : J, with stepsize h = 1/J, using second-order finite differences. For details of the second-order explicit centred differenced schemes see Ivers and Phillips (2003). At the boundary of the conducting sphere the magnetic field is matched to a potential field in the exterior, which vanishes at infinity, there is no flow across the rigid boundary, n · v = 0, and the boundary is a perfect thermal conductor. These conditions imply that the perturbation field coefficients satisfy

∂nn−1 Snm = Tnm = 0,

snm = 0,


m n = 0, at r = 1.

(B.3)

These conditions apply if the flow is inviscid, in which case there are no constraints on the tangential components of the velocity v. If the flow is viscous and there

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

is no tangential slip along the boundary, n × v = 0, the further conditions ∂r snm = tnm = 0,

at

r=1

hold. The following explicit one-sided boundary formulas are used to difference the third- and fourth-order derivatives at rJ−1 = 1 − h, i.e. one grid-point in from the boundary, (3)

fJ−1 =

fJ−4 − 6fJ−3 + 12fJ−2 − 10fJ−1 + 3fJ 2h3 1 + f (5) (ξ)h2 (B.4) 4 (1)

(3)

fJ−1 =

−fJ−3 + 9fJ−1 − 8fJ + 6hfJ 3h3 3 − f (5) (ξ)h2 20

References

(1)

+192fJ−1 − 113fJ + 60hfJ 12h4 1 − f (7) (ξ)h3 . (B.6) 21 Derivatives up to third-order are differenced at the boundary rJ = 1 by the explicit one-sided boundary schemes; (4)

fJ−1 =

(1)

(2)

fJ =

(3)

fJ =

fJ−2 − 4fJ−1 + 3fJ 1 + f (3) (ξ)h2 2h 3 −fJ−3 + 4fJ−2 − 5fJ−1 + 2fJ h2 11 + f (4) (ξ)h2 12

(B.7)

(B.8)

3fJ−4 − 14fJ−3 + 24fJ−2 − 18fJ−1 + 5fJ 2h3 7 + f (5) (ξ)h2 (B.9) 4 (1)

(2)

fJ =

−fJ−2 + 8fJ−1 − 7fJ + 6hfJ 2h2

1 + f (4) (ξ)h2 6 (B.10) (1)

(3)

fJ =

In general, there is one equation for each of the coefficients tnm , snm , Tnm , Snm ,
m n , at each grid-point rj , j = 1 : J − 1. If the value of any coefficient at rJ is not fixed by the boundary conditions, then there is an additional equation for that coefficient at rJ . Finally, there is a further equation for
00 at r0 . The discretisation produces a banded matrix with the narrowest band, in general, if all coefficients and their equations at each grid point are blocked together. The resulting matrix is notionally block-pentadiagonal due to the centred difference formulas. The one-sided formulas (B.3)–(B.11) distort this slightly with a block of J − 4 coefficient values and J − 1 equations, and a block of J − 3 coefficient values and J equations.

(B.5)

−3fJ−4 + 32fJ−3 − 108fJ−2

fJ =

99

fJ−3 − 6fJ−2 + 15fJ−1 − 10fJ + 6hfJ h3 11 + f (5) (ξ)h2 . (B.11) 20

Bullard, E., Gellman, H., 1954. Homogeneous dynamos and terrestrial magnetism. Phil. Trans. R. Soc. Lond. A 247, 213–278. Dudley, M.L., James, R.W., 1989. Time-dependent kinematic dynamos with stationary flows. Proc. R. Soc. Lond. A 425, 407– 429. Fearn, D.R., 1979. Thermal and magnetic instabilities in a rapidly rotating fluid sphere. Geophys. Astrophys. Fluid Mech. 14, 103– 126. Fearn, D.R., Proctor, M.R.E., 1983. Hydromagnetic waves in a differentially rotating sphere. J. Fluid Mech. 128, 1–20. Fearn, D.R., Weiglhofer, W.S., 1991a. Magnetic instabilities in rapidly rotating spherical geometries I. From cylinders to spheres. Geophys. Astrophys. Fluid Dynam. 56, 159–181. Fearn, D.R., Weiglhofer, W.S., 1991b. Magnetic instabilities in rapidly rotating spherical geometries II. More realistic fields and resistive instabilities. Geophys. Astrophys. Fluid Dynam. 60, 257– 294. Fearn, D.R., Weiglhofer, W.S., 1992. Magnetic instabilities in rapidly rotating spherical geometries III. The effect of differential rotation. Geophys. Astrophys. Fluid Dynam. 67, 163–184. Ivers, D.J., Phillips, C.G., 2003. A vector spherical harmonic spectral code for linearised magnetohydrodynamics. Anziam J. 44E, C423– C442. Ivers, D.J., Phillips, C.G., 2005. Scalar and vector spherical harmonic spectral equations of rotating non-linear and linearised magnetohydrodynamics, preprint. James, R.W., 1974. The spectral form of the magnetic induction equation. Proc. R. Soc. Lond. A 340, 287–299. James, R.W., 1976. New tensor sperical harmonics, for application to the partial differential equations of mathematical physics. Phil. Trans. R. Soc. Lond. A 281, 195–221. Latter, H., Ivers, D.J., 2004. Kinematic roll dynamo computations at large magnetic Reynolds number. ANZIAM J. 45E, C905– C920. Lehoucq, R.B., and Sorensen, D.C., 1998. Arpack Users Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. ISBN 0898714079, SIAM, Philadelphia, 140pp. Phillips, G.C., Ivers, D.J., 2000. Spherical anisotropic diffusion models for the Earth’s core. Phys. Earth Planet. Int. 117, 209–223.

100

C.G. Phillips, D.J. Ivers / Physics of the Earth and Planetary Interiors 153 (2005) 83–100

Phillips, G.C., Ivers, D.J., 2001. Spectral interactions of rapidlyrotating anisotropic turbulent viscous and thermal diffusion in the Earth’s core. Phys. Earth Planet. Int. 128, 93–107. Roberts, P.H., 1968. On the thermal instability of a self-gravitating fluid sphere containing heat sources. Phil. Trans. R. Soc. Lond. A 263, 93–117.

Zhang, K., 1994. On coupling between the Poincare equation and the heat equation. J. Fluid Mech. 268, 211–229. Zhang, K., Fearn, D.R., 1994. Hydromagnetic waves in rapidly rotating spherical shells generated by magnetic toroidal decay modes. Geophys. Astrophys. Fluid Dynam. 77, 133–157.