Dust Grain Growth and Dusty Supernovae in Low-metallicity Molecular Clouds

Most massive stars reside within young massive star clusters (Lada & Lada 2003; Portegies Zwart et al. 2010). There, hundreds to thousands of massive stars inject large amounts of metal-enriched gas via stellar winds and supernova explosions. Consequently, a localized dust enrichment of galaxies originating from massive star clusters can be expected (Consiglio et al. 2016). An important amount of this dust is likely produced after the efficient condensation of the ejecta of core-collapse supernovae (e.g., Todini & Ferrara 2001; Indebetouw et al. 2014, and references therein) and pair-instability supernovae (PISNe; Nozawa et al. 2003; Cherchneff 2010). A competing mechanism that may massively enhance the amount of dust is the accretion of gas species onto already existing dust grains (Dwek 1998; Calura et al. 2008; Zhukovska et al. 2008; Asano et al. 2013; Calura et al. 2014). Both processes, dust produced by supernovae and dust grain growth, are expected to play an important role in local and high-redshift galaxies. For instance, Gall & Hjorth (2018) concluded that the amount of dust observed in high-redshift galaxies comes from an efficient supernova dust production, followed by a rapid dust grain reformation if dust is destroyed by high-velocity shocks (see also Priestley et al. 2021). Observations of Lyα systems and quasar host galaxies at high redshifts suggest a rapid transition (on the order of a few megayears) from having a dust-poor to a dust-rich interstellar medium (e.g., Michałowski et al. 2010; Mattsson et al. 2015), particularly if star formation took place with a top-heavy initial mass function (Dwek & Cherchneff2011; Gall et al. 2011). It is also worth mentioning that the cold envelopes of the most massive asymptotic giant branch stars could also be important dust producers at high redshifts (e.g., Valiante et al. 2009; Leśniewska & Michałowski 2021). Moreover, dust formation from early-type carbon-rich Wolf–Rayet binaries could also represent an important source of dust in low-metallicity environments (Z ≤ 0.65 Z_⊙; Lau et al. 2021).

The 3D hydrodynamical evolution of dusty supernova remnants originating within young massive stellar clusters was explored in Martínez-González et al. (2018). For that purpose, Cinder (Cooling INduced by Dust & Erosion Rates), a module for the adaptive mesh refinement code Flash (Fryxell et al. 2000), was developed. With Cinder one can follow the survival rate of dust condensed out of a supernova ejecta and subjected to multiple shock processing, which includes the passage of the reverse shock and its bouncing back, the interaction with shocked stellar winds, and the crossing of sequential supernova forward shocks. In the case of off-centered supernova remnants, facing a steep density gradient, a blowout phase is triggered: The supernova remnants elongate in preferential directions, become Rayleigh–Taylor unstable and the supernova ejecta suffers a rapid decline in density and temperature (Tenorio-Tagle et al. 2015; Jiménez et al. 2021). The main result of the "pyroclastic blowout model" is that clustered supernova explosions are likely to cause a net increase in the amount of dust in the free-wind region surrounding their parental stellar clusters.

However, the mechanical feedback provided by the cluster leads to the creation of a wind-blown superbubble (see Weaver et al. 1977; Mac Low & McCray 1988; Koo & McKee 1992; Bisnovatyi-Kogan & Silich 1995), where a supershell of swept-up interstellar matter encompasses the hot and tenuous free/shocked wind region. Thus, one cannot exclude that supernova remnants originating from the star cluster overrun the free and shocked wind regions, to then impact directly the swept-up supershell. This results in a further processing of both the ejecta dust and the dust grains present in the supershell. Furthermore, supershells represent a fertile environment for the growth of dust grains even in high-redshift (z ∼ 6) galaxies at supersolar metallicities (Martínez-González et al. 2021).

