# A three-phase amplification of the cosmic magnetic field in galaxies

Sergio Martin-Alvarez, Julien Devriendt, Adrianne Slyz, Romain Teyssier

► To cite this version:

Sergio Martin-Alvarez, Julien Devriendt, Adrianne Slyz, Romain Teyssier. A three-phase amplification of the cosmic magnetic field in galaxies. Monthly Notices of the Royal Astronomical Society, 2018, 479 (3), pp.3343-3365. 10.1093/mnras/sty1623. hal-01833727

HAL Id: hal-01833727

https://hal.science/hal-01833727v1

Submitted on 4 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A three-phase amplification of the cosmic magnetic field in galaxies

Sergio Martin-Alvarez,<fnref n="1" />* Julien Devriendt,<fnref n="1" />,<fnref n="2" /> Adrianne Slyz<fnref n="1" /> and Romain Teyssier<fnref n="3" />

Accepted 2018 June 14. Received 2018 June 14; in original form 2018 May 9

## ABSTRACT

Arguably the main challenge of galactic magnetism studies is to explain how the interstellar medium of galaxies reaches energetic equipartition despite the extremely weak cosmic primordial magnetic fields that are originally predicted to thread the inter-galactic medium. Previous numerical studies of isolated galaxies suggest that a fast dynamo amplification might suffice to brgalaxyidge the gap spanning many orders of magnitude in strength between the weak early Universe magnetic fields and the ones observed in high-redshift galaxies. To better understand their evolution in the cosmological context of hierarchical galaxy growth, we probe the amplification process undergone by the cosmic magnetic field within a spiral galaxy to unprecedented accuracy by means of a suite of constrained transport magnetohydrodynamical adaptive mesh refinement cosmological zoom simulations with different stellar feedback prescriptions. A galactic turbulent dynamo is found to be naturally excited in this cosmological environment, being responsible for most of the amplification of the magnetic energy. Indeed, we find that the magnetic energy spectra of simulated galaxies display telltale inverse cascades. Overall, the amplification process can be divided in three main phases, which are related to different physical mechanisms driving galaxy evolution: an initial collapse phase, an accretion-driven phase, and a feedback-driven phase. Whilst different feedback models affect the magnetic field amplification differently, all tested models prove to be sub-dominant at early epochs, before the feedback-driven phase is reached. Thus, the three-phase evolution paradigm is found to be quite robust vis-à-vis feedback prescriptions.

**Key words:** MHD – turbulence – methods: numerical – galaxies: formation – galaxies: magnetic fields – galaxies: spiral.

## 1 INTRODUCTION

Magnetic fields are one of the principal components of the Universe, with a presence that can be traced over an incredibly wide range of scales, from planets to clusters of galaxies and beyond. Whilst magnetism is not expected to have a considerable impact on the formation process of the large-scale structure (LSS) of the Universe, at galactic scales magnetic energy is believed to be in equipartition with the thermal energy in the interstellar medium (ISM; Wolfe, Lanzetta & Oren 1992; Bernet et al. 2008; Beck & Wielebinski 2013; Mulcahy et al. 2014). Therefore, magnetic fields on these scales could impact gas dynamics, and thus emerge as an essential component in the modelling of galaxy formation and evolution. Moreover, observations hint that magnetic fields are critical for several other processes in galaxies. For instance, they are crucial ingredients in the propagation of cosmic rays (Dubois & Commerçon

2016; Alves Batista et al. 2017; Wittor, Vazza & Brüggen 2017) and the generation of jets, not to mention key contributors to the processes of star formation (Li et al. 2009; Hennebelle & Iffrig 2014; Hull et al. 2017) and feedback, both from stellar (Chen et al. 2017; Ntormousi et al. 2017) and active galactic nuclei (Bambic, Morsony & Reynolds 2018).

State-of-the-art numerical methods are able to model the overall structure of the LSS and the properties of individual galaxies and galaxy clusters with remarkable accuracy (e.g. Anglés-Alcázar et al. 2014; Hopkins et al. 2014; Vogelsberger et al. 2014; Santos-Santos et al. 2016; Grisdale et al. 2017; Kaviraj et al. 2017). However, amongst the commonly missing ingredients in simulations of galaxies, magnetic fields stand out for their importance. This can be attributed to their intrinsically three-dimensional divergenceless vectorial nature that significantly increases the complexity of their correct numerical implementation and the daunting dynamical range required to capture their amplification. In an attempt to remedy the situation, various studies have recently explored the issue of the role played by magnetism both in galaxies and in the

<footnote n="1">¹Subdepartment of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="2">²Université de Lyon, Université Lyon 1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="3">³Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland</footnote>

<footnote>* E-mail: sergio.martin@physics.ox.ac.uk</footnote>
LSS (Dubois & Teyssier 2010; Beck et al. 2013; Gressel, Elstner & Ziegler 2013; Hennebelle & Iffrig 2014; Pakmor, Marinacci & Springel 2014; Vazza et al. 2014, 2015; Gheller et al. 2016; Marinacci & Vogelsberger 2016; Rieder & Teyssier 2016; Butsky et al. 2017; Su et al. 2017). However, despite these laudable efforts, key questions remain unanswered.

Arguably the most intriguing puzzle about magnetism is the question of its primordial origin, and how its amplitude evolves to reach present-day values. The exact formation time and mechanism of magnetic fields in the Universe is unclear (see Widrow 2002, for a review of different magnetogenesis scenarios). Theoretical models suggest the formation of weak fields in the very early Universe, typically during the inflationary phase or cosmological phase transitions (Hogan 1983; Kolb & Turner 1990; Ratra 1992). Various astrophysical plasma seeding processes have also been proposed (Biermann 1950; Schlickeiser 2012), suggesting the generation of a magnetic field with strength of the order of $B_0 \lesssim 10^{-19}$ Gauss. These field values are extremely weak, differing by many orders of magnitude from the micro-Gauss fields observed in galaxies and galaxy clusters (Beck & Wielebinski 2013; Mulcahy et al. 2014), and arguably even voids. How is such an enormous gap between primordial and current magnetic fields bridged? Present day averaged cosmic magnetic fields are constrained by observations that provide both lower and upper limits ($\sim 10^{-9} - 10^{-16}$ G; Neronov & Vovk 2010; Taylor, Vovk & Neronov 2011; Planck Collaboration XIX 2016), although the validity of the lower limit is still debated (e.g. Broderick, Chang & Pfrommer 2012). In any case, the range of values compatible with these observations remains wide, and as such does not provide stringent constraints on the generation mechanism.

For the specific case of galaxies, the gap between primordial and observed field values is the largest. The strength of ordered magnetic fields is typically found to be of the order of $10^{-5} - 10^{-6}$ G (Beck 2015; Mulcahy, Beck & Heald 2017). Therefore, extremely effective dynamo mechanisms are required to amplify the primordial fields, especially since equipartition fields are already reported at high redshifts ($z \simeq 2.0$ Wolfe et al. 1992; Bernet et al. 2008). As a consequence, any dynamo mechanism proposed as a solution needs to be able to amplify magnetic fields during the early stages of galaxy formation and evolution, as well as sustain them to the present day. An alternative scenario to the primordial origin of cosmic magnetism, stipulates that astrophysical sources be the progenitors of magnetic fields. These can range from stellar dynamos to active galactic nuclei (Rees 1987; Vazza et al. 2017) and supernovae (SNe) (Beck et al. 2013; Butsky et al. 2017). However, this scenario faces the challenge of how to stretch spatially localized fields to galactic scales and transfer them to the intergalactic medium and LSS.

Various fast and powerful amplification mechanisms able to strengthen the magnetic field on short time-scales have been invoked to address this problem. Some are the $\alpha - \Omega$ dynamo (Wang & Abel 2009; Dubois & Teyssier 2010) and the turbulent dynamo (Pakmor et al. 2014; Vazza et al. 2014; Rieder & Teyssier 2016; Pakmor et al. 2017; Rieder & Teyssier 2017a; Vazza et al. 2018). A magneto-rotational instability-driven dynamo is another possible candidate (Gressel et al. 2013; Machida et al. 2013), once a field strong enough to excite it has developed (Kitchatinov & Rüdiger 2004). Due to the complexity of the problem, numerical studies have become an important tool to probe different dynamos and seeding mechanisms. In fact, most have already been tested in plasma simulations (e.g. Haugen, Brandenburg & Dobler 2004a; Haugen, Brandenburg & Mee 2004b; Schekochihin et al. 2004;

Teyssier, Fromang & Dormy 2006; Schober et al. 2018) and simulations of individual, isolated galaxies (Dubois & Teyssier 2010; Pakmor & Springel 2013; Rieder & Teyssier 2016). Plasma simulations suggest that a turbulent dynamo is efficient enough to amplify the magnetic energy up to equipartition (Schekochihin et al. 2002; Federrath et al. 2011a, 2014; Schober et al. 2015), but such simulations are restricted to scales several orders of magnitude below that of a galaxy. Furthermore, even these simulations are restricted to a limited set of scales, being unable to simultaneously fully resolve both the inertial and sub-viscous ranges (see Kinney et al. 2000; Schekochihin et al. 2002). Meanwhile, in isolated galaxies simulations, the amplification exhibits moderate growth rates. The main reason is that successfully modelling turbulent dynamo amplification requires extreme spatial resolution to capture a significant fraction of the inertial range and a considerable level of galactic turbulence. The task becomes more arduous when this picture is studied in a cosmological context, where galaxies are born intrinsically small at high redshifts and there is a further added cost for also simulating the cosmological environment. However, rapid magnetic amplification has recently been reported in simulations of both isolated (Rieder & Teyssier 2016) and cosmological (Rieder & Teyssier 2017b) dwarf galaxies.

In this paper, we follow the evolution of the magnetic energy in a Milky Way-like galaxy forming in an explicit cosmological environment with high spatial resolution (10 pc), sufficient to trigger the turbulent dynamo. As far as we are aware, this is the first time this experiment is performed at such a resolution whilst employing a code that guarantees negligible divergence of the magnetic field at all times. Our work extends the study of magnetic field amplification occurring in isolated spiral galaxies (Rieder & Teyssier 2016) by accounting for processes such as galaxy mergers and gas accretion. We find that a turbulent dynamo also leads to growth of the magnetic energy, but that the amplification is modulated by environmental effects, thus undergoing different phases where magnetic growth occurs at different rates. Repeating our simulation, in turn, without stellar feedback and with a stellar feedback prescription different from our fiducial one, we establish that feedback only dominates the amplification process during the later stages of galaxy evolution, when gas accretion and mergers have significantly subsided.

The structure of the manuscript is as follows: the numerical setup is described in detail in Section 2. Results are described in Section 3, where three main epochs of amplification are identified: collapse (Section 3.1), accretion-driven (Section 3.2), and feedback-driven (Section 3.3) amplifications. Section 4 discusses the robustness of this three-phase amplification and reviews implementation-dependent aspects of our findings. Our concluding remarks are given in Section 5. We add two appendices addressing numerical resolution (Appendix A), and presenting the method we use to calculate the Fourier power spectra of a disc galaxy embedded within an explicit cosmological context (Appendix B).

## 2 NUMERICAL METHODS

To run the cosmological zoom-in simulations presented in this paper, we employ a modified version of the publicly available magnetohydrodynamical (MHD) code RAMSES (Teyssier 2002) that accounts for the magnetic field dilution caused by the expansion of the universe. RAMSES couples a tree-based Adaptive Mesh Refinement (AMR) Eulerian treatment of the gas with an N-body treatment for the dark matter and stellar components. In this section, we briefly introduce the simulation and discuss various numerical aspects relevant for assessing the robustness of our results.
## 2.1 Simulations

### 2.1.1 Initial conditions, resolution, and cooling/heating

The initial conditions for all simulations are those of the NUT suite (Powell, Slyz & Devriendt 2011) and consist of a cosmological box of $9 h^{-1}$ Mpc comoving side, with a $3 h^{-1}$ Mpc comoving diameter zoom-in spherical region. In this sphere, the code is allowed to refine cells in a quasi-Lagragian way (so that any cell contains roughly the same amount of mass) up to a maximum physical spatial resolution of 10 pc. A $M_{vir} \sim 5 \times 10^{11}$ M$_{\odot}$ DM halo is formed in this region by redshift $z = 0$. The dark matter and star particle minimal masses are $M_{DM} \sim 5 \times 10^4$ M$_{\odot}$ and $M_* \sim 5 \times 10^3$ M$_{\odot}$, respectively. Initial conditions are set at redshift $z = 500$, with relevant cosmological parameters chosen according to WMAP5 values (Dunkley et al. 2009). The simulations also include a UV background to take into account inter-galactic medium photoionization effects (Haardt & Madau 1996) (re-ionization of the simulated volume is assumed to happen instantaneously at $z = 10$) and metallicity dependent gas cooling both at temperatures above $10^4$ K (Sutherland & Dopita 1993) and below (Rosen & Bregman 1995). The gas is assumed to be monatomic with specific heat ratio $\gamma = 5/3$ and follows an ideal equation of state.

### 2.1.2 Star formation

Star formation is generally modelled in simulations following a Kennicutt–Schmidt law (Schmidt 1959) with a constant efficiency $\epsilon$. As long as it is employed on large (kpc) scales, the constant efficiency approximation provides a fair description of the observational data. However, at smaller (pc) scales, it is expected that $\epsilon$ will depend on local ISM properties (Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012). This motivates a thermo-turbulent star formation model, where stars form with a high efficiency ($\epsilon \ge 0.1$) in regions where gravity overcomes the turbulent, magnetic, and thermal pressure support. Such a spatially varying star formation efficiency prescription has been introduced in Kimm et al. (2017), Trebitsch et al. (2017), and Mitchell et al. (2018) and is presented in more detail together with an analysis of its impact in Devriendt et al. (in preparation).

All the runs presented in this work spawn star particles in the highest resolution cells (Rasera & Teyssier 2006) in line with this calculation of the local star formation efficiency. The resulting stellar particles can then generate stellar feedback following one of the two different stellar feedback prescriptions outlined below. A Kroupa initial mass function (Kroupa 2001, IMF) is assumed for the stellar populations. Each SN returns a mass fraction of gas $\eta_{SN} = 0.213$ and newly formed metals $\eta_{metals} = 0.075$ back to the ISM that are then advected by the code as passive scalars.

