11institutetext:  Institute of Astronomy, Graduate School of Science, The
University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan
11email<EMAIL_ADDRESS>22institutetext: School of Astronomy and
Space Science, Nanjing University, Nanjing 210023, China 33institutetext:
School of Physics and Technology, Nanjing Normal University, Nanjing 210023,
China 44institutetext: National Astronomical Observatory of Japan, 2-21-1
Osawa, Mitaka, Tokyo 181-8588, Japan
44email<EMAIL_ADDRESS>55institutetext: Department of Astronomy,
School of Science, The Graduate University for Advanced Studies (SOKENDAI),
2-21-1 Osawa, Mitaka Tokyo, 181-1855 Japan 66institutetext: ALMA Project,
National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo
181-8588, Japan 77institutetext: Department of Physics, Faculty of Science and
Technology, Keio University, 3-14-1 Hiyoshi, Yokohama, Kanagawa 223–8522 Japan
88institutetext: European Southern Observatory, Alonso de Córdova, 3107,
Vitacura, Santiago 763-0355, Chile 99institutetext: Joint ALMA Observatory,
Alonso de Córdova, 3107, Vitacura, Santiago 763-0355, Chile 1010institutetext:
National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA
22903-2475, USA 1111institutetext: Institute of Astronomy and Astrophysics,
Academia Sinica, 11F of AS/NTU Astronomy-Mathematics Building, No.1, Sec. 4,
Roosevelt Rd, Taipei 10617, Taiwan 1212institutetext: Department of Space,
Earth and Environment, Chalmers University of Technology, Onsala Space
Observatory, SE-439 92 Onsala, Sweden 1313institutetext: Leiden Observatory,
Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands
1414institutetext: Centro de Astrobiología (CSIC-INTA), Ctra. de Ajalvir Km.
4, Torrejón de Ardoz, 28850 Madrid, Spain 1515institutetext: New Mexico
Institute of Mining and Technology, 801 Leroy Place, Socorro, NM 87801, USA
1616institutetext: National Radio Astronomy Observatory, PO Box O, 1003
Lopezville Road, Socorro, NM 87801, USA 1717institutetext: Max-Planck-Institut
für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany 1818institutetext:
Astronomy Department, Faculty of Science, King Abdulaziz University, P. O. Box
80203, Jeddah 21589, Saudi Arabia 1919institutetext: Departamento de
Astronomia, Instituto de Astronomia, Geofísica e Ciências Atmosféricas da USP,
Cidade Universitária, 05508-900 São Paulo, SP, Brazil

# Physical Properties of the Southwest Outflow Streamer in the Starburst
Galaxy NGC 253 with ALCHEMI

Min Bao 112233 Nanase Harada 4455 Kotaro Kohno 11 Yuki Yoshimura 11 Fumi Egusa
11 Yuri Nishimura 1166 Kunihiko Tanaka 77 Kouichiro Nakanishi 4455 Sergio
Martín 8899 Jeffrey G. Mangum 1010 Kazushi Sakamoto 1111 Sébastien Muller 1212
Mathilde Bouvier 1313 Laura Colzi 1414 Kimberly L. Emig 1010 David S. Meier
15151616 Christian Henkel 17171818 Pedro Humire 19191717 Ko-Yun Huang 1313
Víctor M. Rivilla 1414 Paul van der Werf 1313 Serena Viti 1313

###### Abstract

Aims. The physical properties of galactic molecular outflows are important as
they could constrain outflow formation mechanisms. In this work, we study the
properties of the southwest (SW) outflow streamer including gas kinematics,
optical depth, dense gas fraction, and shock strength through molecular
emission in the central molecular zone of the starburst galaxy NGC 253.

Methods. We image the molecular emission in NGC 253 at a spatial resolution of
1.6′′($\sim$27 pc at D$\sim$3.5 Mpc) based on data from the ALMA Comprehensive
High-resolution Extragalactic Molecular Inventory (ALCHEMI) large program. We
trace the velocity and the velocity dispersion of molecular gas with CO(1-0)
line and study the molecular spectra in the SW streamer region, the brightest
CO streamer in NGC 253. We constrain the optical depth of CO emission with
CO/13CO(1-0) ratio, the dense gas fraction with HCN/CO(1-0), H13CN/13CO(1-0)
and N2H+/13CO(1-0) ratios, as well as the shock strength with
SiO(2-1)/13CO(1-0) and CH3OH(2k-1k)/13CO(1-0) ratios.

Results. The CO/13CO(1-0) integrated intensity ratio is $\sim$21 in the SW
streamer region, which approximates the C/13C isotopic abundance ratio. The
higher integrated intensity ratio compared to the disk can be attributed to
the optically thinner environment for CO(1-0) emission inside the SW streamer.
The HCN/CO(1-0) and SiO(2-1)/13CO(1-0) integrated intensity ratios both
approach $\sim$0.2 in three giant molecular clouds (GMCs) at the base of the
outflow streamers, which implies the higher dense gas fraction and enhanced
strength of fast shocks in those GMCs than in the disk, while the HCN/CO(1-0)
integrated intensity ratio is moderate in the SW streamer region. The contours
of those two integrated intensity ratios are extended towards the directions
of outflow streamers, which connects the enhanced dense gas fraction and shock
strength with molecular outflow. Moreover, the molecular gas with enhanced
dense gas fraction and shock strength located at the base of the SW streamer
shares the same velocity with the outflow.

Conclusions. The enhanced dense gas fraction and shock strength at the base of
the outflow streamers suggest that the star formation inside the GMCs can
trigger the shocks and further drive the molecular outflow. The increased
CO/13CO(1-0) integrated intensity ratio coupled with the moderate HCN/CO(1-0)
integrated intensity ratio in the SW streamer region are consistent with the
picture that the gas velocity gradient inside the streamer may decrease the
optical depth of CO(1-0) emission, as well as the dense gas fraction in the
extended streamer region.

###### Key Words.:

Galaxies: individual: NGC 253 – Galaxies: evolution – Galaxies: starburst –
Galaxies: kinematics and dynamics

##  1 Introduction

Outflows on the galactic scale are ubiquitous in the local and distant
Universe. Together with gas accretion, star formation, and black hole growth,
outflows can govern the cycle of material between a galaxy and the
circumgalactic medium (Veilleux et al., 2005). Cool outflows, which include
atomic and molecular gas, as well as dust, dominate the mass and energy of the
outflowing material. The study of extragalactic cool outflows only dates back
to $\sim$20 years, which is a relatively new field in the research of galactic
outflows (Veilleux et al., 2020). The physical properties of cool outflows are
gradually uncovered thanks to the development of high-capability radio
telescopes, including ALMA (ALMA Partnership et al., 2015), NOEMA (Lefèvre et
al., 2020), and VLA (Greisen, 2003).

As one of the key phases of cool gas, molecular gas is the raw material for
star formation, an important process in galaxy evolution. More than 300
molecular species have been detected in the circumstellar envelopes of evolved
stars and the interstellar medium (ISM) of the Galaxy
111https://cdms.astro.uni-koeln.de/classic/molecules. One-third of those
species have also been detected in external galaxies (Takano et al., 2019;
Martín et al., 2021). Among them, carbon monoxide (CO) is the second most
abundant molecule after H2 and serves as the main observational tracer of
molecular gas. Moreover, other molecules with different chemical formation
pathways and excitation requirements provide supplementary information to
constrain the evolution of external galaxies (Aalto, 2015).

The formation scenarios for molecular outflows can be roughly divided into
three categories. The outflowing molecular gas can be directly driven by
radiation (Murray et al., 2011) and/or pressure gradients (Socrates et al.,
2008; Uhlig et al., 2012). Alternatively, hot winds can entrain molecular
clouds (Banda-Barragán et al., 2019; Fielding & Bryan, 2022), while the
lifetime of the molecular clouds depends on the balance between radiation,
conduction, and turbulence (Orlando et al., 2005). The third scenario is that
molecular gas forms in situ from hot winds through cooling and/or thermal
instability. The first step in the formation of molecules consists of hydrogen
atoms (H) combining to form hydrogen molecules (H2). Given that H atoms can
combine on the surface of dust grains, those grains are catalysts for the
formation of H2 (Hollenbach & Salpeter, 1971; Pantaleone et al., 2021).
However, along with temperature rises, dust grains undergo destructive
sputtering via gas-grain collisions in the hot winds, which reduces the
efficiency of H2 formation (Le Bourlot et al., 2012). Hence, a cooling
timescale shorter than a dynamical timescale is a premise for the in situ
formation (Efstathiou, 2000; Silich et al., 2003).

There is no agreement on the influence of outflows on the star formation
activity inside host galaxies. On the one hand, the hot ionized winds are one
of the leading processes that can quench star formation (Tacchella et al.,
2015, 2016; Spilker et al., 2019), namely, negative feedback. The quenching is
owing to the energy that is injected into molecular clouds, or to the
entrainment/depletion of molecular clouds. On the other hand, the shocks,
which form through the encounter between winds and molecular clouds, can
compress the molecular clouds and trigger star formation (Klein et al., 1994),
namely, positive feedback.

The engines for outflows can be central starbursts and/or active galactic
nuclei (AGN). Different physical processes including thermal energy,
radiation, cosmic rays, and/or radio jets can work together in either driving
mechanism. AGN-driven molecular outflows have been well studied in a few
nearby cases such as the Seyfert 2 galaxy NGC 1068 (Kormendy & Ho, 2013) and
the Quasar Mrk 231 (Veilleux et al., 2009). However, to study how molecular
outflows relate to star formation, we need to target sources without
contamination from AGN.

