Effects of Forcing on Shocks and Energy Dissipation in Interstellar and Intracluster Turbulences

The dynamics of isothermal, compressible, magnetized gas is described by the following MHD equations:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}+{\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{v}}+\displaystyle \frac{1}{\rho }{\boldsymbol{\nabla }}P-\displaystyle \frac{1}{\rho }({\boldsymbol{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}-{\boldsymbol{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}})=0,\end{eqnarray} \tag{ 3 }$$

where ρ is the gas density, v and B are the velocity and magnetic field vectors. The gas pressure is given as $P=\rho {c}_{s}^{2}$ with a constant sound speed c_s . Note that the unit of B is chosen so that 4π does not appear in Equation (2).

For simulations of turbulence, we solve the above equations numerically using a MHD code based on the WENO scheme (Jiang & Shu 1996; Jiang & Wu 1999). The WENO code used in this work has fifth-order spatial and fourth-order temporal accuracies, respectively, while the Total Variation Diminishing (TVD) code used previously in PR2019 is second-order accurate in both space and time. The Appendix briefly describes the WENO code, including a comparison with the TVD code. Viscous and resistive dissipations are not explicitly included. Simulations are performed in a three-dimensional (3D) cubic box of size L₀ = 1 with 256³ and 512³ uniform Cartesian grid zones. The background medium is initialized with ρ₀ = 1, P₀ = 1 (i.e., c_s = 1), and v ₀ = 0. The initial magnetic field is placed along the x-axis with an uniform strength B₀, which is specified by the plasma beta, ${\beta }_{0}={P}_{0}/{P}_{{\rm{B}},0}=2{c}_{s}^{2}{\rho }_{0}/{B}_{0}^{2}$. We adopt β₀ = 0.1 for the ISM turbulence and β₀ = 10⁶ for the ICM turbulence.

Turbulence is driven by adding a small velocity perturbation, δ v, at each grid zone at each time step; δ v is drawn from a Gaussian random distribution with a Fourier power spectrum, $| \delta {{\boldsymbol{v}}}_{k}{| }^{2}\propto {k}^{6}\exp (-8k/{k}_{\exp })$, where ${k}_{\exp }=2{k}_{0}$ with k₀ = 2π/L₀ (Stone et al. 1998; Mac Low 1999). The injection scale is regarded as the peak of ∣δ v _k ∣² k², ${k}_{\mathrm{inj}}={k}_{\exp }$, i.e., L_inj = L₀/2. The perturbations have random phases, hence the driving is temporally uncorrelated. The amplitude of δ v is adjusted, so that turbulence saturates with the rms velocity of flow motions, v_rms = 〈v²〉^1/2 ≈ M_turb c_s . We aim to obtain M_turb = 10 for supersonic ISM turbulence in molecular clouds, and M_turb = 0.5 for subsonic ICM turbulence.

Forcing with δ v generally leads to a combination of solenoidal and compressive components. By separating the two components in Fourier space, we construct three types of forcing: (1) fully solenoidal forcing with ∇ · δ v = 0, (2) fully compressive forcing with ∇ × δ v = 0, and (3) a mixture of the two. It has been argued that the driving of turbulence in ISM molecular clouds could be specified by a mixture of solenoidal and compressive modes with a 2:1 ratio (Federrath et al. 2010). Hence, we consider the case of 2:1 mixed forcing.

Table 1 summarizes the basic parameters of the turbulence models considered in this paper. The nomenclature for the models have three elements, as listed in the first column, which are self-explanatory: (1) the first element indicates either the "ISM" or "ICM" turbulence; (2) the second element denotes the forcing mode; and (3) the third element shows the grid resolution, Nn_g , where n_g is the number of grid zones in one side of the simulation box. The high-resolution models specified with N512 have 512³ grid zones, while the low-resolution models specified with N256 have 256³ grid zones. The models are defined by the two parameters, M_turb (column 2) and β₀ (column 3), and the forcing mode (column 4). Simulations for the ISM models were run up to t_end = 5 t_cross for N256 models and t_end = 2.2 − 3.5 t_cross for N512 models, whereas those for the ICM models run up to t_end = 30 t_cross regardless of the resolution. Here, t_cross = L_inj/v_rms is the crossing time. In the weakly magnetized ICM turbulence, the magnetic field amplification via small-scale dynamo is much slower in units of t_cross, and hence the timescale to reach saturation is much longer than in the ISM turbulence (see Figure 1 and the discussion in the next section). The values of t_end/t_cross are listed in the fifth column.

Zoom In Zoom Out Reset image size

Table 1. Model Parameters of Simulations and Statistics of Turbulence ^a

Model	Mturb	β0 b	Forcing	tend/tcross c	inj d	βsat b	$({N}_{\mathrm{fa}}/{n}_{g}^{2})$ e	$({N}_{\mathrm{sl}}/{n}_{g}^{2})$ e	fa d	sl d	shock/inj
ISM-Sol-N256	10	0.1	Solenoidal	5	1650	0.0407	4.69	0.640	241	2.63	0.148
ISM-Sol-N512	10	0.1	Solenoidal	2.2	1590	0.0388	5.35	1.40	207	5.62	0.134
ISM-Mix-N256	10	0.1	Mixed	5	1550	0.0429	4.66	0.612	234	2.70	0.153
ISM-Mix-N512	10	0.1	Mixed	3.5	1520	0.0419	5.35	1.42	201	5.79	0.136
ISM-Comp-N256	10	0.1	Compressive	5	1300	0.0564	3.43	0.438	181	2.37	0.141
ISM-Comp-N512	10	0.1	Compressive	3.5	1200	0.0550	4.93	1.01	185	4.58	0.158
ICM-Sol-N256	0.5	106	Solenoidal	30	0.140	5.31 × 101	6.03	0.0	0.0139	0.0	0.100
ICM-Sol-N512	0.5	106	Solenoidal	30	0.140	4.71 × 101	6.43	0.0	0.0150	0.0	0.107
ICM-Mix-N256	0.5	106	Mixed	30	0.150	5.79 × 101	8.18	0.0	0.0281	0.0	0.188
ICM-Mix-N512	0.5	106	Mixed	30	0.146	5.12 × 101	8.38	0.0	0.0281	0.0	0.193
ICM-Comp-N256	0.5	106	Compressive	30	0.340	8.76 × 102	10.7	0.0	0.119	0.0	0.349
ICM-Comp-N512	0.5	106	Compressive	30	0.340	5.59 × 102	11.7	0.0	0.113	0.0	0.331

Notes.

^a The statistics of turbulence, β_sat, N_fa, N_sl, _fa, _sl, and _shock, are the mean values at saturation, which are calculated over 2 t_cross ≤ t ≤ t_end for the ISM models and over 15 t_cross ≤ t ≤ t_end for the ICM models. Here, the subscripts "fa" and "sl" stand for fast and slow shocks, respectively, and _shock = _fa + _sl. ^b The initial plasma beta, β₀, and the plasma beta at saturation, β_sat. ^c The end time of simulations in units of the crossing time, t_cross = L_inj/v_rms (see the main text). ^d The energy-injection rate and the energy-dissipation rate at shocks in computational units of ρ₀ = 1, c_s =1, and L₀ = 1. ^e The numbers of shock zones normalized to n_g ².

Download table as:  ASCIITypeset image

In MHDs, there are two kinds of shocks, i.e., fast-mode and slow-mode shocks. The shock-identification scheme and the formulas used for the calculations of the shock Mach number and the energy dissipation at shocks are basically the same as those in PR2019, except that we here identify shocks with M_s ≥ 1.05, rather than M_s ≥ 1.06, taking the advantage of the high-order accuracy of the WENO code.

"Shock zones" are found by a dimension-by-dimension identification scheme, as follows. (1) Along each coordinate direction, grid zones are marked as "shocked", if ∇ · v < 0 and $\max ({\rho }_{i+1}/{\rho }_{i-1},{\rho }_{i-1}/{\rho }_{i+1})\geqslant {1.02}^{2}$ around the zone i. (2) Considering that shocks spread typically over two to three grid zones, the zone with minimum ∇ · v among the attached shocked zones is tagged as a "shock zone". (3) Around the shock zone, the preshock or postshock zones are chosen, depending on the density jump. (4) The Mach number of the shock zones, M_s , is calculated using the preshock and postshock quantities (see below). If a zone is identified as shocked more than once along the coordinate directions, M_s is determined as ${M}_{s}=\max ({M}_{s,x},{M}_{s,y},{M}_{s,z})$.

The formula for M_s can be derived using the shock jump condition for isothermal MHD flows from Equations (1)–(3). Following PR2019, we use

$$\begin{eqnarray}&&{M}_{s}^{2}=\chi +\displaystyle \frac{\chi }{\chi -1}\displaystyle \frac{{B}_{2}^{2}-{B}_{1}^{2}}{2{c}_{s}^{2}{\rho }_{1}},\end{eqnarray} \tag{ 4 }$$

where χ = ρ₂/ρ₁. Hereafter, the subscripts "1" and "2" denote the preshock and postshock states, respectively. The second term in the right-hand side represents the magnetic field contribution to M_s . We point out that both ${B}_{2}^{2}-{B}_{1}^{2}$ in the numerator and χ − 1 in the denominator goes to zero as M_s approaches to unity in weak shocks. In the ISM turbulence with M_turb ≈ 10, most of the shocks are strong, with M_s ≫ 1, and hence this does not pose a problem. On the other hand, in the ICM turbulence with M_turb ≈ 0.5, the shock population is dominated by weak shocks, and uncertainties may be introduced in the PDF of M_s and also the estimate of the energy dissipation at shocks. Equation (4) can be rewritten as

$$\begin{eqnarray}&&{M}_{s}^{2}=\chi +\displaystyle \frac{\chi {B}_{\perp 1}\left({B}_{\perp 2}+{B}_{\perp 1}\right)}{2{c}_{s}^{2}{\rho }_{1}}\displaystyle \frac{\eta }{\zeta },\end{eqnarray} \tag{ 5 }$$

with η = B_⊥2/B_⊥1 − 1 and ζ = ρ₂/ρ₁ − 1. Hereafter, the subscripts "⊥"and "∥" denote the magnetic field components perpendicular and parallel to the shock normal, respectively. It is the ratio, η/ζ, that may not be accurately reproduced from the numerical solutions for weak shocks. We note that η/ζ is given as ${M}_{s}^{2}/({M}_{s}^{2}-\chi {v}_{A\parallel 1}^{2}/{c}_{s}^{2})$, where ${v}_{A\parallel 1}={B}_{\parallel }/\sqrt{{\rho }_{1}}$. Hence, in the ICM turbulence with β ≫ 1, we use Equation (5), assuming η/ζ ≈ 1, for the calculation of M_s , if the magnetic field around the shock zones is weak, specifically, if ${v}_{A\parallel 1}^{2}/{c}_{s}^{2}\leqslant 0.1;$ if ${v}_{A\parallel 1}^{2}/{c}_{s}^{2}\gt 0.1$, we use Equation (4).

Identified shock zones are classified into either fast or slow shocks, according to the criterion of B_⊥2 > B_⊥1 or B_⊥2 < B_⊥1. In addition, fast shocks should satisfy ${M}_{s}^{2}\gt \chi {v}_{A\parallel 1}^{2}/{c}_{s}^{2}$, while slow shocks satisfy ${M}_{s}^{2}\lt {v}_{A\parallel 1}^{2}/{c}_{s}^{2}$. In our simulated turbulence, about 95% of identified shock zones satisfy these conditions for either fast or slow shocks. The rest have ${v}_{A\parallel 1}^{2}/{c}_{s}^{2}\lt {M}_{s}^{2}\lt \chi {v}_{A\parallel 1}^{2}/{c}_{s}^{2}$, and may be classified as "intermediate shocks", which are known to be nonphysical (e.g., Landau & Lifshitz 1960). The presence of these intermediate shocks could be partly due to our dimension-by-dimension approach for shock identification and also partly due to possible numerical errors in capturing shocks in simulations. In the next section, we present the results excluding intermediate shocks. With a fraction of about 5% or so, the exclusion should not affect the main conclusions of this work.

Although the results are presented with the sonic Mach number of shocks, M_s , in the next section the fast or slow Mach numbers can be calculated as M_fa = M_s c_s /c_fa or M_sl = M_s c_s /c_sl, using the fast and slow wave speeds in the preshock region,

$$\begin{eqnarray}\begin{array}{rcl}{c}_{\mathrm{fa},\mathrm{sl}}^{2} & = & \displaystyle \frac{1}{2}\left({c}_{s}^{2}+{v}_{A\parallel 1}^{2}+{v}_{A\perp 1}^{2}\right)\\ & & \pm \displaystyle \frac{1}{2}\sqrt{{\left({c}_{s}^{2}+{v}_{A\parallel 1}^{2}+{v}_{A\perp 1}^{2}\right)}^{2}-4{v}_{A\parallel 1}^{2}{c}_{s}^{2}},\end{array}\end{eqnarray} \tag{ 6 }$$

where ${v}_{A\perp 1}={B}_{\perp 1}/\sqrt{{\rho }_{1}}$. Hereafter, the subscripts "fa" and "sl" denote fast and slow shocks, respectively. Considering that weak shocks could be confused with waves, we examine only fast and slow shocks with M_fa ≥ 1.05 and M_sl ≥ 1.05, respectively, and also with M_s ≥ 1.05. In addition, another constraint, c_sl/c_s ≥ 0.3, is imposed for slow shocks, since slow shocks with very small c_sl are not clearly distinguished from fluctuations.

We also calculate the energy dissipation at shocks, as in PR2019. With the "heat energy" or the "effective internal energy", $P\mathrm{ln}P$, the equation for the "effective total energy" can be written as (Mouschovias 1974)

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{1}{2}\rho {v}^{2}+P\mathrm{ln}P+\displaystyle \frac{1}{2}{B}^{2}\right)\\ \quad +\,{\boldsymbol{\nabla }}\cdot \left[\left(\displaystyle \frac{1}{2}\rho {v}^{2}+P\mathrm{ln}P+P\right){\boldsymbol{v}}\right.\\ \quad \left.+\,\left({\boldsymbol{B}}\times {\boldsymbol{v}}\right)\times {\boldsymbol{B}}\right]=0.\end{array}\end{eqnarray} \tag{ 7 }$$

