Cosmological Magnetogenesis: The Biermann Battery during the Epoch of Reionization
Omar Attia, Romain Teyssier, Harley Katz, Taysun Kimm, Sergio Martin-Alvarez, Pierre Ocvirk, Joakim Rosdahl
MMNRAS , 1–14 (2021) Preprint 19 February 2021 Compiled using MNRAS L A TEX style file v3.0
Cosmological Magnetogenesis: The Biermann Battery during the Epoch ofReionization
Omar Attia, , ★ Romain Teyssier, Harley Katz, Taysun Kimm, Sergio Martin-Alvarez, Pierre Ocvirk, and Joakim Rosdahl Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 27, CH-8093 Zurich, Switzerland Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland Sub-department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK Department of Astronomy, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea Kavli Institute for Cosmology and Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg (ObAS), UMR 7550, F-67000 Strasbourg, France Université Lyon, Univ Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We investigate the effect of the Biermann battery during the Epoch of Reionization (EoR) using cosmological Adaptive MeshRefinement simulations within the framework of the Sphinx project. We develop a novel numerical technique to solve for theBiermann battery term in the Constrained Transport method, preserving both the zero divergence of the magnetic field and theabsence of Biermann battery for isothermal flows. The structure-preserving nature of our numerical method turns out to be veryimportant to minimise numerical errors during validation tests of the propagation of a Strömgren sphere and of a Sedov blastwave. We then use this new method to model the evolution of a 2.5 and 5 co-moving Mpc cosmological box with a state-of-the-artgalaxy formation model within the Ramses code. Contrary to previous findings, we show that three different Biermann batterychannels emerge: the first one is associated with linear perturbations before the EoR, the second one is the classical Biermannbattery associated with reionization fronts during the EoR, and the third one is associated with strong, supernova-driven outflows.While the two former channels generate spontaneously volume-filling magnetic fields with a strength on the order or below10 − G, the latter, owing to the higher plasma temperature and a marginally-resolved turbulent dynamo, reaches a field strengthas high as 10 − G in the intergalactic medium around massive haloes.
Key words: magnetic fields – galaxies: magnetic fields – galaxies: high-redshift – (cosmology:) dark ages, reionization, firststars – methods: numerical
Magnetic fields are an important component of the conductive plasmafound in many astrophysical objects, such as planets, stars, galaxies,but also in the intergalactic medium (IGM) between galaxies. Theirexact origin is however still debated: they could have formed in theearly Universe (Durrer & Neronov 2013) or at a later epoch, duringstructure formation through microscopic processes such as the Bier-mann battery (Biermann 1950) or the Weibel instability (Lazar et al.2009). Magnetic fields could also be generated inside stars throughstellar dynamos or within galaxies through galactic dynamos (Bran-denburg & Subramanian 2005) and then released into the interstellarmedium (ISM) by supernova (SN) explosions (Bisnovatyi-Koganet al. 1973) or by active galactic nuclei (AGN) jets (Rees 2005), andthen released into the IGM by galactic winds (Bertone et al. 2006).For the stellar or galactic dynamos to operate, a small seed mag-netic field is required. Among many alternative scenarios for suchseed fields, the Biermann battery, also known as the ‘thermo-electric ★ E-mail: [email protected] effect’, is probably the most robust and therefore the most populartheory. It has been observed in laboratory experiments (Fryxell et al.2010) and it has been simulated in the cosmological context withinaccretion shocks (Kulsrud et al. 1997) or ionization fronts (Gnedinet al. 2000) with fields generated by thermally-driven currents at thelevel of 10 − G.Galaxy formation models have evolved quite substantially over thelast two decades, with the increasing role of SN and AGN feedbackprocesses in regulating star formation and determining the morpho-logical type of galaxies (see Naab & Ostriker 2017, and referencestherein). Compared to the earlier galaxy formation models, strongergalactic outflows with correspondingly more intense ionizing radia-tion fields could quantitatively change the conclusions of these earlierworks. It has been also pointed out by Naoz & Narayan (2013) thatthe Biermann battery could be already effective during the linearevolution of the coupled dark matter (DM)/gas fluid during the darkages, leading to small seed fields as large as 10 − G, even beforeany star could have formed.Moreover, Graziani et al. (2015) recently showed that numericallymodelling the Biermann battery can be particularly challenging, andcan lead, in the case of sharp discontinuities in the flow, to a spurious © a r X i v : . [ a s t r o - ph . C O ] F e b O. Attia et al. numerical effect that they called ‘the Biermann catastrophe’. Thesevarious considerations have led us to revisit the early work of Kulsrudet al. (1997) and Gnedin et al. (2000), and study the Biermann batteryduring the Epoch of Reionization (EoR).The process of reionizing the Universe after recombination likelybegan after 𝑧 =
30 with the formation of the first stars (Wise et al.2012) and ended in the redshift range 𝑧 ∼ − We have implemented a novel numerical scheme for the Biermannbattery in the Adaptive Mesh Refinement (AMR) code Ramses(Teyssier 2002), which solves the ideal magneto-hydrodynamics(MHD) equations (Teyssier et al. 2006; Fromang et al. 2006). Thetemporal evolution of the magnetic field B is obtained by solving theinduction equation 𝜕 B 𝜕𝑡 = ∇ × E , (1)where the electromotive force (EMF) is given by E = u × B , with u thegas velocity. Equation (1) is at the core of many important physicalprocesses in cosmic magnetism such as large-scale and small-scalegalactic dynamos (Brandenburg & Subramanian 2005). In the idealMHD limit, however, we cannot create magnetic fields in a previ-ously unmagnetised fluid. Our interest precisely being to create seedmagnetic fields for subsequent dynamo processes, we need an addi-tional, non-ideal mechanism such as the Biermann battery (Biermann Figure 1.
Illustration of the naive implementation of the Biermann battery.The orange cells are used to estimate 𝜕 𝑥 𝑝 e at the centre of the purple cellthey frame. 𝐸 B ,𝑥 is computed at the golden edge by averaging the 4 centredvalues of the purple cells. 𝜕 B 𝜕𝑡 = ∇ × E B ≡ ∇ × (cid:18) 𝑐𝑒𝑛 e ∇ 𝑝 e (cid:19) , (2)where 𝑛 e and 𝑝 e are respectively the electron number density and theelectron pressure, 𝑒 is the elementary charge of the electron, and 𝑐 isthe speed of light. For the latter, note that we use the actual speed oflight, not the reduced speed of light in the radiative transfer module(see Sect. 4.1). The Biermann term stems from microscopic electriccurrents induced by the large difference in the electron and ion massand is able to generate magnetic fields ex nihilo as long as ∇ 𝑛 e and ∇ 𝑝 e are misaligned. In the Ramses code, we just add the non-idealBiermann EMF, E B , to the ideal EMF, E , naturally preserving thedivergence-free nature of our scheme.The Ramses code involves a tree-based data structure allowing re-cursive grid refinements on a cell-by-cell basis. The MHD equationsare solved using the second-order unsplit Godunov method basedon the monotonic upstream-centred scheme for conservation laws(MUSCL-Hancock method, Van Leer 1997; Toro 2013). The Con-strained Transport (CT) approach (Yee 1966; Evans & Hawley 1988)is used to advance Eq. (1) in time, which satisfies the solenoidalconstraint ∇ · B = 𝜕𝜕𝑡 ∬ B · d S = ∮ E · d l , (3)where d S and d l are respectively an elementary surface and length.As a consequence, while all the hydrodynamic variables (density,velocities, total energy) are located at cell centres, this approachrequires the magnetic field components to lie on the cell faces. In thisframework, the EMF, E , is defined at cell edges, so that its circulationdefines the B flux across face surfaces. Consequently, the Biermannbattery term E B must also be defined at cell edges.The most straightforward way to proceed is to first calculate theBiermann term in cell centres, where the quantities involved in itsdefinition are located. To this purpose, ∇ 𝑝 e is computed using the twoadjacent cells in the gradient direction. As an example, the expressionof the 𝑥 -directed Biermann battery term at the centre of the ( 𝑖, 𝑗, 𝑘 ) cell is 𝐸 𝑖, 𝑗,𝑘 B ,𝑥 = 𝑐𝑒𝑛 𝑖, 𝑗,𝑘 e 𝑝 𝑖 + , 𝑗,𝑘 e − 𝑝 𝑖 − , 𝑗,𝑘 e Δ 𝑥 , (4) MNRAS000
Illustration of the naive implementation of the Biermann battery.The orange cells are used to estimate 𝜕 𝑥 𝑝 e at the centre of the purple cellthey frame. 𝐸 B ,𝑥 is computed at the golden edge by averaging the 4 centredvalues of the purple cells. 𝜕 B 𝜕𝑡 = ∇ × E B ≡ ∇ × (cid:18) 𝑐𝑒𝑛 e ∇ 𝑝 e (cid:19) , (2)where 𝑛 e and 𝑝 e are respectively the electron number density and theelectron pressure, 𝑒 is the elementary charge of the electron, and 𝑐 isthe speed of light. For the latter, note that we use the actual speed oflight, not the reduced speed of light in the radiative transfer module(see Sect. 4.1). The Biermann term stems from microscopic electriccurrents induced by the large difference in the electron and ion massand is able to generate magnetic fields ex nihilo as long as ∇ 𝑛 e and ∇ 𝑝 e are misaligned. In the Ramses code, we just add the non-idealBiermann EMF, E B , to the ideal EMF, E , naturally preserving thedivergence-free nature of our scheme.The Ramses code involves a tree-based data structure allowing re-cursive grid refinements on a cell-by-cell basis. The MHD equationsare solved using the second-order unsplit Godunov method basedon the monotonic upstream-centred scheme for conservation laws(MUSCL-Hancock method, Van Leer 1997; Toro 2013). The Con-strained Transport (CT) approach (Yee 1966; Evans & Hawley 1988)is used to advance Eq. (1) in time, which satisfies the solenoidalconstraint ∇ · B = 𝜕𝜕𝑡 ∬ B · d S = ∮ E · d l , (3)where d S and d l are respectively an elementary surface and length.As a consequence, while all the hydrodynamic variables (density,velocities, total energy) are located at cell centres, this approachrequires the magnetic field components to lie on the cell faces. In thisframework, the EMF, E , is defined at cell edges, so that its circulationdefines the B flux across face surfaces. Consequently, the Biermannbattery term E B must also be defined at cell edges.The most straightforward way to proceed is to first calculate theBiermann term in cell centres, where the quantities involved in itsdefinition are located. To this purpose, ∇ 𝑝 e is computed using the twoadjacent cells in the gradient direction. As an example, the expressionof the 𝑥 -directed Biermann battery term at the centre of the ( 𝑖, 𝑗, 𝑘 ) cell is 𝐸 𝑖, 𝑗,𝑘 B ,𝑥 = 𝑐𝑒𝑛 𝑖, 𝑗,𝑘 e 𝑝 𝑖 + , 𝑗,𝑘 e − 𝑝 𝑖 − , 𝑗,𝑘 e Δ 𝑥 , (4) MNRAS000 , 1–14 (2021) he Biermann Battery during the EoR Figure 2.
Illustration of the correct implementation of the Biermann battery.The orange cells are used to calculate 𝑝 e and 𝑇 e at the purple vertex. 𝐸 B ,𝑥 iscomputed at the golden edge by using the two purple and golden end vertices. where Δ 𝑥 is the cell width in the 𝑥 -direction. E B is then derived oneach edge by averaging the values at the centre of the 4 adjacent cells(Fig. 1). While this naive implementation has the advantage of beingsimple, it is none the less vulnerable to spurious numerical effectsthat emerge at sharp discontinuities in the fluid variables. Indeed,the numerical Biermann magnetic field can diverge with increasingresolution at a discontinuity, resulting into what was called the Bier-mann catastrophe by Graziani et al. (2015). In addition, this naivediscretisation of the Biermann current can generate spurious fieldsin regions where ∇ 𝑛 e and ∇ 𝑝 e are aligned because of unavoidabletruncation errors in the gradient calculations.The alternative numerical implementation of the Biermann termwe propose is based on the following simple manipulation, using theelectron fluid equation of state 𝑝 e = 𝑛 e 𝑘 B 𝑇 e , E B = 𝑐𝑘 B 𝑒 ( 𝑇 e ∇ ln 𝑛 e + ∇ 𝑇 e ) . (5)As opposed to the first method, we compute the Biermann currenton the cell edges, using for the gradients the two end vertices, where 𝑛 e and 𝑇 e are calculated beforehand by averaging the cell-centredvalues of the 8 surrounding cells (Fig. 2). The corresponding EMFcomponent in the 𝑥 -direction is thus defined as 𝐸 𝑖, 𝑗 + / ,𝑘 + / B ,𝑥 = 𝑐𝑘 B 𝑒 (cid:18) 𝑇 e [ ln 𝑛 e ] Δ 𝑥 + [ 𝑇 e ] Δ 𝑥 (cid:19) , (6)where the edge temperature is defined as the average of the two endvertex temperatures as 𝑇 e = (cid:16) 𝑇 𝑖 + / , 𝑗 + / ,𝑘 + / e + 𝑇 𝑖 − / , 𝑗 + / ,𝑘 + / e (cid:17) , (7)and the two undivided differences are also computed using the twoend vertex log densities and temperatures as [ ln 𝑛 e ] = ln 𝑛 𝑖 + / , 𝑗 + / ,𝑘 + / e − ln 𝑛 𝑖 − / , 𝑗 + / ,𝑘 + / e , (8)and [ 𝑇 e ] = 𝑇 𝑖 + / , 𝑗 + / ,𝑘 + / e − 𝑇 𝑖 − / , 𝑗 + / ,𝑘 + / e . (9)This approach has the benefit of ensuring a zero field at machineprecision in the isothermal case. In this case, the Biermann electricfield is indeed proportional to the gradient of a scalar function, so thatits curl exactly vanishes. Thanks to the Stokes theorem, we can exactlypreserve this property at the discrete level. This results into whatis usually called a ‘structure-preserving scheme’. This alternative Figure 3.
Map of the temporal derivative of the magnetic field after onecomputational step in the smooth flow test. Both methods yield the sameresult. method, called the ‘correct’ method in this paper, is motivated bythe fact that the photo-ionized IGM tends to be globally isothermalin cosmological simulations. It should be noted that our methodslightly differs from the method proposed by Graziani et al. (2015),who decided to preserve the Rankine-Hugoniot relations at shocksinstead. Note also that the refinement of the edge-defined EMFs isautomatically and suitably carried out by the code (Teyssier et al.2006; Fromang et al. 2006).
In this section, we first present the numerical tests that we carry outto validate the newly-implemented Biermann battery module. Then,two complete suites of astrophysical tests are detailed, namely thepropagation of a Strömgren sphere and a Sedov blast wave. They actas a sandbox where the most important mechanisms in the ultimatecosmological simulations of the EoR are separately investigated, thatis to say propagating ionization fronts and shock waves.
We perform here simple tests of the Biermann battery solver alone,in a static configuration and for a single time step. This allows us tovalidate our new scheme with known analytical solutions. We focuson the differences between the naive and the correct method.
This first test problem proposed originally by Tóth et al. (2012)and Graziani et al. (2016) considers smooth electron density andpressure profiles with orthogonal gradients so that a strong Biermann-induced magnetic field is created with a known analytical solution.For simplicity, we use dimensionless units where 𝑐 = 𝑒 = 𝑘 B = B = , u = , 𝑛 e = 𝑛 + 𝑛 cos ( 𝑘 𝑥 𝑥 ) , and 𝑝 e = 𝑝 + 𝑝 cos (cid:0) 𝑘 𝑦 𝑦 (cid:1) , where 𝑘 𝑥 = 𝑘 𝑦 = 𝜋 / 𝑛 = 𝑝 =
1, and 𝑛 = 𝑝 = .
1. Consequently, 𝑛 e and 𝑝 e gradientsare indeed orthogonal, which maximises the effect of the Biermannbattery. The computational box spans a square domain from −
10 to
MNRAS , 1–14 (2021)
O. Attia et al.
Figure 4. L norm of the relative error to the analytical solution of the smoothflow test as a function of decreasing resolution ( 𝐿 represents the box size).Both methods (blue: naive, orange: correct) yield the same result and followthe dotted line of slope 2. Figure 5.
Map of the temporal derivative of the magnetic field after onecomputational step in the isothermal flow test for the naive (left) and correct(right) methods. The correct map only features machine-round-off errorsaround 0, which are not visible on this scale.
10 in both 𝑥 - and 𝑦 -directions. The resulting magnetic field in the 𝑧 -direction after one time step is given by the analytical form 𝐵 𝑧 Δ 𝑡 = − 𝑘 𝑥 𝑘 𝑦 𝑛 𝑝 sin ( 𝑘 𝑥 𝑥 ) sin (cid:0) 𝑘 𝑦 𝑦 (cid:1) ( 𝑛 + 𝑛 cos ( 𝑘 𝑥 𝑥 )) . (10)We use a resolution of 128 ×
128 points and a single time step torecover the analytic solution. Our two methods (naive and correct)yield the same result (Fig. 3), and both agree with Eq. (10) within arelative error of ∼ × − for this resolution.Exploiting the smooth character of this problem, we are able toperform a convergence study by varying the resolution of the grid andby computing the L norm of the error compared to the analyticalsolution. Figure 4 confirms that the results of the two methods arequite similar in this case, and also shows that the error decreases withthe mesh size as Δ 𝑥 , typical of standard second-order schemes. The second test is similar to the previous one with the difference thatthe initial density and pressure profiles are now parallel to each otherwith 𝑛 e = 𝑛 + 𝑛 (cid:16) sin ( 𝑘 𝑥 𝑥 ) + sin (cid:0) 𝑘 𝑦 𝑦 (cid:1) (cid:17) , (11) 𝑝 e = 𝑝 + 𝑝 (cid:16) sin ( 𝑘 𝑥 𝑥 ) + sin (cid:0) 𝑘 𝑦 𝑦 (cid:1) (cid:17) . (12)This corresponds to strict isothermal conditions. Because both den-sity and pressure gradients are now parallel, we do not expect anyBiermann-battery-generated magnetic field. The resulting maps for Figure 6.
Maps of the electron density (left panel), electron temperature(middle panel) and vertical (out of plane) magnetic field (right panel) at 𝑡 =
15 Myr for the str-20%-c simulation. 𝐵 𝑧 / Δ 𝑡 after one computational step are shown in Fig. 5. Note thatwe still use dimensionless units here. The naive method generatesa spurious magnetic field at the level of the truncation errors in thecalculation of the curl of the pressure gradient, around 10 − in theseunits. The correct method, on the other hand, generates vanishinglysmall values at the level of double precision numerical round-offerrors, around 10 − in these units. This proves the benefit of thecorrect method in preventing the formation of a spurious, residualmagnetic field in isothermal conditions where no battery is expected. We now investigate the effect of the Biermann battery in a situationthat is relevant for cosmological simulations, namely the propaga-tion of a Strömgren sphere. This test and the one presented in thefollowing section are particularly challenging for any Biermann bat-tery implementation, because they feature strong discontinuities thatcould result, as explained in Graziani et al. (2015), in the so-calledBiermann catastrophe. The propagation of a Strömgren sphere is alsothe basic process governing the Biermann battery within ionizationfronts during the EoR (e.g. Subramanian et al. 1994; Gnedin et al.2000). Note that from now on we switch to physical units.The Strömgren sphere simulations are performed in a two-dimensional (2D) square box of size 𝐿 × 𝐿 without AMR refinements,where 𝐿 = Δ 𝑥 = 𝐿 × − ℓ , where ℓ isthe refinement level. A source emitting ionizing photons of energy29.6 eV at a rate of (cid:164) 𝑁 = × s − cm − is placed at the origin,namely the lower left corner of the box. In order to create favourableconditions for the Biermann battery, the initial gas density profile isspatially modulated using 𝑛 = 𝑛 + 𝑛 cos ( 𝑘 𝑥 𝑥 ) , with 𝑘 𝑥 = 𝜋 / 𝐿 ,while the initial temperature profile is kept uniform with 𝑇 = 𝑇 . Thegas is composed of pure hydrogen of cross-section 1 . × − cm .Our initial conditions are similar to the physical conditions in theISM in a typical galaxy, namely 𝑛 = . − and 𝑇 = K.We perform a series of simulations in which we vary the amplitudeof the density modulation and the numerical method used (Table 1).The simulation time is fixed to 15 Myr.Results of this test for the simulation with 20% modulation andthe correct method (called here str-20%-c) are shown in Fig. 6. TheStrömgren sphere is clearly visible as the fully-ionized red bubblein the electron density map (left panel of Fig. 6). The initial densitymodulation can also be seen in the same panel as a blue-green sinewave in the background. This modulation causes a clear distortion ofthe sphere in the 𝑥 -direction, transforming it into an ellipsoid. Two MNRAS , 1–14 (2021) he Biermann Battery during the EoR Figure 7.
Vertical magnetic field profiles at the end of the Strömgren sphere test simulations. Each colour represents a different resolution, as labelled in thebottom left panel.Simulation name 𝑛 Methodstr-0%-n 0.0 × 𝑛 Naivestr-10%-n 0.1 × 𝑛 Naivestr-20%-n 0.2 × 𝑛 Naivestr-0%-c 0.0 × 𝑛 Correctstr-10%-c 0.1 × 𝑛 Correctstr-20%-c 0.2 × 𝑛 Correct
Table 1.
Table of the simulation parameters used in the Strömgren sphere testruns. discontinuities are present in this flow: the ionization front itself anda dynamical shock launched into the surrounding ISM. This rathercomplex structure is more visible in the temperature map (middlepanel of Fig. 6). Looking at the range of the colour map, though, wecan convince ourselves that the temperature evolution is quite closeto being isothermal.The propagation of the Strömgren sphere and its accompanyingshock generate a magnetic field of a few 10 − G in the vertical, 𝑧 -direction. We can also see a clear modulation of the magneticfield inside the Strömgren sphere, whose sign reflects the initialmodulation of the density field. The magnitude of this generatedfield can be estimated using a simple dimensional analysis of Eq. (2),with 𝐵 S = 𝑐𝑒 𝑘 B 𝑇 S 𝑟 𝑡 S , (13)where 𝑇 S , 𝑟 S and 𝑡 S are characteristic temperature, length and timescales. Using 𝑡 S =
15 Myr as the simulation time, 𝑟 S = . 𝑇 S = K, weget 𝐵 S (cid:39) − G, in very good agreement with our results.In Fig. 7, we show the radial profiles of the vertical magnetic fieldgenerated in our 6 different Strömgren sphere simulations. The up-per panels are for the naive method, while the lower panels show theresults of the correct method. The left column corresponds to thecase without any density modulation, for which density and pressuregradients are supposed to be both strictly radial with no Biermann battery. Unfortunately, both methods produce a spurious residualfield in this case because of numerical truncation errors, in particulardue to the mapping of a spherical profile onto a Cartesian grid. Notehowever that the correct method produces a spurious field that ismore than 50 × smaller than the naive method, which is a signifi-cant improvement. In both cases, the spurious field decreases withincreasing resolution, except close to the source at the origin and atthe ionization front. This is a problem similar to the one reported byGraziani et al. (2015), where such singularities and discontinuities inthe flow are particularly difficult to handle by the Biermann batteryterm. Note that the correct method performs significantly better thanthe naive method close to these singularities.In the middle and right columns of Fig. 7, we show the resultsobtained for the 10% and 20% modulation. We can clearly see thatthe correct method is now converging towards a unique profile, withmore difficulties at the ionization front. The naive method, on theother hand, is far from converging. It is clear that, in this case,the spurious field strength is always comparable to the correct fieldstrength, making the naive method totally unfit for this test. Only atthe highest resolution ℓ =
10 and for the largest modulation does thenaive method show signs of convergence towards what seems to bethe correct solution.
In this section, we perform a second series of tests relevant for theBiermann battery in shock waves (e.g. Kulsrud et al. 1997; Davies& Widrow 2000). We consider the well-known case of a 2D Sedovblast wave. While the Strömgren sphere maintains a quasi-isothermalregime, the Sedov blast wave exhibits large temperature variationswith sharp discontinuities, which cause numerical difficulties, es-pecially across shock fronts (Fatenejad et al. 2013; Graziani et al.2015).We proceed in a very similar manner to the previous section. Ourcomputational domain is a 2D square box with a width of 1 kpc.Here again we do not use any AMR refinements. A modulated initialdensity profile is used in order to generate a non-spherical shockwith a proper Biermann battery: 𝑛 = 𝑛 + 𝑛 cos ( 𝑘 𝑥 𝑥 ) and 𝑇 = 𝑇 MNRAS , 1–14 (2021)
O. Attia et al.
Figure 8.
Electron density (blue) and temperature (orange) profiles at 𝑡 = ℓ =
10. The exact analytical Sedov solutionsare over-plotted as dotted curves.
Figure 9.
Maps of the electron density (left panel), electron temperature(middle panel) and vertical magnetic field (right panel) at 𝑡 = ℓ = ( 𝑛 and 𝑇 have the same values as in the previous section). Aninitial energy pulse of 3 . × erg/pc is placed at the origin soas to have a blast wave similar to a cluster of 1000 combined SNexplosions. The adiabatic index of the gas is set to 𝛾 = .
4. We usethe same terminology to identify our different simulations as for theStrömgren sphere tests (i.e. as shown in Table 1). Contrary to theprevious section, the gas is assumed to be always fully ionized. Thesimulation time is set to 5 Myr.We first check that our numerical solution without modulationand with resolution ℓ =
10 reproduces the analytical Sedov solution(Fig. 8). The density jump at the shock front is close to 6, as predictedby the Rankine-Hugoniot relations, and the density and temperatureprofiles nicely follow the analytical solution. Note however that asignificant deviation is observed in the centre. This is due to the initialsingularity only captured by one cell at 𝑡 =
0. After the blast wave hasexpanded to its final radius at 𝑡 = 𝑟 = .
35 kpc.When we apply the initial 10% or 20% density modulation, thespherical symmetry of the propagating shock wave is broken into anelliptical shock (Fig. 9). A vertical magnetic field of strength up to10 − G is generated in the post-shock region. We can motivate sucha value using again a dimensional analysis with Eq. (13), yielding apredicted value of 𝐵 S (cid:39) × − G, using 𝑡 S = 𝑟 S = . 𝑇 S = K.Figure 10 compares our 6 different simulations. The left columnshows the case without modulation, for which no magnetic fieldshould be generated as the profiles are purely radial. Here again, the correct method performs better than the naive one, although not asspectacularly as for the Strömgren sphere test. A noticeable differenceis however visible right after the shock front. The naive methodclearly suffers from the phenomenon described in Graziani et al.(2015) as the Biermann catastrophe, for which the spurious magneticfield generated within the shock discontinuity increases when oneincreases the resolution. The correct method, on the contrary, showsa much better behaviour, with the spurious field decreasing withincreasing resolution, similar to the other correct method proposedby Graziani et al. (2015).The middle and right panels of Fig. 10 show the profile of thevertical magnetic field for the modulated cases. The two methodsseem to converge towards what we believe is the correct solution,with the correct method doing a slightly better job than the naiveone. The central region within 0 .
35 kpc is however not convergingat all in both cases. We believe this is due to the initial singularitygenerating stronger and stronger spurious fields in the centre, as thetemperature becomes higher, ultimately diverging in the very centre.As the resolution is increased, the problematic central region shrinks,but this occurs at a relatively slow rate. This also leads to a spuriousmagnetic field morphology in the centre of the explosion (Fig. 9).To conclude this section, we would like to stress that these tests areparticularly difficult as they feature many singularities and disconti-nuities, as well as a very modest deviation from spherical symmetry.In a more realistic three-dimensional environment, we expect shockwaves and ionization fronts to exhibit stronger non-parallel densityand pressure gradients and a less ambiguous Biermann battery pro-cess. Note however that our correct method produces significantlybetter results that the naive method, especially across shock fronts orin nearly isothermal conditions. It allows us to see converged resultsfor the two test suites we have proposed, although this convergenceappears to be surprisingly difficult to achieve.
In this section, we present the cosmological simulations we run tostudy the Biermann battery at play during the EoR. First, we de-scribe the global framework in which these simulations are carriedout, namely the Sphinx project, including the employed initial con-ditions. Then, we focus on the results of the multiple cosmologicalsimulations.
For our cosmological simulations we use the same setup and sub-gridmodels as in the Sphinx suite of simulations presented in Rosdahlet al. (2018). We refer to that paper for details, but briefly recap themain ingredients here.We run cosmological volumes of size 2.5 and 5.0 co-moving Mpc(cMpc) starting at 𝑧 =
150 and ending at 𝑧 =
6. The Λ CDM initialconditions are generated with Music (Hahn & Abel 2011) usingcosmological parameters consistent with the Planck 2013 results(Planck Collaboration 2014). We use hydrogen and helium massfractions of 𝑋 = .
76 and 𝑌 = .
24, respectively, and the gas has aninitial metal mass fraction of 𝑍 init = . × − in order to have highenough cooling rates to form the first stars at 𝑧 (cid:39)
15 in the absenceof molecular hydrogen formation, which is not modelled.The number of DM particles in the 2.5 and 5 cMpc volumesare 128 and 256 respectively, giving a mass resolution of 𝑚 DM = . × 𝑀 (cid:12) in both volumes. We consider a DM halo to be resolvedif it contains 300 particles or more, meaning we resolve halo masses MNRAS000
15 in the absenceof molecular hydrogen formation, which is not modelled.The number of DM particles in the 2.5 and 5 cMpc volumesare 128 and 256 respectively, giving a mass resolution of 𝑚 DM = . × 𝑀 (cid:12) in both volumes. We consider a DM halo to be resolvedif it contains 300 particles or more, meaning we resolve halo masses MNRAS000 , 1–14 (2021) he Biermann Battery during the EoR Figure 10.
Vertical magnetic field profiles at 𝑡 = (cid:38) × 𝑀 (cid:12) . The largest halo masses at 𝑧 = ∼ × 𝑀 (cid:12) and ∼ 𝑀 (cid:12) respectively .The gas cells (which are also used for the ionizing radiation fieldand gravitational potential from all matter) are adaptively refined,splitting a cell in 8 if it contains 8 DM particles or an equivalentmass in gas and stars, or if the local Jeans length is smaller than 4cell widths. The maximum number of such hierarchical splits fromthe coarse resolution is 8, so the coarse cell width is 19.6 ckpc (2.8kpc at 𝑧 =
6) and the minimum cell width reached in the ISM ofgalaxies is 76.3 cpc (10.9 pc at 𝑧 = 𝑀 (cid:12) , are formed out of the gas if 1) the hydrogen density 𝑛 H >
10 cm − , 2) the turbulent Jeans length is smaller than one cellwidth, and 3) the gas is locally convergent. The gas is converted intostellar particles with an efficiency that varies with the local virialparameter and turbulence (Kimm et al. 2017; Trebitsch et al. 2017;Martin-Alvarez et al. 2020b). The star formation is typically verybursty, as the local star formation efficiency can reach values evenhigher than unity under the right conditions.Each stellar particle produces a series of individual 10 erg typeII SN explosions that are stochastically sampled between 3 and 50Myr of its lifetime. An average mass fraction 𝜂 SN = . 𝑦 = . 𝑚 SN = 𝑀 (cid:12) , meaning each particle produces an averageof 40 SN explosions, which is about 4 times more than for a Kroupa(2001) IMF. This boost in the number of SN explosions is necessaryto suppress galaxy growth and reproduce the UV luminosity func-tion of galaxies in larger Sphinx volumes (see detailed discussionin Rosdahl et al. 2018). SN explosions are implemented with a me-chanical feedback model that injects the predicted final momentumof the snowplow phase if the Sedov–Taylor phase is not resolved, buta pre-snowplow kick otherwise (Kimm & Cen 2014). The halo mass, 𝑀 vir , is defined as the total DM mass within an ellipsoidalvirialised region and identified with the AdaptaHop halo finder (Tweed et al.2009) Ionizing radiation is injected into the grid cells from stellar parti-cles, using age- and metallicity-dependent luminosities derived fromthe Binary Population and Spectral Synthesis code (Bpass) version2.2.1 (Eldridge et al. 2007; Stanway & Eldridge 2018), assuming aninitial mass function with slopes of − . . . 𝑀 (cid:12) and − .
35 from 0 . 𝑀 (cid:12) . We use the M1 moment method to advectthe radiation, using the variable speed of light approximation (Katzet al. 2017) to reduce the computational cost. The speed of light atthe highest AMR resolution, corresponding to the ISM of galaxies, is0 . 𝑐 , but the radiation speeds up to 20% of the real speed of lightat the lowest AMR resolution, corresponding to cosmological voids.The Biermann battery being sensitive to the electron number den-sity, we are not affected by post-reionization artefacts in the residualneutral fraction due to the reduced speed of light framework (Ocvirket al. 2019). Moreover, Deparis et al. (2019) showed that adopting avalue of 10% of the speed of light delayed complete reionization byno more than Δ 𝑧 = . 𝑇 > K, metal-line cooling is computed withCloudy-generated tables (Ferland et al. 1998, version 6.02). Finestructure cooling rates from Rosen & Bregman (1995) are used inthe 𝑇 ≤ K regime, allowing the gas to radiatively cool to adensity-independent temperature floor of 15 K.The simulations also solve in parallel the ideal MHD equationsusing the Godunov scheme and a CT method (Teyssier et al. 2006;Fromang et al. 2006). Only the induction equation (Eq. 1) is modifiedas we add the Biermann source term to the EMF. The system ofequations is formulated using ‘supercomoving variables’ to account
MNRAS , 1–14 (2021)
O. Attia et al.
Simulation name
𝐿 ℓ min ℓ max Method Comment(cMpc)spx2.5-n 2.5 7 14 Naivespx2.5-c 2.5 7 14 Correctspx2.5-c-nosn 2.5 7 14 Correct SNe offspx5.0 5.0 8 15 Correct
Table 2.
Parameters of our various cosmological simulations. 𝐿 is the boxsize. ℓ min and ℓ max respectively refer to the minimum and to the maximumrefinement level. for the expansion of the Universe in the FLRW metric (Martel &Shapiro 1998). In doing so, only one additional source term, involvingthe scale factor 𝑎 , is introduced in the induction equation, whichconveniently leaves all the other MHD equations unchanged (moredetails are provided in Rieder & Teyssier 2017). The initial magneticfield is identically zero everywhere. Note that these computations aredone in double precision. In this section, we present the results of our cosmological simulationsusing radiation-magneto-hydrodynamics and our new implementa-tion of the Biermann battery. Table 2 summarises the parametersused for the various runs. We explore the effect of both the naive,less accurate Biermann battery scheme and the correct one. We alsocheck the influence of the box size on our results. We also performone simulation with SN feedback turned off.
We report here on the results of our fiducial spx2.5-c run, featuringthe correct Biermann battery scheme, SN feedback, and a volume of2.5 cMpc on the side. Figure 11 presents the final snapshot of thesimulation at 𝑧 =
6. The gas distribution shows a typical filamentarystructure surrounding large voids. The reionization of the IGM isalmost complete by that time. The whole box is nearly fully ionizedexcept for a small lower-density patch in which ionization fronts havenot had time yet to percolate. The temperature of the gas in ionizedregions is about 10 K. Note however the high temperature bubbles(10 − K) in the vicinity of the largest haloes due to feedbackfrom SN-driven explosions.The generated magnetic field is presented in Fig. 12. The upperpanel shows the 𝑧 -component of the magnetic field. All three com-ponents follow a similar pattern, namely a complicated morphologyalternating between positive and negative values imprinted by thecomplex shape of the density and temperature gradients. The mag-netic field is generally oriented towards the cross-product of thesegradients, broadly encircling H ii regions (Gnedin et al. 2000).In the lower panel, the magnetic field strength map reveals theco-existence of three different regimes: 1) a very low magnetic fieldstrength with 𝐵 (cid:39) − G in voids filled with neutral hydrogen, 2) alow, intermediate regime with 𝐵 (cid:39) − G in ionized regions, andfinally 3) higher values 𝐵 (cid:39) − G in a regime associated withSN-driven hot bubbles.The first mode is associated with the Biermann battery processduring the linear evolution of gas density perturbations, subject tothe complex thermal history after recombination (Naoz & Narayan2013).The second regime is associated with ionized regions and the clas-sical Biermann battery channel, in which the IGM is progressively magnetised as ionization fronts spread. The process can be depictedas propagating Strömgren spheres gradually filling the whole volumewith 10 − G fields, as illustrated in Sect. 3.2 with the Strömgrensphere test.The third mode can be seen as large bubbles hosting magneticfields two orders of magnitude higher than the previous case. Thesebubbles are localised in the vicinity of star-forming regions. In fact,they precisely lie within the above-mentioned high-temperature bub-bles, suggesting that this third regime is concentrated around massivehaloes where SN explosions occur. At this redshift, these SN-drivenbubbles are far from having percolated. This third magnetic fieldgeneration channel is not associated with ionization fronts but withthe propagation of galactic winds, in which magnetic fields are gen-erated within SN-driven shocks, in a similar way to our Sedov blastwave test in Sect. 3.3. It owes its higher amplitude (10 − G) to thehigher temperature of the gas. We will discuss this point later below.We have also compared the results of our two Biermann battery im-plementations, the naive and the correct one, in the density-magneticfield phase-space diagram shown in Fig. 13. We can clearly see hereagain the three modes we have just introduced. The lower branch,where the magnetic field strength is lower than 10 − G, correspondsto the linear regime prior to the EoR (Naoz & Narayan 2013). Wewill discuss more the origin of these very weak fields in Sect. 4.2.2.The upper branch shows the classical behaviour of ideal MHDwith the field strength scaling as 𝑛 / (shown as the dashed lines).Noticeably, we observe the presence of two distinct and parallel sub-branches that both scale as 𝑛 / . This bi-modal distribution is dueto the co-existence of the ionization generation channel and of theshock generation channel.When comparing the naive method (left panel of Fig. 13) to thecorrect method (right panel of Fig. 13), we see a very good qualitativeagreement, supporting the robustness of our conclusions. Quantita-tively, however, the field strength in the lower branch (associated withthe linear regime of the Biermann battery) is overestimated by thenaive method, owing to its largest truncation errors, as revealed inthe tests section.To support this interpretation, we would like to analytically derivethe amplitude of the Biermann battery generated fields using hereagain Eq. (13). For the ionization front branch, the length scale 𝑟 S can be estimated as the radius of the Strömgren sphere emanatingfrom a typical simulated galaxy with 𝑟 S = (cid:18) (cid:164) 𝑁 𝜋𝛽 rec ¯ 𝑛 (cid:19) / , (14)where 𝛽 rec (cid:39) . × − cm s − is the recombination rate at 𝑇 = K and ¯ 𝑛 is the mean baryon density at 𝑧 =
6, namely¯ 𝑛 (cid:39) − cm − . We consider for the ionizing photon rate (cid:164) 𝑁 = s − , corresponding to the luminosity of our largest simulated halo,after accounting for the proper escape fraction (Rosdahl et al. 2018).The resulting Strömgren sphere size is 𝑟 S (cid:39)
70 kpc, consistentlywith Fig. 12. For 𝑡 S , we use the age of the Universe at the end ofthe simulation with 𝑡 S = 𝑇 S , we use the ionized IGMtemperature with 𝑇 S = K. Eq. (13) hence yields an amplitude of 𝐵 low (cid:39) − G for the mean density of the Universe ¯ 𝑛 .As for the high- 𝐵 branch, the spatial extent of the hot magnetisedbubbles is dictated by galactic winds that spread the fields generatedin SN-driven shocks. To fully describe analytically such a processwould be beyond the scope of the current paper (see Bertone et al.2006, for example). We consider that a good approximation is to usethe observed characteristic radius of the hot bubbles 𝑟 S =
20 kpc (seeFig. 12). The most substantial change between the two regimes is theelectron temperature. An examination of the temperature distribution
MNRAS000
MNRAS000 , 1–14 (2021) he Biermann Battery during the EoR Figure 11.
Mass-weighted maps projected along the 𝑧 -direction for the density (left), temperature (centre), and ionization fraction (right) at the end ( 𝑧 =
6) ofthe spx2.5-c simulation.
Figure 12.
Mass-weighted maps projected along the 𝑧 -direction for 𝐵 𝑧 (top)and magnetic field strength (bottom) at the end ( 𝑧 =
6) of the spx2.5-csimulation. (Fig. 14) shows that the typical temperature in the shocked regionsis 𝑇 S = K at the mean density in the Universe, which, using hereagain Eq. (13), translates into an amplitude of 𝐵 high (cid:39) − G at¯ 𝑛 (cid:39) − cm − . Figure 13.
Mass-weighted 2D histograms of the magnetic field versus densityat the end ( 𝑧 =
6) of the spx2.5-n (left) and spx2.5-c (right) simulations. Thetwo analytical 𝑛 / relations are present in red dotted lines. The previous two analytical estimates are then used to normalisethe two branches as 𝐵 = 𝐵 low ( 𝑛 / ¯ 𝑛 ) / and 𝐵 = 𝐵 high ( 𝑛 / ¯ 𝑛 ) / .We show the resulting predictions as the two red dashed lines inFig. 13. We see that the 𝑛 / scaling is not strictly reproduced byour simulation results. This is usually attributed to turbulent pro-cesses providing additional magnetic field amplification, on top ofbare gravitational contraction (Dubois & Teyssier 2008; Quilis et al.2020). This relatively weak effect can be interpreted as the very firstsign of a galactic dynamo (see Discussion below).In order to demonstrate that the high- 𝐵 branch is indeed generatedby SN-driven galactic winds, we run the same simulation but dis-abling SN feedback (spx2.5-c-nosn). The hot magnetised bubblesthat are previously seen around massive haloes completely disap-pear (Fig. 15). We confirm this by inspecting the phase-space dia-gram, from which the high- 𝐵 branch also disappears, only leavingthe ionization branch and the linear regime. This confirms that thisnew magnetic field generation channel directly originates from SN-driven galactic winds. Note that in the early simulations of Gnedinet al. (2000), galactic winds were not as strong as they are now in ournewest cosmological models, since they did not include explosions MNRAS , 1–14 (2021) O. Attia et al.
Figure 14.
Mass-weighted 2D histogram of the temperature versus density atthe end ( 𝑧 =
6) of the spx2.5-c simulation. of massive stars. As a consequence, this magnetic field generationchannel was missing.
A well-known limitation of radiation hydrodynamics simulations ofthe EoR is their sensitivity to the finite size of the computationalbox (e.g. Iliev et al. 2006). The presence (or absence) of large andrare haloes can significantly affect the reionization history of thesimulated Universe. To understand the sensitivity of the Biermannbattery to these finite volume effects, we have performed the spx5.0simulation, whose volume is 8 times larger than our fiducial spx2.5simulations. At 𝑧 =
6, the 5 cMpc box is more efficiently ionizedthan the 2.5 cMpc box (volume-averaged ionized hydrogen fractionof 𝑥 H ii = .
56 for the former versus 𝑥 H ii = .
19 for the latter). Thereionization of the 5 cMpc volume is significantly slower than inRosdahl et al. (2018). This is primarily due to the updated versionof the Bpass SED model we use, which assumes a smaller fractionof stars living in binary systems and leads to lower escape fractionsof ionizing radiation from galaxies. This will be addressed in moredetail in a forthcoming paper (Rosdahl et al., in prep ).The spatial distribution of the magnetic field intensity (Fig. 16,top panel) qualitatively shows the same properties as the spx2.5-csimulation. A few neutral regions remain with a weak magnetic fieldcorresponding to the linear regime of the Biermann battery. Most ofthe volume is filled with ionized gas at 𝑇 (cid:39) K and correspondsto the Biermann battery driven by ionization fronts, resulting in afield strength 𝐵 (cid:39) − G. Hot gas bubbles are also visible, with amagnetic field strength between 100 to 1000 times higher than themean value at that epoch. Given the larger number of galaxies inthe larger box, they are much more numerous than for our fiducialrun. Some of the bubbles are however significantly larger than in thesmaller box. They correspond to more massive galaxies with higherstar formation rates. In the right side of the map in Fig. 16, we caneven see some of the hot bubbles percolating.We thus recover the three modes of magnetic field generation bythe Biermann battery, as discussed in the previous section. This isclearly confirmed by the 𝐵 − 𝑛 histogram (bottom panel of Fig. 16). Figure 15.
Map of the mass-weighted magnetic field strength projected alongthe 𝑧 -direction (top) and 2D mass histogram of the magnetic field versusdensity (bottom) at 𝑧 = The higher branch that we attribute to SN-driven galactic winds ishowever one order of magnitude higher than for the 2.5 cMpc volume.We believe that we witness the onset of a turbulent dynamo insidethe largest galaxies of our simulated volume. Our increased spatialresolution allows us to start resolving a small-scale dynamo within theturbulent ISM (as in Martin-Alvarez et al. 2018). Indeed, the biggesthaloes tend to be better resolved in the bigger volume, because of theirmore massive DM component, and the highly-resolved regions areof bigger spatial extent there. Even if the resolution is not sufficientlygood to see a clear exponential growth of the field (see Rieder &Teyssier 2016, 2017, for a more detailed discussion), we start to seea weak amplification of the field within the galaxy, which is thenre-distributed on larger scales by the galactic winds. These windspollute the IGM with metals, as well as amplified magnetic fields(Bertone et al. 2006; Butsky et al. 2016; Pakmor et al. 2020). Thisexplains why the upper branch follows a similar 𝐵 ∝ 𝑛 / relation,but with a higher normalisation than the 2.5 cMpc box. It shouldbe noted that a higher temperature in the larger volume could alsoexplain the showcased magnetic field amplification (as per Eq. 13),but we have ruled out this possibility by checking that the spx5.0 𝑇 − 𝑛 histogram is similar to the spx2.5-c one (Fig. 14).In order to better understand the difference in the field generation MNRAS000
Map of the mass-weighted magnetic field strength projected alongthe 𝑧 -direction (top) and 2D mass histogram of the magnetic field versusdensity (bottom) at 𝑧 = The higher branch that we attribute to SN-driven galactic winds ishowever one order of magnitude higher than for the 2.5 cMpc volume.We believe that we witness the onset of a turbulent dynamo insidethe largest galaxies of our simulated volume. Our increased spatialresolution allows us to start resolving a small-scale dynamo within theturbulent ISM (as in Martin-Alvarez et al. 2018). Indeed, the biggesthaloes tend to be better resolved in the bigger volume, because of theirmore massive DM component, and the highly-resolved regions areof bigger spatial extent there. Even if the resolution is not sufficientlygood to see a clear exponential growth of the field (see Rieder &Teyssier 2016, 2017, for a more detailed discussion), we start to seea weak amplification of the field within the galaxy, which is thenre-distributed on larger scales by the galactic winds. These windspollute the IGM with metals, as well as amplified magnetic fields(Bertone et al. 2006; Butsky et al. 2016; Pakmor et al. 2020). Thisexplains why the upper branch follows a similar 𝐵 ∝ 𝑛 / relation,but with a higher normalisation than the 2.5 cMpc box. It shouldbe noted that a higher temperature in the larger volume could alsoexplain the showcased magnetic field amplification (as per Eq. 13),but we have ruled out this possibility by checking that the spx5.0 𝑇 − 𝑛 histogram is similar to the spx2.5-c one (Fig. 14).In order to better understand the difference in the field generation MNRAS000 , 1–14 (2021) he Biermann Battery during the EoR Figure 16.
Map of the mass-weighted magnetic field strength projected alongthe 𝑧 -direction (top) and 2D mass histogram of the magnetic field versusdensity (bottom) at 𝑧 = and amplification processes between the two volumes, we look atthe temporal evolution of the volume-weighted average magneticfield strength and the volume-weighted average temperature for ourdifferent simulations (Fig. 17). This allows us to define three differentepochs of field generation and amplification.The first epoch corresponds to the linear regime of the Biermannbattery. In our simulations, it ends around 𝑧 =
18 (resp. 𝑧 = 𝐵 (cid:39) − G on scales of ∼
20 ckpc (Naoz & Narayan2013). In our simulations, these scales are resolved and the thermo-dynamics of the gas are properly captured, which explains why wecan observe this pre-EoR magnetic field generation as predicted byNaoz & Narayan (2013). As it can also be noticed in Naoz & Narayan(2013), the average field strength remains roughly constant duringthe linear regime. Indeed, this equilibrium value stems from a com-petition between the Biermann process and the cosmic expansion, for
Figure 17.
Evolution of the volume-averaged magnetic field strength (top)and volume-averaged temperature (bottom) as a function of the scale factor(or redshift) for the different cosmological simulations. which 𝐵 would otherwise decay as 𝑎 − / . We can also see in Fig. 17that the naive method creates too much spurious field to properlycapture this process.The second epoch corresponds to the onset of cosmic reionization.Our simulation without SN feedback only captures the reionizationprocess, with its corresponding thermal and magnetic history (seeFig. 17). The mean temperature in the box steadily rises as Strömgrenspheres progressively percolate. The average magnetic field quicklyrises after reionization starts and stabilises to an average field strengthof 𝐵 (cid:39) − G. The third epoch corresponds to the onset of galacticwinds, when the first SNe explode in the star-forming galaxies. It canbe seen as a series of discrete events in the temperature evolution. Theresulting average magnetic field is two orders of magnitude higherthan in the no SNe case. In these two regimes, both the naive andthe correct methods predict very similar average field strengths. Theeffect of the box size, although relatively weak (if present at all) forthe average temperature evolution, is quite important for the averagemagnetic field strength evolution, with an increase close to one orderof magnitude in strength. As explained earlier, we speculate thatthe large box simulation is marginally resolving a turbulent dynamoat small scales, within the turbulent ISM of the largest galaxies,amplifying the field 10 × more than the small box.This conclusion is further supported by the comparison of thevolume-weighted PDF of the magnetic field for the two different boxsizes at 𝑧 = 𝐵 ≤ − Gand a larger volume fraction corresponding to the ionization channelof the Biermann battery with 𝐵 (cid:39) − G. Second, both PDFsshow a smaller peak at higher field strength with 𝐵 ≥ − G. Thissecondary peak corresponds to the galactic winds with a smallervolume fraction. The small box features for this secondary peak amaximum value around 𝐵 (cid:39) − G, while the large box peak isclearly shifted to 𝐵 (cid:39) − G. This corresponds to a more efficient
MNRAS , 1–14 (2021) O. Attia et al.
Figure 18.
Volume-weighted PDF of the magnetic field strength at 𝑧 = dynamo process in the largest galaxies and perfectly explains theaverage magnetic field strength temporal evolution in Fig. 17. An important aspect of this paper is the design of a robust numer-ical method for the Biermann battery. We have implemented andtested both a simple but naive approach and a more elaborate butaccurate approach that preserves the absence of Biermann batteryfor isothermal flows. This structure preserving property proves to beefficient in two important test cases, namely the Strömgren sphereand the Sedov blast wave, in reducing spurious magnetic fields andaccelerating the convergence of the numerical solution. In both tests,however, the presence of sharp discontinuities constitutes an impor-tant difficulty for how fast our results are converging, even for thecorrect method. One could certainly argue that the Biermann battery,as described by Eq. (2), is not a well-posed problem in the presenceof such density and temperature discontinuities. We know that infact the electron temperature 𝑇 e is actually smooth in shocks (Spitzer1962; Zel’dovich & Raizer 1967). Indeed, because of electronic heatconduction, 𝑇 e exhibits a thermal precursor ahead of the shock thatmakes the temperature profile continuous. This mechanism operateson characteristic scales close to the electron mean free path, with(e.g. Smith et al. 2013) 𝜆 e (cid:39) ( 𝑘 B 𝑇 ) 𝑒 𝑛 e ln Λ , (15)where ln Λ is the Coulomb logarithm. In our simulation, using for 𝑛 e the background density at 𝑧 = 𝑇 e a temperature of 10 K,(see Sect. 4.2), we find 𝜆 e (cid:39)
70 pc. These scales are much smallerthan the resolution of our base grid (a few kpc at 𝑧 = 𝜆 e is not resolved (see asimilar discussion in Graziani et al. 2015). The analysis of our testsuite and our cosmological simulations suggest that the accuracy ofthe correct method is high enough to properly capture the Biermannbattery process during the EoR.Our cosmological simulations indeed confirm previous findings,namely that the Biermann battery is already active during the darkages, with the linear regime of density and temperature fluctuationsleading to a field strength of 𝐵 (cid:39) − G. We also recover theclassical regime of the Biermann battery induced by the propagationof ionization fronts (Subramanian et al. 1994; Gnedin et al. 2000)with 𝐵 (cid:39) − G. Our simulations also revealed a third magnetic field generation channel within SN-driven galactic winds. In thischannel, we observe a combination of the Biermann battery withinshocks and a turbulent dynamo inside galactic discs. The resultingmagnetic field is between two and three orders of magnitude strongerthan the field generated in ionization fronts. Note that one expectsto better and better resolve the turbulent dynamo process as theresolution is increased, so that the magnetic field strength injectedinto the IGM by galactic winds is likely to increase until the turbulentdynamo reaches saturation within the star-forming galaxies (Rieder& Teyssier 2017; Aramburo Garcia et al. 2021).It is worth mentioning that the Biermann battery within shockwaves has already been proposed by Kulsrud et al. (1997) andDavies & Widrow (2000). In these earlier simulations, the effectsare triggered by non-radiative accretion shocks around haloes andare observed at relatively low redshifts of 𝑧 <
2. The correspondingBiermann generated field strength is only 𝐵 (cid:39) − G. Owing tothe much stronger SN feedback recipe we used here, the shocks weobserve in galactic winds are much stronger than accretion shocks,and therefore generate stronger magnetic fields. In fact, the virialtemperature is on the order of 10 K for the biggest haloes, as con-firmed by the 𝑇 − 𝑛 diagram (Fig. 14). SN-driven galactic winds aretypically hotter, which explains that our shock channel is strongerthan classical accretion shocks.Our results qualitatively match the recent study of Garaldi et al.(2021), in particular the simulations they have performed with theBiermann battery. At 𝑧 =
6, the field strengths they obtain are all (cid:46) − G (see their Fig. 8), which corresponds to the linear regimeof the Biermann battery. Surprisingly, no fields higher than 10 − Gappear in their simulations during the EoR but only after, once theturbulent dynamos inside galaxies saturate. The time evolution oftheir volume-weighted magnetic field (see their Fig. 9) is consistentwith this picture, with an average field strength quite stable around10 − G at high redshift. We interpret this striking difference asbeing due to the vastly different box sizes and corresponding massresolution between our study and theirs. It is indeed of primaryimportance to resolve in this context all the star-forming mini-haloesand the associated radiation and SN feedback.We would like to stress again that we do not accurately model theturbulent dynamo process within the star-forming galaxies becausewe lack the required resolution (Rieder & Teyssier 2017). At thecurrent resolution ( ∼
10 pc), this would probably require a suitablesub-grid model to capture the turbulent dynamo at unresolved scales.This would result in a fast exponential amplification of the magneticfield inside galaxies, which would then be expelled and re-distributedon larger scales by galactic winds (Rieder & Teyssier 2016; Martin-Alvarez et al. 2020a; Steinwandel et al. 2020). This exponentiationwould then consequently destroy all memory of the Biermann batterygenerated field in the regions it pollutes. However, this phenomenonis limited by the spatial extension of the outflows. Hence, it might stillbe possible to observe imprints of the primeval Biermann-generatedfields in voids where pristine gas has not been affected by galacticwinds (Bertone et al. 2006; Aramburo Garcia et al. 2021).The resulting fields from our simulations are not high enough tohave an influence on gas dynamics (negligible magnetic energy infront of the total energy). In particular, they do not affect the reion-ization process as they have no impact on star formation for suchlow values (Sethi & Subramanian 2005; Marinacci & Vogelsberger2016; Katz et al. 2021). They however provide a natural origin forseed magnetic fields subsequently amplified by the galactic dynamo.Interestingly, this process occurs before (Naoz & Narayan 2013),during, and after the EoR. Other seeding mechanisms have been pro-posed like the Durrive battery (Durrive & Langer 2015; Durrive et al.
MNRAS000
MNRAS000 , 1–14 (2021) he Biermann Battery during the EoR We introduce in this paper a novel numerical implementation of theBiermann battery within the AMR cosmological simulation codeRamses. The structure-preserving nature of our scheme has the ad-vantage of ensuring that no spurious magnetic field is generated for anisothermal flow. This property, together with the strict divergence-free nature of the CT-based ideal MHD solver in Ramses, provesparticularly accurate for several test cases and a suite of cosmologi-cal simulations of the EoR.Our simulation suite allowed us to study the impact of various keyparameters such as the choice of the Biermann battery scheme, thepresence of SN feedback and box size. In particular, we simulate twodifferent large-scale volumes: 2.5 cMpc and 5 cMpc from 𝑧 =
150 to 𝑧 =
6. The main results can be summarised as follows. • Prior to the EoR, a small magnetic field with strength 𝐵 (cid:39) − G field is generated by temperature and density fluctuationsinduced by linear perturbations. The average field strength remainsconstant over time due to a competition between the Biermann batteryand the expansion of the Universe. This process, predicted by Naoz& Narayan (2013), is properly captured only by the correct method. • During the EoR, propagating ionization fronts are able to mag-netise a substantial fraction of the volume with an average fieldstrength close to 𝐵 (cid:39) − G (as in Gnedin et al. 2000). • On top of this volume-filling field, SN-driven galactic windsaround star-forming galaxies generate hot gas bubbles hostingstronger magnetic fields with 𝐵 (cid:39) − G around massive haloes.They are due to SN-driven shocks that are subsequently expelled intothe surrounding IGM by galactic winds. Their higher strength resultsfrom the higher temperature and the more efficient Biermann battery.This channel completely disappears if SN feedback is disabled. • The resulting bi-modal distribution of the magnetic field at theend of the EoR ( 𝑧 =
6) is confirmed by a phase-space analysis of 𝐵 − 𝑛 histograms with two distinct branches both scaling as 𝑛 / butdifferently normalised with approximately a two orders of magnitudeoffset. • The higher branch that corresponds to SN-driven galactic windshas a significantly higher normalisation in the larger box. We interpretthis finite volume effect as a better resolved turbulent dynamo inthe larger galaxies. The corresponding amplified field is then re-distributed on larger scales by the galactic winds. • The morphology of the generated fields is complex, regularly al-ternating between positive and negative values. For ionization frontsas well as for galactic winds, the magnetic field tends to be mostlytoroidal, circling around the regions where it is produced. • The time evolution of the magnetic field is roughly consistent forall simulations. After the linear regime during the dark ages wherethe average field is almost constant, its time evolution reveals twomajor phases: the first one corresponds to the beginning of the EoRand the second one to the onset of galactic winds. Again, the secondepoch is absent if SN feedback is disabled. • The volume-integrated magnetic field at 𝑧 = ACKNOWLEDGEMENTS
This work was performed with the support of the Swiss NationalScience Foundation (SNSF). TK was supported by the NationalResearch Foundation of Korea (NRF-2019K2A9A1A06091377 andNRF-2020R1C1C1007079). SMA acknowledges support by the Eu-ropean Research Council (ERC) Starting Grant 638707 ‘Black holesand their host galaxies: co-evolution across cosmic time’. The simu-lations in this work were carried out on Piz Daint at the Swiss Super-computing Center (CSCS) in Lugano, and the analysis was conductedwith equipment maintained by the Service and Support for Science IT(S3IT), University of Zurich. We made use of the Pynbody package(Pontzen et al. 2013). We acknowledge the PRACE Research Infras-tructure (PRACE project ID 2016153539) for providing us accessto the SuperMUC computing resource based in Garching, Germany,to tune and perform the original Sphinx simulations on which thiswork expands.
DATA AVAILABILITY
The data underlying this article will be shared upon reasonable re-quest to the corresponding author.
REFERENCES
Aramburo Garcia A., Bondarenko K., Boyarsky A., Nelson D., Pillepich A.,Sokolenko A., 2021, arXiv e-prints, p. arXiv:2101.07207Bertone S., Vogt C., Enßlin T., 2006, MNRAS, 370, 319Biermann L., 1950, Zeitschrift Naturforschung Teil A, 5, 65Bisnovatyi-Kogan G. S., Ruzmaikin A. A., Syunyaev R. A., 1973, Soviet Ast.,17, 137Brandenburg A., Subramanian K., 2005, Phys. Rep., 417, 1Butsky I., et al., 2016, MNRAS, 462, 663Davies G., Widrow L. M., 2000, ApJ, 540, 755Deparis N., Aubert D., Ocvirk P., Chardin J., Lewis J., 2019, A&A, 622,A142Dubois Y., Teyssier R., 2008, A&A, 482, L13Durrer R., Neronov A., 2013, A&ARv, 21, 62Durrive J. B., Langer M., 2015, MNRAS, 453, 345Durrive J.-B., Tashiro H., Langer M., Sugiyama N., 2017, MNRAS, 472,1649Eldridge J. J., Izzard R. G., Tout C. A., 2007, MNRAS, 384, 1109Evans C. R., Hawley J. F., 1988, ApJ, 332, 659Fan X., Carilli C. L., Keating B., 2006, ARA&A, 44, 415Fatenejad M., et al., 2013, High Energy Density Physics, 9, 172Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B.,Verner E. M., 1998, PASP, 110, 761Finkelstein S. L., et al., 2019, ApJ, 879, 36Fromang S., Hennebelle P., Teyssier R., 2006, A&A, 457, 371Fryxell B., et al., 2010, High Energy Density Physics, 6, 162MNRAS , 1–14 (2021) O. Attia et al.
Garaldi E., Pakmor R., Springel V., 2021, MNRAS,Gnedin N. Y., Ferrara A., Zweibel E. G., 2000, ApJ, 539, 505Graziani C., Tzeferacos P., Lee D., Lamb D. Q., Weide K., Fatenejad M.,Miller J., 2015, ApJ, 802, 43Graziani C., Tzeferacos P., Lee D., Lamb D. Q., Weide K., Fatenejad M.,Miller J., 2016, in Journal of Physics Conference Series. p. 012018,doi:10.1088/1742-6596/719/1/012018Hahn O., Abel T., 2011, MNRAS, 415, 2101Iliev I. T., Mellema G., Pen U. L., Merz H., Shapiro P. R., Alvarez M. A.,2006, MNRAS, 369, 1625Katz H., Kimm T., Sijacki D., Haehnelt M. G., 2017, MNRAS, 468, 4831Katz H., et al., 2021, arXiv e-prints, p. arXiv:2101.11624Kimm T., Cen R., 2014, ApJ, 788, 121Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017,MNRAS, 466, 4826Kroupa P., 2001, MNRAS, 322, 231Kulkarni G., Keating L. C., Haehnelt M. G., Bosman S. E. I., Puchwein E.,Chardin J., Aubert D., 2019, MNRAS, 485, L24Kulsrud R. M., Cen R., Ostriker J. P., Ryu D., 1997, ApJ, 480, 481Lazar M., Smolyakov A., Schlickeiser R., Shukla P. K., 2009, Journal ofPlasma Physics, 75, 19Marinacci F., Vogelsberger M., 2016, MNRAS, 456, L69Martel H., Shapiro P. R., 1998, MNRAS, 297, 467Martin-Alvarez S., Devriendt J., Slyz A., Teyssier R., 2018, MNRAS, 479,3343Martin-Alvarez S., Katz H., Sijacki D., Devriendt J., Slyz A., 2020a, arXive-prints, p. arXiv:2011.11648Martin-Alvarez S., Slyz A., Devriendt J., Gómez-Guijarro C., 2020b, MN-RAS, 495, 4475Naab T., Ostriker J. P., 2017, ARA&A, 55, 59Naoz S., Barkana R., 2005, MNRAS, 362, 1047Naoz S., Narayan R., 2013, Phys. Rev. Lett., 111, 051303Ocvirk P., Aubert D., Chardin J., Deparis N., Lewis J., 2019, A&A, 626, A77Pakmor R., et al., 2020, MNRAS, 498, 3125Planck Collaboration 2014, A&A, 571, A1Pontzen A., Roškar R., Stinson G., Woods R., 2013, pynbody: N-Body/SPHanalysis for python (ascl:1305.002)Quilis V., Martí J.-M., Planelles S., 2020, MNRAS, 494, 2706Rees M. J., 2005, Magnetic Fields in the Early Universe. Wielebinski, Richardand Beck, Rainer, doi:10.1007/11369875_1Rieder M., Teyssier R., 2016, MNRAS, 457, 1722Rieder M., Teyssier R., 2017, MNRAS, 472, 4368Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436,2188Rosdahl J., et al., 2018, MNRAS, 479, 994Rosen A., Bregman J. N., 1995, ApJ, 440, 634Sethi S. K., Subramanian K., 2005, MNRAS, 356, 778Smith B., O’Shea B. W., Voit G. M., Ventimiglia D., Skillman S. W., 2013,ApJ, 778, 152Spitzer L., 1962, Physics of Fully Ionized Gases. New York: IntersciencePublishersStanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75Steinwandel U. P., Dolag K., Lesch H., Moster B. P., Burkert A., Prieto A.,2020, MNRAS, 494, 4393Subramanian K., Narasimha D., Chitre S. M., 1994, MNRAS, 271, L15Teyssier R., 2002, A&A, 385, 337Teyssier R., Fromang S., Dormy E., 2006, Journal of Computational Physics,218, 44Toro E. F., 2013, Riemann solvers and numerical methods for fluid dynamics:a practical introduction. Springer Science & Business MediaTóth G., et al., 2012, Journal of Computational Physics, 231, 870Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS,470, 224Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506,647van Leer B., 1997, Journal of Computational Physics, 135, 229Wise J. H., Turk M. J., Norman M. L., Abel T., 2012, ApJ, 745, 50Yee K., 1966, IEEE Transactions on Antennas and Propagation, 14, 302 Zel’dovich Y. B., Raizer Y. P., 1967, Physics of shock waves and high-temperature hydrodynamic phenomena. DoverThis paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000