Plasmoid Instability in the Multiphase Interstellar Medium

In Figures 1 and 2, we show 3D volume renderings and 2D slices of our fiducial simulation after five eddy turnover times of the warm gas (5L/v_rms,warm). The top panels of both figures show the full volume, while the surrounding panels show zoom-ins on three regions. These regions (the same in both figures) exemplify the three distinct environments in which we find current sheets undergoing plasmoid-mediated magnetic reconnection: (i) within cold clumps, (ii) at the asymmetric interface of the cold and warm phases, and (iii) within the warm, volume-filling phase.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the whole, the volume is characterized by a bistable configuration with the vast majority of the gas residing near one of the two thermal equilibria. The cold gas morphology is clearly clumpy with structures on scales ranging from nearly the size of the whole box (∼tens of parsecs) down to the grid scale (∼10⁻² pc). Many of these cold gas clumps show an elongated, filamentary structure with large aspect ratios (discussed in Section 3.3).

The field lines shown in Figure 1 give a sense of the large- and small-scale magnetic structure. While some of the volume contains relatively ordered field lines, they also generically exhibit highly disturbed regions manifesting as small-scale helical knots or "flux tubes," which trace plasmoids and are hallmarks of tearing instability–induced fast magnetic reconnection. This highlights that plasmoids form throughout the volume and can be found in the volume-filling warm and clumpy cold phases and at the interface of the two. This is further supported by the 2D slice of the electric current, j ≡ ∇ × B , in the top left panel of Figure 2. It is clear that there are numerous large current sheets with L_CS ≳ L/10 throughout the volume (a well-known phenomenon in MHD turbulence; see, e.g., Zhdankin et al. 2013), many of which are actively undergoing plasmoid-mediated magnetic reconnection. In Appendix A, we show that similar large current sheets form in lower-resolution simulations, but they are much more laminar and do not show signs of plasmoid instability.

The top panels of Figures 1 and 2 further demonstrate that the thermodynamic and magnetic properties interact nontrivially. Regions with cold gas tend to be more magnetically dominated with β < 1 (blue-green colors). This is likely due to magnetic field enhancement during cooling-induced contraction. Moreover, the regions of highly disordered magnetic fields and large j tend to be associated with thermally unstable or cold gas.

Turning now to the zoom-in panels of Figures 1 and 2, we get a better sense of the nature of reconnection in these environments and the morphology and properties of the plasmoids. Each region has one or more plasmoids ¹¹ and several clear X-points where the reconnection is actively occurring. ¹² Each of these regions exemplify the three distinct environments in which we find evidence for the plasmoid instability. These plasmoids closely resemble the structures seen in local simulations of individual plasmoid unstable current sheets (e.g., Daughton et al. 2011; Huang & Bhattacharjee 2016; Yang et al. 2020), thereby bolstering the evidence that these structures arise from the plasmoid instability in our marginally resolved global simulations. Furthermore, the plasmoid unstable current sheets have a thickness of a few cells (consistent with numerical resistivity and the estimate of Grete et al. 2023) and a length of >200 cells. ¹³ This corresponds to an axis ratio of ≳100 and thus a Lundquist number S ≳ 10⁴, which is comparable to or greater than the threshold for plasmoid instability (Loureiro et al. 2007). We can, therefore, be confident that these structures arise from the plasmoid instability. We encourage readers to view the animated versions of our plots at https://dfielding14.github.io/movies/ to get a better sense of the morphology and different evolutionary stages of these structures. In general, the plasmoid instability begins in long, relatively undisturbed current sheets that progress through the linear evolution and develop a series of X-points and plasmoids. The plasmoids eventually evolve into a turbulent nonlinear state.

The zoom region shown on the left is an example of a current sheet that forms in the middle of a cold region. Current sheets like this arise in large part because of the compressive motions that result from the thermally unstable region rapidly contracting in on itself. The energy released during reconnection in these cold current sheets is quickly radiated away, and the gas inside the plasmoid remains cold. This is demonstrated by a temperature profile that remains cold along the ray passing through an X-point in the current sheet.