Then, the jump of the total energy flux in the shock-rest frame is given as

$$\begin{eqnarray}&&\begin{array}{l}\left[\left(\displaystyle \frac{1}{2}\rho {v}^{2}+P\mathrm{ln}P+P+{B}^{2}\right){v}_{\parallel }\right.\\ \quad {\left.-{B}_{\parallel }\left({B}_{\parallel }{v}_{\parallel }+{B}_{\perp }{v}_{\perp }\right)\right]}_{2}^{1}\equiv Q,\end{array}\end{eqnarray} \tag{ 8 }$$

where ${\left[f\right]}_{2}^{1}={f}_{1}-{f}_{2}$ denotes the difference between the preshock and postshock quantities. Here, Q is the energy-dissipation rate per unit area at the shock surface, which can be expressed as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{Q}{{\rho }_{1}{M}_{s}{c}_{s}^{3}}=\displaystyle \frac{1}{2}{M}_{s}^{2}\left[1-\displaystyle \frac{1}{{\chi }^{2}}\right.\\ \quad \left.+\displaystyle \frac{{v}_{{\rm{A}}\perp 1}^{2}\left(\chi -1\right)\left\{{v}_{{\rm{A}}\parallel }^{2}\left(\chi +1\right)-2{M}_{s}^{2}{c}_{s}^{2}\right\}}{{\left({v}_{{\rm{A}}\parallel }^{2}\chi -{M}_{s}^{2}{c}_{s}^{2}\right)}^{2}}\right]-\mathrm{ln}\chi .\end{array}\end{eqnarray} \tag{ 9 }$$

