††thanks: Theres authors contributed equally.††thanks: Theres authors
contributed equally.

# Observation of superradiance in a phase fluctuating dipolar Bose-Einstein
condensate

Bojeong Seo  Department of Physics, The Hong Kong University of Science and
Technology, Clear Water Bay, Kowloon, Hong Kong, China  Mingchen Huang
Department of Physics, The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Ziting Chen  Department of
Physics, The Hong Kong University of Science and Technology, Clear Water Bay,
Kowloon, Hong Kong, China  Mithilesh K. Parit  Department of Physics, The Hong
Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong
Kong, China  Yifei He  Department of Physics, The Hong Kong University of
Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China  Peng Chen
Department of Physics, The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Gyu-Boong Jo<EMAIL_ADDRESS>Department of Physics, The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  IAS Center for Quantum
Technologies, The Hong Kong University of Science and Technology, Kowloon,
Hong Kong, China  Microelectronics Thrust, The Hong Kong University of Science
and Technology (Guangzhou),
Guangzhou, China

###### Abstract

Despite the extensive study of matter-wave superradiance in a Bose-Einstein
condensate (BEC) using its unique coherence property, the controllability of
superradiant process has remained limited in the previous studies exploiting a
phase-coherent condensate with isotropic contact interactions. Here, we
combine tunable s-wave scattering with dipolar interactions in a BEC of 168Er
atoms wherein the asymmetry and threshold of superradiance are independently
controlled. By changing the s-wave scattering length near the Feshbach
resonance, we tune the superradiance threshold with increasing phase
fluctuations. In contrast to collective light scattering from a condensate
only with contact interactions, we observe an asymmetric superradiant peak in
a dipolar BEC by changing the direction of external magnetic field. This
results from the anisotropic excitation spectrum induced by the dipole-dipole
interaction. Our observation is expected to bring forth unprecedented
application of matter-wave optics leading to controlled emission of matter
wave.

Collective light scattering Dicke.1954je7 , a cooperative emission process
inducing directional scattered atoms, has been observed in various atomic
systems ranging from thermal atoms Gross.1982l7 ; Yoshikawa.2005 , degenerate
Bose gases Inouye.1999 ; Stenger.1999rlo ; Inouye.19990n6 ; Fallani.2005 ;
Deng.2010gfl ; Lu.2011je5 ; Kampel.2011 ; Lopes.2014 ; Dimitrova.2017 , free
fermions Wang.2010f8xf  to atoms coupled to the cavity mode Slama.2007 ;
Baumann.2010 ; Kessler.2014 . Among them, a Bose-Einstein condensate (BEC) has
served as a promising platform for exploring a matter-wave superradiant
process owing to its unique coherence property with Inouye.1999 ;
Stenger.1999rlo ; Inouye.19990n6 ; Fallani.2005 ; Deng.2010gfl ; Lu.2011je5 ;
Kampel.2011 ; Lopes.2014 ; Dimitrova.2017  and without external light fields
Clark.2017 ; Fu.2018i3f ; Wu.2019p5j ; Kim.2021 . When the external light
shines atoms in the condensate, collective scattering of light creates a
quasiparticle in the form of recoiling atoms that interfere with condensate
atoms at rest, leading to the generation of matter-wave grating that is
further enhanced by subsequent light scattering. So far, however, it has been
assumed in the previous studies that the condensate is phase-coherent and the
superradiant gain of the process entails only external parameters such as the
sample size and geometry. Furthermore, previous studies have remained in a
weakly interacting regime, and neither the effect of strong s-wave
interactions nor anisotropic dipolar interactions Lahaye.2009 ; Bismut.2012 ;
Wenzel.2018  has been addressed Gross.1982 .

Here, we investigate matter wave superradiance of a dipolar BEC in an
elongated trap, which combines tunable s-wave scattering with dipolar
interactions, such that the superradiance process is significantly modified by
phase fluctuations Petrov.20014x ; Dettmer.2001l97  and anisotropic dispersion
of the condensate Bismut.2012 ; Wenzel.2018 . Such dipolar effects in quantum
gases Chomaz.2022  have recently opened up a new regime where anisotropic
dipole-dipole interactions play a crucial role for realizing new phases of
matter, such as quantum droplets and supersolids Chomaz.2022 . In this work,
dipolar interactions result in anisotropic dispersion of dipolar superfluid
Bismut.2012 ; Wenzel.2018 , which allows us to control the asymmetry of
superradiant peaks by changing the direction of external magnetic field.
Furthermore, by controlling the s-wave scattering length near the Feshbach
resonance of 168Er, we tune the superradiance threshold with increasing phase
fluctuations Petrov.20014x ; Dettmer.2001l97 . Our work demonstrates how the
superradiance process can be controlled by the strength of s-wave scattering
and the dipolar interaction.

When pump photons are scattered from a condensate, long-lived quasi-particles
are excited leading to directional amplified Rayleigh scattering. The onset of
the superradiant process is accounted for with the number of recoiling atoms
$N$, following the rate equation $\dot{N}=(G-L)N$ where $G$ is the gain of the
superradiance and $L$ is the loss term. A dipolar BEC of 168Er offers a new
opportunity for the unprecedented control of superradiance. First, with
increasing s-wave scattering length, phase fluctuation shortens the coherence
length Petrov.20014x ; Dettmer.2001l97  and enhances the decay of the matter-
wave grating resulting in the increase of $L$. Secondly, in contrast to alkali
atoms with contact interactions, a dipolar BEC exhibits an anisotropic
excitation spectrum due to the dipole-dipole interaction Lahaye.2009 ;
Bismut.2012 ; Wenzel.2018  resulting in anisotropic superradiant gain $G$
being sensitive to the external magnetic field direction.

Figure 1:  Superradiance from a phase fluctuating dipolar BEC of 168Er (a) An
elongated condensate is exposed to the superradiance (SR) beam propagating
along the $+$y axis. The magnitude ($B_{0}$) and direction ($\hat{r}(\theta)$)
of external magnetic field is independently controlled leading to tunable
s-wave scattering and dipolar interactions. The pump beam is linearly
polarized along the z-direction while the magnetic field is always on the xy-
plane during the superradiance process. (b) Energy level diagram for relevant
transitions are shown. We use a narrow-line transition
($4f^{12}6s^{2}~{}({}^{3}H_{6})\rightarrow
4f^{12}({}^{3}H_{6})6s6p({}^{3}P_{1})$) with $\Gamma_{583}$=2$\pi\times$190
kHz for superradiance. Orange arrow represents the
$\ket{J=6;m=-6}\rightarrow\ket{J^{\prime}=7;m^{\prime}=-7}$ transition for a
pump beam red-detuned by $\Delta_{SR}$=7.767 MHz (6.6 MHz) at 0 mG (600 mG) of
external magnetic field. The external magnetic field changes between 300 mG
and 912 mG during the experiment. (c) Time-of-flight image of superradiant
dipolar BEC for $\theta=\pi/2$. After applying the pump beam, bi-directional
recoiled atom clouds are generated. Here, $k_{R}=2\pi/\lambda_{583}$ where
$\lambda_{583}$=583 nm.

Our experiments begin with 168Er BEC of approximately 1.8(2)$\times 10^{4}$
atoms in the $\ket{M_{J}=-6}$ state Seo.2020 ; qrl . The atoms are trapped in
the crossed optical dipole trap (ODT) consisting of two 1064 nm laser beams
propagating along x-direction and y-direction with the beam waist of
$w_{y-z}$=20 $\mu m$ and $w_{x-z}$=45 $\mu m$, respectively. To observe
matter-wave superraidance, we gradually change the trap geometry over 45 ms,
resulting in a quasi-1D trap with the trap frequency of
$(\omega_{x},\omega_{y},\omega_{z})=2\pi\times(37,445,443)~{}$Hz and the
maximum chemical potential of $\mu\simeq h\times$ 3.5 kHz at $T\simeq$ 140(10)
nK. During evaporative cooling, the magnetic field remains at 400 mG along the
y-direction. Before pulsing a pump light, we adjust the magnetic field to the
target value at variable angle over $t_{ramp}$ and switch off all ODT’s within
60 $\mu$s. At this stage, phase fluctuation is controlled by the magnitude of
bias magnetic field $B_{0}$ (near the Feshbach resonance) while the
anisotropic dispersion relation is set by the field direction $\theta$ as
described in Fig. 1.

To induce superradiant Rayleigh Scattering, an elongated condensate is
illuminated with a single off-resonant 583 nm laser beam for 100 $~{}\mu$s.
The pump beam, linearly polarized along the z-direction (see Fig. 1), is about
100 times larger than the size of sample - with the Thomas-Fermi radius of 15
$\mu$m - to maintain the same optical intensity along the condensate.
Additionally, the incident angle of the pump beam is perpendicular to the
elongated axis to circumvent the unintended asymmetric superradiance
Lu.2011je5 . Atomic momentum distribution after collective light scattering is
recorded after a 16 ms of time-of-flight expansion using absorption imaging
with a circularly polarized ($\sigma^{-}$) 401 nm broad-linewidth transition
($4f^{1}26s^{2}~{}{}^{3}H_{6}~{}\rightarrow~{}4f^{12}({}^{3}H_{6})6s6p({}^{1}P_{1})$)
along the z-direction. With data from the absorption imaging, we quantify the
number of recoiled atoms. The total number of recoiled atoms is denoted as
${N=N_{L}+N_{R}}$ where $N_{L}$ and $N_{R}$ are the number of atoms in left
and right recoiled atom clouds, respectively (see Fig. 1).

Figure 2: Control of superradiance threshold with the phase fluctuation (a)
Observation of portion of recoiled atoms versus the Rayleigh scattering rate
$R_{sc}$ for different phase fluctuation being proportional to
$T_{\phi}^{-1}$. To extract the threshold of Rayleigh scattering rate,
corresponding to the loss term $L$, the data curve is fitted to the equation
$N=Ae^{(R_{sc}-L)t_{pulse}}$ where $t_{pulse}$=100$\mu$s and $A$ is arbitrary
coefficient. Error bars indicate the standard error of 20 consecutive
measurements. (b) $L$ is plotted against the $T_{\phi}^{-1}$ and $a_{s}$. The
threshold of superradiance linearly increases with the amount of phase
fluctuation in the condensate. Inset illustrated the decrease of the $L$
against the increase of the coherence length $l_{\phi}$. Error bars indicate
90$\%$ confidence interval of loss for the fitting in (a). Dashed lines are
guides for eyes.  Figure 3:  Observation of asymmetric superradiant peaks. (a)
When the magnetic field brakes the left-right symmetry, the atomic excitations
propagating towards left and right superradiant peaks experience anisotropic
dispersion resulting in the generation of asymmetric matter-wave gratings
during the collective light scattering. By recording the number of recoiling
atoms, we obtain the distribution of asymmetry ratio ($R_{asym}$) from the
170-300 consecutive measurements. (b) Histograms display the distribution of
$R_{asym}$ at a different angle $\theta$ of the magnetic field. The atomic
momentum distribution after 16 ms time-of-flight shows left (right) asymmetric
recoiling atoms corresponding to the negative (positive) mean value of
$R_{asym}$. Yellow circles represent the mean values of the histograms. The
red solid line was a guideline based on sinusoidal function.

