[1,3,4]Maarit J. Korpi-Lagg

1]Astroinformatics, Department of Computer Science, Aalto University, P.O. Box
15400, FI-00076, Espoo, Finland 2]Department of Astrophysics, American Museum
of Natural History, 200 Central Park West, New York, NY 10024, USA 3]Nordic
Institute for Theoretical Physics, Hannes Alfvéns väg 12, SE-106 91 Stockholm,
Sweden 4]Max Planck Institute for Solar System Research, Justus-von-Liebig-Weg
3, D-37707 Göttingen, Germany 5]School of Mathematics, Statistics and Physics,
Newcastle University, NE1 7RU, Newcastle upon Tyne, UK

# Computational approaches to modeling dynamos in galaxies

<EMAIL_ADDRESS>Mordecai-Mark Mac Low<EMAIL_ADDRESS>Frederick
A. Gent<EMAIL_ADDRESS>[  [  [  [  [

###### Abstract

Galaxies are observed to host magnetic fields with a typical total strength of
around 15$\mu$G. A coherent large-scale field constitutes up to a few
microgauss of the total, while the rest is built from strong magnetic
fluctuations over a wide range of spatial scales. This represents sufficient
magnetic energy for it to be dynamically significant. Several questions
immediately arise: What is the physical mechanism that gives rise to such
magnetic fields? How do these magnetic fields affect the formation and
evolution of galaxies? In which physical processes do magnetic fields play a
role, and how can that role be characterized? Numerical modelling of
magnetized flows in galaxies is playing an ever-increasing role in finding
those answers. We review major techniques used for these models. Current
results strongly support the conclusion that field growth occurs during the
formation of the first galaxies on timescales shorter than their accretion
timescales due to small-scale turbulent dynamos. The saturated small-scale
dynamo maintains field strengths at a few percent of equipartition with
turbulence. The subsequent action of large-scale dynamos in differentially
rotating discs produces observed modern field strengths in equipartition with
the turbulence and having power at large scales. The field structure resulting
appears consistent with observations including Faraday rotation and
polarisation from synchrotron and dust thermal emission. Major remaining
challenges include scaling numerical models toward realistic scale separations
and Prandtl and Reynolds numbers.

###### keywords:

Dynamo, Galaxy, Numerical methods, Magnetic fields, Interstellar medium

###### Contents

  1. 1 Introduction
    1. 1.1 Fundamental observations of galactic magnetic fields
    2. 1.2 Governing physical equations
    3. 1.3 Basics and challenges of dynamo theory
      1. 1.3.1 Large-scale dynamo instability
      2. 1.3.2 Small-scale dynamo instability
    4. 1.4 Seed fields
  2. 2 Numerical approaches
    1. 2.1 Mean-field modelling
    2. 2.2 Direct numerical simulation versus large eddy schemes
    3. 2.3 Spectral methods
    4. 2.4 Grid methods
      1. 2.4.1 Uniform grid methods
      2. 2.4.2 Hierarchical grid methods
    5. 2.5 Particle methods
    6. 2.6 Hybrid and kinetic methods
  3. 3 Results from numerical approaches
    1. 3.1 Plasma scale
    2. 3.2 Star formation scale
    3. 3.3 Kiloparsec scale
      1. 3.3.1 SN driven SSD
      2. 3.3.2 Large-scale fields without LSD
      3. 3.3.3 SN driven LSD
      4. 3.3.4 Magnetic buoyancy instability
      5. 3.3.5 Magnetorotational instability
    4. 3.4 Galactic scale
      1. 3.4.1 Mean-field (MF) models
      2. 3.4.2 MHD models
    5. 3.5 Cosmological scale
    6. 3.6 Observational comparisons
  4. 4 Properties of turbulence
    1. 4.1 Retrieval of turbulent transport coefficients
    2. 4.2 Measuring helicity fluxes
  5. 5 Conclusions
    1. 5.1 Comparisons with analytic theory
    2. 5.2 Consistency with observations
    3. 5.3 Challenges

##  1 Introduction

Galaxies are observed to host magnetic fields with a typical total strength of
around 15$\mu$G. A coherent large-scale field constitutes up to a few
microgauss of the total, while the rest is built from strong magnetic
fluctuations over a wide range of spatial scales (Beck, 2015; Beck et al,
2019). Until recently, the importance of magnetic fields has often been
overlooked when considering the dynamics and evolution of galaxies. However,
magnetic field energy is sufficiently large to have a significance comparable
to kinetic energy. Several questions immediately arise. What are the physical
mechanisms that give rise to such strong magnetic fields? How do different
physical processes in galaxies contribute to magnetic field growth? What are
the small and large-scale structures of these fields? It is becoming
increasingly clear that both small-scale dynamo (SSD) and large-scale dynamo
(LSD) processes play a fundamental role in determining the answers to these
questions. Progress in understanding dynamo behavior relies ever more strongly
on numerical modelling of the magnetised plasma in galaxies. This review
summarises the current state-of-the-art of these numerical simulations and
their results. See also the recent review by Brandenburg and Ntormousi (2022).

###  1.1 Fundamental observations of galactic magnetic fields

In galaxies, both the magnetic field and the velocity have structure across a
wide spectrum of characteristic scales. First, there are contributions from
the largest scales, such as the differential angular velocity that arises from
the typical flat rotation curves with nearly constant rotation velocity as a
function of radius observed in spiral galaxies embedded in dark matter halos.
At slightly smaller scales, galactic outflows from galaxy centers or from
powerful superbubbles produced by correlated SN explosions from young star
clusters produce ordered velocities in the vertical direction outward from the
galactic disc plane. Similarly, the magnetic field of spiral galaxies has been
verified to have in most cases a coherent large-scale structure that
sometimes, but by no means always, follows the optical spiral arms formed by
the stars.

We have a variety of indirect observations through which we can identify the
properties of the magnetic field in the interstellar medium (ISM),
circumgalactic medium, intracluster medium (ICM) and, more locally, in
molecular clouds, within and between spiral arms and around supernova (SN)
remnants or superbubbles (see, e.g., Table 1 of Beck, 2015). The total
magnetic field $B_{\rm tot_{\perp}}$ perpendicular to the line of sight is
measured through total synchrotron intensity. This is the sum of the regular
field $B_{\rm reg_{\perp}}$, which is defined as consistent within the beam
width of the telescope, and the locally turbulent field $B_{\rm
turb_{\perp}}$, which are the fluctuations on top of the regular field within
the width of the beam, whether due to SSD, tangling or both. The regular field
is discriminated from the total field through polarised synchrotron emission.
The sum of isotropic turbulence $B_{\rm iso_{\perp}}$, with dispersions
matching in all three dimensions, and anisotropic turbulence $B_{\rm
aniso_{\perp}}$, where dispersions impacted by compressions and shear differ,
comprises $B_{\rm turb_{\perp}}$.

The isotropic turbulent field can be distinguished by combining unpolarised
synchrotron intensity, beam depolarisation, and Faraday depolarisation (see,
e.g., Ferrière, 2020). The anistropic turbulent field can be difficult to
distinguish from the regular field, so the ordered field energy $B_{\rm
ord_{\perp}}^{2}=B_{\rm reg_{\perp}}^{2}+B_{\rm aniso_{\perp}}^{2}$. This is
determined via intensity and radio, optical, infrared (IR), and submillimeter
polarisation. The Goldreich-Kylafis (Goldreich and Kylafis, 1981) and
longitudinal Zeeman effects provide $B_{\rm reg_{\perp}}$ and $B_{\rm
reg_{\parallel}}$, the latter being parallel to the line of sight.

When interpreting observations, there is an underlying assumption that total
magnetic field and total cosmic ray energy are in equipartition in regions
emitting synchrotron radiation. There are various complications to deriving
the field strength from the synchrotron intensity, which include the increased
energy losses in cosmic rays (CRs), particularly electrons, in higher density
environments, such as starburst regions and massive spiral arms, and the
differences in the energy distribution of CRs produced by various sources,
from which the integrated CR energy is obtained. Notwithstanding such
uncertainties, a survey of 21 nearby galaxies obtained $B_{\rm tot}=17\pm
13\,\mu$G, including $B_{\rm reg}\simeq=5\pm 3\,\mu$G (Fletcher, 2010). Fields
as strong as 50–100 $\mu$G are observed in starburst galaxies and some barred
galaxies (Beck et al, 2005; Heesen et al, 2011; Adebahr et al, 2013). Such
field strengths clearly indicate a magnetic field capable of exerting
significant influences on the evolving dynamics and structure of the galaxies,
including the processes leading to star formation. The relatively high
strength of the turbulent field relative to the regular field is a common
feature across observations on the scale of the Galaxy (see, e.g., Table 3
Beck et al, 2019), who determine $0.04\leq B_{\rm reg}/B_{\rm tot}\leq 0.43$
across 12 galaxies and at a variety of galactic radii.

Even setting aside the uncertainties in the measurements of the magnetic
field, there is no observational diagnostic beyond its overall magnitude to
distinguish the turbulent field that arises through an SSD from that which
results from tangling of the regular field by the turbulent flow. Recent high
resolution observations of extragalactic magnetic fields with SOFIA (Borlaff
et al, 2021) identify polarisation fractions in the dense star forming regions
traced by far-infrared (FIR) anticorrelated with the star formation rate, and
generally lower than the polarisation of 3 and 6 cm radio emission. This is
suggested to imply complex magnetised structures in the star-forming regions
characteristic of SSDs acting there. An accurate model of the turbulent
structure of the magnetic field in galaxies is required, not only to measure
and understand the magnetic field, but also to correctly interpret its effect
on how observables from distant sources reach the observer.

Tsouros et al (2023) use a 3D map of observed dust density distribution in the
local bubble of the Solar neighbourhood (Leike and Enßlin, 2019; Leike et al,
2020) and simulated cosmic ray detection to infer the structure of the local
magnetic field, assuming weak turbulence. Such observational inference
provides a suitable benchmark against which to test numerical models. Insight
derived from numerical modelling into the correlation between the field and
velocities, gas and dust distributions will improve observational inference.
From extragalactic sources Brown et al (2007) use RM to determine the mean
magnetic field in the Milky Way. It aligns along spiral arms and has a
strength up to about $2\mu$G. Within alternating spiral arms it is observed to
switch sign.

Such field reversals are not common in observations of external galaxies.
Determinations of the magnetic field in M31 (Fletcher et al, 2004) and M51
(Fletcher et al, 2011) show the mean and fluctuating fields to be of roughly
equivalent energy, with strength $5\mu$G each across a radial extent between 6
and 14 kpc. Non-axisymmetric configurations are certainly observed, including
misalignment along the spiral arms (Fletcher et al, 2004; Fletcher, 2010;
Fletcher et al, 2011) and field reversals (Brown et al, 2007; Haverkorn,
2015). Pitch angles in M31 steepen from $-17\deg$ at $6<r<10$ kpc to $-8\deg$
at $12<r<14$ kpc. At high redshift, galaxies have already been observed to
have magnetic fields of order $1\mu$G (Bernet et al, 2008).

###  1.2 Governing physical equations

Galactic dynamos are most often studied using the magnetohydrodynamics (MHD)
framework, which describes the interactions of magnetic fields with a plasma
that is treated as a continuum. Further, the plasma is usually assumed to be
resistive, that is, its magnetic diffusivity
$\eta=\left(\mu_{0}\sigma\right)^{-1}\neq 0$, where $\mu_{0}$ is the vacuum
permeability and $\sigma$ the electrical conductivity of the medium. The
justification for these approximations rests on the tiny Larmor radius of
charged particles in the ISM

$r_{\rm L}=\frac{mcv}{qB}=(36\mathrm{km})\left(\frac{m}{m_{\rm
H}}\right)\left(\frac{v}{10\mathrm{kms}^{-1}}\right)\left(\frac{e}{q}\right)\left(\frac{3\,\mu\mathrm{G}}{B}\right)$
(1)

where $m$ and $q$ are the mass and charge of the charged particles, $v$ and
$B$ are the magnitudes of their perpendicular velocity and the ambient
magnetic field, $m_{\rm H}=1.67\times 10^{-24}\mathrm{g}$ and $e=4.8\times
10^{-10}\mathrm{esu}$ are the mass and charge of a proton, and $c=2.99\times
10^{10}\mathrm{cms}^{-1}$ is the speed of light. The details of reconnection
have been argued to depend on kinetic effects (Yamada et al, 2010), but MHD
turbulence may inevitably drive small-scale reconnecting structures regardless
of kinetic scale effects (Eyink et al, 2011; Oishi et al, 2015; Lazarian et
al, 2020). For modelling dynamo behavior the turbulent reconnection process
can be usefully abstracted with Laplacian Ohmic resistivity.

The evolution of magnetic fields in an electrically conducting fluid is
governed by the induction equation

$\frac{\partial{\bm{B}}}{\partial
t}=\nabla\times\left({\bm{u}}\times{\bm{B}}-\nabla\times\eta\nabla\times{\bm{B}}\right),$
(2)

where ${\bm{B}}$ is the magnetic field and ${\bm{u}}$ is the velocity field.
The governing equation for the dynamics of viscous fluids is the Navier–Stokes
equation, which can be written in the following non-conservative form for
rotating, magnetised galaxies in a gravitational potential:

$\rho\frac{D{\bm{u}}}{Dt}=\rho{\bm{g}}-\bm{\nabla}p-2\bm{\Omega}\times\bm{u}+\bm{J}\times\bm{B}+\bm{\nabla}\cdot\left(2\rho\nu{\sf
S}\right),$ (3)

where ${\bm{g}}$ is the gravitational acceleration field, $p$ is the thermal
pressure, $\bm{\Omega}$ is the rotation vector, $\bm{J}$ is the current
density, $\nu$ is the kinematic viscosity, and ${\sf S}$ is the traceless rate
of strain tensor, ${\sf
S}_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right)-\frac{1}{3}\delta_{ij}u_{k,k}$.
Additional terms such as the cosmic ray pressure gradient, coupled to a model
of cosmic ray production, destruction, advection, and energy loss and gain can
also play a role. When the dynamical viscosity, $\rho\nu$, is homogeneous,
then the viscous force can be simplified as

$\nabla\cdot\left(2\rho\nu{\sf
S}\right)=\rho\nu\left(\bm{\nabla}^{2}\bm{u}+\frac{1}{3}\bm{\nabla}\left(\bm{\nabla}\cdot\bm{u}\right)\right).$
(4)

The turbulent component of the velocity field can have its origin in the
cascade of energy from large scales towards small scales thanks to the
nonlinear interactions in the Navier–Stokes equation (Kolmogorov, 1941).
However in galaxies, deviating from the standard picture of forcing at largest
scales driven, e.g., by large-scale shear, energy injection occurs at multiple
scales from distributed stellar feedback in the form of ionising radiation,
stellar winds, and SNe; gravitational instability of the gas disc; and,
especially for young galaxies, gas accretion onto the galaxy Klessen and
Hennebelle (2010). Ionising radiation contributes an order of magnitude more
energy than the SNe (Abbott, 1982), while SNe in turn contribute an order of
magnitude more energy than stellar winds (Shull and Saken, 1995). The smallest
contribution has been estimated to originate from the large-scale shear, an
order of magnitude less than the least powerful stellar source (Abbott, 1982).
Stellar feedback causes mixing of the ISM at much smaller scales than those
associated with systemic large-scale flows and ordered magnetic fields.

Conservation of mass implies

$\frac{\partial\rho}{\partial t}+\bm{\nabla}\cdot(\rho\bm{u})=0,$ (5)

while conservation of thermal energy then implies

$\rho\frac{De}{Dt}=-p\bm{\nabla}\cdot{\bm{u}}+\bm{\nabla}\cdot\bm{q}+2\rho\nu|{\sf
S}|^{2}+\eta\mu_{0}\bm{J}^{2}+\rho\Gamma-\rho^{2}\Lambda,$ (6)

where $\bm{q}$ is the conductive heat flux and the last four terms describe
the viscous and resistive heating and radiative heating and cooling. The
radiative heating and cooling depend sensitively on the temperature,
ionisation state, chemical abundances, metallicity, and dust properties (see
Klessen and Glover 2016 for a comprehensive review of interstellar processes).
The simplest approximation is to treat heating as a constant rate due to far-
ultraviolet ejection of photoelectrons from dust and follow the cooling as a
function of temperature. More complete models include descriptions of
radiation transport and chemistry including both ionisation and molecule
formation to determine both the heating and the cooling. The temperature-
dependent structure of the cooling function allows gas in thermal equilibrium
in an isobaric medium to exist at two temperatures for typical interstellar
pressures, which are described as thermal phases (Field et al, 1969; Wolfire
et al, 1995, 2003). Hot gas at low density takes sufficiently long to cool
that it acts as a third quasi-stable phase (McKee and Ostriker, 1977).

The most common way of closing the equations is to use the ideal gas law

$p=(\gamma-1)\rho e,$ (7)

where $\gamma=c_{\rm P}/c_{\rm V}$, $c_{\rm P}$ is the heat capacity at
constant pressure, and $c_{\rm V}$ that at constant volume, providing the
relation between density and pressure required to solve this set of equations.
However, note that $\gamma$ also varies with the ionisation and chemical state
of the gas, so in the general case this equation is position dependent.

In addition to the stellar feedback, the ISM can undergo instabilities of
various kinds. Especially the magnetorotational instability (MRI),
gravitational instability, and the instabilities caused by cosmic rays are
considered to be important for galactic dynamo action.

The MRI (Balbus and Hawley, 1991; Piontek and Ostriker, 2004, 2005, 2007) can
alter both large- and small-scale dynamics. It has been semi-analytically
argued, however, that this instability would be suppressed by intense SN
activity up to a galactocentric radius of 15 kpc with Milky Way parameters
(Korpi et al, 2010), but it could still affect the dynamics significantly in
the outer regions (Sellwood and Balbus, 1999) explaining the anomalously high
velocity dispersions there (Tamburro et al, 2009; Beck, 2015).

At smaller scales, SN remnant shells undergo fragmentation, e.g., due to the
Richtmeyer-Meshkov instability (Rayleigh–Taylor instability due to
acceleration rather than gravity) or either of the Vishniac instabilities
(Vishniac, 1983, 1994). These instabilities have been argued to be important
in ’breaking up’ the shells of superbubbles and enable them to overcome the
flux-freezing constraint and act favourably for a magneto-buoyantly driven
dynamo (see, e.g., Kulsrud, 2015). The theory for this instability, dubbed as
the spike instability, has been developed in the early linear stage. It seems
unlikely that the assumptions applied in this stage, such as treating the
escape of magnetic flux as a diffusive process, are valid in the non-linear
stage, however.

Another way of triggering magneto-buoyancy in galactic disks is the Parker
instability powered by cosmic rays released in SN explosions (Parker, 1992).
The significant pressure exterted by the cosmic rays enables the frozen-in
magnetic field, even without fragmentation, to buoyantly rise up to the
galactic halo. Meanwhile a portion of the gas will sink along the magnetic
field lines of the so-formed loops, and finally form clumps at the connected
footpoint of the flux tubes. These magnetic loops can be sheared by
differential rotation, twist, and amplify the magnetic field (Hanasz and
Lesch, 1997).

Another instability induced by interstellar cosmic rays is the streaming
instability. According to the analytic theory (Bell, 2004), positively charged
cosmic rays, streaming along a mean magnetic field, e.g., of LSD origin,
induce constant currents along supernova-driven shock jumps in the density of
the cosmic rays and background plasma. In the presence of a magnetic
perturbation perpendicular to the shock, the action of the Lorentz force
causes the background plasma to move transversely, while the cosmic ray plasma
is much less affected. If the magnetic perturbations are helical, these
motions can amplify the magnetic field by orders of magnitude (see, e.g.,
Riquelme and Spitkovsky, 2009). Even though the scale of the plasma is much
smaller than the Larmor radius of the background magnetic fields, this
instability has recently been argued to cause significant effect on the scale
of the small-scale magnetic field and interstellar turbulence (Riquelme and
Spitkovsky, 2009).

Gas-rich galactic disks characteristic of high-redshift galaxies have strong
gravitational instabilities of the type described by Toomre (1964) and
Goldreich and Lynden-Bell (1965). These drive turbulence from large scales
independent of stellar feedback (e.g. Krumholz et al, 2018) that drive prompt
SSD action early in galactic history (see Sect. 3.5 below). Furthermore,
Klessen and Hennebelle (2010) showed that accretion flows, whether on to
galaxies or individual star-forming cores, can drive substantial amounts of
turbulence. These flows have been shown by numerical models to be dynamo
active in the case of first star formation (e.g., Sur et al, 2010) and during
galaxy formation (Pfrommer et al, 2022, e.g.). For further review, see Sects.
3.2 and 3.4 below.

###  1.3 Basics and challenges of dynamo theory

In order to model dynamos, both the velocity and magnetic field can naturally
be decomposed into coherent and fluctuating parts, so that we can write
${\bm{u}}=\overline{{\bm{u}}}+{\bm{u}^{\prime}}$ and
${\bm{B}}=\overline{{\bm{B}}}+{\bm{b}^{\prime}}$, and follow their evolution
separately. Ideally, the averaging should be over many independent ensembles
of the turbulent velocity field, but in practice, spatial averages are used.
If the averaging operation satisfies the Reynolds rules (see, e.g.,
Brandenburg and Subramanian, 2005), then no strict requirement for clear
separation between the average and fluctuating scales exists, and the
decomposition of the equations into these two parts is generally valid.

As the largest scale galactic magnetic fields are predominantly axisymmetric,
azimuthal averaging obeying the Reynolds rules has traditionally been used.
This choice can be criticised on at least two grounds. Firstly, many magnetic
field configurations in galaxies are not axisymmetric, so a Fourier filtering
based decomposition would be a more general choice for these. Secondly,
observational diagnostics are hardly ever derived basing on such averaging.
Gaussian filtering was successfully applied to decompose into the mean and
fluctuating fields by Gent et al (2013b), which would be a more consistent
choice when comparing with observations (see, e.g., the discussion in Hollins
et al, 2022). The major difficulty with these alternative averaging techniques
is that neither of them obeys the Reynolds rules, although this is no
hindrance to deriving alternative, mathematically sound formalisms for the
mean quantities (e.g. Germano, 1992). The bulk of analytic and mean-field
galactic dynamo models continue to use azimuthal averaging and, hence, we
employ it in this review of the theoretical concepts.

####  1.3.1 Large-scale dynamo instability

After taking the azimuthal average of the induction equation, we obtain an
evolution equation for the mean field

$\frac{\partial\overline{{\bm{B}}}}{\partial
t}=\bm{\nabla}\times\left(\overline{{\bm{u}}}\times\overline{{\bm{B}}}+\mbox{\boldmath${\cal
E}$}{}-\eta\bm{\nabla}\times\overline{{\bm{B}}}\right),$ (8)

where $\mbox{\boldmath${\cal
E}$}{}=\overline{{\bm{u}^{\prime}}+{\bm{b}^{\prime}}}$ is the turbulent
electromotive force, which depends on the correlations of the velocity and
magnetic fluctuations. The evolution equation for the magnetic fluctuations is
obtained by subtracting the mean-field equation from the original induction
equation:

$\frac{\partial{\bm{b}^{\prime}}}{\partial
t}=\bm{\nabla}\times\left(\overline{{\bm{u}}}\times{\bm{b}^{\prime}}+{\bm{u}^{\prime}}\times\overline{{\bm{B}}}+{\bm{u}^{\prime}}\times{\bm{b}^{\prime}}-\mbox{\boldmath${\cal
E}$}{}-\eta\bm{\nabla}\times{\bm{b}^{\prime}}\right).$ (9)

Two regimes are usually distinguished, where the Lorentz force due to the mean
magnetic field (1) is still weak, so the magnetic field will not back react on
the flow, called the kinematic regime; (2) is significant, so the magnetic
field influences the flow, called the nonlinear regime. Further distinction
must be made in SSD-active systems, where also magnetic background turbulence
will exist; under these circumstances, there exist magnetic fluctuations,
${\bm{b}}_{0}^{{}^{\prime}}\neq 0$, unrelated to the presence of
$\overline{{\bm{B}}}$, and the turbulent transport coefficients can depend on
these background flunctuations (see, e.g., Rheinhardt and Brandenburg, 2010).

In the kinematic regime, further restricting to the case with
${\bm{b}}_{0}^{{}^{\prime}}=0$, the presence of a mean magnetic field will
induce a fluctuating magnetic field through the term
${\bm{u}^{\prime}}\times\overline{{\bm{B}}}$. Therefore, even if
${\bm{u}^{\prime}}$ and ${\bm{b}^{\prime}}$ were uncorrelated without the mean
magnetic field, in its presence they will be correlated so that
${\bm{u}^{\prime}}\times{\bm{b}^{\prime}}$ is no longer zero, and in addition
linearly correlated with $\overline{{\bm{B}}}$. If $\overline{{\bm{B}}}$
varies over scales much larger than those of the turbulent fluctuations, then
one can expand ${\cal E}$ in a Taylor series expansion of the form

$\mbox{${\cal
E}_{i}$}{}=a_{ij}\overline{B}_{j}+b_{ijk}\frac{\partial\overline{B}_{j}}{\partial
x_{k}}+...$ (10)

and truncate after the first-order spatial derivatives of
$\overline{{\bm{B}}}$. Here, $a_{ij}$ and $b_{ijk}$ are tensors that describe
the influence of turbulence on the evolution of the mean magnetic field,
commonly described as turbulent transport coefficients. Particularly, the
symmetric part of the $\bm{a}$ tensor, usually denoted with $\bm{\alpha}$,
describes the collective inductive effect resulting from turbulent motions,
and the symmetric part of the rank two tensor acting upon
${\bm{\nabla}}\times\overline{{\bm{B}}}$, usually denoted with $\bm{\beta}$,
describes the enhanced magnetic diffusivity resulting from them. In the
kinematic regime, the transport coefficients may depend on the key
hydrodynamic and thermodynamic quantities, but not on $\overline{{\bm{B}}}$
itself. With this simplification, only the evolution of the magnetic fields at
the largest scales can be captured. Simplifying even further, as done in the
oldest and most traditional first-order smoothing approximation (FOSA) or
second-order correlation approximation (SOCA) (see, e.g. Brandenburg and
Subramanian, 2005) the term
$\bm{\nabla}\times({\bm{u}^{\prime}}\times{\bm{b}^{\prime}}-\overline{{\bm{u}^{\prime}}\times{\bm{b}^{\prime}}})$
is neglected in Eq. (9) for the fluctuating field. This is valid only when the
magnetic Reynolds number, ${\rm Re}_{\rm M}\ll 1$, or when the Strouhal number

$\mathrm{St}=\frac{\tau u}{l}\ll 1,$ (11)

where $\tau$ is the correlation time of the turbulence, and $u$ and $l$ are
typical velocity and length scales, respectively. Under one of these
assumptions, the equation for the fluctuating field, Eq. (9), simplifies to
the extent that it can be solved analytically for a given ${\bm{u}^{\prime}}$,
after which ${\cal E}$ can be calculated with relative ease (Steenbeck et al,
1966).

Particularly, in the case of isotropic and homogeneous turbulence, the
turbulent transport coefficients then reduce to scalars,

$\displaystyle\alpha\approx-\frac{1}{3}\tau\overline{{\bm{u}^{\prime}}\cdot\omega^{\prime}},$
(12)
$\displaystyle\beta\approx\frac{1}{3}\tau\overline{{\bm{u}^{\prime}}^{2}},$
(13)

where $\omega^{\prime}$ is the fluctuating vorticity, and
${\bm{u}^{\prime}}\cdot\omega^{\prime}$ is the small-scale kinetic helicity of
the turbulent flow. It can be seen that $\alpha$ relates to the kinetic
helicity of turbulence, and $\beta$ to the intensity of the turbulent motions.
A more generic form of the turbulent emf can be obtained by relaxing isotropy
and homogeneity, which certainly can be appropriate in galaxies, where
rotation and large-scale shear are present:

$\mbox{\boldmath${\cal
E}$}{}={\bm{\alpha}}\cdot\overline{{\bm{B}}}+{\bm{\gamma}}\times\overline{{\bm{B}}}-{\bm{\beta}}\cdot(\mbox{\boldmath$\nabla$}{}\times\overline{{\bm{B}}})-{\bm{\delta}}\times(\mbox{\boldmath$\nabla$}{}\times\overline{{\bm{B}}})-{\bm{\kappa}}\cdot(\mbox{\boldmath$\nabla$}{}\overline{{\bm{B}}})^{(s)},$
(14)

$\alpha$ and $\beta$ tensors now acquire off-diagonal elements. Completely new
effects arise, as discussed, e.g., by Krause and Rädler (1980): $\gamma$
describes turbulent pumping of the mean magnetic field (e.g., Ossendrijver et
al, 2002); $\delta$ encompasses new dynamo effects arising from current and
rotation ($\Omega\times\bm{\overline{J}}$ or Rädler effect Rädler, 1969) or
shear (shear-current effect Rogachevskii and Kleeorin, 2004);
$(\mbox{\boldmath$\nabla$}{}\overline{{\bm{B}}})^{(s)}$ is the symmetric part
of the derivative tensor of $\overline{{\bm{B}}}$; and $\kappa$ is a rank-
three tensor, which has also been shown to be important in the dynamo process
by Warnecke et al (2021), although it is neither fully diffusive nor
inductive.

In the nonlinear regime, still restricting to the case
${\bm{b}}_{0}^{{}^{\prime}}=0$, the turbulent transport coefficients may
depend on both $\overline{{\bm{B}}}$ and on ${\bm{b}^{\prime}}$. The outcome
of taking into account these nonlinearities is strongly dependent on the
closure adopted for Eq. (9) (see, e.g., Shukurov and Subramanian, 2021). Until
the turn of the millennium, most nonlinear galactic dynamo models adopted a
heuristic closure, where mainly the effect of $\overline{{\bm{B}}}$ on the
$\alpha$ effect is postulated to be important, so that the inductive effect is
suppressed when the mean field is approaching equipartition with the
turbulence:

$\alpha=\frac{\alpha_{0}}{1+|\overline{{\bm{B}}}|/|B_{\rm eq}|}.$ (15)

Here, $\alpha_{0}$ is the kinematic value of the $\alpha$ effect and $B_{\rm
eq}\equiv\rho\overline{{\bm{u}^{\prime}}^{2}}$ is the field strength at
equipartition with the turbulent kinetic energy. Evidence from numerical
simulations in fully periodic domains in the beginning of the 1990’s (e.g.
Cattaneo and Vainshtein, 1991; Vainshtein and Cattaneo, 1992) forced the
dynamo community to seriously consider other, more complicated and complete,
closures, some of which had already been formulated years ago, (e.g. Pouquet
et al, 1976; Vainshtein and Kichatinov, 1983) and derive better nonlinear
models based on them. Numerical simulations of forced turbulence, also
exhibiting exponential amplification of magnetic fluctuations through the SSD
(case of ${\bm{b}}_{0}^{{}^{\prime}}\neq 0$), showed that the $\alpha$ effect
is not quenched according to Eq. (15), but in a catastrophic fashion, where
also the magnetic Reynolds number appears in the denominator. Practically,
this was interpreted to imply that galactic (or any astrophysical)
$\alpha$-effect based dynamos should not exist.

However, this quenching behavior was later understood not to result from the
SSD-generated fluctuations, but to be a consequence of magnetic helicity
conservation (summarised e.g. by Brandenburg and Subramanian, 2005): under the
mean-field approximation, a dynamo produces helicities of equal magnitude but
of opposite signs to satisfy helicity conservation. If no helicity fluxes
escape the dynamo active region, then the small-scale magnetic helicity will
quench the turbulent electromotive force parallel to the mean field, leading
to the catastrophic quenching of the LSD. This is because the $\alpha$ effect
acquires a magnetic modification, so that

$\alpha\equiv\alpha_{\rm K}+\alpha_{\rm m}.$ (16)

One path into deriving the magnetic contribution to the $\alpha$ effect is the
so-called minimal $\tau$ approximation, where the triple correlations
appearing in the emf closure problem are modelled as a relexation term (for
details, see, e.g., Brandenburg and Subramanian, 2005). The turbulent
transport coefficients then read for incompressible flows

$\displaystyle\alpha=-\frac{\tau}{3}\left(\overline{\bm{\omega}^{\prime}\cdot{\bm{u}^{\prime}}}-\overline{\bm{j}^{\prime}\cdot{\bm{b}^{\prime}}}/\overline{\rho}\right),$
(17) $\displaystyle\beta=\frac{\tau}{3}\overline{{\bm{u}^{\prime}}^{2}},$ (18)

where
${\bm{j}}^{{}^{\prime}}={\bm{\nabla}}\times{\bm{b}}^{{}^{\prime}}/\mu_{0}$ is
the fluctuating current density, and $\overline{\rho}$ is the mean density.
According to the magnetic helicity conservation law, the magnetic part has to
satisfy (see, e.g., Brandenburg and Subramanian, 2005) leading to an evolution
equation

$\frac{\partial\alpha_{\rm m}}{\partial t}=-\eta_{\rm
t}k_{0}^{2}\left(\frac{2\mbox{\boldmath${\cal
E}$}{}\cdot\overline{{\bm{B}}}+\nabla\cdot\overline{{\bm{F}}}}{B_{\rm
eq}^{2}}+\frac{2\alpha_{\rm m}}{\widetilde{{\rm Re}_{\rm M}}}\right).$ (19)

Here $\eta_{\rm t}=\eta+\beta$ is the sum of microscopic resistivity, $\eta$,
and the enhanced turbulent resistivity, $\beta$, $\overline{\bm{F}}$ is the
large-scale helicity flux of the fluctuating field, $k_{0}$ is the wavenumber
corresponding to the integral scale of turbulence, and $\widetilde{{\rm
Re}_{\rm M}}\equiv\eta_{\rm t}/\eta$. The most advanced mean-field (MF) models
for galaxies solve this equation together with the induction equation. With
some plausible estimates from observations or theory for the coefficients of
the models, galactic MF dynamos remain useful tools (see, e.g., Shukurov and
Subramanian, 2021), but verifying their applicability rigorously remains an
open question (see the discussion in Sect. 3.4.1).

Figure 1: Schematic illustration of different types of SSD spectra. The solid
line represents the kinetic energy spectrum, with a cut-off scale $k_{\nu}$,
$\nu$ being the kinematic viscosity. The dashed lines represent magnetic
energy spectra in the kinematic stage, with cut-off scales $k_{\eta}$, $\eta$
being the magnetic resistivity. Two different scenarios, differentiated by the
magnetic Prandtl number, ${\rm Pr}_{\rm M}=\nu/\eta$, are shown, corresponding
for low (such as the Sun) and high (such as the ISM). The red long-dashed line
illustrates the contribution of a LSD at the low wavenumber range. Image
adapted from Rempel et al (2023).

####  1.3.2 Small-scale dynamo instability

The ISM of galaxies is, in general, a high magnetic Prandtl number fluid,
meaning that the ratio of magnetic to fluid Reynolds numbers

${\rm Pr}_{\rm M}\equiv\frac{\nu}{\eta}=\frac{{\rm Re}_{\rm M}}{\rm Re},$ (20)

is large; Brandenburg and Subramanian (2005) estimate values of ${\rm Pr}_{\rm
M}\sim\mathcal{O}(10^{11})$. The resistivity, $\eta$, is then much smaller
than the fluid viscosity, $\nu$, meaning that the dissipation scale of the
fluid motions far exceeds that of the magnetic field. This implies that the
magnetic fields can dissipate into much smaller-scale structures than the
fluid motions themselves, and hence a typical magnetic field structure sees
the velocity field as a smooth, larger-scale shearing flow acting on it,
(compare spectra plotted with solid and dashed lines plotted in Fig. 1).

In such a setting, a rapid amplification of magnetic fluctuations by the
smooth fluid motions can occur, independent of the flow possessing helicity or
not. This dynamo instability is called the SSD or fluctuation dynamo. The
scenario was first discussed by Batchelor (1950), and the first rigorous
treatment of the non-helical case was presented in Kazantsev (1968). The
amplification occurs at the timescale of the turbulent eddy turnover time,
which is far shorter than the timescale required for the amplification of the
large-scale field by an LSD. In the kinematic growth phase, the structures
generated have a thickness around the resistive scale, but are curved up to
the scale of the turbulent eddies. The only requirement for this instability
to set in is that ${\rm Re}_{\rm M}$ must exceed a critical value,
analytically derived to be 30–60, a condition that is fulfilled from the
integral to the dissipation scale in galaxies (see e.g. Subramanian, 1998).

Presumably, the setting in which the SSD starts to operate in galaxies is one
where the LSD has not yet produced a significant large-scale field, but
battery mechanisms and plasma instabilities can produce seed fields at small
scales in the regime of dissipation scales rather than within the inertial
range (see Sects. 1.4 and 3.1). Therefore, the SSD excitation likely starts in
a diffusion-dominated setting, in contrast to the other possible, but
unlikely, scenario of SSD starting with larger-scale tangled fields, in the
diffusion-free regime (Rincon, 2019). In the kinematic stage, when the
magnetic fluctuations are still weak, the peak in the magnetic power spectrum
is at the resistive scale, but all scales grow with the same growth rate. At
scales larger than the resistive scale, the magnetic power spectrum is
predicted to follow a positive power law of $k^{3/2}$ (Kazantsev, 1968),
called the Kazantsev spectrum. At smaller scales, the spectrum falls off
steeply, following the Macdonald function (for details, see e.g. Brandenburg
and Subramanian, 2005; Rincon, 2019). For a schematic representation of the
magnetic spectrum, see Fig. 1 dashed line.

How the SSD nonlinearly saturates remains an open problem. By extending the
Kazantsev model, it has been proposed that the SSD grows magnetic fluctuations
at the resistive scale up to and exceeding equipartition with the kinetic
energy of the turbulence, but that the generated field would be concentrated
into resistive-scale ropes, hence not being volume filling, and therefore
being energetically insignificant (e.g. Subramanian, 1998). Whether
reconnection plays a role in the nonlinear regime of SSD is a relatively fresh
research question arising from such numerical research that is discussed in
this review.

In galaxies that have already excited LSD, another source of magnetic
fluctuations comes through tangling of the coherent component by the action of
the turbulent eddies. Currently, the fraction of fluctuations caused by the
SSD versus tangling is not known observationally, and calls for numerical
approaches. Ideally, models where both SSD and LSD instabilities and their
non-linear interactions and saturation regimes are realistically captured are
required to study the relation of the regular and fluctuating components. Such
modelling requires large scale separation ratios, and thereby high resolution,
and the models need to be integrated over timescales of gigayears. Currently,
numerical models reach only modest magnetic Prandtl number regimes, while
extending to the extreme regimes in galaxies is yet impossible and likely to
remain that way in the near future. Even when working in parameter regimes
far-removed from the ISM, numerical simulations with modest Prandtl numbers
have provided some answers, and here we review the current state-of-the-art.

###  1.4 Seed fields

Both SSD and LSD need initial magnetic fields to grow, as their evolution
equations are homogeneous in $\overline{{\bm{B}}}$ and ${\bm{b}^{\prime}}$,
respectively. There have been a number of mechanisms proposed to generate such
seeds. These have included seeds of cosmological origin (see, e.g., the review
by Durrer and Neronov, 2013); battery mechanisms (see, e.g. Subramanian et al,
1994); and stars picking up the battery fields during star formation,
amplifying them through a stellar dynamo, and then ejecting the fields as
seeds in SN explosions (Ruzmaikin et al, 1988; Martin-Alvarez et al, 2021).
However, recent work has demonstrated conclusively that kinetic plasma
instabilities can generate magnetic fields from field-free initial conditions
that in the presence of turbulence grow rapidly with an SSD, as we describe
below in Sect. 3.1. These novel findings alleviate the previous concerns of
the origin of strong enough seed fields to be grown by LSD alone in galaxies.
Therefore, we leave further detailed discussion of the origin of seed fields
out of the current review.

##  2 Numerical approaches

###  2.1 Mean-field modelling

The simplest numerical models use the approximations discussed in Sect 1.3.1
to study the evolution of the large-scale component of the galactic magnetic
field. Furthermore, as the rotation laws of galaxies are well-constrained
observationally, the models most often work in the limit where the large-scale
velocity field is steady instead of being evolved according to the mean-field
Navier–Stokes equation. As such dynamo models do not solve for the full
dynamical problem, but only for the evolution of the mean magnetic field (Eq.
8), they need ad hoc parameterisations of the turbulent transport
coefficients. A comprehensive review of the MF approach and the results with
given differential rotation profiles has been recently presented in Shukurov
and Subramanian (2021). We summarise their theoretical discussion in the
following, while for details and the results in the limit of given
differential rotation laws we refer to their Chapters 11 and 12.

The dominant turbulent effects in galaxies are assumed to be the $\alpha$
effect and the enhanced turbulent diffusion, while turbulent pumping and the
more exotic dynamo effects discussed in Sect. 1.3.1, are usually regarded to
be sub-dominant. Furthermore, the scalar expressions for isotropic and
homogeneous turbulence, Eqs. (12) and (13), are often assumed to be valid. As
discussed earlier, it is observationally well-established that strong magnetic
fluctuations dominate over the mean component in galaxies. The effect of the
fluctuations, most likely implying strong pre-existing magnetic background
turbulence ($b_{0}^{{}^{\prime}}\neq$ 0), are currently not taken into account
through the turbulent transport coefficients as such, but rather through their
non-linear saturation formulae, Eqs. (16)–(19).

Observations allow for estimates of the magnitudes of $\alpha$ and $\beta$,
and also that of the amplifying effect of differential rotation, making
galactic dynamos better constrained than many other astrophysical dynamos,
such as the solar one. Nevertheless, confirmations from numerical models, not
working under the MF model assumptions, are required to test the validity of
the various assumptions used, and to fill in details that the observations can
not reveal due to the limited accuracy and projection effects. Such attempts,
tools used, and the results obtained are discussed in Sect. 4.1.

Cylindrical coordinates ($r,\phi,z$), with the $z$-axis aligned with the
rotation vector, and $(r,\phi)$ being coordinates in the horizontal plane, are
the most often used coordinates to simulate the MF equation. Under the
assumption of axisymmetry, $\partial/\partial\phi$ vanishes. If we can further
assume that the disc is very thin, that is, $\epsilon=h_{0}/R\ll 1$, where
$h_{0}$ is the scale height and $R$ is the galactocentric radius, one may also
neglect all the terms with radial derivatives, so the evolution equation for
$B_{z}$ decouples from that of the horizontal components. One is then left
with the problem of solving the evolution equations for the radial and
azimuthal components. In this approximation, the MF equations contain only
$z$-derivatives, but the turbulent transport coefficients and the mean field
itself can depend parametrically on the radius, and are thus local. The non-
dimensional MF equations contain two control parameters for the LSD, called
dynamo numbers:

$R_{\Omega}=\frac{S_{0}h_{0}^{2}}{\beta},{\ \
}R_{\alpha}=\frac{\alpha_{0}h_{0}}{\beta},$ (21)

where ’0’ subscripts now refer to characteristic values determined, e.g., from
observations, and $S$ describes the shear, for flat galactic rotation curves
$S=-V_{0}/R$ being a good approximation, with $V_{0}$ being the roughly
constant rotational velocity. Their product $D=R_{\alpha}R_{\Omega}$
characterises the overall dynamo efficiency. If we adopt the so-called
$\alpha\Omega$ assumption, where the azimuthal magnetic field is solely
generated by differential rotation, neglecting a term arising from the
$\alpha$ effect, that can be estimated to be much smaller, then the parameter
$D$ is the only non-dimensional control parameter describing the full problem.
Using qualitative arguments, that is, not solving any dynamo equations yet,
the growth rate $\lambda$ is expected to scale as

$\lambda^{-1}\approx\frac{h_{0}^{2}}{\beta}\mbox{{}\rm\sf\leavevmode\hbox{}\hbox{}
D \/}^{-1/2}.
\end{equation}Fornumericalsolutions,vacuumboundaryconditionsareoftenused,thejustificationcomingfromtheestimatethatthemagneticdiffusivitybeyondthediscisonetotwoordersofmagnitudehigherthaninit,andthelarge-
scalemagneticfieldislargelyconfinedtothethindisc.\par
Theassumptionofaxisymmetryisastrongrestriction,asmanygalaxieshavebeenobservedtohostnon-
axisymmetricmagneticfieldconfigurations(Sect.~{}\ref{sect:obs}).Theapproachoutlinedabovecan,however,berelativelyeasilyappliedtoanon-
axisymmetricsystemaswell.NowtheparametricdependenceofthemeanmagneticfieldnotonlyincludesR,butalsotheazimuthalangle\phi.Theturbulenttransportcoefficientscanbeeitheraxi-
ornon-
axisymmetric.Thereisonemajordifferencetobeaccountedfor:differentialrotationactsdifferentlyonaxi-
versusnon-axisymmetricsolutions,efficientlydestroyingthelatter,whilethe\alpha
effectdoesnotmakesuchastrongdistinction.Hence,the\alpha\Omega
approximationisnotapplicableinthisregime,andbothdynamonumbersmustbekeptintheproblem.Again,forverythindiscs,anasymptoticsolutioncanbeconstructed,whereradialderivativesvanish,andtheevolutionequationforthez-
componentdecouples.Qualitativeconsiderationsshowthatthemaincharacteristicsofthedynamosolutions,suchasthegrowthrate,shouldremaininvariantwithrespecttotheaxisymmetricsolutions.Onlynumericalstudiesrevealdifferences,\cite[citep]{(see,
e.g., \@@bibref{AuthorsPhrase1Year}{SS22}{\@@citephrase{, }}{})}.\par
Thefundamentalprocessthroughwhichgalacticdynamosareexpectedtosaturateismagnetichelicityconservation,leadingtothedynamical\alpha
effect,Eq.~{}\eqref{eq:alpdyn}.Galacticdiscs,althoughoftentreatedinMFmodelsasthindiscs,interactstronglywiththeirenvironments,withstellarfeedbackdrivinggalacticfountainsandwindsandmaterialaccretingfromthesurroundingcircumgalacticmedium\cite[citep]{(e.g.
\@@bibref{AuthorsPhrase1Year}{tumlinson2017}{\@@citephrase{,
}}{})}.Undersuchconditions,magneticfieldisexpectedtobetransportedoutwardswhentiedtothehot,ionisedgasinfountainsandwinds,withsomeofitfallingbackwithcooleranddenserstructuresinthefountainsandduringaccretion.Also,aturbulentpumpingeffectmaybeimportanthereintransportingmeanfieldinwardsoroutwards.Alltheseprocesseswillresultinre-
distributionofmagnetichelicityatdifferentscales,andhelicityfluxesbetweendifferentregions.Parameterisingtheseeffectsintermsofhelicityfluxesfromtheexistingnumericalmodelsisatremendouschallenge.WhatiscurrentlyassumedintheMFmodelsisdescribedinSect.~{}\ref{MF:results},andwhatisknownfromnumericalsimulationsisdiscussedinSect.~{}\ref{HelicityFluxes:results}.\par$
(22)

###  2.2 Direct numerical simulation versus large eddy schemes

Nonlinear numerical simulation of dynamo action requires that resistivity be
added, either explicitly or by discretisation error, to the equations of ideal
MHD. The discretisation error present in any numerical method does provide a
nonlinear resistivity that can enable or even suppress dynamo activity.
Theoretically, an LSD can be excited whenever ${\rm Re}_{\rm M}>1$, while an
SSD is known to require higher threshold values of order of tens to hundreds
(see Sects. 1.3.1 and 1.3.2). Therefore, an LSD should be easier to model
numerically than an SSD, but an LSD requires many more pre-requisites from
turbulence, such as the standard expectation of small-scale flows being
helical, to be satisfied than an SSD.

Including explicit diffusivities in numerical models at the extremely small
values characteristic of the ISM, while retaining their correct ratio as
reflected by the extremely large ${\rm Pr}_{\rm M}$ in the ISM, is currently
an unreachable challenge. Often models including explicit diffusivity are
called direct numerical simulations (DNS), but this denotion is not strictly
valid due to the orders of magnitude elevated diffusivities required to exceed
numerical diffusivity. A more proper term of DNS-like models has recently been
proposed (Rempel et al, 2023). To be useful, such implementations must be run
at high enough resolution that the numerical resistivity does not dominate
over the physical values chosen. A parameter study to find the resolution
where the physical resistivity starts to determine the solution must be done
to demonstrate that numerical resistivity is not dominant.

Using explicit diffusivities has the disadvantage that it provides constant
diffusion everywhere, while the necessity for using diffusion is the largest
at the grid scale to ensure numerical stability and allow magnetic field lines
to reconnect, while no such need exists for the resolved scales. Hence,
implicit large-eddy simulations (ILES) are often used, where no explicit
physical diffusion terms are included, but numerical diffusion remains active.
In astrophysical modeling, the most usual sources of numerical diffusion are
provided by truncation error in Riemann solvers (Balsara et al, 2004, e.g.) or
hyperdiffusivities introduced for numerical stability (Brandenburg and Dobler,
2002a; Brandenburg et al, 2020). Grete et al (2023) recently showed that such
ILES schemes produce results that compare surprisingly well to the results of
DNS-like models with well-resolved dissipation scales. The disadvantage of
these sorts of ILES schemes is that control parameters such as $\rm Re$, ${\rm
Re}_{\rm M}$, and ${\rm Pr}_{\rm M}$ become ill-defined. To counter-act this
shortcoming, measurements have been performed of numerical resistivity
(Rembiasz et al, 2017; McKee et al, 2020; Bian et al, 2021) and of numerical
viscosity (Schranner et al, 2015; Obergaulinger and Aloy, 2020; McKee et al,
2020; Bian et al, 2021) by comparing results of models without physical
diffusivity to analytical or numerical solutions including physical
diffusivity. Despite these efforts, the ILES approach remains unsatisfactory,
as the diffusivity has a steep, unphysical scale dependence, with all
diffusivity focused at scales close to the grid.

Explicit large-eddy simulations (ELES), where the solved equations are
filtered at an intermediate scale, and the smaller-scale effects are described
through a (ideally) physically motivated sub-grid-scale (SGS) model, have been
developed for idealized MHD turbulence (Smagorinsky, 1963; Grete et al, 2017).
Although SGS models to model processes as SNe and star formation are popular,
numerical SGS models attempting to describe the dynamo processes in the ISM
remain scarce. One recent example is the model by Liu et al (2022) that
includes an unresolved turbulent dynamo. The most complicated case for the
ELES schemes is presented by an LSD, in which small-scale, unresolved, effects
on the large-scale dynamics, usually referred to as back-scatter in the LES
framework, should be captured by the SGS model. Steps in the direction of a
true ELES MHD model for supersonic turbulence have been taken by Grete et al
(2015) and Grete et al (2017), but such an SGS model has not yet been
incorporated into galactic dynamo simulations.

###  2.3 Spectral methods

The first published numerical simulation of a turbulent dynamo was performed
with a spectral method (Meneguzzi et al, 1981). Spectral methods have the
highest accuracy derivative per grid point in smooth flows (e.g. Maron et al,
2008), but struggle to simulate the super-Alfvenic flows characteristic of the
interstellar medium because of Gibbs ringing at shocks, so they have not
featured prominently in studies of galactic dynamos.

###  2.4 Grid methods

####  2.4.1 Uniform grid methods

Discretisation of the MHD equations onto a uniform Cartesian grid provides the
most straightforward numerical approach. For the dynamo problem, a prime
advantage of uniform grid methods available in codes such as the Pencil Code
(Brandenburg and Dobler, 2002b; Brandenburg et al, 2020), Athena (Stone et al,
2008), Athena++ (Stone et al, 2020), ZEUS-MP (Hayes et al, 2006), NIRVANA-III
(Ziegler, 2004, 2011), or Piernik (Hanasz et al, 2010a, b) is the uniform
numerical dissipation properties across the computational domain. This
advantage is unfortunately lost in the non-Cartesian geometries available in
these codes, where different grid zones have different sizes or aspect ratios.
Kritsuk et al (2011) compares single-grid algorithms for isothermal,
magnetised, supersonic turbulence showing that algorithmic numerical
diffusivity and resolution can be traded against each other to get equally
good solutions.

####  2.4.2 Hierarchical grid methods

Adaptive mesh refinement (AMR) methods focus resolution in subgrids into
regions of interest, typically where strong gradients occur. Kritsuk et al
(2006) demonstrated in isothermal, supersonic, hydrodynamic turbulence without
an explicit viscosity using Enzo (Collins et al, 2011) that if the driving and
dissipation scales are well separated, and the base grid has sufficient
resolution to resolve the larger scales of the turbulent cascade, then AMR of
dissipation regions can reduce the cost of the simulation without changing the
statistical properties of the turbulence. Li et al (2012) similarly
demonstrated using the ORION2 AMR MHD code that the statistics of magnetised
turbulence can equally well be modeled with AMR methods, though with an
emphasis on the properties of the highest density, highest resolution regions
at the cost of the low density regions.

AMR methods have been applied to SSDs in dwarf galaxies by Rieder and Teyssier
(2016) and Rieder and Teyssier (2017), and placed in a cosmological context in
Rieder and Teyssier (2017). Martin-Alvarez et al (2022) argue that to study
SSDs in cosmological simulations of galaxies, refinement down to a uniform
grid encompassing the galactic disc yields better results than a refinement
strategy strictly tied to density gradients that results in lower resolution
in low-density warm regions.

###  2.5 Particle methods

There are multiple methods tied to particles rather than a structured grid
that have been used for dynamo simulations, including smoothed particle
hydrodynamics (SPH), meshless finite mesh (MFM) and meshless finite volume
(MFV) methods, and Riemann solvers on a Voronoi mesh.

Dolag and Stasyszyn (2009) implements MHD in the GADGET SPH code (Springel,
2005) using a Powell et al (1999) eight-wave divergence cleaning algorithm,
which diffuses magnetic divergence across the grid to constrain it at a low
level. This was used by subsequent publications for study of galactic field
evolution in isolation or cosmological context (e.g. Beck et al, 2012, 2013;
Steinwandel et al, 2019, 2020). An Euler potential approach Price and Bate
(2007) is implemented by Kotarba et al (2009) into the VINE SPH code
(Wetzstein et al, 2009), providing excellent conservation of field divergence
but constraining magnetic helicity unphysically. Wissing and Shen (2020)
implements MHD within the geometric density averaged force algorithm used by
the Gasoline2 SPH code (Wadsley et al, 2017). This was used by Wissing and
Shen (2022) to model dynamo activity in an isolated disc galaxy. However,
these methods remain problematic because of the lack of zero-order convergence
(Morris, 1996; Dilts, 1999; Read et al, 2010) and the large kernel sizes
needed to run modern versions stably.

Pakmor and Springel (2013) describes the implementation of MHD in the Voronoi
mesh code AREPO (Springel, 2010) using an approximate Riemann solver to derive
fluxes across the mesh edges. They used a Powell et al (1999) scheme to
maintain field divergence small. This was used by Pakmor et al (2014) and
Pakmor et al (2017) to simulate galactic dynamos in cosmological context. Mocz
et al (2016) showed how to use constrained transport (Evans and Hawley, 1988)
with AREPO, to maintain zero divergence to machine precision even on the
unstructured mesh.

An alternative approach to modelling fluid flow with particles is offered by
weighted meshless schemes (Gaburov and Nitadori, 2011), as implemented in the
MFV and MFM methods used in GIZMO (Hopkins, 2015). These methods use values at
a particle and its neighbors to fit a polynomial to the solution and its
derivatives, rather than averaging over the values and then taking derivatives
using differences of the averages as SPH does. The difference between the two
GIZMO methods lies in the decision of whether the Riemann problem at zone
faces should be solved in the frame of motion of the fluid at the face (MFM)
or at the mean cell velocity (MFV). MHD is implemented with a Godunov HLLD
solver and a Dedner et al (2002) magnetic divergence damping scheme.

###  2.6 Hybrid and kinetic methods

Dynamos in collisionless plasmas, such as those found in the intracluster
medium surrounding galaxies or in the hot early universe, can not be properly
modeled with the MHD approximation because of the presence of fast kinetic
instabilities. Rincon et al (2016) used a hybrid model that treats protons as
kinetic particles and electrons as an isothermal fluid to study this problem.
Pusztai et al (2020) used a fully kinetic model implemented with a discrete
Galerkin solution of the kinetic Boltzmann equation limited to subviscous
scales for computational tractability. Sironi et al (2023) and Zhou et al
(2023) describe fluctuation dynamos modeled with a particle-in-cell treatment
of an electron-positron pair plasma. These methods are vital for understanding
the origin of seed fields, but are less relevant for most flows within the
bulk ISM of galaxies. Extending them to include electron-ion mass ratios and
greater dynamic range in scale will be essential for making contact with the
physical parameter ranges occurring during galaxy formation.

##  3 Results from numerical approaches

Simulations of dynamo action in galaxies have demonstrated its importance at
all scales. Many expectations from analytical theory are supported, and the
detailed field structure and evolution can be traced in this highly nonlinear
problem. However, limitations of numerical resolution and the different
numerical approximations used can result in qualitatively different outcomes
in various situations, necessitating careful evaluation of reported results.

###  3.1 Plasma scale

The generation of magnetic fields by cosmic accretion shocks from unmagnetised
plasma was hypothesised by Gruzinov (2001) based on consideration of the
properties of gamma ray burst shocks. Schlickeiser and Shukla (2003) and Lazar
et al (2009) made analytic arguments for the Weibel (1959) instability in
shocks being the physical mechanism that could generate these fields when the
shocks drive streams of electrons and ions moving in different directions.
These ideas were supported by initial models of homogeneous Weibel instability
performed with particle-in-cell codes by Medvedev et al (2004) and Sakai et al
(2004).

Rincon et al (2016) demonstrated that magnetic field formation from an
unmagnetised turbulent plasma and subsequent dynamo growth can occur using a
hybrid model that includes an electron fluid and proton particles computed
using a full six-dimensional Vlasov equation. This model suppresses both the
Biermann (1950) battery and Weibel (1959) instability because of its
assumption of isothermal electrons. However, Pusztai et al (2020) used their
fully kinetic model to demonstrate that Landau damping can limit the growth of
kinetic dynamos in some regimes, an effect that cannot be seen in the hybrid
approximation.

By including both positively and negatively charged particles in a fully
kinetic simulation of a hot pair plasma that excludes the Biermann (1950)
battery effect, Sironi et al (2023) and Zhou et al (2023) show that magnetic
fields generated by the Weibel (1959) instability from an unmagnetised plasma
can grow by SSD in driven turbulence, as shown in Fig. 2.

Figure 2: Kinetic models of an initially unmagnetised turbulent flow. Top
Magnetic energy $\epsilon_{B}$ (solid lines), kinetic energy $\epsilon_{K}$
(dotted lines), and magnetic energy at the turbulent driving scale $k_{\rm
int}$. Time is given in units of the crossing time $t_{L}$ at $k_{\rm int}$
for models with domain sizes $L$ in terms of the electron skin depth $d_{e}$
of $L/d_{e}$ ranging from 250 (orange) to 2000 (blue). The dashed black line
shows exponential growth of the rms field $B_{\rm rms}\propto\exp(0.4v_{\rm
rms}t/L)$. Bottom Magnetic power spectra immediately after magnetisation by
Weibel instability (dot-dashed lines) and in the quasisteady state after
dynamo saturation at equipartition (solid lines), with the Kazantsev
($M(k)\propto k^{3/2}$ and Kolmogorov $(M(k)\propto k^{-5/3}$ spectra shown in
black. Image reproduced with permission from Sironi et al (2023), copyright
APS.

The separation between driving and dissipation scales in these kinetic models
remains far smaller than in galactic or intergalactic gas, and models still
need to be done with more appropriate ratios of electron to ion masses.
Nevertheless, the mechanisms described seem robust enough that we should
conclude that such an SSD will occur already behind the accretion shocks that
occur during the assembly of gas into galaxies. As a result, magnetic fields
of at least a few percent of equipartition with turbulence likely provide the
initial condition for any galactic dynamo. Simulations starting with extremely
small seed fields may unphysically delay the onset of magnetohydrodynamic
effects in galactic evolution.

###  3.2 Star formation scale

During the formation of the first stars, dynamo action can already occur in
the turbulent accretion flows assembling the gas into haloes. Kinetic
simulations of this process were discussed in the previous section. Schleicher
et al (2010) presented analytic estimates of the resulting field strengths in
halos, showing that near equipartition fields could already be expected even
prior to the formation of the first stars.

After the accreted gas cools sufficiently to produce an accretion disc
surrounding a protostar, magnetorotational instability in this disc can drive
an SSD (Brandenburg et al, 1995). A semi-analytic model of the SSD in such a
disc by Schober et al (2012) established that fast growth would be expected.
This treatment was extended to quickly accreting young galaxies by Schober et
al (2013).

After the action of an SSD during accretion treated in the previous section,
dynamos in protostellar accretion discs are another source of macroscopic
fields on short time scales during galaxy formation. Ideal AMR MHD models with
Flash (Fryxell et al, 2000) by Sur et al (2010) and Federrath et al (2011b)
demonstrated that an SSD can only be resolved in an accretion disc if the
Truelove et al (1997) criterion for resolving self-gravitating discs is
satisfied with at least 32 cells, as shown in Fig. 3.

Figure 3: The number of cells with which the Jeans length is resolved in a
primordial star formation region determines whether growth of the magnetic
field from an SSD can be resolved. Amplification by spherical gravitational
collapse produces field growth proportional to $\rho^{2/3}$, so further field
growth is evidence of the action of the SSD. Image reproduced with permission
from Sur et al (2010), copyright vy AAS.

In this case, exponential growth can occur with timescales shorter than the
orbital time of the disc. McKee et al (2020) and Stacy et al (2022) show
growth of fields in primordial discs to values characteristic of modern star
forming regions in simulations starting with cosmological structure formation
and zooming in to disc scale. By characterising the numerical diffusivity,
they found that the growth rates were consistent with the non-ideal MHD SSD
theory of Xu and Lazarian (2016).

Whether modern protostellar accretion discs maintain their magnetic fields
through the dynamo action of the MRI depends on how ionised they are. Dead
zones where ambipolar or Ohmic diffusion dominate have fields in balance with
magnetocentrifugally-driven winds, as suggested by local shearing-box models
with Athena (Bai and Stone, 2013), and then confirmed in a global disc model
on a spherical-polar mesh using NIRVANA-III (Gressel et al, 2015). Active
zones sufficiently ionised for the gas to fully couple to the magnetic field
indeed have MRI-driven SSD action, as shown (Brandenburg et al, 1995) using
uniform-grid, shearing-box, open-boundary simulations with the Stagger code
(Nordlund and Stein, 1990; Brandenburg et al, 1996), and in stratified,
shearing-box simulations (Stone et al, 1996) with the ZEUS code (Stone and
Norman, 1992a, b).

###  3.3 Kiloparsec scale

The structure of the magnetic field within the ISM on the kiloparsec scale can
be separated into turbulent scales and the large-scale or MF. The separation
scale is given by the effective 75–200 pc forcing scale of the turbulence
driven by SN blast waves (Joung and Mac Low, 2006; de Avillez and
Breitschwerdt, 2007) and stellar winds from massive stars (van Buren, 1985).
While the blastwave of a single SN evolves over a range of scales from its
stellar radius to tens of parsecs, the effective forcing scale for the ISM
turbulence is controlled by the collective action of multiple SNe combining to
blow superbubbles (Korpi et al, 1999a; Joung and Mac Low, 2006; de Avillez and
Breitschwerdt, 2007; Gent et al, 2013a; Hollins et al, 2017). Small-scale
drivers from low-mass stellar winds, protostellar jets, planetary nebulae, and
cloud dynamics also exist, but these are negligible compared to SNe and, to a
lesser extent, the highest mass stars (Mac Low and Klessen, 2004). The
resulting turbulence can support a turbulent dynamo (or SSD), which amplifies
the magnetic field exponentially.

The MF is organised by processes such as the differential rotation of galactic
discs, spiral arms, and other metagalactic scale processes (Brandenburg and
Subramanian, 2005; Shukurov and Subramanian, 2021). ISM turbulence, can,
however, also be driven by large scale processes such as gravitational
instability (Rafikov, 2001), Kelvin-Helmholtz and other gas-dynamical
instabilities, disc-halo circulation, and Parker instability.

####  3.3.1 SN driven SSD

Theoretical approaches to modelling the magnetic field in the ISM commonly
presume approximate equipartion between the energy density contained within
each of the thermal, kinetic, cosmic ray and magnetic field energy.

However, numerical models of SSD action in the ISM of galaxies suggest that it
produces and maintains magnetic fields of only a few percent of equipartition
with the kinetic energy. This low saturation of the kinematic dynamo is
evident for isothermal simulations at high Mach number (e.g., Haugen et al,
2004; Federrath et al, 2011a, 2014). Models of turbulence simulated at the
value of ${\rm Pr}_{\rm M}>1$ relevant in the ISM show that the level of
saturation increases with Pm, but still appears to reach an asymptotic limit
well below equipartition (e.g., Schekochihin et al, 2002; Schober et al, 2012;
Seta et al, 2020).

With isothermal turbulence driven by pure solenoidal forcing (Beattie et al,
2022) it is possible to reach around one third equipartition for ${\rm
Pr}_{\rm M}=4$ at $\rm Re=500$. Seta and Federrath (2020) and Beattie et al
(2022) show that the SSD grows independently of the properties of the seed
field in individual galaxies, while Garaldi et al (2021) comes to the same
conclusion using self-consistent cosmological simulations. SSD in a two-phase
medium grows more slowly than in isothermal turbulence of comparable Mach
number and Rm (Seta and Federrath, 2022). It does appear that the SSD alone is
insufficient to account for the observed strength of the turbulent magnetic
field in galaxies.

Figure 4: Magnetic energy scaled by time-averaged kinetic energy in a periodic
box driven by SN explosions and radiative heating and cooling to produce a
turbulent multiphase medium including explicit Laplacian resistivity and
viscosity. Numerical resolution in parsecs is given in the legend. The
timescale of exponential growth drops dramatically as resolution improves and
ReM increases. Image reproduced with permission from Gent et al (2021),
copyright by AAS.

With more realistic forcing in a multiphase ISM by SNe, at rates more relevant
to spiral galaxies, but without any drivers for LSD, Balsara et al (2004);
Balsara and Kim (2005) and Gent et al (2021) also find that the SSD saturates
at a few percent of equipartition strength, but much faster than in isothermal
models or the cooler two phase medium, as shown in Fig. 4. Because of the
stochastic nature of the driving, this dynamo acts intermittently (Gent et al,
2021). Field growth is fastest in the hot phase where high sound speed
facilitates high vorticity (Gent et al, 2023), but the strength at which the
SSD saturates is lowest in the hot gas and increases in cooler phases. In
particular, the turbulent flows of diffuse, inhomogeneous, hot gas in SN
remnants and suberbubbles host the most rapid SSD (Gent et al, 2023) and
(Kirchschlager et al, 2023, Fig. 1, third row).

In elliptical galaxies lacking differential rotation there is no driver for
LSD, but SSD remains active, even though SN rates are much lower than in disc
galaxies. The theoretical model of Moss and Shukurov (1996) anticipates a
turbulent magnetic field similar in magnitude to those of spiral galaxies,
while Seta et al (2021) conclude from simulated observations of numerical
turbulence that turbulent magnetic field must be an order of magnitude weaker.

Due to the random distribution of SNe and the multiphase structure of the ISM,
the turbulent properties are highly anisotropic and intermittent. A thorough
investigation of dependence of the turbulent magnetic properties on
temperature phase and SN rate, and on galactic radial and vertical anisotropy
remains outstanding. This would form the basis for modelling simulated
observations of the magnetic structure of the ISM as well as constructing
effective subgrid-scale models of MHD turbulence with which to explore large-
scale magnetic effects.

####  3.3.2 Large-scale fields without LSD

To achieve the field strength observed in spiral galaxies, and also to
comprehensively understand the anisotopy, structure and strength of the
turbulent magnetic field, requires the presence of an LSD.

Some investigations have explored the evolution of a large scale field in a
stratified ISM (e.g., de Avillez and Breitschwerdt, 2005; Hill et al, 2012;
Walch et al, 2015; Girichidis et al, 2016), but with an imposed initial
uniform MF, which cannot be sustained in the absence of an MF dynamo
mechanism. The turbulent field in these models arises primarily through
tangling of the initial field. Comparisons between Balsara et al (2004) and
Gent et al (2021) indicate that at the resolution applied in such models and
using these codes, the SSD would take some 100 Myr to initiate and up to a
gigayear to saturate, so most of these models are likely still in the
unsaturated growth stage of the SSD.

####  3.3.3 SN driven LSD

To solve directly the evolution of the MF dynamo with turbulence driven by SN
explosions Korpi et al (1999) and Korpi (1999) model a stratified disc subject
to rotation using a shearing box. Neither the SSD nor the LSD could be
sustained because the model did not have the required resolution to reach
magnetic Reynolds numbers sufficient to excite the SSD. With open vertical
boundaries and a relatively low vertical extent the gas in these models is
rapidly exhausted and overheated within a few hundred megayears. This is too
soon to amplify the MF dynamo, which it is reasonable to expect be primarily
hosted in the continuous warm phase of the ISM, rather than the isolated cold
filaments or hot bubbles.

With the Nirvana code and double the grid resolution Gressel et al (2008) and
Gressel (2008) model an area 800 pc on a side extending $\pm 2.1$ kpc either
side of the galactic disc. Exploring angular velocities ranging from
$1~{}\Omega_{\rm Sn}$ to $8~{}\Omega_{\rm Sn}$, where the Sn subscript denotes
the Solar neighborhood value, with and without shear and with SN explosion
rates between $0.25~{}\sigma_{\rm Sn}$ and $1~{}\sigma_{\rm Sn}$, they find
the LSD to be effective for $\Omega\gtrsim 2~{}\Omega_{\rm Sn}$. In the
absence of Coriolis forces the shear produces no LSD. The efficiency of the
LSD increases slightly as SN rate reduces. The LSD does not saturate in these
results, but it is clear that the timescale would exceed 1.5 Gyr for
$\Omega\leq 4~{}\Omega_{\rm Sn}$.

Kim and Ostriker (2015) using a similar setup, but with the inclusion of a
shearing radial boundary at $1~{}\Omega_{\rm Sn}$, evolve their models beyond
850 Myr. Where the initial field is strong enough, the turbulent tangling and
any residual SSD are sufficient to sustain an equilibrium total magnetic
energy throughout, but for the model with a weak initial field the final state
remains far below equipartition. This setup has the essential ingredients to
model the LSD, but given the absence of any clear exponential growth it is
very difficult to differentiate SSD from tangling and to exclude the effects
of the imposed field on the results.

With the Pencil Code, Gent et al (2013a) almost double the resolution again to
model a domain of horizontal area $(1.024\mathrm{kpc})^{2}$ and extending
$1.088$ kpc either side of the disc. Applying differential rotation via the
shear-periodic boundary, they confirm an LSD for $\Omega\geq 1~{}\Omega_{\rm
Sn}$. For a model with $\Omega=2~{}\Omega_{\rm Sn}$ they obtain an MF dynamo
with a growth rate of approximately 5 Gyr-1, about twice that for the growth
rate of the turbulent field. This growth, however, measures only the period
800 Myr to 1.05 Gyr, well after the SSD has saturated, but while the LSD
remains an order of magnitude less than equipartition.

From comparison of the evolving mean and fluctuation magnetic energy, it is
evident that the SSD is absent for Gressel et al (2008) while it is very
active for Gent et al (2013a). It is not clear why this would impede the
progress of the LSD at low $\Omega$ for the former, while not in the latter.
The LSD should be quite insensitive to the Rm in these ranges. It poses the
question as to the extent the SSD may hinder or assist the action of the LSD.
Gent et al (2013a) do not obtain magnetic energy spectra with an obvious
separation of scales between the mean and fluctuating magnetic fields, which
instead show a single parabola. They do, however, show that a meaningful
separation of scales can be defined with the fluctuating field occupying
scales typically smaller than 250 pc and the large scale field having scales
in excess of 700 pc.

Figure 5: For all models of Gressel et al (2008) (a) mean and (b) turbulent
magnetic energy evolution. Labels Q, H and F denote SN rates of 0.25, 0.5 and
$1~{}\sigma_{\rm Sn}$, respectively, while the numbers give the angular
velocity scaled by $~{}\Omega_{\rm Sn}$. (c) Evolution of the horizontally
averaged and the turbulent magnetic field strength for model H4 and (d)
evolution of total magnetic energy for all models from Gent et al (2023),
where L denotes resolution of $\delta x=4$ pc and H $\delta x=1$ pc. The
suffix ‘cl’ denotes use of SN clustering and ‘B0’ the ongoing exclusion of an
MF to isolate SSD activity. In both sets the numbers in the model labels
denote $\Omega$ scaled by $~{}\Omega_{\rm Sn}$. (e) Evolution of the total
magnetic energy in model L2 and the MF energy calculated from horizontal
averages and convolution with a Gaussian kernel. Figure adapted from Gressel
et al (2008) and Gent et al (2023).

Figure 5 compares the models of the galactic LSD by Gressel et al (2008) and
Gent et al (2023). Figure 5(a) shows that for SN rate $\sigma=~{}\sigma_{\rm
Sn}$ the growth rate of the MF is proportional to $\Omega$ for $\Omega\geq
2~{}\Omega_{\rm Sn}$. Gressel et al (2008) do not find an MF dynamo for
$\Omega=~{}\Omega_{\rm Sn}$. However, Gent (2012) finds an MF dynamo is
present at $\Omega=~{}\Omega_{\rm Sn}$. Model F1 from Gressel et al (2008)
extends only to about 150 Myr, though, which may have been too short to
dissipate transients from the initial azimuthal seed field. The linear
dependence of the growth rate on $\Omega$ is consistent with the theoretical
expectation of Brandenburg and Subramanian (2005) for a supercritical MF
dynamo, which suggests the MF dynamo will be well supported in fast and even
moderately rotating galaxies where the rotation curve is flat. In the range
$0.25~{}\sigma_{\rm Sn}\leq\sigma\leq~{}\sigma_{\rm Sn}$ the LSD is
anticorrelated with $\sigma$.

The interaction of the SSD and LSD is explored by Gent et al (2023). With the
growth and saturation properties of the SSD better defined by Gent et al
(2021, 2023), the effect of an active LSD on the SSD is isolated by
systematically removing any horizontal MF from the components of the magnetic
field in a stratified ISM with and without differential rotation. In the case
with additional energy injected by the large scale shear from the differential
rotation, the SSD grows slightly faster and saturates proportionately higher
than in the non-rotating system (compare the growth rates of L0-B0 and L2-B0
in Fig. 5(d)). When the LSD is not excluded in the differentially rotating
model (L2 in Fig. 5(d), which does start with a higher seed field), the SSD
has a slightly lower growth rate during the kinematic phase. The LSD extracts
energy from the SSD equivalent to the additional energy injected by the large
scale shear. The turbulent field, however, does not stop growing after the SSD
would saturate in the absence of the LSD, though growth slows as the LSD
transitions into the dominant mode of the dynamo. The effect of the SSD on the
LSD requires further study.

From Fig. 5(c) it can be seen that, even though Gressel et al (2008) does not
show SSD activity, the turbulent field rapidly increases to around ten times
the strength of the seed field. This is supported by the analysis of Gent et
al (2023) during the period where the SSD is present, but saturated. The
saturation state of the SSD is confirmed by models L0-B0 and L2-B0 shown in
Fig. 5(d), where the MF is continuously removed. When the SSD saturates for
model L2 (see Fig. 5(e)), the turbulent field continues to grow due to the
presence of the LSD even though it is over ten times stronger than the growing
MF which it is continuing to generate. Tangling by SN turbulence therefore can
produce a turbulent field an order of magnitude stronger than the MF.

(a)(b)(c)(d)

Figure 6: Evolution of the horizontally averaged mean magnetic field
components, normalised by the maximum field strength at each time. From model
H4 of Gressel et al (2008) the field components (a) $B_{r}$ and (b) $B_{\phi}$
are shown. From model L2 of Gent et al (2023) the field components (c) $B_{x}$
and (d) $B_{y}$ are shown. Both models are seeded by fields with maximum rms
strength at the midplane $<1$ nG, and grow to strengths of several microGauss.

To understand how we might expect the dynamo generated mean field to appear
near the galactic plane, we compare in Fig. 6 time-latitude diagrams of
horizontal averages for model H4 from Gressel et al (2008) and L2 from Gent et
al (2023). Model H4, in panels (a) and (b), exhibits a quite regular field,
subject to reversals with a period below 1 Gyr, and has quadrupolar structure.
By contrast, model L2, in panels (c) and (d), initially has a rapidly
fluctuating structure both in time and latitude due to the presence of a
vigorous SSD. Well after the saturation of the LSD, the field becomes more
regular. Early on there is more asymmetry with numerous reversals between
different latitudes. As the field becomes more quadrupolar there is then a
reversal around 2 Gyr, after which the field saturates to a steady azimuthally
dominant field. The comparsion between the models is most relevant up to about
2 Gyr, when both models are still subequipartition. It is noteworthy that the
scale height of the magnetic field expands in the late stages when the
strength of the midplane magnetic field becomes dynamically significant with
respect to the kinetic energy. Whether the presence of reversals in the field
is impacted by the rate of rotation, which differs by a factor of two here,
the scale height of the disc, the SN rate, or the horizontal size of the
domain requires further investigation.

####  3.3.4 Magnetic buoyancy instability

Another source of dynamo instabilities explored at these scales is the
magnetic buoyancy instability (MBI) or when attenuated by cosmic rays, the
Parker instability. A suitably stratified disc, with sufficiently strong
parallel magnetic fields, has been initialized to examine the linear stage of
the instability in an isothermal galaxy by Rodrigues et al (2016). In such
cases the mean magnetic field is imposed, presumably acquired by prior
galactic dynamo, and the exponential growth of the velocity purturbations
driven by the MBI drives a secondary dynamo. The resulting magnetic field has
an expanded scale height and can attain superequipartition strength (Rodrigues
et al, 2016). The growth rates during the linear stage of the instability
obtained from the analytic linearised solutions of Giz and Shu (1993) are of
order 20–30 Gyr-1. Numerical solutions yield somewhat lower growth rates. When
the numerical ratio of magnetic to gas pressure is high the growth rate of the
linear instability is enhanced.

Including solid body rotation is anticipated analytically to dampen MBI growth
(Ryu et al, 2003) but Kim et al (2001) find it makes little difference to the
growth of the linear instability, though it can enhance the formation of
turbulent structure in the magneitc field. The nonlinear state of MBI has been
examined in the context of disc galaxies without (Tharakkal et al, 2023b) and
with (Tharakkal et al, 2023a) rotation. MBI is characterised by an
$\alpha$-effect with opposite sign of helicity to that typically driving the
$\alpha\Omega$-dynamo (Tharakkal et al, 2023a; Qasi et al, 2023). When
rotation is included (Qasi et al, 2023) MBI induces periodic reversals in the
large-scale field (Tharakkal et al, 2023a) and the cosmic ray energy density
is higher. Where the large-scale magnetic field is generated via a dynamo, the
MBI has little effect on the stratified structure of the gas and cosmic rays,
while still expanding the vertical scale height of the magnetic field, whereas
when the field is imposed, it expands the scale height of the field and cosmic
rays so far as to diminish their pressure support, resulting in a thinner gas
disc during the nonlinear stage (Qasi et al, 2023). In Qasi et al (2023) an
$\alpha^{2}$-dynamo is applied to a nonrotating disc, so the effect of
differential rotation and the resulting $\alpha\Omega$-dynamo remains to be
explored, as does the interaction between SN-driven turbulence and the MBI.

####  3.3.5 Magnetorotational instability

Investigating MRI-driven turbulence in a two-phase stratified ISM, Piontek and
Ostriker (2005) find that, even without SNe, MRI is sufficient to excite
strong turbulence in the diffuse ISM away from the midplane or in the outer
galaxy where the ISN in the disc becomes more diffuse. MRI inhibits self-
gravitating and thermal instabilites, likely retarding star formation in the
outer galaxy (Piontek and Ostriker, 2007). While turbulence can increase the
proportion of the ISM in thermally unstable temperatures, the bimodal
character of the ISM persists, consistent with the findings of Gent et al
(2013b) for the trimodal character in the presence of SN-driven turbulence.
The strength of the magnetic field becomes somewhat independent of gas number
density (Piontek and Ostriker, 2005, 2007). However, decomposing the magnetic
field between mean and turbulent contribution Evirgen et al (2017) shows the
turbulent magnetic field to correlate more with the hot phase and the MF with
the warm and cold phases.

###  3.4 Galactic scale

Schober et al (2013) estimate the growth rate for an accretion-driven dynamo
in a spherical protogalaxy to be

$\Gamma\simeq kv_{k}\mathrm{Re}^{1/2}$ (23)

for wavenumber $k$, velocity at that scale $v_{k}$, and Reynolds number Re,
where the exponent is appropriate for Kolmogorov turbulence. They use this to
analytically estimate that an accretion-driven small-scale dynamo can amplify
a weak seed field to near-equipartition with the turbulent kinetic energy in
tens to a few hundred megayears. These are almost certainly overestimates of
the time required for these fields to grow, since equipartition level small
scale fields already grow in the first protostellar accretion discs, as
described in Sect. 3.2.

Models of field growth in modern galaxies that start with small seed fields
can therefore only be interpreted as physics experiments exploring the
behavior of dynamos. However, they can not be intepreted fruitfully to
understand the astrophysical evolution of galaxies in the presence of magnetic
fields, because of the extended period of weak field present in these models
that does not occur in any scenario including cosmological growth of galaxies.

####  3.4.1 Mean-field (MF) models

There exist proof-of-concept non-linear MF models employing the full dynamic
quenching formula with various types of helicity fluxes. For example, the
studies of Shukurov et al (2006) and Sur et al (2007) investigate galactic
fountain-driven fluxes and the latter study also investigates helicity fluxes
resulting from anisotropic turbulence due to the presence of shear, a scenario
proposed by Vishniac and Cho (2001). Also diffusive helicity fluxes have been
included by Chamandy et al (2014). All these studies were able to demonstrate
that including helicity fluxes results in LSD saturating at the equipartition
value with turbulence being excited and maintained. Also, evidence for the
heuristic nonlinearity (Eq. (15)) producing results consistent with the more
complete quenching models was found (e.g. Chamandy et al, 2014). This is
reassuring, as both old and new MF models still often use the heuristic
nonlinearity (e.g. Liu et al, 2022).

####  3.4.2 MHD models

The interrelation between the SSD and the LSD is a key question at the
galactic scale. The physical timescale for the SSD is far shorter than for the
LSD. However, insufficient numerical resolution can dramatically slow the
action of the SSD (Gent et al, 2023), causing the two timescales to apparently
overlap in numerical simulations. Nevertheless, numerical simulations do
produce toroidal fields that generally reproduce the observed topology, as
shown, for example, in Fig. 7.

Figure 7: Magnetic field strength projected edge-on and face-on in an isolated
disc galaxy. Directed stremalines of magnetic field are shown in white. The
length scale is given in kiloparsecs on the left axis. Image reproduced with
permission from Butsky et al (2017), copyright by AAS.

##### Small-scale dynamo

The presence of SSD in global galaxy models has been demonstrated using two
primary diagnostics: the exponential growth of magnetic energy and the
agreement of the power spectrum of the field with the $k^{3/2}$ prediction of
Kazantsev (1968). The first global MHD models of field development in a disc
galaxy were performed by Wang and Abel (2009) and Kotarba et al (2009),
neglecting turbulence driven by stellar feedback. These works already
demonstrated that the turbulence generated by gravity and differential
rotation in galaxies is sufficient to produce exponential field growth to a
few percent of equipartition, as shown in Fig. 8, although at a rate limited
by the low numerical resolution and turbulence strength. Kotarba et al (2009)
further demonstrated that constraint of magnetic helicity by evolution of the
field using Euler potentials (Price and Bate, 2007) suppresses dynamo
activity.

Figure 8: Magnetic energy (asterisks) in disc for a dwarf galaxy simulation at
redshift two without stellar feedback at 26 pc resolution compared to kinetic
energy (squares) and thermal energy (diamonds). including SN feedback at the
finest resolutions given in the legend. Exponential growth with the given time
scales is shown by the black lines. Image reproduced with permission from Wang
and Abel (2009), copyright by AAS.

The first diagnostic, exponential growth of magnetic energy to a few percent
of kinetic energy in global galaxy simulations has been confirmed by multiple
groups in the last decade. Models without explicitly modeled feedback-driven
turbulence include Pakmor and Springel (2013); Rieder and Teyssier (2016);
Pakmor et al (2017); Rieder and Teyssier (2017); Steinwandel et al (2020) and
Pfrommer et al (2022). In most of these models, feedback is assumed to
contribute to a turbulent pressure variable but the only flows contributing to
an SSD are produced by gravitational dynamics of the galactic gas or large-
scale galactic winds. A number of studies without explicit feedback did not
explicitly report the exponential growth of magnetic energy, but rather that
of related but not identical quantities such as the average magnetic field
strength, including Pakmor et al (2014); Steinwandel et al (2019) or the
median magnetic energy density Pakmor et al (2017). Models that do include
stellar feedback, driving turbulence at smaller scales and resulting in faster
magnetic energy growth, include Su et al (2017); Butsky et al (2017);
Ntormousi et al (2020) and Martin-Alvarez et al (2022).

Figure 9: Magnetic energy in domain for a dwarf galaxy simulation including SN
feedback at the finest resolutions given in the legend. Exponential growth
with the given time scales is shown by the black lines. Image reproduced with
permission from Rieder and Teyssier (2016), copyright by the author(s).

Rieder and Teyssier (2016) clearly demonstrated the importance of numerical
resolution in full galaxy models, as shown in Fig. 9. Note that these are
still low resolutions and thus slow growth rates compared to those reached in
kiloparsec-scale models (see Sect. 3.3.1), where growth time scales drop to
under a megayear (Fig. 4 at 0.5 pc resolution. Martin-Alvarez et al (2022)
emphasised the importance of uniformly resolving galactic discs rather than
allowing resolution to follow density to correctly capture the growth of
magnetic fields in lower-density regions where they grow faster.

Figure 10: Power spectrum of magnetic energy at listed times in a model of an
isolated galactic disc. These spectra show the Kazantsev slope at large
scales, and a peak that occurs at increasing scale over time, as expected from
a saturating SSD. The data and the guiding lines (dotted lines) have been
regularised by $k$. The kinetic power spectrum at the final time is also shown
(dashed line), with a slope steeper than expected for incompressible
Kolmogorov turbulence, consistent with the compressible turbulence in this
galaxy. Image reproduced with permission from Butsky et al (2017), copyright
by AAS.

The second diagnostic of SSD action is the development of Kazentsev spectra in
models followed by the shift of the peak in magnetic energy from the
dissipation scale to longer wavelengths as saturation sets in at the
dissipation scale, as shown, for example, in Fig. 10. This result has been
found in full-galaxy models by Rieder and Teyssier (2016, 2017, 2017); Butsky
et al (2017); Pakmor et al (2017); Martin-Alvarez et al (2018); Steinwandel et
al (2019); Martin-Alvarez et al (2022) and Pfrommer et al (2022). Steinwandel
et al (2019) and Ntormousi et al (2020) noted that after saturation, an
Iroshnikov spectrum with magnetic energy proportional to $k^{-3/2}$ appears,
while Pfrommer et al (2022) argues that the saturated spectrum is closer to
Kolmogorov with a $k^{-5/3}$ dependence, although in an accretion-dominated
model of a primordial galaxy.

Martin-Alvarez et al (2018, 2022), and Ntormousi et al (2020) using AMR and
Steinwandel et al (2020) using SPH explicitly analysed SSD and LSD in global
galaxy simulations including a subgrid model of SN feedback (Springel and
Hernquist, 2003; Dubois and Teyssier, 2008) They all compared models with and
without feedback and found that both showed the same exponential growth during
the first 0.5 Gyr. At the resolution of these models, this is the expected
timescale for the action of the SSD, suggesting that the primary driver of the
turbulence required for SSD may be accretion (Klessen and Hennebelle, 2010)
rather than SNe in these models that include circumgalactic gas.

##### Large-scale dynamo

Evidence for the presence of the LSD in global models including SSD has only
started to be presented within the last several years. An early paper that
likely showed LSD action, but did not discuss it explicitly, is Pakmor and
Springel (2013). They compared the growth of volume-averaged magnetic pressure
to that of total and thermal pressure in galaxies of halo mass varying from
$10^{9}$ to $10^{12}\,M_{\odot}$, so from dwarf to Milky Way sizes. As shown
in Fig. 11, all of these galaxies showed at least modest further growth in
magnetic pressure after saturation of the SSD. However, the saturation of this
apparent LSD depended strongly on the halo mass, with the smallest, slowest
rotating halos with the shallowest potentials having magnetic pressures more
than an order of magnitude below the thermal pressure, while the Milky Way
mass halo had magnetic pressure exceed thermal pressure.

Figure 11: Magnetic (dotted lines), thermal (dashed lines), and total (solid
lines) pressure for halos with masses of $10^{9}$ to $10^{12}\,M_{\odot}$. All
pressures are volume-averaged over the galactic discs. These models show LSD
behavior whose relative and absolute saturation strength depends strongly on
potential well depth. Figure from Pakmor and Springel (2013).

As helical turbulence is believed to be one of the most important turbulent
ingredients for an LSD, Ntormousi et al (2020) directly measured the turbulent
kinetic and current helicities (see Eqs. (17) and (18)). A median filter with
390 pc (8 cell) size was used to separate turbulent and mean components of the
velocity and magnetic field. In these AMR models, the magnetic energy remained
below 0.1% of the turbulent kinetic energy, so these are representative of an
unsaturated dynamo. This is supported by the multiple sign reversals across
the midplane seen in Fig. 12.

Figure 12: Integrated turbulent kinetic and current helicities over time above
and below the galactic midplane in a global AMR model embedded in a halo. The
symmetric sign reversals suggest an unsaturated dynamo. Image reproduced with
permission from Ntormousi et al (2020), copyright by ESO.

Figure 13: Vertical dependence of radial (left) and toroidal (right) magnetic
fields at times given in the legend for an SPH global model. Antisymmetry
around the midplane suggests dipole structure, while symmetry suggests
quadrupolar structure. Some indications of both structures are seen at
different times, but the full structure is more complex, and influenced by the
buildup of a magnetic tower flow. Image reproduced with permission from
Steinwandel et al (2020), copyright by the author(s).

Steinwandel et al (2020) and Ntormousi et al (2020) both found that the
classical dipolar and quadrupolar field structures do not fully describe the
large-scale field resulting from their models of LSD in global galaxies, with
a rather more confused pattern probably consisting of multiple eigenfunctions
growing simultaneously instead appearing, as shown for example in Fig. 13.
Similar results were found by Pakmor et al (2018) and Reissl et al (2023)
while computing Faraday rotation measures, as discussed in Sect. 3.6 below.

###  3.5 Cosmological scale

The inclusion of the cosmological environment in models of dynamos in global
galaxies allows a more self-consistent treatment of accretion on to galaxies
and the influence of interactions with neighbors. The cost is a low numerical
resolution compared to the models of single galaxies discussed in the previous
section. However, even with SPH mass resolutions as high as $2\times
10^{7}\,M_{\odot}$, Beck et al (2012, 2013) demonstrated that an SSD produces
exponential field growth during cosmological halo formation, accretion, and
merging even though disc formation was not resolved. Radiative cooling
balanced by star formation driven feedback (implemented with a star formation
threshold of only 0.13 cm-3 in these low resolution models) produced the
fastest growth rates, with e-folding times of order $10^{7}$ yr, and
equipartition being acheived with the turbulence. Saturation occurred with
field strengths of $10^{-5}$ G at the centers of halos.

More recent simulations have followed the formation and evolution of star-
forming galaxies by zooming in to cosmological simulations, improving the
numerical resolution within galaxies in order to resolve dynamo action there.
Pakmor et al (2014) used the Voronoi-mesh code AREPO to model a Milky-Way
sized galaxy with mass resolution of $5\times 10^{4}\,M_{\odot}$ and a
gravitational softening length of 340 pc, while Rieder and Teyssier (2017)
used the AMR code RAMSES to model a dwarf galaxy down to redshift $z=4$ with
most refined resolution of 22.5 physical pc (requiring successive refinements
during expansion of the universe). These workers found that an SSD drives
initial exponential field growth. The small-scale turbulence is initially
produced by chaotic accretion, and then subsequently by star formation. In
these works, increasing the resolution increases the growth rate of the SSD as
expected. After saturation of the SSD, they find linear growth of the LSD in
the disc, while fields in the galactic center remain at the value produced by
saturation of the SSD absent significant rotation.

Larger scale efforts include the Auriga simulations (Grand et al, 2017), a set
of 30 zoom-in models using AREPO of Milky Way mass haloes drawn from the dark
matter only version of the Eagle model (Schaye et al, 2015) at the same mass
resolution of $5\times 10^{4}\,M_{\odot}$ used by Pakmor et al (2014).
Although these models still do not capture the physical growth rate of the
interstellar SSD (Sect. 3.3), they do follow the galactic LSD acting on the
fields generated by halo-scale SSD action during the galaxy formation process.
As the galactic LSD appears insensitive to resolution in this regime (Sect.
3.3), the structures modeled in this way may reflect physical processes and
not just numerical dissipation.

Figure 14: Median energy densities over time drawn from the 30 zoom-in models
of the cosmological Auriga simulations. Measurements are performed in rings at
the given radii. At 10 kpc galactocentric radius, the magnetic field only
approaches equipartition with the turbulent kinetic energy after 10 Gyr of
evolution. Image reproduced with permission from Pakmor et al (2017),
copyright by the author(s).

Perhaps the highest resolution cosmological zoom-in model to date is described
by Martin-Alvarez et al (2023), who have simulated a dwarf galaxy in its
cosmological context with 7 pc finest grid resolution using the RAMSES AMR
code. The MHD runs were terminated at redshift $z=3.5$ for reasons of
computational cost. At that time, fields of 2–5 $\mu$G had built up in their
discs, consistent with observed values of $z=0$ dwarf galaxies (Chyży et al,
2011). Unfortunately, these models have not yet been analysed in any detail
for their dynamo properties.

###  3.6 Observational comparisons

The correlation between spiral structure in magnetic fields and spiral
structure in stars and gas is observed to be weak (e.g. Beck, 2015; Bittner et
al, 2017), with the strongest fields often appearing in the interarm regions.
This behavior is reproduced by the global model of a Milky-Way mass galaxy
including a circumgalactic medium by Steinwandel et al (2019). They
hypothesise that the stronger turbulence in the lower-density, hotter interarm
regions leads to faster field growth there.

Figure 15: Radial profiles of average pitch angle of magnetic field over time
for a Milky Way mass galaxy, at times given in the legend. Negative pitch
angles consistent with observations predominate, although positive pitch
angles occur in regions heavily perturbed by magnetic tower outflows. The
purple triangles show observations of M31 by Fletcher et al (2000); these are
consistent with the Milky Way value of $-15^{\circ}$ at the Solar circle. The
black line is from a dynamo model of M31 by Moss et al (1998). Image
reproduced with permission from Steinwandel et al (2020), copyright by the
author(s).

The pitch angle of the spiral structure in magnetic fields is another point of
comparison between global simulations and observations. Models show a great
deal of scatter in the pitch angle, as shown in Fig. 15, with shallow negative
pitch angles consistent with observations being frequent, but by no means
universal. The situation is further confused by accretion, outflows, and
interactions with other galaxies, all of which can distort the magnetic field
sufficiently for the pitch angle to become positive or strongly negative.

Figure 16: Average depth along the line of sight of contributions to the
product of the magnitude of the magnetic field and the density of thermal
electrons measured by rotation measure. The simulated observations are
performed for one of the Auriga zoom-in models of galaxies in a cosmological
simulation. Particularly at high latitudes, the dominant contribution comes
from the local neighborhood rather than the disc as a whole. Figure adapted
from Pakmor et al (2018).  Figure 17: Comparison between the all-sky rotation
measure maps derived by (a) Oppermann et al (2012) and (b) Hutschenreuter et
al (2022) to the rotation measure map from (c) an unmodified Auriga zoom-in
model of a galaxy, and (d) from the same model and viewpoint but with a
subgrid model of ionising radiation from star clusters determining the thermal
electron distribution. The additional small-scale structure reproduces the
statistical small-scale structure distribution of the observations far better.
Image reproduced with permission from Reissl et al (2023), copyright by the
author(s).

Carteret et al (2023) used a shearing box model of a section 0.8 kpc on a side
of a galactic disk to measure the structure of Faraday rotation and
polarisation due to dynamo-derived fields. They found that different
assumptions about the thermal electron density distribution, not explicitly
modeled in their simulation absent radiative transfer of ionising radiation,
drive leading order variation in the magnitude and scale distribution of the
rotation measure for the same magnetic field structures. Anisotropic structure
in the rotation map occurred under both assumptions, though.

Pakmor et al (2018) and Reissl et al (2023) modeled observed Faraday rotation
on the sky as seen from within one of the galaxies from the Auriga
cosmological zoom-in simulations (Grand et al, 2017) at roughly the Solar
circle. These models were compared to the observed distributions derived by
Oppermann et al (2012) and Hutschenreuter et al (2022) from measurements
towards extragalactic sources. Pakmor et al (2018) found agreement on the
strength and qualitative distribution of intermediate to large scale structure
in the distribution, but found much less small-scale structure than in the
observations. They found that the structure within a few kiloparsecs dominates
the Faraday rotation, particularly at high latitudes, rather than any global
dynamo-formed dipole or quadrupole field, as shown in Fig. 16. Reissl et al
(2023) used cluster population synthesis and explicit transfer of ionising
radiation (Pellegrini et al, 2020) to model the distribution of thermal
electrons down to a much smaller scale, allowing polarised radiative transfer
calculations from a vantage point within a superbubble at roughly the Solar
radius. Consistent with the results of Carteret et al (2023), they find that
the addition of the subgrid model of cluster ionisation dramatically improves
the agreement of the model with the observations at smaller scales as shown in
Fig. 17 and quantified in the paper.

Figure 18: Slices of density and magnetic field magnitude for Milky-Way mass
models including circumgalactic medium structure after 2.4 Gyr of evolution.
The magnetic tower outflow driven by magnetic pressure produces X-shaped
structures in the circumgalactic medium similar to those observed. Image
reproduced with permission from Steinwandel et al (2019), copyright by the
author(s).

A strong magnetic tower flow (Fig. 18) with a velocity of a 400–500 km s-1
driven by the dynamo-driven increase of magnetic pressure in the center of the
modeled galaxy occurs in the model of Steinwandel et al (2019). This structure
appears consistent with the X-shaped magnetic field structures seen emerging
from disc galaxies in observations of radio polarisation (e.g. Veilleux et al,
2005). However, this is a transient effect that lasts only a few hundred
megayears in the model, raising the question of how such a structure can be
commonly observed in modern galaxies.

A detailed comparison of the interaction between fields generated by the SSD
and cosmic rays in an accretion dominated galaxy to the correlation between
far-infrared and radio synchrotron luminosities was presented by Pfrommer et
al (2022). Fig. 19 shows that the correlation can be reproduced by fully
saturated SSDs in these galaxies, a result that is robust to the details of
the modelling of cosmic ray travel and the orientation of the galaxies.
Carteret et al (2023) examined polarisation of synchrotron emission in their
shearing box model with an assumed distribution of cosmic rays. They showed
that the cosmic ray distribution and the field both matter for the synchrotron
emission, and that the observed structure further depends on the wavelength of
observation.

Figure 19: Comparison of observations of the correlation between radio
synchrotron luminosity at 1.4 GHz and far infrared luminosity at 8–1000 $\mu$m
(open circles) to luminosities computed from simulations of accretion
dominated galaxies with the Voronoi mesh code AREPO. The semi-transparent
symbols and tracks show results from unsaturated dynamos and their evolution,
while the solid colors represent fully saturated dynamos. Models with only
cosmic ray advection are labeled “CR adv” in the legend, while models also
modelling cosmic ray diffusion are labeled “CR diff”. The legend also gives
halo masses and orientation for the different galaxies. Image reproduced with
permission from Pfrommer et al (2022), copyright by the author(s).

Another comparison between simulations including cosmic rays and observables
was done by Ponnada et al (2022). They compared synthetic rotation measures to
observational constraints from the circumgalactic medium around galaxies. They
found rotation measures consistent with observed upper limits, but low enough
that they would not be detected with current instruments, suggesting that
improved observations will be needed to fully constrain these models.

##  4 Properties of turbulence

Numerical models have been used to extract dynamo-relevant properties of
turbulence. Such studies have two purposes. Firstly, comparisons with
analytical theories can be made. Secondly, some ill-known parameters and
properties can be deduced, which can be vital to constructing simpler and less
resource-demanding numerical approaches.

Figure 20: Time-latitude diagrams of horizontally averaged $\alpha$-effect for
a SN-driven model with twice the galactic rotation and shear rate from (a) a
lower-resolution model integrated up to the saturation of the dynamo and (b) a
higher-resolution model with a far shorter time extent, but illustrating the
detailed time evolution. Adapted from Gent et al (2023). On the right we show
(c) the kinetic $\alpha$ value and (d)  the magnetic $\alpha$ value from a
Parker-unstable system without SNe, but with galactic rotation rate. Adapted
from Tharakkal et al (2023a).

###  4.1 Retrieval of turbulent transport coefficients

Semi-analytic approaches (see, e.g. Ferrière, 1998, and references therein)
showed early on that isolated SNe cause only weak turbulent effects relevant
for the dynamo, while superbubbles could produce more significant effects in
the required range for producing a dynamo number above the critical one needed
for LSD action. These models, however, excluded SN interactions, such as
colliding shock fronts, where non-linear interactions, such as generation of
vorticity through the baroclinic effect (arising when temperature and entropy
gradients are misaligned) or vortex stretching, have been reported to generate
significant amounts of vorticity (see, e.g., Korpi et al, 1999b; Käpylä et al,
2018) To study such effects, full MHD models are needed.

The most straightforward way of estimating the turbulent inductive and
diffusive effects from numerical simulations is to use kinetic and magnetic
helicity and turbulent intensity as a proxy for them (see Eqs. (17) and (18)).
To date, this method has only rarely been implemented to diagnose the full MHD
simulations. Recently, Gent et al (2023) computed the proxy for the $\alpha$
effect this way in their SN-forced, stratified, rotating and shearing local
model, exciting both SSD and LSD, and ran their lower-resolution models until
the saturation of the LSD (see Fig. 20 (a) and (b)). In the growth stage of
the LSD the kinetic $\alpha$ dominates, is smooth and of the expected sign
($+$ in the upper and $-$ in the lower part of the disc) in regions of more
laminar inflows caused by cooling. On the other hand, it is strong and
spatially incoherent with undetermined sign in the strong outflows caused by
clustered SN activity. The unexpected sign of $\alpha$ in the latter case
seems to originate from the opposite disc plane, which indicates that some
non-trivial helicity fluxes might occur across the midplane of the disc. The
magnetic $\alpha$ appears to have two clear roles: to revert the sign of the
$\alpha$ effect during a reversal of the azimuthal magnetic field at
approximately 2 Gyr, and to saturate LSD growth in the warm gas concentrated
around the midplane.

Tharakkal et al (2023a) derived the $\alpha$ contributions in the same way
from a Parker unstable stratified disc otherwise similar, but without SNe and
Milky Way solar neighborhood rotation and shear rates ((Fig. 20 (c) and (d)),
while Gent et al (2023) had twice the Solar neighborhood values. A noteworthy
feature, in comparison to the SN-driven system in panels (a) and (b) is that
the kinetic $\alpha$ effect is systematically of the opposite sign (up to 2
Gyr). The smooth parts of the $\alpha$ effect in the different cases are
roughly of the same magnitude, while significantly stronger only in the SN
bursts and the resulting outflows. The higher rotation and shear rates in Gent
et al (2023) most likely also enhance the SN-driven $\alpha$ magnitudes, so
they might, in fact, be rather equal in a consistently similar system. The
early change in sign of the magnetic part of the $\alpha$ effect is similarly
relatable to a reversal of the azimuthal magnetic field in the Parker-unstable
case. Also, in this case the LSD is saturated by the growing magnetic
contribution of the opposite sign. Comparing these two studies it is evident
that the Parker-instability induced dynamo is of significance and operates in
a markedly different way than an SN-driven dynamo, while both are expected to
occur in the SN-active part of galactic discs. Hence, future numerical studies
including both these mechanisms is of great importance. Such studies, however,
require large enough domains also in the SN-forced models to capture the
Parker unstable magnetic loops.

Another numerical method to retrieve a limited set of the turbulent transport
coefficients is the imposed field method. It relies on the Taylor expansion of
the mean emf according to Eq. (10) under the assumption that all higher order
terms than the linear term containing the $\alpha$ coefficients vanish; this
requires the imposed magnetic field to be uniform. One imposes a field in some
direction, and then measures the mean electromotive force components and the
mean magnetic field, and finally solves for the relevant $\alpha$ tensor
component(s). The problem is that non-uniform magnetic fields can be induced,
leading to the formation of currents relatively quickly. These currents cause
the breakdown of the method unless the magnetic field is reset at suitable
time intervals (see, e.g., Ossendrijver et al, 2002; Hubbard et al, 2009).
This method has not been applied very extensively in the context of galactic
dynamos, although early attempts (Fig. 21) resulted in a successful
measurement of the main $\alpha$ effect component generating the mean radial
field from the azimuthal one (Korpi, 1999). There was no explicit imposed
field applied, as the system itself was dynamo active, and resulted in the
amplification of mean fields. Due to poor resolution, the dynamo was near its
critical ReM for being active. Therefore, the field never reached high
strength due to the limited duration of the simulation runs, and thus remained
in the kinematic regime. Indications of an $\alpha\Omega$ dynamo with an
$\alpha_{y}$ coefficient of roughly $\pm 6\,\mathrm{km\ s}^{-1}$ were found.
This was a factor of a few in excess of the nearly contemporary semi-analytic
findings that neglected SN interactions, thus indicating that the SNe are,
indeed, important.

Figure 21: The mean azimuthal magnetic field vs. the mean electromotive force.
Averages are taken over space and over a 10 Myr time interval, separately for
the (left) upper and (right) lower disc plane of the Cartesian shearing box
simulation. This time period roughly corresponds to the turnover time in the
local simulation domain. The blue dashed lines show $\alpha_{y}=\pm 6$ km
s$-1$ correlation coefficient from the relation
$\overline{\mbox{\boldmath${\cal E}$}{}}_{y}=\alpha_{y}\overline{{\bm{B}}}$.
Adapted from Korpi (1999).

Brandenburg et al (1995) studied accretion disc dynamos driven by MRI in a
Keplerian shearing box model. They found a self-sustained dynamo, and applied
a similar method to measure the $\alpha_{y}$ component as Korpi (1999). As a
result, a clear correlation between the mean emf and theazimuthal magnetic
field was again found, but the sign of $\alpha_{y}$ was the opposite: negative
in the upper and positive in the lower hemisphere, similar to the Parker-
unstable system of Tharakkal et al (2023a). This suggests that dynamos driven
by magnetic buoyancy might all have this feature in common. Also galactic
discs are susceptible to MRI (Piontek and Ostriker, 2005, 2007) although their
rotation laws are not Keplerian, but typically flat. Therefore, although the
Brandenburg et al (1995) results are not directly applicable, they could
indicate that multiple different dynamo mechanisms could be at play in
galaxies, if SN-driven and MRI-driven turbulence co-exist. This, however,
seems unlikely, as SN activity has analytically been shown to suppress MRI in
the inner 15 kpc region in a Milky Way like galaxy Korpi et al (2010). MRI-
driven turbulence and dynamos could, however, explain the unexpectedly large
velocity dispersions in the extended HI discs of galaxies Sellwood and Balbus
(1999); Tamburro et al (2009).

A method suitable for numerically extracting the turbulent transport
coefficients from anisotropic systems, such as galactic discs, which
definitely have strong anisotropies arising from stratification, (non-uniform)
rotation, and the presence of large-scale magnetic fields, was developed by
Brandenburg and Sokoloff (2002). In such systems, $\alpha$ and $\beta$ are
tensors instead of scalars. Thus, to extract all the tensor components, one
needs many equations. Brandenburg and Sokoloff (2002) proposed that enough
equations can be collected by forming various moments of the mean field
components and their derivatives. Next, correlations with turbulent emf, again
utilising the ansatz Eq. (10), are sought for. All of these quantities are
measurable from a full MHD simulation. Their formulation of the resulting
multidimensional regression problem was either spatially local, considering
the value of the magnetic quantities and the turbulent emf at one single point
only, or else non-local, considering a convolution kernel around the location.
This method is referred to as the method of moments or the correlation method,
and it has been extended to further improve the fitting procedure by using SVD
(Simard et al, 2016; Bendre et al, 2020; Bendre and Subramanian, 2022), which
usually prevents unphysical solutions such as negative values of $\beta$. The
method is known to work best for time-variable magnetic fields, and to
underestimate the magnitude of $\beta$ especially for fields stable over time
(Warnecke et al, 2018).

An illustration of local SVD results from a local Cartesian MHD simulation
with the NIRVANA code are shown in the leftmost set of Fig. 22, where
components of the $\alpha$ tensor and the most relevant components of the
turbulent resistivity tensor are shown as a function of the vertical height
from the midplane. The overall sign convention and peak values of the
$\alpha_{yy}$ component are rather similar to those found by Korpi (1999),
although the profiles are quite complex and there are even sign reversals as
functions of height, and strong gradients near the boundaries. The profiles of
the diagonal turbulent resistivity tensor components are also complex, and
unexpected negative values occur throughout the domain, especially near the
top of it.

Recently, Bendre et al (2023) developed a novel method, called iterative
removal of sources (IROS), for extracting the turbulent transport
coefficients. This method shares many similarities with the SVD method, while
the main differences are the iterative fitting algorithm, which incrementally
refines the estimates of the turbulent transport coefficients, and the
possibility to define prior physical constraints (mIROS scheme). The latter
property significantly decreases the errors, and the processing of extremely
noisy data becomes possible. The first results of the IROS method are
illustrated in Fig. 22 leftmost set, and compared to the SVD method. As is
evident, these two methods give nearly identical results.

An important approach to measuring the turbulent transport coefficients is the
TFM (Schrinner et al 2005, Schrinner et al 2007, recently reviewed in detail
by Brandenburg 2018). In short, this method is based on the exertion of
linearly independent test fields to the system that are passive in the sense
of not interfering with it, but act as a diagnostic tool. The action of the
turbulent flow on the test fields is monitored, and the equation for the
fluctuating fields is directly solved for each test field. Turbulent emf is
approximated as local and instantaneous and its expansion is truncated after
first order of the derivatives of $\overline{{\bm{B}}}$(Eq. (10)). This
approximation is again used to solve for the turbulent transport coefficients.
By directly solving for the fluctuating fields, many detrimental assumptions
of the FOSA/SOCA are avoided. The local and instantaneous formulation imposes
the restriction of capturing only effects of large-scale, slowly varying in
time MFs. This restriction can, however, be relaxed by considering a
formulation non-local in space and time (see, e.g. Rheinhardt and Brandenburg,
2012). Typical results from QKTFM in the kinematic stage of a large-scale,
dynamo-active, MHD model are shown in Fig. 22, rightmost set. Similar
kinematic analysis on MHD models with varying parameters has been performed by
Gressel et al (2008, 2008), and Bendre et al (2015).

In comparison to SVD, the profiles from QKTFM are overall somewhat smoother,
while the correspondence of the magnitudes and vertical profiles of the
$\alpha$ coefficients between the methods is satisfactory. The most notable
difference is in the magnitude of the turbulent resistivity components. Both
diagonal elements are roughly an order of magnitude larger with the QKTFM than
with SVD. The off-diagonal component $\eta_{xy}$ is consistent with zero from
SVD, while it has a clearly positive value from the QKTFM. On the other hand,
$\eta_{yx}$ is consistently non-zero and positive from SVD, while consistent
with zero from QKTFM. This does rule out the existence of the shear-current
dynamo effect in galactic discs by both methods, the existence of which would
require $\eta_{yx}$ being negative(see, e.g., Rogachevskii and Kleeorin,
2003).

The quenching of the turbulent transport coefficients in the non-linear
regime, due to the presence of strong large-scale magnetic field, has also
been studied with the QKTFM by Gressel et al (2013a). Finally, Gressel and
Elstner (2020) studied the importance of non-local and non-instantaneous
effects in similar models, and found out that especially the former are
important for galactic LSD. The $\alpha$ and $\eta$ coefficients were found to
reduce when the scale of the test fields is made smaller, agreeing with
results from simpler systems (e.g. Brandenburg et al, 2008). A surprising
finding was that the turbulent pumping, generally directed towards the disc
plane and aiding LSD, was becoming stronger with decreasing scale. This
behavior is yet to be fully understood.

Figure 22: Derived turbulent transport coefficients from three different
methods: (left) local SVD (red, errors in orange) and IROS (blue dashed,
errors in green) results from Bendre et al (2023); (right) QKTFM results from
Bendre et al (2020). The DNS-like MHD model from which these coefficients are
derived is Model Q in Bendre et al (2015), a local, Cartesian, SN-driven
simulation with $\Omega=4~{}\Omega_{\rm Sn}$ and $\sigma=0.25~{}\sigma_{\rm
Sn}$. The model develops large-scale dynamo action with oscillatory behavior
in the kinematic stage, and dynamically significant, but stable in time,
magnetic fields in the non-linear stage. The turbulent transport coefficients
are derived from the kinematic stage.

Until recently, DNS-like models to study galactic LSD have hovered at
resolutions too coarse to properly capture SSD simultaneously. Such models are
now emerging (Gent et al, 2023). In these models dynamically significant LSD
is excited despite the co-existence of fast-growing SSD. This might not be so
surprising, as an SSD alone in a multiphase, highly compressible fluid is
weak, producing fluctuating fields with an energy of only a few percent of the
kinetic energy in the turbulence (Gent et al, 2023). Also, these systems might
exhibit much larger helicity fluxes than the simpler systems studied
previously (see Sect. 4.2).

In this situation, determining the turbulent transport coefficients with the
QKTFM is no longer straightforwardly possible, as the SSD generates background
magnetic fluctuations that are not accounted for in the method. Viable
extensions (Rheinhardt and Brandenburg, 2010; Käpylä et al, 2020, 2022) have
appeared, but even further extensions, such as relaxing isothermality, are
required to deal with the galactic models. Furthermore, only the kinematic
version of the compressible test-field method could be validated for very
strong magnetic fields (Käpylä et al, 2022). Also, the required computational
resources for the extraction of the turbulent transport coefficients
significantly exceeds that for QKTFM. Evidently, SVD and IROS methods would
not suffer from any of the discussed limitations, although the better noise-
handling capacity of mIROS may turn out to be necessary.

Although the results of the different methods seem to be in rough agreement,
it is still necessary to ask the question of which of the methods is the most
reliable? The fitting-based methods SVD and IROS seem to agree well with each
other, but they are both based on the use of the MF of the MHD simulation
itself in the fitting procedure. In contrast, the test-field-based methods use
a set of linearly independent and well defined test fields. In the fitting
approach, the MF generated by the system could be degenerate, while the test-
field approach guarantees that such a situation cannot occur, and the
coefficients can be unambiguously determined. Also, the TFMs are rigorously
tested against analytical solutions such as various incarnations of the
Roberts flow (see, e.g., Rheinhardt and Brandenburg, 2010), while such
verifications have not been performed on the fitting-based methods. In Gressel
and Elstner (2020) the authors propose that the scale-dependence of the
transport coefficients could provide an explanation to the discrepancies
between SVD and QKTFM, but this appears unlikely. Namely, the SVD approach
uses the DNS-generated MF scale for the determination of both the $\alpha$ and
$\eta$ coefficients, which both would become reduced, if the proposed
explanation would be valid. Only the $\eta$ coefficients, however, show marked
reduction.

The best way to test these discrepancies would be to compare the evolution of
MF models using the measured coefficients to the actual MHD models. Curiously,
this has been done with the SVD coefficients by Bendre et al (2020), showing
the two solutions to be in excellent agreement in the oscillatory kinematic
regime. The period of dynamo cycles is known to directly depend on the
turbulent resistivity (see, e.g., Shukurov and Subramanian, 2021). Then one
would expect to see a discrepancy between the QKTFM and MF results. This
remains to be confirmed, and conclusions drawn accordingly.

The ultimate goal of the efforts to measure and describe the turbulent dynamo
coefficients relevant for the LSD in galaxies is to build computationally more
efficient global modelling frameworks. In essence, the goal is to build
explicit large eddy simulations, where also the backscatter (or inverse
cascade from small scales to larger ones) is captured in as much detail as
possible. This has been attempted either in the MF framework (Gressel et al,
2013b) or in a simpler context for cosmological modelling (Liu et al, 2022).
The most promising explicit large eddy frameworks exist only for simpler
turbulence models (see, e.g., Grete et al, 2017).

###  4.2 Measuring helicity fluxes

Helicity fluxes have also been measured from various types of MHD models. Such
work is vital, as the final proof that helicity fluxes actually alleviate
catastrophic quenching must come from full MHD models. The measurements
necessarily involve many simulations with varying ${\rm Re}_{\rm M}$, reaching
out to as high ${\rm Re}_{\rm M}$ as possible. This is demanding, so models
remain limited to simplified, isotropically forced systems.

The rate of change of the mean helicity density of the small-scale field reads

$\frac{\partial}{\partial t}\overline{{\bf
a^{\prime}}\cdot{\bm{b}^{\prime}}}=-2\mbox{\boldmath${\cal
E}$}{}\cdot\overline{{\bm{B}}}-2\eta\mu_{0}\overline{{\bf