The dissipation rate of the turbulent energy at all shocks inside the entire simulation box is estimated as

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{fa}(\mathrm{sl})}=\sum _{\mathrm{fast}(\mathrm{slow})\ \mathrm{shocks},j}{Q}_{j}{\left({\rm{\Delta }}x\right)}^{2},\end{eqnarray} \tag{ 10 }$$

for either the fast- or slow-shock populations. Here, Δx = L₀/n_g is the size of the grid zones, and the summation goes over all the identified shock zones.

This rate is compared to the injection rate of the energy deposited in the simulation box with the forcing described above:

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{inj}}=\displaystyle \frac{1}{{\rm{\Delta }}t}\int \displaystyle \frac{1}{2}\rho \left[{\left({\boldsymbol{v}}+\delta {\boldsymbol{v}}\right)}^{2}-{{\boldsymbol{v}}}^{2}\right]{dV},\end{eqnarray} \tag{ 11 }$$

where Δt is the simulation time step. As mentioned above, the amplitude of δ v is the adjustable parameter, and the resulting values of _inj are given in the sixth column of Table 1.

We first examine the overall behavior of turbulence in our simulations. Figure 1 shows the rms velocity of turbulent flow motions, v_rms, the volume-integrated kinetic energy, E_K =∫(1/2)ρ v² dV, and the volume-integrated turbulent magnetic energy, E_B = ∫(1/2)δ B² dV, as a function of time, for all the models listed in Table 1. Here, δ B = B − B ₀. Throughout the paper, the turbulence quantities are given in computational units of ρ₀ = 1, c_s = 1, and L₀ = 1, unless otherwise specified. Regardless of forcings, v_rms ≈ 10 and 0.5 is attained at saturation in the ISM and ICM turbulence models, respectively.