In the first set of experiments, we test for the control of the loss term $L$
by adjusting the phase coherence length of the condensate $l_{\phi}$. When the
condensate begins to interact with the superradiance pump photons, recoiling
atoms form a matter-wave grating with the lifetime $\tau_{c}$ being
proportional to $\tau_{c}\propto\frac{l_{\phi}}{v}$ where $v$ and $l_{\phi}$
are the velocity of the recoiled atoms and the characteristic coherence length
of the condesnate, respectively. Here, the loss term increases with
$l_{\phi}^{-1}$ as $L\propto 1/\tau_{c}\propto\frac{v}{l_{\phi}}$. For a
condensate without phase fluctuations, the coherence length can be considered
as the Thomas-Fermi radius $l$ of the condensate. With phase fluctuations,
however, the coherence length is given by $l_{\phi}\simeq l(T_{\phi}/T)$ where
$T$ is the temperature and $T_{\phi}$ is the characteristic temperature of
phase fluctuation defined as
$T_{\phi}=15(\hslash\omega_{z})^{2}N_{tot}/32\mu\propto a_{s}^{-2/5}$ with the
reduced Planck constant $\hbar$. Here, $\mu$ is the chemical potential and the
$a_{s}$ is the scattering interaction length Petrov.20014x ; Dettmer.2001l97 .
Therefore, the loss term can be controlled by the scattering length as
$L\propto T_{\phi}^{-1}\propto a_{s}^{2/5}$. In experiment, we exploit the
Feshbach resonance at $B_{res}\simeq$ 912 mG and control $a_{s}$ from
128$a_{0}$ to 279$a_{0}$ with $T_{\phi}$ ranging from 100 nK to 50 nK where
$a_{0}$ is the Bohr radius.

To control the scattering length $a_{s}$, the magnetic field is changed from
400 mG to the target value within $t_{ramp}$=2.7 ms right before the pump beam
is switched on, while the field direction is kept along the x-axis (i.e.
$\theta$=0). For different magnetic fields, we tune the Rayleigh scattering
rate upto 1.2 ms-1 based on the independent calibration (see supplementary
note), which effectively control the gain term $G$. Subsequently, we record
the portion of recoiled atoms against the change of Rayleigh scattering rate
while we adjust the characteristic temperature of phase fluctuation
($T_{\phi}^{-1}$) as shown in Fig. 2(a). Total recoiled atom number $N$
exponentially grows as the Rayleigh scattering rate (or $G$) increases. Then,
we extract the loss term $L$ by fitting the initial rise of recoiling atom
number with an exponential growth. We deduce the loss term from the measured
threshold value. Fig. 2 (b) shows the change in the loss ($L$) for variable
phase fluctuation characterized by $T_{\Phi}^{-1}$. The inset shows the decay
of loss against the increase of coherence length.

In contrast to a BEC only with isotropic contact interactions, a dipolar BEC
exhibits an anisotropic dispersion relation due to the dipole-dipole
interaction Lahaye.2009 ; Bismut.2012 ; Wenzel.2018 . The superradient
Rayleigh scattering is proportional to the square of the contrast of the
matter wave interference between the condensate and the recoiling atoms. In a
condensate without dipolar interactions, two matter-wave gratings,
corresponding to opposite end-fire modes, are formed in a symmetric manner
since the dispersion relation is isotropic with the constant sound velocity
being insensitive to the angle between the external magnetic field and the
propagation direction of atomic excitation. In other word, for the elongated
condensate the gain $G$ is the same for opposite directional emissions.

In a dipolar BEC, however, the number of recoiling atoms into opposite
directions is sensitive to the orientation of magnetic field. By quantifying
the contribution of dipole-dipole interaction $a_{dd}$ with respect to the
contact interactions as $\epsilon_{dd}=a_{dd}/a_{s}$ where
$a_{dd}=66.3~{}a_{0}$ in our system, the elementary Bogoliubov excitation
spectrum $\hbar\omega(k)$ for a uniform density $n$ is given as

