# Magnetic energy cascade in spherical geometry. I. the stellar convective dynamo case

A. Strugarek, A.S. Brun, S. Mathis, Y. Sarazin

► To cite this version:

A. Strugarek, A.S. Brun, S. Mathis, Y. Sarazin. Magnetic energy cascade in spherical geometry. I. the stellar convective dynamo case. The Astrophysical Journal, 2013, 764 (2), pp.189. 10.1088/0004-637X/764/2/189 . cea-00827693

HAL Id: cea-00827693

https://cea.hal.science/cea-00827693

Submitted on 12 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
MAGNETIC ENERGY CASCADE IN SPHERICAL GEOMETRY. I.
THE STELLAR CONVECTIVE DYNAMO CASE

A. STRUGAREK<fnref n="1" />, A. S. BRUN<fnref n="1" />, S. MATHIS<fnref n="1" />, AND Y. SARAZIN<fnref n="2" />

Received 2012 May 10; accepted 2013 January 8; published 2013 February 6

ABSTRACT

We present a method to characterize the spectral transfers of magnetic energy between scales in simulations of stellar convective dynamos. The full triadic transfer functions are computed thanks to analytical coupling relations of spherical harmonics based on the Clebsch–Gordan coefficients. The method is applied to mean field αΩ dynamo models as benchmark tests. From a physical standpoint, the decomposition of the dynamo field into primary and secondary dynamo families proves very instructive in the αΩ case. The same method is then applied to a fully turbulent dynamo in a solar convection zone, modeled with the three-dimensional MHD Anelastic Spherical Harmonics code. The initial growth of the magnetic energy spectrum is shown to be non-local. It mainly reproduces the kinetic energy spectrum of convection at intermediate scales. During the saturation phase, two kinds of direct magnetic energy cascades are observed in regions encompassing the smallest scales involved in the simulation. The first cascade is obtained through the shearing of the magnetic field by the large-scale differential rotation that effectively cascades magnetic energy. The second is a generalized cascade that involves a range of local magnetic and velocity scales. Non-local transfers appear to be significant, such that the net transfers cannot be reduced to the dynamics of a small set of modes. The saturation of the large-scale axisymmetric dipole and quadrupole is detailed. In particular, the dipole is saturated by a non-local interaction involving the most energetic scale of the magnetic energy spectrum, which points to the importance of the magnetic Prandtl number for large-scale dynamos.

Key words: dynamo – magnetohydrodynamics (MHD) – stars: magnetic field – turbulence

Online-only material: color figures

# 1. INTRODUCTION

Magnetic fields are observed in astrophysical bodies in a broad range of scales, from the object scale to the smallest dissipative scales (Donati & Landstreet 2009). The origin of such fields is, in most cases, due to a hydromagnetic dynamo process. Recent developments in dynamo theory led to a distinction between large-scale and small-scale dynamos (Cattaneo & Hughes 2001). The large-scale dynamos produce magnetic fields at a larger scale than the largest velocity scale (or, the largest driving scale), while small-scale dynamos generate magnetic fields at all scales smaller or equal to the driving scales (Tobias et al. 2011). Large-scale dynamos also sometimes refer to dynamos that develop a large-scale magnetic field in super-equipartition with large-scale kinetic energy (Olson et al. 1999). In that case, small-scale dynamos refer to those that develop a spectrum peaked at small scales. Both dynamos generally act together, like in the Sun, where we observe both large-scale, intense, global magnetic fields (Schrijver & DeRosa 2003; DeRosa et al. 2011, 2012) and small-scale magnetic fields (Hagenaar et al. 2003; Centeno et al. 2007). In the case where multiple scales also coexist in the velocity field, special care is needed to isolate “small”- and “large”-scale tendencies (Tobias & Cattaneo 2008a).

Recent developments in numerical simulations in three-
dimensional (3D) spherical geometry allow us to model fully
nonlinear dynamos in stars, involving a broad range of scales
(Brun et al. 2004; Browning 2008; Brown et al. 2010; Racine
et al. 2011; Käpylä et al. 2012). The flow scales of a stellar
convection zone extend from the large-scale differential rotation
down to the smallest convective scales. In order to properly
characterize such dynamos, one may use specific methods to

tackle the multi-scale aspect of the problem. The principal tools
that have been used in the literature are spectral decomposition
(Frick & Sokoloff 1998; Dar et al. 2001) and wavelet analysis
(Farge 1992). In this paper, we choose to use spherical harmonics
decomposition (which is adapted to the spherical geometry of
stars; Bullard & Gellman 1954) to develop a spectral analysis
of energy transfers in the frame of dynamo theory.

Mainly used to study turbulence (Frisch 1995; Debliquy et al.
2005; Lesieur 2008; Alexakis et al. 2005), spectral analysis is
also a useful tool to characterize magnetohydrodynamic (MHD)
processes such as dynamos (Biskamp 1993; Blackman & Bran-
denburg 2002; Mininni et al. 2005; Livermore et al. 2010) or the
magneto-rotational instability (see Lesur & Longaretti 2011).
Understanding spectral energy transfers between scales in such
processes may reinforce our ability to characterize nonlinear
MHD phenomena. The shell-to-shell or mode-to-mode meth-
ods have been recently and extensively used in the context of
MHD turbulence. Indeed, the classical Kolmogorov approach
to turbulence must be adapted to the MHD case, since the mag-
netic field induces an anisotropy that has to be taken into account
(Iroshnikov 1964; Kraichnan 1965; Biskamp 1993; Goldreich
& Sridhar 1995). Depending on the dimensionality of the prob-
lem, spectrum slopes are often understood to result from local
(direct or indirect) transfers of energy, referred to as *cascades*
(Biskamp 1993; Maron et al. 2004). However, it was found that
non-local interactions in MHD turbulence may also contribute
importantly to the build-up of the spectrum (Schilling & Zhou
2002; Aluie & Eyink 2010). The directions and localizations of
energy transfers are then less obvious to identify, and studies
dedicated to transfer processes in spectral space are essential
to properly understand spectrum slopes in MHD (Politano &
Pouquet 1998; Boldyrev et al. 2009; Pouquet et al. 2011).

<footnote n="1">¹ Laboratoire AIM Paris-Saclay, CEA/Irfu Université Paris-Diderot CNRS/INSU, F-91191 Gif-sur-Yvette, France; antoine.strugarek@cea.fr</footnote>

<footnote n="2">² CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France</footnote>
In the past, spectral analyses have mainly been used with Fourier spectral decomposition, generally in Cartesian coordinates and periodic parallelepipedic boxes. The Fourier decomposition is indeed a natural way to understand spectra, since the Fourier wave numbers represent the inverse of a spatial scale. More recently, Hughes & Proctor (2012) used a Fourier spectral analysis to study the influence of large-scale sheared flows on local convective dynamos. They show that the dynamo process depends on a broad range of scales in this case. We also study the dynamo process in the present paper, though in spherical geometry with a self-consistently generated large-scale sheared flow (the differential rotation). In the case of stars or planets, the spherical geometry of the object makes the spherical harmonics basis much more adapted to the spectral analysis (Bullard & Gellman 1954). For example, Ivers & Phillips (2008) wrote the decomposition of the MHD equations onto spherical harmonics in the framework of geodynamics. They were able to analytically express the nonlinear terms by calculating the coupling between the spherical harmonics with the Clebsch-Gordan coefficients (see also Mathis & Zahn 2005, in a stellar context). The spherical harmonics decomposition was also used by Livermore et al. (2010) to develop a spectral analysis that is similar to the one we present in this paper. They used it to identify the spectral interactions leading to a different saturation level of the large-scale magnetic field in kinematic and nonlinear forced dynamos. In the latter case, they observe a significant reorganization of the magnetic field such that a strong large-scale magnetic field can emerge. In addition to a particular geometry, the choice of a certain basis for the spectral analysis may be motivated by the presence of anisotropy (e.g., between the vertical and horizontal directions), which is poorly described by the classical Fourier decomposition (e.g., see Rincon 2006, in the case of turbulent convection).

The spectral interactions in MHD involve triads coupling, meaning that two modes interact to impact a third one through a triangulation rule. Depending on the ideal MHD invariant considered, these kinds of interactions involve couplings between the velocity and the magnetic field, impacting the magnetic or the velocity field. Shell-to-shell methods generally only consider dual interactions, creating an ambiguity over the medium (third component) of the triadic interaction (Verma et al. 2005). In order to cope with this ambiguity, other studies (e.g., Schilling & Zhou 2002) made use of the eddy-damped quasi-neutral Markovian two-point statistical closure (Frisch et al. 1975; Pouquet et al. 1976) to get an analytical expression of the triadic interactions. We point out that such methods are relevant for, e.g., developing subgrid-scale models for large eddy simulations. Tobias & Cattaneo (2008b) demonstrated that such truncated methods often badly describe coherent turbulent structures in flows, which are thought to be responsible for the generation of large-scale fields in dynamos. This limitation is relevant for reduced spectral models, which aim to reproduce the full turbulent behavior with a reduced number of modes. However, here we directly calculate the full triadic shell-to-shell interactions of all of the scales included in our simulations, (i.e., we do not use any specific closure in spectral space to compute the full triadic interactions).

We here applied our spherical harmonic based method to a numerical simulation of a stellar convection by considering the three spectral components of the triadic interactions. The originality of the method we develop in the present work resides in the facts that (1) we explicitly decompose the spectral interactions for both of the magnetic and velocity fields, (2) we

calculate explicitly all of the coupling coefficients between those fields, and (3) we use it to study dynamo action in a solar-like turbulent convection zone that possesses self-consistent large-scale flows (differential rotation, meridional circulation, ...) as well as a broad range of turbulent scales.

In Section 2, we present the set of MHD equations we will use, derive from them the spectral evolution equation for the magnetic energy in the spherical harmonics formalism, and analytically validate our method. A toy model of an axisymmetric $\alpha - \Omega$ dynamo is analyzed with our spectral method in Section 3. In Section 4, we apply our method to study nonlinear dynamo action in a numerical simulation of a solar convective zone. Finally, conclusions and perspectives are given in Section 5.

## 2. MAGNETIC ENERGY EVOLUTION EQUATION

### 2.1. Main Equations in Physical Space

We use the well-tested Anelastic Spherical Harmonics (ASH) code, which models turbulent stellar convection zones (Clune et al. 1999; Jones et al. 2011). It solves the following 3D MHD set of equations (see Brun et al. 2004) in the anelastic approximation, in a reference frame rotating at the angular velocity $\Omega_0 = \Omega_0 e_z$ (where $e_z$ is the Cartesian vertical axis):

$$ \nabla \cdot (\bar{\rho} \mathbf{U}) = 0, \qquad (1) $$

$$ \nabla \cdot \mathbf{B} = 0, \qquad (2) $$

$$ \bar{\rho} [\partial_t \mathbf{U} + (\mathbf{U} \cdot \nabla) \mathbf{U} + 2\Omega_0 \times \mathbf{U}] = -\nabla P + \rho \mathbf{g} + \frac{1}{4\pi} (\nabla \times \mathbf{B}) \times \mathbf{B} - \nabla \cdot \mathbf{D} - [\nabla \bar{P} - \bar{\rho} \mathbf{g}], \quad (3) $$

$$ \bar{\rho} \bar{T} [\partial_t S + \mathbf{U} \cdot \nabla (\bar{S} + S)] = \nabla \cdot [\kappa_r \bar{\rho} c_p \nabla (\bar{T} + T) + \kappa_0 \bar{\rho} \bar{T} \nabla \bar{S} + \kappa \bar{\rho} \bar{T} \nabla S] + \frac{4\pi \eta}{c^2} \mathbf{J}^2 + 2\bar{\rho} v \left[ e_{ij} e_{ij} - \frac{1}{3} (\nabla \cdot \mathbf{U})^2 \right], \quad (4) $$

$$ \partial_t \mathbf{B} = \nabla \times (\mathbf{U} \times \mathbf{B}) - \nabla \times (\eta \nabla \times \mathbf{B}), \quad (5) $$

where the spherically symmetric background thermodynamical state is denoted by bars (fluctuations with respect to the background state are denoted without bars), $\mathbf{v}$ is the local velocity, $\kappa_r$ is the radiative diffusivity, and $\kappa$, $v$, and $\eta$ are, respectively, the effective thermal diffusivity, the eddy viscosity, and the magnetic diffusivity. The thermal diffusion coefficient $\kappa_0$ is applied at the top of the modeled convective zone (where convective motions vanish) to ensure the heat transport through the upper surface. $\mathbf{J} = (c/4\pi) \nabla \times \mathbf{B}$ is the current density, and the viscous stress tensor $\mathbf{D}$ is defined by

$$ D_{ij} = -2\bar{\rho}v \left[ e_{ij} - \frac{1}{3} (\nabla \cdot \mathbf{U}) \delta_{ij} \right], \quad (6) $$

where $e_{ij}$ is the strain rate tensor, and $\delta_{ij}$ is the Kronecker symbol. The system is closed by using the linearized ideal gas law:

$$ \frac{\rho}{\bar{\rho}} = \frac{P}{\bar{P}} - \frac{T}{\bar{T}} = \frac{P}{\gamma \bar{P}} - \frac{S}{c_p} \quad (7) $$

with $c_p$ as the specific heat at constant pressure and $\gamma$ the adiabatic exponent. The vectorial fields are decomposed into
poloidal and toroidal components:

$$ \mathbf{B}(r, \theta, \varphi) = \nabla \times \nabla \times [C(r, \theta, \varphi)\mathbf{e}_r] + \nabla \times [A(r, \theta, \varphi)\mathbf{e}_r], \quad (8) $$

$$ \bar{\rho}(r)\mathbf{U}(r, \theta, \varphi) = \nabla \times \nabla \times [W(r, \theta, \varphi)\mathbf{e}_r] + \nabla \times [Z(r, \theta, \varphi)\mathbf{e}_r], \quad (9) $$