The collision of supernova remnants with a wind-driven shell has been previously explored (Tenorio-Tagle et al. 1990, 1991; Franco et al. 1991; Różyczka et al. 1993; Dwarkadas 2005, 2007; van Marle et al. 2015; Haid et al. 2016). In Martínez-González et al. (2019), it has been shown that the preexistent dust locked up in a wind-driven shell is mostly not in danger of being significantly destroyed after the collision of blast waves with the surrounding shell. This is expected as blast waves would transmit weakly into the much denser encompassing shell.

In the present work, we extend the models presented in Martínez-González et al. (2018, 2019), hereafter referred to as Papers I and II, and follow the wind-blown superbubble evolution in the context of low-metallicity clumpy molecular clouds. Those low-metallicity environments may be representative of the Green Pea galaxies, a local analog class of high-redshift Lyα emitting galaxies with a high star formation rate (Cardamone et al. 2009; Micheva et al. 2017; Svoboda et al. 2019; Franeck et al. 2022). They are also interesting since an increased presupernova feedback (harder ionizing radiation and increased photon fluxes) is expected in low-metallicity environments (e.g., McLeod et al. 2021). We will consider the potential well of the molecular cloud and the star cluster, the growth of grains originally residing within clumps, as well as dust grains produced in multiple pair-instability and core-collapse supernova explosions.

The Paper is organized as follows: In Section 2 we describe the characteristics of the host molecular cloud (2.1), the central star cluster and the wind-driven superbubble (2.2), the growth of dust grains (2.3), and the occurrence of supernovae (2.4). In Sections 3–5, we present the results of the hydrodynamical simulations regarding the early evolution of the superbubbble, and the dust mass evolution. In Section 6 we highlight the limitations of our model. Finally, in Section 7 we summarize our results and present the main conclusions.

2.1. Clumpy Molecular Cloud

We have modeled a hierarchical molecular cloud consisting of a collection of small clumps and a tenuous interclump envelope (Tenorio-Tagle et al. 2006; Draine 2011). It is assumed that the cloud's metallicity is Z = 0.02 Z_⊙, and that the cloud's clumpy density field is determined by Reynolds et al. (2019)

$$\begin{eqnarray}&&{\rho }_{\mathrm{cl}}({\boldsymbol{x}})={\rho }_{0}\left(1+\sum _{i=1}^{{n}_{\mathrm{cl}}}\displaystyle \frac{{\rho }_{i}}{{\rho }_{0}}{e}^{-2{\left(\parallel {\boldsymbol{x}}-{{\boldsymbol{x}}}_{i}\parallel /{r}_{i}\right)}^{2}}\right),\end{eqnarray} \tag{ 1 }$$

where x = (x, y, z), ρ₀ is the interclump gas mass density, n_cl is the number of clumps, ρ_i and r_i are the mass density and radius of the i th spherical clump centered at x _i = (x_i , y_i , z_i ). Hereafter, we assume that n_cl = 1000 and adopt for the interclump gas density ρ₀ = 2.128 × 10⁻²² g cm⁻³, with a mean mass per ion 1.22m_H , and a mean mass per molecular gas 2.33m_H , where m_H is the proton mass.

Each clump is initially assumed to have a turbulent velocity field, as implemented by Taylor et al. (2018). The clumps' kinetic energy spectra take the form E(k) ∝ k^α , with wavenumbers k between 2 × 2π/L and 32 × 2π/L, where L is the length of the computational domain, and α = −5/3. The virial ratio of the clumps is assumed to be 0.9. We have taken normally distributed values of x_i , y_i , z_i , ρ_i /ρ₀, and r_i . The mean and standard deviation of x_i , y_i , and z_i are 0 pc and 20 pc. For ρ_i /ρ₀ they are 66 and 50, respectively. For r_i these are 1.0 pc and 0.76 pc, respectively. The clumps amount to M_cl ∼ 3.73 ×10⁶ M_⊙ of molecular gas (roughly one third of the total gas mass ∼1.03 × 10⁷ M_⊙), with ∼68% of the clumps located within 13 pc.

