# The ALCHEMI atlas: principal component analysis reveals starburst evolution
in NGC 253

Nanase Harada National Astronomical Observatory of Japan, 2-21-1 Osawa,
Mitaka, Tokyo 181-8588, Japan  Institute of Astronomy and Astrophysics,
Academia Sinica, 11F of AS/NTU Astronomy-Mathematics Building, No.1, Sec. 4,
Roosevelt Rd, Taipei 10617, Taiwan  Astronomical Science Program, Graduate
Institute for Advanced Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo
181-1855, Japan  David S. Meier New Mexico Institute of Mining and Technology,
801 Leroy Place, Socorro, NM 87801, USA  National Radio Astronomy Observatory,
PO Box O, 1003 Lopezville Road, Socorro, NM 87801, USA  Sergio Martín European
Southern Observatory, Alonso de Córdova, 3107, Vitacura, Santiago 763-0355,
Chile  Joint ALMA Observatory, Alonso de Córdova, 3107, Vitacura, Santiago
763-0355, Chile  Sebastien Muller Department of Space, Earth and Environment,
Chalmers University of Technology, Onsala Space Observatory, SE-439 92 Onsala,
Sweden  Kazushi Sakamoto Institute of Astronomy and Astrophysics, Academia
Sinica, 11F of AS/NTU Astronomy-Mathematics Building, No.1, Sec. 4, Roosevelt
Rd, Taipei 10617, Taiwan  Toshiki Saito National Astronomical Observatory of
Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Mark D. Gorski Department
of Space, Earth and Environment, Chalmers University of Technology, Onsala
Space Observatory, SE-439 92 Onsala, Sweden  Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA) and Department of Physics and
Astronomy, Northwestern University, Evanston, IL 60208, USA  Christian Henkel
Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
Astronomy Department, Faculty of Science, King Abdulaziz University, P. O. Box
80203, Jeddah 21589, Saudi Arabia  Kunihiko Tanaka Department of Physics,
Faculty of Science and Technology, Keio University, 3-14-1 Hiyoshi, Yokohama,
Kanagawa 223–8522 Japan  Jeffrey G. Mangum National Radio Astronomy
Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475, USA  Susanne
Aalto Department of Space, Earth and Environment, Chalmers University of
Technology, SE-412 96 Göteborg, Sweden  Rebeca Aladro Max-Planck-Institut für
Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany  Mathilde Bouvier
Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The
Netherlands  Laura Colzi Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de
Ajalvir Km. 4, Torrejón de Ardoz, 28850 Madrid, Spain  Kimberly L. Emig Jansky
Fellow of the National Radio Astronomy Observatory  National Radio Astronomy
Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475, USA  Rubén
Herrero-Illana European Southern Observatory, Alonso de Córdova, 3107,
Vitacura, Santiago 763-0355, Chile  Institute of Space Sciences (ICE, CSIC),
Campus UAB, Carrer de Magrans, E-08193 Barcelona, Spain  Ko-Yun Huang Leiden
Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands
Kotaro Kohno Institute of Astronomy, Graduate School of Science, The
University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan  Sabine König
Department of Space, Earth and Environment, Chalmers University of Technology,
Onsala Space Observatory, SE-439 92 Onsala, Sweden  Kouichiro Nakanishi
National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  Astronomical Science Program, Graduate Institute for Advanced
Studies, SOKENDAI, 2-21-1 Osawa, Mitaka, Tokyo 181-1855, Japan  Yuri Nishimura
Department of Astronomy, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan  Institute of Astronomy, Graduate School of Science, The
University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan  ALMA
Project, National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka,
Tokyo 181-8588, Japan  Shuro Takano Department of Physics, General Studies,
College of Engineering, Nihon University, Tamura-machi, Koriyama, Fukushima
963-8642, Japan  Víctor M. Rivilla Centro de Astrobiología (CAB), CSIC-INTA,
Ctra. de Ajalvir Km. 4, Torrejón de Ardoz, 28850 Madrid, Spain  Serena Viti
Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The
Netherlands  Yoshimasa Watanabe Materials Science and Engineering, College of
Engineering, Shibaura Institute of Technology, 3-7-5 Toyosu, Koto-ku, Tokyo
135-8548, Japan  Paul P. van der Werf Leiden Observatory, Leiden University,
P.O. Box 9513, 2300 RA Leiden, The Netherlands  Yuki Yoshimura Institute of
Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa,
Mitaka, Tokyo 181-0015, Japan

###### Abstract

Molecular lines are powerful diagnostics of the physical and chemical
properties of the interstellar medium (ISM). These ISM properties, which
affect future star formation, are expected to differ in starburst galaxies
from those of more quiescent galaxies. We investigate the ISM properties in
the central molecular zone of the nearby starburst galaxy NGC 253 using the
ultra-wide millimeter spectral scan survey from the ALMA Large Program
ALCHEMI. We present an atlas of velocity-integrated images at a 1$\farcs$6
resolution of 148 unblended transitions from 44 species, including the first
extragalactic detection of HCNH+ and the first interferometric images of C3H+,
NO, HCS+. We conduct a principal component analysis (PCA) on these images to
extract correlated chemical species and to identify key groups of diagnostic
transitions. To the best of our knowledge, our dataset is currently the
largest astronomical set of molecular lines to which PCA has been applied. The
PCA can categorize transitions coming from different physical components in
NGC 253 such as i) young starburst tracers characterized by high-excitation
transitions of HC3N and complex organic molecules (COMs) versus tracers of on-
going star formation (radio recombination lines) and high-excitation
transitions of CCH and CN tracing PDRs, ii) tracers of cloud-collision-induced
shocks (low-excitation transitions of CH3OH, HNCO, HOCO+, and OCS) versus
shocks from star-formation-induced outflows (high-excitation transitions of
SiO), as well as iii) outflows showing emission from HOC+, CCH, H3O+, CO
isotopologues, HCN, HCO+, CS, and CN. Our findings show these intensities vary
with galactic dynamics, star formation activities, and stellar feedback.

Astrochemistry (75) — Starburst galaxies (1570) —

††facilities: ALMA††software: astropy (Astropy Collaboration et al., 2013,
2018), scikit-learn (Pedregosa et al., 2011), CASA (McMullin et al., 2007)

##  1 Introduction

The physical properties of molecular gas determine where and how future stars
will be forming. A powerful way to investigate the gas properties is via the
study of molecular line emission; the physical properties and processes taking
place during star formation affect rates of various chemical reactions and
thus leave an imprint on the chemical properties of molecular clouds, star-
forming regions, and protoplanetary systems.

While the Milky Way provides a good example of a galaxy with a relatively
quiescent star formation activity, it is of great interest to extend molecular
inventories to extragalactic sources and to probe the effects of extreme
starbursts and active galactic nuclei (AGNs), which cannot be studied in our
Galaxy due to the absence of such extreme activities (Meier & Turner, 2005;
Martín et al., 2006, 2011; Aladro et al., 2015; Nakajima et al., 2018; Takano
et al., 2019). The Atacama Large Millimeter/submillimeter Array (ALMA) has
enabled high angular resolution astrochemical observations targeting regions
with extreme activity in the AGN-host NGC 1068 (Takano et al., 2014) and NGC
1097 (Martín et al., 2015), the starburst galaxies NGC 253 (Meier et al.,
2015) and M83 (Harada et al., 2019), AGN-starburst composite NGC 4945 (Henkel
et al., 2018; Bellocchi et al., 2023), the ultra-/luminous infrared galaxies
(U/LIRGs) NGC 4418 (Costagliola et al., 2015), NGC 3256 (Harada et al., 2018),
and Arp 220 (Sakamoto et al., 2021).111See Table 1 of Martín et al. (2021) for
the references of extragalactic spectral line scan studies.

Previous observations found that the center of NGC 253 is particularly
chemically rich (e.g., Martín et al., 2006; Aladro et al., 2015; Meier et al.,
2015; Ando et al., 2017; Martín et al., 2021). It is a barred spiral galaxy at
a distance of 3.5 Mpc (Rekola et al., 2005). It hosts $\sim 2\,M_{\odot}$ yr-1
of star formation within the central molecular zone (CMZ) extending to several
hundred parsecs in size (Leroy et al., 2015; Bendo et al., 2015). The CMZ of
NGC 253 is rich in molecular gas (Sakamoto et al., 2011; Leroy et al., 2015,
2018). It also hosts outflows launched from the energetic starburst (Turner,
1985; Bolatto et al., 2013; Walter et al., 2017; Krieger et al., 2019; Levy et
al., 2021).

To investigate further the physical and chemical properties of the CMZ in NGC
253, and build upon its molecular richness, we conducted wide-frequency
coverage observations as an ALMA large program, ALCHEMI (ALMA Comprehensive
High-resolution Extragalactic Molecular Inventory; Martín et al., 2021,
hereafter MMH21). Thanks to ALCHEMI, a strong influence of cosmic rays has
been found by Holdship et al. (2021) from the CCH fractional abundance, by
Harada et al. (2021) from the HOC+ fractional abundance and HCO+/HOC+
abundance ratio, by Holdship et al. (2022) using the H3O+/SO ratio, and by
Behrens et al. (2022) to explain the low HCN/HNC ratio even in high-
temperature regions. Further findings from ALCHEMI show evidence for strong
PDRs in the star-forming regions (Harada et al., 2021), as well as the
presence of shocks at the orbital intersections by Humire et al. (2022) from
Class I methanol masers, by Harada et al. (2022) from HOCO+, and by Huang et
al. (2023) with HNCO and SiO. This survey also made the first extragalactic
detection of a phosphorus-bearing molecule, PN (Haasler et al., 2022). Tanaka
et al. (2023) make use of the multi-transition and multi-molecule data to
accurately map the physical conditions (gas densities, temperatures, and
molecular column densities) using hierarchical Bayesian analysis and found
high densities compared to those in the center of our own Galaxy.

In addition to the chemical richness, chemical differentiation (i.e.,
variation of chemistry) within galactic centers provides useful information.
Such chemical differentiation illuminates the changing physical conditions,
variations in dynamics, energetics, and evolution of the ISM. Correlations
between the chemistries and these conditions establish their relationships.
Such data of spatially-resolved chemistry tend to have very large datasets.
For their interpretation, statistical methods such as principal component
analysis (PCA) are useful in reducing the data size into a handful of
components that still retain most of the information. PCA has been used widely
in astronomical studies including the interpretation of molecular line
intensities from multiple species, although they are relatively limited in
terms of the number of transitions observed and included in their analysis
(e.g., Ungerechts et al., 1997; Meier & Turner, 2005; Costagliola et al.,
2011, 2015; Gratier et al., 2017; Saito et al., 2022; Gorski et al., 2023).

In this paper, we conduct a PCA on the ALCHEMI data to derive dominant
components within the plethora of lines detectable in the interstellar medium
(ISM) of NGC 253. We aim to find a set of molecular transitions that trace
physical conditions such as column densities, excitation conditions, star
formation activities, and shocks through this analysis. To the best of our
knowledge, this is currently the largest collection of molecular lines in a
given astrophysical object to which PCA has been applied. Not least, the
dataset is also characterized by a high degree of homogeneity, in terms of
observations and calibration procedures, angular and velocity resolutions, and
sensitivity, all of which make the ALCHEMI dataset a perfect statistical
sample to run PCA. This type of analysis allows for a comprehensive
understanding of the ALCHEMI dataset while previous ALCHEMI papers mostly
focused on a few particular species.

This paper is organized as follows. In Section 2, we describe the
observations, line identification, and data products. In Section 3, we report
features in the velocity-integrated images. The results of PCA are discussed
in Section 4. The implications of results are discussed in Section 5 and
Section 6 summarizes the main outcome of the work.

##  2 Observations and data products

The ALCHEMI survey (Project codes: 2017.1.00161.L, 2018.1.00162.S) includes
all the ALMA Bands 3–7 frequency coverage (84–375 GHz; $\lambda=3.6$–$0.8$ mm)
except for the small spectral regions affected by telluric absorption lines.
The observed area extends over $50\arcsec\times 20\arcsec$ ($830\times
330\mbox{ pc}^{2}$ on the sky), which covers most of the CMZ. The data
products are convolved into a common 1$\farcs$6 beam (=27 pc), and the maximum
recoverable scale is equal to or greater than $15\arcsec$. The rms noise
levels for a 10 km s-1 channel are $\sim 20$ mK, but exact values for each
spectral window can be found in MMH21 where a complete description of this
survey, as well as the details on calibration and imaging, are provided.

###  2.1 Line Identification

In order to identify spectral transitions in the dataset, we extracted the
spectra from a single pixel at position $\alpha_{J2000}=0^{\rm h}47^{\rm
m}33.33^{\rm s}$, $\delta_{J2000}=-25^{\circ}17^{\prime}15\farcs 73$. This
corresponds to GMC 6 in Leroy et al. (2015) (see Figure 1 for the location)
and is the position with the brightest molecular emission of most detected
species (Haasler et al., 2022). This position has the most complex emission
spectra and therefore is used as the reference position for the purposes of
line blending evaluation. Line identification was done by fitting simulated
spectra under the assumption of local thermodynamic equilibrium (LTE) from all
the commonly observed species, based on previous extragalactic and Galactic
surveys, using MADCUBA222Madrid Data Cube Analysis on ImageJ is a software
program developed at the Center of Astrobiology (CAB) in Madrid;
https://cab.inta-csic.es/madcuba/  (Martín et al., 2019). The model on which
this line identification is based includes a total of 330 entries, including
the most common rare isotopologues (13C, 18O, 17O, 15N, 34S, 33S, 36S, 29Si,
30Si, 37Cl), vibrationally excited modes for relevant species, as well as
hydrogen and helium recombination lines. In total, 134 species are detected. A
line identification using the ACA resolution was already presented by MMH21,
but a full presentation of the line identification including the 12-meter
array and new detections in the extragalactic ISM will be provided in
dedicated papers.

###  2.2 Velocity-integrated images

From the list of all transitions included in the modeled spectra used for the
line identification, we produced velocity-integrated images of isolated lines.
A transition is considered isolated if it meets one of the following two
criteria. 1) A line has less than 10 % contamination compared with its
intensity in fitted models assuming LTE. 2) A line has its line center
separated by more than 210 km s-1, three times the maximum line width, from
other detected lines close in frequency. This level of contamination is
negligible considering the $\sim 15\%$ calibration uncertainty (MMH21). We
also considered a transition as isolated if blended by the same molecule
(e.g., CCH, CH3CCH, CN), but used masks to include all the transitions as we
describe in the next paragraph (see Appendix A for the list of transitions).
Our criterion for detection was having $>3\sigma_{\rm ch}$ intensity at their
peak, where $\sigma_{\rm ch}$ is the rms value of a single channel. Our line
identification was done at the most molecule-rich position, and the
contamination is likely less than 10% in other positions. We also eliminated
lines with heavy contamination from unidentified lines through the visual
inspection of spectra. Spectroscopic information that is used in this paper is
listed in Appendix B. We note that these criteria for isolated lines may be
too strict in some cases, and some transitions not included in this study may
contain only very little contamination. It is the reason this study omits some
transitions used in prior ALCHEMI papers, in which the individual transitions
were more thoroughly checked for contamination.

To create velocity-integrated (moment 0) images, we applied a 3D mask created
from a reference transition with threshold cutoffs to a cube of each isolated
transition. It is to exclude contamination from other transitions close in
frequency in a similar way as in Harada et al. (2021). We choose a reference
transition with strong emission and determine the threshold cutoff based on
the reference transition, then apply this cutoff to cubes of other
transitions. We use the 12CO($J=$1–0) image cube as a reference image to
create masks for transitions that are not 12CO but have relatively strong
emission (HCN, HCO+, HNC, CS, CCH, and CN). For other species, we used
HCN($J=$1–0) as a reference transition. The 12CO transitions have the highest
S/N ratio in our dataset, and emission is detected in larger volumes of their
cubes than for any other line. Therefore, a 12CO mask is suitable for
collecting weak extended emission in transitions such as HCN. On the other
hand, the HCN(1–0) mask is more suitable for weak lines to eliminate
contamination from neighbors, because there is one location (GMC6) where the
CO line is very broad.

To create the binary 0/1 masks, we first convolved the reference cubes to a
lower angular resolution of 3$\farcs$2 to capture faint extended emission in
the masks. We then set the mask value to 1 only at locations in the cubes with
a $>15\sigma_{\rm ch}$ detection for the CO(1–0) mask or $>5\sigma_{\rm ch}$
detection for the HCN(1–0) mask, where $\sigma_{\rm ch}$ is the single-channel
rms value in the 3$\farcs$2 resolution reference cube. A 1$\farcs$6 resolution
integrated intensity map is created as the input line cube multiplied by the
mask cube on the same grid and integrated over velocity. The above cutoffs for
the masks may seem unnecessarily high, but these masks should still include
all regions with significant emission, because all the transitions except
those of 12CO are more than 10 times weaker than CO(1–0) in extended regions,
and HCN(1–0) has more than 3 times higher S/N compared to the transitions to
which the HCN mask was applied. We confirmed that our masking did not exclude
target emission by visually inspecting the original and masked cubes. This
high cutoff for CO(1–0) also helps to exclude nitrogen sulfide (NS)
transitions at the rest frequency $f_{\rm rest}=115.154$ GHz near CO(1–0),
whose rest frequency is at 115.271 GHz. We applied our masked integration to
images within 400 km s-1 from the systemic velocity of 258 km s-1. The
effective range of integration is narrower than 800 km s-1 because of the
masking.

For 12CO transitions, we created the masks from themselves (e.g., using CO 3–2
as the reference cube for CO 3–2). We again convolved the reference cubes to
3$\farcs$2, twice the original beam size, and kept the positions with emission
brighter than $2\sigma_{\rm ch}$. These masks were applied to the original
1$\farcs$6 resolution cubes.

