SHOCKFIND  An algorithm to identify magnetohydrodynamic shock waves in turbulent clouds
Abstract
The formation of stars occurs in the dense molecular cloud phase of the interstellar medium. Observations and numerical simulations of molecular clouds have shown that supersonic magnetised turbulence plays a key role for the formation of stars. Simulations have also shown that a large fraction of the turbulent energy dissipates in shock waves. The three families of MHD shocks — fast, intermediate and slow — distinctly compress and heat up the molecular gas, and so provide an important probe of the physical conditions within a turbulent cloud. Here we introduce the publicly available algorithm, shockfind, to extract and characterise the mixture of shock families in MHD turbulence. The algorithm is applied to a 3dimensional simulation of a magnetised turbulent molecular cloud, and we find that both fast and slow MHD shocks are present in the simulation. We give the first prediction of the mixture of turbulencedriven MHD shock families in this molecular cloud, and present their distinct distributions of sonic and Alfvénic Mach numbers. Using subgrid onedimensional models of MHD shocks we estimate that % of the volume of a typical molecular cloud in the Milky Way will be shock heated above 50 K, at any time during the lifetime of the cloud. We discuss the impact of this shock heating on the dynamical evolution of molecular clouds.
keywords:
(magnetohydrodynamics) MHD – shock waves – turbulence – ISM: kinematics and dynamics – ISM: clouds1 Introduction
The formation of stars is sensitive to the underlying physics of the supersonic magnetohydrodynamic (MHD) turbulence in molecular clouds (Mac Low & Klessen, 2004; Elmegreen & Scalo, 2004; Krumholz & McKee, 2005; McKee & Ostriker, 2007; Hennebelle & Chabrier, 2008; Padoan & Nordlund, 2011; Hennebelle & Falgarone, 2012). High resolution threedimensional simulations of molecular clouds have shown that the star formation rate and efficiency depend on whether the turbulence is solenoidally or compressively driven (Federrath & Klessen, 2012, 2013), and that the stellar initial mass function is sensitive both to nonideal MHD effects such as ambipolar diffusion (McKee et al., 2010) and to the driving and Mach number of the turbulence (Hennebelle & Chabrier, 2009, 2013; Hopkins, 2013). The importance of stellar feedback (Krumholz et al., 2007; Cunningham et al., 2011; Myers et al., 2014; Nakamura & Li, 2007; Wang et al., 2010; Price & Bate, 2008, 2009; Offner & Arce, 2014; Federrath et al., 2014; Federrath, 2015; Padoan et al., 2015), whether gravity drives turbulent motions (Elmegreen & Burkert, 2010; Klessen & Hennebelle, 2010; VÃ¡zquezSemadeni et al., 2010; Federrath et al., 2011; Robertson & Goldreich, 2012) and the role that turbulence plays in producing the ubiquitously observed filaments (Arzoumanian et al., 2011; André et al., 2014; Smith et al., 2014, 2016; Federrath, 2016; Kainulainen et al., 2016; Hacar et al., 2016) are big questions that continue to be studied. Rigorous observational effects distinctly revealing the presence or dominance of the various physical processes are strongly sought after.
Observed nonthermal linewidths of molecular lines reveal a turbulence in molecular clouds that is highly supersonic, with Mach numbers typically in the range of 3 to 30 (Larson, 1981; Solomon et al., 1987; Ossenkopf & Mac Low, 2002; Heyer & Brunt, 2004; RomanDuval et al., 2011; Schneider et al., 2013; Henshaw et al., 2016). The supersonic flows in these clouds will inevitably form shock waves. MHD simulations by Stone et al. (1998) found that of the turbulent energy is dissipated to shocks. When turbulence is allowed to decay, a large range of weak shocks is responsible for the majority of the dissipation (Smith et al., 2000), whereas a small range of stronger shocks dissipates the turbulence while it is continuously driven (Smith et al., 2000). If other turbulence parameters — such as magnetic field strength, driving mechanism, Mach number, inclusion of other physical effects like selfgravity, stellar feedback, cosmic rays and others — could be linked to distributions of shock waves, then the radiative signatures of shocks become observational diagnostics for these parameters in molecular clouds. The goal of this work is to determine this link, from turbulence parameters to distributions of shock waves, in simulations of turbulent clouds in order to provide observational tests of the various physical processes.
MHD fluids can support three families of shocks: fast, intermediate and slow (De Hoffmann & Teller, 1950; Kennel et al., 1989). Note that these names refer to the associated MHD linear wave modes, and not to the speed of a given shock. In Section 2 we use the MHD jump conditions to detail the fundamental differences between the shock families. Pon et al. (2012) argue that, for molecular clouds, the turbulent cascade leads to lowvelocity shocks (a few km/s) doing the majority of the dissipation. They show that even at low velocities, fast MHD shocks will radiate more strongly than photodissociation regions in mid rotational transitions of CO. Lehmann & Wardle (2016) use twofluid MHD models of shocks in molecular cloud conditions to show that lowvelocity fast and slow MHD shocks distinctly compress and heat the molecular gas. At velocities less than 4 km/s slow shocks can reach compression ratios of up to 500 whereas fast shocks compress the gas less than 10 times the preshock values. In addition, slow shocks can reach peak temperatures up to K whereas fast shocks reach up to K. This is because in a weakly ionized gas — such as in molecular clouds — the ionneutral collision timescale, which determines the heating timescale in fast shocks, is slower than the cooling timescale, which in turn is slower than the neutralion collision timescale that controls the heating in slow shocks. These higher peak temperatures in slow shocks result in stronger CO rotational lines (above ) and lowlying pure rotational lines of H () than in fast shocks of the same velocity.
Molecular clouds are usually assumed to be in chemical equilibrium for the average conditions of the cloud. However, the inhomogeneous structure introduced by turbulence allows for local reaction rates to significantly differ from global average conditions (Hollenbach et al., 1971; Wolfire et al., 1995; Glover et al., 2010). Shocks are an important local mechanism for driving chemical reactions. Kumar & Fisher (2013) followed the chemical evolution of a parcel of gas through a turbulence simulation and found the abundances of molecules, like CH and HCO, to be highly sensitive to shocked regions. Their work considers only hydrodynamic shocks and so cannot capture the qualitatively distinct behaviours of differing MHD shock families. For example, the higher temperatures of lowvelocity slow shocks produce vastly different chemical abundances than fast shocks of the same velocities (Lehmann & Wardle, 2016). Using a simple oxygen chemical network Lehmann & Wardle show that while fast shocks leave preshock abundances mostly untouched, slow shocks can increase the abundances of molecules like OH, O and HO by several orders of magnitude within the hot shock front. To understand the chemical makeup of turbulent molecular clouds it is therefore important to understand which families of MHD shocks are present and how much volume they affect.
While there is some observational evidence for the dissipation of molecular cloud turbulence in fast MHD shocks (Lesaffre et al., 2013; Pon et al., 2014; Pon et al., 2016; Larson et al., 2015) and in slow MHD shocks (Lehmann & Wardle, 2016), thus far no one has determined which kinds of MHD shocks are present in simulations of magnetised turbulent molecular clouds. In Section 3 we present an algorithm, shockfind, to detect and characterise the mixture of MHD shock types in such simulations. We apply our new algorithm to an MHD simulation of molecular cloud turbulence in Section 4. Finally, we discuss these results and conclude our study in Sections 5 and 6 respectively. The shockfind algorithm, written in python, is publicly available and can be found on BitBucket (https://bitbucket.org/shockfind/shockfind) and the python Package Index (https://pypi.python.org/pypi/shockfind). It is released under the Apache license version 2.0, and comes with documentation including a user’s guide.
2 Magnetohydrodynamic Shocks
Here we derive some fundamental differences between the shock families in the ideal limit of MHD. We summarise the relevant theoretical implications of the MHD jump conditions that can be found in Lehmann & Wardle (2016), and illustrate how the various families of MHD shocks characteristically affect the ambient magnetic field.
In the frame of reference comoving with a shock wave — the shock frame — the preshock fluid has a speed, , greater than a linear wave speed in the fluid. The fluid then transitions inside a discontinuity to a fluid velocity less than a wave speed in the post shock fluid. In ideal MHD, the three linear waves supported, the fast, intermediate and slow waves have phase velocities
where is the Alfvén velocity, is the isothermal sound speed with Boltzmann constant and mean mass per particle , and is the angle between the magnetic field and the direction of propagation of the wave. These speeds are plotted as functions of in Fig. 1 for , as is usually the case in molecular clouds.
The three wave speeds demarcate four regions of fluid velocities marked 1 to 4 in Fig. 1. There are six ways to transition across at least one wave speed within the shock front, resulting in three families of MHD shocks: fast, intermediate and slow shocks. Fast shocks cross the fast wave speed only (1–2), intermediate shocks cross the intermediate wave speed (1–3, 1–4, 2–3, and 2–4), and slow shocks cross the slow wave speed only (3–4). This means that for a fast shock the Alfvénic Mach number — defined by — is necessarily greater than unity, whereas for a slow shock the shock speed is subAlfvénic, i.e., . We will use this criterion to distinguish fast from slow shocks in the next section.
For a timeindependent, planeparallel shock wave, the pre and postshock mass density, velocity, gas pressure and magnetic field are related by the RankineHugoniot jump conditions (Kennel et al., 1989). A consequence of these jump conditions are changes in magnetic field geometry across a shock front characteristic of each shock family. In fast shocks, the component of the magnetic field perpendicular to the direction of propagation, , must increase from the preshock to the postshock value. In intermediate shocks, must switch sign. This switch is due to a rotation of the magnetic field within the shock front. Finally, in slow shocks must decrease. For all shocks, the planar symmetry implies that the magnetic field parallel to the direction of propagation is conserved across the shock front. These three characteristic changes of the magnetic field direction are shown schematically in Fig. 2.
The magnetic field strength is proportional to the separation of field lines, so one can see from Fig. 2 that the field strength increases across fast shocks and decreases across slow shocks. We will use this fact as a signature of these two classes of shocks in Section 3. This also means that some of the kinetic energy of a fast shock is converted into magnetic field energy. Hence, for slow shocks at the same velocity as fast shocks, a greater portion of the energy budget is available to heat the gas. As we do not account for the magnetic field geometry, but rather use only the field strength as an indicator of shock type we do not explicitly capture switchon or switchoff shocks. However, these shocks are special cases of fast and slow shocks and so are counted amongst these categories. Finally, the magnetic field strength in intermediate shocks can either increase or decrease across the shock front depending on the initial conditions of a particular shock. Hence there is no simple signature of intermediate shocks in the magnetic field strength. In addition, there has been debate over whether intermediate shocks are physically admissable (e.g., Wu, 1987; Falle & Komissarov, 2001), and so this class of MHD shocks is not considered further in this work.
The structure of the magnetic field across the different shock waves only depends on the ideal MHD jump conditions. The heating inside the shock front, however, can depend on nonideal effects. For example, in molecular clouds the gas is weakly ionized and so ion and neutral species can be decoupled. In this case, a multifluid approach is necessary to model the structure of the shock front. In fast shocks, the strong magnetic pressure behind the shock front drives ion species ahead of the neutrals in a magnetic precursor (Mullan, 1971; Draine, 1980). Collisional heating in fast shocks is thus controlled by the ionneutral collision timescale, which is generally larger than the cooling timescale for lowvelocity shocks (a few km/s) in molecular clouds. If the cooling keeps the temperature low within the shock, the neutral velocity may remain supersonic throughout and the fluid variables will smoothly transition from preshock to postshock values in what is called a Ctype shock. In fast shocks with shock velocity exceeding the magnetosonic speed of the charged fluid, defined by
where is the Alfvén velocity in the charged fluid, a thin jump will form in which the heating is determined by molecular viscosity. Such a fast shock is called Jtype. We see then, that we need to determine the local preshock conditions in turbulent molecular clouds in order to predict the shock heating from fast MHD shocks.
In the twofluid models of Lehmann & Wardle (2016) slow shocks produce a temperature structure distinct from fast shocks. Shocks of this family are driven by the gas pressure of the neutrals, and so heating is determined by the neutralion collision timescale. This timescale is shorter than the cooling timecale and so high peak temperatures are reached in a thin shock front resembling an ordinary hydrodynamic shock. The jump in temperature across this shock front is determined by the sonic Mach number :
(1) 
which reaches 200 K at low shock velocities (2 km/s for km/s) and adiabatic index . This is hot enough to produce chemical abundances and molecular cooling significantly different to ambient molecular cloud conditions.
Fast and slow shocks therefore present distinct temperature structures within molecular clouds. To fully understand the heating and chemistry driven by turbulent dissipation in shock waves it is therefore critical to determine the relative fraction of shock families. We will do this in the following by using numerical MHD simulations of turbulent molecular clouds.
3 Shock Detection Algorithm
In this section, we present a new algorithm, shockfind, to detect and characterise the shocks in MHD simulations of molecular clouds. We then test this algorithm by considering a case study simulation of colliding MHD shocks, which should be a common occurence in supersonic MHD turbulence.
3.1 Algorithm Summary
Here we summarise the sevenstep algorithm, shockfind, for detecting fast or slow shocks in an MHD simulation. For a simulation with mass density , three velocity components , , and , and three magnetic field components , , and , the algorithm will:

Identify shock candidates as computational cells with large convergence:
(2) or large magnitude of the density gradient:
(3) where large is above a userdefined threshold appropriate to the particular simulation;

Compute the shock direction at the location of each candidate cell, , using the gradient of the density:
(4) 
Extract averaged fluid variables along a cylinder perpendicular to the shock front. The cells, at coordinates , on the central axis of the cylinder are defined by
(5) where is the location of the candidate cell and parametrises the line and ranges from a few shock thicknesses, , that the simulation spreads the shock over. The effect of varying the radius of this cylinder is shown in Appendix A;

The cellaveraged variables for are the preshock values, while the cellaveraged variables for represent the postshock values;

Compute the shock speed using the pre and postshock parallel velocities and densities obtained in step (iv):
(6) where ;

Compare the shock speed to the preshock Alfvén velocity to form the Alfvénic Mach number . Fast shocks must have and slow shocks must have . In addition, we compare the pre and postshock magnetic field strengths. Fast shocks must have and slow shocks must have . If a shock candidate consistently satisfies both of these inequalities then we have detected a fast or slow MHD shock, otherwise it is only a candidate and is not included in the further analysis;

Finally, we filter the detected shock cells by ignoring those detections that do not occur at local maxima in the convergence along the extracted line of step (iii). This step avoids extracting multiple cells as individual shocks that actually belong to the same shock. This process is illustrated in Appendix B.
3.2 Test simulation of colliding MHD shocks
Variable  Slow Postshock (1)  Common Preshock (2)  Fast Postshock (3) 

In this section we illustrate the steps of our shock detection algorithm outlined in Section 3.1, by applying it to a simulation of colliding shock waves. We use the code FLASH (Fryxell et al., 2000; Dubey et al., 2008) in version 4 to integrate the ideal, threedimensional MHD equations. The equations are solved on a grid with a total of grid cells using the HLL3R positivedefinite Riemann solver (Waagan et al., 2011). The equations are closed with an ideal gas equation of state with adiabatic index .
In this simulation we initialise a box with a slow MHD shock with shock speed km/s travelling in the direction , and a fast MHD shock with shock speed km/s travelling in the direction . We choose preshock variables (common to both shocks) of density , pressure , and a magnetic field strength of 35 G oriented at to the slow shock front. These preshock values are used to compute postshock values using the MHD jump conditions (Kennel et al., 1989), and we choose a frame of reference such that the fast shock is stationary. Fig. 3 shows the initial density configuration of the three regions — slow postshock, common preshock, and fast postshock — and the fluid variables are listed in Table 1.
3.2.1 Convergence and density gradient
In Fig. 4 we plot slices of mass density of the colliding shock simulation at three times in the simulation. The rightmost column shows a slice after the two shocks have collided and a significant interaction region has developed. The coloured contours show the density with 4 regions marked A–D. Region A is the postshock region of the initial 5 km/s fast shock, region B is the common preshock region, region C is the postshock region of the initial 1 km/s slow shock, and finally region D is the postinteraction region.
Overplotting the density contours in the upper row of Fig. 4 are white line contours of strong convergence, defined by equation (2). This quantity distinctly picks out candidate shock fronts labelled S1S4. S1 and S2 are the initial fast and slow shocks we set up to collide. S3 and S4 are the results of the collision, which have geometries suggestive of a refractive (S3) and reflective (S4) process.
In the lower row of Fig. 4 we plot projected vectors of the gradient of the density over the contours of density. This vector points in the direction of increasing density and so it always points towards the plane of a shock front. This allows us to define a line through a shock at which we extract the fluid variables. It also allows us to compute the fluid velocity in the direction of shock propagation, , by projecting onto the direction given by the gradient.
3.2.2 Shock family criteria
Using the convergence and gradient as described above, we plot in Fig. 5 extracted fluid variables through the shock candidates S1S4 from Fig. 4. The axis is zeroed at the convergence peak, with the postshock region at negative values of and preshock region at positive values of . The convergence and gradient are normalised to their peak values. To compare pre and postshock fluid variables we take the average of the variable over a few cells either side of the convergence peak.
The S1 and S2 shocks (which were the initial fast and slow shocks, respectively) show strong density contrasts though the simulation has spread out the S2 shock (red lines in Fig. 5) over a wider range than the S1 shock. This is also reflected in a much lower and wider convergence peak in the slow shock. As discussed in Section 2 the magnetic field, the lower panel of Fig. 5, in the S1 shock is stronger in its postshock region compared to its preshock region because it is a fast MHD shock. Conversely, in the S2 shock is stronger in the preshock region because it is a slow MHD shock.
Once we have computed pre and postshock densities, and , respectively, and the pre and postshock parallel velocities, and , respectively, we can determine the shock velocity using equation (6). By comparing the shock speed to the Alfvén velocity we form the Alfvénic Mach number . As discussed in Section 2, for fast shocks and for slow shocks . For the S1 and S2 candidates, we find that and respectively, confirming their status as fast and slow MHD shocks.
The S3 and S4 shock candidates both have magnetic field stronger in the postshock region than in the preshock region, indicating that they are fast MHD shocks. Indeed, for S3 the Alfvénic Mach number and so it satisfies both criteria. However, for the S4 candidate , ruling it out as a fast shock. As it is a planar converging structure it could be a fast linear wave or even an intermediate shock, but our simple criteria cannot distinguish these cases. We will call converging structures that show magnetic field ratios inconsistent with Alfvénic Mach number criteria fastlike and slowlike disturbances and exclude them from further analyses.
4 Molecular Cloud Turbulence
Here we present an application of our shockdetection algorithm to a 3D simulation of molecular cloud turbulence. The ultimate goal is to determine the fraction of slow and fast shocks and their respective effects for the heating and evolution of molecular clouds.
We use the turbulent initial conditions of molecular cloud simulation model 21 (GT256mM10B10) from Federrath & Klessen (2012). It is a 3D simulation of an isothermal, ideal MHD, turbulent molecular cloud with mixed compressive and solenoidal driving. The turbulence is driven to maintain a velocity dispersion km/s. As the sound speed , the turbulence contains a large fraction of supersonic gas and is therefore expected to drive shocks with Mach numbers . The initial magnetic field strength is 10 G, and the time step that we analyse here occurs after the turbulence has been fully developed but before selfgravity has been turned on to study star formation. The details of the integration scheme can be found in that paper.
Fig. 6 shows threedimensional renderings of the mass density (, left panel) and convergence (, right panel) of the simulation cube. The mass density is cut off at the average value , so that we plot only the high density regions. These regions are highly filamentary and fill only a small volume of the cloud. The convergence has been normalised by where is the 3D velocity dispersion of the cloud and is the cell size. There is a rough correlation between regions of strong convergence and regions of high density. This would be expected if the highest densities in turbulent clouds are postshock layers.
4.1 Search Thresholds
While shockfind could check for shocks at every cell in the simulation, it would be computationally expensive to check all cells of this simulation. Considering that nonconverging cells can be ruled out as shock candidates without further analysis, we develop search criteria to speed up the process.
Fig. 7 shows the probability distribution functions (PDFs) of the convergence (upper) and the magnitude of the gradient of density (lower). The magnitude of the gradient of density has been normalised by where is the average mass density. We can estimate the convergence and gradient across a shock wave with shock velocity propagating into a gas with preshock mass density as
where is the postshock velocity, is the postshock mass density and is the number of cells the simulation typically spreads a discontinuity over. We use N=3 in this work. Using the relation and normalising as above, these estimates become
where is the compression ratio. We use the compression ratio to control the search thresholds of step (i) of the algorithm (Section 3.1). In Fig. 7 the dashed vertical lines show the thresholds for shock velocity km/s, compression ratio , and preshock density . As these thresholds are treated independently (Step (i) of Section 3.1), cells that do not satisfy one threshold may still be identified as a shock candidate if they satisfy the other threhold. This conservative approach means we look at more cells than if we applied both thresholds simultaneously. At a velocity of 1 km/s, a slow shock will reach a peak temperature of 55 K (equation (1)). Molecular cooling is more efficient at high temperatures and high densities, and so running the algorithm at cells above these thresholds ensures we extract the most observationally relevant shocks.
4.2 Shock family statistics
Using the thresholds defined in the previous section, the spatial distribution of fast and slow shocked cells is shown in Fig. 8. Red points refer to slow shocks and blue points refer to fast shocks. Many of the detected cells form connected shock front sheets and others form long filamentary structures. Around 40% of searched cells fail to consistently satisfy the Alfvénic Mach number and magnetic field ratio criteria (Step (vi) of Section 3.1). A further half of the detected shocked cells are filtered out because they do not lie on local maxima in convergence (Step (vii) of Section 3.1). We analyse the results of this search in the following sections.
Fig. 9 shows the distributions of sonic Mach numbers for fast shocks (blue) and slow shocks (red). The slow shock distribution steeply and monotonically decreases, with a larger number of slow shocks than fast below ( km/s). The fast shock distribution peaks around before slowly decreasing. We estimate the area occupied by the shock fronts by treating each detected shocked cells as having an area of one of its faces: pc. We perform a convergence test shown in Appendix C. While there are more very low velocity shocks () to be found below our thresholds defined in Section 4.1, the Mach number distributions are well converged for the most observationally relevant shocks. Fig. 10 shows the distributions of Alfvénic Mach numbers for fast shocks (blue) and slow shocks (red). By definition slow shocks are subAlfvénic and fast shocks are superAlfvénic and so the distributions distinctly lie on either side of unity. In the twofluid MHD shocks of Lehmann & Wardle (2016) the peak temperature of slow shocks is determined by the sonic Mach number, whereas for fast shocks it is determined by the competition of molecular cooling and ionneutral collisional heating. In Section 5.2 we use these two distributions (Figs. 9 and 10) to make an estimate of turbulencedriven shock heating.
In Fig. 11 we plot the distributions of preshock magnetic field strengths and mass densities for the search described above. These distributions give us the typical preshock variables that should be used to model shocks relevent to molecular cloud turbulence. They show that the typical preshock conditions are different for fast and slow shocks. For example, for fast shocks the average preshock mass density g/cm whereas for slow shocks g/cm. This corresponds to total hydrogen densities, , of cm and cm, respectively (using ). In addition, for fast shocks the average preshock magnetic field strength G, whereas for slow shocks G. The average preshock density and magnetic field strength for fast shocks are within the ranges used by Pon et al. (2012) to model twofluid Ctype fast shocks. While the average density jump is much higher in slow shocks, , than in fast shocks, , the average postshock densities are remarkably similiar: g/cm. This implies that both kinds of shocks are equally important with respect to star formation because it is the postshock gas that sets the initial conditions for densecore and star formation (Padoan & Nordlund, 2011; Federrath & Klessen, 2012; Padoan et al., 2014).
These distributions suggest that to understand the impact that MHD shocks have on molecular clouds we need to understand both fast and slow MHD shocks. In a future paper we will apply this algorithm to other MHD simulations of molecular cloud turbulence (e.g., Federrath & Klessen, 2012; Federrath, 2015). We will investigate how the mixture of shock families may depend on the parameters of turbulence, e.g. the initial magnetic field strength, inclusion of extra physics such as selfgravity or protostellar jet feedback. If the mixture of shock types proves to be sensitive to these parameters, then differences in observational signatures between the shock types become signatures of these parameters. In Section 5.2 we discuss how one can use the shock mixture to compute the filling factor of hot, dense shocked gas.
4.3 Energetics
We may also consider the energetics of the shocks by comparing the kinetic energy dissipated in the shocks to the energy available in the turbulent motions. The kinetic flux through a unit area of shockfront is
and the turbulent energy density is
where is the velocity dispersion. Thus the timescale for dissipation in turbulencedriven shocks is
where is the volume of the simulation ( pc) and the summation is over the detected shocked cells.
For the shocks shown in Fig. 8, fast shocks dissipate times more energy than slow shocks due to the high velocity tail seen in Fig. 9. Note that not all of the kinetic energy from fast shocks dissipates by cooling, however, as some fraction goes into strengthening the magnetic field (as discussed in Section 2). Applying the MHD jump conditions for fast shocks with velocities ranging from 1–2 km/s, propagating into gas with the average fast preshock density g/cm and average preshock magnetic field strength G, we find that 35–70% of the kinetic flux is lost to the magnetic field, depending on the orientation of the shock direction with respect to the magnetic field. As energy stored in the magnetic field is free to further dynamically impact the turbulence, we reduce the energy dissipated by fast shocks by a factor of 2 in order to capture the energy being dissipated as heat.
Considering all shocks together, the shock dissipation rate is . Compared to the turbulent kinetic energy, this gives a dissipation timescale of Myrs. This is times the eddy turnover time where is the size of the diagonal of the simulation and is the onedimensional velocity dispersion. From observations of turbulent dissipation regions, Pon et al. (2014) estimated this ratio to be in the Perseus molecular cloud and Larson et al. (2015) made estimates of 0.94 and 0.65 for two shock models in the Taurus molecular cloud. Previous simulations have suggested that shocks dissipate around 50% of the turbulent kinetic energy (Stone et al., 1998). If we take this into account, our predicted ratio would be reduced by a factor of 2, placing it in the middle of these obversational results.
5 Discussion
We have presented an algorithm to detect and characterise fast and slow MHD shock waves in simulations of turbulent molecular clouds. While there is some observational evidence for the presence of fast (Lesaffre et al., 2013; Pon et al., 2014; Larson et al., 2015) and slow (Lehmann & Wardle, 2016) MHD shocks in molecular clouds, we present in this work the first prediction of the relative fraction of fast and slow shocks in molecular clouds. We characterised the shocks and provide the typical preshock conditions that should be used in future shock models that wish to model turbulencedriven MHD shocks. In the following we compare our work with other shockfinding algorithms, and then discuss how the results of the algorithm can be used to obtain an estimate of the volume of shock heated gas.
5.1 Comparison to previous work
Smith et al. (2000) and Smith et al. (2000) developed a method for counting shocks in MHD simulations of decaying and driven turbulence, respectively. Their method computes the velocity jump across converging regions. They found that the number distribution of these jumps did not substantially differ with the addition of a magnetic field. Our results are qualitatively similar, with weaker shocks dominating the number distribution, though they find much lower Mach numbers in general. This could be because they do not consider the shock reference frame, which introduces a correction to the velocity jump (see step (v) of section 3.1). Their method does not explicitly disentangle the MHD shock types, and so cannot quantify the relative importance of fast or slow shocks.
The importance of shock heating on the chemical evolution of turbulent molecular clouds is highlighted by Kumar & Fisher (2013). Like our work they capture the effects of shock heating on subgrid scales. They do this by postprocessing Lagrangian tracer particles in a simulation of hydrodynamic turbulence. Their subgrid model is a onedimensional integration of the fluid equations including a vast chemical network and molecular cooling. They were able to distinguish between chemicals that trace the mean physical state of cloud, and those that trace the nonequilibrium shockheated gas. Their method also accounts for solenoidal heating and so they can measure the relative importance of these two heating mechanisms. However, as they only consider hydrodynamic turbulence, their results are not sensitive to the distinct effects of MHD shock types. Extending their work to the MHD case is in principle simple. The subgrid model would need to include MHD effects like the shock models of Flower & Pineau Des ForÃªts (2010), Pon et al. (2012) or Lehmann & Wardle (2016). Some difficulty lies in obtaining the preshock state, which requires knowledge of the magnetic field direction with respect to the shock propagation direction. In addition, the preshock state does not uniquely determine the MHD shock type (cf Kennel et al., 1989). We have addressed this problem by using both pre and postshock information in order to ascertain two of the three possible shock types.
The shock finding algorithm most similar to shockfind is that outlined in Schaal & Springel (2015). They look for shocks in cosmological hydrodynamic simulations. Their method flags cells of converging flow, and defines the shock direction using the gradient of the temperature. They then use a series of criteria to filter spurious shock detections. While their work does not include magnetic fields, and thus does not consider different shock families, it would be simple to extend their algorithm to do so. It would only take the addition of further filtering criteria such as we presented in step (vi) of section 3.1. This extension would allow for a comparison of our work to MHD simulations using movingmesh codes such as arepo (Springel, 2010; Pakmor et al., 2011).
5.2 Shock heating
Lehmann & Wardle (2016) showed that Ctype fast MHD shocks and Jtype slow MHD shocks distinctly heat the gas they propagate through. The fast shocks modeled in Lehmann & Wardle reach peak temperatures of K, whereas slow shocks could reach temperatures of K. This is because in the weakly ionized gas that makes up molecular clouds, the heating timescale in fast shocks is determined by the ionneutral collision timescale, which is slower than the cooling timescale. In contrast, in slow shocks, the ionneutral collision timescale is shorter than the cooling timescale, such that heating in slow shocks is more significant. Even though these high temperatures only occupy a thin shock layer, the heating is important because of the rich chemistry it activates. The chemical signatures of shocks may persist in regions of unshocked gas. In this section we use the preshock conditions (Fig. 11) and shock front area (Fig. 9) obtained by shockfind, combined with representative subgrid twofluid shock models from Lehmann & Wardle (2016) to estimate the volume of shocked gas in a turbulent cloud.
5.2.1 Subgrid twofluid shock models
Ionneutral collisions determine shock thickness in twofluid shocks, and so the ionization fraction is a key variable in determining the volume filling fraction of shocks. Ionization sources, such as cosmicrays and ultraviolet photons, are density dependent and so the ionization fraction spatially varies in a turbulent cloud. For simplicity we adopt the ionization fraction of Bergin & Tafalla (2007): where we use the density in the preshock gas. This leads to preshock ionization fractions ranging between and .
For slow MHD shocks, the shock thickness is independent of the preshock magnetic field strength. So, for a given preshock density, we choose the magnetic field such that the Alfvén velocity km/s. This allows us to compute slow shocks with speeds up to , where is the angle between the direction of propagation and preshock magnetic field. The dashed line in Fig. 11 traces the preshock magnetic field and mass density at this Alfvén velocity, which roughly follows the distribution of slow shocks in the simulation. We model slow MHD shocks for allowing shock velocities up to 2.5 km/s. The number of slow shocks above this velocity in the simulation is only % of all slow shocks (see Fig. 9), so our results are only negligibly affected by this limit. Fig. 12 shows the thicknesses of slow MHD shocks heated above 50 K and 100 K, for preshock total hydrogen densities ranging between 10 and cm. The hot shock front is largest for models with the lowest preshock density and lowest velocity, peaking at cm which is of the order of the size of a cell . This is important because it means that the substructure within the shock front would not dynamically affect the scales that the simulation captures.
For fast MHD shocks with a given preshock density, we choose the magnetic field such that the Alfvén velocity km/s. The solid line in Fig. 11 traces this Alfvén velocity, which roughly follows the distribution of fast shocks in the simulation. The shock thickness, , of twofluid fast Ctype shocks is estimated as
where is the number density of ions and cm/s is the rate coefficient for ionneutral scattering. This shock thickness estimate can exceed the simulation box size at low densities and large shock velocities, which means, of course, that some fast shock models are innapropriate as subgrid models in this ideal MHD simulation. Thus we consider models with thickness as small enough to not significantly affect the simulation results. Note that this estimate gives the thickness of steadystate fast shocks. Factoring in the structure of nonsteady shocks is beyond the scope of this work. Finally, we find that models of fast shocks for these Mach numbers and densities in the range – cm are all Ctype shocks.
We bin the detected shocks into preshock total hydrogen densities centred on , , , and cm. These values of density almost cover the entire range of detected shocks, with of slow shocks and of fast shocks falling outside. With this binning, we can use the shock thicknesses from Fig. 12 to estimate the volume of warm shocked gas. Using the area computed for Fig. 9 multiplied by the shock thicknesses we find that of the volume is filled with shocked gas greater than 50 K. This is of the order of the shocked volume filling factor measurement of Pon et al. (2014) from turbulent dissipation regions in the Perseus molecular cloud. In addition, and of the volume is filled with shocked gas greater than 100 K and 150 K, respectively. This warm gas occurs entirely in slow shocks, because the fast shocks at these conditions do not reach peak temperatures above 50 K. Hence, if no distinction of MHD shock families is made and all shocks in an MHD simulation were assumed to be fast Ctype shocks, there would be no warm component of gas with temperature K at all.
5.2.2 Rotational line emission
While this predicted filling factor of gas hotter than 50 K, , is very small, only of the volume is filled with gas at densities higher than the average postshock density g/cm. As an example observational impact of this warm gas, we estimate the intensity of a CO rotational line using the non local thermodynamic equilibrium (LTE) radiative transfer code radex (van der Tak et al., 2007).
For a given radiating molecule, radex requires as input the density of H as the collisional partner, the column density of the radiating molecule, the temperature and the linewidth. In order to avoid geometrical effects, we consider only optically thin lines. By doing this, we can use radex in slab mode and treat each column through the simulation as consisting of a simple addition of slabs. We use the shock thicknesses computed above to define at each cell (on one face) the column density excited by gas at 75, 125 and 175 K. We use a CO abundance of to derive the CO column density. We then take the density of H to be equal to the average postshock density. We assume that the rest of the gas makes up a column of CO excited by 10 K gas at an H density equal to the average density in the simulation. Finally, we use the velocity dispersion of the cloud as the input linewidth.
Fig. 13 shows the synthetic map of radex estimated CO =98 intensities with contours of total hydrogen column density overlaid. We chose the =98 line because it was the lowest CO line that was optically thin at the column densities reached here and because another source of high CO lines, photodissociation regions, may have difficulties producing significant emission at this line and above (Pon et al., 2012). The emission in Fig. 13 is entirely due to shocks—as the background 10 K gas negligibly emits in this line—and is strongest in filamentary structures. This is because an edgeon shock, with respect to the line of sight, presents a larger column of heated gas than a faceon shock. These filamentary emission regions also tend to occur in regions with large hydrogen column densities, suggesting that the preshocked gas is already at a high density. The correlation between large hydrogen column density and emission is not perfect, however, as there is some significant CO emission at regions of below average hydrogen column density. If we add up the emission from over the whole face of the cloud, then the total cloud luminosity at this line is . Notably, if we ignored the distinction of MHD shock families and assumed that all shocks were the wellstudied Ctype fast shocks, we would predict that the CO =98 emission would be negligible.
Estimates of high CO lines like this could provide distinct observational predictions between different simulations of turbulent clouds. The accuracy of this estimate depends on the accuracy of the estimate of the volume of warm gas. This was estimated using the shock thicknesses derived from the twofluid shock models of Lehmann & Wardle (2016). It also only included the gas heated by slow shocks, because some twofluid Ctype fast shocks have thicknesses too large to be applicable to this simulation and the remaining fast shocks do not reach peak temperature of 50 K. This implies that a large proportion of the fast shocks detected here would not have the steadystate structure of twofluid fast shocks. We have also ignored the possibility of intermediate MHD shocks, because they do not have the predictable impact on magnetic fields that this algorithm exploits. In addition, the shock models of Lehmann & Wardle (2016) are highly simplified in order to highlight the differences between fast and slow shocks. Improvements in models of shocks, such as using an expanded chemical network and including the effects of dust grains, would improve the accuracy of the shockheated volume estimate.
6 Conclusion
The publicly available algorithm shockfind^{1}^{1}1Found on BitBucket (https://bitbucket.org/shockfind/shockfind) and the python Package Index (https://pypi.python.org/pypi/shockfind) was developed, which extracts and characterises the shock waves in MHD simulations. This algorithm was applied to a highresolution simulation of a magnetised, turbulent molecular cloud. We presented the first prediction of the relative fraction of fast and slow MHD shocks in this turbulent molecular cloud. The sonic and Alfvénic Mach number distributions for these two families of shocks are distinct and confirm that lowvelocities, below km/s, dominate the population of shocks. By considering the energetics of the detected shocks, we found that the ratio of the shock dissipation timescale to cloud crossing time is comparable to observed values from turbulent dissipation regions in molecular clouds. We have also used simple subgrid models of twofluid MHD shocks from Lehmann & Wardle (2016) to estimate the heating that would occur within the thin shock front of these shocks. Slow MHD shocks were found to produce a low volume filling factor, , component of the cloud heated above 50 K with a small portion of this component reaching temperatures above 150 K.
We used the nonLTE radiative transfer code radex to estimate the intensity of a high CO rotational transition and found that the shockheated gas radiates far above the background cloud intensity. High CO line emission may therefore be an important observational diagnostic of shocks in molecular clouds.
Our shockdetection algorithm is general enough to be applied widely to MHD simulations of other astrophysical phenomena. It would be interesting to see the mixture of shock families that might be present in simulations of supernova shocks, protostellar jets interacting with the interstellar medium, colliding flows, cloudcloud collisions, etc. In a future paper we plan to extract and characterise the MHD shocks in a variety of simulations of turbulent molecular clouds in order to search for correlations between the parameters of turbulence and possible observational effects of MHD shock waves.
Acknowledgements
The authors gratefully acknowledge discussions with James Tocknell and Birendra Pandey. This research made use of Astropy, a communitydeveloped core Python package for Astronomy (Astropy Collaboration et al., 2013). A.L. was supported by an Australian Postgraduate Award. C.F. acknowledges funding provided by the Australian Research Council’s Discovery Projects funding scheme (grants DP130102078 and DP150104329) and M.W. acknowledges funding provided by the Australian Research Council’s Discovery Projects funding scheme (grant DP120101792). We gratefully acknowledge the Jülich Supercomputing Centre (grant hhd20), the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Largescale 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. The software used in this work was in part developed by the DOEsupported Flash Center for Computational Science at the University of Chicago.
References
 André et al. (2014) André P., Di Francesco J., WardThompson D., Inutsuka S.I., Pudritz R. E., Pineda J. E., 2014, Protostars and Planets VI, pp 27–51
 Arzoumanian et al. (2011) Arzoumanian D. et al., 2011, A&A, 529, L6
 Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
 Bergin & Tafalla (2007) Bergin E. A., Tafalla M., 2007, ARA&A, 45, 339
 Cunningham et al. (2011) Cunningham A. J., Klein R. I., Krumholz M. R., McKee C. F., 2011, ApJ, 740, 107
 De Hoffmann & Teller (1950) De Hoffmann F., Teller E., 1950, Phys. Rev., 80, 692
 Draine (1980) Draine B. T., 1980, ApJ, 241, 1021
 Dubey et al. (2008) Dubey A. et al., 2008, in Numerical Modeling of Space Plasma Flows. p. 145
 Elmegreen & Burkert (2010) Elmegreen B. G., Burkert A., 2010, ApJ, 712, 294
 Elmegreen & Scalo (2004) Elmegreen B. G., Scalo J., 2004, ARA&A, 42, 211
 Falle & Komissarov (2001) Falle S. A. E. G., Komissarov S. S., 2001, J. Plas. Phys., 65, 29
 Federrath (2015) Federrath C., 2015, MNRAS, 450, 4035
 Federrath (2016) Federrath C., 2016, MNRAS, 457, 375
 Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761, 156
 Federrath & Klessen (2013) Federrath C., Klessen R. S., 2013, ApJ, 763, 51
 Federrath et al. (2014) Federrath C., Schrön M., Banerjee R., Klessen R. S., 2014, ApJ, 790, 128
 Federrath et al. (2011) Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011, ApJ, 731, 62
 Flower & Pineau Des ForÃªts (2010) Flower D. R., Pineau Des ForÃªts G., 2010, MNRAS, 406, 1745
 Fryxell et al. (2000) Fryxell B. et al., 2000, ApJS, 131, 273
 Glover et al. (2010) Glover S. C. O., Federrath C., Mac Low M.M., Klessen R. S., 2010, MNRAS, 404, 2
 Hacar et al. (2016) Hacar A., Kainulainen J., Tafalla M., Beuther H., Alves J., 2016, A&A, 587, A97
 Hennebelle & Chabrier (2008) Hennebelle P., Chabrier G., 2008, ApJ, 684, 395
 Hennebelle & Chabrier (2009) Hennebelle P., Chabrier G., 2009, ApJ, 702, 1428
 Hennebelle & Chabrier (2013) Hennebelle P., Chabrier G., 2013, ApJ, 770, 150
 Hennebelle & Falgarone (2012) Hennebelle P., Falgarone E., 2012, A&ARv, 20, 55
 Henshaw et al. (2016) Henshaw J. D. et al., 2016, MRNAS, 457, 2675
 Heyer & Brunt (2004) Heyer M. H., Brunt C. M., 2004, ApJ, 615, L45
 Hollenbach et al. (1971) Hollenbach D. J., Werner M. W., Salpeter E. E., 1971, ApJ, 163, 165
 Hopkins (2013) Hopkins P. F., 2013, MRNAS, 430, 1653
 Kainulainen et al. (2016) Kainulainen J., Hacar A., Alves J., Beuther H., Bouy H., Tafalla M., 2016, A&A, 586, A27
 Kennel et al. (1989) Kennel C. F., Blandford R. D., Coppi P., 1989, J. Plas. Phys., 42, 299
 Klessen & Hennebelle (2010) Klessen R. S., Hennebelle P., 2010, A&A, 520, A17
 Krumholz et al. (2007) Krumholz M. R., Klein R. I., McKee C. F., 2007, ApJ, 656, 959
 Krumholz & McKee (2005) Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250
 Kumar & Fisher (2013) Kumar A., Fisher R. T., 2013, MNRAS, 431, 455
 Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
 Larson et al. (2015) Larson R. L., Evans II N. J., Green J. D., Yang Y.L., 2015, ApJ, 806, 70
 Lehmann & Wardle (2016) Lehmann A., Wardle M., 2016, MNRAS, 455, 2066
 Lesaffre et al. (2013) Lesaffre P., Pineau des ForÃªts G., Godard B., Guillard P., Boulanger F., Falgarone E., 2013, A&A, 550, 106
 Mac Low & Klessen (2004) Mac Low M.M., Klessen R. S., 2004, Rev. Mod. Phys., 76, 125
 McKee et al. (2010) McKee C. F., Li P. S., Klein R. I., 2010, ApJ, 720, 1612
 McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565
 Mullan (1971) Mullan D. J., 1971, MNRAS, 153, 145
 Myers et al. (2014) Myers A. T., Klein R. I., Krumholz M. R., McKee C. F., 2014, MNRAS, 439, 3420
 Nakamura & Li (2007) Nakamura F., Li Z.Y., 2007, ApJ, 662, 395
 Offner & Arce (2014) Offner S. S. R., Arce H. G., 2014, ApJ, 784, 61
 Ossenkopf & Mac Low (2002) Ossenkopf V., Mac Low M.M., 2002, Astronomy and Astrophysics, 390, 307
 Padoan et al. (2014) Padoan P., HaugbÃ¸lle T., Nordlund Å., 2014, ApJ, 797, 32
 Padoan & Nordlund (2011) Padoan P., Nordlund Å., 2011, ApJ, 730, 40
 Padoan et al. (2015) Padoan P., Pan L., Haugboelle T., Nordlund A., 2015, ArXiv eprints, 1509, arXiv:1509.04663
 Pakmor et al. (2011) Pakmor R., Bauer A., Springel V., 2011, MNRAS, 418, 1392
 Pon et al. (2014) Pon A., Johnstone D., J. Kaufman M., Caselli P., Plume R., 2014, MNRAS, 445, 1508
 Pon et al. (2012) Pon A., Johnstone D., Kaufman M. J., 2012, ApJ, 748, 25
 Pon et al. (2016) Pon A. et al., 2016, ArXiv eprints, 1606, arXiv:1606.03089
 Price & Bate (2008) Price D. J., Bate M. R., 2008, MNRAS, 385, 1820
 Price & Bate (2009) Price D. J., Bate M. R., 2009, MRNAS, 398, 33
 Robertson & Goldreich (2012) Robertson B., Goldreich P., 2012, ApJ, 750, L31
 RomanDuval et al. (2011) RomanDuval J., Federrath C., Brunt C., Heyer M., Jackson J., Klessen R. S., 2011, ApJ, 740, 120
 Schaal & Springel (2015) Schaal K., Springel V., 2015, MNRAS, 446, 3992
 Schneider et al. (2013) Schneider N. et al., 2013, ApJ, 766, L17
 Smith et al. (2000) Smith M. D., Mac Low M.M., Heitsch F., 2000, A&A, 362, 333
 Smith et al. (2000) Smith M. D., Mac Low M.M., Zuev J. M., 2000, A&A, 356, 287
 Smith et al. (2014) Smith R. J., Glover S. C. O., Klessen R. S., 2014, MRNAS, 445, 2900
 Smith et al. (2016) Smith R. J., Glover S. C. O., Klessen R. S., Fuller G. A., 2016, MRNAS, 455, 3640
 Solomon et al. (1987) Solomon P. M., Rivolo A. R., Barrett J., Yahil A., 1987, ApJ, 319, 730
 Springel (2010) Springel V., 2010, MNRAS, 401, 791
 Stone et al. (1998) Stone J., Ostriker E., Gammie C., 1998, ApJ, 508, L99
 van der Tak et al. (2007) van der Tak F. F. S., Black J. H., Schöier F. L., Jansen D. J., van Dishoeck E. F., 2007, A&A, 468, 627
 VÃ¡zquezSemadeni et al. (2010) VÃ¡zquezSemadeni E., ColÃn P., GÃ³mez G. C., BallesterosParedes J., Watson A. W., 2010, ApJ, 715, 1302
 Waagan et al. (2011) Waagan K., Federrath C., Klingenberg C., 2011, J. Comput. Phys., 230, 3331
 Wang et al. (2010) Wang P., Li Z.Y., Abel T., Nakamura F., 2010, ApJ, 709, 27
 Wolfire et al. (1995) Wolfire M. G., Hollenbach D., McKee C. F., Tielens A. G. G. M., Bakes E. L. O., 1995, ApJ, 443, 152
 Wu (1987) Wu C. C., 1987, GeoRL, 14, 668
Appendix A Cylinder Radius
Here we check the effect of varying the radius of the cylinder used to define average pre and postshock variables (step (iii) in Section 3.1). Fig. A shows the Mach number distributions for fast and slow shocks for runs of shockfind on the same locations, but with cylinder radii of , 3 and 5 in units of cell size. A radius of defines the minimum cylinder size to represent a line that does not suffer aliasing effects. This size, however, is still affected by small scale numerical noise. Averaging over larger radii, and , removes the effect of the numerical noise. There is little difference between these two larger radii, and so we choose as the compromise between noise and lowering computational time.
Appendix B Overcounting the shocked cells
In Fig. 15 we show a close up of the detected shocked cells of a curved shock front. In this figure, it can be seen that the thickness of the shock front is 3 or 4 rows of detected cells. As the shock front area is required to link our work to the results of onedimensional shock models, we are overcounting if we consider every detected cell as a unique unit of shocked area. In order to avoid this over counting we only use (as step (viii) of Section 3.1) cells that occur at the local convergence maxima (filled squares in Fig. 15) along the extracted line (step (iii) of Section 3.1).
Appendix C Convergence Test
Here we check that the search thresholds capture converged distributions of shocks. The convergence threshold defined in Sec. 4.1 searches down to of cells with negative divergence (), i.e., of all possible shocked cells. Fig. 16 shows the effect of doubling the number of cells searched on the Mach number distributions of fast and slow shocks. The distributions are converged for the most observationally important shocks ().