© The Authors 2026

1 Introduction

Understanding dust formation and destruction is key for planetary, stellar, and (extra-)galactic astronomy alike (cf. Draine 2003; Micelotta et al. 2018). Asymptotic giant branch (AGB) stars and supernova remnants (SNRs) are known to provide favorable conditions for dust production. Formed dust grains are also hypothesized to be able to grow efficiently in the interstellar medium (ISM; e.g., Draine 2009; Hirashita 2012; Mattsson & Andersen 2012; Zhukovska et al. 2016, 2018), significantly increasing the ISM dust budget over time. Even so, the relative importance of these dust formation and growth mechanisms across redshifts and the detailed underlying physical processes remain uncertain (cf. Schneider & Maiolino 2024).

Asymptotic giant branch stars play an important role in the dust budget of local galaxies (e.g., Matsuura et al. 2009; Olofsson et al. 2010; Schneider et al. 2014; Dell'Agli et al. 2015). However, at high redshifts, where several galaxies have been observed to contain high dust masses (e.g., Pettini et al. 1994; Bertoldi et al. 2003; Watson et al. 2015; Bowler et al. 2022; Sommovigo et al. 2022; Akins et al. 2023; Witstok et al. 2023), the dust needs to be produced at a much faster rate than typical AGB dust production channels can provide. Asymptotic giant branch stars are the end evolutionary stages of low- and intermediate-mass stars with masses of less than 8 M_⊙, and thus take several tens of million years to start contributing to the dust budget (e.g., Vassiliadis & Wood 1993; Valiante et al. 2009). If the initial mass function is bottom-heavy then this contribution will take even longer with solar mass stars taking billions of years to start their dust production. Therefore, the majority of the stars could not aid in dust production at high redshifts. Past studies show that AGB stars in early galaxies need at least 300 Myr to dominate the dust production (Di Criscienzo et al. 2013; Valiante et al. 2017), which is too long to account for the earliest dusty galaxies.

In recent decades, core-collapse supernovae (CCSNe) have been found to be a promising candidate for dust production in the early universe (e.g., Morgan & Edmunds 2003; Dwek et al. 2007; Gall et al. 2011a,b,c; Leśniewska & Michałowski 2019; Pozzi et al. 2021; Dayal et al. 2022; Palla et al. 2024) because they originate from stars more massive than 8 M_⊙, which have lifetimes of only a few to several tens of million years. Observations in the local universe have shown that, within the first few decades after the SN explosion, 1 M_⊙ of dust can be formed in the SNR ejecta (e.g., Lucy et al. 1989; Dunne et al. 2009; Matsuura et al. 2011; Gomez et al. 2012; Temim et al. 2012; Bevan & Barlow 2016; De Looze et al. 2017, 2019; Niculescu-Duvaz et al. 2022). However, SNRs do not only produce dust; they are also able to destroy a large amount of pre-existent dust, due to their high-velocity forward shock (e.g., Lakićević et al. 2015; Priestley et al. 2021, 2022). This shock induces gas-grain collisions (sputtering) and grain–grain collisions (shattering), thus reducing the ISM dust grain sizes or destroying them completely (e.g., Barlow 1978; Draine & Salpeter 1979; Jones et al. 1994, 1996). Furthermore, the freshly produced ejecta dust needs to survive the reverse shock that forms in a subsample of SNRs with heavy pre-SN mass-loss or high-density environment as a result of pressure build-up behind the forward shock. This could result in a significantly lower dust contribution from SNRs to the ISM than the commonly assumed 1 M_⊙ (e.g., Kirchschlager et al. 2024b). When forward and reverse shock destruction are taken into account, it becomes unclear whether SNe can indeed make a net contribution to the dust budget.

To resolve this uncertainty, it is crucial to constrain the dust destruction rate of the SNR forward shock, which encounters several tens of solar masses of dust during its evolution (e.g., Kirchschlager et al. 2022, 2024a). Many analytical and numerical forward shock dust destruction studies have been carried out (e.g., Jones et al. 1994, 1996; Nozawa et al. 2006; Bocchio et al. 2014; Slavin et al. 2015; Hu et al. 2019; Martínez-González et al. 2019; Kirchschlager et al. 2022, 2024a; Vasiliev & Shchekinov 2024; Dedikov & Vasiliev 2025) and predict a wide range of outcomes between almost total survival (Martínez-González et al. 2019, 0.03% destruction rate) to total destruction (e.g., Bocchio et al. 2014, up to 100% destruction rate). Thus, SNRs can act as dust producers or destroyers depending on the considered forward shock study and the similarly uncertain reverse shock destruction rate (e.g., Bianchi & Schneider 2007; Nozawa et al. 2010; Silvia et al. 2010, 2012; Bocchio et al. 2016; Micelotta et al. 2016; Slavin et al. 2020; Kirchschlager et al. 2019, 2023, 2024b; Vasiliev & Shchekinov 2024; Otaki et al. 2025).

This large disparity originates from the different initial conditions, and also from the various ways of implementing dust into codes (Slavin 2024) and the considered dust transport and destruction mechanisms. In particular, grain–grain collisions were expected to only have a small impact on the dust destruction, and thus were often neglected (e.g., Nozawa et al. 2006; Bianchi & Schneider 2007; Nozawa et al. 2010; Silvia et al. 2010, 2012; Micelotta et al. 2016; Hu et al. 2019; Martínez-González et al. 2019; Slavin et al. 2020; Vasiliev & Shchekinov 2024; Dedikov & Vasiliev 2025). However, Slavin et al. (2015) and Kirchschlager et al. (2019) found that the redistribution of grain sizes can work synergistically with other dust destruction mechanisms, and thus lead to large amounts of destroyed dust masses. Furthermore, a magnetic field was also often not taken into account (e.g., Nozawa et al. 2006; Bianchi & Schneider 2007; Nozawa et al. 2010; Silvia et al. 2010, Silvia et al. 2012; Martínez-González et al. 2019; Kirchschlager et al. 2020; Slavin et al. 2020; Kirchschlager et al. 2022; Vasiliev & Shchekinov 2024; Dedikov & Vasiliev 2025), although dust grains are charged and thus influenced by it. It was found that the effects of the Lorentz force can significantly change the dust destruction efficiency in SNRs provided that high enough magnetic field strengths are present (e.g., Shull 1978; Draine & Salpeter 1979; Fry et al. 2020; Kirchschlager et al. 2023, 2024a).

In addition to the dust destruction mechanisms considered, the properties of the ambient medium surrounding the SN are expected to affect how much dust is destroyed by the forward shock. Most works consider a simple homogeneous ISM in which the SNR forward shock evolves spherically. There are just a few studies we know of that include an inhomogeneous ISM (Hu et al. 2019; Kirchschlager et al. 2022, 2024a; Vasiliev & Shchekinov 2024; Dedikov & Vasiliev 2025), in which the SNR evolution becomes non-spherical because the forward shock is slowed down significantly in the dense cold filaments and can travel faster through voids. Since kinetic sputtering mainly depends on the shock velocity and thermal sputtering on the temperature, it is theorized that the dense and dust-rich filaments can shield dust from destruction. Simple small density fluctuations, as considered by Kirchschlager et al. (2022) and Vasiliev & Shchekinov (2024), show no agreement as the former finds a slightly higher dust destruction compared to the homogeneous case and the latter does not provide such a comparison. A more comprehensive study of different density inhomogeneities of Dedikov & Vasiliev (2025) draws a clearer picture. They find that a clumpy ISM can reduce the dust destruction by up to ∼20%. More complex turbulence, as studied by Hu et al. (2019) and Kirchschlager et al. (2024a), shows the same trend. These studies drive their turbulence by setting off SNe every 1–2 Myr at random locations in an initially homogeneous medium and letting the system evolve for 0.4 Gyr or to a statistical steady state, respectively. Since the turbulence is only driven by SNe, it lacks small-scale structures that influence the SNR evolution, particularly in the early stages, and are therefore expected to affect the dust destruction by the forward shock. Additionally, other feedback mechanisms such as winds and photoionization, or gravity, which would affect the turbulent energy cascade, are neglected. Thus, although this type of turbulence ensures a purely physical origin, it often has different statistical properties compared to the observed multi-scale turbulence (e.g., Federrath et al. 2010) whose main driving mechanisms are still under debate (e.g., Miesch & Bally 1994; Mac Low & Klessen 2004; Brunt et al. 2009; Burkhart & Mocz 2019; Sartorio et al. 2021).

Turbulence can also be driven in simulations by a forcing term which can be employed by an Ornstein–Uhlenbeck process to insert Fourier waves with random propagation directions, and at random locations and scales that go down to the Jeans length. These shocks then interact until the density and velocity probability density functions (PDFs) of the fluid quantities converge. Since small and large energy waves are injected, different kinds of modes will occur in the resulting turbulence. With this method, the initial conditions can be fine-tuned to closely resemble the observed turbulence with simulations.

We used this approach to drive purely solenoidal turbulence in a highly resolved 3D simulation box, in order to resemble the density distribution seen in regions of the ISM, including small-scale structures. Afterward, an SN was injected into the driven inhomogeneities and evolved for 10 kyr. During this period, the SNR forward shock is the most destructive due to its high initial kinetic and thermal energy. Furthermore, it already encounters a significant amount of dust due to its rapid expansion. After 10 kyr, the SNR will continue expanding and destroying dust, but with a lower destruction rate. By including thermal and kinetic sputtering, grain–grain collisions, and a magnetic field, we ensured that all known relevant dust transport and destruction mechanisms were considered. In this way, we aim to improve constraints on the forward shock dust destruction efficiency during the most destructive phase of SNR evolution.

We describe the methods in greater detail in Sect. 2, followed by Sect. 3, which shows and analyzes the results of our magne-tohydrodynamic (MHD) and dust simulations. The findings are discussed and compared to previous studies in Sect. 4, and the paper concludes in Sect. 5.

2 Methods

For this study we performed 3D hydrodynamic simulations with the widely used quasi-Lagrangian code Arepo (Springel 2010), which uses a moving mesh based on Voronoi tessellation. We drove turbulence using an Ornstein–Uhlenbeck process (Sect. 2.1), and considered different average ISM densities of 1 cm⁻³ (low density), 10 cm⁻³ (intermediate density), or 100 cm⁻³ (high density) and a turbulent Mach number of 0 (homogeneous), 0.3 (subsonic), 1 (transonic), or 3 (supersonic). Afterward, an SN explosion was injected as a large amount of thermal energy, distributed over a few cells at the center of our simulation domain. The SNR evolution was followed for 10 kyr in 3D MHD simulations with the same code (Sect. 2.2). Since Arepo does not include detailed dust transport and destruction mechanisms, we used the dust post-processing code Paperboats (Kirchschlager et al. 2019) that assumes an initial dust distribution in the turbulent box and then calculates the dust transport and destruction based on gas and plasma drag, Lorentz acceleration, thermal and kinetic sputtering, grain–grain collisions, and accretion (Sect. 2.3). Using this workflow, we are able to estimate the dust destruction efficiency of SNR forward shocks in different ISM environments.

2.1 Turbulence driving

The 3D turbulence simulations were carried out purely hydrody-namically by using the turbulence driving within Arepo (Mocz et al. 2017; Mocz & Burkhart 2018), which makes use of the Ornstein–Uhlenbeck process. We started the simulations with a homogeneous medium in which we drove the turbulence until it was in a statistically stationary state, and thus until the density and velocity PDF, and the power spectra converged. This method of driving let us generate a turbulence that closely resembles the observed ISM, as a visual comparison in Fig. 1 shows. There, the dust column density of one of our driven turbulence snapshots is qualitatively compared to a dust map of the Polaris flare from Herschel/SPIRE observations (Miville-Deschênes et al. 2010)¹ to show that the morphological features and scales are similar.

In the turbulence simulations the mass of every Voronoi cell was kept approximately the same, which induced refinement but did not change the total number of cells. We did not include additional refinement criteria for the turbulence driving. The box itself was 20 or 40 pc large in each coordinate direction, depending on the used ISM density and turbulence strength. These box sizes were large enough so that the SNR could evolve for up to 10 kyr without reaching the boundaries in every simulation. The boundaries were chosen to be periodic to ensure that no dust escapes from the box in post-processing. To be able to resolve filaments in the ISM that are observed to have a width of ~0.1 pc (Arzoumanian et al. 2011), we used 512³ cells in total. Thus, the resolution in filaments was up to 0.008 pc (0.016 pc) whereas in voids it was 0.15 pc (0.3 pc) for the 20 pc (40 pc) box.

The studied parameter sets include different turbulence strengths, defined by the turbulent Mach number ℳ = σ_v/c that is dependent on the velocity dispersion σ_v and theMsound speed c (Beattie et al. 2019). This Mach number determines the statistical distribution of the density and velocity, as can be seen in Fig. 2 with the density PDF. Lower turbulent Mach numbers only show slight density variations (e.g., ℳ = 0.3, blue line), whereas higher turbulent Mach numbers result in dense filaments and empty voids (e.g., ℳ = 3, pink line). Because the turbulent Mach number of the interstellar and circumstellar medium around stars can differ significantly depending on the exact location of the star, we took different Mach numbers into account – from subsonic (ℳ = 0.3) to transonic (ℳ = 1) and supersonic (ℳ = 1), and compared them to the homogeneous case (ℳ = 0). This Mach number range reflects typical conditions in the ISM, with observed values ranging from subsonic (a few tenths) to supersonic (ℳ = 0.17–8, e.g., Heyer et al. 2006; Ginsburg et al. 2013; Schneider et al. 2013; Marchal & Miville-Deschênes 2021).

We also included different ISM densities since SNe can occur in media of various densities. Stars can remain in their birth molecular cloud throughout their lifetime, where they can explode close to dense filaments or in a medium that is emptied by previous SNRs, or can drift away from the cloud to lower density media (Mayker Chen et al. 2023; Sarbadhicary et al. 2023). Thus, we chose to study different average ISM number densities of n = 1, 10 and 100 cm⁻³ in addition to the set of studied turbulent Mach numbers to comprehensively cover various realistic SNR surroundings. The column density of such a turbulence simulation box can be seen in Fig. 1 (left). Furthermore, slices of all simulations that were carried out in this study are visu-alized in Appendix A, in Figs. A.1, A.2, and A.3 for the low-, intermediate-, and high-density environments, respectively.

The initial ISM temperature in these simulations was chosen to be 10⁴ K. In the interest of maintaining a smaller computational requirement, we did not include magnetic fields, cooling or self-gravity in the turbulence driving. We did include the first two a posteriori in the SN explosion simulations, as laid out later in this section.

Since a magnetic field can play an important role in dust transport and destruction, we added it to the converged turbulent system with the same direction as the velocity vectors and a strength proportional to the number density n with B = bn^1/2 (Crutcher 1999), including the constant factor b. This factor was estimated to b = 0.93 ± 0.87 µG cm^3/2 by averaging several observations of the magnetic field strength in dense filaments (Curran & Chrysostomou 2007; Chapman et al. 2011; Pattle et al. 2018; Coudé et al. 2019; Liu et al. 2019; Wang et al. 2019; Arzoumanian et al. 2021; Eswaraiah et al. 2021; Könyves et al. 2021; Lyo et al. 2021; Ngoc et al. 2021; Pattle et al. 2021; Ching et al. 2022; Kwon et al. 2022; Karoly et al. 2023; Ward-Thompson et al. 2023; Wang et al. 2024). We tested the impact of the direction of the magnetic field by comparing to a simulation with a 90^◦ rotated magnetic field. The dust destruction results only differed by less than 0.1% so that the magnetic field strength seems to have a larger impact than its direction for our models.

Fig. 1 Structural comparison of the dust distribution in a turbulence snapshot of our intermediate-density (10 cm−3) supersonic turbulence (ℳ = 3) simulation (left) and an observed SPIRE 250 µm dust map from Herschel observations of the Polaris Flare (right, originally observed by Miville-Deschênes et al. 2010)1, which has a Mach number of ℳ = 2–10 (Beattie et al. 2019). The physical length scale of 20 × 20 pc of the observed ISM is based on the distance of 350 pc to the Polaris Flare (Schlafly et al. 2014; Panopoulou et al. 2022).

Fig. 2 Density PDFs of the different turbulent Mach numbers ℳ in the converged low-density (1 cm−3) turbulence simulations. With a higher turbulent Mach number, cell densities (ρ) can have a larger deviation from the average density <ρ>.

2.2 SN explosion and SNR evolution

The simulation output of Sect. 2.1 was used as the starting point for our MHD SNR simulations. In these simulations, we set off an SN with an explosion energy of 10⁵¹ erg using the SN feedback method of Arepo that injects thermal energy into 40 central neighboring Voronoi cells, corresponding to a radius of 0.08–0.21 pc in our simulations. The SNR evolution was followed for 10 kyr. We show the final snapshots in Fig. 3 (see also Appendix A for the full evolution of the simulations, including movies). In addition to the mass refinement of Arepo, we implemented refinement with respect to the pressure gradient to follow the SNR forward shock in greater detail, since it is the key region where dust destruction happens. Cell splitting and merging were enabled, so that the total number of cells almost doubled throughout the simulation because the SNR grew. Thus, the average distance between neighboring cell centers at the end of the simulations was of the order of 0.030 and 0.060 pc for the 20 and 40 pc box size, respectively.

In the SNR simulations, radiative cooling is considered. We assumed collisional ionization equilibrium and electron-ion equipartition as we do not follow electron and ion temperatures separately, or ionization states in detail. The used cooling implementation in the code depends on whether the gas has been shocked, which was determined by its temperature. Unshocked ISM gas had an initial temperature of 10⁴ K and cooled according to detailed molecular and atomic cooling functions, assuming solar metallicity and a simplified treatment of gas-phase molecular chemistry (Glover & Clark 2012; Bertram et al. 2016). High-density filaments cool down efficiently before they are shocked, whereas voids remain at higher temperatures, due to their low density. We identified shocked gas in cells where the temperature exceeded 10⁴ K. We assumed that this gas mixed with SN ejecta, with the latter dominating the cooling process. Since our simulations did not include explicitly injected ejecta, the shocked ISM effectively took on the role of the ejecta. To emulate the conditions in the ejecta of Cassiopeia A (Cas A), which we used as a reference, we adopted oxygen-dominated cooling based on the CHIANTI cooling function (Dere et al. 1997; Del Zanna et al. 2015), following the approach of Kirchschlager et al. (2019). This method was chosen because applying the ISM cooling prescription to the ejecta resulted in a faster expansion of the SNR, due to insufficient cooling, leading to a divergence from the observed evolution of SNRs similar to Cas A. This approximation is considered reliable only during the early stages of the SN, when the shocked mass remains comparable to the expected ejecta mass. Because the swept-up ISM material is heated beyond the ISM cooling threshold, it was cooled according to the ejecta prescription, resulting in a slight overestimation of cooling at the shock front, particularly at late times in the high-density simulations.

We did not consider thermal conduction. Its timescale can be calculated by (e.g., Orlando et al. 2005) $$t_{\text{cond}}=\frac{7}{2(\gamma -1)}\frac{n{k}_B{l}^2}{\kappa (T)},$$(1) where l is the length scale of interest and $\kappa \left(T\right)=5.6\times {10}^{-7}{T}^{-\frac{5}{2}}$ erg s⁻¹ K⁻¹ cm⁻¹ is the thermal conductivity. The short est timescale of ∼13 kyr is reached at the contact of the thin shell and the SNR interior of the high ISM density simulations. However, the forward shock is expanding on shorter timescales in our simulations so that thermal conduction can be neglected. Although at very early times of a few hundred years, the thermal conduction timescale between the forward shock and the ISM is on the order of only a few thousand years, the ISM is heated more rapidly by the expansion of the forward shock than by conduction at this early stage.

Having set all initial parameters, the SN was set off in the driven turbulent ISM box. Because stars usually excavate their surroundings with stellar winds and radiative feedback, most of them explode in a low-density cavity (e.g., Lamers & Cassinelli 1999; Rogers & Pittard 2013; Fichtner et al. 2024). Modeling the stellar feedback mechanisms in detail is beyond the scope of this paper; however, we mimicked the impact of pre-SN feedback by always placing the SN in a void of the turbulence, also ensuring that all simulations are easily comparable. The void was manually chosen by shifting the 3D turbulence box to a position with a large void in the center. Therefore, the forward shock encountered filaments only after several hundred to a few thousand years of evolution for the higher Mach number simulations as is also the case for SNe that explode in a wind-driven bubble. The voids in our strong turbulence simulations were roughly 10–20 pc large, which is comparable to the observed and predicted wind-driven cavity sizes of 10–100 pc of massive stars, strongly dependent on the stellar mass and ISM density (McCray 1983; Georgy et al. 2013). If the SNe were to explode closer to dense filamentary structures instead of voids, the filaments would be encountered at higher velocity, resulting in more efficient dust destruction in the early SNR stages. This is discussed in Appendix B.

After evolving the SNR for 10 kyr in MHD simulations, we studied the dust dynamics and destruction in post-processing simulations using Paperboats since Arepo does not include detailed dust processing. These dust simulations are run from snapshots of the SNR simulation, meaning that their time steps are larger than the hydrodynamic time steps. Since the forward shock is faster in the beginning and slows down significantly (from v > 10⁴ km s⁻¹ to v < 400 km s⁻¹), the first few hundred years of evolution must have a higher temporal resolution than the last few hundred years. Thus, we took 512 snapshots of the MHD simulations at logarithmically spaced time steps to feed the post-processing simulations, resulting in time steps of 0.15–120 yr for the dust simulations.

Fig. 3 Number densities of the z = 0 slice of all simulations carried out in this paper at 10 kyr. This includes different average ISM number densities (n) and turbulent Mach numbers (ℳ). The colorbar scales are consistent in each of the average ISM number density cases but the box sizes can vary between 20 or 40 pc.

2.3 Dust dynamics

We post-processed the 512 snapshots from each MHD simulation with the Paperboats code (Kirchschlager et al. 2019) in order to determine the dust dynamics. This code was developed to study the dust destruction of the reverse shock in dense ejecta clumps (Kirchschlager et al. 2019, 2020, 2023, 2024b) and was later also used for forward shock dust destruction (Kirchschlager et al. 2022, 2024a). It simulates the dust by introducing a certain dust mass according to an initially fixed gas-to-dust mass ratio parameter that we adopted to be 100, which is close to the value assumed in dust models characterizing the diffuse ISM in the Milky Way (e.g., Draine & Li 2007; Compiègne et al. 2011; Jones et al. 2013). Based on the amount of dust, its properties, and the fluid conditions from the MHD simulations, Paperboats can determine how the dust will be affected by the gas and other dust grains via gas and plasma drag, Lorentz acceleration, thermal and kinetic sputtering, effects of grain–grain collisions (fragmentation and vaporization), and accretion in the following MHD snapshots. The rotational disruption of dust grains by mechanical torques (METD) or radiation (RATD, Hoang & Lee 2020) is not included in Paperboats. These effects were recently evaluated in the context of SNR dust destruction (Martínez-González 2025), and were found to only have a minor impact on the overall dust destruction. However, this could change for our model setup and will be part of analysis in the near future.

Initially, the dust was given the same velocity as the gas. As such, dust in the central cells would not be given an initial velocity so that it would not move out of the center. Therefore, we did not start with the first snapshot of the MHD simulation but from snapshot 10, which relates to only a few years after the explosion. The dust was then added with the dust grain-size distribution of Mathis, Rumpl and Nordsieck (MRN; Mathis et al. 1977) from 5 to 250 nm, and thus a power-law distribution a^−γ with γ = 3.5. To account for grain destruction and accretion, the grain sizes in our simulations ranged from 0.6 to 350 nm, which were separated into 20 logarithmically spaced bins. In total, we performed 4 dust destruction simulations for each SNR simulation, considering only sputtering, or sputtering and grain–grain collisions to assess the impact of grain–grain collisions, and either silicate or carbonaceous dust to study how the dust destruction efficiency varies for different grain species. More detailed information about the code implementation of Paperboats can be found in Kirchschlager et al. (2019, 2020, 2023).

The gas-dust interactions simulated by Paperboats assume that the ISM gas consists of a single atomic species. We modeled the gas as being dominated by hydrogen atoms and therefore used a mean molecular weight of µ = 1 for our dust simulations. As a result, hydrogen atoms must also dominate in the region where ISM gas and ejecta gas mix. This approximation becomes increasingly accurate as the SNR evolves since the mass of ISM gas in the SNR interior increases over time. In our simulations, the remnant encountered at least 150 M_⊙ of ISM material during 10 kyr, which is significantly more than typical SNR ejecta masses (e.g., ∼4 M_⊙ for Cas A). Thus, hydrogen dominated the total mass in the SNR interior for most of the evolution, although the cooling was still dominated by oxygen. For consistency, we used this simplified molecular weight for our density calculations in both the MHD and dust simulations. However, solar abundances were generally assumed for the rest of the MHD calculations.

Since Paperboats does not support a moving mesh yet, we converted the Voronoi grid of Arepo to a Cartesian grid that could be used by Paperboats. The Cartesian grid had 900 cells in each coordinate direction, almost twice as high as the initial average distance of cell centers in the MHD simulations. The dust simulation resolution was found to be optimal for our setup and to not induce significant errors during the simulations, as also discussed in Appendix C. During the conversion from Voronoi to Cartesian grid, the total mass of the system changed due to the rearrangement of the cells. However, our chosen resolution was high enough to only change the mass by 0.2%, which had a minor impact on the dust destruction.

Furthermore, we only studied the central 2D layer because a full 3D dust calculation would have been numerically too expensive. The final amount of Cartesian grid cells was therefore 900 × 900 × 1 cells in x-, y-, and z-direction. To ensure that the investigated layer was representative of the whole 3D SNR, we performed the dust simulations on three different perpendicular planes. The dust destruction of the forward shock in these planes can differ by several per cent as shown in Appendix D, depending on the different detailed density distribution in the slice. Thus, dust destruction values in this paper always refer to the average of these three planes, unless otherwise stated.

To evaluate the dust destruction efficiency of only the encountered ISM dust instead of the dust in the full simulation box, we had to trace the region which was shocked by the SNR in all snapshots. Since the SNR is still hot and fast after 10 kyr, it was possible to follow it by using a temperature or velocity gradient threshold, which were implemented into Paperboats. The temperature threshold worked best for the subsonic and homogeneous cases and traced all cells which are or were above the temperature of 1.1 × 10⁴ K in any of the MHD snapshots. This threshold was chosen to correctly follow the SNR in all time steps, even in the first few where the ISM is at 10⁴ K but cools down rapidly. For the trans- and supersonic simulations, voids did not cool quickly but can be compressed and thus heated by the slight diffusion of the filaments. Therefore, we increased the temperature threshold to 3 × 10⁴ K for the trans- and supersonic simulations with a low ISM density. For the intermediate- and high-density cases, however, it is possible that, at a given snapshot, shocked gas cells cool down to temperatures lower than the ones retained in the unshocked gas in voids. Hence, we used a v∆v_i threshold that traced all cells that are or were above it at any given time step for the trans- and supersonic simulations of these cases, where ∆v_i is the velocity difference of neighboring cells in the coordinate direction i ∈ {x, y, z}. Because this threshold is resolution dependent and the forward shock evolves differently depending on the ISM density and Mach number, we adjusted the v∆v_i threshold for each of the four simulations for which it was used. This was done by carefully fine-tuning the threshold until no ISM gas cells but all SNR forward shock cells were traced. Since we tracked all cells that are or once were above this threshold, we finally followed the whole SNR interior over time. We manually checked that the SNR evolution was followed correctly for all of our snapshots. This meant that we could follow the dust destruction and grain-size distributions of the SNR shocked region more easily, noting that mixed cells at the shock front may be traced too early or too late. However, cells at the boundary were only a small fraction of all traced SNR cells so that the impact on the total dust destruction rate is negligible.

3 Results

In this section we discuss the impact of the different initial ISM densities and turbulent Mach numbers on the SNR evolution, and analyze how they influence the dust destruction efficiency, considering different dust species and destruction mechanisms. The gas density distribution of all SNR simulations at 10 kyr is shown in Fig. 3 (see Appendix A for the full evolution, including movies). Depending on the considered parameter set, the total dust destruction lies between 27–92%, corresponding to 0.85–11.0 M_⊙ of destroyed dust in the first 10 kyr of SNR evolution. The given destroyed dust masses are not the dust masses in a 2D plane but an extrapolation of the 2D destruction rates multiplied with the 3D encountered dust mass. The dust destruction results, averaged over the three investigated orthogonal planes, are compared in Table 1 and the results of the z = 0 plane dust simulations are visualized in Fig. 4 where the rows show the different Mach numbers and the columns show the average ISM densities (n = 1, 10, 100 cm⁻³).

3.1 Low-density ISM

In a low-density environment (n = 1 cm⁻³, left column of Fig. 3), the forward shock is not efficiently slowed down by ISM gas (247–388 km s⁻¹ at 10 kyr, Table 1) and can thus expand to large sizes, reaching up to 24–40 pc, depending on the studied turbulent Mach number. Due to the low amount of shocked gas in this ISM and the fixed initial gas-to-dust mass ratio, small amounts of dust are encountered, resulting in low destroyed dust masses of 1.4–1.5 M_⊙ for silicate dust, and 0.85–0.97 M_⊙ for carbonaceous dust after 10 kyr in our simulations. The destruction efficiencies are also rather low with 43–57% for silicate and 27–36% for carbonaceous dust (also shown in Table 1 and Fig. 4), which may be unexpected given that the shock front is fast and hot throughout the whole simulation, implying efficient thermal and kinetic sputtering. However, sputtering of dust grains can only occur when they collide with the hot and fast ions of the shock front. Since these ions have a rather low number density and move through a low-density ISM, which has an even lower dust grain number density of 0.01 cm⁻³ on average, the energetic ions of the forward shock do not collide often with dust grains. Thus, they cannot transfer their energy to the ISM grains (i.e., sputter them) efficiently. Due to the low grain number density, grain–grain collisions rarely occur, leading to differences of maximum 2.6% compared to the sputtering-only simulations. Therefore, the processed grain-size distributions with and without grain–grain collisions are similar for grains larger than 5 nm, showing a survival of most grains above ~40 nm and significant destruction fractions for smaller grains. Many grains smaller than 5 nm are generated, increasing even more with included grain–grain collisions that produce many small fragments (see Appendix E).

The filaments in the high Mach, n = 1 cm⁻³ simulations are much denser than the average ISM, but their density is mostly below 100 cm⁻³, which is less than the average density in our high ISM density cases. Thus, they are not dense enough to stop the forward shock but get penetrated by it. This can be seen in Fig. 5, which shows the interaction of the forward shock with a filament in the n = 1 cm⁻³ and ℳ = 1 case. The filament does not preserve its integrity but is processed by the forward shock. Fingers form from the destroyed filament into the ejecta, not only at the shown filament, but also at other forward shock positions (see Fig. 3).

A few high-density cloudlet structures can be seen that are a few pixels wide and have higher densities than the initial filaments (n > 100 cm⁻³, e.g., the ones circled in Fig. 5). These cloudlets result from compressed material of the processed filament and can shield some dust from destruction as they take longer to be disrupted than the lower density surroundings.

Although the filaments do not withstand the energetic forward shock but get processed to lower density fingers and eventually destroyed, they are able to shield some dust in comparison to the homogeneous case in the most destructive, first 10 kyr because it takes longer for them to be penetrated and processed by the shock than in the homogeneous case. Roughly 14% of silicates and 8% of carbonaceous dust grains are shielded by the filaments in the supersonic turbulent Mach number case, decreasing with turbulence strength to 4% of silicate and 3% of carbonaceous dust for a subsonic Mach number as can be seen in Fig. 6. This supports the conclusions of previous papers with a longer evolution time that find that filaments or inhomogeneities in general are able to shield dust from destruction (Hu et al. 2019; Kirchschlager et al. 2024a; Dedikov & Vasiliev 2025). However, the destroyed dust mass remains on the same order of magnitude for all Mach number simulations since the forward shock encounters more dust in the high Mach number cases.

This happens due to the turbulent ISM structures in the high Mach number cases that lead to an irregular expansion of the SNR. While it can rapidly move through the voids that have a very low density, the forward shock is slowed down more in higher density regions. Since filaments are spatially thin in comparison to the voids, the shock can evolve through lower density areas around them (cf. Korolev et al. 2015; Slavin et al. 2017; Dedikov & Vasiliev 2025). Thus, after 10 kyr, the forward shock undergoes greater expansion in the high Mach number ISM than in the homogeneous case and is also more efficiently slowed down on average across the entire shock, due to its low velocity in the filaments.

The velocity of the forward shock, which can be estimated by averaging the gas velocity of the shock cells traced by a momentum threshold, is therefore still 388 km s⁻¹ in the homogeneous case, whereas it slows down to 251 km s⁻¹ in the supersonic turbulence case after 10 kyr, as given in Table 1. Since kinetic sputtering is directly dependent on the forward shock velocity, a comparison of the simulations at a similar end shock velocity instead of similar end time would also be a valuable comparison. This was done in Appendix F, where all simulations are compared at different evolution times but similar average shock velocities of 387–390 km s⁻¹. For the low-density simulations, the spread in dust destruction efficiency between the different Mach numbers increases to 26% for the silicate grains and 16% for the carbonaceous grains; however, the high Mach number simulations show significantly less destroyed dust mass in this comparison, reflecting the trend that the filaments can shield the dust from destruction.

Regardless of the comparison criterion of the simulations (Table 1 or F.1) and the studied parameter set, silicates are always at least 12% more efficiently destroyed than carbonaceous grains in a low-density environment. This difference can generally be explained by the higher sputtering yields of silicates that stem from their higher bulk density². Thus, the silicate grains are also less efficiently accelerated by the forward shock, resulting in higher relative dust-to-gas velocities and thus more efficient kinetic sputtering at the shock front. However, other mechanisms are also affected. Due to their higher bulk density, silicates are less numerous than carbonaceous grains since we consider the same dust sizes and total dust mass for both species, resulting in fewer grain–grain collisions and a lower sputtering efficiency due to a lower total grain surface area. Furthermore, silicates have lower relative grain-to-grain velocities when they are accelerated by the forward shock so that fragmentation and vaporization in general are less efficient for silicates than for carbonaceous grains. Moreover, the size-dependency of sputtering is considered in Paperboats (see Kirchschlager et al. 2019), which arises since smaller grain sizes are used in comparison to the experimental sputtering yields which considered a semi-infinite target (Jurac et al. 1998; Serra Díaz-Cano & Jones 2008). This correction was updated in our study to consider hydrogen ions of the ISM to impinge on dust grains as projectiles instead of oxygen ions as used previously by Kirchschlager et al. (2022, 2024a). It generally leads to fewer silicate grains but more small carbonaceous grains being destroyed, supporting a better survival of silicates. For the low-density simulations, however, the generally higher sputtering yields and more efficient kinetic sputtering during the acceleration of the grains seem to be more important than the higher sputtering efficiency due to a higher total grain surface area, the corrections of the size-dependent sputtering, and the higher fragmentation and vaporization efficiencies.

In total, more than 1.3M_⊙ of silicate dust is destroyed in each low-density simulation in the first 10 kyr, which is higher than the maximum amount of dust (~1 M_⊙) produced in SNRs as inferred from observations. However, this threshold is not reached in the case of carbonaceous grains yet: slightly less than 1 M_⊙ is destroyed in the first 10 kyr. Although this seems to imply that SNRs are net dust destroyers for an ISM with predominantly silicate grains but can be net producers for a carbonaceous dust-dominated environment, the results of our simulations are only lower limits of the dust destruction, since we only investigate the first 10 kyr of SNR evolution. At 10 kyr, the dust destruction efficiency of all low-density simulations has not yet plateaued, meaning that the dust is still efficiently destroyed. Furthermore, the forward shock is still fast and hot at 10 kyr (>240 km s⁻¹) so that it is reasonable to expect efficient dust destruction by kinetic and thermal sputtering beyond 10 kyr. If the so far neglected reverse shock dust destruction would also be taken into account, it is probable that SNRs destroy more dust than they produce in this low-density environment, even if the ISM is dominated by carbonaceous dust.

We also note that, although we see locally reflected shocks that are generated at filaments and reverse through the SNR, a global reverse shock does not occur in any of our simulations. Since we do not include a mass injection with the SN explosion, which would incur a proper ejecta evolution, there is no ejecta mass that could compress the gas behind the forward shock. Thus, there is no pressure build-up leading to a global reverse shock, which could destroy ejecta dust. However, studying the reverse shock destruction of ejecta dust is beyond the scope of this paper.

Fig. 4 Dust destruction efficiencies (left y-axis) and masses (right y-axis) in the z = 0 slice. This destruction efficiency refers to the ratio of the destroyed dust mass at each time step to the total encountered dust mass at the last considered snapshot. Different Mach numbers of ℳ = 0, 0.3, 1, and 3, and average ISM densities of n = 1, 10, and 100 cm−3 are shown. For each Mach number and ISM density parameter set, four different dust simulations were carried out, including either silicate (Sil) or carbonaceous (Car) dust, and either grain–grain collisions and sputtering (gg+sp) or only sputtering (sp).

Fig. 5 16 × 16 pc zoomed-in image of a filament of the n = 1 cm−3 and ℳ = 1 simulation at t = 0 yr and t = 10 kyr. The filament cannot withstand the shock efficiently, but is destroyed leaving behind finger-like structures that arise from the forward shock front into the SNR interior. Furthermore, some high-density cloudlets are formed (circled in black). The arrow indicates the approximate direction of the forward shock.

Fig. 6 Dust destruction efficiencies in the z = 0 slice in a direct comparison between the different Mach numbers (ℳ) used for each ISM density (n) case, considering either carbonaceous (Car) or silicate (Sil) grains.

3.2 Intermediate-density ISM

The SNR forward shock in the intermediate-density environment (n = 10 cm⁻³, central column of Fig. 3) is slowed down more efficiently than in the low-density environment, ending at 198–226 km s⁻¹ (Table 1) and only reaching sizes of 16–28 pc. In total, 58–92% or 4.4–5.6 M_⊙ of silicate dust, and 46–77% or 3.5–4.7 M_⊙ of carbonaceous dust is destroyed. Thus, carbonaceous grains are again more robust than silicate grains with similar differences of at least 12% more destruction of silicates as in the low-density simulations.

With more than 3.5 M_⊙ of destroyed dust for every intermediate-density simulation, significantly more dust is destroyed than in the low-density cases. The destruction efficiency is also higher since collisions of the hot and high-velocity ions with the ISM dust grains occur more often, resulting in more efficient thermal and kinetic sputtering. Grain-grain collisions are still not very efficient in the homogeneous and subsonic Mach number cases, leading to differences of ~1%. They make a bigger impact of 2.9–12% for the more turbulent cases of this ISM density as they are typically more efficient in higher density media such as filaments. For all cases, similar to the low-density environment, the grain-size distribution with and without grain–grain collisions is similar for grains larger than 5 nm but is altered for smaller grains (see Appendix E). In general, the final grain-size distribution is altered in a similar way as in the low-density cases, having more smaller grains destroyed than larger grains. However, as grain–grain collisions become effective in the high Mach number simulations, more and more large grains are destroyed.

Compared to the low-density ISM, the filaments in the intermediate-density transonic and supersonic cases have a higher number density, mostly on the order of 1000 cm⁻³. Thus, they can better withstand the forward shock that itself is generally more decelerated due to the higher average ISM density. When the filaments are encountered by the forward shock, they are not immediately destroyed, but partly swept up and dragged along with the shock, so that a thin and dense shell-like structure can form, which can cool down more efficiently (see Fig. 3). However, most of the filaments are still eventually destroyed, leaving behind finger-like structures that expand into the ejecta, similar to Fig. 5, and high-density cloudlets.

Because of the longer filament survival times and the lower penetration velocity, the dust is shielded much better in the dense regions of the higher Mach number cases, independent of the dust species, as shown in Fig. 6. Thus, we can observe a large difference in dust destruction efficiencies between homogeneous and strongly turbulent ISM of up to 34% for the silicate grain simulations and 31% for the carbonaceous grain simulations, decreasing with Mach number to 0.9% dust destruction difference with silicates and 0.6% difference with carbonaceous grains. As a result, also less dust mass is destroyed by the SNR in the higher Mach number simulations, although it expands to larger sizes.

At similar end shock velocities of 387–390 km s⁻¹ instead of similar end times, this trend remains visible and accounts for up to 39% and 34% of dust destruction difference for the silicate and carbonaceous grains, respectively (Appendix F). The amount of destroyed dust mass also has a larger spread between the Mach number cases, and remains higher than in the low-density ISM; however, the dust destruction difference between intermediate-and low-density ISM cases decreases slightly.

Although filaments are able to shield the dust from destruction, the SNR forward shock is strong enough to always destroy much more dust than 1 M_⊙ in the first 10 kyr. Furthermore, the destruction efficiencies do not plateau yet, indicating that the dust is still rather efficiently destroyed at the end of our simulations. Thus, SNRs in this environment are net dust destroyers, even if the reverse shock would not destroy any produced dust and only the first 10 kyr of evolution would be considered.

3.3 High-density ISM

With the highest studied average ISM density (100 cm⁻³), we aim to replicate a molecular cloud environment with filament densities of the order of 10⁴ cm⁻³. The SNR forward shock is strongly slowed down to 88–128 km s⁻¹ (Table 1), thereby only reaching a maximum radius of 10–17 pc in size after 10 kyr (right column of Fig. 3). After 10 kyr, 76–87% or 9.1–11.0 M_⊙ of silicate dust, and 73–87% or 9.1–10.6 M_⊙ of carbonaceous dust is destroyed (Fig. 4).

In contrast to the intermediate- and low-density cases, the forward shock in the high-density case is compressed to a thin and dense shell from ∼3–4 kyr onward. This happens because larger amounts of gas are swept up than in the lower density cases, the shock is decelerated more to v < 130 km s⁻¹, and the SNR interior cools down enough. The formation of the shell can be seen well in the high-density, homogeneous ISM simulation in Fig. 3, in which the SNR interior is empty whereas the shell has a very high density and starts to exhibit Vishniac thin-shell modes (Vishniac 1983, 1994), which are inhibited by the irregular expansion of the forward shock in the higher Mach turbulence cases.

As soon as the dense shell forms, grain–grain collisions become very efficient. They account for a maximum difference of 36% in dust destruction at 10 kyr, which is significantly more than in our low- and intermediate-density models. Due to the high density of the shell, the forward shock starts to cool down more efficiently, and thus slows down more so that the effectiveness of both thermal and kinetic sputtering decreases. At a certain critical velocity, which is dependent on the grain composition, the kinetic sputtering yield drops steeply, making it even less efficient (Tielens et al. 1994; Nozawa et al. 2006). Hence, after 3–4 kyr, the sputtering-only simulations start to have a shallow slope in dust destruction (Fig. 4). The simulations with all destruction mechanisms included continue to have high destruction efficiencies. Thus, with the onset of the shell, the significance of the destruction mechanisms alters since sputtering becomes inefficient whereas grain shattering gains in effectiveness.

Due to this change, the destruction efficiencies of carbonaceous and silicate grains in the high-density ISM are only 3% different at most. For the weakly turbulent and homogeneous cases, carbonaceous grains are even slightly more efficiently destroyed than silicates. This partly results from the higher shattering efficiency of carbonaceous grains in comparison to silicate grains, as discussed in Sect. 3.1, and becomes evident from the larger differences of at least 6% between the carbonaceous and silicate grain simulations with only sputtering considered, shown in Table 1. However, as described in Sect. 3.1, several mechanisms play a role in destroying more carbonaceous than silicate grains in these particular simulations.

The grain-size distribution also changes more significantly in these high-density cases because larger dust grains are more efficiently destroyed by shattering than smaller grains. The processed grain-size distribution is therefore shifted to smaller grain sizes compared to the initial MRN distribution, as discussed in Appendix E.

Since sputtering becomes significantly less efficient with the formation of the shell, the trend of a higher dust destruction efficiency with higher ISM density also discontinues with the n = 100 cm⁻³ cases although more dust mass is destroyed because more dust is encountered. While the forward shock dust destruction efficiency in a high-density and low Mach number ISM is still higher than in the intermediate-density medium at 3 kyr, it is roughly 4–5% lower than in the intermediate-density cases at 10 kyr. If only sputtering is considered, almost all other Mach number cases also decrease in destruction efficiency in the high-density ISM compared to the intermediate-density cases after 10 kyr, by up to 21%.

However, a comparison at the same end shock velocities as in Appendix F still shows the trend of more dust destruction with a higher average ISM density because the shell has not yet formed during the modeled period. Thus, a higher dust fraction and more dust mass is destroyed in the high-density ISM than in the intermediate-density cases. We note that this comparison has to be judged carefully since the SNR shock does not encounter many filamentary structures in the higher Mach number simulations yet, and thus, mostly evolves in a void with a density n ≲10 cm⁻³, rather resembling an intermediate- or low-density case. Thus, the dust destruction in the higher Mach number simulations is also much lower than in the homogeneous ISM with large differences of up to 36%.

Smaller differences are found after 10 kyr of evolution time. The destruction efficiency from the homogeneous to the strongly turbulent ISM case is only reduced by 11% for the silicate grains and 14% for the carbonaceous grains, as seen in Fig. 6. This is less than in the intermediate-density cases and similar to the low-density cases although the filaments are able to slow down the shock more efficiently and almost stop it at high densities. The forward shock shape at the end of our simulations is similar to the initial filament structures, as shown in Fig. 7, and almost no finger-like shapes are formed but some dust-shielding cloudlets occur. However, due to the shell formation in the high-density cases and the subsequent shift of significant destruction mechanisms, the shielding ability of the filaments decreases. Since shattering is more efficient in higher density media, it shrinks the difference in dust destruction between the homogeneous and strongly turbulent cases by destroying dust in dense filaments more efficiently than in the homogeneous and weakly turbulent runs. Thus, at roughly 4 kyr, the dust destruction efficiency difference between weak and strong turbulence is still on the same order of the intermediate-density cases with over 30%, but decreases again after the shell forms.

In total, SNRs always destroy more than 9 M_⊙ of dust in their first 10 kyr of evolution in a high-density ISM, which is a magnitude more than what SNe have been observed to produce at most, making them significant net dust destroyers in this environment. An even higher dust destruction would be found if the considered ejecta cooling would be less efficient than our rather strong oxygen cooling so that the SNR would expand more and the shell would form later on, resulting in more dust being encountered and thus destroyed and a longer period of efficient sputtering.

Fig. 7 Zoomed-in 9 × 9 pc image of a filament of the n = 100 cm−3 and ℳ = 3 simulation at t = 0 yr and t = 10 kyr. The dense filament almost stops the shock and dictates where it will evolve. The initial filament shape looks nearly identical to the forward shock shape at the end. The arrow indicates the approximate direction of the forward shock.

3.4 Summary

To summarize, the SNR forward shock dust destruction depends on the properties of the surrounding medium (ISM density and turbulent Mach number), the dust grain species, and the considered destruction mechanisms. If only one of these parameters is varied, the dust destruction efficiency can change significantly. Generally, we find that the dust destruction efficiency decreases with lower ISM densities and higher turbulent Mach numbers. It is further reduced if carbonaceous instead of silicate dust is considered. Moreover, grain–grain collisions always generate many grains with sizes lower than 5 nm and can have a large impact on the destruction rate in intermediate- and high-density environments (up to 36%).

The lowest forward shock dust destruction in the first 10 kyr is found in the low-density (n = 1 cm⁻³) and strongly turbulent (ℳ = 3) environment, where only 27% of carbonaceous dust is destroyed. With 0.85 M_⊙ of carbonaceous dust, the least dust mass is destroyed in the n = 1 cm⁻³ and ℳ = 0.3 case. However, the forward shock still has a high velocity of v > 250 km s⁻¹ and is still hot at the end of these simulations, so that especially kinetic and thermal sputtering will process much more dust beyond the considered period. Compared to the average dust mass of 0.1–1 M_⊙ that is observed to be produced by SNRs and that still needs to survive the high-velocity reverse shock, we conclude that SNRs are rather net dust destroyers, in particular if higher gas densities are considered.

If grain growth in the ISM were to be efficient, the total grain surface area would become more important than the destroyed dust mass. The many small grains that are produced by grain– grain collisions in our simulations could then even dominate the surface area of all grains and thus enable more efficient grain growth in the ISM, as discussed in Appendix G.

4 Discussion

Our study provides a comprehensive overview of the dust destruction efficiency by the forward shock during the first 10 kyr of SNR evolution. It is the first to examine dust destruction in a turbulent medium that is specifically modeled to replicate ISM observations and has one of the highest resolutions of forward shock dust destruction studies because it focuses on the critical early stages of SNR evolution. A comparison to other studies is very valuable for this paper as they consider longer evolution times – often even up to 1 Myr.

4.1 Previous studies

By comparing our results to previous studies with longer evolution times, we can roughly estimate the amount of dust that would be destroyed by the forward shock if our simulations were longer, since we already concluded that most likely more dust will be destroyed beyond the considered first 10 kyr of SNR evolution. Additionally, this comparison lets us assess how well previous models predicted their first few thousand years of dust destruction because they did not focus on this most destructive period. Usually, a lower spatial (and often temporal) resolution was used and some dust destruction mechanisms were often neglected. Table 2 compares the main parameters and results of our own and other works after 10 kyr of SNR evolution.

Dedikov & Vasiliev (2025) conduct such a study, using a similar approach to ours by simulating the forward shock dust destruction over 100 kyr in a clumpy instead of turbulent medium, which is generated by fractal cubes. They use different inhomogeneity strengths and their densest clumps have roughly the same density contrasts between high- and low-density regions as our supersonic Mach number simulations. They do not consider magnetic fields nor grain–grain collisions, only kinetic and thermal sputtering. In their paper, they conclude that dense clumps can shield the dust efficiently, resulting in up to ∼20% less destruction over 100 kyr for an ISM with n = 1 cm⁻³ on average. After 10 kyr, they find 1.17 M_⊙ (∼39%) of destroyed dust when considering a homogeneous medium, decreasing with increasing inhomogeneity down to 0.67 M_⊙ (∼32%) of destroyed dust with the densest clumps so that these clumps already shield ∼7% of the dust at 10 kyr. Their considered grain species is not mentioned but the results match reasonably well to our carbonaceous (0.79–0.91 M_⊙ or 25–33%) grain simulations with n = 1 cm⁻³ that include only sputtering, for which we also find significant filament shielding of 8% at 10 kyr. In an intermediate-density environment (n = 10 cm⁻³) with rather dense ISM clumps, ∼4.5 M_⊙ (∼50%) of dust is destroyed in the first 10 kyr, also fitting reasonably well to our results of the intermediate-density case with strong turbulence and carbonaceous grains (3.5–4.2 M_⊙ or 46–63% of dust destroyed). This supports the conclusion that generally more dust is destroyed at higher ISM densities. At 100 kyr, however, this trend can reverse according to their study where 8.0–11.4 M_⊙ dust is destroyed in the low-density case, but only 7.5 M_⊙ in the inhomogeneous intermediate-density case.

If our sputtering results for carbonaceous grains would also align quite well with the results of Dedikov & Vasiliev (2025) beyond 10 kyr, over 8 M_⊙ of dust would be destroyed at 100 kyr, which is almost an order of magnitude more than what SNRs are observed to produce. Nevertheless, we note that the filaments in our study could behave differently than the clumpy medium of Dedikov & Vasiliev (2025) at lower shock velocities, and that we cannot estimate the dust evolution in our own simulations that include grain–grain collisions beyond 10 kyr based on this comparison, as this destruction mechanism could become more efficient at later SNR stages (see for example our n = 100 cm⁻³ case).

Kirchschlager et al. (2024a) also investigate single SNRs in a turbulent ISM, similarly to us. However, they use SN-driven turbulence, which brings the caveats mentioned in Sect. 1 with it. After turbulence driving, they inject an SN in either a low- or intermediate-density area of n = 0.03 or 0.7 cm⁻³, respectively. Thus, although their average ISM density is always lower than our lowest density case (n = 1 cm⁻³), the SN injection in their moderately dense region might be comparable to our simulations because their SNR does not encounter many voids, whereas we set off the SN in a void. Their dust properties and considered dust destruction mechanisms are similar to our study, but they only consider silicate grains.

After 10 kyr, they find a destroyed dust mass of 2.8 M_⊙ in the intermediate-density simulation, which is almost two times more than the destroyed silicate grain mass in our high Mach turbulence study (1.5 M_⊙). This disparity can come from multiple factors. Most importantly, the turbulence driving is different such that the filaments they produce appear much wider than ours. Furthermore, the average ISM density in the box is lower with the SN set off in a higher density environment. Thus, the dust destruction efficiency can change significantly (Kirchschlager et al. 2024a) and the encountered dust mass can be different. Moreover, their simulation is not fine-tuned to the first 10 kyr, meaning that the temporal and spatial resolution is rather coarse (250 yr and 0.5 pc) compared to our study (logarithmic time steps between 0.15–120 yr depending on the evolution time and 0.044 pc in the n = 1 cm⁻³ cases), which can cause inaccuracies in the gas drag. With a lower spatial resolution and thus bigger grid cells, the dust velocity in a single cell is averaged over more grains, resulting in larger deviations compared to the real grain velocities. If the temporal resolution is too low, grains that are dragged with the forward shock may not be as efficiently decelerated by encountered ISM filaments as the gas so that the grains could be able to overtake the shock. This would leave the high-velocity grains exposed in a rather stationary ISM, leading to high relative gas-grain and grain–grain velocities.

Thus, the higher destruction efficiency of Kirchschlager et al. (2024a) could be explained by an insufficient resolution. Other differences can add to these inaccuracies. First of all, the extrapolation from the 2D to the 3D dust mass is different and induces uncertainties. We multiply the averaged destruction efficiency of the three investigated 2D planes by the 3D SNR-encountered dust mass, whereas Kirchschlager et al. (2024a) fully extrapolate the dust mass from one 2D plane, using a correction factor they derive from comparing the SNRs in a homogeneous and an inhomogeneous medium. Thus, the real 3D destruction could be different. Secondly, they do not include partial vaporization, which was implemented in Paperboats by Kirchschlager et al. (2023). This slightly reduces the dust destruction for the largest grains in our simulations. Thirdly, the penetration depth of gas ions into dust grains was calculated considering oxygen ions in their study, but was updated in our study to consider hydrogen instead of oxygen ions, which have a different penetration depth. This value alters the size-dependent sputtering, so that significantly more small silicate grains are destroyed. Lastly, it should be mentioned that the dust destruction in a turbulent ISM also depends on the detailed ISM environment, as we show with the three different orthogonal planes in this paper. The differences in the destruction efficiencies could therefore also partly come from the different surrounding environment and positions of the filaments.

Similar to our study, however, the derived dust destruction efficiencies of Kirchschlager et al. (2024a) do not plateau after 10 kyr, but only start to plateau after 50 kyr in their study. Nevertheless, they find that some dust continues to be destroyed even when their simulations finish at 1 Myr, with a total of 37 M_⊙ of destroyed dust, which is ∼13.5 times more than what they find at 10 kyr. If this trend is applied to our high Mach number and low-density ISM silicate simulation, it would mean that even in one of the lowest dust destruction cases, the SNR would destroy 1.5 M_⊙ ×13.5 ≈ 20 M_⊙ of dust after 1 Myr, being 20 times higher than the observed dust production in SNRs. This implies, once again, that SNRs are rather net dust destroyers.

The study of Kirchschlager et al. (2022) is similar to Kirchschlager et al. (2024a), but only for a homogeneous medium and gas that is perturbed by weak Vishniac–Ostriker–Bertschinger overstabilities, which does not aim to resemble ISM turbulence. This medium can be best compared to our homogeneous and subsonic cases. They use the same setup and resolution as Kirchschlager et al. (2024a), but consider an n = 1 cm⁻³ case, which is more comparable to our study. They find a destroyed dust mass of 7.0 M_⊙ in the first 10 kyr which is reduced to 1.3 M_⊙ if only sputtering is considered, consistent across their inhomogeneous and homogeneous cases, which show almost identical destroyed dust masses. Their resulting dust destruction in the first SNR stages is, similarly to Kirchschlager et al. (2024a), much higher than our results for this low-density environment, where we find ∼5 times less destroyed dust mass in the homogeneous and subsonic cases (1.4–1.5 M_⊙) but roughly the same destroyed dust masses with only sputtering considered (1.3–1.4 M_⊙). In addition to the differences mentioned in the comparison between our study and that of Kirchschlager et al. (2024a), the study of Kirchschlager et al. (2022) also does not include a magnetic field, resulting in higher dust destruction due to weaker gas-grain coupling, as Kirchschlager et al. (2024a) discuss. Given these disparities between our own study and that of Kirchschlager et al. (2022), the similarity in the sputtering-only simulations is most likely coincidental. Their simulations finish at 1 Myr with ~70 M_⊙ of destroyed dust, which is ~10 times more than at 10 kyr. If the SNRs in our study would follow the same trend, they would therefore again always destroy much more dust after 1 Myr than they had originally produced.

Another useful comparison can be made to works that consider pre-SN feedback. Martínez-González et al. (2019) let a wind-driven bubble evolve in a homogeneous ISM with a density of 1 cm⁻³ before they inject the SN. The bubble’s shocked region is compressed and has a density that is roughly three to four orders of magnitude higher than the interior of the bubble in which the SN explodes, resulting in roughly the same difference between filaments and voids as in our low-density and supersonic turbulence simulations. It expands up to 23 pc in radius, meaning that it sweeps up ∼12.6 M_⊙ of dust during the stellar wind evolution. They find that ∼0.36 M_⊙ of dust (2.9% of the shell’s dust mass) is destroyed in the first 10 kyr after SN explosion if only thermal sputtering is considered and a mixture of silicate and carbonaceous grains is used. Although we also find that dust can be shielded by filaments, significantly more dust is shielded by their stellar wind shell than in our sputtering-only results with high Mach turbulence for either grain species (0.88–1.5 M_⊙ and 25–45%). However, their bubble is larger and its shell is roughly 3–5 pc, which is much thicker than our filaments, meaning that it can slow down the shock more efficiently. As a result, most of the dust in the shell has not yet been encountered after 10 kyr so that this dust cannot be destroyed. Furthermore, they do not include kinetic sputtering, which could be efficient and increase the destroyed dust fraction when the SNR forward shock starts to encounter the wind-driven bubble. At this point, grain–grain collisions could also become efficient due to the high dust densities and velocities, which would further increase the dust destruction.

Other studies that investigate the dust destruction of SNR forward shocks in a slightly inhomogeneous or turbulent ISM are Vasiliev & Shchekinov (2024) and Hu et al. (2019), respectively. The former uses the same code as Dedikov & Vasiliev (2025) but only introduces small density perturbations of 10% and does not compare their results to a homogeneous case so that the filament shielding efficiency cannot be assessed. Since they only study a monosized dust population, we cannot easily compare their study with ours because we use an MRN-distributed dust population, which leads to significantly different dust behavior. Hu et al. (2019) drive their turbulence with SNe, similar to Kirchschlager et al. (2024a). However, they study the lifetime of silicate and carbonaceous dust grains over a total of 400 Myr and several SNe. Thus, their destruction efficiencies are also not easily comparable to our study that only covers the first 10 kyr. However, similarly to us, they find that silicates are better destroyed than carbonaceous grains.

4.2 Caveats

Although we include many mechanisms that are often neglected by other dust destruction studies, some uncertainties remain:

• The pre-SN effects of stellar winds and radiation were only mimicked simplistically in our study by setting off the SNe in voids because a full pre-SN feedback model is out of the scope of this paper. In particular the weakly turbulent and homogeneous runs did not incorporate realistic pre-SN structures. For a forthcoming model, we plan to improve this by allowing stellar winds and radiation to evolve naturally in separate simulations, and letting the SN explode inside of the wind-driven bubble;

• Although we consistently set off all SNe in voids, the detailed surrounding ISM structures are not the same, since the turbulence-driving results in different filamentary shapes for each simulation. This affects the 3D SNR evolution in the first 10 kyr and therefore also the dust destruction, resulting in an uncertainty, which could be quantified best by studying many different 3D ISM environments for every parameter set. Since the needed resources for this are too high, we assess this uncertainty with our investigation of the three orthogonal 2D planes, as discussed in Appendix D. We find a maximum uncertainty of 9.4% and include the standard deviation between those planes in our dust destruction results of Tables 1 and F.1. This implies that the SNR dust destruction cannot be precisely predicted based on the initial parameters like ISM density and turbulent Mach number, but needs to take the exact SN ejection location (see also Appendix B) and the surrounding environment structures into account;

• Our post-processing approach does not allow the dust to feed heavy elements coming from destroyed dust grains back to the gas, which physically should happen and would induce for example a faster cooling of the gas. The exact impact of the dust feedback is so far unknown but might be significant. Since on-the-fly approaches are more expensive than post-processing simulations, we used the latter for our calculations to include a more extensive dust destruction model. However, this brings the drawbacks of no dust-to-gas feedback and coarser time steps, which gives inaccuracies in the destruction, transport, and charge processes that we tried to minimize by taking as many logarithmically spaced snapshots as computationally feasible. A comprehensive on-the-fly approach is currently being implemented, tested, and compared against our post-processing method by our group (Sartorio et al., in prep.). Various on-the-fly models were used by other groups before (e.g., Slavin et al. 2015; Martínez-González et al. 2019; Vasiliev & Shchekinov 2024; Dedikov & Vasiliev 2025); however, their dust destruction models were not as complete as in our post-processing code, so we chose Paperboats;

• We did not model the ejecta-dominated phase of the SNR evolution. Instead, we assumed that the ejecta mass can always be neglected in comparison to the encountered ISM mass and we inserted only thermal energy and no ejecta mass. Thus, the SNR evolution started immediately in the Sedov-Taylor phase instead of the ejecta-dominated phase. Furthermore, the pressure support of the ejecta was not modeled in this stage. The maximum amount of dust that would be encountered during the ejecta-dominated phase is 2–6% of the encountered mass after 10 kyr in our low ISM density simulations, depending on the considered ejecta mass (5–15 M_⊙). For our high- and intermediate-density ISM simulations, this fraction is ∼0.5–2.5% and thus not significant. At low densities, the ejecta mass would be significant only if the ejecta-dominated phase were skipped entirely. Instead, we replaced it by an immediate Sedov-Taylor phase, which has a similar evolution to the late ejecta-dominated phase. Thus, only a fraction of the 2–6% of the total dust mass is uncertain due to the missing ejecta mass and the dust destruction in our study can be expected to be only slightly underestimated in our low ISM density cases. However, since we started all of our simulations with an immediate Sedov-Taylor phase, concluded trends would not change qualitatively by considering ejecta mass;

• Because we did not include ejecta mass, we did not aim to resemble any specific SNR or SN type. Instead, we modeled SN type IIP and IIL in general by using a thermal energy of 10⁵¹ erg and having no pre-SN mass loss. Hence, our model only applies to this general approach and no observed SN;

• Chemistry was not fully traced. With an assumed collisional ionization equilibrium and an electron-ion equilibration, cooling times, dust charges, and dust transport could be influenced. Non-equilibrium ionization could for example decrease the temperatures at the shock front (Slavin et al. 2015), resulting in less efficient thermal sputtering. However, tracing ionization states in detail, and ion and electron temperatures separately, is beyond the scope of this study;

• Although Arepo uses time step subcycling to handle unresolved cooling times t_cool (Glover & Mac Low 2007), the cooling length was not fully resolved during the whole simulation. The cooling length, l_cool = c_st_cool with the speed of sound c_s, can be as short as 3 × 10⁻⁷ pc in the shell of our high ISM density simulations because the ejecta cooling function peaks at ∼2 × 10⁵ K. This results in undercooled gas at the shock front, overestimating gas temperatures and therefore thermal sputtering. This undercooling seems to be a general problem of SNR dust destruction studies due to the large difference of the cooling length and the used numerical resolution. Thus, the significance of the unresolved cooling length cannot be easily assessed.

We emphasize that, while several assumptions may influence the absolute values of dust destruction in different ways, our model was applied consistently across all simulations. Therefore, the relative trends identified between simulations remain robust and meaningful.

5 Conclusion

We performed 3D MHD Arepo simulations of SNRs evolving in a turbulent ISM that was driven to resemble observations more closely than in previous studies. Different turbulent ISM environments were considered with average ISM densities of n = 1, 10, and 100 cm⁻³ and turbulent Mach numbers of ℳ =0, 0.3, 1, and 3. We let SNRs evolve for the most destructive first 10 kyr and post-processed these simulations with Paperboats to study the dust destruction by the SNR forward shock, considering either purely carbonaceous or purely silicate dust. In detail, we find the following:

1. The dust destruction by SNR forward shocks in their first 10 kyr of evolution varies significantly, depending on the ISM environment. For n = 1–100 cm⁻³ and ℳ = 0–3, the destruction ranges from 27–92%, which refers to 0.85–11 M_⊙ of dust. The largest dust mass is destroyed in the high-density environment. The least dust is destroyed in the low-density ISM where the shock still has end velocities of 250–400 km s⁻¹. Since the total destroyed dust mass has not plateaued yet, it can be expected that more dust will be destroyed in the following years of SNR evolution, as is also derived from comparisons to other studies. Thus, we find that SNRs destroy more ISM dust than they produce, and thus are net dust destroyers;

2. Filaments that closely resemble observations slow down the forward shock, and therefore shield some dust from destruction, similarly to the conclusions of other studies with different inhomogeneities or turbulence driving methods. While filaments in the low-density ISM only shield 8–14% of the dust in the first 10 kyr of SNR evolution, the intermediate-density cases result in 31–34% less destruction in the highest Mach number simulations compared to the homogeneous ISM. In the high-density medium, the forward shock evolves into a thin and dense shell so that grain–grain collisions become significantly more important, resulting in a lower shielding ability of 11–14%;

3. Silicates are almost always more efficiently destroyed than carbonaceous grains, which is primarily caused by the higher sputtering yields. In a high-density ISM, the destruction difference between these two species decreases, partly caused by more effective grain–grain collisions, which are more efficient for carbonaceous grains than for silicates. In the weakly turbulent and homogeneous high-density cases, silicates are even more robust than carbonaceous grains;

4. The dust destruction rate can vary even for similar turbulent environments with the same statistical parameters (average ISM density and turbulent Mach number) because the detailed structure of the turbulence changes, which affects the evolution time and shock velocity at which the first filaments are encountered. We can estimate this difference from the investigated three orthogonal planes, which show a difference in dust destruction rate of up to 9.4% in the first 10 kyr.

Altogether, filaments in the turbulent ISM cannot shield dust enough for SNRs to act as net dust producers. Future work will incorporate pre-SN feedback and dust cooling into the workflow, which may further reduce dust destruction.

Data availability

The movies associated to Figs. A.1–A.3 are available at https://www.aanda.org

Movies

Movie 1 associated with Fig. A.1 Access here

Movie 2 associated with Fig. A.2 Access here

Movie 3 associated with Fig. A.3 Access here

Acknowledgements

TS, FK, and IDL were supported by funding of the Euro-pean Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program DustOrigin (ERC-2019-StG-851622). TS and IDL acknowledge funding from the Flemish Fund for Scientific Research (FWO-Vlaanderen) through the research project G0A1523N. TS also acknowledges funding from FWO-Vlaanderen through the doctoral research project 11A5O26N. NSS acknowledges the support from the FWO-Vlaanderen in the form of a postdoctoral fellowship (1290123N). The resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the FWO-Vlaanderen and the Flemish Government.

References

Appendix A SNR evolution

Figures A.1, A.2, and A.3 show the temporal SNR evolution of all simulations in this study for the low-, intermediate-, and high-density cases, respectively. Their rows refer to the homogeneous, subsonic, transonic, and supersonic simulations and their columns to three different time steps. Links to movies of the SNR evolution are given in the captions of the respective evolution figures.

Fig. A.1 Number density of the z = 0 slice before SN injection (left), during the SNR evolution at 2.5 kyr (center), and at 10 kyr (right) for an average ISM density of n = 1 cm−3. The rows show the different Mach numbers of ℳ = 0, 0.3, 1, and 3 from top to bottom. The box size is always 40 pc in each direction. Movies of these simulations can be found online.

Fig. A.2 Same as Fig. A.1, but for the n = 10 cm−3 cases. Here, only the ℳ = 3 case has a box length of 40 pc, whereas the other Mach number cases have a box length of 20 pc in each direction. Movies of these simulations can be found online.

Fig. A.3 Same as Figs. A.1 and A.2 but for the n = 100 cm−3 cases. Here, all boxes have a box length of 20 pc in each coordinate direction. Movies of these simulations can be found online.

Appendix B SN explosion close to filaments

The position of the SN in the ISM has a significant impact on its dust destruction. To make all simulations comparable and to mimic pre-SN feedback, we always placed the SN in a large void of the driven turbulence. To quantify the impact of the explosion position, we ran a simulation of the n = 100 cm⁻³ and ℳ = 3 case in which the SN explodes close to filaments, as shown in Fig. B.1, and compared it to the fiducial n = 100 cm⁻³ and ℳ = 3 case used in this paper. In this simulation, the SNR only expands up to ∼5 pc in radius, whereas it expands to ~8 pc in the case where it evolves in a void (bottom panel of Fig. A.3). The average shock velocity at 10 kyr (see Table B.1) is roughly the same between both cases. However, the SNR in a void sweeps up more gas during its evolution because it only encounters dense filaments at the end of the simulation, after it already expanded through a substantial amount of material. The dense filaments can then almost stop the shock since it was continuously slowed down by lower density regions at earlier times. The SNR close to filaments is slowed down early due to the high gas densities around it. At such early stage, its forward shock is energetic enough to penetrate the filaments and keep slowly expanding.

Fig. B.1 Number density of the z = 0 slice before SN injection (left), during the SNR evolution at 2.5 kyr (center), and at 10 kyr (right) for a case in which the SN explodes close to filaments, and ISM parameters of n = 100 cm−3 and ℳ = 3.

Fig. B.2 Dust destruction of the z = 0 slice for the case in which the SN explodes close to filaments (left) and for the case in which the SN explodes in a void (right, also shown in Fig. 4). Both simulations have an ISM density of n = 100 cm−3 and Mach number of ℳ = 3, and include either silicate (Sil) or carbonaceous (Car) grains, and grain–grain collisions with sputtering (gg+sp) or only sputtering (sp).

This difference in SNR evolution affects the dust destruction efficiency. The SNR that explodes close to filaments encounters more dust in its early evolutionary stage where its forward shock is highly energetic. Thus, the destroyed dust mass increases rapidly at the beginning of the simulation, as shown in Fig. B.2 but decreases toward the end where the shock is much slower and does not encounter as much dust anymore. The silicate and carbonaceous grains simulations show similar total dust destruction efficiencies of ∼86% (Table B.1), which is more than 10% more than for the SN explosion in a void (Table 1, 76% and 73% for the silicate and carbonaceous case, respectively). However, the total destroyed dust mass after 10 kyr is similar because the SNR in a void encounters more material.

When the SNR evolves closer to filaments, it furthermore forms a shell at earlier times since radiative cooling quickly becomes effective due to the high masses of swept-up gas at early times. Therefore, grain–grain collisions become important early on (from ∼2–3 kyr onward), resulting in a final difference to the sputtering-only simulations of 25% and 33% for silicate and carbonaceous grains, respectively, as shown in Table B.1. Grain-grain collisions are also the main reason why carbonaceous grains are similarly well destroyed as silicates in this simulation.

All in all, the dust destruction efficiency of the SNR close to filaments is higher than when it explodes in a void, but similar to the homogeneous and weakly turbulent n = 100 cm⁻³ cases after 10 kyr of evolution (87% for the silicate and carbonaceous grain simulations). Hence, the filaments cannot shield the dust from destruction at early times if the SN explodes closer to the dense filamentary structures of ISM turbulence. Instead, they then lead to one more solar mass of destroyed dust than in the homogeneous ISM simulation. However, the SNR at filaments is also slowed down quicker than the homogeneous case SNR so that the destruction efficiency and destroyed dust mass are smaller than in the homogeneous case at similar end shock velocities (see Appendix F). After 1 Myr of SNR evolution, the SNR close to filaments might therefore result in lower dust destruction than for the homogeneous high-density case.

Appendix C Resolution of the dust simulations

Different Cartesian grid resolutions were investigated for the Paperboats simulations and 900 × 900 × 1 was found to be optimal for our setup as it did not introduce significant errors while keeping the computational time feasible. However, at the last few time steps (∼9–10 kyr), this resolution is not sufficient anymore to resolve the forward shock, which gets thinner due to the expansion of the shock. Thus, the dust transport is not sufficiently resolved anymore, leading to dust grains that overtake the forward shock when filaments are encountered since the gas in the shock front is more efficiently stopped than the dust. We carefully investigated this phenomenon and found that it is not of physical origin but disappears with a higher considered resolution of 4500 × 4500 × 1 cells, which would require on the order of 25 times more computational resources, which is beyond our capabilities. The maximum possible resolution with our setup was 1800 cells per box side, which we tested against the 900 cells simulation for the high-density supersonic Mach number simulation. With the higher resolution, the majority of the dust remains inside of the SNR; however, at some parts with an especially dense and thin shell, the dust still overtakes the forward shock. Comparing these two resolutions, the dust destruction only changes by 0.28%, meaning that this phenomenon introduces a negligible error compared to other more significant uncertainties in this study. Since this phenomenon only occurs in the last few time steps and therefore only hardly changes the dust destruction over the full 10 kyr, we chose to carry out the analysis with 900 cells per side, and we note that the last one thousand years should be evaluated more carefully for this study.

Appendix D Influence of the surrounding ISM structure

The SNe in this study were consistently injected in large void regions. The exact distance to nearby filaments and the detailed structure of the ambient ISM, however, still influence the dust destruction by a significant factor because they determine the velocity at which the forward shock encounters high-density and dust-rich regions. Our dust simulations with Paperboats only consider three orthogonal 2D slices of the 3D SNR simulations from which we then take the average value for our results. Since the turbulence in these slices has the same statistical properties and the injection region is always a void, the dust destruction difference between these slices can be used to estimate the influence of the detailed structures found in the ambient medium. The difference in dust destruction efficiency between these slices ranges from 0.01% up to 9.4% (Fig. D.1), proving that this can have a large impact, especially for the high Mach number simulations since the slice structure can differ more significantly there. The dust destruction difference in one snapshot can increase even more if the SNe are injected in regions with higher densities than voids, as shown by Kirchschlager et al. (2024a).

Fig. D.1 Difference in dust destruction between three orthogonal planes of the same 3D MHD simulation (n = 10 cm−3 and ℳ = 3). Although the turbulence always has the same statistical properties, the dust destruction varies by roughly 9.4%, depending on the exact density distribution and close-by filamentary shapes.

Appendix E Grain-size distributions

The grain-size distributions of all simulations are shown in Fig. E.1. The black lines show the initial MRN distribution for the carbonaceous and silicate dust, which is altered due to the destruction by the forward shock. The final distributions are similar for the different Mach number cases but differ significantly between the ISM densities. In a low-density environment, grain–grain collisions are inefficient. Thus, large grains above ~40 nm mostly survive. However, a large fraction of intermediate-sized grains down to 5 nm is sputtered so that the final distribution is altered. Grains that are smaller than 5 nm are not included in the initial distribution but a significant number is generated by the different destruction mechanisms. The final grain-size distributions of the intermediate-density cases are similar but with more destroyed dust. Thus, the difference between initial and final distribution increases. Furthermore, the carbonaceous grain-size distributions show a variation between the different Mach numbers, as large carbonaceous grains are better destroyed in a strongly turbulent medium than in low Mach number turbulence. Only the high-density ISM simulations deviate significantly from this distribution. There, grain–grain collisions become important because of the formation of the thin and dense shell. Therefore, many large grains are fragmented to smaller sizes so that intermediate-sized grains are less efficiently destroyed than larger grains. In addition, many more grains smaller than 5 nm are generated than in the lower density cases. However, if only sputtering is considered, the final grain-sizes are distributed similarly to the intermediate-density cases.

Fig. E.1 Grain-size distributions at 10 kyr compared to the initial MRN distribution (black) for the carbonaceous and silicate grains. The y-axis represents the dust mass density per particle radius, which is the number of dust grains per volume and grain radius dn/da multiplied with the particle mass m. All z = 0 simulation parameter sets are shown, including the different considered turbulent Mach numbers (ℳ) and average ISM densities (n). The plots differentiate between carbonaceous (blue) and silicate (orange) dust, and the grain–grain collisions and sputtering (solid), and the sputtering-only (dotted) simulations.

Appendix F Dust destruction results with similar forward shock end velocities

Table F.1 shows the dust destruction results for all simulations of this paper, similar to Table 1, but ending at the same forward shock velocity averaged over the whole forward shock, instead of the same evolution time. This is a useful comparison between the various Mach number and ISM density cases investigated because the SNR evolves to different evolutionary stages with different shock velocities in our different models. However, most conclusions gathered from Table 1 still hold true in this comparison so that we only attach it in the appendix. The most significant change is that silicates are now always more efficiently destroyed than carbonaceous grains and that the dust destruction increases with higher ISM densities for all simulations without exceptions. Furthermore, the destroyed dust masses decrease with higher Mach numbers in this comparison, contrary to Table 1 where they remain roughly the same. Other found trends stay consistent.

Appendix G Total surface area of dust grains

If grain growth of dust grains in the ISM is found to be efficient, the total grain surface area in the post-shock region is an important parameter to evaluate as these grains could grow in the ISM once they are formed, dependent on their surface area. In all of our parameter sets, we observe that the total number of dust grains significantly increases because grain–grain collisions generate many small grains (see Appendix E). The total grain surface area, however, has a more complex evolution and can decrease if too much dust is destroyed but also increase in the case that grain–grain collisions are fragmenting enough larger dust grains to smaller sizes with a higher total grain surface area but sputtering is not efficient enough to destroy these grains. Thus, the ratio between the total grain surface area at the end to the initial value at the beginning is dependent on the considered parameter set, as can be seen in Table G.1. The total grain surface area of silicate grains is always found to decrease; however, it decreases much less for the high-density simulations in which grain–grain collisions are efficient. For carbonaceous grains, the total grain surface area increases in the high-density cases and the intermediate-density case with a strongly turbulent ISM since grain–grain collisions are more efficient for carbonaceous grains than for silicates (see Sect. 3.1). All other parameter sets show a decrease of carbonaceous grain surface area after 10 kyr compared to the initial dust distribution because grain–grain collisions are not that efficient in these simulations. Together with efficient grain growth in the ISM, this could therefore indicate a significantly lower destruction fraction of the SNR forward shock after the most destructive 10 kyr in intermediate- and high-density cases where grain–grain collisions are important. For some environments, accretion could even be up to almost three times more efficient after the SNR shocks the ISM than before the SN event.

All Tables

All Figures

Fig. 1 Structural comparison of the dust distribution in a turbulence snapshot of our intermediate-density (10 cm−3) supersonic turbulence (ℳ = 3) simulation (left) and an observed SPIRE 250 µm dust map from Herschel observations of the Polaris Flare (right, originally observed by Miville-Deschênes et al. 2010)1, which has a Mach number of ℳ = 2–10 (Beattie et al. 2019). The physical length scale of 20 × 20 pc of the observed ISM is based on the distance of 350 pc to the Polaris Flare (Schlafly et al. 2014; Panopoulou et al. 2022).

In the text

	Fig. 2 Density PDFs of the different turbulent Mach numbers ℳ in the converged low-density (1 cm−3) turbulence simulations. With a higher turbulent Mach number, cell densities (ρ) can have a larger deviation from the average density <ρ>.
In the text

	Fig. 3 Number densities of the z = 0 slice of all simulations carried out in this paper at 10 kyr. This includes different average ISM number densities (n) and turbulent Mach numbers (ℳ). The colorbar scales are consistent in each of the average ISM number density cases but the box sizes can vary between 20 or 40 pc.
In the text

Fig. 4 Dust destruction efficiencies (left y-axis) and masses (right y-axis) in the z = 0 slice. This destruction efficiency refers to the ratio of the destroyed dust mass at each time step to the total encountered dust mass at the last considered snapshot. Different Mach numbers of ℳ = 0, 0.3, 1, and 3, and average ISM densities of n = 1, 10, and 100 cm−3 are shown. For each Mach number and ISM density parameter set, four different dust simulations were carried out, including either silicate (Sil) or carbonaceous (Car) dust, and either grain–grain collisions and sputtering (gg+sp) or only sputtering (sp).

In the text

Fig. 5 16 × 16 pc zoomed-in image of a filament of the n = 1 cm−3 and ℳ = 1 simulation at t = 0 yr and t = 10 kyr. The filament cannot withstand the shock efficiently, but is destroyed leaving behind finger-like structures that arise from the forward shock front into the SNR interior. Furthermore, some high-density cloudlets are formed (circled in black). The arrow indicates the approximate direction of the forward shock.

In the text

	Fig. 6 Dust destruction efficiencies in the z = 0 slice in a direct comparison between the different Mach numbers (ℳ) used for each ISM density (n) case, considering either carbonaceous (Car) or silicate (Sil) grains.
In the text

	Fig. 7 Zoomed-in 9 × 9 pc image of a filament of the n = 100 cm−3 and ℳ = 3 simulation at t = 0 yr and t = 10 kyr. The dense filament almost stops the shock and dictates where it will evolve. The initial filament shape looks nearly identical to the forward shock shape at the end. The arrow indicates the approximate direction of the forward shock.
In the text

	Fig. A.1 Number density of the z = 0 slice before SN injection (left), during the SNR evolution at 2.5 kyr (center), and at 10 kyr (right) for an average ISM density of n = 1 cm−3. The rows show the different Mach numbers of ℳ = 0, 0.3, 1, and 3 from top to bottom. The box size is always 40 pc in each direction. Movies of these simulations can be found online.
In the text

	Fig. A.2 Same as Fig. A.1, but for the n = 10 cm−3 cases. Here, only the ℳ = 3 case has a box length of 40 pc, whereas the other Mach number cases have a box length of 20 pc in each direction. Movies of these simulations can be found online.
In the text

	Fig. A.3 Same as Figs. A.1 and A.2 but for the n = 100 cm−3 cases. Here, all boxes have a box length of 20 pc in each coordinate direction. Movies of these simulations can be found online.
In the text

	Fig. B.1 Number density of the z = 0 slice before SN injection (left), during the SNR evolution at 2.5 kyr (center), and at 10 kyr (right) for a case in which the SN explodes close to filaments, and ISM parameters of n = 100 cm−3 and ℳ = 3.
In the text

Fig. B.2 Dust destruction of the z = 0 slice for the case in which the SN explodes close to filaments (left) and for the case in which the SN explodes in a void (right, also shown in Fig. 4). Both simulations have an ISM density of n = 100 cm−3 and Mach number of ℳ = 3, and include either silicate (Sil) or carbonaceous (Car) grains, and grain–grain collisions with sputtering (gg+sp) or only sputtering (sp).

In the text

	Fig. D.1 Difference in dust destruction between three orthogonal planes of the same 3D MHD simulation (n = 10 cm−3 and ℳ = 3). Although the turbulence always has the same statistical properties, the dust destruction varies by roughly 9.4%, depending on the exact density distribution and close-by filamentary shapes.
In the text

Fig. E.1 Grain-size distributions at 10 kyr compared to the initial MRN distribution (black) for the carbonaceous and silicate grains. The y-axis represents the dust mass density per particle radius, which is the number of dust grains per volume and grain radius dn/da multiplied with the particle mass m. All z = 0 simulation parameter sets are shown, including the different considered turbulent Mach numbers (ℳ) and average ISM densities (n). The plots differentiate between carbonaceous (blue) and silicate (orange) dust, and the grain–grain collisions and sputtering (solid), and the sputtering-only (dotted) simulations.

In the text