EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 519 (September 2010) A&A, 519 (2010) A80 Full HTML
Free Access
Issue
A&A
Volume 519, September 2010
Article Number A80
Number of page(s) 7
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/200913702
Published online 16 September 2010
A&A 519, A80 (2010)

An efficient method for computing the eigenfunctions of the dynamo equation

M. Schrinner1,2 - D. Schmitt1 - J. Jiang1 - P. Hoyng3

1 - Max-Planck-Institut für Sonnensystemforschung, Max-Planck-Str. 2, 37191 Katlenburg-Lindau, Germany
2 - MAG (ENS/IPGP), LRA, École Normale Supérieure, 24 rue Lhomond, 75252 Paris Cedex 05, France
3 - SRON Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands

Received 19 November 2009 / Accepted 19 May 2010

Abstract
Aims. We present an elegant method of determining the eigensolutions of the induction and dynamo equations in a fluid embedded in a vacuum.
Methods. The magnetic field is expanded in a complete set of functions. The new method is based on the biorthogonality of the adjoint electric current and the vector potential with an inner product defined by a volume integral over the fluid domain. The advantage of this method is that the velocity and the dynamo coefficients of the induction and the dynamo equation do not have to be differentiated and thus even numerically determined tabulated values of the coefficients produce reasonable results.
Results. We provide test calculations and compare with published results obtained by the classical treatment based on the biorthogonality of the magnetic field and its adjoint. We especially consider dynamos with mean-field coefficients determined from direct numerical simulations of the geodynamo and compare with initial value calculations and the full MHD simulations.

Key words: magnetohydrodynamics (MHD) - magnetic fields - methods: numerical

1 Introduction

The generation and evolution of magnetic fields in cosmic bodies like the planets and stars is generally thought to be governed by induction processes due to motions in their electrically conducting fluid interior. The magnetic field $\vec{B}$ is described by the induction equation

\begin{displaymath}\frac{\partial\vec{B}}{\partial t} = \vec{\nabla}~\!\times\!~D\vec{B} \end{displaymath} (1)

where

\begin{displaymath}D\vec{B}= \vec{u}~\!\times\!~\vec{B}- \eta\vec{\nabla}~\!\times\!~\vec{B}. \end{displaymath} (2)

Here $\vec{u}$ represents the velocity and $\eta$ the magnetic diffusivity.

In the framework of mean-field theory (e.g. Krause & Rädler 1980; Moffatt 1978), $\vec{u}$ and $\vec{B}$ are considered as means, e.g. ensemble averaged, quantities, whereas the action of the small-scale turbulent flow on the mean magnetic field is parametrised by the so-called dynamo coefficients, $\vec{\alpha}$ and $\vec{\beta}$. They are, in general, tensors of second and third ranks, respectively. We use the following compact notation of the mean field coefficients, which include the so-called $\gamma$, $\delta$, and $\kappa$-effects, see e.g. Rädler (1980). Then, the operator D reads

\begin{displaymath}D\vec{B}= \vec{u}~\!\times\!~\vec{B}+ \vec{\alpha}\!\cdot\!\v... ...a}:(\vec{\nabla}\vec{B}) - \eta\vec{\nabla}~\!\times\!~\vec{B} \end{displaymath} (3)

instead of (2), and acts on the mean magnetic field. Except for the additional $\vec{\alpha}$ and $\vec{\beta}$ terms in the D operator, the induction and the dynamo equation are formally equivalent. Thus, the new method presented here equally applies to both.

The dynamo region is located in a flow domain V with exterior vacuum E. In this work we assume V to be either a sphere or a spherical shell. The magnetic field $\vec{B}$ is continuous through the boundary $\partial V$ and potential in E.

In kinematic dynamo theory, all coefficients ($\vec{u}$, $\vec{\alpha}$, $\vec{\beta}$, and $\eta$) are assumed given and independent of the magnetic field. Thus the dynamo equation is linear in the magnetic field and can be solved by considering an eigenvalue problem

\begin{displaymath}\lambda\vec{B}= \vec{\nabla}\!\times\!D\vec{B} \end{displaymath} (4)

with eigenvalues $\lambda$ describing the time evolution proportional to $\exp(\lambda t)$ of the magnetic field $\vec{B}$.

Many studies have been made of the eigenvalues of the dynamo operator for various celestial bodies and with many forms of the dynamo coefficients (e.g.Bullard & Gellman 1954; Steenbeck & Krause 1969; Kumar & Roberts 1975; Gubbins et al. 2000; Deinzer et al. 1993; Gubbins 1973; Jiang & Wang 2007; Schubert & Zhang 2001; Deinzer & Stix 1971; Livermore & Jackson 2004; Dudley & James 1989; Roberts 1960; Jiang & Wang 2006; Roberts & Stix 1972; Roberts 1972; Schmitt & Schüssler 1989; Livermore & Jackson 2005). Often the coefficients are approximated by simple analytical functions of position, and their tensorial character is disregarded. Recently, the test-field method, developed by Schrinner et al. (2005,2007) (see also Ossendrijver et al. 2002,2001), allows one to determine all tensorial components of $\vec{\alpha}$ and $\vec{\beta}$ directly from self-consistent numerical simulations (Käpylä et al. 2009; Brandenburg et al. 2008). These coefficients are sometimes strongly varying functions of position. This may introduce large errors because the dynamo operator $\vec{\nabla}\!\times\!D$ involves differentiation of the dynamo coefficients, and these are only available as numerically determined tabulated values.