NGC 253 (D$\sim$3.5 Mpc, Rekola et al. 2005) is one of the nearest starburst
galaxies. It possesses a strong bar in the center (Das et al., 2001; Paglione
et al., 2004), but shows no sign of AGN activity (Brunthaler et al., 2009). It
is an edge-on (i$\sim$78.5∘) spiral galaxy and hosts a bipolar outflow
(Turner, 1985; Westmoquette et al., 2011). The outflow was proven to be driven
by the starburst that has been occurring in the central $\sim$500 pc region
during the last $\sim$20-30 Myr (Rieke et al., 1980; Engelbracht et al.,
1998), where the star formation rate is $\sim$3 M⊙ yr-1 (Ott et al., 2005;
Bendo et al., 2015). The outflow has been detected in multi-phases, including
molecular gas phase (Turner, 1985; Sturm et al., 2011; Heesen et al., 2011;
Bolatto et al., 2013; Walter et al., 2017; Zschaechner et al., 2018; Krieger
et al., 2019), ionized gas phase (Heckman et al., 2000; Westmoquette et al.,
2011; Cohen et al., 2020), X-ray emitting gas phase (Pietsch et al., 2000;
Strickland et al., 2000, 2002; Lopez et al., 2023), and dust phase (Levy et
al., 2022). The ionized outflow shows a wide opening angle of $\sim$60∘, a
deprojected velocity of a few 100 km s-1 (Westmoquette et al., 2011), and an
extension of $\sim$10 kpc (Strickland et al., 2002).

There have been several ALMA-based studies of the molecular outflow of NGC
253. Krieger et al. (2019) presented ALMA observation of CO(3-2) in the
central 30′′ starburst region. They obtained the non-disk component by
modeling the disk, deducted the disk from position-velocity diagrams (PVDs),
and found that $\sim$7%-16% of CO luminosity is emitted by the non-disk
component. Bolatto et al. (2013) presented ALMA observation of CO(1-0) in the
central arcminute, and found the extraplanar molecular gas closely tracking
the H$\alpha$ filaments; they estimated the mass outflow rate assuming an
optically thin conversion factor ($\alpha_{\rm CO}$). By comparing the
molecular mass outflow rate (9 $M_{\odot}$ yr-1) with the star formation rate
(2.8 $M_{\odot}$ yr-1, Ott et al., 2005), Bolatto et al. (2013) concluded that
the starburst-driven outflow limits the star-formation activity in NGC 253.
Zschaechner et al. (2018) presented ALMA observations of CO(1-0) and CO(2-1)
in the central 40′′ region. They located the molecular outflow via the
extended structure in the velocity-integrated intensity map of CO(2-1). Based
on the velocity-integrated intensity ratio between CO(2-1) and CO(1-0), they
constrained the optical depth of the CO emission in the outflow region. Walter
et al. (2017) presented a detailed study of one streamer of the molecular
outflow on the southwestern side (SW streamer). In addition to the CO
emission, they found many bright tracers of dense gas such as HCN, CN, HCO+,
and CS in the SW streamer.

This paper is based on the data from the ALMA Comprehensive High-resolution
Extragalactic Molecular Inventory (ALCHEMI) project (Martín et al., 2021). The
ALCHEMI data targeting the central molecular zone (CMZ) of NGC 253 with high
resolution, provide insight into the physical condition in the center of this
galaxy (Haasler et al., 2022). Those data have been used to study the
abundance and excitation of different molecular species (Holdship et al.,
2021; Harada et al., 2022) in the CMZ of NGC 253, and how they vary with
different heating environments (Harada et al., 2021; Holdship et al., 2022;
Behrens et al., 2022). In this paper, we uncover the physical properties of
the molecular outflow in NGC 253, focusing on the properties of the SW
streamer. Taking advantage of the numerous molecular species in the ALCHEMI
survey, we study the properties through various integrated intensity ratios.
In detail, we use the CO/13CO(1-0) ratio to probe the optical depth (Israel,
2020), HCN/CO(1-0), H13CN/13CO(1-0) and N2H+/13CO(1-0) ratios to probe the
dense gas fraction (Gao & Solomon, 2004a; Barnes et al., 2020). Moreover, we
use SiO(2-1)/13CO(1-0) and CH3OH(2k-1k)/13CO(1-0) ratios to probe the shock
strength (García-Burillo et al., 2010).

This paper is organized as follows. The data analysis is presented in Section
2. The physical properties of the molecular outflow, including gas kinematics,
optical depth, dense gas fraction, and shock strength, are studied in Section
3. In Section 4, we analyse the correlation between the dense gas fraction and
the strength of fast shocks, and discuss how the molecular outflow relates to
the star formation in NGC 253. Finally, the results are summarized in Section
5.

##  2 Data analysis

The data used in this study have been obtained as part of the ALCHEMI survey,
which is an ALMA Cycle 5 large program (2017.1.00161.L). It consists of a wide
and unbiased spectral survey covering the CMZ of NGC 253 in Bands 3, 4, 6, and
7. The phase center of observation is $\alpha=00^{\rm h}47^{\rm m}33.26^{\rm
s}$, $\delta=-25^{\circ}17^{\prime}17.7^{\prime\prime}$. In this paper, we
only make use of the Band 3 data, which covers the ground-state transitions of
common molecular lines that can be emitted even from the coldest ($\sim$10 K)
medium. We refer readers to Martín et al. (2021) for more details about the
data analysis, and only enumerate here the relevant information to this work.

The data cubes from the ALCHEMI survey are imaged to a spatial resolution of
1.6′′ ($\sim 27$ pc) and a spectral resolution of $\Delta$v$\sim$10 km s-1. We
uniformly use the continuum emission subtracted cubes, which are created by
the Python-based tool STATCONT (Sánchez-Monge et al., 2018). Those cubes are
primary beam corrected in the ALCHEMI imaging process. The Högbom deconvolver
function is used for all the cubes from the ALCHEMI survey. Meanwhile, we
produce the self-calibrated data cubes with a multi-scale deconvolver function
for the integrated intensities of CO and 13CO in the J = 1-0 transition in
Figs. 5(a) and 10(a). The use of self-calibration allows for a higher signal-
to-noise ratio and is better for images with complicated spatial structures,
which helps us in precisely obtaining the optical depth of CO emission in both
disk and outflow regions (see Section 3.3).

The data cubes from the ALCHEMI survey are subject to the missing flux problem
due to the lack of short spacing in interferometric observations. However, the
band 3 observations here probe scales up to $\sim$30′′, while spatial
filtering should be relevant to larger scales (Martín et al., 2021). In order
to quantify the missing fluxes, we collect the results of IRAM 30m single-dish
observations from Aladro et al. (2015). We unify the beam size of the ALCHEMI
data cubes to that of the IRAM 30m observations using the CASA command
imsmooth (The CASA Team et al., 2022) and compare the integrated intensities
between them. All the missing fluxes of the molecular lines in this work are
less than a few percent.

To ensure sufficient signal-to-noise ratio (S/N), we mask the region without
3$\sigma$-detection. Details of our procedures are as follows. We take the
data cubes before primary beam correction as references, which share the same
dimensions as those after primary beam correction. For each targeted molecular
line, we mask the region without 3$\sigma$-detection for each channel in the
corresponding data cube. The integrated intensity, velocity, and velocity
dispersion maps, i.e., moments 0, 1, and 2 maps, in the CMZ (inner $\sim$500
pc) of NGC 253 are generated by the CASA command immoments (The CASA Team et
al., 2022). The PVDs are generated by the Python-based tool pvextractor
(Ginsburg et al., 2016). Given the 3$\sigma$-detection limits prior to
analyses, all the features presented in the following maps and PVDs will be
robust.

##  3 Physical properties of the molecular outflow

###  3.1 Molecular emission

Fig. 1 shows the integrated intensity maps of CO, 13CO, HCN, H13CN, and N2H+
in the J = 1-0 transition and CH3OH (methanol) in the Jk = 2k-1k transition
series. We extract spectral channels covering $1000\,\rm km\,s^{-1}$ velocity
range around each transition from the continuum-subtracted cubes, where the
rest-frame frequencies are 115.271, 110.201, 88.632, 86.340, 93.174, and
$\sim$96.741 GHz, respectively. The black solid line in Fig. 1(a) marks the
galactic major axis (hereafter major axis) with a position angle of 55∘, and
is defined following Krieger et al. (2019). The strong CO(1-0) emission in
Fig. 1(a), which is marked by the blue contour at 4300 K km s-1, concentrates
in the galactic center and distributes along the major axis. Similar
structures also exist in other panels of Fig. 1. Das et al. (2001) found that
a bar structure along the major axis is necessary for modeling the ionized gas
velocity field in the central $\sim$100 pc region of NGC 253. A bar structure
is also applicable in the CO map in the central $\sim$300 pc region of NGC 253
(Paglione et al., 2004). A stellar bar, which is known as an efficient
mechanism for gas accretion, can trigger strong molecular emission along the
major axis as is seen in Fig. 1. One different point in Fig. 1(f) from the
other panels in Fig. 1 is that the strongest CH3OH(2k-1k) emission exists on
the outskirts of the gas disk, suggesting quasi-thermal emission in agreement
with Humire et al. (2022).

Figure 1: Integrated intensity maps of the transitions used in this work. (a)
Integrated intensity map of CO(1-0). The white contour is drawn at 180 K km
s-1, the grey contours are drawn at [360, 720, 1440, 2880] K km s-1, and the
blue contour is drawn at 4300 K km s-1. The horizontal lines in different
colors inside the colorbar mark the values for corresponding contours, and
remain the same in the following figures. black solid line marks the major
axis, and the black dashed line marks the slice across the SW streamer, on
which the black arrows show the positive directions for the PVDs in the
following figures. The white arrow points to the SW streamer, and remains the
same in the following panels. The black cross marks the phase center of the
observation, and remains the same in the following figures. The beam size of
1.6′′ is shown as a black empty circle at the bottom-left corner of each
panel, and remains the same in the following figures. (b) Integrated intensity
map of 13CO(1-0). The grey contours are drawn at [12, 40, 120, 230, 460] K km
s-1. (c) Integrated intensity map of HCN(1-0). The grey contours are drawn at
[12, 40, 120, 300, 670] K km s-1. (d) Integrated intensity map of H13CN(1-0).
The grey contours are drawn at [2, 8, 20, 50, 100] K km s-1. (e) Integrated
intensity map of N2H+(1-0). The grey contours are drawn at [5, 20, 50, 100,
200] K km s-1. (f) Integrated intensity map of CH3OH(2k-1k). The grey contours
are drawn at [5, 20, 50, 100, 200] K km s-1.