In the ISM turbulence, v_rms as well as E_K and E_B grow until ∼ t_cross, and then, after some adjustments, reach saturation by ∼ 2 t_cross. For the models with the same M_turb, both E_K and E_B are smaller in the compressive forcing models than in the solenoidal forcing models. This is because the density distribution is more intermittent (e.g., Federrath et al. 2008, 2009) and the small-scale dynamo is less efficient (e.g., Federrath et al. 2011; Lim et al. 2020) with compressive forcing. The averaged values of β at saturation (2 t_cross ≤t ≤ t_end), β_sat, is listed in the seventh column of Table 1. It indicates that the magnetic energy grows by more than a factor of two in the solenoidal forcing models, while the growth factor is somewhat less than two in the compressive forcing models. The mixed forcing models show the behaviors in between.

In the ICM turbulence, reaching saturation takes longer in terms of t_cross, as mentioned in Section 2.2. In particular, initially with a very weak seed, the growth of the magnetic field needs a number of eddy turnovers, and hence E_B approaches saturation after t ∼ 15 t_cross. Moreover, the saturated E_B is much smaller in the compressive forcing models (green lines) than in the solenoidal forcing models (red lines), since the solenoidal component of flow motions, which is responsible for most of the magnetic field amplification, is smaller with compressive forcing, as shown before in Porter et al. (2015). The saturated E_B for the mixed forcing models (blue lines) is a bit smaller than, yet similar to, that for the solenoidal forcing models. As a consequence, β_sat at saturation (15 t_cross ≤ t ≤ t_end) is close to ∼50 in the solenoidal and mixed forcing models, while it is an order of magnitude larger in the compressive forcing models (see the seventh column of Table 1). Considering that β in the ICM is estimated to be ∼50–100 (e.g., Ryu et al. 2008; Brunetti & Jones 2014), the solenoidal and mixed forcing models would represent more realistic ICM turbulence. Another notable point is that E_B in the N512 models is somewhat larger than E_B in the N256 models. This tells us that the amplification of the magnetic field by small-scale dynamo is sensitive to the effective Reynolds and Prandtl numbers, which are controlled by numerical resolution in our simulations (see also, e.g., Schober et al. 2012; Roh et al. 2019). In our ICM-Sol-N512 model, E_B /E_K approaches ∼20% at saturation, while E_B /E_K ∼ 30% was obtained in a 2048³ simulation using the TVD code in Porter et al. (2015).

For the analysis of turbulent flows in the following subsections, we examine the mean quantities at saturation, which are calculated over 2 t_cross ≤ t ≤ t_end for the ISM models and 15 t_cross ≤ t ≤ t_end for the ICM models.

Figure 2 shows two-dimensional (2D) slice maps of the density (upper panels) and 3D distributions of the shock zones color-coded by M_s (lower panels) in the N512 ISM models with three different forcings at t = t_end. The density images clearly exhibit the characteristic morphologies with different forcings; the density has a larger contrast and a higher intermittency in the compressive forcing model (right panel) than in the solenoidal forcing model (left panel), which is consistent with previous studies using hydrodynamic simulations (e.g., Federrath et al. 2008, 2009). A complex network of shocks appears, regardless of forcings. While shocks are distributed relatively uniformly throughout the entire simulation box in the solenoidal forcing model, the distribution is more concentrated in the compressive forcing model. Again the mixed forcing model (middle panel) shows the behaviors in between. We note that the shock zones include both fast and slow shocks, while strong shocks with high M_s are mostly fast shocks (see the discussion below, and also Lehmann et al. 2016 and PR2019).

Zoom In Zoom Out Reset image size

Figure 3 shows the power spectra of the flow velocity and density, averaged over the saturated period, for the three N512 ISM models. In the lower panel, P_ρ (k) exhibits a clear dependence on forcing. P_ρ (k) is several to an order of magnitude larger in the compressive forcing model than in the solenoidal forcing model. P_ρ (k) is a bit larger in the mixed forcing model than in the solenoidal forcing model. The difference in P_ρ (k) should reflect the visual impression of the density-slice images in Figure 2. In addition, P_ρ (k) has slopes flatter than the Kolmogorov slope, −5/3, in all the models with different forcings, and this is a characteristic property of supersonic turbulence (e.g., Kim & Ryu 2005; Federrath et al. 2009).

