Magnetic field amplification in turbulent astrophysical plasmas

Magnetic field amplification in turbulent astrophysical plasmas

Christoph Federrath1 Email address for correspondence: christoph.federrath@anu.edu.au
Abstract

Magnetic fields play an important role in astrophysical accretion discs, and in the interstellar and intergalactic medium. They drive jets, suppress fragmentation in star-forming clouds and can have a significant impact on the accretion rate of stars. However, the exact amplification mechanisms of cosmic magnetic fields remain relatively poorly understood. Here I start by reviewing recent advances in the numerical and theoretical modelling of the turbulent dynamo, which may explain the origin of galactic and inter-galactic magnetic fields. While dynamo action was previously investigated in great detail for incompressible plasmas, I here place particular emphasis on highly compressible astrophysical plasmas, which are characterised by strong density fluctuations and shocks, such as the interstellar medium. I find that dynamo action works not only in subsonic plasmas, but also in highly supersonic, compressible plasmas, as well as for low and high magnetic Prandtl numbers. I further present new numerical simulations from which I determine the growth of the turbulent (un-ordered) magnetic field component () in the presence of weak and strong guide fields (). I vary over 5 orders of magnitude and find that the dependence of on is relatively weak, and can be explained with a simple theoretical model in which the turbulence provides the energy to amplify . Finally, I discuss some important implications of magnetic fields for the structure of accretion discs, the launching of jets, and the star formation rate of interstellar clouds.

1 Introduction

Magnetic fields are ubiquitous in the Universe. Examples include astrophysical accretion discs around young stars and active galactic nuclei, the interstellar and intergalactic medium, and even the early Universe when the first stars and galaxies formed.

The turbulent dynamo mechanism is believed to be the main cause of cosmic magnetism (Brandenburg & Subramanian 2005). Magnetic fields are amplified exponentially via the turbulent dynamo, leading to dynamically important magnetic forces and energies on relatively short time scales. The amplification of the magnetic field arises from sequences of “stretching, twisting, and folding” of the field, until the magnetic field lines are so tightly packed that the magnetic energy density becomes comparable to the kinetic energy density provided by turbulent motions.

Dynamo action ranges from the Earth and the Sun (Cattaneo & Hughes 2001), over the interstellar medium to whole galaxies (Beck et al. 1996; Beck 2016). The turbulent dynamo is important for the formation of the large-scale structure of the Universe (Ryu et al. 2008; Miniati & Bell 2011; Iapichino & Brüggen 2012; Vazza et al. 2014; Miniati & Beresnyak 2015; Beresnyak & Miniati 2016), in clusters of galaxies (Subramanian et al. 2006) and in the formation of the first cosmological objects in dark matter haloes (Schleicher et al. 2010). It determines the growth of magnetic energy in solar convection (Cattaneo & Hughes 2001; Moll et al. 2011; Pietarila Graham et al. 2010), in the interior of planets (Roberts & Glatzmaier 2000) and is relevant for liquid metal experiments on Earth (Monchaux et al. 2007). It may further explain the far-infrared–radio correlation in spiral galaxies (Schleicher & Beck 2013). After the turbulent dynamo has amplified tiny seeds of the magnetic field, which can be generated during inflation, the electroweak or the QCD phase transition (Grasso & Rubinstein 2001), the large-scale dynamo kicks in and generates the large-scale magnetic fields observed in planets, stars and galaxies today (Beck et al. 1996; Brandenburg & Subramanian 2005).

One of the most important distinctions that we have to make when considering turbulent gases, fluids and plasmas, is the level of compressibility of the medium. For instance, in the Earth and the Sun, the dynamo is driven by subsonic, nearly incompressible flows. By contrast, interstellar clouds and galaxies are dominated by highly compressible turbulence (Mac Low & Klessen 2004; Elmegreen & Scalo 2004; McKee & Ostriker 2007; Hennebelle & Falgarone 2012; Padoan et al. 2014). Indeed, the gas densities in the interstellar medium range from particle per (the average gas density in a Milky-Way type galaxy; see Ferrière 2001) to (where the gas becomes optically thick and proceeds to the formation of a star).

Numerical studies of non-driven turbulence demonstrated that supersonic turbulence decays quickly, in about a crossing time (Scalo & Pumphrey 1982; Mac Low et al. 1998; Stone et al. 1998; Mac Low 1999). Since we observe highly compressible interstellar turbulence, that turbulence must be driven by some physical forcing mechanisms. Those mechanisms include supernova explosions, ionizing shells from high-mass stellar feedback (McKee 1989; Krumholz et al. 2006; Balsara et al. 2004; de Avillez & Breitschwerdt 2005; Breitschwerdt et al. 2009; Gritschneder et al. 2009; Peters et al. 2010, 2011; Goldbaum et al. 2011; Lee et al. 2012), gravitational collapse and accretion of material (Hoyle 1953; Vazquez-Semadeni et al. 1998; Klessen & Hennebelle 2010; Elmegreen & Burkert 2010; Vázquez-Semadeni et al. 2010; Federrath et al. 2011b; Robertson & Goldreich 2012), and galactic spiral-arm compression and cloud-cloud collisions (Dobbs & Bonnell 2008; Dobbs et al. 2008; Tasker & Tan 2009; Benincasa et al. 2013), as well as magneto-rotational instability (Piontek & Ostriker 2004, 2007; Tamburro et al. 2009). Jets and outflows from young stellar objects (Norman & Silk 1980; Banerjee et al. 2007; Nakamura & Li 2008; Cunningham et al. 2009; Carroll et al. 2010; Wang et al. 2010) and active galactic nuclei (Mukherjee et al. 2016) also drive turbulence. Turbulence in high-redshift galaxies and their low-redshift analogues may be driven by feedback (Green et al. 2010; Fisher et al. 2014). A discussion and comparison of turbulence driving mechanisms is published in Mac Low & Klessen (2004), Elmegreen (2009), and Federrath et al. (2016a).