where $(\mathbf{e}_r, \mathbf{e}_\theta, \mathbf{e}_\varphi)$ are the unit vectors in spherical coordinates. All of the quantities are time-dependent. This decomposition ensures numerically that both the magnetic field and the mass flux remain divergenceless up to the machine precision.

A potential match of the magnetic field ($\nabla \times \mathbf{B} = 0$) is applied at both the bottom and top radial boundaries. For the convective dynamo case (Section 4), the boundary conditions for the velocity are impenetrable and stress-free. A latitudinal entropy gradient is imposed at the bottom (as in Miesch et al. 2006), and we fix a constant entropy gradient at the top of the domain.

## 2.2. Magnetic Energy Transfer Functions

### 2.2.1. The Formalism

In this section, we present a method to obtain a spectral (in the sense of the spherical harmonics) evolution equation for the magnetic energy, starting from the induction equation (5). In order to deal with vectorial fields and spherical harmonics (see Equation (A2)), it is practical to define the vectorial spherical harmonics basis (Rieutord 1987; Mathis & Zahn 2005):

$$ \left\{ \begin{aligned} \mathbf{R}_l^m(\theta, \varphi) &= Y_l^m(\theta, \varphi)\mathbf{e}_r \\ \mathbf{S}_l^m(\theta, \varphi) &= \nabla_{\perp} Y_l^m = \partial_{\theta} Y_l^m \mathbf{e}_{\theta} + \frac{1}{\sin\theta} \partial_{\varphi} Y_l^m \mathbf{e}_{\varphi} \\ \mathbf{T}_l^m(\theta, \varphi) &= \nabla_{\perp} \times \mathbf{R}_l^m = \frac{1}{\sin\theta} \partial_{\varphi} Y_l^m \mathbf{e}_{\theta} - \partial_{\theta} Y_l^m \mathbf{e}_{\varphi} \end{aligned} \right. . \quad (10) $$

It is an orthogonal basis for the scalar product $\int_S d\Omega$, where $S$ is a spherical surface and $d\Omega = \sin\theta d\theta d\varphi$ is the associated infinitesimal solid angle. The mode numbers $m$ and $l$ are the azimuthal wave number and the spherical harmonic degree (which characterize their latitudinal variations). The general properties of this basis may be found in Appendix A.1. The two main vectorial fields that appear in the induction equation (5) are the magnetic and the velocity fields. We want to project those fields onto the vectorial basis (Equation (10)), using the decompositions (Equations (8)and (9)). Fortunately, the curl of a vector is a linear operation that can be expressed very easily in the vectorial spherical harmonics basis (see Equation (A9)). We obtain

$$ \mathbf{B}(r, \theta, \varphi) = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \left\{ \frac{l(l+1)}{r^2} C_m^l(r) \mathbf{R}_l^m + \frac{1}{r} \partial_r C_m^l(r) \mathbf{S}_l^m + \frac{A_m^l(r)}{r} \mathbf{T}_l^m \right\}, \quad (11) $$

$$ \bar{\rho}(r)\mathbf{U}(r, \theta, \varphi) = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \left\{ \frac{l(l+1)}{r^2} W_m^l(r) \mathbf{R}_l^m + \frac{1}{r} \partial_r W_m^l(r) \mathbf{S}_l^m + \frac{Z_m^l(r)}{r} \mathbf{T}_l^m \right\}, \quad (12) $$

where we have projected the toroidal and poloidal components of the fields ($A, C, \bar{Z}$, and $W$) onto the scalar spherical harmonics

basis. We see that in Equation (11), the poloidal $C$ and toroidal $A$ components of $\mathbf{B}$ are, respectively, projected onto $(\mathbf{R}, \mathbf{S})$ and $\mathbf{T}$. Consequently, in the remainder of this paper, the projection of any vectorial field onto $(\mathbf{R}, \mathbf{S})$ will be referred to as poloidal, and the projection onto $\mathbf{T}$ as toroidal.

### 2.2.2. Shell to Shell Analysis

To study the transfers of energy between scales on a spherical surface, we distinguish between the different scales of the axisymmetric ($m=0$) and non-axisymmetric physical fields by defining shells $L^0$ and $L^*$ as follows:

$$ \mathbf{X}_L^0 = \mathcal{A}_0^l \mathbf{R}_l^0 + \mathcal{B}_0^l \mathbf{S}_l^0 + C_0^l \mathbf{T}_l^0, \quad (13) $$

$$ \mathbf{X}_L^* = \sum_{\substack{-l \le m \le l \\ m \ne 0}} \{ \mathcal{A}_m^l \mathbf{R}_m^l + \mathcal{B}_m^l \mathbf{S}_m^l + C_m^l \mathbf{T}_m^l \}. \quad (14) $$

This distinction is natural when studying the generation of a large-scale axisymmetric field. Another choice of shells based on dynamo families will also be used in this paper (see Appendix A.7 and the end of Section 2.2.3). Note that the defined shells are orthogonal, i.e., any scalar product of strictly different shells is zero. In order to simplify the notations, the shell $L$ may represent either axisymmetric or non-axisymmetric shells. Exponents $^0$ and $^*$ denote axisymmetric and non-axisymmetric components in the remainder of this paper.

The shells of magnetic energy in spectral space are then defined by

$$ E_L^{\text{mag}} = \frac{1}{2} \int_S \mathbf{B}_L \cdot \mathbf{B}_L \, d\Omega. \quad (15) $$

This spectrum may be defined at any radial location of the spherical domain (e.g., in a stellar interior). Our spectral analysis is intended to characterize the horizontal scales of velocity and magnetic fields, and can be easily applied at any and as many depths as necessary, depending on the radial locations one wants to focus on. In addition, the different terms of the evolution equation of the magnetic energy (see the next section) explicitly depend upon the radial gradients of the different quantities. The horizontal couplings produced by the vertical interactions are thus taken into account by our description. Finally, it is worth noting that the spectra may be different at various depths in a convective dynamo model. In this work, we will only study pure convection zone dynamos; consequently, we will focus on the spectral interaction in the middle of the convection zone.

### 2.2.3. Spectral Magnetic Energy Equation

In order to obtain the spectral magnetic energy evolution equation, we multiply Equation (5) by $\mathbf{B}_L$ and integrate it over the spherical surface so that

$$ \partial_t E_L^{\text{mag}} = \mathcal{D}_L + \sum_{L_1, L_2} \{ \mathcal{P}_L(L_1, L_2) + \mathcal{F}_L(L_1, L_2) \}, \quad (16) $$

where $\mathcal{D}$ regroups the diffusion terms, $\mathcal{P}$ represents the volumetric production of magnetic energy, and $\mathcal{F}$ is the divergence of the flux of magnetic energy through a spherical surface. Note that the production $\mathcal{P}$ term has to be understood as a general production term that can either be positive (real production) or negative (destruction). The sum over $L_1, L_2$ involves the triangular selection rule $\{|l_1 - l_2| \le l \le l_1 + l_2, m_1 + m_2 = m\}$ that comes naturally from the spherical harmonics coupling (see
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic coupling between spherical harmonics of **U** and **B** for the simple test case. The **B** modes are represented by blue squares, and the **U** modes by red triangles. The black arrows represent the coupling between the modes, the green circles on them represent the resulting modes obtained from the coupling via the triangulation rule of the vectorial product **U** × **B**. The three highlighted green circles in the center correspond to the modes calculated in Appendix B.</figcaption></figure>

(A color version of this figure is available in the online journal.)

Appendix A.2). The expressions of the three contributions $\mathcal{D}$, $\mathcal{P}$, and $\mathcal{F}$ are given by

$$ \mathcal{D}_L(r) = \int_S \{\eta \mathbf{B}_L \cdot \nabla^2 \mathbf{B}_L + \partial_r \eta \mathbf{e}_r \cdot (\mathbf{B}_L \times \mathbf{\nabla} \times \mathbf{B}_L)\} \, d\Omega, \quad (17) $$

$$ \mathcal{P}_L(r, L_1, L_2) = \int_S (\mathbf{U}_{L_1} \times \mathbf{B}_{L_2}) \cdot \mathbf{\nabla} \times \mathbf{B}_L \, d\Omega, \quad (18) $$

$$ \mathcal{F}_L(r, L_1, L_2) = \int_S \mathbf{\nabla} \cdot [(\mathbf{U}_{L_1} \times \mathbf{B}_{L_2}) \times \mathbf{B}_L] \, d\Omega. \quad (19) $$

We have split the diffusive term by considering a diffusivity $\eta$ that only depends on $r$. The production and flux of magnetic energy are discretized over scales so that we compute which scale of the velocity field ($L_1$) is interacting with which scale of the magnetic field ($L_2$) toward a studied scale $L$.

Although expressions (17)–(19) are formally written, they include the evaluation of vectorial products decomposed on the vectorial spherical harmonics basis. This operation is not easily calculated on the (R, S, T) basis, and thus we use an alternative basis (Varshalovich et al. 1975) to compute it. For the sake of simplicity, these details are given in Appendix A.5.

As an alternative to considering axisymmetric and non-axisymmetric spectra, it is instructive to decompose the flow and field into the so-called *primary* (dipolar, antisymmetric) and *secondary* (quadrupolar, symmetric) families (McFadden et al. 1991; Roberts & Stix 1972). These families were proven to be very insightful when characterizing geophysical and astrophysical dynamos (Gubbins & Zhang 1993; DeRosa et al. 2011, 2012). Furthermore, the primary/secondary distinction greatly simplifies the transfer maps of $\mathcal{P}$ and $\mathcal{F}$. Indeed, coupling between fields of the same family always gives *secondary* fields, while coupling between fields of different families always gives a *primary* field (see Appendix A.7). Both of these approaches (i.e., axisymmetric/non-axisymmetric and primary/secondary distinctions) will be used hereafter.

## 2.3. Validation and Illustration of the Method

In this section, we illustrate the coupling calculations for two simple fields. The reader only interested in physical discussions may skip this part and go directly to Section 3.

Since the ASH code is a spectral code, it solves the MHD equations for the spherical harmonics coefficients of the fields. Although it does not explicitly compute the decomposition on the vectorial spherical harmonics basis (Equation (10)), it is straightforward to use this basis in the code by using the transformation relations (Equations (11) and (12)). We have added into the code the ability to compute the different terms (Equations (17)–(19)) of the spectral magnetic energy equation (16).

In order to illustrate and validate both the coupling coefficients (Equation (A17)) for the vectorial product and the general method, we numerically computed a simple analytical test case. We initialize the magnetic and velocity fields in the following way:

$$ \mathbf{B} = a \left( \mathbf{R}_1^0 + \frac{1}{2} \mathbf{R}_1^1 - \frac{1}{2} \mathbf{R}_1^{-1} \right) + b \left( \mathbf{S}_1^0 + \frac{1}{2} \mathbf{S}_1^1 - \frac{1}{2} \mathbf{S}_1^{-1} \right) \\ + c \left( \mathbf{T}_1^0 + \frac{1}{2} \mathbf{T}_1^1 - \frac{1}{2} \mathbf{T}_1^{-1} \right), \\ \bar{\rho} \mathbf{U} = d (\mathbf{R}_2^1 - \mathbf{R}_2^{-1}) + e (\mathbf{S}_2^1 - \mathbf{S}_2^{-1}) + f (\mathbf{T}_2^1 - \mathbf{T}_2^{-1}) $$

where $a, b, c, d, e$ and $f$ are functions of $r$ only. This initialization allows us to test at the same time the axisymmetric/non-axisymmetric and non-axisymmetric/non-axisymmetric coupling schemes between the velocity and the magnetic fields. The low-order harmonics ($l \in \{0, 1, 2\}$) that are involved make the analytic calculation easy. We display in Figure 1 the possible couplings (via vectorial product) between the **U** and **B** fields that we initialized. This is in fact a schematic representation of the triangulation rule that appears in the summation of Equation (A17). The analytical calculation of the values of the three large green circles is given in Appendix B. The resulting vectorial products calculated by the code using the Wigner coefficients show very good agreement with the coefficients calculated analytically (Table 1 in Appendix B).

We stress here that this test has been done for low $m$ and $l$ values. The numerical accuracy of the algorithms used to calculate the Clebsch-Gordan coefficients (and thus $3j$, $6j$, and $9j$ Wigner coefficients) is known to decrease with increasing $l$ and $m$. The calculation routines we use are accurate up to values of $l$ of the order of 500. To do so, we used a multiple-precision package<fnref n="3" /> to simulate large-precision numbers that are needed to compute the ratios of factorials and binomial coefficients that are involved in the Wigner coefficient calculations. However, the calculation time of the transfer functions $\mathcal{P}$ and $\mathcal{F}$ increases dramatically with $l$ and $m$. For practical reasons, when computing fully non-linear dynamos (see Section 4 below), we have chosen to limit the computation of the coupling coefficients to $l_{max} = 70$, even if the effective resolution of such simulations reaches $l_{max} = 340$. From time to time, we do calculate the transfer terms for high $l$s so that we have an indication of how energy is transferred at the smallest scales (see Section 4). Nevertheless, the magnetic-energy-carrying scales in the spectrum are dominated by $l \le 70$ in this case. We thus capture the essential part of the dynamics.

<footnote n="3">³ http://crd-legacy.lbl.gov/~dhbailey/mpdist/</footnote>
<table><caption>**Table 1**

Analytical and Numerical Values of the Vectorial Product **U** × **B** (Validation Case, See Section 2.3)</caption><thead><tr><th>SH Mode</th><th>Analytical Expression</th><th>Analytical Value</th><th>Code Output</th></tr></thead><tbody><tr><td>(1, 0)</td><td><table><tr><td>-<sup>3γ<sub>l</sub></sup>&frasl;<sub>5</sub> &radic;<span style="text-decoration:overline;">4&pi;/3</span></td></tr><tr><td>-<sup>(γ<sub>θ</sub><sup>2</sup> + γ<sub>φ</sub><sup>1</sup>)</sup>&frasl;<sub>5</sub> &radic;<span style="text-decoration:overline;">4&pi;/3</span></td></tr><tr><td><sup>(γ<sub>φ</sub><sup>2</sup> + γ<sub>φ</sub><sup>1</sup>)</sup>&frasl;<sub>5</sub> &radic;<span style="text-decoration:overline;">4&pi;/3</span></td></tr></table></td><td><table><tr><td>-12949906405.7406</td></tr><tr><td>-12949960422.6225</td></tr><tr><td>-0.00337401832383</td></tr></table></td><td><table><tr><td>-12949906405.7402</td></tr><tr><td>-12949960422.6225</td></tr><tr><td>-0.00337401832383</td></tr></table></td></tr><tr><td>(2, 0)</td><td><table><tr><td>0</td></tr><tr><td>0</td></tr><tr><td>0</td></tr></table></td><td><table><tr><td>0.0</td></tr><tr><td>0.0</td></tr><tr><td>0.0</td></tr></table></td><td><table><tr><td>0.0</td></tr><tr><td>0.0</td></tr><tr><td>0.0</td></tr></table></td></tr><tr><td>(3, 0)</td><td><table><tr><td>-<sup>2γ<sub>l</sub></sup>&frasl;<sub>5</sub> &radic;<span style="text-decoration:overline;">4&pi;/7</span></td></tr><tr><td>-<sup>2γ<sub>φ</sub><sup>2</sup></sup>&frasl;<sub>15</sub> &radic;<span style="text-decoration:overline;">4&pi;/7</span></td></tr><tr><td><sup>2γ<sub>φ</sub><sup>2</sup></sup>&frasl;<sub>15</sub> &radic;<span style="text-decoration:overline;">4&pi;/7</span></td></tr></table></td><td><table><tr><td>-5651802509.22854</td></tr><tr><td>3767844764.58569</td></tr><tr><td>-0.00147254232049</td></tr></table></td><td><table><tr><td>-5651802509.22835</td></tr><tr><td>3767844764.58568</td></tr><tr><td>-0.00147254232049</td></tr></table></td></tr></tbody></table>

**Notes.** The values are evaluated at $r = 0.84 R_{\odot}$. The expressions for the $\gamma$ coefficients are given in Appendix B. Numerical results are given with 15 significant digits, i.e., up to the numerical accuracy.

### 3. AXISYMMETRIC $\alpha\Omega$ DYNAMO

In this section, we use the spectral method we developed in Section 2 on two academic cases. First, we explain how the classical $\Omega$ effect (Moffatt 1978) is represented by our formalism (Section 3.1). Then, we calculate the spectral transfers for a mean field $\alpha\Omega$ model (Section 3.2).

#### 3.1. Omega Effect

The complexity of the two spherical harmonics bases may be confusing when it comes to interpreting simple and classical dynamo processes. Thus, hereafter we give a step-by-step explanation of the $\Omega$-effect in the two vectorial spherical harmonics bases formalism.

We start with a purely dipolar poloidal magnetic field that reads (using Equation (A9))

$$ \mathbf{B}_p(r, \theta) = b_r(r)\mathbf{R}_1^0 + b_{\theta}(r)\mathbf{S}_1^0. \quad (20) $$

Then, we want to calculate the effect of a differential rotation that reads

$$ \Omega = A + B \cos^2 \theta. \quad (21) $$

Such differential rotation is usually seen as a “$l = 2$” field. However, it projects on an $l = 3$ component when considering the azimuthal component of the velocity $\mathbf{U}_{\varphi} = r \sin\theta \Omega \mathbf{e}_{\varphi}$ (see Roberts & Stix 1972), which reads

$$ \mathbf{U}_{\varphi}(r, \theta) = U_{\varphi}(r)(A \sin\theta + B \sin\theta \cos^2\theta)\mathbf{e}_{\varphi} \quad (22) $$

$$ \sim U_{\varphi}(r) \partial_{\theta} Y_{3}^{0} \mathbf{e}_{\varphi} = U_{\varphi}(r) \mathbf{T}_{3}^{0}. \quad (23) $$

In general, Equation (22) should project on both $\mathbf{T}_1^0$ and $\mathbf{T}_3^0$. For the sake of simplicity, here we select a profile of differential rotation that is purely described by a ($l = 3, m = 0$) harmonic, which corresponds to $B = -5A$ (Equation (21)). We simply apply the curl operator (Equation (A9)) and make use of the coupling relations (Equation (A17)) to obtain the production of **B** in the induction equation,

$$ \nabla \times (\mathbf{U}_{\varphi} \times \mathbf{B}_p) = f(b_r, b_{\theta}, U_{\varphi})\mathbf{T}_2^0. \quad (24) $$

We determined that the action of differential rotation on a purely axisymmetric poloidal field creates a toroidal field $B_{\varphi} \propto \sin\theta \cos\theta$. With our notations, this type of field will be labeled as an “$l = 2$” field.

An additional feature of the differential rotation can also be observed from this little analysis. We immediately note that for axisymmetric fields, the first Wigner coefficient involved in the coupling between two shells $L_1$ and $L_2$ is zero if $3+l_1+l_2$ is odd (Equation (A18)), i.e., if $l_1$ and $l_2$ are of the opposite parity. The shearing effect of differential rotation will then always couple axisymmetric scales of the magnetic field that are of opposite parity, which will be observed in the transfer maps in more complex cases (e.g., Figures 3(b) and 13).

This simple example strikingly highlights how the vectorial product formula (Equation (A17)) couples together two simple fields. This description of the $\Omega$ effect will guide our analysis in Sections 3.2 and 4.

#### 3.2. Case of a Cyclic Mean Field Dynamo

We use the ASH code (Clune et al. 1999; Brun et al. 2004) to simulate an axisymmetric mean field dynamo (Charbonneau 2010; Jouve et al. 2008). To do so, we solve only the induction equation, uniquely considering a prescribed differential rotation profile (see also Section 3.8 of Jouve et al. 2008, for a similar use of a 3D spherical code to model $\alpha - \Omega$ dynamos).

Our radial domain is defined between $r_b = 0.6 R_{\odot}$ and $r_t = 0.966 R_{\odot}$. We use a resolution of $N_r \times N_{\theta} \times N_{\varphi} = 64 \times 128 \times 256$. We choose a solar differential rotation profile $\Omega_{\text{DR}}(r, \theta) = -\partial_{\theta} Z_{\text{DR}}/(r^2 \sin\theta)$ through the toroidal component $Z_{\text{DR}}$ of the momentum, which is, in the frame rotating at $\Omega_0 = 2.6 \times 10^{-6} \text{s}^{-1}$:

$$ Z_{\text{DR}}(r, \theta) = Z_t(r) \left( A \cos\theta + \frac{B}{3} \cos^3\theta + \frac{C}{5} \cos^5\theta \right). \quad (25) $$

From Schou et al. (1998), we take $A = 257$ nHz, $B = 321$ nHz, and $C = 529$ nHz. The differential rotation then naturally projects on $\mathbf{U}_1^0$, $\mathbf{U}_3^0$, and $\mathbf{U}_5^0$. The radial profile $Z_t(r)$ is chosen such that it simulates a stable region at the base of the domain and is defined by

$$ Z_t(r) = \bar{\rho} \frac{r^2}{2} \left[ 1 + \tanh \left( \frac{r - 4.87 \times 10^{10}}{2 \times 10^9} \right) \right]. \quad (26) $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Solar-like differential rotation profile computed from Schou et al. (1998). The tachocline is located near $r = 0.7 R_\odot$, and the base of the tachocline rotates at the solar rotation rate $\Omega_0 = 2.6 \times 10^{-6} \text{ s}^{-1}$. (b) Butterfly diagram in the axisymmetric $\alpha - \Omega$ dynamo (only $B_r$ as a function of time and latitude is shown in the upper convection zone). The vertical black dotted line represents the time at which we display the spectral interactions in Figure 3.</figcaption></figure>

(A color version of this figure is available in the online journal.)

We initialize our magnetic field with a seed $l = 3$ poloidal (antisymmetric and axisymmetric) field. Finally, we add an $\alpha$ effect to the induction equation such that

$$ \partial_t \mathbf{B} = \nabla \times (\mathbf{U} \times \mathbf{B} + \alpha B_\varphi \mathbf{e}_\varphi) - \nabla \times (\eta \nabla \times \mathbf{B}). \quad (27) $$

Since we do not take into account in this simple case the feedback of the Lorentz force on the flow via the Navier-Stokes equations, since we only solve the induction equation, we need to quench the $\alpha$ effect. Hence, $\alpha$ is defined by

$$ \alpha(r, \theta) = \alpha_0 e^{\left(-\frac{r-0.75 R_\odot}{0.05 R_\odot}\right)^2} \frac{\cos \theta}{1 + \left(\frac{|\mathbf{B}|}{B_q}\right)^2}. \quad (28) $$

This is the simplest $\alpha$ that is needed to trigger an oscillating solar-like dynamo (Charbonneau 2010); it is anti-symmetric with respect to the equator. The radial profile of $\alpha$ is localized near the base of the convection zone and the quenching value is given by $B_q = 10^3$ G. We have deliberately chosen an $\alpha$-effect that operates only on the poloidal component of the induction equation, therefore computing an $\alpha\Omega$ mean field dynamo (Moffatt 1978). This $\alpha\Omega$ dynamo exhibits the characteristic butterfly diagram showed in Figure 2(b) (at $r = 0.92 R_\odot$). Although this $\alpha$ profile is ad hoc and only one among the many profiles that were tested in the literature (e.g., Roberts & Stix 1972; Charbonneau & MacGregor 1997; Bonanno et al. 2002; Zhang et al. 2003; Jouve et al. 2008), we chose this form because it easily triggers an oscillatory dynamo and its effect in spectral space can be easily calculated. It is consequently a good choice to illustrate our new spectral method. With the parameters we chose, the cycle period is of the order of 400 days (see Figure 2(b)).

The extra $\alpha$ effect adds a new term in the spectral energy equation (16) that can lead to complex formula in spectral space. We rewrite the energy equation

$$ \partial_t E_L^{\text{mag}} = \mathcal{D}_L + \int_S \nabla \times (\alpha B_\varphi \mathbf{e}_\varphi)_L \cdot \mathbf{B}_L \, d\Omega + \sum_{L_1, L_2} \{\mathcal{P}_L(L_1, L_2) + \mathcal{F}_L(L_1, L_2)\}. \quad (29) $$

The interested reader may read Appendix A.4 for a complete spectral description of this $\alpha$ effect.

Wherever $|\mathbf{B}|$ is not too large, the quenching part of the $\alpha$ effect is negligible. In that case, the $\alpha$ effect (which restores poloidal field from toroidal field) simply couples an $L$ shell of toroidal field to its neighboring shells of the poloidal field, namely $L-1$ and $L+1$. When $|\mathbf{B}|/B_q$ becomes large, the $\alpha$ effect is quenched and the poloidal magnetic field stops being restored. When it is sufficiently low, $\alpha$ stops being quenched and the poloidal field grows again. This sets up a simple feedback mechanism and a cycle is established.

The magnetic energy spectrum is dominated by an $L=2$ component that sets the phase of the total cycle. Various shells' energy oscillate with roughly the same period, but are generally out of phase. This phase shift is a natural ingredient that allows the reversal of the overall field polarity (Knobloch et al. 1998; Tobias 2002).

We stress here that the magnetic field created in this experiment is of the primary family $((\mathbf{R}_{2l+1}^0, \mathbf{S}_{2l+1}^0) \text{ and } \mathbf{T}_{2l}^0$, see Appendix A.7). Our initial magnetic field is a poloidal primary field $(\mathbf{R}_3^0, \mathbf{S}_3^0)$. As a result, the toroidal field created through the $\Omega$-effect is also a primary field $(\mathbf{T}_2^0$, see Section 3.1). Then, our $\alpha$-effect, which creates the poloidal field from the toroidal field, transforms the primary toroidal field into a primary poloidal field $(\mathbf{R}_{1,3}^0, \mathbf{S}_{1,3}^0)$, which is of the same type as our initial magnetic field. Hence, no secondary field can be created in the simulation (which is confirmed by our results), and the $\alpha$-effect can only act on the primary toroidal field to create a primary poloidal field. This is a direct consequence of the well-known separability property of the induction equation between the dipolar and quadrupolar families when symmetric flows and antisymmetric $\alpha$ effect are chosen (Gubbins & Zhang 1993).

We display in Figure 3(a) the evolution of the different terms of the magnetic energy equation (29) for $L=2$ at $r = 0.92 R_\odot$ during the same time period as the butterfly diagram in Figure 2(b). The primary toroidal field energy clearly evolves due to the production $\mathcal{P}$ (dashed blue line) and flux $\mathcal{F}$ (dashed red line) terms that account for the effect of differential rotation on the magnetic field (as expected, since there are no other production nor advection terms). The two terms cancel each other out with a small time lag, and their sum combines with the ohmic diffusion (dotted line) to produce oscillations (solid line) of the total $L=2$ energy. Note that the $\alpha$ effect plays no role at $r = 0.92 R_\odot$ since it is concentrated at the base of the “convection zone” (Equation (28)).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Evolution of the different terms of Equation (29) for the shell $L = 2$. Production $P$ is the dashed blue line, flux $F$ is the dashed red line, diffusion $D_1 + D_2$ is the dotted black line, and the $\alpha$ effect is the dash-dotted green line. The black plain line is the total of all of the contributions. (b) Flux $F$ contribution to the $L = 2$ shell (see Equation (29)) at $t = 1500$ days (vertical dotted line in panel (a)). The 2D color maps are the $\mathbf{B}^0 - \mathbf{U}^0$ transfer functions (dark red is the maximum value, black the minimum). The oblique black lines represent the boundaries of the triangular selection rule.</figcaption></figure>

(A color version of this figure is available in the online journal.)

Our new method allows us to characterize how scales interact to produce this behavior. We display in Figure 3(b) the transfer map for the flux $F$ term of Equation (29) for the $L = 2$ shell at its maximum. The differential rotation is composed of the $\mathbf{U}_1^0$, $\mathbf{U}_3^0$, and $\mathbf{U}_5^0$ shells (Equation (25)). The transfer maps during minima (not shown here) are qualitatively opposite, which means that all of the couplings between the shells reverse sign during the cycle. This reversal of all shells is a simple, direct consequence of the reversal of the whole magnetic field. At this point, the poloidal magnetic energy (not shown here) evolves because of the ohmic diffusion of the $\alpha$-driven poloidal field at deeper radii. Here, the poloidal magnetic field couples with the differential rotation to transfer energy to the toroidal $L = 2$ magnetic shell. The $\mathbf{U}_3^0 - \mathbf{B}_1^0$ appears to be the dominant interaction that sets the $L = 2$ cycle. Interestingly, we will recover this feature in the turbulent (convective) dynamo described in Section 4 (see Figure 13(b)).

This $\alpha\Omega$ dynamo provides a simple example of how our diagnostic may be interpreted in the context of a stellar dynamo. Based on how our diagnostic highlights the saturating properties of the solar differential rotation in an $\alpha - \Omega$ case, we now apply it to a turbulent dynamo triggered in a stellar convection zone that also exhibits a solar-like differential rotation profile.

## 4. NONLINEAR CONVECTIVE DYNAMO

We use the general method described in Section 2 and validated in Section 3 to study dynamo action in a global (spherical) nonlinear convection zone. Contrary to Section 3, we now solve the full set of MHD equations and do not introduce any $\alpha$ effect. We model a turbulent solar convection zone (Brun et al. 2004; Jouve & Brun 2009; Pinto & Brun 2012) that develops a solar-like differential rotation profile (Figure 4(a)), with a fast equator and slow poles. We display the convective patterns we obtain in Figure 4(b). We recover the well-known “banana”-shaped cells at the equator and more patchy patterns at higher latitudes. Our choice of parameters yields a mildly turbulent state (based on the maximum amplitude of the velocity, the Reynolds number in the middle of the convection zone is of the order of 800).

We display in Figure 4(c) the kinetic energy spectra in the rotating frame at the center of the convection zone as a function of the shell $L$. We separate the axisymmetric component (the plain blue line) from the non-axisymmetric component (the dashed red line), and the dotted black line is the total spectrum. Note that two peaks at $L = 3, 5$ dominate the kinetic energy spectrum. They represent the differential rotation of the azimuthal component of the toroidal velocity (see Section 3.1).

We initialize a peaked $(l, m) = (9, 5)$ non-axisymmetric magnetic field (Figure 5(a)) throughout the convection zone by setting

$$ \mathbf{B} = \frac{10 B_0 R_{\odot}^2}{r^2} \left(\frac{R_b}{r}\right)^9 \mathbf{R}_9^5 - \frac{B_0 R_{\odot}^2}{r^2} \left(\frac{R_b}{r}\right)^9 \mathbf{S}_9^5. \quad (30) $$

We set $B_0 = 100$ G so that the initial magnetic energy contained in the $L = 9$ shell is comparable to the kinetic energy at that scale (see Figure 5(b)).

The magnetic Prandtl number throughout the convection zone is set to $P_m = (v/\eta) = 4$, leading to a magnetic Reynolds number of the order of 3200 at the mid-convection zone (based on the maximum amplitude of the velocity). Such a set of parameters triggers a dynamo instability and the growth of magnetic energy (see Figure 6(c) in the following).

The initialization we chose allows us to directly observe how a significant amount of energy can be transferred to large scales. We also did the same numerical experiment varying the initial conditions. By initializing roughly the same amount of energy distributed over the whole scales, we obtained the same statistical saturated state. Hence, this proves that in this case, the initial scale is forgotten when the dynamo saturates.

The complex interactions between the convective motions and the initially peaked magnetic field lead to the construction of the magnetic energy spectrum. The saturated magnetic energy after 600 days of evolution is displayed in Figure 5(c) in physical space, in the middle of the convection zone. In the remainder of this section, we characterize how such a state is obtained and maintained. We distinguish two regimes: the development of the spectrum shape (the kinematic regime; Section 4.1), and its saturation and sustainment (the nonlinear regime; Sections 4.2 and 4.3). The results of this section will be summarized in
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Progenitor hydrodynamical state. (a) Mean differential rotation profile averaged over 3 months. (b) Time-dependent convective patterns at the top of the convection zone with dark tones representing downflows. (c) Kinetic energy spectra in the middle of the convection zone. The axisymmetric spectra are in plain blue, the non-axisymmetric spectra in dashed red, and the total spectra in dotted black.</figcaption></figure>

(A color version of this figure is available in the online journal.)

<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Initial magnetic energy (the seed poloidal streamfunction is a $l = 9, m = 5$ spherical harmonic) in cgs, in the middle of the convection zone. (b) Magnetic energy spectra one timestep after the introduction of the peaked magnetic field in our turbulent convection zone. The axisymmetric spectra are in plain blue, the non-axisymmetric spectra in dashed red, and the total spectra in dotted black. (c) Saturated magnetic energy in cgs, 600 days after the introduction of the seed magnetic field. The color scale is logarithmic.</figcaption></figure>

(A color version of this figure is available in the online journal.)

Figure 14. We recall here that no $\alpha$ effect has been added to the induction equation (5); dynamo action is naturally achieved since convection is 3D and $R_m > R_m^{crit}$ (Brun et al. 2004).

## 4.1. Creation of Magnetic Energy Spectrum: Kinematic Phase

We plot in Figure 6(a) the evolution of the non-axisymmetric magnetic energy spectrum. The initial spectrum is plotted in blue, and the saturated spectrum in red. In addition, we display in Figure 6(c) the evolution of the magnetic energy for six different $L$ shells. The total energy evolution is also shown (plain thick line). The initial evolution ($t < 100$ days) is shown in logarithmic scale. The saturation of magnetic energy is reached at $t \sim 300$ days.

We also ran another numerical experiment where we artificially suppressed the Lorentz force and the ohmic heating in the momentum and energy equations (i.e., effectively running a kinematic dynamo). On average, the relative difference from the fully nonlinear case starts to become significantly different (departure of order one) roughly 10 days after the introduction of the magnetic field (the exact length of the kinematic phase depends on the scale considered). Hereafter, we detail how the non-axisymmetric (Section 4.1.1) and the axisymmetric (Section 4.1.2) spectra are created during these first days, which we will refer to as the kinematic phase.

### 4.1.1. Creation of the Non-axisymmetric Spectrum

At first, we observe that all of the $L$ shells gain energy (Figure 6(a)), except the $L = 9$ shell, which loses energy

because it is redistributed throughout the entire domain by the convective flows (Figure 6(c)). It stops decaying at $t \sim 17$ days. We identify four regions in the non-axisymmetric spectrum that exhibit different behaviors. We define (I) the large-scale zone by $1 \le L \le 4$, (II) the neighborhood zone by $4 \le L \le 13$, and (III) the plateau zone by $13 \le L \le 60$. (IV) The small-scale zone ($L \gtrsim 60$) starts at the highest diffusive scale, which is the highest viscous scale $l \sim 60$ based on the first scale at which the local Reynolds number is lower than 1. It also includes the magnetic dissipative scales ($L \gtrsim 120$). The four zones are separated by the three dotted vertical lines in Figure 6(a). In order to understand how the spectrum is built, we display the contributions from the different terms of Equations (A25)–(A28) in Figure 6(b). Those contributions are taken shortly after the introduction of the magnetic field. They correspond to the spectrum plotted with a dotted line in Figure 6(a). We recall that we fully calculate all of the coupling terms up to $L = 70$.

The energy transfers around the $L = 9$ shell (the neighborhood of zone II) are dominated by $U^0-B^*$ interactions from both the production $\mathcal{P}$ and the flux $\mathcal{F}$ terms. We recall here that both $\mathcal{P}$ and $\mathcal{F}$ represent generic transfer functions, which can either be positive or negative. Dissipation is negligible in zone (II), even for the $L = 9$ shell that initially contains the energy. The $L = 9$ energy then decreases through the interaction of $B_9^*$ and the differential rotation $U_3^0$ that shears the magnetic field (see Section 3.1). The energy is preferentially redistributed to $E_{11}^*$ and $E_7^*$. For those two shells, the production $\mathcal{P}$ and flux $\mathcal{F}$ terms contribute positively to the creation of the spectrum (Figure 6(b)). We display in Figure 7 the detailed contribution
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a) Evolution of the non-axisymmetric magnetic energy spectra. The initial spectra is blue, the saturated spectrum red. Two slopes ($L^{-6}$ and $L^{-3}$) are given as references for the small scales. (b) Contributions to the non-axisymmetric magnetic energy evolution in the initial phase (dotted line in panel (a)). The total is in plain black, *D* in dotted black, *P* in blue, and *F* in red. (c) Evolution of non-axisymmetric energy of shells 1, 3, 7, 9, 11, and 42. Note the fast early evolution of $E_7^*$ and $E_{11}^*$ due to shearing of $B_0^*$ by the differential rotation. The thick plain line is the total non-axisymmetric energy. The abscissa is in log scale for $t < 100$ days and in linear scale for $t > 100$ days. Even though the initial growth is exponential, we chose to represent it in log-log scale to make it appear clearly in the evolution plot.</figcaption></figure>