For transitions blended with the same species (e.g., CCH), we first made cubes
with their velocity coordinates defined from one of the transitions (a
velocity reference line), with enough width to cover all the detectable
transitions of interest. Then, we created masks to include emissions from
multiple transitions from the same species, using the known velocity shifts
from the velocity reference line. The list of such transitions can be found in
Appendix A.

In addition to the 3D masking applied above, we also applied a threshold
cutoff of $3\sigma$ to the integrated intensity images. The rms $\sigma$ of
the integrated intensities is estimated as

$\sigma=\sigma_{\rm ch}\Delta v\sqrt{N},$ (1)

where $\sigma_{\rm ch}$ is the rms for one channel in the cube, $\Delta v$ is
the channel width (10 km s-1), and $N$ is the number of channels used for
integration.

After removing transitions with peak intensities less than 10 K km s-1 and
Band 5 transitions with channels affected by severe artifacts (SiO 4–3, N2H+
2–1, see also Section 3), we use 148 transitions of molecules or radio
recombination lines (RRLs) from 44 species and two continuum images (one
representative of synchrotron/free-free emission at the lower frequency range
and the other representative of dust thermal emission at the higher frequency
range) for PCA. The velocity-integrated images, original cubes, and masked
cubes will be made public as a part of the high level data products to be made
available through the ALMA science archive. In addition, a matrix of
standardized intensities is also publicly
available333http://github.com/nanaseharada/alchemi_pca.

##  3 Distribution of transition intensities

We first describe locations in the CMZ of NGC 253 using the two color-
composite (red-green-blue; RGB) images of selected transitions in Figure 1.
They show the contrast between images of CO(1–0), H39$\alpha$,
CH3OH($2_{K}-1_{K}$), HC3N(25–24), and CN(3–2). This choice of lines for the
RGB channels is justified by the results of the PCA, from which we select
these five representative lines as highlighting the extreme variance in the
dataset, hence giving the RGB images high contrast. As we discuss later in
Section 4.3, these transitions represent the first three principal components
(PC1, positive and negative PC2, positive and negative PC3). In this figure,
we show the positions of giant molecular clouds (GMCs) that we refer to
throughout the paper. They are close to the GMCs discussed in Leroy et al.
(2015), but slightly modified to better match the molecular emission peaks in
most species. Because the GMC peaks are hard to pinpoint from bright,
ubiquitous, and extended transitions such as CO and HCN, we identified peaks
from species that trace more specific locations such as SiO (GMCs 1, 7, 8, and
9), HC3N(GMCs 3, 4, 6), and H39$\alpha$ (GMC5). These positions are listed in
Table 1.

Some of these clouds reside at orbital intersections of bar orbits (the
$x_{1}$ orbits) and central orbits (the $x_{2}$ orbits) resulting from the
stellar bar potential in NGC 253 (Sorai et al., 2000; Das et al., 2001). The
gas on the bar orbits flows in from the northeast and southwest directions.
GMCs 1, 2, 7, 8, and 9 have signs of shocks detected from Class I methanol
masers (Humire et al., 2022). Other clouds are locations of young massive
star-cluster formation. An RRL, H40$\alpha$, has been detected in GMCs 3–6
(Bendo et al., 2015; Mills et al., 2021). We note that the RRL in GMC3 is not
obvious in our 1$\farcs$6 beam. Only high-angular resolution observations
($\sim$0$\farcs$2) by Mills et al. (2021) found weaker and narrower
H40$\alpha$ in GMC 3 compared with GMCs 4–6.

Table 1: GMC positions and their properties

GMC | RA(ICRS) | DEC(ICRS) | Remarks
---|---|---|---
| 0h47m |  $-25^{\circ}$17′ |
1 | 31$\fs$93 | 29$\farcs$0 | Class I methanol maser$\dagger$
2 | 32$\fs$36 | 18$\farcs$8 | Class I methanol maser$\dagger$
3 | 32$\fs$81 | 21$\farcs$2 | Clumps 1-3$\ddagger$
4 | 32$\fs$95 | 19$\farcs$8 | Clumps 4-7$\ddagger$
5 | 33$\fs$16 | 17$\farcs$3 | Clumps 8-13$\ddagger$
6 | 33$\fs$33 | 15$\farcs$7 | Clump 14$\ddagger$
|  |  | Line ID position∗
7 | 33$\fs$65 | 13$\farcs$1 | Class I methanol maser$\dagger$
8 | 33$\fs$94 | 10$\farcs$9 | Class I methanol maser$\dagger$
9 | 34$\fs$14 | 12$\farcs$0 | Class I methanol maser$\dagger$

Note. — $\dagger$Humire et al. (2022) $\ddagger$Leroy et al. (2018) ∗ Position
of the brightest molecular emission

Figure 1: (Top) Velocity-integrated images of CO(1–0) in blue, H39$\alpha$ in
red, and CH3OH($2_{K}-1_{K}$) around the rest frequency of 96.7 GHz in green.
Rough positions of parts of the $x_{1}$ orbits and the full $x_{2}$ orbits
further inside are shown with dash-dotted and dashed lines, respectively. Note
that these are just some examples of the presumably large families of possible
$x_{1}$ and $x_{2}$ orbits. There likely exists another $x_{2}$ orbit
connecting GMCs 3-6 or 3-7 (see Levy et al., 2022) almost fully edge-on. The
synthetic beam size is shown at the left bottom corner as a white circle.
(Bottom) The same as the top figure, but with HC3N(25–24) in red, and CN(3–2)
in green. Intensities are scaled so that faint lines have similar dynamic
ranges as other lines. We employed the algorithm by Lupton et al. (2004) to
produce RGB images.

Figures 2–22 show the integrated intensity images created using the procedure
described in the previous section. Continuum images at 95 GHz (3.2 mm) and
361.5 GHz (0.83 mm) are also included in the last two panels of Figure 22. As
already illustrated by Meier et al. (2015), morphological differences are
obvious from simple visual inspection. While distributions of some transitions
are more extended, some are more concentrated around the central parts, and
others are enhanced in the outskirts of the CMZ (GMCs 1, 2, 8, and 9) in NGC
253.

In general, CO isotopologues (12C16O, 13C16O, 12C18O, and 12C17O, hereafter
12CO, 13CO, C18O, and C17O) exhibit more extended emission than other species.
Higher-$J$ transitions are slightly more compact than $J=1-0$ transitions,
likely because the gas is more excited near active star-forming regions. Rarer
isotopologues of CO show more compact emission than the main isotopologue
(12C16O) both due to the lower optical depths and lower S/N.

Transitions of species with high critical densities ($n_{\rm
crit}>10^{4}\,$cm-3 in optically thin cases) such as HCN, HCO+, CS, N2H+, and
CCH also show significantly extended emission in their $J=1-0$ transition, but
higher-$J$ transitions are more compact, some of which have already been shown
in our previous literature (Harada et al., 2021; Holdship et al., 2021;
Behrens et al., 2022). We see this trend in CO isotopologues as well, but it
is more prominent in these transitions with high critical densities. Similar
to CO isotopologues, emission from rarer isotopologues of other molecules is
more compact than from their main counterparts.444 Sakamoto et al. (2021)
reported similar size variations among various sub/millimeter emissions from
their ALMA imaging spectral scans toward the nuclei of luminous infrared
galaxies NGC 4418 and Arp 220. Taken together, these observations caution
against the assumption of a common beam-filling factor in the excitation
analysis of molecular lines for unresolved galactic nuclei.

We find somewhat compact distributions in the 3-mm continuum emission, most of
which should originate from free-free emission, and RRLs compared with other
transitions mentioned above. They are expected to arise from the ionized gas,
tracing the ongoing star formation (e.g., Bendo et al., 2015). These types of
emission and transitions such as HCN, HCO+, CCH, and CN with higher upper-
state energies emit strongly in similar locations (see also Holdship et al.,
2021, for CCH).

On the other hand, CH3OH, HNCO, and HOCO+ transitions with low upper-state
energies are enhanced in the outskirts of the NGC 253 CMZ (Figures 11-15). Out
of these species, CH3OH and HNCO are known as tracers of slow shocks555Fast
shocks could enhance methanol as well according to some modeling results,
although there is also an observational result that suggests methanol
destruction in fast shocks (Suutarinen et al., 2014)., and their enhancements
are likely attributed to widespread shocks in those regions (Meier et al.,
2015; Huang et al., 2023). Harada et al. (2022) argued that HOCO+ should be
tracing shock-evaporated CO2, and its relation with shocks is suggested by the
similarity with these species. A tracer of fast shocks, SiO, also shows some
emission at the outskirts of the CMZ. Unlike weaker shock tracers, SiO
transitions emit strongly near the starburst regions as well (Figures 10-11).

To the best of our knowledge, we also present the first extragalactic
detection and/or first extragalactic interferometric images for certain
species. For example, HCNH+ (Figure 21), protonated HCN or HNC (see Figures 3
and 5), is detected for the first time outside the Milky Way. It shows a
similar distribution as that of HCN (Figure 3). One production pathway of
HCNH+ is via protonation of HCN or HNC by proton exchange reactions with
H${}_{3}^{+}$ or HCO+. Another formation route is through HNC+ or HCN+
reacting with molecular hydrogen. A survey of high-mass star-forming regions
suggested that HCNH+ is more abundant in cold starless cores compared with
more evolved sources (Fontani et al., 2021). HCNH+(2–1) in the CMZ of NGC 253
has a similar distribution as that of HCN(1–0), but with a lower S/N (compare
Figures 3 and 21). It is distributed widely in the entire CMZ rather than
being concentrated around the starburst region. This is consistent with the
picture that HCNH+ is more abundant in cold clouds. We also obtained the first
extragalactic images of C3H+, NO, and HCS+. C3H+ (Figure 20) is a tracer of
PDRs, and was first detected by Pety et al. (2012) in a PDR of the Horsehead
nebula. With a lower fractional abundance, this species was also detected in
TMC-1 (Cernicharo et al., 2022). The first extragalactic detection of this
species was made towards a molecular absorber at $z=0.89$ (Tercero et al.,
2020). Interstellar NO has been detected already several decades ago (Liszt &
Turner, 1978). Its first extragalactic detection was made in NGC 253 by Martín
et al. (2003). Protonated CS, HCS+ was detected for the first time by Thaddeus
et al. (1981) rather ubiquitously in hot cores and cold clouds. Its first
extragalactic detection was also in a molecular absorber at $z=0.89$ (Muller
et al., 2013).

There are some moment 0 images that are used in PCA, but not shown in this
paper because they will be included in dedicated publications, as indicated in
Table 3. This is the case of the maps of sulfur-bearing species since M.
Bouvier et al. (in preparation) will provide a comprehensive presentation to
investigate the origin of the emission of various sulfur-bearing species and
discuss physical processes dominate their release/presence onto the gas phase.
Other transitions of CO, HCN, and HCO+ contain rarer isotopologues such as
13C, 18O, 17O, and 15N. J. Butterworth et al. (in preparation) will discuss
isotopic ratios through the lens of these isotopologues in order to
investigate the possible relation between these ratios and the ages of SSCs
observed in NGC 253. These images that are currently not included in this
version of paper will be included later in the arXiv preprints after above
papers are accepted.

We note that some of the transitions in Band 5 have poor image qualities
influenced by artifacts (e.g., H13CO+ 2–1, HOC+ 2–1). These transitions come
from the datasets that were not fully observed and did not meet the
homogeneous sensitivity requested in this survey (see Section 3.2 of MMH21, ).
This band has a water absorption line near the middle of the band, and
transitions nearby are prone to increased noise levels. PCA results for these
transitions are less certain.

Figure 2:  Velocity-integrated images integrated from within 400 km s-1 of the
systemic velocity with masking as described in Sect. 2.1. Contour levels are
set at 10, 30, 50, 70, and 90 % of their peak values. Maximum values of color
scales are set for 120 % of their peak values. We show more commonly observed
species first, and transitions are in ascending order of frequency within the
same species. The transition names are in the upper left corner, and the line
rest frequency and the energy level in the upper state are also in each panel.
The 1$\farcs$6 beam is in the bottom left corner. The plus signs indicate GMC
locations. The images are corrected for the primary beam response patterns of
the ALMA antennas.  Figure 3: Same as Figure 2.

Figure 4: Same as Figure 2.

Figure 5: Same as Figure 2.

Figure 6: Same as Figure 2.

Figure 7: Same as Figure 2.

Figure 8: Same as Figure 2.

Figure 9: Same as Figure 2.

Figure 10: Same as Figure 2.

Figure 11: Same as Figure 2.

Figure 12: Same as Figure 2.

Figure 13: Same as Figure 2.

Figure 14: Same as Figure 2.

Figure 15: Same as Figure 2.

Figure 16: Same as Figure 2.

Figure 17: Same as Figure 2.

Figure 18: Same as Figure 2.

Figure 19: Same as Figure 2.

Figure 20: Same as Figure 2.

Figure 21: Same as Figure 2.

Figure 22: Same as Figure 2.

##  4 Principal component analysis

Astronomical datasets often have a large number of variables that are
correlated to some degree. It is a challenge to extract information from such
data and there is a need for somewhat automatized statistical investigation.
Dimensionality reduction is especially a useful technique in this regard,
deriving a small number of components to reproduce most of the variations
within the data. Multiple methods of dimensionality reduction have been
developed so far, including PCA, non-negative matrix factorization (NMF), and
independent component analysis (ICA). PCA is one of the oldest and most
commonly used with more than a century-old history (Pearson, 1901). As
mentioned in Section 1, PCA has been most widely used for molecular line
studies (see references in Section 1). Moreover, this work aims to push the
PCA analysis to the next step by applying it to a much larger set of molecular
lines.

###  4.1 The basics of PCA

We begin by briefly describing the basics of PCA relevant to this work,
leaving details and proofs to textbooks (e.g., Jolliffe, 2002). A data set of
$p$ variables measured at $n$ targets (e.g., intensities of $p$ emissions at
each of the $n$ sky positions) can be viewed as $n$ data points in a
$p$-dimensional data space. PCA sets up an orthogonal coordinate system in the
data space with the origin at the centroid of the data points and the first
coordinate axis (PC1 axis) in the direction of the largest variance of the
data point positions along the axis. The PC2, PC3, …, axes are successively
defined as orthogonal to all previous axes and having the largest possible
variance of the data positions along the axis. Coordinate values on a PC axis
are called PC scores. The PC scores of the data points have a mean of zero for
each PC and smaller variances for later PCs. Accordingly, the relative
positions of the data points can be described reasonably well by using only
their first $q\;(\leq p)$ PC scores, since all the data points have
approximately the same scores ($\approx 0$) for the rest of the PCs. This
description reduces the dimensionality of the data set from $p$ to $q$.

Individual PC axes are found in the following way. Let $x_{ij}$ be the value
of the $j$-th variable measured at the $i$-th target. The data matrix
$\bm{X}=(x_{ij})_{1\leq i\leq n,1\leq j\leq p}$ has $x_{ij}$ in the $i$-th
row, $j$-th column. We use the data matrix that is standardized for each
variable; i.e., each column has its mean subtracted and then divided by its
standard deviation. This subtraction aligns the coordinate origin with the
centroid of the data points. The correlation matrix for the $p$ variables is
then $\bm{\Sigma}=\frac{1}{n}\bm{X}^{T}\bm{X}$, where T denotes the transpose.
The $p\times p$ real matrix is symmetric and non-negative definite, hence
having $p$ eigenvalues that are real and non-negative. We denote them as
$\lambda_{1}\geq\cdots\geq\lambda_{p}\geq 0$. The direction vector of the
PC$j$ axis is the unit eigenvector of $\bm{\Sigma}$ associated with the $j$-th
largest eigenvalue $\lambda_{j}$, and we denote it as
$\bm{a}_{j}=(a_{1j},\cdots,a_{pj})^{T}$. They are orthonormal, i.e.,
$\bm{a}_{i}\cdot\bm{a}_{j}=\delta_{ij}$. The variance of the PC$j$ scores is
$\lambda_{j}$.

The transformation between the data coordinates and the PC coordinates is as
follows. Let $z_{ij}$ be the PC$j$ score of the $i$-th target. The target is
located at $\bm{z}_{i}=(z_{i1},\cdots,z_{iq})$ in the $q$-dimensional PC
coordinates and $\bm{x}_{i}=(x_{i1},\cdots,x_{ip})$ in the standardized data
coordinates. Since $z_{ij}$ is the projection of $\bm{x}_{i}$ onto the PC$j$
axis along the unit vector $\bm{a}_{j}$, we have $z_{ij}=\bm{x}_{i}\bm{a}_{j}$
and $\bm{z}_{i}=\bm{x}_{i}\bm{A}$, where $\bm{A}$ is the matrix
$\bm{A}=(\bm{a}_{1},\cdots,\bm{a}_{q})=(a_{ij})_{1\leq i\leq p,1\leq j\leq
q}$. The element $a_{ij}$ is the direction cosine between the $i$-th data axis
and the PC$j$ axis, and is usually referred to as a ‘coefficient’ or ‘loading’
for PC$j$. If it is positive, a positive value of the $i$-th variable makes a
positive contribution to the PC$j$ score, while a negative $a_{ij}$ indicates
a negative contribution.

For the simplest case of two variables with a strong correlation, the PC1 axis
would lie in the direction along the linear fit of the data points, while the
PC2 axis would be in the direction perpendicular to it. Therefore, PC1 scores
indicate the strength of the common properties in the two variables, while PC2
scores reflect the strength of the differential properties of the two
variables.

PC scores are uncorrelated between different PCs because $\bm{\Sigma}$ is
diagonal in the PC coordinates. Therefore, if variables are determined by
underlying parameters with some correlations, then there would not be a one-
to-one relationship between those parameters and the principal components. We
discuss this caveat when we apply PCA to the ALCHEMI data.

###  4.2 Application of PCA to the ALCHEMI data

We binned the velocity-integrated images into hexagonal pixels with a
horizontal length of 0$\farcs$8 using Python’s hexbin function. The number of
hex-pixels is $n=2723$. The choice of half the beam size is made for the
Nyquist-like sampling. We filled the masked out channels/pixels with zero.