### 2.1.3 Mechanical feedback

The first feedback prescription (Mech) is a mechanical SN feedback model presented in Kimm & Cen (2014) and Kimm et al. (2015). For a given numerical resolution, this model calculates the proper amount of momentum injection by a Sedov–Taylor blast wave that depends on the stage of the explosion. This momentum – along with energy and mass – depends on various local parameters (e.g. number of occurring SNe, metallicity, and surrounding density of the ISM) and is anisotropically deposited in the cell containing the SN and its neighbours. This approach considerably reduces artificial cooling (Marri & White 2003; Slyz et al. 2005). See Kimm

& Cen (2014) and Kimm et al. (2015) for a detailed explanation of the implementation and Rosdahl et al. (2017) for an extended comparison to other popular SN feedback implementations.

Note that due to its better intrinsic capability to capture the anisotropic expansion of shock waves through local inhomogeneities, this model is expected to better describe shock turbulence. Indeed, it should drive baroclinic vorticity more efficiently and result in higher solenoidal turbulence at small scales. Solenoidal forcing in shocks, in turn, fosters magnetic amplification (Ji et al. 2016), without the need for large amounts of energy injection.

### 2.1.4 Thermal and radiative feedback

The second feedback prescription (RdTh) includes both SNe feedback and stellar radiation. Energy injection from SNe alone might not suffice to generate galactic outflows. Stellar radiation has been invoked to supplement the driving power (Murray, Quataert & Thompson 2010). This model implements such a radiatively boosted stellar feedback, which couples to a standard thermal energy injection by SNe, rendering this latter more efficient.

More specifically, all star particles in the simulation emit UV photons. Dust in the ISM absorbs these photons and re-emits the absorbed energy at infrared wavelengths. Infrared photons then transfer momentum to the gas and drive galactic outflows (Murray et al. 2010; Roškar et al. 2014). The efficiency of such a momentum transfer is controlled by the infrared dust opacity $\kappa_{IR}$ in our implementation. We set this free parameter to a value of $\kappa_{IR} = 20 \text{ cm}^2 \text{ g}^{-1}$, which mimics an opaque ISM, very efficient at absorbing infrared energy (Semenov et al. 2003) and generates large outflows. Detail of the implementation can be found in Roškar et al. (2014). SN explosions are modelled as a simple injection of mass, metals, and thermal energy into their host cell.

### 2.1.5 The simulation set

The simulation set used in this work consists of three principal runs, B-Mech, B-NoFb, and B-RdTh that are distinguished by their stellar feedback prescriptions, and four supplementary runs performed at lower resolution: B-MR20, B-MR40, B-MR80, and B-MR160. The B-Mech and B-MRX runs all employ the mechanical SN feedback implementation (Section 2.1.3). Runs labelled MR20, MR40, MR80, and MR160 have a maximal spatial resolution lower than that of their fiducial counterpart B-Mech by factors of 2, 4, 8, and 16, respectively. The B-NoFb run has no stellar feedback, and the B-RdTh run employs the thermal SN prescription and radiative stellar feedback presented in Section 2.1.4. In any figure where they are inter-compared, the B-Mech, B-NoFb, and B-RdTh simulations are colour-coded blue, red, and green, respectively. The main characteristics of the simulations are summarized in Table 1. Supplementary runs performed with a view to estimate the impact of resolution on our results are mostly discussed in Appendix A. All simulations have been evolved down to $z = 2$.

## 2.2 Magnetic field specifics

### 2.2.1 Constrained transport

The MHD solver of RAMSES is based on the constrained transport (CT; Evans & Hawley 1988) scheme (Fromang, Hennebelle &
<table><caption>**Table 1.** Simulations used in the manuscript. Co-moving seed magnetic field $B_0$, maximum resolution $\Delta x_m$, and stellar feedback model are indicated for each run.</caption><thead><tr><th>Simulation</th><th><i>B</i><sub>0</sub> (G)</th><th>Δ<i>x</i><sub>m</sub> (pc)</th><th>Feedback</th></tr></thead><tbody><tr><td>B-Mech</td><td>3 × 10<sup>-20</sup></td><td>10</td><td>Mech<br>(Section 2.1.3)</td></tr><tr><td>B-NoFb</td><td>3 × 10<sup>-20</sup></td><td>10</td><td><i>x</i></td></tr><tr><td>B-RdTh</td><td>3 × 10<sup>-20</sup></td><td>10</td><td>RdTh<br>(Section 2.1.4)</td></tr><tr><td>B-MR20</td><td>3 × 10<sup>-20</sup></td><td>20</td><td>Mech</td></tr><tr><td>B-MR40</td><td>3 × 10<sup>-20</sup></td><td>40</td><td>Mech</td></tr><tr><td>B-MR80</td><td>3 × 10<sup>-20</sup></td><td>80</td><td>Mech</td></tr><tr><td>B-MR160</td><td>3 × 10<sup>-20</sup></td><td>160</td><td>Mech</td></tr></tbody></table>

Teyssier 2006; Teyssier et al. 2006). It solves the induction equation

$$ \frac{\partial \vec{B}}{\partial t} = \vec{\nabla} \times (\vec{v} \times \vec{B}) + \eta \vec{\nabla}^2 \vec{B}, \quad (1) $$

(with $\vec{v}$ gas velocity; $\eta$ magnetic diffusivity) to evolve the magnetic field $\vec{B}$ in time, guaranteeing that the solenoidal constraint

$$ \vec{\nabla} \cdot \vec{B} = 0, \quad (2) $$

holds to numerical precision. As magnetic diffusivity is negligible on galactic scales, we set $\eta = 0$, therefore all magnetic diffusivity is exclusively due to numerical effects. The divergence free aspect of the scheme is especially relevant at shocks and in regions where turbulent motions dominate, as magnetic field tangling may generate large artificial divergences (Balsara & Kim 2004). Non-constrained methods or divergence cleaning algorithms can typically display a divergence per cell ($\vec{\nabla} \cdot \vec{B}$) $dx_{cell}$ of the order of the strength of the local field $|\vec{B}|$, with $dx_{cell}$ the cell size in such regions (Wang & Abel 2009; Kotarba et al. 2011; Pakmor & Springel 2013). The main problem is that these highly turbulent regions are also the places where most of the magnetic field amplification is expected to occur if a small-scale, fast dynamo is the culprit. For non-constrained methods, the magnetic mono-polar contribution is preserved once generated, casting serious doubts on their ability to reliably capture any physical amplification mechanism. Divergence cleaning codes base their solenoidal constraint fulfilment on a rapid diffusion of the spurious numerical contribution to the divergence (see Powell et al. 1999; Dedner et al. 2002). However, by the time such a diffusion has brought the divergence to an acceptably low level, the mono-polar contribution has already likely had an impact on both the amplification and the structure of the field [see Mocz et al. (2016), where divergence cleaning and CT schemes are compared for an isolated spiral galaxy simulation]. We refer the reader to Hopkins & Raives (2016) for a comparison of how different numerical methods perform in a variety of standard MHD test cases.

To quantify the accuracy of our CT scheme in the context of cosmological AMR simulations of galaxy formation, we show in Fig. 1 the values of the divergence to magnetic field ratio ($\vec{\nabla} \cdot \vec{B}$) $dx_{cell}/|\vec{B}|$ as a function of redshift. All three main runs presented in Section 2.1 are shown. They display insignificant median values (dark solid lines) of this ratio, at a level consistent with numerical precision. Interquartile ranges are narrow and always found within ±1 dex above and below the median. Dashed lines indicate the 3$\sigma$ width of the entire distribution of AMR cells: they are always below a ratio of $10^{-10}$. Finally, large circles, triangles, and squares represent the maximal values ever reached by a single cell. These are labelled with the level of refinement of the cell, with 1 corresponding to maximal resolution and higher numbers indicating ever coarser

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Magnetic field divergence to total magnetic field strength ratio ($\vec{\nabla} \cdot \vec{B}$) $dx_{cell}/|\vec{B}|$ per cell with length $dx_{cell}$, as a function of redshift. Median (solid lines) and 3$\sigma$ (dashed lines) values of the AMR cell distribution are shown for the entire simulation volume. Squares (B-Mech), circles (B-NoFb), and triangles (B-RdTh) indicate maximal values reached in any individual cell across the whole simulation volume, and are labelled with the level of resolution of that cell. Number 1 corresponds to the highest resolution level, and larger numbers indicate ever lower resolution levels. Note that even those individual maxima always lie below 2–3 orders of magnitude below the field strength (grey line).</figcaption></figure>

refinement levels. It can be seen that at all times, and for all cells, the divergence to magnetic field strength ratio occurring remains below per cent levels.

### 2.2.2 Resolution requirements

It has also been argued by Federrath et al. (2011b) that self-gravitating MHD simulations require a minimum of 30 resolution elements per Jeans length to achieve significant dynamo action. Due to the multiphase nature of the ISM and the different impact of each phase on the magnetic energy evolution (Evirgen et al. 2017), attaining this resolution across as much of the galaxy as possible is of key importance. For spiral galaxies, such a feat becomes rapidly challenging given how thin their gas disc become. Fig. 2 displays the typical resolution reached in the disc of our simulated galaxies. The half-thickness $h_s$ estimated from fitting an exponential density profile $\rho(z) \propto e^{-|z|/h_s}$ to the galaxy gas disc in the range of redshift probed by the simulations ($2 < z < 4$), is typically $h_s \sim 0.6$ kpc for the runs with feedback, and $h_s \sim 0.4$ kpc for the run without feedback. The vertical dimension of the galaxy disc thus contains approximately 70 resolution elements already for the second highest level of resolution. For a thickness of the disc comparable to its Jeans length, this quantity is well above the threshold suggested by Federrath et al. (2011b). In Appendix A, we also discuss how the magnetic Reynolds number requirements to capture a turbulent dynamo are met in our simulations. As far as we are aware, this work is the first to combine this kind of resolution with an accurate CT scheme for MHD in cosmological simulations of disc galaxies.

### 2.2.3 Initial conditions

Finally, due to the absence of battery terms in equation (1), magnetic fields have to be initially seeded in the simulations. Given the magnetic seed strengths employed, the spatial distribution of the initial magnetic seed has no dynamical effects in collapsing proto-galaxies. Marinacci et al. (2015) have also shown that the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Edge-on projection view of a (8 kpc)³ region centred on the B-Mech galaxy at z = 2. Colours display density-weighted cell sizes. Ne-glecting the maximal level of resolution (10 pc, in black) only found in the densest regions, the galactic disc is resolved with approximately ~70 resolution elements (20 pc, dark blue), well above the necessary threshold advocated by Federrath et al. (2011b) to achieve significant dynamo action.</figcaption></figure>

orientation of the initial seed has no impact on the field structure and its average properties within haloes. In virialized regions such as galaxy clusters and galaxies, all memory of the initial seed is also expected to be lost once enough turbulence and amplification have occurred (Dubois & Teyssier 2010). As a consequence, the use of a uniform initial field, both in strength and in orientation has become a standard amongst astrophysical simulations, and we follow suit by using a uniform field oriented along the z-axis of the box. All runs are seeded with an identical, weak magnetic field, assuming a typical Biermann battery seed strength (Biermann 1950). This corresponds to a co-moving magnetic field $B_{co} = 3 \times 10^{-20}$ G, which is the natural quantity evolved in the simulations. Indeed, we introduce a super-comoving magnetic field in RAMSES, closely following the implementation suggested by Martel & Shapiro (1998). Our value of $B_{co} = 3 \times 10^{-20}$ G corresponds to a physical field strength of $B_{ph} \sim 10^{-17}-10^{-18}$ G at $z \sim 40$, in keeping with classical Biermann battery estimates quoted by Pudritz & Silk (1989) and Widrow (2002). The selection of such low magnetic seeds ensure that our simulations remain in the kinematic regime of dynamo amplification even after the galaxy undergoes a significant initial compression. Therefore, we can focus on measuring the amount of magnetic turbulent dynamo amplification and whether it leads to saturation.

We stress, however, that implementing a different seeding (structure of the field) even with our weak primordial magnetic field strength might lead to a different amplification level during the early collapse phase. For instance, Rieder & Teyssier (2016) measure such an effect in their simulations of isolated galaxies. We defer a detailed investigation of this issue to future work.

## 2.3 Halo and galaxy identification

Any study of galaxies in cosmological simulations requires a clear discrimination between these objects and their environment. In our simulations, this is achieved by running the HaloMaker software described in Tweed et al. (2009) on the dark matter particle distribution. The resulting structures, obtained via a thresholding of the dark matter density field, are identified as dark matter haloes. Note that sub-haloes are also identified in this way using the AdaptaHOP algorithm (Aubert, Pichon & Colombi 2004). Within these dark matter haloes, the most massive bound baryonic objects are then detected and identified as galaxies.

Various parts of the analysis presented in this paper require an accurate identification of galaxy centres (e.g. measuring rotational velocities or generating density/velocity/energy profiles). The centres of our bound baryonic structures are calculated according to the shrinking sphere method proposed by Power et al. (2003) applied to gas and stars (when these latter are present). Finally, unless otherwise stated, all physical quantities are measured within a sphere centred on each galaxy and extending out to $0.2 r_{halo}$, where $r_{halo}$ stands for the virial radius of the halo.

## 3 RESULTS

In this section, we evaluate magnetic energy growth in simulated galaxies. As a preamble, we note that, by the end of the simulated period, i.e. $z=2$, magnetic fields in observed galaxies have already reached equipartition with other forms of energy (Bernet et al. 2008). Thus, we expect that the most significant part of magnetic energy amplification occurs sometime between the seeding of the primordial fields and this redshift.

In contrast to simulations of isolated galaxies, galaxies formed within an explicit cosmological context evolve following a hierarchical formation scenario. More specifically, Fig. 3 presents the gas density, temperature, magnetic energy density, and optical mock images for our galaxies at $z=2$. In reaching this epoch, galaxy properties have been sculpted by various environmental effects such as cold filamentary and hot diffuse modes of gas accretion and mergers. As galaxies grow, the evolution of their magnetic fields is influenced by such changing environment. We identify three different phases in the evolution of the galaxy magnetic field, from the beginning of the simulations to redshift $z=2$. These three phases are well correlated with different physical processes governing the evolution of the galaxy. They are the initial collapse of the proto-galaxy, an accretion phase during which the main factor affecting galaxy properties is the fuelling of its growth by strong cold flow accretion, and a feedback-dominated phase where stellar feedback finally overcomes the dwindling accretion. During each phase, we assess the effects of stellar feedback and cosmological environment processes.