Zoom In Zoom Out Reset image size

In the upper panel, the power spectra for the solenoidal and compressive components of the velocity, ${P}_{v}^{\mathrm{sol}}(k)$ and ${P}_{v}^{\mathrm{comp}}(k)$, are separately presented with solid and dotted lines. In contrast to P_ρ (k), the differences in P_v (k) with different forcings are not large. In particular, on small scales of k/k_inj ≳ a few, ${P}_{v}^{\mathrm{comp}}(k)$ as well as ${P}_{v}^{\mathrm{sol}}(k)$ are almost identical for the three cases of different forcings. This should be because the magnetic tension quickly converts compressive motions to Alfvén modes, and hence the solenoidal component of the velocity is efficiently generated even if the forcing is compressive. With the spatial extension of shock surfaces typically being much smaller than L₀ (see Figure 2), the frequency of shocks should be reflected mostly to ${P}_{v}^{\mathrm{comp}}(k)$ in the range of k/k_inj ≳ a few. With similar ${P}_{v}^{\mathrm{comp}}(k)$, below we see that the number of shock zones is similar in the models with different forcings. On the other hand, in k/k_inj ≲ a few, ${P}_{v}^{\mathrm{comp}}(k)$ is larger for the compressive (green dotted) forcing model than for the solenoidal (red dotted) and mixed (blue dotted) forcing models; this is consistent with the large-scale distributions of shocks shown in Figure 2.

Another point to note is that while ${P}_{v}^{\mathrm{comp}}(k)$ has the slope close to −2, which is the slope of the Burgers spectrum in shock-dominated flows (e.g., Kim & Ryu 2005; Federrath 2013), ${P}_{v}^{\mathrm{sol}}(k)$ is a bit flatter than the Kolmogorov spectrum. According to the Goldreich–Sridhar scaling, the Kolmogorov slope of −5/3 is expected for ${P}_{v}^{\mathrm{sol}}(k)$ in MHD turbulence, if the Alfvénic mode is dominant (Goldreich & Sridhar 1995). However, some simulations of "incompressible" MHD turbulence exhibited a slope close to −3/2 (e.g., Müller & Grappin 2005), although the slope obtained in numerical simulations is controversial (e.g., Beresnyak 2011). Our simulations of "compressible" MHD turbulence produce ${P}_{v}^{\mathrm{sol}}(k)$ with a slope close to ∼–1.2 in the inertial range, indicating that the compressiblility is possibly involved.

The statistics of shocks are presented in Figure 4. The upper panels plot the time-averaged PDFs of the sonic Mach number for fast and slow shocks, $({{dN}}_{\mathrm{fa}}/{{dM}}_{s})/{N}_{\mathrm{fa}}$ and $({{dN}}_{\mathrm{sl}}/{{dM}}_{s})/{N}_{\mathrm{sl}}$, averaged over the saturated period, for the six ISM models in Table 1. In all the cases of different forcings, the PDFs look similar: they peak at M_peak < M_turb and decrease more or less exponentially at higher Mach numbers. The Mach numbers for fast shocks are much higher than those for slow shocks, as expected. In the N512 models the PDFs for fast shocks have M_peak ≈ 3.5–4.5, and the characteristic Mach number M_char ∼ 6, if they are fitted to $\exp [-({M}_{s}-1)/({M}_{\mathrm{char}}-1)]$. By contrast, the PDFs for slow shocks have M_peak ≈ 1.1, close to the lowest Mach number we identify, and the characteristic Mach number M_char ∼ 1.7, although the distributions deviate from the exponential form at M_s ≳ 6. The majority, ∼80%–85%, of fast shocks have the sonic Mach number less than the turbulent Mach number, while virtually all slow shocks have M_s < M_turb, indicating that strong shocks with M_s > M_turb are relatively rare.

Zoom In Zoom Out Reset image size

The most interesting point in the PDFs is that the difference due to different forcings is not significant, in contrast to the ICM turbulence (see the next subsection). As noted above, with the total magnetic energy being comparable to the kinetic energy, ⁴ it should be the consequence of strong magnetic tension; the incompressible, solenoidal mode dominantly appears in the flow velocity, in particular, for k/k_inj ≳ a few, regardless of forcings. Hence, the compressive mode, which is responsible for the formation of shocks, is subdominant and about the same for the three different forcing models.

Accordingly, the energy dissipation at shocks also shows only a weak dependence on forcing. The lower panels of Figure 4 plot the dissipation rates of the turbulent energy at fast and slow shocks as a function of the sonic Mach number, normalized to the energy-injection rate, (d _fa/dM_s )/_inj and (d _sl/dM_s )/_inj, which are averaged over the saturated period. Similar to the Mach number PDF, the energy dissipation is dominated by shocks with small Mach numbers, while the contribution by high-Mach-number shocks decreases more or less exponentially. In the N512 models, (d _fa/dM_s )/_inj for fast shocks has peaks at M_peak ≈ 5.5 ∼ 7, and the characteristic Mach number M_char ∼ 7.5, if it is fitted to $\exp [-({M}_{s}-1)/({M}_{\mathrm{char}}-1)];$ (d _sl/dM_s )/_inj for slow shocks has peaks at M_peak ≈ 1.5, and the characteristic Mach number M_char ∼ 2.

