Biases in turbulence simulations

As a matter of force – Systematic biases in idealized turbulence simulations


Many astrophysical systems encompass very large dynamical ranges in space and time, which are not accessible by direct numerical simulations. Thus, idealized subvolumes are often used to study small-scale effects including the dynamics of turbulence. These turbulent boxes require an artificial driving in order to mimic energy injection from large-scale processes. In this Letter, we show and quantify how the autocorrelation time of the driving and its normalization systematically change properties of an isothermal compressible magnetohydrodynamic flow in the sub- and supersonic regime and affect astrophysical observations such as Faraday rotation. For example, we find that -in-time forcing with a constant energy injection leads to a steeper slope in kinetic energy spectrum and less efficient small-scale dynamo action. In general, we show that shorter autocorrelation times require more power in the acceleration field, which results in more power in compressive modes that weaken the anticorrelation between density and magnetic field strength. Thus, derived observables, such as the line-of-sight magnetic field from rotation measures, are systematically biased by the driving mechanism. We argue that -in-time forcing is unrealistic and numerically unresolved, and conclude that special care needs to be taken in interpreting observational results based on the use of idealized simulations.

MHD — methods: numerical — turbulence

Philipp Grete

0000-0003-3555-9886]Philipp Grete \affiliation Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA \affiliation Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA \affiliation Department of Computational Mathematics, Science and Engineering, Michigan State University, East Lansing, MI 48824, USA \affiliation National Superconducting Cyclotron Laboratory, Michigan State University, East Lansing, MI 48824, USA


Sandia National Laboratories, Albuquerque, NM 87185-1189, USA


SAND No: SAND2018-2657 J

1 Introduction

Many astrophysical systems are governed by compressible magnetohydrodynamic (MHD) dynamics on macroscopic scales (Galtier, 2016). Moreover, given the large scales involved astrophysical systems are often in a turbulent state (Brandenburg & Lazarian, 2013). Compressible MHD turbulence itself is expected to be a major factor in many processes such as magnetic field amplification via the turbulent dynamo (Brandenburg & Subramanian, 2005; Tobias et al., 2013) and particle acceleration in shock fronts which can eventually be observed as cosmic rays (Brunetti & Jones, 2015).

Astrophysical observations of distant systems are typically a two dimensional map that measure quantities along the third, integrated dimension. In order to interpret those observations numerical simulations are often used as they provide detailed data in four (3 spatial and a temporal) dimensions (Burkhart & Lazarian, 2012; Burkhart et al., 2017; Orkisz et al., 2017). However, the large dynamical range in space and time often prohibit direct numerical simulation of an entire system. Thus, idealized subvolumes are used to study specific effects and processes including small-scale turbulent dynamics in so called turbulent boxes. Turbulent boxes are the workhorse in both astrophysical and general turbulence research, and are used to study a variety of aspects of turbulence including energy transfers in the hydrodynamic (HD) (Kida & Orszag, 1990) and MHD (Yang et al., 2016; Grete et al., 2017) case, or HD (Clark et al., 1979; Germano et al., 1991) and MHD (Chernyshov et al., 2012; Grete et al., 2016) turbulence models. To reach a state of turbulence in simulations a large-scale driving field1 is commonly used.

One fundamental assumption in these simulations is that the dynamics on the large scales (where energy in injected) are decoupled from the dynamics on smaller scales. In other words, large-scale features are lost in the energy cascade towards small scales. We show that this assumption is wrong even for simple, purely solenoidal driving.