The zoom region shown in the middle is closely related to the first, with the important distinction that the current sheet forms at the interface of the warm and cold phases. The pressure and magnetic field strength on either side of the current sheet are, in this example, nearly constant; thus, the asymmetry is in both the temperature and density. ¹⁴ Perpendicular to the current sheet, the temperature changes by nearly a factor of 100 (lower middle panel of Figure 2). These asymmetric current sheets are formed as a result of the flow of warm gas as it cools onto a cold clump with an oppositely oriented magnetic field. In the ISM, these strong cooling flows, which occur throughout the volume, can lead to the development of plasmoid unstable current sheets. Here reconnection may modify the phase balance, as the heating is occurring directly in the thermally unstable phase, where it can be most impactful. The plasmoids that form in these asymmetric current sheets tend to be highly elongated, often have cold interiors, and follow the surface of the larger cold clumps.

These first two cases are unique to the thermally unstable ISM-like environments we are considering and represent a new channel for magnetic reconnection. The reconnection is brought on by thermal instability–induced tearing instabilities and leads to the formation of what we refer to as cold plasmoids (i.e., plasmoids that form as a result of radiative cooling and remain cold because of their short cooling times). As we discuss in Section 4, cold plasmoids may be readily observable as H i absorption and/or emission and as sources for strong radio scattering. Furthermore, because of their partial magnetic insulation/shielding, these cold plasmoids (as shown in Figure 1) can survive for an appreciable fraction of the global dynamical time. In particular, their magnetic structure can make their lifetimes significantly longer than cold clumps of a similar size, which tend be mixed away on an eddy turnover time (t_mix ∼ r_cl/v_turb(r_cl) ≲ r_cl/c_s,cold ∼ 0.1 Myr for r_cl = 0.1 pc).

The cold and asymmetric plasmoid unstable current sheets can be contrasted to the final zoom region shown in the right panels, in which a current sheet forms in the volume-filling warm phase as a result of the intermittency in the turbulent flow. The current sheet in this case has β ≳ 100, and the energy release from the reconnection has heated the immediate surrounding medium. This hot plasmoid formation, as we refer to it because it takes place in the hotter, volume-filling phase, has been seen in simulations of single-phase highly magnetized turbulent environments (e.g., Chernoglazov et al. 2021; Dong et al. 2022; Galishnikova et al. 2022) and can be expected to be common throughout the volume-filling region of the ISM.

Looking at the profiles along the rays that pass through the X-point of the most prominent current sheets in each of the zoom regions gives us a more detailed view into the properties of the current sheets. The temperature profiles clearly demonstrate our distinction between the cold, asymmetric, and warm current sheets. Each of the profiles exhibits a large reversal in the reconnecting component of the magnetic field (B_rec) over a width of ∼L/500. The behavior of the component of the magnetic field aligned with the guide field (B_guide; i.e., the component aligned with the field direction in the middle of the current sheet) is noticeably different between the three profiles. Near the cold current sheet (left), B_guide ∼ B_rec; in the asymmetric current sheet, B_guide ≲ B_rec/4 and remains coherent above and below the current sheet; and in the warm current sheet, B_guide ∼ 0. Therefore, in addition to representing three different thermal properties, these three current sheets also represent three guide field properties, namely, strong, weak, and zero guide field, respectively.

Also shown in the profiles are the inflow velocity v_rec (velocity toward the current sheet that is driving reconnection) and the shear velocity v_shear (velocity along the current sheet) normalized by one-tenth of the average local Alfvén speed 0.1v_A. In all cases, v_rec ∼ 0.01–0.03v_A. This reconnection rate is comparable to the fast reconnection expected for tearing-mediated reconnection (e.g., Loureiro et al. 2007; Bhattacharjee et al. 2009; Uzdensky et al. 2010). Furthermore, in all cases, v_shear < v_A, which implies that the current sheets are stable to the Kelvin–Helmholtz instability (Loureiro et al. 2013; Chernoglazov et al. 2021). That said, it is likely that some of the plasmoid/flux tube structures in our simulation arise as a result of instabilities other than the plasmoid instability. This only bolsters our result that these structures are likely present in the ISM.

We now shift our attention to the statistical properties of the simulation with an eye toward understanding the combined effects of turbulence, thermal instabilities, and magnetic reconnection.

3.2.1. Magnetothermal Phase Structure

Figure 3 provides an overview of the thermal and magnetic properties and their correlations. We unpack these results incrementally, starting with a discussion of the middle left panel, which shows the temperature distribution of the mass (solid lines) and volume (dashed lines). The total mass of the domain is dominated by a broad distribution of cold gas with T < T₀/10. There is a peak of mass at the warm-phase equilibrium near 10T₀. Connecting these two peaks and populating the thermally unstable regime is a nearly flat distribution that is continually depleted and repopulated by the combined effects of radiative cooling and heating, turbulent mixing, and reconnection heating. The volume, on the other hand, is dominated by the warm phase. This is a reflection of the multiphase nature shown in the previous figures.

