© ESO 2019

1. Introduction

Recent increased interest in the large-scale magnetic fields of the Milky Way (MW) and other spiral galaxies is driven by a number of factors. Their role in the dynamics of the interstellar medium (ISM) has been widely appreciated although not completely understood. Their significance for the feedback processes in evolving galaxies has also been recognised and is being actively explored (e.g. Hopkins et al. 2018). Furthermore, the separation of the Galactic foreground, including its magnetic field and associated emission, from extragalactic contributions is essential to identify the sources of ultra-high energy cosmic rays (UHECR; Kotera & Olinto 2011; Mollerach & Roulet 2018) and for cosmological interpretations of sensitive CMB observations (e.g. Planck Collaboration XI 2019).

Our understanding of the Galactic magnetic field (GMF) is based on observations of Faraday rotation of polarised radio emission, synchrotron emission of energetic electrons, and polarised emission and absorption by dust (e.g. Haverkorn 2015). Observations of external spiral galaxies have provided rich polarisation data for their discs and haloes (Beck 2015; Wiegert et al. 2015; Mao et al. 2015). Interpretation of such data in terms of three-dimensional magnetic field structures requires parametrised models of the magnetic field based on the understanding of their nature and origin. It can be expected that models that rely less on specific theoretical models would have a larger number of free parameters whose physical meaning would be less clear. On the contrary, physically motivated models can be more flexible and lead to physically transparent interpretations of observations.

A range of heuristic models for the structure of the GMF have been proposed (e.g. Sun et al. 2008; Jaffe et al. 2010; Van Eck et al. 2011; Jansson & Farrar 2012a,b; Ferrière & Terral 2014; Terral & Ferrière 2017). Magnetic fields in most of these models are superpositions of various ad hoc parts whose parameters are selected by fitting to a range of observables. The heuristic nature of the models makes them rather inflexible. Furthermore, their parameters are not necessarily related to the ISM properties (often lacking a transparent physical meaning), and some such models fail to satisfy even such fundamental constraints as ∇ ⋅ B = 0 or to allow for the global helical nature of galactic magnetic fields that imprints on its structure and symmetries. The ambiguities, problems and limitations of the current approaches to GMF modelling have been discussed in Planck Collaboration XLII (2016) and by the IMAGINE consortium (Boulanger et al. 2018).

A possible way to overcome heuristics in GMF modelling would be to extract the field structure from physical simulations of galaxy formation and evolution, which include all relevant processes of magnetic field generation and are specific to the MW. Such simulations of generic galaxies have led to important insights into various magnetic structures in spiral galaxies (e.g. Pakmor et al. 2017; Pakmor & Springel 2013) but their resolution remains insufficient to capture galactic dynamo action, as it is controlled by turbulent processes on scales less than 100 pc. Moreover, the simulations would have to be constrained in a way flexible enough to reproduce an ever increasing set of observational data.

Here we propose a simple (that is, flexible, adjustable, analytic) and yet realistic (i.e., based on relevant equations and observations) approach to model the large-scale (mean) magnetic field in the disc and halo of the Milky Way and other spiral galaxies. The model has been implemented as a publicly available software package GALMAG (Rodrigues 2018), which can be used in the framework of Bayesian optimisers (e.g. Steininger et al. 2018; Steininger 2018).

The text is organised as follows. In Sect. 2, we lay out the modelling strategy, basic equations and fundamental assumptions. In Sect. 3, we describe the solutions for magnetic field in the disc and in Sect. 4, for the halo. In Sect. 5 we discuss possible applications of this model to the interpretation of observations of synchrotron emission and Faraday rotation (Sect. 5.1) as well some of the ways in which it could be used to model the evolution of galactic magnetic fields (Sect. 5.2); possible extensions to this model are also discussed (Sect. 5.3), whilst Sect. 5.4 introduces the publicly available GALMAG software package that implements the model. Our results are summarised in Sect. 6.

2. Basic equations

A physically meaningful model of a galactic magnetic field has to rely on a clear physical picture of its origin and maintenance, as well as on a specific galaxy model. Mean-field dynamo action is the most plausible mechanism of generation and maintenance of large-scale magnetic fields in spiral galaxies (Ruzmaikin et al. 1988; Beck et al. 1996; Brandenburg & Subramanian 2005; Beck 2015). Therefore, we first explore the nature of magnetic structures compatible with dynamo action. However, magnetic fields observed in galaxies have many features that emerge independently of the dynamo process. To allow for such features, we make our model rather independent of the specific properties of dynamo-generated magnetic fields and use solutions of the dynamo equations just as a convenient functional basis to parametrise a wide class of magnetic structures. Furthermore, we discuss how an even more general type of the governing equations can be used for these purposes.

The dynamo converts kinetic energy of random (turbulent) flows into magnetic energy. The mean helicity of the random flow, that emerges because of the overall rotation and stratification, leads to the generation of magnetic fields at scales much larger than the correlation scale of the random flow (this is described as the α-effect). Differential rotation can accelerate the energy conversion by stretching the large-scale magnetic fields in the direction of the flow (the ω-effect). The spatial scale and structure of the magnetic field are controlled by the mean-field transport coefficients, quantifying the averaged induction effects of the random flows, and the large-scale velocity shear rate. The magnetic field structure also depends upon the geometric shape of the dynamo region (e.g. spherical, toroidal or flat). When the magnetic field is weak, so that its effect on the velocity field is negligible, the dynamo leads to an exponential amplification of the large-scale magnetic field. As the Lorentz force becomes stronger, the field growth slows down and the system gradually settles into a (statistically) steady state: the dynamo action saturates.

The mean-field dynamo equation has the form

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times (\alpha \mathit{B})+\mathrm{\nabla }\times (\mathit{V}\times \mathit{B})+\beta {\mathrm{\nabla }}^2\mathit{B},\end{array}$$(1)

where B is the large-scale magnetic field, α and β are the turbulent transport coefficients representing the mean induction effects of the helical interstellar turbulence (the α-effect) and turbulent magnetic diffusion, respectively, and V is the large-scale velocity field. The latter is dominated by differential rotation but can also include galactic outflows (fountain or wind) and accretion.

Detailed reviews of the galactic dynamo and comprehensive references can be found in Ruzmaikin et al. (1988), Beck et al. (1996), Shukurov et al. (2007), and Shukurov & Subramanian (2018). We present here a very short overview of the theory with the number of specific references reduced to a minimum.

The mean-field galactic dynamo equation has been solved under various approximations that cover a wide range of galactic environments. Our goal is to present a general class of physically-motivated magnetic field models that can be used to fit observations without the need to delve too deeply into the theory. Thus, we present a parametrised model with the large-scale magnetic field in the form of an expansion over appropriate basis functions, specifically the modes of free decay which solve the diffusion equation in the disc and spherical geometries. The form of the magnetic structures obtained is controlled by the choice of the expansion coefficients. The large-scale magnetic field in the model consists of a superposition of approximate solutions of the kinematic mean-field dynamo equations for a thin disc (Sect. 3) and a spherical gaseous halo (Sect. 4). Our use of the dynamo equations is less restrictive than it might seem since their solutions can be used as a functional basis to represent a wide class of complex magnetic configurations, not necessarily produced by dynamo action, in terms of a small number of the expansion coefficients. Unlike representations in terms of the Euler potentials (Ferrière & Terral 2014), magnetic configurations presented here can be helical. Another advantage of our approach is that all variables and parameters of the model are either directly observable or related to observable quantities.

The approximate nature of the solutions that we use is due to the approximations adopted to solve the equations as well as the simple superposition of separate disc and halo solutions of the dynamo equation. Such a superposition is not quite consistent with the presumably non-linear nature of galactic dynamos. However, the non-linear effects do not, plausibly, affect the spatial distribution of the large-scale magnetic field too strongly and the marginally stable dynamo eigenfunctions often provide a satisfactory approximation for the non-linear solutions (Chamandy et al. 2014).

The analytic solutions of the mean-field dynamo equations for the galactic discs and halos presented here are obtained assuming that the disc is thin and the halo is spherical. The thin-disc solutions are applicable at those distances to the galactic centre s where h/s ≲ 0.1, where h is the scale height of the warm ionised gas that is assumed to host the mean magnetic field.

2.1. Symmetries of galactic magnetic fields

It is convenient to introduce cylindrical polar coordinates (s, ϕ, z) with the origin at the galactic centre and the z-axis parallel to the angular velocity Ω. Hence, Ω = (0, 0, Ω) and z = 0 at the galactic mid-plane. However, spherical coordinates (r, θ, ϕ) are natural for the halo, with the polar axis θ = 0 aligned with the z-axis of the cylindrical frame and the mid-plane (equator) at θ = π/2.

Solutions of the dynamo equation are sensitive to the geometry of the dynamo region. In a thin disc, large-scale magnetic fields of even parity strongly dominate, that is, B_s, ϕ(−z) = B_s, ϕ(z) and B_z(−z) = −B_z(z); this configuration corresponds to quadrupolar symmetry. Without dynamo action, quadrupolar magnetic fields in a thin disc decay slower than dipolar ones. As a result, for realistic values of parameters, dipolar magnetic fields can be supported by the dynamo only in the central parts of the discs of spiral galaxies, ≲1 kpc. In a quasi-spherical halo, however, both odd (dipolar) and even (quadrupolar) magnetic fields can be maintained with almost equal ease, with B_r, ϕ(−z)= −B_r, ϕ(z) and B_θ(−z) = B_θ(z) for the dipolar symmetry and B_r, ϕ(−z) = B_r, ϕ(z) and B_θ(−z) = −B_θ(z) in the quadrupolar field (with z = rcosθ). Moreover, magnetic fields in the two halves of the halo, z >  0 and z <  0, can be disconnected by the disc, so that the symmetry of the magnetic field in the halo is only weakly constrained. The model proposed here provides freedom in selecting any symmetry of the magnetic field in the disc and the halo independently.

Despite deviations from axial symmetry, mainly associated with the spiral pattern, galactic discs are sufficiently symmetric in azimuth that the axially-symmetric modes dominate the dynamo. Therefore, deviations from axial symmetry in the large-scale magnetic field, however strong they might be, can be included as distortions of a background axially symmetric magnetic structure. In this paper, we mainly consider axially symmetric galaxies and axially symmetric magnetic fields and discuss extensions to more general configurations in Sect. 5.3.

2.2. Boundary conditions

The simplest and best explored solutions of the mean-field dynamo equations are obtained with the so-called vacuum boundary conditions, that is, under the assumption that the electric current density outside the dynamo region is negligible in comparison with any electric currents within it. Equivalently, the magnetic diffusivity (inversely proportional to electric conductivity) outside the dynamo region is assumed to be much larger than within it. With vanishing electric current density, ∇ × B = 0, the magnetic field is potential.

Neglecting external electric currents in the gaseous halo appears to be reasonable given the low density of intergalactic plasma. Regarding galactic discs surrounded by a gaseous halo, it is important to note that the magnetic diffusivity relevant to a large-scale magnetic field is the turbulent diffusivity. Poezd et al. (1993) argue that the turbulent magnetic diffusivity in galactic haloes is about 50 times larger than in the disc. This justifies the application of vacuum boundary conditions to the large-scale magnetic field at the disc surface as well. In other words, we assume that the extension of the disc’s magnetic field into the halo is a potential magnetic field that adds to the magnetic field produced in situ in the halo.