The total numbers of fast- and slow-shock zones, N_fa and N_sl, normalized to the grid resolution, n_g ², are listed in the eighth and ninth columns of Table 1, while the total energy-dissipation rates at fast and slow shocks, integrated over the Mach number, _fa and _sl, are listed in the tenth and eleventh columns. N_fa is several times larger than N_sl in our results. This is partly because only the shock zones with M_s ≥ 1.05 are counted; then, while most of the fast shocks should be included, a substantial fraction of slow shocks might be missed since slow shocks could have even M_s < 1. If shocks with M_s < 1.05 are included, N_sl would be larger (see Lehmann et al. 2016 and PR2019). On the other hand, the contribution of slow shocks with M_s < 1.05 to the energy dissipation should not be significant. As a matter of fact, _fa is much larger than _sl, as also noted in PR2019.

The spatial frequency of shocks can be expressed in terms of the mean distance between shock surfaces, $\langle {d}_{\mathrm{shock}}\rangle \,={L}_{0}/({N}_{\mathrm{shock}}/{n}_{g}^{2})\propto {N}_{\mathrm{shock}}^{-1}$, where N_shock = N_fa + N_sl. For fast and slow shocks of M_s ≥ 1.05 altogether, the mean distance is estimated to be 〈d_shock〉 ∼ 0.3L_inj with L_inj = L₀/2 in the N512 models. The ratio of the energy dissipated at shocks and the injected energy, _shock/_inj, where _shock = _fa + _sl, is listed in the twelfth column of Table 1. It is ∼15% for the three different forcing models. This is roughly the fraction of the turbulent energy dissipated at shocks, while the rest, ∼85%, of the turbulent energy should dissipate through the turbulent cascade.

Although the dependence on forcing is not large, there are still some differences in the shock statistics with different forcings. For instance, both N_shock and _shock are larger in the solenoidal forcing model than in the compressive forcing model. This is partly because larger _inj should be adopted for the solenoidal forcing model to achieve the same M_turb ≈ 10 (see Table 1). We point out that smaller _inj with compressive forcing would be an unexpected result, since compressive motions are expected to dissipate faster and hence compressive forcing would require larger _inj. Indeed, we see larger _inj for the compressive case in the ICM turbulence with weak magnetic fields (Table 1). Again, this should be due to strong magnetic fields in the ISM turbulence. The magnetic tension seems to efficiently convert compressive motions to incompressive Alfvén modes, and hence _inj is not necessarily larger with compressive forcing.

Also, the shock statistics depend somewhat on the numerical resolution; N_shock is larger in the N512 models, while _shock is larger in the N256 models. However, considering uncertainties in the identification of shock zones and the calculations of M_s and Q, we regard the statistics as being reasonably resolution converged. Comparing the shock statistics of the ISM-Sol-N512 model to those of the similar model in PR2019, 1024M7-b0.1 (M_turb ≈ 7 and β₀ = 0.1), ISM-Sol-N512 has larger ${N}_{\mathrm{shock}}/{n}_{g}^{2}$ and _shock/_inj. This should be partly owing to the higher M_turb of ISM-Sol-N512, and also because the current WENO code with a higher order of accuracy seems to capture shocks better than the TVD code used in PR2019, particularly in complex flows with strong magnetic fields.

An interesting consequence of shocks is the density enhancement, which could have implications on the evolution of molecular clouds including the star formation rate (SFR). As shown in Figure 2, the density fluctuations are larger in the compressive forcing model; hence with more frequent occurrence of high-density regions, the SFR would be enhanced. As a matter of fact, Federrath & Klessen (2012) showed that the SFR would be about 10 times larger with compressive forcing than with solenoidal forcing in the turbulent ISM. We here examine how the forcing affects the density fluctuations at shocks as well as in the entire computational volume with the density PDF.

In isothermal turbulence, the density PDF is often approximated as the lognormal distribution (e.g., Vazquez-Semadeni 1994; Padoan et al. 1997; Passot & Vázquez-Semadeni 1998). The standard deviation of the density distribution, σ, is expected to be larger with compressive forcing than with solenoidal forcing. Federrath et al. (2008), for instance, showed that the radio of the compressive to solenoidal forcing cases is σ_comp/σ_sol ∼ 3, for hydrodynamic turbulence with M_turb ∼ 5. The value of σ depends on M_turb and the magnetic field strength, or β₀ (e.g., Federrath et al. 2008; Molina et al. 2012); so does this ratio. The top panel of Figure 5 shows the PDFs of $\mathrm{log}\rho$ in the entire computational volume, averaged over the saturated period, for the N512 ISM models with three different forcing modes. ⁵ In our simulations, σ_comp/σ_sol ≈ 2.1; σ_mix for mixed forcing is similar to σ_sol. With a larger M_turb ≈ 10, the smaller ratio should be caused by strong magnetic fields.

Zoom In Zoom Out Reset image size

The middle and bottom panels of Figure 5 show the PDFs of $\mathrm{log}\rho$ at the preshock (dotted lines) and postshock (solid lines) regions, for fast and slow shocks, respectively. Again, compressive forcing leads to broader density distributions with larger σ's in both the preshock and postshock regions than solenoidal forcing. In all the cases, σ_comp/σ_sol is in the of range ≈2.2–2.9 at both the preshock and postshock regions, and again σ_mix is similar to σ_sol. These values around the shocks are comparable to, or slightly larger than, σ_comp/σ_sol for the entire volume.