Zoom In Zoom Out Reset image size

What determines the exact balance between the amount of mass in each of the two phases remains an open question. The answer depends on a nontrivial combination of the primary physical processes at play. The strength of conduction has been shown to regulate the balance in more simplified systems (e.g., Zel'Dovich & Pikel'Ner 1969; Choi & Stone 2012; Kim & Kim 2013; Jennings & Li 2021). Furthermore, we found that in a simulation with a longer cooling time relative to the eddy turnover time (e.g., a smaller box), the fraction in the cold phase decreased. Likewise, increasing the strength of turbulence can significantly modify the phase distribution and, in the extreme limits, erase the bimodal structure altogether (e.g., Vázquez-Semadeni et al. 2000; Gronke et al. 2022). How magnetic reconnection may impact this balance is an important future consideration given the difficulty of adequately resolving this process in ISM simulations.

To partially illuminate the role of both magnetic reconnection on the thermal structure and thermal instabilities on the magnetic structure, we show with purple lines in the middle left panel of Figure 3 the mass and volume temperature distributions of gas with large magnetic curvatures, ${K}_{B}\,\equiv \,(\hat{b}\,\cdot \,{\boldsymbol{\nabla }})\,\hat{b}$, where $\hat{b}={\boldsymbol{B}}/| {\boldsymbol{B}}|$. We adopt a threshold of K_B > 512/L, which corresponds to the field changing direction in four cells or less (the middle right panel shows the full K_B distribution, and the small dashed line indicates our choice of threshold). High K_B is a good indicator of magnetic reconnection because of the large curvature inherent to plasmoids and X-points (as can be seen in the bottom right panel, which shows the curvature distribution of the plasmoid subvolumes). The temperature distributions of the high-K_B gas are skewed toward more cold gas than the volume as a whole. This indicates that current sheet formation and reconnection are enhanced in the cold phase, likely as a result of the compression of the field due thermal instability–induced contraction.

The curvature, K_B , alone is an imperfect tracer of plasmoids, so in the bottom left panel of Figure 3, we show the temperature distribution of the high-K_B gas in the three example subvolumes that contain clear plasmoids shown in the previous figures. All three subvolumes exhibit a multiphase structure; however, the cold current sheet has a negligible amount of warm, highly magnetically curved material. This demonstrates a key point that the heating by magnetic reconnection can be rapidly lost due to radiative cooling and is thus insufficient to drive the material to the hotter phase. How this heating impacts the global phase structure is an interesting discussion for future work.

The volume distribution of K_B in the whole volume (black), cold gas (blue), warm gas (yellow), and example plasmoid subvolumes is shown in the middle and bottom right panels of Figure 3. As expected, the plasmoid subvolumes have a large high-curvature volume fraction owing to the fact that plasmoids are, by definition, comprised of highly curved field lines. At the high-curvature end, the total curvature distribution scales proportional to $\propto \,{K}_{B}^{-8/5}$. On average, magnetic fields in cold gas (T < T₀/10) tend to be more curved than in warm gas (T > 8T₀), but the highest K_B values are equally likely in either phase. The curvature distribution of all of the gas is remarkably similar to the isothermal driven turbulence simulation we ran for comparison.

To better understand the magnetothermal phase distributions, in the top left panel of Figure 3, we show the joint distribution of pressure and temperature (the contours trace the mass distribution, and the colors show the average β). As expected, the majority of the mass is found near the two stable thermal equilibria (i.e., along the negatively sloping portions of the equilibrium curve). The gas near the warm equilibrium is tightly clustered around P ≈ P₀ and T ≈ 10T₀ and has an average β ∼ 1. The mass near the cold equilibrium is spread out along the stable equilibrium curve with 10^−2.5 ≲ T/T₀ < 0.1 and 0.1 ≲ P/P₀ < 3. The β also varies in the cold phase, with the lower-pressure cold gas having β ≲ 1/10, indicating a strong compensation of the thermal pressure deficit by enhanced magnetic pressure (B²/2). The enhancement of magnetic pressure in regions with lower thermal pressure means that, although the gas near the cold and warm thermal equilibria are not in thermal pressure equilibrium, they are close to being in total pressure equilibrium with $P+{B}^{2}/2\approx \mathrm{constant}$. There is some mass with very low temperatures, T ≲ T₀/100, that has both high thermal pressure and low β, which is found in rapidly contracting thermally unstable regions that often host cold current sheets.