We standardized our data for each variable (i.e., transition or emission) as
described in Section 4.1. The division by the standard deviation is a common
procedure for a data set whose variables differ greatly in their variances. It
also makes our PCA independent of the choice of data units between Jy and K.

The variables in our analysis are transitions and the targets are sky
positions, as in Ungerechts et al. (1997) or Meier & Turner (2005). With 148
lines and 2 continuum images, we have $p=150$ variables. Instead, we could
have used different transitions as targets and their intensities at the pixels
as variables, as in facial recognition, where eigenvectors are eigenfaces.
However, we did not use the approach to avoid the problem of high dimension
and low sample size (i.e., $p>n$)666See Takeuchi et al. (2022) for how to
properly perform PCA on such a dataset.. It is also possible to conduct a PCA
on image cubes, instead of integrated intensities. Nonetheless, we use
integrated intensity images first, because integrated intensities tend to have
higher signal-to-noise ratios. A PCA on image cubes, possibly with a lower
number of transitions, is planned as a future work. Consequently, there are
PC$k$ scores for all hexagonal pixels, $\\{z_{ik}\\}_{i=1,\cdots,n}$. We used
the PCA package from the scikit-learn project (Pedregosa et al., 2011).

The number of principal components $q$ must be chosen to retain most of the
information in the data set. The information in the context of PCA is the
configuration of the data points and is evaluated with the variances. For
example, PC$k$ has the variance of $\lambda_{k}$ in its scores and therefore
has less information for larger $k$. Figure 23, called a scree plot, shows the
fraction of the variance of each PC in the total variance,
$f(k)\equiv\lambda_{k}/\sum_{j=1}^{p}\lambda_{j}$. In this plot, the PCs that
exceed the least-squares linear fit to the high-order PCs (PCs 6–20) are
considered more significant. The first three PCs have relatively large
fractions of information (i.e., variance) $f(k)=$ 79.4, 9.0, and 5.0 %,
respectively, for $k=$1, 2, 3, with a cumulative contribution of 93% of the
total. PCs 4 and 5 also show some significance ($f(k)=$ 1.8 and 1.2 %)
exceeding the fit at PCs 6–20, each of which has $f(k)\leq$ 0.6%. As shown in
Section 4.3, the PC score maps for PC6 and higher are noisier than for PCs
1–5. Therefore, we only discuss PCs 1–5 in this paper; they have
$\sum_{j=1}^{5}\lambda_{j}/\sum_{j=1}^{p}\lambda_{j}=96.4$% of the total
variance.

PCA is one of the unsupervised machine learning methods, which means that any
known properties of data, or labels, are not used. This means that PCA results
are derived independent of what kind of species the molecules are, or energy
levels of the transitions. We will see below (in Section 4.4) that our PCA
could recover such information to some extent from the data alone.

Figure 23: (Top) A scree plot for our PCA. (Bottom) The same as the top
figure, but for a narrower range of the y-axis.

###  4.3 PC score maps

Figure 24: PC scores of the 1st to 6th principal components (i.e.,
$\\{z_{ij}\\}_{1\leq i\leq n}$ for $j=1-6$). GMC positions listed in Table 1
are shown with black crosses. The origin is at the phase center of the
observation, $\alpha=00^{h}47^{m}33.26^{s}$,
$\delta=-25^{\circ}17^{\prime}17.7^{\prime\prime}$ (ICRS).

Figure 25: The same as 24, but for PCs7-10.

Figures 24 and 25 show the projection of data points onto the PC axes
($\\{z_{ik}\\}_{i=1,\cdots,n}$ for PC$k$ in the notation of Section 4.1), or
the maps of PC scores for the first ten PCs. Here we provide a brief overview
of the relationship between the individual PCs and the emission of certain
molecular species. Note that the PC scores are dimensionless and have little
meaning in their absolute values. Only their relative magnitudes between
different PCs and positions are relevant. The PC1 score map (Figure 24 top
left) shows a rough shape of the CMZ, similar to the distribution of CO
isotopologues. The positive PC2 scores (Figure 24 top right) appear near the
central starburst regions (GMC4-6), where the high-$J$ transitions of HCN
emission and RRLs have high intensities. On the other hand, the negative PC2
scores are at the outskirts of the CMZ, where low-$J$ transitions of CH3OH and
HNCO are enhanced (GMC1-3,7-9). Note that the green (CH3OH $2_{K}$–$1_{K}$)
and red (H39$\alpha$) colors in Figure 1 (top) show remarkable similarity with
blue (negative) and red (positive) colors in the PC2 score map (Figure 24).
High absolute values of the PC3 scores (Figure 24 middle left) appear around
positions with positive values of the PC2 score map. Positive PC3 scores
appear at the peak location of HC3N vibrationally excited lines (GMC6) while
the negative PC3 scores are at the peak of the RRL (GMC5). We again note that
the green (CN 3–2) and red (HC3N 25–24) colors in Figure 1 (bottom) resemble
the blue and red colors in the PC3 map.

As mentioned earlier, contributions of PCs 4 and 5 (Figure 24 middle right and
bottom left) are less than those of PCs 1-3. They may contain multiple
physical components because PCs (eigenvectors) are by definition orthogonal to
each other while each physical component does not necessarily influence
intensities orthogonally in dimensions of transition intensities. For example,
PC4 has positive PC scores in most of the GMCs except for GMCs 3 and 7. The
PC4 score map appears similar to that of PC2, except for GMCs in the outskirts
(GMCs 1, 2, 7, 8, and 9). The PC4 score map also has an extended structure out
of the galactic plane with the negative component. PC5 shows negative PC
scores in most of the GMCs, but it has positive PC scores in GMC6 and diffuse
extended structure.

In principle, PC6 and higher-order maps still provide some information about
this galaxy, not merely contributions due to noise. However, PC score maps
from PC6 to PC10 cannot be uniquely linked to the emission of individual
species and therefore cannot easily be linked to specific physical processes
(Figures 24 and 25).

###  4.4 Relationship with physical parameters

Here we examine whether each PC has any clear relationship with physical
parameters such as column densities, kinetic temperatures, and volume
densities. Figure 26 shows the relationship between the PC scores in our
results and physical conditions derived by Tanaka et al. (2023) from the
ALCHEMI data. These physical parameters were obtained from the hierarchical
Bayesian analysis method developed by Tanaka et al. (2018). Tanaka et al.
(2023) present results from different sets of transitions. Among their
results, we use their “high-HB” model derived from species with higher
critical densities such as HCN, HCO+, and their isotopologues, HC3N, and SiO,
instead of results derived from CO isotopologues. This is because our PCA used
more species with high critical densities than CO isotopologues. Pixels from
$N_{\rm H2}$-weighted averaged images along the velocity axis are binned to
the same hexagonal pixels as our PCA. We omitted hexagonal pixels that do not
have physical parameters derived for all the original pixels within the
hexagon. We show pixels within one standard deviation of all the data points
of PC1 scores from zero in grey because they represent pixels with low signal-
to-noise ratios. The standard deviation of PC1 scores is 21. We note, however,
that not all the grey points are noise. We only show relations with PCs 1 and
2 because there are no clear trends for PC3 and higher.

PC1 scores show a clear increasing trend with increasing column densities.
They also have a similar, but weaker correlation with volume densities, and
even to a lesser degree, with temperatures. These results suggest that PC1
primarily represents the overall molecular gas content. The volume densities
tend to be higher in high column-density regions in Tanaka et al. (2023), and
this is likely why PC1 and the volume densities have a positive correlation.
The high-density regions are also regions with high star formation rates,
which explains the trend of increasing PC1 scores with increasing
temperatures.

PC2 scores increase with the increasing values of all of these physical
parameters (column densities, volume densities, and temperature) if we focus
on high S/N pixels (blue points in Figure 26). This correlation between PC2
scores and physical parameters is tightest in the temperature, while it gets
slightly weaker in the density. The scatter becomes larger between PC2 scores
and the column density. To consider the relationship between PC2 and physical
parameters in a different perspective, we consider PC2 coefficients ($a_{i2}$
in the notation of Section 4.1), instead of PC scores discussed above, of
multiple transitions of HCN, CS, C34S, and HC3N. Figure 27 shows the
dependence of PC2 coefficients on upper-state energies (left figure) and
critical densities (right figure). For all the species, PC2 coefficients tend
to be higher with higher upper-state energies and critical densities.
Transitions with high PC2 coefficients clearly emit preferentially from
increasingly higher-excitation regions.

Figure 26: Physical parameters (top: column densities, middle: volume
densities, bottom: temperatures) and PC scores in individual hexagonal pixels.
Physical parameters are taken from “high-HB” results from Tanaka et al.
(2023), which are the ones derived from species with relatively high critical
densities. Grey points represent data points that lie within one standard
deviation of PC1 scores from zero, pixels with low signal-to-noise ratios.

Figure 27: Comparisons between PC2 coefficients and (left) upper state
energies and (right) critical densities from the bright linear rotors, where
the reasonable connection between $J$, $E_{\rm up}$, and $n_{\rm cr}$ is
expected. For HCN and CS, values of critical densities are taken from Shirley
(2015) considering multiple levels. We assume the same critical densities as
CS for C34S. We take values of critical densities from Wernli et al. (2007)
for HC3N.

###  4.5 PC coefficients of each species

Figures 28 – 32 show the coefficients of various transitions for PCs 1–5,
i.e., $a_{ij}$ for a transition $i$ and PC$j$ for $j=$ 1–5. A table of these
coefficients can be found in Appendix C. Figure 33 also shows the same
information, but for transitions that are not shown in Figures 28 \- 32.
Errors of the coefficients are evaluated in Appendix D and are negligible.
Figures 28 \- 32 are separated into four panels. The first two panels are for
dense gas tracers (HCN, HNC, HCO+, CN, and their isotopologues, as well as
N2H+, HC3N) and shock tracers (CH3OH, SiO, HNCO and HOCO+)777While HOCO+ is
not a conventional shock tracer, we consider it a shock tracer here because of
the results obtained by Harada et al. (2022) who found an association between
HOCO+ and shocks in NGC 253.. Colors indicate the upper state energies of
transitions from 0 to 100 K in cyan to blue for dense gas tracers, and lime to
green for shock gas tracers. The third panel shows CO isotopologues,
continuuum, RRLs, vibrationally-excited HC3N, and CH2NH, and the fourth panel
is for sulfur-bearing molecules. The fourth panel also shows upper-state
energies with colors (0 to 100 K with pink to purple).

Figures 28 and 33 show that all the transitions have positive coefficients in
PC1. The variation of PC1 coefficients among transitions is not large. This
result can easily be explained by the fact that PC1 correlates with the
emission brightness of these transitions. Consequently, PC1 in general traces
overall gas column density as described in Section 4.4 because regions with
higher column densities tend to have higher brightness. Star formation tracers
such as RRLs and the 3-mm continuum also have positive PC1 values because the
star formation rate is generally higher towards the high column-density
regions.

High absolute values of PC3 scores (around GMCs 5 and 6) appear in regions
with positive values of PC2 scores (Figure 24). Negative PC3 coefficients are
seen for RRLs, while high-excitation dense gas tracers have positive PC3
coefficients (Figure 30). It is interesting to note that most of the shock
tracers and sulfur-bearing molecules have positive PC3 coefficients, unlike
dense gas tracers. Meanwhile, CO isotopologues and continua all have negative
PC3 coefficients. PC3 separates high-excitation transitions with high PC2
coefficients into different groups as illustrated in Figure 34. High
excitation transitions of CN, HNC, HCN, and RRLs have negative PC3
coefficients (the fourth quadrant in Figure 34)888We number quadrants counter-
clockwise from the top right. while HC3N transitions have positive ones (the
first quadrant in Figure 34).

As mentioned previously, there are some similarities between the PC2 and PC4
score maps near the central starbursts. However, there are also interesting
exceptions in the outskirts (Figure 24). Therefore, if a transition has
positive coefficients for PC2 and PC4 (the first quadrant in Figure 35), it
has stronger emission in GMCs 4, 5, and 6 than that can be accounted for by
PC2 alone. On the other hand, if a transition has negative PC2 and positive
PC4 coefficients (the second quadrant in Figure 35), contributions from the
CMZ outskirts are significant (GMCs 1, 2, 7, 8, 9).

Positive PC5 scores are either at GMC6 or the component extended towards
higher galactic latitudes likely associated with outflows, while the area of
negative PC4 scores contain these outflow-related extended structures and GMC3
(Figure 24). Because the vertically extended structure is likely caused by
outflows and is of interest, we picked out transitions tracing this structure
by comparing PC4 and PC5 coefficients in Figure 36. The transitions of our
interest should have negative PC4 and positive PC5 coefficients and in the
second quadrant of the Figure. Because these outflow features also appear
weakly as negatives in the PC3 score map, we also compare PC3 and PC5
coefficients in Figure 37, where the outflow tracers should be in the second
quadrant. It turns out that most transitions that appear in the second
quadrant in Figure 36 also appear in the second quadrant of Figure 37.
Transitions having this outflow feature are CO isotopologues, CCH(1–0),
HOC+(1–0), and H3O+. Although they are not shown in the figure, dense gas
tracers seen in an outflow (Southwest streamer) by Walter et al. (2017) such
as HCN, HCO+, CS in J=1–0, and CN(N=1–0) also have negative PC4 and positive
PC5 coefficients. We note that CCH(1–0) has been detected also in the outflow
in Maffei 2 (Meier & Turner, 2012) and from an AGN in NGC 1068 (García-Burillo
et al., 2017; Saito et al., 2022).

It is important to note that line intensities depend on column densities (or
abundances), temperatures, and densities. Consequently, our results show very
different PC coefficients, which means different distributions, even for the
same species depending on transitions. This illustrates that caution is needed
when using a particular species as a tracer of certain phenomenon simply by
observing a single transition.

Figure 28: PC1 coefficients for selected transitions are shown in a descending
order. For dense gas tracers, shock tracers, and sulfur species, upper-state
energies are color-coded. Values from 0 K to 100 K are shown from cyan to dark
blue (dense gas tracers), from lime to dark green (shock tracers), and pink to
purple (sulfur species). Colors of CO, continuum, RRL, vibrationally excited
lines, and COMs indicate categories, not energies.

Figure 29: The same as Figure 28, but for PC2.

Figure 30: The same as Figure 28, but for PC3.

Figure 31: The same as Figure 28, but for PC4.

Figure 32: The same as Figure 28, but for PC5.  Figure 33: The same as Figure
28, but for PC1-5 of the rest of transitions. Dense gas tracers are shown with
blue, reactive ions with violet, carbon chains with pink, and others with
grey.  Figure 34: PC2 and PC3 coefficients are shown in abscissa and ordinate,
respectively, for dense gas tracers and RRLs. Other transitions are also shown
as grey arrows. Arrow coloring is the same as in Figure 28. Transitions in the
first quadrant (filled with light yellow) are distributed near young
starbursts ($<10^{5}\,$yr) while transitions in the fourth quadrant (filled
with light green) are near well-developed starbursts ($>10^{5}\,$yr) as
discussed in Section 5.3.  Figure 35: Same as Figure 34, but showing PC2 and
PC4 for shock tracers.  Figure 36: Same as Figure 34, but showing PC4 and PC5
for CO isotopologues, reactive ions, and carbon chains. Arrows show CO
isotopologues in black, ions in violet, and carbon chains in pink. Transitions
with outflow-like features should lie in the quadrant marked with light blue
as discussed in Section 4.5.  Figure 37: Same as Figure 36, but showing PC3
and PC5. Transitions with outflow-like features should lie in the quadrant
marked with light blue as discussed in Section 4.5. These transitions are the
same as ones in the second quadrant in Figure 36 except for C3H+(7-6).

###  4.6 Similarities between species

The similarity between a pair of transitions can be evaluated by the
correlation coefficient of their integrated intensity maps, or the relevant
element of the matrix $\bm{\Sigma}$ used for the PCA999We also make the
correlation matrix among all the transitions and continua available at
http://github.com/nanaseharada/alchemi_pca.. Figure 38 shows the correlation
coefficients (or similarities) between a reference transition and other
transitions, for reference transitions of 12CO(1–0), HC3N(37–36),
CH3OH($2_{K}$–$1_{K}$) at 96.741 GHz, and H39$\alpha$. Transitions highly
correlated with 12CO(1–0) are other transitions of CO isotopologues, $J=1-0$
or $N=1-0$ transitions of HCN, HCO+, CN, and CCH, and CS(2–1) (Figure 38).
These transitions are all ubiquitous in translucent ($n\sim 10^{3}$ cm-3)
molecular gas. On the other hand, transitions that have similar distributions
as that of HC3N(37–36) are highly excited. The examples are high-excitation
transitions of other HC3N, sulfur-bearing species (rarer isotopologues of CS
and SO, H2S), SiO, and CH2NH. Some of these transitions are known to trace
shocks (SiO), but others are considered as hot core tracers. Transitions that
are strongly correlated with CH3OH(2–1) are low-excitation transitions of
HNCO, HOCO+, and OCS(8–7). Methanol and HNCO are both known slow shock
tracers, and HOCO+ is shown to be enhanced in shocks in NGC 253 due to the ice
sputtering of CO2 (Harada et al., 2022). OCS may also be enhanced in shocks if
OCS is a major constituent of ice. H39$\alpha$, other RRLs, and the 3-mm
continuum are all very strongly correlated, all of which are known to be
related to star formation. High-excitation radicals such as HOC+(3–2),
CN(3–2), CCH(4–3) are also highly correlated with RRLs. These species are
radicals, and are abundant in PDRs or cosmic-ray dominated regions.

Figure 38: Covariance between reference transitions and other transitions. The
reference transitions are 12CO(1–0) (top left panel), HC3N(37–36) (top right
panel), CH3OH($2_{K}-1_{K}$) (bottom left panels), and H39$\alpha$ (right
panels). Transitions that appear close to 1 have similar distributions as that
of the reference transition.

