L Chen

W Herreman

K Li

P W Livermore

J W Luo

A Jackson

The optimal kinematic dynamo driven by steady flows in a sphere

Keywords:

This is a repository copy of The optimal kinematic dynamo driven by steady flows in a sphere.

Reuse

Items deposited in White Rose Research Online are protected by copyright, with all rights reserved unless indicated otherwise. They may be downloaded and/or printed for private study, or other acts as permitted by national copyright laws. The publisher or other rights holders may allow further reproduction and re-use of the full text version. This is indicated by the licence information on the White Rose Research Online record for the item.

Takedown

If you consider content in White Rose Research Online to be in breach of UK law, please notify us by emailing eprints@whiterose.ac.uk including the URL of the record and the reason for the withdrawal request.

Introduction

The Earth, as well as certain other planets, moons and stars, is known to possess a self-sustaining magnetic field. Such magnetic fields are generated via fluid motion in a confined conductive interior through a process that is termed dynamo action [START_REF] Moffatt | Magnetic field generation in electrically conducting fluids[END_REF]). Mathematically such dynamics can be modelled by a coupled system of differential equations referred to as magnetohydrodynamics (MHD). If we prescribe the flow field U and ignore the backreaction from the magnetic field B on U , we obtain a simple system called the kinematic dynamo problem. In this reduced system, one can study how a flow field amplifies a seed magnetic field. Early numericists found flow fields that are capable of dynamo action in a sphere with electrically insulating boundary conditions [START_REF] Backus | A class of self-sustaining dissipative spherical dynamos[END_REF][START_REF] Pekeris | Kinematic dynamos and the earth's magnetic field[END_REF][START_REF] Kumar | A spectral solution of the magneto-convection equations in spherical geometry[END_REF][START_REF] Gubbins | Numerical solutions of the kinematic dynamo problem[END_REF][START_REF] Dudley | Time-dependent kinematic dynamos with stationary flows[END_REF]. However, there is no universal recipe on how to obtain such flow fields; known dynamo solutions do not necessarily share similar spatial structures [START_REF] Gubbins | Implication of kinematic dynamo studies for the geodynamo[END_REF].

Optimisation offers a tool to explore the many possible ways of getting dynamo action. Historically optimisation of dynamo models has taken place predominantly with a restricted type of velocity field: for example, working in a sphere, Love & Gubbins (1996b) optimized the flow proposed by [START_REF] Kumar | A spectral solution of the magneto-convection equations in spherical geometry[END_REF] (abbreviated as KR) while [START_REF] Holme | Optimised axially-symmetric kinematic dynamos[END_REF] optimized flows proposed by [START_REF] Dudley | Time-dependent kinematic dynamos with stationary flows[END_REF] (abbreviated as DJ); later studies by Love & Gubbins (1996a); [START_REF] Holme | Three-dimensional kinematic dynamos with equatorial symmetry: Application to the magnetic fields of uranus and neptune[END_REF]; Gubbins et al. (2000a) and Gubbins et al. (2000b) yielded more efficient solutions using only poloidal flows of KR type. Allied work also includes the search for optimal dynamos driven by ABC flows in a periodic domain [START_REF] Alexakis | Searching for the fastest dynamo: Laminar abc flows[END_REF], maximization of helicity for the cylindrical Riga experiment [START_REF] Stefani | Velocity profile optimization for the riga dynamo experiment[END_REF], optimisation of boundary driven flows for the spherical Madison plasma experiment [START_REF] Khalzov | Optimized boundary driven flows for dynamos in a sphere[END_REF], optimisation of flow structures and conducting layers for the cylindrical VKS experiment [START_REF] Marie | Numerical study of homogeneous dynamo based on experimental von karman type flows[END_REF][START_REF] Ravelet | Towards an experimental von karman dynamo: numerical studies for an optimized design[END_REF], and most recently, the study of optimal forcing for chaotic dynamos [START_REF] Sadek | Optimal length scale for a turbulent dynamo[END_REF]. From a more fundamental point of view, we want to ask: "Imposing no constraints other than the spherical boundaries, what is the most efficient dynamo, and what are its characteristics?". Optimisation that potentially covers an infinitely large parameter space of flow structures, is needed to answer these questions.

A variational approach has successfully optimized kinematic dynamo models in a cube [START_REF] Willis | Optimization of the magnetic dynamo[END_REF][START_REF] Chen | Optimal dynamo action by steady flows confined in a cube[END_REF][START_REF] Herreman | Minimal flow perturbations that trigger kinematic dynamo in shear flows[END_REF]. This method sets up a Lagrangian consisting of an objective functional that usually measures the magnetic energy at later times, subject to a number of constraints, then searches for the optimal state, i.e. a stationary point of the Lagrangian in phase space with respect to the variation of each variable, including that of the velocity field and the seed magnetic field. The important control parameter in optimisation is the magnetic Reynolds number Rm. It measures the effects of advection and stretching of magnetic field compared to the effects of Ohmic decay. In order to sustain dynamo action, Rm must be sufficiently large such that field creation outweighs field destruction through dissipation, and indeed Rm must exceed certain known lower bounds [START_REF] Backus | A class of self-sustaining dissipative spherical dynamos[END_REF][START_REF] Busse | A necessary condition for the geodynamo[END_REF][START_REF] Childress | Thèorie magnétohydrodynamique de l'effet dynamo[END_REF][START_REF] Proctor | On Backus' necessary condition for dynamo action in a conducting sphere[END_REF][START_REF] Proctor | Necessary conditions for the magnetohydrodynamic dynamo[END_REF]. Historically, different definitions of Rm have been used, and one needs to be careful which measure to use. It has been proven that the kinetic energy-based magnetic Reynolds number (denoted by Rm u ) has no lower bound [START_REF] Proctor | Energy requirement for a working dynamo[END_REF], since a given flow field can be shrunk to a smaller region with a lower value of Rm u while dynamo action is retained through the increased shear. As a consequence we use an enstrophybased normalization to define Rm in our optimisation which is equivalent to normalizing the global magnitude of shear.

In this work, we extend these previous studies to find the optimal kinematic dynamo in a sphere. As a first step, our aim is to look for optimal dynamo action among general steady and incompressible flows in a sphere with widely-used electrically insulating and no-slip boundary conditions. This requires the use of an inner product defined over infinite space, which we refer to as an all-space norm, and a primitive formulation of the Lagrangian. The all-space norm removes all surface integrals during the derivation of the variations of the Lagrangian and is key to the successful application of our method. We also developed an adjoint model for this particular optimisation problem, akin to those first written down by [START_REF] Namikawa | Kinematic dynamo problem[END_REF] and [START_REF] Li | Variational data assimilation for the initial value dynamo problem[END_REF]. Lastly, a specially constructed orthogonal Galerkin basis [START_REF] Livermore | Quasi-L p norm orthogonal Galerkin expansions in sums of Jacobi polynomials[END_REF][START_REF] Livermore | Galerkin orthogonal polynomials[END_REF][START_REF] Li | An optimal galerkin scheme to solve the kinematic dynamo eigenvalue problem in a full sphere[END_REF] for the vector fields ensures a rapid convergence in spectral space towards the optimal state.

The structure of this article is as follows: we introduce the methodology in §2 and present the results in §3, wherein we describe the optimal flow field, together with its magnetic eigenvector at the minimal critical magnetic Reynolds number for the onset of dynamo action Rm c,min = 64.45. We discuss the characteristics of this optimal solution, contrasting them to theoretical bounds and other spherical dynamo models; we also compare Rm c,min directly with the optimized dynamo models in a cube [START_REF] Chen | Optimal dynamo action by steady flows confined in a cube[END_REF] by rescaling typical length scales and volumes in §3.8. We close with conclusions in §4.

Methods

Definition of the optimisation objective

Consider a sphere V * filled with an electrically conducting fluid. Outside the sphere, we suppose a current free region V * that extends to infinity. If L * denotes dimensional values for the spherical radius and η * denotes the fluid's magnetic diffusivity, we then non-dimensionalize the flow field using units [U ] = ω * L * for velocity, [x] = L * for length and [t] = L * 2 /η * for time. Here ω * is the dimensional root mean enstrophy,

ω * = 1 V * V * (∇ × U * ) 2 dV * , (2.1)
V * is the dimensional volume of the sphere, and U * is the dimensional velocity field; we set η * = 1/(µ * 0 σ * ), where µ * 0 is permeability of free space, σ * is the electrical conductivity. The magnetic Reynolds number is then defined as

Rm = ω * L * 2 η * . (2.2)
Note here it is necessary to use the enstrophy measure of flow strength for reasons explained in the introduction. In this model, we impose no-slip boundary conditions on the flow field U * . This then leads to the condition that the root mean enstrophy ω * = S * , where S * = 1 V * V * 2 e * ij e * ij dV * is the dimensional global shear magnitude, and e * ij is the dimensional strain rate tensor.