The top right panel of Figure 3 shows the joint mass distribution of thermal and magnetic pressure and is colored by the average temperature in each bin. Most of the mass is found at a pressure of P ≈ 0.3P₀ and a magnetic pressure of B²/2 ≈ 2P₀, corresponding to β ≈ 10^−0.75. The mass distribution has a distinct shape with a tail toward extremely low magnetic pressures. This tail follows a line of nearly constant total pressure. Gas moves both ways along lines of nearly constant total pressure. On the one hand, as warm gas cools and contracts, it loses thermal pressure and gains a nearly equal proportion of magnetic pressure. On the other hand, highly magnetized gas that undergoes reconnection moves in the other direction in this phase space; the conversion of magnetic energy into thermal energy raises the thermal pressure of the gas as it loses magnetic pressure. In some cases, however, the magnetic reconnection heating is not sufficient to overcome the strong radiative cooling, so the end result is cold gas with low thermal and magnetic pressure, as can be seen by the lobe in the mass distribution of cold gas with P ∼ 0.3P₀ and B²/2 ≲ 10⁻¹. This is the case in the cold plasmoid formation channel exemplified by the cold and asymmetric current sheets shown in the previous two figures and the bottom panels of Figure 3. Note that in addition to heating, reconnection also accelerates gas, so some appreciable fraction of the magnetic energy is also being converted to small-scale motions.

The detailed analyses presented in Figure 3 highlight the rich interplay between the thermodynamics, turbulence, and magnetic fields that combine to regulate the magnetothermal phase structure of the ISM. While there is clearly still much to be understood (particularly in the computationally daunting limit of resolved explicit resistivity, conduction, and viscosity), we have shown that there are important correlated magnetic and thermal properties that, when simulated at high enough resolution, can lead to the formation of cold and hot plasmoids.

3.2.2. Scale-dependent Magnetic and Kinetic Structure

We now move on to an investigation of the structure of the magnetic and velocity fields as a function of scale, which is summarized in Figure 4. The top row shows the magnetic (M(k) ≡ ∫dΩ_k k²〈∣ B ( k )∣²〉/2) and velocity (E(k) ≡∫dΩ_k k²〈∣ v ( k )∣²〉/2) power spectra compensated by k^5/3 for the turbulent thermal instability (left; E(k) of a nonmagnetized simulation shown in green) and isothermal driven turbulence simulations (right). On large scales (small k) in the turbulent thermal instability simulation, E(k) > M(k), but M(k) is nearly flat, while E(k) ∝ k^−5/3, which leads to a reversal at intermediate scales (4 ≲ kL/2π ≲ 32), in which E(k) < M(k). Starting on intermediate scales, there is a broad range of k over which the M(k) steepens to approximately follow ∝ k^−9/5, and the E(k) flattens to follow ∝ k^−4/3. This transition and change in slopes occurs at kL/2π ≈ 12. On small scales (kL/2π ≳ 32), E(k) again dominates M(k). These general trends and slopes are well converged with numerical resolution down to 16 times lower resolution (Δx = L/128) than our fiducial simulation (Δx = L/2048; see Appendix A). Most of these same trends were found in the similar (albeit significantly lower-resolution) turbulent thermally unstable MHD simulations presented by Kritsuk et al. (2017), with the notable exception of the dominance of E(k) at large k, which we ascribe to slower cooling relative to turbulent mixing in their simulations. The slope of the power spectra in simulations with additional physics, such as those presented in Padoan et al. (2016) that include self-gravity and SNe, exhibit somewhat steeper E(k) and flatter M(k). We remind the reader that these simulations are highly idealized and are not attempting to reproduce specific observations, which tend to show somewhat different slopes than what arises in our simulations (e.g., Chepurnov 1998; Lee & Lee 2020).

Zoom In Zoom Out Reset image size

In contrast, in the isothermal simulation, M(k) > E(k) on all scales, although they do approach equipartition on the largest (kL/2π ∼ 1) and smallest (kL/2π ≲ L/16Δx) scales. The velocity spectrum follows a shallower ∝ k^−4/3 from kL/2π ∼ 1 down to ∼100, which is in agreement with other similar simulations (e.g., Grete et al. 2021). The magnetic spectrum, on the other hand, peaks around kL/2π ∼ 4 and then follows a relatively steep ∝ k^−9/5 slope down to the dissipative scale. These slopes, which are robust to changes in resolution, agree very well with the slopes of the thermally unstable simulations at kL/2π ≳ 12.