###  4.7 GMCs in principal components

Figure 39 shows the nine GMCs in their scores of PC1 through PC5. We measured
the PC scores from the hexagonal pixels at the GMC positions. Each panel has
dashed lines at the PC score of zero to indicate the mean score of the entire
CMZ. All GMCs have positive PC1 scores since the GMCs are emission peaks by
definition, and the PC1 score reflects the overall strength of the emissions,
approximately the ISM column density, at each position. PCs 2–5 are defined to
be orthogonal to PC1. Therefore, they disregard the effect of overall emission
strengths or ISM quantity and should only reflect the qualities of the ISM.
The nine GMCs are clearly different in PCs 2–5. For example, they spread in
the PC2–PC3 plane with the range of PC scores matching that in the entire CMZ.
Thus, the nine GMCs are distinct in their characteristics unrelated to their
bulk emission strengths.

We categorized the GMCs into three groups on the basis of their PC scores
using the k-means clustering method (KMeans class in the scikit-learn
package). The k-means method is one of the simplest methods of categorization
of data points into a desired number of clusters. This algorithm minimizes the
following quantity: $\sum^{k}_{i=1}\sum_{z\in S_{i}}|z-\mu_{i}|^{2}$ for a set
of clusters S = {$S_{1}$, $S_{2}$, …, $S_{k}$} where $\mu_{i}$ is a mean of a
set $S_{i}$ (MacQueen, 1967), that is, the sum of a square of a distance
between a data point and a cluster mean for each data point. Computationally,
this minimum is found by 1) randomly assigning data points into
clusters101010We use k-means++, the default method in the KMeans class, where
the initial assignment is not totally at random, but chosen so that the
centers of mass are not close to each other., 2) calculating the center of
mass within each cluster, 3) re-assigning each data point into a cluster
having the closest center of mass, and 4) iterating Step 2 and 3 until it
converges. In our case, the distances are Euclidean in the five-dimensional PC
space.

Using the cluster number of three based on the inspection of Figure 39, we
categorized the nine GMCs into the groups of GMCs [1,2,7,8,9], [3,6], and
[4,5] and refer to them as categories A, B, and C, respectively. This
categorization indeed appears to correspond to physical differences between
the GMCs. The first GMC category (A) contains the locations where Class I
methanol masers have been detected, which indicates shocked locations. The
second category (B) contains GMCs with relatively high intensities of complex
molecules. The third category (C) has locations where RRLs are strong. We
discuss in Section 5.3 that the difference between the second and third
categories is likely attributed to the evolutionary stages of starbursts.

Figure 39: PC scores in GMCs. Different colors show three different groups
categorized with the k-means clustering method. Each category is referred to
as A (red), B (green), and C (blue) in the text.

##  5 Discussion

###  5.1 Overall gas content and star formation

We suggest in Section 4.4 that PC1 is mostly associated with the molecular
hydrogen column density. Figure 24 and Figure 28 show that all transitions
have positive PC1 coefficients (in the range of 0.042–0.091, see Table 4),
which means that all transitions are correlated with the column density.
Understandably, molecular transitions are correlated with the total column
densities because the emission is expected to be stronger where there is more
molecular gas. RRLs should also somewhat correlate with the total column
density because when there is molecular gas, star formation is also expected
in general (e.g., Kenicutt-Schmidt Law, Kennicutt, 1998). Although star
formation and molecular gas content are not always co-spatial depending on the
GMC evolutionary stage and spatial resolution (e.g., Schinnerer et al., 2019),
their correlation is relatively strong in the CMZ of NGC 253 where its star
formation is relatively young. We note that this strong correlation may also
be attributed to the edge-on inclination of this target.

Although PC1 generally represents the overall molecular gas column density,
transitions with the highest coefficients for PC1 are not what we expect given
this relation. It is not 12CO(1–0) that has the highest coefficients for PC1
among CO isotopologues; instead these are the $J=3-2$ transitions of C18O and
C17O. This may be partly because 12CO(1–0) is optically thick, and does not
trace the overall gas column density properly. If the optical depth is the
major effect, PC1 appears to trace the overall molecular gas better than
12CO(1–0). On the other hand, there is a factor that could cause deviation of
PC1 from overall molecular gas. Many high-excitation or optically-thin
transitions or transitions from complex molecules show concentration on the
central starburst regions with fewer transitions of optically thick, low-
excitation, and extended transitions. That results in PC1 being biased towards
central starburst regions, not purely indicating the total column density.
Given the above factors, PC1 in general is a good tracer of overall molecular
gas because many transitions have similar PC1 coefficients as C18O and C17O
(3–2).

###  5.2 Dense gas tracers and star formation

Our PCA shows that PC2 is related to excitation (Figure 27). PC2 coefficients
are positive and large for RRLs (Figure 29), which is also evident from the
resemblance between the positive part of the PC2 score map (Figure 24) and
moment 0 maps of RRLs. PC2 scores are the highest in the GMCs of category C
(GMCs 4 & 5), and relatively high in category B (GMCs 3 & 6), as shown in
Figure 39. They are GMCs with massive-star formation judging from the RRLs. In
addition to RRLs, dense gas tracers with high upper-state energies ($\gtrsim
50\,$K) also have large positive coefficients for PC2 (Figure 29). Meanwhile,
their lower excitation transitions (e.g., HCN $J=1-0$) do not have such a
strong correlation to RRLs (e.g., see their PC2 coefficients in Figure 29),
and have negative PC2 coefficients comparable to the PC2 coefficients of CO
isotopologues.

The correlation between star formation and low-excitation transitions with
moderately high critical densities ($n_{\rm crit}\sim 10^{5}\,$cm-3 in an
optically-thin case, e.g. HCN $J=1-0$) is not very tight in our study,
although such a correlation has been reported from extragalactic observations
in single-dish resolution (e.g., Gao & Solomon, 2004). In fact, this
correlation appears unlikely to be tight on $<100\,$pc scale in general.
Galactic observations on the 1–10 parsec scale (in Orion A, Orion B, W51,
W33(OH), and Perseus) have already shown that HCN(1–0) does not trace only
high-density regions, but can also originate from translucent clouds
($n\gtrsim 10^{3}\,$cm-3) (Pety et al., 2017; Watanabe et al., 2017; Kauffmann
et al., 2017; Nishimura et al., 2017; Tafalla et al., 2021). This may be in
conflict with the results by Wu et al. (2010) who claim that the correlation
between the HCN luminosity and star formation rate does not break down within
$<100\,$pc. However, this difference may come from the fact that they focused
on dense clumps where stars are already forming in the field of view of 4 pc
on average. Although the mean mass-weighted density of the molecular gas in
NGC 253 CMZ is $\sim 10^{3.5}-10^{4}\,$cm-3 (Tanaka et al., 2023), higher than
in the GMCs in the Galactic disk ($10^{2}-10^{3}\,$cm-3; Pety et al., 2017),
gas detected with low-excitation “dense-gas tracers” should be low enough in
density not to immediately form stars within the timescale for the gas to
orbit to different locations in our observations. For example, it takes $\sim
10^{6}\,$yrs for a cloud to move 100 pc if it is moving with a velocity of 100
km s-1, a few times the free-fall time of a cloud with $n=10^{4}\,$cm-3. Note,
though, that the star formation threshold is likely higher in the CMZ of NGC
253 than that in the Galactic center and the Galactic disk (Tanaka et al.,
2023).

The similarity between the high-$J$ dense gas tracers and star formation has
at least two possible interpretations. One is that dense gas is intimately
related to star formation. Because star formation takes place in dense gas,
this high-density gas traced with high-$J$ transitions may co-exist with star
formation. Alternatively, the high-excitation dense gas may have resulted from
star formation. Star formation heats the gas, and molecules around young stars
may be more excited owing to the high temperature. The distinction between
these two scenarios is discussed further in Section 5.3.

###  5.3 Evolutionary stages of proto super star clusters

There are some works reported on the evolutionary different stages among GMCs
in the center of NGC 253, and it is of interest to see whether it has an
impact on the physical and chemical properties that appear in our data. Rico-
Villas et al. (2020) showed different ages among these GMCs using the ratios
between luminosities of protostars and zero-age main-sequence stars
(ZAMS)111111This age estimation by Rico-Villas et al. (2020) applies to young
clusters up to the ZAMS stage., concluding that GMCs 3 (their clusters 1, 2,
and 3) and 6 (their cluster 14) have proto-super star clusters that are
younger (proto-super-star clusters $<10^{5}$ yr) than most of the stars in
GMC5 (their clusters 8-13; mostly ZAMS with $t_{\rm age}=10^{5}-10^{6}\,$yr)
although there are a couple of clusters that are young in GMC5 (their clusters
8 and 13). Mills et al. (2021) also concluded that their clusters 1, 2, 3, 8,
13, and 14 (our GMCs 3 and 6, a part of GMC 5) are younger than other clusters
with the assumption that the stellar mass to dust mass ratios increases with
age.

The literature above indicates that GMCs in category B defined in Section 4.7
host younger proto-super-star clusters, while GMCs in category C are more
developed clusters (GMCs 4 and 5). It is therefore natural to suspect that the
evolutionary stages of SSCs alter chemical characteristics in GMCs. The
difference that appears in PC3 likely contributes to these characteristics in
the PC score categorization. Positive PC3 scores seem to trace GMCs that have
relatively young star formation (because its tracers have positive
coefficients for PC3 as seen in the 1st quadrant of Figure 34). There,
relatively large molecules like cyanopolyynes (HC2n+1N) especially
vibrationally-excited lines are abundant. PC3 is also positive in GMC3, where
Rico-Villas et al. (2020) show the youngest clusters. Meanwhile, negative PC3
values likely trace more well-developed star formation (the 4th quadrant of
Figure 34). GMC5 and a part of GMC4 have negative PC3 scores, and they contain
more developed star formation ($>10^{5}$ yr).

We note that the age sequence proposed by Rico-Villas et al. (2020) and Mills
et al. (2021) is not consistent with that reported by Krieger et al. (2020),
who assumed that the HCN/HC3N ratios decrease with evolutionary stage owing to
the increasing density. Instead, their HCN/HC3N intensity ratios in the 1-mm
band tends to be higher in GMCs 5 and 6. We argue that the HCN/HC3N ratios may
not follow the same stages as those measured by Rico-Villas et al. (2020).
Krieger et al. (2020) used the HCN/HC3N ratios to estimate the ratios between
the very dense gas and the dense gas. This ratio likely increases as a
molecular cloud evolves to form stars, but it is unclear if the density
continues to increase after stars form. In addition, the age differences among
all the super-star clusters in NGC 253 CMZ may not be enough to cause chemical
differentiation as already noted by Krieger et al. (2020).

These young and more developed starbursts embedded in the gas are
characterized by different chemistry. In GMC6, high-excitation transitions of
cyanopolyynes (HC3N and HC5N), CH2NH, and shock tracers have high intensities.
On the other hand, in GMC5, there are high intensities of high excitation
transitions of CN and CCH. These features can be attributed to differences in
shock strengths/frequencies and/or the presence of hot molecular cores or
PDRs. Embedded young starbursts may encounter a stronger influence from the
interaction between the starburst and its surrounding medium, which causes
shocks. Cyanopolyynes and CH2NH may also be enhanced from the ice sublimation
from shocks. Dense and hot molecular cores heated from star formation can also
increase abundances of cyanopolyynes, CH2NH, and methanol. These GMCs are hot
and dense, having a similar environment as that of protostellar cores (“hot
cores”), although individual protostellar cores are too small to resolve in
our beams. If shocks or hot molecular cores are the only cause of this
chemical differentiation, it does not fully explain the strong correlation
between high-excitation radicals (CCH 4–3, CN 3–2) and RRLs. Dense PDRs around
GMC5 may explain this trend as they can both increase the abundances of
radicals and cause high excitation. Galactic PDRs around high-mass star-
forming regions (e.g., Taniguchi et al., 2022) show decreased abundances of
cyanopolyynes likely due to photodissociation. While low-excitation CN and CCH
can be emitted from extended low-density gas or cosmic-ray dominated regions,
we suggest that higher-excitation transitions such as CCH(4-3), CN(3-2),
HOC+(3-2) may be relatively uncontaminated tracers of PDRs in the vicinity of
starburst. Higher stellar-to-dust mass ratios in some clumps in GMC5 seen by
Mills et al. (2021) may support this scenario because GMCs with more ionizing
sources (massive stars) per unit mass are capable of creating more PDRs within
them.

There are possible arguments against the strong influence of PDRs. For
example, Holdship et al. (2021, 2022) concluded that the total column
densities of GMCs are extremely high ($N_{\rm H}\sim 10^{24}\,$cm-2), and PDRs
should not have significant contributions assuming that the media are largely
uniform, and not clumpy; instead, radicals can be enhanced with high cosmic-
ray ionization rates ($\zeta\gtrsim 10^{-13}\,$s-1). We note that the presence
of PDRs, as previously claimed based on the detection of significant
abundances of HOC+ and CO+ (Martín et al., 2009), does not necessarily exclude
high cosmic-ray ionization rates. In fact, they are expected to be high in
starburst galaxies, and we do not dispute that cosmic-ray ionization rates are
in general much higher than that in the Galactic disk. Among the cosmic rays,
lower-energy ones are attenuated even with low column densities while higher-
energy cosmic rays are likely to influence a larger extent than individual GMC
scale. If there are cosmic-ray sources within GMCs, cosmic rays can cause
differences between individual GMCs. However, known supernovae lie outside of
GMCs (Ulvestad & Antonucci, 1997). Therefore, cosmic rays are unlikely to
account for chemical differences among GMCs shown in this work. Moreover,
cosmic rays have stronger influences in lower-density regions because the
chemistry roughly scales as $\zeta/n$, and molecules can be dissociated with
extremely high cosmic-ray ionization rates ($\zeta\gtrsim 10^{-13}\,$s-1). One
needs a reasonable explanation for the wide-spread presence of CO in low-
density regions ($n\sim 10^{3}\,$cm-3) if the cosmic-ray ionization rates are
that high over a large extent. Therefore, we suggest that PDRs can be a main
driver of chemistry in some regions even if cosmic-ray ionization rates are
also high.

How feedback is acting in these GMCs is an interesting question. Simulations
have modeled quenching mechanisms of star formation in GMCs. Stellar feedback
may include UV, optical, and infrared radiation pressure, stellar winds,
supernovae, photoionization, and photoelectric heating (Grudić et al., 2018).
Photoionization feedback may be important in low-column-density clouds, but
not in high-column-density clouds (Fukushima & Yajima, 2021). Radiation
pressure may instead be more important in many cases (Grudić et al., 2018).
Therefore, the incident radiation field producing the PDR-like emission
characteristics that we see in our results may not directly measure the total
feedback strength. In other words, a larger degree of feedback may be acting
in young embedded proto-super-star clusters (e.g., Levy et al., 2021) than in
clusters showing PDR-like features. However, we expect such PDRs when
molecular clouds are already disrupted by feedback and UV photons can travel
spatially further away than when they are embedded (e.g., the ISM is made more
porous).

We note that a strong continuum in GMC5 (so-called “TH2” from Turner & Ho,
1985) causes absorption features for some transitions such as H13CO+(1-0),
HC13N(1-0), and SiO(2-1). The absorption features are not obvious for other
transitions, and our PCA results are unlikely to be significantly affected.
But this effect should be checked against higher angular resolution data.

###  5.4 Origins of shocks and relationship with star formation

The clustering analysis in the PC space found GMCs 1, 2, 7, 8, and 9 in the
same group (category A, Figure 39). These are locations with detected Class I
methanol masers (Humire et al., 2022), indicating the presence of shocks.
Shocks are likely to be at intersections between orbits of bar ($x_{1}$
orbits) and circumnuclear ring ($x_{2}$ orbits) (Meier et al., 2015; Harada et
al., 2022; Humire et al., 2022). Transitions that are enhanced in these
regions are low-$J$ transitions of CH3OH, HNCO, and HOCO+ , tracers of weak
shocks. Although star formation could also cause shocks, it is unlikely to be
causing shocks here because there is no obvious evidence of star formation.
Meanwhile, these shocks or cloud collisions could promote future star
formation through the compression of molecular clouds (e.g., Böker et al.,
2008; Harada et al., 2019).

A well-known tracer of shocks, SiO, does not show strong enhancement near
orbital intersections according to our result (e.g., PC2; Figure 29). High
excitation transitions of SiO (e.g., J=8-7, 6-5) correlate relatively well
with RRLs when PC2 and PC4 are considered (Figure 35). In addition to these
SiO lines, some HNCO transitions are relatively well-correlated with RRLs.
They are the $K_{a}\geq 1$ transitions of
HNCO$(J_{K_{a},K_{b}}-J^{\prime}_{K_{a}^{\prime},K_{b}^{\prime}})$, which
Churchwell et al. (1986) suggested to be radiation-pumped, not collisionally
pumped, in Galactic hot cores. For them to be radiation-pumped, a strong far-
infrared radiation field needs to be present, which may come from embedded
protostars. If these HNCO $K_{a}=1$ transitions originate from protostars,
they do not necessarily need shocks. It has been suggested that HNCO can be
abundant near young protostars (Quan et al., 2010). We note here that these
SiO and HNCO transitions are similar to RRLs in contributions to PC2 and PC4,
but not to PC3. This trend implies that these transitions are associated with
younger starbursts, not developed ones (Section 5.3).

