As a Matter of Force—Systematic Biases in Idealized Turbulence Simulations

We use a modified version⁶ 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 512³ cells. We start with uniform initial conditions (all in code units) ρ = 1, ${\boldsymbol{u}}={\boldsymbol{0}}$, and ${\boldsymbol{B}}={(1/6,0,0)}^{T}$ (subsonic) or ${\boldsymbol{B}}\,={(2/3,0,0)}^{T}$ (supersonic) corresponding to initial plasma betas of 72 and 4.5, respectively. The two different initial ${\boldsymbol{B}}$ lead to the the same Alfvénic Mach number in both regimes. We evolve the system for five dynamical times T = V/0.5L. Here, V is the characteristic velocity in the stationary phase, which corresponds to the rms sonic Mach number ${M}_{s}$ as we fix the isothermal sound speed to 1. In the subsonic simulations V = 0.5, and in the supersonic simulations V = 2. The characteristic length 0.5L is half of the box size (L = 1), because our acceleration spectrum is parabolic (Schmidt et al. 2009) and peaks at k = 2 (using normalized wavenumbers). The acceleration field is purely solenoidal, i.e., ${\rm{\nabla }}\cdot {\boldsymbol{a}}=0$. We store 20 equidistant snapshots per dynamical time. The stationary phase is reached after approximately 2.5T, and we calculate statistical properties for this phase based on 50 snapshots between 2.5T < t < 5T.

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 simulations in the subsonic regime that vary only in their autocorrelation time (and their rms acceleration value $\langle a\rangle$ as explained in Section 3.2). We identify these simulations with ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{4}}{\mathtt{T}}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{16}}{\mathtt{T}}}$, and ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$ corresponding to correlation times of 1T, 0.25T, 0.0625T, and 10⁻⁹T, and $\langle a\rangle$ of 1, $\sqrt{2}$, 2, and 40, respectively. The last simulation is effectively δ-in-time correlated, as the smallest timestep in the simulation is  >10⁻⁵.

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 $\dot{E}=0.1$ in order to reach the same rms sonic Mach number as the other simulations. We refer to this simulation as ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ throughout the Letter.

To verify our findings in a more compressive regime we conduct supersonic simulations of the three corner cases: δ-in-time normalized to $\dot{E}(=6.4)$, δ-in-time normalized to $\langle a\rangle \,(=1000)$, and 1T-correlated normalized to $\langle a\rangle \,(=16)$. These simulations are referred to as ${{\mathtt{F}}}_{\dot{E}}^{\delta }{\mathtt{M}}{\mathtt{2}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }{\mathtt{M}}{\mathtt{2}}$, and ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}{\mathtt{M}}{\mathtt{2}}$, respectively, and are separately presented in Section 3.4.

All subsonic simulations reach a stationary regime with a sonic Mach number ${M}_{s}\approx 0.5$ after ≈2.5T, as illustrated in Figure 1(a). In contrast to this, the Alfvénic Mach number ${M}_{a}=\sqrt{\rho }u/B$ varies substantially between different normalizations. While ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ is ≈3 in the stationary regime, it reaches only ≈2 for all other simulations, shown in Figure 1(b). Given that ${M}_{a}$ is a proxy for the ratio of kinetic to magnetic energy, a higher value for ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ corresponds to a lower saturation value of the magnetic energy.

Zoom In Zoom Out Reset image size

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 ${\boldsymbol{u}}\,=\,{{\boldsymbol{u}}}_{{\rm{c}}}+{{\boldsymbol{u}}}_{{\rm{s}}}+{{\boldsymbol{u}}}_{0}$ with ${\rm{\nabla }}\cdot {{\boldsymbol{u}}}_{{\rm{s}}}=0$ and ${\rm{\nabla }}\,\times {{\boldsymbol{u}}}_{{\rm{c}}}={\boldsymbol{0}}$. The simulation ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ has significantly less power (≈50% 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 6 ≲ k ≲ 20. It is steeper, as is expected from a strongly compressive turbulence phenomenology (Burgers 1948).

Zoom In Zoom Out Reset image size

The rotational and total energy spectra for the simulations ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{16}}{\mathtt{T}}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{4}}{\mathtt{T}}}$, and ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$ 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 $\langle a\rangle$), there is more power in the compressive modes. This can be explained by a more detailed analysis of the link between $\langle a\rangle$ and ${T}_{\mathrm{corr}}$ in the following subsection.