(A color version of this figure is available in the online journal.)

of *P* and *F* to $E_7^*$. We only display contributions from **B**\* because the axisymmetric magnetic energy is very small initially. The **U**\*-\**B**\* interactions are displayed in panel (a), and the **U**\*\*-\**B**\* interactions in panel (b). We sum over the velocity shells to plot

the production term against **B**\* in panel (c). We note that the summed contribution is dominated by $B_9^* - U_3^0$ interactions, as expected. Also, we note that energy is directly transferred from $E_9$ to $E_7^*$, such that the $L=8$ shell is not involved in the transfer. This is true for all of the shells in zone (II) and implies that the transfer of energy is non-local, even for shells close to the initial energetic shell.

Due to the triangular selection rule, the $U_3^0 - B_9^*$ interaction can only act in zone (II). Indeed, $L$ must be strictly greater than 12 in zones (III–IV) and strictly lower than 6 in zone (I). $U_3^0$ and $B_9^*$ initially dominate, respectively, the kinetic and magnetic energy spectra. Their interaction was consequently dominant in zone (II), and we expect different kinds of spectral transfers in the other zones. This zone exists because of our choice of initial condition. The very early evolution would have been changed if we had chosen a different initial shell. However, as stated before, this initial scale is forgotten when the saturated state is reached (Figure 6(a)).