All quantities in this paper are expressed in physical units (by opposition to co-moving) unless otherwise indicated.

### 3.1 Formation of the galaxy: the collapse phase

The growth of a galaxy in a cosmological cold dark matter simulation follows a hierarchical collapse scenario (e.g. White & Frenk 1991). Initial perturbations expand until they reach the turnaround point, when their magnetic energy reaches a minimum. At this stage, they decouple from the expansion and their collapse begins. The deepening gravitational potential builds up the galaxy by accumulating baryonic mass until stars form. In our simulations, we identify the end of the collapse phase ($t_{coll}$ or $z_{coll}$) with the moment
<figure><img src="image_5.png" /><figcaption>Figure 4. Gas (red), dark matter (yellow), and stellar (blue) masses within the galactic region for the B-Mech (top), B-NoFb (middle), and B-RdTh (bottom) runs. Vertical lines identify the end of the collapse (dotted) and the accretion (dot-dashed) phases.</figcaption></figure>

specific angular momentum to the dark matter enclosed within the virial radius of its host halo. Indeed, this yields $r_{gal} \sim 4 \times \lambda r_{halo}$ (equation 14 of Mo, Mao & White 1998) given the measured spin parameter of our halo, as defined in Bullock et al. (2001) remains close to $\lambda \sim 0.05$ at all times Kimm et al. (2011). The central region thus defined will be referred to as the galactic region in the remainder of the paper. Finally, at very early epochs, i.e. between the turnaround time when the density perturbation decouples from the expansion of the universe and the formation of the first virialized dark matter halo we are able to detect given our resolution, we use a fixed galactic region physical size $r_{gal} = 0.7$ kpc $\simeq 0.2 r_{halo}(z_{coll})$, centred on the co-moving position of the collapsed proto-galaxy at the end of collapse. We have checked that reasonable alterations to this somewhat arbitrary choices for size and centre of the galactic region do not significantly impact our results, as all density fields tend to become very rapidly uniform at these early epochs.

Changes in galaxy mass due to accretion, mergers, and outflows lead to ambiguities when comparing absolute energies within the galactic region. Similarly, energy densities are sensitive to the shape and volume of the galaxy, as magnetic lines are frozen in the fluid in ideal MHD. This is especially relevant for the study of disc galaxies. To circumvent this issue, we therefore resort to using specific energies

$$ \epsilon_x = \frac{E_x}{M_{\text{gas}}}, \qquad (3) $$

obtained by taking the ratio between total energies $E_x$ (not energy densities, $e_x$, although the two quantities are, of course, related by a simple volume factor) and gas mass $M_{gas}$ within the given region. This allows us to disentangle the amplification of the magnetic energy caused by the growth of the galaxy from that of an intrinsic dynamo effect.

During the initial phase of formation, the growth of the magnetic energy density in an idealized, isolated halo closely follows an isotropic adiabatic approximation as the gas collapses, so that $e_{mag} \propto \rho_g^{4/3}$ where $\rho_g$ is the gas density (e.g. Wang & Abel 2009; Dubois & Teyssier 2010; Rieder & Teyssier 2016). However, the cosmological formation of a galaxy involves a more complex accretion/collapse history that fosters higher turbulence levels than the isolated case (Elmegreen & Burkert 2010; Klessen & Hennebelle 2010). The relevance of accretion-driven turbulence increases when the growth of the gravitational potential is taken into account (Elmegreen & Burkert 2010). Fig. 5 displays the evolution of the specific magnetic energy $\epsilon_{mag}$ within the galactic region during the collapse phase and the beginning of the accretion phase. Note that specific magnetic energies are normalized to their respective (but only different by a few percent) values, $\epsilon_{mag, ta}$, within the galactic region at the time of turnaround for each of the runs. As expected, the collapse times of haloes and thus the magnetic amplification levels reached during collapse, are quite robust vis-à-vis different stellar feedback prescriptions, at the resolution reached in our simulations.

More quantitatively, the estimate of the collapse time

$$ t_{\text{coll}} = \frac{2}{3H_0(1 + z_{\text{coll}})^{3/2}}, \qquad (4) $$

where $z_{coll}$ is defined as the redshift at which a spherical top hat density perturbation reaches the linearly extrapolated critical overdensity contrast in an Einstein-de Sitter universe (e.g. Peebles 1994) can be compared with the onset of star formation in simulated galaxies. We find the value of the collapse time defined in this manner, $t_{coll} \sim 0.34$ Gyr ($z_{coll} \sim 13$) (and by extension the turnaround time $t_{ta} \simeq 0.5$ $t_{coll}$; $z_{ta} \sim 22$), to be in fair agreement with the measured galaxy formation time (compare the vertical dotted line on Fig. 5 to the earliest peak of the dashed curves). Specific magnetic energies measured in the galactic region at exactly $t_{coll}$ are used as reference in the rest of the paper and denoted $\epsilon_{mag, c}$. They provide an upper limit of the magnetic amplification expected from the collapse phase because the expansion of the universe ensures that magnetization of any gas subsequently fed to the halo will be lower. They are indicated by horizontal lines on Fig. 6 and spread out by a factor of 4 between runs (with values of $55 \epsilon_{mag, c}$ for B-Mech, 220 for B-NoFb, and B-RdTh in between), hinting that stellar feedback might play a role in reducing magnetic field amplification in the later stages of the collapse. Note that if we relax the collapse time criterion to some degree, allowing the halo in the different runs to collapse at slightly different times (which is bound to happen given the chaotic nature of the system) and identify $\epsilon_{mag, c}$ with the peak
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Specific magnetic energy growth during the collapse phase and the beginning of the accretion phase within the galactic region (inner $0.2r_{\text{halo}}$ of the dark matter halo). The three main runs B-Mech (blue), B-NoFb (red), and B-RdTh (green) are displayed. Dashed orange lines represent the expected growth of the specific magnetic energy if the collapse were to proceed adiabatically and isotropically. All runs are normalized to their value at the epoch of turnaround $t_{\text{ta}}$ (see the text for detail). The vertical dotted line marks the end of the collapse phase.</figcaption></figure>

in specific magnetic energy for each run (solid curves), the gap between runs is still present, even though it is reduced to a factor of 1.7 (130 $\epsilon_{\text{mag,ta}}$ for B-Mech, 220 for B-NoFb and B-RdTh in between, see Fig. 5).

From Fig. 5, it is clear that in the simulations, specific magnetic energies $\epsilon_{\text{mag}}$ grow by a factor of ~100–200 between the time of turnaround (beginning of the time axis) and the end of the collapse phase (vertical dotted lines). These values are about an order of magnitude above those derived in the analytical study of Lesch & Chiba (1995), who use the spherical top hat model to estimate that the magnetic field strength is amplified by $f_c \sim [r_{\text{ta}}/(0.2r_{\text{halo}})]^2 \sim 100$ from turnaround until disc formation. Indeed, this corresponds to a specific magnetic energy increase of $f_c^2 [(0.2r_{\text{halo}})/r_{\text{ta}}]^3 \sim 10$

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Specific magnetic energies growth (solid lines) in the galactic regions for the B-Mech (blue), B-NoFb (red), and B-RdTh (green) runs. Dashed lines show the specific magnetic energies at the end of the collapse phase $\epsilon_{\text{mag, C}}$ for each run. Vertical lines correspond to the beginning of the accretion (dotted) and feedback (dot-dashed) phases. All energies are normalized to their turnaround value.</figcaption></figure>

in the galactic region – the gas mass within the sphere of radius $r_{\text{ta}}$ is assumed to be conserved exactly in the analytic estimate – which matches quite well our own estimate of the evolution of the specific magnetic energy expected from pure isotropic adiabatic collapse (orange dashed curve in Fig. 5). This latter estimate is obtained assuming that each individual cell in the galactic region has undergone such a collapse from the time of turnaround, i.e., that its magnetic energy is given by

$$E_{\text{mag}}^{\text{cell}}(t) = \frac{dx_{\text{cell}}^3}{\frac{4}{3}\pi r_{\text{gal}}^3(t_{\text{ta}})} \times \langle E_{\text{mag}}(t_{\text{ta}}) \rangle \left( \frac{\rho_{\text{g}}^{\text{cell}}(t)}{\langle \rho_{\text{g}}(t_{\text{ta}}) \rangle} \right)^{4/3}, \quad (5)$$

where $\langle \rangle$ stands for averages over the whole galactic region at turnaround. One then obtains the adiabatic isotropic estimate for the specific magnetic energy of the galactic region at time $t$ by adding the contributions of all the cells enclosed within the region at that time and dividing it by the total amount of gas present in the region.

We emphasize that such an estimate assumes not only that the collapse is adiabatic and isotropic but also that no decrease of the field occurs due to cosmological expansion. In reality, the strength of the magnetic field should be lower post-collapse because (a) the expanding Universe reduces the magnetization of material accreted later ($B \propto a^{-2}$), and (b) even though the initial collapse might be close to isotropic, material will subsequently be fed to the galaxy through filamentary accretion and mergers. Filamentary accretion is expected to bring in material with lower magnetization due to the simple spatial trajectories followed by the infalling gas, given the low amount of cold streams turbulence expected at the mass range of our halo (Mandelker et al. 2017). The impact of mergers on magnetic amplification is more complex and its treatment is deferred until Section 3.2. In any case, shortly before the initial collapse phase ends ($t \gtrsim 0.25$ Gyr), this crude estimate of isotropic and adiabatic amplification starts to systematically underestimate the measured specific magnetic energies by a significant amount (solid curves in Fig. 5). In practice, it is found to oscillate at a level between $\epsilon_{\text{mag}} \sim 10-20 \epsilon_{\text{mag,ta}}$ throughout the simulations (i.e. until $z=2$), and is thus more representative of a configuration where the galaxy accretes material with constant specific magnetic energy
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time evolution of specific magnetic energies in a thin shell around the galaxy halo ($r_{halo} < r < 1.05 r_{halo}$, dashed lines) for the three main runs compared to the same quantities within the galactic region (solid lines, identical to Fig. 6). Specific magnetic energies within the outer halo shell decrease over time as gas accretes. Note the contribution of outflows driven by stellar feedback to $\epsilon_{mag}$ in the shell.</figcaption></figure>

throughout its life time, i.e. without any kind of amplification or decay of the magnetic field beyond the initial compression.

In summary, our simulations clearly show (Fig. 5) that by the end of the collapse phase, the level of turbulence induced in the proto-galaxy has already amplified the specific magnetic energy within the galactic region by approximately an order of magnitude above the isotropic adiabatic collapse prediction. Such a significant deviation has also been reported for larger virialized objects such as galaxy clusters (e.g. Dolag et al. 2005; Dubois & Teyssier 2008; Marinacci et al. 2015). Smaller structures, however, generally still seem to follow the $\epsilon_{mag} \propto \rho_g^{4/3}$ scaling, even in the most recent work of Pakmor et al. (2017), which can probably be attributed to their lower spatial and mass resolution (based on the resolution study provided in Appendix A). Although not easily gathered by looking at Fig. 5, towards the end and immediately after the collapse phase ($0.25 \lesssim t \lesssim 0.5$ Gyr), the deviation from the isotropic case is mainly driven by turbulent accretion. Indeed, Fig. 6, where we plot the three runs together on the same panel, presents evidence that the more powerful the feedback prescription (the B-RdTh run being the most powerful), the lower the difference with the isotropic collapse estimate.

## 3.2 Turbulent amplification: the accretion phase

As previously mentioned, accretion of pristine gas from $t_{coll}$ onwards contributes to reducing specific magnetic energies due to its lower magnetization. As shown in Fig. 7, the specific magnetic energy measured within a thin spherical shell located at $r_{halo} < r < 1.05 r_{halo}$ is significantly lower than the specific magnetic energy within the galactic region and exhibits a tendency to slowly decline with time. The presence of stellar feedback somewhat alleviates this decay by pushing a fraction of the magnetic field amplified in the galactic region out in the halo via outflows. This explains why both the B-Mech and B-RdTh runs show specific magnetic energies in the outer halo shell that are systematically larger than those of the B-NoFb run. Note that the specific magnetic energy of the virialized intra-halo region (i.e. $r < r_{halo}$) would there-

fore decrease if the halo were only to accrete this pristine, weakly magnetized fresh gas.

In the accretion-dominated phase, extending from $z_{coll} \sim 13$ until $z_{acc} \sim 4$, the galaxy continuously grows in mass and size, principally through accretion of cold filamentary flows (Tillson et al. 2015). Once the collapse phase has ended and the average gas density of the galaxy is more or less fixed, growth of specific magnetic energy beyond $\epsilon_{mag}$ requires stretching the magnetic field through dynamo processes. In particular, a turbulent ISM can generate exponential amplification of the magnetic energy through a turbulent dynamo (Rieder & Teyssier 2016). The level of ISM turbulence is highly dependent on gas accretion, either directly (Elmegreen et al. 2009a; Klessen & Hennebelle 2010) or through accretion-fuelled processes such as gravitational instabilities (Elmegreen & Burkert 2010; Krumholz & Burkhart 2016). Even when turbulence is mainly attributed to stellar feedback processes, accretion is still credited with playing an increasingly relevant role, the higher the redshift becomes (Hopkins, Kereš & Murray 2013). Klessen & Hennebelle (2010) find that it mainly is the total mass of the galaxy that determines the importance of accretion-fuelled turbulence. The halo mass also sets whether cold accretion operates (Birnboim & Dekel 2003; Dekel & Birnboim 2006; Ocvirk, Pichon & Teyssier 2008). As galaxies and haloes accumulate mass, cold accretion progressively shuts down and the resulting turbulence rapidly decays.

In the absence of feedback, virialized hot halo gas can cool and fall onto the galaxy, but at the redshifts and halo masses discussed in this work, feedback is required to maintain the presence of a significant hot gas phase (Agertz et al. 2009a) as accretion onto the galaxy proceeds mainly along cold filaments (Tillson et al. 2015). However, feedback could also be reducing the budget of both absolute and specific energies available to drive turbulence through accretion by launching massive and disruptive galactic winds.

