Supernova 2009ip outbursts in 2012: From scenario to model

© 2022 . N. N. Chugai111email<EMAIL_ADDRESS>
1Institute of astronomy, Russian academy of science, Moscow

Keywords:  stars — evolution; stars — supernovae — SN 2009ip

PACS codes:

###### Abstract

Spherical and aspherical models are presented for two outbursts in 2012 of
supernova 2009ip. Models are based on a scenario which suggests that the
August 2012 outburst is caused by the explosive shell ejection from LBV-
presupernova. The model predicts an emergence of an unobserved outburst in
late July 2012 related to a shock breakout and a subsequent diffusive
radiative cooling of the ejected envelope. The luminosity of the first
observed outburst in August 2012 was presumably powered by the central source,
whereas the second, more powerful outburst in late September 2012, was caused
by the ejecta interaction with the circumstellar envelope. Models provide
estimates of the ejecta energy and mass along with the mass of the
circumstellar shell.

##  1 Introduction

Supernova SN 2009ip in NGC7259 galaxy is a complex phenomenon following the
20/08/2009 outburst with the absolute magnitude of -14.5 (Maza et al. 2009).
Presupernova is found in archive image of NGC7259 taken by the Hubble
telescope (HST hereafter) with the absolute magnitude of -9.5 and was
identified as LBV (Smith et al. 2010). Two years after the discovery SN 2009ip
experienced several outbursts in the range of -10.5… -13.5 in the span of 160
days (Pastorello et al. 2013).

Key events happened after 2012 August 8, when SN 2009ip experienced outburst
up to -15 mag (aka 2012a outburst) followed by the late September, 2012
powerful outburst (2012b outburst) up to -18 mag (Pastorello et al. 2013) with
the subsequent monotonic flux decrease. In December 2021 the absolute
magnitude in F606W filter of HST is found to be -8.5, one magnitude fainter
compared to -9.5 mag in 1999, so authors concluded that the LBV-presupernova
has gone (Smith et al. 2022).

Starting with the discovery spectra have been showing strong narrow hydrogen
emission lines with the FWHM of 800-1000 km s-1 (Pastorello et al. 2013,
Margutti et al. 2014) indicating heavy mass loss. The first spectrum
09/08/2012 of the 2012a outburst revealed broad hydrogen P Cyg line profiles
with wing velocities up to 14000 km s-1 (Pastorello et al. 2013). This
observation argued a scenario in which the 2012a outburst was the supernova
explosion, whereas the subsequent more powerful outburst 2012b was caused by
the interaction of the supernova ejecta with the circumstellar envelope
(Mauerhan et al. 2013, Pastorello et al. 2013). An alternative possibility
suggests two explosive events related respectively with 2012a and 2012b
outbursts with the subsequent CS interaction after the second explosion that
presumably occured 13/09/2012 (Margutti et al. 2014). Details of this scenario
are bounded by analytical estimates. For the first, single explosion scenario,
the light curve modelling and estimates of supernova energy and ejected mass
are absent at all. Yet this scenario deserves detailed study since it requires
only one explosion.

The identification of the SN 2009ip presupernova with the massive LBV implies
that the observed phenomena in this case are related to nuclear burning
flashes in oxygen core initiated by the instability cased by pair creation
(Woosley 2017, Smith et al. 2022). Supernovae related to such phenomena are
called pulsational pair-instability supernovae (PPISN for short). Yet one has
to emphasise that full resemblance between SN 2009ip and PPISN models (Woosley
2017) is lacking. The crucial mismatch is the point that SN 2009ip spectra
during 2012a outburst show high expansion velocity up to 13000-14000 km s-1,
which is 5-10 larger compared to velocities of ejected shell in PPISN models.

The paper addresses the scenario of a single explosion of a LBV related to the
2012a outburst of SN 2009ip. The presented model suggests explosive shell
ejection (2012a outburst) and its subsequent interaction with the CS envelope
(2012b outburst). Given significant polarization of SN 2009ip radiation
(Mauerhan et al. 2014) we explore models with spherical and aspherical
distribution of the CS matter.

##  2 Light curve modelling

###  2.1 Preliminary consideration

Based on the premis that 2012a outburst is caused by the explosive mass
ejection by the LBV one should adopt the presupernova radius of 70 $R_{\odot}$
in line with the evolutionary models of LBV with masses of 80-90
$M_{\odot}$(Woosley 2017). With such a radius the initial luminosity related
to the shock breakout (call it "2012a0" outburst) can exceed the luminosity of
2012a outburst. Moreover, since the initial luminosity peak of 2012a0 outburst
and the subsequent light curve has to differ from the light curve of 2012a,
one should admit that the initial 2012a0 outburst was missing, so the
supernova explosion occured somewhat earlier than turn on of the 2012a
outburst. Yet the explosion should be later than the flash of 24/07/2012 at
the level of $V\approx 17.4$ (Drake et al. 2012).

Along with the explosive ejection of the supernova (2012a outburst), as we see
below, one has to admit additional thermal energy injection into the expanding
envelope of the order of $\sim 2\times 10^{48}$ erg. The supernova light curve
will be calculated based on the analytical model of the radiation diffusion
with the adiabatic cooling and the additional source of the thermal energy
(Arnett 1982). The subsequent CS interaction is treated based on the thin
shell approximation (Giuliani 1982, Chevalier 1982). The detailed realization
of the CS interaction light curve model is described earlier (Chugai 2001,
2021). We therefore omit model description, but remind that the supernova
shell deceleration in the dense CS gas proceeds via formation of outer
(direct) and reverse shock waves with the formation of a cold dense shell
(CDS) inbetween. The CDS can be optically thick and in this case the supernova
shows smooth continuum without broad absorption lines; this phenomenon took
place in SN 1998S first month after explosion (Chugai 2001).

We consider spherical and aspherical options for the CS mater distribution.
Non-sphericity is supported by the significant polarization of the SN 2009ip
radiation during 2012b outburst (Mauerhan et al. 2014). Non-spherical model
will be calculated based on the spherical model with the CS interaction
luminosity reduced by an asphericity factor being equal to the ratio of a
visible area of the aspherical-to-spherical photosphere: $\xi=A_{ns}/A_{s}<1$.
To a first approximation $\xi$ is equal to the ratio of radii squared of the
photometric photosphere (Margutti et al. 2014) and the photosphere of a
spherical model. In fact, the photometric radius determined based on the
blackbody approximation is lower than the geometric value because the escaping
flux from the atmosphere with the absorption and scattering is lower than
blackbody flux by a factor of $\sqrt{\epsilon}/(1+\sqrt{\epsilon})<1$, where
$\epsilon$ is the photon loss probability per one act of extinction (cf.
Mihalas 1978). It should be emphasised that in aspherical model the supernova
envelope is assumed to be spherical; asphericity refers to the CS shell only.

Preliminary modelling of the 2012b outburst indicates that a steep luminosity
rise requires the steep density drop in outer layers ofthe supernova shell
$\rho\propto v^{-\omega}$ with $\omega\geq 20$. The density in the inner
layers are assumed to be uniform. The CS shell should have a sharp inner
boundary as well to provide the steep luminosity rise. The inner zone of the
CS medium is presumably produced by a steady wind of a moderate density with
the density parameter $w=\dot{M}/u=10^{14}$ g cm-1 that corresponds to the
mass loss rate of $1.5\times 10^{-4}$ $M_{\odot}$ yr-1 for the speed of 1000
km s-1. The dense CS shell presumably flows with the same speed of $u=1000$ km
s-1. For the CS shell mass of $\sim 1$ $M_{\odot}$ its Thomson optical depth
in the model is less than unity. which permits us to neglect by the radiation
diffusion delay compared to the light travel time $r/c$ and adopt the
approximation of the instant radiation escape. The modelling shows that the
light travel effect does not modify significantly the light curve of 2012b
outburst as well.

Model parameters (kinetic energy and mass of supernova ejecta, radius and mass
of the CS envelope) are constrained by the light curve and maximal velocity of
the undisturbed supernova ejecta that are determined from a blue wing of the
absorption component of H$\alpha$, H$\beta$ at the 2012a-stage and He I 10830
Å at 2012b-stage.

Table 1 Parameters of supernova and CS envelope

Model |  $E_{51}$ a | $M/M_{\odot}$ |  $t_{inj}$ (d) |  $L_{0,40}$ b | $M_{cs}/M_{\odot}$ |  $r_{1,15}$ c | $r_{2,15}$
---|---|---|---|---|---|---|---
sm | 0.34 | 0.2 | 24 | 90 | 1.3 | 8.4 | 13
nsm | 2 | 0.8 | 22 | 110 | 3.2 | 8.8 | 13.4

Units: a $10^{51}$ erg; b $10^{40}$ erg s-1; $10^{15}$ cm

###  2.2 Modelling results

The supernova light curve description (2012a outburst), as noted earlier,
requires a slow post-explosion injection of the thermal energy. The temporal
behavior of the power released by a central source is adopted as

$L=L_{0}(t/t_{inj})^{p}\exp{[-(t/t_{inj})^{2})]}\,,$ (1)

where the power index $p\sim 1.5-2$, while $t$ is the time lapsed from the
adopted explosion moment JD 2456133.5 (25/07/2012).

First month after the light maximum of the 2012b outburst the CDS is optically
thick as indicated by the smooth continuum. Since the ejecta mass and the CDS
optical depth decreases as the kinetic energy decreases, one expects the
minimal acceptable energy. The modelling indicates that the kinetic energy
cannot be lower than $4\times 10^{50}$ erg. For the spherical fiducial model
(model sm, Table 1, Figure 1) we adopt $E=3.4\times 10^{50}$ erg. For the
aspherical model (nsm, Table 1, Figure 2) with the asphericity parameter
$\xi=0.05$ the minimal energy is significantly higher, $E=2\times 10^{51}$
erg. The Table 1 includes supernova ejecta energy and mass, the timescale and
power of the additional energy injection, CS mass, the inner and outer radius
of the CS shell.

Both models describe the bolometric light curves of 2012a and 2012b outbursts
and are consistent with the observational maximal velocities of unshocked
supernova ejecta at about 100 days when CDS becomes transparent in the
continuum. Noteworthy, the model expansion velocity of the undisturbed
supernova envelope can be larger, but cannot be lower, than the maximum
velocity recovered from broad absorption lines.

At the very late phase (400 days) in the Keck-2 spectrum (Graham et al. 2017)
the velocity estimated from wings of emission lines of H$\beta$, He I 7065 Å,
and Fe II 5018 Å, are equal to 1500-2000 km s-1 that is lower than the CDS
velocity and maximal velocity of the undisturbed supernova shell in the
spherical model (Figure 1). In a context of the spherical model this suggests
that the contribution of the CDS and unshocked ejecta in these lines is
relatively small, whereas the dominant source of emission lines might be
shocked dense CS clouds following a scenario for the intermadiate component of
emission lines in SN 1988Z (Chugai & Danziger 1994). In the aspherical model,
however, the CDS velocity on day 400 ($\approx 1400$ km s-1, Figure 2) is
close to the emission lines FWHM, so the CDS could be the emitting site for
these lines.

An easy way to realise an aspherical model would be the supernova interaction
with a CS shell arranged as a disk (torus). This structure is invoked, e.g.,
to account for the co-existence of broad emission lines and the intermediate
component in SN 1988Z (Chugai & Danziger 1994). A similar structure is
considered to account for the polarization of the continuum originated
presumably from the shocked interaction at the inner edge toroidal CS shell
(Mauerhan et al. 2014). In the aspherical model CS shell asphericity parameter
$\xi=0.05$ corresponds to the equatorial belt in the latitude range $\pm
3^{\circ}$ for the inclination angle of $90^{\circ}$. The value $\xi=0.05$
approximately corresponds to the ratio of radii squared of the photometric
photosphere (Margutti et al. 2014) and the photosphere of spherical model with
a possible correction on the lower surface brightness compared to the
blackbody intensity.

A high density contrast of the toroidal CS shell with respect to the supernova
envelope result in a situation when almost unaccelerated CS torus submerges
into the expanding supernova envelope. Eventually we would obtain a picture in
which the interaction region at the inner torus edge becomes an internal
source for the continuum radiation on which background broad absorption lines
could form. On the other hand, the presence of absorption lines in the
spectrum at the 2012b-stage is a serious problem for the spherical model which
makes the aspherical model more advantageous.