Most importantly, the majority of turbulence drivers (e.g., supernova explosions, high-mass stellar winds, and accretion) are expected to drive compressible modes, so we refer to these as “compressive drivers”. On the other hand, solenoidal modes can be generated directly by shear (Federrath et al. 2016b) and the MRI (so we call them “solenoidal drivers”), and indirectly by nonlinear interactions of multiple colliding shock fronts (Vishniac 1994; Sun & Takayama 2003; Kritsuk et al. 2007; Federrath et al. 2010), baroclinity (Padoan et al. 2016), rotation and shear (Del Sordo & Brandenburg 2011), as well as viscosity (Mee & Brandenburg 2006; Federrath et al. 2011a). Thus, turbulence driven by purely compressive drivers can still contain up to half of its kinetic power in solenoidal modes (Federrath et al. 2010, Figure 14).

In order to understand the dependence of the dynamo on the compressibility of the plasma, we present in §2 a systematic study in which the turbulent driving and Mach number are varied, covering the subsonic, nearly incompressible regime as well as the high-compressible, supersonic regime of turbulence. We determine the growth rate and saturation level of the magnetic field in both regimes. In §3 we determine the dependence of the compressible turbulent dynamo on the magnetic Prandtl number and the kinematic Reynolds number.

Many astrophysical systems are characterised by weak or strong magnetic guide fields, i.e., fields that have an ordered mean component along a specified direction. This is in contrast to the turbulent dynamo, where primarily the turbulent un-ordered magnetic field component is amplified on small scales. In order to measure and understand magnetic field amplification in the presence of a strong guide field, we present new simulations in §4, where we systematically measure the strength of the turbulent magnetic field component as a function of increasing magnetic guide field.

Finally, in §5 we discuss the implications of magnetic fields for astrophysical accretion discs, the structure of the interstellar medium, and for star formation.

2 The Mach number dependence of the turbulent dynamo

Here we investigate fundamental properties of the turbulent dynamo—its growth rate and saturation level—in simulations with extremely different driving of the turbulence (solenoidal versus compressive) and with a range of Mach numbers from to . Some of the results in this section are published in Federrath et al. (2011a). While these simulations are highly idealised, they serve as a systematic numerical experiment in which the driving and Mach number are controlled. Turbulent astrophysical systems can often be characterised by a few fundamental numbers that cover the basic physical behaviour. For example, the turbulence in astrophysical accretion discs around young stars is subsonic () and is primarily driven by shear or magneto-rotational instability, which drive solenoidal motions (solenoidal driving). By contrast, molecular cloud turbulence is highly compressible and supersonic (), and is driven by a range of physical processes, most of which induce compression (compressive driving), such as supernova explosions, expanding radiation fronts from high-mass stellar feedback and/or galactic spiral shocks. For a review of potential drivers of the turbulence in galaxies and the interstellar medium, please see the articles by Elmegreen (2009), Federrath et al. (2010), Federrath & Klessen (2012), and Padoan et al. (2014).

2.1 Background and open questions

Most studies of the turbulent dynamo concentrate on incompressible plasmas, with only few studies approaching the effects of compressibility. For example, Haugen et al. (2004b) obtained critical Reynolds numbers for dynamo action, but did not investigate growth rates or saturation levels, and only studied Mach numbers in the range . The energy released by e.g. supernova explosions, however, drives interstellar and galactic turbulence with Mach numbers of 10–100 (Mac Low & Klessen 2004). Thus, much higher Mach numbers have to be investigated in order to understand dynamo action in interstellar clouds. It is furthermore tempting to associate such supernova blast waves with compressive driving of turbulence (Mee & Brandenburg 2006; Schmidt et al. 2008; Federrath et al. 2010). Mee & Brandenburg (2006) concluded that it is very hard to excite the turbulent dynamo with this curl-free driving, because it does not directly inject vorticity, .

Here we show that the turbulent dynamo can be driven by curl-free driving mechanisms, and we quantify the amplification as a function of compressibility of the plasma. The main questions addressed are: How does the turbulent dynamo depend on the Mach number and driving of the turbulence? What is the field geometry and amplification mechanism? What are the growth rates and saturation levels in the subsonic and supersonic regime?

2.2 Methods

2.2.1 Magnetohydrodynamical equations

We use a modified version of the FLASH code (Fryxell et al. 2000; Dubey et al. 2008) to solve the three-dimensional (3D), compressible, magnetohydrodynamical (MHD) equations on uniform, periodic grids, including viscous and resistive dissipation terms,

(2.0)
(2.0)
(2.0)
(2.0)
(2.0)

Here, , , , , and denote the gas density, velocity, pressure (thermal plus magnetic), magnetic field, and energy density (internal plus kinetic, plus magnetic), respectively. Physical shear viscosity is included with the traceless rate of strain tensor, in the momentum equation (2.2.1), and controlled by the kinematic viscosity (). Physical diffusion of is controlled by the magnetic diffusivity (the inverse of the electrical conductivity ) in the induction equation (2.2.1). The MHD equations are closed with an isothermal equation of state, , where is the sound speed.

We note that equations (2.2.1) and (2.2.1) only contain the shear viscosity , while the bulk viscosity is assumed to be zero. For monatomic ideal gases, the bulk viscosity is indeed identically zero, which can be derived from kinetic theory (Mihalas & Mihalas 1984). For polyatomic molecules, this does not need to be the case, but only, if a relaxation process takes place that is of the same order or slower than a typical fluid timescale (Mihalas & Mihalas 1984). The value of strongly depends on the composition of the polyatomic gas and measurements of from different experiments give different results (Tisza 1942). For simplicity, we set the bulk viscosity and only consider the well-understood shear viscosity . We further note that using an equation of state for the gas that relates pressure with density and temperature (which is the standard approach in fluid dynamics) implies (Truesdell 1952).

2.2.2 Turbulence driving