Comparing the spectra of the thermally unstable and isothermal simulations highlights the enhanced kinetic energy at small scales in the thermally unstable simulations. This cannot be due to magnetic reconnection alone, since both exhibit plasmoid-mediated magnetic reconnection. Furthermore, with our resolution, Δx = L/2048, only the largest current sheets are resolved, and the scales on which plasmoids are forming are somewhat impacted by numerical dissipation. Therefore, we cannot say that at yet higher resolution, plasmoid-mediated magnetic reconnection would or would not change the power spectra, as is predicted for turbulence with a strong mean field (Boldyrev & Loureiro 2017; Dong et al. 2022). Comparing to the spectrum of a nonmagnetized turbulent thermally unstable simulation, which also exhibits a similar large k peak in E(k) (green line in top left panel of Figure 4), supports the notion that strong radiative cooling is driving small-scale motions in the thermally unstable simulations. Further evidence supporting this picture comes from the fact that as we reduced the strength of cooling in our simulations (not shown), the magnitude of the small-scale (large k) peak in the velocity spectra also decreased. The spectra in both the magnetized and nonmagnetized turbulent thermally unstable simulations change slopes at k ∼ 10–20, which is a reflection of the imprint of cooling on the flows.

To further understand the structure of the magnetic and turbulent velocity fields, we turn to the structure functions shown in the middle row of Figure 4. They encode the scale-dependent structure and anisotropy with respect to the local magnetic field direction. We adopt the two-point second-order structure function that is defined for the magnetic field to be SF₂[ B ]( ℓ ) ≡ 〈∣ B ( x ) − B ( x + ℓ )∣²〉 and equivalently for the velocity field. In a sufficiently magnetized plasma, these structure functions can show a dependence on the degree of alignment between the separation vector ℓ and the local mean magnetic field, which we define to be ${\overline{{\boldsymbol{B}}}}_{{\ell }}\equiv ({\boldsymbol{B}}({\boldsymbol{x}})\,+{\boldsymbol{B}}({\boldsymbol{x}}+{\boldsymbol{\ell }}))/2$ (Cho & Vishniac 2000). We therefore consider the parallel and perpendicular structure functions, which are defined to be ${\mathrm{SF}}_{2}[{\boldsymbol{B}}]({{\ell }}_{\parallel })\equiv \langle | {\boldsymbol{B}}({\boldsymbol{x}})-{\boldsymbol{B}}({\boldsymbol{x}}+{\boldsymbol{\ell }}){| }^{2};0\leqslant {\theta }_{{\ell }}\lt \tfrac{\pi }{36}\rangle$ and ${\mathrm{SF}}_{2}[{\boldsymbol{B}}]({{\ell }}_{\perp })\equiv \langle | {\boldsymbol{B}}({\boldsymbol{x}})-{\boldsymbol{B}}({\boldsymbol{x}}+{\boldsymbol{\ell }}){| }^{2};\tfrac{17\pi }{36}\leqslant {\theta }_{{\ell }}\lt \tfrac{\pi }{2}\rangle$, where ${\theta }_{{\ell }}\,\equiv \arccos | {\overline{{\boldsymbol{B}}}}_{{\ell }}\cdot {\boldsymbol{\ell }}/{\overline{B}}_{{\ell }}{\ell }|$ is the angle between the displacement vector and the mean magnetic field direction (cf. St-Onge et al. 2020). The average difference in the magnetic field and velocity at a given distance perpendicular to the local mean magnetic field is significantly smaller than the difference in the same distance parallel to the local mean magnetic field in both the thermally unstable and isothermal simulations. This is in line with the notion that it is easier to perturb the field and the flow perpendicular to the local mean magnetic field than it is parallel.