In a thin disc, the vacuum boundary conditions have the form B_ϕ(±h) = 0 and B_s(±h)≈0, where z = ±h(s) is the disc surface. The boundary condition for B_ϕ is exact whereas the accuracy of the boundary condition for B_s is higher when the disc is thinner (Priklonsky et al. 2000; Willis et al. 2004). The potential magnetic field around the disc, that satisfies these boundary conditions, is purely vertical. The vacuum boundary conditions for the halo are B_ϕ(r_h) = 0 and ∇ × B = 0 at r >  r_h, where r_h is the halo radius.

3. Magnetic field in the disc

3.1. Rotation curve and disc thickness

While the model can be applied to any galaxy, our choice of fiducial parameters is motivated by the Milky Way, with the disc rotation curve of Clemens (1985). To illustrate the impact of the rotation curve on the magnetic field, we also consider a flat rotation curve (i.e., the rotational speed is nearly independent of s at large distances from the disc axis),

$$\begin{array}{c}\hfill V(s)=V_0\frac{1-exp(-s/s_{\ast })}{1-exp(-s_0/s_{\ast })},\end{array}$$(2)

where s₀ is a reference galactocentric distance defined below, V₀ is the rotation speed at s = s₀ and we take s_* = 250 pc. The two rotation curves and the corresponding velocity shear rates, S = sdΩ/ds, are shown in Fig. 1.

The disc scale height is assumed to increase exponentially with the cylindrical radius (a flared disc),

$$\begin{array}{c}\hfill h(s)=h_{0}exp\left(\frac{s-s_0}{s_{\mathrm{h}}}\right),\end{array}$$(3)

where we adopt a flaring length scale of s_h = 5 kpc, similar to that of the MW H I disc, and s₀ is the Galactocentric distance of the Sun (Kalberla & Kerp 2009). In our fiducial model, we adopted the characteristic height h₀ = 0.5 kpc, which is the scale height of the Lockman layer (Lockman 1984; Dickey & Lockman 1990) near the Sun.

The radius of the dynamo-active part of the disc is chosen to be s_d = 17 kpc, similar to the radius of the supernova distribution in the MW (Case & Bhattacharya 1998). The fiducial values of the parameters that appear in the model (which are also input parameters for the GALMAG code) are shown in Table 1. Parameters of the halo are introduced in Sect. 4.2.

3.2. Thin-disc dynamos

In terms of cylindrical coordinates (s, ϕ, z), the radial and azimuthal components of the axisymmetric α²ω-dynamo Eq. (1) can be written as

$$\begin{array}{cc}\hfill \frac{\partial B_s}{\partial t}& =-\frac{\partial (\alpha B_{\varphi })}{\partial z}+\beta \frac{{\partial }^2{B}_s}{\partial z^2}+\beta \frac{\partial }{\partial s}\left[\frac{1}{s}\frac{\partial (s{B}_s)}{\partial s}\right],\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \frac{\partial B_{\varphi }}{\partial t}& =S{B}_s+\frac{\partial (\alpha B_s)}{\partial z}+\beta \frac{{\partial }^2{B}_{\varphi }}{\partial z^2}+\beta \frac{\partial }{\partial s}\left[\frac{1}{s}\frac{\partial (s{B}_{\varphi })}{\partial s}\right],\hfill \end{array}$$(5)

where S = s ∂Ω/∂s is the velocity shear due to differential rotation. We do not exhibit the equation for B_z since, in a thin disc, it decouples from the equations shown and can be solved separately (equivalently, B_z can be derived from ∇ ⋅ B = 0).

It is convenient to use dimensionless variables, denoted with tilde,

$$\begin{array}{c}\hfill \stackrel{\sim }{s}=s/s_0\text{and}\stackrel{\sim }{z}=z/h_0,\text{with}{h}_0=h(s_0),\end{array}$$(6)

where s₀ is the reference cylindrical radius within the disc – for example, s₀ = s_⊙ ≈ 8.5 kpc is a convenient choice in the MW. Using different length units across and along the disc allows us to make the disc thinness explicit and quantified with the (small) aspect ratio

$$\begin{array}{c}\hfill ϵ=h_0/s_0.\end{array}$$(7)

The large-scale velocity is measured in the units of a characteristic rotational speed V₀ = V(s₀),

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathit{V}}=\frac{\mathit{V}}{V_0},\stackrel{\sim }{\mathrm{\Omega }}=\mathrm{\Omega }\frac{s_0}{V_0}·\end{array}$$(8)

Velocity shear due to differential rotation is non-dimensionalised similarly,

$$\begin{array}{c}\hfill \stackrel{\sim }{S}=\frac{S}{S_0},\text{with}{S}_0=S(s_0).\end{array}$$(9)

The unit of time is the turbulent magnetic diffusion time across the disc. With β_d the turbulent magnetic diffusivity in the disc, the dimensionless time is

$$\begin{array}{cc}\hfill \stackrel{\sim }{t}& =t\beta {}_{\mathrm{d}}/h_0^2.\hfill \end{array}$$(10)

The magnitude of the α-effect can be estimated as

$$\begin{array}{c}\hfill \alpha \simeq min(l^2\mathrm{\Omega }/h,v),\end{array}$$(11)

where l and v are the turbulent scale and speed, and the corresponding fiducial value is used to non-dimensionalise α:

$$\begin{array}{c}\hfill \stackrel{\sim }{\alpha }=\frac{\alpha }{{\alpha }_0},{\alpha }_0=\frac{l^2{V}_0}{h_0{s}_0}·\end{array}$$(12)

The magnitude of α cannot exceed α = v because it is a measure of the helical part of the turbulent flow speed, hence α/v cannot exceed unity. This limit is usually important only in the central parts of galaxies (where the thin-disc approximation does not apply anyway). Equation (11) gives the magnitude of α and its dependence on s through the variations of Ω and h with s. It is expected that α is an odd function of z. Gressel et al. (2008a,b) and Bendre et al. (2015) confirm this and discuss the dependence of α on z in numerical simulations of the supernova-driven interstellar medium. We adopt a factorised form for α(s, z),

$$\begin{array}{c}\hfill \stackrel{\sim }{\alpha }(s,z)=\frac{\stackrel{\sim }{\mathrm{\Omega }}(s)}{\stackrel{\sim }{h}(s)}a(z),\text{where}a(-z)=-a(z),\end{array}$$(13)

where we assume that l, in Eqs. (11) and (12), is independent of s. The model can be generalised straightforwardly to more general forms of α(s, z).

When galactic outflow is neglected, the dynamo is fully characterised by two dimensionless control parameters that quantify the intensity of the mean magnetic induction due to helical turbulence and differential rotation, respectively:

$$\begin{array}{c}\hfill R_{\alpha \text{d}}=\frac{h_0{\alpha }_0}{\beta {}_{\mathrm{d}}}=\frac{l^2{V}_0}{s_0\beta {}_{\mathrm{d}}},R_{\omega \text{d}}=\frac{h_0^2{S}_0}{\beta {}_{\mathrm{d}}},\end{array}$$(14)

where the subscript “d” refers to the disc (similar parameters are defined slightly differently in the halo – see Sect. 4).

In terms of dimensionless variables, Eqs. (4) and (5) reduce to

$$\begin{array}{cc}\hfill \frac{\partial B_s}{\partial \stackrel{\sim }{t}}& =-R_{\alpha \text{d}}\frac{\partial (\stackrel{\sim }{\alpha }{B}_{\varphi })}{\partial \stackrel{\sim }{z}}+\frac{{\partial }^2{B}_s}{\partial {\stackrel{\sim }{z}}^2}+{ϵ}^2\frac{\partial }{\partial \stackrel{\sim }{s}}\left[\frac{1}{\stackrel{\sim }{s}}\frac{\partial (\stackrel{\sim }{s}{B}_s)}{\partial \stackrel{\sim }{s}}\right],\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill \frac{\partial B_{\varphi }}{\partial \stackrel{\sim }{t}}& =R_{\omega \text{d}}\stackrel{\sim }{S}{B}_s+R_{\alpha \text{d}}\frac{\partial (\stackrel{\sim }{\alpha }{B}_s)}{\partial \stackrel{\sim }{z}}+\frac{{\partial }^2{B}_{\varphi }}{\partial {\stackrel{\sim }{z}}^2}+{ϵ}^2\frac{\partial }{\partial \stackrel{\sim }{s}}\left[\frac{1}{\stackrel{\sim }{s}}\frac{\partial (\stackrel{\sim }{s}{B}_{\varphi })}{\partial \stackrel{\sim }{s}}\right]·\hfill \end{array}$$(16)

It is now clear that the solutions are fully determined by R_αd, R_ωd, $\stackrel{\sim }{S}(s)$, $\stackrel{\sim }{\alpha }(s)$ and ϵ.

In a thin disc, the magnetic field distribution along z is established over a time scale h²/β_d ≃ 5 × 10⁸ yr which is ϵ⁻² = (s₀/h₀)² times shorter than the time scale at which the radial distribution evolves, $s_0^2/{\beta }_{\text{d}}$. Because of this difference, the radial derivatives in Eqs. (15) and (16) have ϵ² as a factor. Therefore, the magnetic field distribution in a thin disc can be represented as a local solution (at a given galactocentric distance s), b(z; s), modulated by an envelope Q(s). The local solution also depends on s since α, S and h vary with s, but this variation is parametric. It is convenient to normalise the local solution to unit surface magnetic energy density, ${\displaystyle {\int }_{-h}^h|}\mathit{b}{|}^2\text{d}z=1$, at all values of s: then Q(s) represents magnetic field strength at the galactocentric radius s (an envelope of the local solutions). Thus, asymptotic solutions of Eqs. (15) and (16) for ϵ ≪ 1 have the form

$$\begin{array}{c}\hfill \mathit{B}(s,z,t)=exp(\mathrm{\Gamma }t)Q(s)\mathit{b}(z;s).\end{array}$$(17)

The magnetic field varies exponentially with time at a rate Γ in the kinematic dynamo stage. In a saturated thin-disc dynamo, the solution has a similar form but with Γ = 0 (Poezd et al. 1993). The local solution is discussed in Sect. 3.3, whereas Sect. 3.4 presents the radial solution Q(s).

To simplify the notation, we use exclusively the dimensionless variables with the tilde suppressed in the remaining part of this section unless otherwise stated.

3.3. Local solutions

Governing equations for the local solution, b = exp(γt) (b_s, b_z), follow from Eqs. (15) and (16) when we put ϵ = 0:

$$\begin{array}{cc}\hfill \gamma (s)b_s& =-R_{\alpha \text{d}}\frac{\partial }{\partial z}\left[\alpha (s,z)b_{\varphi }\right]+\frac{{\partial }^2{b}_s}{\partial z^2},\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill \gamma (s)b_{\varphi }& =R_{\omega \text{d}}S(s)b_s+\frac{{\partial }^2{b}_{\varphi }}{\partial z^2}+R_{\alpha \text{d}}\frac{\partial }{\partial z}\left[\alpha (s,z)b_s\right].\hfill \end{array}$$(19)

To allow for the disc flaring, we introduce the following new variables:

$$\begin{array}{c}\hfill \widehat{z}=\frac{z}{h(s)},{\widehat{b}}_s=\frac{b_s}{R_{\alpha \text{d}}\mathrm{\Omega }(s)}\text{and}\widehat{\gamma }=\gamma (s)h^2(s).\end{array}$$(20)

Since R_αd ≪ R_ωd at s ≳ 1 kpc in most spiral galaxies, we can omit the term proportional to R_αd in Eq. (19), thus obtaining the αω-dynamo approximation:

$$\begin{array}{cc}\hfill \widehat{\gamma }(s){\widehat{b}}_s& =-\frac{\partial }{\partial \widehat{z}}\left[a(z)b_{\varphi }\right]+\frac{{\partial }^2{\widehat{b}}_s}{\partial {\widehat{z}}^2},\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill \widehat{\gamma }(s)b_{\varphi }& =\mathcal{D}(s){\widehat{b}}_s+\frac{{\partial }^2{b}_{\varphi }}{\partial {\widehat{z}}^2},\hfill \end{array}$$(22)

where a(z), the z-dependent part of α(s, z), is defined in Eq. (13) and 𝒟(s), the local dynamo number, includes the radial variation of the dynamo parameters:

$$\begin{array}{c}\hfill \mathcal{D}(s)=R_{\alpha \text{d}}{R}_{\omega \text{d}}\mathrm{\Omega }(s)S(s)h^2(s).\end{array}$$(23)

It is also convenient to introduce the local (in galactocentric radius) values of R_αd and R_ωd,

$$\begin{array}{c}\hfill {\mathcal{R}}_{\alpha }(s)=R_{\alpha \text{d}}\mathrm{\Omega }(s),{\mathcal{R}}_{\omega }(s)=R_{\omega \text{d}}S(s)h^2(s),\end{array}$$(24)

so that 𝒟 = ℛ_αℛ_ω. The normalisation condition for the local solution reduces to

$$\begin{array}{c}\hfill {\displaystyle {\int }_{-1}^1(|{\widehat{b}}_s{|}^2+{|b_{\varphi }|}^2)}\mathrm{d}\widehat{z}=1,\end{array}$$(25)

where we have neglected ${\widehat{b}}_z$ because this component of magnetic field is, on average, weaker than the other two. To this order in ϵ, b_z does not enter the equations for b_s and b_ϕ and can be solved for separately. The result can be shown to be identical to that obtained from ∇ ⋅ B = 0.

The vacuum boundary conditions are

$$\begin{array}{c}\hfill {\widehat{b}}_s=b_{\varphi }=0\text{at}\widehat{z}=\pm 1.\end{array}$$(26)

Since a(z) is an odd function of z, solutions of Eqs. (21) and (22) split into two independent classes, even and odd in z (or quadrupolar and dipolar, respectively). These can be distinguished using the symmetry conditions at the galactic mid-plane,

$$\begin{array}{cc}\hfill \frac{\partial {\widehat{b}}_s}{\partial \widehat{z}}& =\frac{\partial b_{\varphi }}{\partial \widehat{z}}=b_z=0\text{at}z=0\text{(even parity)},\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill {\widehat{b}}_s& =b_{\varphi }=\frac{\partial b_z}{\partial \widehat{z}}=0\text{at}z=0\text{(odd parity)}.\hfill \end{array}$$(28)

Approximate solutions for both parity families can be obtained in the form of an expansion over the free-decay modes obtained as solutions of Eqs. (21) and (22) for a(z) = 𝒟(s)=0. The procedure is discussed in detail by Shukurov et al. (2008) and Chamandy et al. (2014), and here we only provide the results.

From Eqs. (21)–(27), we obtain the following approximate solutions of quadrupolar symmetry, for $a(\widehat{z})=sin(\pi \widehat{z}/2)$:

$$\begin{array}{cc}\hfill b_s(z;s)& \approx {\mathcal{R}}_{\alpha }(s)K_0(s)[cos\frac{\pi z}{2h(s)}+\frac{3\sqrt{-\mathcal{D}(s)}}{4{\pi }^{3/2}}cos\frac{3\pi z}{2h(s)}],\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill b_{\varphi }(z;s)& \approx -2K_0(s)\sqrt{-\frac{\mathcal{D}(s)}{\pi }}cos\frac{\pi z}{2h(s)},\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill \gamma (s)& \approx \frac{1}{h^2(s)}[-\frac{{\pi }^2}{4}+\frac{1}{2}\sqrt{-\pi \mathcal{D}(s)}],\hfill \end{array}$$(31)

where

$$\begin{array}{c}\hfill K_0(s)={[1-\frac{4\mathcal{D}(s)}{\pi }-\frac{9\mathcal{D}(s)}{16{\pi }^3}]}^{-1/2}\end{array}$$(32)

is the normalization factor obtained using Eq. (25). These solutions are shown in Fig. 2. Equation (31) provides the critical value of the local dynamo number required for the local amplification and maintenance of the magnetic field: γ ≥ 0 for 𝒟 ≤ 𝒟_c ≈ −π³/4 ≈ −8. Other choices of the functional form of $a(\widehat{z})$ lead to slightly different values of the critical dynamo number (e.g. −11 for a = z), but the difference hardly has any practical consequences.

In the odd-parity solutions, the term proportional to (−𝒟)^1/2 in b_s vanishes for $a(\widehat{z})=sin(\pi \widehat{z}/2)$. Therefore, it is more convenient to use a similar solution with $a(z)=\widehat{z}$ that satisfies the symmetry condition Eq. (28), here written to the lowest order in (−𝒟)^1/2:

$$\begin{array}{cc}\hfill b_s(z;s)& \approx K_1(s){\mathcal{R}}_{\alpha }(s)\sqrt{2}sin\frac{\pi z}{h(s)},\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill b_{\varphi }(z;s)& \approx -2K_1(s)\sqrt{-\mathcal{D}(s)}sin\frac{\pi z}{h(s)},\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill \gamma (s)& \approx \frac{1}{h^2(s)},\hfill \end{array}$$(35)

with

$$\begin{array}{c}\hfill K_1(s)={[1-4\mathcal{D}(s)]}^{-1/2}.\end{array}$$(36)

The dipolar modes can be excited, γ ≥ 0, for 𝒟 ≤ −2π⁴ ≈ −195, a threshold much higher than for the quadrupolar modes. This is true for any plausible form of $a(\widehat{z})$ and explains the predominance of quadrupolar magnetic fields in thin discs.

The local solutions are derived, formally, for |𝒟|≪1 but they remain reasonably accurate for |𝒟| as large as about 50 or more (Ji et al. 2014). We compare in Fig. 2 the local quadrupolar eigenfunctions obtained from numerical solution of Eqs. (21) and (22) with the approximate solutions Eqs. (29)–(31) for 𝒟 = −20 and −50, values typical of the main parts of spiral galaxies (see Fig. 3).

The local eigenfunctions presented above are approximate solutions of the mean-field dynamo equations. However, the functional basis of these solutions, the free-decay modes, can be used as a complete functional basis to represent any magnetic field configuration. The free-decay modes are solutions of Eqs. (21) and (22) with a(z) = 0 and 𝒟(s) = 0, so that the equations decouple and can easily be solved. Normalised as in Eq. (25), the dipolar and quadrupolar free-decay eigenfunctions and eigenvalues (identified with superscripts d and q, respectively) have the following respective forms with n = 1, 2, 3, …:

$$\begin{array}{c}\hfill {\mathit{b}}_n^{(\text{d})}=\left(\begin{array}{c}sin(\pi n\widehat{z})\\ 0\end{array}\right),{\mathit{b}}_n^{(\text{d})\prime }=\left(\begin{array}{c}0\\ sin(\pi n\widehat{z})\end{array}\right),{\widehat{\gamma }}_n^{(\text{d})}=-{\pi }^2{n}^2,\end{array}$$(37)

$${\mathit{b}}_n^{(\text{q})}=\left(\begin{array}{c}\cos \left[\pi \left(n-\frac{1}{2}\right)\widehat{z}\right]\\ 0\end{array}\right),{\mathit{b}}_n^{(\text{q})\prime }=\left(\begin{array}{c}0\\ \cos \left[\pi \left(n-\frac{1}{2}\right)\widehat{z}\right]\end{array}\right),$$(38)

$${\widehat{\gamma }}_n^{(\text{q})}=-{\pi }^2{\left(n-\frac{1}{2}\right)}^2$$(39)

The free-decay disc modes are double degenerate as two orthogonal modes of the same parity (distinguished by prime) correspond to each eigenvalue.

3.4. Radial solution

When Eq. (17) is substituted into Eqs. (15) and (16), and Eqs. (18) and (19) are allowed for, equations for both B_s and B_ϕ reduce to the same equation for Q(s), that is,

$$\begin{array}{c}\hfill {ϵ}^2\frac{\partial }{\partial s}\left[\frac{1}{s}\frac{\partial }{\partial s}(sQ)\right]+[\gamma (s)-\mathrm{\Gamma }]Q=0,\end{array}$$(40)

where γ(s) is the local growth rate obtained as a part of the local solution, Eq. (31) or (35). This equation applies to both even and odd local solutions, and to both α²ω and αω-dynamos, that is, with and without the term proportional to R_αd in Eq. (19). The boundary conditions adopted are

$$\begin{array}{c}\hfill Q(0)=Q(s{}_{\mathrm{d}})=0.\end{array}$$(41)

The condition at s = 0 follows from axial symmetry, whereas that at the outer boundary of the dynamo-active region s = s_d is adopted for the sake of simplicity.

When Ω(s) and h(s) are known, the global growth rate Γ and the radial distribution of the magnetic field strength Q(s) can be obtained by solving Eq. (40) numerically. However, for our present purposes it is more useful to obtain an approximate analytical solution for Q(s). This allows faster computations at the cost of an additional approximation. When the magnetic field model is used in Bayesian analyses, the computation speed is of primary importance.

To obtain a simple analytical solution of Eq. (40), consider γ to be a constant, γ = γ₀, within the disc radius, s <  s_d, and zero outside,