In order to drive turbulence with a target Mach number, we apply the driving field as a source term in the momentum equation (2.2.1). The driving field is constructed with a stochastic Ornstein-Uhlenbeck (OU) process (Eswaran & Pope 1988; Schmidt et al. 2009; Price & Federrath 2010), implemented by Federrath et al. (2010) in the public version of the FLASH code (http://www.flash.uchicago.edu/site/flashcode/). The OU process yields a driving pattern that varies smoothly in space and time with an auto-correlation time equal to the eddy-turnover time (also called turbulent box-crossing time), on the largest scales () in the periodic simulation domain of side length . The turbulent sonic Mach number is defined as , which is the ratio of the turbulent velocity dispersion on scale and the sound speed . The driving is constructed in Fourier space such that the kinetic energy is injected at the smallest wave numbers, . The peak of energy injection is on scale , i.e., , and falls off parabolically towards smaller and higher wave numbers such that the driving power reaches zero exactly at and . This procedure limits the direct effect of the driving to a narrow wave number band () and allows the turbulence to develop self-consistently on smaller scales ().

By applying a projection in Fourier space, we can decompose the driving field into its solenoidal and compressive parts. In index notation, the projection operator reads

(2.0)

where and are the solenoidal and compressive projection operators, respectively. This yields a ratio of compressive driving power to total driving power as a function of the parameter ,

(2.0)

The projection operation allows us to construct a solenoidal (divergence-free) driving field () or a compressive (curl-free) driving field () by setting (sol) or (comp). Mixed driving ratios are obtained by picking values in the range . Setting yields a “natural mixture”, i.e., a driving field with . The latter is equivalent to the result of randomly-chosen driving modes in 3D, where on average 1 spatial dimension is occupied by longitudinal (compressive) modes and 2 spatial dimensions are occupied by transverse (solenoidal) modes (Federrath et al. 2008).

2.2.3 Initial conditions and numerical methods

We start our numerical experiments by setting , zero initial velocities , uniform density , , and with in -direction, corresponding to an extremely high initial plasma . These values are motivated by dynamo studies of primordial clouds (Schleicher et al. 2010; Sur et al. 2010; Federrath et al. 2011b; Sur et al. 2012), but in the following, we scale all quantities to dimensionless units to address fundamental questions of magnetic field amplification in compressible plasmas. The physics of these turbulent systems is fully determined by the dimensionless and plasma , as well as the driving parameter (sol vs. comp). The actual plasma densities, velocities, and magnetic fields can be scaled to astrophysical systems that are described by the same set of basic dimensionless numbers (, , ).

For most of the simulations, we set the kinematic viscosity and the magnetic diffusivity to zero, and thus solve the ideal MHD equations. In this case, the dissipation of kinetic and magnetic energy is caused by the discretisation of the MHD equations. However, we do not add artificial viscosity. We use Riemann solvers, which capture shocks well in the absence of artificial viscosity. In addition to the ideal MHD simulations, we also solve the full, non-ideal MHD equations (2.2.1)–(2.2.1) for 4 representative models to show that our results are physical and robust against changes in the numerical scheme. For the ideal MHD simulations, we use the positive-definite, split HLL3R Riemann scheme (Waagan et al. 2011) in FLASH v2.5, while our non-ideal MHD simulations were performed with the un-split staggered-mesh scheme in FLASH v4 (Lee & Deane 2009), which uses a third-order reconstruction, constrained transport to maintain to machine precision, and the HLLD Riemann solver (Miyoshi & Kusano 2005). We run simulations with , , and grid cells, in order to test numerical convergence of our results.

2.3 Results

2.3.1 Time evolution

Figure 1: Turbulent sonic Mach number () as a function of the turbulent crossing time () for all runs with solenoidal (sol) and compressive (comp) driving. These simulations cover compressible plasmas from subsonic turbulence () up into the highly compressible, supersonic regime ().

The turbulence is fully developed after an initial transient phase that lasts for (e.g., Kitsionas et al. 2009; Federrath et al. 2009; Price & Federrath 2010), and the Mach number reaches its target value, fluctuating on a 10% level. Figure 1 shows the time evolution of in all runs with solenoidal and compressive driving of the turbulence. Some simulations had to be run for a few hundred crossing times in order to reach saturation of the magnetic field. This is the case in all the compressively driven simulations, which—as we will see below—have significantly lower dynamo growth rates than in the case of solenoidal driving. Interesting to note is the drop in for the solenoidally driven runs with towards late times, . This slight drop in Mach number indicates that the magnetic field has reached very high saturation levels, such that the back reaction of the magnetic field via the Lorentz force has become significant in these models. However, the magnetic field has little dynamical impact on the turbulent flow in all supersonic runs and in all runs with compressive driving. Although the kinematics of the turbulent flow is not strongly affected in those cases, the structure and fragmentation of the gas is still influenced by the presence of a turbulent magnetic field (Hennebelle & Teyssier 2008; Federrath & Klessen 2012; Padoan et al. 2014).

Figure 2: Magnetic-to-kinetic energy ratio () as a function of the turbulent crossing time () for all runs with solenoidal (sol) and compressive (comp) driving. The time-averaged sonic Mach number () of each model is indicated in the legend (see figure 1 for the time evolution of ). The thin solid lines show exponential fits in the regime of turbulent dynamo amplification, followed by constant fits in the saturation phase. The evolution of reveals extremely different efficiencies of the dynamo, depending on the Mach number and driving of the turbulence.

Figure 2 shows the ratio of magnetic-to-kinetic energy () as a function of time. We immediately see that the magnetic energy grows exponentially over several orders of magnitude in any of the simulation models and finally reaches saturation at different levels (indicated by the fitted constant line towards late times when the curves saturate).

2.3.2 Density and magnetic field structure

Figure 3: 3D renderings of the gas density contrast on a logarithmic scale () (from white to blue), and magnetic field lines (orange) for solenoidal driving at (top left) and (bottom left), and compressive driving at (top right) and (bottom right). The stretch-twist-fold mechanism of the dynamo (Brandenburg & Subramanian 2005) is evident in all models, but operates with different efficiency due to the varying compressibility, flow structure, and formation of shocks in the supersonic plasmas. From Federrath et al. (2011a). An animation is available at http://www.mso.anu.edu.au/~chfeder/pubs/dynamo_prl/dynamo_prl.html.

Figure 3 shows 3D volume renderings of some the extreme cases (solenoidal driving on the left and compressive driving on the right, each for in the top panels and in the bottom panels). We see that the supersonic simulation runs (bottom panels) are dominated by shocks. Compressive driving yields stronger density enhancements for similar Mach numbers compared to solenoidal driving (Nordlund & Padoan 1999; Federrath et al. 2008, 2010; Price et al. 2011; Konstandin et al. 2012; Federrath & Klessen 2013; Ginsburg et al. 2013; Kainulainen et al. 2013; Federrath 2013). The magnetic field occupies large volume fractions with relatively straight field lines in the compressively driven cases, while solenoidal driving produces more space-filling, tangled field configurations. This suggests that the turbulent dynamo operates more efficiently in the case of solenoidal driving, which is quantified and explained in the following.

2.3.3 Dynamo growth rate and saturation level

Figure 4: Growth rate (top), saturation level (middle), and solenoidal ratio (bottom) as a function of Mach number, for all simulations with solenoidal (crosses) and compressive driving (diamonds). The thin solid lines show empirical fits with equation (2.3.3), the parameters of which are listed in table 1. The arrows point to 4 simulations (, for sol. and comp. driving), which used ideal MHD on grid cells (a), non-ideal MHD on (b), and grid cells (c), demonstrating convergence for the given magnetic Prandtl number, , and kinematic Reynolds number, . The theoretical predictions for the saturation level from Schober et al. (2015) are added as grey lines (middle panel) in the limit (Kolmogorov scaling exponent: ) and (Burgers scaling exponent: ).
Growth rate Saturation level Solenoidal ratio
(sol) (comp) (sol) (comp) (sol) (comp)
-18.71 -2.251 -0.020 -0.037 -0.808 -0.423
-0.051 -0.119 -2.340 -1.982 -2.850 -1.970
-1.059 -0.802 -23.33 -0.027 -1.238 0
-2.921 -25.53 -2.340 -3.601 -2.850 -1.970
-1.350 -1.686 1 -0.395 1 -0.535
-0.313 -0.139 0 -0.003 0 0
1/3 1/3 0 0 0 0
Table 1: Parameters in equation (2.3.3) for the fits in figure 4.

Figure 4 (top and middle panels) show the dynamo growth rate () in the relation and the saturation level as a function of Mach number for all simulation models. Both and depend on and on the driving of the turbulence. Solenoidal driving yields growth rates and saturation levels higher than for compressive driving. Both and change significantly at the transition from subsonic to supersonic turbulence. Within a turbulent medium, the transition from supersonic to subsonic turbulence occurs on the sonic scale (Vázquez-Semadeni et al. 2003; Federrath & Klessen 2012; Hopkins 2013; Federrath 2016). We conclude that the formation of shocks at is responsible for destroying some of the coherent vortex motions necessary to drive the dynamo (Haugen et al. 2004b). However, as is increased further, vorticity generation in oblique, colliding shocks (Sun & Takayama 2003; Kritsuk et al. 2007) starts to dominate over the destruction of vorticity.

The data points labelled with ‘a’, ‘b’, and ‘c’ in figure 4 show simulations with different numerical resolution (, and cells, respectively) and with different numerical solvers/schemes (‘a’: HLL3R with and set to zero, i.e., ideal MHD; ‘b’ and ‘c’: HLLD with non-ideal MHD; cf. §2.2). These convergence studies demonstrate that our results do not depend on the choice of numerical solver and, in particular, do not depend on the choice of dissipation mechanism (see also Benzi et al. 1993). The reason for this is that the physical dissipation was chosen to be similar to the numerical one in the low-resolution simulations, while in the high-resolution simulations, the physical dissipation dominates the numerical dissipation. This is why our simulations are converged and the results do not depend on the choice of numerical solver.

The very small growth rates in the subsonic, compressively driven models are due to the fact that hardly any vorticity is generated. In order to quantify this, we show the solenoidal ratio, i.e., the specific kinetic energy in solenoidal modes of the turbulent velocity field divided by the total specific kinetic energy, in the bottom panel of figure 4. We see that the solenoidal ratio drops sharply towards lower Mach numbers in the case of compressive driving. In the absence of the baroclinic term, , the only way to generate vorticity, , with compressive (curl-free) driving is via viscous interactions in the vorticity equation (Mee & Brandenburg 2006),

(2.0)

While the second term on the right hand side is diffusive and dampens , the last term can actually generate vorticity through viscous interactions in the presence of logarithmic density gradients. Thus, even in the absence of initial vorticity, small seeds of can be generated. Those are then exponentially amplified by the non-linear term, , if the Reynolds numbers are sufficiently large (Frisch 1995)The general mechanism for vorticity amplification is thus analogous to the dynamo amplification of small seed fields, because of the similar form of equations (2.3.3) and (2.2.1).. For very low Mach numbers, however, density gradients start to vanish, thus explaining the steep drop of the dynamo growth rate for compressively driven turbulence with low Mach numbers.

Analytic estimates by Moss & Shukurov (1996) suggest that the dynamo growth rate in compressively driven, acoustic turbulence, shown as the dotted line in the top panel of figure 4 (we note that we defined , while Moss & Shukurov (1996) defined , which differs by a factor of from our definition). Our simulations are consistent with this, though simulations at even lower Mach number for compressive driving are required to confirm the scaling.

The solid lines in figure 4 show fits with the empirical function,

(2.0)

The fit parameters are given in table 1. We caution that these fits do not necessarily reflect the true asymptotic behaviour of and , but they provide a reasonably good fit for the range of Mach numbers tested in our simulations.

The subsonic, solenoidally driven models show very high saturation levels, –60%, in excellent agreement with the theoretical predictions from Schober et al. (2015) for Kolmogorov (1941) turbulence (velocity scaling exponent ). On the other hand, for the highly compressible, supersonic limit, Schober et al. (2015) predict times lower saturation levels, of the order of , which is qualitatively consistent with our MHD simulations. The prediction by Schober et al. (2015) for is based on the scaling exponent of supersonic Burgers (1948) turbulence, , which is indeed measured in numerical simulations (Kritsuk et al. 2007; Federrath et al. 2010; Federrath 2013).

The high saturation levels in the case of subsonic, solenoidally driven turbulence cause a strong back reaction of the magnetic field, explaining the Mach number drop in the saturation regime that we saw in figure 1 for these models. This is consistent with the simulations in Haugen et al. (2003, 2004a). For the growth rate, we fix such that for , in good agreement with our models up to . However, simulations with even higher are required to see if holds in this limit. We find that depends much less on in the case of solenoidal driving. Nevertheless, a drop of the growth rate at is present in both driving cases.

Dynamo theories based on the phenomenological model of incompressible, purely solenoidal turbulence by Kolmogorov (1941) predict no dependence of on . For instance, Subramanian (1997) derived based on the Kolmogorov-Fokker-Planck equation, in the limit of large magnetic Prandtl number, with the kinetic and magnetic Reynolds numbers and . For (which is typically the result of numerical dissipation in the ideal MHD approximation; see Lesaffre & Balbus 2007), and , corresponding to our simulations, we find slightly smaller growth rates than predicted in Subramanian (1997), but in agreement with analytic considerations (Boldyrev & Cattaneo 2004), and with numerical simulations of incompressible turbulence for (Schekochihin et al. 2007; Cho et al. 2009). Thus, an extension of the dynamo theory to small is needed and will be presented in §3 below. Moreover, extending the theory from Kolmogorov to Burgers, shock-dominated turbulence is a necessary next step in order to develop a more generalised theory of turbulent dynamos, potentially with predictive power for the supersonic, highly compressible regime.

In summary, we find that the growth rate and saturation level of the turbulent dynamo strongly depend on the Mach number and the driving of the turbulence. We conclude that the compressibility of the plasma plays a crucial role for the amplification of turbulent magnetic fields.

3 Dependence of the dynamo on the magnetic Prandtl number

Previous studies of incompressible turbulence have demonstrated that the turbulent dynamo depends on the magnetic Prandtl number, , defined as the ratio of kinematic viscosity to magnetic diffusivity (Schekochihin et al. 2004). The magnetic Prandtl number also plays a role for the magneto-rotational instability and the transport of angular momentum in accretion disks (Fromang 2010; Fromang et al. 2010).

Here we extend this study of into the compressible, supersonic regime of turbulence. Some of the results discussed in this section are published in Federrath et al. (2014a).

On large cosmological scales and in the interstellar medium, we typically have (for details on how to estimate the viscosity and magnetic diffusivity, see Wardle & Ng 1999; Pinto & Galli 2008; Schober et al. 2012b; Krumholz 2014, for example, in molecular clouds, we have and , hence ), while for the interior of stars and planets, the case with is more relevant (Schekochihin et al. 2007). However, numerical simulations are typically restricted to , because of limited scale separation that can be achieved with currently accessible numerical resolutions. Simulations by Iskakov et al. (2007) have demonstrated that the turbulent dynamo operates for in incompressible gases, even though an asymptotic scaling relation could not be confirmed. While the bulk of previous work was dedicated to exploring the turbulent dynamo for incompressible plasmas (Brandenburg et al. 2012), most astrophysical systems show signs of high compressibility. This is especially true for the early Universe when the first stars and galaxies formed (Latif et al. 2014), in the interstellar medium of present-day galaxies (Larson 1981), as well as in the intergalactic medium (Vazza et al. 2009; Iapichino et al. 2013; Miniati & Beresnyak 2015). As in §2, the compressibility of the plasma is characterised in terms of the sonic Mach number .

Based on the Kazantsev model (Kazantsev 1968; Kazantsev et al. 1985) (see also recent work by Xu & Lazarian (2016)), Schober et al. (2012a) derived an analytical dynamo solution for the limiting cases and , considering different scaling relations of the turbulence. Bovino et al. (2013) later computed numerical solutions of the Kazantsev equation for intermediate values of . These theoretical studies suggested that the turbulent dynamo operates for a range of , as long as the magnetic Reynolds number is sufficiently high. A central restriction of the Kazantsev model is the assumption of an incompressible velocity field, for which a separation into solenoidal and compressible parts is not necessary. But the distinction between solenoidal and compressible modes is essential for highly compressible, supersonic turbulence (see §2). Moreover, the Kazantsev framework assumes that the turbulence is -correlated in time, but the resulting uncertainties introduced by this assumption are only a few per cent (Schekochihin & Kulsrud 2001; Kleeorin et al. 2002; Bhat & Subramanian 2014). However, the assumption of incompressibility may introduce significantly higher uncertainties in the analytic and semi-analytic estimates based on the Kazantsev model. Ultimately one needs full 3D simulations to determine the dependence of the dynamo growth rates and saturation levels on .

Here we present a systematic study of the turbulent dynamo and its dependence on the magnetic Prandtl number in the highly compressible regime. For this purpose, we consider supersonic turbulence with , and magnetic Prandtl numbers between and . The simulation results are compared with analytic and semi-analytic predictions based on the Kazantsev model.

3.1 Numerical simulations

The simulations follow the same basic approach as explained in detail in §2.2. We run most of the simulations until saturation of the magnetic field is reached. We determine the growth rate by fitting an exponential function to the time evolution of the magnetic energy. The saturation level is determined by measuring the magnetic-to-kinetic energy ratio in the saturated phase. We note that the turbulent dynamo has also been studied with “shell models” (Frick et al. 2006, and references therein). Shell models can also provide theoretical predictions for the magnetic energy growth in the exponential and saturated regimes of the dynamo, and they are complementary to the 3D numerical simulations presented here.

Simulation Model
Dyn_512_Pm0.1_Re1600 n/a
Dyn_1024_Pm0.1_Re1600 n/a
Dyn_512_Pm0.2_Re1600
Dyn_512_Pm0.5_Re1600
Dyn_512_Pm2_Re1600
Dyn_256_Pm5_Re1600
Dyn_512_Pm5_Re1600
Dyn_1024_Pm5_Re1600
Dyn_256_Pm10_Re1600
Dyn_512_Pm10_Re1600
Dyn_1024_Pm10_Re1600
Dyn_128_Pm10_Re4.7 n/a
Dyn_256_Pm10_Re4.6 n/a
Dyn_128_Pm10_Re15 n/a
Dyn_256_Pm10_Re15 n/a
Dyn_128_Pm10_Re26
Dyn_256_Pm10_Re26
Dyn_128_Pm10_Re39
Dyn_256_Pm10_Re38
Dyn_512_Pm10_Re38 n/a
Dyn_512_Pm10_Re88
Dyn_512_Pm10_Re190
Dyn_512_Pm10_Re390
Dyn_256_Pm10_Re790
Dyn_512_Pm10_Re790
Table 2: Turbulent dynamo simulations with different magnetic Prandtl number ().

Here we study the dependence of the turbulent dynamo on , which is accomplished by varying the physical viscosity and magnetic diffusivity. Table 2 provides a complete list of all simulations and key parameters. The MHD equations (2.2.1)–(2.2.1) are solved with the positive-definite second-order accurate HLL3R Riemann scheme (Waagan et al. 2011). As a numerical convergence test, we run simulations with grid points.

3.2 Dynamo theory

Most theoretical models of the turbulent dynamo are based on the Kazantsev framework (Kazantsev 1968; Subramanian 1997; Rogachevskii & Kleeorin 1997; Brandenburg & Subramanian 2005), The Kazantsev equation assumes zero helicity, -correlation in time, and does not take into account the mixture of solenoidal-to-compressible modes in the turbulent velocity field. These limitations are related to the fact that the Kazantsev equation was historically only applied to incompressible turbulence. Here we present an extension of the Kazantsev model to compressible, supersonic regime of turbulence (the details of which are published in Schober et al. 2012c, a; Bovino et al. 2013; Schleicher et al. 2013; Schober et al. 2015).

The form of the Kazantsev equation is similar to the Schrödinger equation. This allows us to solve the equation both numerically and analytically with standard methods from quantum mechanics such as the Wentzel-Kramers-Brillouin (WKB) approximation. In order find a solution, we have to assume a scaling of the turbulent velocity fluctuations in the inertial range (),

(3.0)

where and are the viscous and integral scales, respectively. The exponent varies from for incompressible, non-intermittent Kolmogorov (1941) turbulence up to for supersonic, shock-dominated turbulence (Burgers 1948). Based on numerical simulations of mildly supersonic turbulence with Mach numbers , scaling exponents were found in numerical simulations (Boldyrev et al. 2002; Kowal & Lazarian 2010; Federrath et al. 2010). Highly supersonic turbulence with asymptotically approaches the Burgers limit, (Federrath 2013). Observations of interstellar clouds suggest similar velocity scaling exponents with (Larson 1981; Solomon et al. 1987; Ossenkopf & Mac Low 2002; Heyer & Brunt 2004; Roman-Duval et al. 2011). Given this range of turbulence scaling exponents from numerical simulations and observations, we here compute theoretical estimates of the dynamo growth rate for , , and .

3.3 Results and discussion

3.3.1 Magnetic field structure in low- and high- supersonic plasmas

Figure 5: Magnetic energy slices through the mid plane of our dynamo simulations with grid resolutions of points. The magnetic field grows more slowly for low magnetic Prandtl number (left-hand panel) compared to (right-hand panel). However, dynamo action occurs in both cases, and for the first time shown in highly compressible, supersonic plasmas (Federrath et al. 2014a). An animation is available at http://www.mso.anu.edu.au/~chfeder/pubs/dynamo_pm/dynamo_pm.html.

Figure 5 provides a visual comparison of the differences in magnetic energy between low- and high- supersonic plasmas. Magnetic dissipation is stronger in low- compared to high- plasmas (for ). Here we find that the dynamo operates in both cases, but with extremely different efficiency. This is the first time that dynamo action was confirmed in , highly compressible, supersonic plasmas (Federrath et al. 2014a).

3.3.2 Growth rate and saturation level as a function of and

Figure 6: Left panels: dynamo growth rate (top) and saturation level (bottom) as a function of for fixed . Resolution studies with , and grid cells demonstrate convergence, tested for the extreme cases and . Theoretical predictions for by Schober et al. (2012a) and Bovino et al. (2013) and for by Schober et al. (2015) are plotted with different line styles for a typical range of the turbulence scaling exponent, (dotted), (solid) and (dashed). Right panels: same as left panels, but and are shown as a function of for fixed . The dot-dashed line is a fit to the simulations, yielding a constant saturation level of for and the triple-dot-dashed line shows the result of our modified model for the saturation level originally proposed by Subramanian (1999).

In order to compare the analytical and numerical solutions of the Kazantsev equation from §3.2 with the MHD simulations, we now determine the dynamo growth rate as a function of for fixed and as a function of for fixed . Depending on and we find exponential magnetic energy growth over more than 6 orders of magnitude for simulations in which the dynamo operates. We determine both the exponential growth rate and the saturation level . The measurements are listed in table 2 and shown in figure 6.

Dependence on magnetic Prandtl number:

In the left-hand panel of figure 6 we see that first increases strongly with for . For it keeps increasing, but more slowly. The analytic and semi-analytic models by Schober et al. (2012a) and Bovino et al. (2013) both predict an increasing growth rate with . These models are shown for three different turbulence scaling exponents, (dotted lines), (solid lines), and (dashed lines) in equation (3.2). We note that the case with a turbulence scaling exponent (the dashed lines) is the most appropriate here, because the turbulence in all our simulations is supersonic with Mach numbers in the range 4–11, for which previous high-resolution numerical simulations are consistent with (Kritsuk et al. 2007; Federrath et al. 2010). The purely analytical solution of the Kazantsev by Schober et al. (2012a) yields power laws for . However, the semi-analytic solution by Bovino et al. (2013) yields a sharp cutoff for , closer to the results of our 3D MHD simulations.

The theoretical predictions are qualitatively consistent with the MHD simulations. Quantitative discrepancies arise because the theoretical models assume zero helicity, -correlation of the turbulence in time, and do not distinguish solenoidal and compressible modes in the turbulent velocity field. Finite time correlations, however, do not seem to affect the Kazantsev model significantly (Bhat & Subramanian 2014). We therefore conclude that one needs to distinguish between the solenoidal and compressible modes in future dynamo theories, because the dynamo is primarily driven by solenoidal modes and the amount of vorticity strongly depends on the driving and Mach number of the turbulence (Mee & Brandenburg 2006; Federrath et al. 2011a, see §2).

The saturation level as a function of is shown in the bottom left-hand panel of figure 6. It increases with similar to the growth rate and is also well converged with increasing numerical resolution. For , seems to become independent of (and thus independent of , because we have constant here). The theoretical predictions for from Schober et al. (2015) are shown in the low- and high- limits—as for the growth rate, is the most relevant case for this comparison. The analytic prediction agrees qualitatively with the results of the MHD simulations, but similar to the limitations of the theories for the growth rate, more work is needed to incorporate the mode mixture (solenoidal vs. compressive) in the saturation models.

Brandenburg (2014) performed similar simulations and investigated the saturation level, but in the subsonic regime of turbulence (with ), while here we study Mach 10 MHD turbulence. Brandenburg (2014) found that is independent of for fixed . Here instead, we show that increases with for fixed , which implies that grows with for . This is qualitatively consistent with the simulations (e.g., model X3 versus Y7) in Brandenburg (2014).

Dependence on kinematic Reynolds number:

The right-hand panels of figure 6 show the growth rate and saturation level as a function of kinematic Reynolds number . Similar to the dependence on , we find a non-linear increase in with , which is qualitatively reproduced by the semi-analytic solution of the Kazantsev equation in Bovino et al. (2013). However, the critical Reynolds number for dynamo action is much lower in the MHD simulations than predicted by the theoretical models, which may have the same reasons as the discrepancy found for the dependence on , i.e., the lack of distinction of mixed turbulence modes in the Kazantsev model.

The top panels of figure 6 show that the growth rate () depends on both and . To take both dependences into account and to determine the critical magnetic Reynolds number for dynamo action, we perform fits with the empirical model function,

(3.0)

using the fit parameters and , which are related to the critical magnetic Reynolds number . From the fits with equation (3.3.2) to all our simulations we find that dynamo action is suppressed for in highly compressible, supersonic MHD turbulence. Our result for in the supersonic simulations is significantly higher than measured in subsonic, incompressible turbulence by Haugen et al. (2004a), who found for . It is also higher than in mildly compressible turbulence, where for and (Haugen et al. 2004b). The reason for the higher compared to incompressible turbulence could be the sheet-like than vortex-like structure of supersonic turbulence (Boldyrev 2002; Schmidt et al. 2008; Federrath et al. 2009) and the reduced fraction of solenoidal modes (Mee & Brandenburg 2006; Federrath et al. 2010, 2011a, see §2). The difference of the MHD simulations compared to the Kazantsev models is primarily in . Bovino et al. (2013) predicted a much higher for . However, fits to their theoretical model yield , which is in agreement with the range found in our MHD simulations (). This demonstrates that the discrepancy between the simulations and the Kazantsev model is primarily in the critical magnetic Reynolds number, while the qualitative behaviour (determined by the parameter) is reproduced in the Kazantsev model.

The saturation level shown in the bottom right-hand panel of figure 6 is consistent with a constant level of for . Given our measurement of , we can compute the theoretical prediction by Subramanian (1999) for the saturation level, . This is significantly smaller than our simulation result, assuming that is the turbulent crossing time on the largest scales of the system. However, Subramanian noted that the timescale is an “unknown model parameter”. A more appropriate timescale for saturation may be the eddy-turnover timescale on the viscous scale, , for a given turbulent velocity scaling following equation (3.2), because this is where the field saturates first. We find and with , we obtain for a typical range of the velocity scaling exponent from molecular cloud observations (e.g., Larson 1981; Ossenkopf & Mac Low 2002; Heyer & Brunt 2004; Roman-Duval et al. 2011) and simulations of supersonic turbulence (Kritsuk et al. 2007; Schmidt et al. 2009; Federrath et al. 2010; Federrath 2013). The saturation level of our MHD simulations agrees within the uncertainties with our modified version of Subramanian’s saturation model.

3.3.3 Evolution of magnetic energy spectra

Figure 7: Time evolution of the magnetic energy power spectra for simulations with (dotted lines; from bottom to top: , , , , ) and (dashed lines; from bottom to top: , , , , ). Note that for , the last magnetic energy spectrum () has just reached saturation on small scales—the runs did not reach saturation within the limited computing time available, because the growth rates are extremely small for this model; cf. figure 6). The Kazantsev spectrum () is shown as a dash-dotted line for comparison. The solid lines show the time-averaged kinetic energy spectra.