The dynamics of zones (I), (III), and (IV) are dominated by two initially competing effects: a direct non-local $B_9^* - U^*$ transfer of energy and an effective shearing of neighbor shells by the large-scale differential rotation ($B^* - U_3^0$ interactions). These two effects are demonstrated in Figure 8 for the $L=42$ (zone III).

In the case of the large-scale zone (I) (not shown here), the evolution is dominated by both *P* and *F*. The interactions between $B_9^*$ and $U^*$ alternate signs depending on the $U^*$ shell considered. We also stress that the interactions involving other $B^*$ shells are not negligible. The differential rotation action is completely negligible compared to $U^* - B^*$ interactions in zone (I).

In the case of the plateau zone (III), an almost flat profile in the log-log plot is observed in Figure 6(a) (hence its name). This plateau is characteristic of convective flows that usually exhibit a broad spectrum between the injection and inertial ranges (Figure 4(c)). The evolution of the spectrum is dominated only by the $P_L(L_1, L_2)$ contributions ($F$ is negligible), and in particular by the coupling between (non-axisymmetric) $U^*$ and $B^*$ (Figure 6(b)). Hence, it is a non-local transfer of magnetic energy that creates the spectrum. All of the shells in zone (III) receive energy mainly through this non-local mechanism. As a result, the energy transfer is very sensitive to the kinetic energy contained in the $U_L^*$ shells involved in the coupling. This explains why the magnetic energy spectrum reflects the kinetic energy spectrum in this region.

Although the $U^* - B^*$ interactions dominate (Figure 8), we stress that the $U^0 - B^*$ interactions exhibit a direct cascade pattern. $E_{42}^*$ receives energy from $E_{40}^*$ through the $U_3^0 - B_{40}^*$ interactions, and gives energy to $E_{44}^*$ through the $U_3^0 - B_{44}^*$ interactions (see Figure 8(b)). Even if the triadic interaction involves the large-scale velocity $U_3^0$, we nonetheless refer to this effect as a *cascade*. The velocity field only acts here as a *mediator*, and the scales of the magnetic field involved in the magnetic energy transfer are at the same scale. It is consequently a cascade when considering the scales of the magnetic field.

The energy transfers in zone (IV) (not shown here) are very similar to zone (III). A noticeable difference is that the cascade of energy triggered by the shear of the differential rotation $U_3^0$ is much less efficient since the smallest scales hardly feel the large-scale rotation profile. Finally, ohmic diffusion acts over the entire zone (IV) and tends to dissipate energy. It has a sufficiently lower amplitude than the non-local transfers so that it does not initially dictate the spectrum shape. It will nevertheless contribute to the saturation process (Section 4.2).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Production $\mathcal{P}$ and flux $\mathcal{F}$ contributions to the non-axisymmetric shell $L = 7$. The 2D color maps are the $B^* - U^*$ and $B^* - U^0$ transfer functions and the 1D plot is the sum of the transfer functions over the $U$ shells. The horizontal dotted line labels the $B_7$ shell. The oblique black lines represent the boundaries of the triangular selection rule.</figcaption></figure>