2.2. The Star Cluster and the Wind-driven Superbubble

We have set up the 3D hydrodynamical evolution of a superbubble within a molecular cloud with the characteristics described in Section 2.1 (see Table 1), using a similar setup to that described in Papers I and II. The superbubble is driven by the mechanical feedback of a coeval young massive stellar cluster located at the cloud's center. To this end, we have used the adaptive mesh refinement code FLASH v4.3 (Fryxell et al. 2000). The employed hydrodynamical solver is a modified version of the Piecewise Parabolic Method introduced by Colella & Woodward (1984). Similarly to Paper I, the simulations take into account the star cluster's gravitational potential and the self-gravity of the gas calculated by the tree-based solver developed by Wünsch et al. (2018), and the optically thin cooling function for gas at temperatures T ≥ 10⁴K (Schure et al. 2009), and for gas at temperatures T < 10⁴ K—gas cooling due to radiation emitted from hydrogen-deuteride ground-state rotational transition (Dalgarno & McCray 1972). We have used Cinder (Martínez-González et al. 2018, 2019) to follow the injection and destruction of dust grains due to thermal sputtering, as well as the cooling induced by gas–grain collisions in hot plasmas (T ≳ 10⁵ K). As in Papers I and II, we have not included other grain destruction mechanisms such as kinetic sputtering and grain shattering (see Section 6). Additionally, the process of dust grain growth is incorporated into Cinder for the first time.

Table 1. The Host Molecular Cloud and the Star Cluster

Model	MMC	Z	ncl	ρ0	MSC	Rch	RSC		No. Stars	Grain Growth	Ejecta Dust
	(M⊙)	(Z⊙)	−	(g cm−3)	(M⊙)	(pc)	(pc)	(%)	(≥ 71 M⊙)	−	−
GrainGrowth	107	0.02	1000	2.26 × 10−22	5.6 × 105	1.0	3.0	3.7	186	Yes	No
SNDust	107	0.02	1000	2.26 × 10−22	5.6 × 105	1.0	3.0	3.7	186	No	Yes