Many different driving schemes are used in numerical turbulence research. Here, we focus on schemes that are popular in the astrophysical turbulence community and often correspond to the default turbulence in a box setup in many codes. These schemes can be differentiated by two main properties: the autocorrelation time and the normalization applied to the driving field. The autocorrelation time determines on which timescale the driving field evolves. Most commonly two extreme cases are used. On the one hand, a -in-time forcing calculates a new random driving field on each timestep that is completely uncorrelated to the driving field of the previous timestep. On the other hand, smoothly evolving driving fields are used with a given autocorrelation time (often set to the dynamical time of the simulation) realized, e.g., by a stochastic process (Eswaran & Pope, 1988). In both cases the driving field undergoes a random change on each timestep, which leads to a random change in its power. Hence, the second differentiating property is how the amplitude of the acceleration field is normalized on each timestep. Again, two possibilities are commonly used. First, the driving field is normalized to have constant power over time, i.e., the root mean square (RMS) value is constant. Second, the driving field is normalized to have a constant energy injection rate throughout the simulation. Examples for -in-time forcing with constant energy injection include Stone et al. (1998), Lemaster & Stone (2009), or Kim & Ryu (2005). Examples for -in-time forcing with constant power include Brandenburg & Dobler (2002). Examples using a driving field that evolves on a dynamical timescale with constant power include Cho et al. (2009), Federrath et al. (2008), and Schmidt et al. (2009).

A previous study by Yoon et al. (2016) analyzed simulations with two different driving mechanisms: a -in-time forcing normalized to , and a driving field with a finite correlation time and constant power. They find differences in the correlation between density and magnetic field strength and in statistical moments of density related fields including the probability density function, the dispersion measure and Faraday rotation measure. Yoon et al. (2016) attribute those differences to a link between the autocorrelation time of the driving and the ability of the system to reach pressure equilibrium. Here, we go one step further and show that the autocorrelation time is only a secondary parameter. The primary driver in the observed differences is the power in the compressive modes of the resulting flow.

In particular, we show that the power in the acceleration field, which varies with the autocorrelation time for similar stationary regimes, directly affects compressive modes in the simulation. Varying compressive power, in turn, results in different statistical properties such as the slope in the kinetic energy spectrum or the correlation between density and magnetic field strength. Moreover, these differences manifest in changing observable quantities, e.g., Faraday rotation measures. Thus, a systematic bias is introduced when simulations are used as a basis to interpret observations such as the line-of-sight magnetic field strength.

This Letter is organized as follows. In Section 2, we introduce the simulations and the implementation of the different driving mechanisms. In Section 3, we present the key results and differences in the simulations, and discuss the implications in Section 4. Finally, we conclude in Section 5 and provide future directions.

2 Method

In this work we are dealing with the compressible, ideal MHD equations


that are closed by an isothermal equation of state. The symbols have their usual meaning, i.e., density , velocity , thermal pressure , and magnetic field , which includes a factor . Vector quantities that are not in boldface refer to the norm of the vector and denotes the outer product. The details of the acceleration field  that we use to mechanically drive our simulations are described below in Section 2.2.

2.1 Simulations

We use a modified version2 of the astrophysical MHD code Athena 4.2 (Stone et al., 2008). All simulations use second order reconstruction with slope-limiting in the primitive variables, the HLLD Riemann solver, constrained transport for the magnetic field, and the MUSCL-Hancock integrator (Stone & Gardiner, 2009) on a uniform, static grid with cells. We start with uniform initial conditions (all in code units) , , and (subsonic) or (supersonic) corresponding to initial plasma betas of 72 and 4.5, respectively. The two different initial lead to the the same Alfvénic Mach number in both regimes. We evolve the system for five dynamical times . Here, is the characteristic velocity in the stationary phase, which corresponds to the root mean square (RMS) sonic Mach number as we fix the isothermal sound speed to 1. In the subsonic simulations and in the supersonic simulations . The characteristic length L is half the box size (), because our acceleration spectrum is parabolic (Schmidt et al., 2009) and peaks at (using normalized wavenumbers). The acceleration field is purely solenoidal, i.e., . We store 20 equidistant snapshots per dynamical time. The stationary phase is reached after approximately T, and we calculate statistical properties for this phase based on 50 snapshots between .