The CO(1-0) contours (especially the white one at 180 K km s-1) in Fig. 1(a)
are extended towards the SW streamer defined by Walter et al. (2017), which is
indicated by a white arrow. Similar extent can be found in the HCN(1-0)
contours (Fig. 1c), while the 13CO(1-0) contours are not so extended (Fig.
1b). Based on the ratios between the main and rarer isotopologues, Meier et
al. (2015) found that CO(1-0) and HCN(1-0) emission is on average optically
thick in the central kpc-scale region of NGC 253. As 13C-bearing molecular
species, 13CO(1-0) emission should be moderately optically thin (Martín et
al., 2019). The different extents towards the SW streamer of CO(1-0) and
HCN(1-0) emission from the 13CO(1-0) emission imply that the optical depths of
CO(1-0) and HCN(1-0) emission in the SW streamer region may be lower than that
in the gas disk. Moreover, the fainter emission of H13CN(1-0), N2H+(1-0) and
CH3OH(2k-1k) in Figs. 1(d), 1(e) and 1(f) share similar extents towards the SW
streamer with 13CO(1-0).

Figure 2: The intensities of molecular lines in the PVDs along the major axis.
The major axis of NGC 253 is marked by the black solid line in Fig. 1(a), on
which the black arrow shows the positive direction of position. (a) CO(1-0)
PVD. The grey contours are drawn at [2, 10] K, and the blue contours are drawn
at [20, 30] K. The white arrow points to the SW streamer, and remains the same
in the following panels. (b) 13CO(1-0) PVD. The grey contours are drawn at
[0.3, 1.2, 2.0, 3.5] K. (c) HCN(1-0) PVD. The grey contours are drawn at [0.3,
1.2, 2.0, 3.5] K. (d) H13CN(1-0) PVD. The grey contours are drawn at [0.06,
0.25, 0.45] K, and the blue contour is drawn at 0.62 K. (e) SiO(2-1) PVD. The
grey contours are drawn at [0.06, 0.25, 0.45] K, and the blue contour is drawn
at 0.62 K. (f) CH3OH(2k-1k) PVD. The grey contours are drawn at [0.12, 0.50,
1.00, 1.50] K.

Fig. 2 shows the intensities of five molecular lines from Fig. 1 and SiO in
the J = 2-1 transition ($\sim$86.847 GHz, see Huang et al. 2023 for its
integrated intensity map) in the PVDs along the major axis with systemic
velocity reserved and distance referenced from the phase center. The CO(1-0)
PVD in Fig. 2(a) mainly follows a rotating pattern, which is consistent with
Figure 5 from Krieger et al. (2019). Another common point is that the CO(1-0)
PVD shows non-disk features towards both redshifted and blueshifted sides. The
strong CO(1-0) emission that is marked by blue contours at [20, 30] K is
located in the region within a galactocentric distance range of [-10, 10]
arcsec and follows a rotating pattern. The white arrow in Fig. 2(a) points to
the SW streamer, and the strongest CO(1-0) emission is located near the base
of this streamer.

Figure 3: The intensities of molecular lines in the PVDs along the SW slice.
The SW slice is marked by the black dashed line in Fig. 1(a), on which the
black arrow shows the positive direction of the position. (a) CO(1-0) PVD. The
grey contours are drawn at $\sim$[1.5, 10, 25] K. The white dashed profiles in
panels (a), (b) and (c) outline the outflow in the SW streamer region. (b)
13CO(1-0) PVD. The grey contours are drawn at $\sim$[0.1, 0.9, 2.3] K. (c)
HCN(1-0) PVD. The grey contours are drawn at $\sim$[0.1, 0.9, 3.2] K. (d)
H13CN(1-0) PVD. The grey contours are drawn at $\sim$[0.03, 0.25] K. (e)
SiO(2-1) PVD. The grey contours are drawn at $\sim$[0.03, 0.25] K. (f)
CH3OH(2k-1k) PVD. The grey contours are drawn at $\sim$[0.03, 0.25] K.

The 13CO(1-0) PVD in Fig. 2(b) is less extended than the CO(1-0) and HCN(1-0)
PVDs, but all have similar strong emission patterns. The H13CN(1-0), SiO(2-1),
and CH3OH(2k-1k) PVDs in the second row of Fig. 2 are the least extended. For
H13CN(1-0) and SiO(2-1) PVDs in Figs. 2(d) and 2(e), the strong emission
marked by the blue contours at 0.62 K is distributed in a few clumpy regions,
among which the strongest two clumps are located near the base of the SW
streamer (indicated by the white arrow). For the CH3OH(2k-1k) PVD in Fig.
2(f), the strongest emission originates from the outskirts of the gas disk,
which agrees with the integrated intensity map in Fig. 1(f).

Fig. 3 shows the intensities in the PVDs along a slice across the SW streamer
(SW slice, black dashed line in Fig. 1a) with systemic velocity reserved and
distance taken zero in the slice center. The slice center is the same as
Figure 1 from Walter et al. (2017), and the position angle of the slice equals
150∘. In Fig. 3(a), the CO(1-0) emission in the PVD is dominated by the disk
component, and is extended towards the SW streamer in the region with
D$\sim$[-15, 0] arcsec and V$\sim$200 km s-1. Such extension outlined by the
white dashed profiles, is consistent with Figure 2 from Walter et al. (2017).
The 13CO(1-0) and HCN(1-0) PVDs in Figs. 3(b) and 3(c) are less extended, but
have similar patterns as the CO(1-0) PVD. The H13CN(1-0), SiO(2-1) and
CH3OH(2k-1k) PVDs in the second row of Fig. 3 are least extended, where the
molecular emission only extends at the base of the SW streamer with
D$\sim$[-5, 0] arcsec and V$\sim$200 km s-1. The intensities of N2H+(1-0) in
the PVDs along the major axis and along the SW slice are also checked, which
show similar patterns to the H13CN(1-0) PVDs.

###  3.2 Gas kinematics

The CO(1-0) line is generally the best tracer of total molecular gas content
thanks to its high abundance and low critical density for excitation. Fig. 4
shows the kinematics of the molecular gas in NGC 253 traced with CO(1-0).
There are two interesting kinematic features in Figs. 4(a) and 4(b): (1) the
velocity field showing gradients along both major and minor axes; (2) the
blueshifted velocity and high velocity dispersion in the SW streamer region.

The white contours showing the velocity gradient in the center within a radius
of $\sim$100 pc in Fig. 4(a) have the same pattern as the CO velocity field
from Das et al. (2001) and Krieger et al. (2019). Their studies showed that
the molecular bar shares the same direction as the stellar bar in the near-
infrared band, and tilts $\sim$18∘ with respect to the major axis (Peng et
al., 1996; Iodice et al., 2014). The direction of molecular and stellar bars
is shown as a white solid line in Fig. 4(a). There is a velocity gradient
(white contours) in Fig. 4(a) along the direction of the molecular bar (white
solid line). Yet the ionized bar traced with the H92$\alpha$ recombination
lines is almost parallel to the major axis (Cohen et al., 2020). Das et al.
(2001) explained this phenomenon as different velocities for the gas moving in
different bar orbits and suggested that the perturbation in the gas velocity
field in NGC 253 is due to an accretion event occurred $\sim$10 Myr ago.

The gas velocity in the SW streamer region (pointed out by the black arrow) is
blueshifted ($\sim\rm 200\,km\,s^{-1}$) as shown in Fig. 4(a), which means the
outflowing gas is on the approaching side of NGC 253. The blueshifted outflow
is visibly distinguishable from the redshifted ($\sim\rm 380\,km\,s^{-1}$)
disk rotation. Meanwhile, the other outflow streamers are indistinguishable in
Fig. 4(a), because they are fainter than SW streamer and/or locate behind the
galactic disk (Bolatto et al., 2013). In Fig 4(b), the gas in the SW streamer
region presents high velocity dispersion, which can be contributed by
blueshifted outflow and redshifted rotation.

To compare the kinematics between the outflow and disk, we decompose two
components from the CO(1-0) spectra in the SW streamer region. From a circular
region in the SW streamer region having the same area as the beam size (the
black empty circle in Fig. 4b), we extract the averaged CO(1-0) spectrum,
which is displayed by the black solid profile in Fig. 4(c). The averaged
CO(1-0) spectrum shows a double-peak structure, where we take the velocity
$\sim$300 km s-1 as a reference velocity for the emission near the base of the
SW streamer. Combining the velocity field of CO(1-0) in Fig. 4(a), the
blueshifted component relative to the reference velocity is emitted by the SW
streamer, and the redshifted component is emitted by the gas disk. We also
extract the averaged CO(1-0) spectrum from a region with high-velocity-
dispersion (high-$\sigma$) in the SW streamer region (the black empty ellipse
in Fig. 4b), which is displayed by the black solid profile in Fig. 4(d) also
showing a double peak structure. Using Python-based tool curve_fit, we fit the
CO(1-0) spectra in Figs. 4(c) and 4(d) with double Gaussian function. The blue
dashed profiles show the blueshifted Gaussian components, the red dashed
profiles show the redshifted Gaussian components and the green solid profiles
show the superpositions. The velocity and full width at half maximum (FWHM) in
the velocity space of each Gaussian component are listed in Table 1.

Figure 4: Kinematic features of CO(1-0). (a) Velocity field of CO(1-0). The grey contours are the same as in Fig. 1(a). The white contours mark the [175, 300] km s-1 velocity range with a step of 25 km s-1 in the central $\sim$100 pc. The white solid line shows the direction of the molecular bar. The black solid and dashed lines are the same as Fig. 1(a). The black arrow is the same as the white one in Fig. 1(a), and indicates the same region in panel (b) and the following figures. (b) Velocity dispersion field of CO(1-0). The grey contours are the same as in panel (a). The black empty ellipse outlines the high-velocity-dispersion region in the SW streamer region. The black empty circle inside the ellipse marks a beam-size region in the SW streamer region. (c) The averaged CO(1-0) spectrum from the black empty circle in panel (b). The black profile shows the CO(1-0) line. The blue dashed profile shows the blueshifted Gaussian component. The red dashed profile shows the redshifted Gaussian component. The green solid profile is the superposition. (d) The averaged CO(1-0) spectrum from the black empty ellipse in panel (b). The color codes are the same as panel (c). Table 1: Results of double Gaussian fits for the CO(1-0) spectra. Region | VStreamer (km s-1) | FWHMStreamer (km s-1) | VDisk (km s-1) | FWHMDisk (km s-1) | FWHMStreamer/FWHMDisk
---|---|---|---|---|---
Beam-size | 235.5 $\pm$ 0.2 | 81.8 $\pm$ 0.6 | 386.9 $\pm$ 0.2 | 87.5 $\pm$ 0.5 | 0.93
High-$\sigma$ | 219.7 $\pm$ 0.5 | 118.4 $\pm$ 1.3 | 373.9 $\pm$ 0.2 | 82.4 $\pm$ 0.5 | 1.44

