Compression of turbulent magnetized gas in Giant Molecular Clouds

TL;DR

This study uses 3D MHD simulations to explore how compression in magnetized, turbulent gas affects energy dissipation and velocity dispersion in molecular clouds, revealing reduced dissipation and higher turbulence levels.

Contribution

It introduces a novel simulation-based analysis of compressing, magnetized turbulence and develops an analytic model for the effective equation of state in such conditions.

Findings

Reduced dissipation rate in compressing, magnetized turbulence.

Higher equilibrium velocity dispersion compared to non-magnetized or non-compressing cases.

Implications for turbulence driving and chemical processes in molecular clouds.

Abstract

Interstellar gas clouds are often both highly magnetized and supersonically turbulent, with velocity dispersions set by a competition between driving and dissipation. This balance has been studied extensively in the context of gases with constant mean density. However, many astrophysical systems are contracting under the influence of external pressure or gravity, and the balance between driving and dissipation in a contracting, magnetized medium has yet to be studied. In this paper we present three-dimensional (3D) magnetohydrodynamic (MHD) simulations of compression in a turbulent, magnetized medium that resembles the physical conditions inside molecular clouds. We find that in some circumstances the combination of compression and magnetic fields leads to a rate of turbulent dissipation far less than that observed in non-magnetized gas, or in non-compressing magnetized gas. As a result, a compressing, magnetized gas reaches an equilibrium velocity dispersion much greater than would be expected for either the hydrodynamic or the non-compressing case. We use the simulation results to construct an analytic model that gives an effective equation of state for a coarse-grained parcel of the gas, in the form of an ideal equation of state with a polytropic index that depends on the dissipation and energy transfer rates between the magnetic and turbulent components. We argue that the reduced dissipation rate and larger equilibrium velocity dispersion produced by compressing, magnetized turbulence has important implications for the driving and maintenance of turbulence in molecular clouds, and for the rates of chemical and radiative processes that are sensitive to shocks and dissipation.

Tables

Table 1. Table 1: List of simulation parameters at the beginning of the contraction phase.

Simulation model	$⟨B_z⟩$	$B_{\mathrm{turb}}$	$H$	$ℳ$	${\sigma }_{\mathrm{sat}}(\%)$	$N_{\mathrm{res}}^3$
HD	$0$	$0$	$-200$	$9.4$	$0$	${512}^3$
HD-H1	$0$	$0$	$-1$	$9.4$	$0$	${512}^3$
noGF-Medium	$0$	$9.2$	$-200$	$9.8$	$8$	${512}^3$
noGF-Strong	$0$	$24$	$-200$	$9.1$	$56$	${512}^3$
noGF-Strong-LR	$0$	$21$	$-200$	$9.5$	$40$	${256}^3$
noGF-Strong-H1	$0$	$24$	$-1$	$9.1$	$56$	${512}^3$
GF-Weak	$0.35$	$5.5$	$-200$	$9.3$	$3$	${512}^3$
GF-Medium	$3.5$	$18$	$-200$	$9.5$	$30$	${512}^3$
GF-Medium-H1	$3.5$	$18$	$-1$	$9.5$	$30$	${512}^3$
GF-Strong	$35$	$21$	$-200$	$9.8$	$140$	${512}^3$

Equations

$$\frac{\partial }{\partial t}\rho +\nabla \cdot \left(\rho \mathbf{v}\right)=0,$$

$$\frac{\partial }{\partial t}\left(\rho \mathbf{v}\right)+\nabla \cdot \left(\rho \mathbf{v}\otimes \mathbf{v}-\frac{1}{4\pi }\mathbf{B}\otimes \mathbf{B}\right)+\nabla P_{\mathrm{tot}}=0,$$

$$\frac{\partial }{\partial t}e+\nabla \cdot \left[\left(e+P_{\mathrm{tot}}\right)\mathbf{v}-\frac{1}{4\pi }\left(\mathbf{B}\cdot \mathbf{v}\right)\mathbf{B}\right]=0,$$

$$\frac{\partial }{\partial t}\mathbf{B}-\nabla \times \left(\mathbf{v}\times \mathbf{B}\right)=0,\nabla \cdot \mathbf{B}=0.$$

$$\rho$$

$$\mathbf{B}$$

$$P_{\mathrm{tot}}$$

$$e$$

$${ϵ}_{\mathrm{int}}$$

$$\frac{\partial }{\partial t}\rho +\nabla \cdot \left(\rho \mathbf{v}\right)=0,$$

$$\displaystyle\!\begin{multlined}\frac{\partial{\!}}{\partial t}\left(\rho\mathbf{v}\right)+\nabla\cdot\left(\rho\mathbf{v}\!\otimes\!\mathbf{v}-\frac{1}{4\pi}\mathbf{B}\!\otimes\!\mathbf{B}\right)+\nabla P_{\mathrm{tot}}={}\\ -2H\rho\mathbf{v},\end{multlined}\frac{\partial{\!}}{\partial t}\left(\rho\mathbf{v}\right)+\nabla\cdot\left(\rho\mathbf{v}\!\otimes\!\mathbf{v}-\frac{1}{4\pi}\mathbf{B}\!\otimes\!\mathbf{B}\right)+\nabla P_{\mathrm{tot}}={}\\ -2H\rho\mathbf{v},$$

$$\displaystyle\!\begin{multlined}\frac{\partial{\,}}{\partial t}e+\nabla\cdot\left[\left(e+P_{\mathrm{tot}}\right)\mathbf{v}-\frac{1}{4\pi}\left(\mathbf{B}\cdot\mathbf{v}\right)\mathbf{B}\right]={}\\ \phantom{..........................}-H\left[(3\gamma-1)\rho\epsilon_{\mathrm{int}}+2\rho\mathbf{v}\cdot\mathbf{v}\right],\end{multlined}\frac{\partial{\,}}{\partial t}e+\nabla\cdot\left[\left(e+P_{\mathrm{tot}}\right)\mathbf{v}-\frac{1}{4\pi}\left(\mathbf{B}\cdot\mathbf{v}\right)\mathbf{B}\right]={}\\ \phantom{..........................}-H\left[(3\gamma-1)\rho\epsilon_{\mathrm{int}}+2\rho\mathbf{v}\cdot\mathbf{v}\right],$$

$$\frac{\partial }{\partial t}\mathbf{B}-\nabla \times \left(\mathbf{v}\times \mathbf{B}\right)=-\frac{3}{2}H\mathbf{B},\nabla \cdot \mathbf{B}=0.$$

$${\eta }^{\ast }\frac{v_{\mathrm{rms}}(a)}{a\lambda }=\eta \frac{v_{\mathrm{rms}}}{aL}=|H|=-\frac{\dot{a}}{a},$$

$$P_{\mathrm{kin}}$$

$$P_{\mathrm{B}}$$

$$\frac{P}{P_{\mathrm{th}}}\propto {\rho }^{\gamma -1}\propto a^{-3(\gamma -1)}.$$

$$\overline{\varphi }$$

$${\varphi }^{\prime }$$

$$\tilde{\varphi }$$

$$0$$

$${\mathbf{\tau }}_{\mathrm{R},\mathrm{SFS}}$$

$${\mathbf{\tau }}_{\mathrm{M},\mathrm{SFS}}$$

$$P_{\mathrm{kin}}$$

$${\mathbf{\tau }}_{\mathrm{R},\mathrm{SFS}}$$

$$P_{\mathrm{B}}$$

$${\mathbf{\tau }}_{\mathrm{M},\mathrm{SFS}}$$

$$0$$

$$e_{\mathrm{kin},\mathrm{SFS}}$$

$$e_{\mathrm{B},\mathrm{SFS}}$$

$$P_{\mathrm{kin}}$$

$$P_{\mathrm{B}}$$

$$P=\left(\gamma -1\right)\rho ϵ,$$

$$P_{\mathrm{tot}}$$

$$=$$

$$\gamma ={\left(\frac{\partial \ln P}{\partial \ln \rho }\right)}_s,$$

$${\gamma }_{\mathrm{eff}}=\frac{d\ln P_{\mathrm{tot}}}{d\ln \rho }=\frac{\rho }{P_{\mathrm{tot}}}\frac{{\dot{P}}_{\mathrm{tot}}}{\dot{\rho }},$$

$$\dot{P}=\left(\gamma -1\right)\left(\dot{ϵ}\rho +ϵ\dot{\rho }\right),$$

$$\dot{ϵ}=-P\dot{V}-q,$$

Full text

††pagerange: Compression of turbulent magnetized gas in Giant Molecular Clouds–B††pubyear: 2017

1 Introduction

Magnetized plasma is ubiquitous in astrophysical systems. Particularly, gas in the interstellar medium (ISM) is observed to be magnetized, and a large fraction of its energy content is in the form of magnetic fields (see Beck et al., 1996; Ferrière, 2001; Federrath, 2016, and references within). The magnetic field in the ISM of disk galaxies consists of an ordered rotating component on galactic disk scales that is consistent with slow winding of the magnetic field via macroscopic dynamo processes, and small-scale magnetic fields that are generated by the winding of the magnetic field via turbulent dynamo processes (Brandenburg & Subramanian, 2005; Brandenburg et al., 2012). For the Galaxy, the values of the two components are comparable, with typical values of $2$–$5\,\mu\mathrm{G}$.