Bearing this possibly complex picture in mind, the evolution of turbulence within the galactic region is analysed by measuring in turn specific turbulent energies (Fig. 8), radial velocity profiles (Fig. 9), and turbulent kinetic energy spectra (Fig. 13).

All specific energies displayed in Fig. 8 have been normalized to $\epsilon_{ta} = 10^{10} \text{ cm}^2 \text{ s}^{-2}$ that represents the specific thermal energy of the galactic region at turnaround. The turbulent energy $E_{turb}$ is calculated by summing over the contribution of each cell enclosed in the galactic region:

$$E_{\text{turb}} = \sum_{i=1}^{N_{\text{cells}}} \frac{1}{2} m_i v_{\text{turb},i}^2, \qquad (6)$$

where $m_i$ is the gas mass, and the turbulent velocity of the gas cell is defined as $v_{\text{turb},i} = \sqrt{v_{r,i}^2 + (v_{t,i} - v_c(r))^2}$, with $v_{r,i}$ and $v_{t,i}$ the radial and tangential components of the gas velocity, respectively, and the bulk motion of the galactic region has been subtracted. To remove the circular velocity $v_c(r) = \sqrt{GM(r)/r}$, where $M(r)$ is the total mass within the sphere of radius $r$, two distinct procedures are employed depending on whether a disc with substantial rotational support is present. If the measured spin parameter of the galactic region, $\lambda = L/(\sqrt{2} M(r_{gal}) v_c(r_{gal}) r_{gal}) < 0.5$ [with $L$ the angular momentum as in Bullock et al. (2001)], $v_c(r)$ is subtracted from the tangential velocity (i.e. the vector formed by the two non-radial components of the velocity field) in each cell. After the formation of the disc, when $\lambda > 0.5$, $v_c$ is only subtracted from the toroidal velocity in a cylindrical coordinate system, where the z-axis is aligned with the gas angular momentum of the disc and goes through the centre of the galactic region. Turbulent velocity dispersions profiles, $\sigma_{\text{turb}}(r)$, built on this definition of $v_{\text{turb},i}$, are presented in
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of the specific magnetic (red), thermal (yellow), and turbulent (blue) energies within the galactic region for the B-Mech (top), B-NoFb (middle), and B-RdTh (bottom) runs. All energies have been normalized to a unique $\epsilon_{\text{ta}}$ value (see the text). To facilitate the comparison, $\epsilon_{\text{mag}}$ curves have been boosted by 15 orders of magnitude.</figcaption></figure>

Fig. 9. The formation of a stable gas disc occurs around $z \sim 4.5–4.0$ in all three simulations, and is confirmed by visual inspection. This transition correlates with the end of the accretion phase and occurs shortly after the gravitational potential of the galactic region becomes dominated by the stellar component. For simplicity of the analysis, the transition is taken to occur at the same $t_{\text{acc/fb}}$ ($z_{\text{acc/fb}} = 4$) in all runs.

Two features are identified during the accretion-dominated phase in Fig. 6: (i) sharp peaks are present throughout the entire phase, and (ii) the underlying continuum upon which these peaks rest grows steadily with time. We attribute the presence of peaks to individual mergers, whilst the growth of the smooth underlying component corresponds to *in situ* turbulent amplification. They are consequently addressed separately.

Monotonic turbulent amplification of the specific magnetic energy results in an exponential growth $\epsilon_{\text{mag}} \propto e^{\Gamma_{\epsilon}t}$, with characteristic time-scale $\Gamma_{\epsilon}$. Such an exponential growth typically indicates a turbulent dynamo process (Federrath 2016). Exact values of characteristic time-scales obtained by fitting exponential curves to $\epsilon_{\text{mag}}$ (see Fig. 6) are given in Table 2 for both accretion- and feedback-dominated epochs. During the accretion-dominated epoch, they are approximately of the order of $\Gamma_{\epsilon, \text{acc}} \sim 2 \text{ Gyr}^{-1}$, in fair agreement with the analytic estimate $\Gamma_{\epsilon} \sim \sigma_{\text{turb}}/l \sim 1.5–3 \text{ Gyr}^{-1}$, provided one picks a dissipation scale, $l \sim 10–20$ pc, typical of the maximal resolution achieved in our AMR simulations. From the numbers listed in Table 2, we see that during the accretion phase, $\Gamma_{\epsilon}$ is larger in the absence of feedback. The amplification is also reduced as the strength of the sub-grid feedback implementation increases (the B-RdTh run yields a lower value than the B-Mech run). Stellar feedback thus has a detrimental effect on magnetic energy amplification during the accretion phase. Potential explanations of this behaviour are an increase in turbulent diffusion caused by stellar feedback (Gressel et al. 2008), a suppression of gravitational instabilities, feedback diluting magnetic fields into the hot phase, or a reduction of small scale solenoidal turbulence relative importance (Grisdale et al. 2017). More specifically, stellar feedback could be hampering clump formation triggered by larger scale gravitational instabilities in the galaxy. As solenoidal modes are expected to produce higher growth rates (Federrath et al. 2011a) and stellar feedback can efficiently drive solenoidal turbulence (Korpi et al. 1999), one would expect that in the presence of stellar feedback, the fraction of specific turbulent energy in solenoidal modes within the galactic region should increase. However, we find the opposite during the accretion phase: stellar feedback reduces the fraction of solenoidal specific turbulent energy in our runs on scales $L < 1$ kpc by $\sim 10$ per cent. A detailed study of this complex interaction between feedback, turbulence modes, and magnetic fields amplification is beyond the scope of this paper. However, we note that, despite the anisotropic nature of the stellar feedback implementation in the B-Mech run, which has the potential to yield a higher amplification of the magnetic energy than the B-NoFb run because of an intrinsically more efficient capability to drive solenoidal turbulence than its isotropic B-RdTh counterpart, this does not happen in our simulations. This interpretation is underpinned by the measurement of the specific turbulent energies presented in Fig. 8. It is clear from Fig. 8 that, at least for $z > 4$, the level of turbulent energy (blue curves) measured in the B-NoFb run is at least comparable to, if not higher than that measured in both feedback runs. We point out that these findings corroborate both the observational standpoint that stellar feedback does not rank as the favoured mechanism to drive turbulence in high-redshift clumpy galaxies (e.g. Elmegreen et al. 2009b) and cosmological simulations of magnetized Milky Way-like galaxies performed with a different (moving mesh + divergence cleaning) technique (Pakmor et al. 2014).

As feedback heats the ISM gas, another significant difference between the B-NoFb run and its feedback counterparts, is that this former displays a larger difference between thermal and turbulent specific energies (yellow and blue curves in Fig. 8, respectively). This supports the claim by Evirgen et al. (2017) that the contribution from the hot phase of the ISM to the amplification of magnetic energy is negligible. A possible explanation may be that flow resides in this phase for a short time compared to the dynamo time-scale. By contrast, these two time-scales become comparable for the warm and cold phases, indicating that these latter are responsible for most of the amplification. For these reasons, we expect the cold phase to be the most efficient at amplifying the field in our simulations.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Galactic region radially averaged velocity profiles for B-Mech (blue), B-NoFb (red), and B-RdTh (green). Cylindrical coordinates are employed for $z = 2$, with the poloidal component pointing in the direction of the galactic region angular momentum vector. Solid lines display the tangential $v_{\text{tan}}$ (or toroidal $v_{\phi}$) component and dashed lines the turbulent velocity dispersions $\sigma_{\text{turb}}$. Vertical dashed lines represent $0.1 r_{\text{halo}}$. As redshift decreases, galaxies become less dominated by turbulence and a clear rotation curve is established.</figcaption></figure>

<table><caption>**Table 2.** $\Gamma_{\epsilon}$ time-scales derived from exponential fits to specific magnetic energies in the accretion phase, $\Gamma_{\epsilon, \text{acc}}$, and the feedback phase, $\Gamma_{\epsilon, \text{fb}}$. All $\Gamma_{\epsilon}$ time-scales are measured in Gyr$^{-1}$. See the text (Section 4) for details of the fitting procedure.</caption><thead><tr><th>Run</th><th>Γ<sub>ε, acc</sub></th><th>Γ<sub>ε, fb</sub></th></tr></thead><tbody><tr><td>B-NoFb</td><td>2.7 ± 0.2</td><td>-2.0 ± 0.5</td></tr><tr><td>B-Mech</td><td>2.3 ± 0.1</td><td>0.3 ± 0.2</td></tr><tr><td>B-RdTh</td><td>2.0 ± 0.1</td><td>1.7 ± 0.2</td></tr></tbody></table>

In an effort to quantify this, Fig. 10 displays density–temperature ($\rho_{\text{gas}}-T_{\text{gas}}$) phase diagrams, colour coded by specific magnetic energy, for the galactic region at $z = 6$. In the three runs, we have divided the ISM following Gent (2012) and Evirgen et al. (2017). We distinguish three different ISM phases: cold and dense, warm, and hot and diffuse. Hereafter, we will simply refer to them as cold (blue), warm (orange), and hot (purple). These phases are separated by constant specific entropy lines in Fig. 10. The fraction of gas mass contained in each phase is shown in Fig. 11. We find that the fraction of gas in the hot phase is negligible (less than 10 per cent) in all three runs and at all times (see the bottom panel of Fig. 11). During the very early stages of evolution ($z \ge 7$), the amount of hot gas in the B-NoFb run is significantly lower than that of both feedback runs (Agertz, Teyssier & Moore 2009b), but it progressively grows as the cooling time-scale overtakes the compression time-scale and the stable standing shock is pushed outside of the galactic region (Dekel & Birnboim 2006). Indeed, at later times ($z \le 4.5$), the B-NoFb run even has comparatively more hot gas than both feedback runs due to its less efficient, quasi-metal-free cooling.<fnref n="1" /> In any case, as is clear from Figs 10 and 11, most of the specific magnetic energy concentrates in the mass dominant cold phase as well as, to a lesser extent, the colder part of the warm phase. This trend is especially marked for the no B-NoFb run (the middle

panel of Fig. 10), where $\epsilon_{\text{mag}}$ is strongly concentrated in the density and temperature ranges $10^{-23}$ g cm$^{-3} < \rho_{\text{gas}} < 10^{-20}$ g cm$^{-3}$ and $10^3$ K $< T < 10^4$ K, respectively. In the feedback runs, the specific magnetic energy occupies a wider locus in the phase diagram, and is more homogeneously distributed across the cold and warm phases. As previously mentioned, this reinforces the idea that SN feedback is detrimental to magnetic energy amplification during the accretion phase: because feedback heats up the gas, during most of the accretion-dominated phase ($12 < z > 4.5$), the B-NoFb run contains on average 10–20 per cent more gas in the cold phase than both feedback runs (Fig. 11).

Another striking feature of the phase diagrams presented in Fig. 10 is the absence of correlation between specific magnetic energy and gas density in all three runs, which would have, had it been present, indicated an adiabatic-like amplification of the magnetic field. Instead, Fig. 9 shows that the B-NoFb run has the largest absolute velocity dispersion $\sigma_{\text{turb}}$ of the three simulations throughout a large fraction of the galactic region ($r \lesssim 0.1 r_{\text{halo}}$) at high redshift ($z \ge 4.5$). This is not a very large effect, as these values are generally comprised between 40 and 80 km s$^{-1}$ in the B-NoFb run, compared to 30 and 60 km s$^{-1}$ for the B-Mech and B-RdTh runs, but it is nonetheless systematic. In all three runs, the turbulent velocity profiles subtly increase with $r$, so that in the external galactic region ($r \gtrsim 0.1 r_{\text{halo}}$), they all reach typical values of the order of $\sigma_{\text{turb}} \sim 50–100$ km s$^{-1}$, mainly dominated by radial velocity dispersion. This is not surprising as turbulence is expected to be higher in the external region of galaxies when significant accretion flows are present (Klessen & Hennebelle 2010; Hopkins et al. 2013). As redshift decreases, the velocity dispersion in the B-RdTh and B-Mech runs increases monotonically, to reach a similar (or even slightly higher) level that in the B-NoFb run by $z = 4$. It is interesting to note it is at this redshift when the disc forms, even though there are hints of an earlier presence in the B-NoFb run.

Throughout the accretion period and particularly at the end ($z \sim 4$), all three runs display sharp peaks in their specific magnetic energies (see Fig. 6). These peaks are correlated with mergers, which stretch the magnetic field out of galaxies in tidal streams during close encounters, as reported in previous studies (see Roettiger,

<footnote n="1">¹Note however that, at least for $z > 2.5$, this hot halo gas is unable to sever the filaments of cold gas directly supplying the galactic region (Tillson et al. 2015).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Gas density ($\rho_{\text{gas}}$)–temperature ($T_{\text{gas}}$) phase diagrams colour coded by specific magnetic energy ($\epsilon_{\text{mag}}$). Diagrams correspond to the galactic regions at $z = 6$. Simulations from left to right are B-Mech, B-NoFb, and B-RdTh. Following Gent (2012), three ISM phases are identified by shaded regions: hot (purple, top left), warm (red, central), and cold phases (blue, bottom right). Divisions between ISM phases correspond to lines of constant specific entropy $s \propto [\log T - (\gamma - 1)\log \rho_{\text{gas}}]$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Gas mass fraction within the different ISM phases of the galactic region (hot – purple; warm – orange; cold – blue). Runs are from top to bottom B-Mech, B-NoFb, and B-RdTh. For clarity, the bottom panel presents a zoom version of the hot phase mass fractions for the three runs.</figcaption></figure>

Stone & Burns 1999; Kotarba et al. 2011), and are followed by a subsequent decrease in magnetic energy generally attributed to numerical dissipation (Kotarba et al. 2011; Geng et al. 2012; Pakmor et al. 2014). Note that the main difference between our work and these simulations (besides our order of magnitude better spatial and mass resolution that makes the merger peaks more prominent) is that they typically reach magnetic field values close to saturation before the merger happens, regardless of the initial value of the magnetic field they use. This does not happen in our simulations with the magnetic seed employed here. We do reach energy equipartition for higher initial values of the magnetic field (Martin-Alvarez et al., in preparation). Such a saturation of the magnetic field is generally interpreted in these papers as evidence of numerical convergence. However, given the high level of residual magnetic divergence inherent to some of the divergence cleaning methods used in these works, one is left to wonder about the role played by numerics in establishing the mean level of field amplification. This is especially worrisome given the results reported in the comparative study presented by Mocz et al. (2016), where the magnetic pressure in two idealized isolated disc simulations that differ only by their induction equation solver (CT versus divergence cleaning) ‘converges’ to values different by more than an order of magnitude (see fig. 6 of Mocz et al. 2016).