In this paper we present a new method that does not require differentiation, so it is also applicable to numerically determined dynamo coefficients. The method is based on the biorthogonality of the electric current and the vector potential with an inner product defined by a volume integral over the fluid. This property has already been noted by Rädler & Bräuer (1987), Hoyng (1988), Fuchs et al. (1993), Hoyng & Schutgens (1995), and Rädler et al. (2002).

The method is described in detail in Sect. 2. Extensive test calculations have been performed and compared with published results by other eigenvalue methods. Some of these tests are presented in Sect. 3. In Sect. 4 we apply the new method to eigenmodes of the dynamo operator with coefficients obtained from geodynamo models. Our conclusion are drawn in Sect. 5.

2 Eigenvalue problem

We expand the field $\vec{B}$ of the dynamo in a complete set of functions $\vec{b}_i(\vec{r})$:

\begin{displaymath}\vec{B}= \sum_i e_i\vec{b}_i(\vec{r}) \equiv e_i\vec{b}_i. \end{displaymath} (5)

Here and in the following we make use of the summation convention for two identical indices. The expansion functions are often eigenfunctions of some differential operator. Since this operator is, in general, not self-adjoint, the functions are not orthogonal. This problem is handled by using the adjoint set $\vec{\hat{b}}_k(\vec{r})$, with the following inner product

\begin{displaymath}(\vec{\hat{b}}_k,\vec{b}_i)_X \equiv \int_X \vec{\hat{b}}_k\!\cdot\!\vec{b}_i ~ {\rm d}^3\vec{r}= \delta_{ki}. \end{displaymath} (6)

The integration volume X can be either the whole space V+E or the fluid domain V alone. The base functions $\vec{\hat{b}}_k(\vec{r})$ and $\vec{b}_i(\vec{r})$ constitute a biorthogonal set. For a given set of functions $\vec{b}_i$, the adjoint set $\vec{\hat{b}}_k$ depends on the choice of the integration domain X, so in principle we have two different sets $\vec{\hat{b}}_k$, one for X=V and one for X=V+E.

Later we adopt the free magnetic decay modes, for which the base functions $\vec{b}_i$and their adjoints $\vec{\hat{b}}_k$ are known. But at this point there is no need to specify which set $\vec{b}_i$ we actually use.

2.1 Biorthogonal sets

Starting from a set $\vec{\hat{b}}_k$ and $\vec{b}_i$ that is biorthogonal on V+E, a very useful biorthogonal set on V is provided by the associated electric current $\vec{\hat{\jmath}}=\vec{\nabla}\times \vec{\hat{b}}$ and the vector potential $\vec{a}$ where $\vec{\nabla}\times\vec{a}=\vec{b}$, with the inner product

\begin{displaymath}(\vec{\hat{b}}_k,\vec{b}_i)_{V+E} = (\vec{\hat{\jmath}}_k,\vec{a}_i)_V = (\vec{\hat{a}}_k,\vec{j}_i)_V =\delta_{ki} . \end{displaymath} (7)