2.2 Forcing mechanisms

In order to quantify the influence of the autocorrelation time of the driving field we implemented a stochastic forcing mechanism as presented by Schmidt et al. (2009) in Athena, and conduct four identical simulation in the subsonic regime that vary only in their autocorrelation time (and their RMS acceleration value as explained in subsection 3.2). We identify these simulations with , , , and  corresponding to correlation times of T, T, T, and T, and of , , and , respectively. The last simulation is effectively -in-time correlated as the smallest timestep in the simulation is .

This allows a direct comparison to the existing forcing implementation in Athena. It produces -correlated realizations that are normalized by the energy input rate. We conduct one simulation in the subsonic regime with this mechanism and set in order to reach the same RMS sonic Mach number as the other simulations. We refer to this simulation as throughout the paper.

To verify our findings in a more compressive regime we conduct supersonic simulations of the three corner cases: -in-time normalized to  , -in-time normalized to  , and 1T-correlated normalized to  . These simulations are referred to as , , and , respectively, and are separately presented in Section 3.4.

3 Results

3.1 Temporal evolution

Figure 1: Temporal evolution of the spatial root mean square (RMS) sonic Mach number (a) and Alfvénic Mach number (b). The gray area between indicates the temporal range we use as stationary regime throughout the paper.

All subsonic simulations reach a stationary regime with a sonic Mach number after T, as illustrated in Fig. 1(a). In contrast to this, the Alfvénic Mach numbers varies substantially between different normalizations. While is in the stationary regime, it reaches only for all other simulations, shown in Fig. 1(b). Given that is a proxy for the ratio of kinetic to magnetic energy, a higher value for corresponds to a lower saturation value of the magnetic energy.

Figure 2: Compressive (), rotational , and total (i.e., compressive plus rotational) kinetic energy spectra compensated by . The total energy spectra are shifted vertically by a factor of 100 for clarity. All lines correspond to the temporal mean during the stationary phase and the shaded areas indicate the standard deviation over time. The energy spectra are calculated based on the Fourier transforms of and are virtually identical to the ones based on .

We attribute this to a less efficient small-scale dynamo, which is driven by rotational motion. Figure 2 shows the mean kinetic energy spectra (rotational, compressive, and total) in the stationary regime based on the Helmholtz decomposition of the velocity field with and . The simulation has significantly less power (% compared to the other simulations) in rotational modes on scales smaller than the injection scale. Similarly, there is much more power in compressive motions in that simulation, to the degree that the total (compressive plus rotational) kinetic energy spectrum exhibits a different slope in the power law regime . It is steeper as expected from a strongly compressive turbulence phenomenology (Burgers, 1948).

The rotational and total energy spectra for the simulations , , , and are all virtually identical given that the compressive modes are weaker by at least one order of magnitude. However, there are clear differences in the compressive power. With decreasing correlation time (and, thus, increasing ), there is more power in the compressive modes. This can be explained by a more detailed analysis of the link between and in the following subsection.

3.2 Linking , and compressive modes

In all simulations, the driving field is purely solenoidal, i.e., it carries no compressive power itself. Nevertheless, a strong (high ) driving is expected to seed compressive modes, c.f., the canonical idea of wave steepening.

Figure 3: Temporal evolution of the spatial root mean square (RMS) acceleration (a), and the temporal mean PDF of the cosine of the angle between the velocity/momentum field and the acceleration field (b). All lines in panel (b) correspond to the temporal mean during the stationary phase (gray area) and the shaded colored areas indicate the standard deviation over time.

Figure 3(a) shows for all subsonic simulations. By construction is constant for , , , and , while it varies strongly over more than two orders of magnitude for , which is normalized for a constant energy injection rate . To keep at a given value can take arbitrarily large (positive) values for the following reason. The driving field consists of new (large-scale) random vectors at every single timestep, making it possible that it is locally perpendicular to the large-scale velocity field in the entire box. In this case, the resulting energy injection (based on the scalar product of and ) would be negligible for “small” . Thus, during the normalization step is increased (decreased) to match a desired for an predominantly perpendicular (aligned) combination of flow configuration and acceleration field.