We find that the blueshifted and redshifted components of the averaged CO(1-0)
spectrum from the beam-size region (the black empty circle in Fig. 4b) have
similar FWHMs in the velocity space (line widths). The redshifted component of
the averaged CO(1-0) spectrum from the high-$\sigma$ region (the black empty
ellipse in Fig. 4b) has similar FWHM to the components from the beam-size
region, while the blueshifted component is $\sim$50% wider than the redshifted
component. Moreover, there is a gradual velocity shift along the SW streamer
in the CO(1-0) PVD outlined by white dashed profiles in Fig. 3(a). The
velocity gradient inside the SW streamer can be evidence for an inside-out
acceleration on the gas velocity, i.e., the gas inside the SW streamer is
accelerated as it outflows. It can also be attributed to the ability for the
fast ejecta to reach farther away than the slow ejecta. Those two
possibilities were also discussed in Walter et al. (2017). If we take the
velocity center of the blueshifted component from the high-$\sigma$ region
($\sim$219.7 km s-1) as the averaged velocity for the SW streamer, then the
projected local velocity of the SW streamer equals $\sim$80 km s-1.
Considering the inclination angle of the gas disk ($\sim$78.5∘) and the
outflow being perpendicular to the disk, the deprojected local velocity of the
SW streamer approaches $\sim$400 km s -1, which is comparable with the outflow
velocity in the ionized gas phase (Westmoquette et al., 2011).

###  3.3 Optical depth

Figs. 1(a) and 1(b) show the contours of CO(1-0) emission being more extended
towards the SW streamer region than 13CO(1-0) emission, which can be a result
of the different optical depths. Fig. 5(a) shows the integrated intensity
ratio map of CO/13CO(1-0) based on self-calibration, where we masked the
regions with ratio S/N being less than $3\sigma$. The CO/13CO(1-0) ratio
increases in the SW, southeastern (SE) and northwestern (NW) directions (named
after Bolatto et al. 2013), which are pointed out by the black/grey arrows and
are perpendicular to the gas disk (Fig. 5). The black contour with
CO/13CO(1-0) ratio of 15 defines the boundary between the outflow streamers
and the disk. The CO/13CO(1-0) ratio further increases in part of the outflow
region and is outlined by the blue contour at 21, where the CO/13CO(1-0) ratio
is highest in the SW streamer region. CO is generally optically thick, while
its isotopologue 13CO is optically thinner. The molecular gas in the CMZ of
NGC 253 shows high kinetic temperature in the previous studies (Paglione et
al., 2004; Sakamoto et al., 2011), which weakens the fractionation effects on
C-bearing species (Colzi et al., 2020). Moreover, 13C preferentially forms in
the center, while the C and 13C can be well mixed in the outflow region. The
CO/13CO(1-0) ratio is expected to be similar to the C/13C isotope ratio in the
optically thin limit, and decreases as the optical depth increases.

Figure 5: Integrated intensity ratio maps. (a) Integrated intensity ratio map
of CO/13CO(1-0). The ratios are taken in $\rm K\,km\,s^{-1}$ units, and remain
the same in the following integrated intensity ratio plots. The black contour
is drawn at 15, and the blue contour is drawn at 21. The grey arrows point to
the outflow in the SE and NW directions. (b) Integrated intensity ratio map of
HCN/13CO(1-0). The black contour is drawn at 1.4, and the blue contour is
drawn 2. (c) Integrated intensity ratio map of HCN/CO(1-0). The black solid
line marks the major axis. The black contour is drawn at 0.09, and the blue
contour is drawn at 0.15. The white contour is drawn at CO/13CO(1-0) ratio of
15. The orange, purple, and yellow stars mark the positions at -1.5′′, -4′′
and -6.5′′ offsets along the SW slice. The labels ‘3’, ‘6’, and ‘7’ mark three
GMCs. The grey arrows point to the outflow in the SE and NW directions.

Martín et al. (2019) estimated a value of the isotopic ratio
C/13C$\sim$21$\pm$6 via the integrated intensity ratio between C18O(1-0) and
13C18O(1-0). In Fig. 5(a), the observed CO/13CO(1-0) ratio in the gas disk is
lower than the C/13C ratio, which can be contributed by the optically thick CO
emission in the disk. The CO/13CO(1-0) ratio in the SW, SE and NW outflow
regions approaches or exceeds the C/13C ratio, where the fluctuation in
isotopic abundance is weak. Part of the outflow region (outlined by the blue
contour at 21), including the SW streamer region, shows the CO/13CO(1-0) ratio
consistent with the C/13C ratio. Except for NGC 253, Weiß et al. (2005) found
the integrated intensity ratio CO/13CO(1-0) of the prominent molecular
streamers being comparable to that of the starburst disk in M82, which is
different from the increasing CO/13CO(1-0) ratio towards the directions of
outflow in NGC 253.

Figure 6: Intensity ratio in the PVDs. (a) Intensity ratio of CO/13CO(1-0) in
the PVD along the major axis. The black arrow points to the SW streamer, and
remains the same in the following figures. (b) Intensity ratio of CO/13CO(1-0)
in the PVD along the SW slice. The black dashed profiles outline the outflow
in the SW streamer region.

To decompose different velocity components, we plot in Figs. 6(a) and 6(b) the
intensity ratio of CO/13CO(1-0) in the PVDs along the major axis and the SW
slice, which are marked by the black solid and dashed lines in Fig. 1(a). In
Fig. 6(a), the CO/13CO(1-0) ratio increases in the non-disk components and is
highest on the blueshifted side, which implies the decrease in the optical
depth of CO emission in the molecular gas on the approaching side of NGC 253.
This phenomenon is more obvious in the PVD along the SW slice in Fig. 6(b),
where the CO/13CO(1-0) ratio is low along the gas disk and increases inside
the SW streamer with D$\sim$[-15, 0] arcsec and V$\sim$200 km s-1 outlined by
the two black dashed profiles. Combining Figs. 5(a), 6(a), and 6(b), we infer
that the decrease in the optical depth of CO emission happens inside the
molecular outflow including the SW streamer, which can be attributed to the
gas velocity gradient inside the outflow.

Measuring the molecular mass outflow rate is important because it may indicate
the rate at which the fuel for star formation is expelled, which may
significantly suppress star formation activity. The optical depth of CO
emission is a key factor in estimating the molecular mass outflow rate.
Zschaechner et al. (2018) derived the optical depth with the integrated
intensity ratio of CO(2-1)/CO(1-0), and suggested that the majority of the CO
emission is optically thick in the outflow region of NGC 253. However, the
13C-bearing molecular species being optically thinner than the 12C-bearing
ones implies that the CO/13CO(1-0) ratio is a more reliable indicator of CO
optical depth. Hence, the agreement between the CO/13CO(1-0) ratio (Fig. 5a)
and the C/13C ratio (Martín et al., 2019) suggests that the CO emission in a
considerable portion of the outflow region is optically thin in NGC 253. Such
phenomenon supports the estimation from Bolatto et al. (2013) on the total
molecular mass outflow rate $\sim 9\ M_{\odot}\ yr^{-1}$, which is three times
the star formation rate $\sim 3\ M_{\odot}\ yr^{-1}$ (Ott et al., 2005) in NGC
253.

HCN is one of the most abundant high dipole-moment molecules that trace dense
gas. The permanent electric dipole-moment of HCN ($\mu_{e}\sim$ 2.99 D) is
much higher than CO ($\mu_{e}\sim$ 0.11 D) (Ebenstein & Muenter, 1984;
Goorvitch, 1994), which makes the critical density of HCN three orders of
magnitude higher than that of CO. Fig. 5(b) shows the integrated intensity
ratio map of HCN/13CO(1-0). The HCN/13CO(1-0) ratio is lowest in the gas disk,
increases towards the outflow directions, and becomes the highest in parts of
the outflow region (outlined by the blue contour at a ratio of 2) including
the SW streamer region. Such a trend reminds us of the increasing CO/13CO(1-0)
ratio in the outflow region with an increasing distance from the major axis in
Fig. 5(a). Given that HCN(1-0) emission is optically thick in the central
$\sim$kpc region of NGC 253 (Meier et al., 2015), the increased line width
attributed to gas velocity gradient in the SW streamer region not only can
decrease the optical depth of CO(1-0) emission tracing molecular gas, but also
can decrease the optical depth of HCN(1-0) emission tracing dense gas.

###  3.4 Dense gas fraction

The HCN/CO(1-0) ratio is widely used as an indicator of the fraction of dense
gas, which is immediately responsible for the star formation inside galaxies
(Gao & Solomon, 2004a, b; Lada et al., 2012). Tanaka et al. (2024) presented
non-LTE analyses with ALCHEMI data for the CMZ of NGC 253, and compared the
dense gas fraction in this region with the center of the Milky Way. The
difference turned out to be consistent with that in the HCN/CO(1-0) ratio
between the CMZ of NGC 253 and the Galactic center, which confirms the
HCN/CO(1-0) ratio as a good measurement of the mass fraction of the dense gas
to the entire molecular gas.