In figure 7 we show the time evolution of the magnetic energy power spectra in our simulations with and and numerical resolution of grid points. They are qualitatively consistent with incompressible dynamo studies (Brandenburg & Subramanian 2005; Mason et al. 2011; Bhat & Subramanian 2013). We see that the power spectra for dissipate on larger scales (lower ) than the spectra, consistent with the theoretical expectation, by a factor of for our . Even for , we see the characteristic increase in the magnetic energy on all scales. The magnetic energy spectra roughly follow the Kazantsev spectrum () on large scales (Kazantsev 1968; Bhat & Subramanian 2014) in the case, but we would expect the same to hold in the case, if our simulations had larger scale separation, i.e., higher numerical resolution. The final spectrum for has just reached saturation on small scales (approaching the kinetic energy spectrum at high ) and continues to grow on larger scales during the non-linear dynamo phase. The runs did not have enough time to reach saturation within the limited available compute time (cf. figure 6), but we expect a non-linear dynamo phase for that is qualitatively similar to the case.

3.4 Summary

In this section we presented a quantitative comparison of the turbulent dynamo in the analytical and semi-analytical Kazantsev models by Subramanian (1999), Schober et al. (2012c, a, 2015) and Bovino et al. (2013) with 3D simulations of supersonic MHD turbulence. We found that the dynamo operates at low and high magnetic Prandtl numbers, but is significantly more efficient for than for . The Kazantsev models agree qualitatively with MHD simulations, but not quantitatively. We attribute the quantitative differences to the fact that the current dynamo theories do not take into account the varying mixture of solenoidal and compressible modes in the velocity field in the case of compressible, supersonic plasmas. An extension in this direction is of high priority, and would lead to theoretical dynamo models with predictive power for the realistic high-Reynolds number regime of many astrophysical plasmas, currently inaccessible to 3D numerical simulations.