In all simulations, the driving field is purely solenoidal; i.e., it carries no compressive power itself. Nevertheless, a strong (high $\langle a\rangle$) driving is expected to seed compressive modes; compare this to the canonical idea of wave steepening.

Figure 3(a) shows $\langle a\rangle$ for all subsonic simulations. By construction $\langle a\rangle$ is constant for ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{16}}{\mathtt{T}}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{4}}{\mathtt{T}}}$, and ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$, while it varies strongly over more than two orders of magnitude for ${{\mathtt{F}}}_{\dot{E}}^{\delta }$, which is normalized for a constant energy injection rate $\dot{E}$. To keep $\dot{E}$ at a given value $\langle a\rangle$ 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 ${\boldsymbol{a}}$ and ${\boldsymbol{u}}$) would be negligible for "small" $\langle a\rangle$. Thus, during the normalization step $\langle a\rangle$ is increased (decreased) to match a desired $\dot{E}$ for a predominantly perpendicular (aligned) combination of flow configuration and acceleration field.

Zoom In Zoom Out Reset image size

The same mechanism, i.e., the alignment of ${\boldsymbol{a}}$ and ${\boldsymbol{u}}$, is responsible for requiring a higher $\langle a\rangle$ value for smaller ${T}_{\mathrm{corr}}$ to reach the same ${{M}}_{s}$ in the other simulations. Figure 3(b) illustrates the mean PDF of the angle between ${\boldsymbol{a}}$ and ${\boldsymbol{u}}$. A larger correlation time results in a distribution for which there is a tendency for ${\boldsymbol{a}}$ and ${\boldsymbol{u}}$ to be 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 ${\boldsymbol{a}}$ and ${\boldsymbol{u}}$) and require larger $\langle a\rangle$ for lower ${T}_{\mathrm{corr}}$ (Eswaran & Pope 1988). For the extreme cases, ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ and ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$, 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 ${T}_{\mathrm{corr}}$, another feature is introduced to the flow by larger $\langle a\rangle$: compressive modes. At locations where ${\boldsymbol{u}}$ and ${\boldsymbol{a}}$ are aligned, which is always the case given the non-zero PDF around $\cos (\sphericalangle ({\boldsymbol{u}},{\boldsymbol{a}}))=1$ in Figure 3(b), large $\langle a\rangle$ leads to a strong acceleration, resulting in immediate downstream compression. A careful examination of the compressive power spectra in Figure 2 reveals signatures of this effect. The compressive power of ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ peaks at k ≈ 3, 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.

A link between the Pearson correlation coefficient between density and magnetic field strength, $\mathrm{Corr}[\rho ,B]$, 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. Figure 4(a) shows the correlation coefficient over time for all subsonic simulations. ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{4}}{\mathtt{T}}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{16}}{\mathtt{T}}}$, and ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$ build one family for which the correlation coefficient in the stationary regime increases from −0.803(12), −0.70(3), −0.56(3), to −0.50(4), respectively. In other words, there exists a strong anticorrelation for a ${T}_{\mathrm{corr}}=1{\text{}}T$ that decreases with smaller ${T}_{\mathrm{corr}}$ toward to a still-significant non-zero value of −0.5 for δ-in-time forcing (${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$). In contrast, ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ exhibits virtually no correlation between the density field and magnetic field strength $\mathrm{Corr}[\rho ,B]=-0.064(18)$. This is in agreement with the results of Yoon et al. (2016), who analyzed two forcing configurations corresponding to our ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$ and ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ cases. Given our additional simulations with varying ${T}_{\mathrm{corr}}$, we argue that the power in the acceleration field $\langle a\rangle$ and not ${T}_{\mathrm{corr}}$ is the primary driver of changes in $\mathrm{Corr}[\rho ,B]$. Otherwise, ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$ and ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ should yield identical results, which is not observed here. Nevertheless, $\langle a\rangle$ and ${T}_{\mathrm{corr}}$ are tightly linked, as shown in the previous subsection.

Zoom In Zoom Out Reset image size

Clear differences between the simulations are also observed in the logarithmic density PDFs, as illustrated in Figure 4(b). ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$ 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 ${T}_{\mathrm{corr}}$, it it still present in the ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$ simulation. Again, ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ contrasts with these results, showing an almost symmetrical distribution. In general, the high-density tails (with ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ allowing for the most and ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$ for the least extreme value) nicely illustrate how larger $\langle a\rangle$ lead to immediate compression.

Finally, in order to illustrate how these differences translate to biases in interpretations of astrophysical observations, we calculate LOS magnetic field strength derived from rotation measures⁷ as

$$\begin{eqnarray}&&\mathrm{LOS}{B}_{i}=\displaystyle \frac{{\int }_{L}\rho (l){B}_{\parallel }(l){\rm{d}}l}{{\int }_{L}\rho (l){\rm{d}}l}\end{eqnarray} \tag{ 4 }$$

with B_∥ being the LOS component of the magnetic field. Beck et al. (2003) derived how an (anti)correlation between ρ and B changes this measurement of the mean magnetic field strength, which is exact for uncorrelated fields. For anticorrelated fields Equation (4) underestimates the true LOS B. This effect can be observed in Figure 4(c), where the mean PDFs of the LOS B_x (over the 512² available LOSs) are shown. With increasing anticorrelation from ${{\mathtt{F}}}_{\dot{E}}^{\delta }$, ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{16}}{\mathtt{T}}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{4}}{\mathtt{T}}}$, to ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$ the derived values of 0.1640(8), 0.1589(12), 0.1584(12), 0.1558(8), and 0.1533(11), respectively, deviate further from the real value 0.1667. More strikingly, the PDF of ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ is much more peaked with a standard deviation of 0.058(3) compared to ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{16}}{\mathtt{T}}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}/{\mathtt{4}}{\mathtt{T}}}$, and ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$ with 0.103(5), 0.104(4), 0.125(12), and 0.116(5), respectively.