When portions of this magnetized fluid are subject to rapid radiative cooling, for example in molecular clouds, the result is a highly supersonic, strongly magnetized flow. Within such a flow, the velocity dispersion is dictated by the balance between driving and dissipation processes. This balance, particularly the dissipation part of it, has been studied extensively for both non-magnetized and magnetized flows in the context of periodic boxes with constant mean density (e.g. Mac Low et al., 1998; Stone et al., 1998; Ostriker et al., 1999; Padoan & Nordlund, 1999; Mac Low, 1999; Kritsuk et al., 2007; Lemaster & Stone, 2009). The general result from these simulations is that the turbulence decays on a timescale comparable to a large eddy turnaround time, and that the rates of decay are not substantially altered by the presence or absence of a magnetic field.

The problem of the balance between driving and decay for magnetized turbulence is most acute in molecular clouds. Since these have linewidths indicating the presence of supersonic flow, the fast dissipation of turbulence found by these simulations necessitates a mechanism to reinject the energy equally quickly. A number of candidates have been proposed, including internal feedback from H ii regions (Matzner, 2002; Krumholz et al., 2006; Goldbaum et al., 2011) or protostellar outflows (Li & Nakamura, 2006; Nakamura & Li, 2007; Wang et al., 2010; Federrath et al., 2014a), driving of turbulence by ongoing accretion (Klessen & Hennebelle, 2010; Goldbaum et al., 2011; Lee & Hennebelle, 2016) or gravitational contraction on small scales (Federrath et al., 2011b; Sur et al., 2012), thermal instability driving (Koyama & Inutsuka, 2002; Hennebelle & Inutsuka, 2006) and injection of energy from external supernova shocks (Mac Low & Klessen, 2004; Padoan et al., 2016b, a; Pan et al., 2016). Alternately, it is possible that the linewidths do not reflect turbulent motion at all, and instead indicate global gravitational collapse (Ballesteros-Paredes et al., 2011; Zamora-Avilés & Vázquez-Semadeni, 2014). Each of these proposals, however, faces challenges – internal feedback must maintain large linewidths without destroying the clouds in which they occur, driving by accretion faces the problem of what happens when the accretion eventually ends, thermal instability seems unlikely to be a viable mechanism in molecule-dominated galaxies that lack a significant warm phase, and external driving requires efficient coupling between the low density external medium and the dense clouds. The view that clouds are in global collapse is hard to reconcile with the observed very low rates of star formation found even in gas at densities $\gtrsim 10^{5}$ cm*-3* (Krumholz & Tan, 2007; Krumholz et al., 2012; Federrath & Klessen, 2012; Padoan et al., 2014; Evans et al., 2014; Usero et al., 2015; Salim et al., 2015; Vutisalchavakul et al., 2016; Heyer et al., 2016).

The problem of the persistence of turbulence in molecular clouds is significantly eased if gravitational compression is able to pump energy into turbulent motion, since this would provide a mechanism to both power the turbulence and slow the collapse. The phenomenon has been explored for non-magnetized flows by Robertson & Goldreich (2012). In their work, an initially-turbulent gas is compressed in a scale-free manner by renormalizing the thermodynamic variables according to the expected values from a uniform collapse. As gas is compressed, the amplitude of the velocity field increases because the compression does $P\,dV$ work against the kinetic pressure. On the other hand, the typical size of the eddies is reduced by the compression, and this accelerates the decay of turbulence. Depending on the compression rate, one process or the other dominates, and the turbulence either increases or decays. Qualitatively, the results are consistent with what one would have derived by naively equating the rate of $P\,dV$ work with a decay timescale of $\sim 1$ eddy turnover time derived from non-compressing driven turbulence simulations: the turbulence is amplified when the box compression time is short compared to the eddy turnover time, and decays if the converse holds. However, Robertson & Goldreich (2012) did not include magnetic fields in their simulations, and we know that all clouds in the ISM are magnetized to a level that corresponds to a near equipartition between turbulent and magnetic energy densities (Crutcher, 2012, and references therein).

In this work we seek to determine whether the result by Robertson & Goldreich (2012) is altered in the presence of a magnetic field. We have already noted that, in driven turbulence simulations, magnetic fields make no qualitative difference. However, driving turbulence by global compression is qualitatively different than direct driving of the gas. In the first case, the scaling relations of velocity and distance for global compression enhance all modes similarly, and some modes, for which the dissipation is faster and that are not replenished quickly enough by a turbulent cascade, can decay and disappear. In the second, the turbulence forcing arbitrarily sets the geometry of the flow and phases of the various modes, preventing the flow from achieving a more relaxed state.

A magnetic field might change the situation for a compressing flow in two ways. First, a magnetic field and gas motions can exchange energy via a turbulent dynamo (Kazantsev, 1968; Subramanian, 1997, 1999; Brandenburg & Subramanian, 2005; Schleicher et al., 2010; Schober et al., 2012b; Schober et al., 2012c; Schober et al., 2012a; Bovino et al., 2013; Schober et al., 2015). For driven turbulence, the amount of energy stored in the dynamo is limited by the back reaction of the Lorentz forces on the gas, (Federrath et al., 2011a; Federrath, 2016), and as a result the energy stored in the magnetic field is always subdominant compared to the turbulence. However, gravitational compression will amplify magnetic fields differently than gas motions (Sur et al., 2010; Federrath et al., 2011b; Sur et al., 2012), potentially leading to magnetic-turbulent interactions not found in driven, non-compressing boxes. Second, magnetic fields will impose anisotropy on the flow, and anisotropic turbulence shows a different cascade pattern and a different decay rate than isotropic turbulence (Cho et al., 2002; Cho & Lazarian, 2003; Hansen et al., 2011).

In this paper we examine the effect of magnetic fields on global compression of turbulent gas in idealized 3D MHD simulations. We distinguish between cases of zero net magnetic flux and cases with non-zero net flux of various amplitudes. This is of theoretical importance because a magnetic field with finite flux increases monotonically as gas contracts, and of practical interest because fields with non-zero net flux are likely present in proto-GMCs (Li & Henning, 2011; Li et al., 2011; Pillai et al., 2015). Rather than introducing cooling, we assume that the gas is isothermal, which is a reasonable approximation for GMCs over a wide range of densities. We describe our setup in §2 and our simulated results in §3. We then construct an analytic prediction for the effective equation of state of a system with mixed thermal, kinetic and magnetic pressure components (§4), use some of the physical insights to further analyze the dissipation in the simulations and compare our predictions to the simulations (§4.4). In §5 we discuss possible implication of our results to ISM and GMCs, and in §6 we summarize and conclude.

2 Simulations

2.1 The FLASH code