4 Turbulent magnetic fields in the presence of ordered guide fields

In this section we present new simulations of turbulent plasmas in the presence of a guide field—an ordered magnetic field component. This is relevant to a number of astrophysical applications such as accretion discs where an ordered magnetic field component is often observed along the rotation axis of the disc at a few scale heights above and below the disc mid plane (Frank et al. 2014). Ordered guide fields are also observed towards molecular clouds in the Milky Way arms (Li et al. 2014) and in the Galactic Centre cloud G0.253+0.016 (Pillai et al. 2015; Federrath et al. 2016b). We note that the ordered field on a particular scale may actually be part of a turbulent field on much larger scales. However, here we are interested in understanding the amplification and evolution of the turbulent magnetic field (pressure) as a function of the guide-field strength, on the same scale.

4.1 Simulation parameters and initial conditions

The simulations follow the same equations (2.2.1)–(2.2.1) as in the dynamo studies from §2 and §3, but here we add stronger ordered guide fields along the axis and systematically vary the strength of . We fix the driving of the turbulence to solenoidal driving ( in equations 2.2.2 and 2.2.2) and we keep the average energy injection rate of the turbulence constant, resulting in a roughly constant Mach number for all simulations. We use normalised, dimensionless values of all basic variables, i.e., a mean density , sound speed , and box length . This gives a constant total mass and a turbulent box crossing time . We set the kinematic viscosity to and the magnetic diffusivity to , which gives a kinematic Reynolds number of and magnetic Reynolds number of . Thus, the magnetic Prandtl number is fixed to . Each simulation was run with a numerical resolution of grid points, but we have also tested a few cases with grid cells and found convergence of our results.