Figure 6 shows 2D slice maps of the density (upper panels) and the color-coded M_s (lower panels) in the N512 ICM models with different forcings at t = t_end. The density images exhibit expected distributions with different forcings, for instance a higher intermittency in the compressive forcing model (right panel) than in the solenoidal forcing model (left panel). With β_sat ≫ 1, i.e., subdominant magnetic fields at saturation, the turbulence should be almost hydrodynamic; hence the density distributions are similar to those of previous hydrodynamic studies, such as Federrath et al. (2008, 2009), although the details would be different if M_turb was different. Our ICM turbulence model is subsonic, with M_turb ≈ 0.5, yet shocks are ubiquitous (see also Porter et al. 2015). Similar to the density distributions, the shock distributions exhibit differences with different forcings; the distribution looks more organized with stronger shocks in the compressive forcing model than in the solenoidal forcing model. In the mixed forcing model (middle panel), the density and shock distributions show a slightly larger intermittency and more shocks, compared to the solenoidal forcing model.

Zoom In Zoom Out Reset image size

Figure 7 shows the power spectra for the solenoidal and compressive components of the flow velocity, ${P}_{v}^{\mathrm{sol}}(k)$ and ${P}_{v}^{\mathrm{comp}}(k)$ (upper panel), and the density, P_ρ (k) (lower panel), averaged over the saturated period, for the three N512 ICM models. As in the ISM case, P_ρ (k) is several to an order of magnitude larger in the compressive forcing model than in the other forcing models; P_ρ (k) is a bit larger in the mixed forcing model than in the solenoidal forcing model. Even though there are discontinuities in the density distribution formed by shocks, in the ICM turbulence models with small M_turb, P_ρ (k) is steeper than the Kolmogorov spectrum in all the models with different forcings. This is consistent with the previous finding of Kim & Ryu (2005), that in hydrodynamic turbulence, while P_ρ (k) flattens as M_turb increases, the slope is less than −5/3 for M_turb ≲ 1.

Zoom In Zoom Out Reset image size

The upper panel shows that P_v behaves differently from that of the ISM turbulence. While ${P}_{v}^{\mathrm{sol}}(k)$ dominates over ${P}_{v}^{\mathrm{comp}}(k)$ in all the wavenumbers in the solenoidal and mixed forcing models, ${P}_{v}^{\mathrm{comp}}(k)$ is much larger than ${P}_{v}^{\mathrm{sol}}(k)$ in the compressive forcing model, as was previously shown, for instance, in Federrath et al. (2011) and Porter et al. (2015). With negligible magnetic tension in the high-β ICM, the memory of forcing is persistent in the flows of fully developed turbulence. ${P}_{v}^{\mathrm{comp}}(k)$ is several times larger in the compressive forcing model than in the other forcing models, and P_comp(k) is a bit larger in the mixed forcing model than in the solenoidal forcing model. With larger ${P}_{v}^{\mathrm{comp}}(k)$, shocks would be more abundant in the compressive forcing model, as noted with the spatial distribution of shocks in Figure 6. For all the models with different forcings, ${P}_{v}^{\mathrm{comp}}(k)$ is steeper than the Kolmogorov spectrum and has the slope close to −2, which is the slope of the shock-dominated Burgers spectrum. By contrast, ${P}_{v}^{\mathrm{sol}}(k)$ shows a slight concavity around k ∼ 10–30 in the solenoidal and mixed forcing models. It has been argued that in ICM turbulence, the power spectrum for the kinetic energy also has a concavity but, being compensated by the magnetic power, the power spectrum for the total energy roughly follows the Kolmogorov spectrum (see, e.g., Porter et al. 2015). Hence, while the magnetic field is subdominant, it would still affect the flow motions and P_v (k), especially in those solenoidal and mixed forcing models.

Unlike in the ISM turbulence models, virtually all the shocks formed in the ICM turbulence models are fast shocks, and the fraction of slow shocks is very small; for instance, ∼1% in the ICM-Sol-512 model and even smaller in the compressive and mixed forcing models. This is because the magnetic field strength decreases across slow shocks, and hence they can appear only when the preshock magnetic field is sufficiently strong to satisfy the condition ${v}_{A\parallel 1}^{2}/{c}_{s}^{2}\gt {M}_{s}^{2}$. With the small fraction and almost no contribution to the energy dissipation, below we do not further consider slow shocks, and present only the statistics of fast shocks for the ICM models.

Figure 8 presents the statistics of the fast shocks. The upper panel shows the PDFs of the sonic Mach number, $({{dN}}_{\mathrm{fa}}/{{dM}}_{s})/{N}_{\mathrm{fa}}$, averaged over the saturated period, for the six ICM models in Table 1. As in the ISM cases, the PDFs peak at small Mach numbers and decrease more or less exponentially at high Mach numbers. However, unlike in the ISM cases, the PDFs differ substantially with different forcings; there are more shocks with higher Mach numbers in the compressive forcing model (green lines). With M_turb < 1, the peak occurs at M_peak ≈ 1.05, the lowest Mach number we identify, in all the forcing models. But when the PDFs are fitted to $\exp [-({M}_{s}-1)/({M}_{\mathrm{char}}-1)]$, the characteristic Mach numbers are M_char ∼ 1.07, 1.08, and 1.13 for the solenoidal, mixed, and compressive forcing models, respectively. This reveals that with compressive forcing, shocks with higher compression and higher M_s are more abundant. We note that these characteristic Mach numbers agree well with those of Porter et al. (2015), who quoted M_char ∼ 1.08 and 1.125 for the solenoidal and compressive forcing cases, although different numerical codes and different shock-identification schemes are used.