Fig. 5(c) shows the integrated intensity ratio map of HCN/CO(1-0). The
HCN/CO(1-0) ratio is the highest in three clumpy regions along the major axis
(black solid line), which are marked by the blue contour at a ratio of 0.15.
The physical sizes of the three regions in Fig. 5(c) are $\lesssim$100 pc,
which fit the scales of giant molecular clouds (GMCs), and can be further
resolved into star-forming clumps ($\sim$10 pc) (Ando et al., 2017). The
positions of those GMCs are consistent with the observation in Leroy et al.
(2015), which were numbered 3, 6, and 7. The bar of NGC 253 traced with
ionized gas is also along the major axis (Das et al., 2001). The inflowing gas
along the bar interacting with the gas in the disk can trigger star formation
inside those GMCs. The existence of supernova remnants and Hii region in the
inner 200 pc of NGC 253 has been revealed about thirty years ago (Ulvestad &
Antonucci, 1997). The forming super star clusters and the winds from young
clusters were detected in high spatial resolution (Leroy et al., 2018; Levy et
al., 2021, 2022). The GMCs 3, 6, and 7 in Fig. 5(c) are located at the base of
outflow streamers outlined by the white contour, and the black contour at 0.09
is extended towards the directions of outflow streamers. Given that the
dynamical age of the SW streamer equaling $\sim$1 Myr is short (Walter et al.,
2017), the star formation inside the GMCs can be the engine for the outflow
streamers (Krieger et al., 2019).

Walter et al. (2017) measured the ratio of peak intensities between HCN(1-0)
and CO(1-0), which equals $\sim$1/10 both in the SW streamer region and the
starburst center of NGC 253. The orange, purple, and yellow stars in Fig. 5(c)
mark the corresponding positions at -1.5′′, -4′′ and -6.5′′ offsets from the
slice center along the SW slice following Walter et al. (2017). Based on the
spatially resolved map, we can estimate the integrated intensity ratios at the
three stars to be $\sim$0.19, 0.13, and 0.08. It is understandable that the
dense gas fraction is the highest inside the GMC, and monotonously decreases
away from the GMC.

Figure 7: First row: PVDs along the major axis. Second row: PVDs along the SW
slice. (a)(d) HCN/CO(1-0) ratio in PVDs. The black dashed profiles in panel
(d) outline the outflow in the SW streamer region. (b)(e) H13CN/13CO(1-0)
ratio in PVDs. In panel (e), the black ellipse named region A outlines the gas
with D$\sim$0 arcsec and V$\sim$300 km s -1, the black ellipse named region B
outlines the gas with D$\sim$[-5, 0] arcsec and V$\sim$200 km s-1, and keep
the same in the following figures. (c)(f) N2H+/13CO(1-0) ratio in PVDs.

Figs. 7(a) and 7(d) show the intensity ratio of HCN/CO(1-0) in the PVDs along
the major axis and the SW slice. In Fig. 7(a), the HCN/CO(1-0) intensity ratio
increases at the bottom of the PVD where the blueshifted non-disk component is
located. In Fig. 7(d), the HCN/CO(1-0) ratio is highest in the SW streamer
region with D$\sim$[-5, 0] arcsec and V$\sim$200 km s -1. Figs. 4(c) and 4(d)
show that the averaged projected velocity of the outflowing gas in the SW
streamer region is around 200 km s -1. Hence, the HCN/CO(1-0) ratio increases
in the gas at the base of the SW streamer that shares the same velocity with
the molecular outflow. The HCN/CO(1-0) ratio in the extended streamer region
with D$\sim$[-15, -5] is moderate as in the gas disk (Fig. 7d), which agrees
with the tendency at -6.5′′ offset in Figure 7 from Walter et al. (2017).

It is possible to detect strong molecular emission at densities below the
critical density. The effective excitation density is defined as the density
at which the integrated intensity of a molecular line equals 1 K km s-1 with
reasonable assumptions about the column density and kinetic temperature.
Shirley (2015) calculated the optically thin critical densities and effective
excitation densities at assumed column density and kinetic temperature for the
common dense gas tracers. The kinetic temperatures are assumed to be $\sim$100
K (Mangum et al., 2019) in the following comparisons on the effective
excitation densities between different tracers. Shirley (2015) calculated the
critical density and effective excitation density of HCN(1-0) at a column
density of 1014 cm-2, to be $1.1\times 10^{5}$ cm-3 and $1.7\times 10^{3}$
cm-3, respectively.

To confirm the distribution of the dense gas in the CMZ of NGC 253, we turn to
tracers of dense gas with lower opacity. The critical density and effective
excitation density of H13CN(1-0) at a column density of 1012.3 cm-2, equal
$9.7\times 10^{4}$ cm-3 and $6.5\times 10^{4}$ cm-3. Even though the higher
effective excitation density originates from a lower reference column density
compared with that of HCN (Shirley, 2015), H13CN(1-0) as an isotopologue can
trace the dense gas with less limitation of optical depth. Fig. 7(b) shows the
intensity ratio of H13CN/13CO(1-0) in the PVD along the major axis, where the
ratio also increases in the blueshifted non-disk component. Although the
H13CN/13CO(1-0) ratio in the PVD along the SW slice in Fig. 7(e) is not as
extended as HCN/CO(1-0), we still observe high ratios in two regions. One is
the ellipse named region A with D$\sim$0 arcsec and V$\sim$300 km s -1
(reference velocity) that represents the gas inside GMC 3. The other is the
ellipse named region B with D$\sim$[-5, 0] arcsec and V$\sim$200 km s-1
(outflow velocity) that represents the gas at the base of the SW streamer.

The critical density and effective excitation density at a column density of
1013 cm-2 of another dense gas tracer N2H+(1-0) equal $2\times 10^{4}$ cm-3
and $2.6\times 10^{3}$ cm-3. Even though the critical density of N2H+(1-0) is
lower than HCN(1-0) and H13CN(1-0), Galactic parsec-scale observations show
that the N2H+(1-0) emission is exclusively associated with rather dense gas
(Kauffmann et al., 2017; Tafalla et al., 2021). The intensity ratio of
N2H+/13CO(1-0) in the PVD along the major axis (Fig. 7c) has an accordant
pattern with H13CN/13CO(1-0) PVD (Fig. 7b). In the PVD along the SW slice
(Fig. 7f), the N2H+/13CO(1-0) ratio increases at the base of the SW streamer
(region B). On one hand, the increased patterns in Figs. 7(d) and 7(f) are
highly consistent. Given that N2H+(1-0) is not safely optically thin and the
optical depths of HCN(1-0) and CO(1-0) can be different in the central region,
the N2H+/13CO(1-0) and HCN/CO(1-0) ratios inside GMC 3 may be affected by
optical depths. On the other hand, all three ratios in Figs. 7(d), 7(e) and
7(f) are higher in region B, which implies a negligible influence of optical
depths at the base of the SW streamer, not to mention in the extended streamer
region (Fig. 6b).

In summary, we find that the dense gas fraction is high inside GMC 3 (with
reference velocity, region A) and at the base of the SW streamer (with outflow
velocity, region B). The dense gas fraction in the extended streamer region
(outlined by the black dashed profiles) is moderate in Fig. 7(d) and doesn’t
show visible signs of accumulation of dense gas traced with HCN/CO(1-0).
Combining the low optical depths of CO(1-0) and HCN(1-0) emission in the
extended streamer region (Figs. 5a and 5b), we suggest the existence of gas
velocity gradient prevents the gas from accumulation inside the SW streamer of
NGC 253.

###  3.5 Shock strength

Si-bearing material and solid-phase methanol are located in different parts of
the dust grains. Silicon tends to reside in the core of grains and can be
sputtered to the gas phase by fast shocks (v${}_{\rm shock}\gtrsim$ 15-20 km
s-1, Kelly et al. 2017). Gas-phase silicon reacts with molecular oxygen or a
hydroxyl radical and forms SiO (Schilke et al., 1997). On the other hand,
solid-phase methanol is in the icy mantle of grains. Slow shocks (v${}_{\rm
shock}\lesssim$ 15-20 km s-1, Nesterenok 2022) are able to impact the icy
mantle and inject gas-phase methanol into the ISM without destroying the grain
core nor the methanol molecule (Millar et al., 1991; Charnley et al., 1995).
Huang et al. (2023) studied the shock tracers SiO and HNCO with the ALCHEMI
data, and revealed a picture that most of the GMCs are subjected to shocks.
They declared that HNCO may not be a unique tracer of slow shocks, while a
high abundance of silicon in the gas phase can only be explained by fast
shocks. In this section, we will explore the methanol emission to verify the
state of slow shocks, and inspect the relation between fast shocks traced by
SiO emission and outflow streamers.

Figure 8: Integrated intensity ratio maps. (a) Integrated intensity ratio map
of SiO(2-1)/CO(1-0). The black contour is drawn at 0.003, and the blue contour
is drawn at 0.01. The labels ‘3’, ‘6’, and ‘7’ mark the three GMCs. (b)
Integrated intensity ratio map of SiO(2-1)/13CO(1-0). The black contour is
drawn at 0.04, and the blue contour is drawn at 0.1. The grey contour is drawn
at CO/13CO(1-0) ratio of 15. (c) Integrated intensity ratio map of
CH3OH(2k-1k)/13CO(1-0). The black contour is drawn at 0.15, and the blue
contour is drawn at 0.4. (d) Integrated intensity ratio map of
SiO(2-1)/CH3OH(2k-1k). The black contour is drawn at 0.5, and the blue contour
is drawn at 1.

Fig. 8(a) shows the integrated intensity ratio map of SiO(2-1)/CO(1-0). The
three dense GMCs (labeled by ‘3’, ‘6’, ‘7’ following Fig. 5c) also show
increased SiO(2-1)/CO(1-0) ratio that is marked by the blue contour at 0.01.
Moreover, two SiO(2-1)/CO(1-0) ratio peaks appear around GMC 3 in Fig. 8(a),
with the second one consistent with the position of GMC 4 in Leroy et al.
(2015). To examine the influence of optical depth, we also plot the integrated
intensity ratio map of SiO(2-1)/13CO(1-0) in Fig. 8(b). The increased
SiO(2-1)/13CO(1-0) ratio exists in the four GMCs (marked by blue contour at
0.1). The increased ratios in the GMCs imply the enhanced fast shocks there.
The grey contour in Fig. 8(b) outlines the outflow streamers that originate
from the four GMCs. Both the distributions of SiO(2-1)/CO(1-0) (marked by
black contour at 0.003) and SiO(2-1)/13CO(1-0) (marked by black contour at
0.04) are extended towards the outflow streamers, which connect the fast
shocks with the formation of the molecular outflow.