We run a total of 14 simulations. The ordered guide field () is varied over about 5 orders of magnitude from to , which yields Alfvén Mach numbers with respect to the guide field in the range , i.e., covering the full range from sub-Alfvénic to super-Alfvénic turbulence. Since the simulation domain is periodic, the guide-field strength remains globally constant throughout the whole simulation box because of magnetic-flux conservation. The basic simulation parameters and the derived turbulent magnetic field strengths () are listed in table 3.

Simulation Model
256sMS10MA3000
256sMS10MA1000
256sMS10MA300
256sMS10MA100
256sMS10MA50
256sMS10MA20
256sMS10MA10
256sMS10MA5
256sMS10MA2
256sMS10MA1
256sMS10MA05
256sMS10MA02
256sMS10MA01
256sMS10MA005
Table 3: List of turbulence simulations with different guide-field strength ().

4.2 Results and discussion

Figure 8: Turbulent (un-ordered) magnetic field () as a function of (ordered) magnetic guide field () from the simulations (crosses with 1 error bars) listed in table 3. We can distinguish three regimes: i) the dynamo regime for weak , ii) the intermediate regime, and iii) the strong guide-field regime. The dashed and dotted lines indicate fits of for the dynamo regime and the intermediate regime, respectively. The boxes show the theoretical estimate of in the strong-field regime (equation 4.3.2). The long-dashed line shows , which separates the intermediate from the strong-field regime.

