Shock identification and classification in 2D magnetohydrodynamiccompressible turbulence—Orszag–Tang vortex

Abstract

Compressible magnetohydrodynamic (MHD) turbulence is a common feature of astrophysical systems such as the solar atmosphere and interstellar medium. Such systems are rife with shock waves that can redistribute and dissipate energy. For an MHD system, three broad categories of shocks exist (slow, fast, and intermediate); however, the occurrence rates of each shock type are not known for turbulent systems. Here, we present a method for detecting and classifying the full range of MHD shocks applied to the Orszag–Tang vortex. Our results show that the system is dominated by fast and slow shocks, with far less-frequent intermediate shocks appearing most readily near magnetic reconnection sites. We present a potential mechanism that could lead to the formation of intermediate shocks in MHD systems, and study the coherency and abundances of shocks in compressible MHD turbulence.

1. Introduction

Many astrophysical systems are both highly compressible and turbulent, for example, the solar atmosphere (Carlsson & Stein, 1992; Houston et al., 2020; Reardon et al., 2008; Ulmschneider, 1970) and interstellar medium (Draine et al., 1983; Elmegreen & Scalo, 2004). 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, 1971). 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, 2012; Orta et al., 2003; Park & Ryu, 2019); however, intermediate shocks are far less studied, and their existence has been somewhat controversial in the past (Karimabadi, 1995; Wu, 1988). These shock jumps are fully permitted by MHD equations (Hau & Sonnerup, 1989) and have recently been observed in the solar chromosphere (Houston et al., 2020). 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, 1989). In two-fluid simulations, it was shown that a slight upstream acceleration of the plasma can result in an intermediate shock forming (Snow & Hillier, 2019). Previous studies have suggested that intermediate shocks may be related to magnetic reconnection (La Belle-Hamer et al., 1994; Ugai & Shimizu, 1994).

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, 1950), 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, 1989).

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. 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. 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. 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. 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, 1979). 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, 1989; Jiang & Wu, 1999; Parashar et al., 2009; Picone & Dahlburg, 1991; Uritsky et al., 2010). 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., 2016), which has been used previously to study shocks (e.g., Snow & Hillier, 2021). 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, 1989; Picone & Dahlburg, 1991).

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, 1995) may result in far fewer slow-mode shocks being detectable in turbulent simulations (Park & Ryu, 2019). 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, 2021).

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, 1964), 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, 2011). 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, 1989). 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 (2011); 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,384²) 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,024² simulation; see Figure 1). Intermediate shocks are proportionally less-frequent; intermediate 2–4 shocks are $ \approx $3% as frequent as slow shocks at 1,024², 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.