Figure 9: First row: PVDs along the major axis. Second row: PVDs along the SW
slice. (a)(d) SiO(2-1)/13CO(1-0) ratio in PVDs. (b)(e) CH3OH(2k-1k)/13CO(1-0)
ratio in PVDs. (c)(f) SiO(2-1)/CH3OH(2k-1k) ratio in PVDs.

Fig. 8(c) shows the integrated intensity ratio map of CH3OH(2k-1k)/13CO(1-0).
The increased CH3OH(2k-1k)/13CO(1-0) ratio that is marked by black contour at
0.15 located on the outskirts of the gas disk. Locations with high
CH3OH(2k-1k)/13CO(1-0) ratio are similar to those with high CH3OH(2k-1k)
integrated intensity shown in Fig. 1(f). Such a tendency is in agreement with
HNCO 40,4-30,3 emission in Huang et al. (2023) that shows high integrated
intensity in the outermost CMZ of NGC 253. The strong methanol emission on the
outskirts, with positions in consistent with the methanol masers found by
Gorski et al. (2017), implies the frequent occurrence of slow shocks.
Meanwhile, the weak methanol emission in the center could be a result of the
infrequent slow shocks, or the depletion of methanol molecules (Hartquist et
al., 1995; Ellingsen et al., 2017) by fast shocks or photo-dissociation due to
intense star formation (Humire et al., 2022). The non-cospatial distributions
of fast and slow shocks were also found in another nearby galaxy NGC 1068
(Kelly et al., 2017) and imply that fast shocks do not necessarily occur with
slow shocks. To obtain the relative strength between the fast and slow shocks,
we plot the integrated intensity ratio map of SiO(2-1)/CH3OH(2k-1k) in Fig.
8(d), where the central region turns out to be dominated by the fast shocks.

Fig. 9 shows intensity ratios of the shock tracers in PVDs, where the top
panels show the PVDs along the major axis and the bottom panels show the PVDs
along the SW slice. The asymmetric and increased pattern (including regions A
and B) of SiO(2-1)/13CO(1-0) ratio in PVDs (Figs. 9a and 9d) are highly
similar to the H13CN/13CO(1-0) ratio in PVDs (Figs. 7b and 7e). We will study
in detail the correlation between the dense gas fraction and the strength of
fast shocks in the next section. The increased intensity ratio of
CH3OH(2k-1k)/13CO(1-0) on the outskirts of the gas disk in Fig. 9(b) is
coherent with the distribution of the increased integrated intensity ratio in
Fig. 8(c). It is worth noting that a few red pixels located in region B of
Fig. 9(e) indicate an enhancement of the slow shocks at the base of the SW
streamer. In Figs. 9(c) and 9(f), we plot the SiO(2-1)/CH3OH(2k-1k) ratio in
PVDs. The region A in Fig. 9(f) is dominated by the fast shocks, while the
SiO(2-1)/CH3OH(2k-1k) ratio in region B is not so high as region A. We infer
that the fast shocks are triggered by the star formation inside GMC, which
supports the second scenario for the origin of shocks in Huang et al. (2023).
Meanwhile, the fast and slow shocks co-exist in the SW streamer, and can
further be related to the formation of the molecular outflow.

##  4 Formation of the molecular outflow

Walter et al. (2017) presented molecular spectra in the SW streamer region and
quantified the relationship between the different molecular transitions,
including CO(1-0) and HCN(1-0). Based on the data from the ALCHEMI survey, we
can add several weaker molecular lines. In Fig. 10(a), we plot orange, purple,
and yellow stars the same as Fig. 5(c) onto the self-calibrated integrated
intensity ratio map of CO/13CO(1-0), where the increased ratio indicated by
red color represents the outflow region. The purple star is located at -4′′
offset on the SW slice following Walter et al. (2017). The molecular emission
shown in Fig. 10(b) are averaged intensity profiles from a beam-size region
centering on this purple star. The black profile shows the averaged intensity
of the CO(1-0) spectrum (scaled down by a factor of 100), the red and grey
profiles show the averaged intensities of HCN(1-0) and 13CO(1-0) spectra
(scaled down by a factor of 10), as well as pink, orange, blue and green
profiles, show the averaged intensities of H13CN(1-0), N2H+(1-0), SiO(2-1) and
CH3OH(2k-1k) spectra.

Figure 10: Molecular spectra in the SW streamer region. (a) Integrated
intensity ratio map of CO/13CO(1-0). The black contour is drawn at 15, and the
blue contour is drawn at 21. The orange, purple, and yellow stars are the same
as Fig. 5(c). (b) Averaged molecular spectra extracted from a beam size region
centering at the purple star in panel (a). The color codes are labeled in the
top left. The averaged intensity of the CO(1-0) line is scaled down by a
factor of 100. The averaged intensities of 13CO(1-0) and HCN(1-0) lines are
scaled down by a factor of 10. Table 2: Peak intensities from double Gaussian
fits for molecular spectra from the SW streamer region (Fig. 10b).

Line | Streamer (mJy beam-1) | Disk (mJy beam-1) | Streamer-to-disk ratio
---|---|---|---
CO(1-0) | 161.75 $\pm$ 1.43 | 367.44 $\pm$ 1.47 | 0.440 $\pm$ 0.004
13CO(1-0) | 7.54 $\pm$ 0.36 | 25.17 $\pm$ 0.38 | 0.300 $\pm$ 0.015
HCN(1-0)* | 12.60 $\pm$ 1.29 | 22.02 $\pm$ 1.43 | 0.572 $\pm$ 0.069
H13CN(1-0)* | 0.64 $\pm$ 0.10 | 1.22 $\pm$ 0.09 | 0.524 $\pm$ 0.091
N2H+(1-0)* | 1.07 $\pm$ 0.05 | 2.14 $\pm$ 0.05 | 0.500 $\pm$ 0.028
SiO(2-1) | 0.27 $\pm$ 0.42 | 0.93 $\pm$ 0.41 | 0.290 $\pm$ 0.471
CH3OH(2k-1k) | 0.76 $\pm$ 0.09 | 2.66 $\pm$ 0.10 | 0.286 $\pm$ 0.034

  * *

The star symbols in the first column mark the dense gas tracers.

All the molecular lines in Fig. 10(b) show double peak structures, which are
in agreement with Figure 7 from Walter et al. (2017) and CO(3-2) lines in
Figure 4 from Levy et al. (2022). The relatively weak blueshifted component is
emitted by the SW streamer, and the relatively strong redshifted component is
emitted by the gas disk. We use the Python-based tool curve_fit to conduct a
double Gaussian fit on each molecular line. The peak intensities together with
errors from the Gaussian fit of the streamer and disk components for each
molecular line are listed in the second and third columns of Table 2, and the
corresponding line name is listed in the first column. Although the streamer
component is weak in the SiO(2-1) line, it exists in all the molecular lines
that are extracted from the beam-size region (purple star in Fig. 10a) in the
SW streamer.

As described in the introduction, there are three main formation scenarios
previously proposed for the molecular outflow. One is the molecular outflow
directly driven by radiation and/or pressure, one is the molecular cloud
entrained by the hot wind, and the other is the molecular outflow in situ
forming from the hot wind. The last scenario assumes the cooling timescale to
be shorter than the dynamical timescale (Efstathiou, 2000; Silich et al.,
2003). Walter et al. (2017) estimated that the molecular gas inside the SW
streamer of NGC 253 is ejected from the disk $\sim$1 Myr ago, which is too
short to create new CO molecules that can emit (Clark et al., 2012). Moreover,
the phenomenon that all the molecular lines in Table 2 show streamer features
indicates that the SW streamer starts as dense, shocked, and chemically rich
outflowing gas. The existence of dense gas tracers, e.g., HCN, H13CN and N2H+,
inside the outflow further excludes the possibility of in situ formation.
However, it is still difficult to distinguish whether the outflowing molecular
gas is directly driven by radiation/pressure or is entrained by the hot wind.

###  4.1 Correlation between the dense gas fraction and the strength of fast
shocks

Fig. 11(a) shows the integrated intensity ratio map of H13CN/13CO(1-0), where
the GMCs have the highest ratios marked by the blue contour at 0.12, and the
black contour at 0.06 is extended towards the outflow streamers. Such
increased patterns are similar to that of SiO(2-1)/13CO(1-0) in Fig. 8(b).
Moreover, the increased H13CN/13CO(1-0) and SiO(2-1)/13CO(1-0) ratios in PVDs
in Figs. 7 and 9 also share similar patterns. The consistency in the enhanced
emissions of H13CN(1-0) and SiO(2-1) implies a positive correlation between
the dense gas fraction and the strength of fast shocks in NGC 253. We regrid
the H13CN(1-0), SiO(2-1) and 13CO(1-0) data cubes to a 1.6′′ beam-size to
avoid oversampling, then statistically quantify the linear relationship
between H13CN/13CO(1-0) and SiO(2-1)/13CO(1-0).

Figure 11: Correlation between the dense gas fraction and the strength of fast
shocks. (a) Integrated intensity ratio map of H13CN/13CO(1-0). The black
contour is drawn at 0.06, and the blue contour is drawn at 0.12. The labels
‘3’, ‘6’, and ‘7’ mark the same GMCs as Fig. 5(c). (b) Correlation between
H13CN/13CO(1-0) and SiO(2-1)/13CO(1-0). The blue points mark the ratios from
pixels with a spatial resolution equaling the beam size. The grey solid line
shows the result of a linear fit to the blue points.