$$\begin{array}{c}\hfill \gamma (s)=\{\begin{array}{cc}{\gamma }_0,\hfill & 0\le s<s{}_{\mathrm{d}},\hfill \\ 0,\hfill & s\ge s{}_{\mathrm{d}}.\hfill \end{array}\end{array}$$(42)

A suitably averaged value of γ(s) within the disc can be adopted for γ₀. The relevance of this approximation depends on the specific case, particularly the rotation curve and the rate of the disc flaring, as shown in Fig. 3. The inner parts of the disc, s ≲ 1 kpc, should be disregarded since the thin-disc approximation is not applicable there. At s ≳ 1 kpc, the approximation γ(s) = const does not appear unreasonable, especially for the Milky Way rotation curve, as shown in the lower panel of Fig. 3.

For γ(s) of Eq. (42), Eq. (40) with boundary conditions Eq. (41) can be solved to yield the eigensolutions

$$\begin{array}{cc}\hfill Q_n(s)& =J_1(k_{n}s/s{}_{\mathrm{d}}),n=1,2,3,\dots ,\hfill \end{array}$$(43)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_n& ={\gamma }_0-{ϵ}^2{k}_n^2,\hfill \end{array}$$(44)

where J₁(x) is the standard Bessel function (of order one) and k_n ≈  3.83, 7.02, 10.17, … are its zeros, thus, J₁(k_n) = 0. The first four modes of Q_n(s) are shown in Fig. 4. Independently of the form of γ(s), the lowest radial mode, Q₁(s), is sign-constant but Q_n(s) has n − 1 zeros, and thus the scale of variation of Q_n(s) decreases with n. This feature of the solution is responsible for the reversals of the large-scale magnetic field discussed in Sect. 3.9.

The evolving radial distribution of magnetic field is obtained as a superposition of the eigensolutions of Eq. (40):

$$\begin{array}{c}\hfill Q(s,t)={\displaystyle \sum _{n=1}^{\infty }{C}_n{\mathrm{e}}^{{\mathrm{\Gamma }}_{n}t}{Q}_n(s)},\end{array}$$(45)

where, in the context of dynamo models, the coefficients C_n are determined by the initial conditions. Alternatively, the eigenfunctions Q_n(s)b(z; s) can be used as the basis functions to represent a given magnetic field with exp(Γ_nt) absorbed into C_n. The set of radial eigenfunctions Q_n is complete, so that any radial distribution of magnetic field can be represented in this form. The only constraint on the form of the solution is due to the fact that the set of the local solutions b is incomplete. Nevertheless, a wide class of magnetic field distributions along z can be represented as a superposition of various quadrupolar and dipolar local modes, so that the lack of the functional completeness of the local solutions is not likely to be restrictive in practice. Otherwise, the complete set of local free-decay modes ${\mathit{b}}_m^{\text{(q,d})}(z;s)$ of Eqs. (37)–(39) can be used instead of the local dynamo solutions to construct a more general expansion over a complete set of basis functions of the form

$$\begin{array}{c}\hfill \mathit{B}(s,z)={\displaystyle \sum _{m=1}^{\infty }}{\displaystyle \sum _{n=1}^{\infty }{C}_{\mathit{mn}}{Q}_n(s){\mathit{b}}_m^{\text{(q,d})}(z;s)},\end{array}$$

which can be used to parametrise (in terms of the coefficients C_mn) an arbitrary magnetic field that does not need to be a solution of the dynamo equations.

3.5. Vertical magnetic field

The axially symmetric vertical component of the magnetic field, B_z, can be obtained from

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathit{B}=\frac{1}{s}\frac{\partial }{\partial s}(s{B}_s)+\frac{\partial B_z}{\partial z}=0,\end{array}$$(46)

using B_s(s, z) = b_s(z; s)Q(s) from Eq. (29) or (33), and Eq. (43). We have assumed that γ₀ = const to derive the simple form Eq. (43), implying 𝒟 = const. It is therefore justifiable to neglect the dependence of 𝒟, Ω and h on s when differentiating B_s in Eq. (46) and only retain the dependence of Q on s (this simplification can easily be relaxed if required). By virtue of linearity, Eq. (46) can be solved for each radial mode n separately:

$$\begin{array}{cc}\hfill -\frac{\partial B_z^{(n)}}{\partial z}& =\frac{1}{s}\frac{\partial }{\partial s}(s{B}_s^{(n)})={\mathcal{R}}_{\alpha }{K}_0{C}_n\frac{1}{s}\frac{\mathrm{d}}{\mathrm{d}s}\left[s{J}_1(k_{n}s/s_{}\text{d})\right]\hfill \\ \hfill & \times (cos\frac{\pi z}{2h(s)}+\frac{3}{4{\pi }^{3/2}}\sqrt{-\mathcal{D}(s)}cos\frac{3\pi z}{2h(s)})·\hfill \end{array}$$(47)

Then, for the quadrupolar parity,

$$\begin{array}{cc}\hfill B_z^{(n)}& =-{\mathcal{R}}_{\alpha }{K}_0{C}_n\frac{k_n}{s_{}\text{d}}{J}_0\left(\frac{k_{n}s}{s_{}\text{d}}\right)\hfill \\ \hfill & \times {\displaystyle {\int }_0^z}(cos\frac{\pi z{}^{\prime }}{2h}+\frac{3\sqrt{-\mathcal{D}}}{4{\pi }^{3/2}}cos\frac{3\pi z{}^{\prime }}{2h})\mathrm{d}z{}^{\prime }\hfill \end{array}$$(48)

where the dependence of ℛ_α 𝒟, h and K₀ on s should be allowed for.

3.6. Magnetic field in the disc

We can now collect the solutions from Sects. 3.3–3.5 to write the approximate solution of the dynamo equation (evolving or at a fixed time) for a thin disc as a sum of N_r radial eigenmodes,

$$\begin{array}{c}\hfill \left(\begin{array}{c}{B}_s\\ B_{\varphi }\\ B_z\end{array}\right)={\displaystyle \sum _{n=1}^{N_{\mathrm{r}}}}{C}_n\left(\begin{array}{c}{B}_s^{(n)}\\ B_{\varphi }^{(n)}\\ B_z^{(n)}\end{array}\right),\end{array}$$(49)

where, for the quadrupolar symmetry,

$$\begin{array}{cc}\hfill B_s^{(n)}& =K_0(s){\mathcal{R}}_{\alpha }(s)J_1(k_{n}s/s_{}\text{d})\hfill \\ \hfill & \times [cos\frac{\pi z}{2h(s)}+\frac{3\sqrt{-\mathcal{D}(s)}}{4{\pi }^{3/2}}cos\frac{3\pi z}{2h(s)}],\hfill \end{array}$$(50)

$$\begin{array}{cc}\hfill B_{\varphi }^{(n)}& =-2\sqrt{-\frac{\mathcal{D}(s)}{\pi }}{K}_0(s)J_1(k_{n}s/s_{}\text{d})cos\frac{\pi z}{2h(s)},\hfill \end{array}$$(51)

$$\begin{array}{cc}\hfill B_z^{(n)}& =-\frac{2k_{n}h(s)}{\pi s_{}\text{d}}{K}_0(s){\mathcal{R}}_{\alpha }(s)J_0(k_{n}s/s_{}\text{d})\hfill \\ \hfill & \times [sin\frac{\pi z}{2h(s)}+\frac{\sqrt{-\mathcal{D}(s)}}{4{\pi }^{3/2}}sin\frac{3\pi z}{2h(s)}],\hfill \end{array}$$(52)

and similarly for the dipolar symmetry. The local dynamo parameters 𝒟(s) and ℛ_α(s) are defined in Eqs. (23) and (24), whereas an example of the form of h(s) is given by Eq. (3). In what follows (and in the GALMAG code), we normalise each eigenmode so that each expansion coefficient represents the strength of the corresponding part of the magnetic field at the reference radius s₀, that is, |C_n|=|B⁽ⁿ⁾(s₀, 0)|.

The coefficients C_n can be chosen to fix the strength and to reproduce any radial distribution of the magnetic field. For instance, when used to approximate the magnetic field observed in the disc, C_n are used to fit the observations to a desired accuracy. When used to simulate a growing magnetic field, its evolution is introduced through $C_n=C_n^{(0)}{\text{e}}^{{\Gamma }_{n^t}}$, where the initial values $C_n^{(0)}$ at t = 0 are obtained from the similar expansion for the seed magnetic field. Galactic evolution can be included by an appropriate time variation of 𝒟(s), Ω(s), h(s) and s_d.

3.7. Dynamo parameters

In order to use the solution presented above, the dimensionless dynamo control parameters R_αd and R_ωd have to be specified. We adopt the reference radius s₀ = s_⊙ ≈ 8.5 kpc and, correspondingly, V₀ = 220 km s⁻¹ and S₀ = −35 km s⁻¹ kpc⁻¹ as obtained from the rotation curve. An estimate of the turbulent magnetic diffusivity widely used in various applications derives from the mixing length theory,

$$\begin{array}{c}\hfill \beta =\frac{1}{3}lv,\end{array}$$(53)

where l and v are the turbulent scale and speed, respectively. With l ≃ 50 pc (see Hollins et al. 2017, and references therein) and v ≃ 10 km s⁻¹ (e.g. Mac Low & Klessen 2004), we have β ≃ 5 × 10²⁵ cm² s⁻¹. Equations (53) and (14) then yield

$$\begin{array}{c}\hfill R_{\alpha \text{d}}=0.39,R_{\omega \text{d}}=-53.\end{array}$$(54)

A summary of parameters used by GALMAG and their fiducial values can be found in Table 1.

3.8. The role of the disc flaring and details of the rotation curve

Figure 5 shows the vertical cross section of magnetic field in the disc, constructed using the first two radial eigenmodes with (C₁, C₂) = (4.6 μG, −1.5 μG). These distributions have an absolute maximum of |B| at large s. This happens because the local dynamo number of Eq. (23) increases with galactocentric distance in an exponentially flared disc, 𝒟 ∝ ΩSh² ∝ s⁻² exp(2s/s_h) for a flat rotation curve, Ω ∝ S ∝ s⁻¹. If h indeed increases with s faster than s², the outer parts of galactic discs can have relatively strong magnetic fields at early stages of magnetic field growth when these kinematic solutions apply. Such distributions may occur in young and evolving galaxies. Non-linear dynamo effects eventually limit the local magnetic field strength to a value related to equipartition between magnetic and turbulent kinetic energies, B² ≈ 4πρv²(𝒟/𝒟_cr − 1)∝s⁻² exp(2s/s_h − s/s_ρ) assuming that ρ ∝ exp(−s/s_ρ). The radial profile of magnetic field strength then depends on the relation between the radial length scales of the gas density and disc thickness. The number density of H I in the MW has s_ρ ≃ 3 kpc (Kalberla & Kerp 2009) which is close to s_h/2 ≃ 2.5 kpc. Therefore, we cannot exclude the possibility that the magnetic field remains strong in the outer MW. The effective boundary of the dynamo active region is then determined by the rapid increase of the local dynamo time scale γ⁻¹(s)≃h²(s)/β_d with s in a flared disc: this time scale exceeds 10¹⁰ yr where h ≳ 1 kpc for β_d = 5 × 10²⁵ cm² s⁻¹, and the growth of magnetic field becomes practically negligible. For h(s) given by Eq. (3), this happens at s ≳ 11 kpc provided β_d is independent of s. When the magnetic field is stronger in the outer parts of the disc, the boundary condition Q(s_d) = 0 may be too restrictive. We have considered solutions with ∂(sQ)/∂s|_s = sd = 0 to confirm that the magnetic field distribution in the main part of the disc is not significantly affected but the outer field maximum becomes more pronounced.

The increase of the local dynamo number with distance from the galactic centre enhances the large-scale magnetic field in the outer parts of a galactic disc in either the kinematic or saturated dynamo. As a result, the magnetic field energy density may decrease with radius slower than other energy densities in the interstellar medium, as suggested by Beck (2007).

3.9. Field reversals along the galactocentric radius

Magnetic field reversals can be reproduced naturally in the model because the radial eigenfunction Q_n(s) has n − 1 zeros. With an appropriate selection of the expansion coefficients C_n, any desired number of reversals located at any prescribed positions can be produced. While the strength of the magnetic field is controlled by the magnitudes of the expansion coefficients C_n, the number of reversals and their positions along the radius are controlled by the ratios of the coefficients, for instance, C_n/C₁. Since Q₂(s) has one zero while Q₃(s) has two zeros, retaining only the two leading terms in the expansion Eq. (49) allows us to obtain a magnetic field with one radial reversal, whereas in order to have two reversals, Q₃(s) needs to be included. To help selecting the coefficients as required to obtain a magnetic field that has a given number of reversals at desired positions, we present Fig. 6 where the ranges of C₂/C₁ and C₃/C₁ that produce one or two reversals can be read off the left-hand panel. The desired positions of the reversals can be converted into the coefficient ratios using the middle and right-hand panels. Similar diagrams can be constructed for any number of reversals if required.

Figure 7 illustrates the structure of axisymmetric magnetic fields that have radial reversals. Model A, with a reversal at s = 7 kpc, represents the same field as in the bottom right panel of Fig. 5. Model B, constructed using the three leading eigenfunctions, has reversals at s = 7 and 12 kpc. In both models, the magnetic field is normalised so that the azimuthal magnetic field, shown in Fig. 7c, has the strength $B_{\varphi }^{(\mathrm{d})}{|}_{s=s_0,z=0}=-3\mu \mathrm{G}$ in the mid-plane z = 0 at the reference radius s = s₀ = 8.5 kpc.

Figure 7d shows the vertical cross-section (a meridional plane) of the magnetic structure in Model B to demonstrate that, in a thin disc, the horizontal magnetic field components dominate over the vertical field, |B_ϕ| > |B_s|≫|B_z|, only on average. Locally, and especially near the reversals and the disc axis, the vertical magnetic field dominates. This is a direct consequence of the solenoidality of magnetic field: if B_ϕ and B_s are weak, B_z needs to be stronger to ensure that ∇ ⋅ B = 0. The dominance of B_ϕ over B_s is less general in origin: this is a consequence of the stretching of the radial magnetic field by differential rotation, and the larger is the velocity shear the larger is the ratio |B_ϕ/B_s| and the smaller is the magnitude of the magnetic pitch angle p = arctan(B_s/B_ϕ).

4. Magnetic field in the halo

Magnetic field in the spherical halo is obtained as the perturbation solution of the mean-field dynamo equation with free-decay eigenfunctions as the unperturbed solutions. This approach is similar to that employed to obtain the local disc solution in Sect. 3.3.

In the spherical halo, it is convenient to use spherical coordinates (r, θ, ϕ). As for the disc, we define convenient dimensionless variables distinguished by the tilde: spherical radius and time are measured in the units of the halo radius r_h and the corresponding magnetic diffusion time, respectively,

$$\begin{array}{c}\hfill \stackrel{\sim }{r}=r/r_{\mathrm{h}}\text{and}\stackrel{\sim }{t}=t{\beta }_{\mathrm{h}}/r_{\mathrm{h}}^2,\end{array}$$(55)

with β_h the turbulent magnetic diffusivity in the halo. The velocity field and the α-coefficient are normalised as

$$\begin{array}{c}\hfill \stackrel{\sim }{\alpha }=\alpha /{\alpha }_{\mathrm{h}},\stackrel{\sim }{V}=V/V_{\mathrm{h}},\end{array}$$(56)

where α_h is the α-coefficient at the north pole, (r, θ) = (r_h, 0), and V_h is the equatorial rotation velocity at the boundary (r, θ) = (r_h, π/2).

In terms of the dimensionless variables, the mean-field dynamo Eq. (1) reduces to

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial \stackrel{\sim }{t}}=R_{\alpha \text{h}}\stackrel{\sim }{\mathrm{\nabla }}\times (\stackrel{\sim }{\alpha }\mathit{B})+R_{\omega \text{h}}\stackrel{\sim }{\mathrm{\nabla }}\times (\stackrel{\sim }{\mathit{V}}\times \mathit{B})+{\stackrel{\sim }{\mathrm{\nabla }}}^2\mathit{B},\end{array}$$(57)

