Dust Grain Growth&Dusty Supernovae in Low-Metallicity Molecular Clouds

We present 3-D hydrodynamical models of the evolution of superbubbles powered by stellar winds and supernovae from young coeval massive star clusters within low metallicity ($Z = 0.02$Z$_{\odot}$), clumpy molecular clouds. We explore the initial stages of the superbubble evolution, including the occurrence of pair-instability and core-collapse supernovae. Our aim is to study the occurrence of dust grain growth within orbiting dusty clumps, and in the superbubble's swept-up supershell. We also aim to address the survival of dust grains produced by sequential supernovae. The model accounts for the star cluster gravitational potential and self-gravity of the parent cloud. It also considers radiative cooling (including that induced by dust) and a state-of-the-art population synthesis model for the coeval cluster. As shown before, a superbubble embedded into a clumpy medium becomes highly distorted, expanding mostly due to the hot gas streaming through low density channels. Our results indicate that in the case of massive ($\sim10^7$M$_{\odot}$) molecular clouds, hosting a super star cluster ($\sim5.6\times10^5$M$_{\odot}$), grain growth increments the dust mass at a rate $\sim4.8\times10^{-5}$M$_{\odot}$ yr$^{-1}$ during the first $2.5$Myr of the superbubble's evolution, while the net contribution of pair-instability and core-collapse supernovae to the superbubble's dust budget is $\sim1200$M$_{\odot} (M_{SC}/5.6\times10^{5}$M$_{\odot})$, where $M_{SC}$ is the stellar mass of the starburst. Therefore, dust grain growth and dust injection by supernovae lead to create, without invoking a top-heavy initial mass function, massive amounts of dust within low-metallicity star-forming molecular clouds, in accordance with the large dust mass present in galaxies soon after the onset of cosmic reionization.


INTRODUCTION
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;Zhukovska et al. 2008;Asano et al. 2013;Calura et al. 2008Calura 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 Myr) 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 (Gall et al. 2011;Dwek & Cherchneff 2011). It is also worth mentioning that the cold envelopes of the most massive AGB stars could also be important dust producers at high-redshifts (e.g. Valiante et al. 2009;Leśniewska & Michałowski 2019). Moreover, dust formation from early-type carbon-rich Wolf-Rayet binaries could also represent an important source of dust in lowmetallicity environments (Z 0.65 Z Lau et al. 2021). The three-dimensional 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, that 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 superno-va 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 . The collision of supernova remnants with a wind-driven shell has been previously explored (Tenorio-Tagle et al. 1990Franco et al. 1991;Różyczka et al. 1993;Dwarkadas 2005Dwarkadas , 2007van Marle et al. 2015;Haid et al. 2016). In Martínez-González et al. (2019), it has been shown that the pre-existent dust locked-up in a winddriven 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 as Papers I andII), 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. 2021). 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.

Clumpy Molecular Cloud
We have modelled 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) where x = (x, y, z), ρ 0 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 ith 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 ρ 0 = 2.128 × 10 −22 g cm −3 , 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 /ρ 0 and r i . The mean and standard deviation of x i , y i , and z i are 0 pc and 20 pc. For ρ i /ρ 0 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 6 M of molecular gas (roughly one third of the total gas mass ∼ 1.03 × 10 7 M ), with ∼ 68 per cent of the clumps located within 13 pc.

The Star Cluster & the Wind-driven Superbubble
We have set-up the three-dimensional 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 treebased solver developed by Wünsch et al. (2018), and the optically thin cooling function for gas at temperatures T 10 4 K (Schure et al. 2009), and for gas at temperatures T < 10 4 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 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 5 K). As in Papers I and II, we have not included other grain destruction mechanisms such as kinetic sputtering and grain shattering (see subsection 6). Additionally, the process of dust grain growth is incorporated into Cinder for the first time. 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 cut-off 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. Figure 1. The stream of hot gas through low-density channels. The upper and bottom panels display slides in the z-plane of the distribution of gas mass density and gas temperature, respectively, at three different times (left: 100 kyr; center: 2.0 Myr; right: 2.65 Myr). The flow of thermalized matter between low-temperature overdensities leads first to the distortion of the starburst-driven shell, and then to the development of multiple ragged shells at later times. Note the left and central panels correspond to model GrainGrowth, while the right panels depict model SNDust, once PISNe are occurring and the superbubble's interior is filled with wind and supernova matter. At 2.65 Myr, a young supernova remnant can be identified as a ring-like hot structure exhibiting its forward and reverse shocks, and surrounded by the free-wind region with a density that falls as ∼ r −2 . Further out, there is a global reverse shock thermalizing the free cluster wind and causing a temperature larger than 10 7 K (appearing as a red zone at ∼ 10 pc from the center). Much further out there is a leading shock that has managed to display the original fragments from the central regions, while advancing with different speeds in different directions given the uneven density left in the surrounding medium.
The simulations are inscribed into cubes (40 pc) 3 and (140 pc) 3 (when indicated) in a grid 256 3 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. 2020;Franeck et al. 2021). 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(Vink et al. ( , 2001 type mass loss in the OB phase and Nieuwen-huijzen & 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 3-D hydrodynamical simulations. We apply a Kroupa initial stellar mass function (IMF) in the mass interval [0.01-120] M . The adopted stellar The mass of the ejecta corresponds to the progenitor's final mass, and the ejecta-condensed dust mass is assumed to be 3.2 per cent of the ejecta mass (Cherchneff 2010).

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ṁ gr = 3m grȧ /a, where m gr is the grain mass andȧ 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 ∼ a −1 exp{−0.5[log(a/a 0 )/σ] 2 }; with a 0 = 0.1 µm and σ = 0.7 and minimum and maximum grain sizes a min = 5 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, D (i) , in a bin with representative grain size, a (i) m , at time t + ∆t is (Valentini & Brighenti 2015) whereȧ is given bẏ 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 1 This value is inspired by the derived mass of NGC 604 (Relaño et al. 2016). As we shall see, the simulated star cluster, superbubble and molecular cloud resemble the complex morphology of massive HII regions such as NGC 604.
f gr is set to 0.5 for both graphite and silicate grains), f (i) M d is the initial dust mass fraction of the ith grain size, k B and µ r = 18.11m H are the Boltzmann constant and the mean mass per refractory atom (i.e. excluding non-condensable gas species). D max = 2.1 × 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 −6 , a quantity that is well approximated by the relation D max ≈ 10 −2 × 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 5 cm −3 have an equilibrium temperature 30 K. Consequently, the grain sticking coefficient S is approximated as 1.0 (Ferrara et al. 2016).

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-metaillicity stellar tracks, that occurs 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 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 SN = 5 × 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. 2020). We have thus explored the occurrence of the first 50 core-collapse supernovae starting after the last PISNe. Our assumption is that each core-collapse supernova releases M ej = 10 M of gas and E SN = 10 51 erg in the form of kinetic energy. For each supernova, the initial ejecta gas mass density distribution is assumed to be given by where r 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 v ej = 2 (5 − ω) (3 − ω)