The non-dimensionalized flow field U = U (x) lives within the parameter space E u of all steady and incompressible fields with no-slip boundary conditions, i.e. ∀ U ∈ E u ,

∇ • U = 0, x ∈ V, (2.3) U | Σ-= 0, (2.4)
where we denote by Σ ± either the outer surface (+) or the inner surface (-) of the sphere at the interface between non-dimensionalized domains V and V. In a kinematic approach, the magnetic field B and the electric field E need to satisfy the non-dimensional equations,

E = ∇ × B -Rm U × B, x ∈ V, (2.5)
∂ t B = -∇ × E, x ∈ V ∪ V, (2.6) ∇ × B = 0, x ∈ V, (2.7) ∇ • B = 0, x ∈ V ∪ V. (2.8)
We have used an arbitrary scale [B] to non-dimensionalize the magnetic field, and [E] = ω * L * [B] to non-dimensionalize the electric field. Unlike in [START_REF] Willis | Optimization of the magnetic dynamo[END_REF] and [START_REF] Chen | Optimal dynamo action by steady flows confined in a cube[END_REF], we introduce Ampère's law (2.5) and Faraday's law (2.6) separately as a primitive formalism rather than collectively as the induction equation (2.24), for reasons that we will discuss later in §2.3. This step is needed to impose a truly current-free exterior region V in the Lagrangian formalism discussed in §2.2. Equation (2.7) is the current-free condition, equation (2.8) is Gauss' law for magnetism, and we also impose the continuity conditions,

B| Σ+ -B| Σ-= 0, r × E| Σ+ -E| Σ-= 0, (2.9)
where r is the radial unit vector. Since the flow is steady, the kinematic dynamo problem admits exponential solutions for the magnetic field, (2.10) where b(x) is an eigenvector, γ is the growth rate and Ω is an oscillation frequency. As t → ∞, one can expect that the eigenvector with the largest growth rate will eventually dominate the solution, assuming an arbitrary noisy initial magnetic field B 0 = B(x, 0). In this optimisation study, we want to identify the flow field U ∈ E u that maximizes the growth rate γ for a given value of Rm, and subsequently find the lower bound on the critical magnetic Reynolds number, denoted by Rm c,min , below which no kinematic dynamo is possible ∀ U ∈ E u .

B(x, t) ∼ b(x) e (γ+iΩ)t ,

Euler-Lagrange formalism in primitive format

Since we want to optimize the growth rate γ, it seems most obvious to propose an optimisation method that directly maximizes γ. This can indeed be done in low-dimensional optimisation problems, see for example [START_REF] Holme | Optimised axially-symmetric kinematic dynamos[END_REF]. However, considering the infinitely large dimension of the functional space E u we have, optimising the growth rate directly is not practical. Instead, we adapt the optimisation strategy of [START_REF] Willis | Optimization of the magnetic dynamo[END_REF]; [START_REF] Chen | Optimal dynamo action by steady flows confined in a cube[END_REF]; [START_REF] Herreman | Minimal flow perturbations that trigger kinematic dynamo in shear flows[END_REF] that proved successful in a cube. These methods are designed to find the optimal flow field U and seed magnetic field B 0 that maximize a specified global measure of B T , where B T = B(x, T ) is the magnetic field at some finite time T that needs to be long enough to overcome initial transient growth phases. Optimising B 0 then allows us to reach this exponential regime faster. It is worth remarking that, in general, the initial field B 0 that has greatest projection onto the most rapidly growing eigenvector B T is the corresponding eigenvector of the related adjoint eigenvalue problem [START_REF] Farrell | Generalized stability theory. part 1: Autonomous operators[END_REF] (distinct from the adjoint formalism we will use). However, we opt to discover the structure of B 0 directly by an optimising scheme as we will set out below, rather than introducing the additional complexity of solving the adjoint eigenvalue problem.

The variational method we use is built on the following Lagrangian:

L = ln (B T ) 2 -λ 1 1 V (∇ × U ) 2 -1 -λ 2 (B 0 ) 2 -1 -Π∇ • U - T 0 ψ † ∇ • B dt - T 0 B † • [∂ t B + ∇ × E] dt - T 0 E † • [σ r E + Rm U × B -∇ × B] dt. (2.11)
We denote the volume of domain V by V , the all-space integral:

. . .

= ∞ 0 2π 0 π 0 • • • r 2 sin θ dθdφdr, (2.12)
and we suppose that all functions and fields are regular at r = 0 and decay sufficiently quickly as r → ∞ for this measure to make sense. The first term in (2.11) is the objective functional that we want to maximize, namely the logarithm of the energy of magnetic field at time T . From the second to the last term in (2.11), we impose a series of constraints in this model, namely the normalization constraint for U and B 0 , the solenoidal conditions on U and B, Faraday's law and Ampère's law. Note that the flow field U is entirely confined to the sphere V, so there is no contribution to the enstrophy from V. The variables λ 1 , λ 2 , Π(x), ψ † (x, t), B † (x, t) and E † (x, t) are Lagrange multipliers. Π, ψ † , B † , E † are often referred to subsequently as adjoint fields. The number σ r is a relative electrical conductivity:

σ r = 1, x ∈ V, 0, x ∈ V. (2.13) At the optimum the Lagrangian L(U , B, B 0 , B T , E, λ 1 , λ 2 , Π, ψ † , B † , E †
) must be stationary, meaning that δL = 0 for arbitrary variations in each of its variables. Each variational derivative needs to disappear separately, which defines eleven Euler-Lagrange (EL) equations in our optimisation problem. The physical constraints define six equations: the solenoidal conditions as in (2.3) and (2.8), the normalization constraints on B 0 and U , Ampère's law and Faraday's law. The remaining five non-trivial variational derivatives are

δL δB T = 2B T (B T ) 2 -B † T , (2.14) δL δE = -σ r E † -∇ × B † , (2.15) δL δB = ∂ t B † + Rm U × E † + ∇ × E † + ∇ψ † , (2.16
) (2.18) where B † T = B † (x, T ) and B † 0 = B † (x, 0). A number of boundary terms arise during the derivation of these non-trivial variational derivatives by partial integration, and their sum needs to cancel at the optimum in order to have a consistent theory. In this model, the magnetic field B decays at least as fast as r -3 in the insulating region V, where r is the distance to the origin, so all boundary variations on the surface of V as r → ∞ are negligible. We then need all boundary terms to cancel on Σ ± :

δL δU = Rm T 0 B × (-σ r E † ) dt -2λ 1 ∇ × ∇ × U + ∇Π, (2.17)
δL δB 0 = B † 0 -2λ 2 B 0 ,
0 = -Π(r • δU ) Σ-dS -2λ 1 δU × (∇ × U ) • r Σ-dS + T 0 (r × B † ) • δE Σ--(r × B † ) • δE Σ+ dS dt + T 0 (r × E † ) • δB Σ--(r × E † ) • δB Σ+ dS dt - T 0 ψ † (r • δB) Σ--ψ † (r • δB) Σ+ dS dt.
(2.19)

If we suppose that the field variations δU satisfy no-slip conditions (2.4), the first two boundary terms are zero. Then, the inclusion of the external region V allows us to cancel the remaining boundary terms by matching them from Σ -to Σ + . The primitive formalism also allows us to simplify the treatment of boundary terms on Σ ± : using continuity of the magnetic field variations δB and the tangential electric field variations δE as in (2.9), we immediately obtain the boundary conditions for the adjoint fields,

ψ † | Σ+ -ψ † | Σ-= 0, r × (B † | Σ+ -B † | Σ-) = 0, r × (E † | Σ+ -E † | Σ-) = 0. (2.20)
Because of (2.15) and the presence of ∇ψ † in (2.16), the adjoint magnetic field has a gauge degree of freedom B † → B † + ∇φ. Therefore, we can add a supplementary requirement that the adjoint magnetic field is also solenoidal,

∇ • B † = 0, x ∈ V ∪ V. (2.21)
along with the boundary condition,

r • (B † | Σ+ -B † | Σ-) = 0. (2.22)
By requiring the sum of all boundary terms to vanish, we just derived the continuity conditions for ψ † , B † and E † . This brings us to the conclusion that, with a gauge degree of freedom, adjoint fields B † and E † satisfy exactly the same boundary conditions as the direct fields B and E. This is an eminently desirable property, since it allows us to use the same solenoidal Galerkin expansions for both direct field B and adjoint field B † later in §2.4.

Reduced set of Euler-Lagrange equations

Now we have dealt with the boundary terms, we can eliminate the fields E and E † from the previous "primitive" set of EL equations and present the reduced set of EL equations that we will actually solve. We build an optimisation loop that will locate the stationary point of the Lagrangian (2.11) iteratively. Upon initialization (possibly combined with an update) and at each iteration in the optimisation loop, the flow field U and the seed magnetic field B 0 satisfy the normalization constraints,

