As a matter of force – Systematic biases in idealized turbulence simulations
Abstract
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 smallscale effects including the dynamics of turbulence. These turbulent boxes require an artificial driving in order to mimic energy injection from largescale 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 intime forcing with a constant energy injection leads to a steeper slope in kinetic energy spectrum and less efficient smallscale 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 lineofsight magnetic field from rotation measures, are systematically biased by the driving mechanism. We argue that intime 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.
Philipp Grete
0000000335559886]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 871851189, USA
SAND No: SAND20182657 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
smallscale 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 largescale driving field
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, largescale 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 intime 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 intime forcing with constant energy injection include Stone et al. (1998), Lemaster & Stone (2009), or Kim & Ryu (2005). Examples for intime 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 intime 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 lineofsight 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
(1)  
(2)  
(3) 
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 version
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 intime 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: intime normalized to , intime normalized to , and 1Tcorrelated normalized to . These simulations are referred to as , , and , respectively, and are separately presented in Section 3.4.
3 Results
3.1 Temporal evolution
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.
We attribute this to a less efficient smallscale 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(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 (largescale) random vectors at every single timestep, making it possible that it is locally perpendicular to the largescale 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 largescale forcing patterns, which evolve (or exist) for a reasonable fraction of a dynamical time, leave an imprint on the largescale 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 largescale pattern on the flow, which is demonstrated by flat PDFs as expected.
While the largescale 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 nonzero 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
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 intime 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 lineofsight (LOS) magnetic field strength derived
from rotation measures
(4) 
with being the lineofsight 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 linesofsight) 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 lineofsight 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 largescale intime forcing
While the idea of a intime forcing is appealing on first sight due to its random,
uncorrelated nature, we argue that it is unrealistic for two reasons.
First, no largescale process (on some length scale ) in nature evolves instantaneously
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 intime 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, selfconsistent 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)
(5) 
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. intime 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 lineofsight 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, superAlfvénic regime we are probing is potentially dominated by slow modes, which is inferred from the observed  anticorrelation (Passot, T. & VázquezSemadeni, E., 2003). This suggests a review of the DavisChandrasekharFermi 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 superAlfvénic, i.e., on average kinetic motions dominate magnetic field dynamics. The superAlfvé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 intime 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 intime 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 smallscale 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 lineofsight magnetic field derived from rotation measures in the sub and supersonic regime.
Overall, we argue that largescale intime forcing is neither realistic nor numerically resolved, and conclude that results from simulations with a intime 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.
Footnotes
 In this Letter, we use forcing, driving and acceleration interchangeably to describe a mechanical energy injection process.
 The implementation of the stochastic forcing used in this paper is available at https://github.com/pgrete/AthenaCversion.
 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.
 Smallscale 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 smallscale forcing.
References
 Beck, R., Shukurov, A., Sokoloff, D., & Wielebinski, R. 2003, A&A, 411, 99
 Brandenburg, A., & Dobler, W. 2002, Computer Physics Communications, 147, 471
 Brandenburg, A., & Lazarian, A. 2013, Space Science Reviews, 178, 163
 Brandenburg, A., & Subramanian, K. 2005, Physics Reports, 417, 1
 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
 Burgers, J. 1948, in Advances in Applied Mechanics, ed. R. V. Mises & T. V. KÃ¡rmÃ¡n, Vol. 1 (Elsevier), 171
 Burkhart, B., & Lazarian, A. 2012, The Astrophysical Journal Letters, 755, L19
 Burkhart, B., Stalpes, K., & Collins, D. C. 2017, The Astrophysical Journal Letters, 834, L1
 Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113
 Chernyshov, A., Karelsky, K., & Petrosyan, A. 2012, Flow, Turbulence and Combustion, 89, 563
 Cho, J., Vishniac, E. T., Beresnyak, A., Lazarian, A., & Ryu, D. 2009, The Astrophysical Journal, 693, 1449
 Clark, R. A., Ferziger, J. H., & Reynolds, W. C. 1979, Journal of Fluid Mechanics, 91, 1
 Davis, L. 1951, Phys. Rev., 81, 890
 Eswaran, V., & Pope, S. 1988, Computers & Fluids, 16, 257
 Federrath, C. 2016, Journal of Plasma Physics, 82
 Federrath, C., Klessen, R. S., & Schmidt, W. 2008, The Astrophysical Journal Letters, 688, L79
 Galtier, S. 2016, Introduction to Modern Magnetohydrodynamics (Cambridge University Press)
 Gaspari, M., Ruszkowski, M., & Sharma, P. 2012, The Astrophysical Journal, 746, 94
 Germano, M., Piomelli, U., Moin, P., & Cabot, W. H. 1991, Physics of Fluids A, 3, 1760
 Grete, P., O’Shea, B. W., Beckwith, K., Schmidt, W., & Christlieb, A. 2017, Physics of Plasmas, 24, 092311
 Grete, P., Vlaykov, D. G., Schmidt, W., & Schleicher, D. R. G. 2016, Physics of Plasmas, 23
 Kida, S., & Orszag, S. A. 1990, Journal of Scientific Computing, 5, 85
 Kim, J., & Ryu, D. 2005, The Astrophysical Journal Letters, 630, L45
 Kroupp, E., Stambulchik, E., Starobinets, A., et al. 2018, Phys. Rev. E, 97, 013202
 Lemaster, M. N., & Stone, J. M. 2009, The Astrophysical Journal, 691, 1092
 Li, Y., Bryan, G. L., Ruszkowski, M., et al. 2015, The Astrophysical Journal, 811, 73
 Meece, G. R., Voit, G. M., & OâShea, B. W. 2017, The Astrophysical Journal, 841, 133
 Nolan, C. A., Federrath, C., & Sutherland, R. S. 2015, Monthly Notices of the Royal Astronomical Society, 451, 1380
 Orkisz, J. H., Pety, J., Gerin, M., et al. 2017, A&A, 599, A99
 Passot, T., & VázquezSemadeni, E. 2003, A&A, 398, 845
 Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, Astronomy & Astrophysics, 494, 127
 Stone, J. M., & Gardiner, T. 2009, New Astronomy, 14, 139
 Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, The Astrophysical Journal Supplement Series, 178, 137
 Stone, J. M., Ostriker, E. C., & Gammie, C. F. 1998, The Astrophysical Journal Letters, 508, L99
 Tobias, S. M., Cattaneo, F., & Boldyrev, S. 2013, in Ten Chapters in Turbulence, ed. P. Davidson, Y. Kaneda, & K. Sreenivasan (Cambridge University Press), 351
 Towns, J., Cockerill, T., Dahan, M., et al. 2014, Computing in Science & Engineering, 16, 62
 Wu, Q., Kim, J., & Ryu, D. 2015, New Astronomy, 34, 21
 Yang, Y., Shi, Y., Wan, M., Matthaeus, W. H., & Chen, S. 2016, Phys. Rev. E, 93, 061102
 Yoon, H., Cho, J., & Kim, J. 2016, The Astrophysical Journal, 831, 85