The blue points in Fig. 11(b) show the H13CN/13CO(1-0) and SiO(2-1)/13CO(1-0)
ratios in each beam size region from Figs 8(b) and 11(a). To check the linear
relationship between two ratios, which is visibly tight, we calculate the
Pearson correlation coefficient using the Python-based package pearsonr. The
correlation coefficient ($r$) equals 0.85, which reveals a tight positive
correlation between H13CN/13CO(1-0) and SiO(2-1)/13CO(1-0) ratios. We perform
a linear fit for the data set, which is shown by the grey solid line in Fig.
11(b). We also calculate the Pearson correlation coefficient between
H13CN(1-0) and SiO(2-1) integrated intensities, which approaches 0.96 and
reveals a tight correlation between the intensities of two lines. Therefore,
the correlation between H13CN/13CO(1-0) and SiO(2-1)/13CO(1-0) ratios is real,
it is not a result from the common denominator. Meanwhile, we keep in mind
that the high critical densities of H13CN(1-0) and SiO(2-1) may contribute to
the tight positive correlation between two ratios. Theoretically, the star
formation that happens in the region with a high dense gas fraction can
trigger fast shocks. Combining the similar extensions of the SiO(2-1) and
H13CN(1-0) emissions towards the outflow in Figs 8(b) and 11(a), the star
formation inside the GMCs can trigger fast shocks and contribute to the
formation of the molecular outflow in NGC 253.

The difference between the tracers of dense gas and shocks presents when we
further analyse the molecular spectra in Fig. 10(b). Given that all the lines
show double peak structures, we can distinguish the contributions of steamer
and disk via the double Gaussian fits. In the last column of Table 2, we list
the streamer-to-disk peak intensity ratios (streamer-to-disk ratios for short)
together with errors for different molecular lines. Each ratio is calculated
by dividing the peak intensity of the streamer component by that of the disk
component.

The streamer-to-disk ratio of CO(1-0) is higher than its isotopologue
13CO(1-0), which can be explained by the lower opacity environment inside the
SW streamer than in the gas disk for CO emission as discussed in Section 3.3.
It is interesting to find that all the dense gas tracers (marked by the star
symbols in the first column of Table 2) have higher streamer-to-disk ratios
than the others. The optical depth decreasing inside the SW streamer for HCN
emission can explain the highest streamer-to-disk ratio of HCN(1-0) among the
dense gas tracers. In addition, the SiO(2-1) and methanol, as the tracers of
fast and slow shocks, respectively, have similar streamer-to-disk ratios to
13CO(1-0). It seems that the enhancement of dense gas fraction is stronger in
the SW streamer than in the gas disk, while the shock strength is equivalently
enhanced in the SW streamer and the gas disk. More detailed clues are needed
to explain such a difference.

We suggest the physical pictures that are related to the SW streamer in NGC
253 as follows: (i) The GMCs (Fig. 5c) in the direction of the major axis are
related to the gas accretion along the bar structure that further provides
material for the star formation. Meanwhile, the star formation inside the GMCs
that are located at the base of the outflow streamers contributes to driving
the molecular outflow, which is in agreement with the model presented by Levy
et al. (2022). (ii) There can be intense star formation inside the GMCs also
presenting as higher dense gas fraction, and results in enhanced shock
strength at the base of the outflow (Figs. 7 and 9) including the SW streamer.
(iii) The optical depths of CO and HCN emission decrease in the SW streamer
(Figs. 5a and 5b), and the dense gas is diluted in the extended SW streamer
region (Fig. 7d), which can be attributed to the gas velocity gradient inside
the molecular outflow.

##  5 Summary

In this work, we analyse data from the ALCHEMI survey and study the physical
properties of the molecular outflow in the starburst galaxy NGC 253.

  * (1)

The emission of CO(1-0), 13CO(1-0), HCN(1-0), H13CN(1-0), N2H+(1-0) and
CH3OH(2k-1k) is extended towards the SW streamer. The CO(1-0) and HCN(1-0)
emission is the most extended, which can be attributed to the optically thick
environments in the gas disk and decreased optical depths towards the SW
streamer.

  * (2)

The molecular outflow in the SW streamer region is blueshifted with a
deprojected local velocity of $\sim$400 km s-1, which is consistent with
previous studies. The wider blueshifted component than the disk component
suggests an inside-out acceleration of molecular outflow or the fast ejecta
getting farther away than the slow ejecta. All the molecular spectra from the
SW streamer region show double peak structures, which indicate that the SW
streamer starts as dense, shocked and chemically rich outflowing gas, rather
than in situ formation from hot wind.

  * (3)

The integrated intensity ratio maps of CO/13CO(1-0) and HCN/13CO(1-0) show
similar patterns. Both are lowest in the gas disk, increase towards outflow
directions that are perpendicular to the gas disk, and become highest in the
SW streamer region. Combining the isotopic ratio of C/13C from a previous
study, we suggest that the CO(1-0) emission is optically thin in the SW
streamer region, which may be a result of the gas velocity gradient.

  * (4)

Three GMCs present in the integrated intensity ratio map of HCN/CO(1-0), which
are aligned with the direction of the major axis and can be caused by the gas
accretion along the bar structure. Those GMCs are located at the base of the
outflow streamers including the SW streamer, where the star formation can
drive molecular outflow. The HCN/CO(1-0), H13CN/13CO(1-0) and N2H+/13CO(1-0)
intensity ratios in PVDs show high dense gas fraction at the base of the SW
streamer. The HCN/CO(1-0) intensity ratio in PVDs suggests moderate dense gas
fraction in the extended streamer region without visible signs of accumulation
of dense gas, which may also be a result of gas velocity gradient.

  * (5)

The SiO(2-1)/CO(1-0) and SiO(2-1)/13CO(1-0) integrated intensity ratios show
enhanced fast shocks in the GMCs. The SiO(2-1)/13CO(1-0) and
CH3OH(2k-1k)/13CO(1-0) intensity ratios in PVDs indicate that fast shocks can
be triggered by the star formation inside GMC, while fast and slow shocks co-
exist at the base of the SW streamer and can be related to the formation of
molecular outflow.

  * (6)

There is a tight positive correlation between the dense gas fraction traced
with H13CN/13CO(1-0) and the strength of fast shocks traced with
SiO(2-1)/13CO(1-0). The dense gas fraction is high in GMCs, where the star
formation can trigger fast shocks and contribute to the formation of molecular
outflow. One difference is that the enhancement of the dense gas fraction is
more tightly related to the SW streamer than the gas disk, while the shock
strength is equivalently enhanced in the SW streamer and the gas disk.

###### Acknowledgements.

We thank the referee for constructive comments that greatly improved this
paper. This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2017.1.00161.L. 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. This work was supported by JSPS KAKENHI Grant Number JP17H06130 and the
NAOJ ALMA Scientific Research Grant Number 2017-06B. M.B. acknowledges support
from the National Natural Science Foundation of China (No. 12303009) and the
China Scholarship Council (No. 202006860042). N.H. acknowledges support from
JSPS KAKENHI Grant Number JP21K03634. F.E. is financially supported by JSPS
KAKENHI Grant Numbers JP17K14259 and JP20H00172. L.C. acknowledges financial
support through the Spanish grant PID2019-105552RB-C41 funded by
MCIN/AEI/10.13039/501100011033. V.M.R. has received support from the project
RYC2020-029387-I funded by MCIN/AEI /10.13039/501100011033. S.V., M.B., and
K-Y. Huang acknowledge support from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme MOPPEX
833460. M.B. acknowledges Prof. Yu Gao for his remarkable discovery of the HCN
J = 1-0 as an indicator of dense molecular gas, as well as for his gentle
clarification to her regarding the ALMA instrument in 2018.