As for the simulation models in §2 and §3, we evolve the simulations in table 3 until they reach a converged state in the turbulent magnetic field component (). Figure 8 shows as a function of . We find the remarkable result that although we varied over 5 orders of magnitude, only varies by a factor of . Our measurements of are listed in the last column of table 3.

We identify three different regions in figure 8. First, for low guide fields, we find . This is the ‘dynamo regime’. Second, in the ‘intermediate regime’ we find that increases with as . Finally, in the third regime—the ‘strong-field regime’— decreases with increasing and follows the relation shown by the boxes. The intermediate and strong-field regimes are separated by the line . We now derive the theoretical model, , for the strong guide-field regime.

4.3 Theoretical model for the turbulent magnetic field

We provide a simple analytical model for in the limits of weak and strong guide field, respectively ( and ). We start by separating the magnetic field into an ordered component () and an un-ordered, turbulent component (),

(4.0)
(4.0)

Note that the mean field and , while the magnetic energy density is proportional to .

4.3.1 Weak magnetic guide field

The limit of a weak or vanishing guide field is the dynamo limit that we explored in detail in §2 and §3. In the limit , equation (4.3) becomes and thus the turbulent magnetic energy density . From the models of dynamo saturation (c.f. §2 and §3), we know that the magnetic energy can only reach a fraction of the turbulent kinetic energy because of the back reaction of the field via the Lorentz force. Thus, in the dynamo limit we make the ansatz,

(4.0)
(4.0)
(4.0)
(4.0)

To evaluate this equation for the present case ( and ; see §4.1), we take the measured saturation level from the middle panel of figure 4 for solenoidal driving at Mach 10 and find . This is in excellent agreement with the simulations in the dynamo limit shown in figure 8.

4.3.2 Strong magnetic guide field

In the limit of a strong guide field (), equation (4.3) can be approximated as and the un-ordered (turbulent) magnetic energy density is . We now make the ansatz that this energy is provided by the turbulent kinetic energy via tangling of the ordered field component,

(4.0)
(4.0)
(4.0)