1 V (∇ × U ) 2 = 1, (B 0 ) 2 = 1, (2.23)
the solenoidal conditions (2.3) and (2.8), and the boundary conditions (2.4) and (2.9).

In the optimisation loop itself, there are three parts: a forward integration, a backward integration and an update, similar to the procedure described by [START_REF] Kerswell | An optimization approach for analysing nonlinear stability with transition to turbulence in fluids as an exemplar[END_REF]. First we solve the forward/direct problem. In region V, the magnetic field satisfies (2.7). In region V, by combining Ampère's law and Faraday's law, we solve the induction equation forward in time for t : 0 → T , (2.24) subject to the solenoidal condition (2.8) and continuity conditions (2.9). With the field B T known, we then initialize the adjoint magnetic field at time t = T as

∂ t B = Rm∇ × (U × B) -∇ × (∇ × B), x ∈ V,
B † T = 2B T (B T ) 2 . (2.25)
Elimination of E † in the adjoint Ampère's and Faraday's laws derived from (2.15) and (2.16), yields the adjoint problem:

∂ t B † = Rm U × (∇ × B † ) + ∇ × ∇ × B † -∇ψ † , x ∈ V, (2.26) ∇ × B † = 0, x ∈ V.
(2.27)

In the second part, we solve (2.26) backward in time for t : T → 0 while respecting the solenoidal condition (2.21) and the boundary conditions (2.20) and (2.22). In the third part we consider the update using variational derivatives (2.17) and (2.18) which lead to the only two EL equations that are not automatically satisfied. We iterate through the optimisation loop so that δL/δU → 0 and δL/δB 0 → 0. To calculate the integral term in δL/δU , we must know B and B † at all times. Lagrange multipliers λ 1 , λ 2 , Π are still undetermined at this stage. As in [START_REF] Chen | Optimal dynamo action by steady flows confined in a cube[END_REF] and [START_REF] Pringle | Minimal seeds for shear flow turbulence: using nonlinear transient growth to touch the edge of chaos[END_REF], we use a preconditioned ascent method in which we use only part of the second variation,

δ 2 L ≈ 2λ 1 δU • ∇ 2 δU -2λ 2 δB 0 • δB 0 + . . . (2.28)
with respect to variables U and B 0 to allow faster convergence. To be more specific, let us denote the update as

U : = U + α 1 ∆U , B 0 : = B 0 + α 2 ∆B 0 , (2.29)
with α 1 , α 2 ≪ 1 being the two relaxation parameters. We calculate the increments ∆U and ∆B as if (2.28) would be exact such that α 1 = α 2 = 1 would correspond to a Newton step,

δL δU + 2λ 1 ∇ 2 ∆U = 0, δL δB 0 -2λ 2 ∆B 0 = 0. (2.30)
The magnetic field update is always rather trivial, since we can explicitly evaluate ∆B 0 , to arrive at

B 0 := B 0 + α 2 2λ 2 δL δB 0 . (2.31)
The value of λ 2 is finally set by requiring that the updated B 0 remains normalized as in (2.23). The supplementary requirement (2.21) is necessary for the updated B 0 to remain within the same parameter space as B with the correct boundary conditions. The velocity field update is much less trivial. Due to the presence of the Laplacian operator ∇ 2 in front of ∆U , we cannot write at this point an explicit update formula for the flow. There is also a second difficulty: owing to formula (2.17) for δL/δU , there is no guarantee that the updated velocity field will automatically satisfy the no-slip boundary conditions. We tackle both issues by projecting U on a very particular vector field basis. The explicit equation for the velocity field update is given in §2.5.

2.4.

Galerkin expansions for U , B and B † We solve this optimisation problem using a Galerkin method. This method builds the solenoidal condition and boundary conditions into the field expansions of U , B and B † . We denote

      U B B †       = lmax l=1 l m=-l nmax n=1       t nml U t nml + s nml U p nml T nml B t nml + S nml B p nml T † nml B t nml + S † nml B p nml       .
(2.32)

Here t nml , s nml , T nml (t), S nml (t), T † nml (t), and S † nml (t) are the spectral coefficients. We can use the same basis vector fields of B to represent B † since the adjoint magnetic field is solenoidal in the present gauge and has the same boundary conditions as B. We use the following as the toroidal (with a superscript t) and poloidal (with a superscript p) basis vector fields for U and B: (2.34) where Y m l (θ, φ) are the spherical harmonics, and f l n (r), g l n (r), h l n (r), and k l n (r) are the radial basis functions. In Appendix A we provide more technical details on how to use this basis. Here we only present the essential properties. All basis fields are infinitely differentiable at the origin r = 0. This gives the form of the polynomials as r l+1 (a 0 + a 1 r 2 +a 2 r 4 +• • • ) for each l. All basis fields satisfy boundary conditions at r = 1. These are no-slip boundary conditions for the flow field and continuity conditions for the magnetic field. This requires

U t nml (r, θ, φ) =∇ × f l n (r) Y m l (θ, φ)r , r ∈ V, U p nml (r, θ, φ) =∇ × ∇ × g l n (r) Y m l (θ, φ)r , r ∈ V, (2.33) B t nml (r, θ, φ) = ∇ × h l n (r) Y m l (θ, φ)r , r ∈ V, 0, r ∈ V, B p nml (r, θ, φ) = ∇ × ∇ × k l n (r) Y m l (θ, φ)r , r ∈ V, -lk l n (1)∇ r -(l+1) Y m l (θ, φ) , r ∈ V,
f l n (1) = 0, g l n (1) = ∂g l n ∂r (1) = 0, h l n (1) = 0, ∂k l n (1) ∂r + l k l n (1) = 0. (2.35)
The vector field basis is orthonormal with respect to a predefined inner product. The radial functions for the velocity field and for the magnetic field are chosen such that

1 V (∇ × U ρ nml ) • (∇ × U ρ n ′ m ′ l ′ ) = δ nn ′ δ mm ′ δ ll ′ , B ρ nml • B ρ n ′ m ′ l ′ = δ nn ′ δ mm ′ δ ll ′ , ρ = t, p.
(2.36)

Using partial integration, the no-slip boundary conditions, and the solenoidal property of the basis functions, (2.36) is equivalent to

1 V -U ρ nml • ∇ 2 U ρ n ′ m ′ l ′ = δ nn ′ δ mm ′ δ ll ′ , ρ = t, p.
(2.37)

The orthonormal condition (2.37) simplifies the update step for the flow, which we will describe later in §2.5. For the velocity field U the radial functions are

f l n (r) = N f n,l r l+1 P (0,l+1/2) n (2r 2 -1) -P (0,l+1/2) n-1 (2r 2 -1) , (2.38) g l n (r) = N g n,l r l+1 3 i=1 c i P (0,l+1/2) n+2-i (2r 2 -1), c 1 = 2l + 4n + 1, c 2 = -2(2l + 4n + 3), c 3 = 2l + 4n + 5, (2.39)
and for the magnetic field B the radial functions are

h l n (r) = N h n,l r l+1 (1 -r 2 )P (2,l+1/2) n-1 (2r 2 -1), (2.40) k l n (r) = N k n,l r l+1 c 1 P (0,l+1/2) n (2r 2 -1) + c 2 P (0,l+1/2) n-1 (2r 2 -1) , c 1 = n(2l + 2n -1), c 2 = -(n + 1)(2n + 2l + 1).
(2.41)

where P (α,β) λ are Jacobi polynomials. The numbers N j n,l , j = f, g, h, k are normalization factors and they are chosen such that the basis is orthonormal as in (2.36). Without this factor the basis is only orthogonal. The procedure of constructing this polynomial basis is described by [START_REF] Livermore | A compendium of galerkin orthogonal polynomials[END_REF], [START_REF] Livermore | Quasi-L p norm orthogonal Galerkin expansions in sums of Jacobi polynomials[END_REF], and [START_REF] Livermore | Galerkin orthogonal polynomials[END_REF]. The scalar products are positive definite and symmetric and they have particular benefits when it comes to calculating norms:

1 V (∇ × U ) 2 = l,m,n t nml 2 + s nml 2 , B 2 = l,m,n T nml 2 + S nml 2 .
(2.42) 2.5. The optimisation algorithm We initialize the algorithm with randomly distributed spectral coefficients or we restart from previously stored fields. In all cases the initial coefficients are normalized such that l,m,n

t nml 2 + s nml 2 = 1, l,m,n T nml (0) 2 + S nml (0) 2 = 1. (2.43)
We then enter the optimisation loop. We inject the field expansion into the forward problem (2.24) and project it onto the magnetic field basis using the orthonormality condition. This gives a system of ordinary differential equations

∂ t T nml =Rm B t nml • [∇ × (U × B)] + n ′ T l nn ′ T n ′ ml , (2.44) ∂ t S nml =Rm B p nml • [∇ × (U × B)] + n ′ P l nn ′ S n ′ ml , (2.45)
to be integrated for t : 0 → T . On the right hand side, the matrix elements