LARGE-SCALE EVOLUTION
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 modelled the central (40 pc) 3 at three spatial resolutions, 0.15 pc and 0.31 pc, respectively. Additionally, we tested the results with a larger computational domain, (140 pc) 3 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 dustto-gas mass ratio is D = 1 2 D max (and zero elsewhere). Thus the initial dust mass within molecular clumps is 1 2 D max M cl ≈ 373 M . 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 PISNe. 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,b, see Figure 1). 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. 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  In model GrainGrowth, dust grain growth is promoted when the molecular clumps become filamentary and mix with the interclump material. It is also favoured 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 −5 M yr −1 , 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 per cent, 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.
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 CIN-DER), 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 . Higher resolution simulations are required in such a case.

DUST INJECTED BY SUPERNOVAE
Model SNDust explores the injection of dust by sequential PISN and core-collpase 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 2 in the cavity is very large due to its low density and high temperature. It is assumed that each PISNe leads to the condensation of 3.2 per cent of the progenitor's final mass (around 3.25 M for a 120 M progenitor star) onto dust grains (see Figure 2). Such 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 normallydistributed, with a mean of 0.6 M and standard deviation of 0.1 M . Similarly to the grains originally lockedup 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 ejectacondensed 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 −4 M yr −1 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 per cent (∼ 61 M ) remains in the superbubble's cavity ∼ 320 kyr after the last PISN. This value is in agreement with the conclusions raised in Paper I, that 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 pre-existing 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 per cent (∼ 9.5 M ) survives after 175, 000 years from the first core-collapse supernova explosion. If we extrapolate the rate of net dust injection by corecollapse supernovae until ∼ 21 Myr of the cluster's evolution (the time needed to explode a 10 M star in the BoOST low-metaillicity stellar tracks), we can expect a net increase in the dust mass by ∼ 1200 M (M SC /5.6 × 10 5 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 9 M , we obtain a net dust contribution by supernovae ∼ 4.3 × 10 6 M . This value fits pretty well the dust mass (6 × 10 6 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 topheavy IMF was not necessary in order to explain such a large dust mass produced by pair-instability and corecollapse supernovae. We have also tested the occurrence of the first five PISNe at a resolution 0.15 pc in a domain (40 pc) 3 (see the right panels in Figure 1), and the results agree within 15 per cent.

RADIATIVE HEATING, INTERSTELLAR MAGNETIC FIELD, & GAS-DUST COUPLING
Radiative heating provided by the central stellar cluster is not included into the simulations. We have, however, set a standard heating rate 2 × 10 −26 erg s −1 (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 which 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 modelled other grain destructive processes that require a shock to kinematicallydecouple 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 Paper II). 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 5 K), dust grains tend to be neutral (McKee et al. 1987).

CONCLUDING REMARKS
By conducting 3-D hydrodynamical simulations, we have explored the early (∼ 3.2 Myr) evolution of a starburstdriven superbubble and the onset of PISNe in a lowmetallicity (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 corecollapse supernovae. The star cluster's mechanical feedback has been modeled using the state-of-the-art Bo-OST stellar model grids (Szécsi et al. 2020) and the Syn-Stars stellar population synthesis code (Franeck et al. 2021). 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 −5 M yr −1 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 evolution, the most massive stars in the central star cluster start to explode as PISNe. Their forward shocks move into the shockheated, 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 non-thermal, shock-induced dust grain disruptive processes. Overall, ∼ 13 and ∼ 32 per cent 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 −4 M yr −1 and 5.4 × 10 −5 M yr −1 , 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.