Noteworthy, the postmaximum fast luminosity decline of 2012b outburst can be
well described by the power law $L\propto t^{-6}$ that follows from a model of
the shock interaction of homologously expanding envelope with the spherical
shell of the radius $R$ at rest. Indeed, the incoming speed at the radius $R$
decreases with time as $v\propto R/t$, whereas the density drops as
$\rho\propto t^{-3}$, so the luminosity of the reverse shock is $L\propto
R^{2}\rho v^{3}\propto t^{-6}$. This simple model can reflect the supernova
interaction with the dense toroidal CS shell that is weeakly accelerated by
the rarefied flow of supernova ejecta.

##  3 Discussion and conclusion

The paper has been aimed at the modelling of two outbursts of SN 2009ip in
2012 based on a scenario of explosive shell ejection by LBV presupernova. The
scrnario suggests that the supernova explosion corresponds to the relatively
weak outburst in August 2012 (2012a outburst), whereas powerful outburst of
late September 2012 (2012b outburst) is caused by the shock interaction of
supernova shell with the CS shell (Pastorello et al. 2013, Mauerhan et al.
2013). According to the model the explosive ejection of the envelope by the
LBV presupernova brings about strong initial flash 2012a0, preceded 2012a,
with the luminosity of $\sim 10^{42}$ erg s-1 and the timescale of one week
that is missing in observations.

The light curve of both outbursts 2012a and 2012b and the expansion velocity
of the supenova envelope are reproducible in both spherical model and models
with toroidal (disk) CS shell. However, if for the spherical model the energy
$\sim 3\times 10^{50}$ erg is sufficient, the aspherical model, on the other
hand, requires substantially higher value, of $\sim 2\times 10^{51}$ erg that
may be in discomfort with PPISN models (Woosley 2017). Summarizing,
phenomenologically the aspherical model is preferred, although it poses a
serious problem for the explosion theory.

Two features of the presented models should be emphasised. First, the model
supernova light curve (2012a outburst, Figures 1, 2) requires the prolong
additional thermal energy injection. This energy is presumably supplied by the
continued nuclear burning flash. Second, the explosive ejection of the shell
related to the 2012a outburst is preceeded by the initial week-long outburst
after the shock breakout (2012a0 outburst). The existence of this outburst
around 25/07/2012 cannot be verified by available observations.

This initial outburst is a generic feature of the LBV scenario with the
presupernova radius of 70 $R_{\odot}$ and is appropriate for the verification.
In this regards one should be emphasise that for the presupernova with the low
radius, say 10 $R_{\odot}$, the initial luminosity peak 2012a0 turns out to be
rather weak, at the level of initial luminosity of 2012a outburst. Since the
full resemblance between SN 2009ip and PPISN models doesn’t take place, the
case of low presupernova radius cannot be rulled out. This remark is directly
related to the test of the LBV scenario based on the observational
presence/absence of the initial 2012a0 peak.

Following the termination 2012a outburst the LBV with the presupernova mass
reduced by $\sim 4$ $M_{\odot}$ due to supernova and CS shells probably
experienced collapse into a black hole, since the HST photometry in December
2021 suggests that the LBV-presupernova has gone (Smith et al. 2022). Since
the presupernova mass according to the HST photometry in 1999 was estimated as
50-80 $M_{\odot}$ (Smith et al. 2010) the black hole mass should be in the
range of 50-80 $M_{\odot}$.

Noteworthy, SN 2009ip is only partially unique, primarily due to a long set of
observational data. As to the phenomenon of 2012a and 2012b outburst, the
similar photometric and spectral behavior was observed for SN 2010mc (Ofek et
al., Smith et al. 2014). Yet it is premature to conclude on the full
resemblance, since SN 2010mc spectra at the stage corresponding to 2012a
outburst are lacking, so one cannot say, whether the supernova explosion
occurred long before the major outburst.

This study was partially supported by RFBR and DFG, project number
21-52-12032.

References

Arnett W. D., Astrophys. J. 253, 785 (1982)

Chevalier R. A., Astrophys. J. 259, 302 (1982)

Chugai N. N., Mon. Not. R. Astron. Soc. 508, 6023 (2021)

Chugai N. N,, Mon. Not. R. Astron. Soc. 326, 1448 (2001)

Chugai N. N., Danziger I. J., Mon. Not. R. Astron. Soc., 268, 173 (1994) Drake
A. J., ATel 4334 (2012)

Graham M. L., Bigley A., Mauerhan J. C. et al., Mon. Not. R. Astron. Soc. 469,
1559 (2014)

Giuliani J. L., Astrophys. J. 256, 624 (1982)

Margutti R., Milisavljevic D., Soderberg A. M. et al., Astrophys. J. 780, 21
(2014)

Mauerhan J., Williams G., Smith N. et al., Mon. Not. R. Astron. Soc. 442, 1166
(2014)

Mauerhan J., Smith N., Filippenko A. V. et al., Mon. Not. R. Astron. Soc. 430,
1801 (2013)

Maza J. et al., Cent. Bur. Electron. Telegrams, 1928, 1 (2009)

Mihalas D. Stellar armospheres (San Francisco: Freeman and company, 1978)

Ofek E. O., Sullivan M., Cenko S. B. et al., Nature, 494, 66 (2013) Pastorello
A., Cappellaro E., Inserra C., et al., Astrophys. J. 767, 1 (2013)

Smith N., Andrews J., Filippenko A. V. et al., Mon. Not. R. Astron. Soc. 515,
715 (2022)

Smith N., Mauerhan J. C., Prieto J. L., Mon. Not. R. Astron. Soc. 438, 1191
(2014)

Smith N., Miller A., Li W., Filippenko A. V. et al., Astron. J. 139, 1452
(2010)

Woosley S. E., Astrophys. J. 836, 244 (2017)

Figure 1: Left. Bolometric light curve (thick line) in the spherical model
(model sm, Table 1) compared to the observational light curve (crosses from
Margutti et al. 2014; Y-symbols from Graham et al. 2017). First unobserved
maximum corresponds to the diffusion radiation cooling of supernova after the
shock breakout, the second is the 2012a outburst that reflects the supernova
emission powered by the additional energy release by the central source, the
third maximum is the result of the supernova interaction with the CS shell.
Right. Model velocity of the CDS (thick line) and the boundary velocity of the
undisturbed supernova ejecta compared to the maximal velocity in the blue wing
of absorption in H$\alpha$ and H$\beta$ (empty diamonds) and He I 10830 Å
(filled diamonds). Inset shows CS density distribution.

Figure 2: The same as in Figure 1, but for the apherical model (nsm, Table 1).


# Scattering of relativistic electron beams by the anode mesh in high-current
vircators

Anishchenko S.V<EMAIL_ADDRESS>Research Institute for Nuclear Problems
of Belarusian State University
Bobruiskaya Str. 11, 220030 Minsk, Belarus  Baryshevsky V.G<EMAIL_ADDRESS>Research Institute for Nuclear Problems of Belarusian State University
Bobruiskaya Str. 11, 220030 Minsk, Belarus  Gurinovich A.A<EMAIL_ADDRESS>Research Institute for Nuclear Problems of Belarusian State University
Bobruiskaya Str. 11, 220030 Minsk, Belarus

###### Abstract

In a virtual cathode oscillator, the scattering of a high-current relativistic
electron beam by the anode mesh leads to formation of an electron cloud near
the anode. The cloud particles possess low energy and large spread in
velocities caused by multiple scattering and ionization losses. The electrons
captured by the cloud do not participate in the oscillations of the virtual
cathode and partially block a vircator. As a result, the amplitude of the
electric field oscillations is reduced. In order to increase the oscillation
amplitude, the thickness of the anode mesh should be equal to the mean free
path for electrons in the mesh material.

vircator, multiple scattering, ionization losses, radiation losses

###### pacs:

41.75.-i, 34.50.Bw

##  I Introduction

The interaction of charged particles with matter plays an important role in
various fields of science and technology. Radiation therapy, design of modern
particle detectors, and radiation protection of spacecrafts are impossible
without a thorough analysis of the passage of charged particles through matter
and quantitative treatment of multiple scattering and energy losses.

The interaction of charged particles with matter is also important in high-
current electronics dealing with the propagation of high-power ion and
electron beams in electrodynamic structures Creedon1990 . Indeed, some of the
beam particles, when they inevitably hit one of the structural elements (anode
mesh, collector, drift tube, etc.), are reflected due to multiple scattering,
return back into the interaction area and continue to interact with
electromagnetic fields inside the system. For example, about a quarter of
particles normally incident on a steel element of the structure are reflected
off. If the incidence angle significantly differs from normal then the
fraction of reflected particles is much higher. The spectrum of kinetic energy
for reflected particles extends from zero to the initial energy. Since the
number of incident electrons is comparable with the number of reflected ones,
the latter can significantly affect the operation of high-current devices.

A large number of theoretical results devoted to multiple scattering of
particles are unsuitable for studying the interaction of high-current electron
beams with matter. First, many studies use a small-angle approximation
Williams1939 ; GS1940 ; Moliere1948 ; Bethe1953 ; Nigam1959 ; Marion1967 ;
Voinov1969 . Secondly, the approximate values of elastic electron-atom
scattering cross sections differ significantly from those obtained by
numerical calculation at energies < 1 MeV Mayol1997 . Thirdly, multiple
scattering cannot be considered separately from ionization Bohr1913 ;
Bethe1932 ; Bloch1932 ; BetheAshkin1953 ; Rohrlich1953  and radiation
Bethe1934 ; Pratt1977  losses. These losses are significant in the case of
<<thick>> electrodynamic structures.

A consistent quantitative description of particle-matter interaction can be
obtained using Monte-Carlo simulations Hebbard1955 ; Berger1963a ; Berger1963b
; Berger1964 ; Perkins1962 ; Seltzer1974 ; Seltzer1991 ; Jenkins1987 ;
Fernandez1993  based on the most rigorous Goudsmit-Sanderson multiple
scattering theory GS1940 . The algorithms used in numerical calculations
should be, on the one hand, sufficiently accurate, and, on the other hand,
fast.

In this paper, we describe an approach to modeling the interaction of electron
beams with electrodynamic structures in high-current electronic devices. The
approach will be validated by calculating electron transmission and reflection
coefficients and comparing the obtained results with the numerical simulations
Berger1963b ; Seltzer1974  and experimental data Tabata1971 ; Ito1993
published in the literature.

After validation, the numerical method will be integrated into a one-
dimensional program for modeling high-current devices with an oscillating
virtual cathode (VC). We will demonstrate that the scattering of relativistic
particles by the anode mesh leads to the formation of an electron cloud near
the anode. The cloud particles possess large energy spread and cause
significant decrease in the amplitude of field oscillations in vircators. It
will be shown that the use of an anode mesh with the thickness approximately
equal to the electron mean free path in the anode material leads to decrease
in the number of cloud particles and could contribute to increase of the
oscillation amplitude.

##  II Passage of relativistic electrons through matter

###  II.1 Monte-Carlo simulation

The basis of numerical simulation of electron-matter interaction is the Monte-
Carlo method Berger1963a ; Seltzer1974 ; Seltzer1991 ; Jenkins1987 , which can
be described as follows. The trajectory of each electron in matter is divided
into many small segments. For each segment, the electron energy is assumed to
be constant. When passing from one segment to another, the particle changes
its energy in accordance with the theory of ionization and radiation losses
Rohrlich1953 . Angular distribution of scattered particle is described by the
Goudsmit-Sanderson distribution Berger1963a . The calculation is carried out
until the particle either leaves the substance or loses a significant part of
initial energy (in practice, calculations stop when the electron energy
approaches 10 keV).

A significant drawback of the Monte-Carlo method is the computational
complexity. To calculate the Goudsmit-Sanderson distribution, it is necessary
to sum a large number of terms. Each term contains an integral of a rapidly
oscillating function. This circumstance significantly complicates the use of
the standard Monte-Carlo method in high-current electronics due to the huge
number of particles used in the simulation of high-current devices.

A way out in this situation could be a method, in which the distribution
function of the scattering angles would be found using a simple formula. And
such a method was proposed in Fernandez1993 . It is based on a rigorous
theoretical results obtained by Lewis Lewis1950 . According to Lewis1950 , the
most significant mean values depend only on the transport cross-sections: the
first one

$\sigma_{1}=2\pi\int_{0}^{\pi}(1-\cos\chi)\frac{d\sigma(\chi)}{d\Omega}\sin\chi
d\chi$ (1)

and second that