$\hbar\omega(k)=\sqrt{E(k)(E(k)+2gn(1+\epsilon_{dd}(3\cos^{2}{\theta_{k}}-1))}$

where $E(k)=\frac{\hbar^{2}k^{2}}{2m}$ and $\theta_{k}$ is the angle between
the external field and the atomic excitation propagating along $\vec{k}$,
which reveals the anisotropy of the dispersion relation Lahaye.2009 . Here,
$g=4\pi\hbar^{2}a_{s}/m$ with atomic mass $m$. With the sound velocity given
as $c(\theta_{k})=c_{0}\sqrt{1+\epsilon_{dd}(3\cos^{2}{\theta_{k}}-1)}$ where
$c_{0}=\sqrt{2gn}$, excitations propagating along the magnetic field direction
(i.e. $\theta_{k}$=0) are stiffer with a larger sound velocity of
$c(0)=c_{0}\sqrt{1+2\epsilon_{dd}}$ than the one propagating in the
perpendicular direction (i.e. $\theta_{k}$=$\pi/2$) with
$c(\pi/2)=c_{0}\sqrt{1-\epsilon_{dd}}$ . In a confined system with
inhomogeneous density, however, roton softening at finite momentum $k$ is
expected to change the collective excitation, but it is expected this effect
is not significant in the current experiment with the small dipolar strength
of $\epsilon_{dd}\lesssim 0.5$, revealing an anisotropic dispersion
Wilson.2010 ; Ticknor.2011 .

Figure 4:  Insensitivity of asymmetric superradiance to the scattering length
(a) The asymmetry of superradiance peaks is insensitive to the scattering
length $a_{s}$. The error bar is the standard error and the solid line is a
guide for eyes. (b) The full distribution of asymmetry does not change with
$a_{s}$ within our experimental resolution. In this measurement, the magnetic
field is changed within $t_{ramp}$=0.5 ms with phase fluctuations being
negligible. (c) The distribution becomes narrower with the increasing Rayleigh
scattering rate. In (b,c), the error bar represents the 1 sigma confidence
interval of standard deviation of $R_{asym}$.

In our experiment, the phase-matching condition in a given trap geometry
allows scattered atoms to be only at 45∘ with respect to the x axis. When the
magnetic field is oriented at $\theta$=0∘ or 90∘, atomic excitations that
involves opposite matter-wave gratings occur with the same strength as
$\theta_{k}$ is $\pi/4$ for both cases. The effect of anisotropic dispersion
becomes prominent when magnetic field is aligned to one of excitation
directions (e.g. $\theta$=45∘ or 135∘). At $\theta$=45∘, for example, the
atomic excitation involving the right superradiant peak ($N_{R}$) undergoes
stiffer dispersion with a larger excitation resonance frequency (for
$\theta_{k}$=0) compared to the left peak ($N_{L}$). According to the f-sum
rule, the excitation involving the increased resonance frequency (i.e. right
peak) is relatively more suppressed than the left peak with a smaller
resonance frequency Nozieres.1990  (see Fig. 3(a)).

To elucidate the effect of magnetic field direction on superradiance, we keep
the Rayleigh scattering rate constant and therefore effectively the gain term
of superradiance the same. In our experiment, the end-fire mode always occurs
along the x direction regardless of magnetic field direction. In this case,
the spontaneous emission strength along the x direction is sensitive to the
quantization direction (i.e. external magnetic field) following the radiation
pattern (see supplementary note). Therefore, we calibrate the intensity of the
pump beam such that spontaneous emission along the longitudinal direction of
the atom cloud remains constant for the different quantization directions.
Furthermore, to avoid the shape change of condensate induced by the field
direction we switch on a pump beam immediately after rapidly rotating the
magnetic field in less than 0.5 ms.

In Fig. 3, we characterize a superradiant sample by recording an asymmetry of
superradiant peaks. For each condition, we obtain around 170-300 data and
plotted a histogram for an asymmetry ratio
($R_{asym}=\frac{N_{R}-N_{L}}{N_{R}+N_{L}}$) where $N_{R}$ ($N_{L}$)
represents the atom number of right (left) superradiant peak. When the left
(right) atom cloud is dominant, the mean value of the histogram has a negative
(positive) value (see yellow circles in Fig. 3).

It is interesting to notice that the anisotropic dispersion relation depends
only on the absolute magnitude of dipolar interaction $a_{dd}$ \- not the
relative strength $\epsilon_{dd}$ \- as shown in sound velocity
$c(\theta_{k})=\sqrt{gn/m+\frac{4\pi\hbar^{2}na_{dd}}{m^{2}}(3\cos^{2}{\theta_{k}}-1)}$.
To test for the insensitivity of asymmetric superradiance to the chemical
potential, we monitor the asymmetry of superradiant peaks at different
scattering lengths. In Fig. 4(a), the mean value of $R_{asym}$ does not change
with increasing $a_{s}$ confirming our understanding. In this case, the full
distribution of $R_{asym}$ is also insensitive to $a_{s}$ as shown by solid
green circles in Fig. 4(b).

Finally, for experimentally observed asymmetry $R_{asym}$ of around 100-150
consecutive measurements we find a finite width of the $R_{asym}$ distribution
even when the superradiance process is symmetric with the mean value of
$R_{asym}$ close to zero. We attribute this partly to the finite size of the
condensate with an aspect ratio of 12, which results in possible tilting of
the end fire mode by $\pm$5∘. This may introduce the uncertainty of the
direction of recoiling atoms, $\delta\theta_{k}\sim$0.055$\pi$. We also find
that the distribution becomes narrower with the increasing Rayleigh scattering
rate in Fig. 4(c), which suggests asymmetric superradiance occurs at the early
stage of the process.

In conclusion, we have explored a superradiant process in a dipolar BEC with
phase fluctuations highlighting versatile tools for controlling a matter-wave
superradiance process. The anisotropic dispersion allows us to control
asymmetric superradiant peaks by changing the direction of external magnetic
field. Furthermore, we demonstrate how to tune the superradiance threshold by
adjusting coherence property of BEC. In the future, it will be interesting to
investigate a superradiant process in an oblate dipolar condensate
Ticknor.2011 . It is then anticipated that the direction of recoiling atoms
may be controlled by the external magnetic field direction. Another
interesting future direction is to explore collective light scattering in
dipolar fermions Lu.2012 ; Aikawa.2014  or dipolar molecules Ni.2008 .

## Acknowledgements

We thank Hui Zhai for fruitful discussions. We acknowledge support from the
RGC through 16304918, 16308118, 16306119, 16302420, 16302821, 16306321
C6005-17G, C6009-20G and RFS2122-6S04. This works is also supported by the
Guangdong-Hong Kong Joint Laboratory and the Hari Harilela foundation.

## References

  * (1) R. H. Dicke, Physical Review 93, 99 (1954).
  * (2) M. Gross and S. Haroche, Physics Reports 93, 301 (1982).
  * (3) Y. Yoshikawa, Y. Torii and T. Kuga, Physical Review Letters 94, 083602 (2005).
  * (4) S. Inouye, A. Chikkatur, D. M. Stamper-Kurn, J. Stenger, D. Pritchard and W. Ketterle, Science 285, 571 (1999).
  * (5) J. Stenger, S. Inouye, D. Stamper-Kurn, A. Chikkatur, D. Pritchard and W. Ketterle, Applied Physics B 69, 347 (1999).
  * (6) S. Inouye, T. Pfau, S. Gupta, A. P. Chikkatur, A. Gorlitz, D. E. Pritchard and W. Ketterle, Nature 402, 641 (1999).
  * (7) L. Fallani, C. Fort, N. Piovella, M. Cola, F. S. Cataliotti, M. Inguscio and R. Bonifacio, Physical Review A 71, 033612 (2005).
  * (8) L. Deng, E. W. Hagley, Q. Cao, X. Wang, X. Luo, R. Wang, M. G. Payne, F. Yang, X. Zhou, X. Chen et al., Physical Review Letters 105, 220404 (2010).
  * (9) B. Lu, X. Zhou, T. Vogt, Z. Fang and X. Chen, Physical Review A 83, 033620 (2011).
  * (10) N. S. Kampel, A. Griesmaier, M. P. H. Steenstrup, F. Kaminski, E. S. Polzik and J. H. Muller, Physical review letters 108, 090401 (2011).
  * (11) R. Lopes, A. Imanaliev, M. Bonneau, J. Ruaudel, M. Cheneau, D. Boiron and C. I. Westbrook, Physical Review A 90, 013615 (2014).
  * (12) I. Dimitrova, W. Lunden, J. Amato-Grill, N. Jepsen, Y. Yu, M. Messer, T. Rigaldo, G. Puentes, D. Weld and W. Ketterle, Physical Review A 96, 051603 (2017).
  * (13) P. Wang, L. Deng, E. W. Hagley, Z. Fu, S. Chai and J. Zhang, Physical Review Letters 106, 210401 (2010).
  * (14) S. Slama, S. Bux, G. Krenz, C. Zimmermann and P. W. Courteille, Physical Review Letters 98, 053603 (2007).
  * (15) K. Baumann, C. Guerlin, F. Brennecke and T. Esslinger, Nature 464, 1301 (2010).
  * (16) H. Kessler, J. Klinder, M. Wolke and A. Hemmerich, Physical Review Letters 113, 070404 (2014).
  * (17) L. W. Clark, A. Gaj, L. Feng and C. Chin, Nature 551, 356 (2017).
  * (18) H. Fu, L. Feng, B. M. Anderson, L. W. Clark, J. Hu, J. W. Andrade, C. Chin and K. Levin, Physical Review Letters 121, 243001 (2018).
  * (19) Z. Wu and H. Zhai, Physical Review A 99, 063624 (2019).
  * (20) K. Kim, J. Hur, S. Huh, S. Choi and J.-y. Choi, Physical Review Letters 127, 043401 (2021).
  * (21) T. Lahaye, C. Menotti, L. Santos, M. Lewenstein and T. Pfau, Reports on Progress in Physics 72, 126401 (2009).
  * (22) G. Bismut, B. Laburthe-Tolra, E. Marechal, P. Pedri, O. Gorceix and L. Vernac, Physical Review Letters 109, 155302 (2012).
  * (23) M. Wenzel, F. Bottcher, J.-N. Schmidt, M. Eisenmann, T. Langen, T. Pfau and I. Ferrier-Barbut, Physical Review Letters 121, 030401 (2018).
  * (24) M. Gross and S. Haroche, Phys. Rep. 93, 301 (1982).
  * (25) D. S. Petrov, G. V. Shlyapnikov and J. T. M. Walraven, Physical Review Letters 87, 050404 (2001).
  * (26) S. Dettmer, D. Hellweg, P. Ryytty, J. J. Arlt, W. Ertmer, K. Sengstock, D. S. Petrov, G. V. Shlyapnikov, H. Kreutzmann, L. Santos et al., Physical Review Letters 87, 160406 (2001).
  * (27) L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev and T. Pfau, arXiv:2201.02672 , (2022).
  * (28) B. Seo, P. Chen, Z. Chen, W. Yuan, M. Huang, S. Du and G.-B. Jo, Physical Review A 102, 013319 (2020-07).
  * (29) B. Seo, Z. Chen, M. Huang, M. K. Parit, Y. He, P. Chen and G.-B. Jo, Journal of Korean Physical Society In press, (2022).
  * (30) R. M. Wilson, S. Ronen and J. L. Bohn, Physical Review Letters 104, 094501 (2010).
  * (31) C. Ticknor, R. M. Wilson and J. L. Bohn, Physical Review Letters 106, 065301 (2011).
  * (32) P. Nozieres, The Theory of Quantum Liquids (1990).
  * (33) M. Lu, N. Q. Burdick and B. L. Lev, Physical Review Letters 108, 215301 (2012).
  * (34) K. Aikawa, A. Frisch, M. Mark, S. Baier, R. Grimm and F. Ferlaino, Physical Review Letters 112, 010404 (2014).
  * (35) K. K. Ni, S. Ospelkaus, M. H. G. d. Miranda, A. Pe’er, B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne, D. S. Jin and J. Ye, Science 322, 231 (2008-10).

Supplementary Note

## Calibration of Rayleigh scattering rate

In our experiment, the Rayleigh scattering rate has been independently
calibrated at different magnetic field. Considering Lande g-factors for
$\ket{J^{\prime}=7;m^{\prime}=-7}$ and $\ket{J=6;m=-6}$ are $g_{e}=-11.71$
MHz/G and $g_{g}=-9.77$ MHz/G, respectively, we compensate the change of
detuning $\delta$ induced by the Zeeman shift by adjusting the intensity $I$
of the pump beam following the relation

$R_{scatt}=\frac{\Gamma_{583}}{2}\frac{\Omega_{R}^{2}/2}{\delta^{2}+\Omega_{R}^{2}/2+\Gamma_{583}^{2}/4}$

with the Rabi frequency $\Omega_{R}=\Gamma_{583}\sqrt{\frac{I}{2I_{sat}}}$. As
an example, we select 300 mG and 400 mG in Fig. S1 where the phase fluctuation
can be avoided with the small scattering length, and monitor the number of
rocoiling atoms as a function of Rayleigh scattering rate.

Figure S1:  Rayleigh scattering rate calibration. The Rayleigh scattering rate
is independently controlled at different magnetic fields. For example, we
calibrate the superradiance at 300 mG and 400 mG, where $a_{s}$ are almost the
same. The dashed lines are a guideline. The error bar is the standard
deviation.

## Calibration of effective pumping intensity

Radiation from an atom has an angular dependence on a quantization axis.
Therefore, the Rayleigh scattering rate changes while we rotate the external
magnetic field during experiments. Additionally, only radiation along the
longitudinal direction of the atom cloud is amplified by the gain medium and
contributes the superradiance. Therefore, we calibrate the intensity of the
pump beam such that spontaneous emission along the longitudinal direction of
the atom cloud remains constant for the different quantization directions.
Since the superradiance process involves $\sigma$-transitions, the angular
distribution of emission spectrum $I_{eff}$ is given as
$I_{eff}(\alpha)=I(1+\cos^{2}{\alpha})/2$ where $I$ is the intensity of the
pump beam and $\alpha$ is an angle of the emission direction with respect to
the quantization axis.

To observe asymmetric superradiance induced by the dipole-dipole interaction,
we keep the constant scattering rate into the end- fire mode (i.e. $I_{eff}$)
by using variable pump intensity $I(\theta)$ as
$I(\theta)=2I_{eff}/(1+\cos^{2}{\theta})$ where the magnetic field has an
angle $\theta$ with respect to the x axis. Fig. S2 shows the angular
dependence of the required pump beam intensity at different angle of magnetic
field.

Figure S2:  Calibration of effective pumping intensity When collective lights
scattering occurs, we keep the emission rate along the x-direction constant as
the magnetic field rotates. The required pump beam intensity scales as
$1/(1+\cos^{2}{\theta})$ as described by a solid curve. Ellipses and dumbbells
in insets are atom condensates and radiation patterns, respectively at 0∘,
45∘, and 90∘. Here, $\theta$ is an angle of magnetic field (green arrow) with
respect to the longitudinal direction of the atom cloud.

## Rotation of magnetic field

Bias magnetic field was generated by three pairs of Helmholtz coils
perpendicular to each other in Fig. 1. To rotate the direction of the magnetic
field, we tuned currents in coils within 500 us. Due to the dipole-dipole
interaction, the shape of dipolar condensates could be changed when the
direction of the magnetic field was tuned. Moreover, the cloud shape change
affects the gain of the superradiance. Hence, we prevented the cloud shape
change by sending a pump beam immediately after rapidly rotating the magnetic
field in less than 0.5 ms.


# Dust Grain Growth at High Redshift: Starburst-driven CMB-Dark Supershells

Sergio Martínez-González ${}^{1},$ Sergiy Silich ${}^{2},$ Guillermo Tenorio-
Tagle${}^{2},$
1CONACYT-Instituto Nacional de Astrofísica, Óptica y Electrónica, AP 51, 72000
Puebla, México
2Instituto Nacional de Astrofísica Óptica y Electrónica, AP 51, 72000 Puebla,
México
E-mail<EMAIL_ADDRESS>
(Accepted July 27th, 2021.)

###### Abstract

We present a novel scenario for the growth of dust grains in galaxies at high-
redshift ($z\sim 6$). In our model, the mechanical feedback from massive star
clusters evolving within high-density pre-enriched media allows to pile-up a
large amount of matter into massive supershells. If the gas metallicity
($\geq\thinspace\hbox{$\hbox{Z}_{\odot}$}$), number density ($\geq 10^{6}$
cm-3) and dust-to-gas mass ratio ($\sim 1/150\times Z$) within the supershell
are sufficiently large, such supershells may become optically thick to the
starlight emerging from their host star clusters and even to radiation from
the Cosmic Microwave Background (CMB). Based on semi-analytic models, we argue
that this mechanism, occurring in the case of massive ($\geq
10^{7}\thinspace\hbox{$\hbox{M}_{\odot}$}$) molecular clouds hosting $\geq
10^{6}\thinspace\hbox{$\hbox{M}_{\odot}$}$ star clusters, allows a large mass
of gas and dust to acquire a temperature below that of the CMB, whereupon dust
grain growth may occur with ease. In galaxies with total stellar mass $M_{*}$,
grain growth within supershells may increase the dust mass by $\sim
10^{6}\thinspace\hbox{$\hbox{M}_{\odot}$}(M_{*}/10^{8}\thinspace\hbox{$\hbox{M}_{\odot}$})$.

###### keywords:

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

††pubyear: 2021

##  1 Introduction

Violent star formation within the first $\sim 10^{9}$ years of cosmic time
($z\gtrsim 6$) induced a super-solar metal enrichment as observed in quasar
host galaxies (e.g. Jiang et al., 2007; Calura et al., 2014). Furthermore, the
detection of massive amounts of dust in the gravitationally-lensed galaxy
A2744_YD4 ($z\sim 8.38$) reveals that an early dust enrichment occurs during
the epoch of cosmic reionisation (Laporte et al., 2017), when the universe was
$\sim 200$ Myr old. It has been argued that galaxies at high-redshifts must
experience a rapid ($\sim$ a few million years long) transition between
relatively dust-free and dust-rich (Mattsson et al., 2015), and that star
formation may have proceeded with a top-heavy initial mass function (Gall et
al., 2011; Dwek & Cherchneff, 2011). Several possible dust producers at high
redshift have been discussed in the literature. For example, the cold
envelopes of the most massive AGB stars have been regarded as possible dust
producers at high-redshift (e.g. Valiante et al., 2009; Leśniewska &
Michałowski, 2019), particularly in the case of a top-heavy initial mass
function (Chiosi et al., 1998), while Todini & Ferrara (2001, and references
therein); Indebetouw et al. (2014, and references therein) have claimed that a
large fraction of the early dust enrichment follows from the very fast and
highly efficient dust grain condensation in supernova ejecta.

The relative importance of dust formation in supernova ejecta, dust
destruction by supernova-driven shocks and dust grain growth (by accretion of
gas-phase species) in dense and cold molecular clouds (MCs) is still matter of
intense debate. For instance, Calura et al. (2008); Calura et al. (2014),
confronted the dust destruction rate obtained by McKee (1989) for SNe
exploding in a three-phase interstellar medium to the dust production rates
derived from galactic chemical evolutionary models (see also Dwek, 1998;
Zhukovska et al., 2008; Asano et al., 2013). They asserted that dust grain
growth is a necessary mechanism to account for the presence of dust in
galaxies, and to effectively counteract dust destruction in supernova-driven
shocks. The processes of grain coagulation and grain growth are tightly
related, but the former has been shown to be unimportant in calculating the
evolution of the total dust mass in galaxies (Hirashita, 2012). More recently,
e.g. Gall & Hjorth (2018), have provided supporting evidence in that the dust
content in local and high-redshift galaxies is consistent with an efficient
dust production by supernovae and dust reformation shortly after destruction.

Notwithstanding, Ferrara et al. (2016) has noted that the process of grain
growth must overcome several complications under the conditions prevailing in
galaxies at high redshift. First, at $z\gtrsim 6$, the Cosmic Microwave
Background (CMB) sets the minimum possible grain temperature as
$T_{cmb}\gtrsim 20$ K (da Cunha et al., 2013), this translates into warmer
grain surfaces and reduced accretion rates. Second, after the host molecular
clouds are disrupted (typically in $\sim 10$ Myr), any material accreted onto
the grain surfaces will almost immediately be photo-desorbed as the grains are
exposed to the radiative feedback. A third complication arises from the
Coulomb repulsive forces between positively-charged ions and grains in the
presence of a strong UV radiation field. Given the aforementioned problems,
their view favours the notion that the large quantities of dust observed at
high redshift must have been readily condensed in supernova ejecta. In this
paper we will reexamine these objections in the context of starburst-driven
superbubbles.

Figure 1: Schematic representation of the structure of a starburst-driven
superbubble, and the central star cluster (not to scale). Each region of the
superbubble is indicated: the region filled with shock-heated gas, the
dense/warm ionised supershell, the very dense/cold neutral supershell, and the
unshocked turbulent molecular cloud. The bath of CMB photons reaching the
neutral supershell is represented by wavy arrows.

In the “standard bubble model”, the stellar winds of massive stars,
interacting with the surrounding ISM, create a four-distinct-zone structure,
namely: (i) the central hot and tenuous free wind region, (ii) a (reverse-)
shocked wind region separated by a contact discontinuity from (iii) a massive,
rapidly-cooling, supershell of swept-up interstellar matter that is able to
absorb a large part of the impinging ionising and non-ionising radiation, and
(iv) the ambient medium, an MC in our case (see Weaver et al., 1977; Mac Low &
McCray, 1988; Koo & McKee, 1992; Bisnovatyi-Kogan & Silich, 1995).

As the gas in the supershell recombines efficiently and dust absorbs large
amounts of radiation, it may eventually occur that the ionisation front
becomes trapped within the supershell and thus a fifth zone is formed (see
Figure 1): (v) an outer dusty neutral layer in the supershell. For example,
for a supershell driven by the mechanical feedback of a
$10^{6}\thinspace\hbox{$\hbox{M}_{\odot}$}$ stellar cluster, embedded into a
medium with constant density $\sim 10^{3}$ cm-3, the ionisation front becomes
trapped in $\sim 1.5$ Myr (Martínez-González et al., 2014), and earlier at
higher densities. Following Weaver et al. (1977), we will not distinguish
between regions (i) and (ii), because the radius of the eventually-formed
reverse shock is usually much smaller than the radius of the forward shock.
Thus, in the following we will refer to them as the shock-heated gas region.

Moreover, radiative cooling induced by dust grains (Whitworth, 2016) and
hydrogen-deuteride (HD) molecules, with the latter efficiently formed in post-
shocked gas, may contribute to decrease the neutral supershell temperature to
the floor set by the CMB (Johnson & Bromm, 2006).

Furthermore, if the dust column density and the average grain cross section
within the neutral supershell are sufficiently large, there is an opportunity
to trap all incoming CMB radiation before reaching the supershell’s ionised
part. Consequently, some part of the neutral supershell may become self-
shielded from the inner starlight and from the outer CMB radiation.
Interestingly, self-shielding and cooling via adiabatic expansion have been
invoked to explain that the outflowing matter in the Boomerang Nebula is
colder than the current CMB temperature (Sahai & Nyman, 1997; Sahai et al.,
2017). Our investigation tests the feasibility of dust grain growth within the
neutral parts of very dense ($\geq 10^{6}$ cm-3), dusty (with dust-to-gas mass
ratio $\sim 1/150\times Z$), and heavily pre-enriched
($\geq\thinspace\hbox{$\hbox{Z}_{\odot}$}$), starburst-driven supershells.

The paper is organised as follows: Section 2 describes the superbubble model,
including the equations governing the supershell’s ionised and neutral zones,
and the main parameters and relevant time-scales. In Section 3 we change
several parameters to explore the evolution of supershells under different
ambient conditions and star cluster masses. In Section 4 we study the
feasibility of dust grain growth in supershells at high-redshift. Finally, in
Section 5 we summarise our conclusions.

##  2 The large-scale evolution

###  2.1 The Host Molecular Cloud

We consider a young massive stellar cluster located at the centre of a pre-
enriched molecular cloud supported by turbulent pressure (Elmegreen & Efremov,
1997; Johnson et al., 2015; Calura et al., 2015; Elmegreen, 2017, see Appendix
A). We assume that the molecular cloud is formed out of material with super-
solar metallicity ($Z_{ISM}\geq\thinspace\hbox{$\hbox{Z}_{\odot}$}$). The
turbulent molecular cloud’s assumed gas density distribution is a power-law of
the form (e.g. Lee et al., 2015; Raskutti et al., 2017)

$\displaystyle\rho=\rho_{c}\left(\frac{r}{R_{c}}\right)^{-\omega}\,~{},\,\,\,\,\,\mbox{for
}r\geq R_{c}~{},$ (1)

where $r$ is the radial distance, $R_{c}$ is a characteristic length scale or
“core radius”, and $\rho_{c}$ is the gas density at $r=R_{c}$. Within the
central zone, we assume that the gas density distribution is homogeneous,

$\displaystyle\rho=\rho_{c}\,~{},\,\,\,\,\,\mbox{for }r<R_{c}~{}.$ (2)

The molecular cloud is truncated at a radius $R_{cut}$$\gg$$R_{c}$.

####  2.1.1 Superbubble Evolution

We assume that the kinetic energy injected by massive stars is thermalized in
neighboring stellar wind collisions. Once the thermalized gas pressure
overcomes the turbulent pressure in the ambient molecular cloud, the shock-
heated gas begins to expand. In the spherically-symmetric case, this occurs at
$t_{0}=4\pi R^{3}_{c}P_{t}(R_{c})/(3(\gamma-1)L_{SC}$, where $L_{SC}$ is the
star cluster mechanical luminosity, $P_{t}(R_{c})$ is the molecular gas
turbulent pressure at the star cluster edge (see Appendix A) and $\gamma=5/3$
is the ratio of specific heats. The dynamics of the expanding bubble is
determined by the set of mass, momentum and energy conservation equations
(Weaver et al., 1977; Mac Low & McCray, 1988; Bisnovatyi-Kogan & Silich,
1995). In the case of a homogeneous ambient gas density distribution these
equations have a well-known power-law solution. This solution, however, is
asymptotic as it does not account for the initial conditions (at $t=0$ years,
the shell radius is zero while the shell velocity is infinite) and also it
does not take into consideration the thermalized gas radiative cooling. We add
to the Weaver’s et al. equations the re-inserted gas radiative cooling and
adopt that the ambient gas density distribution is determined by equation (1).
The conservation equations that determine the starburst-driven bubble
evolution then are

$\displaystyle
M_{sh}(r)=\frac{4\pi\rho_{c}R^{3}_{c}}{3-\omega}\left[\left(\frac{r}{R_{c}}\right)^{(3-\omega)}-1\right],$
(3) $\displaystyle\frac{{\rm d}u}{{\rm d}t}=\frac{4\pi
r^{2}}{M_{sh}(r)}\left[P_{b}-\rho_{c}u^{2}\left(\frac{R_{c}}{r}\right)^{\omega}\right],$
(4) $\displaystyle\frac{{\rm d}E_{b}}{{\rm d}t}=L_{SC}-4\pi
r^{2}P_{b}u-\frac{3(\dot{M}t)^{2}\Lambda(T_{b},Z_{ISM})}{4\pi
r^{3}\eta^{2}_{ion}},$ (5) $\displaystyle\frac{{\rm d}r}{{\rm d}t}=u,$ (6)

where $M_{sh}(r)$ is the supershell’s gas mass, the shock-heated gas thermal
pressure $P_{b}=3(\gamma-1)E_{b}/(4\pi r^{3})$, $\mu_{ion}=14/11m_{H}$ is the
mean mass per ion in the mass-loaded, shock-heated plasma with 1 helium atom
per each 10 hydrogen atoms, $m_{H}$ is the mass of the hydrogen atom,
$\Lambda(T_{b},Z_{ISM})$ is the Raymond et al. (1976) cooling function and
$\dot{M}=2L_{SC}/V_{\infty}^{2}$ is the mass input rate within the star
cluster volume that may include the re-inserted gas mass loading by the
residual gas.

The inside bubble gas mass density and temperature are

$\displaystyle\rho_{b}=\frac{3\dot{M}t}{4\pi r^{3}},$ (7) $\displaystyle
T_{b}=\frac{\mu_{a}P_{b}}{k_{B}\rho_{b}},$ (8)

where $\mu_{a}=14/23m_{H}$ is the mean mass per particle in the completely
ionised gas with 1 helium atom per each 10 hydrogen atoms and $k_{B}$ is the
Boltzmann constant. We solve equations (3) - (6) numerically upon the
assumptions that, at $t=t_{0}$, the shell radius is slightly larger than the
star cluster radius, $r=\alpha R_{c}$, the shell velocity $u=0$ km s-1 and
$\alpha=1.05$.

In the case of a homogeneous ambient gas density distribution and negligible
thermalized gas cooling, the numeric solution rapidly approaches the classic
power-law solution.

If the residual gas mass is evenly incorporated into the outflow during a time
$t_{ml}$, then the mass-loading factor is

$\displaystyle\eta_{ml}=\frac{M_{MC}(r<R_{c})}{\dot{M}_{w}t_{ml}}~{},$ (9)

where $M_{MC}(r<R_{c})$ is

$\displaystyle M_{MC}(r<R_{c})=\frac{4}{3}\pi\rho_{c}R_{c}^{3}~{}.$ (10)

As the outflow is mass-loaded, the adiabatic terminal speed is changed to

$\displaystyle V_{\eta,\infty}=V_{\infty}(1+\eta_{ml})^{-1/2}~{}.$ (11)

###  2.2 The Micro-Physics and the Shell Inner Structure

The supershell’s inner structure is calculated by integrating the following
set of equations, starting from the supershell’s inner edge, $R_{b}$, (Draine,
2011; Martínez-González et al., 2014)

$\displaystyle\frac{{\rm d}}{{\rm
d}r}\left(\frac{\mu_{ion}}{\mu_{a}}nk_{B}T\right)$ $\displaystyle=$
$\displaystyle n\sigma_{d}\frac{\left[L_{n}e^{-\tau}+L_{i}\phi\right]}{4\pi
r^{2}c}+\frac{n^{2}\beta_{2}L_{i}}{Q_{0}c}~{},$ (12)
$\displaystyle\frac{Q_{0}}{4\pi r^{2}}\frac{{\rm d}\phi}{{\rm d}r}$
$\displaystyle=$
$\displaystyle-\beta_{2}n^{2}-\frac{n\sigma_{d}Q_{0}\phi}{4\pi r^{2}}~{},$
(13) $\displaystyle\frac{{\rm d}\tau}{{\rm d}r}$ $\displaystyle=$
$\displaystyle n\sigma_{d}~{}.$ (14)

In the above equations $L_{i}$, $L_{n}$ and $Q_{0}$ are the cluster’s ionising
and non-ionising luminosity, and the number of ionising photons that emerge
from the cluster per unit time; $n$ and $T$ are the supershell’s gas number
density and temperature, respectively; $\sigma_{d}$ is the grain absorption
cross section per particle, $\phi$ is the fraction of the ionising power that
gets to a surface with radius $r$, while $\tau$ is the absorption optical
depth. Finally, $c$ and $\beta_{2}=2.59\times 10^{-13}$ cm3 s-1, are the speed
of light, and the recombination coefficient to the excited states of H
(Osterbrock, 1989), respectively. The set of equations (12-14) allows to
calculate the thickness of the supershell’s ionised section.

The grain cross section per particle, $\sigma_{d}$, averaged by the Planck
function, $B_{\lambda}$, at blackbody temperature $T_{\lambda}$ and wavelength
$\lambda$, and grain size distribution is calculated as (e.g. Ferrara et al.,
2017)

Figure 2: The grain absorption cross section per particle averaged by a log-
normal grain size distribution as a function of the blackbody temperature for
carbonaceous grains and metallicity $2.5\thinspace\hbox{$\hbox{Z}_{\odot}$}$.
The blue solid and red dashed lines depict $\sigma_{d}$ for a log-normal grain
size distribution (see equation 17 below) with $a_{0}=0.1\mu$m, and
$\delta=0.01$ and $\delta=0.05$, respectively. The vertical line indicates the
value of $T_{cmb}(z=6)=19.1$ K.

$\displaystyle\sigma_{d}$ $\displaystyle=$
$\displaystyle\int_{a_{min}}^{a_{max}}\int_{0}^{\infty}\frac{\pi^{2}a^{2}Q_{abs}(a,\lambda)B_{\lambda}(T_{\lambda})}{\sigma_{SB}T_{\lambda}^{4}}\frac{\partial
n}{\partial a}\,{\rm d}{a}\,{\rm d}{\lambda}~{},$ (15)

where $a$ denotes the grain radius, $Q_{abs}(a,\lambda)$ is the dust
absorption efficiency (Draine, 2003; Compiègne et al., 2011) and $\sigma_{SB}$
is the Stefan-Boltzmann constant. $\sigma_{d}$ in equations (12-14), is
evaluated from equation (15) at the blackbody temperature of the integrated
stellar spectrum (see Figure 2).

At the supershell’s inner edge, the boundary conditions are that $\tau$ is
equal to zero and $\phi$ equal to unity. Given sufficient density and dust
content, the ionisation front can be trapped within the supershell. This
occurs when $\phi$ goes to zero. Thus a neutral skin is formed in the
supershell’s outer part (Martínez-González et al., 2014; Rahner et al., 2017),
whose thickness is determined by the thermal pressure and the difference
between the swept-up mass (see equation (26)) and the mass of ionised gas.

###  2.3 Neutral Shell

The optical depth to CMB photons is obtained by considering that
$N_{H}(r^{\prime})$ is the gas column density measured from the supershell’s
outer edge inwards:

$\displaystyle\tau_{cmb}(r^{\prime})=\sigma_{d}N_{H}(r^{\prime})~{},$ (16)

where $\sigma_{d}$ is evaluated at $T_{\lambda}=T_{cmb}(z)$.

We have assumed a log-normal grain size distribution of the form (e.g. Hensley
& Draine, 2017)

$\displaystyle\displaystyle\frac{\partial n}{\partial a}\sim
a^{-1}\exp\left\\{-\frac{1}{2\delta^{2}}\left[\displaystyle\ln\left(\frac{a}{a_{0}}\right)\right]^{2}\right\\}~{},$
(17)

where $a_{0}$ and $\delta$ set the peak grain size and width of the size
distribution, and $a_{min}$ and $a_{max}$ are the lower and upper size cut-
offs, respectively. The grain size distribution is normalised to the dust-to-
gas mass ratio $\mathcal{D}$, which is assumed to scale with $Z_{ISM}$, as

$\displaystyle\mathcal{D}=\mathcal{\thinspace\hbox{$\hbox{D}_{\odot}$}}Z_{ISM}=\frac{4\pi}{3\mu_{mol}}\rho_{gr}\int_{a_{min}}^{a_{max}}a^{3}\frac{\partial
n}{\partial a}{\rm d}{a},$ (18)

where $\mu_{mol}=2.33m_{H}$ is the mean mass of molecular gas, $\rho_{gr}$ is
the grains’ bulk density (equal to $2.26$ g cm-3 for carbonaceous grains),
$Z_{ISM}$ is measured in solar units and $\hbox{D}_{\odot}$ is set to $1/150$
as in the solar vicinity.

Figure 3: The supershell’s gas mass density structure for model A3500 at $0.5$
Myr. The supershell’s ionised and neutral parts, and the cloud’s unshocked
part are indicated. The colour scale depicts the fraction of the CMB radiation
that is trapped as it paths through the supershell.

####  2.3.1 Hydrogen-Deuteride Molecules (HD) and Dust Cooling

Hydrogen-deuteride molecules, efficiently formed in post-shocked gas (Vasiliev
& Shchekinov, 2005), allow the supershell’s neutral part to cool down as a
result of radiation emitted from their ground state rotational transition. As
long as the energy levels of HD molecules are populated by CMB photons, the
gas cannot cool below $T_{cmb}(z)$. The time-scale to reach the CMB
temperature via HD cooling is approximated by (Johnson & Bromm, 2006)

$\displaystyle
t_{cmb}\simeq\frac{1}{2A_{10}X_{HD}}\left(\frac{k_{B}T_{cmb}(z)}{h\nu_{10}}\right)^{2}\exp\left(\frac{h\nu_{10}}{k_{B}T_{cmb}(z)}\right),$
(19)

where $A_{10}=5\times 10^{8}$ s-1 is the HD Einstein coefficient for
spontaneous emission, $X_{HD}$ is the fractional abundance of HD,
$\nu_{10}\approx 2.7$ THz and $h$ is the Planck constant. At $z=6$,
$t_{cmb}\approx 6.21X_{HD}^{-1}$ years.

At densities $\gtrsim 10^{6}$ cm-3, however, dust grains are generally the
dominant coolants (Klessen & Glover, 2016) and the time required to cool the
neutral supershell to $T_{cmb}(z)$ can be significantly reduced (Meece et al.,
2014). The gas temperature evolution as a consequence of cooling via gas-grain
collisions is given by

$\displaystyle\frac{{\rm d}T}{{\rm d}t}=\frac{2\Lambda_{gas\rightarrow
dust}}{5nk_{B}}~{},$ (20)

where $\Lambda_{gas\rightarrow dust}\sim[1.2\times 10^{2}$ cm2 g-1]
$(\mu_{mol}n)^{2}c_{s}^{3}$ (Whitworth, 2016), and
$c_{s}=(k_{B}T\mu_{mol}^{-1})^{1/2}$ is the isothermal sound speed. The above
equation can be integrated to obtain the time required to cool the neutral
supershell down to $T_{cmb}$. The strong radiative cooling in the supershell
leads to its rapid collapse, making it extremely thin (see Figure 3).

We calculate the supershell’s inner structure using time-dependent parameters,
$Q_{0}$, $L_{i}$ and $L_{n}$, obtained from the Starburst99 synthesis
evolutionary code (Leitherer et al., 1999). We have considered in all cases a
Kroupa initial mass function sampled in the interval ($0.1-100$)
$\hbox{M}_{\odot}$ and Padova evolutionary models including AGB stars.

##  3 Supershell Evolution Models

Model | $M_{SC}$ | $L_{SC}$ | $V_{\eta,\infty}$ | $R_{c}$ | $R_{cut}$ | $n_{c}$ | $\eta_{ml}$ | $t_{ml}$ | $M_{MC}$ | $\epsilon$ | $\omega$
---|---|---|---|---|---|---|---|---|---|---|---
|  $\hbox{M}_{\odot}$ | 1041 erg s-1 | km s-1 | pc | pc | cm-3 | $-$ | Myr |  $\hbox{M}_{\odot}$ | $-$ | $-$
A3500 | $2\times 10^{6}$ | $3.12$ | $3500$ | $5.0$ | $100$ | $3.0\times 10^{4}$ | $0.0$ | $\infty$ | $5.3\times 10^{7}$ | $3.7\%$ | $2.0$
A1750 | $2\times 10^{6}$ | $3.12$ | $1750$ | $5.0$ | $100$ | $3.0\times 10^{4}$ | $3.0$ | $3.8$ | $5.3\times 10^{7}$ | $3.7\%$ | $2.0$
B3500 | $2\times 10^{6}$ | $3.12$ | $3500$ | $5.0$ | $48$ | $3.0\times 10^{4}$ | $0.0$ | $\infty$ | $5.3\times 10^{7}$ | $3.7\%$ | $1.5$
B1750 | $2\times 10^{6}$ | $3.12$ | $1750$ | $5.0$ | $48$ | $3.0\times 10^{4}$ | $3.0$ | $3.8$ | $5.3\times 10^{7}$ | $3.7\%$ | $1.5$
C3500 | $2\times 10^{6}$ | $3.12$ | $3500$ | $5.0$ | $100$ | $3.0\times 10^{4}$ | $0.0$ | $\infty$ | $2.0\times 10^{6}$ | $10\%$ | $2.5$
C1750 | $2\times 10^{6}$ | $3.12$ | $1750$ | $5.0$ | $100$ | $3.0\times 10^{4}$ | $3.0$ | $3.8$ | $2.0\times 10^{6}$ | $10\%$ | $2.5$
D3500 | $1\times 10^{6}$ | $1.56$ | $3500$ | $3.0$ | $100$ | $5.0\times 10^{4}$ | $0.0$ | $\infty$ | $3.2\times 10^{6}$ | $3.1\%$ | $2.0$
D1750 | $1\times 10^{6}$ | $1.56$ | $1750$ | $3.0$ | $100$ | $5.0\times 10^{4}$ | $3.0$ | $2.7$ | $3.2\times 10^{6}$ | $3.1\%$ | $2.0$
Table 1: The Table presents the main parameters of our six models: the star
cluster mass, the mechanical luminosity, the adiabatic terminal speed, the
cloud’s core and cut-off radii, the cloud’s gas number density at the core
radius, the mass-loading factor, the mass-loading time-scale, the molecular
cloud’s gas mass, the star formation efficiency, and the cloud’s gas density
distribution power-law index.

Our reference case A3500 (see Table 1), considers a star cluster with mass
$M_{SC}=2\times 10^{6}$ $\hbox{M}_{\odot}$, located within the central
($r<R_{c}=5$ pc), dense ($n_{c}=\rho_{c}\mu_{mol}^{-1}=3.05\times 10^{4}$
cm-3) zone of a molecular cloud (see equations 1 and 2), and a total mass
$M_{MC}=5.3\times 10^{7}\thinspace\hbox{$\hbox{M}_{\odot}$}$. This assumes
that the molecular cloud follows a power-law gas density distribution with
$\omega=2$ at $r\geq R_{c}$. This value is compatible with an isothermal,
self-gravitating molecular cloud (e.g. Rahner et al., 2019, and references
therein). The molecular cloud is truncated at $R_{cut}=100$ pc and thus the
star formation event occurs with a global efficiency $\epsilon\equiv
M_{SC}/M_{MC}\sim 3.7\%$. The metallicity in the host molecular cloud is
chosen as $2.5\thinspace\hbox{$\hbox{Z}_{\odot}$}$. This choice is motivated
by the rapid build-up of metals in starbursting quasar-host galaxies where the
gas metallicity reaches super-solar values already at $\sim 0.1$ Gyr (e.g.
Calura et al., 2014). According to the outputs from Starburst99, the values of
$L_{SC}$ and $V_{\infty}$ at the start of the cluster’s evolution are $\sim
3.12\times 10^{41}$ erg s-1 and $\sim 3500$ km s-1, respectively, implying a
star cluster’s mass loss rate of $\dot{M}_{w}\sim 8.0\times 10^{-2}$
$\hbox{M}_{\odot}$ yr-1.

The fractional HD abundance is selected as $X_{HD}=6.8\times 10^{-5}$ (Meece
et al., 2014). This implies, from evaluation of equation (19), that the
neutral supershell reaches the CMB temperature in $\sim 9.13\times 10^{4}$
years via HD cooling alone. However, if one accounts for dust cooling
(equation 20), the time required to reach $T_{cmb}=19.1$ K would be $\sim 350$
years if the supershell’s neutral section has an initial temperature $\sim
10^{4}$ K (right after it ceases to be ionised) and a density $\sim 9.4\times
10^{5}$ cm-3. Moreover, dust-induced cooling is largely capable of
counteracting the heating induced by H2 formation at densities $\geq 10^{6}$
cm-3 and metallicities $\geq 10^{-2}\thinspace\hbox{$\hbox{Z}_{\odot}$}$
(Meece et al., 2014).

Regarding the grain size distribution defined in equation (17), we have fixed
$\delta=0.05$, $a_{0}=0.1$ $\mu$m, $a_{min}=0.001$ $\mu$m and $a_{max}=0.5$
$\mu$m. For simplicity we consider only carbonaceous grains. We have imposed
that supershells effectively self-shield from the CMB only when
$[1-\exp(-\tau_{cmb})]\gtrsim 99.0\%$.

For the sake of comparison, seven additional cases (A1750, B3500, B1750,
C3500, C1750, D3500 and D1750) were also considered. The chosen parameters in
all of them fall within the observational constraints derived for the case of
a molecular cloud ($\sim 9\times 10^{7}\thinspace\hbox{$\hbox{M}_{\odot}$}$)
at $z\sim 6$, hosting a compact proto-globular cluster candidate with a
stellar mass of a few $\sim 10^{6}$ $\hbox{M}_{\odot}$(Vanzella et al., 2019;
Calura et al., 2021).

Model A1750 ($\epsilon\sim 3.7\%$) differs from the reference case A3500 as it
includes mass-loading that results in a terminal speed $V_{\eta}=1750$ km s-1.
Models B3500 and B1750 ($\epsilon\sim 3.7\%$) assume a flatter gas density
distribution ($\omega=1.5$). This choice is supported by the results of Lee et
al. (2015) and Raskutti et al. (2017). Models C3500 and C1750 ($\epsilon\sim
10.0\%$) consider a steeper gas density distribution ($\omega=2.5$). Finally,
models D3500 and D1750 ($\epsilon\sim 3.1\%$) explore the case of a central
star cluster with $10^{6}\thinspace\hbox{$\hbox{M}_{\odot}$}$. All models
denoted with 3500 have $V_{\eta}=3500$ km s-1. Models with with 1750
correspond to mass-loaded outflows with $V_{\eta}=1750$ km s-1,
$\eta_{ml}=3.0$, and $t_{ml}=3.8$ Myr for models A1750, B1750 and C1750, and
$t_{ml}=2.7$ Myr for model D1750. The central zone gas mass is
$M_{MC}(r<R_{c})\sim 9.2\times 10^{5}$ $\hbox{M}_{\odot}$ in all but models
D3500 and D1750. In those models, $M_{MC}(r<R_{c})$ is $\sim 3.3\times 10^{5}$
$\hbox{M}_{\odot}$.

In models A3500, A1750, C3500, C1750, D3500 and D1750, we have truncated the
clouds at $R_{cut}=100$ pc. For models B3500 and B1750 the clouds are
truncated at $R_{cut}=48$ pc so they have the same mass than in the reference
case. In all models we have selected a redshift $z=6$. Our calculations are
stopped once the supershells reach the edge of their host molecular cloud. The
input parameters are summarised in Table 1. In all the presented cases, the
outward-directed force resulting from the shock-heated gas thermal pressure
overcomes the inward-directed gravitational force acting on the supershell
during the whole supershell evolution:

$\displaystyle P_{b}>\frac{1}{4\pi
r^{2}}\left(\frac{GM_{sh}(r)M_{0}}{r^{2}}+\frac{1}{2}\frac{GM_{sh}^{2}}{r^{2}}\right)~{},$
(21)

where $M_{0}=M_{MC}(r<R_{c})+M_{SC}$. This warrants that the star cluster’s
mechanical energy is sufficient to overcome the binding energy of the star
cluster (Baumgardt et al., 2008) and form an expanding supershell.

In our reference model A3500, the supershell self-shields from the star
cluster ionising and non-ionising radiation as soon as it starts to grow (i.e.
both ($1-\phi$) and [$1-\exp(-\tau)$] go rapidly to unity).

Figure 3 presents the density structure of the starburst-driven supershell at
$0.5$ Myr of evolutionary time in the reference model A3500. In this case, the
fraction of CMB radiation that is absorbed within the supershell’s neutral
part is remarkably high, with $[1-\exp(-\tau_{cmb})]$ reaching $\gtrsim
99.9\%$, during $\sim 0.65$ Myr, and $\gtrsim 99.0\%$ during $\sim 1.0$ Myr.
In all models with mass-loading (decreased terminal speed), the supershell
self-shields from the CMB radiation for a longer time. Interestingly, if we
consider a flatter gas density distribution, as in cases B3500 and B1750, the
self-shielded gas mass increments monotonically, and the supershells in those
cases maintain $[1-\exp(-\tau_{cmb})]\gtrsim 99.9\%$ during the whole time
they evolve within their host molecular cloud (see Section 4).

Models C3500, C1750, D3500 and D1750 show similar trends to the ones observed
in models A3500 and A1750, but the time in which $[1-\exp(-\tau_{cmb})]\gtrsim
99\%$ is reduced to $\sim 0.36$ Myr, $\sim 0.4$ Myr, $\sim 0.6$ Myr, and $\sim
0.75$ Myr, respectively.

Figure 4: The supershell expansion velocity in different molecular gas density
distributions for models A3500 (black solid line), B3500 (red dashed line),
and C3500 (blue dotted line). Note that the calculations are stopped when the
supershells reach the edge of their host molecular clouds.

The supershell expansion velocity in different ambient gas density
distributions is shown in Figure 4, where solid, dashed and dotted lines
correspond to models A3500, B3500 and C3500, respectively. The supershells in
all our models do not stall within their host molecular clouds. However, in
clouds with a power-law index $\omega>2$, the supershells accelerate and,
inevitably, become Rayleigh-Taylor unstable (Krause et al., 2012). Note also
that in simulations with larger mass-loading rates (terminal speeds
$V_{\eta}<1500$ km s-1), catastrophic gas cooling sets in the shock-heated gas
zone that prevents the formation of supershells.

We have restricted our models to the spherically-symmetric case. However, if
the host molecular cloud is clumpy, the starburst-driven superbubble would
expand preferentially through channels in between clumps (Lucas et al., 2020).
We will explore that scenario in a forthcoming communication by means of
three-dimensional hydrodynamical simulations.

##  4 The Growth of Dust Grains

As a consequence of the supershell’s self-shielding against the incoming
CMB/starlight radiation, the restrictions for dust grain growth at high-
redshift (Ferrara et al., 2016) are evaded. Indeed in this case:

  * •

Dust grains in supershells can cool below $T_{cmb}$. For instance, a dust
grain would take just a few days to cool down to $15$ K or less (see Appendix
B).

  * •

The grains will not photo-desorb their accreted species, nor will they become
Coulomb repulsive.

The grain growth time-scale is given by (Spitzer, 1978)

$\displaystyle
t_{growth}=\left[\frac{3S\mathcal{D}\mu_{mol}n}{4\rho_{gr}a}\left(\frac{8k_{B}T}{\pi
m_{s}}\right)^{1/2}\right]^{-1},$ (22)

where $S=[0,1]$ is the sticking coefficient, and $m_{s}$ is the mass of the
accreting species. The redshift-dependent CMB photon number density is given
by (e.g. Grupen, 2005)

$\displaystyle
n_{\gamma}=16\pi\zeta(3)\left(\frac{kT_{CMB}(z)}{ch}\right)^{3}~{},$ (23)

where $\zeta$ is the Riemann zeta function. The CMB photon number density
scales with redshift as $\approx 411(1+z)^{3}$ cm-3, so its value (attenuated
by dust absorptions) at $z=6$ is $\approx 1.41\times
10^{5}[1-\exp(-\tau_{cmb}(r^{\prime}))]$ cm-3. The dust absorption rate per
unit volume for the case of CMB photons, is

$\displaystyle\kappa_{cmb}=n\sigma_{d}F_{\gamma}[1-\exp(-\tau_{cmb}(r^{\prime}))]\,~{},$
(24)

where $F_{\gamma}=n_{\gamma}\times c/\sqrt{3}$ is the mean CMB photon flux
(Heacox, 2017; Kim et al., 2019).

Figure 5, panel a, shows $\kappa_{cmb}$ calculated at the inner edge of the
neutral supershell. In all cases $\kappa_{cmb}$ is tiny during much of the
supershell’s evolution. This implies that the CMB radiation and the gas and
dust are thermally-decoupled. Panel b in the same Figure shows that,
initially, the self-shielded gas mass increases as it collects more mass.
However, depending on the cloud’s radial gas density distribution, the fast
supershell expansion eventually leads to a decreased supershell’s density.
This provokes a decline in the amount of self-shielded gas mass in the case of
a steep ambient gas density distribution (in cases A3500 and A1750 this occurs
after $\sim 0.5$ Myr). This decline does not occur for our models with a
flatter ($\rho\sim r^{-1.5}$) gas density distribution.

The grain growth time-scale within the neutral supershell is so fast in our
models (see panel c in Figure 5) and the dust-induced cooling is so efficient,
that one can expect that all the available refractory elements will inevitably
be accreted onto dust grains. For instance, for $a=0.01$ $\mu$m, $S=1$, and
$T=T_{cmb}(z=6)$, $t_{growth}$ is $\lesssim 100$ years during the first $\sim
2$ Myr, while if $S$ acquires a modest value, e.g. 0.1, $t_{growth}$ would be
$\lesssim 1000$ years during that period. Panel c also shows that adiabatic
cooling (lines above the axis break, see equation 30 in Appendix B) is at
play, but with a reduced effect due to its significantly longer characteristic
time-scale. Hence, grains can grow very rapidly in the extremely dense,
heavily-enriched and cold neutral part of the dusty supershell.

Figure 5: The starburst-driven supershell properties. Panel a presents he dust
absorption rate per unit volume. The amount of dust in the self-shielded
region as a function of time is presented in panel b. Lines below and above
the axis break on panel c display the adiabatic cooling and the dust grain
growth (for $S=1$ and $a=0.01$ $\mu$m) time-scales, respectively. In each
panel, we present the results for the A3500 (thick black solid lines), A1750
(thin black solid lines), B3500 (thick red dashed lines), B1750 (thin red
dashed lines), C3500 (thick blue dotted lines), C1750 (thin blue dotted
lines), D3500 (thick dash-dotted lines) and D1750 (thin dash-dotted lines)
cases, respectively. Figure 6: Bar graph of the growth of dust mass in the
self-shielded region for each studied case. The bars in purple, cyan and green
indicate the pre-existent swept-up dust mass in the self-shielded region (in
units $10^{3}\thinspace\hbox{$\hbox{M}_{\odot}$}$), the growth of the dust
mass due to accretion of refractory elements, and the estimated growth of dust
mass due to icy mantles, respectively. Figure 7: $T_{peak}$ as a function of
grain radius for the case of grains colliding with an HD-emitted photon (panel
a). The time required to cool from $T_{peak}$ to $T_{d}=10$ K as a function of
grain radius (panel b). In both panels the solid lines assume the initial
grain temperature to be $19.1$ K, while the dashed lines depict the case for
an initial grain temperature equal to $15$ K.

For an assumed maximum dust-to-gas mass ratio $\sim 1/100$ (Valentini &
Brighenti, 2015) at solar metallicity, scaling with $Z_{ISM}$, the dust mass
will have a net increase of $\sim 8.8\times 10^{4}$ $\hbox{M}_{\odot}$ for
model A3500. and $\sim 8.6\times 10^{4}$ $\hbox{M}_{\odot}$ for model A1750.
In models B3500 and B1750, the dust mass increase is $\sim 3.26\times 10^{5}$
$\hbox{M}_{\odot}$ and $\sim 3.31\times 10^{5}$ $\hbox{M}_{\odot}$,
respectively. The accelerating supershells in models C3500 and C1750 rupture
as they become Rayleigh-Taylor unstable, and consequently, it is uncertain
whether dust grain growth could occur. Nevertheless, we report the results of
our calculations for these models in Figure 6. Finally, the dust mass
increments are $\sim 3.1\times 10^{4}$ $\hbox{M}_{\odot}$ for both models
D3500 and D1750. The above quantities disregard the dust mass in the form of
ice mantles, which can further boost the grains’ mass. The bars in Figure 6
represent the amount of dust that was originally swept-up (purple), the amount
of dust that is obtained after all refractory elements are locked-up onto the
grains (cyan), and the grain mass in the form of ice mantles, which is assumed
to increase the total dust mass by a factor of two (Kruegel, 2003). As a
result, the supershells’ self-shielded sections will have, within $\sim 1$
Myr, three times more dust mass than the original self-shielded dust mass.

As expected, the most successful cases are those in which the stellar masses
are larger. If one assumes a cluster mass function
${\rm{dN}}/{\rm{d}}M_{\rm{*}}\propto M_{\rm{*}}^{\beta}$, with
$\beta\approx-2$ (e.g. Cook et al., 2019), in a galaxy with lower and upper
cluster masses $10^{3}\thinspace\hbox{$\hbox{M}_{\odot}$}$ and
$10^{7}\thinspace\hbox{$\hbox{M}_{\odot}$}$, respectively, then the number of
clusters with masses $\gtrsim 10^{6}\thinspace\hbox{$\hbox{M}_{\odot}$}$
scales as $\sim 10(M_{*}/10^{8}\thinspace\hbox{$\hbox{M}_{\odot}$})$. In the
above $M_{*}$ stands for the total stellar mass. If the typical net increment
in dust mass within CMB-dark supershells in clusters with $\gtrsim
10^{6}\thinspace\hbox{$\hbox{M}_{\odot}$}$ is $\sim
10^{5}\thinspace\hbox{$\hbox{M}_{\odot}$}$ (and all of them reside within
dense molecular clouds), then the proposed mechanism can account for the
build-up of $\sim
10^{6}\thinspace\hbox{$\hbox{M}_{\odot}$}(M_{*}/10^{8}\thinspace\hbox{$\hbox{M}_{\odot}$})$
of dust.

This rough estimate takes into account that the fraction of stars that form in
bound clusters grows with redshift, and at $z\sim 6$ they likely made the
larger proportion (e. g. Vanzella et al., 2019, and references therein).

We have seen so far that dust grain growth within starburst-driven supershells
can be both very rapid and very efficient in time-scales of order $\lesssim 1$
Myr, before any supernova goes off in the central starburst. However, between
$\sim 3-40$ Myr, the massive stars in the clusters will explode as core-
collapse supernovae. In that scenario, supernova blast waves will ram through
the hot gas cavity and collide with the encompassing supershell. Such
situation was studied by Martínez-González et al. (2018, 2019) by means of
three-dimensional hydrodynamical simulations. Upon the collision, the
supernova remnants become strongly radiative and do not experience the Sedov-
Taylor phase (see Figure 5 in Martínez-González et al., 2019). Moreover, blast
waves are mostly reflected and soon catch up with the reverse shock as they
only penetrate a very thin layer into the supershell (Tenorio-Tagle et al.,
1990). Thus, the dust grains locked-up in the supershell will remain largely
unaffected.

Supershells may eventually become pressure-confined, perhaps outside their
host molecular clouds, and the time required for them to stall is shorter than
the life-time of massive stars and their continuous winds (several tens of
Myr). This implies that the central cluster will continue to produce a high
velocity outflow ($\gtrsim 1000$ km s-1). The dense ($\gtrsim 10^{6}$ cm-3)
stalling supershell rests then above a lower density ($\lesssim 10^{-1}$
cm-3), hot ($\gtrsim 10^{7}$ K) gas. In that case the standing supershell
feels an inward gravity and becomes Rayleigh-Taylor unstable. Consequently,
the blast waves will collide with supershell fragments and the hot ejecta will
pass through channels in between them, establishing pressure equilibrium. The
blast waves will then attempt to ram through the fragments, but one can expect
that they will only be weakly transmitted. Thus the large majority of the
supershell’s dust grains will survive to eventually mix with the unshocked
ISM.

This situation differs from that described by Weaver et al. (1977) (our models
with $\omega\leq 2.0$), where decelerating supershells are stable because the
hot gas rests on top of the dense expanding supershell, and then a co-moving
parcel of fluid feels an effective outward gravity.

##  5 Concluding Remarks

Based on semi-analytic models, we have analysed the large-scale evolution of
starburst-driven supershells and the micro-physics of the dust grain growth
process within supershells at high-redshift ($z\sim 6$).

The main model predictions are summarised as follows:

  * •

Supershells can self-shield from the stellar radiation field and from the
Cosmic Microwave Background. The latter requires them to be sufficiently
metal-rich ($Z\geq\thinspace\hbox{$\hbox{Z}_{\odot}$}$), dense ($n\geq 10^{6}$
cm-3) and dusty ($\mathcal{D}\sim 1/150\times Z$), and that the gas density
distribution in the host molecular cloud is not too steep ($\omega\leq 2$).

  * •

As a result of dust-induced and HD cooling (and less importantly adiabatic
expansion), the supershell’s dust grains can cool below the temperature set by
the CMB ($T_{cmb}(z=6)=19.1$ K), e.g. to a temperature $\sim 15$ K or less
within a few days.

  * •

The grain growth time-scale in the self-shielded section of supershells is
sufficiently short ($\lesssim 100$ years) to allow all refractory elements
within the supershell’s self-shielded section to be accreted onto dust grains.

  * •

These grains do not photo-adsorb the accreted species, nor become Coulomb
repulsive, and thus can survive being self-shielded from the central cluster
starlight.

  * •

The total amount of dust produced by this mechanism may reach $\sim
10^{6}\thinspace\hbox{$\hbox{M}_{\odot}$}(M_{*}/10^{8}\thinspace\hbox{$\hbox{M}_{\odot}$})$.

This mechanism takes places before the first supernova explosion in the
central star cluster ($\sim 3$ Myr). Nevertheless, supernovae do not pose a
significant threat to the survival of the dust grains locked-up within the
encompassing supershell because the latter mostly reflects the colliding blast
waves, leaving its dust content largely unaffected (Martínez-González et al.,
2019). We have thus shown that dust grain growth at high redshift may take
place rapidly and efficiently.

##  6 Acknowledgements

This study was supported by CONACYT-México research grant A1-S-28458. S.M.G.
also acknowledges support by CONACYT through Cátedra n.482. 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. The authors
thank the anonymous Referee for a careful reading and helpful suggestions
which greatly improved the paper.

##  7 Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Asano et al. (2013) Asano R. S., Takeuchi T. T., Hirashita H., Nozawa T., 2013, MNRAS, 432, 637
  * Badjin et al. (2016) Badjin D. A., Glazyrin S. I., Manukovskiy K. V., Blinnikov S. I., 2016, MNRAS, 459, 2188
  * Baumgardt et al. (2008) Baumgardt H., Kroupa P., Parmentier G., 2008, MNRAS, 384, 1231
  * Bisnovatyi-Kogan & Silich (1995) Bisnovatyi-Kogan G. S., Silich S. A., 1995, Reviews of Modern Physics, 67, 661
  * Calura et al. (2008) Calura F., Pipino A., Matteucci F., 2008, A&A, 479, 669
  * Calura et al. (2014) Calura F., Gilli R., Vignali C., Pozzi F., Pipino A., Matteucci F., 2014, MNRAS, 438, 2765
  * Calura et al. (2015) Calura F., Few C. G., Romano D., D’Ercole A., 2015, ApJ Let, 814, L14
  * Calura et al. (2021) Calura F., et al., 2021, MNRAS, 500, 3083
  * Chiosi et al. (1998) Chiosi C., Bressan A., Portinari L., Tantalo R., 1998, A&A, 339, 355
  * Compiègne et al. (2011) Compiègne M., et al., 2011, A&A, 525, A103
  * Cook et al. (2019) Cook D. O., et al., 2019, MNRAS, 484, 4897
  * Draine (2003) Draine B. T., 2003, ApJ, 598, 1026
  * Draine (2011) Draine B. T., 2011, ApJ, 732, 100
  * Dwek (1986) Dwek E., 1986, ApJ, 302, 363
  * Dwek (1998) Dwek E., 1998, ApJ, 501, 643
  * Dwek & Cherchneff (2011) Dwek E., Cherchneff I., 2011, ApJ, 727, 63
  * Elmegreen (2017) Elmegreen B. G., 2017, ApJ, 836, 80
  * Elmegreen & Efremov (1997) Elmegreen B. G., Efremov Y. N., 1997, ApJ, 480, 235
  * Ferrara et al. (2016) Ferrara A., Viti S., Ceccarelli C., 2016, MNRAS, 463, L112
  * Ferrara et al. (2017) Ferrara A., Hirashita H., Ouchi M., Fujimoto S., 2017, MNRAS, 471, 5018
  * Gall & Hjorth (2018) Gall C., Hjorth J., 2018, ApJ, 868, 62
  * Gall et al. (2011) Gall C., Hjorth J., Andersen A. C., 2011, A&ARv, 19, 43
  * Grupen (2005) Grupen C., 2005, Astroparticle Physics, doi:10.1007/3-540-27670-X.
  * Heacox (2017) Heacox W. D., 2017, Contemporary Physics, 58, 190
  * Hensley & Draine (2017) Hensley B. S., Draine B. T., 2017, ApJ, 836, 179
  * Hirashita (2012) Hirashita H., 2012, MNRAS, 422, 1263
  * Indebetouw et al. (2014) Indebetouw R., et al., 2014, ApJ Let, 782, L2
  * Jiang et al. (2007) Jiang L., Fan X., Vestergaard M., Kurk J. D., Walter F., Kelly B. C., Strauss M. A., 2007, AJ, 134, 1150
  * Johnson & Bromm (2006) Johnson J. L., Bromm V., 2006, MNRAS, 366, 247
  * Johnson et al. (2015) Johnson K. E., Leroy A. K., Indebetouw R., Brogan C. L., Whitmore B. C., Hibbard J., Sheth K., Evans A. S., 2015, ApJ, 806, 35
  * Kim et al. (2019) Kim J.-G., Kim W.-T., Ostriker E. C., 2019, ApJ, 883, 102
  * Klessen & Glover (2016) Klessen R. S., Glover S. C. O., 2016, Saas-Fee Advanced Course, 43, 85
  * Koo & McKee (1992) Koo B.-C., McKee C. F., 1992, ApJ, 388, 93
  * Krause et al. (2012) Krause M., Charbonnel C., Decressin T., Meynet G., Prantzos N., Diehl R., 2012, A&A, 546, L5
  * Kruegel (2003) Kruegel E., 2003, The physics of interstellar dust
  * Laporte et al. (2017) Laporte N., et al., 2017, ApJ Let, 837, L21
  * Lee et al. (2015) Lee E. J., Chang P., Murray N., 2015, ApJ, 800, 49
  * Leitherer et al. (1999) Leitherer C., et al., 1999, ApJ, 123, 3
  * 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. (2014) Martínez-González S., Silich S., Tenorio-Tagle G., 2014, ApJ, 785, 164
  * 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., Muñoz-Tuñón C., Silich S., Tenorio-Tagle G., 2018, ApJ, 866, 40
  * Martínez-González et al. (2019) Martínez-González S., Wünsch R., Silich S., Tenorio-Tagle G., Palouš J., Ferrara A., 2019, ApJ, 887, 198
  * Mattsson et al. (2015) Mattsson L., Gomez H. L., Andersen A. C., Matsuura M., 2015, MNRAS, 449, 4079
  * McKee (1989) McKee C., 1989, in Allamandola L. J., Tielens A. G. G. M., eds, Vol. 135, Interstellar Dust. p. 431
  * Meece et al. (2014) Meece G. R., Smith B. D., O’Shea B. W., 2014, ApJ, 783, 75
  * Osterbrock (1989) Osterbrock D. E., 1989, Astrophysics of gaseous nebulae and active galactic nuclei. Mill Valley, CA: University Science Books
  * Rahner et al. (2017) Rahner D., Pellegrini E. W., Glover S. C. O., Klessen R. S., 2017, MNRAS, 470, 4453
  * Rahner et al. (2019) Rahner D., Pellegrini E. W., Glover S. C. O., Klessen R. S., 2019, MNRAS, 483, 2547
  * Raskutti et al. (2017) Raskutti S., Ostriker E. C., Skinner M. A., 2017, ApJ, 850, 112
  * Raymond et al. (1976) Raymond J. C., Cox D. P., Smith B. W., 1976, ApJ, 204, 290
  * Sahai & Nyman (1997) Sahai R., Nyman L.-Å., 1997, ApJ Let, 487, L155
  * Sahai et al. (2017) Sahai R., Vlemmings W. H. T., Nyman L. Å., 2017, ApJ, 841, 110
  * Spitzer (1978) Spitzer L., 1978, Physical processes in the interstellar medium, doi:10.1002/9783527617722.
  * Tenorio-Tagle et al. (1990) Tenorio-Tagle G., Bodenheimer P., Franco J., Rozyczka M., 1990, MNRAS, 244, 563
  * 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
  * Vanzella et al. (2019) Vanzella E., et al., 2019, MNRAS, 483, 3618
  * Vasiliev & Shchekinov (2005) Vasiliev E. O., Shchekinov Y. A., 2005, Astrophysics, 48, 491
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Whitworth (2016) Whitworth A. P., 2016, MNRAS, 458, 1815
  * Zhukovska et al. (2008) Zhukovska S., Gail H. P., Trieloff M., 2008, A&A, 479, 453
  * da Cunha et al. (2013) da Cunha E., et al., 2013, ApJ, 766, 13

##  Appendix A Turbulent Molecular Cloud

The pressure gradient in the static initial cloud is determined by the
equation (e.g. Calura et al., 2015)

$\displaystyle\frac{{\rm d}P_{t}}{{\rm d}r}=\frac{-GM(r)\rho}{r^{2}}~{},$ (25)

where $G$ is the gravitational constant. In the case of a power-law density
distribution (see equation 1), the mass enclosed within a radius $r$,
$M_{MC}(r)$, is

$\displaystyle
M_{MC}(r)=M_{0}+\frac{4\pi\rho_{c}R_{c}^{3}}{(3-\omega)}\left[\left(\frac{r}{R_{c}}\right)^{3-\omega}-1\right]~{},$
(26)

where $M_{0}$ is the total (the star cluster and the residual gas) mass within
the central zone with radius $r<R_{c}$. The second term in equation (26) gives
the mass contained in the swept-up supershell, $M_{sh}$.

In the case of a power-law density distribution, equation (25) is easily
integrated:

$\displaystyle P_{t}$ $\displaystyle=$
$\displaystyle\displaystyle\frac{GM_{0}\rho_{c}}{(\omega+1)R_{c}}\left(\frac{r}{R_{c}}\right)^{-(\omega+1)}+\frac{2\pi
G\rho_{c}^{2}R_{c}^{2}}{(3-\omega)(\omega-1)}$ (27) $\displaystyle\times$
$\displaystyle\left[\left(\frac{r}{R_{c}}\right)^{2(1-\omega)}-\frac{2(\omega-1)}{\omega+1}\left(\frac{r}{R_{c}}\right)^{-(\omega+1)}\right]~{},$

if $\omega>1$ and $P_{t}(\infty)=0$.

##  Appendix B Heating and Cooling in the Self-Shielded Neutral Shell

Dust grains residing within the supershell’s self-shielded, neutral part do
not experience stochastic temperature fluctuations due to collisions with CMB
photons, however, they still have temperature fluctuations as they collide
with neutral species and/or absorb HD-emitted ($h\nu_{10}\sim 1.1\times
10^{-2}$ eV) and infrared re-emitted photons. After each collision/absorption
event, a dust grain with initial temperature $T_{0}$ will be heated to a peak
temperature, $T_{peak}$, given by

$\displaystyle E=\int_{T_{0}}^{T_{peak}}C(a,T_{d})\mbox{ d}T_{d}.$ (28)

where $C(a,T_{d})$ is the grain heat capacity as a function of the grain’s
radius, chemical composition and temperature, $T_{d}$. Upon reaching
$T_{peak}$, the grain starts to cool down in a time-scale given by (Dwek,
1986; Martínez-González et al., 2016, 2017)

$\displaystyle t_{d,cool}=\int_{T_{d}}^{T_{peak}}\frac{C(a,T_{d})\mbox{
d}T_{d}}{|4\pi a^{2}\sigma_{SB}\langle Q_{abs}\rangle T_{d}^{4}|},$ (29)

where $\langle Q_{abs}\rangle$ is the Planck-averaged grain absorption
efficiency. Figure 7 shows $T_{peak}$ and $t_{d,cool}$ as a function of grain
radius for the case of grains colliding with an HD-emitted photon, whose
energy roughly corresponds with the energy of an infrared photon with
wavelength $\sim 100$ $\mu$m. We have assumed values of $T_{0}$ equal to $15$
K and $19.1$ K. A significant temperature increase will be expected only for
grains $a\leq 0.003$ $\mu$m, while $t_{d,cool}$ between $T_{peak}$ and
$T_{d}=10$ K takes only some tens of days.

On the other hand, the supershell cools as it expands and the characteristic
time-scale for adiabatic cooling is (Badjin et al., 2016)

$\displaystyle t_{adiab}=\frac{(\gamma+1)}{6(\gamma-1)}\frac{R_{S}}{V_{S}},$
(30)

where $R_{S}$ and $V_{S}$ are the forward shock radius and velocity.