We can take this one step further by investigating the scale-dependent anisotropy of both the magnetic field and velocity fluctuations, which are shown in the bottom row of Figure 4. This is calculated by equating the parallel and perpendicular structure functions, which yields the parallel length scale ℓ_∥ as a function of the perpendicular length scale ℓ_⊥ (or vice versa) for both the magnetic fields (ℓ_{∥,δ B} ) and the velocity (ℓ_{∥,δ v} ) fluctuations (e.g., Cho & Vishniac 2000; Maron & Goldreich 2001). In the isothermal simulation, when ℓ_⊥ ≳ 8Δx, the magnetic field and velocity fluctuations are increasingly elongated along the local mean magnetic field with decreasing scale following the predicted ${{\ell }}_{\parallel }\propto {{\ell }}_{\perp }^{2/3}$ (Goldreich & Sridhar 1995; Lazarian & Vishniac 1999). On the other hand, in the thermally unstable simulation, only the velocity fluctuations follow this scaling, while the magnetic fluctuations instead have ℓ_{∥,δ B} scaling linearly or even superlinearly with ℓ_{⊥,δ B} . This change in the anisotropy scaling must come about due to the cooling-inducing flows that are sufficiently strong to rearrange the field. Note that there is a transition to isotropy at small scales due to numerical dissipation.

We end this section with a brief look at the properties of the cold clumps that form in our fiducial turbulent thermally unstable simulation. We identify cold clumps in the simulation using a pressure–temperature cut that selects gas that is close to the cold thermal equilibrium ¹⁵ (see top left panel of Figure 3). The top left panel of Figure 5 shows the clump mass, M_cl, distribution, and the bottom left panel shows the amount of mass contained in the cold clumps per logarithmic bin in clump mass. Over nearly 4 orders of magnitude in clump masses, the distribution closely follows a Zipf's law–like distribution with an index of −1, which has commonly been observed in other turbulent multiphase systems (see Gronke et al. 2022, for a recent investigation of the origin of this distribution). This distribution yields nearly equal mass per logarithmic bin in clump mass. In terms of total mass, however, the distribution is dominated by a single large clump that contains >30% of all of the mass in the domain. The amount of mass contained in this single large cloud is dependent on the strength of cooling. In similar simulations with weaker cooling relative to turbulent mixing, the amount of mass in the largest cloud is reduced dramatically, while the slope of the main distribution is nearly unchanged.

Zoom In Zoom Out Reset image size

In the top right panel of Figure 5, we show the intra- and interclump velocity distributions as a function of clump size and separation, respectively. The intraclump velocity (defined internal velocity standard deviation v_cl,std) versus clump size distribution is shown in pink along with the median relation. The distribution closely follows a Larson's law–like scaling with ${v}_{\mathrm{cl},\mathrm{std}}\propto {r}_{\mathrm{cl}}^{1/2}$ (Larson 1981). This is also consistent with the scaling for supersonic Burgers turbulence, as expected, since v_cl,std > c_s,cl ≈ 0.1v₀ for all clumps. The interclump velocity (difference between two clumps, the mean of which is the first-order clump velocity structure function) follows a somewhat flatter relation with 〈∣δ v_cl(ℓ)∣〉 ∝ ℓ^1/3, which follows the prediction for Kolmogorov-like turbulence and is consistent with recent observations (Ha et al. 2022). The difference between these two scalings points to the cold clouds being entrained in a background sub-/trans-sonic turbulent flow while having supersonic internal motions.

In the bottom right panel of Figure 5, we show the elongation distribution as a function of clump size. To quantify the elongation, we adopt ${(2{I}_{3}/{I}_{1}-1)}^{1/2}$ as a metric, where I₃ and I₁ are the largest and smallest moments of inertia of the clump, which, for a prolate (cigar-shaped) ellipsoid, equals the ratio of the length of the longest axis to the shortest. For a sphere, this takes a value of 1, and for a thin disk, it reaches 5^1/2. Although this quantity is imperfect for distinguishing between marginal cases, it clearly demonstrates that the cold clumps in our simulation tend to be highly elongated structures.

We present the highest-resolution (to our knowledge) MHD simulation of a thermally unstable turbulent medium to date. We show that large current sheets unstable to plasmoid-mediated reconnection form regularly throughout the volume (the simulations in previous related studies, such as Padoan et al. 2016 and Kritsuk et al. 2017, had significantly lower resolution, which prevented them from resolving the critical Lundquist number, S_crit, and, therefore, inhibited the plasmoid instability). The current sheets that reconnect and give rise to plasmoids develop in three distinct environments: within the coldest phase, at the asymmetric interface of the cold and warm phases, and within the volume-filling warm phase.