Fig. 12 displays projections of the gas density, temperature, and pressure, along with the specific magnetic energy in a volume centered on the galactic region and coincident with a merger occurring around $z \sim 6$ in all three runs. The peaks in the specific magnetic energy coincident with the merger event are also marked as black circles in Fig. 6. From the combination of Figs 6 and 12 is clear that mergers boost magnetic energy within the galactic region (solid circles on the gas density panels of Fig. 12) by stretching field lines. Moreover, $\epsilon_{\text{mag}}$ peaks in the B-NoFb run achieve considerably higher values than their feedback run counterparts, as field lines smoothly follow merger stream flows in the absence of feedback. In the presence of feedback these very same field lines are strongly disrupted, and even though the level of turbulence may rise within the galactic region and amplify the field, a fraction of the resulting magnetic energy is expelled by outflows into the halo
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Projections of a (25 kpc)³ cubic volume centred on each galactic region for the three runs during a merger even ($z \sim 6$). Density, temperature, magnetic energy density, and thermal pressure are displayed in columns from left to right, respectively (B-Mech run: top panels; B-NoFB run: centre panels; B-RdTh run: bottom panels). In both feedback runs, magnetic field lines around the central galaxy are disrupted by the strong outflows of hot gas. These outflows reduce the efficiency of magnetic stretching resulting from the merger, generating shallower peaks in the specific magnetic energy compared to the no feedback run. Galactic regions are indicated as circles on the gas density images (left column) and intensive variables are mass weighted. The exact redshifts for each image correspond to the points circled in black on Fig. 6.</figcaption></figure>

(see also Dubois & Teyssier 2010), generating the radial filamentary magnetic structures seen in the magnetic energy density maps of Fig. 12. More quantitatively, a higher specific magnetic energy is measured in the outer regions of haloes in the feedback runs (see Fig. 7). As a consequence, even though a significant amount of magnetic energy is generated by bulk stretching during mergers, a large fraction of this energy is lost through outflows in the presence of feedback. Once the merger concludes and the galaxy settles, significant numerical magnetic reconnection occurs as the gas (and field) falls onto the galaxy and further of the magnetic energy gained in the merger is lost (see Figs A1 and A2 where the amplitude of the merger peaks is greatly reduced in lower spatial resolution, but otherwise identical runs). Since these losses are strongly correlated with resolution, better resolving the circumgalactic medium seems required in order to properly capture the amplification of magnetic energy by mergers. At fixed resolution, the use of fixed nested grids as implemented in Vazza et al. (2014) for LSS simulations, which set a unique turbulent and magnetic energy dissipation scale, also alleviate the problem. Bearing in mind these resolution caveats, we find, in agreement with Rieder & Teyssier (2017b) for cosmological simulations of dwarf galaxies, that mergers have a negligible

contribution to the evolution of the magnetic energy within the galactic halo.

To further assess the presence of a turbulent dynamo, we measure turbulent kinetic and magnetic energy spectra for each simulation (Fig. 13). These are calculated by performing a Fast Fourier Transform (FFT) of cubic boxes with 1024³ cells, centred on the galactic region in each run, and with physical box length of $l = 0.4 r_{\text{halo}}$, so as to encompass the whole galactic region. For the sake of numerical performance, our FFT calculation assumes periodic boundary conditions, which, although obviously unrealistic at these scales, still provides a very good approximation because of the high gas density and magnetic field contrast between the galaxy and its surroundings. In other words, the tenuous circumgalactic medium acts as a zero-padding region that isolates the galaxy, ensuring that the power in Fourier modes at any scale but the largest is dominated by that of the galaxy ISM. The impact of periodicity and shape of the galaxy disc on the measured spectra are discussed in more detail in Appendix B, so we only mention at this stage that the disc shape results in a power-law decay of the density spectra as a function of the spatial angular wavenumber $\propto k^{-7/4}$ for $k_s \ge \frac{\pi}{h_s}$. As previously discussed, kinetic energies were calculated from velocities with
<figure><img src="image_14.png" /><figcaption>Figure 13. Comoving turbulent kinetic (left) and magnetic (right) energy power spectra for B-Mech (top row), B-NoFb (middle row), and B-RdTh (bottom row). Dashed lines represent Kolmogorov (black), and Burgers (grey) turbulent scalings for the kinetic energy; and Kazantsev scaling for the magnetic energy. Spikes in the turbulent kinetic energy spectra correspond to resonances with nested AMR grid levels.</figcaption></figure>

bulk motion and circular velocity $v_c(r)$ removed. All spectra are displayed in comoving units: dark matter halo physical radii grow $\propto (1+z)^{-1}$, so their comoving size remains roughly constant with time, which facilitates the analysis and comparison of the different spectra across the entire redshift range.

The turbulent kinetic energy spectra for B-NoFb from $z \sim 8$ to $z \sim 4$ exhibit a cascade with power-law scaling between Kolmogorov $\mathcal{E}_{\text{turb}} \propto k^{-5/3}$ (Kolmogorov 1941, indicated by black dashed lines in the left column panels of Fig. 13) and Burgers $\mathcal{E}_{\text{turb}} \propto k^{-2}$ (Burgers 1948, indicated by grey dashed lines on the same panels) at intermediate scales $0.3 \text{ kpc}^{-1} < k_{\text{co}} < 3 \text{ kpc}^{-1}$, i.e. in the inertial range, where $k_{\text{co}}$ represents the co-moving spatial angular wavenumber. These cascade spectra for the turbulent kinetic energy indicate a transfer of energy from large scales towards smaller eddies, consistent with fully developed turbulence resulting from continuous accretion. On the other hand, the feedback runs, and especially B-Mech, inject energy and/or momentum on small scales (see e.g. Kimm et al. 2015, for detail), as expected from SNe explosions. This generates extra power on scales ranging typically from $3 \text{ kpc}^{-1}$ <

$k_{\text{co}} < 30 \text{ kpc}^{-1}$, flattening the power-law index of the spectra as the scale height of the disc increases (see Appendix B). Note that the energy injection scale into the galaxy can be identified by the break present in all spectra around comoving wavenumber $k_{\text{co}} \sim 0.2 \text{ kpc}^{-1}$ (corresponding to physical wavelengths of $\sim 3$ and $\sim 10 \text{ kpc}$ at redshifts 8 and 2, respectively; see the left column panels of Fig. 13).

Because small-scale turbulence more rapidly amplifies the magnetic field due to its short e-folding times, the magnetic energy generated by the dynamo is fed back towards large scales, as evidenced by an inverse cascade in the magnetic energy spectra (the right column panels of Fig. 13). In agreement with simulations by Bhat & Subramanian (2013) and Federrath (2016), these inverse cascade spectra typically have a shallower slope than the Kazantsev spectrum ($\mathcal{E}_{\text{mag}} \propto k^{3/2}$ Kazantsev 1968) on large scales ($k_{\text{co}} \le 1 \text{ kpc}^{-1}$), characteristic of compressible flows with magnetic Prandtl numbers $Pm > 1$ in the kinematic stage (the magnetic energy is much smaller than the kinetic energy see Fig. 8). We estimate the numerical $2 < Pm < 10$ for our ideal MHD scheme, see Rembiasz et al. (2017) for detail. We argue that part of this shallowness arises from the
<figure><img src="image_15.png" /><figcaption>**Figure 14.** B-Mech comoving magnetic energy spectrum evolution throughout the merger occurring around $z = 4$. Redshift evolves from high (red) to low values (blue). As the merger proceeds, power is gained on small scales due to the window function set by the disc being blurred. Once the merger has finished, the field is confined in the galaxy again and the decay due to the shape of the disc progressively propagates to larger scales.</figcaption></figure>

disc-like shape of our galaxies and the complexity of the system (bursty injection of energy and momentum by SN in a multiphase ISM subject to accretion and mergers). Indeed, when going from the largest (box size) to the smallest scales (single grid cell), the power in magnetic energy spectra of isothermal turbulent boxes first increases as steeply as Kazantsev, then flattens as it reaches a maximum on scales corresponding to multiple elements of resolution, before decaying (e.g. fig. 5 of Haugen et al. 2004a). Our simulations display the same global qualitative behaviour, but the inverse cascade spans a range of wave numbers that comprises both the scale length and scale height of the galaxy, so the galaxy shape naturally influences the energy spectra slopes. A perhaps clearer imprint of this effect can be found on the smallest scales (large $k$s) where our spectra exhibit a steep 'disc-like' decay following $E_{mag} \propto k^{-7/4}$ (see Fig. B2 and discussion in the Appendix B) instead of the somewhat shallower decay observed in simulations of turbulent boxes (both for compressible and incompressible fluids; Haugen et al. 2004a; Schekochihin et al. 2004; Bhat & Subramanian 2013; Federrath 2016). Isolated galaxy simulations are expected to have less power on large scales (greater than the galaxy length scale) by construction. However, similar magnetic spectra to ours are also observed in such simulations (Rieder & Teyssier 2016; Butsky et al. 2017), and in a cosmological simulation of a dwarf (Rieder & Teyssier 2017b). A difference worth noting is that Butsky et al. (2017, see their fig. 5) obtain magnetic energy spectra in the kinematic regime steeper than Kazantsev on intermediate to large scales, which we speculate might be the result of their SN seeding of magnetic field.

For a more detailed look at the time evolution of the magnetic energy spectra when a merger occurs, Fig. 14 displays the evolution of $E_{mag}$ for B-Mech, around a merger event at $z = 4$. As redshift decreases, the stretching resulting from material and magnetic field being pulled out of the galaxy amplifies the field. Furthermore, as the magnetic field is stretched, the window defined by the disc of the galaxy is temporarily blurred, eliminating the $k^{-7/4}$ decay, and the spectrum displays an inverse cascade all the way to $k_{co} \sim 10 \text{ kpc}^{-1}$ ($k_{ph} \sim 50 \text{ kpc}^{-1}$). However, once the merger completes, tidal gas streams fall back onto the galaxy and magnetic energy is lost due to numerical reconnection. Present-day simulations seem to the lack resolution to properly capture the amplification of turbulence and magnetic energy associated with the merger process. As a result,

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Ratio of the inflowing (top) and outflowing (bottom) gas across the galactic region to the total baryonic mass contained within it during a dynamical time for all three studied runs. Dynamical time-scales of the galactic regions are computed for each run as $\tau_{dyn}(z) = 2\pi 0.2 r_{gal}(z) v_c^{-1} (r_{gal}(z))$.</figcaption></figure>

mergers only act as temporary enhancers of the field, both in the galaxy medium and in the circumgalactic medium (see also Rieder & Teyssier 2017b).

The presence of a turbulent cascade in the kinetic turbulent spectra and an inverse cascade in the magnetic energy spectra between $z = 8$ and $z = 4$, regardless of whether SN feedback is implemented or not, proves that during the accretion phase, a turbulent galactic dynamo is triggered, amplifying the magnetic energy on small scales and feeding it back to larger galactic scales.

## 3.3 Feedback-driven amplification

The final period of amplification is dominated by local galactic effects instead of accretion. Galaxies have grown to a size where gas accretion represents a small fraction of their baryonic mass content, as shown in Fig. 15. The progressive shutdown of cold accretion flows and the increased stability of the gas disc as the stellar component grows also reduce the impact of accretion (Dekel & Birnboim 2006; Klessen & Hennebelle 2010; Cacciato, Dekel & Genel 2012). Thus, inflow-induced turbulence declines, especially in the inner regions of the discs. Outer regions still seem to display accretion-driven turbulence (see Fig. 9 and Klessen & Hennebelle 2010; Hopkins et al. 2013) that could maintain or even amplify the field locally (Pakmor et al. 2017). Transfer of magnetic energy from the outer parts of galaxies towards the centre has been reported in other studies (Machida et al. 2013), significant enough to strengthen the inner field once a dynamo powerful enough is activated in the external region. However, this does not happen in the simulations presented in this work. Instead, as there are no sub-grid processes forcing turbulence nor injecting energy or momentum in the B-NoFb run, and accretion loses importance, both specific
and absolute turbulent energies, required to maintain the growth of the field, decay (Table 2). Numerical resistivity also contributes to the decay of the magnetic energy, but the decline of turbulence eventually leads to the onset of the final field topology as argued in Rieder & Teyssier (2017a). Note however that, contrary to what we measure for the B-NoFb run, even without an explicit feedback model (that of Springel & Hernquist 2003), Pakmor et al. (2014) report that after the formation of a Milky Way-like galaxy they still find linear growth of the magnetic field, which they attribute to differential rotation.

On the other hand, even after accretion has significantly diminished, runs with stellar feedback still exhibit considerable amounts of turbulent energy throughout the galaxy (at least an order of magnitude more that B-NoFb, see Fig. 13 and Federrath 2016 for a discussion of stellar feedback-induced turbulence). It is now well established that compressive forcing due to SNNe can drive solenoidal turbulence (Korpi et al. 1999; Mee & Brandenburg 2006; Federrath et al. 2010), and considerable amplification is expected in trailing SN shock waves (Ji et al. 2016). Having said that, because of considerable numerical turbulent and magnetic dissipation due to limited resolution, and the neglect of magnetic field injection by SN that might constitute a considerable source of magnetic energy (see e.g. Butsky et al. 2017), our results should be viewed as only providing a lower limit on the importance of SN during the feedback-dominated phase.