The supersonic simulations ${{\mathtt{F}}}_{\dot{E}}^{\delta }{\mathtt{M}}{\mathtt{2}}$, ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }{\mathtt{M}}{\mathtt{2}}$, and ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}{\mathtt{M}}{\mathtt{2}}$ exhibit a very similar behavior to their subsonic counterparts. All of the simulations reach a stationary regime after ≈2.5T with the same ${M}_{s}\approx 2.1$, but different ${M}_{a}$ of 2.84(13), 1.96(19), and 1.62(7), respectively. Again, a higher ${M}_{a}$ for similar ${M}_{s}$ in the saturated regime implies less-effective magnetic field amplification in the ${{\mathtt{F}}}_{\dot{E}}^{\delta }{\mathtt{M}}{\mathtt{2}}$ case. Similarly, the ratio of compressive versus rotational power in the kinetic energy spectrum decreases from ${{\mathtt{F}}}_{\dot{E}}^{\delta }{\mathtt{M}}{\mathtt{2}}$ to ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }{\mathtt{M}}{\mathtt{2}}$ to ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}{\mathtt{M}}{\mathtt{2}}$. 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 Figure 3(b).

Quantitative differences between both regimes are first observed in the ρ–B correlation, as illustrated in Figure 4(d). ${{\mathtt{F}}}_{\dot{E}}^{\delta }{\mathtt{M}}{\mathtt{2}}$ and ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }{\mathtt{M}}{\mathtt{2}}$ exhibit virtually no correlation with correlation coefficients of 0.04(5) and −0.09(7) whereas ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}{\mathtt{M}}{\mathtt{2}}$ still shows a weak anticorrelation with a coefficient of −0.23(3). The logarithmic density PDFs of the supersonic simulations in Figure 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 LOS magnetic field strength measurements are again systematically affected as shown in Figure 4(f). The derived value in the ${{\mathtt{F}}}_{\dot{E}}^{\delta }{\mathtt{M}}{\mathtt{2}}$ simulation of 0.630(12) is closest the real value of 2/3 and the PDF is most peaked with a standard deviation of 0.229(17) whereas ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }{\mathtt{M}}{\mathtt{2}}$ and ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}{\mathtt{M}}{\mathtt{2}}$ underestimate the magnetic field strength with 0.562(19) and 0.544(15), respectively, and generally broader PDFs with deviations of 0.39(2) and 0.46(4), respectively.

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 L_p) in nature evolves instantaneously.⁸ For a δ-in-time evolution with ${T}_{p}\to 0$, the characteristic velocity of that process is ${U}_{p}\to \infty$. This leads to the second, numerical argument. By construction, a δ-in-time forcing in a numerical simulation is not resolving the physical timescale. The timestep Δ_t 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 ${U}_{p}\to \infty$ the required timestep ${{\rm{\Delta }}}_{t}\to 0$. This restriction can never be satisfied. Therefore, a large-scale δ-in-time forcing is never numerically resolved.