The same mechanism, i.e., the alignment of and , is responsible for requiring a higher value for smaller to reach the same in the other simulations. Figure 3(b) illustrates the mean probability density function (PDF) of the angle between and . A larger correlation time results in a distribution for which there is a tendency of and being more aligned. This illustrates how large-scale forcing patterns, which evolve (or exist) for a reasonable fraction of a dynamical time, leave an imprint on the large-scale flow pattern. Despite its clear signal in the PDFs, this alignment should not be overrated, because for our chosen binning a perfect alignment in the entire box would correspond to a -peak with a value of 64. Nevertheless, it is enough to influence the energy injection efficiency (via the local scalar product between and ) and requiring larger for lower (Eswaran & Pope, 1988). For the extreme cases, and , there is no imprint of the large-scale pattern on the flow, which is demonstrated by flat PDFs as expected.

While the large-scale imprint vanishes for smaller , another feature is introduced to the flow by larger : compressive modes. At locations where and are aligned, which is always the case given the non-zero PDF around in Fig. 3(b), large lead to a strong acceleration, resulting in immediate downstream compression. A careful examination of the compressive power spectra in Fig. 2 reveals signatures of this effect. The compressive power of peaks at , which corresponds to the smallest scales in the power spectrum of the acceleration field. Additional signatures of this compression effect are also visible in statistics of the density field, as illustrated in the following subsection.

3.3 Density field dynamics

Figure 4: Temporal evolution of the correlation between density and magnetic field magnitude (a,d), and temporal mean PDFs of the logarithmic density (b,e) and line of sight magnetic field (c,f). Subsonic simulations are shown in the left panels and supersonic simulations in the panels on the right. All lines in the bottom two rows correspond to the temporal mean during the stationary phase (gray area) and the shaded colored areas indicate the standard deviation over time. The true line-of-sight magnetic field strength is illustrated by the vertical dashed line in the two bottom panels.

A link between the Pearson correlation coefficient between density and magnetic field strength, , and the correlation time of the forcing has already been recognized by Yoon et al. (2016). However, the correlation time is only one of two integral parts, see Fig. 4(a), which shows the correlation coefficient over time for all subsonic simulations. , , , and build one family for which the correlation coefficient in the stationary regime increases from , , , to , respectively. In other words, there exists a strong anticorrelation for a T that decreases with smaller towards to a still significant, non zero value of for -in-time forcing (). In contrast, exhibits virtually no correlation between the density field and magnetic field strength . This is in agreement with the results of Yoon et al. (2016), who analyzed two forcing configurations corresponding to our and cases. Given our additional simulations with varying , we argue that the power in the acceleration field and not is the primary driver of changes in . Otherwise, and should yield identical results, which is not observed here. Nevertheless, and are tightly linked as shown in the previous subsection.

Clear differences between the simulations are also observed in the logarithmic density PDFs, as illustrated in Fig. 4(b). exhibits a pronounced negative skew, i.e., the low density tail is longer so that lower than average density values are more likely than higher values. While the skewness decreases with decreasing it it still present in the simulation. Again, contrasts with these results showing an almost symmetrical distribution. In general, the high density tails (with allowing for the most and for the least extreme value) nicely illustrate how larger lead to immediate compression.

Finally, in order to illustrate how these differences translate to biases in interpretations of astrophysical observations, we calculate line-of-sight (LOS) magnetic field strength derived from rotation measures3 as