T l nn ′ = B t nml , ∇ 2 B t n ′ ml , P l nn ′ = B p nml , ∇ 2 B p n ′ ml , (2.46)
capture the effect of magnetic diffusion; they are precomputed in Mathematica. To calculate the non-linear induction term in practice, we choose to introduce a physical space grid that is composed of Fourier grid points in φ and Gaussian quadrature points in r and θ, more details are given in Appendix A.3. The numerical scheme for time-integration is Crank-Nicolson for the diffusive terms and second order Adams-Bashforth for the induction terms. The initialization condition for the adjoint magnetic field is given by the magnetic field at time T as in (2.25). In terms of spectral coefficients, we have

T † nml (T ) = 2 T nml (T ) l,m,n T nml (T ) 2 + S nml (T ) 2 , S † nml (T ) = 2 S nml (T ) l,m,n T nml (T ) 2 + S nml (T ) 2 . (2.47)
A similar projection of the adjoint induction equation (2.26) onto the basis of toroidal magnetic field B t yields

∂ t T † nml = Rm B t nml • [U × (∇ × B † )] - n ′ T l nn ′ T † n ′ ml , (2.48)
which we integrate backwards in time for t : T → 0, similarly for the poloidal magnetic field B p , see Appendix A.3. We also use a similar time-discretisation scheme as in the forward model. The diffusion matrices T l nn ′ and P l nn ′ are unchanged from (2.46). The field ψ † is never used within our numerical scheme: when we use this Galerkin basis, the projection of ψ † onto the basis of B gives:

B t nml • ∇ψ † = B p nml • ∇ψ † = 0.
The seed magnetic field update is simple, since we can directly project the update formula (2.31) onto the vector field basis. In terms of spectral coefficients, we have

T nml (0) := T nml (0) + α 2 2λ 2 T † nml (0) -2λ 2 T nml (0) = T nml (0) + δT nml (0), (2.49)
for the toroidal magnetic field where δT nml denotes the incremental change of spectral coefficient between consecutive iterations. We use the same update formula for the poloidal field. By requiring that the updated seed magnetic field remains normalized, we find a quadratic polynomial that sets λ 2 :

l,m,n 4λ 2 2 (α 2 -2) + 4λ 2 (1 -α 2 ) T nml (0) T † nml (0) + S nml (0) S † nml (0) +α 2 T † nml (0) 2 + S † nml (0) 2 = 0. (2.50)
Provided α 2 < 2, there is always just one positive root of (2.50) which we choose.

To update the velocity field U , we must first calculate the integral term that appears in δL/δU . We store the spectral coefficients T nml (t i ) and S nml (t i ) of B at all discrete times t i . Since the code is parallelized for a distributed memory architecture, a checkpointing strategy is not required. During the backward time integration, we then calculate the integral term projected on the vector basis of U using the trapezoidal rule,

I t,p nml = N i=0 β i ∆t U t,p nml • B(x, t i ) × (∇ × B † (x, t i )) , (2.51)
with the superscript indicating either toroidal t or poloidal p components, ∆t is the timestep size, N is the total number of timesteps and

β 0 = β N = 0.5, β i = 1, i = 1, • • • , N -1
are the integration weights. In a similar vein to the calculation of the induction terms, we need a transformation to physical space to calculate the non-linear terms in (2.51), see Appendix A.3. We expand the field ∆U onto the vector field basis of U and project (2.30) onto the vector field basis. Using the property (2.37), we then find an explicit formula for the spectral coefficients of ∆U that finally sets the update formula for the spectral coefficients of U . The toroidal coefficients are given by

t nml := t nml + α 1 2λ 1 (I p nml -2λ 1 t nml ) = t nml + δt nml , (2.52)
where δt nml denotes the incremental change; a similar formula applies for the update of poloidal coefficients. The value of λ 1 is calculated in the same way as λ 2 . Once this update is performed, we return to the forward integration and iterate through the whole loop until convergence. At each iteration in the optimisation loop, we measure progress through the incremental change in the fields

r B0 2 = l,m,n δT nml (0) 2 + δS nml (0) 2 , r U 2 = l,m,n δt nml 2 + δs nml 2 , (2.53)
and the total residue is

r t = r 2 U + r 2 B0 .
(2.54)

Besides tracking the residue, we also follow how (B T ) 2 evolves. We terminate the optimisation when r t is smaller than a fixed value and when (B T ) 2 no longer changes significantly. The relaxation parameters α 1 , α 2 are prescribed and adjusted depending on the value of r t . We typically choose 0.1 < α j < 0.5, j = 1, 2 for the first 5 iterations if starting from random initial conditions, then sharply reduce α j to about 0.01. This reduction in α j is to prevent the optimizer from oscillating around a fixed point in the parameter space of U . If the residue r t < 0.005, which indicates we are close to the optimum, we increase the relaxation parameters again to 0.04 < α j < 0.1 to speed up (lmax, nmax) (16, 8) (24, 12) (32, 32) γ -6.89873355 -6.92808061 -6.92884872

Table 1. Benchmark of the growth rate for the most unstable magnetic field eigenmode (m=0) for the axisymmetric MDJ t1s2 flow at Rm = 545.9, compared with a previously reported value of γ = -6.92884871 [START_REF] Livermore | A comparison of numerical schemes to solve the magnetic induction eigenvalue problem in a spherical geometry[END_REF].

convergence. We directly use step size 0.04 < α j < 0.1 if restarting from previously stored fields, which in general are moderately converged with a relatively small r t . This is close to the value of 0.05 used in [START_REF] Duguet | Minimal transition thresholds in plane couette flow[END_REF]. We also repeated the optimisations with different choices of α j but nevertheless reached the same optimum.

Degeneracy of the optimum

This optimisation problem does not have a unique solution due to an infinite but trivial degeneracy. The symmetry group of the sphere is the orthogonal group O(3), which includes rotational symmetries and reflection symmetries. Any particular transform in the group can be represented using an orthogonal matrix R. For any given flow U or magnetic field B, we can define transformed fields as

U (x) = R T U (Rx), B(x, t) = R T B(Rx, t).
(2.55)

Under the proper transformation rules for the adjoint variables and multipliers, it is not difficult to show that our Lagrangian is invariant under rotations and reflections, meaning that

L(U , B, • • • ) = L( U , B, • • • ). (2.56)
If U is an optimal dynamo that drives the field B, then any transformed field U is also an optimal dynamo that drives the transformed field B. Note that degeneracy should not be confused with symmetries within a particular solution-degeneracy is a general property of any optimal soltution in our setup.

Results

Preliminary tests

Before running the optimisation, we first benchmark our forward model using an axisymmetric flow termed the MDJ t 1 s 1 flow, as studied by [START_REF] Livermore | A comparison of numerical schemes to solve the magnetic induction eigenvalue problem in a spherical geometry[END_REF], [START_REF] Li | An optimal galerkin scheme to solve the kinematic dynamo eigenvalue problem in a full sphere[END_REF]; see Appendix B. We use three different spectral resolutions and the results are shown in Table 1. We see that the highest resolution (l max , n max ) = (32, 32) gives the best match with the known value, but a smaller resolution (l max , n max ) = (24, 12) is good enough to give 3 decimal places of accuracy and has much less computational cost.

For most of the optimisation, we use either (l max , n max ) = (16, 8) or (24, 12), then for verification purposes we add one more optimisation run with (l max , n max ) = (24, 24). We also verify the projection of the variational derivative (2.17) onto the basis of U . This quantity contains information about the adjoint model at all time steps and is important for the update. We set λ 1 = 0 in (2.17) for this test due to the convenient choice of orthonormal properties of the vector bases in (2.37), also to remove the constraint on unity enstrophy in (2.23). We want to compare how the Lagrangian varies with respect to a small perturbation on the flow field U using two methods. First, we compute the variations using a forward model. Since we have all other constraints satisfied, each component of the variation δL is given by

δL ρ nml = ln (B T per ) 2 -ln (B T ) 2 ε , (3.1)
where B T per is the final magnetic field at time T generated by a given random flow field plus a small perturbation: U + εU ρ nml , ρ = t, p. Each component with a different index in n, m, l needs to be calculated separately using a forward model, while B T is computed once without adding perturbations in the random flow field U using a forward model. Second, we calculate the variational derivative δL/δU projected onto the vector space of U as in (2.51) by launching the adjoint model only once. The two variations are related by taking the limit of the perturbation amplitude ε to zero:

lim ε→0 δL ρ nml = δL δU • U ρ nml , ρ = t, p. (3.2)
We then compare the RHS of (3.1) and (3.2) in spectral space by setting ε to a very small value. All realizations used here have the same initial B 0 and U . The test result shows good agreement, see Figure 1.

A systematic survey identifies the Rm c,min