Similar to the autocorrelation time discussion, choosing a normalization to a constant energy injection rate $\dot{E}$ over a constant rms acceleration $\langle a\rangle$ is appealing at first glance. From a turbulence analysis point of view $\dot{E}$ 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 $\dot{E}$ and $\langle a\rangle$ is less important. Normalizing by $\langle a\rangle$ for the stochastic forcing used here naturally leads to a statistical constant energy injection rate if ${T}_{\mathrm{corr}}$ 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 $\langle a\rangle$ 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 $\dot{E}$ can be seen as a process that depends on the interaction with the environment, for example, cold-mode active galactic nuclei (AGN) accretion (Gaspari et al. 2012; Li et al. 2015; Meece et al. 2017).

The strong anticorrelation between the density and the magnetic field strength observed in the ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$ run is consistent with the expectation of a (statistical) total pressure equilibrium (Beck et al. 2003)

$$\begin{eqnarray}&&{B}^{2}/2+{p}_{\mathrm{th}}={p}_{\mathrm{tot}}\approx \mathrm{const}.\end{eqnarray} \tag{ 5 }$$

In addition, the isothermal equation of state used in the simulations mandates ${p}_{\mathrm{th}}\propto \rho$. Hence, in order 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 ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ equivalent simulation by the forcing timescale. δ-in-time forcing evolves so fast that the system is not able to reach pressure equilibrium. Thus, ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$ 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 $\langle a\rangle$ (and, thus, smaller ${T}_{\mathrm{corr}}$) is the main driver behind a decreasing ρ–B anticorrelation. Given Alfvén's theorem in MHD, compression naturally leads to a positive correlation between ρ and B. 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 ρ–B anticorrelation is successively weakened by increasing compressive modes associated with a positive ρ–B correlation. This is also in agreement with the supersonic simulations, which naturally have more power in compressive modes.

The disparity of the ρ–B correlations for the different forcing parameters is less pronounced with an increasing sonic Mach number. We expect that it becomes negligible in the hypersonic (${M}_{s}\gtrsim 5$) regime, as, for example, found in molecular clouds, for example, because compressive modes become dynamically important independent of the forcing scheme.

On the other hand, there is no indication that the density and LOS B PDFs of ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ and ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$ or ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$ generally converge with an increasing Mach number. While we expect that the autocorrelation time (and, thus, the power in the acceleration field) becomes less important for the statistics if the acceleration field is normalized to $\langle a\rangle$ (i.e., ${{\mathtt{F}}}_{\langle a\rangle }^{\delta }$ and ${{\mathtt{F}}}_{\langle a\rangle }^{{\mathtt{1}}{\mathtt{T}}}$), the bias for ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ is likely to remain. We base this expectation on the increasing extreme variations in $\langle a\rangle$ in comparison to the velocity dispersion when the acceleration field is normalized to $\dot{E}$.

Overall, this suggests that in the hypersonic regime a δ-in-time forcing with constant power (while still being unrealistic) can probably be used without major implications on turbulence statistics. Nevertheless, a more detailed study is required to verify this statement.

An empirical relation between the derived LOS 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 ${{\mathtt{F}}}_{\dot{E}}^{\delta }$ forcing and assumes no ρ–B correlation, as reported by Wu et al. (2015). While there are differences in our $\langle a\rangle$-normalized simulations, the differences are overall much less pronounced compared to what is observed in ${{\mathtt{F}}}_{\dot{E}}^{\delta }$. Thus, a modified LOS B 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 ρ–B anticorrelation. Their estimate assumes underlying isotropic Alfvénic perturbations of the flow. However, the subsonic, super-Alfvénic regime that we are probing is potentially dominated by slow modes, which is inferred from the observed ρ–B anticorrelation (Passot & Vázquez-Semadeni 2003). This suggests a review of the Davis–Chandrasekhar–Fermi method for different regimes.

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 of 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 ${M}_{a}$. While Yoon et al. (2016) conducted and analyzed simulations with varying ${M}_{a}$, 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 such as the link between $\langle a\rangle$ and compressive modes, or the increasing alignment of velocity and acceleration field with increasing ${T}_{\mathrm{corr}}$, also in the pure hydrodynamic case.

Finally, we conducted all of the simulations at a resolution of 512³, which is nowadays commonly used for turbulent boxes. While a higher resolution (and, thus, a larger dynamical range) is preferable, we observe the same effects in higher-resolution simulations at individual points in the probed parameter space (Grete et al. 2017). Thus, there is no indication that the described processes are affected by the resolution of the simulation.

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.