High-$J$ SiO transitions appear to be the only transitions of shock tracers
that is centrally concentrated, if HNCO $K_{a}=1$ transitions are not related
to shocks. This trend suggests the prevalence of strong shocks in the central
regions. One has to be careful making this conclusion because we did not image
high-$J$ HNCO transitions ($J>10$) due to the contamination criteria discussed
in Section 2.2, but the contamination to these transitions from neighboring
transitions does not seem very severe except for $J=12$ (Huang et al., 2023).
However, it seems also true that SiO is more concentrated around the starburst
regions, GMCs 4-6, compared with CH3OH, HNCO, and HOCO+ of similar upper-state
energy levels. The only exception to it is when the low-$J$ SiO line is
affected by absorption against a strong continuum source. An LVG analysis by
Huang et al. (2023) found that SiO/HNCO column density ratios are higher in
GMCs 4 and 6 compared with those in GMCs 3 and 7 although this difference is
not significant considering errors and there are some in the outskirts (GMCs 2
and 9) that have high SiO/HNCO ratios as well. Note that they did not model
GMC5 because of the self-absorption. It is not clear whether fast shocks are
more dominant than slow shocks in this central starburst region as Huang et
al. (2023) concluded, but there are indeed strong shocks near the starburst
regions.

What are the origins of the shocks that caused the high-$J$ SiO line emission
near the central starbursts? One scenario is again cloud collisions. Levy et
al. (2022) suggested another $x_{2}$ orbits connecting GMCs 3-6 in addition to
the larger one shown in Figure 1. If this orbit crosses with some $x_{1}$
orbits, cloud collisions could occur. However, features from such cloud
collisions should appear in GMCs 3 and 6 where $x_{1}$ orbits likely crosses
with the inner $x_{2}$ orbits instead of GMCs 4-6. Another explanation is that
the medium in NGC 253 CMZ is in general turbulent, and high-excitation
transitions are seen simply because the gas is excited near the starbursts.
This scenario is hard to exclude, but it does not solely explain variations
among SiO and other shock tracers. The other cause of shocks is the starburst
itself. Cluster-scale outflows have been found in high-angular-resolution
observations through P-Cygni profiles in GMCs 4 and 6 (Levy et al., 2021) and
through broad line wings of H40$\alpha$ in GMCs 4 and 5 (Mills et al., 2021).
Therefore, outflows from young stellar objects may cause strong shocks, which
yield high-excitation SiO emission. In addition to such cluster-scale
outflows, fast shocks are common around young massive stars from stellar
winds. We consider this starburst-induced shock scenario is most plausible for
now because of the signatures of outflows, but this claim is still
speculative.

###  5.5 Vertically extended structures likely associated with outflows

Outflows on the galactic scale have been detected in the center of NGC 253 in
previous studies, and some of them contribute to the vertically extended
structures or emission at high galactic altitudes seen in PC4 and PC5 (Turner,
1985; Bolatto et al., 2013; Walter et al., 2017; Krieger et al., 2019). In
Section 4.5, we show that CO isotopologues, $J$ or $N$=1–0 transitions of HCN,
HCO+, CN, CS, CCH, HOC+, and H3O+($3_{2,0}-2_{2,1}$) have emission in this
structure. Among them, HOC+ (Figure 10) and H3O+ (Figure 20) have relatively
weak overall intensities but still show this extended structure. This is
likely because these species are reactive ions that are enhanced by strong UV
radiation or cosmic rays (Harada et al., 2021; Holdship et al., 2022). This
result suggests a high ionization degree in these outflows.

In general, there are multiple possible mechanisms creating starburst-driven
and AGN-driven outflows such as thermal energy, radiation, cosmic rays, and
radio jets (Veilleux et al., 2020). While our observations do not suggest
which mechanisms are most effective in producing energy, they do suggest an
influence of cosmic rays because H3O+ is preferentially enhanced with cosmic
rays compared to UV photons (Holdship et al., 2022).

###  5.6 Comparisons with other PCA studies

We compare our study with previous image-based PCA studies conducted for
nearby galaxies (Meier & Turner, 2005, 2012; Chidiac, 2020; Saito et al.,
2022) and Galactic molecular clouds (Ungerechts et al., 1997; Gratier et al.,
2017). Comparisons are not necessarily straightforward because the number of
molecules and transitions in this work is unprecedented. The results of PCAs
are affected by the transitions used, their intensities, and the normalization
method. In addition, the field of view and angular resolutions are different
in all studies. Here, our discussion takes these factors into consideration.
Similar to our results, most studies found positive PC1 coefficients for all
input variables used in PCA, except for Saito et al. (2022). This is because
they used not only molecular transitions or continuum but also other
quantities such as the [Siii]/[Sii] ratio and [Ci] intensity that do not
necessarily trace overall molecular content. Note also that Gratier et al.
(2017) showed that the interstellar radiation field anticorrelates with all
the molecular emission intensities.

Meier & Turner (2005) observed the center of the modest starburst galaxy IC
342, with a similar field of view and twice larger beam size (50 pc) in
comparison with our study. They used 12CO(1-0), C18O(1-0), CCH, C34S(2-1),
N2H+(1-0), CH3OH, HNCO, HNC, HC3N(10-9), SO, and continuum in the 3-mm band.
Their study shows that CH3OH and HNCO emission also peaks at the outskirts of
the CMZ, and has a high absolute value in their PC2. This is in agreement with
our study; high-excitation transitions have high absolute values of PC2 while
CH3OH and HNCO have the opposite sign. In the 3-mm band observations of Meier
& Turner (2005), there were no such high-excitation transitions. Instead, C34S
and CCH showed different degrees of contribution in PC2 from CH3OH and HNCO
coming from high column density regions (C34S) or PDRs (CCH). In our study,
these transitions have somewhat different coefficients for PC2 compared to
CH3OH and HNCO, but the difference is not as large as for the high-excitation
transitions. Similar results are found by Meier & Turner (2012) in Maffei 2.

In the study of M82 by Chidiac (2020), they used CCH, CN, CS, HCN, HCO+, HNC,
and 13CO in the 3-mm band. Because of the limited number of transitions, 13CO
and other transitions had different contributions to PC2 likely because of the
difference in the critical densities.

The nature of PCA in Saito et al. (2022) is quite different from ours in the
sense that they focused on the region around the AGN in NGC 1068. They used
sixteen transitions in the 3-mm band line survey together with [Ci] and
[Siii]/[Sii] ratio from archival data. They separated out the circumnuclear
disk in PC1 and the AGN-driven outflow in PC2. In their result, CN, HNC, and
CCH showed high contribution from the outflow while N2H+, HC3N, and CO
isotopologues show the opposite trend according to the categorization from
PC2. In our result, CN, HNC, and CCH also have outflow features, but so do the
CO isotopologues, HCN, and HCO+.

Ungerechts et al. (1997) mapped the Orion BN-KL regions with 12CO, CS, HCN,
HNC, HCO+ and their 13C isotopologues, C18O, C34S, N2H+, CN, CCH, HC3N,
CH3CCH, C3H2, SO, and CH3OH. Because their line widths are much narrower than
ours, their observations spectrally resolved hyperfine structure lines of N2H+
and various $K$-ladders of CH3CCH. Their field of view is at the pc scale with
the beam size of 0.1 pc and so is qualitatively different from the
extragalactic studies. Their results separate among transitions enhanced near
the high-mass star-forming regions of BN-KL such as HC3N and CH3OH, ones
enhanced around the Ridge such as CS and CN, and ones enhanced near a
relatively cold cloud such as N2H+. This distinction is not obvious in our
result. One reason may be the difference in spatial resolution. Our 27-pc beam
is not small enough to separately resolve dense clouds and star-forming
regions within a GMC. Another reason may be that the difference between dense
clouds and star-forming regions is better highlighted in our study by the
difference between higher-$J$ HC3N and N2H+ rather than HC3N(10-9) and N2H+ in
the 3-mm band. It is also important to note that NGC 253 contains more extreme
star formation compared with Orion-KL.

Gratier et al. (2017) conducted a PCA study using maps of 12CO, 13CO, C18O,
CS, HCN, HCO+, SO, CN, HNC, CCH, N2H+, and CH3OH in the Orion B molecular
cloud with the linear resolution of $0.05-0.07\,$pc covering the area of $\sim
6$ pc$\times 8\,$pc. In addition to these line intensities, they also examined
the correlation between these molecular line emission intensities and
logarithms of the H2 column density $N_{\rm H2}$ obtained from dust continuum
flux, volume density $n_{\rm H}$ from the dust mass and core sizes, and UV
radiation field normalized by the mean value $U/\bar{U}$ from PAH emissivity.
Their coefficients for PC2 show a large difference between N2H+ or CH3OH vs.
CO isotopologues. This is likely due to the difference in the density and
temperature as N2H+ and CH3OH should be more abundant in cold and dense
regions. Our result also shows a similarity between N2H+ and CH3OH. They also
showed a contrast between CCH, CN and other species. CCH and CN are known as
PDR tracers, and their PC3 should have contributions from the UV radiation
field. In fact, they suggested a tight positive correlation between PC1 and
$N_{\rm H2}$, a positive but weak correlation between PC2 and $n_{H}$, and a
moderate positive correlation between PC3 and the UV radiation field.

The above studies show that PCA is a useful way to identify coupled
chemistries, reveal the underlying dynamics driving them, and summarize rich
molecular line datasets. Extragalactic PCA studies seem to be consistent in
showing distinct contributions in low-excitation transitions of methanol and
HNCO, especially near the intersection of the bar orbit and the central
region. Galactic and extragalactic PCAs may show different categorizations of
transitions, likely due to the difference in spatial scale.

##  6 Summary

We have presented integrated-intensity images of isolated line emission from
the ALCHEMI survey, a wide-frequency imaging spectral scan of the CMZ of the
starburst galaxy NGC 253. With these images, we have performed a principal
component analysis (PCA) as a statistical approach to reduce the dimensions of
the dataset to investigate underlying physical conditions that cause
differences in line intensities, and to relate these physical conditions to
starburst activities.

We succeeded in reducing the dimension of the data space from 150 variables to
5 principal components through PCA while retaining more than 96% of the
original information. It reflects the fact that spatial distributions of
various sub/millimeter emissions are highly correlated in the CMZ. One needs
only a handful of uncorrelated parameters to explain the most of the
standardized emission maps. Our PCA found that the nine GMCs in the CMZ are
distinct from each other in their qualities after subtracting the effect of
their overall emission strengths. We interpreted the PCs in physical terms and
assessed various emissions with their contributions to the PCs. We list our
findings below.

  * •

As expected, the largest factor that changes line intensities is the amount of
molecular gas. There is an overall correlation among all the molecular
emission lines, radio recombination lines (RRLs), and continuum images. This
is indicated by the first principal component (PC1). Locations with larger
molecular gas column densities in general emit stronger molecular transitions.
These locations are also likely to have higher star formation rates (Kennicut-
Schmidt Law), which causes the higher intensity of RRLs. We also note that
there are some biases that cause PC1 to deviate from an accurate measure of
the total column density as discussed in the text. The deviation from this
correlation appears in PCs 2 or higher.

  * •

High-excitation dense gas tracers and RRLs have similar spatial distributions
and large positive contributions to PC2. Meanwhile, low-excitation ($J=1-0$)
dense gas tracers have negative coefficients for PC2; indeed, we found strong
correlations between contributions to PC2 and critical densities as well as
upper-state energies of the transitions of dense gas tracers. The difference
between low-excitation dense gas and star formation tracers suggests that
these low-excitation dense gas tracers also emit from less dense regions than
indicated by their critical densities. This result is consistent with Galactic
large-scale cloud observations in the literature.

  * •

Evolutionary stages of giant molecular clouds cause differences in the
chemistry. Young starbursts (GMCs 3 and 6) tend to have peaks of high-
excitation transitions of HC3N and COMs, while more developed starbursts (GMCs
4, 5) have peaks of RRLs and high-excitation transitions of radicals (CN,
CCH). This difference appears as positive and negative contributions to PC3,
respectively. The lack of HC3N and COMs in developed starbursts may be caused
by photodissociation, which shows the effect of stellar feedback. These PDRs
near the starbursts can be traced by high-excitation transitions of CN and
CCH. These transitions should be observable for some sources at $z=2-3$ in the
3-mm band, and can be useful for studying the ISM variation during $z=0-2$.

  * •

Comparisons between PC2 and PC4 show that low-excitation shock tracers (CH3OH,
HNCO, HOCO+, and OCS) do not correlate with star formation, and likely trace
shocks caused by cloud collisions of $x_{1}$ and $x_{2}$ orbits. On the other
hand, high-excitation transitions of SiO may originate from star formation as
their emissions are from locations of cluster-scale outflows in the
literature.

  * •

Vertically-extended structures, some of which must originate from outflows,
have emissions of CO isotopologues, CCH, HOC+, H3O+, HCN, HCO+, CS, and CN.
The enhanced emission of HOC+, H3O+ in outflows indicates that these outflows
are highly ionized, possibly due to cosmic rays.

  * •

We have also made the first extragalactic detection of HCNH+, and the first
extragalactic images for C3H+ NO, and HCS+. From the characteristics of these
species and their distributions, we suggest that C3H+ emission preferentially
comes from PDRs near the starburst while HCNH+ emits in cold clouds.

We have identified some categories of molecular transitions that are
characteristics of physical conditions through the extensive spectral scan
ALCHEMI. While not many observations can afford a wide-band scan like ALCHEMI,
overall chemical structures can be found by observing sets of transitions such
as CO isotopologues, shock tracers (SiO, CH3OH, HNCO), high-excitation dense
gas tracers (HC3N) or COMs, high-excitation radicals (CN, CCH) as starburst
tracers. Reactive ions such as HOC+ and H3O+ are also useful in tracing highly
ionized regions as our previous studies found.

##  7 acknowledgements

We appreciate the constructive comments by the anonymous referee. We thank
Francesco Costagliola for setting the path to the ALCHEMI survey. We are also
grateful to the ALMA staff at JAO and ARCs for their dedicated work to keep
ALMA running, which made this survey possible. This paper makes use of the
following ALMA data: ADS/JAO.ALMA#2017.1.00161.L, ADS/JAO.ALMA#2018.1.00162.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 National
Radio Astronomy Observatory is a facility of the National Science Foundation
operated under cooperative agreement by Associated Universities, Inc. The
Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Data analysis
was in part carried out on the Multi-wavelength Data Analysis System operated
by the Astronomy Data Center (ADC), National Astronomical Observatory of
Japan. N.H. acknowledges support from JSPS KAKENHI grant No. JP21K03634, and
the ALMA Japan Research Grant of NAOJ ALMA Project, NAOJ-ALMA-304. K.S.
acknowledges the support from the Ministry of Science and Technology (MOST) of
Taiwan through the grant MOST 111-2112-M-001-039. This work was supported by
NAOJ ALMA Scientific Research grant No. 2021-18A. L.C. acknowledges financial
support through the Spanish grant PID2019-105552RB-C41 funded by
MCIN/AEI/10.13039/501100011033. S.V., M.B., and K.-Y.H. acknowledge support
from the European Research Council (ERC) Advanced Grant MOPPEX 833460. K.K.
acknowledges the support by JSPS KAKENHI Grant Number JP17H06130 and the NAOJ
ALMA Scientific Research Grant Number 2017-06B.

##  Appendix A Transitions properties of blended with the same species

Transitions of species such as CN, CCH are blended in our observations due to
large line widths. We combined these blended neighboring transitions to make
velocity-integrated images of these species. Table 2 lists these transitions.

Table 2: Spectroscopic Properties for transitions blended with the same
species

