Computation of kinematic and magnetic $\alpha$-effect and eddy diffusivity tensors by Pad\'e approximation

TL;DR

This paper develops Padé approximation methods to compute the $oldsymbol{ ext{α}}$-effect and eddy diffusivity tensors in fluid flows, providing accurate estimates of magnetic field growth rates and flow stability in dynamo theory.

Contribution

It introduces Padé approximation techniques for tensors in kinematic dynamo problems, enabling better predictions of magnetic field growth and flow stability.

Findings

Padé approximants accurately estimate growth rates for small magnetic diffusivity.

Eddy diffusivity can tend to negative infinity near dynamo onset.

First 49 coefficients are insufficient to capture certain behaviors.

Abstract

We present examples of Pad\'e approximation of the $\alpha$-effect and eddy viscosity/diffusivity tensors in various flows. Expressions for the tensors derived in the framework of the standard multiscale formalism are employed. Algebraically the simplest case is that of a two-dimensional parity-invariant six-fold rotation-symmetric flow, where eddy viscosity is negative, indicating intervals of large-scale instability of the flow. Turning to the kinematic dynamo problem for three-dimensional flows of an incompressible fluid, we explore application of Pad\'e approximants for computation of tensors of magnetic $\alpha$-effect and, for parity-invariant flows, of magnetic eddy diffusivity. We construct Pad\'e approximants of the tensors expanded in power series in the inverse molecular diffusivity $1/\eta$ around $1/\eta=0$. This yields the values of the dominant growth rate due to the action of the $\alpha$-effect or eddy diffusivity to satisfactory accuracy for $\eta$, several dozen times smaller than the threshold, above which the power series is convergent. For one sample flow, we observe eddy diffusivity tending to negative infinity when $\eta$ tends from above to the point of the onset of small-scale dynamo action in a symmetry-invariant subspace where a neutral small-scale magnetic mode resides. However, 49 first coefficients in the power series in $1/\eta$ prove insufficient for Pad\'e approximants to reproduce this behaviour. We do computations in Fortran in the standard `double' (real*8) and extended `quadruple' (real*16) precision, as well as perform symbolic calculations in Mathematica.

Tables

Table 1. Table 1: First 39 non-zero coefficients of ( 6 ) for the DHF (prime decomposition, where presented). The first 5 significant figures of these coefficients are given in VGF . ` ​ ` ​ a ⋘ p ⋙ c ​ " very-much-less-than ` ` 𝑎 𝑝 very-much-greater-than 𝑐 " ``a\lll p\ggg c" denotes a natural number containing p 𝑝 p decimal digits between a 𝑎 a and b 𝑏 b .

$n$	Coefficient ${\nu }_E^{\left(n\right)}$ (exact rational number)
1	$\frac{3}{2^2}$
3	$-\frac{3\times 5\times 11\times 1931\times 80491}{2^9\times 7^4\times {13}^2\times {19}^2}$
5	$-\frac{3\times {53}^2\times 222967\times 1994517983033813651288306079222192539}{2^{19}\times 5^2\times 7^9\times {13}^6\times {19}^7\times {31}^3\times {37}^2\times {43}^3\times {61}^2}$
7	$\frac{3^3\times 23\times 17401\times 11608063\times 570396658307516795186040829874710499\times 595146062519802577066082838776447096016784218965582671080441286999}{2^{25}\times 5^2\times 7^{17}\times {13}^{10}\times {19}^{11}\times {31}^5\times {37}^7\times {43}^5\times {61}^6\times {67}^3\times {73}^3\times {79}^3\times {97}^3\times {103}^2\times {109}^2}$
9	$-\frac{9606359879⋘188⋙\mathrm{\hspace{0.25em}5777697637}}{2^{33}\times 5^2\times 7^{24}\times {11}^2\times {13}^{14}\times {19}^{15}\times {31}^7\times {37}^{11}\times {43}^7\times {61}^{10}\times {67}^5\times {73}^5\times {79}^5\times {97}^5\times {103}^3\times {109}^6\times {127}^3\times {139}^3\times {151}^3\times {157}^2\times {163}^2}$
11	$-\frac{7129561983⋘324⋙\mathrm{\hspace{0.25em}7108258721}}{8493879641⋘326⋙\mathrm{\hspace{0.25em}4312960000}}$
$⋮$	$⋮$
39	$-\frac{1648936106⋘9785⋙\mathrm{\hspace{0.25em}2091564067}}{3775138782⋘9788⋙\mathrm{\hspace{0.25em}0000000000}}$

Table 2. Table 2: Order parameter L 𝐿 L of the Padé approximants [ 2 ​ L − 1 / 2 ​ L ] delimited-[] 2 𝐿 1 2 𝐿 [2L-1/2L] (ratios of polynomials in 1 / η 1 𝜂 1/\eta ) constructed by the algorithm GGT for six independent entries of the symmetrized α 𝛼 \alpha -effect tensor 𝕬 s superscript 𝕬 s {{}^{\rm s}\boldsymbol{\mathfrak{A}}} .

tol$\sum$	${}_{}{}^s𝔄_1^1$	${}_{}{}^s𝔄_1^2$	${}_{}{}^s𝔄_1^3$	${}_{}{}^s𝔄_2^2$	${}_{}{}^s𝔄_2^3$	${}_{}{}^s𝔄_3^3$
${10}^{-{10}^{}}$	5	5	6	5	5	5
${10}^{-12}$	6	6	6	6	6	6
${10}^{-14}$	7	7	8	7	7	7
${10}^{-16}$	8	8	9	8	8	8
${10}^{-18}$	9	8	10	10	9	9
${10}^{-20}$	10	10	11	10	10	10

Equations

$$[0,L_1]\times [0,L_2]\ni \mathbf{x}=(x_1,x_2),L_1=\sqrt{3}{L}_2=2\pi .$$

$${\partial }_t{\nabla }^2\Psi +J({\nabla }^2\Psi ,\Psi )=\nu {\nabla }^2{\nabla }^2\Psi +\partial f_1/\partial x_2-\partial f_2/\partial x_1.$$

$${\nu }_{\mathrm{E}}(\Psi ,\nu )=\nu -⟨\left(Q+2\frac{\partial S}{\partial x_1}\right)\frac{\partial \Psi }{\partial x_2}⟩.$$

$$⟨g(x_1,x_2)⟩=\frac{1}{L_1{L}_2}{\int }_{x_1=0}^{L_1}{\int }_{x_2=0}^{L_2}g(x_1,x_2)\mathrm{d}{x}_2\mathrm{d}{x}_1,$$

$$\mathfrak{G}Q=\frac{\partial {\nabla }^2\Psi }{\partial x_2},\mathfrak{G}S=Q\frac{\partial {\nabla }^2\Psi }{\partial x_2}+2J\left(\Psi ,\frac{\partial Q}{\partial x_1}\right)-{\nabla }^{2}Q\frac{\partial \Psi }{\partial x_2}+4\nu \frac{\partial {\nabla }^{2}Q}{\partial x_1},$$

$$\mathfrak{G}:\psi \mapsto J({\nabla }^2\psi ,\Psi )+J({\nabla }^2\Psi ,\psi )-\nu {\nabla }^2{\nabla }^2\psi .$$

$$G(x_1,x_2)=\sum _{p,q\in \mathbb{Z}}{\hat{G}}_{p,q}\exp (2\pi \mathrm{i}(p{x}_1/L_1+q{x}_2/L_2)),$$

$${\nabla }^{-2}G=-(2\pi {)}^{-2}\sum _{p,q\in \mathbb{Z}\backslash \{0\}}\frac{{\hat{G}}_{p,q}}{(p/L_1{)}^2+(q/L_2{)}^2}\exp (2\pi \mathrm{i}(p{x}_1/L_1+q{x}_2/L_2)).$$

$$\begin{array}[]{rl}\Psi(x_{1},x_{2})=\big{[}\!\!\!&-\cos 2x_{1}-\cos(x_{1}+\sqrt{3}x_{2})-\cos(x_{1}-\sqrt{3}x_{2})\\ &+\cos(4x_{1}+2\sqrt{3}x_{2})+\cos(5x_{1}-\sqrt{3}x_{2})+\cos(x_{1}-3\sqrt{3}x_{2})\\ &-\cos 4x_{1}-\cos(2x_{1}+2\sqrt{3}x_{2})-\cos(2x_{1}-2\sqrt{3}x_{2})\\ &+\cos(4x_{1}-2\sqrt{3}x_{2})+\cos(5x_{1}+\sqrt{3}x_{2})+\cos(x_{1}+3\sqrt{3}x_{2})\big{]}/2.\end{array}$$

$${\nu }_{\mathrm{E}}(\Psi ,\nu )=\nu +\sum _{n=1}^{\infty }{\nu }_{\mathrm{E}}^{(n)}(\Psi ){\nu }^{-n}.$$

$$Q=\sum _{n=1}^{\infty }{Q}^{(n)}{\nu }^{-n},S=\sum _{n=1}^{\infty }{S}^{(n)}{\nu }^{-n}.$$

$${\nu }_{\mathrm{E}}^{(n)}(\Psi )=⟨-Q^{(n)}\frac{\partial \Psi }{\partial x_2}+2S^{(n)}\frac{{\partial }^2\Psi }{\partial x_1\partial x_2}⟩.$$

$$Q^{(n)}$$

$$S^{(n)}$$

$$\mathfrak{B}:f\mapsto {\nabla }^{-2}{\nabla }^{-2}\left[J({\nabla }^{2}f,\Psi )+J({\nabla }^2\Psi ,f)\right].$$

$$\text{L}\mathbf{b}$$

$$\text{L}:\mathbf{b}$$

$$\mathbf{b}$$

$$\lambda$$

$$⟨\mathbf{f}⟩=(2\pi {)}^{-3}{\int }_{{\mathbb{T}}^3}\mathbf{f}(\mathbf{X},\mathbf{x})\mathrm{d}\mathbf{x}=\sum _{k=1}^3{⟨\mathbf{f}⟩}_k{\mathbf{e}}_k,\{\mathbf{f}\}=\mathbf{f}-⟨\mathbf{f}⟩.$$

$$\text{L}{\mathbf{s}}_k=0,⟨{\mathbf{s}}_k⟩={\mathbf{e}}_k.$$

$$\text{L}{\mathbf{b}}_1+2\eta (\nabla \cdot {\nabla }_{\mathbf{X}}){\mathbf{b}}_0+{\nabla }_{\mathbf{X}}\times (\mathbf{v}\times {\mathbf{b}}_0)={\lambda }_1{\mathbf{b}}_0,$$

$${\nabla }_{\mathbf{X}}\times (\mathfrak{A}⟨{\mathbf{b}}_0⟩)={\lambda }_1⟨{\mathbf{b}}_0⟩,{\nabla }_{\mathbf{X}}\cdot ⟨{\mathbf{b}}_0⟩=0$$

$${\mathfrak{A}}_k=⟨\mathbf{v}\times {\mathbf{s}}_k⟩,$$

$$⟨{\mathbf{b}}_0⟩=\mathbf{B}{\mathrm{e}}^{\mathrm{i}\mathbf{q}\cdot \mathbf{X}},\mathbf{B}\cdot \mathbf{q}=0$$

$${\lambda }_{1_{\pm }}(\mathbf{q})=\frac{\mathrm{i}}{2}\left(({\mathfrak{A}}_3^2-{\mathfrak{A}}_2^3)q_1+({\mathfrak{A}}_1^3-{\mathfrak{A}}_3^1)q_2+({\mathfrak{A}}_2^1-{\mathfrak{A}}_1^2)q_3\right)\pm \sqrt{a},a=\mathbf{q}\cdot (\det {}^{\mathrm{s}}\mathfrak{A}){{}^{\mathrm{s}}\mathfrak{A}}^{-1}\mathbf{q}.$$

$$a={\alpha }_1{\alpha }_2(q_3^{\prime }{)}^2+{\alpha }_2{\alpha }_3(q_1^{\prime }{)}^2+{\alpha }_1{\alpha }_3(q_2^{\prime }{)}^2,$$

$${\gamma }_{\alpha }\equiv \underset{|\mathbf{q}|=1}{\max }\mathrm{Re}{\lambda }_{1_{\pm }}(\mathbf{q})=\sqrt{\max ({\alpha }_1{\alpha }_2,{\alpha }_2{\alpha }_3,{\alpha }_1{\alpha }_3)}$$

$${\nabla }^2{\mathbf{s}}_k=-{\eta }^{-1}\nabla \times (\mathbf{v}\times {\mathbf{s}}_k).$$

$${\mathbf{s}}_k=\sum _{n=0}^{\infty }{\mathbf{s}}_k^{(n)}{\eta }^{-n}$$

$${\mathbf{s}}_k^{(0)}={\mathbf{e}}_k,{\mathbf{s}}_k^{(n)}=-{\nabla }^{-2}\left(\nabla \times (\mathbf{v}\times {\mathbf{s}}_k^{(n-1)})\right)\text{mbox}forn\ge 1$$

$$\mathfrak{M}:\mathbf{s}\mapsto -{\nabla }^{-2}\left(\nabla \times (\mathbf{v}\times \mathbf{s})\right)$$

Full text

1 Introduction

Power series expansion of analytic functions is perhaps the most powerful tool of numerical analysis. Let us just note that most algorithms for numerical integration of ordinary differential equations, such as the Runge–Kutta methods, rely on Taylor series expansions for derivation. The truncated series — i.e., polynomials — are easy to compute, and thus provide an important basic algorithm for evaluation of many analytic functions.

However, there are two caveats. One is related to a finite precision of computations, stemming from the hardware architecture and employed in the overwhelming majority of computer codes. A well-known example of the resultant failure of a computational procedure is a straightforward attempt to compute by this technique an exponent of a real large negative number FMM . Mathematically, this does not present any difficulty — the large individual terms in the Taylor expansion around zero are guaranteed to mutually cancel out and yield the final result which is less than unity. For finite-precision computations, however, the cancellation is not any more guaranteed, and the initial growth of individual terms can result in ultimate loss of accuracy.

The other one stems from finiteness of the radius of convergence of most power series encountered in computational practice. A complementary technique is then needed to continue analytically a function defined by the power series outside its circle of convergence. This can be achieved by constructing the so-called Padé approximants Gi ; B90 ; B96 , i.e., an implementation of the continuation in the form of the ratio of two polynomials. Let us cite the words of appraisal in nr : “Padé approximation has the uncanny knack of picking the function you had in mind from among all the possibilities. Except when it doesn’t! That is the downside of Padé approximation: it is uncontrolled. There is, in general, no way to tell how accurate it is, or how far out in $x$ it can usefully be extended. It is a powerful, but in the end still mysterious, technique.”

A not less mysterious notion is that of eddy diffusivity St , also known as eddy (or turbulent) viscosity when fluid viscosity, the source of diffusion in hydrodynamics, is considered. Like the magnetic $\alpha$-effect and anisotropic kinetic alpha- (aka AKA) effect, eddy diffusivities are often encountered in magnetohydrodynamics when generation of large-scale magnetic fields by flows of electrically conductive fluids is considered. At first sight, it is in direct contradiction with the second principle of thermodynamics. Of course, in fact no physical laws are violated. Both notions just describe the mean influence of the small scales on large-scale structures.

According to the modern paradigm, cosmic magnetic fields (such as the solar or geomagnetic field) exist due to the dynamo processes in the moving electrically conductive medium (such as melted rocks in the outer Earth’s core) Mb . The generating flows are typically turbulent and feature a vast hierarchy of spatial and temporal scales. Small-scale fluctuations of the flow (called “cyclonic events” by E. Parker) give rise to small fluctuations of magnetic field. The interaction of the small-scale components of the magnetic field and flow velocity produces a mean electromotive force (e.m.f.) that may have a non-zero component parallel to the mean magnetic field, and this can be beneficial for magnetic field generation P55 . The part of the mean e.m.f. linear in the mean field gives rise to the so-called magnetic $\alpha$-effect. If the flow is parity-invariant, the $\alpha$-effect disappears and the impact of yet smaller spatial scales becomes apparent; the mean e.m.f. is then a linear combination of the first-order spatial derivatives of the mean field, giving rise to eddy (turbulent) diffusivity. These physical ideas are treated under various simplifying assumptions in the mean-field electrodynamics SKR ; KR , and, relying only on the first principles, by asymptotic methods of homogenization of elliptic operators in the magnetohydrodynamic multiscale stability theory GO1 ; GO2 ; V86 ; V87 ; DF ; La ; VZ . Weakly nonlinear stability problems are also amenable to these methods CZ .

Analysis of equations makes it evident that eddy viscosity/diffusivity acquires the unusual properties, because it acts on mean fields only, i.e., essentially an open physical system is considered. In fact, in this class of MHD systems the inverse energy cascade is important, energy proliferating from small scales towards large ones; the source of energy for the developing large-scale magnetic, hydrodynamic or combined MHD perturbation is the forcing applied to maintain the perturbed (also sometimes called basic) fluid flow.

Evaluation of the $\alpha$-effect and eddy diffusivity tensors involves solving the so-called auxiliary problems, which are linear problems for the respective elliptic operators of linearization. In the large-scale dynamo problem, computing the magnetic $\alpha$-effect tensor requires considering three such problems; the number increases to 12, when the tensor of eddy diffusivity is sought (unless auxiliary problems for the adjoint operator come into play, decreasing the number of auxiliary problems to be treated to just 6, see, e.g., ABNZ ; RCZ ). Interesting results (e.g., instability to large-scale perturbation or dynamos) are typically obtained for relatively small molecular viscosity or magnetic diffusivity. Consequently, high spatial resolution is needed when solving the auxiliary problems, which makes the problems computationally intensive. However, the tensors can be easily expanded in the respective Reynolds number (i.e., in the inverse viscosity or diffusivity provided the size of the flow periodicity box and the flow velocity are order unity), when it is small, i.e., for large viscosities and diffusivities. When the parameter tends to the critical value for the onset of the small-scale instability (i.e., in the dynamo context, to the value for which generation of small-scale magnetic field starts), the tensors usually exhibit singular behaviour ZPF ; RCZ ; ACZ ; ACZ2 of a simple pole type, bounding from above the radius of convergence of the series. This suggests to try Padé approximants for computing the tensors and the respective large-scale magnetic field / instability growth rates.

We report here numerical experiments exploring these ideas. The paper is organized as follows. In section 2 we apply the Padé approximants techniques for evaluation of the eddy viscosity in a two-dimensional flow with two symmetries, in whose presence the eddy viscosity tensor reduces to a scalar. In view of the first caveat discussed in the beginning of this introduction, we have chosen to perform the calculations in precise arithmetics allowing an arbitrary number of correct digits; for this purpose, we have used the programming language Mathematica, giving an opportunity to make symbolic computations. In section 3 we revert to the standard “double precision” computations (real8, in Fortran speak) of the magnetic $\alpha$-effect tensor, using the “quadruple precision” (real16) computations for comparison. In section 4 we again employ the symbolic capabilities of Mathematica to evaluate the magnetic eddy diffusivity tensor. Our findings are summarized in section 5.

2 Calculation of eddy viscosity

No truly two-dimensional flows exist in nature, but they mimick properties of natural objects, such as the atmosphere or ocean Tay ; Fb ; Lin . We analyze the eddy viscosity tensor, $\varepsilon_{ijk\ell}$, DF of a two-dimensional flow of incompressible fluid that has two symmetries: parity invariance ($S1$) and the six-fold rotation symmetry ($S2$).

Since an $S1$-symmetric flow has a center of symmetry, it cannot possess the large-scale anisotropic kinetic $\alpha$-effect FSS . This is important, because in the presence of the AKA effect the large-scale dynamics is essentially dispersive and non-diffusive, concealing the impact of the eddy viscosity. The symmetry $S2$ implies the isotropy of fourth-order tensors (see, for instance, LL ), in particular, $\varepsilon_{ijk\ell}=\nu_{\rm E}\delta_{ij}\delta_{k\ell}$, where the scalar $\nu_{\rm E}$ is the (standard) eddy viscosity and $\delta_{mn}$ is the Kronecker symbol. Although the assumption that a flow features the two symmetries is mathematically convenient, it may be non-realistic when considering natural or engineering problems CH .

Two-dimensional flows endowed with the symmetries $S1$ and $S2$ can be constructed as follows. A space-periodic flow is assumed, the periodicity cell being the rectangle

[TABLE]

Its stream-function $\Psi(t,{\bf x})$ is then a sum of Fourier modes, whose wave vectors are $p(2,0)+q(1,\sqrt{3})$, where $p$ and $q$ are integer. Any two such modes, that have wave vectors mutually related by rotations by $\pi/3$, must both be involved in the sum with the same real amplitude.

We begin this section by recalling the expression for the scalar eddy viscosity in terms of the solutions to two auxiliary problems and the analytical framework for evaluating the expansion of the eddy viscosity tensor in powers of the inverse of the molecular viscosity. We carry on by recalling the standard terminology and definitions of Padé approximants to a series. Finally, we discuss our results and conclusions.

2.1 Eddy viscosities and multiscale techniques

For parity-invariant and six-fold rotation-symmetric flows of incompressible fluid, eddy viscosities were calculated in DF by multiscale techniques. They can be expressed in terms of solutions to two auxiliary problems, which can be solved analytically only in special cases DF (e.g., if the flow depends only on a single spatial coordinate).

Briefly, the eddy viscosity in an isotropically forced two-dimensional flow is calculated as follows DF ; GVF . In terms of the stream-function $\Psi(t,{\bf x})$, the two-dimensional forced Navier–Stokes equation for incompressible fluid takes the form

[TABLE]

Its solution, $\Psi$, defines a basic flow. Here, $J(g_{1},g_{2})=(\partial g_{1}/\partial x_{1})(\partial g_{2}/\partial x_{2})-(\partial g_{1}/\partial x_{2})(\partial g_{2}/\partial x_{1})$ is the Jacobian determinant of functions $g_{1}(x_{1},x_{2})$ and $g_{2}(x_{1},x_{2})$, the operator $\nabla^{2}=\partial^{2}\!/\partial x^{2}_{1}+\partial^{2}\!/\partial x^{2}_{2}$ is the Laplacian, $\nu$ the kinematic molecular viscosity, and ${\bf f}=(f_{1},f_{2})$ the external force. (In our numerical examples, the flow norm and the size of the small-scale periodicity cell are order unity; thus, the inverse molecular viscosity can be regarded as the local Reynolds number, which is the key dimensionless parameter of the problem.) Now assume that the basic flow possesses the symmetries $S1$ and $S2$. Then the (scalar) eddy viscosity coefficient, $\nu_{\rm E}=\nu_{\rm E}(\Psi,\nu)$, that depends only on the basic flow and molecular viscosity, is GVF

[TABLE]

Here, the angle brackets denote the average over the periodicity cell:

[TABLE]

the scalar functions $Q(x_{1},x_{2})$ and $S(x_{1},x_{2})$ are solutions to two auxiliary problems

[TABLE]

and $\boldsymbol{\mathfrak{G}}$ denotes the linearization of the Navier–Stokes equation around $\Psi$,

[TABLE]

We restrict it to zero-mean functions of the same space periodicity as the basic flow. In this functional space we can define the inverse Laplace operator, that we denote $\nabla^{-2}$. A field $G$ from this space can be expressed as a Fourier series

[TABLE]

where $\left\langle G\right\rangle=\widehat{G}_{0,0}=0$. The equation $\nabla^{2}F=G\ \Leftrightarrow\ F=\nabla^{-2}G$ can now be readily solved: given $\left\langle F\right\rangle=0$, we find

[TABLE]

Existence of a deterministic time-independent space-periodic flow, which has an isotropic negative eddy viscosity when the molecular viscosity is below a critical value, was established in VGF . The so-called decorated hexagonal flow (DHF), on which we will also focus here, is

[TABLE]

The phenomenon of negative eddy viscosity is quite common among two-dimensional divergenceless space-periodic basic flows with the symmetries $S1$ and $S2$: about 1/3 of such flows feature negative eddy viscosity for sufficiently low molecular viscosity GVF . The auxiliary problems (2) can be solved either numerically by spectral methods, or by expanding in powers of $\nu^{-1}$ to high orders and afterwards extending analytically (relying on their meromorphy GVF ) beyond the disk of convergence. We examine the latter method enabling us to perform all calculations exactly.

2.2 Eddy viscosity expansion in powers of $\nu^{-1}$

Let us expand (1):

[TABLE]

To calculate $\nu_{\rm E}^{(n)}$, we expand the solutions $Q$ and $S$ to (2) in Maclaurin series in $\nu^{-1}$:

[TABLE]

Substituting the series into (1) and integrating by parts yields

[TABLE]

Here $Q^{(1)}=-\nabla^{-2}(\partial\Psi/\partial x_{2})$ and $S^{(1)}=-4\nabla^{-2}(\partial Q^{(1)}/\partial x_{1})$, and the subsequent terms satisfy the recurrence relations

[TABLE]

where the operator $\boldsymbol{\mathfrak{B}}$ is defined as

[TABLE]

Since we consider here parity-invariant flows, their stream-functions being even, i.e., $\Psi(-{\bf x})=\Psi({\bf x})$, the series (6) involves only odd powers of $\nu^{-1}$ GVF , i.e., $\nu_{\rm E}^{(n)}=0$ for all even $n$.

By definition, the $[L/M]$ Padé approximant to a series, whose first $m\geq L+M+1$ terms are known, is the ratio of a polynomial of degree $\leq L$ to a polynomial of degree $\leq M$, such that the first $L+M+1$ terms of the expansion of the ratio coincide with the respective terms of the series.

We use Padé approximants to reconstruct the dependence of the eddy viscosity on the inverse molecular viscosity employing the expansion (6). The Padé approximation techniques can also be naturally applied for exploring the poles of the eddy viscosity. A pole on the real axis can usually be linked to the onset of linear instability to small-scale perturbations, or sometimes (if it appers again on decreasing $\nu>0$) to its cessation.

2.3 Results of calculations

The performance of present-day computers and efficiency of symbolic programming software gives an opportunity to calculate the coefficients (7) using recurrence relations (8), and construct the Padé approximants exactly. All calculations reported in this section are performed with full precision by Mathematica HMM . Table 1 shows some of the first twenty non-zero coefficients $\nu_{\rm E}^{(n)}$ ($n=1,3,5,...,39$) for the DHF. One of the reasons to perform full precision symbolic computations by Mathematica has been a hope to discover useful relations between the coefficients of the approximants. Unfortunately, none are visible in the data of Table 1.

We truncate the series (6) at orders up to 39, even terms missing. Roots of their Padé approximants quickly stabilize near $\nu=\nu^{\star}\approx 0.58$ (see Fig. 1), indicating a transition to negative eddy viscosity at lower molecular viscosities. The root $\nu^{\star}$ is simple; a sharpened estimate is $1/\nu^{\star}=1.72144\pm 10^{-5}$.

These results indicate that Padé approximants is a reliable alternative to other methods for calculation of the point of the onset of large-scale instability in this type of flows. Not only they furnish stable estimates (provided the approximated function is meromorphic and the series is long enough), but also serve as precursors to future work. An illustration is Fig. 1 (right), where the $[14/14]$ approximant (as many others) exhibits a singular behaviour near $1/\nu\approx 2.81$ apparently related to the onset of linear instability to small-scale perturbations. The non-monotonicity of the eddy viscosity as a function of $\nu$ can be regarded as a manifestation of the complexity of the two-dimensional turbulent flow.

3 Computation of the magnetic $\alpha$-effect tensor

As we have seen in the previous section, the use of Padé expansions for reconstructing the dependence of eddy viscosity on the molecular one is possible and does not require very high degrees of the polynomials involved at least for moderate (not very small) molecular viscosities. However, arbitrary-precision symbolic calculations cannot be regarded as a practical realization of the approximation algorithm. Here we construct Padé approximants for evaluation of the magnetic $\alpha$-effect tensor using arithmetics of floating point numbers of the conventional double (real8) and the extended quadruple (real16) precision. The problem now at hand is to construct approximations applicable for small magnetic molecular diffusivities.

3.1 The multiscale formalism revealing the magnetic $\alpha$-effect

We review here the multiscale expansions VZ describing the kinematic generation of large-scale magnetic field by small-scale zero-mean space-periodic steady flows. Our magnetic modes depend on two three-dimensional spatial variables, the fast, $\bf x$, and slow, ${\bf X}=\varepsilon\bf x$, one (the flow $\bf v$ depends exclusively on $\bf x$). A magnetic mode $\bf b$ is an eigenfield of the magnetic induction operator $\L$:

[TABLE]

Here $\eta$ is the magnetic molecular diffusivity and Re $\lambda$ the growth rate of the mode $\bf b$. Both the mode and flow are solenoidal.

The scale ratio $\varepsilon$ is a small parameter of the problem, in which the magnetic mode $\bf b$ and its growth rate are expanded:

[TABLE]

By substituting the expansions into (\theparentequation.1) and the solenoidality conditions, we derive a hierarchy of equations for the coefficients in (10).

Like in the previous section, we denote by angle brackets the mean over the periodicity cell ${\mathbb{T}}^{3}$ in the fast variables and by the braces the fluctuating part:

[TABLE]

Here ${\bf e}_{k}$ are unit Cartesian coordinate vectors.

The relevant solution to the first (order $\varepsilon^{0}$) equation $\L{\bf b}_{0}=\lambda_{0}{\bf b}_{0}$ is $\lambda_{0}=0$ and a linear combination ${\bf b}_{0}=\sum_{k=1}^{3}\left\langle{\bf b}_{0}\right\rangle_{k}{\bf s}_{k}$ of small-scale solenoidal neutral magnetic modes ${\bf s}_{k}({\bf x})$ that are solutions to the three so-called auxiliary problems of type I:

[TABLE]

Averaging the second (order $\varepsilon^{1}$) equation,

[TABLE]

we obtain an eigenvalue problem

[TABLE]

(the subscript $\bf X$ denotes differentiation in the slow variables). Here $\boldsymbol{\mathfrak{A}}$ is the tensor of magnetic $\alpha$-effect. The $k$th column of this $3\times 3$ matrix is

[TABLE]

in agreement with the Parker’s P55 concept of the interaction of fine structures of the flow, $\bf v$, and magnetic field, $\sum_{k=1}^{3}\left\langle{\bf b}_{0}\right\rangle_{k}\{{\bf s}_{k}\}$, giving rise to a mean e.m.f., $\boldsymbol{\mathfrak{A}}\left\langle{\bf b}_{0}\right\rangle$, linear in the large-scale magnetic field $\left\langle{\bf b}_{0}\right\rangle$. For space-periodic mean fields

[TABLE]

where $\bf q$ is an arbitrary unit vector, straightforward algebra RCZ yields solutions to the eigenvalue problem (13):

[TABLE]

Here ${{}^{\rm s}\mathfrak{A}}_{k}^{l}=(\mathfrak{A}_{k}^{l}+\mathfrak{A}_{l}^{k})/2$ are entries of the symmetrized $\alpha$-tensor ${{}^{\rm s}\boldsymbol{\mathfrak{A}}}=(\boldsymbol{\mathfrak{A}}+\boldsymbol{\mathfrak{A}}^{*})/2$.

For $a\leq 0$, the $\alpha$-effect just sustains harmonic oscillations of the mean magnetic field in the slow time $T_{1}=\varepsilon t$. When $a>0$, the slow-time growth rate Re$\,\lambda_{1}({\bf q})=\sqrt{a}$ of the large-scale magnetic mode depends only on the symmetrized tensor ${{}^{\rm s}\boldsymbol{\mathfrak{A}}}$, whose eigenvalues $\alpha_{i}$ are real. In the Cartesian coordinate system, whose axes coincide with eigenvectors of ${{}^{\rm s}\boldsymbol{\mathfrak{A}}}$, (16) takes the form

[TABLE]

where $q^{\prime}_{i}$ are components of $\bf q$ in this basis. Thus,

[TABLE]

is the maximum slow-time growth rate of large-scale magnetic modes generated by the $\alpha$-effect. While the entries of the $\alpha$-effect tensor, $\mathfrak{A}^{k}_{l}$, are smooth functions of $\eta$, the graph of $\gamma_{\alpha}$ (17) has cusps at the points $\eta$, where $\alpha_{1}<\alpha_{2}=0<\alpha_{3}$ RCZ (see, e.g., two such cusps in Fig. 2).

When $a\neq 0$ and the kernel of the magnetic induction operator $\L$ does not involve small-scale zero-mean modes (generically both conditions are satisfied), all terms in the expansions (10) can be determined from the hierarchy of equations obtained by substituting the series into the eigenvalue equation (\theparentequation.1). If the symmetrized tensor ${{}^{\rm s}\boldsymbol{\mathfrak{A}}}$ is positively or negatively defined (and if the spatial periodicity of the eigenfunction is compatible with that of the flow), then the series (10) are summable V87 for sufficiently small $\varepsilon$ and constitute an analytical in $\varepsilon$ eigensolution for the large-scale magnetic induction operator; a unique $\varepsilon$-parameterized branch of the eigenvalues (\theparentequation.2) originates from any simple eigenvalue $\lambda_{1}$ of the $\alpha$-effect operator.

3.2 Padé approximation

The modes ${\bf s}_{k}$ (and, consequently, elements of the magnetic $\alpha$-effect tensor) are functions of molecular eddy diffusivity $\eta$, meromorphic in this parameter. (By contrast, the slow-time growth rates of modes generated by the $\alpha$-effect are not, because the square root present in (16) gives rise to branch points.) A power series expansion of ${\bf s}_{k}$ in $\eta^{-1}$ for large $\eta$ can be constructed like in the hydrodynamic problem considered in the previous section. We divide (11) by $\eta$ to obtain

[TABLE]

Consequently, the coefficients in the expansion

[TABLE]

satisfy the recurrence relations

[TABLE]

(cf. (5.14)–(5.15) in GO1 ). Here $\nabla^{-2}$ denotes the inverse Laplacian in the fast variables acting in the functional space of zero-mean fields. (Actually, we consider the problem for a flow, whose r.m.s. velocity is unity; the size of the periodicity box also being order unity, the inverse molecular diffusivity $\eta^{-1}$ can be regarded, like in the hydrodynamic case, to be equal to the local magnetic Reynolds number, the dimensionless parameter characterizing the mathematical properties of the problem; (18) can thus be understood as an expansion in a small Reynolds number.) These recurrence relations can be used to compute the coefficients in (18) by pseudospectral methods. It must be noted that mathematically they are perfectly suitable for numerical work: indeed, the presence of the inverse Laplacian in the second relation (19) suggests, that on increasing the number $n$ of the coefficient ${\bf s}^{(n)}_{k}$ they become smoother and their energy spectrum decay is steeper. This is in sharp contrast, for instance, with the recurrence relations for the Lagrangian time-Taylor coefficients in the expansions of solutions to the Euler equation for incompressible fluid flow PZF .

It is straightforward to determine the radius of convergence of the series (18), $\rho$, regarded as a function of $1/\eta$. The recurrence relations (19) involve the operator

[TABLE]

acting in the functional space of solenoidal zero-mean space-periodic fields. Since it is compact, its spectrum consists of a countable set of eigenvalues $\mu_{i}$, tending to zero. Consequently, $\rho\geq 1/\max_{i}|\mu_{i}|$; generically the equality holds, but $\rho>1/\max_{i}|\mu_{i}|$, if the expansion of ${\bf s}^{(1)}_{k}$ in the basis of eigenfunctions of the operator $\boldsymbol{\mathfrak{M}}$ does not involve eigenfunctions associated with any eigenvalue $\mu_{i}$ such that $|\mu_{i}|=\max_{i}|\mu_{i}|$. Therefore, generically the radii of convergence of the series for all the three ${\bf s}_{k}$ are the same.

Clearly, the radius of convergence of the ensuing series for the $\alpha$-effect tensor (14),

[TABLE]

is not smaller than that of the series (18) for ${\bf s}_{k}$. Let us note a symmetry property of (21). Denote by the superscript minus objects pertinent to the reverse flow $-\bf v$:

[TABLE]

The $\alpha$-effect tensor $\boldsymbol{\mathfrak{A}}^{-}$ for the reverse flow $-\bf v$ is obtained from the tensor $\boldsymbol{\mathfrak{A}}$ for $\bf v$ by transposition, i.e. ${(\boldsymbol{\mathfrak{A}}^{-}_{k})}_{l}=(\boldsymbol{\mathfrak{A}}_{l})_{k}$ RCZ . Since the recurrence relations (19) are linear in $\bf v$,

[TABLE]

where $({\bf s}^{-}_{k})^{(n)}$ denote the coefficients in the expansion of ${\bf s}^{-}_{k}$ in power series (18). This implies

[TABLE]

Therefore, the coefficients in the series (21) are symmetric matrices for odd $n$, and antisymmetric ones for even $n$. In other words, the symmetrized $\alpha$-effect tensor ${{}^{\rm s}\boldsymbol{\mathfrak{A}}}$ involved in computation of the discriminant $a$ in (16) is expanded in odd powers of $\eta^{-1}$, and of the common imaginary part of $\lambda_{1_{\pm}}({\bf q})$ in even powers; the latter expansion is not needed when computing just the growth rates.

3.3 Numerical results

Here, we construct Padé approximants for the $\alpha$-effect tensor (21) for a sample solenoidal flow, and compare the maximum growth rate values $\gamma_{\alpha}$ (17) obtained for the approximated tensor to those computed directly at individual values of $\eta$ by spectral methods.

For this purpose, a sample solenoidal flow has been synthesized as a Fourier series with pseudo-random coefficients, corrected to make it solenoidal and zero-mean. It involves Fourier harmonics for wave numbers not exceeding 10. The coefficients are scaled so that the energy spectrum decays exponentially by 10 orders of magnitude and the r.m.s. velocity is unity.

Solutions to auxiliary problems have been computed by the code Zh employing standard pseudo-spectral methods. For $\eta>0.05$, the problem was preconditioned by the operator $(-\nabla^{2})^{-1/2}$, readily available in the Fourier space. The resolution of $64^{3}$ Fourier harmonics was used. Energy spectra of the neutral modes ${\bf s}_{k}$ decay for this flow by at least 9 orders of magnitude for the smallest considered $\eta=0.035$. The Lebesgue space $L_{2}$ norms of ${\bf s}^{(n)}_{k}$ from $n=0$ to $n=128$ decay by 38 orders of magnitude, and hence the power series (18) and (21) converge for $\eta\apprge 0.50$ .

3.3.1 Approximation by the algorithm GGT

We have tried two approaches for construction of Padé approximants of the entries of the $\alpha$-effect tensor. Here we discuss the results obtained by the algorithm proposed in GGT .

Padé approximant $[M/L]_{f}$ of a function $f(y)=\sum_{n=0}^{\infty}f^{(n)}y^{n}$ is the ratio of two polynomials of degrees $M$ (numerator) and $L$ (denominator), whose $M+L+1$ first Taylor expansion coefficients coincide with those of $f$. For the sake of argument, let us assume $M\geq L$. Then $L+1$ coefficients of the denominator $d(y)=\sum_{n=0}^{L+1}d_{n}y^{n}$ satisfy the linear system of equations of the form

[TABLE]

(finding the coefficients of the numerator upon solving (24) is straightforward, see Gi ; B96 for details).

In our case, the coefficients ${\bf s}^{(n)}$ are obtained by applying iteratively the operator $\boldsymbol{\mathfrak{M}}$ (20). It is compact, and its eigenvalues, except for a finite number of them, are below unity in absolute value. The respective spectral components decay during the iterations according to the power law and sooner or later reduce in magnitude below the accuracy of computations. Consequently, for large $L$ the Toeplitz matrix of size $L\times(L+1)$ in the l.h.s. of (24) becomes numerically degenerate, i.e., its rank effectively falls below $L$. To construct an approximant under such adverse numerical conditions, it was proposed in GGT to compute the singular value decomposition of this matrix and to regard its effective rank as equal to the number of singular values whose absolute value exceeds the given relative tolerance tol (i.e., is not smaller than ${\tt tol}\|(f^{(1)},f^{(2)},...,f^{(M+L-1)},f^{(M+L)})\|$, where $\|\cdot\|$ is the standard Lebesgue space $L_{2}$ norm), decreasing the degrees of the polynomials involved in the Padé approximation, $M$ and $L$, by the number of the “missing” dimensions. To counter noise due to rounding errors, ${\tt tol}=10^{-14}$ was often used in GGT .

Beyond poor spectral properties of the system of equations for the Padé coefficients, there exist two other reasons for amplification of the numerical noise originally due to round-off errors:

$\bullet$ Pseudospectral methods used in computation of space-periodic solutions ${\bf s}_{k}$ to the auxiliary problems (11) and their coefficients ${\bf s}^{(n)}_{k}$ (19) involve fast Fourier transforms. These algorithms are very efficient. However, they operate by computing various linear combinations of the Fourier coefficients. Typically, at least for moderate molecular diffusivities, the energy spectra of these fields decay fast. In a sum of a large coefficient with a small (in absolute values) one, a significant part of the accuracy of the smaller coefficient is lost.

$\bullet$ Insufficiency of the spatial resolution can result in significant numerical errors. We may note that while increasing the resolution improves solutions, it aggravates the FFT accuracy problems.

We have tried the algorithm GGT for a set of $\tt tol$ values ranging between $10^{-10}$ and $10^{-20}$ using the MATLAB procedure provided by the authors of GGT . We have computed 65 first coefficients of the power series expansions of the symmetrized $\alpha$-effect tensor entries $({{}^{\rm s}\boldsymbol{\mathfrak{A}}}_{k})^{(n)}_{l}$ (which involve only odd powers of $1/\eta$, see section 3.2) up to order $\eta^{-129}$ terms with the spatial resolution of $64^{3}$ Fourier harmonics. The MATLAB procedure has been requested to construct the [63/64] Padé approximants for each entry. The results are shown in Fig. 2. We observe that the approximations of the maximum growth rates $\gamma_{\alpha}$ are relatively accurate for ${\tt tol}=10^{-12}$ and $10^{-16}$ for $\eta\apprge 0.1$. This bound is roughly 5 times smaller than the minimal $\eta$ for which the power series for ${{}^{\rm s}\mathfrak{A}}^{l}_{k}$ converge. Table 2 sheds light on the reasons why the gain is unsatisfactory (for $\eta\geq 0.1$ spectral computations of the $\alpha$-effect growth rates for individual $\eta$’s are efficient): since the rank decreases, when small in absolute value singular values are discarded, the algorithm ends up with very moderate orders $[2L/2L-1]$.

The four remaining panels in Fig. 2 reveal the presence of the so-called Froissart doublets in the approximants of some entries. Froissart doublets is a factor of the form $(1/\eta-a)/(1/\eta-\tilde{a})$ in the approximant, where the two constants $a$ and $\tilde{a}$ are close but distinct. Such a factor implies a singular behaviour of the approximant for $1/\eta$ close to $\tilde{a}$, not altering much its behaviour at distances from $\tilde{a}$ significantly larger than $|a-\tilde{a}|$. Often such factors are artifacts emerging due to noise in the data. Almost vertical segments of the plots are signatures of the Froissart doublets (see Fig. 2). They extend to both positive and negative infinity in the graphs of the approximants, but since (17) are nonlinear functions of the tensor entries, the respective segments of graphs of $\gamma_{\alpha}$ may be bounded from below and/or above. (Their detection has proved unexpectedly difficult; in order to reliably show their full range in the vertical direction, we have plotted the approximation step $10^{-9}$ along the abscissa.) We thus see that although the algorithm GGT is supposed to be robust, it is prone to yield approximants involving Froissart doublets.

3.3.2 Approximation by the algorithm nr

Since we have failed to obtain satisfactory resuts with the use of the algorithm GGT , we have also tested the algorithm nr . Quoting from nr , although the equations for the coefficients of the approximant involve a matrix in the Toeplitz form, “experience shows that the equations are frequently close to singular, so that one should not solve them by the methods” relying on this form, “but rather by full LU decomposition. Additionally, it is a good idea to refine the solution by iterative improvement (routine mprove in §2.5)”. This is implemented in their pade procedure. We have used it with one alteration: the routine mprove stops when the discrepancy increases; instead, it has been allowed to make up to 1000 improvement iterations permitting the discrepancy to temporarily grow and storing the minimum-discrepancy solution obtained in the course of these iterations (however, it has often been forced to stop before the allowed number of iterations has been performed, the iterative process blowing up with an overflow).

Approximate maximum slow-time growth rates (17) of large-scale magnetic modes generated by the $\alpha$-effect have been computed for the same flow using again $[2L-1/2L]$ Padé approximants for the entries ${{}^{\rm s}\mathfrak{A}}_{l}^{k}(\eta)$. They are compared in Fig. 3 for increasing orders $L$ with the actual maximum growth rates obtained by direct computation of the fields ${\bf s}_{k}$ using spectral methods for individual molecular diffusivities $\eta$.

Four Padé approximants $[2L-1/2L]$ of the entries of the $\alpha$-effect tensor have been constructed for each considered $L$, using the resolution of $64^{3}$ or $512^{3}$ Fourier harmonics, and running our code with the double (real8) or extended quadruple (real16) precision of the floating-point number arithmetics. In computers built around Intel and compatible processors, the former is standard, and the latter is not supported by hardware, but is software-emulated; however, many compilers do not require modifying the Fortran source code to use it, all the floating-point data and computations can be readily promoted to the real*16 precision by using the appropriate compiler option such as -r16. Higher precision and resolution has been expected to improve the accuracy of the coefficients of the approximants, to augment the orders of the approximants beyond those produced by the algorithm GGT and to increase the $\eta$ interval, where the growth rate values determined for the approximated $\alpha$-effect tensor are close to the actual growth rates.

We show in Fig. 3 the resultant approximations of $\gamma_{\alpha}$ for $L=4$, 8, 14, 20 and 28 to 31. The four graphs for $L=4$ are visually indistinguishable and we show only one of them; the same holds true for $L=8$. For $L=14$ and 20, the plots of the approximated $\gamma_{\alpha}$, computed with the quadruple precision for the two spatial resolutions, also visually coincide (Fig. 3(d),(f) and Fig. 3(h),(j)), but this is wrong for the respective double precision approximations. For higher $L$, all four plots are visually distinct. The quadruple precision approximations are plagued much less by the Froissart doublets (and never involve multiple occurrences of the doublets) than the double precision ones. Three high-$L$ quadruple precision approximations are reasonably accurate: for $L=28$ and the $512^{3}$ resolution for $\eta\geq 0.055$ (Fig. 3(n)); for $L=29$ and the $64^{3}$ resolution for $\eta\geq 0.05$ (Fig. 3(p)); and for $L=30$ and the $64^{3}$ resolution for $\eta\geq 0.05$ (Fig. 3(t)). Thus, the left end of the interval of validity of Padé approximations has decreased roughly twice compared to that obtained by the algorithm GGT . All other quadruple precision approximations for $L\geq 29$ (Fig. 3(r), (v), (x) and (z)) can also give reasonable accuracy for $\eta\apprge 0.06\div 0.07$ upon removal of Froissart doublets from the affected approximants of ${{}^{\rm s}\mathfrak{A}}_{l}^{k}$.

These results suggest, that Padé approximants are useful for representing the functional dependence of the slow-time growth rates (17) of large-scale magnetic modes generated by the $\alpha$-effect for fairly low magnetic molecular diffusivities. However, for construction of Padé approximants, accurate enough to serve small $\eta$, the quadruple precision arithmetics must be used, and hence run times become comparable to those of direct computation of the growth rates at individual $\eta$’s (note that real16 computations are typically ten times slower than real8 ones). Therefore another strategy is perhaps also sensible: in computations for an individual $\eta$, to use relatively low-order not-very-precise Padé approximants for neutral modes ${\bf s}_{k}$ (constructed, for instance, at each point in space or for each Fourier harmonics within the employed resolution) as the initial data for further refinement by the usual iterative methods.

4 Computation of the magnetic eddy diffusivity tensor

As discussed in section 2, an important class are parity-invariant flows. This symmetry is compatible with the equations of fluid dynamics (the Navier–Stokes or Euler equations) provided the forcing has the same property. For such a flow, the domain of the operator of magnetic induction $\L$ splits into the subspaces of parity-invariant fields (such that ${\bf f}(-{\bf x})=-{\bf f}({\bf x})$) and of parity-antiinvariant ones (such that ${\bf f}(-{\bf x})={\bf f}({\bf x})$). Solutions to the auxiliary problems (11), ${\bf s}_{k}({\bf x})$, are therefore parity-antiinvariant. This implies $\boldsymbol{\mathfrak{A}}=0$ (see (14)), i.e., no $\alpha$-effect acts in such flows, and hence $\lambda_{1}=0$.

4.1 The multiscale formalism revealing the magnetic eddy diffusivity

By (12),

[TABLE]

where the small-scale zero-mean (non-solenoidal!) fields ${\bf g}_{mk}({\bf x})$ solve nine auxiliary problems of type II:

[TABLE]

For parity-invariant $\bf v$, ${\bf g}_{mk}({\bf x})$ are also parity-invariant; moreover, ${\bf b}_{n}$ are parity-antiinvariant for all even $n$ and parity-invariant for odd $n$ VZ in the expansion (\theparentequation.1), and no odd powers of $\varepsilon$ enter the series (\theparentequation.2) for the eigenvalue $\lambda$.

Averaging the third (order $\varepsilon^{2}$) equation in the hierarchy yields

[TABLE]

where

[TABLE]

is the so-called tensor of magnetic eddy diffusivity correction. Again assuming that the mean field $\left\langle{\bf b}_{0}\right\rangle$ is a Fourier harmonics (15), we find RCZ

[TABLE]

The minimum

[TABLE]

is called the minimum magnetic eddy diffusivity; when it is negative, the interaction of the fluctuating small-scale velocity and magnetic field is capable of generating large-scale magnetic fields. For this reason, this quantity is of prime interest in the large-scale dynamo theory.

Expression (27) can be transformed into

[TABLE]

where ${\bf Z}_{l}$ are zero-mean solutions to three auxiliary problems for the adjoint operator:

[TABLE]

and $\L^{*}:{\bf z}\mapsto\eta\nabla^{2}{\bf z}-\{{\bf v}\times(\nabla\times{\bf z})\}$ is the operator adjoint to $\L$ acting in the space of zero-mean space-periodic fields.

4.2 Padé approximation

Relation (27) suggests to construct expansions of the solutions to the auxiliary problems in the inverse molecular diffusivity, (18) and

[TABLE]

Dividing (25) by $\eta$ yields

[TABLE]

whereby

[TABLE]

Clearly, the flow being parity-invariant, all ${\bf s}_{k}^{(n)}$ are parity-antiinvariant and all ${\bf g}_{mk}^{(n)}$ parity-invariant. By (27) and (32),

[TABLE]

Thus, algorithm I consists of the following steps:

$\bullet$ find the fields ${\bf s}_{k}^{(n)}$ employing (19);

$\bullet$ find the fields ${\bf g}_{mk}^{(n)}$ employing (33);

$\bullet$ calculate the coefficients $\boldsymbol{\mathfrak{D}}_{mk}^{(n)}$ employing (34).

Expressions (34) reveal symmetry properties of coefficients in the eddy diffusivity tensor expansion (similar to those of the coefficients of the $\alpha$-effect tensor expansion). Eddy diffusivity tensor for the reverse flow $-\bf v$ is related to that of the flow $\bf v$ by the relations ${(\boldsymbol{\mathfrak{D}}^{-}_{ml})}_{k}=-\mathfrak{D}_{mk}^{l}$ ABNZ . By (33), ${\bf g}_{mk}^{(n)}=(-1)^{n}({\bf g}_{mk}^{-})^{(n)}$. Thus, identities, analogous to those used in the case of the $\alpha$-effect tensor, reveal that for each fixed $m$ the coefficients $(\mathfrak{D}_{mk}^{l})^{(n)}$ in the series (34) are symmetric $3\times 3$ matrices for even $n$, and antisymmetric ones for odd $n$. This implies that the symmetrized matrix ${{}^{\rm s}\bf D}$ (\theparentequation.3), determining the discriminant $d$ (\theparentequation.2), is expanded in even powers of $1/\eta$, and the antisymmetric one $\bf D-{{}^{\rm s}\bf D}$, determining the common part of $\lambda_{2_{\pm}}({\bf q})$ (\theparentequation.1), in odd powers.

It is simple to show that, like in the case of $\alpha$-effect dynamos, for a given flow $\bf v$ the radius of convergence of all the series (32) and (34) (regarded as functions of $1/\eta$) is generically equal to $1/\max_{i}|\mu_{i}|$, where $\mu_{i}$ are eigenvalues of the compact operator $\boldsymbol{\mathfrak{M}}$ (20); convergence of the series is guaranteed for $\eta>\max_{i}|\mu_{i}|$.

An alternative form of the eddy diffusivity tensor can be exploited. Comparing (11) and (31), we find

[TABLE]

where the superscript minus denotes objects pertinent to the reverse flow $-\bf v$ (see (22)). Using (35) to eliminate ${\bf s}_{k}$ in (30) yields ABNZ

[TABLE]

By (35), all ${\bf Z}_{l}$ are parity-invariant, and hence (31) is equivalent to

[TABLE]

which implies a power series expansion

[TABLE]

where the coefficients satisfy recurrence relations

[TABLE]

By linearity of (38) in $\bf v$, the coefficients of such an expansion for the reverse flow are linked:

[TABLE]

This implies algorithm II for calculation of $(\mathfrak{D}^{l}_{mk})^{(n)}$ based on (36):

$\bullet$ find the fields ${\bf Z}_{l}^{(n)}$ applying (38);

$\bullet$ for $n\geq 1$ calculate

[TABLE]

Using (35), it is straightforward albeit tedious to transform (36) into

[TABLE]

Relations (35) imply

[TABLE]

Thus the coefficients in the expansion (34) have the entries

[TABLE]

Algorithm III consists of the following steps:

$\bullet$ determine coefficients in the expansion of the neutral magnetic modes ${\bf s}_{k}$ applying recurrence relations (19);

$\bullet$ in the course of these calculations, determine coefficients in the expansion of $\nabla\cdot{\bf Z}_{l}$ using (41);

$\bullet$ calculate the coefficients $(\mathfrak{D}^{l}_{mk})^{(n)}$ applying (42).

If the flow $\bf v$ involves a small number of Fourier harmonics, it is unclear a priori, which of the three algorithms is more efficient. Algorithm I involves computation of two sets of coefficients, for ${\bf s}_{k}$ and ${\bf g}_{mk}$; algorithm II only of the set of coefficients for ${\bf Z}_{l}$; algorithm III is an intermediate case involving two sets of coefficients, for ${\bf s}_{k}$ and $\nabla\cdot{\bf Z}_{l}$. However, calculation of the $n$-term sums (39) and (42) in algorithms II and III, respectively, is computer-intensive.

4.3 Numerical results

For numerical experimentation we have applied two types of flows: the so-called cosine flows introduced in RCZ , and their curls considered in ACZ2 . They are of interest in that the latter have a pointwise zero vorticity (kinematic) helicity, and the former have a pointwise zero velocity helicity, and nevertheless they are capable of both small- and large-scale magnetic field generation (see ibid). Involving a small number of trigonometric functions, they are particularly useful for calculating the Taylor series coefficients (34).

The cosine flows are defined as follows:

[TABLE]

Here ${\bf a}=(a_{1},a_{2},0)$ and ${\bf b}=(b_{1},b_{2},0)$ are constant horizontal vectors, and

[TABLE]

so that the r.m.s. flow velocity is unity.

Because of many symmetries of the cosine flows, all entries of the eddy diffusivity tensor $\boldsymbol{\mathfrak{D}}$ vanish, except for five pairs (see RCZ ):

[TABLE]

Consequently, the minimum eddy diffusivity (29) takes a simple form:

[TABLE]

We have considered the particular sample flow for

[TABLE]

and used Mathematica again to implement algorithm I: we have calculated exactly the coefficients of the expansions (18) and (32) of solutions to the auxiliary problems (11) and (25) using the recurrence relations (19) and (33), and of the coefficients $\boldsymbol{\mathfrak{D}}_{mk}^{(n)}$ of the series (34) up to order $\eta^{-49}$. The precise coefficients ${\bf s}_{k}^{(49)}$ and ${\bf g}_{mk}^{(49)}$ require about 2 Gbytes of memory for storage (in the ASCII form). For the flow (43), the coefficients $\boldsymbol{\mathfrak{D}}_{mk}^{(n)}$ turn out to be rational; in the ASCII form, the vectors $\boldsymbol{\mathfrak{D}}_{mk}^{(49)}$ occupy 10 to 20 Kbyte of memory.

Fig. 4 shows the sequence of the ratios $|(\mathfrak{D}^{l}_{mk})^{(2n-1)}/(\mathfrak{D}^{l}_{mk})^{(2n+1)}|^{1/2}$ for five independent entries. The limit of this sequence for $n\to\infty$ is equal to the radius of convergence of the series (34) (regarded as a function of $1/\eta$) for the respective entry (note that due to the antisymmetry in $l$ and $k$, the entries involve only odd powers of $1/\eta$). The figure demonstrates that the series for the five entries have the same radius of convergence and converge for $\eta\apprge 0.5$.

Because precise Mathematica calculations require considerable computer resources, we have not considered high-order Padé approximants. The amount of calculations reduces if meromorphic functions

[TABLE]

are approximated only, in terms of which (see (44))

[TABLE]

The poor quality of the resultant approximation (see Fig. 6(a)) is due to the presence of two Froissart doublets in the approximants of $q_{2}$ and $q_{3}$ (Fig. 6(b)). Gaps are present in the plot in Fig. 6(a), where the approximant of $q_{3}$ becomes negative due to its singular behaviour and thus the square root in (47) can not be extracted. Upon factoring the doublets out (which is simple in Mathematica) in the two plagued approximants, the quality of the approximation becomes very similar to that obtained by using [18/18] approximants of $q_{i}$ (see Fig. 6(c)). Increasing the orders to [24/24] does not significantly improve the approximated $\eta_{\rm eddy}$ (see Fig. 6(d)). A better approximation is obtained if the elements of the eddy diffusivity tensor $\boldsymbol{\mathfrak{D}}$ are Padé-approximated individually (see Fig. 6); this gives reasonably accurate values of $\eta_{\rm eddy}$ for $\eta\apprge 0.02$, which is roughly 25 times larger than the minimum $\eta$, for which the power series in $1/\eta$ for the fields ${\bf s}_{k}$ and ${\bf g}_{mk}$, as well as for the elements of the tensor $\boldsymbol{\mathfrak{D}}$ are convergent.

Upon shifting by a quarter of the period in the vertical coordinate $x_{3}$, the curl of (43) takes the form

[TABLE]

where we now assume the normalizing factor

[TABLE]

for which the r.m.s. flow velocity is again 1. Since this flow possesses all the symmetries of (43), the expression (44) for the minimum eddy diffusivity still applies.

Following algorithm II, for a sample flow (48) for

[TABLE]

we have calculated by Mathematica 49 coefficients $\boldsymbol{\mathfrak{D}}_{mk}^{(n)}$ of the series (34) up to order $\eta^{-49}$. The graph of the ratios $|(\mathfrak{D}^{l}_{mk})^{(2n-1)}/(\mathfrak{D}^{l}_{mk})^{(2n+1)}|^{1/2}$ for the five independent entries (see Fig. 7(a)) of the eddy diffusivity tensor shows that the series (34) converge for $\eta\apprge 2.2$. The highest-order (for this set of coefficients) [25/24] approximants of the entries are free of Froissart doublets (see Fig. 7(b)). They yield a satisfactory approximation of dependence on $\eta$ of the minimum magnetic eddy diffusivity for $\eta\apprge 0.03$, which is roughly 70 times smaller than the bound obtained for convergence of the Taylor series (34) for $\mathfrak{D}^{l}_{mk}$ (see Fig. 7(c),(d)). However, their fidelity is insufficient to reproduce the singularity of the minimum eddy diffusivity observed in Fig. 7(c),(d).

The symmetries of the flow (48), (49) imply, that the neutral modes ${\bf s}_{k}$ reside in invariant subspaces of the magnetic induction operator, that can be categorized in terms of the Fourier harmonics ${\bf b_{n}}{\rm e}^{{\rm i}\bf n\cdot x}$, comprising the Fourier series for ${\bf s}_{k}$ (by virtue of the mode periodicity, the wave vectors $\bf n$ have integer components). Only the harmonics that have the following properties enter the Fourier series for ${\bf s}_{k}$:

$\bullet$ $\bf b_{n}=b_{-n}$ are real;

$\bullet$ the numbers $n_{1}$ and $n_{1}/2+n_{2}+n_{3}$ are even;

$\bullet$ ${\bf s}_{1}$ and ${\bf s}_{2}$ are symmetric in $x_{3}$ (i.e., $b_{\bf n}^{1}=b_{{\bf n}^{*}}^{1},~{}b_{\bf n}^{2}=b_{{\bf n}^{*}}^{2},~{}b_{\bf n}^{3}=-b_{{\bf n}^{*}}^{3}$) and ${\bf s}_{3}$ is antisymmetric in $x_{3}$ (i.e., $b_{\bf n}^{1}=-b_{{\bf n}^{*}}^{1},~{}b_{\bf n}^{2}=-b_{{\bf n}^{*}}^{2},~{}b_{\bf n}^{3}=b_{{\bf n}^{*}}^{3}$), where ${\bf n}^{*}=(n_{1},n_{2},-n_{3})$.

We have checked that for $\eta>\eta_{\rm cr}=0.00420516$ small-scale modes from the subspace, where ${\bf s}_{1}$ and ${\bf s}_{2}$ are located, are not generated, but at $\eta=\eta_{\rm cr}$ the small-scale generation starts in the subspace, where ${\bf s}_{3}$ resides. It is known ZPF ; RCZ ; ACZ2 that the point of the onset of the small-scale generation is typically associated with a singularity of the $\alpha$-effect or eddy diffusivity tensors; this is the case for the flow under consideration. We observe in Fig. 7(d) that the least-squares fit by a hyperbola through 20 computed values of minimum eddy diffusivity at equispaced points in the interval $0.0045\leq\eta\leq 0.014$ is very accurate (actually, only 19 points out of the 20 are shown; for the smallest $\eta=0.0045$, eddy diffusivity -20.217501, also well approximated by the hyperbola, is out of the vertical range of Fig. 7(d)). The hyperbolic fit yields the location of the singularity at $\eta=0.0041988$ (the vertical asymptote is shown by a dashed line in Fig. 7(d)) which is very close to the point of the onset of the small-scale generation $\eta_{\rm cr}=0.00420516$ computed by spectral methods; the hyperbola through the three smallest $\eta$ from this interval yields a closer value 0.00420233 .

5 Conclusions

We have tested Padé approximants of the $\alpha$-effect and eddy diffusivity tensors, responsible for generation of large-scale fields, as functions of the respective molecular diffusivity: the viscosity $\nu$ when hydrodynamic perturbations are studied, and the magnetic diffusivity $\eta$ when kinematic dynamo problem is under scrutiny. We have tried different computational tools: Fortran codes relying on the floating point arithmetics, and Mathematica for symbolic and arbitrary precision calculations. A relatively high (several dozens) order of Padé approximants is needed to obtain a reasonable accuracy of approximation of the tensor entries. For this, high precision of computations (in particular, the quadruple precision in Fortran) has proved indispensable. For our sample flows the Padé-approximated tensors yield large-scale magnetic field growth rates to satisfactory accuracy for $\eta$, several dozen times smaller then those, for which power series in the inverse molecular diffusivity converge, for both large-scale generating mechanisms (the $\alpha$-effect and negative magnetic eddy diffusivity).

Application of these techniques in computational fluid dynamics and magnetohydrodynamics seems natural for estimating transport coefficients quantifying the influence of small scales on the evolution of large-scale fields in the spirit of Large Eddy Simulation methods. Our findings, while promising, suggest that to achieve this goal additional algorithms are needed for determination

$\bullet$ of Froissart doublets in approximants of tensor entries and their elimination (the approach of BLM may prove useful for monitoring the absence of the doublets);

$\bullet$ of the interval in molecular diffusivity, where the approximation is sufficiently accurate;

$\bullet$ of the realistic orders of a Padé approximant, for which the length of such interval is close to the maximum.

It is relatively easy to perform these tasks manually by trial-and-error methods — the difficulty lies in performing them automatically.

Acknowledgements

This work was supported by CMUP (Centro de Matemática da Universidade do Porto, UID/ MAT/00144/2019) (SG$+$VZ) and SYSTEC (Centro de Investigação em Sistemas e Tecnologias, POCI-01-0145-FEDER-006933/SYSTEC) (RC), which are funded by FCT with national (MCTES) and European structural funds through the programs FEDER (Fundo Europeu de Desenvolvimento Regional / European Regional Development Fund) under the partnership agreement PT2020, and projects STRIDE [NORTE-01-0145-FEDER-000033] funded by FEDER – NORTE 2020 (SG$+$VZ $+$RC) and MAGIC [POCI-01-0145-FEDER-032485] funded by FEDER via COMPETE 2020 – POCI (SG$+$VZ).