We show the rich magnetothermal phase structure that develops, in which most of the mass is in the highly magnetized cold phase with β < 1 (albeit super-Alfvénic, ${{ \mathcal M }}_{A}\gt 1$), and that the regions of highest magnetic curvature span a broad range of temperatures. Furthermore, we show that thermal instabilities and magnetic fields change the scale-dependent magnetic and velocity structure. This is achieved by comparing our fiducial simulation to an isothermal simulation and a simulation without magnetic fields. The comparison highlights that the velocity spectrum dominates on both large and small scales with cooling likely being responsible for the increase in small-scale kinetic energy. Additionally, in the isothermal simulation, both the magnetic and velocity perturbations follow the predicted increase in anisotropy relative to the local mean magnetic field with decreasing scale. In the thermally unstable simulation, however, only the velocity perturbations follow this scaling, whereas the magnetic perturbations are more isotropic and maintain the same small degree of anisotropy on all scales.

Lastly, we show that the clumps of cold gas in the thermally unstable simulation mostly follow a scale-free mass distribution with equal mass per logarithmic bin in clump mass, with the exception of the largest clump, which contains most of the mass in the entire simulation. These clumps tend to be highly elongated and exhibit a size versus internal velocity relation consistent with supersonic turbulence and a relative clump distance–velocity scaling consistent with subsonic motion.

While our simulations presents a major step forward in resolving the small-scale magnetothermal properties of the ISM, they are clearly and intentionally highly idealized, which means there are important caveats that should be taken into account. Although our piecewise power-law cooling and heating source terms capture the important general trends applicable to a range of systems, they miss important subtleties essential for detailed comparison to specific systems. Likewise, our simulations do not track changes in the chemistry and ionization state of the system, which can have important observational and physical impacts. For our 3D simulations, we rely on numerical dissipation to break the magnetic flux freezing, when in reality, explicit nonideal MHD effects, such as ohmic dissipation, should be mediating the reconnection. Resolving these effects in the plasmoid-mediated regime in 3D is extremely computationally expensive and has only been achieved once in a single-phase, decaying turbulence simulation (Dong et al. 2022). We stress, however, that in Appendix C, we use 2.5D simulations to test the impact of using resolved explicit dissipation and find strong indications that the physical processes we see are robust. Similarly, we do not include explicit conduction and viscosity in our 3D simulations, which are expected to be relevant (e.g., Cowie & McKee 1977) and highly anisotropic (e.g., Choi & Stone 2012) and may alter the phase structure and flows in important ways. The SNe, which we have not included in our simulations, generate a third, hottest phase (McKee & Ostriker 1977) and drive more realistic turbulence than what we assumed. It remains to be seen how reconnection proceeds in the hot phase, which will be much more subsonic, and in the presence of more compressive flows. Finally, our simulations do not include the impact of gravity, which may modify the dynamics, particularly in the densest regions, and may lead to the development of converging flows that can form current sheets. Despite these caveats, our findings are relevant to many outstanding ISM observation puzzles and may provide a useful guide to better understanding the physical processes at play in the ISM.

In the ISM, it is well known that structures exist across a vast range of scales spanning more than 10 orders of magnitude from hundreds of parsecs down to hundreds of kilometers (Armstrong et al. 1995). On the few to ∼10⁴ au scale, there is ample observational evidence for tiny-scale atomic and ionized structures (TSAS and TSIS, respectively) that probe the poorly understood dissipation processes (for a useful overview, see Stanimirović & Zweibel 2018). These structures manifest as extremely small or fast spatial or temporal variations in H i absorption (TSAS) and as extreme scattering events and interstellar scintillation in radio observations of both quasars and pulsars (TSIS). Observational models favor either very high density spherical structures or highly elongated structures viewed lengthwise. The plasmoids in our simulation are both dense and elongated, and they are found on the smallest scales we are able to probe, which is comparable to the larger end of the observed sizes (Δx ∼ 5 × 10³ au). Plasmoids will only become more common at smaller scales if we are able to probe higher resolutions. The thermal instability–induced cold plasmoids are thus good candidates for TSAS, and the hot plasmoids are good candidates for TSIS. A careful comparison will need to be made to show this connection more concretely, but the qualitative agreement is compelling. Furthermore, the magnetic topology of the plasmoids provides a natural confinement and insulation that will tend to extend the lifetime of these small structures. Our scale separation between the outer and dissipative scales is, however, limited, so a better understanding of how these results scale is needed.