with being the line-of-sight component of the magnetic field. Beck, R. et al. (2003) derived how an (anti)correlation between and changes this measurement of the mean magnetic field strength, which is exact for uncorrelated fields. For anticorrelated fields Eq. (4) underestimates the true LOS . This effect can be observed in Fig. 4(c) where the mean PDFs of the LOS  (over the available lines-of-sight) are shown. With increasing anticorrelation from , , , , to the derived values of , , , , and , respectively, deviate further from the real value . More strikingly, the PDF of is much more peaked with a standard deviation of compared to , , , and with , , , and , respectively.

3.4 Supersonic simulations

The supersonic simulations , , and exhibit a very similar behavior to their subsonic counterparts. All simulations reach a stationary regime after 2.5T with the same , but different of , , and , respectively. Again, a higher for similar in the saturated regime implies less effective magnetic field amplification in the case. Similarly, the ratio of compressive versus rotational power in the kinetic energy spectrum decreases from to to . However, the differences are less pronounced compared to the subsonic regime given that there is overall more power in compressive modes as expected from the supersonic regime. The dynamical alignment between velocity and acceleration field in the supersonic regime is virtually identical to the corresponding subsonic simulations shown in Fig. 3(b).

Quantitative differences between both regimes are first observed in the correlation as ilustrated in Fig. 4(d). and exhibit virtually no correlation with correlation coefficients of and whereas still shows a weak anticorrelation with a coefficient of . The logarithmic density PDFs of the supersonic simulations in Fig. 4(e) follow the same trend observed in the subsonic simulations. With decreasing power in the acceleration field the PDFs get a more pronounced negative skew. Finally, derived line-of-sight magnetic field strength measurements are again systematically affected as shown in Fig. 4(f). The derived value in the simulation of is closest the real value of and the PDF is most peaked with a standard deviation of whereas and underestimate the magnetic field strength with and , respectively, and generally broader PDFs with deviations of and , respectively.

4 Discussion

4.1 Unrealistic large-scale -in-time forcing

While the idea of a -in-time forcing is appealing on first sight due to its random, uncorrelated nature, we argue that it is unrealistic for two reasons. First, no large-scale process (on some length scale ) in nature evolves instantaneously4. For a -in-time evolution with , the characteristic velocity of that process is . This leads to the second, numerical argument. A -in-time forcing in a numerical simulation, by construction, is not resolving the physical timescale. The timestep in a simulation (of the type discussed in this Letter) is restricted so that information locally travels no further than to adjacent cells. Thus, with the required timestep . This restriction can never be satisfied. Therefore, a large-scale -in-time forcing is never numerically resolved.

4.2 Forcing normalizations

Similar to the autocorrelation time discussion, choosing a normalization to a constant energy injection rate over a constant RMS acceleration is appealing at first glance. From a turbulence analysis point of view automatically fixes the energy dissipation rate in the simulation. However, in allowing the flow to reach a stationary state that has no realistic counterpart, it also masks the effects of using an unresolved -in-time forcing.

For finite autocorrelation times, the practical choice between normalizing by and is less important. Normalizing by for the stochastic forcing used here naturally leads to a statistical constant energy injection rate if is adjusted appropriately (Eswaran & Pope, 1988). In fact, both approaches can mimic realistic processes depending on the feedback mechanism. On the one hand, normalizing by can be seen as an external, self-consistent process that regulates itself without significant feedback from the environment, for example, energy injection from a jet. On the other hand, normalizing by can be seen as a process that depends on the interaction with the environment, for example, cold mode AGN accretion (Gaspari et al., 2012; Li et al., 2015; Meece et al., 2017).

4.3 Total pressure equilibrium and - correlations

The strong anticorrelation between the density and the magnetic field strength observed in the run is consistent with the expectation of a (statistical) total pressure equilibrium (Beck, R. et al., 2003)