(A color version of this figure is available in the online journal.)

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Production $\mathcal{P}$ contribution to the non-axisymmetric shell $L = 42$.</figcaption></figure>

(A color version of this figure is available in the online journal.)

### 4.1.2. Creation of the Axisymmetric Spectrum

We now characterize the creation of the axisymmetric spectrum. We display in Figure 9(a) the evolution of the axisymmetric component of the magnetic energy. We recall that since we initialized the dynamo with a purely non-axisymmetric field, the initial axisymmetric spectrum is null. After one timestep, the axisymmetric magnetic energy is orders of magnitude lower than the non-axisymmetric spectrum (Figure 5(b)). The global shape of the axisymmetric spectrum is created very rapidly; all

of the shells gain energy at about the same rate until they saturate. The initial exponential growth rate is the same for both the axisymmetric and non-axisymmetric spectrum, approximately $0.6 \text{ days}^{-1}$ (which corresponds to a timescale that is approximately 17 times lower than the convective turnover time). This can also be observed in Figure 9(c), where we plot the evolution of few shells against time. They all gain energy at about the same rate initially and then slowly tend to a saturated state. The axisymmetric shells considered have comparable energy since the spectrum is essentially flat at scales $L \le 30$ (Figure 9(a)),
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Same as for Figure 6, but for the axisymmetric part of the spectrum. In panel (a), the initial spectrum is zero.</figcaption></figure>

(A color version of this figure is available in the online journal.)

which was not the case for the non-axisymmetric spectrum (Figure 6(c)). We observe in Figure 9(b) that the flux term $\mathcal{F}$ plays a major role between $L = 2$ and $L = 13$. This means that the creation of the spectrum is dominated by the radial interactions at those scales. The two flux curves exhibit a *saw-tooth* pattern that is again reminiscent of the differential rotation energy shells (see Sections 3.1 and 4.1.1). At higher $L$, the evolution of the spectrum is the result of a complex interplay between the production and flux terms.

More interestingly, the dipole ($L = 1$) evolution is dominated by the production term through the interaction between the non-

axisymmetric magnetic field and velocity field. We display the detailed transfer maps for this scale in Figure 10. We observe that the large-scale magnetic field is mainly created by the interplay between $\mathbf{B}_7^*$ and $\mathbf{U}_6^{*-8}$. The transfers involving $\mathbf{B}_9^*$ (where the energy is originally mainly contained) act negatively and do not dominate the transfer of magnetic energy. This is consistent with the fact that the whole axisymmetric spectrum shape is rapidly created and only gains energy globally afterward. It does not depend on the scale at which we initially put the non-axisymmetric magnetic energy. Since the energy is not transferred directly from the initial reservoir of energy $E_0^*$, we already see preferred transfers toward the large-scale dipole involving $\mathbf{B}_7^*$, which is one of the highest energy scales of the non-axisymmetric spectrum at this time. This effect shall be confirmed during the saturation phase (Section 4.3). The creation of the axisymmetric magnetic energy spectrum seems to depend essentially on the initial hydrodynamic convective spectrum (as expected in such kinematic phase).

## 4.2. Nonlinear Saturation of the Smallest Scales

Following Section 4.1, we now detail the saturation and sustainment of the magnetic energy spectrum at small scales. By 500 days, the axisymmetric and non-axisymmetric spectra are saturated (Figures 9(c) and 6(c)).

The flux $\mathcal{F}$ contribution is likely to never be null at the largest scales since it represents the flux of magnetic energy through the horizontal surface at the middle of the turbulent convection zone. In order to saturate the magnetic energy (i.e., to get $dE_L^{\text{mag}}/dt = 0$), $\mathcal{D}$ and/or $\mathcal{P}$ have to compensate for $\mathcal{F}$. In the first three zones, diffusion is negligible. Hence, $\mathcal{P}$ naturally tends to cancel out $\mathcal{F}$ in those zones (see Section 3.2 for a simple version of this cancellation effect). The cancellation effect is such that $\mathcal{F}(\mathbf{U}^*, \mathbf{B}^*)$ tends to cancel out $\mathcal{P}(\mathbf{U}^*, \mathbf{B}^*)$. This is also the case for $\mathcal{F}(\mathbf{U}^0, \mathbf{B}^*)$, $\mathcal{F}(\mathbf{U}^*, \mathbf{B}^0)$, and $\mathcal{F}(\mathbf{U}^0, \mathbf{B}^0)$.

In spite of the cancellation of the different contributions, characteristic patterns can still be identified. The more distinctive pattern we identified in Section 4.1 was the direct cascade of magnetic energy in zone (III). It turns out that we still observe it and that it slightly dominates the transfer terms during the saturation phase. We display in Figure 11 the production contribution to the non-axisymmetric magnetic energy evolution 600 days after the magnetic field introduction. We recover the direct cascade of energy in the production contribution, which was already present in Figure 8. This direct cascade of energy is associated with an inverse cascade of energy carried by the flux contribution, which opposes the production term during the saturation phase. Both cascades are of the same order of magnitude and tend to cancel each other out. They are associated with the axisymmetric component $\mathbf{U}_3^0$ (the differential rotation) and the non-axisymmetric components of $\mathbf{B}$. The contributions of the non-axisymmetric components of $\mathbf{U}$ involve more shells, but their net effect is a bit lower than the shear from differential rotation (panel (a) in Figure 11). On this panel, no particular global pattern can be identified.

The transfers of magnetic energy appear to be very interesting in zone (IV), where diffusion acts significantly. In order to saturate, $\mathcal{P}$ and $\mathcal{F}$ have to combine to cancel $\mathcal{D}$. For the non-axisymmetric spectrum, it is the production term that dominates over the flux term to compensate for diffusion. In addition, the production term in zone (IV) exhibits a very particular generalized cascade shape. This cascade could not be identified during the early evolution because it was dominated by the non-local transfer from $\mathbf{B}_0^*$. We display in Figure 12(a) the
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Production *P* for the axisymmetric shell *L* = 1 during the initial state. Both **U**<sup>0</sup>–**B**<sup>0</sup> and **U**\*–**B**\* couplings are shown. (A color version of this figure is available in the online journal.)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11.** Production *P* contribution to the non-axisymmetric shell *L* = 42 during the saturation phase. Interactions between **B**<sup>0</sup> and **U**\* are negligible. (A color version of this figure is available in the online journal.)</figcaption></figure>

**U**\*-**B**\* production map toward *E*<sub>152</sub><sup>*</sup>. The other interactions are negligible. We note that the map is dominated by a positive contribution (red) under the horizontal dashed line (*L* = 152) and by a negative contribution (blue) above. This is confirmed by the plot in panel (b) where the transfers have been summed over the velocity shells. This cascade is of a different kind than that observed in zone (III) (Figure 11). Here, no clear velocity shell dominates the transfer map (Figure 12(a)). It is a *generalized* *cascade* that results from the coupling between many magnetic shells (around *L* = 152) and all of the largest velocity scales.

Hence, the velocity scales involved in the cascade are not local compared to the magnetic field scale considered.

Trying to simplify the complex two-dimensional (2D) transfer maps, one may isolate the main contributing couplings to the different evolution terms. By constantly doing this for the non-axisymmetric spectrum at small scales, we find that the percentage of couplings that account for 90% of the contributing terms typically varies from nearly 1% to 70% of the calculated couples. As a result, we demonstrate here that the complex dynamo process occurring in a 3D turbulent convection zone
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Production *P* contribution to the non-axisymmetric shell *L* = 152 during the saturation phase. The red line is the smoothed contribution that exhibits the characteristic S-shape of the generalized cascade.</figcaption></figure>

(A color version of this figure is available in the online journal.)

involves many modes that interact through non-trivial triadic interactions. Then, the dynamics of the smallest scales can hardly be reduced to the evolution of a small set of modes.

Finally, the analysis of the axisymmetric αΩ dynamo in Section 3.2 sheds light on the importance of the families of symmetry (with respect to the equator) of the fields. The instantaneous convective motions do not exhibit any particular symmetry at any scale and the kinetic energy spectrum is a mixture of both primary and secondary velocities. The differential rotation is the only velocity feature that has a clear symmetry (secondary family, see Section 3.1) and that has a large influence on the magnetic energy spectrum. It is involved in the magnetic energy cascade in zone (III) and shears both primary and secondary magnetic fields to cascade primary and secondary magnetic energy. Thus, it does not select a particular symmetry. Indeed, the ratio of primary (antisymmetric) to secondary (symmetric) magnetic energy varies with time for all of the shells and does not settle even during the saturation phase. The presence of complex flows, often breaking the equatorial symmetry, yields a strong coupling of both dynamo families (as in the Sun; see DeRosa et al. 2012), contrary to the simpler mean field dynamo models (see Section 3).

### 4.3. Sustainment of the Mean Large-scale Magnetic Field

Given their key role in setting the overall magnetic polarity in the Sun (DeRosa et al. 2012), we now detail the main contributions to the saturation and sustainment of the large-scale axisymmetric dipole ($l = 1, m = 0$) and quadrupole ($l = 2, m = 0$) fields.

At the late phase of the simulation, the large-scale axisymmetric spectrum is fully saturated (Figure 9(c)). The saturation is obtained thanks to the compensation of the production and flux terms, which is similar to the saturation of the mid-scales (see previous section). The large-scale dipole ($l = 1, m = 0$) saturation process differs significantly from its creation. We display in Figure 13 the production maps for the axisymmetric dipole and quadrupole averaged over 150 days during the saturated state. The transfers maps of $\mathcal{F}$ (not shown here) are the exact opposite of maps (a) and (c) for $\mathcal{P}$. We see that both the axisymmetric and non-axisymmetric fields significantly contribute to the saturation and sustainment of the large-scale dipole. In particular, two main contributors emerge. First (panel (a)), the coupling of the differential rotation $\mathbf{U}_3^0$ with the large-scale $\mathbf{B}_0^0$ field dominates the axisymmetric contributions. This effect is more likely

to represent the shearing of the large-scale poloidal multipole $\mathbf{B}_p^0$ by the large-scale toroidal differential rotation.

Second, the non-axisymmetric contributions (panel (c)) are at least equally important for the saturation of the dipole. In particular, the interaction $\mathbf{U}_{23}^* - \mathbf{B}_{23}^*$ dominates the non-axisymmetric contributions. Thus, it is a *non-local* interaction that saturates the large-scale magnetic dipole. Furthermore, $\mathbf{B}_{23}^*$ is one of the most energetic shells of the magnetic energy spectrum (Figure 6(a)). This directly indicates the importance of the mid-scale part of both the kinetic and magnetic energy spectra for the saturation level of the large-scale magnetic dipole.

We can note here that the major contributions of $\mathcal{P}$ for the saturation of the dipole are all positive. They are balanced by negative contributions from $\mathcal{F}$. Consequently, if the differential rotation was more efficient, or if the $\mathbf{U}_{23}^* - \mathbf{B}_{23}^*$ interaction possessed more energy, then the saturation level of the large-scale dipole would be much higher.

Since our magnetic Prandtl number is 4, the peaks of the kinetic and magnetic energy spectra are likely to be shifted. At saturation, the couplings are nonetheless dominated by the peak of the magnetic energy spectrum that occurs at a smaller scale than the peak of the kinetic energy spectrum. Changing the magnetic Prandtl number will cause the separation of the peaks to change. If the peaks separate more, our results suggest that the saturating interaction will involve smaller scale velocity and magnetic fields. The velocity field involved is likely to be less energetic, which could trigger a smaller saturating interaction, and in turn a lower energy state for the large-scale dipole. If the peaks are closer (or eventually switch), then the picture becomes more complicated and we cannot predict if the saturating interaction will remain fixed by the peak of the magnetic energy spectrum. The exploration of this parameter space is left for future work.

The large-scale quadrupole also saturates thanks to both the axisymmetric and non-axisymmetric fields (panels (b) and (d)). The axisymmetric contributions (panel (b)) are very similar to the dipole case and are dominated by the differential rotation. The differential rotation shears both $\mathbf{B}_1^0$ and $\mathbf{B}_5^0$, which is opposed by the flux term to saturate the quadrupole. Again, this effect accounts for the saturation of the toroidal quadrupolar field. Hence, the saturated level of the poloidal dipolar field (panel (c)) plays a major role in the saturation of the toroidal quadrupolar field.

The poloidal quadrupolar field is then saturated through non-axisymmetric interactions (panel (d)). The contributions are
<figure><img src="image_18.png" /><figcaption>**Figure 13.** Saturation of the large-scale axisymmetric dipole (panels (a) and (c), production $\mathcal{P}$ for $E_1^0$) and quadrupole (panels (b) and (d), production $\mathcal{P}$ for $E_2^0$). The interaction maps are time averaged over a period of ~150 days. Both $U^0-B^0$ and $U^*-B^*$ couplings are shown.
(A color version of this figure is available in the online journal.)</figcaption></figure>

again very non-local, though in this case no particular scale dominates the saturation process. Hence, we may expect that the saturation process of the axisymmetric quadrupole will have a very different dependency on the magnetic Prandtl number than the axisymmetric dipole.

## 5. CONCLUSIONS AND PERSPECTIVES

In this paper, we developed and validated a new spectral analysis method suited for spherical objects. Using a two vectorial spherical harmonics basis, we were able to calculate the transfer functions of magnetic energy in spectral space. We can calculate the coupling coefficients up to $l \sim 500$. For the first time in such studies, the complete 2D transfer maps have been calculated to characterize the full triadic interactions.