Zoom In Zoom Out Reset image size

The lower panel of Figure 8 plots the dissipation rate of the turbulent energy at fast shocks normalized to the energy-injection rate, (d _fa/dM_s )/_inj, which are averaged over the saturated period. Again the distributions peak at small Mach numbers and decrease more or less exponentially at high Mach numbers in all the cases of different forcings, whereas it shifts to higher Mach numbers in the compressive forcing model. The peak and characteristic Mach numbers are M_peak ≈ 1.10, 1.14, and 1.33, and M_char ∼ 1.11, 1.12, and 1.16, for the solenoidal, mixed, and compressive forcing models, respectively. ⁶

The total number of fast-shock zones, N_fa, normalized to the grid resolution, n_g ², and the total energy-dissipation rate at fast shocks, integrated over the Mach number, _fa, are given in Table 1. As expected from the PDFs in Figure 8, N_fa is substantially larger in the compressive forcing model than in the other forcing models, and also N_fa is noticeably larger in the mixed forcing model than in the solenoidal model. Specifically, the N_fa of ICM-Comp-N512 and ICM-Mix-N512 is ∼1.8 and ∼1.3 times the N_fa of ICM-Sol-N512. The mean distance between shock surfaces is 〈d_shock〉 ∼ 0.31, 0.24, and 0.17 L_inj in the N512 ICM models with solenoidal, mixed, and compressive forcings, respectively. An interesting point is that fast shocks are more frequent in our ICM turbulence models than in the ISM models, although M_turb is 20 times smaller. While strong magnetic fields in the ISM models limit the gas compression and hence inhibit the formation of shocks, subdominant magnetic fields in the ICM models do not significantly affect the occurrence of shocks.

Accordingly, _fa is much larger in the compressive forcing model than in the other forcing models. For instance, the _fa of ICM-Comp-N512 is ∼7.5 times the _fa of ICM-Sol-N512. On the other hand, the energy-injection rate, _inj, differs substantially with different forcings, and the _inj of ICM-Comp-N512 is ∼2.4 times the _inj of ICM-Sol-N512. As a result, the normalized energy-dissipation rate, _fa/_inj, is only ∼3.1 times larger in ICM-Comp-N512 than in ICM-Sol-N512, and _fa/_inj is ∼1.8 times larger in ICM-Mix-N512 than in ICM-Sol-N512. The fraction of the turbulent energy dissipated at shocks is estimated to be _shock/_inj ∼ 11%, ∼ 19%, and ∼ 33% in the N512 ICM models with solenoidal, mixed, and compressive forcings, respectively. Compared to the fractions in the ISM turbulence, the fraction in the solenoidal forcing model is slightly smaller, but the fractions in the compressive and mixed forcing models are larger.

The dependence of the shock statistics on numerical resolution is rather weak in the ICM turbulence models, as can be seen in Figure 8 and Table 1. For instance, ${N}_{\mathrm{fa}}/{n}_{g}^{2}$ and _shock/_inj are ∼ 6% larger in ICM-Sol-N512 than in ICM-Sol-N256, and the differences are even smaller in the other forcing models. Again, considering uncertainties in the estimation of those quantities, the statistics would be regarded as being reasonably resolution converged.

The shock statistics of the ICM-Sol-N512 model may be compared to those of 1024M0.5-b10 in PR2019 (M_turb ≈ 0.5 and β₀ = 10). The shock frequency, ${N}_{\mathrm{fa}}/{n}_{g}^{2}$, is larger in ICM-Sol-N512 than in 1024M0.5-b10, partly because shocks with M_s ≥ 1.05 are counted in this work, while those with M_s ≥ 1.06 are included in PR2019. Also in ICM-Sol-N512 with a much weaker initial magnetic field (β₀ = 10⁶), β_sat is about 10 times larger than in 1024M0.5-b10, that is, the magnetic energy at saturation is about 10 times smaller, and hence more shocks form. However, the shock dissipation fraction, _shock/_inj, is actually smaller in ICM-Sol-N512. This is because the term involving v_A in Equation (9), which represents the energy dissipation through the magnetic field, makes a relatively small contribution to _shock in ICM-Sol-N512, while it is substantial in 1024M0.5-b10, especially at weak shocks.

The upper panel of Figure 9 displays the PDFs of $\mathrm{log}\rho$ in the entire computational volume, averaged over the saturated period, for the three N512 ICM models with different forcings. As in the ISM cases, compressive forcing leads to a larger density contrast with larger σ than solenoidal forcing. In our simulations, σ_comp/σ_sol ≈ 3.1; this value is roughly the same as that of Federrath et al. (2008) for hydrodynamic turbulence with M_turb ∼ 5, indicating σ_comp/σ_sol would not be very sensitive to M_turb, as long as the magnetic field is subdominant. The width of the PDF in the mixed forcing model is somewhat larger than in the solenoidal forcing model, with σ_mix/σ_sol ≈ 1.4. The bottom panel plots the PDFs of $\mathrm{log}\rho$ at the preshock and postshock regions of fast shocks. While the peaks of the PDFs are located at low and high densities for the preshock and postshock regions, respectively, the widths look similar to those for the whole computational volume. Quantitatively, σ_comp/σ_sol ≈ 3.0 and 2.9 for the preshock and postshock regions, respectively, and σ_mix/σ_sol ≈ 1.3 for both the preshock and postshock regions.

Zoom In Zoom Out Reset image size