In addition, the isothermal equation of state used in the simulations mandates . Hence, to maintain a constant total pressure low density regions correspond to regions with higher magnetic field strength and vice versa. Yoon et al. (2016) also followed this reasoning and explained the lack of anticorrelation in their equivalent simulation by the forcing timescale. -in-time forcing evolves so fast that the system is not able to reach pressure equilibrium. Thus, should also exhibit no significant anticorrelation. However, we observe a moderate anticorrelation in that simulation, arguing that the autocorrelation timescale alone is not sufficient to explain the results.

We argue that the increasing power in compressive modes injected by larger (and, thus, smaller ) is the main driver behind a decreasing - anticorrelation. Given Alfvén’s theorem in MHD, compression naturally leads to a positive correlation between and . Any compression that is locally not exactly aligned with the magnetic field direction compresses the magnetic field in the other two directions, which results in an increased magnetic flux. Thus, the total pressure equilibrium induced strong - anticorrelation is successively weakened by increasing compressive modes associated with a positive - correlation. This is also in agreement with the supersonic simulations, which naturally have more power in compressive modes.

4.4 Observational consequences

An empirical relation between the derived line-of-sight magnetic field strength, the sonic Mach number, and the widths of the rotation measure distribution was suggested by Wu et al. (2015). This relation targets rapid estimates of the LOS magnetic field in the turbulent warm ionized medium in our Galaxy. However, the relation is based on simulations employing forcing and assumes no - correlation as reported by Wu et al. (2015). While there are differences in our -normalized simulations, the differences are overall much less pronounced compared to what is observed in . Thus, a modified LOS estimate could still be obtained and will be explored in future work.

Similarly, magnetic field estimates in the plane of the sky as proposed by Davis (1951) and Chandrasekhar & Fermi (1953) are potentially affected by the observed - anticorrelation. Their estimate assumes underlying isotropic Alfvénic perturbations of the flow. However, the subsonic, super-Alfvénic regime we are probing is potentially dominated by slow modes, which is inferred from the observed - anticorrelation (Passot, T. & Vázquez-Semadeni, E., 2003). This suggests a review of the Davis-Chandrasekhar-Fermi method for different regimes.

4.5 Limitations

The main purpose of the present study is to highlight and explain observed differences in turbulence simulations employing one of the most commonly used idealized setups: isothermal, solenoidally driven, stationary turbulence. Independent of the driving scheme analysis, our results indicate that compressibility and pressure dynamics leave a clear imprint on the flow and derived observables. Thus, other factors are also expected to be dynamically relevant, in particular compressive modes in the acceleration field (Federrath, 2016), and an adiabatic equation of state (Nolan et al., 2015).

Moreover, all our simulations are super-Alfvénic, i.e., on average kinetic motions dominate magnetic field dynamics. The super-Alfvénic regime could be relaxed in two directions. On the one hand, increasing the background magnetic field strength decreases the Alfvénic mach number . While Yoon et al. (2016) conducted and analyzed simulations with varying , the interplay between compressive modes and a dynamically important background field remains open. On the other hand, we expect to see magnetic field unrelated features, for example, the link between and compressive modes, or the increasing alignment of velocity and acceleration field with increasing , also in the pure hydrodynamic case.

We leave a more detailed analysis of effects pertaining to compressibility, the equation of state, and the strength of a background magnetic field to future work.

5 Conclusions

In this Letter, we studied the effects of driving parameters on statistical quantities and observables in stationary, isothermal, compressible MHD turbulence simulations. All simulation were driven to reach the same subsonic (supersonic) Mach number of () with varying autocorrelation time and normalization of the acceleration field. We varied the autocorrelation time between T dynamical times, i.e., between an effective -in-time forcing and a forcing field that evolves on the dynamical timescale of the flow. In these simulations (identified with , , , and ) the acceleration field was normalized so that the power in the acceleration field is constant over time where shorter necessitate higher . For the -in-time case, we also varied to normalization in order to keep the energy injection rate exactly constant over time (). This allowed for varying power in the acceleration field given the current flow configuration.