Bearing in mind these caveats, SN-driven turbulence still causes an exponential amplification of the magnetic energy in runs B-Mech and B-RdTh through a turbulent dynamo (as indicated by the turbulent kinetic and magnetic energy spectra corresponding to this phase in Fig. 13). Table 2 shows that B-RdTh displays a larger growth rate, $\Gamma_{\epsilon_{mag}}$, than B-Mech. However, both runs have lower $\Gamma_{\epsilon_{mag}}$ than during their accretion-driven phase. We suspect this is related to the fact that accretion and gravitational instabilities in a differentially rotating disc preferentially drive solenoidal turbulence that is more efficient at amplifying the magnetic field. Admittedly, this latter property will depend both on the exact sub-grid implementation of stellar feedback and on resolution, as vorticity generated behind shocks is difficult to capture, but compressible turbulence conversion into solenoidal modes within the small volume covered by a 2D-like galaxy is limited. As B-RdTh is more efficient at perturbing the gaseous disc of the galaxy and driving a fountain, it is therefore not surprising that this run displays a higher $\Gamma_{\epsilon_{mag}}$ than its B-Mech counterpart. B-Mech is, in turn, more efficient driving outflows (see Fig. 15), which are responsible for injecting magnetic energy in the rest of the halo, causing the higher specific magnetic energies seen in Fig. 7. In summary, whilst in the accretion driven phase, feedback, if anything, is detrimental to amplification ($\Gamma_{\epsilon_{mag}}$ values in Table 2 are slightly lower for B-Mech and B-RdTh than for B-NoFb), in the feedback-driven phase it is the main driver of the amplification the field, which results in the values of $\Gamma_{\epsilon_{mag}}$ being both very model dependent during this phase, and of opposite sign to B-NoFb (see Table 2).

Finally, in comparison with isolated simulations of MW galaxies (Rieder & Teyssier 2016), where there is no accretion phase, cosmological simulations seem to display a lower growth rate $\Gamma_{\epsilon_{mag}}$. This is especially interesting in the case of B-RdTh, as this run features a similar sub-grid feedback as one of the isolated galaxy runs by Rieder & Teyssier (2016). Small disparities could arise from a combination of a slightly diminished feedback specific energy injection and the sub-dominant cosmic accretion of poorly magnetized gas, which will act to somewhat reduce our $\Gamma_{\epsilon_{mag}}$. However, the most likely cause of significant discrepancy is that their fixed gravitational

potential, where no fresh gas is accreted after the initial collapse phase, facilitates the escape of outflows from the galaxy far into the halo, thus promoting the establishment of a galactic fountain able to stretch (and therefore amplify) the field more efficiently than in our cosmological set-up.

## 4 DISCUSSION

The main result of this paper is the identification of three distinct phases of magnetic energy amplification during a galaxy lifetime. In this section, we discuss the robustness of this result. Whilst there are physical arguments that support the idea that such an amplification behaviour in three periods might be universal, the employed numerical prescriptions affect the duration and amplification occurring during the phases in various manners.

First, a collapse-driven phase of amplification is identified separately from the later evolutionary phases of the galaxy. It is constrained to occur between the turnaround time of the collapsing perturbation and the moment when the forming galaxy has generated its first stars. We measure these two moments as $z_{ta} \sim 22$ and $z_{coll} \sim 13$. The latter redshift strongly correlates with a change of behaviour in the magnetic energy amplification, which during the collapse phase, is mostly the result of frozen magnetic field lines being dragged by the collapsing gas. Motions are dominated by gravity, and magnetization levels are comparable across all three simulations, with the stellar feedback sub-grid implementation having a noticeable, but sub-dominant effect. Whilst the existence of such a phase is difficult to fault, this does not mean that the amplification is unaffected by numerical considerations. In fact, two aspects have considerable impact on the amplification found during the collapse phase. One is numerical resolution, which is discussed in more detail in Appendix A. It mainly affects the amount of turbulent amplification measured above the expected level for isotropic collapse. The second one is the orientation and structure (i.e. spatial variation) of the initial seed. This aspect is not explored here, and will be further analysed in Martin-Alvarez et al. (in preparation), but we note that such variation of the initial seed field could alter the amplification factor of the specific magnetic energy by about a dex.

After this initial collapse, galaxies rapidly develop a multiphase ISM and are mostly fed by cold filamentary gas flows (see Tillson et al. 2015). The second phase is identified as the accretion phase of amplification. Turbulent motions in the multiple phase ISM are deemed responsible for the growth of the magnetic field. This implies the presence of turbulent dynamo amplification, for which $E_{mag}$ follows an exponential growth. During this epoch, we only measure a weak rotational support of the gas and any differential rotation amplification should therefore be negligible. Per contra, a turbulent dynamo fed by resolved vortical motions should exponentially amplify the magnetic energy. Filamentary accretion flows are expected to naturally drive solenoidal turbulence, but whether they can transfer a significant amount of turbulent kinetic energy to the entire ISM is still a matter of debate. However, we find that whenever implemented, stellar feedback, whilst not preventing amplification of magnetic energy to occur, appears to be detrimental to it. As a consequence, the turbulent amplification found during this phase is primarily attributed to accretion flow-driven gravitational instabilities. As amplification is linked to turbulence, a higher numerical resolution will yield considerably higher growth rates by pushing down the numerical dissipation scale and resolving a larger span of the turbulent cascade. The typical growth rates of specific magnetic energy that we measure during this phase are of the order
of $\Gamma_{\varepsilon\text{mag}} \sim 2 \text{ Gyr}^{-1}$ (see Table 2), and are relatively independent of the feedback prescription employed.

The beginning of the accretion phase coincides with the end of the collapse phase $z_{\text{coll}} \sim 13$, itself determined by cosmological (and star formation) arguments, but turbulent amplification inside the proto-galaxy is likely to have already started before this redshift. As time progresses, the mass growth of the galaxy and the decrease in gas accretion rate by the galaxy steadily reduce the relevance of cold accretion flows, eventually leading to the termination of accretion-driven magnetic energy amplification. In the absence of other turbulence driving mechanisms (B-NoFb run), any amplification is likely confined to the outer regions of the galaxy and an overall decay of the specific magnetic energy is observed. This decay is, by-and-large, the result of numerical resistivity. However, galaxies are expected to host energetic phenomena, such as stellar feedback, stirring up the ISM. The diminishing importance of accretion-driven turbulence thus leads to a smooth transition to a feedback-dominated regime. In our simulations, this latter is observed to take place after the galaxy mass becomes star dominated and shortly after the formation of a thin gas disc. In fact, one could even argue for the presence of an intermediate transition phase once the galaxies become star dominated (especially in the absence of feedback), when the stability of the gas disc is enhanced by the presence of a dominant stellar component, and the growth rate is slightly reduced. This should be interpreted with caution though, due to the complication of having to disentangle merger/interaction effects with external objects. The transition from accretion to feedback-driven amplification is estimated to occur at $t_{\text{acc/fb}} \sim 1.7 \text{ Gyr}$ (or equivalently, $z_{\text{acc/fb}} \sim 4$).

The final phase of amplification, which we dub the feedback phase of amplification, is characterized by lower amplification rates. Tables are now turned and stellar feedback is the main mechanism responsible for the amplification. As a result, the B-NoFb run that displayed the highest growth rate now exhibits a magnetic energy decay, whilst the B-RdTh run shifts from the lowest growth rate to the highest. As for the accretion-dominated phase, the turbulent nature of the amplification implies that higher numerical resolution will better capture the flow turbulence and lead to higher growth rates. A large spread of growth rates are found for this phase, which depend sensitively on feedback implementation $-2.5 \text{ Gyr}^{-1} < \Gamma_{\varepsilon\text{mag}} < 1.5 \text{ Gyr}^{-1}$.

The distinction between the accretion and feedback dominated phases is thus underpinned by the robustness of our $\Gamma_{\varepsilon\text{mag}}$ measurements. In Table 2, best fits of this parameter to magnetic growth time-scales are shown. These are calculated employing seven bins to fit data spanning either the entire duration of the accretion phase or that of the feedback phase, with an exponential function $A \exp[\Gamma_{\varepsilon\text{mag}}t]$. For each phase, data are smoothed using a 0.2 Gyr window function to minimize the contribution of sharp features and better recover the smooth background behaviour of the curve. As the end of the collapse phase is considered to be clearly defined, A is fixed to its value $\varepsilon_{\text{mag},C}$ for the accretion phase. It is left as a free parameter for the feedback phase. To probe the robustness of the transition between the two phases, in Fig. 16 we display the different values of $\Gamma_{\varepsilon\text{mag}}$ we obtain when varying $t_{\text{acc/fb}}$. We find that modifying the binning of the data or the size of the window function has a negligible impact on the fitted values of $\Gamma_{\varepsilon\text{mag}}$, but that they are more sensitive to the choice of the epoch at which the phase transition takes place, $t_{\text{acc/fb}}$. This is not very surprising, as varying $t_{\text{acc/fb}}$ results in the inclusion of a larger fraction of the 'wrong' phase

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Measured exponential growth parameter $\Gamma_{\varepsilon\text{mag}}$ versus time of phase transition $t_{\text{acc/fb}}$ for the three runs B-NoFb (red), B-Mech (blue), and B-RdTh (green). Solid coloured lines indicate accretion phase values $\Gamma_{\varepsilon\text{mag,acc}}$, whilst dashed lines represent feedback phase values $\Gamma_{\varepsilon\text{mag,fb}}$. Note that the distinction between phases is robust regardless of data binning or specific choice for the time at which the phase transition takes place.</figcaption></figure>

in the estimate of the value of $\Gamma_{\varepsilon\text{mag}}$. However, as shown by Fig. 16, the amount of uncertainty is quite limited: values of $\Gamma_{\varepsilon\text{mag}}$ remains fairly constant over a reasonable range of $t_{\text{acc/fb}}$ for all phases and all runs. Even for the phase where the estimate $\Gamma_{\varepsilon\text{mag}}$ varies the most (feedback phase of the B-NoFb run, dashed red line of Fig. 16), we measure a clear decay of the magnetic energy, in stark contrast with a clear amplification in the matching accretion phase (solid red line on the same plot).

The amplification of the specific magnetic energy reached in our galaxies by $z=2$ becomes all the more remarkable when considering that these latter mainly grow via low-magnetized gaseous inflows, and lose magnetic energy via outflows. Contrary to isolated simulations of galaxies (e.g. Dubois & Teyssier 2008; Rieder & Teyssier 2016), stellar feedback need not necessarily be the dominant component that drives magnetic amplification in cosmological simulations. However, once accretion becomes unimportant, these idealized simulations resemble more closely to the cosmic scenario, and stellar feedback probably plays an important role in preventing the decay of the magnetic energy or generating further amplification. In any case, whatever the driving process, numerical resolution is the decisive factor to obtain the correct rate of amplification (Rieder & Teyssier 2017b). Extrapolating our numerical results to physical ISM magnetic resistivities (and assuming non-ideal MHD effects on small scales such as ambipolar diffusion do not play a major role) would lead to saturation of the magnetic energy in time-scales shorter than a Gyr, i.e. within the accretion phase.

# 5 CONCLUSION

In this work, we performed high-resolution cosmological zoom-in simulations of a Milky Way-like galaxy (extending the NUT suite to MHD, see Powell et al. 2011 for detail about the original simulations) with the RAMSES-MHD code (Teyssier 2002; Fromang et al. 2006; Teyssier et al. 2006), which uses a CT scheme to solve the magnetic field induction equation and therefore guarantees as low a divergence of the magnetic field as numerically possible.

We studied the evolution of the magnetic energy within this galaxy starting from an extremely weak primordial seed field and us-
ing three simulations dubbed B-NoFb, B-Mech, and B-RdTh where the only difference between simulations is the stellar feedback sub-grid model implemented. We find that:

(i) the evolution of the magnetic energy can be decomposed as a smooth exponential amplification/decay component on top of which sharp spikes, corresponding to merger events, superimpose.

(ii) the amplification of the magnetic energy occurs in three main phases: an initial collapse phase, an accretion-driven phase, and a stellar feedback-driven phase. The two latter amplification phases occur because of the presence of a turbulent dynamo.

(iii) during the collapse phase, magnetic energy amplification closely follows an isotropic adiabatic approximation, but is already accompanied by turbulent motions that drive a larger growth of the magnetic field.

(iv) the accretion-driven phase accounts for a significant fraction of the amplification, and does not require the presence of stellar feedback to sustain the turbulent dynamo. As a matter of fact, during this phase, feedback is found to be *detrimental* to amplification, as the B-NoFb run reaches the highest amplification level.

(v) the contribution of mergers to the overall amplification is negligible. Mergers trigger temporary enhancement of the field, primarily by stretching magnetic field lines. In the absence of stellar feedback, these lines are not disrupted and generate sharp spikes of amplification, but once merging events conclude, the field lines fall back onto the galaxy and most of the corresponding magnetic energy is lost through numerical reconnection.

(vi) as the galaxies grow in size and cosmic accretion dwindles, their discs stabilize and, in the absence of stellar feedback, the level of turbulence decreases and the specific magnetic energy decays. This typically occurs around $z \sim 4$.

(vii) in the runs where stellar feedback is present, it is able to drive a significant turbulent dynamo *after* the accretion phase, which further amplifies the field. The most efficient feedback implementation (B-RdTh run), which was the most detrimental during the accretion phase, sustains the highest growth of specific magnetic energy.

Although we have performed quantitative measurements of the specific magnetic energy growth rates, $\Gamma_{\epsilon, \text{acc}}$, the numerical values reported in this work should only be considered as lower limits. This is because the numerical resolution, which determines the magnetic resistivity (as well as the shear and bulk viscosities) of the gas flow is necessarily limited. We therefore expect that, as the resolution is increased, the value of this latter will drop and the growth rate of magnetic energy will rise accordingly. We provided a first attempt at quantify this resolution effect in Appendix A, which we believe is somewhat encouraging. However, further work is needed to properly establish the level of convergence we have reached, given the complex interplay between the numerics, multiphase ISM and implemented sub-grid physics.

This difficulty notwithstanding, we can compare our results to other, similar numerical approaches. By the end of the collapse phase, the specific magnetic energy of our galaxies is found to lie a factor ~5–10 above the level predicted by isolated galaxy simulations or analytical estimates. Using a similar set-up to our B-RdTh run, albeit to simulate a smaller dwarf galaxy, Rieder & Teyssier (2017b) measure, in the entire zoomed region, a growth rate for the magnetic energy density $\Gamma_{\epsilon, \text{acc}}$ that is an order of magnitude larger than that in our galactic region during the accretion phase. On the other hand, the magnetic energy growth during the feedback phase in our simulation only is slightly weaker than in isolated spirals simulations by Rieder & Teyssier (2016).