We use a modified version of the grid-based code FLASH (Fryxell et al., 2000; Dubey et al., 2008) (http://www.flash.uchicago.edu/site/flashcode/) to solve the three-dimensional (3D), compressible, ideal magnetohydrodynamical (MHD) equations,

[TABLE]

Here, $\rho$, $\mathbf{v}$, $P_{\mathrm{tot}}=P_{\mathrm{th}}+(1/8\pi)\left|\mathbf{B}\right|^{2}$, $\mathbf{B}$, and $e=\rho\epsilon_{\mathrm{int}}+(1/2)\rho\left|\mathbf{v}\right|^{2}+(1/8\pi)\left|\mathbf{B}\right|^{2}$ denote the gas density, velocity, pressure (thermal plus magnetic), magnetic field, and total energy density (internal, plus kinetic, plus magnetic), respectively. The MHD equations are closed with a quasi-isothermal equation of state (EoS), $P_{\mathrm{th}}=(\gamma-1)\rho\epsilon_{\mathrm{int}}$, where we set $\gamma=1.00001$. Using this setting we model a gas with an extremely high number of degrees of freedom, $f=2/(\gamma-1)\sim 2\times 10^{5}$, effectively resulting in a gas that is isothermal. This is a standard procedure to obtain a quasi-isothermal EoS and results in the same thermodynamic response of the gas as a polytropic EoS, $P_{\mathrm{th}}\propto\rho^{\Gamma}$ with $\Gamma=1$ (Federrath et al., 2014a). The practical reason for choosing the simple ideal gas EoS with $\gamma=1.00001$ is to keep track of how much energy is dissipated by the turbulence, i.e., we solve the energy equation (3) and record the change in energy every timestep.

The system of ideal MHD equations (1–4) are solved with the robust HLL3R Riemann scheme by Waagan et al. (2011a), based on previous developments in applied mathematics to preserve positive density and pressure by construction (Bouchut et al., 2007; Klingenberg et al., 2007; Waagan, 2009; Bouchut et al., 2010). For our particular simulations, the magnetic field is shown to remain divergence free to a reasonable degree (Appendix §A.1)

2.2 Numerical scheme for solving the MHD equations in an expanding or contracting coordinate system

The cosmology unit in FLASH allows one to use any hydrodynamics solver written for a non-expanding universe to work unmodified in a cosmological context. This is achieved by solving the MHD equations in the co-moving reference frame and accounting for the additional terms in the MHD equations that appear due to the expansion or contraction of the system. All calculations are assumed to take place in co-moving coordinates $\mathbf{x}$ = $\mathbf{r}/a$, where $\mathbf{r}$ is the physical (proper) position vector and $a(t)$ is the time-dependent cosmological scale factor. When transforming the MHD equations to the co-moving frame, the spatial derivative transforms as $\nabla_{\mathbf{x}}=a\nabla_{\mathbf{r}}$ and the time derivative transforms as $(\partial/\partial t)_{\mathbf{x}}=(\partial/\partial t)_{\mathbf{r}}+H\mathbf{r}\cdot\nabla_{\mathbf{r}}$, where the Hubble constant is defined as $H=\dot{a}/a$. The physical (proper) velocity is given as $\tilde{\mathbf{v}}=H\mathbf{r}+a\dot{\mathbf{x}}$, where the first term is the Hubble flow and the second term contains the co-moving velocity $\mathbf{v}=\dot{\mathbf{x}}$.

Using these relations between the physical and co-moving derivatives in addition to the following transformations from physical (with tilde) to co-moving hydrodynamical quantities (without tilde),

[TABLE]

the MHD equations in co-moving coordinates have exactly the same form as Equations (1–4) with additional Hubble source terms on the right-hand sides of the momentum, energy and induction equations:

[TABLE]

Note that we have changed all time and space derivatives in these equations to the co-moving frame, i.e., $\partial/\partial t\equiv(\partial/\partial t)_{\mathbf{x}}$ and $\nabla\equiv\nabla_{\mathbf{x}}$.

Since the form of these equations is identical to the conservation Equations (1–4) without the Hubble source terms, we can use any existing hydrodynamical scheme to solve this set of equations in the co-moving frame. In order to account for the Hubble source terms on the right-hand side of these equations, we use an operator-splitting approach, where the co-moving hydrodynamical variables are modified in each time step (after the hydro step) to account for the source terms.

First we note that the mass continuity equation is unchanged between physical and co-moving coordinates. The momentum equation has the Hubble source term $-2H\rho\mathbf{v}$. Expanding the co-moving momentum equation (13) with respect to the change in $a$, we find $\cancel{\dot{\rho}v}+\rho\dot{v}+\cancel{\nabla(\dots)}=-2H\rho v$, where $\dot{\rho}=(d\rho/da)(da/dt)=0$ because $d\rho/da=0$, and any spatial derivatives cancel, because $a$ does not depend on space. This leaves us with the simple differential equation, $\dot{v}/v=-2\dot{a}/a$, for which the solution is $v^{\prime}=v(a/a^{\prime})^{2}$, where $v$ and $a$ are the velocity and scale factor before accounting for the Hubble term (i.e., before the hydro step) and $v^{\prime}$ and $a^{\prime}=a(t+\Delta t)$ are the velocity and scale factor after the current time step $\Delta t$. An analogous correction has to be made in the co-moving energy equation to account for the Hubble source term, i.e., $\epsilon_{\mathrm{int}}^{\prime}=\epsilon_{\mathrm{int}}(a/a^{\prime})^{3\gamma-1}$.

These procedures to account for the Hubble flow in pure hydrodynamics (without magnetic fields) were already implemented in the cosmology module of the public version of FLASH. However, MHD was not supported. Here we implemented the necessary modifications of the induction equation with the Hubble source term $-(3/2)H\mathbf{B}$ in Equation (17), which requires a modification of the co-moving magnetic field with $B^{\prime}=B(a/a^{\prime})^{3/2}$, analogous to the operator-split corrections for the velocity and energy explained in the previous paragraph.

2.3 Initial driving of turbulence

In order to establish a fully-developed turbulent state, we first drive turbulence for a few crossing times. The state after this initial driving phase serves as the initial condition for our numerical experiments on the statistics of MHD turbulence in a contracting reference frame. Since we are focussing on MHD turbulence in molecular clouds, we drive turbulence to a target mass-weighted (MW) Mach number $\mathcal{M}=\langle v_{\mathrm{rms}}/c_{\mathrm{s}}\rangle_{\mathrm{MW}}=9$–$10$ (i.e., supersonic turbulence; see Larson, 1981; Solomon et al., 1987; Ossenkopf & Mac Low, 2002; Heyer & Brunt, 2004; Heyer et al., 2009; Roman-Duval et al., 2011; Schneider et al., 2013a) by applying a driving field $\rho{\bf F}$ as a source term in the momentum equation (2). The sound speed is chosen as $c_{\mathrm{s}}=1$ in normalised units.

The turbulence driving field is constructed with a stochastic Ornstein-Uhlenbeck (OU) process (Eswaran & Pope, 1988; Schmidt et al., 2009; Price & Federrath, 2010), implemented by Federrath et al. (2010) and available in the public version of the FLASH code. The OU process creates a spatial and temporal driving pattern that varies smoothly in space and time with an auto-correlation timescale equal to the turbulent turnover time (also called turbulent box-crossing time), $t_{\mathrm{turb}}=L/(2\mathcal{M}c_{\mathrm{s}})=0.05$ for $\mathcal{M}=10$ on the largest scales ($L/2$) in our periodic simulation domain of side length $L=1$ (normalised units). The driving field ${\bf F}$ is constructed in Fourier space such that most power is injected at the smallest wave numbers, $1<\left|\mathbf{k}\right|L/2\pi<3$. The peak of energy injection is on scale $L/2$, i.e., $k=2$, and falls off as a parabola towards smaller and higher wave numbers, such that the driving power is identically zero at $k=1$ and $k=3$, as in our previous studies of driven turbulence (e.g., Federrath et al., 2010; Federrath, 2013, 2016). This procedure confines the effect of the driving to a narrow wave number range and allows the turbulence to develop self-consistently on smaller scales ($k\geq 3$).

In constructing the driving field, we apply a Helmholtz decomposition in Fourier space, in order to separate the driving field into its solenoidal and compressive parts. This allows us to construct a solenoidal (divergence-free) driving field ($\nabla\cdot\mathbf{F}=0$) or a compressive (curl-free) driving field ($\nabla\times\mathbf{F}=0$). The influence of different driving on the statistics of turbulence, the amplification of magnetic fields, and on the star formation rate has been determined in Federrath et al. (2008, 2009); Federrath et al. (2010); Federrath et al. (2011a), Federrath & Klessen (2012, 2013), Federrath (2013, 2016), and Federrath et al. (2016a, 2017). For simplicity and since here we simply want to seed a fully-developed initial turbulent state before starting the contraction, we chose to use purely solenoidal (divergence-free) driving.

2.4 Initial conditions and list of simulations

We start from gas with uniform density $\rho_{0}=1$ (normalised units) at rest and drive turbulence for $t_{0}=4\,t_{\mathrm{turb}}=0.2$ in a fixed (non-contracting) reference frame ($a=1$), which establishes fully-developed turbulence. After this, we begin the contraction phase, $a(t)<1$, at which point the driving is deactivated and turbulence as well as magnetic-field dynamics are solely determined by the contraction of the gas in the co-moving reference frame given by Equations (10–17).

Table 1 provides a list of all the simulations performed. We distinguish between three main cases: two purely hydrodynamical (HD) runs, four MHD runs without magnetic guide field (noGF), and four MHD runs that include a constant guide field (GF) $\langle B_{z}\rangle$ in the $z$-direction of the simulation domain; for each of the MHD cases we consider multiple field strengths in order to determine the sensitivity of the results to this parameter. The simulations without guide field use an initial turbulent field generated with a flat power spectrum in the range $k/(2\pi)=2$–$20$, which produces initial turbulent fields (after the driving phase) of $B_{\mathrm{turb}}=9.2$, $24$, and $21$, respectively (Table 1, middle section). The simulations with guide field were initialised with $\langle B_{z}\rangle=0.35$, $3.5$, and $35$, respectively, giving rise to initial turbulent fields (after the driving phase) of $B_{\mathrm{turb}}=5.5$, $18$, and $21$, respectively (Table 1, bottom section). We note that the field strength is in normalised units, so the Alfvén speed $v_{\mathrm{A}}=B/(4\pi\rho_{0})^{1/2}$ in normalised units. This means that $B\sim 3.5$ corresponds to an Alfvén speed of one and $B\sim 35$ to an Alfvén speed of 10, comparable to the turbulent velocity dispersion. Thus we can think of our three guide field cases as representing three regimes of plasma $\beta$ and Alfvén Mach number $\mathcal{M}_{A}$ (computed with respect to the guide field): GF-Weak has $\beta\gg 1$, $\mathcal{M}_{A}\ll 1$, GF-Medium has $\beta\sim 1$, $\mathcal{M}_{A}\ll 1$, and GF-Strong has $\beta\ll 1$, $\mathcal{M}_{A}\sim 1$.

All our simulations use the same resolution of $N_{\mathrm{res}}^{3}=512^{3}$ grid cells (except MHD-noGF-Strong-LR with $N_{\mathrm{res}}=256$, used to investigate numerical convergence in Appendix §A.2).

Finally, our simulations (Tab. 1) use the same time evolution for the scale factor $a(t)=\exp{[H(t-t_{0})]}$ for $t\geq t_{0}$ with $H=-(t_{\mathrm{turb}}/10)^{-1}=-200$, i.e., fast contraction on a time scale ten times shorter than the initial turbulent crossing time. However, we also run an HD-H1, MHD-noGF-Strong-H1 and MHD-GF-Medium-H1 simulation with $H=-1$, in order to demonstrate that our main conclusions do not depend on the choice of $H$ (see Appendix §B). Robertson & Goldreich (2012) discussed three different cases for the contraction law in the pure HD limit, while here we are primarily interested in the case of fast contraction (compression), focussing on the effect of the magnetic field (MHD runs).

3 Simulation Results

3.1 Evolution of the Mach number and energy

Fig. 1 and 2 show the time evolution of the rms Mach number and the kinetic and thermal energies, respectively, in all simulations. We show these quantities as a function of time during the initial driving phase, and as a function scale factor once compression begins, with the two phases separated by the solid vertical lines in the plots. Since the scale factor is exponential in time ($H=\dot{a}/a=-200$) and the plots use a logarithmic scale, position on the $x$-axis is proportional to time during both phases, albeit with different scalings. The top panels show our noGF simulations (without magnetic guide field) and the bottom panel shows our GF simulations. We show the pure hydrodynamic simulation (HD) in both panels to help guide the eye.

First examine Fig. 1. We see that, starting from the fully-developed turbulent state at $a=1$, all simulations start compression (decreasing $a$), which drives turbulence, i.e., increasing $v$ and $\mathcal{M}$. The turbulence is initially supersonic with $\mathcal{M}\sim 9$–$10$ (Tab. 1), and increases to a peak of $\mathcal{M}\approx 20$. However, at $a\sim 0.2$, in all simulations except noGF-Strong, the evolution reverses and turbulence begins to decay. This change is also apparent in the kinetic energy density evolution shown in Fig. 2, which increases sharply from $a=1$ to $a\approx 0.2$, but then shows an inflection point and increases less steeply thereafter.

This qualitative change from increasing to decaying turbulence is not related to the turbulence becoming sonic or subsonic, which only happens much later. Instead, it can be explained by the change of dissipation with $a$. The dissipation timescale is proportional to the largest eddy turnover time (Mac Low, 1999; Robertson & Goldreich, 2012) and the dissipation rate becomes comparable to the compression rate when

[TABLE]

with $v_{\mathrm{rms}}$ the root mean square of the velocity field, $\lambda=L/2$ the largest eddy size and $H=-200$, the compression rate in our simulations (§2). The coefficient $\eta^{*}=\eta/2$ is a dimensionless dissipation efficiency (see §4) and has been calibrated to be $\eta^{*}\approx 0.9$ (see §4.4). We can solve eq. 18 numerically for $a$ using the value of $v_{\rm rms}(a)$ measured from the simulations, and the result is $a\approx 0.2$; we show the exact solution for the HD run as the vertical dashed line in Fig. 1. It is evident that this typical timescale for equality between compression and dissipation successfully predicts the onset of efficient dissipation and roughly coincides with the beginning of the decaying stage of turbulence.

From $a=1$ to $a\approx 0.2$, the magnetic field has only minor effects on the evolution in all runs except GF-Strong. Compared to the HD case, in the noGF models the magnetic field stores additional energy which replenishes some of the kinetic energy that is dissipated. This slightly delays the onset of the decaying stage, and allows higher maximum velocities or Mach numbers by about $10$–$20\%$, but this is clearly a modest effect. However, at later times the MHD and HD runs show profound differences. In all the MHD runs the Mach number (Fig. 1) eventually stops decreasing and begins to increase again. Corresponding to this, the slope of the kinetic energy versus $a$ curve (Fig. 2) steepens again.

The value of $a$ at which the switch from decaying to increasing turbulence happens appears to depend both on whether there is a guide field, and on the saturation level of the turbulent dynamo, as parameterized by $\sigma_{\rm sat}\equiv e_{\rm B}/e_{\rm kin}$, at the onset of compression; we report this quantity in Tab. 1. The noGF-medium run has $\sigma_{\rm sat}=8\%$, and does not switch from decaying to increasing until $a\approx 0.001$, while the noGF-Strong ($\sigma_{\rm sat}=56\%$), GF-Weak ($\sigma_{\rm sat}=3\%$), and GF-Medium ($\sigma_{\rm sat}=30\%$) all reverse at $a\approx 0.01$. The GF-Strong case ($\sigma_{\rm sat}=130\%$) never goes through a decaying phase at all, and instead has a Mach number that increases almost monotonically. In call cases, the difference between the MHD and HD cases is large and growing with time. Even the noGF-Medium case, with $\sigma_{\rm sat}=8\%$, has $\sim 10$ times as much kinetic energy as the pure HD case by $a=0.001$. The GF-Strong case has 10 times the kinetic energy of the HD case even at $a\approx 0.1$, and by $a=0.01$ this gap has grown to more than two orders of magnitude.

3.2 Dissipationless flows

3.2.1 The transition to dissipationless flow

Having seen that the presence of a magnetic field causes a major change in the behavior of compressive turbulence, we now investigate in more detail the origin of this behavior. We shall show that this change the result of a shift in the flow pattern to one that is nearly dissipationless. As a first step in this direction, we note that the switch from decaying to increasing Mach number is associated the ratio of magnetic to kinetic pressure. In our dimensionless units, the volume-averaged thermal pressure is $1/V$, where $V$ is the box volume, and we define the volume-averaged kinetic and magnetic pressures by

[TABLE]

Note that the factor of $1/3$ in the definition of $P_{\rm B}$ might at first seem surprising, but we shall see the justification for it in §4.

We plot the time evolution of $P_{\rm kin}$ and $P_{\rm B}$ in all our runs in Fig. 3. As in Figs. 1 and 2, the $x$-axis is separated (by the vertical solid line) into the initial driving stage on the left, and the contraction phase on the right. At the beginning of compression, dissipation is comparatively unimportant because the compression timescale is small compared to the eddy turnover timescale. Thus the flow is nearly dissipationless. We show below that, for adiabatic contraction, kinetic pressure acts as a gas with $\gamma=5/3$ and turbulent magnetic pressure acts as a gas with $\gamma=4/3$, and we expect these pressures to scale as

[TABLE]

Thus we expect the kinetic to thermal ratio to scale as $a^{-2}$, and the magnetic to thermal ratio to scale as $a^{-1}$ for the case without a guide field. With a guide field, flux conservation requires that the mean magnetic field rise as $B_{\rm mean}\propto a^{-2}$, and thus the scaling is the same, though for a somewhat different reason. We show lines with slopes of $-1$ and $-2$ in Fig. 3, and they are indeed good descriptions of the slope at early times.

Unsurprisingly, the kinetic and magnetic pressures begin to drop when the dissipation rate becomes comparable to the compression rate. However, the kinetic term drops more steeply than the magnetic term. This effect is due to the fact that the only true dissipation channel in the system is via the kinetic term, and that the dissipation of the magnetic component is bottlenecked by the rate at which the now over-magnetized gas can transfer energy back into the kinetic component.

As a result of the difference in the dissipation rates, the magnetic pressure ultimately exceeds the kinetic pressure in all cases except GF-Weak. In all the other cases, the transition from decreasing to increasing Mach number occurs almost exactly when this crossover happens, although if one closely compares GF-Medium to noGF-Medium, it is clear that at equal field strength the transition occurs earlier, in terms of both $a$ and in terms of ratio of $P_{\rm kin}$ to $P_{\rm B}$, in the presence of a guide field. Whether the flow is subsonic or supersonic appears to make little difference to the transition, consistent with the findings of Mac Low (1999) that the dissipation rate is not greatly affected by whether the flow is subsonic or supersonic.

3.2.2 The nature of the dissipationless flow

When the magnetic pressure begins to dominate, or even earlier in the presence of a net magnetic flux, the flow re-arranges itself into a fundamentally different topology, characterized by a much lower rate of dissipation. We illustrate this topology in Fig. 4, which shows density field maps along the major three axes, and velocity streamlines colour-coded by their $z$-component of the Mach number. The left column presents the initial state ($a=1$) and the right column a highly compressed stage for which the flow has had time to settle into a self consistent non-driven mode ($a=10^{-3}$). The runs presented here are the HD run (top), noGF-Strong (middle) and GF-Medium (bottom panels), but the other MHD runs are qualitatively the same as the two shown in the figure.

Both the MHD runs exhibit a behaviour such that after compression has taken place, the flow settles into two main sheets sliding across each other at supersonic velocities. Since there is no preferred direction in the noGF run, and in the x-y plane of the GF run, in the simulation setup, the division into the domains is arbitrary (for this specific simulation it roughly coincides with the $x$-axis. It is clear that this flow, which has naturally developed from standard turbulence, is highly non-random, and that the expected dissipation of these flows is greatly reduced compared to the standard flow of the HD run or the initial state.

We can illustrate the reduces dissipation more directly by examining the ratio of flow power in compressible modes to the total power in all modes, which we refer to as the compressive ratio. We show the time evolution of this quantity for all runs in Fig. 5. We compute the energy in solenoidal ($\langle v_{\mathrm{s}}^{2}\rangle$) and compressible ($\langle v_{\mathrm{c}}^{2}\rangle$) modes by performing a Helmholtz decomposition of the velocity field (Federrath et al., 2010; Federrath et al., 2011a; Pan et al., 2016; Jin et al., 2017). The ratio of $\langle v_{\mathrm{c}}^{2}\rangle/(\langle v_{\mathrm{s}}^{2}\rangle\!+\!\langle v_{\mathrm{c}}^{2}\rangle)\sim 0.2$–$0.6$ for supersonic turbulence depends on the driving mode (Federrath et al., 2011a). In the absence of magnetic fields, the flow remains in this range of values even after driving ceases, during the compressive phase (the HD case). However, Fig. 5 (top panel) demonstrates that when the magnetic field begins to dominate ($a\lesssim 0.01$ for noGF-Medium, and $a\lesssim 0.03$ for noGF-Strong; see also Fig. 3, top panel), the compressive ratio drops rapidly. These values are marked by the vertical dashed lines in Fig. 5 (top panel).

The GF-Strong run is particularly noteworthy in that it has a compressive ratio $\lesssim 10\%$ even before the onset of compression, simply as a result of the strong magnetic field that prevents flows across field lines. As a result, it never experiences significant dissipation, and never goes through a phase when the turbulence decays.

4 Theoretical framework for compressible MHD turbulence

Having seen that magnetic fields lead to novel and initially-unexpected effects in MHD turbulence, we now seek to construct a theoretical model that we can use to interpret the results. Our basic approach will be to think of the region we are simulating as a small portion of a much larger cloud. We will then coarse-grain the MHD equations over the scale of our box, allowing us to write down an effective pressure in the box. We will use the results of our numerical experiments, together with some basic physical arguments, to provide an effective equation of state to describe this pressure and its evolution, so that we can interpret our numerical results in thermodynamic terms.

4.1 Coarse-grained pressures

We begin by following the usual method of constructing a set of coarse-grained equations (e.g., Germano, 1992; Kuncic & Bicknell, 2004; Schmidt et al., 2006; Schmidt & Federrath, 2011). We define a spatial filter $F_{\Delta}(\mathbf{x})$ with characteristic scale $\Delta$ with which we can convolve all the fluid variables. For any field $\phi(\mathbf{x})$, we define

[TABLE]

Here $\overline{\phi}$ is the filtered variable, obtained by convolving $\phi$ with the filter, and $\phi^{\prime}$ is the fluctuating part that remains after the filtered part has been removed.

Convolving the MHD equation of momentum conservation, eq. 2, with the filter $F_{\Delta}(\mathbf{x})$ gives

[TABLE]

where $\mathbf{I}$ is the identity tensor. Per the usual approach, we now write the averages over correlated terms as differences of the filtered quantities and the sub-filter-scale (SFS) quantities,

[TABLE]

where

[TABLE]

are the Reynolds stress and Maxwell stress exerted by the SFS components of the fluid velocity and magnetic field, respectively.

As standard for the microphysical stress tensor, we decompose the SFS stresses $\bm{\tau}_{\rm R,SFS}$ and $\bm{\tau}_{\rm M,SFS}$ into on- and off-diagonal components, and identify the former as effective pressures. That is, we define the effective kinetic and magnetic pressures by

[TABLE]

We use the notation $P_{\rm B}$ for the effective magnetic pressure to distinguish it from $P_{\rm mag}$, the true, microphysical magnetic pressure, since we shall see below that they are somewhat different. For homogenous, isotropic turbulence the tensors $\bm{\pi}_{\rm R,SFS}$ and $\bm{\pi}_{\rm M,SFS}$ have zero on their diagonals. In the presence of a large-scale guide field where isotropy is broken, this is not necessarily the case, and in principle the on-diagonal components of $\bm{\pi}_{\rm R,SFS}$ and $\bm{\pi}_{\rm B,SFS}$ can be as large as $P_{\rm kin}$ and $P_{\rm B}$. However, since we are only after a heuristic model, we will ignore this complication. With these definitions, the filtered momentum equation reads

[TABLE]

The final step in defining the coarse-grained pressures via an equation of state is to relate the pressures as we have defined them to the energy content of the gas. The SFS kinetic and magnetic energies per unit volume are simply the differences between the true energies per unit volume and their analogs defined using the filtered quantities, i.e.,

[TABLE]

From the definitions of $\bm{\tau}_{\rm R,SFS}$, $\bm{\tau}_{\rm M,SFS}$, $P_{\rm kin}$, and $P_{\rm B}$, it is immediately clear that we have

[TABLE]

i.e., the kinetic pressure is simply $2/3$ of the sub-filter-scale kinetic energy density, and the magnetic pressure is $1/3$ the scale of the sub-filter-scale magnetic energy density. Note that this relationship between pressure and energy density is different than the ones that obtain between the microscopic pressures and energy densities, for which $P_{\rm th}=(\gamma-1)e_{\rm th}$ and $P_{\rm mag}=e_{B}=B^{2}/8\pi$. This difference is the reason for the factor of $1/3$ we introduce into $P_{B}$ as computed in eq. 20. Interpreted in terms of an adiabatic index $\gamma$, we see that SFS kinetic pressure acts like a fluid with a $\gamma=5/3$ equation of state, while SFS magnetic pressure acts like a fluid with a $\gamma=4/3$ equation of state.

If the thermal pressure also obeys an equation of state

[TABLE]

where $\epsilon$ is the energy per unit mass, then we can write the total coarse-grained pressure as

[TABLE]

where the various $\epsilon$ terms are the thermal, SFS kinetic, or SFS magnetic energy per unit mass.

4.2 An effective EoS for supersonic magnetized gas

We now wish to model the reaction of the full system of isothermal, turbulent and magnetized gas to compression, taking into account dissipation terms and interactions between the various components. Since these additional energy transfers are time dependent, they cannot be modelled as a proper EoS, that is only a function of the thermodynamic state. Instead, we model this behaviour as an effective EoS which also depends on the thermodynamic trajectory of a parcel of gas. A similar approach has been successfully implemented for stability analysis of gravitationally collapsing haloes, filaments and sheets in a cosmological context in Birnboim & Dekel (2003), Dekel & Birnboim (2006), and Birnboim et al. (2016). For the sake of notational simplicity, we shall from this point forward drop the overlines and the SFS notation, and we will understand that, unless otherwise stated, all quantities except energies are unit mass are filtered quantities, while specific energies are SFS quantities.

In analogy to the definition of the adiabatic index (with the subscript $s$ indicating constant entropy)

[TABLE]

we define $\gamma_{\rm eff}$ as the full derivative of the pressure to density of a fixed parcel of gas along a Lagrangian path,

[TABLE]

with the upper dot indicating full time derivative and $P_{\rm tot}$ as defined in eq. (39).

4.2.1 Ideal EoS with dissipation

The time derivative of an ideal EoS (eq. 38) can be separated into its isentropic and non-isentropic part. We first differentiate the pressure of a Lagrangian parcel of gas:

[TABLE]

The time derivative of the specific energy, $\dot{\epsilon}$, can be taken from its thermodynamic definition,

[TABLE]

with $V$ the specific volume ($V=\rho^{-1}$) and $q$ a general non-adiabatic energy sink rate. A negative value of $q$ corresponds to an energy source. Inserting eq. 44 into eq. 43 and using eq. 38 again, we find,

[TABLE]

The interaction between two forms of energy, such as the transfer between kinetic and magnetic components associated with dynamo action, can be incorporated into this framework by introducing a positive $q$ term into one form of energy (for example the kinetic), compensated by a negative contribution of equal magnitude to the other (for example magnetic energy).

4.2.2 Kinetic EoS with dissipation

Following Mac Low et al. (1998) and Robertson & Goldreich (2012) we model the dissipation rate of turbulence as proportional to the largest eddy turnover time (see eq. 18),

[TABLE]

with $\eta$ a dimensionless free parameter (that depends on the numerical scheme and resolution) and $\lambda$ the largest eddy scale ($L/2$). When decay is efficient, we expect $\eta$ to be of order unity, but once a dissipationless flow pattern develops, it will be much smaller. This term dissipates kinetic into thermal energy, and formally should appear with a negative sign in the thermal component. However, by using an isothermal EoS for the gas, the thermal energy of the gas is fixed, and any heating of the gas is assumed to radiate out instantly. We simply introduce this term as cooling, directly from the kinetic pressure.

4.2.3 Energy transfer between the kinetic and magnetic components

Turbulence enhances initially small seeds of the magnetic field via small-scale dynamo processes (Brandenburg & Subramanian, 2005). The amplification rate is initially exponential and eventually decreases to zero as the magnetic energy and the turbulent energy approach their saturation ratio, which is a function both of the Mach number and turbulent driving pattern (Federrath et al., 2011a) and the ratio of turbulent and magnetic dissipation, i.e, the magnetic Prandtl number (Schekochihin et al., 2007; Federrath et al., 2014b; Schober et al., 2015). Inspired by this behavior, we model the energy transfer rate from the kinetic to magnetic as

[TABLE]

where we recall that $\epsilon_{\rm B}$ and $\epsilon_{\rm kin}$ are the sub-filter scale magnetic and kinetic energies. The saturation ratio, $\sigma_{\rm sat}=\epsilon_{\rm B}^{\rm sat}/\epsilon_{\rm kin}^{\rm sat},$ is the ratio of the two components before compression, and will be calibrated for each simulation depending on its setup. By analogy with the dissipation rate, the growth rate of the magnetic field, $\Gamma,$ is also taken as a fraction of the largest eddy turnover rate,

[TABLE]

with the dimensionless coefficient $\eta_{\rm B}$ of order unity.

While this internal energy transfer does not change the total energy content of the gas, it does change the total pressure, because of the difference in $\gamma$ for each component. We thus treat it as a sink term in $\dot{\epsilon}_{\rm kin}$ and as a source term in $\dot{\epsilon}_{\rm B}.$ We also note that this rate can become negative if the magnetic energy exceeds its saturation level with respect to the kinetic energy. While energy in such a physical case will flow from magnetic to kinetic, as the fields are strong enough to rearrange the material into a lower magnetic energy state, there is no justification to assume that the rate at which this happens is related to the eddy turnover rate. Regardless, we use eq. 47 even for that case. As we show later, this modeling, with the same coefficients for $\epsilon_{\rm B}\rightarrow\epsilon_{\rm kin}$ as $\epsilon_{\rm kin}\rightarrow\epsilon_{\rm B},$ leads to reasonable results, although the value of $\eta_{\rm B}$ that we end up requiring is considerably less than unity. We leave further investigation into this point for later studies.

4.3 Calculation of $\gamma_{\rm eff}$

We have now everything in place for the calculation of $\gamma_{\rm eff}.$ It is convenient to write the total pressure (eq. 39) as:

[TABLE]

with

[TABLE]

with $\mathscr{M}$ the rms Mach number of the flow, and $\beta$ the plasma $\beta$ parameter for the coarse-grained case.

If the sound speed is constant as we have assumed, and we use the relation in eq. 43 for the kinetic and magnetic parts, we get:

[TABLE]

By noting that (c.f., eqs. 46 and 47)

[TABLE]

and using the relation

[TABLE]

and eq. 50, we get:

[TABLE]

Plugging eq. 54 into eq. 42, we finally get:

[TABLE]

Eq. (55) manifests a few important physical properties. The equation properly interpolates between various extremes: when the thermal component is dominant over the kinetic and magnetic components, $1\gg\alpha_{k},\beta^{-1}$, $\gamma_{\rm eff}$ reduces to $1$, as expected (for an isothermal gas). Likewise, when the kinetic component dominates, $\alpha_{k}\gg 1,\beta^{-1}$, $\gamma_{\rm eff}$ reduces to $5/3$ (and to $4/3$ when the magnetic component dominates). The sign of the dissipation term and the energy transfer term depends on the sign of $\dot{a}.$ This is consistent with the expected behaviour that the dissipation always acts to reduce the pressure. When gas contracts ($\dot{a}<0$), $\gamma_{\rm eff}$ drops, and the pressure growth is reduced. When gas expands ($\dot{a}>0$), $\gamma_{\rm eff}$ increases, and the pressure drop due to the expansion, $P\propto\rho^{\gamma_{\rm eff}}$ drops even faster because of the dissipation. The magnitude of the dissipation term is not theoretically bound, and $\gamma_{\rm eff}$ can become negative or very large. This does not contradict our physical understanding. When large dissipation is present ($v\gg|\dot{a}|$), it is possible that as gas is compressed its pressure decreases, corresponding to a negative $\gamma_{\rm eff}$.

4.4 Comparison to simulations

In this section we analyze the reaction of multi-component gas to contraction, and compare the simulated runs (§2 and §3) to our analytic predictions. To quantify this reaction we use $\gamma_{\rm eff}$ (eq. 42), and, in order to compare the simulations to our predictions, we derive this quantity in two separate ways. First, we take the logarithmic derivative of the total pressure $P_{\rm tot}$ (eqs. 19, 20, 39) with respect to density directly from the simulation, $\gamma_{\rm sim}.$ Then, we compare it to our analytic prediction, $\gamma_{\rm pred},$ according to eq. 55. It is of pedagogical and practical value to consider first a simplified version of our analytic model, for which zero dissipation is assumed. This is discussed in §4.4.1. A comparison to the full model is presented in §4.4.2.

4.4.1 $\gamma_{\rm eff}$ comparison without dissipation

Dissipation sinks energy from the gas as it is compressed. Therefore, we expect gas to be more compressible (i.e., lower $\gamma_{\rm eff}$) when dissipation cannot be neglected. By contrast, if dissipation is weak enough, $\gamma_{\rm eff}$ approaches the value expected for a gas with a mixture of thermal, turbulent, and magnetic pressure, with the different components weighted according to the relative pressures of each component. In terms of our model, this amounts to setting $\eta=\eta_{\rm B}=0$. As discussed in the previous sections, we expect this approximation to be valid at the initial stage, and then, to some degree, at late stages when a dissipationless flow (DLF) forms.

Fig. 6 compares $\gamma_{\rm sim}$ (solid lines) with $\gamma_{\rm pred}$ (dashed lines). The horizontal dashed lines mark values of $\gamma_{\rm eff}=1$, $4/3$, and $5/3$ that correspond to the expected behavior for purely isothermal, purely magnetic and purely kinetic gas, respectively. For the HD simulation, the gas initially starts close to the value predicted for supersonic turbulence, $\gamma_{\rm eff}\sim 5/3$, and as dissipation becomes dominant $\gamma_{\rm sim}$ drops to $\approx 1/2$. Once the turbulence has decayed significantly (Fig. 1 top panel) and the kinetic pressure drops below the thermal pressure (Fig. 3), no energy is left to dissipate, and the gas behaves like an isothermal gas with $\gamma_{\rm eff}=1$. The MHD runs without magnetic guide fields (noGF runs, top panel) initially behave similarly, with $\gamma_{\rm eff}$ dropping as dissipation becomes important. However, the dissipation stops for a different reason. At late times, after the transition to DLF, the gas behaves as magnetic field-dominated, with $\gamma_{\rm eff}\approx 4/3$. The asymptotic behavior at the initial stage and then again for $a\lesssim 10^{-2}$ is well recovered by our analytic model without dissipation (top panel, dashed lines), for the HD and noGF runs. Note that at $a\approx 1$ the predicted values are slightly higher than the simulated ones, indicating that some dissipation exists even for the very early stages when compression starts.

The bottom panel presents a similar calculation, but for the runs with guide fields (GF runs). In all three GF runs, $\gamma_{\rm sim}$ approaches $\gamma_{\rm eff}\approx 4/3$ at large compressions. This occurs even when the magnetic pressure is sub-dominant to the kinetic one (GF-Weak) or comparable (GF-Medium and GF-Strong), and our model predicts that the gas behaves like a purely kinetic gas (GF-Weak, $\gamma_{\rm pred}\approx 5/3$) or intermediate (GF-Medium and GF-Strong, $4/3<\gamma_{\rm pred}<5/3$). This discrepancy indicates that some dissipation is present even at those late stages. In summary, the simplified non-dissipative model provides a very good prediction of the asymptotic values of $\gamma$ for the HD and noGF runs, and provides a reasonable prediction but with a slight overestimate for the asymptotic $\gamma_{\rm eff}$ of the GF runs.

4.4.2 $\gamma_{\rm eff}$ comparison of full model

Fig. 7 shows a similar plot to Fig. 6, but with $\gamma_{\rm pred}$ that includes dissipation as the dashed lines. First examine the top panel. In this plot, we have used $\eta=1.8$ and $\eta_{\rm B}=0.01$ at compressions $a\gtrsim 10^{-2}$. These are by-eye fits; a more systematic calibration is certainly possible, but since it is likely that $\eta$ depends on the numerical scheme and resolution, and $\eta_{\rm B}$ depends on the specific geometry of the driving, we find little practical value in estimating them more robustly here. To obtain a meaningful physical results, we would need to explicitly include physical dissipation terms, i.e., kinematic viscosity and magnetic resistivity (Federrath et al., 2011a; Federrath et al., 2014b; Federrath, 2016).

For our purpose, it suffices to show that we can fit the numerical results with reasonable values of $\eta$ and $\eta_{\rm B}$. The inclusion of the kinetic dissipation term, $\eta v/aL$ improves the small mismatch at $\rho\approx 1$ that was present in Fig. 6. This illustrates that for our simulation setup, dissipation makes a minor difference from the beginning of the compression phase, even though it was not directly observed in Fig. 3. For the HD and noGF runs (top panel), the predicted $\gamma_{\rm eff}$ curve is a reasonably good fit to the simulated one, even when the dissipation dominates.

However, the fit would fail miserably when flows enter the self-avoiding channel-flow phase. At that stage, our model would predict an ever-increasing dissipation because of the rising turbulent velocity when, in fact, very little dissipation takes place. Similarly, as the magnetic energy exceeds its saturation level, more and more energy would be predicted by our model to be pumped from magnetic to kinetic when, in fact, little energy does. To mimic this drop in transfer terms, we simply set $\eta=\eta_{\rm B}=0$ when $(\alpha_{k}\beta)^{-1}=P_{\rm B}/P_{\rm kin}>3$. In the Figure, the shutdown of the transfer terms produces the discontinuity in the predicted $\gamma_{\rm eff}.$ Without it, $\gamma_{\rm eff}$ would first shoot up because of the large $\dot{\epsilon}_{\rm KB}$, and then shoot down because of the large $q_{\rm dis}.$ Obviously, the transition to the ineffective dissipation regime is not sharp, and a better fit between the calculated and predicted models can be obtained with a more sophisticated transition model. However, at this stage, we value the simplicity of the model over its accuracy. The success of the dissipation model when $P_{\rm B}\lesssim P_{\rm kin},$ combined with its failure beyond that point and the success of the dissipationless model there, is the key physical finding of this paper.

Our model is somewhat less successful in the case when a guide field is present, as shown in the lower panel of Fig. 7. The model is reasonable at first, but, as noted above, the switch to DLF appears to depend not just on the ratio of magnetic to kinetic pressures, but also on the field topology. Our switch at $P_{\rm B}/P_{\rm kin}>3$ is too conservative for this case, and as a result predicts strong dissipation at intermediate values of $a$ when in fact our simulations are already switching to a low-dissipation state. We could improve the fits by hand-tuning when we switch to $\eta=\eta_{\rm B}=0$, but without an understanding of exactly how this switch depends on the magnetic topology, we would have to do so independently for each run, which seems of little value.

Regardless of our ability to predict exactly when the switch to DLF will occur, we can still give a physical interpretation to our results as follows. With our simulation setup, the turbulence at the compression stage is driven by globally re-normalizing the velocity and distance. This type of driving enhances existing modes, and allows the system to settle into “natural modes” which are not forced by an arbitrary external driving. From our results it appears that even a relatively small decay of some of the modes (e.g., even in the case of a sub-dominant magnetic field) is enough to significantly decrease the power in these modes and allows the flow to settle into a DLF.

5 Implications for ISM

While the study of dissipation in scale-free compression of turbulence described above is general, our motivation for conducting this study is not. We seek to find how magnetic fields alter the collapse of GMCs and, in particular, to test whether realistic initial conditions may delay the dissipation and consequent collapse of the GMCs. The ISM, in which GMCs form, is a multi-phase medium, with subsonic turbulence in the diffuse, warm and hot phases and supersonic turbulence in the denser, cold phases such as GMCs. The entire ISM is furthermore interlaced with magnetic fields and immersed in cosmic rays and radiation fields. The turbulence is maintained by various energy sources, including driving by feedback from star formation, gravitational collapse, and galaxy dynamics (Federrath et al., 2016a, 2017). In this work we simplify this system to manageable levels by focusing on its most basic aspects: we start with gas that is supersonic and magnetized and test its response to quick compression. The compression pumps energy into the turbulence and into the magnetic field on all scales, without imposing any particular scale or randomness through external driving. We find that even weak magnetic fields eventually force the gas to settle into a channel-flow pattern with greatly-reduced dissipation. This is in stark contrast to the case where magnetic fields are absent, where we find that (in agreement with previous purely hydrodynamical simulations – Robertson & Goldreich 2012) that the kinetic energy steadily decays and dissipation remains significant throughout the entire evolution of the gas.

The greatly reduced dissipation rate for the flow that we find in the presence of a magnetic field indicates that significantly less energy input may be required to produce the observed linewidths in GMCs than had previously been conjectured. However, a full exploration of this issue will need to take into account many more processes, among them the multiphase nature of the ISM and realistic stellar feedback. The latter process, which will drive random turbulence on a typical scale corresponding to the distance between stars within each GMC, may actually increase the dissipation rate by disturbing the dissipationless flow. Thus, while stellar feedback is highly energetic, it could prove to be an effective cooling agent by increasing the dissipation. Even within our toy model, a more systematic test of various compression models and rates is needed to test the universality (or lack-thereof) of the dissipation model and the calibrated rates. We leave these tests to future work.

At a minimum, however, we note that our results suggest that dissipation inside a forming GMC will be much less than has commonly been assumed. The gas from which GMCs form is observed to be threaded by a significant net magnetic flux (Li & Henning, 2011; Li et al., 2011; Pillai et al., 2015), and in this case the flow can be nearly dissipationless almost immediately after compression begins. Even in the limiting case of zero net flux but fields close to equipartition, as observed, the dissipation rate is substantially reduced once the cloud has compressed by a factor of $\sim 100$ in linear dimension. This is not all that much by interstellar standards: the mean density of the Milky Way’s ISM is $n\sim 1$ cm*-3*, so a factor of 100 linear compression corresponds to a density $n\sim 10^{6}$ cm*-3*, i.e., typical of a prestellar core.

The onset of dissipationless flow will not necessarily halt collapse. Even in the most favorable case, when no dissipation occurs, highly magnetized gas is unstable to scale-free collapse because its effective EoS is $\gamma_{\rm eff}=4/3$, which is right at the critical value for hydrostatic gravitational stability ($\gamma_{\rm crit}=4/3$). This indicates that, as a cloud compresses, the force exerted by the magnetic field outwards grows just as fast as the gravitational forces increases inwards. This general point has been overlooked in the past, because the focus was on comparing timescales rather than using the formal requirement for hydrostatic atmospheres. However, we note that if the external compression is filamentary as is often found, the critical value for stability drops to $\gamma_{\rm eff}=1,$ and, if linear magnetic fields prevent compression perpendicular to the magnetic-field direction, the collapse can only occur along one direction, for which the critical value is $\gamma_{\rm eff}=0$ (Birnboim et al., 2016).

The compression rate of observed GMCs can be characterized in a way that relates it to our ideal simulations. We define the non-dimensional compression rate $H_{\mathrm{ND}}=-t_{\mathrm{sound}}H$ which is our (absolute value of the) Hubble coefficient for the compression in units of sound crossing time. If compression is driven by gravitational collapse, this non-dimensional compression rate is (up to order of unity corrections) $t_{\mathrm{sound}}/t_{\mathrm{ff}},$ with $t_{\mathrm{ff}}$ the gravitational free-fall time. By replacing $t_{\mathrm{sound}}$ of our turbulent cloud by $\mathcal{M}t_{\mathrm{turb}}$ (which is, again, correct up to order of unity corrections) we get $H_{\mathrm{ND}}=\mathcal{M}t_{\mathrm{turb}}/t_{\mathrm{ff}}$. The ratio of the turbulent timescale to the free fall timescale for reasonable GMC’s is $t_{\mathrm{turb}}/t_{\mathrm{ff}}=2/\sqrt{\alpha_{\rm vir}},$ with $\alpha_{\rm vir}\approx 1,$ the virial parameter for virialized GMCs (see sec. 8.3 of Krumholz, 2017), yielding $H_{\mathrm{ND}}\approx 2\mathcal{M}.$ In our numerical simulations $t_{\mathrm{sound}}=1$ and $H_{\mathrm{ND}}=-H$ spanning a range between $1$ and $200$ (see Table 1 and Appendix §B). Since typical Mach numbers for local GMCs (Schneider et al., 2013b) and in the Central Molecular Zone Cloud G0.253+0.016 (Federrath et al., 2016b) are $\mathcal{M}\approx 10$ and for high-z GMCs or ULIRGS can be as high as $\mathcal{M}\approx 100,$ we argue that our simulations bracket the observed range. Furthermore, as gas is compressed, dissipation always dominates eventually over the adiabatic compression. Faster compression rate simply extends the initial stage in which dissipation can be neglected. We note, however, that unlike realistic GMCs, our Mach number is independent of the compression rate and is set by the initial driving. Our Mach numbers ($\mathcal{M}\approx 10$) are smaller than is observed for ULIRGS, and is perhaps comparable to that of local GMCs. Additionally, it does not follow the correlation set by $H_{\mathrm{ND}}\approx 2\mathcal{M}$ that is expected by observations. Since the onset of the dissipation is predicted by comparing the dissipation timescale (ie. turbulent timescale) and compression timescale, we expect that realistic GMCs will always start near the onset of the dissipation phase, without the long adiabatic phase seen for $H=-200$ in our simulations. We leave a more systematic test of the dependence of our conclusion on the compression rate and Mach numbers for future work.

6 Summary and conclusions

In this paper, we study the evolution of magnetized supersonically turbulent gas as it is compressed. The compression is scale-free and corresponds to gravitational compression that operates on all scales, much like the expansion of the universe in cosmology, but with a scale factor that decreases with time (negative Hubble constant). Our simulations of this scale-free compression are performed by using a modified version of the cosmological expansion model in the FLASH MHD code (§2), and explore a range of magnetic field strengths and net fluxes.

The scale-free compression enhances all turbulent and magnetic modes by the same factor, and does not impose any arbitrary scale or randomness of phases. Consequently, the system is allowed to relax into a self-consistent state for which naturally decaying modes decay away while non-decaying modes are enhanced by the compression. We find that this relaxation, combined with a magnetic field, produces a surprising result: after some time the gas re-arranges itself into a self-avoiding channel flow, in which state the dissipation rate is nearly zero. This occurs whether or not there is a net magnetic flux, but the transition happens more readily for non-zero magnetic flux and for stronger fields, and it does not occur at all in the absence of a magnetic field.

We interpret the simulations by comparing them to a theoretical model for coarse-grained MHD turbulence. Our model treats the flow has having three distinct energy reservoirs (thermal, kinetic, and magnetic) that are coupled by dynamo action and dissipation. This model allows us to construct an equation of state with an effective adiabatic index $\gamma_{\rm eff}$, whose value depends on the relative balance between the different energy reservoirs and on the overall rate of dissipation. We calibrate the transfer and dissipation terms from the simulations, and show that, once calibrated, the model provides a good match to our numerical experiments. A key features of the model, and of the numerical results it describes, is that once the flow is sufficiently magnetically-dominated the dissipation rate for the flow is nearly zero, and compression drives a continual increase in the kinetic energy per unit mass and the Mach number.

The existence of this dissipationless state may have significant implications for the ISM, and for giant molecular clouds in particular. These clouds are assembled by compression of regions of atomic ISM with significant magnetic flux. If this compression resembles the idealized gravitational contraction we consider here, then GMCs may be in flow states where the dissipation rate is much less than is commonly assumed based on earlier work in non-compressing regions and on non-magnetized flows. This in turn might significantly ease the problem of how GMCs maintain their large linewidths while still forming stars with very low efficiencies.

Finally, we note that the analytic model derived here constitutes a first step towards a physically-motivated sub-resolution model for the ISM. Given some idea (perhaps calibrated by observations; see Birnboim et al., 2015) about the content of the turbulence, thermodynamics and magnetic fields of the ISM, our model predicts the behavior of the gas in a way that can be implemented into large-scale, low-resolution simulations that only resolve the ISM as a coarse-grained mixture of the components. This too will be investigated in future work.

acknowledgements

We thank Chalence Safranek-Shrader and Romain Teyssier for their help during the implementation of the Hubble source terms for MHD (§2). Y.B. wishes to thank the Research School of Astronomy and Astrophysics at The Australian National University for hosting him on a sabbatical during 2016–2017. C.F. gratefully acknowledges funding provided by the Australian Research Council’s Discovery Projects (grants DP150104329 and DP170100603). M. R. K.’s work was supported under the Australian Research Council’s Discovery Projects funding scheme (project DP160100695). The simulations presented in this work used high performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391), the Partnership for Advanced Computing in Europe (PRACE grant pr89mu), the Australian National Computational Infrastructure (grant ek9), and the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia, in the framework of the National Computational Merit Allocation Scheme and the ANU Allocation Scheme. The simulation software FLASH was in part developed by the DOE-supported Flash Center for Computational Science at the University of Chicago.

Appendix A Numerical tests

A.1 The $\nabla\cdot B=0$ constraint

Physically, no energy should propagate into longitudinal modes of the magnetic field. For highly supersonic, low plasma $\beta$ turbulence simulations in FLASH using the HLL3R and HLL5R solver (Waagan et al., 2011a) this has been demonstrated and compared with alternative schemes in Waagan et al. (2011b). However, the compression in the simulations presented here could, in principle, change this conclusion because the amplitude of the magnetic fields is enhanced due to adiabatic compression. Fig. 8 presents the ratio of magnetic field energy in longitudinal modes to the total energy in magnetic fields. As is evident, this value is never larger than $10^{-4}$, and, except for GF-Weak, smaller than $10^{-5}$. Additionally, this value drops as compression occurs, and is largest near the comencement of compression. We therefore do not expect the numerical errors in $\nabla\cdot B$ to effect the conclusions of this paper.

A.2 Convergence with numerical grid resolution

It is well known that the necessary grid resolution for simulating fully-developed turbulence with an inertial subrange is at least $1024^{3}$ cells (e.g., Klein et al., 2007; Kritsuk et al., 2007; Schmidt et al., 2009; Federrath et al., 2010; Federrath, 2013). However, this was not our goal in this paper. We expect that for our particular focus on coarse-grained values averaged over large scales, such a high resolution is not critical. In the interest of computational efficiency and to allow for many different simulations, we only ran $512^{3}$-cell boxes. Indeed, previous studies have demonstrated that large-scale averages converge even at a resolution of $256^{3}$ grid cells (Kitsionas et al., 2009; Federrath et al., 2010; Price & Federrath, 2010; Kritsuk et al., 2011). In this appendix we briefly demonstrate that our resolution of $512^{3}$ cells is sufficient for our needs.

We check for convergence by comparing two physically similar MHD runs: noGF-Strong and noGF-Strong-LR (see table 1). The two simulations differ only in that the the former has a lower resolution of $256^{3}$, while the latter has our standard $512^{3}$ resolution. We find that, at the end of the initial driving stage, the Mach numbers for the two runs differ by $5\%$, while the saturation levels for the magnetic field strength differ by about $30\%$. While these differences are not negligible, they do not change the qualitative results once compression begins. We demonstrate this in Fig. 9, which shows the Mach number evolution, and Fig. 10, which shows $\gamma_{\rm eff}$; these figures can be compared to Fig. 1 and Fig. 7 in the main text. We see that the overall behaviour of both the Mach number and the adiabatic index are nearly identical in the two runs.

Appendix B Slow ($H=-1$) compression

All the simulations used in the main text have $H=-200$ corresponding to very fast compression. The motivation for this choice is that, given our high Mach numbers ($\mathcal{M}\approx 10$, motivated by the properties of observed molecular clouds), the dissipation timescale is very short. Consequently, values of $H$ near unity would not show a distinct amplification stage before the onset of compression, and instead would proceed directly to the dissipation stage. We find that ensuring the presence of a distinct amplification phase helps elucidate the physics of the problem, which is why we elected to use $H=-200$ as our standard choice. However, it is important to demonstrate that our central result, the onset of dissipationless flow, is independent of this choice. For this reason, in Fig. 11 we show the evolution of the Mach number for three simulations (HD-H1, noGF-Strong-H1 and GF-Medium-H1; see table 1) with a compression rate $H=-1,$ as well as our fiducial HD run to guide the eye. This compression rate corresponds to the gas compression at roughly the sound speed, significantly slower than what would be expected if the gas were to contract at free-fall.

In these runs, the Mach number declines immediately once driving is turned off and compression starts, as we would expect for such slow compression. However, in these runs we still see the characteristic increase in Mach number at late times that occurs due to the onset of dissipationless flow. We demonstrate this more clearly in Fig. 12, which shows $\gamma_{\rm eff}$ measured for the $H=-1$ simulations, as compared to our theoretical model using the same values of $\eta$ and $\eta_{\rm B}$ (including the condition when we set these terms to zero) calibrated from the $H=-200$ simulations, and used for Fig. 7 in the main text. We first note that, as in the $H=-200$ simulations, at late times the MHD runs have $\gamma_{\rm eff}$ substantially above unity, demonstrating the reduced dissipation that is the central result of this work. Moreover, the plot shows that our model continues to provide a very good description of the value of $\gamma_{\rm eff}$. Since dissipation is dominant from the beginning, our model predicts the initial value of $\gamma_{\rm eff}$ to be very low, as is the case. At late times, our model predicts that $\gamma_{\rm eff}$ should increase as energy is pumped into the magnetic component, and that eventually dissipation should turn off, leading to an increase in $\gamma_{\rm eff}$. Both of these predictions are borne out by the simulations. The success of our model’s $\gamma_{\rm pred}$ at reproducing the value $\gamma_{\rm sim}$ we measure from the simulations, using the same calibration of $\eta$ and $\eta_{\rm B}$, suggests that our results are not dependent on the choice of a particular compression rate $H$.