$\sigma_{2}=2\pi\int_{0}^{\pi}\frac{3}{2}(1-\cos^{2}\chi)\frac{d\sigma(\chi)}{d\Omega}\sin\chi
d\chi.$ (2)

Here, $\frac{d\sigma(\chi)}{d\Omega}$ is the differential elastic electron-
atom cross-section. For example, the average longitudinal velocity $v$ and the
root-mean-square deviation of the transverse velocity in a segment of length
$s$ 111The length of $s$ must be much less than the interval over which the
particle energy changes significantly. change in accordance with the formulas
as follows, respectively:

$v<\cos\theta>=v\exp(-n\sigma_{1}s)$ (3)

and

$v^{2}(1-<\cos^{2}\theta>)=v^{2}\Big{(}1-\frac{1+2\exp(-n\sigma_{2}s)}{3}\Big{)}.$
(4)

Symbol $n$ denotes the particle density, $\theta$ is the polar angle between
two vectors of particle velocity. The first vector corresponds to the entering
to the segment $s$ and the second one corresponds to the segment exit.

The idea stated in Fernandez1993  is as follows. If a simple distribution
leads to the formulas (3) and (4), then the computational complexity of the
Monte-Carlo method is significantly reduced. At the same time, the results of
modeling particle passage through matter does not change. This was
convincingly demonstrated by authors Fernandez1993  by calculating the angular
distribution of electrons after passed through gold foils of various
thicknesses.

###  II.2 Reflection and transmission coefficients

To verify the described approach to simulation of the interaction of electrons
with matter, we investigated electron scattering by thick foils made of
various materials (Be, C, Al, Fe, Ag, Au, U) and calculated reflection
coefficients. The energy of particles normally incident on a target varied
from 0.1 to 2 MeV. Comparison of our simulation results and experimental data
Tabata1971 ; Ito1993  demonstrates better agreement for heavy elements (left
plot in figure 1).

Simulation results obtained for light elements (right plot in figure 1)
deviate from experimental data in a greater extent: this is due to ignoring
the electron-electron collisions in the simulation that increases the angle of
multiple scattering by approximately $1+1/Z$ times. The nuclear charge in the
denominator indicates that this phenomenon is more significant for light
elements. Difference observed in simulation results and experiment for heavy
elements (see curves for gold and uranium in left plot in figure 1) could be
explained by rearrangement of the shells of valence electrons in a solid
matter.

Figure 1: Electron reflection coefficients: solid lines depict our simulation
results, points correspond to experimental data Tabata1971 ; Ito1993 .

Figures 2, 3 illustrate comparison of simulation results with calculations
Seltzer1974  for transmission and reflection coefficients for electrons
scattered by aluminum and beryllium foils of different thicknesses; comparison
with experimental data Ito1993  is also provided. Our results are in good
agreement with the calculations Seltzer1974 . However, experimentally obtained
data for reflection coefficient (right plots in Figures 2, 3) differ from both
our simulations and calculations Seltzer1974  that is most likely due to the
neglect of electron-electron collisions.

Figure 2: Transmission (left) and reflection (right) coefficients for
electrons scattered by aluminum foils of various thicknesses: solid lines
correspond to the simulation results; red, green and black dots are the
calculations from the paperSeltzer1974 ; blue dots depict experimental data
Ito1993

Figure 3: Transmission (left) and reflection (right) coefficients for
electrons scattered by beryllium foils of various thicknesses: solid lines
correspond to the simulation results; red, green and black dots are the
calculations from the paperSeltzer1974 ; blue dots depict experimental data
Ito1993

Figure 4 shows the dependence of the reflection coefficient on the kinetic
energy of a particle incident on an aluminum plate. Three solid curves present
the simulation results obtained for different angles of incidence. At energy
values up to $\sim 2$ MeV our simulation results are in good agreement with
both calculation Berger1963b  and experimental data Ito1993 . At the energy
$\sim 2$ MeV, the reflection coefficients given in Berger1963b  demonstrate a
noticeable decrease in contrast to our calculations. Such behavior applies to
the case of normal incidence.

Slight deviation of our simulation results from those presented in Berger’s
paper Berger1963b  could be associated with different approaches to the
consideration of multiple scattering. We used the approach Fernandez1993
based on rigorous results obtained by Lewis Lewis1950  from the Goudsmit-
Sanderson theory, while Berger resorted to the approximate Moli$\grave{e}$re
theory. Note that experimental data obtained at normal incidence
($\theta=0^{o}$) (figure 1) does not show such a decrease in the reflection
coefficient.

Figure 4: Dependencies of the electron reflection coefficient on the kinetic
energy at different angles of incidence on an aluminum plate: solid lines
depict our simulation results and dots show calculation results presented in
Berger1963b .

Thus, the simulation results obtained in this work demonstrate good agreement
with experimental data and numerical calculations published in the literature.
Some discrepancy ($\sim 20\%$) with the experimental data is observed in the
reflection coefficient for electrons with energy $\sim 0.1$ MeV scattered by
beryllium foil.

##  III VC oscillations

Figure 5: One-dimensional model of a high-current device with an oscillating
VC.

The approach developed to simulate the passage of electrons through matter was
integrated into a one-dimensional code designed to simulate high-current
devices with an oscillating virtual cathode. This code calculates motion of
relativistic electrons in a self-consistent longitudinal electric field by the
particle-in-cell method. The computational domain consists of two parts
(figure 5): the cathode-anode gap and the drift space. The potential of the
right wall of the drift space can be either equal to the anode potential (in
the case of a vircator) or to the cathode potential (in the case of a reflex
triode). Injection of particles into the system is carried out under
conditions of unlimited emission capability of the cathode. This condition
corresponds to the regime of explosive electron emission that takes place in
high-current accelerators.

The initial version of the PIC code took no account of the scattering by anode
mesh: the anode was assumed to be semitransparent and was characterized by a
single parameter, the geometric transparency $\eta$ which varied from zero to
one. When a particle passed through the anode, its charge was multiplied by
the transparency coefficient. The procedure corresponded to the partial
absorption of particles by the anode mesh. Neither the reflection of
particles, nor the deceleration of particles in the anode material, nor spread
of particle velocity due to multiple scattering were taken into account by the
code. In what follows, we will refer to the simulation model just described as
the absorption model. Figures 6 and 7 show the electric field oscillations,
spectrum, and the phase portrait of the beam in the vircator in the absorption
model. The geometric transparency is $\eta=0.7$.

In the updated version of PIC code the scattering of particles by the anode
mesh was added via consideration of electron-matter interaction. (We will
refer to the new simulation model as the scattering model.) The anode mesh was
described by three parameters: the geometric transparency of the mesh, the
anode material, and anode thickness $d$. The geometric transparency $\eta$ was
introduced in terms of probability ($1-\eta$) for a particle to hit the anode
mesh, which was supposed to be a steel foil of thickness $d$ with holes. The
ratio of the sum of hole areas to the total foil area defined the geometric
transparency $\eta$.

Figure 6: Electric field in the VC region and electric field spectrum.
Absorption model.

Figure 7: Phase portrait of an electron beam in the absorption model. The
longitudinal coordinates and particle velocities are shown along the abscissa
in the units of the cathode-anode gap and along the ordinate in units of speed
of light, respectively.

Thus, a particle passing through the anode is scattered with probability
$1-\eta$ rather than $\eta$. Particle scattering by the anode is calculated
exactly in the same way as scattering by a solid foil222The possible complex
structure of the anode is not taken into account. For example, the mesh may be
made of round wires, which cannot be described in terms of a one-dimensional
model.. The scattered (unabsorbed) particle enters either the cathode-anode
gap or the drift space with a random component of the transverse momentum and
with a reduced longitudinal velocity due to ionization energy loss. In this
case, time of particle motion inside the anode material is neglected. Since
the characteristic thickness $d$ is much smaller than the characteristic value
of the cathode-anode gap, we neglect time, which particle spends inside the
metal.

Figures 8 and 9 show the electric field oscillations, spectra and the phase
portraits in the case of steel anode of different thicknesses (thicknesses are
indicated on the plot in units of the mean free path $s_{0}$). The geometric
transparency of anode mesh is $\eta=0.7$. Simulation of vircators showed that
electron scattering by the anode mesh leads to the formation of a cloud of
low-energy electrons with a large spread in velocities. These electrons appear
due to loss of longitudinal momentum after scattering by the anode mesh. Due
to deceleration, these electrons do not have enough energy to reach the either
real or virtual cathode. As a result, they oscillate in the potential well
between the cathode and the virtual cathode until complete absorption.
Scattering by the anode leads to a random change in phase of electron
oscillations. This change prevents particles from generating coherent
oscillations. In addition, an electron cloud partially shields the anode and
blocks the vacuum diode. As a consequence, the amount of electrons, which
could participate in collective oscillations, appears smaller. As a result,
the amplitude of oscillations in a vircator decreases.

The vircator operation can be improved if the thickness of the partitions in
the anode mesh is chosen to be close to the electron path length in the anode
material. In this case, the number of scattered electrons passing through the
anode could be significantly reduced. As a consequence, fewer particles in the
cloud would affect oscillation amplitude decreasing.

Figure 8: Electric field and its spectrum in the VC region.

Figure 9: Phase portraits of the electron beam at different thicknesses of the
steel anode. The phase portraits shown on the left and on the right correspond
to the thickness $0.8s_{0}$ and $0.2s_{0}$, respectively. The longitudinal
coordinates and particle velocities are shown along the abscissa in the units
of the cathode-anode gap and along the ordinate in units of speed of light,
respectively

##  IV onclusion

Multiple scattering of a relativistic electron beam by an electrodynamic
structure accompanied by ionization and radiation energy losses significantly
influences on operation of a high-current virtual cathode oscillator. In a
vircator, the scattering of a high-current relativistic electron beam by the
anode mesh leads to formation of an electron cloud near the anode. The cloud
particles possess low energy and large spread in velocities caused by multiple
scattering and ionization losses. The electrons captured by the cloud do not
participate in the oscillations of the virtual cathode and partially block a
vircator. As a result, the amplitude of the electric field oscillations is
reduced. In order to increase the oscillation amplitude, the thickness of
wires forming the anode mesh should be equal to the mean free path for
electrons in the mesh material. Use of anode mesh corresponding the above
requirements enables to minimize the number of scattered electrons passing
through the anode.

We applied the approach for modeling of multiple scattering described by
Fernandez1993  in the simulation code developed by us to calculate the
influence of electron ionization losses and multiple scattering by the anode
mesh on particles’ dynamics in a vircator. The above approach significantly
reduces calculation time as compared with the most popular concept based on
the Goudsmit-Sanderson distribution.