where we defined, analogously to Eq. (14), the dynamo parameters

$$\begin{array}{c}\hfill R_{\alpha \text{h}}=r_{\mathrm{h}}{\alpha }_{\mathrm{h}}/{\beta }_{\mathrm{h}},R_{\omega \text{h}}=-r_{\mathrm{h}}{V}_{\mathrm{h}}/{\beta }_{\mathrm{h}}.\end{array}$$(58)

To avoid excessively heavy notation, we suppress the tilde on the dimensionless variables and work exclusively with dimensionless variables unless otherwise stated.

Solutions of Eq. (57), growing or decaying at a rate Γ, are sought in the form of an expansion

$$\begin{array}{c}\hfill \mathit{B}=exp(\mathrm{\Gamma }t){\displaystyle \sum _{i=1}^N}{a}_i{\mathit{B}}_i(\mathit{r}),\end{array}$$(59)

in the free-decay modes B_i which are obtained as solutions of Eq. (57) with R_αh = R_ωh = 0,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^2{\mathit{B}}_i={\gamma }_i{\mathit{B}}_i.\end{array}$$(60)

where γ_i <  0 is the rate of exponential decay of the mode B_i. Outside the halo, an electromagnetic vacuum is assumed, implying a potential magnetic field, ∇ × B_i = 0. The boundary conditions that ensure a continuous matching, at the halo boundary r = 1, of the interior magnetic field to a potential exterior magnetic field that decays at infinity as the point dipole (the lowest magnetic multipole) are given by (Moffatt 1978)

$$\begin{array}{c}\hfill \left[{\mathit{B}}_i\right]=0\text{at}r=1,{\mathit{B}}_i=\mathcal{O}\left(r^{-3}\right)\text{for}r\to \infty ,\end{array}$$(61)

where the square brackets denote the jump of the corresponding quantity.

The spatial form and decay rates of the spherical modes of free decay are derived in Appendix A; here we briefly discuss their properties. The free decay modes form a complete, orthonormal set of basis functions (related to spherical harmonics), each either purely poloidal (comprising the field components B_r and B_θ) or purely toroidal (consisting of B_ϕ alone). They can be divided into two classes based on their symmetry about the equator θ = π/2: the symmetric modes are quadrupolar (indicated with superscript “q”) whereas the anti-symmetric modes have a dipolar symmetry (superscript “d”). Their analytic forms can be found in Appendices A.1 and A.2, respectively. Figure 8 shows the structure of the four free-decay modes of each symmetry that have the largest γ_i.

4.1. The perturbation solution

Equation (57) can be conveniently written as

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\widehat{\mathcal{W}}\mathit{B}+{\mathrm{\nabla }}^2\mathit{B},\end{array}$$(62)

where the perturbation operator $\widehat{\mathcal{W}}$ corresponding to the α²ω-dynamo is given by

$$\begin{array}{c}\hfill \widehat{\mathcal{W}}\mathit{B}=R_{\alpha \text{h}}\mathrm{\nabla }\times (\alpha \mathit{B})+R_{\omega \text{h}}\mathrm{\nabla }\times (\mathit{V}\times \mathit{B}).\end{array}$$(63)

As discussed in Appendix B, GALMAG has also an option to use the αω-dynamo operator but this approximation may be questionable in the case of the halo.

We substitute Eq. (59) into Eq. (62), take the scalar product of the result with B_i and integrate over the whole space. As a result, we obtain a homogeneous system of algebraic equations for the expansion coefficients a_i of the form

$$\begin{array}{c}\hfill a_j({\gamma }_j-\mathrm{\Gamma })+{\displaystyle \sum _{i=1}^N}{a}_i{W}_{\mathit{ij}}=0,j=1,2,\dots ,N,\end{array}$$(64)

where

$$\begin{array}{c}\hfill W_{\mathit{ij}}={\displaystyle {\int }_V{\mathit{B}}_i·\widehat{\mathcal{W}}{\mathit{B}}_j}{\mathrm{d}}^3\mathit{r}\end{array}$$(65)

are the matrix elements of the perturbation operator, with integration performed over the whole space. Since the operator $\widehat{\mathcal{W}}$ transforms a poloidal field into a toroidal one and vice versa (and the two are orthogonal), it follows that W_ii = 0 and each non-vanishing matrix element involves at least one toroidal and one poloidal free-decay eigenfunction. Because the toroidal eigenfunctions vanish at r >  r_h, the integrals are in fact restricted to the interior of the halo. The solvability condition of the system of equations for a_i, the vanishing of its determinant, yields the growth rate Γ. Once the matrix elements have been computed and the system Eq. (64) has been solved for a_i, Eq. (59) yields the solution of the dynamo equation. One of the coefficients a_i remains arbitrary because the dynamo equation is linear in magnetic field and hence its solution is determined up to an arbitrary factor. This freedom is used to fix the magnetic field strength at any desired value.

4.2. Parameters of galactic haloes

Velocity fields in galactic haloes are poorly known. Random velocities are likely to increase with altitude, and H I observations of Kalberla et al. (1998; see also Kalberla & Kerp 2009) suggest a three-dimensional velocity dispersion of about 100 km s⁻¹, close to the sound speed at a temperature 10⁶ K. The scale of these motions is uncertain. The size of supernova remnants above the galactic disc is expected to be of order 0.3 kpc (McKee & Ostriker 1977). The size of the hot gas bubbles rising from the disc and the scale of the Parker instability are of order 0.5–1 kpc (e.g. Rodrigues et al. 2016). Adopting the random speed and scale as v = 100 km s⁻¹ and l = 0.5 kpc, the turbulent diffusivity is estimated by $\beta {}_{}\text{h}\simeq \frac{1}{3}lv=5\times {10}^{27}{\mathrm{cm}}^2{\mathrm{s}}^{-1}$. The corresponding magnetic diffusion time across the halo radius is $r_{\text{h}}^2/{\beta }_{\text{h}}\simeq 1.4\times {10}^{10}\text{yr}$.

The knowledge of the variation of the rotation speed with position within galactic haloes is rather rudimentary. Both the rotational speed and its radial gradient decrease in spiral galaxies with distance from the mid-plane, with a typical vertical gradient of order ∂V/∂z = −(15–25) km s⁻¹ kpc⁻¹ within a few kiloparsecs from the mid-plane (Zschaechner et al. 2015). In our fiducial model, the halo is assumed to have a rotation curve of the form

$$\begin{array}{c}\hfill \mathit{V}(\mathit{r})=V_{\mathrm{h}}f(r,\theta )\widehat{\mathit{\varphi }}\end{array}$$(66)(expressed in terms of dimensional variables), with $\widehat{\mathit{\varphi }}$ the unit azimuthal vector and

$$\begin{array}{c}\hfill f(r,\theta )=\frac{1-exp\left(-s/s_v\right)}{1-exp(-r_{\mathrm{h}}/s_v)},\text{with}s=rsin\theta ,\end{array}$$(67)

where the turnover radius is chosen to be s_v = 3 kpc, the typical value found in observations and simulations of MW-type galaxies (Reyes et al. 2011; Schaller et al. 2015). For simplicity, the rotation curve of Eqs. (66) and (67) has no z-dependence but it can easily be introduced. The role of the variation of Ω with z is to produce B_ϕ from B_z, arguably a process somewhat less important than the stretching of the radial magnetic field in the azimuthal direction at a rate S = s∂Ω/∂s.

We adopt a simple form for the α-coefficient often used in spherical mean-field dynamo models; in dimensional variables,

$$\begin{array}{c}\hfill \alpha (\mathit{r})={\alpha }_{\mathrm{h}}cos\theta ,\end{array}$$(68)

implying the largest absolute value of α near the poles whilst α also vanishes at the equator, reflecting the fact that the mean helicity of the random flows is produced by the Coriolis force.

We consider an axially symmetric magnetic field in the halo and assume that the dynamo operates within a region of r_h = 15 kpc in radius. We take V_h = 220 km s⁻¹ (similar to that in the disc). With the turbulent magnetic diffusivity β_h = 5 × 10²⁷ cm² s⁻¹, this leads to R_ωh ≃ 200.

Estimating R_αh in the halo is more difficult given the uncertainty of the random flow parameters. The standard estimate of Eq. (11) yields α_h ≃ l²Ω/h ≃ 1 km s⁻¹ and R_αh = α_hr_h/β_h ≃ 1 for Ω = 26 km s⁻¹ kpc⁻¹ and h = 3 kpc, the gas density scale height in the halo. As the fiducial value for R_αh, we select its marginal value corresponding to the vanishing dynamo growth rate (see Sect. 4.3 for details): R_αh^(q) = 4.3 for symmetric solutions and R_αh^(d) = 8.1 for anti-symmetric ones. The symmetric mode is preferred to the anti-symmetric one only slightly, R_αh^(q)/R_αh^(d) ≃ 0.5 for the marginal values. The similarity of the marginal values of R_αh for the dipolar and quadrupolar magnetic structures in the halo reflects the fact that, unlike the disc dynamo, spherical dynamos usually do not exhibit a strong preference for either symmetry.

