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

Sergio Martínez-González11affiliation: CONACYT-Instituto Nacional de
Astrofísica, Óptica y Electrónica, AP 51, 72000 Puebla, México:
<EMAIL_ADDRESS>, Richard Wünsch22affiliation: Astronomical Institute,
Czech Academy of Sciences, Boční II 1401/1, 141 00 Praha 4, Czech Republic ,
Guillermo Tenorio-Tagle33affiliation: Instituto Nacional de Astrofísica,
Óptica y Electrónica, AP 51, 72000 Puebla, México , Sergiy
Silich33affiliation: Instituto Nacional de Astrofísica, Óptica y Electrónica,
AP 51, 72000 Puebla, México , Dorottya Szécsi44affiliation: I. Physikalisches
Institut, Universität zu Köln, Zülpicher Strasse 77, D-50937 Köln, Germany
55affiliation: Institute of Astronomy, Faculty of Physics, Astronomy and
Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Toruń,
Poland , Jan Palouš22affiliation: Astronomical Institute, Czech Academy of
Sciences, Boční II 1401/1, 141 00 Praha 4, Czech Republic

###### Abstract

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\,\thinspace\hbox{$\hbox{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 ($\sim
10^{7}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$) molecular clouds, hosting a
super star cluster ($\sim 5.6\times
10^{5}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$), grain growth increments the
dust mass at a rate $\sim 4.8\times
10^{-5}\,\thinspace\hbox{$\hbox{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 $\sim
1200\,\thinspace\hbox{$\hbox{M}_{\odot}$}(M_{SC}/5.6\times
10^{5}\,\thinspace\hbox{$\hbox{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.

###### Subject headings:

galaxies: star clusters: general — (ISM:) dust, extinction — Physical Data and
Processes: hydrodynamics

††slugcomment: Accepted by ApJ: June 11th, 2022.

footnote

##  1\. 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., 2008, 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$\alpha$ 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 low-metallicity
environments ($Z\leq 0.65\,\thinspace\hbox{$\hbox{Z}_{\odot}$}$ 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
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\sim 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
pre-existent 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 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$\alpha$ 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 pre-supernova 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\. Model

###  2.1. 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\,\thinspace\hbox{$\hbox{Z}_{\odot}$}$, and that the cloud’s clumpy
density field is determined by (Reynolds et al., 2019)

$\rho_{cl}({\bf
x})=\rho_{0}\left(1+\sum_{i=1}^{n_{cl}}\frac{\rho_{i}}{\rho_{0}}e^{-2(\|{\bf
x}-{\bf x}_{i}\|/r_{i})^{2}}\right),$ (1)

where ${\bf x}=(x,y,z)$, $\rho_{0}$ is the interclump gas mass density,
$n_{cl}$ is the number of clumps, $\rho_{i}$ and $r_{i}$ are the mass density
and radius of the $i$th spherical clump centered at ${\bf
x}_{i}=(x_{i},y_{i},z_{i})$. Hereafter, we assume that $n_{cl}=1000$ and adopt
for the interclump gas density $\rho_{0}=2.128\times 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)\propto k^{\alpha}$, with wavenumbers $k$ between $2\times
2\pi/L$ and $32\times 2\pi/L$, where $L$ is the length of the computational
domain, and $\alpha=-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}$,
$\rho_{i}/\rho_{0}$ and $r_{i}$. The mean and standard deviation of $x_{i}$,
$y_{i}$, and $z_{i}$ are $0$ pc and $20$ pc. For $\rho_{i}/\rho_{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}\sim 3.73\times
10^{6}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ of molecular gas (roughly one
third of the total gas mass $\sim 1.03\times
10^{7}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$), with $\sim 68$ per cent of the
clumps located within $13$ pc.

###  2.2. 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 tree-based solver developed by Wünsch et
al. (2018), and the optically thin cooling function for gas at temperatures
$T\geq 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, 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\gtrsim 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.

Table 1 The Host Molecular Cloud & the Star Cluster

Model | $M_{MC}$ | $Z$ | $n_{cl}$ | $\rho_{0}$ | $M_{SC}$ | $R_{ch}$ | $R_{SC}$ | $\epsilon$ | No. Stars | Grain Growth | Ejecta dust
---|---|---|---|---|---|---|---|---|---|---|---
|  $\hbox{M}_{\odot}$ |  $\hbox{Z}_{\odot}$ | $-$ |  g cm-3 |  $\hbox{M}_{\odot}$ |  pc |  pc | per cent | $\geq 71\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ | $-$ | $-$
GrainGrowth | $10^{7}$ | $0.02$ | $1000$ | $2.26\times 10^{-22}$ | $5.6\times 10^{5}$ | $1.0$ | $3.0$ | $3.7$ | $186$ | Yes | No
SNDust | $10^{7}$ | $0.02$ | $1000$ | $2.26\times 10^{-22}$ | $5.6\times 10^{5}$ | $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 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 $\sim 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 $\sim
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\,\thinspace\hbox{$\hbox{Z}_{\odot}$}$ (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 in the OB phase 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 3–D hydrodynamical
simulations.

We apply a Kroupa initial stellar mass function (IMF) in the mass interval
$[0.01$–$120]\,\thinspace\hbox{$\hbox{M}_{\odot}$}$. The adopted stellar
mass111This 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.
is $5.6\times 10^{5}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$, distributed
according to a Schuster stellar density profile
$\rho_{\star}\propto[1+(r/R_{ch})^{2}]^{-\beta}$ with characteristic scale
radius $R_{ch}=1$ pc, cut-off radius $R_{SC}=3$ pc and index $\beta$=1.5
(Palouš et al., 2013; Martínez-González et al., 2016, 2017). The global star
formation efficiency is $\epsilon=M_{SC}/M_{MC}=3.7$ per cent.

###  2.3. Dust Grain Growth

Figure 2.— The plot shows the ejecta mass (solid line) and the assumed ejecta-
condensed dust mass (dashed line) as functions of the progenitor’s mass
(measured at the zero-age main-sequence). 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). Figure 3.— The
wind-blown superbubble’s evolution in an initially clumpy medium. The central
star cluster has a mass $5.6\times
10^{5}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$. The upper, and bottom panels
show the gas column density, and the dust mass surface density (model
GrainGrowth), respectively. The left and central and right panels correspond
to three different evolutionary times (left: $100$ kyr, center: $2.5$ Myr, for
model GrainGrowth; right: $3.2$ Myr for model SNDust). Soon, the host
molecular cloud and the supperbubble acquire complex morphologies due to the
action of gravity and the flow of hot gas through channels between
clumps/filaments. An animated version of this Figure can be found on Github.

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}_{gr}=3m_{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[\log(a/a_{0})/\sigma]^{2}\\}$; with $a_{0}=0.1$ $\mu$m and
$\sigma=0.7$ and minimum and maximum grain sizes $a_{min}=5$ nm and
$a_{max}=0.5$ $\mu$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+\Delta t$ is (Valentini & Brighenti, 2015)

$\displaystyle\mathcal{D}^{(i)}(t+\Delta t)$ $\displaystyle=$
$\displaystyle\mathcal{D}^{(i)}(t)\left(1+\displaystyle\frac{3|\dot{a}|\Delta
t}{a_{m}^{(i)}}\right),$ (2)

where $\dot{a}$ is given by

$\displaystyle\dot{a}=\frac{f_{gr}\rho_{r}S}{4\rho_{gr}}\left[f_{M_{d}}^{(i)}\mathcal{D}^{max}-\mathcal{D}^{(i)}\right]\left(\frac{8k_{B}T}{\pi\mu_{r}}\right)^{1/2}.$
(3)

In equation (3) $\rho_{r}$ is the local refractory gas mass density,
$\rho_{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 $\mu_{r}=18.11m_{H}$
are the Boltzmann constant and the mean mass per refractory atom (i.e.
excluding non-condensable 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\,\thinspace\hbox{$\hbox{Z}_{\odot}$}$),
refractory elements make up a gas mass fraction of $\lesssim 10^{-6}$, 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\leq 1000$ K (Nozawa
et al., 2012). At these temperatures, grains with radius $\geq 5$ nm embedded
into a gas with molecular gas mass density $\leq 10^{5}$ cm-3 have an
equilibrium temperature $\leq 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 ($\gtrsim
65\,\thinspace\hbox{$\hbox{M}_{\odot}$}$) 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 $\gtrsim 71\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ 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)\,\thinspace\hbox{$\hbox{M}_{\odot}$}$,
$[80$–$90)\,\thinspace\hbox{$\hbox{M}_{\odot}$}$,
$[90$–$100)\,\thinspace\hbox{$\hbox{M}_{\odot}$}$,
$[100$–$110)\,\thinspace\hbox{$\hbox{M}_{\odot}$}$, and
$[110$–$120]\,\thinspace\hbox{$\hbox{M}_{\odot}$}$, with one supernova per
$4250$ years, $3300$ years, $3050$ years, $2880$ years, and $2460$ years,
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_{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 $\sim[10$–$65]\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ 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\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ of gas and
$E_{SN}=10^{51}\,\mbox{erg}$ in the form of kinetic energy.

For each supernova, the initial ejecta gas mass density distribution is
assumed to be given by

$\displaystyle\rho_{ej}=\frac{(3-\omega)}{4\pi}\frac{M_{ej}}{R_{SN}^{3}}\left(\frac{R_{SN}}{r^{\prime}}\right)^{\omega},$
(4)

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

$\displaystyle
v_{ej}=\left(2\frac{(5-\omega)}{(3-\omega)}\frac{E_{SN}}{M_{ej}}\right)^{1/2}\left(\frac{r^{\prime}}{R_{SN}}\right).$
(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 dust-to-gas mass ratio is
$\mathcal{D}=\displaystyle\frac{1}{2}\mathcal{D}^{max}$ (and zero elsewhere).
Thus the initial dust mass within molecular clumps is
$\displaystyle\frac{1}{2}\mathcal{D}^{max}M_{cl}\approx
373\,\thinspace\hbox{$\hbox{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 PISNe.

Figure 4.— The dust mass evolution. Panel a presents the increment of dust
mass due to grain growth in the central part ($40^{3}$ pc3) of the molecular
cloud/supershell during the pre-supernova era (model GrainGrowth). The thick
line depicts the result of the run with the $0.15$ pc spatial resolution,
while the dashed line to the $0.31$ pc spatial resolution run. Panel b
displays the cumulative ejecta-condensed dust mass deposited by PISNe. The
vertical dashed lines in panel b mark the onset and finale of the PISN era.
Note the different time intervals depicted in both panels.

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, $\Sigma_{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
$\Sigma_{dust}$ corresponds to supernova-condensed dust only.

##  4\. Dust grain growth at low metallicity

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 $\sim 4.8\times
10^{-5}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ yr-1, leading to a net increase
of $\sim 120\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ 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 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.

##  5\. 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 length222The cooling length can be defined as the product of the
local sound speed and the cooling time, and increases with decreasing gas
density and metallicity (Smith et al., 2017). 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\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ for a
$120\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ 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 normally-distributed, with a mean of
$0.6\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ and standard deviation of
$0.1\,\thinspace\hbox{$\hbox{M}_{\odot}$}$. 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^{-4}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ yr-1 between $2.55$ and $3.18$
Myr. To put it into perspective, we find that from
$470\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ of ejecta dust that were injected
during the era of PISNe, about $13$ per cent ($\sim
61\,\thinspace\hbox{$\hbox{M}_{\odot}$}$) remains in the superbubble’s cavity
$\sim 320$ kyr after the last PISN. This value is in agreement with the
conclusions raised in Paper I, that explored, at a higher resolution ($\sim
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
$\sim[10$–$65]\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ explode as core-collapse
supernovae. We have found that from the $\sim
30\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ of ejecta dust that were injected
into the simulation, $\sim 32$ per cent ($\sim
9.5\,\thinspace\hbox{$\hbox{M}_{\odot}$}$) survives after $175,000$ years from
the first core-collapse supernova explosion.

If we extrapolate the rate of net dust injection by core-collapse supernovae
until $\sim 21$ Myr of the cluster’s evolution (the time needed to explode a
$10\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ star in the BoOST low-metaillicity
stellar tracks), we can expect a net increase in the dust mass by $\sim
1200\,\thinspace\hbox{$\hbox{M}_{\odot}$}(M_{SC}/5.6\times
10^{5}\,\thinspace\hbox{$\hbox{M}_{\odot}$})$. Applying this estimate to the
gravitationally-lensed galaxy A2744_YD4 (at a redshift $z\sim 8.38$, Laporte
et al., 2017), with a total stellar mass $\sim 2\times
10^{9}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$, we obtain a net dust
contribution by supernovae $\sim 4.3\times
10^{6}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$. This value fits pretty well the
dust mass ($6\times 10^{6}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$) 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)3 (see the right panels in Figure 1), and the results agree within $15$ per
cent.

##  6\. 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\times
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 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 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 ($\gtrsim 2\times
10^{5}$ K), dust grains tend to be neutral (McKee et al., 1987).

##  7\. Concluding Remarks

By conducting 3–D hydrodynamical simulations, we have explored the early
($\sim 3.2$ Myr) evolution of a starburst-driven superbubble and the onset of
PISNe in a low-metallicity ($Z=0.02\,\thinspace\hbox{$\hbox{Z}_{\odot}$}$),
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., 2020) and the SynStars 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 $\sim 4.8\times 10^{-5}$ $\hbox{M}_{\odot}$ 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 $\sim 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 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 non-thermal, shock-induced
dust grain disruptive processes. Overall, $\sim 13$ and $\sim 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}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ yr-1 and $5.4\times
10^{-5}\,\thinspace\hbox{$\hbox{M}_{\odot}$}$ 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\sim 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.

##  8\. Acknowledgements

This study was supported by CONACYT-México research grant A1-S-28458. S.M.-G.
also acknowledges support by CONACYT through project n.482 of the “Programa
Investigadoras e Investigadores por México”. The authors thankfully
acknowledge the computer resources, technical expertise 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. & 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. The authors thank the anonymous Referee for a careful
reading and helpful suggestions which greatly improved the paper.
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).

## References

  * Alūzas et al. (2012) Alūzas, R., Pittard, J. M., Hartquist, T. W., Falle, S. A. E. G., & Langton, R. 2012, MNRAS, 425, 2212
  * Asano et al. (2013) Asano, R. S., Takeuchi, T. T., Hirashita, H., & Nozawa, T. 2013, MNRAS, 432, 637
  * Behrens et al. (2018) Behrens, C., Pallottini, A., Ferrara, A., Gallerani, S., & Vallini, L. 2018, MNRAS, 477, 552
  * Bisnovatyi-Kogan & Silich (1995) Bisnovatyi-Kogan, G. S., & Silich, S. A. 1995, Reviews of Modern Physics, 67, 661
  * Calura et al. (2014) Calura, F., Gilli, R., Vignali, C., et al. 2014, MNRAS, 438, 2765
  * Calura et al. (2008) Calura, F., Pipino, A., & Matteucci, F. 2008, A&A, 479, 669
  * Cardamone et al. (2009) Cardamone, C., Schawinski, K., Sarzi, M., et al. 2009, MNRAS, 399, 1191
  * Cherchneff (2010) Cherchneff, I. 2010, in Astronomical Society of the Pacific Conference Series, Vol. 425, Hot and Cool: Bridging Gaps in Massive Star Evolution, ed. C. Leitherer, P. D. Bennett, P. W. Morris, & J. T. Van Loon, 237
  * Colella & Woodward (1984) Colella, P., & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174
  * Consiglio et al. (2016) Consiglio, S. M., Turner, J. L., Beck, S., & Meier, D. S. 2016, ApJ Let, 833, L6
  * Crutcher (2012) Crutcher, R. M. 2012, ARA&A, 50, 29
  * Crutcher & Kemball (2019) Crutcher, R. M., & Kemball, A. J. 2019, Frontiers in Astronomy and Space Sciences, 6, 66
  * Dalgarno & McCray (1972) Dalgarno, A., & McCray, R. A. 1972, ARA&A, 10, 375
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium
  * Dwarkadas (2005) Dwarkadas, V. V. 2005, ApJ, 630, 892
  * Dwarkadas (2007) Dwarkadas, V. V. 2007, ApJ, 667, 226
  * Dwek (1998) Dwek, E. 1998, ApJ, 501, 643
  * Dwek & Cherchneff (2011) Dwek, E., & Cherchneff, I. 2011, ApJ, 727, 63
  * Elmegreen (1988) Elmegreen, B. G. 1988, Astrophysical Letters and Communications, 26, 207
  * Ferrara et al. (2016) Ferrara, A., Viti, S., & Ceccarelli, C. 2016, MNRAS, 463, L112
  * Ferriere et al. (1991) Ferriere, K. M., Mac Low, M.-M., & Zweibel, E. G. 1991, ApJ, 375, 239
  * Fowler & Hoyle (1964) Fowler, W. A., & Hoyle, F. 1964, ApJ, 9, 201
  * Franco et al. (1991) Franco, J., Tenorio-Tagle, G., Bodenheimer, P., & Rozyczka, M. 1991, PASP, 103, 803
  * Franeck et al. (2021) Franeck, A., Wünsch, R., Martínez-González, S., et al. 2021, ApJ, submitted
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJ, 131, 273
  * Gall & Hjorth (2018) Gall, C., & Hjorth, J. 2018, ApJ, 868, 62
  * Gall et al. (2011) Gall, C., Hjorth, J., & Andersen, A. C. 2011, A&A Rev., 19, 43
  * Haid et al. (2016) Haid, S., Walch, S., Naab, T., et al. 2016, MNRAS, 460, 2962
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Indebetouw et al. (2014) Indebetouw, R., Matsuura, M., Dwek, E., et al. 2014, ApJ Let, 782, L2
  * Jiménez et al. (2021) Jiménez, S., Tenorio-Tagle, G., & Silich, S. 2021, MNRAS, 505, 4669
  * Kirchschlager et al. (2020) Kirchschlager, F., Barlow, M. J., & Schmidt, F. D. 2020, ApJ, 893, 70
  * Koo & McKee (1992) Koo, B.-C., & McKee, C. F. 1992, ApJ, 388, 93
  * Koyama & Inutsuka (2002) Koyama, H., & Inutsuka, S.-i. 2002, ApJ Let, 564, L97
  * Kozyreva et al. (2014) Kozyreva, A., Blinnikov, S., Langer, N., & Yoon, S. C. 2014, A&A, 565, A70
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57
  * Lancaster et al. (2021a) Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021a, ApJ, 914, 89
  * Lancaster et al. (2021b) Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021b, ApJ, 914, 90
  * Laporte et al. (2017) Laporte, N., Ellis, R. S., Boone, F., et al. 2017, ApJ Let, 837, L21
  * Lau et al. (2021) Lau, R. M., Hankins, M. J., Kasliwal, M. M., et al. 2021, ApJ, 909, 113
  * Leśniewska & Michałowski (2019) Leśniewska, A., & Michałowski, M. J. 2019, A&A, 624, L13
  * Lucas et al. (2020) Lucas, W. E., Bonnell, I. A., & Dale, J. E. 2020, MNRAS, 493, 4700
  * Mac Low & McCray (1988) Mac Low, M.-M., & McCray, R. 1988, ApJ, 324, 776
  * Martínez-González et al. (2021) Martínez-González, S., Silich, S., & Tenorio-Tagle, G. 2021, MNRAS, 507, 1175
  * Martínez-González et al. (2016) Martínez-González, S., Tenorio-Tagle, G., & Silich, S. 2016, ApJ, 816, 39
  * Martínez-González et al. (2017) Martínez-González, S., Wünsch, R., & Palouš, J. 2017, ApJ, 843, 95
  * Martínez-González et al. (2018) Martínez-González, S., Wünsch, R., Palouš, J., et al. 2018, ApJ, 866, 40
  * Martínez-González et al. (2019) Martínez-González, S., Wünsch, R., Silich, S., et al. 2019, ApJ, 887, 198
  * Mattsson et al. (2015) Mattsson, L., Gomez, H. L., Andersen, A. C., & Matsuura, M. 2015, MNRAS, 449, 4079
  * McKee et al. (1987) McKee, C. F., Hollenbach, D. J., Seab, G. C., & Tielens, A. G. G. M. 1987, ApJ, 318, 674
  * McLeod et al. (2021) McLeod, A. F., Ali, A. A., Chevance, M., et al. 2021, MNRAS, arXiv:2109.08703
  * Michałowski et al. (2010) Michałowski, M. J., Watson, D., & Hjorth, J. 2010, ApJ, 712, 942
  * Micheva et al. (2017) Micheva, G., Oey, M. S., Jaskot, A. E., & James, B. L. 2017, ApJ, 845, 165
  * Nieuwenhuijzen & de Jager (1990) Nieuwenhuijzen, H., & de Jager, C. 1990, A&A, 231, 134
  * Nozawa et al. (2012) Nozawa, T., Kozasa, T., & Nomoto, K. 2012, ApJ Let, 756, L35
  * Nozawa et al. (2003) Nozawa, T., Kozasa, T., Umeda, H., Maeda, K., & Nomoto, K. 2003, ApJ, 598, 785
  * Palouš et al. (2013) Palouš, J., Wünsch, R., Martínez-González, S., et al. 2013, ApJ, 772, 128
  * Portegies Zwart et al. (2010) Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431
  * Priestley et al. (2021) Priestley, F. D., Chawner, H., Matsuura, M., et al. 2021, MNRAS, 500, 2543
  * Relaño et al. (2016) Relaño, M., Kennicutt, R., Lisenfeld, U., et al. 2016, A&A, 595, A43
  * Reynolds et al. (2019) Reynolds, D. R., Gardner, D. J., Balos, C. J., & Woodward, C. S. 2019, arXiv e-prints, arXiv:1909.12966
  * Rogers & Pittard (2013) Rogers, H., & Pittard, J. M. 2013, MNRAS, 431, 1337
  * Różyczka et al. (1993) Różyczka, M., Tenorio-Tagle, G., Franco, J., & Bodenheimer, P. 1993, MNRAS, 261, 674
  * Schure et al. (2009) Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751
  * Shull (1977) Shull, J. M. 1977, ApJ, 215, 805
  * Smith et al. (2017) Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217
  * Spitzer (1978) Spitzer, L. 1978, Physical processes in the interstellar medium, doi:10.1002/9783527617722
  * Svoboda et al. (2019) Svoboda, J., Douna, V., Orlitová, I., & Ehle, M. 2019, ApJ, 880, 144
  * Szécsi et al. (2015) Szécsi, D., Langer, N., Yoon, S.-C., et al. 2015, A&A, 581, A15
  * Szécsi & Wünsch (2019) Szécsi, D., & Wünsch, R. 2019, ApJ, 871, 20
  * Szécsi et al. (2020) Szécsi, D., Wünsch, R., Agrawal, P., & Langer, N. 2020, arXiv e-prints, arXiv:2004.08203
  * Taylor et al. (2018) Taylor, R., Wünsch, R., & Palouš, J. 2018, MNRAS, 479, 377
  * Tenorio-Tagle (2002) Tenorio-Tagle, G. 2002, in Revista Mexicana de Astronomia y Astrofisica Conference Series, Vol. 12, Revista Mexicana de Astronomia y Astrofisica Conference Series, ed. W. J. Henney, J. Franco, & M. Martos, 50–55
  * Tenorio-Tagle & Bodenheimer (1988) Tenorio-Tagle, G., & Bodenheimer, P. 1988, ARA&A, 26, 145
  * Tenorio-Tagle et al. (1990) Tenorio-Tagle, G., Bodenheimer, P., Franco, J., & Rozyczka, M. 1990, MNRAS, 244, 563
  * Tenorio-Tagle et al. (2006) Tenorio-Tagle, G., Muñoz-Tuñón, C., Pérez, E., Silich, S., & Telles, E. 2006, ApJ, 643, 186
  * Tenorio-Tagle et al. (2015) Tenorio-Tagle, G., Muñoz-Tuñón, C., Silich, S., & Cassisi, S. 2015, ApJ Let, 814, L8
  * Tenorio-Tagle et al. (1991) Tenorio-Tagle, G., Rozyczka, M., Franco, J., & Bodenheimer, P. 1991, MNRAS, 251, 318
  * Todini & Ferrara (2001) Todini, P., & Ferrara, A. 2001, MNRAS, 325, 726
  * Valentini & Brighenti (2015) Valentini, M., & Brighenti, F. 2015, MNRAS, 448, 1979
  * Valiante et al. (2009) Valiante, R., Schneider, R., Bianchi, S., & Andersen, A. C. 2009, MNRAS, 397, 1661
  * van Marle et al. (2015) van Marle, A. J., Meliani, Z., & Marcowith, A. 2015, A&A, 584, A49
  * Vink et al. (2000) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295
  * Vink et al. (2001) Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574
  * Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377
  * Wünsch et al. (2017) Wünsch, R., Palouš, J., Tenorio-Tagle, G., & Ehlerová, S. 2017, ApJ, 835, 60
  * Wünsch et al. (2018) Wünsch, R., Walch, S., Dinnbier, F., & Whitworth, A. 2018, MNRAS, 475, 3393
  * Zhukovska et al. (2008) Zhukovska, S., Gail, H. P., & Trieloff, M. 2008, A&A, 479, 453