After a quick numerical validation, we first applied our method to a simplified $\alpha\Omega$ dynamo case. Such axisymmetric models are very well known to trigger cyclic dynamos (Charbonneau 2010) with our choice of a symmetric (with respect to equator) velocity field and an antisymmetric $\alpha$ effect. The clear separation between the dipolar and quadrupolar families was illustrated thanks to our new diagnostic. The production (i.e., on a spherical surface) and flux (i.e., through a spherical surface) contributions were shown to quasi-cancel each other out for all shells.

Our method was then successfully applied to a 3D turbulent convective dynamo case. We initialize a highly non-axisymmetric magnetic field and let the dynamo develop a turbulent spectrum of magnetic energy. We distinguished between the kinematic phase with exponential growth of the magnetic energy spectrum and the nonlinearly saturated phase. The

first phase is dominated by both a non-local transfer of energy from the initial scale of magnetic energy, coupled with the convective scales, toward all of the other magnetic scales, and the shearing by the large-scale differential rotation. A large part of the magnetic energy spectrum is then dictated by the kinetic energy spectrum developed by the convection.

The saturation phase is more subtle and greatly depends on the considered scale in the spectrum. The saturating interactions for the different spectral scales are illustrated in Figure 14 and summarized hereafter.

Our new method allowed us to distinguish between two clear cascades of magnetic energy at the smallest scales of our simulation, for $13 \le L \le 60$ and $L > 60$ (highest $L$'s). In the former case, the differential rotation profile mediates the cascade by shearing the magnetic field. It results in an efficient cascade of magnetic energy.

The latter cascade is also direct and involves all of the highest velocity scales (the large-scale differential rotation does not dominate in this case). It is a generalized cascade over a large range of magnetic scales. The velocity scales involved in the cascade are not local with respect to the magnetic scales. As a consequence, we cannot predict if this generalized cascade would hold at the lowest scales in the case of a real convective dynamo where scale separation is much higher. Besides, the saturation also involves the non-local coupling that can eventually be of the order of the sheared cascade for the intermediate scales. We proved in that case that the transfers cannot be reduced to a limited set of modes.

The saturation of the large-scale axisymmetric dipole and quadrupole appears to be radically different than the small-scale saturation. The toroidal components are mainly saturated
<figure><img src="image_19.png" /><figcaption>**Figure 14.** Saturating interactions of the nonlinear convective dynamo summarizing the results of Section 4. KE and ME stand for kinetic and magnetic energy and DR for differential rotation.</figcaption></figure>

(A color version of this figure is available in the online journal.)

by the balance of the shearing effect of the differential rotation on the large-scale poloidal fields, and the flux transfer through the spherical shell due to the effect of the differential rotation. The poloidal components are mainly saturated by non-local non-axisymmetric interactions. The dipole is saturated by the scale of maximum (highest) magnetic energy, and the quadrupole saturation is not dominated by any particular scale. These two observations point to the two main dependencies of the saturating interactions for the large-scale fields. First, the rotation rate of the star (which is linked to the saturating interaction through the differential rotation) can determine the ability of the dynamo to build wreaths (Brown et al. 2010) and/or to be in a strong or weak regime (Christensen & Aubert 2006; Featherstone et al. 2009; Simitev & Busse 2009). Second, the magnetic Prandtl number $P_m = v/\eta$ determines the position of the peak of magnetic energy and then affects the saturating interaction (e.g., see Schekochihin et al. 2004). We will explore in detail how the saturating interactions depend on those two effects in future work.

Finally, it is worth comparing these results with the previous related work of Livermore et al. (2010). When using forced helical flows and allowing the dynamo field to back-react on the flow, they found that the saturation of the large-scale poloidal dipole was dominated by non-local interaction with a particular magnetic scale ($l = 10, m = 10$ in their case). They showed that magnetic energy was transferred to this scale by the large-scale toroidal magnetic field. In our case, the large-scale dipole is also saturated due to non-local interactions. However, the flow we consider is significantly different because (1) it is obtained from the convective instability and (2) its spectrum is dominated by the large-scale differential rotation that develops self-consistently. Hence, the dynamo process is different and we find that the toroidal large-scale field is saturated by the effect of the large-scale differential rotation, and the large-scale poloidal dipole by the non-local transfer of energy from the magnetic scale of maximal energy.

Our results also suggest that no significant large-scale magnetic field is growing over dissipative timescales in our simulation (the ohmic dissipation timescale for the axisymmetric dipole is typically of the order of $\tau_\eta \sim 700$ days in the simulation). Again, the fast saturation of the dynamo (less than

300 days) may not hold for lower magnetic Prandtl number dynamos.

We developed a diagnostic on the magnetic energy that is an invariant of ideal MHD. In the case of non-ideal MHD, the existence of the selective decay (Taylor 1974; Matthaeus & Montgomery 1980; Mininni & Montgomery 2006) introduces a decoupling between, e.g., the evolution timescales of the magnetic (or total) energy (fast) and the magnetic helicity (slow). As mentioned before, we were interested, in this work, in fast phenomena compared to the ohmic diffusion time. For such processes, the ideal invariants of MHD are still the appropriate quantities to interpret the interactions between scales. The dynamo saturation is necessarily achieved through a modification of the kinetic energy spectrum. As mentioned before, the case we studied in this paper is in the weak branch of the dynamo (i.e., the large-scale poloidal magnetic field does not dominate the magnetic energy spectrum; Christensen & Aubert 2006; Simitev & Busse 2009; Gastine et al. 2012). The detailed modification of the kinetic energy spectrum implied by the saturation of the dynamo process will be discussed for both the strong and weak branches in a future work. Although the magnetic energy is the relevant quantity to characterize nearly kinematic dynamos (where the Lorentz force plays little role), a diagnostic of kinetic energy would be highly valuable for nonlinearly saturated dynamos. In addition to that, the detailed spectral transfers of magnetic helicity are also mandatory to fully address the complexity of the dynamo process. Evolution equations of kinetic energy, magnetic helicity, and cross helicity in the framework introduced in this paper are under development and will be published in a forthcoming paper.

Finally, those diagnostics may also prove very useful for non-dynamo transfer related MHD phenomena. For example, spectral analysis applied to the relaxation and the stability of low-$l$ fossil field (see Braithwaite & Nordlund 2006; Brun 2007; Zahn et al. 2007; Duez & Mathis 2010; Duez et al. 2010) will be studied in a future publication.

We acknowledge inspiring interactions with the participants of the fifth and sixth Festival de Théorie, held in Aix en Provence, 2009 and 2011 July. S. Mathis warmly thanks Dr. N. Tran Minh for fruitful discussions on Racah-Wigner
algebra and its applications. We acknowledge the valuable comments of the anonymous referee that have tightened the focus of this paper. A. Strugarek acknowledges S. Matt for his careful reading of the manuscript. A. S. Brun and A. Strugarek

acknowledge funding by the European Research Council through ERC grant STARS2 207430 (www.stars2.eu). The simulations were performed using HPC resources of GENCI [IDRIS], project 1623.

# APPENDIX A
## DEFINITION AND PROPERTIES OF VECTORIAL SPHERICAL HARMONICS

### A.1. Classical Vectorial Spherical Harmonics Basis

#### A.1.1. Definitions

We define from Rieutord (1987) and Mathis & Zahn (2005):

$$ 
\begin{cases}
\mathbf{R}_l^m(\theta, \varphi) = Y_l^m(\theta, \varphi)\mathbf{e}_r \\
\mathbf{S}_l^m(\theta, \varphi) = \nabla_\perp Y_l^m = \partial_\theta Y_l^m \mathbf{e}_\theta + \frac{1}{\sin\theta} \partial_\varphi Y_l^m \mathbf{e}_\varphi \\
\mathbf{T}_l^m(\theta, \varphi) = \nabla_\perp \times \mathbf{R}_l^m = \frac{1}{\sin\theta} \partial_\varphi Y_l^m \mathbf{e}_\theta - \partial_\theta Y_l^m \mathbf{e}_\varphi
\end{cases}
\quad (A1) $$

where $(\mathbf{e}_r, \mathbf{e}_\theta, \mathbf{e}_\varphi)$ defines the spherical basis and $Y_l^m$ are the Laplace spherical harmonics defined by

$$ Y_l^m(\theta, \varphi) = (-1)^{\frac{m+|m|}{2}} \sqrt{\frac{(2l+1)(l-|m|)!}{4\pi(l+|m|)!}} P_l^{|m|}(\cos\theta)e^{im\varphi}, \quad (A2) $$

where $P_l^m$ are the associated Legendre polynomials. The basis Equations ((A1)) have the following properties:

$$ \int_S \mathbf{R}_{l_1}^{m_1} \cdot (\mathbf{R}_{l_2}^{m_2})^{cc} \, d\Omega = \delta_{l_1, l_2} \delta_{m_1, m_2}, \quad (A3) $$

$$ \int_S \mathbf{S}_{l_1}^{m_1} \cdot (\mathbf{S}_{l_2}^{m_2})^{cc} \, d\Omega = \int_S \mathbf{T}_{l_1}^{m_1} \cdot (\mathbf{T}_{l_2}^{m_2})^{cc} \, d\Omega = l_1(l_1+1)\delta_{l_1, l_2} \delta_{m_1, m_2}, \quad (A4) $$

where $S$ is a spherical surface, $d\Omega = \sin\theta d\theta d\varphi$ the solid angle, $cc$ means complex conjugate, and $\delta$ is the Kronecker symbol. We also have

$$ (\mathbf{S}_l^m)^{cc} = (-1)^m \mathbf{S}_l^{-m}, \quad (A5) $$

and all of the other scalar cross products are 0. We remind the reader that the poloidal fields are described by their projection on $(\mathbf{R}_l^m, \mathbf{S}_l^m)$, and the toroidal fields by their projection on $\mathbf{T}_l^m$.

#### A.1.2. Scalar Fields

Defining $\psi(r, \theta, \varphi) = \sum_{l=1}^\infty \sum_{m=-l}^l \{\psi_m^l(r) Y_l^m(\theta, \varphi)\}$, we get

$$ \nabla \psi = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \left\{ \partial_r \psi_m^l \mathbf{R}_l^m + \frac{\psi_m^l}{r} \mathbf{S}_l^m \right\}, \quad (A6) $$

$$ \nabla \cdot \nabla \psi = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \Delta_l \psi_m^l Y_l^m \quad (A7) $$

where $\Delta_l = \partial_{rr}^2 + (2/r)\partial_r - l(l+1)/r^2$.

#### A.1.3. Vectorial Fields

For a vector $\mathbf{X}(r, \theta, \varphi) = \sum_{l=1}^\infty \sum_{m=-l}^l \{A_m^l(r)\mathbf{R}_l^m + B_m^l(r)\mathbf{S}_l^m + C_m^l(r)\mathbf{T}_l^m\}$, we obtain

$$ \nabla \cdot \mathbf{X} = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \left\{ \left[ \frac{1}{r^2} \partial_r (r^2 A_m^l) - l(l+1) \frac{B_m^l}{r} \right] Y_l^m \right\}, \quad (A8) $$

$$ \nabla \times \mathbf{X} = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \left\{ \left[ l(l+1) \frac{C_m^l}{r} \right] \mathbf{R}_l^m + \left[ \frac{1}{r} \partial_r (r C_m^l) \right] \mathbf{S}_l^m + \left[ \frac{A_m^l}{r} - \frac{1}{r} \partial_r (r B_m^l) \right] \mathbf{T}_l^m \right\}, \quad (A9) $$

$$ \nabla^2 \mathbf{X} = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \left\{ \left[ \Delta_l A_m^l - \frac{2}{r^2} (A_m^l - l(l+1)B_m^l) \right] \mathbf{R}_l^m + \left[ \Delta_l B_m^l + 2 \frac{A_m^l}{r^2} \right] \mathbf{S}_l^m + \left[ \Delta_l C_m^l \right] \mathbf{T}_l^m \right\}. \quad (A10) $$
### A.1.4. Recurrence Relations

In addition to the expression of the different operators, we also give here two useful coupling relations between spherical harmonics. First, according to Varshalovich et al. (1975), the coupling between cos $\theta$ and the spherical harmonic $Y_l^m$ is given by

$$ \cos \theta Y_l^m = \sqrt{\frac{(l-m+1)(l+m+1)}{(2l+1)(2l+3)}} Y_{l+1}^m + \sqrt{\frac{(l-m)(l+m)}{(2l-1)(2l+1)}} Y_{l-1}^m. \quad (A11) $$

Then, one easily deduces the following properties:

$$ \cos \theta \mathbf{S}_l^m = \frac{l}{l+1} \sqrt{\frac{(l-m+1)(l+m+1)}{(2l+1)(2l+3)}} \mathbf{S}_{l+1}^m + \frac{l+1}{l} \sqrt{\frac{(l-m)(l+m)}{(2l-1)(2l+1)}} \mathbf{S}_{l-1}^m - \frac{im}{l(l+1)} \mathbf{T}_l^m, \quad (A12) $$

$$ \cos \theta \mathbf{T}_l^m = \frac{l}{l+1} \sqrt{\frac{(l-m+1)(l+m+1)}{(2l+1)(2l+3)}} \mathbf{T}_{l+1}^m + \frac{l+1}{l} \sqrt{\frac{(l-m)(l+m)}{(2l-1)(2l+1)}} \mathbf{T}_{l-1}^m + \frac{im}{l(l+1)} \mathbf{S}_l^m. \quad (A13) $$

## A.2. An Alternative Vectorial Basis

### A.2.1. Definitions

The vectorial spherical harmonics basis defined in Appendix A.1 is very efficient for calculating scalar products or a linear differential operator on vectors. Nevertheless, it is quite hard to use it to express vectorial products. Instead, we define the following basis (e.g., see Varshalovich et al. 1975):

$$ \mathbf{Y}_{l,l+\nu}^{m}(\theta, \varphi) = \sum_{\mu=-1}^{1} \left\{ (-1)^{l-m} \sqrt{2l+1} \begin{pmatrix} l & l+\nu & 1 \\ m & \mu-m & -\mu \end{pmatrix} Y_{l+\nu}^{m-\mu} \mathbf{e}_{\mu} \right\}, \quad (A14) $$

where (...) is the 3j Wigner coefficient linked to Clebsch-Gordan coefficients, and the vectors $\mathbf{e}_{\mu}$ are

$$ \left\{ \begin{aligned} \mathbf{e}_{-1} &= \frac{1}{\sqrt{2}}(\mathbf{e}_x - i\mathbf{e}_y) \\ \mathbf{e}_0 &= \mathbf{e}_z \\ \mathbf{e}_1 &= -\frac{1}{\sqrt{2}}(\mathbf{e}_x + i\mathbf{e}_y) \end{aligned} \right. , \quad (A15) $$

where $(\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z)$ defines the Cartesian basis. Note that the equivalent of the conjugation rule (Equation (A5)) is then

$$ (\mathbf{Y}_{l,l+\nu}^{m})^{cc} = (-1)^{m+\delta_{0\nu}} \mathbf{Y}_{l,l+\nu}^{-m}. \quad (A16) $$

Again, we recall that the poloidal fields are described by their projection on $(\mathbf{Y}_{l,l+1}^m, \mathbf{Y}_{l,l-1}^m)$ ($\nu \in \{-1; 1\}$), and the toroidal fields are described by their projection on $\mathbf{Y}_{l,l}^m$ ($\nu = 0$).

### A.2.2. Vectorial Product

We decompose a vector $\mathbf{X}$ on this basis in the following way:

$$ \mathbf{X}(r, \theta, \varphi) = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \sum_{\nu=-1}^{1} X_{l,l+\nu}^{m}(r) \mathbf{Y}_{l,l+\nu}^{m}. $$

Evaluating the vectorial product of two vectors $\mathbf{X}_{12} = \mathbf{X}_1 \times \mathbf{X}_2$, one gets

$$ X_{12; l_1 l_2 l_1 l_2}^{m_1 m_2} = \sum_{\substack{l_1, l_2 = 1 \\ l_1 l_2 \ge l_1 - l_2 \\ l_1 l_2 \le l_1 + l_2}}^{\infty} \sum_{\substack{m_1 = -l_1 \\ m_2 = -l_2 \\ m_1 + m_2 = m_2}}^{l_1 l_2} \sum_{\nu_1, \nu_2} X_{1; l_1 l_1 + \nu_1}^{m_1} X_{2; l_2 l_2 + \nu_2}^{m_2} \mathcal{J}_{l_1, m_1; \nu_1, l_2, m_2; \nu_2}^{l_1 l_2, m_1 m_2; \nu_1 \nu_2}, \quad (A17) $$

where

$$ \begin{aligned} \mathcal{J}_{l_1, m_1; \nu_1, l_2, m_2; \nu_2}^{l_1, m_1 + m_2; \nu} &= i(-1)^{\nu_1 - \nu_2 + (m_1 + m_2)} \\ & \sqrt{\frac{3}{2\pi}(2l_1 + 1)(2l_1 + 2\nu_1 + 1)(2l_2 + 1)(2l_2 + 2\nu_2 + 1)(2l + 1)(2l + 2\nu + 1)} \\ & \left\{ \begin{pmatrix} l_1 & l_2 & l \\ l_1 + \nu_1 & l_2 + \nu_2 & l + \nu \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l \\ m_1 & m_2 & -(m_1 + m_2) \end{pmatrix} \begin{pmatrix} l_1 + \nu_1 & l_2 + \nu_2 & l + \nu \\ 0 & 0 & 0 \end{pmatrix}, \right. \end{aligned} \quad (A18) $$

with $\{\cdots\}$ being the 9j Wigner coefficient.
### A.3. Basis Change Relations

For a vector **X** decomposed in the following manner:

$$ 
\begin{aligned}
\mathbf{X} &= \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \{ \mathcal{A}_m^l \mathbf{R}_l^m + \mathcal{B}_m^l \mathbf{S}_l^m + \mathcal{C}_m^l \mathbf{T}_l^m \} \\
&= \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \sum_{\nu=-1}^{1} \{ X_{l,l+\nu}^m \mathbf{Y}_{l,l+\nu}^m \},
\end{aligned}
$$

we have the two following relations to change from one basis to the other:

$$ 
\left\{ 
\begin{aligned}
\mathcal{A}_m^l &= \frac{1}{\sqrt{2l+1}} \left[ \sqrt{l} X_{l,l-1}^m - \sqrt{l+1} X_{l,l+1}^m \right] \\
\mathcal{B}_m^l &= \frac{1}{\sqrt{2l+1}} \left[ \frac{1}{\sqrt{l}} X_{l,l-1}^m + \frac{1}{\sqrt{l+1}} X_{l,l+1}^m \right] \\
\mathcal{C}_m^l &= \frac{i}{\sqrt{l(l+1)}} X_{l,l}^m
\end{aligned}
\right. , \quad 
\left\{ 
\begin{aligned}
X_{l,l-1}^m &= \sqrt{\frac{l}{2l+1}} (\mathcal{A}_m^l + (l+1)\mathcal{B}_m^l) \\
X_{l,l}^m &= -i\sqrt{l(l+1)}\mathcal{C}_m^l \\
X_{l,l+1}^m &= \sqrt{\frac{l+1}{2l+1}} (-\mathcal{A}_m^l + l\mathcal{B}_m^l)
\end{aligned}
\right. . \qquad (A19)
$$

### A.4. Expression of the α Effect

The α effect introduces the spectral coupling of a scalar and a vector, which was not treated before. In the special case of an axisymmetric α and an axisymmetric vectorial field $B_\phi \mathbf{e}_\phi$ (which is the case in this paper, see Section 3), we write the α coefficient

$$ \alpha(r, \theta) = \sum_{l=0}^{+\infty} \alpha_0^l(r) Y_l^0, $$

and we rewrite the magnetic field from Equation (11)

$$ B_{\phi} \mathbf{e}_{\phi} = \sum_{l=0}^{+\infty} \frac{A_0^l}{r} \mathbf{T}_l^0 . $$

By introducing the coefficient

$$ \mathcal{H}_{l_1, l_2}^l = -\sqrt{\frac{1}{4\pi} l(l+1)l_2(l_2+1)(2l+1)(2l_1+1)(2l_2+1)} \begin{pmatrix} l_1 & l_2 & l \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l_1 & l_2 & l \\ 0 & 1 & -1 \end{pmatrix}, \qquad (A20) $$

we can write the α effect as

$$ \alpha B_{\phi} \mathbf{e}_{\phi} = \sum_{l=0}^{+\infty} \left( \sum_{\substack{l_1, l_2 = 0 \\ l \ge |l_1 - l_2| \\ l \le l_1 + l_2}}^{+\infty} \mathcal{H}_{l_1, l_2}^l \alpha_{l_1}^0 \frac{A_{l_2}^0}{r} \right) \mathbf{T}_l^0 = \sum_{l=0}^{+\infty} \frac{\xi_0^l}{r} \mathbf{T}_l^0 . \qquad (A21) $$

Finally, one gets

$$ \nabla \times (\alpha B_{\phi} \mathbf{e}_{\phi}) = \sum_{l=0}^{+\infty} \left\{ \frac{l(l+1)}{r^2} \xi_0^l \mathbf{R}_l^0 + \frac{1}{r} \partial_r (\xi_0^l) \mathbf{S}_l^0 \right\} . \qquad (A22) $$

### A.5. Couplings in the Magnetic Energy Equation

The detailed expressions of the different terms of the magnetic energy equation (Equations (17)–(19)) are given here. We write the magnetic field **B** and the current **J** as

$$ \mathbf{B} = \sum_{l=1}^{+\infty} \sum_{m=-l}^{l} \sum_{\nu=-1}^{1} B_{l,l+\nu}^{m}(r) \mathbf{Y}_{l,l+\nu}^{m}, \qquad (A23) $$

$$ \mathbf{J} = \sum_{l=1}^{+\infty} \sum_{m=-l}^{l} \sum_{\nu=-1}^{1} J_{l,l+\nu}^{m}(r) \mathbf{Y}_{l,l+\nu}^{m}. \qquad (A24) $$
Based on this, the vectorial product may be evaluated thanks to a coupling coefficient $J_{l_1,m_1,v_1,l_2,m_2,v_2}^{l,m,v}$ given in Equation (A18). The transformation rules from one basis to the other are given in Appendix A.3, and they allow us to easily evaluate the integrals (17)–(19). By separating the diffusive terms into two $D_1$ and $D_2$ terms, we get

$$D_1(L, r) = \eta \sum_L (-1)^m l(l+1) \left\{ \frac{l(l+1)}{r^3} \Delta_l \left(\frac{C_m^l}{r}\right) C_{-m}^{l} + \frac{1}{r^2} \partial_r \left[ r \Delta_l \left(\frac{C_m^l}{r}\right) \right] \partial_r C_{-m}^{l} + \Delta_l \left(\frac{A_m^l}{r}\right) \frac{A_{-m}^l}{r} \right\}, \quad (A25)$$

$$D_2(L, r) = - \frac{c}{\sqrt{4\pi}} \sum_L \sum_{v_1, v_2} B_{l,l+v_1}^m J_{l,l+v_2}^{-m} J_{l,m,v_1,l,-m,v_2}^{0,0,1}, \quad (A26)$$

where $\sum_L$ stands for a summation over all of the spherical harmonics contained in the shell L (one element in an axisymmetric shell, and 2l elements in a non-axisymmetric shell). The production and flux terms then read

$$\mathcal{P}_L(L_1, L_2, r) = \frac{c}{4\pi} \sum_L \frac{(-1)^m}{r^2} \frac{l(l+1)}{\sqrt{2l+1}} \left\{ \left[ \sqrt{l} (\mathbf{U}_{L_1} \times \mathbf{B}_{L_2})_{l,l-1}^m - \sqrt{l+1} (\mathbf{U}_{L_1} \times \mathbf{B}_{L_2})_{l,l+1}^m \right] A_{-m}^l \right. \\ \left. + \left[ \frac{1}{\sqrt{l}} (\mathbf{U}_{L_1} \times \mathbf{B}_{L_2})_{l,l-1}^m + \frac{1}{\sqrt{l+1}} (\mathbf{U}_{L_1} \times \mathbf{B}_{L_2})_{l,l+1}^m \right] r \partial_r A_{-m}^l - i \sqrt{\frac{2l+1}{l(l+1)}} (\mathbf{U}_{L_1} \times \mathbf{B}_{L_2})_{l,l}^m r^2 \Delta_l \left(\frac{C_{-m}^l}{r}\right) \right\}, \quad (A27)$$

$$\mathcal{F}_L(L_1, L_2, r) = -\frac{\sqrt{4\pi}}{r^2} \partial_r \left\{ r^2 \sum_L \sum_{v_1=-1}^1 \sum_{v_2=-1}^1 (\mathbf{U}_{L_1} \times \mathbf{B}_{L_2})_{l,l+v_1}^m B_{l,l+v_2}^{-m} J_{l,m,v_1,l,-m,v_2}^{0,0,1} \right\}. \quad (A28)$$

The Laplacian formula used for $D_1$ in the vectorial spherical harmonics basis is given in Equation (A10). In the production term, we simply made use of the basis transformation (Equation (A19)). Finally, the expressions for $D_2$ and the flux term need some intermediate steps to be properly explained. These details are given in Appendix A.6 for the flux of magnetic energy, and the same procedure may be applied in the case of the second diffusive term.

### A.6. Simplification of the Magnetic Energy Flux

The flux of magnetic energy can be simplified, if one notes that it has the general form

$$F = \int_S \nabla \cdot \mathbf{X} \, d\Omega \quad \text{and} \quad \mathbf{X} = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \{ A_m^l \mathbf{R}_l^m + B_m^l \mathbf{S}_l^m + C_m^l \mathbf{T}_l^m \}.$$

Then, one can easily deduce that

$$F = \int_S \nabla \cdot \mathbf{X} \, d\Omega = \frac{\sqrt{4\pi}}{r^2} \partial_r (r^2 A_0^0).$$

Recalling from the system (Equation (A19)) that $\mathbf{R}_0^0 = -\mathbf{Y}_{0,1}^0$, and if one assumes that $\mathbf{X} = \mathbf{X}_1 \times \mathbf{X}_2$, then one obtains, for an integral similar to the magnetic energy flux,

$$\int_S \nabla \cdot (\mathbf{X}_1 \times \mathbf{X}_2) \, d\Omega = -\frac{\sqrt{4\pi}}{r^2} \partial_r (r^2 X_{0,1}^0) = -\frac{\sqrt{4\pi}}{r^2} \partial_r \left\{ r^2 \sum_{l=0}^{+\infty} \sum_{m=-l}^{l} \sum_{v_1, v_2} X_{1;l,l+v_1}^m X_{2;l,l+v_2}^{-m} J_{l,m,v_1,l,-m,v_2}^{0,0,1} \right\}.$$

### A.7. On Primary and Secondary Families

Previous studies on dynamos in stars shed light on the important distinction of primary (or dipolar, antisymmetric with respect to the equator) and secondary (or quadrupolar, symmetric with respect to the equator) families of the magnetic field. For a vector $\mathbf{X} = \sum_{l=1}^{\infty} \sum_{m=-l}^{l} \{A_m^l \mathbf{R}_l^m + B_m^l \mathbf{S}_l^m + C_m^l \mathbf{T}_l^m\}$, Roberts & Stix (1972) define the primary family as

$$\mathbf{X}^p = A_m^{m+1} \mathbf{R}_{m+1}^m + B_m^{m+1} \mathbf{S}_{m+1}^m + C_m^{m+1} \mathbf{T}_{m+1}^m + A_m^{m+3} \mathbf{R}_{m+3}^m + B_m^{m+3} \mathbf{S}_{m+3}^m + C_m^{m+2} \mathbf{T}_{m+2}^m + \dots$$

and the secondary family as

$$\mathbf{X}^s = A_m^m \mathbf{R}_m^m + B_m^m \mathbf{S}_m^m + C_m^{m+1} \mathbf{T}_m^m + A_m^{m+2} \mathbf{R}_{m+2}^m + B_m^{m+2} \mathbf{S}_{m+2}^m + C_m^{m+3} \mathbf{T}_{m+3}^m + \dots$$

It can also be easily shown that in the $(\mathbf{Y}_{l,l+v}^m)_{v=-1,0,1}$ basis, a primary field always satisfies $l+m+v$ as even, and a secondary field always satisfies $l+m+v$ as odd.
Note that the vectorial product depends on the 3j Wigner

$$ \begin{pmatrix} l_1 + \nu_1 & l_2 + \nu_2 & l + \nu \\ 0 & 0 & 0 \end{pmatrix}. $$

Recalling that $m_1 + m_2 = m$, this 3j Wigner is zero if $l_1 + \nu_1 + l_2 + \nu_2 + l + \nu = (l_1 + \nu_1 + m_1) + (l_2 + \nu_2 + m_2) + (l + \nu + m) - 2m$ is odd. Consequently, in order to have a non-zero 3j Wigner, if U and B are from different families, then their vectorial product is a secondary field; and if they are from the same family, then their vectorial product is a primary field. If $C = A \times B$, then this means that

$$ \left. \begin{array}{l} \mathbf{A}^p \times \mathbf{B}^s \\ \mathbf{A}^s \times \mathbf{B}^p \end{array} \right\} \rightarrow \mathbf{C}^s \quad \text{and} \quad \left. \begin{array}{l} \mathbf{A}^p \times \mathbf{B}^p \\ \mathbf{A}^s \times \mathbf{B}^s \end{array} \right\} \rightarrow \mathbf{C}^p, \qquad (A29) $$

where the superscripts p and s stand for primary and secondary. This was already acknowledged by McFadden et al. (1991) and Gubbins & Zhang (1993).

## APPENDIX B

### NUMERICAL VALIDATION

In order to validate the way we implemented the complex interactions between spherical harmonics in the ASH code, we compared an analytic calculation for a simple setup with numerical results. We summarize here those calculations.

We start from a mixed $m = 0, \pm 1$ and $l = 1$ state for the magnetic field, and an $(l, m) = (2, \pm 1)$ state for the velocity field. We initialize the magnetic field in the following way:

$$ \mathbf{B} = \frac{2R_{\odot}^2 R_b}{r^3} \left( \mathbf{R}_1^0 + \frac{1}{2}\mathbf{R}_1^1 - \frac{1}{2}\mathbf{R}_1^{-1} \right) - \frac{R_{\odot}^2 R_b}{r^3} \left( \mathbf{S}_1^0 + \frac{1}{2}\mathbf{S}_1^1 - \frac{1}{2}\mathbf{S}_1^{-1} \right) + \frac{\beta^2}{r R_{\odot}^2} \left( \mathbf{T}_1^0 + \frac{1}{2}\mathbf{T}_1^1 - \frac{1}{2}\mathbf{T}_1^{-1} \right), $$

where $R_{\odot}$ is the solar radius, $R_b$ is our inner boundary radius, $R_t$ is our outer boundary radius, and $\beta = (R_t - r)^2(r - R_b)^2$. The velocity is initialized by

$$ \bar{\rho} \mathbf{U} = \frac{3\beta^2}{r^2 R_{\odot}^4} (\mathbf{R}_2^1 - \mathbf{R}_2^{-1}) + \frac{3\chi r}{2R_{\odot}^2} (\mathbf{S}_2^1 - \mathbf{S}_2^{-1}) + \frac{\beta^2}{2r R_{\odot}^2} (\mathbf{T}_2^1 - \mathbf{T}_2^{-1}), $$

where $\chi = (R_t - r)^2(r - R_b)^2(R_t + R_b - 2r)/(r^2 R_{\odot}^2)$. By rewriting those fields in the conventional spherical harmonics method of writing, we may calculate the axisymmetric components of the vectorial product $\bar{\rho} \mathbf{U} \times \mathbf{B}$ and obtain

$$ \bar{\rho} \mathbf{U} \times \mathbf{B} = \begin{pmatrix} -\frac{2\gamma_t}{5} \sqrt{\frac{4\pi}{7}} Y_3^0 - \frac{3\gamma_t}{5} \sqrt{\frac{4\pi}{3}} Y_1^0 \\ -\frac{2\gamma_t^2}{15} \sqrt{\frac{4\pi}{7}} \partial_{\theta} Y_3^0 - \left(\frac{\gamma_t^2}{5} + \gamma_{\theta}^1\right) \sqrt{\frac{4\pi}{3}} \partial_{\theta} Y_1^0 \\ -\frac{2\gamma_{\phi}^2}{15} \sqrt{\frac{4\pi}{7}} \partial_{\theta} Y_3^0 - \left(\frac{\gamma_{\phi}^2}{5} + \gamma_{\phi}^1\right) \sqrt{\frac{4\pi}{3}} \partial_{\theta} Y_1^0 \end{pmatrix}, \qquad (B1) $$

where the $\gamma$ coefficients are defined by

$$ \gamma_r = \frac{9\sqrt{5}\chi\beta^2}{8\pi R_{\odot}^2} + \gamma_{\theta}^1, $$

$$ \gamma_{\theta}^1 = \frac{3\sqrt{5}\beta^2 R_b}{8\pi r^4}, \quad \gamma_{\theta}^2 = \frac{9\sqrt{5}\beta^5}{8\pi r^3 R_{\odot}^6} - \gamma_{\theta}^1 $$

$$ \gamma_{\phi}^1 = \frac{9\sqrt{5}\chi R_b}{8\pi r^2} \quad \text{and} \quad \gamma_{\phi}^2 = -\frac{9\sqrt{5}\beta^3 R_b}{8\pi r^5 R_{\odot}^2} - \gamma_{\phi}^1. $$

These coefficients match exactly the outputs from the code (Table 1). The production and flux terms in Equation (16) are then simple scalar products involving the vectorial product (Equation (B1)). They also have been checked by comparison with the analytical calculation.

## REFERENCES

* Alexakis, A., Mininni, P. D., & Pouquet, A. 2005, PhRvE, 72, 46301

* Aluie, H., & Eyink, G. L. 2010, PRL, 104, 81101

* Biskamp, D. 1993, Nonlinear Magnetohydrodynamics (Cambridge Monographs on Plasma Physics; New York: Cambridge Univ. Press)

* Blackman, E. G., & Brandenburg, A. 2002, ApJ, 579, 359

* Boldyrev, S., Mason, J., & Cattaneo, F. 2009, ApJL, 699, L39

* Bonanno, A., Elstner, D., Rüdiger, G., & Belvedere, G. 2002, A&A, 390, 673

* Braithwaite, J., & Nordlund, Å. 2006, A&A, 450, 1077

* Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2010, ApJ, 711, 424
Browning, M. K. 2008, ApJ, 676, 1262

Brun, A. S. 2007, AN, 328, 1137

Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073

Bullard, E., & Gellman, H. 1954, RSPTA, 247, 213

Cattaneo, F., & Hughes, D. W. 2001, A&G, 42, 18

Centeno, R., Socas-Navarro, H., Lites, B., et al. 2007, ApJL, 666, L137

Charbonneau, P. 2010, LRSP, 7, 3

Charbonneau, P., & MacGregor, K. B. 1997, ApJ, 486, 502

Christensen, U. R., & Aubert, J. 2006, GeoJI, 166, 97

Clune, T. L., Elliott, J. R., Miesch, M. S., Toomre, J., & Glatzmaier, G. A. 1999, ParC, 25, 361

Dar, G., Verma, M. K., & Eswaran, V. 2001, PhyD, 157, 207

Debliquy, O., Verma, M. K., & Carati, D. 2005, PhPl, 12, 2309

DeRosa, M. L., Brun, A. S., & Hoeksema, J. T. 2011, in IAU Symp. 271, Astrophysical Dynamics: From Stars to Galaxies, ed. N. H. Brummell, A. S. Brun, M. S. Miesch, & Y. Ponty (Cambridge: Cambridge Univ. Press), 94

DeRosa, M. L., Brun, A. S., & Hoeksema, J. T. 2012, ApJ, 757, 96

Donati, J.-F., & Landstreet, J. D. 2009, ARA&A, 47, 333

Duez, V., Braithwaite, J., & Mathis, S. 2010, ApJL, 724, L34

Duez, V., & Mathis, S. 2010, A&A, 517, 58

Farge, M. 1992, AnRFM, 24, 395

Featherstone, N. A., Browning, M. K., Brun, A. S., & Toomre, J. 2009, ApJ, 705, 1000

Frick, P., & Sokoloff, D. 1998, PhRvE, 57, 4155

Frisch, U. 1995, Turbulence: The Legacy of A. N. Kolmogorov (Cambridge University Press, Cambridge)

Frisch, U., Pouquet, A., Leorat, J., & Mazure, A. 1975, JFM, 68, 769

Gastine, T., Duarte, L., & Wicht, J. 2012, A&A, 546, 19

Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763

Gubbins, D., & Zhang, K. 1993, PEPI, 75, 225

Hagenaar, H. J., Schrijver, C. J., & Title, A. M. 2003, ApJ, 584, 1107

Hughes, D. W., & Proctor, M. R. E. 2012, JFM, arXiv:1211.5339

Iroshnikov, P. S. 1964, SvA, 7, 566

Ivers, D. J., & Phillips, C. G. 2008, GeoJI, 175, 955

Jones, C. A., Boronski, P., Brun, A. S., et al. 2011, Icar, 216, 120

Jouve, L., & Brun, A. S. 2009, ApJ, 701, 1300

Jouve, L., Brun, A. S., Arlt, R., et al. 2008, A&A, 483, 949

Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012, ApJL, 755, L22

Knobloch, E., Tobias, S. M., & Weiss, N. O. 1998, MNRAS, 297, 1123

Kraichnan, R. H. 1965, PhFI, 8, 1385

Lesieur, M. 2008, Turbulence in Fluids (Berlin: Springer)

Lesur, G., & Longaretti, P.-Y. 2011, A&A, 528, 17

Livermore, P. W., Hughes, D. W., & Tobias, S. M. 2010, PhFI, 22, 7101

Maron, J., Cowley, S., & McWilliams, J. 2004, ApJ, 603, 569

Mathis, S., & Zahn, J.-P. 2005, A&A, 440, 653

Matthaeus, W. H., & Montgomery, D. 1980, NYASA, 357, 203

McFadden, P. L., Merrill, R. T., McElhinny, M. W., & Lee, S. 1991, JGR, 96, 3923

Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 641, 618

Mininni, P. D., Alexakis, A., & Pouquet, A. 2005, PhRvE, 72, 46302

Mininni, P. D., & Montgomery, D. C. 2006, PhFI, 18, 6602

Moffatt, H. K. 1978, Magnetic Field Generation in Electrically Conducting Fluids (Bristol: Cambridge)

Olson, P., Christensen, U., & Glatzmaier, G. A. 1999, JGR, 104, 10383

Pinto, R. F., & Brun, A. S. 2012, ApJ, submitted

Politano, H., & Pouquet, A. 1998, GeoRL, 25, 273

Pouquet, A., Brachet, M.-E., Lee, E., et al. 2011, in IAU Symp. 271, Astrophysical Dynamics: From Stars to Galaxies, ed. N. H. Brummell, A. S. Brun, M. S. Miesch, & Y. Ponty (Cambridge: Cambridge Univ. Press), 304

Pouquet, A., Frisch, U., & Leorat, J. 1976, JFM, 77, 321

Racine, É., Charbonneau, P., Ghizaru, M., Bouchat, A., & Smolarkiewicz, P. 2011, ApJ, 735, 46

Rieutord, M. 1987, GApFD, 39, 163

Rincon, F. 2006, JFM, 563, 43

Roberts, P. H., & Stix, M. 1972, A&A, 18, 453

Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., & McWilliams, J. C. 2004, ApJ, 612, 276

Schilling, O., & Zhou, Y. 2002, JPIPh, 68, 389

Schou, J., Antia, H. M., Basu, S., et al. 1998, ApJ, 505, 390

Schrijver, C., & DeRosa, M. 2003, SoPh, 212, 165

Simitev, R. D., & Busse, F. H. 2009, EPL, 85, 19001

Taylor, J. B. 1974, PRL, 33, 1139

Tobias, S. M. 2002, RSPTA, 360, 2741

Tobias, S. M., & Cattaneo, F. 2008a, JFM, 601, 101

Tobias, S. M., & Cattaneo, F. 2008b, PRL, 101, 125003

Tobias, S. M., Cattaneo, F., & Brummell, N. H. 2011, ApJ, 728, 153

Varshalovich, D. A., Moskalev, A. N., & Khersonskii, V. K. 1975, Quantum Theory of Angular Momentum (Singapore: World Scientific), 514

Verma, M. K., Ayyer, A., & Chandra, A. V. 2005, PoP, 12, 2307

Zahn, J.-P., Brun, A. S., & Mathis, S. 2007, A&A, 474, 145

Zhang, K., Chan, K. H., Zou, J., Liao, X., & Schubert, G. 2003, ApJ, 596, 663