1. Introduction
Many astrophysical systems are both highly compressible and turbulent, for example, the solar atmosphere (Carlsson & Stein, Reference Carlsson and Stein1992; Houston et al., Reference Houston, Jess, Keppens, Stangalini, Keys, Grant, Jafarzadeh, McFetridge, Murabito, Ermolli and Giorgi2020; Reardon et al., Reference Reardon, Lepreti, Carbone and Vecchio2008; Ulmschneider, Reference Ulmschneider1970) and interstellar medium (Draine et al., Reference Draine, Roberge and Dalgarno1983; Elmegreen & Scalo, Reference Elmegreen and Scalo2004). A common feature of compressible turbulence is shocks, where the fluid properties change drastically over a small area. Such sharp structures are important for dissipation. For systems that are governed by magnetohydrodynamic (MHD) equations, there are three broad categories of stable shocks: slow, fast, and intermediate (e.g., Tidman & Krall, Reference Tidman and Krall1971). Understanding the types of shocks that form in a system and where they are likely to form is critical to understand the energy transfer in a turbulent system.
Significant work has been performed to analyze the fast and slow shocks in turbulent systems (e.g., Komissarov, Reference Komissarov2012; Orta et al., Reference Orta, Huerta and Boynton2003; Park & Ryu, Reference Park and Ryu2019); however, intermediate shocks are far less studied, and their existence has been somewhat controversial in the past (Karimabadi, Reference Karimabadi1995; Wu, Reference Wu1988). These shock jumps are fully permitted by MHD equations (Hau & Sonnerup, Reference Hau and Sonnerup1989) and have recently been observed in the solar chromosphere (Houston et al., Reference Houston, Jess, Keppens, Stangalini, Keys, Grant, Jafarzadeh, McFetridge, Murabito, Ermolli and Giorgi2020). However, when and how intermediate shocks can form, and their relative stability, is far from understood. The “strong” intermediate shock (transition from super-Alfvénic to subslow) is smoothly connected to the slow-mode shocks (Hau & Sonnerup, Reference Hau and Sonnerup1989). In two-fluid simulations, it was shown that a slight upstream acceleration of the plasma can result in an intermediate shock forming (Snow & Hillier, Reference Snow and Hillier2019). Previous studies have suggested that intermediate shocks may be related to magnetic reconnection (La Belle-Hamer et al., Reference La Belle-Hamer, Otto and Lee1994; Ugai & Shimizu, Reference Ugai and Shimizu1994).
Here, we present a method to detect and classify the full range of MHD shocks from a 2D simulation and apply this to the Orszag–Tang (OT) vortex. The results indicate that the vast majority of the shocks are either fast or slow, and hence these shocks are responsible for the majority of the dynamic consequences. Less-frequent intermediate shocks are detected and appear to form in high-current regions, that may arise due to turbulent inflow to a reconnecting region.
2. Methods
2.1. Classification of shocks
Shocks can be classified based on their upstream and downstream velocity relative to the characteristic speeds of the system (see Table 1). For correct classification, the velocity must be in the shock frame, that is, a frame of reference where the shock is stationary, for example, the de Hoffmann–Teller frame (de Hoffmann & Teller, Reference de Hoffmann and Teller1950), where the electric field is zero either side of the shock. In the de Hoffmann–Teller frame, it is possible to derive an equation that gives all possible steady-state transitions as a function of the upstream (u) and downstream (d) Alfvén Mach numbers, and the upstream plasma-β and angle of magnetic field (Hau & Sonnerup, Reference Hau and Sonnerup1989).
Table 1. Possible stable magnetohydrodynamic shock transitions based on the upstream and downstream states
Note: Transitions of the form u > d are entropy-forbidden, so not listed here.
2.2. Identification technique
To automatically identify shocks in our simulations, we use the following procedure:
1. Identify shock candidates based on the divergence of the velocity field (
$ \nabla \cdot \mathbf{v}<0 $ is a necessary condition for a shock).2. Calculate the parallel and perpendicular unit vectors based on the maximum density gradient. The list of possible shocks is then filtered by checking the density gradient along a line perpendicular to the shock front. If the shock candidate is not a local maxima or minima of the perpendicular density gradient, then it is not the center of the feature and is discarded. This prevents the numerical (and physical) finite width of the shock resulting in multiple detections.
3. Estimate the shock frame from the steady-state conservation of mass equation. In the de Hoffman–Teller shock frame, the conservation of mass equation becomes
$ {\left[\rho {v}_{\perp}\right]}_d^u=0 $. To transfer the simulation laboratory frame to the shock frame, we need the shock velocity $ {v}_s $. By rewriting the mass conservation including this constant shock velocity as $ {\left[\rho \left({v}_{\perp }+{v}_s\right)\right]}_d^u=0 $, where $ {v}_{\perp } $ is directly from the simulation, we can estimate the shock velocity by rearranging as $ {v}_s=\frac{\rho^d{v}_{\perp}^d-{\rho}^u{v}_{\perp}^u}{\rho^u-{\rho}^d} $.4. Compare the upstream and downstream Mach numbers to calculate the shock transition.
There are a number of assumptions that go into this identification technique, with the strongest being that the shocks are steady-state. In a highly dynamic simulation, this may prevent transient shocks being detected as readily and could lead to erroneous shock detection. As such, our identification procedure is likely an underestimate of the number of shocks in the system. Furthermore, we include additional checks to confirm the compressibility is greater than unity for all shocks, and that there is a reversal in the magnetic field across the interface for the intermediate shocks.
2.3. MHD simulation
To test this shock detection code, we use the OT shock vortex (Orszag & Tang, Reference Orszag and Tang1979). The OT initial conditions produce a decaying compressible turbulent system starting from the following initial conditions:
$$ {\displaystyle \begin{array}{l}{v}_x=-{v}_0\sin \left(2\pi y\right),\quad (1)\qquad \rho =25/\left(36\pi \right),\quad (4)\quad \\ {v}_y={v}_0\sin \left(2\pi x\right),\qquad (2)\qquad P=5/\left(12\pi \right),\quad (5) \\ {}{B}_x=-{B}_0\sin \left(2\pi y\right),\quad (3)\qquad {v}_0=1, \qquad\qquad(6)\\ {} {B}_y={B}_0\sin \left(4\pi y\right),\qquad(4)\quad {B}_0=1/\sqrt{4\pi },\quad\quad (7)\end{array}} $$for density 
$ \rho $, pressure P, velocity 
$ \mathbf{v}=\left[{v}_x,{v}_y,0\right] $, and magnetic field 
$ \mathbf{B}=\left[{B}_x,{B}_y,0\right] $. The OT vortex has been well studied in the literature (Dahlburg & Picone, Reference Dahlburg and Picone1989; Jiang & Wu, Reference Jiang and Wu1999; Parashar et al., Reference Parashar, Shay, Cassak and Matthaeus2009; Picone & Dahlburg, Reference Picone and Dahlburg1991; Uritsky et al., Reference Uritsky, Pouquet, Rosenberg, Mininni and Donovan2010). The initial conditions are evolved in 2D for ideal MHD equations using the fourth-order central-difference solver in the (PIP) code (Hillier et al., Reference Hillier, Takasao and Nakamura2016), which has been used previously to study shocks (e.g., Snow & Hillier, Reference Snow and Hillier2021). The simulations are performed using 1,024 × 1,024 cells with periodic boundary conditions.
3. Results
Figure 1 shows the evolution of the OT vortex through time, with the identified shock pixels colored according to their classification. The simulation shows the full range of shock transitions exists in our simulation, with the most prevalent being the fast- and slow-mode shocks. Intermediate shocks are detected in the simulation; however, these are far less-frequent and appear later than the fast or slow shocks.
Figure 1. Shocks detected in the Orszag–Tang vortex in magnetohydrodynamic at t = 0.3 (a), 0.6 (b), 0.8 (c), and 1.0 (d). Background shows the out-of-plane current. Overplotted are the detected slow (red), fast (blue), and intermediate (yellow, magenta, cyan, and green) shocks. The occurrence counts of each type of shock through time are shown in the panel (e). The panel (f) shows the coherency of fast and slow shocks through time.
Figure 1e shows the number and type of detected shocks in the simulation through time. At very early times in the simulation (t < 0.15), the detections of shocks are sporadic, because the system is evolving rapidly and the assumption of a steady-state shock is not valid. After t = 0.15, large-scale fast-mode shocks are generated due to the initial conditions of the OT vortex. As the system develops and the shocks interact, slow-mode shocks become increasingly detected in the system, and fast and slow shocks have comparable detection counts near the end of the simulation (t = 1). The relative counts of fast- and slow-mode shocks are likely dependent on the initial conditions of the system (Dahlburg & Picone, Reference Dahlburg and Picone1989; Picone & Dahlburg, Reference Picone and Dahlburg1991).
A metric to determine the coherency of shocks can be calculated by measuring the number of similar shocks within a given radius. Here, we use a radius of three grid cells and classify any shocks that have more than six nearby shocks of the same type as coherent. Figure 1f shows the proportion of coherent fast and slow shocks through time. At early times in the simulation, the system is dominated by long, continuous, fast-mode shocks, and hence the proportion of coherent shocks is very high. Because the simulation evolves and the fast-mode shocks interact and break up, the coherency decreases. For the slow-mode shocks, the coherency is fairly constant through time at around 40%. From Figure 1b–d, one can see that the large-scale slow-mode shocks may exist, but are not detected as continuous features, possibly due to either limitations of the detection method.
It has been suggested that the corrugation instability (e.g., Stone & Edelman, Reference Stone and Edelman1995) may result in far fewer slow-mode shocks being detectable in turbulent simulations (Park & Ryu, Reference Park and Ryu2019). However, here, there are comparable numbers of fast and slow shocks in the system after t = 0.8, and large un-corrugated coherent slow-mode shock structures exist in the simulation (see Figure 1b–d). A recent study has shown that the corrugation instability can actually increase the number of detected shock pixels due to an increased shock length (Snow & Hillier, Reference Snow and Hillier2021).
3.1. Potential mechanism for the formation of intermediate shocks
The intermediate shocks appear to form near and around high-current regions. As such, one may consider the link between magnetic reconnection and intermediate shocks in turbulent systems. Note that not all high current regions feature intermediate shocks, but intermediate shocks are usually located at a high-current region. For fast magnetic reconnection, one may expect switch-off slow-mode shocks (e.g., Petschek, Reference Petschek1964), where the inflow Alfvén Mach number is unity. High-resolution MHD simulations have shown that switch-off slow-mode shocks are integral to plasmoid formation (Zenitani & Miyoshi, Reference Zenitani and Miyoshi2011). Now, since a switch-off slow-mode shock has an inflow Alfvén Mach number of 1, it is on the cusp of being a 
$ 2\to 4 $ intermediate shock (where the inflow is super-Alfvàic and the downstream is subslow, e.g., Figure 1 in Hau & Sonnerup, Reference Hau and Sonnerup1989). As such, for a turbulent reconnection region, one would expect variations in the inflow quantities that could result in a super-Alfvénic inflow and an intermediate shock. This appears to be the formation mechanism occurring here. Specifically isolating a plasmoid in a high-resolution MHD simulation, one sees remarkable similarities to the plasmoid schematic in Zenitani and Miyoshi (Reference Zenitani and Miyoshi2011); however, here, intermediate shocks form on the inflow (see Figure 2). In this example figure, the inflow to the plasmoid is asymmetric, and slow-mode shocks form on the lower inflow, whereas intermediate shocks form on the upper inflow.
Figure 2. Plasmoid with intermediate shocks forming on inflow region. Fundamental structures of the plasmoid are annotated.
In terms of relative stability, in these simulations, the intermediate shocks can be regularly detected across a few time steps, and by eye, one can track certain shocks across outputs; however, others are less stable (see Figure 3). Through time, the intermediate shocks are routinely detected at the high-current regions. Further work is needed to analyze the stability and evolution of intermediate shocks.
Figure 3. Snapshots of a reconnecting region in the time frame t = 0.66–0.73 (from [a]–[h]) at intervals of 0.01. Background color is the out-of-plane current density.
3.2. Consistency of results at higher grid resolutions
A high-resolution (16,3842) simulation of the OT vortex was performed to analyze the consistency of the results. At time t = 1, the number of pixels that satisfy a shock jump condition in the high-resolution simulation are: 160,819 (slow), 316,029 (fast), 350 (intermediate 1–2), 23 (intermediate 1–4), 836 (intermediate 2–3), and 1,394 (intermediate 2–4). The increased grid resolution has led to many more shocks being detected; however, proportionally, the results are similar. There is a factor of 
$ \approx 2 $ faster than slow shocks (as in the 1,0242 simulation; see Figure 1). Intermediate shocks are proportionally less-frequent; intermediate 2–4 shocks are 
$ \approx $3% as frequent as slow shocks at 1,0242, and 
$ \approx $0.1% in the high-resolution simulation. This may be due to the larger scale shocks contributing more to the overall shock count and requires further study. However, fast and slow shocks are most abundant, and the intermediate shocks form at high-current regions in both simulations indicating that these are consistent results.
4. Conclusion
The developed method for identifying shocks in MHD systems detects the full range of MHD shocks in the OT vortex. The MHD simulation is dominated by fast and slow shocks implying that the vast majority of the shock-related dissipation occurs due to fast and slow shocks. Furthermore, the fast and slow shocks occur with roughly equal detected pixel counts at t = 1 implying that the corrugation instability is not preventing slow-mode shocks from being detected. Less-frequent intermediate transitions are also detected in the simulation that appear to be connected to high-current regions. In particular, the 2–4 intermediate shock may be formed due to fluctuations of the inflow velocity from a reconnection-generated switch-off shock. Further work is needed to analyze the evolution and stability of intermediate shocks occurring in turbulent MHD systems.
Acknowledgments
We would like to thank the stimulating discussions with Joanne Mason, Pierre Lesaffre, Andrew Lehmann, and Thibaud Richard.
Funding Statement
B.S. and A.H. are supported by STFC research grant ST/R000891/1 and ST/V000659/1.
Conflict of Interest
The authors have no conflicts of interest to declare.
Authorship Contributions
B.S. designed the study and wrote the manuscript. A.H. and G.J.J.B. provided insight and assisted in understanding the results. G.M. assisted in benchmarking the shock identification code.
Data availability Statement
The simulation data from this study are available from B.S. upon reasonable request. The (PIP) code is available at https://github.com/AstroSnow/PIP. The IDL scripts used for this manuscript are available at https://github.com/AstroSnow/Orszag-Tang-Shock-Detection.
Open access
Comments
Comments to the Author: This is an excellent paper that characterizes fast, slow, and intermediate shocks in the Orszag-Tang vortex problem -- a classic test problem for numerical simulations and an example of supersonic decaying turbulence. I recommend it for publication, following consideration of my comments below.
I found the presentation very clear and the analysis straightforward. The results about which types of shocks are responsible for dissipating energy in this test problem is interesting.
The one comment I have is that the authors’ simulate the ideal MHD equations at a single fixed resolution. Without physical resistivity, the solution is known not to converge, and dissipation and formation of plasmoids is resolution dependent (e.g. one would not even form plasmoids at low-resolution 64x64). So my question is, how robust are the findings of which shock types dominate to increased resolution, e.g. 2048^2, 4096^2, 8192^2? My experience simulating the OT problem is that plasmoid structure can change significantly with increased resolution. It is possible that results/dissipation are converged in 2D because 2D turbulence is an inverse-cascade. Numerical studies of dissipation in 3D problems though, where you have cascade to small-scales, have proved to be more difficult due to effects of numerical viscosity/resistivity.
Typos:
* page 4. "vortex.As" --> "vortex. As"