In a systematic survey, we launch optimisations for various values of the control parameters. The timestep ∆t is set to at most 10 -4 . The total integration time is mainly T = 2 except in one optimisation where T = 2.5 for verification. We have initiated the optimisation loops from either random or previously stored fields. The tolerance level on r t gradually decreases in a refinement approach. For each of these runs, we measure the optimal growth rate γ near the end of the time-interval. Detailed information on all optimisation runs is given in the supplementary materials. In Figure 2, we show the growth rate from multiple optimisation runs, performed at different resolutions in order to finally locate the optimal dynamo threshold at Rm c,min ≈ 64.45.

(3.3)

We refer to this optimal solution as model R1. To ensure that we did not get stuck in a local optimum, we ran an independent optimisation at this Rm c,min , without restarting from a previous optimum, which we labeled as model R2. This second independent optimisation confirms the location of the threshold. In fact, model R2 is a degenerate solution (see §2.6) of model R1. In Figure 3 a), we show the logarithm of the magnetic energy as a function of time for optimized solutions near the minimal onset of dynamos at Rm c,min . We clearly see that the exponential regime is well-established and that the optimal dynamos are non-oscillatory (Ω = 0 in (2.10)). This indicates T = 2 is long enough to overcome the initial transient. The red curve highlights the run at the optimal dynamo threshold. An optimisation run generally needs a few thousand iterations to converge from random initial conditions. It takes about 2 minutes per iteration on 24 cores for (l max , n max ) = (24, 12) resolution, ∆t = 10 -4 , T = 2. In Figure 3 b), we show an example of the convergence as a function of iterations from a random start at Rm = 64.45 using model R2. Clearly most of the improvement happens in the first thousand iterations. 3.3. The spatial structure of the optimal solution at Rm c,min

In order to have a complete picture of the optimal dynamo, we analyse the spatial structure of the optimal fields U and B T at Rm c,min in both physical and spectral space.

Analysis in physical space

In physical space, we show the 3D structure of the streamlines of U and B T from model R1 in Figure 4 a) and b), field lines are coloured by the local field intensity, small amplitude structures are not shown in the plots for clarity. In Figure 4 a), the field lines appear to be twisting in a twofold manner. We see a drastic increase in the flow speed in the centre. Correspondingly, we also see a strong magnetic field close to the origin in Figure 4 b). Overall we find a quite complex structure in both U and B T that apparently lacks reflectional symmetry. This lack of reflectional symmetry was also observed in the optimal dynamos in cubes with perfectly conducting or pseudo-vacuum boundaries (homogeneous cases NNN and TTT in [START_REF] Chen | Optimal dynamo action by steady flows confined in a cube[END_REF]). However, we observe a rotational symmetry of π that the cubic models lack, i.e. only even m modes are present when the symmetry axis is aligned with the ẑ direction.

We also observe that the velocity and vorticity fields are nearly parallel in the centre, which gives a strong helical motion, shown in Figure 5 a). The magnetic field seems to be enhanced strongly next to the spiral. This is evident when we combine Figures 4 a)-b) and Figure 5 b), where we show the interaction of the flow field U and the magnetic field B T , together with the structure of the vorticity field ∇ × U . We found the field structure of B T is mostly dipolar outside the sphere, shown in Figure 6 a) where streamlines of the magnetic field extend to two radii. The Mollweide map of the radial magnetic field at the sphere surface also confirms the dominant component is dipolar in Figure 6 b). All plots mentioned in this paragraph are based on model R1; model R2 has the same spatial structures up to a rotation.

Besides looking directly at the optimal fields, we also want to understand what are the averaged spatial distributions. We introduce a radial distribution of kinetic energy

U (r) 2 = U t (r) 2 + U p (r) 2 , U ρ (r) 2 = 1 4π 2π 0 π 0 (U ρ ) 2 sin θ dθdφ, ρ = t, p, (3.4)
for the optimal flow U , toroidal part U t and poloidal part U p . Similarly we define B T (r) 2 , B t T (r) 2 , B p T (r) 2 for the magnetic eigenvector. We also show the radial distribution of two important quantities: helicity and shear (measured by the absolute maximal strain rate). The toroidal part of the radially varying helicity is given by

H t (r) = 1 4π 2π 0 π 0 U t • ∇ × U p sin θ dθdφ. (3.5)
Similarly H p (r) measures the poloidal part. In total the radial distribution of helicity is H(r) = H t (r) + H p (r). The radial profile of shear is given by

S max (r) = max φ∈[0,2π], θ∈[0,π] eig (∇U + ∇U ⊺ )/2 . (3.6)
The results are shown in Figure 7 again showing the alignment of U and ∇ × U breaks down within a certain range of r. We also observe that the magnetic energy in Figure 7 b) almost follows the shear distribution in Figure 7 d). This may at first seem to point to the Omega effect [START_REF] Moffatt | Magnetic field generation in electrically conducting fluids[END_REF], although in our model the relation between shear and flow is more complex so there is no easy way to relate U and B T through S max . For this reason, we cannot give a definite answer regarding the physics involved in the optimal dynamo. What we can say is that shear (as measured by the absolute maximal strain rate) plays an important role in amplifying the magnetic field.

Analysis in spectral space

In spectral space, we report global measures for U and B T as a function of rotationally invariant spherical harmonic degree l in Table 2. We introduce the partial sums for the flow U as

U l = m,n t nml U t nml + s nml U p nml , (3.7)
and similarly for B T l . We notice the l = 1 and l = 3 modes dominate the enstrophy, and 94% of the enstrophy (squared value) comes from the first three spherical harmonic degrees. The dominant flow field with only these three degrees with l max = 3 has a critical magnetic Reynolds number Rm c = 66.29 when the non-dimensional root mean enstrophy is rescaled to 1. This slight increase in the critical point indicates very little influence from higher spherical harmonic degrees in the optimal U . We also notice that the fields have low rms velocities compared to the root mean enstrophy. This indicates that we have strong shear in our optimal flow. In Figure 8, we show power spectra of the optimal root mean enstrophy and magnetic energy as a function of spherical harmonic degree l at Rm c,min . Both models R1 and R2 give the same power spectra. They are rapidly decaying, which confirms that the spatial resolution we have used is sufficient. Alternatively, we can also analyse each spherical harmonic mode with index l, m by summing over indices n. We define the partial sums as

U t m l = n t nml U t nml , U pm l = n s nml U p nml .
(3.8) Power spectra of the optimal enstrophy (∇ × U l ) 2 /V and the magnetic energy (BT l ) 2 as a function of the spherical harmonic degree l.

l = 1 l = 2 l = 3 l = 4 l = 5 (∇ × U l ) 2 /V 0.442 0.130 0.376 0.042 0.006 (U l ) 2 /V 4.40 × 10 -3 2.12 × 10 -3 4.77 × 10 -3 0.42 × 10 -3 0.04 × 10 -3 (BT l ) 2 1.27 0.62 0.16 0.04 0.02
Table 2. The partial enstrophy, kinetic energy and magnetic energy for the first five spherical harmonic degrees l from models R1 and R2, both giving the same value up to three digits.

U p2 3 U p0 1 U t 0 1 U p-2 2 U t -2 2 (∇ × U m l ) 2 /V 0.285 0.231 0.211 0.054 0.046 BT t 1 1 BT p1 1 BT t 1 2 BT p1 2 BT t -1

2

(BT m l ) 2 0.618 0.593 0.384 0.151 0.054 Table 3. Leading components of enstrophy for U m l and magnetic energy for BT m l from model R1 after rotation and reduction. A summation over n has been performed for the comparison.

Since the amplitude of each mode changes with rotation, we must fix an orientation for the optimal U . We rotate the field in spectral space using Wigner D-matrices such that the rotational symmetry axis is along the ẑ direction and the sum of root mean enstrophy of cos(mφ) modes is maximized. We impose a cut-off of absolute value 0.01 for all spectral coefficients after rotation. The leading spherical harmonic modes are given in Table 3. After a further reduction to keep only l 3, the critical magnetic Reynolds number becomes 66.34 when the root mean enstrophy is rescaled to 1, which shows the reduced flow is a good approximation to the full optimum. Capturing 83% of the enstrophy, the optimal flow can be approximated by five spherical harmonic modes:

U ≈ U t 0 1 + U p0 1 + U t -2 2 + U p-2 2 + U p2 3 . (3.9)
Positive m corresponds to cosine modes and negative m corresponds to sine modes. We see the dominant modes in (3.9) mostly describe large scale structures. We also see here an even m symmetry for the optimal flow that matches with the spatial structures in Figure 4 a). Owing to the selection rules of [START_REF] Bullard | Homogeneous dynamos and terrestrial magnetism[END_REF], the magnetic field can only have either all even or all odd m modes when U has only even m modes. In this model, only the odd m modes of B T are present. 3 using model R1. The cross correlation of the perturbed field and the optimal flow field is given by Cp and the measured growth rate is γp. Each square represents an independent perturbation on U .

