A small-scale dynamo in feedback-dominated galaxies as the origin of cosmic magnetic fields - I. The kinematic phase
The origin and evolution of magnetic fields in the Universe is still an open question. Their observations in galaxies suggest strong magnetic fields already at high redshift as well as at present time. However, neither primordial magnetic fields nor battery processes can account for such high field strengths, which implies the presence of a dynamo process with rapid growth rates in high-redshift galaxies and subsequent maintenance against decay.
We investigate the particular role played by feedback mechanisms in creating strong fluid turbulence, allowing for a magnetic dynamo to emerge. Performing magnetohydrodynamic simulations of isolated cooling gas halos, we compare the magnetic field evolution for various initial field topologies and various stellar feedback mechanisms. We find that feedback can indeed drive strong gas turbulence and dynamo action. We see typical properties of Kolmogorov turbulence with a kinetic energy spectrum, as well as a small-scale dynamo, with a magnetic energy spectrum predicted by Kazantsev dynamo theory. We also investigate simulations with a final quiescent phase. As turbulence decreases, the galactic fountain settles into a thin, rotationally supported disk. The magnetic field develops a large-scale, well-ordered structure with even symmetry, which is in good agreement with magnetic field observations of nearby spirals. Our findings suggest that weak initial seed fields were first amplified by a small-scale dynamo during a violent, feedback-dominated early phase in the galaxy formation history, followed by a more quiescent evolution, where the fields have slowly decayed or were maintained via large-scale dynamo action.
keywords:galaxies: magnetic fields - methods: numerical - MHD - turbulence
Large-scale magnetic fields are measured with strengths of up to several G in nearby galaxies (Beck et al., 1996), and possibly even higher field strength have been detected in earlier galaxies at high redshift (Bernet et al., 2008). The preferred theory to explain their origin is based on the early generation of seed fields at the epoch of cosmic re-ionisation, through the microscopic process known as ”Biermann battery” (Biermann, 1950; Naoz & Narayan, 2013), later amplified during the galaxy formation era, through the large-scale galactic dynamo process (Parker, 1970; Brandenburg & Subramanian, 2005). The Biermann battery is likely to generate fields at the level of G, leaving to the galactic dynamo process more than 14 orders of magnitude of field amplification during the 10 Gyr of cosmic evolution. The challenge for galactic dynamos is even more severe, if one considers that strong fields are already in place at high redshift (Widrow, 2002), and are probably even stronger than they are today (Bernet et al., 2008).
Successful theoretical models for large-scale galactic dynamos report exponential growth rates of the order of to , where is the galactic disk rotation rate (Pariev et al., 2007). For typical, present day spirals, this translates into e-folding amplification time scale of roughly 1 Gyr, making the task of amplifying the field over 14 orders of magnitude virtually impossible. One noticeable exception is the cosmic-ray-driven dynamo proposed by Parker (1992) and simulated by Hanasz et al. (2004), leading to a measured growth rate , although the numerical experiment was performed over only a relatively limited time, since the reported magnetic field amplification was over only 3 orders of magnitude (Hanasz et al., 2004). On the theoretical side, classical mean field dynamos are plagued by the catastrophic -quenching effect, leading to very low saturated values for the large-scale magnetic field (Kulsrud & Anderson, 1992; Vainshtein & Cattaneo, 1992), owing to the strict conservation of magnetic helicity in a closed system. A possible solution to this problem is the effect of galactic winds, that could drag the magnetic field lines outside of the dynamo-active disk, therefore alleviating the aforementioned quenching issue (Del Sordo et al., 2013).
The theory of galaxy formation has significantly evolved over the past decade, with the ever increasing role of feedback processes (Scannapieco et al., 2012) and their associated galactic winds (Oppenheimer & Davé, 2006), together with the dominant accretion mechanism through cold streams (Kereš et al., 2005; Ocvirk et al., 2008; Dekel et al., 2009). On the observational side, galactic winds are indeed ubiquitous in star bursting local galaxies (Martin, 1999), but also in many “normal” high redshift galaxies (Steidel et al., 2010). One of the most spectacular observational constraints on galaxy formation theories was obtained by matching the stellar mass of the central galaxies to their parent halo mass (Behroozi et al., 2013; Moster et al., 2013). This has led theorists to consider much stronger feedback processes, in order to regulate star formation throughout cosmic time, especially at high redshift, when the star formation efficiency was so low (Agertz et al., 2013; Hopkins et al., 2014; Roškar et al., 2014; Wang et al., 2015).
In this rather violent, feedback-dominated scenario, dwarf galaxies play a very important role. They are the dominant galaxy population at high redshift, probably responsible for the cosmic re-ionisation (Kimm & Cen, 2014). They are also the progenitors of the Milky Way satellites, which are useful laboratories to test our current galaxy formation paradigm. For the latter, violent feedback mechanisms have also been invoked to explain the absence of cusp in the dark matter density profile, and the presence of a dark matter core in low surface brightness galaxies (de Blok et al., 2001). Cosmological simulations of dwarf galaxies have been performed with strong feedback recipes, confirming in this case the formation of a dark matter core (Governato et al., 2010, 2012). Recently, we have also performed idealised simulations of an isolated, cooling gaseous dwarf halo, obtaining, in this well-controlled numerical experiment, the formation of a dark matter core (Teyssier et al., 2013). The dark matter core formation mechanism is now well understood (Pontzen & Governato, 2012). It is due to repeated, energetic feedback events due to many supernovae explosions, leading to violent oscillations of the gravitational potential, due to the large gas mass variations within the central kilo parsec of the galaxy. A possible observational signature of this effects is a typical, oscillatory star formation history, mimicking a breathing mode in the gas distribution (Kauffmann, 2014).
In this paper, we want to study the impact of a strong feedback scenario on the growth of magnetic fields in dwarf, as well as in larger galaxies. The velocity field on both small and large scales, resulting from repeated giant feedback events, can have a direct influence on the growth of the magnetic energy. Indeed, supernovae explosions in the Milky Way have been considered for quite a long time as a source of helical gas motions, promoting a large-scale -dynamo in the Galaxy (Ferriere, 1992). The Milky Way is however a rather quiescent galaxy, with moderate supernovae activity. In this paper, we are considering feedback-dominated galaxies, with high star formation rates and violent turbulent motions, together with large-scale galactic fountains or winds.
Several simulations including magnetic fields have been performed recently in the context of galaxy formation (Wang & Abel, 2009; Dubois & Teyssier, 2010). These simulations, based on the popular “cooling halo” numerical set-up, have achieved only moderate magnetic field amplification. The important property of these simulations is the absence of feedback (Wang & Abel, 2009), or the relative weakness of the feedback recipe used at that time (Dubois & Teyssier, 2010).
A first exception is the simulation reported in Beck et al. (2012), based on a MHD version of the SPH code GADGET with divergence cleaning. They observed a fast exponential growth of the magnetic field, which they attribute to a small-scale dynamo. Feedback processes were included through an effective Equation-of-State (EoS), without any explicit source of turbulence in these relatively smooth, thermally-supported flows. These authors however reported very strong growth rates, with e-folding times as small as 10 Myr, although analytical estimates based on small-scale dynamo theory predicted e-folding times closer to 100 Myr.
A second exception is the recent simulation reported in Pakmor & Springel (2013), using the new Magneto-Hydrodynamics (MHD) solver developed for the AREPO code (Pakmor et al., 2011), where strong magnetic field amplification has also been observed, although, here again, stellar feedback effects were not considered explicitly, but only indirectly as a modified thermal EoS, leading to the formation of relatively smooth, two-dimensional flows, in which dynamo amplification is in principle notoriously difficult to obtain.
In the present paper, we will use a similar set-up as in all those previous studies, namely a cooling isolated gaseous halo, considering simulations with (but also without) strong stellar feedback. We will use the Adaptive Mesh Refinement code RAMSES (Teyssier, 2002), adopting the “Constrained Transport”, strictly divergence-free-preserving, MHD solver presented in Teyssier et al. (2006) and in Fromang et al. (2006). The paper is organised as follows: In 2, we will present our numerical methods, both in terms of galaxy formation physics and magnetic field modelling. In 3, we describe our initial conditions for the isolated, magnetised cooling halo. In 4, we present our main results, outlining the difference between the feedback and the no-feedback cases. Finally, in 5, we discuss the implications of our results in the context of galactic dynamo theory, as well as possible further studies to confirm and broaden our findings.
2 Numerical methods
We have performed MHD simulations of isolated, cooling haloes, using the Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier, 2002). These simulations feature a collisionless fluid (for dark matter and stars) and a magnetised gaseous component, coupled through gravity. In this section, we describe the simulation technique used to follow the evolution of our isolated halo. First, we describe in details our AMR implementation for solving the ideal MHD equations, together with simple test cases to show that it works as intended in the context of galactic dynamo. We then describe the adopted galaxy formation physics, such as gas cooling, metal enrichment, star formation and stellar feedback, leading to what we believe to be a realistic model of the interstellar medium (ISM).
2.1 Ideal MHD solver
The equations that we solve are the ideal MHD equations (written here without gravity and cooling source terms for the sake of simplicity)
where is the gas density, is the momentum, is the magnetic field, is the total energy, and is the internal energy. The total pressure is given by where we assume a perfect gas equation of state . The system of equations is completed by the soloinoidal constraint
The code is grid-based with a tree-based adaptively refined mesh. The equations are solved using the second-order unsplit Godunov scheme based on the MUSCL-Hancock method. We chose the HLLD Riemann solver with the MinMod slope limiter for the hydro variables which are cell-centred. The magnetic field on the other hand is treated as a face-centered variable. This allows the use of the Constrained Transport (CT) method to advance the induction equation ((4)) in time, which preserves the divergence of the magnetic field to the numerical precision level (Teyssier et al., 2006). The CT method involves a spatial interpolation of the EMF on the cell edges for the predictor step and solving a 2D Riemann problem for the corrector step. For the 2D problem, we use the HLLD solver as well and for the magnetic field in general, the MonCen slope limiter.
Boundary conditions were chosen to allow for free outflow. For the 5 hydro variables, this is done by imposing a vanishing gradient at the domain boundary (zero-gradient method). The same can be applied to the transverse magnetic field component parallel to the boundary face , but would cause a non-zero divergence of the magnetic field if applied to the normal component which is perpendicular to the face. Instead, we use a linear interpolation for so that . Note that this method can cause an inward Poynting flux which transports magnetic energy from the outside into the computational domain. Since the magnetic field at the border is many orders of magnitude weaker than the average, this does not contribute significantly to the overall magnetic energy evolution (see Dubois & Teyssier, 2010).
Special care needs to be taken also when refining and de-refining cells, in order to enforce the constraint, when interpolating the magnetic field. A solution to this problem within the CT framework has been proposed by Balsara (2001) and Tóth & Roe (2002), and we adopt it here for newly refined cells, but also for temporary ghost cells used to set proper boundary conditions at coarse-fine level boundaries.
In the context of galactic dynamos, it is worth mentioning that our code has been tested extensively against well-known flows triggering fast dynamos, such as in the ABC flow (Galloway & Frisch, 1986; Childress & Gilbert, 1995) or in the Ponomarenko dynamo (Ponomarenko, 1973). We have shown in Teyssier et al. (2006) that our numerical scheme for the ideal MHD equations was in fact slightly resistive, with, for a regular Cartesian grid, a numerical magnetic Reynolds number roughly inversely proportional to the square of the number of grid points per box length. This scaling is to be expected for second-order schemes and smooth flows. In the context of AMR and highly complex, turbulent flows, determining the exact effective numerical Reynolds number of the simulated flow is impossible. Qualitatively, though, it is important to bear in mind that magnetic reconnection and other diffusive processes occur in the simulation at a typical scale probably very close to the grid scale. This scale plays a very important role in dissipating the kinetic energy of the turbulence, and also controls the magnetic energy losses due to reconnection or (numerical) Ohmic dissipation.
2.2 Cooling and star formation
In addition to solving the ideal, self-gravitating MHD equations, we also include many physical processes relevant to galaxy formation. One of the key physical ingredient is gas cooling, which leads the hot, initially hydrostatic halo gas to loose pressure support and to condense in the centre as a centrifugally supported disc. When this atomic gas of K is allowed to cool even more due to fine-structure metal line cooling or molecular cooling, the disc fragments into dense clumps, leading to the formation of a turbulent, multiphase medium. To model gas cooling, we use standard H and He cooling processes with an additional contribution from metals based on Sutherland & Dopita (1993) for temperatures above K and metal fine-structure cooling below K, as in Rosen et al. (1995). The metallicity, denoted as , is modelled as a passive scalar, representing the mass fraction of atoms heavier than Helium in the gas. It is initialised to in the halo, mimicking molecular Hydrogen cooling in a zero metallicity gas. The metallicity is increased further by supernova feedback events, using a metal yield of .
Our refinement strategy is based on a quasi-Lagrangian approach: each cell is refined if it contains more than 8 dark matter particles or if its baryonic mass (including gas and star particle mass) exceeds , where is the adopted mass resolution of the simulation. Refinement are performed recursively, on a cell–by–cell basis, until the adopted maximum level of refinement is reached (noted ). It is crucial for astrophysical simulations to resolve spatially the Jeans length (Truelove et al., 1997). Requiring that the Jeans mass is resolved by at least mass resolution elements, , and adopting a realistic minimum temperature for the gas, noted , one can compute the corresponding Jeans length, and require it to be resolved by 4 cells, . We can then determine the maximum required level of refinement corresponding to the adopted mass resolution . To prevent the gas from accumulating and locally violating the Jeans length criterion, we also add an artificial pressure floor,
so that the gas density will never significantly exceed a typical value given by the relation .
Stars are treated as collisionless particles which are created stochastically from the gas according to a Schmidt law (as in Rasera & Teyssier, 2006)
if the local density is above a threshold density . In this paper, we always choose the star formation threshold density to be equal to the previously defined Jeans density . The star formation efficiency per free-fall time is always set to ; this value is based on observations of nearby molecular clouds (Krumholz & Tan, 2007). Creation of stellar particles is a local random Poisson process with Poisson parameter where is the simulation time-step and
the mass of the resulting stellar particle, which is equal to the smallest cell mass at the density threshold. For each simulation, precise numbers for all these parameters are given in Table 1.
2.3 Stellar feedback
In this paper, we explore the consequences of strong feedback scenarios on the amplification of the magnetic field in dwarf and Milky-Way-sized galaxies. The proper modelling of stellar feedback mechanisms, such as supernovae explosions, photo-ionised bubbles or infrared radiation in dusty environment is subject to intensive research throughout the ISM and galaxy formation literature. Understanding in details these various processes goes far beyond the scope of this paper. Our goal is merely to use various phenomenological recipes to model such feedback mechanisms very crudely, and produce dynamical properties that we believe are relevant for high-redshift galaxies, the most important one being the gas velocity field, highly turbulent, explosive and fountain-like, which could result in a fast magnetic dynamo.
For this purpose, we used a numerical model for supernovae feedback developed in the context of dwarf galaxies evolution, and that turned out to lead to the formation of a dark matter core (Teyssier et al., 2013). The main ingredient is the use of a non-thermal energy variable, and its associated pressure, treated as a passively advected scalar quantity , which represents various small-scale, non-thermal energies released by supernovae (e.g. turbulence, magnetic fields or cosmic rays). The evolution of this non-thermal energy is specified by
where the dissipation time scale is fixed to Myr and the energy injection per supernovae is set by
where the mass fraction in massive stars is set to and the local star formation rate is set by our adopted Schmidt law. For details about the implementation, we refer to Teyssier et al. (2013)
Because supernovae are probably not energetic enough to trigger strong winds in Milky-Way-sized galaxies, it has been proposed recently to consider stellar radiation as an additional feedback mechanism (Murray et al., 2010). Interstellar dust indeed absorbs UV photons, much of it being subsequently re-emitted as thermal radiation in the infrared band. This radiation will transfer momentum to the gas through radiation pressure (Murray et al., 2010; Hopkins et al., 2012; Agertz et al., 2013; Roškar et al., 2014) In this paper, we consider radiation feedback only for Milky-Way-sized galaxy simulations. We use here again a very crude model to capture the energy from the stellar UV radiation, using a simple escape probability model as
with the dust mass density is assumed to be (here, Z denotes the gas metallicity). The dust opacity at is taken to be (Draine & Li, 2007), and the total energy released during the first 10 Myr of a progenitor . The same cell is then assumed to absorb the energy in the infrared
where is the dust opacity in the IR band, which is a free parameter in our feedback implementation, usually around 10 cmg (Draine & Li, 2007; Semenov et al., 2003). The energy is added in the non-thermal energy equation for in the supernova feedback so that it contributes to . Details about the implementation can be found in Roškar et al. (2014). Here again, we would like to stress that our goal is not to study in great details realistic feedback mechanisms, but rather to generate galactic velocity fields in qualitative agreement with high-redshift galaxies and their associated strong outflows, in the context of galactic dynamos.
3 Initial conditions
We have performed a series of non-cosmological simulations of isolated halos in hydrostatic equilibrium, varying the initial set-up with two different halo sizes (a typical dwarf and a typical Milky-Way) and two different initial magnetic field topologies (dipole and quadrupole), in addition to the various options for stellar feedback that we have discussed in the previous section. We will now describe more precisely our initial set-up, and a summary of the various run parameters is given in Table 2. We refer the interested reader to Teyssier et al. (2013) for a more detailed analysis of the non-magnetohydrodynamic properties in the dwarf halo case.
3.1 Initial halo
Our initial halo follows a Navarro et al. (1997) (hereafter NFW) density profile with a concentration parameter . The smaller of the two halos, representative of a typical dwarf galaxy (we use the acronym DW), has a circular velocity of and a virial mass of , both measured at the virial radius . The halo is truncated at , resulting in the total enclosed mass of . This is essentially the same set of parameters we used in the hydrodynamic simulations of Dubois & Teyssier (2010) and Teyssier et al. (2013).
The larger halo is chosen to be a typical Milky Way galaxy (we use the acronym MW), where we increased the circular velocity to , corresponding to a virial radius of and a virial mass of . It is truncated at , resulting in the total enclosed mass of .
In all other aspects, both initial configurations follow the same prescription as in Teyssier et al. (2013). We consider a gas fraction equal to the universal mean value , and the gas density is also following a NFW profile. The gas temperature is initialised by solving the hydrostatic equilibrium equation. The gaseous halo is set in slow rotation around the z-axis, using the angular momentum profile from cosmological simulations and a spin parameter . The dark matter halo is sampled by dark matter particles, whose initial positions and velocities were computed with the density-potential pair approach of Kazantzidis et al. (2004) and Read et al. (2006). The stability of the resulting gas-dark matter equilibrium was shown in Teyssier et al. (2013) to be sufficiently good for our present purpose.
3.2 Initial magnetic field
A fundamental ingredient in any MHD simulation is the adopted initial magnetic field configuration. The simplest possible choice would be a constant field parallel to one direction, e.g. a vertical uniform field
However, we argue here that this choice is not appropriate for initially concentrated mass distributions in general and for cosmological halos in particular. This simple choice results indeed in a completely uniform magnetic energy distribution. As collapse proceeds, because of the initially peaked density distribution, less and less mass is added to the central galaxy, especially at late time. Magnetic energy, however, is still being accreted efficiently, even at late time, and artificially added to the central object. We believe it is more appropriate, in the context of ideal MHD (frozen-in magnetic flux), to consider that the initial magnetic energy follows closely the initial density distribution, in which case the magnitude of the field would scale roughly as
This requirement translates into a more complex field topology, and we need to work harder to initialise the field, compared to the uniform, vertical field case. Another important property of galactic dynamos is the field parity with respect to the system’s mid-plane. To explore possible effects related to the direction of the vector field, with odd or even parity across the mid plane, we consider initially two typical topologies: dipole-like or quadrupole-like. The dipole-like field (we use the acronym D) is defined using the following vector potential
where is the initial gas density given by the NFW profile and is the unit vector along the toroidal direction. Note that here, coordinates are given in a cylindrical coordinate system centered around the halo center and aligned with the rotation axis. The corresponding magnetic field has a vertical component which is symmetric with respect to the mid plane, while its radial component is antisymmetric. The quadrupole field (we use the acronym Q) is defined using the vector potential
and has the reversed symmetries. Its vertical component is antisymmetric and its radial component is symmetric with respect to the mid plane. Figure 1 illustrates the shape of the two vector fields. Note that in both cases, the initial field strength follows a profile peaked around the center with magnitude roughly proportional to , as expected. The magnetic field is then initialised on each cell face of the AMR grid, using a finite-difference approximation of the curl:
This ensures that the divergence of the magnetic field is initially exactly zero (to machine precision), and, thanks to the Constrained Transport method, remains zero during the entire simulation.
In the present paper, we would like to explore the purely kinematic regime of the magnetic field evolution. This corresponds to very small values of the magnetic field, for which there is no back reaction on the flow (the Lorentz force can be ignored). We therefore only solve the induction equation, which is linear with respect to the magnetic field. The exact value of the parameter is therefore irrelevant, and we will always quote magnetic field intensity as a function of the initial intensity, or as a function of the average intensity. We will study the saturation regime, and how the field reaches equipartition with the thermal and kinetic energies of the gas in a companion paper (hereafter Paper II). Close to saturation and equipartition, the exact value of the field matters a lot, and in this case, the initial intensity plays a very important role. Here, however, only the initial spatial distribution and the initial topology of the field are important, but not its overall initial normalisation.
3.3 Summary of additional physics parameters
The feedback mechanism, whose details were explained in the previous section, can be switched on and off and in case of radiative feedback its effective strength can be controlled by changing the surrounding dust opacity parameter . The dwarf halo simulations were run without any radiative feedback. The only option was to have supernova feedback (simulations dubbed ’SN’) or no feedback at all. In the Milky Way case, however, adding to those two options we also tested two more set-ups with radiation feedback. The values tested were one medium scale dust opacity value of (’K-5’) for 70 K dust and a rather opaque gas with (’K-20’) corresponding to a dust temperature of 140 K (cite Semenov et al. 2003). The feedback parameters for the dwarf as well as the Milky Way simulations are given in Table 1 and Table 2. Both configurations have the same star formation and supernova feedback efficiencies. The temperature floor used to prevent the gas from fragmenting below our resolution limit is given by
where the critical temperature T is a cool 100 K for the dwarf halo and warm 2000 K in the Milky Way case.
|Name||Topology||SN Feedback||Opacity [cm/g]|
4 Magnetic fields amplification through feedback processes
We will now present the results of our halo simulations, where we studied the influence of various stellar feedback parameters. This section is organised as follows: First, we present our dwarf galaxy simulations, without feedback, then using supernovae feedback. Second, we present the Milky-Way-sized galaxy simulations. For the latter case, supernovae feedback does not differ strongly from the no-feedback case, although it introduces slightly more turbulence in the gas. Radiation feedback makes however a big difference, and we explore two different dust opacities, resulting into two different scenarios for the galactic outflows.
4.1 Dwarf galaxy
All our simulations begin in a similar way, which is the classical scenario for these cooling halo set-up. The gas, although initially in strict hydrostatic equilibrium, loses thermal energy through radiative cooling. It thus collapses and a centrifugally supported disk forms from the inside out, thanks to the initial angular momentum profile. In the dwarf galaxy case, the disc is relatively thick at first: Atomic cooling sets a natural temperature floor around K. Low temperature radiative processes (here mostly fine-structure cooling of metals) cools the gas further, leading to the formation of a thin disc which quickly fragments into dense gas clumps. The gas density in these clumps reaches the star formation density threshold and the first stars form.
No feedback case
In absence of feedback, the disc remains very thin, and the gas clumps are long-lived. Although our star formation efficiency was set very low (one percent), most of the gas inside the dense clumps is converted into stars, after a few disc orbital time. The resulting galaxy is very efficient at transforming most of its baryons into stars, which is at odd with observed dwarf galaxies in the universe. Moreover, the resulting circular velocity profile is strongly peaked towards the centre, although dwarf galaxies circular velocity profiles are usually declining towards the centre.
The magnetic energy evolution can be seen in Figure 2. Our new simulations confirm the earlier finding of Wang & Abel (2009) and Dubois & Teyssier (2010): During the early magnetic field amplification due to the collapse of the gas in the first few Myr, the magnetic energy increases as , where is the ratio of the gas density after and before the collapse. Note that the magnetic field topology matters a lot in this early evolution. Without feedback, the dipole configuration leads to magnetic reconnection in the mid plane, as expected from the antisymmetry in the poloidal field. In the quadrupole field configuration, because of the symmetry of the poloidal field, the magnetic energy is not affected by field cancellation effects.
After the collapse, one can see in Figure 2 that the magnetic energy still grows, but much more slowly. This can be explained from field lines being twisted by the differential rotation. In this almost perfectly axisymmetric geometry, one can indeed approximate the induction equation as (see e.g. Dubois & Teyssier, 2010)
The toroidal field grows therefore only linearly with time, while the poloidal field (mostly radial) remains constant. This results in quadratic time relation of the magnetic energy
with the magnetic energy after collapse E, which depends on the field topology, and the shearing rate . We illustrate the contribution of this model for shearing amplification to for the no feedback simulations in Figure 2 (left).
On the other hand, one can see directly from (19) that a fast, exponential amplification of the field can be obtained only if the radial component grows as fast as the tangential component. For two dimensional, axisymmetric flows like our smooth rotating disk, this cannot be the case, according to the famous Zel’dovich and Cowling anti-dynamo theorems (Charbonneau & Steiner, 2012). A further inspection of Figure 2 reveals several spikes in the magnetic energy evolution. These are due to collapsing, rotating gas clumps, which trigger short episodes of field amplification. As anticipated by Wang & Abel (2009), these vortex modes do indeed amplify the field locally, but as soon as the clumps dissolve in the large-scale rotating flow, so does their magnetic energy.
These clumps also trigger three-dimensional turbulence in the gas, thanks to clump-clump interactions (Agertz et al., 2009). This could in principle increase the magnitude of the radial component of the magnetic field, but the induced effects remain too weak to affect the large-scale dynamo. Figure 7 shows the velocity dispersion of the dwarf galaxy in the no-feedback case: It barely reaches 10% of the tangential velocity. As a consequence, the corresponding magnetic field remains mostly toroidal, as shown in Figure 3. One can also clearly see in this Figure that the initial field parity (odd for the dipole and even for the quadrupole) has been conserved during the collapse and the subsequent shear amplification, providing a direct dependence of the final field parity on the initial halo field topology.
Supernovae feedback case
We now describe our results for the dwarf galaxy with supernovae feedback enabled. The evolution is drastically different, with violent outflows terminating quickly the life of the dense star-forming clouds. The resulting star formation rate is reduced by one order of magnitude, compared to the no-feedback case. As shown in Teyssier et al. (2013), the galactic circular velocity is now in much better agreement with observed dwarf galaxies, exhibiting a kpc-sized core in the dark matter distribution. Star formation also proceeds in successive starbursts, leading to the ejection of massive quantities of gas into a galactic fountain. The gas falls back after a dynamical time, triggering a new star formation event.
The corresponding magnetic energy evolution can be seen in in Figure 2. We observe, for both dipole and quadrupole initial conditions, a very fast, exponential growth with e-folding time of around 200 Myr. The measured growth rate is therefore quite fast, comparable to the rotation rate . Note that the magnetic energy has been amplified by almost 18 orders of magnitude, which correspond to 9 orders of magnitude in the magnetic field itself. Figure 7 compares the rotational velocity to the vertical velocity dispersion . We have in the feedback case , a clear sign of strong turbulence. As a result, field lines are violently twisted in random directions, allowing the radial and vertical components of the field to reach a similar strength than the toroidal component. Figure 3 illustrates this isotropy in the magnetic field. The field is now highly turbulent, with the largest fluctuation seen on the smallest scales (a few cell in size).
To study further the interplay between the supernovae-driven turbulence and the growth of the field, we have computed the power spectra of both the gas kinetic energy and the gas magnetic energy in Figure 8. For the former, we see the almost immediate onset of a power law power spectrum like the one predicted by Kolmogorov (), quite typical of supersonic turbulence. The kinetic energy power spectrum appears very stable throughout the evolution, maintained at this high level by supernovae explosions and rotational energy. Note that the injection scale of the kinetic energy is very large here: it is the size of the entire galaxy. To get an idea of the nature of the forcing of the turbulence in this feedback-dominated galaxy, we have plotted in Figure 5 a rendering of the gas velocity field. It shows large-scale upward and downward motions, together with a clear overall rotation pattern. The largest scale at which kinetic energy is injected turns out the be roughly the halo scale radius , for which marks the transition between the two power law regimes in the dark matter distribution (from deep inside to in the outskirts). In what follows, this radius will be identified to the kinetic energy injection scale, also noted . Note that in our spectral analysis, the global rotation was not removed from the velocity field before computing its Fourier transform, because turbulence is clearly dominating. As an additional caveat, we also note that both the gas density and the magnetic field are far from being homogeneously distributed in the box where the spectrum is computed, so that the signals are only approximately isotropic and far from periodic.
The magnetic energy power spectrum, on the other hand, is plotted in Figure 8 (bottom). We see here again that its amplitude is exponentially growing, while its shape remains roughly the same, with on the large-scale end. The power spectrum reaches a maximum at a scale corresponding to 5 cells, then slowly declines as its gets to the Nyquist frequency of the grid. This is exactly what is predicted from Kazantsev’s theory of turbulent dynamos (Kazantsev, 1968), and confirmed in the forced-turbulence, periodic box MHD simulations of Haugen et al. (2004). In the present dwarf galaxy simulations, we also obtain a small-scale magnetic dynamo, for which the forcing scale would be the size of the entire galaxy kpc, and for which the magnetic dissipation scale would be set by the adopted numerical resolution.
In the small-scale dynamo theory, a critical ingredient is the magnetic Reynolds number which encode the magnitude of the small-scale magnetic dissipation. In our notation, the velocity dispersion at the forcing scale is and is the magnetic dissipation coefficient. As discussed in Brandenburg & Subramanian (2005), exponential growth of the field is obtained if , where the critical magnetic Reynolds number was observed to be around 30-35. In our case, where no explicit magnetic dissipation has been included in the induction equation, this translates into a critical spatial resolution, beyond which we expect to see exponential amplification of the field.
For this reason, we repeat the dwarf simulations at different resolutions to study their effect. We show in Figure 9 the time evolution of the total magnetic energy in our dwarf galaxy with 3 different maximum resolutions: 36, 18 and 9 pc. We see a tendency for stronger amplification from compression with increasing resolution in all cases because the gas can build structures with higher density. Yet, without supernova feedback the subsequent shearing amplification does not show a significant impact of resolution in the resulting magnetic energy at the end. The feedback-driven dynamo, on the other hand, shows a strong growth rate dependence on the maximum resolution. In the low-resolution run, we obtain a rather slow amplification, with . This means we are close to or slightly better than the critical resolution for small-scale dynamo. For our fiducial resolution, we see a fast exponential growth with . In the high resolution case, we observe an even faster growth with . In Figure 10, images of gas density projections are plotted to show how the gas structure is resolved at different resolutions. With increasing resolution, we see more dense substructures and importantly finer details in the flow and stronger winds.
This behaviour can be explained nicely within the framework of small-scale dynamos, for which the growth rate is determine by the inverse of the eddy turn-over time at the dissipation scale , namely for Kolmogorov’s turbulence. At higher resolution, the eddy turn-over time scale become shorter, so that the growth rate becomes larger. In the kinematic phase, well before the saturation phase, when the magnetic energy will reach equipartition with the kinetic energy, the growth rate of the small-scale dynamo is therefore determined by the smallest resolved scale of the turbulence, , and propagates through an inverse cascade all the way up to the forcing scale , following the Kazantsev’s power law.
Since actual microscopic dissipation processes are occurring on very tiny scales, completely unresolved even by our highest resolution run, one expects the actual growth rate to be even higher than the already fast rates we have measured in our numerical experiments. It is therefore reasonable to assume that field saturation will be reached very quickly. We however defer the study of the saturation phase to a follow-up, companion paper (hereafter Paper II).
4.2 Milky-Way-like galaxy
We now describe our results for the Milky-Way-sized halo. In this case, the gas quickly cools down to K, which leads to the immediate formation of a very thin disk. Our spatial resolution in the Milky-Way case is limited to 84 pc, so we can’t allow the gas to cool much below 3000 K, as summarised in Table 1. Nevertheless, the disk is so massive that it also quickly fragments into massive clumps that actively form stars.
No feedback and supernovae feedback cases
In the no-feedback case, the galactic disk develops a massive central condensation of stars, which results in a central circular velocity close to 500 km/s. With supernovae feedback, however, we manage to reduce significantly the central bulge mass, with a maximum rotational velocity approaching only 250 km/s (see Figure 15). Note that even with supernovae feedback included, we do not reduce the overall star formation efficiency, and after several Gyr, most of the baryons have been converted into stars. This can be explained by the relatively low specific energy released by supernovae, making it very difficult for the gas to reach the escape velocity of the halo (see for example Roškar et al., 2014, for a complete discussion). A small galactic fountain sets in, so that the gas remains mostly bound to a weakly turbulent disc, as can be seen on Figure 15 in the vertical velocity dispersion profile.
Figure 11 shows the magnetic energy evolution for these two cases, and they exhibit the same features as the no-feedback dwarf galaxy case: early magnetic field amplification due to gravitational collapse of the cooling halo gas, followed by a weak shear amplification, with some fluctuations associated with fragmenting-clumps-induced vortex modes. The corresponding magnetic field topology appears as very well organised on large scales, with a dominating toroidal component.
Radiative Feedback case
Following the methodology explored for the first time in Roškar et al. (2014), we now consider a feedback model based on radiation from young stars efficiently absorbed by dust, and converted into kinetic energy through the infrared radiation force. Bear in mind that this model has been designed to be very optimistic, in order to maximise the effect of the radiation pressure on dust grains. Realistic modeling of the physical underlying processes is not our main objective here and we rather want to obtain an efficient feedback mechanism, launching a strong enough galactic wind, and analyse the possible effect of the resulting fountain flow on a magnetic dynamo. Images of gas density projections are shown in Figure 12. We have modelled our Milky-Way-like galaxy using two different values for the dust opacity: and , which span a range of realistic dust temperatures. The latter model gives rise to the strongest galactic fountain, and resembles in many aspects to the dwarf galaxy case with supernovae feedback. The former, lower opacity case appears as less energetic, with a weaker wind and slightly smaller turbulence. These differences can be seen in the gas images in Figure 13, and are expressed quantitatively using the vertical velocity dispersion and the rotational velocity profiles (see Figure 15), the high opacity, more extreme case giving rise to a quasi-spheroidal galaxy, with .
When looking at the corresponding magnetic energy evolution on Figure 11, it is interesting to notice that the higher velocity dispersion corresponds to the faster growth rate. Although both galaxies show an exponential growth of the magnetic energy, only the high opacity simulation reaches a growth rate as high as . The low opacity case only reaches . This is again in line with the theory of small-scale magnetic dynamos, for which the growth rate , is proportional to the amplitude of the velocity fluctuations. Here again, our turbulent forcing scale is the size of the entire galaxy, and the dynamo growth rate is set by the the dissipation scale, which in our case corresponds to the numerical spatial resolution.
We have performed MHD simulations of feedback-dominated galaxies, both dwarf galaxies and Milky-Way-sized galaxies. We have shown that, if feedback processes are strong enough, a small-scale dynamo sets in, with a typical kinetic energy injection scale corresponding to the size of the entire galaxy (which is in our case is close to the halo scale radius ), and a typical magnetic dissipation scale corresponding in our case to the adopted spatial resolution. We have observed an exponential increase of the magnetic energy, with growth rate , higher than the galaxy rotation rate, possibly much higher if one considers realistic microscopic diffusion processes instead of only numerical diffusion.
Three important aspects are missing in order to apply our findings to the origin of galactic magnetic fields: 1- We have considered rather idealised simulations of galaxies in isolation. 2- We have only described the kinematic phase, deferring the discussion of the saturation to a companion paper (Paper II). 3- We have considered feedback-dominated galaxies, which are relevant for the high-redshift universe. What will happen after this feedback-dominated epoch, for quiescent, razor-thin galactic discs ? In this section, we speculate on possible cosmological consequences of our results on the nature of the magnetic field in high redshift galaxies, as well as the magnetic strength and topology in lower redshift galaxies.
5.1 Cosmological implications for high redshift galaxies
Although our numerical simulations were not performed in a realistic cosmological context, we can still draw conclusions for the cosmic evolution of magnetic fields, assuming that the universe is made of a collection of halos of various masses, and generalised our results using simple analytical estimates. For this purpose, we will assume that high-redshift galaxies are all dominated by efficient feedback processes, so that galactic winds can drive a powerful fountain, resulting in a gas rich, turbulence-dominated corona, with a size equal to , the halo scale radius, that sets the turbulence injection scale. As shown in the previous sections, a very efficient small-scale dynamo is likely to develop, with a growth rate larger (possibly much larger) than the rotation rate of the galaxy. Although a detailed study of the saturation phase is required to study how fast the field will reach equipartition (and at what scales, see Paper II), we postulate here that each halo reaches equipartition between magnetic and turbulent kinetic energy almost instantaneously, within a volume set by the halo scale radius . This leads to the equipartition value for the field:
For the gas density, we assume that its average value can be approximate by the baryon fraction of the total mass halo scale density
and for the turbulence velocity dispersion, we adopt the value inferred by our simulations, namely the halo maximum circular velocity (at )
Using standard redshift-dependant functions for these two halo parameters and , assuming an average halo concentration parameter and a baryon fraction , we obtain a prediction for the saturated field
The resulting magnetic pressure scales as a function of mass and redshift exactly like a “viral pressure” in the halo. Note that this equipartition field increases quite fast with increasing redshift, but quite slowly with increasing halo mass. For a feedback-dominated, Milky-Way-sized galaxy, at redshift , one still predicts a rather strong magnetic field, slightly above 100. Although very interesting in setting the foundations for a theory of small-scale dynamo in the cosmological context, our present discussion remains speculative, in the sense that we do not properly model cosmological infall and the associated hierarchical merging of smaller structures. This could lead to a dilution of the dynamo-amplified field and lower the growth rate. It is therefore of primary importance to simulate such feedback-dominated galaxy formation models with both MHD and a realistic cosmological environment.
5.2 Transition to quiescent, low redshift galaxies
The cosmological applications we have derived from a simple analytical extension of our numerical results was within the context of feedback-dominated, high-redshift galaxies. Low redshift galaxies are however quite different. We see in our nearby universe many grand design, quiescent disk galaxies, with very thin, low velocity dispersion disks, and for larger halo masses, we even see red and dead elliptical galaxies. The present day universe is therefore very quiet, and strong feedback effects are absent, except may be in some very intense starbursts, usually triggered by (very rare) merger events. Our present methodology, based on isolated, gas rich, cooling halos, does not allow to explore the low redshift universe self-consistently, unless one artificially switches off feedback processes. This is the strategy we adopt in this section, in order to explore the consequence of evolving our simulated objects from a feedback-dominated state, to a more quiescent state, without strong galactic fountains, resulting in much thinner, rotationally supported disks.
It is indeed very important to estimate how the magnetic field, amplified first through a turbulent-driven small-scale dynamo, could evolve into a large-scale field, mostly driven by rotational shear. The obvious question one might ask is: Does the magnetic energy disappear, as the small-scale turbulent field reconnects on small scales? How intense will the surviving large-scale and mostly toroidal magnetic field be, after the galaxy develops into a thin, centrifugally supposed disk? For this reason, we decide to re-start our dwarf galaxy simulations after 3 Gyr of small-scale dynamo amplification, but without any stellar feedback. In the absence of galactic winds, the turbulent corona rapidly collapses back to a thin disk, and a mostly toroidal field appears after a few rotations. Images for gas density, velocity field and magnetic field of the galaxy after it has collapsed are shown in Figure 16.
We see in Figure 17 the evolution of the magnetic energy versus time for this “suppressed feedback” simulation. We see that overall, after the thin disk appears, the magnetic energy is for the most part conserved. Two competing effects are indeed at work here: gravitational collapse that amplifies the field from a corona-diluted state to a thin disk-concentrated state, on one hand, and magnetic losses due to reconnection of mostly small-scale field lines in the mid plane of the disk. What is interesting and highly non-trivial, is that these two effects basically cancel each other, and that the final magnetic energy in the thin disk is the same than the initial magnetic energy in the large, turbulence-supported corona.
Large-scale, collapse-amplified magnetic fields have replaced small-scale, reconnection-suppressed fluctuations. Large-scale modes arise in the small-scale dynamo picture because the field is amplified on all scales at the same rate, up to the turbulence forcing scale, a well-known property of small-scale dynamos. Enough magnetic energy has been stored on large scales, so it can survive and compensate for reconnection and field cancellation effects during the collapse.
We would like to stress here again that this new type of simulations with suppressed feedback was performed in the weak field, pure kinematic regime as well. Since we believe that after the feedback-driven, small-scale dynamo phase, the field has probably quickly reached equipartition, we need to study this transition from thick-corona to thin-disk using a properly saturated field, which is the purpose of paper II.
Nevertheless, we can also study the topology of the field after the disk has settled in a thin, rotation-supported state. Figure 19 shows the toroidal, radial and vertical field in the disk 2 Gyr after the feedback has been suppressed. The field is mostly tangential, with however stronger radial and vertical components, compared to the no-feedback case. This is in agreement with observations which show pitch angles as high as (Patrikeev et al., 2006). The most interesting results is the topology of the field, which appears quadrupole-like, even if we start with only a dipole in the initial conditions. Dipole-like modes present in the large-scale magnetic field of the corona have odd parity. They will cancel in the mid plane after the gas has collapsed into a thin disk (like in the no-feedback case). Quadrupole-like modes, on the other hand, have even parity and they will be combined non-destructively in the mid plane after the collapse of the turbulent corona. One interesting prediction of the scenario which we consider here is therefore a systematic quadrupole field topology, independently of the initial topology of the primordial field. This result is in very good agreement with observational data of the Milky Way (Taylor et al., 2009; Oppermann et al., 2012), as well as external galaxies (Braun et al., 2010; Mao et al., 2012).
We have performed MHD simulations of cooling halos, for both dwarf and Milky-Way-sized haloes, in the kinematic regime where magnetic field are weak enough so that the effect of the Lorentz force on the turbulent flow is insignificant. Using supernovae feedback for dwarf galaxies, as well as radiation feedback for large galaxies, we have shown that a small-scale dynamo quickly sets in, with the turbulent energy injection scale roughly equal to the halo scale radius , and dissipation scale roughly equal to 4-5 cell sizes. In agreement with small-scale dynamo theory, we observe an exponential amplification of the magnetic energy on all scales, up to the injection scale, with a growth rate at least equal to twice the disk rotation rate. The growth rate in the kinematic phase is set by the adopted numerical resolution. In our simulated halo, we need to have a resolution of 20 pc or better (for kpc) to obtain a significant growth rate (as high as ).
This corresponds roughly to 100 resolution elements per turbulent energy injection scale. We believe this resolution corresponds to the critical magnetic Reynolds number , beyond which magnetic field amplification through small-scale dynamo is possible. Note that for weaker feedback scenarios, the kinetic energy injection scale is likely to be smaller, which translates into more stringent resolution requirements. This analysis should be kept in mind, when one wants to extend this to the cosmological context. One difficulty might arise when one considers infall of pristine, low magnetic field gas as a source of field dilution. This could require a larger small-scale dynamo growth rate for the field to be able to increase, and therefore results in even more demanding resolution requirements.
We have shown that a small-scale dynamo, driven by strong feedback processes in high-redshift galaxies, could be the origin of galactic magnetic fields. This scenario is completely different than the more traditional large-scale dynamo approach, in the sense that the small-scale dynamo acts very quickly at amplifying the field up to the turbulent injection scale. The new ingredient is here the fact that feedback processes at high redshift are probably energetic enough, so that this injection scale is the size of the entire galaxy (more precisely the halo scale radius ). We therefore have a small-scale dynamo, together with a large-scale forcing, hence enabling the fast amplification of the field all the way up to the scale of the entire galaxy.
We have also shown, using simple numerical experiments with suppressed feedback, that the field can evolve into a large-scale toroidal, quadrupole field in a low-redshift quiescent state, although it has been amplified by a small-scale dynamo. Compared to the large-scale dynamo picture, for which the field is amplified in a thin disk over several Gyr, the small-scale dynamo picture completely reverses the point of view, with a very strong field already in place at high redshift (after a feedback-dominated epoch) that slowly evolves until the present epoch. This large-scale field could even slowly decay, on resistive time scales, and still be large enough to account for the observed field strength in nearby galaxies. In this new picture, large-scale dynamos are not required anymore to amplify the field from its primordial value around to levels. Instead, they are needed to maintain the field at a roughly constant level by compensating dissipative losses.
We thank the anonymous referee for their thoughtful comments and helpful suggestions. This work was funded by the Swiss National Science Foundation SNF. All Computations were performed on the local zBox super computer and on the Piz Dora cluster at the Swiss National Supercomputing Center CSCS.
- Agertz O., Lake G., Teyssier R., Moore B., Mayer L., Romeo A. B., 2009, Monthly Notices of the Royal Astronomical Society, 392, 294
- Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, The Astrophysical Journal, 770, 25
- Balsara D. S., 2001, Journal of Computational Physics, 174, 614
- Beck R., Brandenburg A., Moss D., Shukurov A., Sokoloff D., 1996, Annual Review of Astronomy & Astrophysics, 34, 155
- Beck A. M., Lesch H., Dolag K., Kotarba H., Geng A., Stasyszyn F. A., 2012, Monthly Notices of the Royal Astronomical Society, 422, 2152
- Behroozi P. S., Wechsler R. H., Conroy C., 2013, The Astrophysical Journal, 770, 57
- Bernet M. L., Miniati F., Lilly S. J., Kronberg P. P., Dessauges-Zavadsky M., 2008, Nature, 454, 302
- Biermann L., 1950, Zeitschrift Naturforschung Teil A, 5, 65
- Brandenburg A., Subramanian K., 2005, Physics Reports, 417, 1
- Braun R., Heald G., Beck R., 2010, Astronomy & Astrophysics, 514, A42
- Charbonneau P., Steiner O., 2012, Solar and Stellar Dynamos: Saas-Fee Advanced Course 39 Swiss Society for Astrophysics and Astronomy
- Childress S., Gilbert A. D., 1995, The Fast Dynamo, -1
- Dekel A., et al., 2009, Nature, 457, 451
- Del Sordo F., Guerrero G., Brandenburg A., 2013, Monthly Notices of the Royal Astronomical Society, 429, 1686
- Draine B. T., Li A., 2007, The Astrophysical Journal, 657, 810
- Dubois Y., Teyssier R., 2010, Astronomy and Astrophysics, 523, 72
- Ferriere K., 1992, Astrophysical Journal, 391, 188
- Fromang S., Hennebelle P., Teyssier R., 2006, \aap, 457, 371
- Galloway D., Frisch U., 1986, Geophysical and Astrophysical Fluid Dynamics (ISSN 0309-1929), 36, 53
- Governato F., et al., 2010, Nature, 463, 203
- Governato F., et al., 2012, Monthly Notices of the Royal Astronomical Society, 422, 1231
- Hanasz M., Kowal G., Otmianowska-Mazur K., Lesch H., 2004, The Astrophysical Journal, 605, L33
- Haugen N. E. L., Brandenburg A., Dobler W., 2004, Physical Review E, 70, 016308
- Hopkins P. F., Hopkins P. F., Quataert E., Murray N., 2012, Monthly Notices of the Royal Astronomical Society, 421, 2654
- Hopkins P. F., Kereš D., Onorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, Monthly Notices of the Royal Astronomical Society, 445, 581
- Kauffmann G., 2014, Monthly Notices of the Royal Astronomical Society, 441, 2717
- Kazantsev A. P., 1968, Soviet Physics JETP, 26, 1031
- Kazantzidis S., Magorrian J., Moore B., 2004, The Astrophysical Journal, 601, 37
- Kereš D., Katz N., Weinberg D. H., Davé R., 2005, Monthly Notices of the Royal Astronomical Society, 363, 2
- Kimm T., Cen R., 2014, The Astrophysical Journal, 788, 121
- Krumholz M. R., Tan J. C., 2007, The Astrophysical Journal, 654, 304
- Kulsrud R. M., Anderson S. W., 1992, Astrophysical Journal, 396, 606
- Mao S., et al., 2012, The Astrophysical Journal, 759, 25
- Martin C. L., 1999, The Astrophysical Journal, 513, 156
- Moster B. P., Naab T., White S. D. M., 2013, Monthly Notices of the Royal Astronomical Society, 428, 3121
- Murray N., Quataert E., Thompson T. A., 2010, The Astrophysical Journal, 709, 191
- Naoz S., Narayan R., 2013, Physical Review Letters, 111, 051303
- Navarro J. F., Frenk C. S., White S. D. M., 1997, Astrophysical Journal v.490, 490, 493
- Ocvirk P., Pichon C., Teyssier R., 2008, Monthly Notices of the Royal Astronomical Society, 390, 1326
- Oppenheimer B. D., Davé R., 2006, Monthly Notices of the Royal Astronomical Society, 373, 1265
- Oppermann N., et al., 2012, Astronomy & Astrophysics, 542, A93
- Pakmor R., Springel V., 2013, Monthly Notices of the Royal Astronomical Society, 432, 176
- Pakmor R., Bauer A., Springel V., 2011, Monthly Notices of the Royal Astronomical Society, 418, 1392
- Pariev V. I., Colgate S. A., Finn J., 2007, The Astrophysical Journal, 658, 129
- Parker E. N., 1970, Astrophysical Journal, 160, 383
- Parker E. N., 1992, Astrophysical Journal, 401, 137
- Patrikeev I., Fletcher A., Stepanov R., Beck R., Berkhuijsen E., Frick P., Horellou C., 2006, Astronomy & Astrophysics, 458, 441
- Ponomarenko Y. B., 1973, Journal of Applied Mechanics and Technical Physics, 14, 775
- Pontzen A., Governato F., 2012, Monthly Notices of the Royal Astronomical Society, 421, 3464
- Rasera Y., Teyssier R., 2006, Astronomy and Astrophysics, 445, 1
- Read J. I., Wilkinson M. I., Evans N. W., Evans N. W., Gilmore G., Kleyna J. T., Kleyna J. T., 2006, Monthly Notices of the Royal Astronomical Society, 367, 387
- Rosen A., Rosen A., Bregman J. N., 1995, Astrophysical Journal v.440, 440, 634
- Roškar R., Teyssier R., Agertz O., Wetzstein M., Moore B., 2014, Monthly Notices of the Royal Astronomical Society, 444, 2837
- Scannapieco C., et al., 2012, Monthly Notices of the Royal Astronomical Society, p. 2970
- Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, Astronomy and Astrophysics, 410, 611
- Steidel C. C., Erb D. K., Shapley A. E., Pettini M., Reddy N., Bogosavljević M., Rudie G. C., Rakic O., 2010, The Astrophysical Journal, 717, 289
- Sutherland R. S., Dopita M. A., 1993, Astrophysical Journal Supplement Series (ISSN 0067-0049), 88, 253
- Taylor A., Stil J., Sunstrum C., 2009, The Astrophysical Journal, 702, 1230
- Teyssier R., 2002, \aap, 385, 337
- Teyssier R., Fromang S., Dormy E., 2006, Journal of Computational Physics, 218, 44
- Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, Monthly Notices of the Royal Astronomical Society, 429, 3068
- Tóth G., Roe P. L., 2002, Journal of Computational Physics, 180, 736
- Truelove J. K., Klein R. I., McKee C. F., Holliman J. H. I., Howell L. H., Greenough J. A., 1997, Astrophysical Journal Letters v.489, 489, L179
- Vainshtein S. I., Cattaneo F., 1992, Astrophysical Journal, 393, 165
- Wang P., Abel T., 2009, The Astrophysical Journal, 696, 96
- Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, preprint, (arXiv:1503.04818)
- Widrow L. M., 2002, Reviews of Modern Physics, 74, 775
- de Blok W. J. G., McGaugh S. S., Bosma A., Rubin V. C., 2001, The Astrophysical Journal, 552, L23