4.3. Basic magnetic structures

Figure 9 shows two examples of magnetic structures in the halo that are marginally stable with respect to the mean-field dynamo action, ∂B/∂t = 0, one symmetric with respect to the equator and the other anti-symmetric:

$$\begin{array}{cc}\hfill {\mathit{B}}_{}{\text{h}}^{\text{(d)}}& \approx -0.48{\mathit{B}}_1^{\text{(d)}}-0.38{\mathit{B}}_2^{\text{(d)}}-0.70{\mathit{B}}_3^{\text{(d)}}-0.12{\mathit{B}}_4^{\text{(d)}},\hfill \\ \hfill {\mathit{B}}_{}{\text{h}}^{\text{(q)}}& \approx 0.14{\mathit{B}}_1^{\text{(q)}}+0.86{\mathit{B}}_2^{\text{(q)}}+0.10{\mathit{B}}_3^{\text{(q)}}-0.41{\mathit{B}}_4^{\text{(q)}}.\hfill \end{array}$$

The eigenfunctions are normalised to have $B_{\text{h}}^{\text{(d)}}=-0.5\mu \text{G}$ and $B_{\text{h}}^{\text{(q)}}=-0.01\mu \text{G}$ at (s, z) = (8.5, 0.02) kpc in the anti-symmetric and symmetric cases, respectively, so that they have similar maximum magnetic field strengths.

The poloidal magnetic lines (the left-hand panels) have the so-called X shape detected in the halos of some galaxies, especially pronounced in the quadrupolar structure. This is a generic field structure typical of any divergence-free vector field that can be enhanced further by a large-scale velocity shear of the galactic wind or fountain. Unlike the symmetric eigenfunction, the anti-symmetric one has a maximum away from the equator. The position of the maxima depends on the spatial forms of α(r) and V(r); in our case, the rotation speed is independent of z, and α(r) alone controls this feature.

It is not clear which of the two symmetries may dominate in galactic haloes: this depends on the strength of the magnetic coupling between the disc and the halo and between the two hemispheres of the halo. If the disc-halo coupling is strong or the disc disrupts magnetic connection between the two halo hemispheres, the quadrupolar disc field could enforce a symmetric field structure in the halo. Halo fields of mixed parity are also a possibility, but their modelling requires non-linear dynamo solutions rather than superpositions of linear eigenmodes that we use here. The strength of the disc-halo magnetic connection depends on the ratio of the turbulent magnetic diffusivities in the two regions: the larger the value of β_h/β_d, the weaker the coupling. Existing models of the mean-field dynamo action in galactic disc-halo systems only considered the range β_h/β_d ≤ 30.

The top two panels of Fig. 10 illustrate how the growth rates and oscillation frequencies of the symmetric and anti-symmetric modes depend on R_αh when all other parameters are fixed to the fiducial values of Table 1. These solutions involve the first four symmetric or anti-symmetric free-decay modes with expansion coefficients shown in the lower half of Fig. 10. As shown in Fig. 10b, both the symmetric and anti-symmetric eigenmodes are typically oscillatory (ImΓ ≠ 0). The fact that ImΓ = 0 for the anti-symmetric mode at R_αh ≳ 10 appears to be an artefact of including only a small number of the free-decay modes into the perturbation series. The series Eq. (59) converges rather slowly (Rädler & Wiedemann 1989; Rädler et al. 1990) and adding a few more terms does not always improve the accuracy (Sokoloff et al. 2008). Therefore, the model for the halo magnetic field that involves only a modest number of modes can reproduce only relatively simple magnetic configurations (and yet quite non-trivial – see Fig. 11). This does not appear to be a serious problem, though, since the scale of the mean magnetic field in galactic haloes is unlikely to be smaller than a few kiloparsecs.

The dependencies of the growth rate and oscillation frequency of the magnetic field on the turnover radius of the rotation curve s_v are shown in Fig. 12. Larger values of s_v correspond to weaker differential rotation and, therefore, lower growth rates and oscillation frequencies.

5. Discussion

Magnetic fields obtained above for the disc and halo are combined by a simple superposition with arbitrary weights. For illustration, we show in Fig. 11 a magnetic structure that has two radial field reversals in the disc (Model B of Fig. 7) and a symmetric (quadrupolar) field in the halo shown in the bottom row of Fig. 9. The complexity of the resulting magnetic structure clearly illustrates the possibilities of the model. The most important limitations of the model in its current form is that it is axially symmetric and does not include galactic outflows. Both can be addressed rather straightforwardly within the framework of this approach. In the following we discuss some applications and extensions of our approach.

5.1. Synthetic radio maps

The model can be used to interpret observations of synchrotron emission and Faraday rotation as soon as the distributions of cosmic ray and thermal electrons have been specified. The total and polarised synchrotron intensities, polarisation angle and Faraday rotation measure can be derived as described in Appendix C. As a simple illustrative example, we assume a uniform distribution for the cosmic-ray electrons, n_γ = const in terms of the number density, and adopt exponential profiles for the number density of thermal electrons,

$$\begin{array}{c}\hfill n_{\mathrm{e}}(s,z,\varphi )=n_{0}exp[-\frac{z}{h(s)}-\frac{s}{s_{}\text{e}}],\end{array}$$(69)

where h(s) is given by Eq. (3) and s_e = 3 kpc, is the scale radius of the disc (chosen to be similar to the case of the Milky Way, Binney & Tremaine 2008). Any other model (e.g. Cordes & Lazio 2002) could be used instead but we prefer to avoid exaggerating the amount of detail in the magnetic field model that may arise from details such as spiral arms in more complicated models for n_e.

In the top panel Fig. 13, we show the synchrotron emission in a galaxy seen edge-on with the magnetic field of Fig. 11. The other two panels show the polarised emission at two wavelengths, λ = 5 cm and λ = 20 cm, as in C- and L-bands of the VLA used, for example, in the CHANG-ES survey (Irwin et al. 2012). At 5 cm, most of the polarisation signal is dominated by the disc component and localised around the mid-plane of the disc. At longer wavelengths, most of the emission from the galactic plane is depolarised and two conical lobes of about 5 kpc in height are prominent in the halo, similar to the so-called X-shaped structures observed in edge-on galaxies (Wiegert et al. 2015).

5.2. Evolution of galactic magnetic fields

Another possible application is a simple, approximate model for the evolution of a large-scale magnetic field in a galaxy. In this application, an initial (seed) magnetic field has to be prescribed and then it can be evolved using the growth rates of the magnetic modes derived above. Suitable initial conditions must then be selected. The simplest approach is to assume that all the modes are equally represented in the initial state, that is, C_n independent of n and chosen to obtain an initial large-scale magnetic field of any given strength. A physically better motivated initial magnetic field represents a random field produced by the fluctuation dynamo in a young galaxy or protogalaxy (Poezd et al. 1993). Because of the finite size of the dynamo region, the projection of such a random field onto the dynamo eigenmodes does not vanish. Ruzmaikin et al. (Sect. VII.14 in 1988) estimate the corresponding initial dimensional magnitudes of the radial disc modes as

$$\begin{array}{c}\hfill C_n^{(0)}=\frac{b}{N_n^{1/2}}\frac{l}{\delta s_n},\end{array}$$(70)

where b is the root-mean square strength of the random magnetic field b, N_n is the number of the correlation cells of b within a cylindrical annulus of an axial extent 2h, radius s and width δs_n, with δs_n the radial scale of Q_n(s), and l is the scale of b. With δs ≃ s_disc/n, N ≃ hss_disc/(nl³), l = 100 pc, s_disc = 20 kpc, we have

$$\begin{array}{c}\hfill C_n^{(0)}\simeq n^{3/2}{\left(\frac{s}{s_{}\text{disc}}\right)}^{-1/2}{\left(\frac{l^5}{h{s}_{}{\text{disc}}^4}\right)}^{1/2}\simeq {10}^{-5}{n}^{3/2}b{\left(\frac{s}{s_{}\text{disc}}\right)}^{-1/2}·\end{array}$$(71)

Thus, the seed for the large-scale dynamo due to the small-scale magnetic field favours higher-order modes being proportional to n^3/2, because they have smaller scale, and decreases with radius as s^−1/2. A plausible estimate is b ≃ 5 μG by analogy with observational estimates for nearby spiral galaxies. Otherwise, if a dependence on the interstellar gas parameters is required, a suitable estimate is

$$\begin{array}{c}\hfill b\simeq {(4\pi \rho )}^{1/2}v,\end{array}$$(72)

where ρ is the gas density in the diffuse warm interstellar gas and v is the turbulent speed. The standard estimates of the latter are ρ ≃ 1, 7 × 10⁻²⁴ g cm⁻³ corresponding to the number density of 1 cm⁻³ and v ≃ 10 km s⁻¹. This yields b ≃ 3 μG.

5.3. Extensions of the model

There are several directions in which the model can be extended. Perhaps most important is to include non-axisymmetric magnetic fields. This is straightforward to implement. In the disc, the local equations of Sect. 3.3 and their solutions remain unchanged but the radial part of the eigenfunction Q_n of Sect. 3.4 becomes a function of both radius and azimuth. Solutions for Q_m(s, ϕ) were obtained in largely the same manner as above by Baryshnikova et al. (1987), Krasheninnikova et al. (1989), and Bykov et al. (1997), and are reviewed by Ruzmaikin et al. (Sect. VII.8 in 1988) and Krasheninnikova et al. (1990). Introducing non-axisymmetric magnetic fields in the halo would only require that non-axisymmetric free-decay modes are included into the perturbation solution. This is straightforward to do and does not require any significant modification of the formalism of Sect. 4.1.

Another physically important generalisation is the inclusion of galactic outflows and accretion flows, that is, large-scale poloidal velocity fields U. The additional velocity components appear in the perturbation operators. Within the disc, U_z enters the local Eqs. (18) and (19) whereas U_s is included in the radial Eq. (40). The modified solutions are discussed by Bardou et al. (2001) and Moss et al. (2000), respectively. In the halo, the poloidal velocity just enters the perturbation operator Eq. (63) without affecting the procedure of perturbation analysis.

The solutions used in the model are kinematic (linear in magnetic field) as they are derived for V, α and β independent of B. The linear nature of the solution is not restrictive in the present context since its aim is to provide a convenient functional basis to parametrise a desired magnetic configuration. On the other hand, Chamandy et al. (2014) show that a wide class of non-linear solutions are well approximated by the marginally stable eigenfunction (i.e., that obtained for ∂B/∂t = 0). Non-linear dynamo effects, leading to solutions sensitive to the gas density and other relevant parameters, can be introduced in the radial thin-disc Eq. (40) as discussed by Poezd et al. (1993), via a non-linear modification (quenching) of the local growth rate which becomes a function of Q:

$$\begin{array}{c}\hfill \gamma (s,Q)=\gamma (s)[1-Q^2/B_0^2(s)],\end{array}$$

where γ(s) is the kinematic local growth rate obtained as discussed in Sect. 3.3. Since the time scale of magnetic field evolution in the halo is comparable to 10¹⁰ yr (Sect. 4), non-linear dynamo effects are likely to be less important in galactic haloes.