Our main findings are the following:

  • With increasing more power is injected in compressible modes even though the acceleration field itself is purely solenoidal. In the most extreme case, , the resulting power in compressive modes becomes dynamically relevant so that the total kinetic energy spectrum exhibits a steeper slope compared to the other simulations. In addition, the relatively weaker rotational modes in result in less efficient small-scale dynamo action and, in turn, a lower magnetic energy saturation value.

  • With increasing energy injection is more efficient (and, thus, require lower ) because the acceleration and velocity field are more aligned.

  • In the subsonic regime, there is a strong anticorrelation (correlation coefficient of -0.81) between density and magnetic field strength for . With decreasing (increasing ) the anticorrelation constantly weakens down to -0.5 for . Effectively no correlation is observed for . The anticorrelation itself in the present regime is expected from a total pressure equilibrium. We attribute the decreasing anticorrelation from increasing compressive modes, which are associated with a positive - correlation due to Alfvén’s frozen in flux theorem.

  • In the supersonic regime, correlation coefficients are generally higher, which supports the argument for the importance of compressive modes.

  • We confirm that the presence of a - anticorrelation leads to bias in observables, for example, an underestimated and more variable line-of-sight magnetic field derived from rotation measures in the sub- and supersonic regime.

Overall, we argue that large-scale -in-time forcing is neither realistic nor numerically resolved, and conclude that results from simulations with a -in-time forcing should be interpreted with care. As such, our findings have implications for observations of magnetized turbulence in both astrophysical and terrestrial environments, in particular those that depend on probability distribution functions of fluid quantities (e.g., Kroupp et al., 2018) or on correlations between turbulent fluctuations (e.g., Wu et al., 2015). A more detailed analysis of the interplay between compressive modes and pressure fluctuations in the context of observations will be presented in future work.