Finally, given the level of resolution we can currently achieve, a much higher amplitude of the initial cosmic seed field, or an injection by stellar sources within the galaxy (as in e.g. Butsky et al. 2017) is needed to reach saturation of specific magnetic energy through turbulent galactic dynamo amplification. We plan to address the potential impact of the back-reaction induced by such a saturated field on the three-phase amplification of cosmic magnetic fields in future work.

## ACKNOWLEDGEMENTS

The authors kindly thank the referee for insightful comments and suggestions that contributed to improve the quality of this manuscript. This work was supported by the Oxford Hintze Centre for Astrophysical Surveys which is funded through generous support from the Hintze Family Charitable Foundation. SMA thanks M. Rieder for useful discussions and hospitality during his visit to the University of Zurich (UZH). This work is part of the Horizon-UK project, which used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure. The authors would like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work. http://dx.doi.org/10.5281/zenodo.22558. The authors also thank Taysun Kimm for sharing his mechanical feedback prescription and acknowledge the usage of the FFTW library: http://www.fttw.org/.

## REFERENCES

Agertz O., Lake G., Teyssier R., Moore B., Mayer L., Romeo A. B., 2009a, MNRAS, 392, 294

Agertz O., Teyssier R., Moore B., 2009b, MNRAS, 397, L64

Alves Batista R., Shin M.-S., Devriendt J., Semikoz D., Sigl G., 2017, Phys. Rev. D, 96, 023010

Anglés-Alcázar D., Davé R., Özel F., Oppenheimer B. D., 2014, ApJ, 782, 84

Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

Balsara D. S., Kim J., 2004, ApJ, 602, 1079

Bambic C. J., Morsony B. J., Reynolds C. S., 2018, ApJ, 857, 84

Beck R., 2015, A&A, 578, A93

Beck R., Wielebinski R., 2013, Magnetic Fields in Galaxies. Springer, Berlin, Germany, p. 641

Beck A. M., Dolag K., Lesch H., Kronberg P. P., 2013, MNRAS, 435, 3575

Bernet M. L., Miniati F., Lilly S. J., Kronberg P. P., Dessauges-Zavadsky M., 2008, Nature, 454, 302

Bhat P., Subramanian K., 2013, MNRAS, 429, 2469

Biermann L., 1950, Zeitschrift Naturforschung Teil A, 5, 65

Birnboim Y., Dekel A., 2003, MNRAS, 345, 349

Brandenburg A., Subramanian K., 2005, Phys. Rep., 417, 1

Broderick A. E., Chang P., Pfrommer C., 2012, ApJ, 752, 22

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240

Burgers J. M., 1948, Advances in Applied Mechanics. Academic Press, New York, p. 171

Butsky I., Zrake J., Kim J.-h., Yang H.-I., Abel T., 2017, ApJ, 843, 113

Cacciato M., Dekel A., Genel S., 2012, MNRAS, 421, 818

Chen Z., Jiang Z., Tamura M., Kwon J., Roman-Lopes A., 2017, ApJ, 838, 80

Dedner A., Kemm F., Kröner D., Munz C.-D., Schnitzer T., Wesenberg M., 2002, J. Comput. Phys., 175, 645
Dekel A., Birnboim Y., 2006, MNRAS, 368, 2

Devriendt J. E. G., Guiderdoni B., Sadat R., 1999, A&A, 350, 381

Dolag K., Grasso D., Springel V., Tkachev I., 2005, J. Cosmol. Astropart. Phys., 1, 009

Dubois Y., Commerçon B., 2016, A&A, 585, A138

Dubois Y., Teyssier R., 2008, A&A, 482, L13

Dubois Y., Teyssier R., 2010, A&A, 523, A72

Dunkley J. et al., 2009, ApJS, 180, 306

Elmegreen B. G., Burkert A., 2010, ApJ, 712, 294

Elmegreen B. G., Elmegreen D. M., Fernandez M. X., Lemonias J. J., 2009a, ApJ, 692, 12

Elmegreen D. M., Elmegreen B. G., Marcus M. T., Shahinyan K., Yau A., Petersen M., 2009b, ApJ, 701, 306

Evans C. R., Hawley J. F., 1988, ApJ, 332, 659

Evirgen C. C., Gent F. A., Shukurov A., Fletcher A., Bushby P., 2017, MNRAS, 464, L105

Federrath C., 2016, J. Plasma Phys., 82, 535820601

Federrath C., Klessen R. S., 2012, ApJ, 761, 156

Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M.-M., 2010, A&A, 512, A81

Federrath C., Chabrier G., Schober J., Banerjee R., Klessen R. S., Schleicher D. R. G., 2011a, Phys. Rev. Lett., 107, 114504

Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011b, ApJ, 731, 62

Federrath C., Schober J., Bovino S., Schleicher D. R. G., 2014, ApJ, 797, L19

Fromang S., Hennebelle P., Teyssier R., 2006, A&A, 457, 371

Geng A., Kotarba H., Bürzle F., Dolag K., Stasyszyn F., Beck A., Nielaba P., 2012, MNRAS, 419, 3571

Gent F., 2012, PhD thesis, Newcastle Univ.

Gheller C., Vazza F., Brüggen M., Alpaslan M., Holwerda B. W., Hopkins A. M., Liske J., 2016, MNRAS, 462, 448

Gressel O., Elstner D., Ziegler U., Rüdiger G., 2008, A&A, 486, L35

Gressel O., Elstner D., Ziegler U., 2013, A&A, 560, A93

Grisdale K., Agertz O., Romeo A. B., Renaud F., Read J. I., 2017, MNRAS, 466, 1093

Haardt F., Madau P., 1996, ApJ, 461, 20

Haugen N. E., Brandenburg A., Dobler W., 2004a, Phys. Rev. E, 70, 016308

Haugen N. E. L., Brandenburg A., Mee A. J., 2004b, MNRAS, 353, 947

Hennebelle P., Chabrier G., 2011, ApJ, 743, L29

Hennebelle P., Iffrig O., 2014, A&A, 570, A81

Hogan C. J., 1983, Phys. Rev. Lett., 51, 1488

Hopkins P. F., Raives M. J., 2016, MNRAS, 455, 51

Hopkins P. F., Kereš D., Murray N., 2013, MNRAS, 432, 2639

Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Hull C. L. H. et al., 2017, ApJ, 842, L9

Ji S., Oh S. P., Ruszkowski M., Markevitch M., 2016, MNRAS, 463, 3989

Kaviraj S. et al., 2017, MNRAS, 467, 4739

Kazantsev A. P., 1968, Sov. J. Exp. Theor. Phys., 26, 1031

Kimm T., Cen R., 2014, ApJ, 788, 121

Kimm T., Devriendt J., Slyz A., Pichon C., Kassin S. A., Dubois Y., 2011, preprint (arXiv:1106.0538)

Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900

Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS,

Kinney R. M., Chandran B., Cowley S., McWilliams J. C., 2000, ApJ, 545, 907

Kitchatinov L. L., Rüdiger G., 2004, A&A, 424, 565

Klessen R. S., Hennebelle P., 2010, A&A, 520, A17

Kolb E. W., Turner M. S., 1990, The Early Universe, Addison-Wesley Publishing Company, Redwood City, California

Kolmogorov A. N., 1941, Akademiia Nauk SSSR Doklady, 32, 16

Korpi M. J., Brandenburg A., Shukurov A., Tuominen I., 1999, A&A, 350, 230

Kotarba H., Lesch H., Dolag K., Naab T., Johansson P. H., Donnert J., Stasyszyn F. A., 2011, MNRAS, 415, 3189

Kroupa P., 2001, MNRAS, 322, 231

Krumholz M. R., Burkhart B., 2016, MNRAS, 458, 1671

Lesch H., Chiba M., 1995, A&A, 297, 305

Li H.-b., Dowell C. D., Goodman A., Hildebrand R., Novak G., 2009, ApJ, 704, 891

Machida M., Nakamura K. E., Kudoh T., Akahori T., Sofue Y., Matsumoto R., 2013, ApJ, 764, 81

Mandelker N., van Dokkum P. G., Brodie J. P., van den Bosch F. C., Ceverino D., 2017, preprint (arXiv:1711.09108)

Marinacci F., Vogelsberger M., 2016, MNRAS, 456, L69

Marinacci F., Vogelsberger M., Mocz P., Pakmor R., 2015, MNRAS, 453, 3999

Marri S., White S. D. M., 2003, MNRAS, 345, 561

Martel H., Shapiro P. R., 1998, MNRAS, 297, 467

Mee A. J., Brandenburg A., 2006, MNRAS, 370, 415

Mitchell P. D., Blaizot J., Devriendt J., Kimm T., Michel-Dansac L., Rosdahl J., Slyz A., 2018, MNRAS, 474, 4279

Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319

Mocz P., Pakmor R., Springel V., Vogelsberger M., Marinacci F., Hernquist L., 2016, MNRAS, 463, 477

Mulcahy D. D. et al., 2014, A&A, 568, A74

Mulcahy D. D., Beck R., Heald G. H., 2017, A&A, 600, A6

Murray N., Quataert E., Thompson T. A., 2010, ApJ, 709, 191

Neronov A., Vovk I., 2010, Science, 328, 73

Ntormousi E., Dawson J. R., Hennebelle P., Fierlinger K., 2017, A&A, 599, A94

Ocvirk P., Pichon C., Teyssier R., 2008, MNRAS, 390, 1326

Padoan P., Nordlund Å., 2011, ApJ, 730, 40

Pakmor R., Springel V., 2013, MNRAS, 432, 176

Pakmor R., Marinacci F., Springel V., 2014, ApJ, 783, L20

Pakmor R. et al., 2017, MNRAS, 469, 3185

Peebles P. J. E., 1994, ApJ, 429, 43

Planck Collaboration et al., 2016, A&A, 594, A19

Powell K. G., Roe P. L., Linde T. J., Gombosi T. I., De Zeeuw D. L., 1999, J. Comput. Phys., 154, 284

Powell L. C., Slyz A., Devriendt J., 2011, MNRAS, 414, 3671

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

Pudritz R. E., Silk J., 1989, ApJ, 342, 650

Rasera Y., Teyssier R., 2006, A&A, 445, 1

Ratra B., 1992, ApJ, 391, L1

Rees M. J., 1987, QJRAS, 28, 197

Rembiesz T., Obergaulinger M., Cerdá-Durán P., Aloy M.-Á., Müller E., 2017, ApJS, 230, 18

Rieder M., Teyssier R., 2016, MNRAS, 457, 1722

Rieder M., Teyssier R., 2017a, MNRAS, 471, 2674

Rieder M., Teyssier R., 2017b, MNRAS, 472, 4368

Roettiger K., Stone J. M., Burns J. O., 1999, ApJ, 518, 594

Rosdahl J., Schaye J., Dubois Y., Kimm T., Teyssier R., 2017, MNRAS, 466, 11

Rosen A., Bregman J. N., 1995, ApJ, 440, 634

Roškar R., Teyssier R., Agertz O., Wetzstein M., Moore B., 2014, MNRAS, 444, 2837

Santos-Santos I. M., Brook C. B., Stinson G., Di Cintio A., Wadsley J., Domínguez-Tenreiro R., Gottlöber S., Yepes G., 2016, MNRAS, 455, 476

Schekochihin A. A., Cowley S. C., Hammett G. W., Maron J. L., McWilliams J. C., 2002, New J. Phys., 4, 84

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

Schlickeiser R., 2012, Phys. Rev. Lett., 109, 261101

Schmidt M., 1959, ApJ, 129, 243

Schober J., Schleicher D. R. G., Federrath C., Bovino S., Klessen R. S., 2015, Phys. Rev. E, 92, 023010
Schober J., Rogachevskii I., Brandenburg A., Boyarsky A., Froehlich J., Ruchayskiy O., Kleeorin N., 2018, ApJ, 858, 124

Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A, 410, 611

Slyz A. D., Devriendt J. E. G., Bryan G., Silk J., 2005, MNRAS, 356, 737

Springel V., Hernquist L., 2003, MNRAS, 339, 289

Su K.-Y., Hopkins P. F., Hayward C. C., Faucher-Giguère C.-A., Kereš D., Ma X., Robles V. H., 2017, MNRAS, 471, 144

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Taylor A. M., Vovk I., Neronov A., 2011, A&A, 529, A144

Teyssier R., 2002, A&A, 385, 337

Teyssier R., Fromang S., Dormy E., 2006, J. Comput. Phys., 218, 44

Tillson H., Devriendt J., Slyz A., Miller L., Pichon C., 2015, MNRAS, 449, 4363

Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

Vazza F., Brüggen M., Gheller C., Wang P., 2014, MNRAS, 445, 3706

Vazza F., Ferrari C., Brüggen M., Bonafede A., Gheller C., Wang P., 2015, A&A, 580, A119

Vazza F., Brüggen M., Gheller C., Hackstein S., Wittor D., Hinz P. M., 2017, Class. Quantum Gravity, 34, 234001

Vazza F., Brunetti G., Brüggen M., Bonafede A., 2018, MNRAS, 474, 1672

Vogelsberger M. et al., 2014, MNRAS, 444, 1518

Wang P., Abel T., 2009, ApJ, 696, 96

White S. D. M., Frenk C. S., 1991, ApJ, 379, 52

Widrow L. M., 2002, Rev. Mod. Phys., 74, 775

Wittor D., Vazza F., Brüggen M., 2017, MNRAS, 464, 4448

Wolfe A. M., Lanzetta K. M., Oren A. L., 1992, ApJ, 388, 17

## APPENDIX A: RESOLUTION TESTS

As mentioned in the main body of the text, one of the factors that is expected to have a significant impact on the amplification of magnetic energy throughout the entire simulation is numerical resolution. The higher the resolution, the smaller the scale that can be captured by the simulation, and the larger the turbulence cascade that can be captured. Smaller scale turbulence has smaller e-folding times and therefore yields faster amplification. At the same time, the numerical viscosity, which also scales with resolution, is reduced and turbulence is thus dissipated at larger wavenumbers, increasing the inertial range. Finally, the numerical magnetic resistivity decreases as well, and the decay of the magnetic field is thus less pronounced.