## References

  * (1) J.M. Creedon, Multiple scattering of electrons in the reflex triode, J. Appl. Phys. 1990. Vol. 68. No. 11. P. 5494.
  * (2) E.J. Williams, Concerning the scattering of fast electrons and cosmic-ray particles, Proc. R. Soc. Lond. A. 1939. Vol. 169, P. 531–572.
  * (3) S. Gouldsmit and J.L. Saunderson, Multiple scattering of electrons, Phys. Rev., 1940. Vol. 57, P. 244–29.
  * (4) G. Moliere, Theorie der Streuung schneller geladener Teilchen II, Mehrfach-und Vielfachstreuung, Z. Naturforsch. Ba, 78 (1948).
  * (5) H.A. Bethe, Moliere’s theory of multiple scattering, Phys. Rev. 1953. Vol. 89. P. 1256–1266.
  * (6) B.P. Nigam, M.K. Sundaresan, and Ta-You Wu, Theory of multiple scattering: second Born approximation and corrections to Moliere’s work, 1959. Vol. 115, P. 491–502.
  * (7) J.B. Marion and B.A. Zimmerman, Multiple scattering of charged particles, Nucl. Instrum. Methods, 1967. Vol. 51., P. 93–101.
  * (8) V.G. Voinov, I.Ya. Chasnikov. Multiple scattering in nuclear emulsions. — Alma-Ata: Nauka, 1969 [in Russian].
  * (9) R. Mayol, F. Salvat, Total and transport cross sections for elastic scattering of electrons in atoms, Atomic Data and Nuclear Data Tables 65 (1997).
  * (10) N. Bohr, On the theory of the decrease of velocity of moving electrified particles on passing through matter, Phil. Mag. 1913. Vol. 25., P. 10–31.
  * (11) H. Bethe, Bremsformel ffir Elektronen relativistischer Geschwindigkeit, Zeits. f. Physik, 1932, 76, 293.
  * (12) F. Bloch, Zur Bremsung rasch bewegter Te2chenbeim Durchgarzg durch Materie, Ann. d Physik 16, 285 (1933).
  * (13) A. Bethe and J. Ashkin, Experimenal Nuclear Physics, L. Segre, ed., Vol. 1, p. 276 (John Wiley and Sons, Inc., New York, 1953)
  * (14) F. Rohrlich and B.C. Carlson, Positron-electron differences in energy loss and multiple scattering, Phys. Rev., 1954. Vol. 93., P. 38.
  * (15) H. Bethe and W. Heitler, Proc. Roy. Soc. (London) A146 (1934) 83.
  * (16) R.H. Pratt et al., Bremsstrahlung energy spectra from electrons of kinetic energy 1 keV $\leq T_{1}\leq 2000$ keV incident on neutral atoms $2\leq Z\leq 92$, Atomic Data and Nuclear Data Tables 20 (1977).
  * (17) D.F. Hebbard and P.R. Wilson, The effect of multiple scattering on electron energy loss distributions, 1955.
  * (18) J.F. Perkins, Monte-Carlo calculation of transport of fast electrons, Phys. Rev, 1962. Vol. 126. P. 1781–1784.
  * (19) M.J. Berger, Monte-Carlo calculation of the penetration and diffusion of fast charged particles, 1963.
  * (20) M.J. Berger, Transmission and reflection of electrons by aluminum foils, 1963. NBS TN187.
  * (21) M.J. Berger and S.M. Seltzer, Tables of energy losses and ranges of electrons and positrons, NASA SP-3012 (1964).
  * (22) S.M. Seltzer and M.J. Berger, Transmission and reflection of electrons by foils, Nucl. Instrum. Methods, 1974. Vol. 119, P. 157–179.
  * (23) S.M. Seltzer, Electron-photon Monte-Carlo calculations: the ETRAN Code, Appl. Radiat. Isot. 1991. Vol. 42, No. 10. P. 917–941, 1991.
  * (24) Monte-Carlo transport of electrons and photons, eds. Th. Jenkins, W.R. Nelson, and A. Rindi, Erice, Italy, 1987.
  * (25) J.M. Fernandez-Varea, R. Mayol, J. Baro, and F. Salvat, On the theory and simulation of multiple elastic scattering of electrons, Nucl. Instrum. Methods in Phys. Res. B73 (1993) 447–473.
  * (26) T. Tabata, S. Okabe, An empirical equation for the backscattering coefficients of electrons, Nucl. Instrum. Methods. 1971. Vol. 94, P. 509-513.
  * (27) R. Ito, P. Andreo, T. Tabata, Reflection ratios of electrons and photons from solids, Bull. of University of Osaka Prefecture, Series A, Vol. 41, Issue 2, 1993, 69–76.
  * (28) H.W. Lewis, Multiple scattering in an infinite medium, Phys. Rev. 1950. Vol. 78. P. 526.


In the description of the properties of the Moore–Penrose pseudoinverse above, we noted that $(F^+)^\dag = (F^\dag)^+$. In fact, this property generalizes even further, motivating the definition of the induced inner product on $\Hom$ spaces.

Suppose we have inner product spaces $\Uprod$, $\Vprod$, and $\Wprod$, together with $\Hom(U,W)$ and $\Hom(V,W)$ with their induced (Frobenius) inner products, and a linear map $F \co U \rightarrow V$. Then $(F^*)^+ = (F^+)^*$. Both are maps $\Hom(U,W) \rightarrow \Hom(V,W)$.

The proof is an exercise in checking that $(F^+)^*$ satisfies the Moore–Penrose conditions for a pseudoinverse of $F^* \co \Hom(V,W) \rightarrow \Hom(U,W)$.

\begin{equation*}
(F^+)^* F^* (F^+)^* = (F^+ F F^+)^* = (F^+)^* \quad\text{and}\quad
F^* (F^+)^* F^* = (F F^+ F)^* = F^*.
\end{equation*}