5.4. The GALMAG software package

The model presented in this paper has been implemented as the Python software package GALMAG¹ (Rodrigues 2018), which is publicly available under the GNU General Public License v3. Further details can be found in the on-line code documentation², which includes a tutorial.

Since GALMAG uses Python objects of the D2O package (Steininger et al. 2016) instead of regular NUMPY arrays, when it is invoked using MPI, all the array operations are automatically performed in parallel. As a stand-alone package, it can synthesise three-dimensional magnetic field structures from a provided set of expansion coefficients and compute synthetic maps of the Stokes parameters of synchrotron emission and Faraday rotation.

A more flexible use of GALMAG is to employ it in modular magnetic field optimisation frameworks like the IMAGINE pipeline (Steininger et al. 2018; Steininger 2018), where it serves as a magnetic field generator and is interfaced to multi-purpose observable generators, such as the HAMMURABI code³ (Waelkens et al. 2009). This allows us to not only compute maps of observables from any point of view and thus to compare with observations, but also provides sophisticated sampling techniques to optimise the GALMAG parameters.

6. Conclusions

We have presented an approach to develop parametrised models for large scale magnetic fields of the Milky Way and other disc galaxies based on fundamental equations of magnetic field generation and evolution. Implemented in the software package GALMAG, it is designed to be used in interpretations of observations of Faraday rotation, synchrotron and dust emission, and other observational tracers. In this paper, we have presented the basic formalism of the approach, and demonstrated its capabilities in illustrative examples.

The model is based on the expansion of the large-scale magnetic field over a basis of eigenfunctions of the mean-field dynamo Eq. (1), and the standard induction equation is its special case obtained for α = 0. As long as the functional basis is complete, any magnetic structure, whether or not produced by the dynamo, can be represented as a superposition of the eigenfunctions. Therefore, an alternative use of the magnetic field model is to represent any magnetic configuration of interest in terms of a relatively small number of parameters. The resulting magnetic field is physically realisable, being a solution of the induction equation or its modification with α ≠ 0, as desired. Furthermore, the fact that the model parameters have clear physical meaning would help to refine it so as to satisfy any additional constraints. Magnetic fields of the model are obtained in the form of series expansions, and the series can be truncated to achieve the desired amount of detail in the resulting solution.

The novelty and strength of our approach lies in advancing both the flexibility and physical plausibility of GMF models. It will be useful in Bayesian optimization machines, which seek among many reasonable morphological approaches to the GMF structure the one that gives the best results both in terms of physical plausibility and the explanation of existing data.

---

Acknowledgments

We have benefited from fruitful discussions with the members of the ISSI International Team 323, Bayesian modeling of the Galactic magnetic field constrained by space- and ground-based radio-millimetre and ultra-high energy cosmic ray data (http://www.issibern.ch/teams/bayesianmodel/), and the IMAGINE Consortium (https://www.astro.ru.nl/imagine/). Special thanks are due to Torsten Enßlin, Marijke Haverkorn, Jens Jasche and Andrew Fletcher for their useful comments, and Theo Steininger for help with the software development and the D2O Python package AS, LFSR and JPR acknowledge financial support and hospitality of the International Space Science Institute (ISSI) in Bern, Switzerland. AS and LFSR are supported by STFC (ST/N000900/1, Project 2), and AS and PB are supported by the Leverhulme Trust (RPG-2014-427). This research has made use of NASA’s Astrophysics Data System.

References

Appendix A: Spherical free-decay modes

For axisymmetric free-decay modes, a solution to Eq. (60) is obtained in terms of scalar potentials, as discussed by Krause & Rädler (1980) and Moffatt (1978). For the reader’s convenience, we present an outline of the solution. Any magnetic field B can be represented as the sum of a poloidal field ∇ × A_P, where A_P is its vector potential, and a toroidal field B_T:

$$\begin{array}{c}\hfill \mathit{B}=\mathrm{\nabla }\times {\mathit{A}}_{}\text{P}+{\mathit{B}}_{}\text{T},\end{array}$$(A.1)

and ∇ ⋅ B = 0 provided

$$\begin{array}{c}\hfill {\mathit{A}}_{\text{P}}=-\mathit{r}\times \mathrm{\nabla }S,{\mathit{B}}_{\text{T}}=-\mathit{r}\times \mathrm{\nabla }T,\end{array}$$(A.2)

where r is the position vector normalised such that r = 1 is the halo surface (r = r_h in dimensional variables) and S and T are known as the scalar potentials. In terms of the scalar potentials and assuming axial symmetry, Eq. (60) reduces in spherical coordinates (r, θ, ϕ) to

$$\begin{array}{c}\hfill \frac{1}{r^2}\frac{\partial }{\partial r}{\displaystyle \left(r^2\frac{\partial S}{\partial r}\right)+\frac{1}{r^{2}sin\theta }\frac{\partial }{\partial \theta }\left(sin\theta \frac{\partial S}{\partial \theta }\right)}=\gamma S,\end{array}$$(A.3)

$$\begin{array}{c}\hfill \frac{1}{r^2}\frac{\partial }{\partial r}{\displaystyle \left(r^2\frac{\partial T}{\partial r}\right)+\frac{1}{r^{2}sin\theta }\frac{\partial }{\partial \theta }\left(sin\theta \frac{\partial T}{\partial \theta }\right)}=\gamma T,\end{array}$$(A.4)

for r <  1 and,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}^{2}S=0T=0\text{for}r>1,\end{array}$$(A.5)

with the vacuum boundary conditions

$$\begin{array}{c}\hfill [S]=[\partial S/\partial r]=[T]=0\text{at}r=1,\end{array}$$(A.6)

where [X] denotes the jump of X, and [X]=0 means continuity. We also require both potentials to be finite at r = 0.

The potentials satisfy identical equations at r <  1, so consider this equation for G equal to either S or T,

$$\begin{array}{c}\hfill \frac{1}{r^2}\frac{\partial }{\partial r}{\displaystyle \left(r^2\frac{\partial G}{\partial r}\right)+\frac{1}{r^{2}sin\theta }\frac{\partial }{\partial \theta }\left(sin\theta \frac{\partial G}{\partial \theta }\right)}-\gamma G=0.\end{array}$$(A.7)

Using separation of variables, G(r, θ)=R(r)Θ(θ), Bessel’s equation is obtained in r and Legendre’s equation in θ, with the separation constant n(n + 1) (n = 1, 2, 3, …):

$$\begin{array}{c}\hfill r^2\frac{{\mathrm{d}}^{2}R}{\mathrm{d}{r}^2}+2r\frac{\mathrm{d}R}{\mathrm{d}r}-[\gamma r^2+n(n+1)]R=0,\end{array}$$(A.8)

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\theta }(sin\theta \frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}\theta })+n(n+1)\mathrm{\Theta }sin\theta =0.\end{array}$$(A.9)

In terms of $x=\sqrt{-\gamma }r$ and Q(x)=x^1/2R(x) in Eq. (A.8) and x = cosθ in Eq. (A.9), we have

$$\begin{array}{c}\hfill x^2\frac{{\mathrm{d}}^{2}Q}{\mathrm{d}{x}^2}+x\frac{\mathrm{d}Q}{\mathrm{d}x}+[x^2-{(n+\frac{1}{2})}^2]Q=0,\end{array}$$(A.10)

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}x}\left[(1-x^2)\frac{\mathrm{d}\mathrm{\Theta }}{\mathrm{d}x}\right]+n(n+1)\mathrm{\Theta }=0.\end{array}$$(A.11)

Non-singular solutions of Eqs. (A.3) and (A.4) then follow as

$$\begin{array}{cc}\hfill T& ={\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \sum _{l=1}^{\infty }}{c}_{\mathit{nl}}{T}_{\mathit{nl}}(r)P_n(cos\theta ),\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill S& =r_{}\text{h}{\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \sum _{l=1}^{\infty }}{c}_{\mathit{nl}}{S}_{\mathit{nl}}(r)P_n(cos\theta ),\hfill \end{array}$$(A.13)

where

$$\begin{array}{c}\hfill T_{\mathit{nl}}(r)=S_{\mathit{nl}}(r)=\frac{1}{{\xi }_{\mathit{nl}}\sqrt{r}}{J}_{n+1/2}({\xi }_{\mathit{nl}}r),\end{array}$$(A.14)

with constants c_nl,

$$\begin{array}{c}\hfill {\xi }_{\mathit{nl}}=\sqrt{-{\gamma }_{\mathit{nl}}},\end{array}$$(A.15)

and ξ_nl are solutions to (A.18). The factor r_h in Eq. (A.13) is introduced to ensure dimensional consistency when obtaining the magnetic field from these potentials. The boundary conditions (A.6) reduce to

$$\begin{array}{c}\hfill T_{\mathit{nl}}=0,S_{\mathit{nl}}=d_n,\frac{\partial S_{\mathit{nl}}}{\partial r}=-(n+1)d_n\text{at}r=1,\end{array}$$(A.16)

where d_n are constants. Eliminating d_n, the boundary conditions for S_nl reduce to the recurrence relation

$$\begin{array}{c}\hfill \frac{\partial S_{\mathit{nl}}}{\partial r}+(n+1)S_{\mathit{nl}}=0\text{at}r=1.\end{array}$$(A.17)

Together with the requirement that T_nl and S_nl do not vanish simultaneously, this gives

$$\begin{array}{c}\hfill J_{n-1/2}({\xi }_{\mathit{nl}})J_{n+1/2}({\xi }_{\mathit{nl}})=0,\end{array}$$(A.18)

which determines the admissible values for ξ_nl, and ${\gamma }_{nl}=-{\xi }_{nl}^2$ yields the decay rates γ_nl given in Table A.1.

For l odd, J_{n − 1/2}(ξ_nl)=0 for all n. Hence, T_nl = 0 for l odd. Conversely, when l is even, J_n + 1/2(ξ_nl)=0 for all n. Hence, S_nl = 0 for l even. The solutions satisfying the boundary conditions can be written as follows:

$$\begin{array}{cc}\hfill T& ={\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \sum _{l\text{even}}}\frac{c_{\mathit{nl}}}{\sqrt{r}}{J}_{n+1/2}({\xi }_{\mathit{nl}}r)P_n(cos\theta ),\hfill \end{array}$$(A.19a)

$$\begin{array}{cc}\hfill S& ={\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \sum _{l\text{odd}}}\frac{d_{\mathit{nl}}}{\sqrt{r}}{J}_{n+1/2}({\xi }_{\mathit{nl}}r)P_n(cos\theta ),\hfill \end{array}$$(A.19b)

where c_nl and d_nl are constants. Individual terms in the sums are arranged according to increasing magnitude of the decay rates, |γ_nl|.

The free-decay modes form two separate families based on their symmetry about the equator θ = π/2, the anti-symmetric (dipolar) and symmetric (quadrupolar) ones. The anti-symmetric modes, denoted with superscript (d), occur when both n and l are either odd or even, whilst the symmetric modes that have superscript (q) occur otherwise.

Explicit forms of a few of the lowest free-decay modes are given in the next section, each normalised to

$$\begin{array}{c}\hfill {\displaystyle {\int }_V}{|{\mathit{B}}_n|}^2{\mathrm{d}}^3\mathit{r}=1,\end{array}$$(A.20)