As a simple way to quantify the importance of this resolution effect, five B-Mech like runs (Table 1) that differ exclusively in their maximum spatial resolution (B-Mech, B-MR20, B-MR40, B-BMR80, and BMR160 with $\Delta x_m = 10, 20, 40, 80$, and $160$ pc, respectively) are compared in Fig. A1. This figure is an analogue of Fig. 5, but for different spatial resolutions instead of feedback prescriptions. The adiabatic isotropic amplification estimate of the specific magnetic energy for the four runs with the highest resolution reaches a similar value (~$10\epsilon_{ta}$), but appears slightly inferior for the lowest resolution (~$7\epsilon_{ta}$). The collapse phase is completed around ~0.35 Gyrs for the B-Mech run and seems to be delayed no more than 0.05 Gyrs for the four highest resolutions, whilst it is delayed ~0.1 Gyrs for the lowest resolution (the gravitational force softening scales with spatial resolution). As expected, both the initial collapse amplification and the later turbulent growth are augmented as resolution increases, although there seems to be a reasonable level of convergence between B-Mech and B-MR20.

<figure><img src="image_18.png" /><figcaption>**Figure A1.** Specific magnetic energy growth during the collapse phase and the early stages of the accretion phase. The five runs plotted are identical B-Mech simulations with decreasing spatial resolution $\Delta x_m = 10, 20, 40, 80$, and $160$ pc from top to bottom, respectively. The approximate moment of collapse is delayed as AMR levels are removed. The amplification rate $\Gamma_{\epsilon_{mag}}$ is reduced as the turbulence is less and less resolved (Table A1).</figcaption></figure>
<table><caption>**Table A1.** Specific magnetic energy exponential fits for $\Gamma_{\epsilon_{mag}}$ in the five resolution runs. Values are for the period displayed in Fig. A1 (early), and the accretion (acc) and feedback (fb) phases. All $\Gamma_{\epsilon_{mag}}$ time-scales are measured in Gyrs$^{-1}$. Calculations of the growth rates are as for Table 2. Last column displays a bulk estimate of the magnetic Reynolds number $Rm$ in the galaxy (see the text). The estimated $Rm$ is above or similar to the critical magnetic Reynolds number required for turbulent amplification $Rm_{crit} \simeq 30$ in all runs but B-MR160.</caption><thead><tr><th>Run</th><th>Γ<sub>εmag,early</sub></th><th>Γ<sub>εmag,acc</sub></th><th>Γ<sub>εmag,fb</sub></th><th>Rm</th></tr></thead><tbody><tr><td>B-Mech</td><td>2.3 ± 0.2</td><td>2.3 ± 0.1</td><td>0.3 ± 0.1</td><td>160</td></tr><tr><td>B-MR20</td><td>2.4 ± 0.1</td><td>1.9 ± 0.1</td><td>0.6 ± 0.3</td><td>120</td></tr><tr><td>B-MR40</td><td>2.0 ± 0.1</td><td>1.9 ± 0.1</td><td>-0.2 ± 0.1</td><td>60</td></tr><tr><td>B-MR80</td><td>1.0 ± 0.2</td><td>1.2 ± 0.1</td><td>0.1 ± 0.2</td><td>30</td></tr><tr><td>B-MR160</td><td>0.7 ± 0.1</td><td>0.9 ± 0.1</td><td>-0.3 ± 0.1</td><td>15</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Figure A2.** Specific magnetic energy growth in the galactic region. The lines represent runs with decreasing resolution from dark blue (B-Mech) to crème (B-MR160). Regardless of the employed resolution, the three-phase evolution is recovered. B-MR160 appears to be a limiting case, displaying almost no growth. The displayed phases and divisions are the same as shown in Fig. 6.</figcaption></figure>

We expect this to be the result of only a fraction of the galaxy being resolved at 10 pc in B-Mech (see Fig. 2). The growth parameters measured for this early accretion phase, $\Gamma_{\epsilon_{mag,early}}$ are shown in Table A1. As expected, increasing spatial resolution leads to more rapid amplification rates (Rieder & Teyssier 2016; Pakmor et al. 2017).

To illustrate the robustness of the accretion and feedback phases with resolution, Fig. A2 presents the evolution of the specific magnetic energy at different spatial resolutions for the whole duration of the B-Mech simulation. The existence of the accretion and feedback phases appears to be quite insensitive to the numerical resolution. Once again, growth rates $\Gamma_{\epsilon_{mag}}$ measured for the different phases are shown in Table A1. These values are calculated in the same way as those reported in Table 2. Note that during the feedback phase only B-Mech and B-MR20 possess enough resolution to capture turbulent dynamo growth. In B-MR40 and B-MR160, the magnetic field decays. Similarly, the growth factor for B-MR80 is small and compatible with null.

Turbulent amplification also depends on the value of magnetic Reynolds numbers. Haugen et al. (2004a) and Brandenburg & Subramanian (2005) indicate that a magnetic Reynolds number $Rm$ above a critical value $Rm_{crit} \sim 30–35$ is required in order to trigger a turbulent dynamo, given $Pm = 1$. The critical Reynolds

number is found to be lower for the range of $Pm$ values expected in our simulations (see equation 10 and fig. 2 Haugen et al. 2004a, $Rm_{crit}(Pm) \sim 35Pm^{-1/2}$).

If we define a local magnetic Reynolds number as a function of the local diffusivity $\eta$, and typical length $L_{sys}$, and velocity $V_{sys}$ of the studied system

$$Rm = \frac{L_{sys} V_{sys}}{\eta}, \qquad (A1)$$

we can compute a bulk estimate of the magnetic Reynolds number in the ISM of the galaxy. In order to do that, we fix the scale length of the system to be equal to the entire thickness of the disc $L_{sys} = 2h_s \sim 1.2$ kpc. Note that $h_s$ is in the range of 0.6–0.4 kpc between redshifts $z=8$ and $z=2$, and that this scale is expected to somewhat increase with decreasing resolution.

Teyssier et al. (2006) argue that the amount of local numerical diffusivity introduced by the Godunov solver in RAMSES is of the order of $\eta \sim 0.5|u|\Delta x$, with $|u|$ the (1D) flow velocity, and $\Delta x$ the size of a resolution element. Assuming that $|u| \sim V_{sys}$, we can compute an approximate magnetic Reynolds number $Rm \sim Rm(\Delta x) = 2L_{sys}/\Delta x$. Under these assumptions, the calculated magnetic Reynolds numbers (presented in Table A1) fulfil $Rm \gtrsim Rm_{crit}$ for all runs but B-MR160. Accordingly, we expect some degree of turbulent amplification for all the different resolution runs but B-MR160, and possibly B-MR80. This seems in good agreement with the growth parameters presented in Table A1.

Growth rates likewise depend on turbulence properties such as Mach number $M$, or Helmholtz mode fraction. For significantly higher magnetic Reynolds numbers, Federrath et al. (2011a) find purely solenoidal turbulence forcing to produce higher growth rates than compressional forcing. Nonetheless, they find amplification for both types of forcing, specially in the case of supersonic turbulence. We thus expect a similar behaviour for our simulations, with $\langle M \rangle \sim 10^{1.5}$ and a typical fraction of solenoidal turbulence $\sim 50–70$ per cent in the galactic region. We also expect that the scaling of their supersonic growth rates $\propto M^{1/3}$ applies to our simulations.

This confirms the three-phase character of the turbulent dynamo amplification during its purely kinematic regime. Further increase of resolution is expected to lead to a higher level and more rapid amplification.

## APPENDIX B: FFTS OF A GALAXY

In this appendix, the characteristics of the spectrum of a disc galaxy embedded within a cubic box are addressed, as this is used extensively to interpret results in the main body of the text and disentangle physical properties from numerical artefacts in a given spectrum. We caution that our approach does not fully differentiate features caused in the spectral analysis by the gas density structure. A more localized analysis based on wavelets would be required to fully characterize which of the spectral features are exclusively generated by turbulence. The aim of this section is to provide a better understanding of the overall impact of the morphology of the galaxy. More specifically, one assumption commonly made in the calculation of an FFT is that of domain periodicity. Whilst it is not adequate for the case of an individual galaxy, an accurate spectrum can still be calculated if precautions are taken. Similarly, it is important to discriminate the galaxy from its environment in a cosmological simulation and to recognize features associated with its overall shape (disc scale height and length) instead of the dynamics of the gas flow within it.
<figure><img src="image_20.png" /><figcaption>**Figure B1.** Spectra for fiducial shapes of a galaxy. The impact of the shape parameters is assessed by plotting different thicknesses and radii for a disc (top) and a sphere (bottom). The three top lines for a disc have been shifted a factor × 1000 for the sake of clarity. Purple lines display the spectrum corresponding to an exponential density profile rather than a uniform one, as expected for a more realistic galaxy. The dashed black lines display the exponent of the power-law decay: $F_{\text{disc}}(k) \propto k^{-7/4}$ for a disc and $F_{\text{sphere}}(k) \propto k^{-2}$ for a sphere.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure B2.** B-Mech magnetic energy spectrum at $z = 2$ calculated as indicated in Section 3.2. The black dashed lines represent a power-law fit to the spectrum in the corresponding $k$ region, with the result of the fit shown on the figure (1 $\sigma$ errors shown in brackets). After $k \sim 10$ kpc$^{-1}$, the spectrum progressively converges to the $\mathcal{E}_{\text{mag}} \propto k^{-7/4}$ regime expected for an exponential disc (see the text for detail). This shape domination of the spectra is observed for studied spectra whenever magnetic fields are well confined to the galaxy.</figcaption></figure>

First, we have to discriminate the studied object from its environment. This occurs by definition in isolated simulations of galaxies, where all the effects only are the result of the presence of a galaxy. Fortunately, physical quantities (i.e. kinetic and magnetic energies) have values in the galaxy several orders of magnitude above the ones found in its environment. This naturally acts as an effective 0-padding: the Fourier modes outside the galaxy are negligible with respect to the amplitude obtained for the ones within the galaxy, as long as the wavelength remains smaller than the size of the region of interest. In physical terms, the energy $E(k)$ contained in a given scale $k$ smaller than the scale of the galaxy $k > k_{\text{gal}}$ is found to be negligible in the environment compared with the energy corresponding to the same mode within the galaxy. Accordingly, the studied scale $k$ has to be smaller than the largest scale corresponding to the galaxy. As the size of the studied galaxy roughly scales with the size of the box (with a side length $l = 0.2 r_{\text{halo}}$), and the radial scale of the galaxy is typically $\sim 0.1 r_{\text{halo}}$, this corresponds to $k_{\text{gal}} \sim 2 k_{\text{min}}$, where $k_{\text{min}}$ stands for the minimum spectrum scale represented. Smaller $k$s should therefore be discarded from the analysis even though we still display them on the figures.

In addition, wavenumber scales $k$ of the order of the overlapping of the galaxy with its periodic occurrences should be neglected. Correlations can be found between modes that bridge the shortest distance between the studied object and its periodic analogue. For the aforementioned size of the box and galaxy, the maximum scales $k$ that can be affected by periodicity are again $k_p \sim 2 k_{\text{min}}$.

Due to the small volume of the box covered by the galaxy, the next consideration is the actual shape of the galaxy. It can be simplified to either a sphere during early stages of formation and a thick disc otherwise. For the case of a disc, a uniform density disc of half-thickness $h = \alpha h_{\text{gal}}$ and polar radius $r_{\text{polar}} = \beta r_{\text{gal}}$ is studied, where $h_{\text{gal}} \sim 200$ pc and $r_{\text{gal}} \sim 3$ kpc are the half-thickness and radius resulting from fitting the gas density profile of the B-Mech galaxy with an exponential form. Fig. B1 displays how the variations of the $\alpha$ and $\beta$ parameters for a disc affect its spectrum. The strong oscillations seen in the spectrum are caused by the disc being implemented as a Heaviside-like window function. The thickness of the disc affects the scale of the first minimum, $k = \pi/h$ (i.e. the beginning of the decline). On the other hand, the radius of the disc affects how sharp the decay to the minimum is and the amount of power contained in the largest scales, as well as the mean value at $k = 0$. As the radius becomes larger, the spectrum gains power on scales $k < 3k_{\text{min}}$. Conversely, smaller radii lead to a lowering at $k < 3k_{\text{min}}$. As expected, the normalization of the curves is proportional to the size of the disc, affected by both parameters ($\alpha, \beta$). Finally, the decay followed by all disc spectra is fitted by $F(k) \propto k^{-7/4}$ (black dashed line, Fig. B1).

A spherical galaxy shape, on the other hand, leads to a more simple spectrum. It contains roughly constant power on the largest scales and starts decaying at the scale of the sphere radius. In Fig. B1, the spectra of a uniform density sphere of radii $r = \alpha r_{\text{gal}}$ for various $\alpha$ are displayed. The decay of the spectrum for the sphere follows $F(k) \propto k^{-1.92} \sim k^{-2}$. Changing the radius of the solid sphere exclusively affects the normalization and once the radius becomes small enough, causes a decay at scales larger than the radius scale due to periodicity effects disappearing.

Finally, to present a more realistic case, both the disc and the sphere are implemented decreasing exponential density profiles

$$v_{\text{cyl}}(r_{\text{polar}}, h) = \exp(-r_{\text{polar}}/r_{\text{gal}}) \cdot \exp(-h/h_{\text{gal}}), \qquad (\text{B1})$$

$$v_{\text{sp}}(r) = \exp(-r/r_{\text{gal}}). \qquad (\text{B2})$$
Their spectra transform in the same way as their uniform counterparts, but no longer present oscillations.

The measured disc decay $\mathcal{F}(k) \propto k^{-7/4}$ occurs on the smallest scales whenever a disc is present in the simulations. As an illustration, Fig. B2 displays the magnetic energy spectrum for B-Mech at $z = 2$. The decay on small scales starts to be dominated by the

shape of the galaxy around $k \sim 10\text{kpc}^{-1}$ and eventually follows the aforementioned $\mathcal{E}_{\text{mag}} \propto k^{-7/4}$, as indicated by the dashed black line.

This paper has been typeset from a TeX/Latex file prepared by the author.