Referring to induced inner product for the Frobenius inner product, and using standard properties of $\dag$ and $*$[Here we've used several properties of $\dag$ and $\star$: $(FG)^\dag = G^\dag F^\dag$ in the second equality; $F^{*\dag} = F^{\dag*}$ and $(F^+F)^\dag = F^+ F$ in the third, and $(FG)^* = G^* F^*$ in the last one. ]

\[
\Frp{F^*(F^+)^* A}{B} = \tr (F^*(F^+)^*A)^\dag B = \tr A^\dag (F^+ F)^{*\dag} B = \tr A^\dag (F^+F)^* B = \Frp{A}{F^* (F^+)^* B},
\]

so $(F^* (F^+)^*)^\dag = (F^* (F^+)^*)$. Similarly,

\begin{multline*}
\Frp{(F^+)^* F^* A}{B} = \tr ((F^+)^* F^* A)^\dag B = \tr A^\dag ((F^+)^* F^*)^\dag B \\
=\tr A^\dag (F F^+)^{*\dag} B = \tr A^\dag (F F^+)^* B
%= \tr A^\dag (F^+)^* F^* B
= \Frp{A}{(F^+)^* F^* B},
\end{multline*}

so $((F^+)^* F^*)^\dag = (F^+)^* F^*$.


# A Low-mass Cold and Quiescent Core Population in a Massive Star Protocluster

Shanghuo Li  Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro,
Yuseong-gu, Daejeon 34055, Republic of Korea  Shanghai Astronomical
Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030,
People’s Republic of China  Center for Astrophysics $|$ Harvard & Smithsonian,
60 Garden Street, Cambridge, MA 02138, USA  University of Chinese Academy of
Sciences, 19A Yuquanlu, Beijing 100049, People’s Republic of China  Xing Lu
National Astronomical Observatory of Japan, National Institutes of Natural
Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Qizhou Zhang  Center
for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA
02138, USA  Chang-Won Lee  Korea Astronomy and Space Science Institute, 776
Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea  University of
Science and Technology, 217 Gajeong-ro, Yuseong-gu, Daejeon 34113, Republic of
Korea  Patricio Sanhueza  National Astronomical Observatory of Japan, National
Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Department of Astronomical Science, SOKENDAI (The Graduate University for
Advanced Studies), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Henrik Beuther
Max Planck Institute for Astronomy, Konigstuhl 17, 69117 Heidelberg, Germany
Izaskun, Jiménez-Serra  Centro de Astrobiología (CSIC-INTA), Ctra. de Torrejón
a Ajalvir, Km. 4, Torrejón de Ardoz, 28850 Madrid, Spain  Keping Qiu  School
of Astronomy and Space Science, Nanjing University, 163 Xianlin Avenue,
Nanjing 210023, People’s Republic of China  Aina Palau  Instituto de
Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, P.O.
Box 3-72, 58090, Morelia, Michoacán, México  Siyi Feng  National Astronomical
Observatories, Chinese Academy of Sciences, Beijing 100101, People’s Republic
of China  Academia Sinica Institute of Astronomy and Astrophysics, No. 1,
Section 4, Roosevelt Road, Taipei 10617, Taiwan, Republic of China  National
Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo, 181-8588,
Japan  Thushara Pillai  Institute for Astrophysical Research, Boston
University, 725 Commonwealth Ave, Boston, MA 02215, USA  Kee-Tae Kim  Korea
Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon
34055, Republic of Korea  University of Science and Technology, 217 Gajeong-
ro, Yuseong-gu, Daejeon 34113, Republic of Korea  Hong-Li Liu  Department of
Astronomy, Yunnan University, Kunming, 650091, People’s Republic of China
Josep Miquel. Girart  Institut de Ciències de l’Espai (IEEC-CSIC), Campus UAB,
Carrer de Can Magrans s/n, 08193 Cerdanyola del Vallès, Catalonia, Spain  Tie
Liu  Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai 200030, People’s Republic of China  Junzhi Wang  Shanghai
Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road,
Shanghai 200030, People’s Republic of China  Ke Wang  Kavli Institute for
Astronomy and Astrophysics, Peking University, 5 Yiheyuan Road, Haidian
District, Beijing 100871, People’s Republic of China  Hauyu Baobab Liu
Academia Sinica Institute of Astronomy and Astrophysics, 11F of AS/NTU
Astronomy-Mathematics Building, No.1, Sec. 4, Roosevelt Rd, Taipei 10617,
Taiwan, Republic of China  Howard A. Smith  Center for Astrophysics $|$
Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Di Li
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, People’s Republic of China  NAOC-UKZN Computational Astrophysics
Centre, University of KwaZulu-Natal, Durban 4000, South Africa  Jeong-Eun Lee
School of Space Research, Kyung Hee University, 1732, Deogyeong-Daero,
Giheung-gu, Yongin-shi, Gyunggi-do 17104, Republic of Korea  Fei Li  School of
Astronomy and Space Science, Nanjing University, 163 Xianlin Avenue, Nanjing
210023, People’s Republic of China  Juan Li  Shanghai Astronomical
Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030,
People’s Republic of China  Shinyoung Kim  Korea Astronomy and Space Science
Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea
University of Science and Technology, 217 Gajeong-ro, Yuseong-gu, Daejeon
34113, Republic of Korea  Nannan Yue  Kavli Institute for Astronomy and
Astrophysics, Peking University, 5 Yiheyuan Road, Haidian District, Beijing
100871, People’s Republic of China  Shaye Strom  Center for Astrophysics $|$
Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

###### Abstract

Pre-stellar cores represent the initial conditions of star formation. Although
these initial conditions in nearby low-mass star-forming regions have been
investigated in detail, such initial conditions remain vastly unexplored for
massive star-forming regions. We report the detection of a cluster of low-mass
starless and pre-stellar core candidates in a massive star protocluster
forming cloud, NGC 6334S. With the ALMA observations at a $\sim$0.02 pc
spatial resolution, we identified 17 low-mass starless core candidates that do
not show any evidence of protostellar activity. These candidates present small
velocity dispersions, high fractional abundances of NH2D, high NH3 deuterium
fractionations, and are completely dark in the infrared wavelengths from 3.6
up to 70 $\mu$m. Turbulence is significantly dissipated and the gas kinematics
are dominated by thermal motions toward these candidates. Nine out of the 17
cores are gravitationally bound, and therefore are identified as pre-stellar
core candidates. The embedded cores of NGC 6334S show a wide diversity in
masses and evolutionary stages.

Unified Astronomy Thesaurus concepts: Early-type stars (430), Infrared dark
clouds (787), Star forming regions (1565), Star formation (1569), Interstellar
medium (847), Interstellar line emission (844), Protoclusters (1297)

††facilities: ALMA, JVLA, Herschel.††software:  CASA (McMullin et al., 2007),
APLpy (Robitaille & Bressert, 2012), Astropy (Astropy Collaboration et al.,
2013), Matplotlib (Hunter, 2007), PySpecKit (Ginsburg & Mirocha, 2011).

##  1 Introduction

Stars form in the high-density parts of molecular clouds, i.e., dense cores. A
dense core prior to the protostellar phase is called a starless core, and a
gravitationally bound/unstable starless core is referred to as a pre-stellar
core (Ward-Thompson et al., 1994). The pre-stellar phase is considered as the
starting point in the star formation process (Bergin & Tafalla, 2007; Caselli
& Ceccarelli, 2012). Pre-stellar cores can form single or multiple stellar
systems under the combined effect of gravity, magnetic fields, and turbulence
(Pineda et al., 2015). Understanding the properties of pre-stellar cores is
critical to characterize the initial conditions of cluster formation. The pre-
stellar cores in low-mass star-forming regions have been intensively studied
toward nearby molecular clouds, e.g., Perseus, Ophiuchus, Chamaeleon, Serpens,
and Taurus (Lee et al., 1999, 2001; Enoch et al., 2008; Schnee et al., 2010;
Dunham et al., 2016; Kirk et al., 2017; Tokuda et al., 2020). However, the
studies of pre-stellar cores in massive cluster-forming regions are still
limited by the low number of statistics available (Nony et al., 2018; Li et
al., 2019; Sanhueza et al., 2019; Sahu et al., 2021).

The massive infrared dark cloud NGC 6334S (also known as IRDC G350.56+0.44) is
located at the southwestern end of the NGC 6334 molecular cloud complex that
is a nearby (1.3 kpc; Chibueze et al., 2014) young and massive ‘mini-
starburst’ star-forming region (Willis et al., 2013). With a mass of 1.3
$\times\,10^{3}\,M_{\odot}$ (Li et al., 2020a), comparable to the clumps with
embedded massive protostars and protoclusters in the complex, NGC 6334S has
the potential to form a cluster encompassing both low- ($<$ 2 $M_{\odot}$) and
high-mass ($>$ 8 $M_{\odot}$) stars. NGC 6334S provides an ideal laboratory to
search for and study the early stages (e.g., pre-stellar phase) of star
formation in the high-mass regime.

In order to investigate massive star and cluster formation, we performed
Atacama Large Millimeter/submillimeter Array (ALMA) and Karl G. Jansky Very
Large Array (JVLA) observations of NGC 6334S. A study of 49 continuum dense
cores (hereafter continuum cores) revealed by the 3 mm wavelength continuum
image was presented in Li et al. (2020a), in which we reveal that the
nonthermal motions are predominantly subsonic and transonic in both clump
scale and continuum core scale. Here we use molecular lines to study the
starless/pre-stellar cores.

##  2 Observations

###  2.1 ALMA Observations

We carried out a 55–pointing mosaic of NGC 6334S using the ALMA 12 m array in
2017 March (ID: 2016.1.00951.S). Two 234.4 MHz width spectral windows with a
spectral resolution of 61 kHz ($\sim$0.21 km s-1 at 86 GHz) were configured to
cover the H13CO+ (1–0, 86.754 GHz) and ortho-NH2D ($1_{11}-1_{01}$ , 85.926
GHz) lines, respectively. Three additional 1.875 GHz wide spectral windows at
88.5 GHz, 98.5 GHz, and 100.3 GHz with a coarse spectral resolution (3.0 – 3.3
km s-1) were employed to obtain broad band continuum and rotational
transitions of several molecular species (e.g., HCO+ 1–0, HCN 1–0, CS 2–1,
HNCO $4_{0,4}-3_{0,3}$, H15NC 1–0, CH3OH $5_{1,4}-4_{1,3}$, SO $2_{2}-1_{1}$,
HC3N 11–10). The maximum recoverable scales (MRS) is $\sim$25′′ in the ALMA
data. The details of the observations can be found in Li et al. (2020a).

Data calibration and imaging were performed using CASA 4.7.0 (McMullin et al.,
2007). We used Briggs’ robust weighting of 0.5 to the visibilities for both
the continuum and lines, which results in a synthesized beam of 3.′′6 $\times$
2.′′4 with a position angle (P.A.) of 81∘ and 4.′′1 $\times$ 2.′′8 (P.A. =
83∘) for continuum and line images, respectively. The achieved 1$\sigma$ rms
noise level is about 6 mJy beam-1 per 0.21 km s-1 for the line images and 30
$\mu$Jy beam-1 for the continuum image. All images shown in the paper are
prior to primary beam correction, while all measured fluxes are corrected for
the primary beam attenuation.

###  2.2 JVLA Observations

We carried out a 4-pointing mosaic of the central region of NGC 6334S using
the JVLA C-configuration in 2014 August (ID:14A-241). The NH3 (1, 1) through
(5, 5) metastable inversion transitions and H2O maser were simultaneously
covered in this observation. CASA versions 4.7.0 and 5.1.1 were used to
perform data calibration and imaging. An elliptical Gaussian with a FWHM of
6′′$\times$ 3′′(P.A. = 0∘) was used to taper the visibilities, in order to
increase the signal-to-noise ratio (S/N). This yields a synthesized beam size
of about 10′′$\times$ 5′′(P.A. = 26∘) with a 1$\sigma$ rms noise of 9 mJy
beam-1 per 0.2 km s-1 for the NH3 lines. More details on the observations are
presented in Li et al. (2020a).

Figure 1:  Panel a: the gray scale in the background shows the NH2D velocity-
integrated intensity ($W_{\rm NH_{2}D}$) image. Cyan open triangles show the
NH2D cores. Yellow open squares present 49 continuum cores identified by ALMA
3 mm continuum image (Li et al., 2020a). Purple cross ‘x’ and blue plus ‘+’
symbols are 25 Class I and 58 Class II YSOs (Willis et al., 2013),
respectively. The beam size of the NH2D image is shown on the bottom left of
the panel. Dashed purple contour shows the area mosaicked with VLA. Panel b:
the core-averaged spectra of NH2D ($1_{11}-1_{01}$) and H13CO+ (1-0) for
continuum core #2 (black solid line) and NH2D core M1 (blue solid line). Panel
c: the core-averaged spectra of three wide (1.875 GHz) spectral windows for
continuum core #2 (black solid line) and NH2D core M1 (blue solid line),
respectively.  Figure 2:  The 3 mm continuum image (black contours) overlaid
on the velocity-integrated intensity of NH2D for the candidates that are
gravitationally bound with $\alpha_{\rm vir}<$ 2\. The integrated velocity
range is -11.5 and 5.6 km s-1. Black contours are (5, 7, 9, 11, 15, 17, 20,
24, 29)$\times\sigma$, where $\sigma$ = 30 $\mu$Jy beam-1 is the rms noise
level of the continuum image. The yellow ellipses show the identified NH2D
cores. The green open stars indicate continuum cores. The beam size of NH2D
image is shown on the bottom right of M10 panel.

##  3 Results and Analysis

###  3.1 Cold and Quiescent Cores

NH2D ($1_{11}-1_{01}$, critical density $n_{\rm cr}\sim$ 105 cm-3) is a good
tracer of cold and dense molecular gas (Crapsi et al., 2007; Sanhueza et al.,
2013), which can survive in the gas phase in the dense interior region of pre-
stellar cores (Harju et al., 2017). In our ALMA data, the NH2D
($1_{11}-1_{01}$) line emission is in general well correlated with the 3 mm
continuum emission, but there are some exceptions; there are 17 bright compact
structures in the NH2D emission that are associated with weak (3$\sigma$ –
11$\sigma$) or no continuum emission ($<3\sigma$) at all (Figure 1 and 2), and
with no young stellar object (YSO) counterparts (e.g., Class I/II; Willis et
al., 2013). There is no reliable Class 0 catalogue in this region. The details
on the identification of NH2D compact structures are summarized in Appendix A.
We refer to these NH2D compact structures as NH2D cores, naming those in M1,
M2, M3 … in order of descending NH2D velocity-integrated intensity. They have
diameters ranging from 0.018 to 0.04 pc (see Table 1).

Among the NH2D cores, M1 is the most prominent pre-stellar candidate with high
signal-to-noise (S/N) ratio and a relatively isolated environment that
eliminates the contamination from the outflows driven by the continuum cores.
The following analysis and discussion will therefore use M1 as a showcase,
while the physical parameters of the remaining NH2D cores will only be
included in statistical analyses.

Using the HCO+, CS, HCN, SO, HNCO, and H2O maser lines (see Figure 1), we have
searched for signatures of protostellar activity, such as
bipolar/monopolar/multipolar outflow activity and H2O maser emission. All the
NH2D cores show no signs of protostellar activity. Although we currently
cannot rule out the existence of unresolved or weak outflows due to the coarse
spectral resolution and the lack of better low-mass outflow tracers such as
CO, these NH2D cores appear to be starless with the available evidence. In
addition, the NH2D cores do not show any appreciable ($<3\sigma$) continuum
emission in the infrared wavelengths from 3.6 up to 70 $\mu{}$m, and therefore
they are completely dark in these infrared images. The upper limit on the
internal luminosity ($L_{\rm int}$) of NH2D cores is estimated to be
$\sim$1.26 $L_{\odot}$ using Herschel/PACS 70 $\mu{}$m data (Tigé et al.,
2017, see Appendix C for detailed derivation of internal luminosity). These
results suggest that these NH2D cores are cold and quiescent in terms of star-
forming activity, though some undetected very faint embedded low-mass
protostars cannot be fully ruled out (Dunham et al., 2008; Evans et al., 2009;
Offner & McKee, 2011; Fischer et al., 2017). An example of the comparison
between NH2D core (M1) and continuum core (#2) is shown in Figure 1, which
includes all of the observed spectral windows. Continuum core #2 is a
protostellar core as evidenced by the clear molecular outflows in the CS, HCN,
and HCO+ line emission, and it is one of the most chemically rich among the
continuum cores. Nonetheless, only a few lines are detected toward continuum
core #2. This indicates that the continuum cores #2 is at a stage prior to the
hot cores/corinos. We note that M1 shows even less chemical complexity than
continuum core #2, with fewer lines, molecular species, and narrower line
profiles. This suggests that the evolutionary stage of the NH2D cores are
likely earlier than the continuum cores. On the other hand, we cannot
completely rule out the possibility that some of continuum cores are
relatively evolved pre-stellar cores, because of the lack of high sensitivity
and high spectral resolution outflow tracers (e.g., CO, SiO) to distinguish
the evolutionary stages (Li et al., 2020b). In this study, we will focus only
on the NH2D cores.

Figure 3: Panel a: NH2D spectrum (black curve) extracted from single pixels
toward M1. Red curve shows the best result from the hyperfine fitting that
includes NH2D 36 hyperfine components (Ginsburg & Mirocha, 2011; Daniel et
al., 2016), with the best-fit parameters displayed in the panel. Panel b: the
NH2D integrated intensity (black contours) overlaid on its $\sigma_{\rm obs}$
toward M1. The contours are (3, 5, 10, 20, 30, 40, 50)$\times\sigma$, with
$\sigma$ = 0.011 Jy beam-1 km s-1. The beam size is shown on the lower left of
the panel. Panel c: the annularly averaged observed velocity dispersions
($\sigma_{\rm obs}$) as a functions of radial distance from the center (given
in Table 1) of M1. The error bars show the statistical standard deviation at
each bin.

###  3.2 Small Velocity Dispersions

The NH2D line cube is fitted pixel by pixel, following the fitting processes
described in Li et al. (2020a). In the hyperfine fits, the uncertainties of
the LSR velocities and line widths are small, while those of the excitation
temperature and optical depth are relatively large (Figure 3). Figure 4 shows
the observed velocity dispersions ($\sigma_{\rm obs}$) distribution for each
NH2D core. All the NH2D cores show small velocity dispersion
($\langle\sigma_{\rm obs}\rangle\sim$ 0.16 km s-1) and the majority of them
are smaller than the sound speed, $c_{\rm s}$(10 K) = 0.19 km s-1. The
observed velocity dispersions may be regarded as upper limits, given the
limited spectral resolution, the blending of velocity components along the
line of sight and opacity broadening effects. For the NH2D cores, the line
widths of detected lines are smaller than those of the majority of continuum
cores (see Figure 1 and 4).

Among all the NH2D cores, M1 shows the strongest NH2D emission and the
smallest velocity dispersion. For M1, pixel by pixel fit to the NH2D line cube
yields a small mean velocity dispersion $\langle\sigma_{\rm
obs}\rangle\approx$ 0.107 km s-1 (Figure 4), corresponding to an intrinsic
observed velocity dispersion of $\sigma_{\rm obs,int}$ = 0.06 km s-1 after
removing the smearing effect due to the channel width using $\sigma_{\rm
obs,int}=\sqrt{\sigma_{\rm obs}^{2}-\sigma_{\rm ch}^{2}}$, where $\sigma_{\rm
ch}=\triangle_{\rm ch}/(2\sqrt{2\;\rm ln2})$ = 0.089 km s-1 is the channel
width. For the extremely narrow line emission, a higher spectral resolution
data is required to accurately characterize the actual intrinsic velocity
dispersion. Considering the $\langle\sigma_{\rm obs}\rangle$ is 1.2 times
larger than the $\sigma_{\rm ch}$, the aforementioned derived $\sigma_{\rm
obs,int}$ can be regarded as a reasonable approximation for the following
analysis (see Appendix C in Li et al., 2020a). Assuming the observed
linewidths of the NH2D line are caused only by thermal broadening, the gas
temperature is found to be about 7 K, $\sigma_{\rm th,NH_{2}D}$(7 K)111 The
molecular thermal velocity dispersion can be calculated by $\sigma_{\rm
th}\;=\;\sqrt{(k_{\rm B}T)/(\mu m_{\rm H})}\;=\;9.08\times
10^{-2}\rm\;km\;s^{-1}(\frac{T}{K})^{0.5}\;\mu^{-0.5}$, where
$\mu\;=\;m/m_{\rm H}$ is the molecular weight, $m$ is the molecular mass and
$m_{\rm H}$ is the proton mass. = $\sigma_{\rm obs,int}$ = 0.06 km s-1. This
indicates that the gas kinematics is dominated by thermal motions with a
kinetic temperature that is 7 K at most. In addition to M1, there are four
NH2D cores (M9, M12, M13 and M15) that also show extremely small velocity
dispersions comparable to the NH2D thermal velocity dispersion in the
temperature range of 7 – 10 K. The velocity dispersions of remaining cores are
relatively higher than the NH2D thermal velocity dispersion at a temperature
of 10 K ($\sigma_{\rm th,NH_{2}D}\sim$ 0.07 km s-1), while they are much
smaller than the sound speed (Figure 4). This indicates that the gas
kinematics are dominated by thermal motions toward these NH2D cores.

Figure 3 shows the annularly averaged observed velocity dispersion as a
function of the radial distance ($R_{\rm dist}$) from the center of M1. The
observed velocity dispersion appears to increase with increasing $R_{\rm
dist}$ (see Figure 3), which prevails in pre-stellar and starless cores
(Crapsi et al., 2005; Auddy et al., 2019). The increase of $\sigma_{\rm obs}$
along Rdist is also detected in the other 9 cores (M4, M6, M7, M8, M9, M10,
M13, M15, M17), although the annularly averaged $\sigma_{\rm obs}$ starts to
suffer from low S/N toward the outer edges of the cores. This increasing trend
of $\sigma_{\rm obs}$ along Rdist indicates that the turbulence is
significantly dissipated toward the center of the core. The results indicate
that the turbulence dissipation is commonly seen toward these NH2D cores.
Alternatively, in the larger-scale collapse scenario, if the pre-stellar cores
undergo accretion, the infall speed decreases towards the center, and thus the
line width will be smaller for more central regions (Gómez et al., 2021).
Although no clear infall signatures are detected in M1 using coarser spectral
resolution (3.0 – 3.3 km s-1) optically thick lines (e.g., HCO+ 1-0, HCN 1-0
and CS 2-1), this possibility cannot be ruled out and should be explored using
higher spectral resolution data and other tracers (e.g., HCN 5-4/4-3/3-2 and
HCO+ 5-4/3-2; Chira et al., 2014).

Overall, the observed small velocity dispersions and their spatial
distributions resemble those observed in the well known pre-stellar cores in
Ophiuchus/H-MM1 and L 1544 (Crapsi et al., 2007; Harju et al., 2017; Punanova
et al., 2018). This suggests that the identified NH2D cores may have not been
affected by the star formation activities yet.

Figure 4:  Panel a to c: violin plots of the $\sigma_{\rm obs}$, $N_{\rm
NH_{2}D}$, and $\chi(\rm NH_{2}D)$ distributions for all NH2D cores. The shape
of each distribution shows the probability density of the data smoothed by a
kernel density estimator. The blue bars from the top to bottom represent the
maximum, mean, and minimum values, respectively. Panel a: the red solid line
is the sound speed at a temperature of 10 K. The blue solid line is the mean
observed velocity dispersion of the continuum cores, $\langle\sigma_{\rm
obs,cont}\rangle$ = 0.23 $\pm$ 0.08 km s-1. The shaded gray region shows the
spectral resolution limit of the NH2D line. Panel b: the red and blue solid
lines are observed $\rm N(NH_{2}D)$ in the pre-stellar core L1544 (1.75
$\times\,10^{14}$ cm-2; Crapsi et al., 2007) and Ophiuchus/H-MM1 (1.1
$\times\,10^{14}$ cm-2; Harju et al., 2017). Panel c: the blue and red solid
lines are observed abundance $X(\rm NH_{2}D)$ in the pre-stellar cores
Ophiuchus/H-MM1 (2 $\times\,10^{-9}$, Harju et al., 2017) and L1544 (2
$\times\,10^{-10}$, Crapsi et al., 2007), respectively.

###  3.3 Gas Masses and Dynamical States

Assuming a typical temperature ($T_{\rm dust}$ = 10 K) of pre-stellar cores
(see also $\S$ 3.4), the gas masses ($M_{\rm gas}$) of the NH2D cores are
estimated from continuum emission in two different approaches depending on
their intensities (see Appendix B for detailed derivation of gas mass). First,
the continuum emission is used to compute the gas mass if the continuum peak
emission is higher than $3\sigma$. Second, a 3$\sigma$ mass sensitivity of
0.13 $M_{\odot}$ is used as an upper limit for the core mass if the continuum
peak emission is $\leqslant 3\sigma$. The derived cores masses are between $<$
0.13 and 0.87 $M_{\odot}$, with a mean value of 0.45 $M_{\odot}$ (Table 1).
Following the same procedure, the core-averaged H2 column densities ($N_{\rm
H_{2}}$) are derived to be between $<$ 1.1 $\times$ 1022 and 3.0 $\times$ 1022
cm-2 (Appendix B), with a mean value of 2.0 $\times$ 1022 cm-2.

In order to evaluate the gravitationally bound of the NH2D cores, we estimated
the virial masses ($M_{\rm vir}$) that range from 0.43 to 2.54 $M_{\odot}$
(see Appendix B for derivation of virial mass). The derived virial parameters,
$\alpha_{\rm vir}$ = $M_{\rm vir}/M_{\rm gas}$, are between 0.93 and 6.03.
Nine out of the 17 NH2D cores are gravitationally bound with virial parameters
$<$ 2, and the remaining eight NH2D cores are gravitationally unbound with
$\alpha_{\rm vir}>$ 2 if external pressure is ignored (Kauffmann et al.,
2013).

###  3.4 NH2D Abundances and NH3 Deuterium Fractionation

Assuming a constant temperature of 10 K and the ortho-to-para ratio of 3
(Harju et al., 2017), the estimated NH2D column densities ($N(\rm NH_{2}D)$)
range from 7.3 $\times$ 1012 cm-2 to 3.4 $\times$ 1015 cm-2, with a mean value
of 5.2 $\times$ 1014 cm-2 (see Figure 4). The derived abundance $X(\rm
NH_{2}D)$ = $N(\rm NH_{2}D)$/$N(\rm H_{2})$ toward the NH2D cores range from
1.8 $\times$ 10-9 to 2.0 $\times$ 10-7, with a mean value of 2.5 $\times$
10-8.

The JVLA NH3 observations only cover the central region of NGC 6334S (see
Figure 1). There are ten NH2D cores located in the field of the NH3
observations. Limited by the relatively low S/N, both NH3 (1,1) and (2,2)
transitions are detected in M4 and M10 (S/N $\sim$ 4–5), while only NH3 (1,1)
is marginally detected in M13 and undetected in the other seven NH2D cores.
Using PySpecKit we fit the core-averaged spectrum, and we derive a kinetic
temperature of $T_{\rm k}$ = 12.7 $\pm$ 3.1 K and a column density of $N_{\rm
NH_{3}}$ = $(4.4\,\pm\,2.5)\times 10^{14}$ cm-2 for M4, and the derived
parameters are $T_{\rm k}$ = 12.1 $\pm$ 3.2 K and $N_{\rm NH_{3}}$ = $(8.9\pm
5.0)\times 10^{14}$ cm-2 for M10, assuming the ortho-to-para ratio of 1 for
NH3 (Harju et al., 2017). The column density ratio NH2D-to-NH3 ($D_{\rm
NH_{3}}$, also known as NH3 deuterium fractionation) is found to be 0.25 and
0.11 for M4 and M10, respectively. A 3$\sigma$ integrated intensity is used to
estimate the upper limit column density (3.8 $\times$ 1014 cm-2) for the other
eight cores assuming an excitation temperature of 10 K. The resulting lower
limits $D_{\rm NH_{3}}$ range from 0.11 to 0.39 (see Table 1), which are close
to those in the pre-stellar core L 1544 (0.5$\pm$0.2; Crapsi et al., 2007) and
Ophiuchus/H-MM1 (0.45$\pm$0.09; Harju et al., 2017).

##  4 Discussion

As shown in Figure 4, the derived NH2D column densities and fractional
abundances are higher than the values of Ophiuchus/H-MM1
($1.1\,\times\,10^{14}$ cm-2 and $2\,\times\,10^{-9}$; Harju et al., 2017) and
L1544 ($1.8\,\times\,10^{14}$ cm-2 and $2\,\times\,10^{-10}$; Crapsi et al.,
2007). The high abundances of NH2D could benefit from the significant
depletion of CO, which occurs faster at cold and dense environments (e.g.,
starless or pre-stellar cores; the detail chemical processes related to the
NH2D formation can be found in Sipilä et al., 2015; Harju et al., 2017).

The small velocity dispersion, small Mach number (Table 1), and positive
$\sigma_{\rm obs}$–$R_{\rm dist}$ relation indicate that the turbulence is
likely dissipated toward these NH2D cores. In addition, it implies that star
formation has not yet started; outflows motions would widen the lines as can
be seen in the continuum cores (Figure 1). Therefore, these results in
conjunction with the high NH3 deuterium fractionation, and the lack of
embedded YSOs and 70 $\mu{}$m emission indicate that the identified NH2D cores
are still in a starless phase. Nine out of 17 starless cores are
gravitationally bound with $\alpha_{\rm vir}<$ 2, therefore are identified as
excellent pre-stellar core candidates. The dust temperature may drop to
$\sim$5.5 K toward the core center (e.g., L1544; Crapsi et al., 2007). The
derived gas mass can increase by a factor of 2.25 if a temperature of 5.5 K is
assumed, so that the resulting virial parameter will decrease by a factor of
2.25. The levels of the missing flux are unknown due to the lack of single
dish data. However, the estimated masses and radius of cores may not
significantly affected by the missing flux, because the observations are
carried out with the most compact 12 m array configuration (C40-1) and the MRS
is $\sim$25′′ that is much larger than the cores sizes. Therefore, the
estimated virial parameters are not seriously affected by the missing flux. In
addition, the ambient pressure from the parental clumps and filaments could
also provide additional confinement to make these NH2D cores bound (e.g., Chen
et al., 2019; Li et al., 2020a). Therefore, we cannot completely rule out the
possibility that the other 8 starless cores could also be pre-stellar core
candidates.

The estimated masses of NH2D cores range from 0.13 to 0.87 $M_{\odot}$,
suggesting that they are in the low-mass core regime. The majority (15) of
NH2D cores are associated with filamentary structures revealed by the H13CO+
line emission (Li et al. 2021, in preparation). Therefore, the NH2D cores have
the potential to accrete a significant amount of material from the parental
clump via filamentary accretion. The masses of continuum cores are found to
range from low- to high-mass regime (0.17 – 14.03 $M_{\odot}$; Li et al.,
2020a). These results suggest that cores with a large range of masses are
simultaneously forming in this cluster-forming cloud. The massive cores are
found toward the center of cloud (see Li et al., 2020a), while the low-mass
cores are wide spread throughout the region. The evolutionary stages vary
significantly over the sample of NH2D cores, continuum cores, and YSOs (Class
I and II). This indicates that there is a temporally continuous formation of
low- to high-mass cores.

In summary, NGC 6334S is forming a stellar protocluster, its embedded cores
show a wide diversity in masses (0.13 – 14.03 $M_{\odot}$) and evolutionary
stages (pre-stellar – Class II phases), and the embedded cores are expected to
grow in mass by gas accretion from parental clump via filaments. This seems to
be consistent with the competitive accretion massive star formation scenario
(Bonnell & Bate, 2006), in which the dense cores located near the center of
the gravitational potential continue accreting material to form massive stars.

##  5 Summary

We present ALMA and JVLA high spatial resolution observations toward the
massive IRDC NGC 6334S. We have identified 17 low-mass starless core
candidates that show bright NH2D emission, but no YSO counterparts and no
signs of protostellar activity, although we cannot completely rule out the
possibility of very weak thus undetected outflows or very faint embedded low-
mass protostars. These candidates show almost-thermal velocity dispersions
(down to $\sigma_{\rm tot}\sim$ 0.06 km s-1), high NH2D abundances (up to
$\sim$10-7), high NH3 deuterium fractionations (up to $>$ 0.39), and are
completely dark in the infrared wavelengths from 3.6 up to 70 $\mu{}$m. In
addition, the $\sigma_{\rm obs}$ appear to decrease toward the center of cores
in 10 of them. These results suggest that turbulence has significantly
dissipated and the NH2D abundance has significantly enhanced toward these
cores. The gas kinematics resemble the well known pre-stellar cores, e.g.,
Ophiuchus/H-MM1 and L 1544, but with one to two orders of magnitude greater
NH2D abundance (Crapsi et al., 2007; Harju et al., 2017). Nine out of the 17
NH2D cores are gravitationally bound, therefore are identified as pre-stellar
core candidates. Therefore, the NH2D line could be a powerful tracer to reveal
the starless and pre-stellar cores that do not show significant dust continuum
emission. This is the first detection of a cluster of low-mass starless and
pre-stellar core candidates in a massive star cluster-forming cloud. Low- to
high-mass cores are simultaneously forming in this cluster-forming cloud.

We thank the anonymous referee for constructive comments that helped improve
this paper. This work was partially supported by National Natural Science
Foundation of China (NSFC)(Grant Nos. 11988101, 11911530226). C.W.L. is
supported by the Basic Science Research Program through the National Research
Foundation of Korea (NRF) funded by the Ministry of Education, Science and
Technology (NRF-2019R1A2C1010851). H.B. acknowledges support from the European
Research Council under the European Community’s Horizon 2020 framework program
(2014-2020) via the ERC Consolidator Grant ‘From Cloud to Star Formation
(CSF)’ (project number 648505). H.B. also acknowledges funding from the
Deutsche Forschungsgemeinschaft (DFG) via the Collaborative Research Center
(SFB 881) ‘The Milky Way System’ (subproject B1). I.J.-S. has received partial
support from the Spanish State Research Agency (AEI; project number
PID2019-105552RB-C41). K.Q. is partially supported by National Key R&D Program
of China No. 2017YFA0402600, and acknowledges the National Natural Science
Foundation of China (NSFC) grant U1731237. A.P. acknowledges financial support
from the UNAM-PAPIIT IN111421 grant, the Sistema Nacional de Investigadores of
CONACyT, and from the CONACyT project number 86372 of the ‘Ciencia de Frontera
2019’ program, entitled ‘Citlalcóatl: A multiscale study at the new frontier
of the formation and early evolution of stars and planetary systems’, México.
J.M.G. acknowledges the support of the grant AYA2017-84390-C2-R (AEI/FEDER,
EU). K.W. acknowledges support by the National Key Research and Development
Program of China (2017YFA0402702, 2019YFA0405100), the National Science
Foundation of China (11973013, 11721303), NSFC grant 11629302, and the
starting grant at the Kavli Institute for Astronomy and Astrophysics, Peking
University (7101502287). H.B.L. is supported by the Ministry of Science and
Technology (MoST) of Taiwan (Grant Nos. 108-2112-M-001-002-MY3). This paper
makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.00951.S. ALMA is a
partnership of ESO (representing its member states), NSF (USA) and NINS
(Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ.

Table 1: Physical parameters of the NH2D cores.

Core ID | R.A. | Decl. | $\rm Maj\times Min$ | P.A. | $R$ | $\sigma_{\rm obs}$ | $v_{\rm LSR}$ | $\mathcal{M}$ | $W_{\rm NH_{2}D}^{\rm peak}$ | $W_{\rm NH_{2}D}$ | $S_{\rm\nu}^{\rm peak}$ | $S_{\rm\nu}$ | $M_{\rm gas}$ | $M_{\rm vir}$ | $\alpha_{\rm vir}$ | $N_{\rm NH_{2}D}$ | $N_{\rm H_{2}}$ | $X(\rm NH_{2}D)$ | $D_{\rm NH_{3}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| (J2000) | (J2000) | (′′$\times$′′) | (deg) | (pc) | (km/s) | (km/s) |  | (Jy/beam km/s) | (Jy/beam km/s) | (mJy/beam) | (mJy) | ($M_{\odot}$) | ($M_{\odot}$) |  | (cm-2) | (cm-2) |  |
M1 | 17:19:09.38 | -36:03:37.96 | 7.1$\times$2.8 | 107.6 | 0.028 | 0.11 | -3.59 |  $<$0.68 | 0.69 | 140.95 | 0.174 | 0.294 | 0.43 |  $<$0.68 |  $<$1.57 | 2.2E+15 | 1.8E+22 | 1.1E-07 | …
M2 | 17:19:04.28 | -36:06:57.17 | 9.3$\times$4.3 | 126.6 | 0.040 | 0.30 | -1.36 | 2.41 | 0.40 | 136.13 | 0.186 | 0.412 | 0.61 | 2.54 | 4.18 | 6.1E+14 | 1.6E+22 | 3.2E-08 |  $>$0.39
M3 | 17:19:05.61 | -36:11:01.29 | 7.5$\times$3.9 | 75.4 | 0.034 | 0.20 | -2.95 | 1.38 | 0.41 | 104.66 | 0.292 | 0.462 | 0.68 | 1.25 | 1.84 | 7.0E+14 | 2.4E+22 | 2.8E-08 | …
M4 | 17:19:11.12 | -36:05:41.37 | 9.4$\times$3.2 | 79.5 | 0.034 | 0.17 | -5.07 | 1.10 | 0.31 | 82.90 | 0.214 | 0.589 | 0.87 | 1.00 | 1.15 | 4.8E+14 | 2.8E+22 | 1.8E-08 | 0.25
M5 | 17:19:07.39 | -36:10:24.52 | 7.5$\times$2.1 | 72.1 | 0.025 | 0.15 | -2.47 | 0.84 | 0.44 | 82.36 | 0.180 | 0.258 | 0.38 | 0.63 | 1.65 | 9.3E+14 | 1.8E+22 | 4.5E-08 | …
M6 | 17:19:12.04 | -36:06:55.76 | 7.3$\times$3.9 | 57.3 | 0.034 | 0.18 | -2.80 | 1.13 | 0.28 | 69.97 | 0.230 | 0.458 | 0.68 | 1.08 | 1.60 | 4.0E+14 | 2.4E+22 | 1.7E-08 |  $>$0.25
M7 | 17:19:09.91 | -36:09:05.62 | 6.6$\times$5.2 | 21.1 | 0.037 | 0.23 | -4.17 | 1.69 | 0.23 | 65.62 | 0.206 | 0.447 | 0.66 | 1.83 | 2.77 | 1.8E+14 | 2.2E+22 | 8.0E-09 | …
M8 | 17:19:10.30 | -36:09:09.49 | 5.8$\times$4.4 | 122.3 | 0.032 | 0.19 | -4.22 | 1.32 | 0.24 | 56.76 | 0.237 | 0.450 | 0.66 | 1.27 | 1.91 | 3.0E+14 | 2.5E+22 | 1.2E-08 | …
M9 | 17:19:11.40 | -36:08:44.12 | 4.9$\times$3.9 | 78.0 | 0.028 | 0.12 | -1.68 | 0.34 | 0.24 | 45.38 | 0.158 | 0.201 | 0.30 | 0.73 | 2.47 | 2.7E+14 | 1.4E+22 | 1.4E-08 | 0.22
M10 | 17:19:05.12 | -36:06:45.00 | 4.2$\times$3.7 | 9.9 | 0.025 | 0.15 | -3.47 | 0.78 | 0.26 | 45.19 | 0.310 | 0.366 | 0.54 | 0.79 | 1.46 | 6.9E+14 | 2.8E+22 | 2.5E-08 |  $>$0.11
M11 | 17:18:58.66 | -36:07:04.76 | 4.1$\times$2.2 | 48.5 | 0.019 | 0.11 | -2.01 |  $<$0.73 | 0.33 | 44.30 | … | … |  $<$0.13 |  $<$0.51 | 3.82 | 7.5E+14 |  $<$1.5E+22 |  $>$7.1E-08 |  $>$0.32
M12 | 17:19:06.86 | -36:06:22.65 | 6.8$\times$3.6 | 124.3 | 0.031 | 0.13 | -3.51 | 0.39 | 0.18 | 41.11 | 0.244 | 0.538 | 0.79 | 0.74 | 0.93 | 2.5E+14 | 3.0E+22 | 8.9E-09 |  $>$0.16
M13 | 17:19:07.39 | -36:05:55.13 | 6.9$\times$3.7 | 116.2 | 0.032 | 0.10 | -5.36 |  $<$0.57 | 0.15 | 36.21 | … | … |  $<$0.13 |  $<$0.80 | 6.03 | 1.3E+14 |  $<$1.5E+22 |  $>$7.8E-09 |  $>$0.14
M14 | 17:19:04.55 | -36:05:31.47 | 3.7$\times$2.2 | 29.7 | 0.018 | 0.16 | -4.20 | 0.93 | 0.26 | 33.59 | … | … |  $<$0.13 | 0.54 |  $>$4.04 | 2.8E+14 |  $<$1.5E+22 |  $>$1.2E-08 | …
M15 | 17:19:08.88 | -36:08:01.55 | 5.3$\times$2.7 | 78.2 | 0.024 | 0.11 | -4.26 |  $<$0.70 | 0.19 | 29.68 | … | … |  $<$0.13 |  $<$0.62 | 4.64 | 2.0E+14 |  $<$1.5E+22 |  $>$2.0E-08 |  $>$0.17
M16 | 17:19:07.21 | -36:08:16.60 | 6.4$\times$2.3 | 125.7 | 0.024 | 0.15 | -4.18 | 0.81 | 0.16 | 29.59 | 0.160 | 0.168 | 0.25 | 0.60 | 2.44 | 1.6E+14 | 1.2E+22 | 1.1E-08 |  $>$0.14
M17 | 17:19:07.78 | -36:10:07.54 | 4.8$\times$1.8 | 57.5 | 0.019 | 0.13 | -2.52 | 0.49 | 0.13 | 18.00 | 0.158 | 0.170 | 0.25 | 0.43 | 1.70 | 1.6E+14 | 1.6E+22 | 7.2E-09 | …
Mean |  |  |  |  | 0.028 | 0.16 |  | 0.96 | 0.29 | 62.50 | 0.211 | 0.370 | 0.45 | 0.94 | 2.60 | 5.1E+14 | 2.0E+22 | 2.6E-08 | 0.22
Median |  |  |  |  | 0.028 | 0.15 |  | 0.81 | 0.26 | 45.38 | 0.206 | 0.412 | 0.43 | 0.74 | 1.91 | 3.0E+14 | 1.8E+22 | 1.7E-08 | 0.20
Minimum |  |  |  |  | 0.018 | 0.10 |  | 0.34 | 0.13 | 18.00 | 0.158 | 0.168 | 0.13 | 0.43 | 0.93 | 1.3E+14 | 1.2E+22 | 7.2E-09 | 0.11
Maximum |  |  |  |  | 0.040 | 0.30 |  | 2.41 | 0.69 | 140.95 | 0.310 | 0.589 | 0.87 | 2.54 | 6.03 | 2.2E+15 | 3.0E+22 | 1.1E-07 | 0.39

††footnotetext: Notes. For M1, M11, M13, and M15, $\sigma_{\rm nt,NH_{2}D}$ =
$\sigma_{\rm obs}$ was used since the $\sigma_{\rm
obs,int}\leqslant\sigma_{\rm nt,NH_{2}D}(\rm 10K)$. (4)-(5) beam-deconvolved
size. (6) beam-deconvolved effective radius. (7)-(8) $\sigma_{\rm obs}$ and
$v_{\rm LSR}$ are derived by the core-averaged spectrum. (9) Mach number
$\mathcal{M}$ = $\sqrt{3}\sigma_{\rm nt,NH_{2}D}/c_{\rm s}$. (10) the peak
value of NH2D velocity-integrated image. (11) the integrated flux density of
the NH2D line emission. (12)-(13) the peak and integrated intensity of
continuum emission. (14) gas mass estimated by continuum emission. (15) virial
mass. (16) virial parameter. (17) mean NH2D column density. (18) mean H2
column density. (19) mean NH2D abundance fraction. (20) the NH3 deuterium
fractionation $D_{\rm NH_{3}}$.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Auddy et al. (2019) Auddy, S., Myers, P. C., Basu, S., et al. 2019, ApJ, 872, 207, doi: 10.3847/1538-4357/ab0086
  * Bergin & Tafalla (2007) Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339, doi: 10.1146/annurev.astro.45.071206.100404
  * Bonnell & Bate (2006) Bonnell, I. A., & Bate, M. R. 2006, MNRAS, 370, 488, doi: 10.1111/j.1365-2966.2006.10495.x
  * Caselli & Ceccarelli (2012) Caselli, P., & Ceccarelli, C. 2012, A&A Rev., 20, 56, doi: 10.1007/s00159-012-0056-x
  * Chen et al. (2019) Chen, H. H.-H., Pineda, J. E., Goodman, A. A., et al. 2019, ApJ, 877, 93, doi: 10.3847/1538-4357/ab1a40
  * Chibueze et al. (2014) Chibueze, J. O., Omodaka, T., Handa, T., et al. 2014, ApJ, 784, 114, doi: 10.1088/0004-637X/784/2/114
  * Chira et al. (2014) Chira, R.-A., Smith, R. J., Klessen, R. S., Stutz, A. M., & Shetty, R. 2014, MNRAS, 444, 874, doi: 10.1093/mnras/stu1497
  * Crapsi et al. (2005) Crapsi, A., Caselli, P., Walmsley, C. M., et al. 2005, ApJ, 619, 379, doi: 10.1086/426472
  * Crapsi et al. (2007) Crapsi, A., Caselli, P., Walmsley, M. C., & Tafalla, M. 2007, A&A, 470, 221, doi: 10.1051/0004-6361:20077613
  * Daniel et al. (2016) Daniel, F., Coudert, L. H., Punanova, A., et al. 2016, A&A, 586, L4, doi: 10.1051/0004-6361/201527768
  * Dunham et al. (2008) Dunham, M. M., Crapsi, A., Evans, Neal J., I., et al. 2008, ApJS, 179, 249, doi: 10.1086/591085
  * Dunham et al. (2016) Dunham, M. M., Offner, S. S. R., Pineda, J. E., et al. 2016, ApJ, 823, 160, doi: 10.3847/0004-637X/823/2/160
  * Enoch et al. (2008) Enoch, M. L., Evans, Neal J., I., Sargent, A. I., et al. 2008, ApJ, 684, 1240, doi: 10.1086/589963
  * Evans et al. (2009) Evans, Neal J., I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321, doi: 10.1088/0067-0049/181/2/321
  * Fischer et al. (2017) Fischer, W. J., Megeath, S. T., Furlan, E., et al. 2017, ApJ, 840, 69, doi: 10.3847/1538-4357/aa6d69
  * Ginsburg & Mirocha (2011) Ginsburg, A., & Mirocha, J. 2011, PySpecKit: Python Spectroscopic Toolkit, Astrophysics Source Code Library.  http://ascl.net/1109.001
  * Gómez et al. (2021) Gómez, G. C., Vázquez-Semadeni, E., & Palau, A. 2021, MNRAS, doi: 10.1093/mnras/stab394
  * Harju et al. (2017) Harju, J., Daniel, F., Sipilä, O., et al. 2017, A&A, 600, A61, doi: 10.1051/0004-6361/201628463
  * Hildebrand (1983) Hildebrand, R. H. 1983, QJRAS, 24, 267
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90
  * Kauffmann et al. (2008) Kauffmann, J., Bertoldi, F., Bourke, T. L., Evans, N. J., I., & Lee, C. W. 2008, A&A, 487, 993, doi: 10.1051/0004-6361:200809481
  * Kauffmann et al. (2013) Kauffmann, J., Pillai, T., & Goldsmith, P. F. 2013, ApJ, 779, 185, doi: 10.1088/0004-637X/779/2/185
  * Kirk et al. (2017) Kirk, H., Dunham, M. M., Di Francesco, J., et al. 2017, ApJ, 838, 114, doi: 10.3847/1538-4357/aa63f8
  * Lee et al. (1999) Lee, C. W., Myers, P. C., & Tafalla, M. 1999, ApJ, 526, 788, doi: 10.1086/308027
  * Lee et al. (2001) —. 2001, ApJS, 136, 703, doi: 10.1086/322534
  * Li et al. (2013) Li, D., Kauffmann, J., Zhang, Q., & Chen, W. 2013, ApJ, 768, L5, doi: 10.1088/2041-8205/768/1/L5
  * Li et al. (2019) Li, S., Zhang, Q., Pillai, T., et al. 2019, ApJ, 886, 130, doi: 10.3847/1538-4357/ab464e
  * Li et al. (2020a) Li, S., Zhang, Q., Liu, H. B., et al. 2020a, ApJ, 896, 110, doi: 10.3847/1538-4357/ab84f1
  * Li et al. (2020b) Li, S., Sanhueza, P., Zhang, Q., et al. 2020b, ApJ, 903, 119, doi: 10.3847/1538-4357/abb81f
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Nony et al. (2018) Nony, T., Louvet, F., Motte, F., et al. 2018, A&A, 618, L5, doi: 10.1051/0004-6361/201833863
  * Offner & McKee (2011) Offner, S. S. R., & McKee, C. F. 2011, ApJ, 736, 53, doi: 10.1088/0004-637X/736/1/53
  * Pineda et al. (2015) Pineda, J. E., Offner, S. S. R., Parker, R. J., et al. 2015, Nature, 518, 213, doi: 10.1038/nature14166
  * Punanova et al. (2018) Punanova, A., Caselli, P., Feng, S., et al. 2018, ApJ, 855, 112, doi: 10.3847/1538-4357/aaad09
  * Robitaille & Bressert (2012) Robitaille, T., & Bressert, E. 2012, APLpy: Astronomical Plotting Library in Python, Astrophysics Source Code Library.  http://ascl.net/1208.017
  * Rosolowsky et al. (2008) Rosolowsky, E. W., Pineda, J. E., Kauffmann, J., & Goodman, A. A. 2008, ApJ, 679, 1338, doi: 10.1086/587685
  * Sahu et al. (2021) Sahu, D., Liu, S.-Y., Liu, T., et al. 2021, ApJ, 907, L15, doi: 10.3847/2041-8213/abd3aa
  * Sanhueza et al. (2013) Sanhueza, P., Jackson, J. M., Foster, J. B., et al. 2013, ApJ, 773, 123, doi: 10.1088/0004-637X/773/2/123
  * Sanhueza et al. (2019) Sanhueza, P., Contreras, Y., Wu, B., et al. 2019, ApJ, 886, 102, doi: 10.3847/1538-4357/ab45e9
  * Schnee et al. (2010) Schnee, S., Enoch, M., Johnstone, D., et al. 2010, ApJ, 718, 306, doi: 10.1088/0004-637X/718/1/306
  * Sipilä et al. (2015) Sipilä, O., Harju, J., Caselli, P., & Schlemmer, S. 2015, A&A, 581, A122, doi: 10.1051/0004-6361/201526468
  * Tigé et al. (2017) Tigé, J., Motte, F., Russeil, D., et al. 2017, A&A, 602, A77, doi: 10.1051/0004-6361/201628989
  * Tokuda et al. (2020) Tokuda, K., Fujishiro, K., Tachihara, K., et al. 2020, ApJ, 899, 10, doi: 10.3847/1538-4357/ab9ca7
  * Ward-Thompson et al. (1994) Ward-Thompson, D., Scott, P. F., Hills, R. E., & Andre, P. 1994, MNRAS, 268, 276, doi: 10.1093/mnras/268.1.276
  * Willis et al. (2013) Willis, S., Marengo, M., Allen, L., et al. 2013, ApJ, 778, 96, doi: 10.1088/0004-637X/778/2/96

##  Appendix A NH2D core identification

The astrodendro222 http://dendrograms.org/ algorithm was used to preselect the
compact NH2D emission structures (i.e., the leaves in the terminology of
astrodendro; Rosolowsky et al., 2008) from the NH2D velocity-integrated
intensity image (integrated between -11.5 and 5.6 km s-1). The following
parameters are used in the computation of astrodendro: the minimum pixel value
is 3$\sigma$ noise level (1$\sigma\sim$ 0.02 Jy beam-1 km s-1 for $W_{\rm
NH_{2}D}$); the peak value is $\geqslant$ 6$\sigma$; the minimum difference in
the peak intensity between neighboring compact structures is 1$\sigma$; and
the minimum number of pixels required for a structure to be considered an
independent entity is 40 that is approximately the synthesized beam area. From
the identified leaves, we further search for the compact structures that are
no associated with continuum cores, YSOs, and showing a single peak and no
extreme elongations (major axis smaller than 3 times minor axis). The physical
parameters (size, flux density, peak intensity) of these selected compact
structures are extracted from the NH2D velocity-integrated intensity image
using the CASA-imfit task (see Table 1). We stress that this identification is
likely to be incomplete. Potential NH2D compact structures could have been
missed if they cannot be distinguished from the clumpy or filamentary
structures. Four NH2D cores (M11, M13, M14, M15) are not associated with
continuum emission at all ($<3\sigma$, 1$\sigma\sim$ 30 $\mu$Jy beam-1 for
continuum image), and the rest of the NH2D cores are partially associated with
weak continuum emission ($3\sigma-7\sigma$), except for M10 that has about 10%
region at the edge of the core associated with $8\sigma-11\sigma$ continuum
emission (Figure 2).

##  Appendix B Gas mass and Virial mass

With the observed 3 mm continuum fluxes, we estimated the gas mass ($M_{\rm
gas}$) of identified cores assuming optically thin modified black body
emission in the Rayleigh-Jeans limit, following

$M_{\rm gas}=\eta\frac{S_{\nu}d^{2}}{B_{\nu}(\rm T)\,\kappa_{\nu}},$ (B1)

where $\eta$ = 100 is the assumed gas-to-dust mass ratio, $d$ is the source
distance, $S_{\nu}$ is the continuum flux at frequency 98.5 GHz, $B_{\nu}(\rm
T)$ is the Planck function at temperature $T$, and $\kappa_{\nu}$ is the dust
opacity at frequency $\nu$. We adopt $\kappa_{\rm 98.5\;GHz}$ = 0.235 cm2g-1,
assuming $\kappa_{\nu}$ = 10(${\nu}$/1.2 THz)β cm2g-1 and $\beta$ = 1.5 for
all of dense cores (Hildebrand, 1983).

The column density ($N_{\rm H_{2}}$) is computed with

$N_{\rm H_{2}}=\eta\frac{S_{\nu}^{\rm beam}}{B_{\nu}({\rm
T})\,\kappa_{\nu}\,\Omega\,\mu\,m_{\rm H}},$ (B2)

where $S_{\nu}^{\rm beam}$ is the continuum flux density, $\Omega$ is the beam
solid angle, $\mu$ = 2.8 is the mean molecular weight of the interstellar
medium (Kauffmann et al., 2008), and $m_{\rm H}$ is the proton mass.

The virial mass is estimated using:

$M_{\rm vir}=\frac{5}{a\beta}\frac{\sigma_{\rm tot}^{2}R}{G},$ (B3)

where $\sigma_{\rm tot}\;=\;\sqrt{\sigma_{\rm nt,NH_{2}D}^{2}+\sigma_{\rm
th,\langle m\rangle}^{2}}\;=\;\sqrt{(\sigma_{\rm obs,int}^{2}-\sigma_{\rm
th,NH_{2}D}^{2})+\sigma_{\rm th,\langle m\rangle}^{2}}$ is the total velocity
dispersion, $R$ is the effective radius, $G$ is the gravitational constant,
parameter $a$ equals to $(1-b/3)/(1-2b/5)$ for a power-law density profile
$\rho\propto r^{-b}$, and $\beta=(\arcsin e)/e$ is the geometry factor (see
also Li et al., 2013, 2020a). Here $\langle m\rangle$ is 2.37$m_{\rm H}$
(Kauffmann et al., 2008), $\sigma_{\rm obs,int}$ is derived by the core-
averaged spectrum, and $b$ is assumed to be 2. For M1, M11, M13, and M15,
$\sigma_{\rm tot}$ = $\sigma_{\rm th,\langle m\rangle}(\rm 10\,K)$ was used
since the $\sigma_{\rm obs,int}\leqslant\sigma_{\rm th,NH_{2}D}(\rm 10\,K)$.

##  Appendix C Internal Luminosity

The internal luminosity of dense object can be estimated by using the
empirical relationship between internal luminosity ($L_{\rm int}$) and
observed 70 $\mu{}$m flux (Dunham et al., 2008):

$L_{\rm int}(L_{\odot})=3.3\times 10^{8}\left(F_{\rm
70}\,\left(\frac{d}{140\rm pc}\right)^{2}\right)^{0.94},$ (C1)

where $F_{\rm 70}$ is the flux at 70 $\mu{}$m in cgs units (erg s-1 cm-2 Hz-1)
and $d$ is the source distance in pc. A 3$\sigma$ sensitivity of $\sim$0.3 Jy
from Herschel/PACS 70 $\mu{}$m observation is used to derived an upper limit
($\sim$1.26 $L_{\odot}$) of internal luminosity for the identified NH2D cores
(Tigé et al., 2017), with the exception of M3, M5, and M17, which are
illuminated by the western YSOs rather than internal heating. The measured 70
$\mu{}$m fluxes of M3, M5, and M7 are about 5 – 6$\sigma$.