Here we have absorbed a factor of $4\pi/c$ in the definition of the current $\vec{j}$. The relation (7) is derived with the help of the vector identity $\vec{\nabla}\cdot(\vec{a}_i\times\vec{\hat{b}}_k)=\vec{\hat{b}}_k\cdot(\vec{\na... ...c{\hat{b}}_k) =\vec{\hat{b}}_k\cdot\vec{b}_i-\vec{a}_i\cdot\vec{\hat{\jmath}}_k$ and a volume integration over V+E; $\vec{a}$and $\vec{b}$ go fast enough to zero at infinity. Surface integrals vanish because the field and the vector potential are continuous through $\partial V$. Volume integrals containing currents are restricted to V since $\vec{j}=0$ in E. The inner product (7) is invariant under a gauge transformation $\vec{a}\to\vec{a}+\vec{\nabla}\psi$ because $\int_V\vec{\hat{\jmath}}\cdot\vec{\nabla}\psi~{\rm d}^3\vec{r}= \int_V\vec{\nab... ...d}^3\vec{r}=\int_{\partial V}\psi\vec{\hat{\jmath}}\cdot{\rm d}^2\vec{\sigma}=0$, as currents and their adjoints run parallel to the boundary.

Electric currents and vector potentials thus form a biorthogonal set on V. This is essential for the new eigenvalue method presented in Sect. 2.3.

2.2 Classical eigenvalue method

Inserting the expansion (5) in the dynamo eigenvalue Eq. (4) yields

\begin{displaymath}\lambda e_i\vec{b}_i = \vec{\nabla}\!\times\!(De_i\vec{b}_i). \end{displaymath} (8)

Subsequently, we take the inner product (6) based on V with the adjoint magnetic field. This leads to

\begin{displaymath}\lambda e_k = M_{ki}e_i \;\;{\rm with}\;\; M_{ki}=(\vec{\hat{b}}_k,\vec{\nabla}\!\times\!D\vec{b}_i)_V. \end{displaymath} (9)

A partial integration to shift the curl from the second to the first term, as done in (7) above and used in the new method below, is not possible because the surface term $\int_{\partial V}(D\vec{b}_i\!\times\!\vec{\hat{b}}_k)\cdot{\rm d}^2\vec{\sigma}$ need not vanish here.

We mention as an aside that the magnetic field is often decomposed in its poloidal and toroidal components (see Appendix A) after which the dynamo equation is formulated in terms of the defining scalars P and T. If the dynamo coefficients possess certain symmetry properties, the solutions can be split into two independent subsets, describing magnetic fields symmetric and antisymmetric with respect to the equator.

2.3 New eigenvalue method

We start again with (8), which we uncurl to obtain

\begin{displaymath}\lambda e_i\vec{a}_i = De_i\vec{b}_i +\vec{\nabla}\psi. \end{displaymath} (10)

Taking now the inner product (7) with the adjoint current results in

\begin{displaymath}\lambda e_k = N_{ki}e_i \;\;{\rm with}\;\; N_{ki}=(\vec{\hat{\jmath}}_k,D\vec{b}_i)_V. \end{displaymath} (11)

The gradient term drops out as discussed in Sect. 2.1 above. The corresponding adjoint functions $\vec{\hat{b}}_k$ here are different from those in Sect. 2.2 as they pertain to a different inner product.

The matrices Mki and Nki have the same eigenvalues $\lambda$. The advantage of the new method using Nki in (11) instead of Mki in (9) is that no differentiation of the operator D is required, so even numerically computed or tabulated values of $\vec{u}$, $\vec{\alpha}$, and $\vec{\beta}$ produce accurate results.

2.4 Choice of b$_{\sf i}$ and numerical handling of Eq. (11)

For the set of base functions, we adopt the free magnetic decay modes whose magnetic fields $\vec{b}_i$ are known analytically in V+E in terms of the defining scalars P and T as described in Appendix A. The decay modes are continuous through $\partial V$ and potential in E, so they satisfy the boundary conditions. They are characterised by three numbers, the radial order n, the latitudinal degree l, and the azimuthal order m.

Another advantage of the decay modes is that they are self-adjoint on V+E so that the adjoint functions are the complex conjugates $\vec{\hat{b}}_k=\vec{b}^*_k$ and likewise $\vec{\hat{\jmath}}_k=\vec{j}^*_k$. Normalisation on V+E, i.e. $(\vec{b}_k^*,\vec{b}_i)_{V+E}=(\vec{j}_k^*,\vec{a}_i)_V=\delta_{ki}$, is thus readily achieved, see Appendix A.

The computation of the matrix elements Nki is now straightforward. Once we know the matrix elements, the eigenvalue problem (11) is solved numerically using LAPACK routines (http://www.netlib.org/lapack), and we obtain the eigenvalues $\lambda_k$ and eigenvectors $\{e_{ki}\}$, such that

\begin{displaymath}\vec{B}_k=e_{ki}\vec{b}_i\end{displaymath} (12)

is eigenfunction of $\vec{\nabla}\times D$ with eigenvalue $\lambda_k$. Each mode k contains, in general, a mixture of n, l, and mvalues.

Table 1:   Eigenvalues of the fundamental dipolar mode and the fifth and tenth overtones of the $\alpha ^2$-sphere.

In the following we consider only velocities and dynamo coefficients that are independent of azimuth $\varphi$, but this is not a necessary constraint. Thus each value of m can be treated separately. Although we present only results for m=0 here, we have tested and applied other values of m as well. We employ the robust Gauss-Legendre quadrature in r and $\cos\th$ to compute the matrix elements since the basis functions are heavily oscillatory in r for high values of n and in $\theta$ for high degree l. For the Gauss-Legendre integration we used 66 quadrature points here in the radial and 80 in the latitudinal direction, respectively. In general, this depends of course on the required resolution.

2.5 Adjoint eigenfunctions

We now show how one may construct the adjoint set of eigenfunctions $\vec{\hat{B}}_p$ of a set of eigenfunctions $\vec{B}_i$ of the dynamo operator $\vec{\nabla}\!\times\!D$. Although these adjoints are not needed in the present paper, they appear in applications. For example, let $\vec{B}$ be the actual magnetic field of the dynamo, then it is often advantageous to expand $\vec{B}$ in dynamo eigenfunctions, i.e., $\vec{B}(\vec{r},t)=\sum_i c_i(t)\vec{B}_i(\vec{r})$. To find the coefficients ci(t) we use the adjoint set, to find $c_i=(\vec{\hat{B}}_i,\vec{B})_{V+E}=(\vec{\hat{J}}_i,\vec{A})_V$.

This illustrates that we need the adjoints $\vec{\hat{B}}_p$, and these may be constructed as follows. Let $\vec{B}_k=e_{ki}\vec{b}_i$ be the representation of $\vec{B}_k$in terms of the self-adjoint magnetic decay modes as above. Then we write $\vec{\hat{B}}_p=f^*_{pi}\vec{b}^*_i$ and $\vec{\hat{J}}_p=f^*_{pi}\vec{j}^*_i$, and we require

\begin{displaymath}\delta_{pk}=(\vec{\hat{J}}_p,\vec{A}_k)_V=(f^*_{pi}\vec{j}^*_i,e_{kj}\vec{a}_j)_V= f^*_{pi}e_{kj}\delta_{ij}=f^*_{pi}e_{ki}~. \end{displaymath} (13)

Therefore ${\rm f}^\dagger={\rm e}^{-1}$ in matrix notation, and we find a unique biorthogonal set. Here, $\dagger$ indicates the Hermitean adjoint, ${\rm f}^\dagger=({\rm f}^{\rm T})^*$, where T indicates the transposed and * complex conjugation.

Three important messages follow from this construction: (i) the adjoint of $\vec{\nabla}\times\vec{B}$ is $\vec{\nabla}\times\vec{\hat{B}}$, that is, the adjoint operation commutes with $\vec{\nabla}$; (ii) to obtain the adjoint eigenfunctions, it is not necessary to know the explicit form of the adjoint dynamo operator $\vec{\nabla}\times\hat{D}$; and (iii) the eigenfunctions and their adjoints have the same boundary conditions because they are a linear combination of the decay modes and their complex conjugates, respectively.

Table 2:   Eigenvalues of the fundamental mode for the $\vec{t_1s_2}$ flow for two magnetic Reynolds numbers $R_{\rm m}$.

3 Test results

3.1 $\alpha ^2$-sphere

We first consider the so-called $\alpha ^2$-sphere of unit radius r0=1, represented by $\vec{u}=\vec{\beta}=\vec{0}$, $\alpha_{ij}=R_\alpha\delta_{ij}$, and $\eta=1$, which also can be treated analytically (Krause & Rädler 1980, Chap. 14). The eigenvalues are independent of azimuth m, and the eigenfunctions decouple in latitudinal quantum number l. For $R_\alpha\neq 0$, the modes couple in radial number n, as they do between the poloidal and toroidal components.

For $R_\alpha=4.493409458$[*], the first mode is a stationary dipole, while the overtones decay with the rates given by (Hoyng & van Geffen 1993, HvG93). We successfully reproduced the fundamental mode and the overtones. In Table 1 we consider the convergence in the eigenvalues as a function of the maximum radial number $n_{\rm max}$ for some dipolar (l=1) modes. Higher l modes behave similarly. We also reproduced the eigenfunction plots as provided by Krause & Rädler (1980).

Table 3:   Eigenvalues of the first and fifth modes of the $\alpha ^2\Omega $-dynamo (see Sect. 3.3).

3.2 Spherical flows

As a next test, we apply the spherical stationary $\vec{t}_1 \vec{s}_2$ (MDJ) flow of Livermore & Jackson (2004) which is given by

                          $\displaystyle \vec{u}$ = $\displaystyle u_0K^{-1}\vec{\nabla}\!\times\!\left(r^2(1-r^2)P_1^0(\cos\th)\vec{\hat{r}}\right)$  
  $\textstyle \quad +$ $\displaystyle u_0 K^{-1}\epsilon~\vec{\nabla}\!\times\!\vec{\nabla}\!\times\!\left(r^3(1-r^2)^2P_2^0(\cos\th)\vec{\hat{r}}\right)$ (14)

with $K^{-1}=\sqrt{9009/572}$ and $\epsilon=0.5\sqrt{143/1008}$ such that the rms poloidal to toroidal energy ratio is 0.5, and the flow has an rms value of u0. Plm and $\vec{\hat{r}}$ are defined in Appendix A.

Like Livermore & Jackson (2005) we consider the axisymmetric (m=0) and equatorially antisymmetric magnetic field solution for a unit sphere (r0=1) embedded in a vacuum. Table 2 shows the convergence in the eigenvalue with the largest real part as a function of truncation $n_{\rm max}$ and $l_{\rm max}$ for two magnetic Reynolds numbers $R_{\rm m} = u_0r_0/\eta=10$ and $R_{\rm m}=100$, together with the converged values given by Livermore & Jackson (2005) (LJ05). There is a difference of about one permille between their value and ours for $R_{\rm m} = 10$.

3.3 $\alpha ^2$and $\alpha ^2\Omega $-dynamos

We reproduced the critical dynamo numbers $R_\alpha$ further for the dipolar (l=1) mode of an isotropic $\alpha ^2$-dynamo with $\alpha_{rr}=\alpha_{\th\th}= \alpha_{\varphi\varphi}=R_\alpha\sin N\pi(r-r_i)$ and N=1,2 in a spherical shell of inner and outer radius ri and r0 with r0-ri=1 and ri/r0=0.35 and 0.8 surrounded by a vacuum and either an insulating or a conducting inner core, as reported in Table 2 of Schubert & Zhang (2001). With this test we treated in particular two different aspect ratios of a spherical shell (a thick and a thin one) and two different molecular diffusivities (insulating or conducting) of the inner core.

Finally we applied our method to an $\alpha ^2\Omega $-dynamo of Jiang & Wang (2006) who employ the classical eigenvalue treatment for the poloidal and toroidal scalars P and T expanded in spherical harmonics in the angular coordinates and in Chebychev polynomials in r-direction. We set $u_r=u_\th=0$, $u_\varphi= R_\Omega r^3\sin^3\th$, $\alpha_{ij}=\beta_{ijk}=0$, except $\alpha_{rr}=\alpha_{\th\th}=\alpha_{\varphi\varphi}=R_\alpha\sin2\pi(r-r_i)/(r_0-r_i)\cos\th$ with ri=0.5, r0=1, embedded in a vacuum inside and outside, and $R_\alpha=\alpha_0r_0/\eta=10$ and $R_\Omega=u_0r_0/\eta=1000$. Some results obtained by the new and the classical methods are compiled in Table 3. Numbers in parentheses $(\dots~,~\dots)$ are the real and imaginary parts of complex eigenvalues. The real part denotes the growth rate, the imaginary part the frequency of the mode in units of $\eta/r_0^2$. The modes are axisymmetric (m=0), the fundamental mode is monotonously growing and symmetric (indicated by S) with respect to the equator, and the fourth overtone is damped, oscillatory and antisymmetric (indicated by A). Modes with higher m are more strongly damped. $n_{\rm max}$ refers to the maximum radial number of the decay modes (spherical Bessel functions) for the new method and to the maximum degree of the Chebychev polynomials for the code of Jiang & Wang (2006), respectively. Since the modes have smaller length scales in latitudinal than in radial direction, higher values of l than of n are required for convergence. We find remarkably similar convergence of the eigenvalues for both methods. This also applies to modes with higher m. Of course we have also verified that the eigenfunctions obtained with the two methods are identical.

4 Geodynamo models

Having proven that the new method works correctly and efficiently, we now apply it to determine the eigensolutions of the dynamo operator with mean-field coefficients obtained from self-consistent numerical simulations of the geodynamo. For a recent review of numerical geodynamo simulations, see Christensen & Wicht (2007). Schrinner et al. (2007) developed an efficient method of calculating all tensorial mean-field coefficients  $\vec{\alpha}$ and $\vec{\beta}$ and compared the results of mean-field and direct numerical simulations of the geodynamo. We plan to use the eigenmodes of the dynamo equation to decompose the magnetic field of the numerical simulations and to determine the statistical properties of the mode coefficients (Hoyng 2009) to analyse the working of the geodynamo.

4.1 Benchmark dynamo

We examine a quasi-steady geodynamo model which has been used before as a numerical benchmark dynamo (Christensen et al. 2001, case 1). The governing parameters are Ekman number E=10-3, Rayleigh number Ra = 100, Prandtl number Pr=1, and magnetic Prandtl number Pm = 5. The convection pattern is columnar with a natural 4-fold azimuthal symmetry and is stationary except for an azimuthal drift. The intensity of the fluid motion is characterised by a magnetic Reynolds number of $R_{\rm m} \simeq 40$, defined with a characteristic flow velocity, the thickness of the convecting shell, and the molecular magnetic diffusivity. The magnetic energy density exceeds the kinetic one by a factor of 20.

In Schrinner et al. (2007), the mean-field coefficients are derived from the numerical simulation. We solved the dynamo equation with these mean-field coefficients by the new method and obtained the eigenvalues and eigenfunctions. Since the coefficients are spatially variable to a considerable degree, converged solutions require high truncation levels in n and l. The eigenvalues of the first two modes are shown in Table 4. Beyond $n_{\rm max} \simeq 20$ and $l_{\rm max} \simeq 20$, the eigensolution of the first mode does not change significantly and is displayed in Fig. 1. The convergence of the second mode requires a larger $l_{\rm max}$ of about 32. The results for high values of $n_{\rm max}$ may be affected by the spatial variation of the mean-field coefficients and would require more than 66 radial quadrature points to compute the matrix elements.

A comparison of Fig. 1 with its counterpart Fig. 10 of Schrinner et al. (2007) shows that the field of the antisymmetric fundamental mode resembles the field of an initial-value mean-field dynamo calculation remarkably well as it does the axisymmetric component of the direct numerical simulation. The mode here grows slightly with a rate around $\lambda_0 \simeq 4.2\eta/L^2$, the field of the initial value calculation decays slightly with a rate of approximately - $0.25\eta/L^2$, while the solution of the direct numerical simulation is stationary[*]. Here L=r0-ri=1 is the thickness of the spherical shell. The difference in these rates between the eigenvalue and initial value calculation comes from the higher numerical diffusivity of the latter at the chosen resolution of 33 radial and 80 latitudinal grid points. The difference is actually small, much less than one effective decay rate, because the relevant turbulent diffusivity, described by the $\vec{\beta}$ coefficient with values up to $33\eta$, is much higher than the molecular one.

Besides the true physical eigenmodes, we find growing unphysical spurious eigenmodes. Their eigenvalues depend strongly on the resolution, and their eigenfunctions are highly structured. We attribute their appearance to a locally confined inappropriate parametrisation of the mean electromotive force by the mean-field coefficients $\vec{\alpha}$ and $\vec{\beta}$ (Schrinner et al. 2007). The spurious modes are present neither in the initial value calculation nor in the following example of a time-dependent dynamo, because of a higher numerical and molecular diffusivity, respectively.

4.2 A time-dependent dynamo in the columnar regime

The next example has stronger forcing with parameters E=10-4, Ra=334, Pr=1, and Pm=2. The numerical simulation by Olson et al. (1999, case 2) shows a highly time-dependent, but still dominantly columnar convection characterised by a magnetic Reynolds number of $R_{\rm m} \simeq 88$. The magnetic energy exceeds the kinetic energy by a factor of three. The magnetic field has a strong axial dipole contribution. Although chaotically time-dependent, the velocity field is symmetric and the magnetic field antisymmetric with respect to the equatorial plane.

Table 4:   Eigenvalues in units of $\eta /L^2$ of the first two eigenmodes of the benchmark dynamo.

\begin{figure} \par\includegraphics[width=8.8cm,clip]{13702fg1.eps} \end{figure} Figure 1:

Magnetic field structure of the fundamental antisymmetric eigenmode for the benchmark dynamo. Compare with Fig. 10 of Schrinner et al. (2007). For each plot the grey scale is separately adjusted to its maximum modulus with white as negative and black as positive. The contour lines correspond to $\pm $0.1, $\pm $0.3, $\pm $0.5, $\pm $0.7, and $\pm $0.9 of the maximum modulus.

Open with DEXTER

The mean-field coefficients are obtained as before by the test-field method of Schrinner et al. (2007). The coefficients are now of course also highly time-dependent. A time average yields coefficients that roughly resemble those for the benchmark dynamo, although there are differences in some profiles and amplitudes.

For the time-averaged dynamo operator the eigenvalues of the first two antisymmetric eigenmodes for various values of $n_{\rm max}$ and $l_{\rm max}$ are shown in Table 5. It seems that a value of $n_{\rm max}\simeq16$ is sufficient for convergence, while $l_{\rm max}\simeq32$ is needed. Figure 2 shows the eigenfunctions of these modes. The eigensolutions for $\lambda_1=(-6.298, 0.0)$and $\lambda_2= (-28.712, \pm5.364)$, values for $n_{\rm max}=16$ and $l_{\rm max}=32$, are symmetric with respect to the equator.

Table 5:   Eigenvalues in units of $\eta /L^2$ of the first two antisymmetric eigenmodes of the temporally averaged dynamo operator obtained from the time-dependent dynamo (case 2, Sect. 4.2).

\begin{figure} \par\resizebox{9cm}{!}{\includegraphics{13702fg2.eps}} \end{figure} Figure 2:

Radial component of the first two antisymmetric eigenmodes for the case2 dynamo. Left: fundamental mode; middle: real part of the first overtone; right: imaginary part of the first overtone. Grey scales and contours as in Fig. 1.

Open with DEXTER

An initial-value, mean-field dynamo calculation with the same mean velocity and dynamo coefficients shows a slighly decaying solution with a decay rate of approximately $-5.9\eta/L^2$ which is to be compared with the eigenvalue $\lambda_0\simeq-3.87\eta/L^2$ of the fundamental mode. Again, the turbulent diffusivity exceeds the molecular one by a factor of up to 23 in this case. The difference in the decay rates is therefore much less than one effective decay rate. As for the benchmark dynamo, the profile of the antisymmetric fundamental mode is again remarkably similar to the solution of the initial value calculation and to the axisymmetric component of the direct numerical simulation.

A decomposition of the actual magnetic field of the simulation by Olson et al. (1999, case 2) in eigenfunctions of the time-averaged dynamo operator, i.e., $\vec{B}(\vec{r},t)=\sum_i c_i(t)\vec{B}_i(\vec{r})$, shows that the antisymmetric fundamental mode contributes to about 75 percent and, together with the first antisymmetric overtone (see Table 5 and Fig. 2), to about 85 percent of the total magnetic energy. The variability in time of the magnetic field of the direct numerical simulation is reflected in the variability of the expansion coefficients. More details are presented in Schrinner et al. (2010).

5 Conclusions and outlook

We presented a new method for computing the eigenvalues and eigenfunctions of the induction and the dynamo equation. The method is based on the biorthogonality of the adjoint electric current and the vector potential with an inner product defined by a volume integral over the fluid domain. The advantage of the method is that the velocity and dynamo coefficients do not have to be differentiated. The method is therefore well-suited for spatially strongly variable dynamo coefficients.

We tested the new method against the classical treatment and proved that it works correctly and efficiently. We applied it to two cases with dynamo coefficients derived from direct numerical simulations of the geodynamo. The obtained dynamo eigenmodes are promising candidates for decomposing the magnetic field of the numerical simulations and for analysing the statistical properties of the mode coefficients as proposed by Hoyng (2009).

Acknowledgements
We thank Ulrich Christensen, Johannes Wicht, and Robert Cameron for many useful discussions and support. We further thank the referee, Matthias Rheinhardt, for his detailed comments that helped to improve the paper.

Appendix A: Free magnetic decay modes in a sphere or spherical shell embedded in vacuum

We decompose the magnetic field in its poloidal and toroidal components

\begin{displaymath}\vec{B}=\vec{\nabla}\!\times\!\vec{\nabla}\!\times\!P\vec{\hat{r}}+\vec{\nabla}\!\times\!T\vec{\hat{r}}\end{displaymath} (A.1)

with defining scalars $P(r,\th,\varphi)$ and $T(r,\th,\varphi)$ and unit vector in radial direction $\vec{\hat{r}}=(1,0,0)$ in spherical coordinates $(r,\th,\varphi)$. The equation of free magnetic decay

\begin{displaymath}\frac{\partial\vec{B}}{\partial t}=-\eta\vec{\nabla}\!\times\!\vec{\nabla}\!\times\!\vec{B}\end{displaymath} (A.2)

with constant magnetic diffusivity $\eta$ then reads

\begin{displaymath}\frac{\partial}{\partial t}(P,T)=\eta\left(\Delta_{\rm H} +\frac{\partial^2}{\partial r^2}\right)(P,T) \end{displaymath} (A.3)

where $\Delta_{\rm H}$ is the horizontal Laplacian

\begin{displaymath}\Delta_{\rm H} =\frac{1}{r^2\sin^2\th}~\frac{\partial}{\parti... ...ac{1}{r^2\sin^2\th}~\frac{\partial^2}{\partial\varphi^2} \cdot \end{displaymath} (A.4)

The solutions are the free magnetic decay modes
$\displaystyle P_{lmn} = \exp\left(\lambda_{ln}^Pt\right)xf_l(p_{ln}x)~Y_l^m(\th,\varphi),$     (A.5)
$\displaystyle T_{lmn} = \exp\left(\lambda_{ln}^Tt\right)xg_l(t_{ln}x)~Y_l^m(\th,\varphi)$     (A.6)

with x=r/r0 where r0 is the radius of the sphere. The growth rates are given by

\begin{displaymath}\lambda_{ln}^P=-\eta p_{ln}^2/r_0^2 \;\;{\rm and}\;\; \lambda_{ln}^{\rm T} = -\eta t_{ln}^2/r_0^2 \end{displaymath} (A.7)

and are independent of the azimuthal degree m. The constants pln and tln are

\begin{displaymath}p_{ln}=j_{l-1,n} \;\; {\rm and} \;\; t_{ln}=j_{l,n} \end{displaymath} (A.8)

where jl,n is the nth zero of jl. The Ylm are the spherical harmonics and normalised to unity by taking

\begin{displaymath}Y_l^m(\th,\varphi)=\left(\frac{4\pi}{2l+1}\frac{(l+m)!}{(l-m)!}\right)^{-1/2} P_l^m(\cos\th)~{\rm e}^{im\varphi} \end{displaymath} (A.9)

using Ferrer's definition of the Legendre functions of first kind Plm with degree l and order m.

For a sphere embedded in vacuum the radial functions are given by

$\displaystyle f_l(p_{ln}x) = \left\{ \begin{array}{ll} a_{ln}~j_l(p_{ln}x) & 0\le x\le 1 \\ [0.5mm] a_{ln}~j_l(p_{ln})x^{-l} & x\ge 1 \end{array} \right.$     (A.10)
$\displaystyle g_l(t_{ln}x) = \left\{ \begin{array}{ll} b_{ln}~j_l(t_{ln}x) & 0\le x\le 1 \\ [0.5mm] 0 & x\ge 1 \end{array} \right.$     (A.11)

with the spherical Bessel functions of first kind jl. This ensures regularity in the origin of the sphere, vanishing toroidal component at its outer boundary and smooth transition of the poloidal component to a potential field in the vacuum outside.

For a spherical shell with inner radius ri (xi=ri/r0) and outer radius r0 (x0=1) embedded in vacuum the radial functions inside the shell are given by

fl(plnx)=jl(plnx)-yl(plnx)jl+1(plnxi)/yl+1(plnxi) (A.12)

and

gl(tlnx)=jl(tlnx)-yl(tlnx)jl(tlnxi)/yl(tlnxi), (A.13)

and the constants in the arguments are the roots of

jl+1(plnxi)yl-1(pln)-jl-1(pln)yl+1(plnxi)=0 (A.14)

for pln and of

jl(tln)yl(tlnxi)-jl(tlnxi)yl(tln)=0 (A.15)

for tln. Here yl are the spherical Bessel functions of second kind.

The magnetic field of the decay modes $\vec{B}_i$ is obtained by inserting the spatial parts of the defining scalars Plmn and Tlmn, respectively, into Eq. (A.1). Here we have comprised the three indices into one. The decay modes are self-adjoint on V+E, so that the adjoint functions are obtained simply by complex conjugation: $\vec{\hat{B}}_k=\vec{B}_k^*$ and likewise $\vec{\hat{J}}_k=\vec{J}_k^*$. Normalisation on V+E, i.e., $(\vec{B}_k^*,\vec{B}_i)_{V+E}=(\vec{J}_k^*,\vec{A}_i)_V=\delta_{ki}$, is thus straightforward. For a unit sphere the radial functions are normalised to unity by scaling the fl with

\begin{displaymath}a_{ln}=\left(\frac{1}{2r_0}l(l+1)j_{l-1,n}^2~j_l^2(j_{l-1,n})\right)^{-1/2}\end{displaymath} (A.16)

and the gl with

\begin{displaymath}b_{ln}=\left(\frac{r_0}{2}l(l+1)j_{l+1}^2(j_{l,n})\right)^{-1/2}\cdot \end{displaymath} (A.17)

For a spherical shell the normalisation constants are more lengthy expressions, which we suppress here.

The free magnetic decay modes form a complete and orthogonal set of functions, and they obey the boundary conditions of the magnetic field between the dynamo volume V and the exterior vacuum E.

We mention for completeness that the poloidal decay modes are not self-adjoint on V, i.e., $(\vec{B}_k^*,\vec{B}_i)_V \ne \delta_{ki}$. If we like to work with an inner product defined on V, the adjoint functions $\vec{\hat{B}}_k$ can be constructed by requiring $(\vec{\hat{B}}_k,\vec{B}_i)_V=\delta_{ki}$, similar to the one described in Sect. 2.5.

References

  1. Brandenburg, A., Rädler, K.-H., Rheinhardt, M., & Subramanian, K. 2008, ApJ, 687, L49 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bullard, E., & Gellman, H. 1954, Phil. Trans. R. Soc. Lond. A, 247, 213 [Google Scholar]
  3. Christensen, U. R., & Wicht, J. 2007, in Treatise of Geophysics, Core Dynamics, ed. G. Schubert, 8, 245 (Amsterdam: Elsevier) [Google Scholar]
  4. Christensen, U. R., Aubert, J., Cardin, P., et al. 2001, Phys. Earth Planet. Inter., 128, 25 [Google Scholar]
  5. Deinzer, W., & Stix, M. 1971, A&A, 12, 11 [NASA ADS] [Google Scholar]
  6. Deinzer, W., Grosser, H., & Schmitt, D. 1993, A&A, 273, 405 [NASA ADS] [Google Scholar]
  7. Dudley, M., & James, R. 1989, Proc. R. Soc. Lond. A, 425, 407 [NASA ADS] [CrossRef] [Google Scholar]
  8. Fuchs, H., Rädler, K.-H., & Schüler, M. 1993, IAU Symp., 157, 129 [NASA ADS] [Google Scholar]
  9. Gubbins, D. 1973, Phil. Trans. R. Soc. Lond. A, 274, 493 [NASA ADS] [CrossRef] [Google Scholar]
  10. Gubbins, D., Barber, C., Gibbons, S., & Love, J. 2000, Proc. R. Soc. Lond. A, 456, 1333 [NASA ADS] [CrossRef] [Google Scholar]
  11. Hoyng, P. 1988, ApJ, 332, 857 [NASA ADS] [CrossRef] [Google Scholar]
  12. Hoyng, P. 2009, Phys. Rev. E, 79, 046320 [NASA ADS] [CrossRef] [Google Scholar]
  13. Hoyng, P., & van Geffen, J. H. G. M. 1993, Geophys. Astrophys. Fluid Dynamics, 68, 203 [NASA ADS] [CrossRef] [Google Scholar]
  14. Hoyng, P., & Schutgens, N. A. J. 1995, A&A, 293, 777 [NASA ADS] [Google Scholar]
  15. Jiang, J., & Wang, J. X. 2006, Chin. J. Astron. Astrophys., 6, 227 [NASA ADS] [CrossRef] [Google Scholar]
  16. Jiang, J., & Wang, J. X. 2007, MNRAS, 377, 711 [NASA ADS] [CrossRef] [Google Scholar]
  17. Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2009, A&A, 500, 633 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Krause, F., & Rädler, K.-H. 1980, Mean-field magnetophydrodynamics and dynamo theory, (Oxford: Pergamon Press) [Google Scholar]
  19. Kumar, S., & Roberts, P. H. 1975, Proc. R. Soc. Lond. A, 344, 235 [NASA ADS] [CrossRef] [Google Scholar]
  20. Livermore, P. W., & Jackson, A. 2004, Proc. R. Soc. Lond. A, 460, 1453 [NASA ADS] [CrossRef] [Google Scholar]
  21. Livermore, P. W., & Jackson, A. 2005, Geophys. Astrophys. Fluid Dynamics, 99, 467 [Google Scholar]
  22. Moffatt, H. K. 1978, Magnetic field generation in electrically conducting fluids (Cambridge: Cambridge University Press) [Google Scholar]
  23. Olson, P., Christensen, U., & Glatzmaier, G. A. 1999, J. Geophys. Res., 104, 10383 [NASA ADS] [CrossRef] [Google Scholar]
  24. Ossendrijver, M., Stix, M., & Brandenburg, A. 2001, A&A, 376, 713 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Ossendrijver, M., Stix, M., Brandenburg, A., & Rüdiger, G. 2002, A&A, 394, 735 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Rädler, K.-H. 1980, Astron. Nachr., 301, 101 [Google Scholar]
  27. Rädler, K.-H., & Bräuer, H.-J. 1987, Astron. Nachr., 308, 101 [NASA ADS] [CrossRef] [Google Scholar]
  28. Rädler, K.-H., Rheinhardt, M., Apstein, E., & Fuchs, H. 2002, Magnetohydrodynamics, 38, 41 [NASA ADS] [EDP Sciences] [Google Scholar]
  29. Roberts, P. H. 1960, J. Math. Analysis Applic., 1, 195 [CrossRef] [Google Scholar]
  30. Roberts, P. H. 1972, Phil. Trans. R. Soc. Lond. A, 272, 663 [Google Scholar]
  31. Roberts, P. H., & Stix, M. 1972, A&A, 18, 453 [NASA ADS] [Google Scholar]
  32. Schmitt, D., & Schüssler, M. 1989, A&A, 223, 343 [NASA ADS] [Google Scholar]
  33. Schrinner, M., Rädler, K.-H., Schmitt, D., Rheinhardt, M., & Christensen, U. R. 2005, Astron. Nachr., 326, 245 [NASA ADS] [CrossRef] [Google Scholar]
  34. Schrinner, M., Rädler, K.-H., Schmitt, D., Rheinhardt, M., & Christensen, U. R. 2007, Geophys. Astrophys. Fluid Dynamics, 101, 81 [CrossRef] [Google Scholar]
  35. Schrinner, M., Schmitt, D., Cameron, R., & Hoyng, P. 2010, Geophys. J. Int., 182, 675 [NASA ADS] [CrossRef] [Google Scholar]
  36. Schubert, G., & Zhang, K. 2001, ApJ, 557, 930 [Google Scholar]
  37. Steenbeck, M., & Krause, F. 1969, Astron. Nachr., 291, 49 and 271 [NASA ADS] [CrossRef] [Google Scholar]

Footnotes

...$R_\alpha=4.493409458$[*]
The numerical value of $R_\alpha$ is equal to the first zero of the spherical Bessel function j1.
... stationary[*]
In Schrinner et al. (2007) the mean flow entered with a sign error into the initial value calculation, leading to a stronger decay of - $3.5\eta/L^2$. We apologise and correct this value here.
Copyright ESO 2010

All Tables

Table 1:   Eigenvalues of the fundamental dipolar mode and the fifth and tenth overtones of the $\alpha ^2$-sphere.

Table 2:   Eigenvalues of the fundamental mode for the $\vec{t_1s_2}$ flow for two magnetic Reynolds numbers $R_{\rm m}$.

Table 3:   Eigenvalues of the first and fifth modes of the $\alpha ^2\Omega $-dynamo (see Sect. 3.3).

Table 4:   Eigenvalues in units of $\eta /L^2$ of the first two eigenmodes of the benchmark dynamo.

Table 5:   Eigenvalues in units of $\eta /L^2$ of the first two antisymmetric eigenmodes of the temporally averaged dynamo operator obtained from the time-dependent dynamo (case 2, Sect. 4.2).

All Figures

  \begin{figure} \par\includegraphics[width=8.8cm,clip]{13702fg1.eps} \end{figure} Figure 1:

Magnetic field structure of the fundamental antisymmetric eigenmode for the benchmark dynamo. Compare with Fig. 10 of Schrinner et al. (2007). For each plot the grey scale is separately adjusted to its maximum modulus with white as negative and black as positive. The contour lines correspond to $\pm $0.1, $\pm $0.3, $\pm $0.5, $\pm $0.7, and $\pm $0.9 of the maximum modulus.

Open with DEXTER
In the text

  \begin{figure} \par\resizebox{9cm}{!}{\includegraphics{13702fg2.eps}} \end{figure} Figure 2:

Radial component of the first two antisymmetric eigenmodes for the case2 dynamo. Left: fundamental mode; middle: real part of the first overtone; right: imaginary part of the first overtone. Grey scales and contours as in Fig. 1.

Open with DEXTER
In the text

Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.