Species | Transition | $\nu_{rest}^{(a)}$ |  $E_{\rm up}$(b) |  $log(A_{\rm ul})$(c)
---|---|---|---|---
|  | (GHz) | (K) | (s-1)
CN | N=1-0, J=1/2-1/2, F=1/2-3/2 | 113.144 | 5.43 | -4.98
| N=1-0, J=1/2-1/2, F=3/2-1/2 | 113.170 | 5.43 | -5.29
| N=1-0, J=1/2-1/2, F=3/2-3/2 | 113.191 | 5.43 | -5.18
CN | N=1-0, J=3/2-1/2, F=3/2-1/2 | 113.488 | 5.45 | -5.17
| N=1-0, J=3/2-1/2, F=5/2-3/2 | 113.491 | 5.45 | -4.92
| N=1-0, J=3/2-1/2, F=1/2-1/2 | 113.500 | 5.45 | -4.97
| N=1-0, J=3/2-1/2, F=3/2-3/2 | 113.509 | 5.45 | -5.28
CN | N= 2-1, J=3/2-3/2, F=1/2-1/2 | 226.287 | 16.3 | -4.99
| N= 2-1, J=3/2-3/2, F=1/2-3/2 | 226.299 | 16.3 | -5.08
| N= 2-1, J=3/2-3/2, F=3/2-1/2 | 226.303 | 16.3 | -5.38
| N= 2-1, J=3/2-3/2, F=3/2-3/2 | 226.315 | 16.3 | -5.00
| N= 2-1, J=3/2-3/2, F=3/2-5/2 | 226.332 | 16.3 | -5.34
| N= 2-1, J=3/2-3/2, F=5/2-3/2 | 226.342 | 16.3 | -5.50
| N= 2-1, J=3/2-3/2, F=5/2-5/2 | 226.360 | 16.3 | -4.79
| N= 2-1, J=3/2-1/2, F=1/2-3/2 | 226.617 | 16.3 | -4.97
| N= 2-1, J=3/2-1/2, F=3/2-3/2 | 226.632 | 16.3 | -4.37
| N= 2-1, J=3/2-1/2, F=5/2-3/2 | 226.660 | 16.3 | -4.02
| N= 2-1, J=3/2-1/2, F=1/2-1/2 | 226.664 | 16.3 | -4.07
| N= 2-1, J=3/2-1/2, F=3/2-1/2 | 226.679 | 16.3 | -4.28
| N= 2-1, J=5/2-3/2, F=5/2-3/2 | 226.874 | 16.3 | -4.02
| N= 2-1, J=5/2-3/2, F=7/2-5/2 | 226.875 | 16.3 | -3.94
| N= 2-1, J=5/2-3/2, F=3/2-1/2 | 226.876 | 16.3 | -4.07
| N= 2-1, J=5/2-3/2, F=3/2-3/2 | 226.887 | 16.3 | -4.56
| N= 2-1, J=5/2-3/2, F=5/2-5/2 | 226.892 | 16.3 | -4.74
| N= 2-1, J=5/2-3/2, F=3/2-5/2 | 226.905 | 16.3 | -5.95
CN | N= 3-2, J=5/2-5/2, F=3/2-3/2 | 339.447 | 32.6 | -4.64
| N= 3-2, J=5/2-5/2, F=3/2-5/2 | 339.460 | 32.6 | -5.36
| N= 3-2, J=5/2-5/2, F=5/2-3/2 | 339.463 | 32.6 | -5.53
| N= 3-2, J=5/2-5/2, F=5/2-5/2 | 339.476 | 32.6 | -4.67
| N= 3-2, J=5/2-5/2, F=5/2-7/2 | 339.493 | 32.6 | -5.52
| N= 3-2, J=5/2-5/2, F=7/2-5/2 | 339.499 | 32.6 | -5.63
| N= 3-2, J=5/2-5/2, F=7/2-7/2 | 339.517 | 32.6 | -4.60
| N= 3-2, J=5/2-3/2, F=3/2-5/2 | 339.992 | 32.6 | -5.41
| N= 3-2, J=5/2-3/2, F=5/2-5/2 | 340.008 | 32.6 | -4.21
| N= 3-2, J=5/2-3/2, F=3/2-3/2 | 340.020 | 32.6 | -4.03
| N= 3-2, J=5/2-3/2, F=7/2-5/2 | 340.032 | 32.6 | -3.41
| N= 3-2, J=5/2-3/2, F=3/2-1/2 | 340.035 | 32.6 | -3.54
| N= 3-2, J=5/2-3/2, F=3/2-1/2 | 340.035 | 32.6 | -3.49
| N= 3-2, J=7/2-5/2, F=7/2-5/2 | 340.248 | 32.7 | -3.42
| N= 3-2, J=7/2-5/2, F=7/2-5/2 | 340.248 | 32.7 | -3.38
| N= 3-2, J=7/2-5/2, F=5/2-3/2 | 340.249 | 32.7 | -3.43
| N= 3-2, J=7/2-5/2, F=5/2-5/2 | 340.262 | 32.7 | -4.35
| N= 3-2, J=7/2-5/2, F=7/2-7/2 | 340.265 | 32.7 | -4.47
CH3CCH | J=5-4, K=4 | 85.431 | 128. | -6.65
| J=5-4, K=3 | 85.443 | 77.3 | -6.40
| J=5-4, K=2 | 85.451 | 41.2 | -6.28
| J=5-4, K=1 | 85.456 | 19.5 | -6.23
| J=5-4, K=0 | 85.457 | 12.3 | -6.21
CH3CCH | J=6-5, K=5 | 102.499 | 198. | -6.48
| J=6-5, K=4 | 102.517 | 133. | -6.22
| J=6-5, K=3 | 102.530 | 82.3 | -6.09
| J=6-5, K=2 | 102.540 | 46.1 | -6.02
| J=6-5, K=1 | 102.546 | 24.5 | -5.98
| J=6-5, K=0 | 102.548 | 17.2 | -5.96
CH3CCH | J=10-9, K=5 | 170.824 | 226. | -5.41
| J=10-9, K=4 | 170.854 | 161. | -5.36
| J=10-9, K=3 | 170.876 | 110. | -5.33
| J=10-9, K=2 | 170.893 | 74.0 | -5.30
| J=10-9, K=1 | 170.903 | 52.3 | -5.29
| J=10-9, K=0 | 170.906 | 45.1 | -5.29
CH3CCH | J=11-10, K=5 | 187.904 | 235. | -5.26
| J=11-10, K=4 | 187.936 | 170. | -5.22
| J=11-10, K=3 | 187.961 | 119. | -5.19
| J=11-10, K=2 | 187.979 | 83.0 | -5.17
| J=11-10, K=1 | 187.990 | 61.4 | -5.16
| J=11-10, K=0 | 187.994 | 54.1 | -5.16
CH3OH | 2( 1, 2)- 1( 1, 1) E | 96.739 | 12.5 | -5.59
| 2( 0, 2)- 1( 0, 1) A | 96.741 | 6.96 | -5.47
| 2(- 0, 2)- 1(- 0, 1) E | 96.745 | 20.1 | -5.47
| 2(- 1, 1)- 1(- 1, 0) E | 96.756 | 28.0 | -5.58
CCH | N= 1-0, J=3/2-1/2, F= 1- 1 | 87.284 | 4.19 | -6.59
| N= 1-0, J=3/2-1/2, F= 2- 1 | 87.317 | 4.19 | -5.82
| N= 1-0, J=3/2-1/2, F= 1- 0 | 87.329 | 4.19 | -5.90
| N= 1-0, J=1/2-1/2, F= 1- 1 | 87.402 | 4.20 | -5.90
| N= 1-0, J=1/2-1/2, F= 0- 1 | 87.407 | 4.20 | -5.81
| N= 1-0, J=1/2-1/2, F= 1- 0 | 87.446 | 4.20 | -6.58
CCH | N= 2-1, J=5/2-3/2, F= 3- 2 | 174.663 | 12.6 | -4.83
| N= 2-1, J=5/2-3/2, F= 2- 1 | 174.668 | 12.6 | -4.87
| N= 2-1, J=3/2-1/2, F= 2- 1 | 174.722 | 12.6 | -4.95
| N= 2-1, J=3/2-1/2, F= 1- 0 | 174.728 | 12.6 | -5.09
| N= 2-1, J=3/2-1/2, F= 1- 1 | 174.733 | 12.6 | -5.29
| N= 2-1, J=3/2-3/2, F= 2- 2 | 174.807 | 12.6 | -5.57
CCH | N= 3-2, J=7/2-5/2, F= 4- 3 | 262.004 | 25.1 | -4.28
| N= 3-2, J=7/2-5/2, F= 3- 2 | 262.006 | 25.1 | -4.29
| N= 3-2, J=5/2-3/2, F= 3- 2 | 262.065 | 25.2 | -4.31
| N= 3-2, J=5/2-3/2, F= 2- 1 | 262.067 | 25.2 | -4.35
| N= 3-2, J=5/2-3/2, F= 2- 2 | 262.079 | 25.2 | -5.22
| N= 3-2, J=5/2-5/2, F= 3- 3 | 262.209 | 25.2 | -5.40
CCH | N= 4-3, J=9/2-7/2, F= 5- 4 | 349.338 | 41.9 | -3.88
| N= 4-3, J=9/2-7/2, F= 4- 3 | 349.339 | 41.9 | -3.89
| N= 4-3, J=7/2-5/2, F= 4- 3 | 349.399 | 41.9 | -3.90
| N= 4-3, J=7/2-5/2, F= 3- 2 | 349.401 | 41.9 | -3.92
C4H | N = 12-11, J = 25/2 - 23/2, F = 12 - 11 | 114.183 | 35.6 | -4.44
| N = 12-11, J = 25/2 - 23/2, F = 13 - 12 | 114.183 | 35.6 | -4.44
| N = 12-11, J = 23/2 - 21/2, F = 11 - 10 | 114.221 | 35.6 | -4.44
| N = 12-11, J = 23/2 - 21/2, F = 12 - 11 | 114.221 | 35.6 | -4.44

Note. — $(a)$ Rest frequency; $(b)$ Upper level energy of the transition;
$(c)$ Einstein $A$ coefficient; All values were taken from the Cologne
Database for Molecular Spectroscopy (CDMS; https://cdms.astro.uni-koeln.de;
Müller et al., 2001, 2005)

##  Appendix B Transition properties

Table 3 shows the spectroscopic properties of transitions used in our PCA in
an ascending order in frequency.

Table 3: A list of molecular transitions and RRLs used in PCA

Transition | $\nu_{\rm rest}^{(a)}$ |  $E_{\rm up}$(b) | Transition | $\nu_{\rm rest}^{(a)}$ |  $E_{\rm up}$(b)
---|---|---|---|---|---
| (GHz) | (K) |  | (GHz) | (K)
CH3CCH($5_{K}-4_{K}$) | 85.457 | 12.3 | HOCO+($4_{0,4}-3_{0,3}$) | 85.531 | 10.26
CCS($7_{6}-6_{5}$)$\dagger$ | 86.181 | 23.35 | H13CN(1-0)$\dagger$ | 86.34 | 4.14
H13CO+(1-0) | 86.754 | 4.16 | SiO(2-1) | 86.847 | 6.25
HN13C(1-0) | 87.091 | 4.18 | CCH($1-0$) | 87.317 | 4.19
HNCO($4_{1,4}-3_{1,3}$) | 87.597 | 53.78 | HNCO($4_{0,4}-3_{0,3}$) | 87.925 | 10.55
H52$\beta$ | 88.406 |  | HCN(1-0) | 88.632 | 4.25
HCO+(1-0) | 89.189 | 4.28 | HOC+(1-0) | 89.487 | 4.29
C3H+(4-3) | 89.958 | 10.79 | HC5N(34-33) | 90.526 | 76.03
HNC(1-0) | 90.664 | 4.35 | HC3N(10-9) | 90.979 | 24.01
13CS(2-1) | 92.494 | 6.65 | N2H+(1-0) | 93.173 | 4.47
H58$\gamma$ | 93.776 |  | CCS($7_{8}-6_{7}$)$\dagger$ | 93.87 | 19.89
C34S(2-1) | 96.413 | 6.25 | CH3OH($2_{K}-1_{K}$) | 96.741 | 6.96
C33S(2-1) | 97.172 | 7.0 | OCS(8-7)$\dagger$ | 97.301 | 21.01
CH3OH($2_{1,1}-1_{1,0}$) | 97.583 | 21.56 |  34SO($2_{3}-1_{2}$)$\dagger$ | 97.715 | 9.09
CS(2-1)$\dagger$ | 97.981 | 7.1 | HC5N(37-36) | 98.513 | 89.83
H40$\alpha$ | 99.023 |  | HC3N(11-10) | 100.076 | 28.82
HC3N,v7=1(111-10-1) | 100.322 | 349.73 | CH3CCH(6K-5K) | 102.548 | 17.23
H2CS(30,3-20,2)$\dagger$ | 103.04 | 9.89 | H49$\beta$ | 105.302 |
CCS(89-78)$\dagger$ | 106.348 | 70.39 | H39$\alpha$ | 106.737 |
HOCO+(50,5-40,4) | 106.914 | 15.39 | CH3OH($3_{1,3}-4_{0,4}$) | 107.014 | 28.35
HC3N(12-11) | 109.174 | 34.06 | C18O(1-0)$\dagger$ | 109.782 | 5.27
HNCO($5_{0,5}-4_{0,4}$) | 109.906 | 15.82 |  13CO(1-0)$\dagger$ | 110.201 | 5.29
C17O(1-0)$\dagger$ | 112.359 | 5.39 | CN($1_{3/2}-0_{3/2}$) | 113.191 | 5.43
CN($1_{5/2}-0_{3/2}$) | 113.491 | 5.45 | C4H(12-11) | 114.221 | 35.64
CO(1-0) | 115.271 | 5.53 | HC3N(14-13) | 127.368 | 45.85
HOCO+($6_{0,6}-5_{0,5}$) | 128.295 | 21.55 | SiO(3-2) | 130.269 | 12.5
HNCO($6_{1,6}-5_{1,5}$) | 131.394 | 65.35 | HNCO($6_{0,6}-5_{0,5}$) | 131.886 | 22.15
HC3N(15-14) | 136.464 | 52.4 | SO($3_{4}-2_{3}$) | 138.179 | 15.86
13CS(3-2) | 138.739 | 13.32 | CH3OH($3_{1,3}-2_{1,2}$) | 143.866 | 28.35
CCS($11_{12}-10_{11}$)$\dagger$ | 144.245 | 43.94 | C34S(3-2) | 144.617 | 11.8
C33S(3-2) | 145.756 | 13.99 | HCNH+(2-1) | 148.221 | 10.67
HOCO+($7_{0,7}-6_{0,6}$) | 149.676 | 28.73 | HNCO(71,6-61,5) | 154.415 | 72.92
HC3N(17-16) | 154.657 | 66.8 | CH3OH($2_{1,2}-3_{0,3}$) | 156.602 | 21.44
C3H+(7-6) | 157.419 | 30.22 | HC3N(18-17) | 163.753 | 74.66
H${}_{2}^{34}$S($1_{1,0}-1_{1,0,1}$)$\dagger$ | 167.911 | 27.83 | H2S($1_{1,0}-1_{0,1}$)$\dagger$ | 168.763 | 27.88
H2CS($5_{1,5}-4_{1,4}$)$\dagger$ | 169.114 | 37.52 | CH3OH($3_{2,2}-2_{1,1}$) | 170.061 | 36.17
HCS+(4-3) | 170.692 | 20.48 | CH3CCH($10_{K}-9_{K}$) | 170.906 | 45.11
H13CN(2-1)$\dagger$ | 172.678 | 12.43 | HC3N(19-18) | 172.849 | 82.96
H13CO+(2-1)$\dagger$ | 173.507 | 12.49 | HN13C(2-1) | 174.179 | 12.54
CCH($2-1$) | 174.663 | 12.57 | HNCO($8_{0,8}-7_{0,7}$) | 175.844 | 37.98
HCN(2-1) | 177.261 | 12.76 | HCO+(2-1) | 178.375 | 12.84
SO(45-34)$\dagger$ | 178.605 | 24.43 | HOC+(2-1) | 178.972 | 12.88
HNC(2-1) | 181.325 | 13.05 | CH3CCH($11_{K}-10_{K}$) | 187.994 | 54.13
C34S(4-3) | 192.818 | 18.98 | CH3OH($4_{1,3}-3_{1,2}$) | 195.147 | 37.95
CS(4-3)$\dagger$ | 195.954 | 23.5 | HNCO($9_{0,9}-8_{0,8}$) | 197.821 | 47.47
CH2NH($3_{1,2}-2_{1,1}$) | 199.823 | 27.07 | HC3N(22-21) | 200.135 | 110.46
CH3CCH($12_{K}-11_{K}$) | 205.081 | 63.98 | c-C3H2($3_{3,0}-2_{2,1}$) | 216.279 | 13.53
H2S($2_{2,0}-2_{1,1}$)$\dagger$ | 216.71 | 83.98 | SiO(5-4) | 217.105 | 31.26
C18O(2-1) | 219.56 | 15.81 |  13CO(2-1) | 220.399 | 15.87
C17O(2-1)$\dagger$ | 224.714 | 16.18 | CN(2-1) | 226.66 | 16.31
HC3N(25-24) | 227.419 | 141.89 | CO(2-1) | 230.538 | 16.6
HC3N(26-25) | 236.513 | 153.25 | C33S(5-4) | 242.914 | 34.97
CS(5-4) | 244.936 | 35.3 | NO($3_{1,3}$-2-1,2) | 250.437 | 19.23
NO($3_{-1,3}-2_{1,2}$) | 250.796 | 19.28 |  34SO($6_{7}-5_{6}$)$\dagger$ | 256.878 | 46.71
H13CN(3-2)$\dagger$ | 259.012 | 24.86 | H13CO+(3-2) | 260.255 | 24.98
SiO(6-5) | 260.518 | 43.76 | CCH($3-2$) | 262.004 | 25.15
HCN(3-2) | 265.886 | 25.52 | CH2NH($4_{1,3}-3_{1,2}$) | 266.27 | 39.84
CH3OH($5_{2,4}-4_{1,3}$) | 266.838 | 57.07 | HCO+(3-2) | 267.558 | 25.68
HOC+(3-2) | 268.451 | 25.77 | HNC(3-2) | 271.981 | 26.11
HC3N(30-29) | 272.885 | 203.01 |  13CS(6-5) | 277.455 | 46.61
H2CS($8_{1,7}-7_{1,6}$)$\dagger$ | 278.888 | 73.41 | N2H+(3-2) | 279.512 | 26.83
CH3OH($6_{1,6}-5_{1,5}$) | 287.671 | 62.87 | C34S(6-5) | 289.209 | 38.19
CS(6-5) | 293.912 | 49.4 | H2S(3${}_{3,0}-3_{2,1}$)$\dagger$ | 300.506 | 168.9
CH3OH($1_{1,0}-1_{0,1}$) | 303.367 | 16.88 | CH3OH($3_{1,2}-3_{0,3}$) | 305.473 | 28.59
CH3OH($6_{1,5}-6_{0,6}$) | 311.853 | 63.71 | CH3OH($7_{1,6}-7_{0,7}$) | 314.86 | 80.09
CH3OH($6_{2,5}-5_{1,4}$) | 315.267 | 71.0 | H27$\alpha$ | 316.415 |
C18O(3-2) | 329.331 | 31.61 |  13CO(3-2) | 330.588 | 31.73
HC3N(37-36) | 336.52 | 306.91 | C17O(3-2)$\dagger$ | 337.061 | 32.35
CN(3-2) | 340.032 | 32.66 | CS(7-6) | 342.883 | 65.8
CO(3-2) | 345.796 | 33.19 | H13CO+(4-3) | 346.998 | 41.63
SiO(8-7) | 347.331 | 75.02 | CCH($4-3$) | 349.338 | 41.91
HNCO($16_{1,16}-15_{1,15}$) | 350.333 | 186.2 | HCN(4-3) | 354.505 | 42.53
HCO+(4-3) | 356.734 | 42.8 | HNC(4-3) | 362.63 | 43.51
H3O+($3_{2,0}-2_{2,1}$) | 364.797 | 97.06 | N2H+(4-3) | 372.672 | 44.71