The authors thank David Collins, Alexei Kritsuk, Jeffrey Oishi, and Wolfram Schmidt for useful discussions. PG and BWO acknowledge funding by NASA Astrophysics Theory Program grant #NNX15AP39G. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. BWO acknowledges additional funding by NSF AAG grant #1514700. The simulations were run on the NASA Pleiades supercomputer through allocation SMD-16-7720 and on the Comet supercomputer as part of the Extreme Science and Engineering Discovery Environment (XSEDE Towns et al., 2014), which is supported by National Science Foundation grant number ACI-1548562, through allocation #TG-AST090040. Athena is developed by a large number of independent researchers from numerous institutions around the world. Their commitment to open science has helped make this work possible.


  1. In this Letter, we use forcing, driving and acceleration interchangeably to describe a mechanical energy injection process.
  2. The implementation of the stochastic forcing used in this paper is available at
  3. Technically, the relation is valid for the number density of thermal electrons, but given the isothermal single fluid MHD approximation employed we use the fluid density instead.
  4. Small-scale processes, for example, energy injection from supernovae within a galaxy, can occur almost instantaneously on the dynamical time of the galaxy. However, this corresponds to small-scale forcing.


  1. Beck, R., Shukurov, A., Sokoloff, D., & Wielebinski, R. 2003, A&A, 411, 99
  2. Brandenburg, A., & Dobler, W. 2002, Computer Physics Communications, 147, 471
  3. Brandenburg, A., & Lazarian, A. 2013, Space Science Reviews, 178, 163
  4. Brandenburg, A., & Subramanian, K. 2005, Physics Reports, 417, 1
  5. Brunetti, G., & Jones, T. W. 2015, Cosmic Rays in Galaxy Clusters and Their Interaction with Magnetic Fields, ed. A. Lazarian, M. E. de Gouveia Dal Pino, & C. Melioli (Berlin, Heidelberg: Springer Berlin Heidelberg), 557
  6. Burgers, J. 1948, in Advances in Applied Mechanics, ed. R. V. Mises & T. V. Kármán, Vol. 1 (Elsevier), 171
  7. Burkhart, B., & Lazarian, A. 2012, The Astrophysical Journal Letters, 755, L19
  8. Burkhart, B., Stalpes, K., & Collins, D. C. 2017, The Astrophysical Journal Letters, 834, L1
  9. Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113
  10. Chernyshov, A., Karelsky, K., & Petrosyan, A. 2012, Flow, Turbulence and Combustion, 89, 563
  11. Cho, J., Vishniac, E. T., Beresnyak, A., Lazarian, A., & Ryu, D. 2009, The Astrophysical Journal, 693, 1449
  12. Clark, R. A., Ferziger, J. H., & Reynolds, W. C. 1979, Journal of Fluid Mechanics, 91, 1
  13. Davis, L. 1951, Phys. Rev., 81, 890
  14. Eswaran, V., & Pope, S. 1988, Computers & Fluids, 16, 257
  15. Federrath, C. 2016, Journal of Plasma Physics, 82
  16. Federrath, C., Klessen, R. S., & Schmidt, W. 2008, The Astrophysical Journal Letters, 688, L79
  17. Galtier, S. 2016, Introduction to Modern Magnetohydrodynamics (Cambridge University Press)
  18. Gaspari, M., Ruszkowski, M., & Sharma, P. 2012, The Astrophysical Journal, 746, 94
  19. Germano, M., Piomelli, U., Moin, P., & Cabot, W. H. 1991, Physics of Fluids A, 3, 1760
  20. Grete, P., O’Shea, B. W., Beckwith, K., Schmidt, W., & Christlieb, A. 2017, Physics of Plasmas, 24, 092311
  21. Grete, P., Vlaykov, D. G., Schmidt, W., & Schleicher, D. R. G. 2016, Physics of Plasmas, 23
  22. Kida, S., & Orszag, S. A. 1990, Journal of Scientific Computing, 5, 85
  23. Kim, J., & Ryu, D. 2005, The Astrophysical Journal Letters, 630, L45
  24. Kroupp, E., Stambulchik, E., Starobinets, A., et al. 2018, Phys. Rev. E, 97, 013202
  25. Lemaster, M. N., & Stone, J. M. 2009, The Astrophysical Journal, 691, 1092
  26. Li, Y., Bryan, G. L., Ruszkowski, M., et al. 2015, The Astrophysical Journal, 811, 73
  27. Meece, G. R., Voit, G. M., & O’Shea, B. W. 2017, The Astrophysical Journal, 841, 133
  28. Nolan, C. A., Federrath, C., & Sutherland, R. S. 2015, Monthly Notices of the Royal Astronomical Society, 451, 1380
  29. Orkisz, J. H., Pety, J., Gerin, M., et al. 2017, A&A, 599, A99
  30. Passot, T., & Vázquez-Semadeni, E. 2003, A&A, 398, 845
  31. Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, Astronomy & Astrophysics, 494, 127
  32. Stone, J. M., & Gardiner, T. 2009, New Astronomy, 14, 139
  33. Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, The Astrophysical Journal Supplement Series, 178, 137
  34. Stone, J. M., Ostriker, E. C., & Gammie, C. F. 1998, The Astrophysical Journal Letters, 508, L99
  35. Tobias, S. M., Cattaneo, F., & Boldyrev, S. 2013, in Ten Chapters in Turbulence, ed. P. Davidson, Y. Kaneda, & K. Sreenivasan (Cambridge University Press), 351
  36. Towns, J., Cockerill, T., Dahan, M., et al. 2014, Computing in Science & Engineering, 16, 62
  37. Wu, Q., Kim, J., & Ryu, D. 2015, New Astronomy, 34, 21
  38. Yang, Y., Shi, Y., Wan, M., Matthaeus, W. H., & Chen, S. 2016, Phys. Rev. E, 93, 061102
  39. Yoon, H., Cho, J., & Kim, J. 2016, The Astrophysical Journal, 831, 85
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description