At the level of individual spherical harmonic modes, we notice from Table 3 that the first spherical harmonic mode (l, m) = (1, 0) after rotation and reduction has approximately the same enstrophy for both poloidal and toroidal components. It turns out this mode has a normalized helicity H 0 1 = 0.95 where

H m l = U m l • ∇ × U m l / U m l 2 (∇ × U m l ) 2 . (3.10)
This corresponds to the strong central flow we have shown in Figures 4 a) and 6 a). Other secondary structures are weaker, so the overall normalized helicity for the optimal flow U is only ∼ 0.65. Helicity is an important measure in fluid dynamics and its significance for dynamo action has been discussed in many places, e.g. review by [START_REF] Moffatt | Magnetic field generation in electrically conducting fluids[END_REF]; [START_REF] Gubbins | Implication of kinematic dynamo studies for the geodynamo[END_REF], as well as specific examples given by [START_REF] Livermore | The role of helicity and stretching in forced kinematic dynamos in a spherical shell[END_REF]. While our optimal dynamo does have large helicity for some dominant components, helicity alone cannot explain all the optimal structures we have. We also find in Table 3 the equatorial dipole component of B T is relatively large.

Perturbation study

In this section, we perform a perturbation study to verify the optimality of the growth rate at Rm c,min . A perturbed field U per is the optimal field U plus a small portion of a randomly generated flow field, subsequently rescaled to have a unit enstrophy. We define a correlation amplitude

C p = U per • U (U per ) 2 U 2 (3.11)
similar to that in [START_REF] Chen | Optimal dynamo action by steady flows confined in a cube[END_REF] but only for the dominant modes l 3, then launch a forward run with U per and measure the perturbed growth rate γ p . The result of multiple random perturbations is shown in Figure 9. The growth rate rapidly decreases as more perturbations are added to U , which indicates that random perturbations quickly reduce the flow's optimality.

Transient growth

In dynamo models, a seed magnetic field can be transiently amplified by a conductive flow field for a short period even when it decays later [START_REF] Livermore | On magnetic energy instability in spherical stationary flows[END_REF], 2006). This can clearly be seen in Figure 3 a) for t < 0.15. Using our optimisation method with very short optimising time windows T , we can maximize the transient growth and then measure the magnetic Reynolds number Rm t,min at which this maximum transient growth is zero. Numerically, using T = 0.001 and ∆t = 10 -5 , we find the minimal transient magnetic Reynolds number Rm t,min = 30.21, (3.12) see Figure 10 a). No flow U ∈ E u is able to admit magnetic solutions with d B 2 /dt 0 at Rm < Rm t,min . Furthermore, this lower bound Rm t,min also applies to any timedependent flow in MHD models that can be approximated by a steady flow U ∈ E u at a given instant t. Thus, the value of Rm t,min provides an important physical measure of the ultimate lower bound for sustaining magnetic energy. Compared with the optimal flow U found with T = 2 in Figure 10 b), we see this magnetic field with maximized transient growth decays more rapidly at later times. This Rm t,min is 53% below the optimal threshold Rm c,min for sustained dynamo action which is quite different to that observed by [START_REF] Willis | Optimization of the magnetic dynamo[END_REF] in the periodic cube. We also find the corresponding optimal flow field U tr at Rm t,min has equatorial reflectional symmetry.

Comparison with theoretical bounds

Our magnetic Reynolds number Rm is not the most conventional in dynamo theory, but is the one that relates to the class of enstrophy normalized fields U ∈ E u . We can define two other magnetic Reynolds numbers using different scales for the flow field:

Rm u = U * L * η * , Rm s = S * max L * 2 η * (3.13)
where U * is the dimensional rms flow velocity, and S * max is the dimensional absolute maximum strain rate. This allows us to convert the results from root mean enstrophy or shear based Rm to kinetic-energy based (Rm u ) and strain rate based (Rm s ) magnetic Reynolds numbers: (3.14) using the non-dimensional optimal flow field U . Both models R1 and R2 give the same value. This yields Rm u = 6.96, Rm s = 215.

Rm u = U 2 /V Rm, Rm s = max V eig (∇U + ∇U ⊺ )/2 Rm = S max Rm,
(3.15) The value of Rm u is very low and close to the values ∼ 12 we have measured in the cube with unit sizes [START_REF] Chen | Optimal dynamo action by steady flows confined in a cube[END_REF]. According to [START_REF] Proctor | Energy requirement for a working dynamo[END_REF], such a low value should be no real surprise since a theoretical lower bound on Rm u does not exist. We can further check how far we are above theoretical lower bounds using definitions other than Rm u . In particular, we choose three bounds to compare with, they are [START_REF] Backus | A class of self-sustaining dissipative spherical dynamos[END_REF] : Rm s > π 2 , [START_REF] Proctor | On Backus' necessary condition for dynamo action in a conducting sphere[END_REF] : Rm s > 12.29, [START_REF] Childress | Thèorie magnétohydrodynamique de l'effet dynamo[END_REF] : Rm U max > π.

(3.16)

The Backus/Proctor bounds use a strain-rate based magnetic Reynolds number, and the Childress bound uses the maximum speed of the flow as a typical scale. Using our conventions, we find Rm s /12.29 ≈ 17 and Rm U max /π ≈ 14, which are clearly above the theoretical limits. This can be attributed to the fact that theoretical bound calculations systematically overestimate the spatial extent of the magnetic fields [START_REF] Herreman | Minimal flow perturbations that trigger kinematic dynamo in shear flows[END_REF].

Besides the bounds listed above, [START_REF] Proctor | Necessary conditions for the magnetohydrodynamic dynamo[END_REF] derived a bound for dynamo action that is immediately applicable to our present results. He derived that Clearly our dynamo operates well above this bound. Proctor speculated that the best lower bound might be ten times larger than the one given here, and indeed we find that the best dynamo is operating at a value of Rm over one hundred times larger.