## References

  * Aalto (2015) Aalto, S. 2015, Revolution in Astronomy with ALMA: The Third Year, 499, 85
  * ALMA Partnership et al. (2015) ALMA Partnership, Fomalont, E. B., Vlahakis, C., et al. 2015, ApJ, 808, L1
  * Aladro et al. (2015) Aladro, R., Martín, S., Riquelme, D., et al. 2015, A&A, 579, A101
  * Ando et al. (2017) Ando, R., Nakanishi, K., Kohno, K., et al. 2017, ApJ, 849, 81
  * Banda-Barragán et al. (2019) Banda-Barragán, W. E., Zertuche, F. J., Federrath, C., et al. 2019, MNRAS, 486, 4526
  * Barnes et al. (2020) Barnes, A. T., Kauffmann, J., Bigiel, F., et al. 2020, MNRAS, 497, 1972
  * Behrens et al. (2022) Behrens, E., Mangum, J. G., Holdship, J., et al. 2022, ApJ, 939, 119
  * Bendo et al. (2015) Bendo, G. J., Beswick, R. J., D’Cruze, M. J., et al. 2015, MNRAS, 450, L80
  * Bolatto et al. (2013) Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013, Nature, 499, 450
  * Brunthaler et al. (2009) Brunthaler, A., Castangia, P., Tarchi, A., et al. 2009, A&A, 497, 103
  * Charnley et al. (1995) Charnley, S. B., Kress, M. E., Tielens, A. G. G. M., et al. 1995, ApJ, 448, 232
  * Clark et al. (2012) Clark, P. C., Glover, S. C. O., Klessen, R. S., et al. 2012, MNRAS, 424, 2599
  * Cohen et al. (2020) Cohen, D. P., Turner, J. L., & Consiglio, S. M. 2020, MNRAS, 493, 627
  * Colzi et al. (2020) Colzi, L., Sipilä, O., Roueff, E., et al. 2020, A&A, 640, A51
  * Das et al. (2001) Das, M., Anantharamaiah, K. R., & Yun, M. S. 2001, ApJ, 549, 896
  * Ebenstein & Muenter (1984) Ebenstein, W. L., & Muenter, J. S. 1984, J. Chem. Phys., 80, 3989
  * Efstathiou (2000) Efstathiou, G. 2000, MNRAS, 317, 697
  * Ellingsen et al. (2017) Ellingsen, S. P., Chen, X., Breen, S. L., et al. 2017, MNRAS, 472, 604
  * Endres et al. (2016) Endres, C. P., Schlemmer, S., Schilke, P., et al. 2016, Journal of Molecular Spectroscopy, 327, 95
  * Engelbracht et al. (1998) Engelbracht, C. W., Rieke, M. J., Rieke, G. H., et al. 1998, ApJ, 505, 639
  * Fielding & Bryan (2022) Fielding, D. B. & Bryan, G. L. 2022, ApJ, 924, 82
  * Gao & Solomon (2004a) Gao, Y. & Solomon, P. M. 2004, ApJ, 606, 271
  * Gao & Solomon (2004b) Gao, Y. & Solomon, P. M. 2004, ApJS, 152, 63
  * García-Burillo et al. (2010) García-Burillo, S., Usero, A., Fuente, A., et al. 2010, A&A, 519, A2
  * Goldsmith & Kauffmann (2017) Goldsmith, P. F. & Kauffmann, J. 2017, ApJ, 841, 25
  * Goorvitch (1994) Goorvitch, D. 1994, ApJS, 95, 535
  * Gorski et al. (2017) Gorski, M., Ott, J., Rand, R., et al. 2017, ApJ, 842, 124
  * Greisen (2003) Greisen, E. W. 2003, Information Handling in Astronomy - Historical Vistas, 285, 109
  * Ginsburg et al. (2016) Ginsburg, A., Robitaille, T., & Beaumont, C. 2016, Astrophysics Source Code Library. ascl:1608.010
  * Haasler et al. (2022) Haasler, D., Rivilla, V. M., Martín, S., et al. 2022, A&A, 659, A158
  * Harada et al. (2021) Harada, N., Martín, S., Mangum, J. G., et al. 2021, ApJ, 923, 24
  * Harada et al. (2022) Harada, N., Martín, S., Mangum, J. G., et al. 2022, ApJ, 938, 80
  * Hartquist et al. (1995) Hartquist, T. W., Menten, K. M., Lepp, S., et al. 1995, MNRAS, 272, 184
  * Heckman et al. (2000) Heckman, T. M., Lehnert, M. D., Strickland, D. K., et al. 2000, ApJS, 129, 493
  * Heesen et al. (2011) Heesen, V., Beck, R., Krause, M., et al. 2011, A&A, 535, A79
  * Hollenbach & Salpeter (1971) Hollenbach, D. & Salpeter, E. E. 1971, ApJ, 163, 155
  * Huang et al. (2022) Huang, K.-Y., Viti, S., Holdship, J., et al. 2022, A&A, 666, A102
  * Huang et al. (2023) Huang, K.-Y., Viti, S., Holdship, J., et al. 2023, A&A, 675, A151
  * Holdship et al. (2021) Holdship, J., Viti, S., Martín, S., et al. 2021, A&A, 654, A55
  * Holdship et al. (2022) Holdship, J., Mangum, J. G., Viti, S., et al. 2022, ApJ, 931, 89
  * Humire et al. (2022) Humire, P. K., Henkel, C., Hernández-Gómez, A., et al. 2022, A&A, 663, A33
  * Iodice et al. (2014) Iodice, E., Arnaboldi, M., Rejkuba, M., et al. 2014, A&A, 567, A86
  * Israel (2020) Israel, F. P. 2020, A&A, 635, A131
  * Kahle et al. (2023) Kahle, K. A., Hernández-Gómez, A., Wyrowski, F., et al. 2023, A&A, 673, A143
  * Kauffmann et al. (2017) Kauffmann, J., Goldsmith, P. F., Melnick, G., et al. 2017, A&A, 605, L5
  * Kelly et al. (2017) Kelly, G., Viti, S., García-Burillo, S., et al. 2017, A&A, 597, A11
  * Klein et al. (1994) Klein, R. I., McKee, C. F., & Colella, P. 1994, ApJ, 420, 213
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511
  * Krieger et al. (2019) Krieger, N., Bolatto, A. D., Walter, F., et al. 2019, ApJ, 881, 43
  * Lada et al. (2012) Lada, C. J., Forbrich, J., Lombardi, M., et al. 2012, ApJ, 745, 190
  * Le Bourlot et al. (2012) Le Bourlot, J., Le Petit, F., Pinto, C., et al. 2012, A&A, 541, A76
  * Lefèvre et al. (2020) Lefèvre, C., Kramer, C., Neri, R., et al. 2020, mm Universe @ NIKA2 - Observing the mm Universe with the NIKA2 Camera, 228, 00014
  * Leroy et al. (2015) Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015, ApJ, 801, 25
  * Leroy et al. (2018) Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2018, ApJ, 869, 126
  * Lopez et al. (2023) Lopez, S., Lopez, L. A., Nguyen, D. D., et al. 2023, ApJ, 942, 108
  * Levy et al. (2021) Levy, R. C., Bolatto, A. D., Leroy, A. K., et al. 2021, ApJ, 912, 4
  * Levy et al. (2022) Levy, R. C., Bolatto, A. D., Leroy, A. K., et al. 2022, ApJ, 935, 19
  * Mangum et al. (2019) Mangum, J. G., Ginsburg, A. G., Henkel, C., et al. 2019, ApJ, 871, 170
  * Martín et al. (2019) Martín, S., Muller, S., Henkel, C., et al. 2019, A&A, 624, A125
  * Martín et al. (2021) Martín, S., Mangum, J. G., Harada, N., et al. 2021, A&A, 656, A46
  * McGuire (2022) McGuire, B. A. 2022, ApJS, 259, 30
  * Meier et al. (2015) Meier, D. S., Walter, F., Bolatto, A. D., et al. 2015, ApJ, 801, 63
  * Michiyama et al. (2018) Michiyama, T., Iono, D., Sliwa, K., et al. 2018, ApJ, 868, 95
  * Millar et al. (1991) Millar, T. J., Herbst, E., & Charnley, S. B. 1991, ApJ, 369, 147
  * Murray et al. (2011) Murray, N., Ménard, B., & Thompson, T. A. 2011, ApJ, 735, 66
  * Nesterenok (2022) Nesterenok, A. V. 2022, MNRAS, 509, 4555
  * Orlando et al. (2005) Orlando, S., Peres, G., Reale, F., et al. 2005, A&A, 444, 505
  * Ott et al. (2005) Ott, J., Weiss, A., Henkel, C., et al. 2005, ApJ, 629, 767
  * Paglione et al. (2004) Paglione, T. A. D., Yam, O., Tosaki, T., et al. 2004, ApJ, 611, 835
  * Pantaleone et al. (2021) Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2021, ApJ, 917, 49
  * Peng et al. (1996) Peng, R., Zhou, S., Whiteoak, J. B., et al. 1996, ApJ, 470, 821
  * Pietsch et al. (2000) Pietsch, W., Vogler, A., Klein, U., et al. 2000, A&A, 360, 24
  * Rekola et al. (2005) Rekola, R., Richer, M. G., McCall, M. L., et al. 2005, MNRAS, 361, 330
  * Rieke et al. (1980) Rieke, G. H., Lebofsky, M. J., Thompson, R. I., et al. 1980, ApJ, 238, 24
  * Sakamoto et al. (2011) Sakamoto, K., Mao, R.-Q., Matsushita, S., et al. 2011, ApJ, 735, 19
  * Sánchez-Monge et al. (2018) Sánchez-Monge, Á., Schilke, P., Ginsburg, A., et al. 2018, A&A, 609, A101
  * Schilke et al. (1997) Schilke, P., Walmsley, C. M., Pineau des Forets, G., et al. 1997, A&A, 321, 293
  * Shirley (2015) Shirley, Y. L. 2015, PASP, 127, 299
  * Silich et al. (2003) Silich, S., Tenorio-Tagle, G., & Muñoz-Tuñón, C. 2003, ApJ, 590, 791
  * Socrates et al. (2008) Socrates, A., Davis, S. W., & Ramirez-Ruiz, E. 2008, ApJ, 687, 202
  * Spilker et al. (2019) Spilker, J. S., Bezanson, R., Weiner, B. J., et al. 2019, ApJ, 883, 81
  * Strickland et al. (2000) Strickland, D. K., Heckman, T. M., Weaver, K. A., et al. 2000, AJ, 120, 2965
  * Strickland et al. (2002) Strickland, D. K., Heckman, T. M., Weaver, K. A., et al. 2002, ApJ, 568, 689
  * Sturm et al. (2011) Sturm, E., González-Alfonso, E., Veilleux, S., et al. 2011, ApJ, 733, L16
  * The CASA Team et al. (2022) CASA Team et al. 2022, PASP, 134, 114501
  * Tacchella et al. (2015) Tacchella, S., Carollo, C. M., Renzini, A., et al. 2015, Science, 348, 314
  * Tacchella et al. (2016) Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 458, 242
  * Tafalla et al. (2021) Tafalla, M., Usero, A., & Hacar, A. 2021, A&A, 646, A97
  * Takano et al. (2019) Takano, S., Nakajima, T., & Kohno, K. 2019, PASJ, 71, S20
  * Tanaka et al. (2024) Tanaka, K., Mangum, J. G., Viti, S., et al. 2024, ApJ, 961, 18
  * Tang et al. (2019) Tang, X. D., Henkel, C., Menten, K. M., et al. 2019, A&A, 629, A6
  * Turner (1985) Turner, B. E. 1985, ApJ, 299, 312
  * Uhlig et al. (2012) Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, MNRAS, 423, 2374
  * Ulvestad & Antonucci (1997) Ulvestad, J. S. & Antonucci, R. R. J. 1997, ApJ, 488, 621
  * Veilleux et al. (2005) Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769
  * Veilleux et al. (2009) Veilleux, S., Kim, D.-C., Rupke, D. S. N., et al. 2009, ApJ, 701, 587
  * Veilleux et al. (2020) Veilleux, S., Maiolino, R., Bolatto, A. D., et al. 2020, A&A Rev., 28, 2
  * Walter et al. (2017) Walter, F., Bolatto, A. D., Leroy, A. K., et al. 2017, ApJ, 835, 265
  * Weiß et al. (2005) Weiß, A., Walter, F., & Scoville, N. Z. 2005, A&A, 438, 533
  * Westmoquette et al. (2011) Westmoquette, M. S., Smith, L. J., & Gallagher, J. S. 2011, MNRAS, 414, 3719
  * Zschaechner et al. (2018) Zschaechner, L. K., Bolatto, A. D., Walter, F., et al. 2018, ApJ, 867, 111