Note. The table summarizes the main properties of the host molecular cloud (mass, metallicity, number of molecular clumps, and the interclump mass density), and the star cluster (stellar mass, the stellar distribution's characteristic scale, and cutoff radii, the global star formation efficiency, and number of progenitor stars that end as PISNe). The last two columns indicate whether the processes of dust grain growth and dust injection by supernovae are considered.

Download table as:  ASCIITypeset image

The simulations are inscribed into cubes (40 pc)³ and (140 pc)³ (when indicated) in a grid 256³ in Cartesian geometry, with outer boundary conditions set to outflow.

We follow the evolution of the star cluster wind mechanical output, mass-loss rate, and terminal speed using the BoOST stellar model grids and the SynStars stellar population synthesis code (Szécsi et al. 2022; Franeck et al. 2022). In particular, we rely on a stellar grid with metallicity 0.02 Z_⊙ (originally presented by Szécsi et al. 2015), which consists of slowly rotating single stars computed with the "Bonn" stellar evolution code with standard wind mass-loss prescriptions (Vink et al. (2000, 2001) type mass loss for O and B stars and Nieuwenhuijzen & de Jager (1990) mass loss in the supergiant phase). As discussed by Szécsi & Wünsch (2019), these models spend their main-sequence lifetimes as hot OB stars, and their post-main-sequence as cool supergiants; none of them form Wolf–Rayet stars, as the metallicity is too low for the self-stripping by the wind. The stellar winds of these massive stars do contribute to the mechanical energy inserted into the cluster, but neither the wind yields, nor the supernova yields from these models are included into the 3D hydrodynamical simulations.

We apply a Kroupa initial stellar mass function (IMF; Kroupa et al. 2001) in the mass interval [0.01–120] M_⊙. The adopted stellar mass ⁶ is 5.6 × 10⁵ M_⊙, distributed according to a Schuster stellar density profile ${\rho }_{\star }\propto {[1+{(r/{R}_{\mathrm{ch}})}^{2}]}^{-\beta }$ with characteristic scale radius R_ch = 1 pc, cutoff radius R_SC = 3 pc, and index β = 1.5 (Palouš et al. 2013; Martínez-González et al.2016, 2017). The global star formation efficiency is = M_SC/M_MC = 3.7%.

2.3. Dust Grain Growth

We have incorporated the process of dust grain growth via the accretion of gas species (Spitzer 1978), so individual dust grains increase their mass at a rate ${\dot{m}}_{\mathrm{gr}}=3{m}_{\mathrm{gr}}\dot{a}/a$, where m_gr is the grain mass and $\dot{a}$ is the rate of increase in their radius a. Similarly to Paper II, we consider graphite and silicate dust grains and a log-normal grain size distribution of the form $\sim {a}^{-1}\exp \{-0.5{[\mathrm{log}(a/{a}_{0})/\sigma ]}^{2}\};$ with a₀ = 0.1 μm and σ =0.7, and minimum and maximum grain sizes ${a}_{\min }=5\,\mathrm{nm}$ and ${a}_{\max }=0.5$ μm, respectively. The grain size distribution is sampled by 10 logarithmically spaced, size bins. The dust-to-gas mass ratio, ${{ \mathcal D }}^{(i)}$, in a bin with representative grain size, ${a}_{m}^{(i)}$, at time t + Δt is (Valentini & Brighenti 2015)

$$\begin{eqnarray}&&{{ \mathcal D }}^{(i)}(t+{\rm{\Delta }}t)={{ \mathcal D }}^{(i)}(t)\left(1+\displaystyle \frac{3| \dot{a}| {\rm{\Delta }}t}{{a}_{m}^{(i)}}\right),\end{eqnarray} \tag{ 2 }$$

where $\dot{a}$ is given by

$$\begin{eqnarray}&&\dot{a}=\displaystyle \frac{{f}_{\mathrm{gr}}{\rho }_{r}S}{4{\rho }_{\mathrm{gr}}}\left[{f}_{{M}_{d}}^{(i)}{{ \mathcal D }}^{\max }-{{ \mathcal D }}^{(i)}\right]{\left(\displaystyle \frac{8{k}_{B}T}{\pi {\mu }_{r}}\right)}^{1/2}.\end{eqnarray} \tag{ 3 }$$

In Equation (3) ρ_r is the local refractory gas mass density, ρ_gr is the grain's bulk density, S is the sticking coefficient, f_gr is the fraction of grain species (for simplicity f_gr is set to 0.5 for both graphite and silicate grains), ${f}_{{M}_{d}}^{(i)}$ is the initial dust mass fraction of the i th grain size, k_B and μ_r = 18.11m_H are the Boltzmann constant and the mean mass per refractory atom (i.e., excluding noncondensable gas species); ${{ \mathcal D }}^{\max }=2.1\times {10}^{-4}$ is the dust-to-gas mass ratio in the case all refractory elements are locked up onto dust grains. This value is taken from the initial mass fractions of refractory elements in the "Bonn" stellar evolution code. At the explored molecular gas metallicity (Z = 0.02 Z_⊙), refractory elements make up a gas mass fraction of ≲10⁻⁶, a quantity that is well approximated by the relation ${{ \mathcal D }}^{\max }\,\approx {10}^{-2}\times Z$ (Martínez-González et al. 2021).

We allow dust grain growth to occur at gas temperatures T ≤ 1000 K (Nozawa et al. 2012). At these temperatures, grains with radius ≥5 nm embedded into a gas with molecular gas mass density ≤10⁵ cm⁻³ have an equilibrium temperature ≤30 K. Consequently, the grain sticking coefficient S is approximated as 1.0 (Ferrara et al. 2016).

2.4. Supernovae

At their oxygen-burning phase, massive CO stellar cores (≳65 M_⊙) undergo the so-called electron–positron pair-creation instability (Fowler & Hoyle 1964). This leads to the disruption of the whole star as a PISN, leaving behind neither a black hole, nor a neutron star (Kozyreva et al. 2014, and references therein). According to the BoOST low-metallicity stellar tracks, this occurs tracks occur for stars with masses ≳71 M_⊙ at the zero-age main sequence. From the assumed IMF, we thus follow the occurrence of 186 PISNe between 2.55 and 3.18 Myr. The supernova rate is a function of the star cluster mass, the IMF, and the metallicity of the parent cloud, and thus evolves with time. We, however, have set a constant supernova rate for progenitors in the mass intervals [71–80) M_⊙, [80–90) M_⊙, [90–100) M_⊙, [100–110) M_⊙, and [110–120] M_⊙, with one supernova per 4250, 3300, 3050, 2880, and 2460 yr, respectively. Supernovae are randomly distributed following the assumed stellar density profile. As PISNe completely obliterate their progenitor stars, their ejecta masses, M_ej , correspond to the final masses obtained from the BoOST stellar model grids. With respect to the total energy released in a single PISN event, we assume ${E}_{\mathrm{SN}}=5\times {10}^{51}$ erg. This energy is inserted into a sphere of radius equivalent to five grid cells. Massive stars with final masses ∼[10–65] M_⊙ explode as core-collapse supernovae (Szécsi et al. 2022). We have thus explored the occurrence of the first 50 core-collapse supernovae starting after the last PISN. Our assumption is that each core-collapse supernova releases M_ej = 10 M_⊙ of gas and ${E}_{\mathrm{SN}}={10}^{51}\,{\rm{erg}}$ in the form of kinetic energy.

$$\begin{eqnarray}&&{\rho }_{\mathrm{ej}}=\displaystyle \frac{(3-\omega )}{4\pi }\displaystyle \frac{{M}_{\mathrm{ej}}}{{R}_{\mathrm{SN}}^{3}}{\left(\displaystyle \frac{{R}_{\mathrm{SN}}}{r^{\prime} }\right)}^{\omega },\end{eqnarray} \tag{ 4 }$$

where $r^{\prime}$ is the radial distance from the center of the explosion and ω = 5/2 was set for all supernovae. The initial ejecta velocity distribution is assumed to be

$$\begin{eqnarray}&&{v}_{\mathrm{ej}}={\left(2\displaystyle \frac{(5-\omega )}{(3-\omega )}\displaystyle \frac{{E}_{\mathrm{SN}}}{{M}_{\mathrm{ej}}}\right)}^{1/2}\left(\displaystyle \frac{r^{\prime} }{{R}_{\mathrm{SN}}}\right).\end{eqnarray} \tag{ 5 }$$

In our first model GrainGrowth, we follow the occurrence of dust grain growth during the evolution of the molecular cloud and the superbubble driven by the central star cluster. Given the low gas temperatures involved in the host molecular cloud, and hence short cooling timescales, and small cooling lengths, we have modeled the central (40 pc)³ at three spatial resolutions, 0.15 pc and 0.31 pc, respectively. Additionally, we tested the results with a larger computational domain, (140 pc)³ at a resolution 0.54 pc, which encompasses all the molecular clumps. In this model we initially assume that half the mass of refractory elements within molecular clumps are already depleted onto dust grains, so that their dust-to-gas mass ratio is ${ \mathcal D }=\tfrac{1}{2}{{ \mathcal D }}^{\max }$ (and zero elsewhere). Thus the initial dust mass within molecular clumps is $\tfrac{1}{2}{{ \mathcal D }}^{\max }{M}_{\mathrm{cl}}\approx 373\,{M}_{\odot }$. The initial population of dust grains is assumed to have been originated in prior (very early) generations of supernovae (Nozawa et al. 2012). We stop our calculations just before the occurrence of the first PISN.

As the molecular cloud evolution proceeds, clumps become elongated and end up forming filamentary structures while they revolve around the gravitational potential well. We note that the potential well is dominated by the molecular gas, rather than the central star cluster. During the process, the clumps interact, collide, and coalesce with other clumps (Elmegreen 1988). The superbubble thus expands preferentially along the paths of least resistance, i.e., through the lower density channels between clumps/filaments (Tenorio-Tagle et al. 2006; Alūzas et al. 2012; Rogers & Pittard 2013; Lucas et al. 2020; Lancaster et al. 2021a, 2021b, see Figure 1).

Zoom In Zoom Out Reset image size

The end result is that the superbubble morphology becomes highly distorted. The flow of thermalized matter between clumps/filaments leads first to the distortion of the star cluster wind-driven shell (displayed in the bottom panels of Figure 1 by the blue/cool contours encompassing the red/hot zones) and then to the development of multiple shells; some of them appearing as nested shells when they are seen projected against the background (see Figures 2 and 3 in Tenorio-Tagle et al. 2006). Those multiple shells are all interconnected, thus forming a large-scale supershell. Eventually, the superbubble, powered mostly by the SN shocks, will breakout/blowout from the host molecular cloud (Tenorio-Tagle 2002); however this moment is not captured in our simulations within our computational domain.

Zoom In Zoom Out Reset image size

The top and bottom panels in Figure 3 present the evolution of the gas column density, N, and dust mass surface density, Σ_dust, integrated along the line of sight at 100 kyr and 2.5 Myr for model GrainGrowth. The right panels in the Figure correspond to 3.2 Myr in the case of model SNDust, so Σ_dust corresponds to supernova-condensed dust only.

Zoom In Zoom Out Reset image size

In the model GrainGrowth, dust grain growth is promoted when the molecular clumps become filamentary and mix with the interclump material. It is also favored when interclump gas is incorporated into the dense shells formed as the hot gas streams through low-density channels. We note that as the shells cool down quickly, a negligible amount of swept-up dust is destroyed. As a result, the dust mass increases at a rate ∼4.8 × 10⁻⁵ M_⊙ yr⁻¹, leading to a net increase of ∼120 M_⊙ within 2.5 Myr (see panel (a) in Figure 4). The results with resolutions 0.31 pc and 0.54 pc show a reduction of 13% and 25%, respectively, on the amount of dust grain growth. This is expected since decreasing the resolution inhibits efficient mixing of dust-free and dust-rich gas via numerical viscosity.

Zoom In Zoom Out Reset image size

As we have set the outer boundary conditions to outflow, there is some fraction of the dust mass that leaves the computational domains. However, the majority of the gas/dust overdensities remain within the central part of the molecular cloud.

We note that if dust-induced cooling of molecular gas were at play (not included in the current version of CINDER), the shells would be much thinner, denser, and much colder, and thus they may promote a much more efficient growth of dust grains at higher gas metallicities (Martínez-González et al. 2021). Higher resolution simulations are required in such a case.

Model SNDust explores the injection of dust by sequential PISNe and core-collapse supernovae, while molecular clumps are set to be dust free (i.e., no grain growth occurs). This model is studied at a resolution 0.54 pc since the cooling length ⁷ in the cavity is very large due to its low density and high temperature.

It is assumed that each PISN leads to the condensation of 3.2% of the progenitor's final mass (around 3.25 M_⊙ for a 120 M_⊙ progenitor star) onto dust grains (see Figure 2). Such a high value is consistent with theoretical expectations for dust condensation in PISNe occurring from zero-metallicity progenitors (Cherchneff 2010). The dust mass condensed out of supernova ejecta for core-collapse supernovae is assumed as normally distributed, with a mean of 0.6 M_⊙ and standard deviation of 0.1 M_⊙. Similarly to the grains originally locked up within molecular clumps, the simulated supernovae are also assumed to inject graphites and silicates in equal mass fractions following a log-normal grain size distribution. We have not included the process of ion trapping onto dust grains (see Kirchschlager et al. 2020), although such a process can produce a further growth of ejecta-condensed dust grains.

Panel (b) in Figure 4 shows that, despite the occurrence of several PISN explosions, the amount of ejecta-condensed dust increments steadily, at a rate 10⁻⁴ M_⊙ yr⁻¹ between 2.55 and 3.18 Myr. To put it into perspective, we find that from 470 M_⊙ of ejecta dust that were injected during the era of PISNe, about 13% (∼ 61 M_⊙) remain in the superbubble's cavity ∼320 kyr after the last PISN. This value is in agreement with the conclusions raised in Paper I, which explored, at a higher resolution (∼0.11 pc), the case of sequential core-collapse supernova explosions, occurring at a similar rate, within a star cluster wind. It also builds on the conclusion that massive shells represent insurmountable barriers that protect the vast majority of the preexisting surrounding dust from being shock processed (Paper II).

Progenitors with final masses ∼[10–65] M_⊙ explode as core-collapse supernovae. We have found that from the ∼30 M_⊙ of ejecta dust that were injected into the simulation, ∼32% (∼9.5 M_⊙) survive after 175,000 yr from the first core-collapse supernova explosion.

If we extrapolate the rate of net dust injection by core-collapse supernovae until ∼21 Myr of the cluster's evolution (the time needed to explode a 10 M_⊙ star in the BoOST low-metallicity stellar tracks), we can expect a net increase in the dust mass by ∼1200 M_⊙ (M_SC/5.6 × 10⁵ M_⊙). Applying this estimate to the gravitationally lensed galaxy A2744_YD4 (at a redshift z ∼ 8.38, Laporte et al. 2017), with a total stellar mass ∼2 ×10⁹ M_⊙, we obtain a net dust contribution by supernovae ∼4.3 × 10⁶ M_⊙. This value fits pretty well with the dust mass (6 × 10⁶ M_⊙) derived by Laporte et al. (2017; see also Behrens et al. 2018). Note that our estimate relies on a standard Kroupa initial stellar mass function, i.e., a top-heavy IMF was not necessary in order to explain such a large dust mass produced by pair-instability and core-collapse supernovae. We have also tested the occurrence of the first five PISNe at a resolution 0.15 pc in a domain (40 pc)³ (see the right panels in Figure 1), and the results agree within 15%.

Radiative heating provided by the central stellar cluster is not included into the simulations. We have, however, set a standard heating rate 2 × 10⁻²⁶ erg s⁻¹ (Koyama & Inutsuka 2002). This does not prevent molecular gas from cooling to low temperatures provided a sufficient density. Additionally, the gas can also cool adiabatically, which occurs in low-density regions behind shocks.

Our model does not incorporate interstellar magnetic fields that may play an important role in molecular cloud dynamics. Nevertheless, recent measurements of the magnetic field strength in molecular clouds show that gravitational forces may exceed those due to magnetic pressure gradients (Crutcher 2012; Crutcher & Kemball 2019). The presence of an interstellar magnetic field could initially increase the supershell's thickness, thus increasing the timescale for dust grain growth within the supershell (Ferriere et al. 1991).

As already stated, we have not modeled other grain destructive processes that require a shock to kinematically decouple the gas particles and dust grains, such as kinetic sputtering and grain shattering. However, the superbubble's interior is very tenuous and thus gas–grain and grain–grain encounters are not very frequent even in the case of large gas–grain and grain–grain relative velocities (see Appendix A.2 in the Paper by Martínez-González et al. 2019). This fact, and the fact that sequential supernova shocks soon decay into sound waves within the thermalized superbubble's interior (Tenorio-Tagle & Bodenheimer 1988), are likely to decrease the relative importance of such processes. In addition, betatron acceleration, which occurs when charged dust grains gyrate along magnetic field lines (Shull 1977), is not likely to play a significant role within the superbubble's interior given that at high temperatures (≳2 × 10⁵ K), dust grains tend to be neutral (McKee et al. 1987).

By conducting 3D hydrodynamical simulations, we have explored the early (∼3.2 Myr) evolution of a starburst-driven superbubble and the onset of PISNe in a low-metallicity (Z = 0.02 Z_⊙), clumpy molecular cloud. Our main purposes were to study the process of dust grain growth at low metallicities and the fate of dust grains condensed from the ejecta of pair-instability and core-collapse supernovae. The star cluster's mechanical feedback has been modeled using the state-of-the-art BoOST stellar model grids (Szécsi et al. 2022) and the SynStars stellar population synthesis code (Franeck et al. 2022).

Dense clumps that randomly move in the gravitational potential well, coalesce, and mix with the interclump matter. At the same time, the hot gas that fills the superbubble's cavity flows through the channels in between overdensities, to then create a net of interconnected shells. It was found that the mixing of clump–interclump material acts as a trigger for an important enhancement (at a rate ∼4.8 × 10⁻⁵ M_⊙ yr⁻¹ during the first 2.5 Myr of the superbubble evolution) of the dust mass in the filamentary/clumpy molecular cloud and the supershell.

Around ∼2.5 Myr of the superbubble's evolution, the most massive stars in the central star cluster start to explode as PISNe. Their forward shocks move into the shock-heated, low-density cavity, and soon they decay to sound waves (Tenorio-Tagle & Bodenheimer 1988). This behavior might inhibit the destruction of the ejecta-condensed dust grains via nonthermal, shock-induced dust grain disruptive processes. Overall, ∼13% and ∼32% of the dust masses injected by PISNe and core-collapse supernovae, respectively, survive even after being processed in multiple supernova collisions, with net dust injection rates of the order 10⁻⁴ M_⊙ yr⁻¹ and 5.4 × 10⁻⁵ M_⊙ yr⁻¹, respectively. The destruction of dust grains locked up in the shell is also largely inhibited as the multiple supernova blast waves are marginally transmitted into the supershell; a result originally presented in Paper II.

The net dust contribution by supernovae alone is sufficient to explain the large dust mass present in the gravitationally lensed galaxy A2744_YD4 (at a redshift z ∼ 8.38) reported by Laporte et al. (2017), without the need to invoke a top-heavy IMF.

We thus have demonstrated that the processes of dust grain growth and dust injection by supernovae are both efficient pathways that lead to massive amounts of dust to be present in low-metallicity environments populated by young stellar clusters.

The authors thank the anonymous referee for a careful reading and helpful suggestions, which greatly improved the paper. This study was supported by CONACYT-México research grant A1-S-28458. S.M.-G. also acknowledges support by CONACYT through project No. 482 of the "Programa Investigadoras e Investigadores por México." The authors thankfully acknowledge the computer resources, technical expertize, and support provided by the Laboratorio Nacional de Supercómputo del Sureste de México, CONACYT member of the network of national laboratories, and by the Laboratorio Nacional de Cómputo de Alto Desempeño (LANCAD), project 13-2021. R.W. and J.P. acknowledge financial support from the Czech Science Foundation project No. 19-15008S and by the Astronomical Institute of the Czech Academy of Sciences through the project RVO:67985815. D.Sz. has been supported by the Alexander von Humboldt Foundation. This research was funded in part by the National Science Center (NCN), Poland under grant number OPUS 2021/41/B/ST9/00757. For the purpose of Open Access, the author has applied a CC-BY public copyright license to any Author Accepted Manuscript (AAM) version arising from this submission.

Software: FLASH v4.3 (Fryxell et al. 2000), Numpy (Harris et al. 2020), Wind (Wünsch et al. 2017), TreeRay (Wünsch et al.2018), Cinder (Martínez-González et al. 2018).