D * > πη * 2 4L * (3.17

Comparison with other flows

In this section we compare several properties of our optimal flow to some selected models in the literature. We do not include self-consistent dynamos in a spherical shell here due to the complexity of these models (but note that in [START_REF] Christensen | Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields[END_REF], the limit for dynamo action is Rm u ∼ 60 when the radius is used as the length scale, which is about 8 times higher than our optimal dynamo). Instead, we focus on kinematic dynamos driven by two families of large scale flows: the KR type inspired by [START_REF] Kumar | A spectral solution of the magneto-convection equations in spherical geometry[END_REF], and the DJ type inspired by [START_REF] Dudley | Time-dependent kinematic dynamos with stationary flows[END_REF]. The KR flow consists of three parts: differential rotation, meridional circulation and a simplified model for convection columns. The DJ flow has axisymmetric structures. The boundary conditions are not limited to no-slip for these models. The precise definition of the flows abbreviated as KR, GKR, STW, DJ t 1 s 1 , DJ t 1 s 2 , DJ t 2 s 2 , MDJ t 1 s 2 is given in Appendix B. We could not include the optimized DJ flow by [START_REF] Holme | Optimised axially-symmetric kinematic dynamos[END_REF] in our comparison table due to the unknown explicit field representation.

We first compare the critical magnetic Reynolds numbers for the onset of dynamo action. All the studies we compare with in spheres have used L * as the dimensional radius, but there is little homogeneity in what the scale for the flow field should be (e.g. maximal dimensional speed, rms dimensional speed). In order to compare critical magnetic Reynolds numbers of different dynamo studies, we need to rescale all results to a consistent definition. Denoting U the originally reported flow field that was nondimensionalized using a scale other than [U ] = ω * L * , we systematically convert Rm to Rm u , Rm and Rm s respectively:

Rm u = U 2 /V Rm, Rm = (∇ × U ) 2 /V Rm, Rm s = max V eig (∇ U + ∇ U ⊺ )/2 Rm = S max Rm. (3.20)
The comparison of critical magnetic Reynolds numbers and the rescaling factors in (3.20) are shown in Table 4. We see our optimal flow U has improved the existing lower bound on the critical Rm of DJ t 2 s 2 by at least a factor of 3. The Rm u for our optimal flow U is also significantly lower than the others. The trade-off of this is that the maximal strain rate remains comparable to DJ flows (we compared S max and S max ( U [START_REF] Livermore | On magnetic energy instability in spherical stationary flows[END_REF], shown in Table 5.

2 /V ) 1/2
Additionally, we also compare important properties of the optimal flow such as helicity, energy ratio and maximum speed in Table 6. We use the poloidal/toroidal kinetic energy and enstrophy given by

E P /E T = (U p ) 2 / (U t ) 2 , ω P 2 /ω T 2 = (∇ × U p ) 2 / (∇ × U t ) 2 , (3.21)
where the superscript p and t denote the poloidal and toroidal component of flow field respectively. The choice of scale [U ] does not influence the energy ratios, so we use the same notation for all models. The total helicity in the sphere and the hemispherical helicity for flows U are given by

H V = 3 4π U • ∇ × U , H V 2 = 3 2π 1 0 2π 0 π 2 0 U • ∇ × U r 2 sin θ dθdφdr. (3.22)
We use the same formula without the hat symbol for the flow field that has been nondimensionalized using scale [U ] = ω * L * . From Table 6 it seems that there is no clear pattern on these measured values across different dynamo models. We observe that for our optimal flow the poloidal enstrophy to toroidal enstrophy ratio is roughly the same as poloidal kinetic energy to toroidal kinetic energy ratio.

Comparison with the cubic model

It is of interest to compare our present results to our previous results in a cube [START_REF] Chen | Optimal dynamo action by steady flows confined in a cube[END_REF]. When mixed boundary conditions were used, the best dynamo had Rm c,min = 7.5π 2 where the length scale was based on the length of the side of the cube. For the sphere the choice of a length scale is clear (it is the diameter), but for the cube there is no obvious analogous scale. If we just consider all chords passing through the centre of a unit cube, it is clear that the lengths l vary as one moves from a chord perpendicular to two faces to one joining opposite corners, such that l is bounded by 1 l √ 3. There is probably no fully satisfactory way of defining the relevant length scale in a cube, but we can derive an objective choice. We let E{x} denote the expected value of x and consider all chords present in both the sphere and the cube. For the sphere of radius L it is clear that E{l 2 } = 4L 2 = (E{l}) 2 , whereas in the cube, we find

E{L 2 } = 5/3 (3.23)
with the ancillary result (not needed) of E{L} ≈ 1.28. Using the value (3.23) for the length scale gives Rm c,min = 123 in the cube with mixed boundary conditions, whereas Rm c,min = 197 with homogeneous (either pseudo-vacuum or superconducting) boundary conditions. These compare to the values in the sphere we report here when converted using E{l 2 } = 4L 2 of Rm c,min = 258 for insulating boundary conditions. A complementary way of comparing is to imagine that we are given a unit volume of electrically conducting fluid and consider the two possible geometrical arrangements at our disposal: a cube and a sphere. Obviously the size of the sphere (using radius R) is 4π/3R 3 = 1, thus R = (3/(4π)) 1/3 ≈ 0.62. We choose a common length scale, which is unity for the cube and for the sphere must be measured in these units. We previously referred our Rm to the radius, but the present lengthscale inflates the length units by (0.62) -1 .When we reinterpret our present value for Rm on this length scale we have a new value for Rm of 64.45/0.62 2 ≈ 168, now above all the cubic results. Probably we have taken this analysis as far as is sensible, but we speculate that the extra roughness of the boundaries present in the cube might be responsible for additional shearing in the optimal flow that then appears to be a better dynamo.

Conclusions

In this study, we have found the optimal kinematic dynamo in a sphere with no-slip and insulating boundary conditions. The enstrophy based minimal critical Rm c,min for the onset of a dynamo is 64.45. Compared to other known dynamo models, our optimal dynamo has lowered the critical Rm at least by a factor of 3. The rms speed is much lower compared to the global shear magnitude or root mean enstrophy in the system. The optimal flow field at Rm c,min has a rotational symmetry of order 2 and a very concentrated helical structure near the centre. This indicates that a localised helix plus secondary twofold spirals are favourable for the onset of a kinematic dynamo in a sphere. The dominant spherical harmonic modes are (l, |m|) = (1, 0), (2, 2), (3, 2), when the rotational symmetry axis is aligned with ẑ axis. The optimal flow structures near the boundary do not play a significant role. Therefore, we expect that the use of other boundary conditions for the flow, such as free-slip boundary conditions, will not lower Rm c,min by much. The magnetic field at the dynamo onset has mainly dipolar structures. With respect to the rotational symmetry axis of U , the fastest growing magnetic field eigenmode has only odd m modes. We also find the minimal critical magnetic Reynolds number for transient growth Rm t,min = 30.21, indicating any time-dependent flows cannot lower the critical magnetic Reynolds number below 30.21.

Our optimisation method is designed to identify the best dynamo, but it does not answer the question "why are these flows so efficient in generating dynamos?" quantitatively. Obtaining and analysing the optimal flow as we have done here is just a first step, while providing understanding what makes them optimal physically is a much more relevant but rather difficult second step. In the supplementary material of this article, we have provided both the spectral coefficients and field values on grids for the optimal flow field and the associated magnetic eigenvector for possible follow-up studies. Since our polynomial basis is not standard, we also include a Mathematica file that allows one to generate the optimal flow in physical space and on arbitrarily fine grids.

We plan to extend this study in two parallel directions. One is to study how different boundary conditions on U and B affect the optimum, the other is to study the optimum within a restricted subclass of the present symmetry group O(3), e.g. to find the optimal kinematic dynamo generated by axisymmetric flows. It is also feasible to study the optimum within geophysically interesting classes of flows, such as Rossby waves or inertial waves.

ǫ 1 = 0.03, ǫ 2 = ǫ 3 = 0.04, K = 4.4 is a normalization factor that has been used in Love & Gubbins (1996a) to get unit rms speed. The factors of the square root of π are related to the normalization of the spherical harmonics.

• GKR flow This is a dynamo studied by Gubbins et al. (2000a) and [START_REF] Li | An optimal galerkin scheme to solve the kinematic dynamo eigenvalue problem in a full sphere[END_REF]. The flow field is a generalized KR type flow given by t 1 0 (r) = ǫ 0 r 2 (1 -r 2 ), s 2 0 (r) = ǫ 1 r 6 (1 -r 2 ) 3 , s 2 2 (r) = ǫ 2 r 4 (1 -r 2 ) 2 sin(3πr), s 2 -2 (r) = ǫ 3 r 4 (1 -r 2 ) 2 cos(3πr), (B 3) where ǫ 0 = 6.421324129676529, ǫ 1 = -4.632495546344649, ǫ 2 = ǫ 3 = 1.414456257660752.

• STW flow This is dynamo solution that satisfies the thermal wind equation found by [START_REF] Sarson | Kinematic dynamos driven by thermal wind flows[END_REF]. The flow field is t 1 0 (r) = K 4π/3 r 2 (1 -r 2 ), s 2 -2 (r) = K 4π/5 ǫ 2 r 4 (1 -r 2 ) 2 cos(3πr), s 2 2 (r) = K 4π/5 ǫ 3 r 4 (1 -r 2 ) 2 sin(3πr), • DJ flows This is a group of three simple axisymmetric flows found by [START_REF] Dudley | Time-dependent kinematic dynamos with stationary flows[END_REF]: DJ t 1 s 1 : t 1 0 (r) = 4π/3 r sin(πr), s 1 0 (r) = 0.17 4π/3 r sin(πr), DJ t 1 s 2 : t 1 0 (r) = 4π/3 r sin(πr), s 2 0 (r) = 0.13 4π/5 r 2 sin(πr), DJ t 2 s 2 : t 2 0 (r) = 4π/5 r 2 sin(πr), s 2 0 (r) = 0.14 4π/5 r 2 sin(πr). (B 5)

• MDJ t 1 s 2 flow A modified Dudley James flow is discussed by [START_REF] Livermore | On magnetic energy instability in spherical stationary flows[END_REF], [START_REF] Li | An optimal galerkin scheme to solve the kinematic dynamo eigenvalue problem in a full sphere[END_REF]. The flow field is given by

t 1 0 (r) = a 1 r 2 (1 -r 2 ), s 2 0 (r) = a 2 r 3 (1 -r 2 ) 2 , ( B 6)
where a 1 = 1.193271237996972 and a 2 = 8.107929179422066. In our Galerkin basis, it can be expressed as a multiple of toroidal basis f l n (r), for n = 1, l = 1 and poloidal basis g l n (r), for n = 1, l = 2.

Figure 1 .

 1 Figure 1. Verification of the adjoint model. Parameters used in the test: Rm = 60, T = 0.1, ∆t = 10 -4 , ε = 10 -4 , and (lmax, nmax) = (24, 12). a) Toroidal flow coefficients. b) Poloidal flow coefficients. There is almost exact overlay of the curves depicting the variations in L obtained through two methods.

Figure 2 .Figure 3

 23 Figure 2. Main figure: the growth rate as a function of Rm for multiple optimisation runs. Inset figure: the near-zero growth rate as a function of Rm coloured by different resolutions. Orange stars: (lmax, nmax) = (16, 8). Green crosses: (lmax, nmax) = (24, 12). Purple diamond: (lmax, nmax) = (24, 24).

Figure 4 .

 4 Figure 4. Spatial structures at Rmc,min from model R1: a) Streamlines of U , b) Streamlines of BT . Field lines are coloured by local field intensity (red = intense, blue = weak) and indicate a localized structure near the centre. Structures with very small field intensity are not shown here. The plotted domain is r ∈ [0, 1].

Figure 5 .

 5 Figure 5. a) The velocity field direction (red) and vorticity field direction (black) are nearly parallel within r 0.1. b) Streamlines of the velocity field (orange), the vorticity field (black) and the magnetic field BT (blue) within r 1. The magnetic field is strongest next to the helical fluid motion.