Note. — $(a)$ Rest frequency; $(b)$ Upper level energy of the transition; All
values were taken from the Cologne Database for Molecular Spectroscopy (CDMS;
https://cdms.astro.uni-koeln.de; Müller et al., 2001, 2005). $\dagger$ Maps
that are used in PCA, but not shown in this paper because they are included in
the future publications.

##  Appendix C PC coefficients

The coefficients of all the 150 emissions for PC1 through PC 10 are in Table
4.

Table 4: PC coefficients for each species  Transition | PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7 | PC8 | PC9 | PC10
---|---|---|---|---|---|---|---|---|---|---
CO(1-0) | 0.069 | -0.127 | -0.129 | -0.046 | 0.163 | 0.029 | 0.161 | -0.001 | -0.116 | -0.065
CO(2-1) | 0.071 | -0.115 | -0.138 | -0.055 | 0.140 | -0.017 | 0.170 | -0.024 | -0.132 | -0.058
CO(3-2) | 0.075 | -0.097 | -0.139 | -0.057 | 0.101 | -0.060 | 0.125 | -0.011 | -0.087 | -0.035
13CO(1-0) | 0.071 | -0.131 | -0.102 | -0.030 | 0.151 | 0.133 | 0.081 | 0.064 | -0.001 | -0.110
13CO(2-1) | 0.082 | -0.080 | -0.102 | -0.068 | 0.079 | 0.053 | 0.034 | 0.038 | 0.015 | -0.059
13CO(3-2) | 0.087 | -0.031 | -0.093 | -0.071 | 0.019 | 0.015 | -0.024 | 0.046 | 0.034 | -0.022
C18O(1-0) | 0.076 | -0.123 | -0.078 | -0.028 | 0.120 | 0.145 | 0.055 | 0.065 | 0.052 | -0.079
C18O(2-1) | 0.086 | -0.064 | -0.071 | -0.075 | 0.054 | 0.058 | 0.014 | 0.043 | 0.066 | -0.015
C18O(3-2) | 0.090 | -0.005 | -0.057 | -0.067 | -0.005 | 0.031 | -0.044 | 0.047 | 0.056 | -0.013
C17O(1-0) | 0.075 | -0.120 | -0.073 | -0.012 | 0.099 | 0.188 | -0.008 | 0.090 | 0.102 | -0.074
C17O(2-1) | 0.087 | -0.057 | -0.060 | -0.068 | 0.036 | 0.067 | -0.027 | 0.056 | 0.072 | -0.057
C17O(3-2) | 0.090 | 0.017 | -0.030 | -0.058 | -0.009 | 0.019 | -0.073 | 0.051 | 0.082 | -0.017
HCN(1-0) | 0.080 | -0.101 | -0.093 | -0.073 | 0.068 | -0.063 | 0.089 | -0.109 | -0.093 | -0.001
HCN(2-1) | 0.082 | -0.073 | -0.111 | -0.045 | 0.063 | -0.085 | 0.125 | -0.053 | -0.094 | 0.037
HCN(3-2) | 0.088 | -0.019 | -0.078 | 0.000 | -0.035 | -0.112 | 0.016 | -0.032 | -0.033 | 0.046
HCN(4-3) | 0.089 | 0.012 | -0.060 | 0.025 | -0.054 | -0.098 | -0.021 | -0.023 | -0.011 | 0.037
HCO+(1-0) | 0.079 | -0.094 | -0.090 | -0.117 | 0.066 | -0.033 | 0.106 | -0.128 | -0.077 | -0.035
HCO+(2-1) | 0.082 | -0.063 | -0.113 | -0.091 | 0.066 | -0.056 | 0.138 | -0.063 | -0.062 | 0.011
HCO+(3-2) | 0.088 | -0.005 | -0.087 | -0.042 | -0.029 | -0.063 | 0.011 | -0.029 | -0.001 | 0.025
HCO+(4-3) | 0.089 | 0.027 | -0.072 | -0.007 | -0.043 | -0.015 | -0.041 | -0.015 | 0.031 | 0.029
CN($1_{3/2}-0_{3/2}$) | 0.083 | -0.053 | -0.101 | -0.128 | 0.061 | -0.047 | -0.008 | -0.073 | -0.025 | 0.038
CN($1_{5/2}-0_{3/2}$) | 0.082 | -0.050 | -0.116 | -0.120 | 0.045 | -0.067 | 0.058 | -0.099 | -0.057 | 0.023
CN(2-1) | 0.085 | 0.016 | -0.123 | -0.047 | -0.025 | -0.070 | -0.019 | -0.015 | -0.009 | 0.039
CN(3-2) | 0.083 | 0.065 | -0.120 | 0.016 | -0.052 | -0.052 | -0.053 | 0.023 | -0.003 | 0.019
HNC(1-0) | 0.086 | -0.043 | -0.098 | -0.063 | 0.041 | -0.008 | -0.014 | -0.067 | -0.015 | 0.040
HNC(2-1) | 0.087 | 0.032 | -0.083 | 0.004 | -0.069 | -0.040 | -0.075 | -0.068 | 0.073 | 0.053
HNC(3-2) | 0.087 | 0.045 | -0.085 | 0.027 | -0.022 | -0.034 | -0.062 | -0.069 | 0.009 | 0.027
HNC(4-3) | 0.086 | 0.072 | -0.065 | 0.048 | -0.007 | -0.023 | -0.089 | -0.094 | 0.007 | 0.015
N2H+(1-0) | 0.085 | -0.091 | -0.040 | 0.053 | 0.026 | -0.017 | -0.015 | -0.136 | -0.033 | -0.032
N2H+(3-2) | 0.088 | 0.015 | -0.055 | 0.101 | -0.031 | -0.043 | -0.029 | -0.162 | 0.017 | -0.051
N2H+(4-3) | 0.086 | 0.058 | -0.045 | 0.105 | -0.030 | -0.083 | -0.053 | -0.150 | 0.003 | -0.060
CS(2-1) | 0.084 | -0.092 | -0.062 | -0.062 | 0.067 | -0.045 | 0.046 | -0.013 | 0.016 | 0.003
CS(4-3) | 0.090 | -0.023 | -0.047 | -0.043 | 0.002 | -0.082 | 0.001 | 0.019 | 0.018 | 0.017
CS(5-4) | 0.091 | 0.014 | -0.034 | -0.029 | -0.012 | -0.075 | -0.011 | 0.032 | 0.024 | 0.016
CS(6-5) | 0.090 | 0.040 | -0.023 | -0.004 | -0.009 | -0.070 | -0.013 | 0.065 | 0.015 | 0.005
CS(7-6) | 0.088 | 0.064 | -0.019 | 0.028 | 0.008 | -0.068 | -0.027 | 0.092 | -0.008 | 0.014
CCH($1-0$) | 0.085 | -0.059 | -0.087 | -0.097 | 0.048 | 0.009 | -0.007 | -0.006 | 0.005 | 0.017
CCH($2-1$) | 0.086 | 0.051 | -0.057 | -0.032 | -0.101 | 0.003 | -0.118 | 0.002 | 0.022 | -0.009
CCH($3-2$) | 0.088 | 0.051 | -0.072 | -0.018 | -0.029 | -0.004 | -0.067 | 0.036 | -0.015 | 0.005
CCH($4-3$) | 0.084 | 0.080 | -0.088 | 0.021 | -0.039 | 0.006 | -0.056 | 0.015 | -0.006 | 0.007
H13CN(1-0) | 0.085 | -0.089 | -0.001 | -0.054 | 0.051 | -0.070 | -0.108 | -0.125 | -0.049 | 0.006
H13CN(2-1) | 0.091 | -0.007 | -0.023 | -0.019 | 0.020 | -0.079 | -0.069 | -0.047 | -0.045 | 0.021
H13CN(3-2) | 0.089 | 0.058 | 0.014 | 0.040 | 0.040 | -0.085 | -0.071 | -0.025 | 0.004 | 0.008
H13CO+(1-0) | 0.086 | -0.065 | -0.010 | -0.126 | 0.050 | 0.020 | -0.118 | -0.085 | -0.006 | -0.009
H13CO+(2-1) | 0.081 | 0.048 | 0.089 | -0.040 | 0.046 | 0.144 | 0.027 | 0.078 | 0.066 | -0.173
H13CO+(3-2) | 0.087 | 0.077 | -0.025 | -0.011 | -0.011 | 0.045 | -0.070 | -0.028 | -0.022 | -0.011
H13CO+(4-3) | 0.085 | 0.093 | 0.009 | 0.014 | 0.032 | 0.067 | -0.071 | -0.063 | -0.021 | -0.013
HN13C(1-0) | 0.089 | -0.036 | -0.003 | -0.049 | 0.022 | 0.009 | -0.171 | -0.100 | 0.042 | 0.005
HN13C(2-1) | 0.083 | 0.088 | 0.006 | 0.010 | 0.001 | 0.124 | -0.089 | -0.170 | -0.010 | -0.004
13CS(2-1) | 0.086 | -0.033 | 0.069 | -0.066 | -0.003 | -0.111 | -0.048 | 0.036 | 0.152 | -0.016
13CS(3-2) | 0.087 | 0.031 | 0.073 | -0.065 | 0.025 | -0.070 | -0.026 | 0.020 | 0.153 | -0.017
13CS(6-5) | 0.070 | 0.117 | 0.078 | 0.126 | 0.243 | -0.063 | 0.059 | 0.039 | 0.131 | 0.093
C34S(2-1) | 0.087 | -0.069 | -0.009 | -0.060 | 0.059 | -0.041 | 0.018 | 0.009 | 0.061 | 0.013
C34S(3-2) | 0.090 | -0.012 | -0.014 | -0.036 | 0.048 | -0.057 | 0.011 | 0.043 | 0.096 | 0.033
C34S(4-3) | 0.089 | 0.043 | 0.027 | -0.014 | 0.047 | -0.084 | -0.025 | 0.039 | 0.087 | 0.019
C34S(6-5) | 0.084 | 0.092 | 0.019 | 0.039 | 0.102 | -0.075 | -0.044 | 0.085 | 0.073 | 0.027
C33S(2-1) | 0.087 | -0.015 | 0.059 | -0.093 | -0.020 | -0.086 | -0.046 | 0.024 | 0.171 | -0.053
C33S(3-2) | 0.086 | 0.057 | 0.078 | -0.021 | 0.061 | -0.045 | -0.015 | -0.011 | 0.079 | 0.003
C33S(5-4) | 0.081 | 0.084 | 0.085 | 0.038 | 0.105 | 0.021 | -0.024 | -0.002 | -0.036 | 0.013
HOC+(1-0) | 0.069 | -0.120 | -0.075 | -0.151 | 0.138 | 0.070 | 0.174 | 0.184 | 0.206 | 0.197
HOC+(2-1) | 0.081 | 0.048 | 0.015 | -0.120 | -0.065 | 0.205 | -0.103 | 0.293 | -0.075 | 0.023
HOC+(3-2) | 0.084 | 0.073 | -0.068 | -0.069 | -0.015 | 0.049 | -0.182 | 0.008 | -0.109 | 0.065
SiO(2-1) | 0.081 | -0.112 | 0.034 | 0.012 | 0.025 | -0.058 | -0.117 | -0.118 | -0.102 | 0.060
SiO(3-2) | 0.087 | -0.070 | 0.033 | 0.039 | -0.004 | -0.057 | -0.063 | -0.071 | -0.087 | 0.060
SiO(5-4) | 0.090 | 0.013 | 0.044 | 0.039 | -0.026 | -0.095 | -0.061 | 0.046 | -0.072 | 0.024
SiO(6-5) | 0.088 | 0.055 | 0.046 | 0.050 | 0.010 | -0.080 | -0.036 | 0.079 | -0.084 | 0.020
SiO(8-7) | 0.072 | 0.073 | 0.166 | 0.034 | 0.186 | -0.072 | 0.034 | 0.083 | 0.007 | 0.018
CH3OH($2_{K}-1_{K}$) | 0.048 | -0.211 | -0.010 | 0.191 | 0.053 | 0.032 | 0.040 | 0.023 | -0.054 | 0.013
CH3OH($2_{1,1}-1_{1,0}$) | 0.060 | -0.163 | 0.090 | 0.178 | -0.120 | -0.103 | 0.049 | -0.066 | 0.058 | 0.057
CH3OH($3_{1,3}-4_{0,4}$) | 0.078 | 0.030 | 0.155 | -0.106 | -0.023 | 0.166 | 0.051 | -0.004 | -0.139 | -0.019
CH3OH($3_{1,3}-2_{1,2}$) | 0.065 | -0.151 | 0.093 | 0.140 | -0.125 | -0.137 | 0.075 | -0.077 | 0.026 | 0.074
CH3OH($2_{1,2}-3_{0,3}$) | 0.084 | 0.010 | 0.082 | -0.065 | -0.073 | 0.084 | 0.115 | 0.119 | 0.114 | 0.133
CH3OH($3_{2,2}-2_{1,1}$) | 0.076 | -0.135 | 0.041 | 0.019 | -0.033 | 0.002 | 0.050 | 0.059 | 0.135 | 0.077
CH3OH($4_{1,3}-3_{1,2}$) | 0.079 | -0.076 | 0.110 | 0.023 | -0.139 | -0.139 | 0.115 | -0.125 | -0.042 | 0.014
CH3OH($5_{2,4}-4_{1,3}$) | 0.086 | -0.029 | 0.087 | -0.057 | -0.094 | -0.048 | 0.006 | 0.049 | -0.048 | -0.013
CH3OH($6_{1,6}-5_{1,5}$) | 0.077 | -0.041 | 0.135 | -0.028 | -0.176 | -0.106 | 0.220 | -0.197 | 0.030 | -0.001
CH3OH($1_{1,0}-1_{0,1}$) | 0.089 | -0.033 | 0.040 | -0.061 | 0.007 | 0.076 | -0.033 | 0.042 | -0.028 | 0.106
CH3OH($3_{1,2}-3_{0,3}$) | 0.089 | -0.012 | 0.049 | -0.057 | -0.045 | 0.084 | -0.049 | -0.011 | -0.045 | 0.067
CH3OH($6_{1,5}-6_{0,6}$) | 0.088 | 0.041 | 0.064 | -0.011 | -0.042 | 0.035 | -0.017 | -0.115 | -0.079 | 0.006
CH3OH($7_{1,6}-7_{0,7}$) | 0.088 | 0.041 | 0.056 | 0.001 | -0.028 | 0.019 | -0.027 | -0.071 | -0.145 | 0.025
CH3OH($6_{2,5}-5_{1,4}$) | 0.085 | -0.009 | 0.102 | -0.033 | -0.082 | 0.031 | 0.002 | -0.015 | -0.214 | 0.078
HNCO($4_{1,4}-3_{1,3}$) | 0.081 | 0.066 | -0.032 | 0.020 | -0.127 | 0.106 | 0.175 | -0.112 | 0.031 | -0.093
HNCO($4_{0,4}-3_{0,3}$) | 0.043 | -0.201 | -0.009 | 0.260 | 0.069 | 0.130 | -0.061 | 0.066 | -0.121 | 0.062
HNCO($5_{0,5}-4_{0,4}$) | 0.053 | -0.197 | -0.002 | 0.208 | 0.042 | 0.097 | -0.023 | 0.052 | -0.015 | 0.045
HNCO($6_{1,6}-5_{1,5}$) | 0.077 | -0.049 | 0.108 | -0.004 | -0.110 | 0.161 | 0.101 | -0.151 | 0.004 | -0.027
HNCO($6_{0,6}-5_{0,5}$) | 0.062 | -0.188 | 0.015 | 0.142 | 0.003 | 0.065 | 0.026 | 0.041 | 0.076 | 0.048
HNCO(71,6-61,5) | 0.078 | 0.056 | 0.124 | -0.027 | -0.009 | 0.181 | 0.213 | 0.047 | -0.091 | -0.018
HNCO($8_{0,8}-7_{0,7}$) | 0.082 | -0.087 | 0.058 | -0.049 | -0.136 | -0.003 | 0.010 | 0.056 | 0.217 | 0.032
HNCO($9_{0,9}-8_{0,8}$) | 0.082 | -0.048 | 0.099 | -0.105 | -0.141 | 0.014 | 0.093 | 0.050 | 0.142 | -0.019
HNCO($16_{1,16}-15_{1,15}$) | 0.080 | 0.109 | 0.036 | 0.059 | 0.061 | 0.163 | 0.058 | -0.166 | -0.067 | 0.023
HOCO+($4_{0,4}-3_{0,3}$) | 0.042 | -0.196 | 0.038 | 0.276 | 0.010 | 0.123 | -0.157 | 0.069 | -0.117 | 0.050
HOCO+(50,5-40,4) | 0.050 | -0.193 | 0.023 | 0.237 | -0.018 | 0.076 | -0.097 | 0.045 | -0.023 | -0.032
HOCO+($6_{0,6}-5_{0,5}$) | 0.076 | -0.116 | 0.027 | 0.124 | -0.108 | 0.002 | -0.026 | -0.000 | 0.173 | -0.036
HOCO+($7_{0,7}-6_{0,6}$) | 0.074 | -0.056 | 0.142 | -0.107 | -0.167 | -0.008 | 0.111 | -0.035 | 0.226 | -0.155
H40$\alpha$ | 0.068 | 0.113 | -0.171 | 0.115 | -0.096 | 0.015 | 0.053 | 0.073 | 0.009 | 0.040
H39$\alpha$ | 0.070 | 0.107 | -0.165 | 0.103 | -0.099 | 0.011 | 0.060 | 0.062 | 0.020 | 0.034
H27$\alpha$ | 0.066 | 0.118 | -0.161 | 0.147 | -0.083 | 0.042 | 0.124 | 0.018 | 0.017 | 0.053
H52$\beta$ | 0.064 | 0.119 | -0.170 | 0.150 | -0.108 | 0.002 | 0.097 | 0.065 | -0.014 | -0.017
H49$\beta$ | 0.065 | 0.119 | -0.171 | 0.141 | -0.106 | 0.026 | 0.075 | 0.027 | -0.011 | 0.029
H58$\gamma$ | 0.057 | 0.116 | -0.192 | 0.168 | -0.139 | 0.039 | 0.166 | 0.031 | 0.025 | 0.002
CH3CCH($5_{K}-4_{K}$) | 0.086 | -0.072 | 0.001 | -0.059 | -0.005 | 0.172 | -0.052 | -0.018 | 0.072 | 0.021
CH3CCH(6K-5K) | 0.088 | -0.050 | -0.015 | -0.061 | -0.040 | 0.146 | -0.035 | 0.001 | 0.096 | 0.044
CH3CCH($10_{K}-9_{K}$) | 0.089 | 0.010 | -0.009 | -0.087 | -0.050 | 0.097 | -0.039 | -0.040 | -0.019 | 0.039
CH3CCH($11_{K}-10_{K}$) | 0.083 | 0.070 | 0.067 | -0.034 | -0.052 | 0.232 | 0.078 | -0.035 | -0.126 | -0.035
CH3CCH($12_{K}-11_{K}$) | 0.089 | 0.033 | -0.001 | -0.062 | -0.061 | 0.097 | -0.020 | -0.027 | -0.022 | 0.011
SO($3_{4}-2_{3}$) | 0.090 | -0.036 | 0.008 | 0.030 | 0.012 | -0.041 | -0.057 | 0.039 | -0.035 | -0.012
SO(45-34) | 0.080 | -0.006 | 0.010 | 0.102 | 0.050 | -0.134 | 0.002 | 0.140 | -0.128 | -0.748
34SO($2_{3}-1_{2}$) | 0.086 | 0.061 | 0.050 | 0.013 | -0.075 | 0.049 | 0.114 | -0.021 | 0.040 | -0.106
34SO($6_{7}-5_{6}$) | 0.070 | 0.100 | 0.125 | 0.061 | 0.234 | 0.042 | -0.031 | -0.240 | -0.008 | 0.017
H2S($1_{1,0}-1_{0,1}$) | 0.088 | -0.010 | -0.040 | 0.003 | 0.149 | -0.090 | -0.057 | -0.030 | 0.013 | -0.010
H2S($2_{2,0}-2_{1,1}$) | 0.082 | 0.103 | 0.048 | 0.053 | 0.113 | 0.034 | -0.020 | -0.063 | 0.029 | 0.042
H2S(3${}_{3,0}-3_{2,1}$) | 0.075 | 0.118 | 0.058 | 0.099 | 0.194 | -0.022 | -0.044 | -0.047 | 0.066 | 0.065
OCS(8-7) | 0.063 | -0.162 | 0.071 | 0.187 | -0.024 | 0.079 | -0.046 | 0.003 | 0.097 | 0.012
H${}_{2}^{34}$S($1_{1,0}-1_{1,0,1}$) | 0.081 | 0.088 | 0.077 | 0.049 | 0.124 | 0.020 | 0.006 | -0.099 | 0.146 | -0.006
HC3N(10-9) | 0.087 | -0.081 | 0.033 | 0.014 | -0.040 | -0.031 | 0.008 | -0.009 | -0.027 | 0.032
HC3N(11-10) | 0.088 | -0.066 | 0.032 | 0.006 | -0.032 | -0.035 | 0.004 | -0.014 | -0.018 | 0.017
HC3N(12-11) | 0.089 | -0.053 | 0.042 | 0.003 | -0.042 | -0.054 | 0.016 | 0.001 | -0.031 | 0.015
HC3N(14-13) | 0.090 | -0.022 | 0.046 | -0.002 | -0.040 | -0.046 | 0.006 | 0.005 | -0.049 | 0.017
HC3N(15-14) | 0.090 | -0.018 | 0.043 | -0.000 | -0.037 | -0.067 | 0.013 | 0.035 | -0.090 | 0.002
HC3N(17-16) | 0.090 | 0.009 | 0.038 | 0.011 | -0.033 | -0.070 | -0.003 | 0.083 | -0.077 | 0.033
HC3N(18-17) | 0.090 | 0.001 | 0.033 | -0.010 | -0.009 | -0.086 | 0.006 | 0.078 | -0.057 | 0.025
HC3N(19-18) | 0.089 | 0.030 | 0.041 | -0.012 | -0.028 | -0.094 | -0.017 | 0.094 | -0.089 | 0.004
HC3N(22-21) | 0.086 | 0.068 | 0.072 | 0.018 | 0.018 | -0.056 | 0.031 | 0.127 | -0.072 | -0.004
HC3N(25-24) | 0.087 | 0.062 | 0.051 | 0.020 | 0.063 | -0.087 | 0.021 | 0.100 | 0.005 | 0.042
HC3N(26-25) | 0.086 | 0.076 | 0.056 | 0.034 | 0.057 | -0.078 | 0.014 | 0.104 | -0.050 | -0.003
HC3N(30-29) | 0.084 | 0.093 | 0.038 | 0.048 | 0.063 | -0.027 | 0.015 | 0.110 | 0.006 | 0.003
HC3N(37-36) | 0.075 | 0.111 | 0.090 | 0.096 | 0.192 | 0.013 | 0.060 | 0.023 | 0.067 | 0.050
HC3N,v7=1(111-10-1) | 0.068 | 0.094 | 0.149 | 0.068 | 0.246 | 0.019 | 0.052 | -0.061 | 0.073 | 0.053
H3O+($3_{2,0}-2_{2,1}$) | 0.087 | 0.046 | -0.077 | -0.023 | 0.025 | 0.004 | -0.091 | 0.065 | -0.004 | 0.090
CCS($7_{6}-6_{5}$) | 0.086 | 0.020 | -0.041 | -0.048 | -0.081 | -0.113 | -0.090 | 0.125 | 0.080 | -0.025
CCS($7_{8}-6_{7}$) | 0.089 | 0.001 | 0.032 | 0.004 | -0.064 | 0.016 | -0.084 | 0.016 | 0.031 | -0.094
CCS(89-78) | 0.088 | -0.000 | 0.045 | -0.017 | -0.036 | -0.051 | -0.041 | -0.024 | 0.072 | -0.075
CCS($11_{12}-10_{11}$) | 0.081 | 0.101 | 0.054 | 0.062 | 0.069 | 0.046 | 0.059 | 0.048 | -0.035 | -0.076
H2CS(30,3-20,2) | 0.081 | -0.085 | 0.076 | -0.001 | -0.093 | -0.076 | -0.110 | 0.018 | 0.046 | 0.010
H2CS($5_{1,5}-4_{1,4}$) | 0.078 | -0.018 | 0.130 | 0.033 | 0.013 | -0.169 | 0.025 | 0.142 | -0.028 | -0.118
H2CS($8_{1,7}-7_{1,6}$) | 0.080 | 0.087 | 0.094 | 0.060 | 0.098 | 0.011 | 0.142 | -0.011 | -0.014 | 0.044
HC5N(34-33) | 0.070 | -0.022 | 0.176 | -0.066 | -0.085 | -0.106 | 0.151 | 0.260 | -0.335 | 0.255
HC5N(37-36) | 0.082 | 0.032 | 0.140 | -0.030 | -0.047 | 0.019 | 0.166 | -0.034 | -0.014 | -0.060
C3H+(4-3) | 0.084 | 0.070 | -0.017 | -0.064 | -0.016 | 0.107 | -0.212 | 0.034 | -0.053 | 0.018
C3H+(7-6) | 0.065 | 0.132 | -0.108 | 0.192 | 0.007 | 0.012 | 0.168 | -0.082 | 0.126 | 0.034
CH2NH($3_{1,2}-2_{1,1}$) | 0.085 | 0.047 | 0.054 | -0.064 | -0.104 | 0.150 | 0.050 | -0.152 | -0.090 | -0.029
CH2NH($4_{1,3}-3_{1,2}$) | 0.085 | 0.069 | 0.090 | 0.006 | 0.029 | 0.053 | 0.103 | -0.007 | -0.001 | 0.013
C4H(12-11) | 0.085 | 0.071 | 0.019 | -0.045 | -0.042 | 0.132 | -0.040 | 0.106 | -0.117 | -0.055
HCNH+(2-1) | 0.080 | -0.054 | 0.088 | 0.047 | -0.067 | 0.130 | -0.104 | -0.038 | 0.085 | -0.043
HCS+(4-3) | 0.086 | 0.033 | 0.063 | -0.068 | -0.050 | -0.042 | -0.059 | -0.025 | -0.046 | -0.083
c-C3H2($3_{3,0}-2_{2,1}$) | 0.088 | 0.057 | -0.035 | -0.039 | -0.061 | 0.032 | -0.096 | -0.086 | 0.005 | -0.016
NO($3_{1,3}$-2-1,2) | 0.091 | -0.011 | -0.015 | 0.048 | -0.030 | 0.008 | -0.091 | 0.024 | -0.053 | -0.043
NO($3_{-1,3}-2_{1,2}$) | 0.091 | 0.004 | -0.012 | 0.043 | -0.018 | 0.026 | -0.102 | 0.048 | -0.047 | -0.037
3mm continuum | 0.068 | 0.101 | -0.182 | 0.101 | -0.096 | -0.008 | 0.096 | 0.112 | 0.026 | 0.036
0.8mm continuum | 0.085 | 0.067 | -0.078 | 0.042 | -0.065 | 0.101 | -0.039 | -0.036 | 0.040 | -0.024

##  Appendix D Estimation of errors in PCA

To estimate errors, we tried the so-called bootstrapping method. We first
created the rms error map, then multiplied it with a random number
distribution that has the same size as the number of pixels, mean of 0 and
standard distribution of 1. We added this noise map to the original image, and
ran PCA. We examined how results changed for 1000 different sets of random
number distributions for each transition tested. We only used selected
transitions for bootstrapping: HCN(1-0) and 13CO(1-0) for high-intensity
transitions, H39$\alpha$ for intermediate, HC3N v7=1 (111-10-1),
CH2NH(41,3–31,2), and OCS(8-7) for low-intensity transitions. Results are
shown in Table 5. For high-intensity transitions, addition of noise did not
change results in any significant way, and results are not shown. Variation
among 1000 different noise map added is overall small. However, we also note
that, for a compact transition such as HC3N v7=1 (111-10-1), an act of adding
noise itself changes the PCs more than the variation among different random
noise sets. Nonetheless, qualitative results do not change even for these
transitions.

Table 5: Bootstrapping results

PC | Original | Mean(Bootstrap) | Standard deviation(Bootstrap)
---|---|---|---
H39$\alpha$
PC1 | 0.070 | 0.070 | 0.000
PC2 | 0.107 | 0.107 | 0.001
PC3 | -0.165 | -0.165 | 0.001
PC4 | 0.103 | 0.102 | 0.001
PC5 | -0.099 | -0.099 | 0.002
PC6 | 0.011 | 0.011 | 0.003
PC7 | 0.060 | 0.060 | 0.003
PC8 | 0.062 | 0.063 | 0.004
PC9 | 0.020 | 0.021 | 0.004
PC10 | 0.034 | 0.034 | 0.005
HC3N v7=1 (111-10-1)
PC1 | 0.068 | 0.066 | 0.000
PC2 | 0.094 | 0.092 | 0.002
PC3 | 0.149 | 0.146 | 0.003
PC4 | 0.068 | 0.067 | 0.005
PC5 | 0.246 | 0.247 | 0.005
PC6 | 0.019 | 0.019 | 0.008
PC7 | 0.052 | 0.057 | 0.012
PC8 | -0.061 | -0.070 | 0.016
PC9 | 0.073 | 0.082 | 0.014
PC10 | 0.053 | 0.072 | 0.029
CH2NH(41,3-31,2)
PC1 | 0.085 | 0.084 | 0.000
PC2 | 0.069 | 0.068 | 0.001
PC3 | 0.090 | 0.089 | 0.002
PC4 | 0.006 | 0.006 | 0.003
PC5 | 0.029 | 0.030 | 0.003
PC6 | 0.053 | 0.054 | 0.005
PC7 | 0.103 | 0.105 | 0.006
PC8 | -0.007 | -0.007 | 0.008
PC9 | -0.001 | -0.002 | 0.007
PC10 | 0.013 | 0.013 | 0.010
OCS(8-7)
PC1 | 0.063 | 0.063 | 0.000
PC2 | -0.162 | -0.162 | 0.000
PC3 | 0.071 | 0.071 | 0.001
PC4 | 0.187 | 0.187 | 0.001
PC5 | -0.024 | -0.023 | 0.001
PC6 | 0.079 | 0.079 | 0.002
PC7 | -0.046 | -0.045 | 0.003
PC8 | 0.003 | 0.002 | 0.004
PC9 | 0.097 | 0.099 | 0.003
PC10 | 0.012 | 0.011 | 0.005

Note. — Transitions with strong emission such as 13CO(1-0) and HCN(1-0) did
not show any significant difference between the original and bootstrapped
cases, with very little standard deviation from the errors. Therefore, these
transitions are not shown.

## References

  * Aladro et al. (2015) Aladro, R., Martín, S., Riquelme, D., et al. 2015, A&A, 579, A101, doi: 10.1051/0004-6361/201424918
  * Ando et al. (2017) Ando, R., Nakanishi, K., Kohno, K., et al. 2017, ApJ, 849, 81, doi: 10.3847/1538-4357/aa8fd4
  * 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
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Behrens et al. (2022) Behrens, E., Mangum, J. G., Holdship, J., et al. 2022, ApJ, 939, 119, doi: 10.3847/1538-4357/ac91ce
  * Bellocchi et al. (2023) Bellocchi, E., Martín-Pintado, J., Rico-Villas, F., Martín, S., & Jiménez-Sierra, I. 2023, MNRAS, 519, L68, doi: 10.1093/mnrasl/slac154
  * Bendo et al. (2015) Bendo, G. J., Beswick, R. J., D’Cruze, M. J., et al. 2015, MNRAS, 450, L80, doi: 10.1093/mnrasl/slv053
  * Böker et al. (2008) Böker, T., Falcón-Barroso, J., Schinnerer, E., Knapen, J. H., & Ryder, S. 2008, AJ, 135, 479, doi: 10.1088/0004-6256/135/2/479
  * Bolatto et al. (2013) Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013, Nature, 499, 450, doi: 10.1038/nature12351
  * Cernicharo et al. (2022) Cernicharo, J., Agúndez, M., Cabezas, C., et al. 2022, A&A, 657, L16, doi: 10.1051/0004-6361/202142992
  * Chidiac (2020) Chidiac, C. 2020, PhD thesis, Rheinische Friedrich Wilhelms University of Bonn, Germany
  * Churchwell et al. (1986) Churchwell, E., Wood, D., Myers, P. C., & Myers, R. V. 1986, ApJ, 305, 405, doi: 10.1086/164256
  * Costagliola et al. (2011) Costagliola, F., Aalto, S., Rodriguez, M. I., et al. 2011, A&A, 528, A30, doi: 10.1051/0004-6361/201015628
  * Costagliola et al. (2015) Costagliola, F., Sakamoto, K., Muller, S., et al. 2015, A&A, 582, A91, doi: 10.1051/0004-6361/201526256
  * Das et al. (2001) Das, M., Anantharamaiah, K. R., & Yun, M. S. 2001, ApJ, 549, 896, doi: 10.1086/319430
  * Fontani et al. (2021) Fontani, F., Colzi, L., Redaelli, E., Sipilä, O., & Caselli, P. 2021, A&A, 651, A94, doi: 10.1051/0004-6361/202140655
  * Fukushima & Yajima (2021) Fukushima, H., & Yajima, H. 2021, MNRAS, 506, 5512, doi: 10.1093/mnras/stab2099
  * Gao & Solomon (2004) Gao, Y., & Solomon, P. M. 2004, ApJ, 606, 271, doi: 10.1086/382999
  * García-Burillo et al. (2017) García-Burillo, S., Viti, S., Combes, F., et al. 2017, A&A, 608, A56, doi: 10.1051/0004-6361/201731862
  * Gorski et al. (2023) Gorski, M. D., Aalto, S., König, S., et al. 2023, A&A, 670, A70, doi: 10.1051/0004-6361/202245166
  * Gratier et al. (2017) Gratier, P., Bron, E., Gerin, M., et al. 2017, A&A, 599, A100, doi: 10.1051/0004-6361/201629847
  * Grudić et al. (2018) Grudić, M. Y., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2018, MNRAS, 475, 3511, doi: 10.1093/mnras/sty035
  * Haasler et al. (2022) Haasler, D., Rivilla, V. M., Martín, S., et al. 2022, A&A, 659, A158, doi: 10.1051/0004-6361/202142032
  * Harada et al. (2018) Harada, N., Sakamoto, K., Martín, S., et al. 2018, ApJ, 855, 49, doi: 10.3847/1538-4357/aaaa70
  * Harada et al. (2019) —. 2019, ApJ, 884, 100, doi: 10.3847/1538-4357/ab41ff
  * Harada et al. (2021) Harada, N., Martín, S., Mangum, J. G., et al. 2021, ApJ, 923, 24, doi: 10.3847/1538-4357/ac26b8
  * Harada et al. (2022) —. 2022, ApJ, 938, 80, doi: 10.3847/1538-4357/ac8dfc
  * Henkel et al. (2018) Henkel, C., Mühle, S., Bendo, G., et al. 2018, A&A, 615, A155, doi: 10.1051/0004-6361/201732174
  * Holdship et al. (2021) Holdship, J., Viti, S., Martín, S., et al. 2021, A&A, 654, A55, doi: 10.1051/0004-6361/202141233
  * Holdship et al. (2022) Holdship, J., Mangum, J. G., Viti, S., et al. 2022, ApJ, 931, 89, doi: 10.3847/1538-4357/ac6753
  * Huang et al. (2023) Huang, K. Y., Viti, S., Holdship, J., et al. 2023, arXiv e-prints, arXiv:2303.12685, doi: 10.48550/arXiv.2303.12685
  * Humire et al. (2022) Humire, P. K., Henkel, C., Hernández-Gómez, A., et al. 2022, A&A, 663, A33, doi: 10.1051/0004-6361/202243384