On somewhat larger scales, recent work has uncovered highly elongated thin H i structures referred to as H i fibers (Clark et al. 2014). These fibers are dense (n ∼ 14 cm⁻³) and have lengths on the order of a few parsecs and widths of ≲0.1 pc (or possibly ≪0.1 pc, since the widths are unresolved in current observations). The properties of the cold plasmoids in our simulations are highly reminiscent of the properties of the H i fibers and may provide a physical explanation for their origin. This supports a physical picture in which TSAS and H i fibers are cold plasmoids seen either lengthwise or edgewise. A detailed look (which we defer to a future work) at the properties of the cold clumps in our simulations, in particular how aligned they are with the local magnetic field, may help shed light on the nature of ISM turbulence (e.g., González-Casanova & Lazarian 2017) and far-ranging observations such as foregrounds in the cosmic microwave background (Clark et al. 2021).

Our findings have potentially major implications for the transport of CRs in the ISM. The CRs play an essential role in the evolution of galaxies (Zweibel 2017), but their impact depends strongly on their transport (Salem & Bryan 2014; Ruszkowski et al. 2017; Hopkins et al. 2020; Quataert et al. 2022a, 2022b). The CR transport is regulated by the scattering with small-scale electromagnetic fluctuations that are either self-excited (Kulsrud & Pearce 1969; Skilling 1971), preexist in the ambient turbulent medium (e.g., Chandran 2000; Yan & Lazarian 2004), or some combination of both (Xu & Lazarian 2022). Recent work has shown that both of these models make predictions that are in significant tension with existing observations (Hopkins et al. 2022; Kempski & Quataert 2022). Two aspects of our simulations represent attractive solutions to the shortcomings of the model in which CRs scatter as a result of turbulence in the ambient medium. First, the magnetic fluctuations in our thermal instability simulation do not exhibit an increase in anisotropy with decreasing scale (Figure 4). This is crucial because anisotropic eddies are very inefficient at scattering CRs (Chandran 2000). Second, plasmoids are regions of extremely high magnetic curvature, so a CR confined to move along a field line that encounters a plasmoid may efficiently scatter if its Larmor radius is comparable to the curvature radius, thus making plasmoids an attractive means to scatter CRs, as suggested by Kempski & Quataert (2022). ¹⁶ Furthermore, the size spectrum of the plasmoids may naturally lead to an energy-dependent CR diffusion coefficient in line with what is needed to reproduce Milky Way observations (Amato & Blasi 2018). We plan to investigate CR scattering in plasmoid-dominated MHD turbulence in a future work.

An additional implication of our ISM plasmoid findings is that the lifetime of the smallest structures may be significantly enhanced. The highly anisotropic nature of thermal conduction in the ISM (e.g., Spitzer 1962) means that conductive evaporation will be reduced, since the interior of a plasmoid will be well insulated due to the fact that, as can be clearly seen in Figure 1, it is almost entirely magnetically disconnected from the surrounding medium. Therefore, we expect to find clumps on scales smaller than the standard isotropic Field length (Begelman & McKee 1990). This is supported by the findings of thermal instability simulations with anisotropic thermal conduction (Sharma et al. 2010; Choi & Stone 2012; Jennings & Li 2021). Furthermore, the enshrouding magnetic fields will also limit hydrodynamic mixing of the plasmoid with the surroundings (Chandrasekhar 1961), thereby extending the lifetime and minimum size that should be expected.

More generally, our findings demonstrate that we are now getting to the stage where we can move beyond simulations of individual current sheets (e.g., Daughton et al. 2011; Huang & Bhattacharjee 2016; Zhang et al. 2021; Beg et al. 2022) and instead study the interaction between turbulence and reconnection (e.g., Schekochihin 2020; Chernoglazov et al. 2021; Nättilä & Beloborodov 2021; Comisso & Sironi 2022; Dong et al. 2022). Such studies will help disentangle the intrinsic turbulence (e.g., Kowal et al. 2017) that arises from reconnection within a layer and the extrinsically driven turbulence. This will help in understanding the longevity of plasmoids because, unlike in 2D, in 3D, they do not have a clear circular structure but do exhibit clear coherence in the third dimension. How these plasmoids/flux tubes survive will depend on the interplay of reconnection and turbulence, which will determine how susceptible they are to kinking and dissolving into the background.

This work represents the beginning of a new line of inquiry into the nature of the magnetized multiphase ISM. There remains much to be done to strengthen and extend our initial findings and connect them to observations and other physical processes. We have clearly demonstrated that the plasmoid instability and associated fast magnetic reconnection is a common occurrence in the ISM and can be induced by thermal instabilities. Furthermore, given the simplicity of our approach, our findings are likely to apply to a broad range of systems beyond the ISM, such as the circumgalactic medium, intracluster medium, and stellar coronae.