Figure 7 .

 7 Figure 7. Radial distributions of the optimal solution. a) Kinetic energy. b) Magnetic energy. c) Helicity. d) Absolute maximal strain rate over all angles θ, φ.

  Figure8. Power spectra of the optimal enstrophy (∇ × U l ) 2 /V and the magnetic energy (BT l ) 2 as a function of the spherical harmonic degree l.

Figure 9 .

 9 Figure9. Perturbation study of the optimal flow U with l 3 using model R1. The cross correlation of the perturbed field and the optimal flow field is given by Cp and the measured growth rate is γp. Each square represents an independent perturbation on U .

Figure 10 .

 10 Figure 10. The optimized transient growth. a) Optimized magnetic energy growth for Rm ∈ [30.20, 30.22] with optimising time T = 0.001, resolution (lmax, nmax) = (24, 12). b) Comparison of optimized transient growth (red solid line) and R1 (blue dashed line) at the same Rmt,min = 30.21.

  ) for a sphere of radius L * , where D * = e * ij e * ij dV * is proportional to viscous dissipation and e * ij is the dimensional strain rate tensor. Translating this to our notation, the no-slip conditions, S * = ω * , and we find immediately that

  0.04 and ǫ 3 = -0.04. Here K = 4.2449 is a normalization factor so the rms speed of the flow is 1.

  ). It is not yet clear how to reduce local shear and global shear at the same time. As for the transient growth, our ultimate lower bound Rm t,min = 30.21 is at least 2 times lower than the critical transient magnetic Reynolds number (Rm t ) for selected flows from Livermore

Table 4 .

 4 A comparison of critical magnetic Reynolds numbers between our optimised flow (top row) and those reported in the literature (bottom section), the rescaling factor rms (• • • ) is equivalent to ( (• • • ) 2 /V ) 1/2 ,see (3.14), (3.20) and Appendix B for more definitions.

			Rmu Rm Rms rms U	Smax
	U		7	64.45 215	0.108	3.33
		Rm Rmu Rm Rms	rms U	Smax	rms ∇ × U
	KR	890 890 4486 5002	1	5.62	5.04
	GKR	131 131 1327 1451	1	11.08	10.13
	STW	479 479 2936 3741	1	7.81	6.13
	DJ t1s1	155 104 476 274	0.671	1.77	3.07
	DJ t1s2	95	60	310 233	0.635	2.45	3.26
	DJ t2s2	54	32	193 153	0.587	2.84	3.57
	MDJ t1s2 55	55	300 296	1	5.38	5.46
				U tr	KR STW MDJ t1s2
		Rmt	30.21	80	83	62

Table 5 .

 5 The critical transient magnetic Reynolds number for the optimal transient flow U tr and selected flows from[START_REF] Livermore | On magnetic energy instability in spherical stationary flows[END_REF], converted using ( (∇ × U ) 2 /V ) 1/2 as the rescaling factor.

		H V H V 2	max V	|U | ωP 2 /ωT 2 EP /ET
	U	0.07 0.08 ‡	0.69	1.885	1.824
		H V H V 2	max V	| U | ωP 2 /ωT 2 EP /ET
	KR	0	0.210	1.96	0.149	0.017
	GKR	0	2.078	2.58	8.117	1
	STW	0	0.362	2.38	0.088	0.015
	DJ t1s1	1.969 1.969	1.12	0.630	0.592
	DJ t1s2	0	1.380	1.07	0.838	0.425
	DJ t2s2	1.992 1.992	0.94	0.791	0.679
	MDJ t1s2	0	2.906	1.64	0.662	0.254

Table 6 .

 6 Properties of the flow fields that can generate a dynamo (see Appendix B for the definition of flows). ‡ : We choose the hemisphere with the larger value of helicity when the symmetry axis is aligned with the ẑ axis.

† Email address for correspondence: ajackson@ethz.ch

Appendix A. Orthonormal bases of toroidal/poloidal vector fields

In this section we give the definitions of our Galerkin bases and show how to project a vector field onto these bases.

A.1. The Poloidal-toroidal expansions

The poloidal-toroidal decomposition for the flow field is

where f l n , g l n are radial functions and Y m l are the spherical harmonics. The field decomposition for B follows similar definitions. There are in total 2 n max l max (l max + 2) spectral coefficients for each vector field. For input and output, these coefficients are stored in memory as a one dimensional array. In the optimisation process, the memory is distributed to at least l max number of cores. In this model, we use fully normalized real spherical harmonics. The spherical harmonics are defined as

where P m l (cos θ) are the associated Legendre functions, defined as

where Pn (x) is a Legendre polynomial of x, and

(l+|m|)! is a normalization factor such that we have the orthonormality condition for the spherical harmonics:

The derivatives of spherical harmonics satisfy

where âm l =

(l+|m|)(l-|m|) (2l+1)(2l-1) . They are chosen in this way such that we have a simple formula for taking the derivatives, an alternative way would be to use only positive m, but then we need to have separate expression for sine and cosine functions.

A.2. Orthogonal Galerkin polynomials

The optimisation scheme relies on a convenient choice of orthogonal polynomials. These polynomials are constructed by combining Jacobi polynomials into a Galerkin basis that satisfies the required boundary conditions using the techniques pioneered by [START_REF] Livermore | Quasi-L p norm orthogonal Galerkin expansions in sums of Jacobi polynomials[END_REF] and [START_REF] Livermore | Galerkin orthogonal polynomials[END_REF]. The key idea is that by using a Gram-Schmidt procedure, a family of orthogonal Galerkin polynomials can be iteratively constructed. However, as these authors showed, it is possible to explicitly write down the structure of Galerkin polynomial of arbitrary index as the terse sum of Jacobi polynomials. This allows us to deduce that asymptotically the recombined Jacobi polynomials have the same optimal properties as the Jacobi polynomials themselves. Using field expansions (A 1), the orthonormality conditions (2.36) in this model then requires

for radial functions, where the radial inner products are defined as

and

and the last term in (A 10) derives from the fact that the inner product of magnetic field is defined over all space; it handles the vacuum contribution from the poloidal magnetic field. For simplicity, we use the same notation for the poloidal flow field under the condition: g l n (1) = 0.

A.3. Projection of non-linear terms and the diffusion term

To solve the induction equation forward in time, we need to project the induction term onto the basis of B:

The cross product

is evaluated on the grid points. This transformation from spectral space to physical space uses an inverse Fourier transform in φ and a matrix multiplication transform (MMT) in r and θ. The MMT requires us to store grid values of each radial expansion function and each associated Legendre functions in memory. Along r and θ, we use N r and N θ Gaussian quadrature points and along φ we use N φ Fourier grid points. The number of spatial grid points for resolution (l max , n max ) are

where ⌈ ⌉ is the ceiling function. We then calculate their projection on the spherical harmonic basis:

Finally, we extract the spectral coefficients by an integration along radial direction. For the poloidal components,

We have used partial integration, no-slip and electrically insulating boundary condition to derive (A 16). For the toroidal components,

The diffusion does not involve non-linear terms and can be directly calculated in spectral space. This operation is represented by diffusion matrices given by

For the backward time integration of the adjoint problem, the adjoint induction term needs to be projected onto the basis of B. In practice, we evaluate ∇ × B † directly on the physical grids using an inverse Fourier transform along φ and MMT along r and θ.

We then calculate the terms

on the gridpoints and project onto (r, l, m) space to get [H µ ] m l (r, l, m) as defined in (A 15) using a Fourier transform and MMT. Using the recurrence rules of the spherical harmonics, we get

which we can evaluate using the Gaussian quadrature rule.

In the velocity field update, we need to calculate the integral term I t nml (t i ) at each discrete time t i . This projection is structurally similar to the ones we have in the adjoint induction term, so we follow almost the same method. Let B ti = B(x, t i ), we calculate

on gridpoints, then project it onto (r, l, m) space as [G µ ] m l (r, l, m). Finally we evaluate

, (A 25) using the Gaussian quadrature rule.

Appendix B. Flow fields for comparison

The description of the flow fields used for comparison in Table 4 and Table 6 is given below where we have used the notation

that is compatible with fully normalized spherical harmonics in Appendix A.

• KR flow [START_REF] Kumar | A spectral solution of the magneto-convection equations in spherical geometry[END_REF] defined a flow field that contains both azimuthal differential rotation and meridional circulation. The flow is given by t 1 0 (r) = K 4π/3 r 2 (1 -r 2 ), s 2 0 (r) = K 4π/5 ǫ 1 r 6 (1 -r 2 ) 3 , s 2 -2 (r) = K 4π/5 ǫ 2 r 4 (1 -r 2 ) 2 cos(3πr), s 2 2 (r) = K 4π/5 ǫ 3 r 4 (1 -r 2 ) 2 sin(3πr), (B 2)