where the integral is taken over the sphere r ≤ 1, to form an orthonormal set. Although each eigenmode is either poloidal or toroidal, their superpositions Eq. (A.20) necessarily contain both poloidal and toroidal parts: purely toroidal and purely poloidal fields cannot sustain Ohmic dissipation and unavoidably decay.

A.1. Symmetric modes

The four leading quadrupolar free-decay modes are shown in the upper row of Fig. 8. The quadrupolar mode of the slowest decay has (n, l)=(2, 1) and is poloidal,

$$\begin{array}{c}\hfill {\mathit{B}}_1^{\text{(q)}}=A_1(\frac{Q_1(r)}{r}(3{cos}^2\theta -1),-\frac{sin\theta cos\theta }{r}\frac{\mathrm{d}}{\mathrm{d}r}[r{Q}_1(r)],0),\end{array}$$(A.21)

where A₁ ≈ 0.662 and

$$\begin{array}{c}\hfill Q_1(r)=\{\begin{array}{cc}{r}^{-1/2}{J}_{5/2}(q_{1}r),\hfill & r\le 1,\hfill \\ r^{-3}{J}_{5/2}(q_1),\hfill & r>1,\hfill \end{array}{q}_1\approx 4.493.\end{array}$$(A.22)

The next mode, (n, l)=(1, 2), is toroidal and has the same eigenvalue,

$$\begin{array}{c}\hfill {\mathit{B}}_2^{\text{(q)}}=A_2(0,0,Q_2(r)sin\theta ),\end{array}$$(A.23)

where A₂ ≈ 1.330 and

$$\begin{array}{c}\hfill Q_2(r)=\{\begin{array}{cc}{r}^{-1/2}{J}_{3/2}(q_{1}r),\hfill & r\le 1,\hfill \\ r^{-2}{J}_{3/2}(q_{1}r),\hfill & r>1.\hfill \end{array}\end{array}$$(A.24)

The modes ${\mathit{B}}_3^{\text{(q)}}$ and ${\mathit{B}}_4^{\text{(q)}}$, poloidal and toroidal respectively, also form a doublet with the common eigenvalue and correspond to (n, l)=(4, 1) and (n, l)=(3, 2), respectively:

$$\begin{array}{c}\hfill {\mathit{B}}_3^{\text{(q)}}=A_3(-20\frac{Q_3(r)}{r}{S}_1(\theta ),-r^{-1}\frac{\mathrm{d}}{\mathrm{d}r}[r{Q}_3(r)]\frac{\mathrm{d}{S}_1(\theta )}{\mathrm{d}\theta },0),\end{array}$$(A.25)

where A₃ ≈ 0.133, S₁(θ)=35cos⁴θ − 30cos²θ + 3 and

$$\begin{array}{c}\hfill Q_3(r)=\{\begin{array}{cc}{r}^{-1/2}{J}_{9/2}(q_{3}r),\hfill & r\le 1,\hfill \\ r^{-5}{J}_{9/2}(q_3),\hfill & r>1,\hfill \end{array}{q}_3\approx 6.988,\end{array}$$(A.26)

and

$$\begin{array}{c}\hfill {\mathit{B}}_4^{\text{(q)}}=A_4(0,0,-Q_4(r)\frac{\mathrm{d}{S}_2(\theta )}{\mathrm{d}\theta }),\end{array}$$(A.27)

where A₄ ≈ 0.763, S₂(θ)=5cos³θ − 3cosθ and

$$\begin{array}{c}\hfill Q_4(r)=\{\begin{array}{cc}{r}^{-1/2}{J}_{7/2}(q_{3}r),\hfill & r\le 1,\hfill \\ r^{-4}{J}_{7/2}(q_3),\hfill & r>1.\hfill \end{array}\end{array}$$(A.28)

A.2. Anti-symmetric modes

The spherical components of magnetic field in a few leading anti-symmetric modes have the following form, illustrated in the bottom row of Fig. 8.

The mode that decays most slowly is poloidal, with (n, l)=(1, 1):

$$\begin{array}{c}\hfill {\mathit{B}}_1^{(\mathrm{d})}=C_1{\displaystyle (\frac{2}{r}{Q}_1(r)cos\theta ,-\frac{sin\theta }{r}\frac{\mathrm{d}}{\mathrm{d}r}[r{Q}_1(r)],0)},\end{array}$$(A.29)

where C₁ ≈ 0.346 and

$$\begin{array}{c}\hfill Q_1(r)=\{\begin{array}{cc}{r}^{-1/2}{J}_{3/2}(k_{1}r),\hfill & r\le 1,\hfill \\ r^{-2}{J}_{3/2}(k_1),\hfill & r>1.\hfill \end{array}{k}_1=\pi ,\end{array}$$(A.30)

The next two modes ${\mathit{B}}_2^{(\text{d})}$ and ${\mathit{B}}_3^{(\text{d})}$, poloidal and toroidal with (n, l)=(3, 1) and (n, l)=(2, 2), respectively, form a degenerate pair:

$$\begin{array}{cc}\hfill {\mathit{B}}_2^{(\mathrm{d})}& =C_2(\frac{2cos\theta }{r}(5cos2\theta -1)Q_2(r),\hfill \\ \hfill & -\frac{sin\theta }{r}(5{cos}^2\theta -1)\frac{\mathrm{d}}{\mathrm{d}r}\left[r{Q}_2(r)\right],0),\hfill \end{array}$$(A.31)

where C₂ ≈ 0.250 and

$$\begin{array}{c}\hfill Q_2(r)=\{\begin{array}{cc}{r}^{-1/2}{J}_{7/2}(k_{2}r),\hfill & r\le 1,\hfill \\ r^{-4}{J}_{7/2}(k_2),\hfill & r>1,\hfill \end{array}{k}_2\approx 5.763.\end{array}$$(A.32)

The toroidal mode of the doublet has the form

$$\begin{array}{c}\hfill {\mathit{B}}_3^{(\mathrm{d})}=C_3(0,0,Q_3(r)sin\theta cos\theta ),\end{array}$$(A.33)

where C₃ ≈ 3.445 and

$$\begin{array}{c}\hfill Q_3(r)=\{\begin{array}{cc}{r}^{-1/2}{J}_{5/2}(k_{2}r),\hfill & r\le 1,\hfill \\ r^{-3}{J}_{5/2}(k_{2}r),\hfill & r>1.\hfill \end{array}\end{array}$$(A.34)

The fourth antisymmetric mode is also poloidal, with (n, l)=(1, 3):

$$\begin{array}{c}\hfill {\mathit{B}}_4^{(\mathrm{d})}=C_4(\frac{2}{r}{Q}_4(r)cos\theta ,-\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}[r{Q}_4(r)]sin\theta ,0),\end{array}$$(A.35)

where C₄ ≈ 0.244 and

$$\begin{array}{c}\hfill Q_4(r)=\{\begin{array}{cc}{r}^{-1/2}{J}_{3/2}(k_{4}r),\hfill & r\le 1,\hfill \\ r^{-2}{J}_{3/2}(k_4),\hfill & r>1,\hfill \end{array}{k}_4=2\pi .\end{array}$$(A.36)

Appendix B: The αω-dynamo approximation for the halo

In the main part of this paper, we use an α²ω-dynamo model for the halo, where the toroidal magnetic field is produced from the poloidal one by both differential rotation and the α-effect. The contribution of the α²-dynamo to the generation of the toroidal field is usually weaker and therefore often neglected to simplify the solutions. To assess the consequences of this approximation, here we provide solutions of the αω-dynamo in the halo which should be compared with results presented in Sect. 4.3.

For the αω-dynamo, the operator defined in Eq. (63) is simplified to

$$\begin{array}{c}\hfill \widehat{\mathcal{W}}\mathit{B}=R_{\alpha \text{h}}(\mathrm{\nabla }\times (\alpha \mathit{B})-{[\mathrm{\nabla }\times (\alpha \mathit{B})]}_{\varphi })+R_{\omega \text{h}}\mathrm{\nabla }\times (\mathit{V}\times \mathit{B}),\end{array}$$(B.1)

thus removing the contribution of the α-effect to the azimuthal magnetic field. In Fig. A.1, similar to Fig. 10, we show the resulting eigenvalues and expansion coefficients of the perturbation solution. The difference between the two solutions is noticeable.

The α²-dynamo is negligible in comparison with the αω mechanism when |R_ω/R_α|≫1 or even $|R_{\omega }/R_{\alpha }^2|\gg 1$ (Ruzmaikin et al. 1980). For the fiducial value of parameters, |R_ω/R_α|≃140 in the disc and 50–100 in the halo (and |R_ωh/R_αh² = 10–50). Thus, unlike the case of galactic discs with strong differential rotation, the αω-approximation cannot be recommended for the halo because it neglects a potentially important part of the dynamo mechanism.

Appendix C: Synchrotron emission and Faraday rotation

For the line of sight along the x-axis of a Cartesian reference frame r = (x, y, z), the synchrotron emissivity at a wavelength λ is derived as

$$\begin{array}{c}\hfill ϵ(\mathit{r},\lambda )\propto {[B_y^2(\mathit{r})+B_z^2(\mathit{r})]}^{(\kappa +1)/4}{\lambda }^{(\kappa -1)/2},\end{array}$$(C.1)

assuming a uniform distribution of cosmic ray electrons and the cosmic ray energy spectrum N(E) dE ∝ E^−κ dE with κ = 3. The Stokes parameters are computed as

$$\begin{array}{cc}\hfill I(y,z,\lambda )& ={\displaystyle {\int }_{-\infty }^{\infty }}ϵ(x^{\prime },y,z,\lambda )\mathrm{d}{x}^{\prime },\hfill \end{array}$$(C.2)

$$\begin{array}{cc}\hfill Q(y,z,\lambda )& =p_0{\displaystyle {\int }_{-\infty }^{\infty }}ϵ(x^{\prime },y,z,\lambda )cos[2\psi (x^{\prime },y,z)]\mathrm{d}{x}^{\prime },\hfill \end{array}$$(C.3)

$$\begin{array}{cc}\hfill U(y,z,\lambda )& =p_0{\displaystyle {\int }_{-\infty }^{\infty }}ϵ(x^{\prime },y,z,\lambda )sin[2\psi (x^{\prime },y,z)]\mathrm{d}{x}^{\prime },\hfill \end{array}$$(C.4)

with the intrinsic polarisation degree p₀ = 0.75, and the local polarisation angle ψ(r) is obtained from

$$\begin{array}{cc}\hfill \psi (\mathit{r})=& \frac{\pi }{2}+arctan\left[\frac{B_z(\mathit{r})}{B_y(\mathit{r})}\right]\hfill \\ \hfill & +0.81R_{\alpha \text{d}}{\left(\frac{\lambda }{1\mathrm{m}}\right)}^2{\displaystyle {\int }_x^{\infty }}\frac{n_{\mathrm{e}}(x^{\prime },y,z)}{1{\mathrm{cm}}^{-3}}\frac{B_x({\mathit{r}}^{\prime })}{1\mu \mathrm{G}}\frac{\mathrm{d}{x}^{\prime }}{1\mathrm{pc}},\hfill \end{array}$$(C.5)

where the thermal electron density n_e(r) is given by Eq. (69) in both the disc and the halo. The polarised intensity, observed polarisation angle and fractional polarisation follow as

$$\begin{array}{c}\hfill P=\sqrt{Q^2+U^2},\mathrm{\Psi }=\frac{1}{2}arctan(U/Q),p=P/I.\end{array}$$(C.6)

The Faraday rotation measure is calculated as RM = ∂Ψ/∂(λ²).

All Tables

All Figures