# Energizing Star Formation: The Cosmic Ray Ionization Rate in NGC 253 Derived
From ALCHEMI Measurements of H3O+ and SO

Jonathan Holdship Leiden Observatory, Leiden University, PO Box 9513, NL-2300
RA Leiden, The Netherlands  Department of Physics and Astronomy, University
College London, Gower Street, London WC1E 6BT  Jeffrey G. Mangum National
Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA
22903-2475, USA  Serena Viti Leiden Observatory, Leiden University, PO Box
9513, NL-2300 RA Leiden, The Netherlands  Department of Physics and Astronomy,
University College London, Gower Street, London WC1E 6BT  Erica Behrens
Department of Astronomy, University of Virginia, P. O. Box 400325, 530
McCormick Road, Charlottesville, VA 22904-4325  Nanase Harada National
Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
Department of Astronomy, School of Science, The Graduate University for
Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo, 181-1855 Japan
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  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
Sebastien Muller Department of Space, Earth and Environment, Chalmers
University of Technology, Onsala Space Observatory, SE-43992 Onsala, Sweden
Kunihiko Tanaka Department of Physics, Faculty of Science and Technology, Keio
University, 3-14-1 Hiyoshi, Yokohama, Kanagawa 223–8522 Japan  Kouichiro
Nakanishi National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka,
Tokyo 181-8588, Japan  Department of Astronomy, School of Science, The
Graduate University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka,
Tokyo, 181-1855 Japan  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  Yuki Yoshimura  Institute of Astronomy, Graduate School of Science, The
University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan  Rebeca
Aladro Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn,
Germany  Laura Colzi Centro de Astrobiología (CSIC-INTA), Ctra. de Ajalvir Km.
4, 28850, Torrejón de Ardoz, Madrid, Spain  INAF-Osservatorio Astrofisico di
Arcetri, Largo E. Fermi 5, I-50125, Florence, Italy 5  Kimberly L. Emig
National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA
22903-2475, 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  Yuri
Nishimura Institute of Astronomy, 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  Víctor M. Rivilla
Centro de Astrobiología (CSIC-INTA), Ctra. de Ajalvir Km. 4, 28850, Torrejón
de Ardoz, Madrid, Spain  INAF-Osservatorio Astrofisico di Arcetri, Largo E.
Fermi 5, I-50125, Florence, Italy 5  Paul P. van der Werf Leiden Observatory,
Leiden University, PO Box 9513, NL - 2300 RA Leiden, The Netherlands J.
Holdship<EMAIL_ADDRESS>
###### Abstract

The cosmic ray ionization rate (CRIR) is a key parameter in understanding the
physical and chemical processes in the interstellar medium. Cosmic rays are a
significant source of energy in star formation regions, which impacts the
physical and chemical processes which drive the formation of stars. Previous
studies of the circum-molecular zone (CMZ) of the starburst galaxy NGC 253
have found evidence for a high CRIR value; $10^{3}-10^{6}$ times the average
cosmic ray ionization rate within the Milky Way. This is a broad constraint
and one goal of this study is to determine this value with much higher
precision. We exploit ALMA observations towards the central molecular zone of
NGC 253 to measure the CRIR. We first demonstrate that the abundance ratio of
H3O+and SO is strongly sensitive to the CRIR . We then combine chemical and
radiative transfer models with nested sampling to infer the gas properties and
CRIR of several star-forming regions in NGC 253 from emission from their
transitions. We find that each of the four regions modelled has a CRIR in the
range $(1-80)\times 10^{-14}$ s-1 and that this result adequately fits the
abundances of other species that are believed to be sensitive to cosmic rays
including C2H, HCO+, HOC+, and CO. From shock and PDR/XDR models, we further
find that neither UV/X-ray driven nor shock dominated chemistry are a viable
single alternative as none of these processes can adequately fit the
abundances of all of these species.

starburst galaxies,Galaxies: ISM, Galaxies: active, Galaxies: abundances

††journal: The Astrophysical Journal

tablenum SIXtablenum

##  1 Introduction

Cosmic rays play an important role in the interstellar medium as a source of
heating and ionization. They drive chemistry by ionizing atoms; many gas-phase
reaction chains begin with ionization followed by a barrierless reaction
(Williams & Viti, 2013). Moreover, interactions with cosmic rays heat the gas
(Goldsmith, 2001), which can spur more complex chemistry and also affects the
dynamics of the system in question.

It is clear that measuring the cosmic ray ionization rate (CRIR) is vital to
understanding the chemical and dynamical evolution of molecular clouds.
Careful measurement of the CRIR is required in order to properly characterize
the energy budget within star formation regions. For this reason, there have
been many efforts to measure the CRIR in the Milky Way (e.g., Padovani et al.,
2009; Indriolo et al., 2015) and in extragalactic environments (e.g.,
González-Alfonso et al., 2013; Muller et al., 2016). These measurements
typically use ratios of OH+, H2O+, and H3O+ to infer the CRIR assuming
chemical equilibrium has been reached.

In this work, we attempt to infer the average CRIR for several regions in the
central molecular zone (CMZ) of the starburst galaxy NGC 253 using other
molecular line ratios. Due to its rich molecular emission, NGC 253 was
selected as the target of the ALMA Comprehensive High-Resolution Extragalactic
Molecular Inventory (ALCHEMI), an ALMA large program (Martín et al., 2021).
The goal of ALCHEMI is to produce the most complete molecular inventory of an
extragalactic object and use that to drive understanding of the CMZ.

The CMZ of NGC 253 contains several large ($\sim 50$ pc), dense
($\sim$${10}^{5}\text{\,}{\mathrm{cm}}^{-3}$), well-studied molecular clouds
(e.g., Sakamoto et al., 2011; Leroy et al., 2018) which are labelled in Fig.
1. Often labelled as giant molecular clouds (GMCs) due to their size, they
have much higher masses and have higher velocity dispersions (Leroy et al.,
2015) than typical GMCs in the Milky Way. We refer to these GMC-like
structures as GMCs for simplicity throughout this work. See Leroy et al.
(2015) for position information associated with these GMCs, and Mangum et al.
(2019) for the association of these GMCs with other identified dense gas
positions in NGC 253. It is for the brightest of these GMCs (3, 4, 5, 6 and
7), closest to the centre of the CMZ, that we intend to infer the CRIR.

Previous work has attempted to measure the CRIR in the GMCs of NGC 253. For
example, Holdship et al. (2021) found that a high C2H abundance in the GMCs
was likely caused by cosmic rays. However, an attempt to use chemical
modelling to measure the CRIR found that it could only be constrained in the
range
${10}^{-14}\text{\,}{\mathrm{s}}^{-1}{10}^{-11}\text{\,}{\mathrm{s}}^{-1}$ due
to the low sensitivity of the C2H abundance to the CRIR. Meanwhile, Harada et
al. (2021) found that HCO+ and HOC+ measurements indicated a cosmic ray
ionization rate >${10}^{-14}\text{\,}{\mathrm{s}}^{-1}$ confirming - but not
further constraining - the previous measurement. Despite these large
uncertainties, it is clear the CRIR in these regions is much higher than the
Galactic value which is typically measured to be in the range
$1\text{\times}{10}^{-17}\text{\,}{\mathrm{s}}^{-1}10\text{\times}{10}^{-17}\text{\,}{\mathrm{s}}^{-1}$
(Padovani et al., 2009; Indriolo et al., 2015). More sensitive tracers are
required to obtain robust and accurate measurements of the CRIR in star
forming regions. In this work, we consider H3O+ and SO since their abundance
ratio has previously been found to be strongly dependent on the CRIR (Bayet et
al., 2011). We first confirm this dependency using a large grid of chemical
models and then use emission from these species to infer the CRIR in the CMZ
of NGC 253.

In Sect. 2, the data reduction steps followed to extract line intensities from
the ALCHEMI image cubes are described. In Sect. 3 we present preliminary
modelling to justify our analysis and the approach used to infer the CRIR for
the targeted GMCs of NGC 253. In Sect. 4 we present the results of the
analysis which is discussed in Sect. 5. We present our summary in Sect. 6.

##  2 Observational Data

###  2.1 ALCHEMI Data

We make use of data acquired as part of the ALCHEMI ALMA large program. A full
description of the ALCHEMI observations and data reduction can be found in
Martín et al. (2021) but important details are given here. ALCHEMI is an
unbiased spectral survey of the central molecular zone (CMZ) of the starburst
galaxy NGC 253 covering ALMA bands 3 through 7 (84 to 373 GHz).

Figure 1: Sample H3O+ (top) and SO (bottom) integrated intensity (moment 0)
images toward NGC 253. For each image the green ellipse and black scale bar in
the lower-left and lower-right corners show the final imaged beam size (1.6
arcsec) and physical scale, respectively. Red numbers indicate the locations
of the dense molecular emission regions identified by Leroy et al. (2015,
Table 4). Star shaped markers locate the positions of the 2 cm radio continuum
emission peaks (Ulvestad et al., 1997), with a square indicating the position
of the strongest radio continuum peak identified by Turner et al. (1985). The
lower integrated intensity limit for each transition is set to 3$\sigma$ in
the integrated intensity (see Section 2.4). Overlain in contours is the
associated continuum emission distribution for each transition. Continuum
contours are in steps of 3, 6, 9, 12, 30, 120, 240, and 900 times the
respective continuum RMS, where the peak continuum intensity dictates the
number of these levels actually used for a given panel. The continuum RMS
values for the H3O+ and SO transitions shown are 1.0 and 0.15 mJy/beam,
respectively.

NGC 253 was observed toward a nominal phase center of
$\alpha=00^{h}47^{m}33.26^{s}$,
$\delta=-25^{\circ}17^{\prime}17.7^{\prime\prime}$ (ICRS). Observations were
configured to cover a common rectangular area of $50\arcsec\times 20\arcsec$
with a position angle of $65^{\circ}$ (East of North). All ALCHEMI image cubes
presented in this article were imaged to a common beam size of 1$\farcs$6\.
The maximum recoverable scale for the ALCHEMI observations is 15″. If the
distance of NGC 253 is taken to be 3.5 Mpc (Rekola et al., 2005), these
angular size scales correspond to linear size scales of 28 and 250 pc
respectively. The spectral range of the ALCHEMI data includes 37 SO and 2 H3O+
rotational transitions. These SO and H3O+ transitions are listed with their
frequencies and quantum numbers in Table 1. The recommended $1\sigma$ absolute
flux calibration uncertainty for measurements from the ALCHEMI survey is 15%
(Martín et al., 2021), a value which we adopt in the analysis presented in
this article. We have also extracted the continuum emission associated with
the ALCHEMI measurements listed in Table 1 using the continuum subtraction and
imaging process described in Martín et al. (2021).

Table 1: SOaaSO energy levels designated using the NJ notation. and H3O+ Transitions and Frequencies Target | Frequency / GHz | EU / K | Interloper | Overlap CorrectionbbOverlap correction key: GMCs (3,4,5,6,7), except if all are 0.0 or 1.0.
---|---|---|---|---
H3O+ $1_{1}-2_{1}$ | 307.192 | 79.5 | CH3OH $4_{1}-4_{0}$ | (0.25,0.23,0.21,0.21,0.24)
H3O+ $3_{2}-2_{2}$ | 364.797 | 139.8 | OCS $30-29$ | (1.0,0.86,0.88,0.83,1.0)
| . |  | HC3N v7=1 40-39 | (1.0,1.0,0.71,1.0,1.0)
H3O+ $3_{0}-2_{0}$ | 396.272 | 169.1 |  | 1.0
SO $2_{2}-1_{1}$ | 86.093 | 19.3 | HC15N $1-0$ | 1.0
SO $2_{3}-1_{2}$ | 99.299 | 9.22 | NH2CN $5_{1}5-4_{1}4$ | 1.0
SO $5_{4}-4_{4}$ | 100.029 | 38.6 | HC3N $11-10$ | 1.0
SO $3_{2}-2_{1}$ | 109.252 | 21.1 | HC3N $12-11$ | 1.0
SO $3_{3}-2_{2}$ | 129.138 | 25.5 |  | 1.0
SO $6_{5}-5_{5}$ | 136.634 | 50.7 |  | 1.0
SO $3_{4}-2_{3}$ | 138.178 | 15.9 |  | 1.0
SO $4_{3}-3_{2}$ | 158.971 | 28.7 |  | 1.0
SO $4_{4}-3_{3}$ | 172.181 | 33.7 | HC15N $2-1$ | 1.0
SO $4_{5}-3_{4}$ | 178.605 | 24.4 |  | 1.0
SO $5_{4}-4_{3}$ | 206.176 | 38.6 | CCS $16_{1}5-15_{1}4$ | 1.0
SO $8_{7}-7_{7}$ | 214.357 | 81.2 | H${}_{2}^{34}$S $2_{2}0-2_{1}1$ | 0.0
SO $5_{5}-4_{4}$ | 215.220 | 44.1 |  | 1.0
SO $5_{6}-4_{5}$ | 219.949 | 35.0 |  | 1.0
SO $2_{1}-1_{2}$ | 236.452 | 15.8 | HC3N $26-25$ | 0.0
SO $3_{2}-2_{3}$ | 246.404 | 21.1 | several | 0.0
SO $6_{5}-5_{4}$ | 251.825 | 50.7 | CH3OH $6_{3}3-6_{2}4$ | (1.0,0.98,0.82,0.97,1.0)
| . |  | CH3OH $4_{3}2-4_{2}3$ | (0.84,0.55,0.87,0.81,0.65)
SO $9_{8}-8_{8}$ | 254.573 | 99.7 |  | 1.0
SO $6_{6}-5_{5}$ | 258.256 | 56.5 | HC15N $3-2$ | 1.0
SO $6_{7}-5_{6}$ | 261.844 | 47.6 | CH3OH $2_{1}1-1_{0}1$ | (0.88,0.76,0.74,0.82,0.89)
| . |  | C2H $3_{4}-2_{3}$ | (1.0,1.0.0.91,1.0,1.0)
SO $4_{3}-3_{4}$ | 267.198 | 28.7 | HCN,v2=1 $3-2$ | 0.0
SO $1_{1}-0_{1}$ | 286.340 | 15.2 | H2C18O $4_{1}3-3_{1}2$ | 1.0
SO $5_{4}-4_{5}$ | 294.768 | 38.6 | NH2CHO $14_{2},13-13_{2},12$ | 0.0
SO $10_{9}-9_{9}$ | 295.356 | 120.2 | H${}_{2}^{34}$S $3_{3}0-3_{2}1$ | 0.0
SO $7_{6}-6_{5}$ | 296.550 | 64.9 | HCNH+ $4-3$ | (1.0,1.0,0.86,1.0,1.0)
SO $7_{7}-6_{6}$ | 301.286 | 71.0 |  | 1.0
SO $7_{8}-6_{7}$ | 304.078 | 61.1 | SiO $7-6$ | 1.0
| . |  | OCS $25-24$ | (1.0,1.0,0.87,1.0,1.0)
| . |  | CH3OH $2_{1}1-2_{0}2$ | (1.0,1.0,0.96,1.0,1.0)
SO $2_{2}-1_{2}$ | 309.502 | 19.3 |  | 1.0
SO $2_{1}-1_{0}$ | 329.385 | 15.8 | C18O $3-2$ | 0.0
SO $11_{10}-10_{10}$ | 336.554 | 142.8 | HC3N $37-36$ | 0.0
SO $3_{3}-2_{3}$ | 339.341 | 25.5 |  | 1.0
SO $8_{7}-7_{6}$ | 340.714 | 81.2 | HC18O+ $4-3$ | (0.92,0.87,0.93,0.85,1.0)
SO $8_{8}-7_{7}$ | 344.311 | 87.5 | HC15N $4-3$ | (1.0,0.99,0.94,0.96,1.0)
SO $3_{2}-1_{2}$ | 345.705 | 21.1 | CO $3-2$ | 0.0
SO $8_{9}-7_{8}$ | 346.528 | 78.8 | HC3N,v7=1 $38-37$ | (1.0,1.0,1.0,0.90,1.0)
| . |  | SO2 $19_{1},19-18_{0},18$ | (1.0,1.0,1.0,0.98,1.0)
SO $7_{6}-6_{7}$ | 361.351 | 64.9 |  | 1.0

###  2.2 Additional ALMA Archival Data

Only two transitions of H3O+ are contained in the ALCHEMI data. This may limit
our ability to constrain our fits and therefore we augmented the ALCHEMI data
analyzed with ALMA archival measurements of the H3O+ $3_{0}-2_{0}$ transition.
This transition has a rest frequency of 396.272 GHz and the data was taken
from ALMA project 2016.1.01285.S (PI: Jesus Martín-Pintado). These data
include imaging of a single primary beam (PB) toward phase center position
$\alpha=00^{h}47^{m}33.134^{s}$,
$\delta=-25^{\circ}17^{\prime}19.68^{\prime\prime}$ (ICRS) with the ACA
($\theta_{PB}=38$″) and 12m Array ($\theta_{PB}=22$″). The observed phase
center position for these measurements is within 2″of the observed phase
center for the ALCHEMI measurements (Section 2.1). In order to directly
compare these H3O+ $3_{0}-2_{0}$ measurements with our other H3O+ and SO
measurements extracted from the ALCHEMI archive, we have imaged the H3O+
$3_{0}-2_{0}$ observations using the ALCHEMI imaging pipeline (Martín et al.,
2021). Starting with the ALMA-calibrated measurement set from the ALMA
archive, we produced continuum-subtracted H3O+ $3_{0}-2_{0}$ image cubes using
a robust parameter of 0.5. To match the spatial and spectral resolution of the
ALCHEMI measurements, the H3O+ $3_{0}-2_{0}$ images were imaged to a spatial
and spectral resolution of 1.6″ and 10 km s-1, respectively. The final
spectral channel RMS of these images is 2.6 mJy beam-1. The ALMA
recommendation for $1\sigma$ absolute flux calibration uncertainty for
measurements from ALMA Band 8 is 20% (ALMA Cycle 4 Proposer’s Guide, Section
A.9.2), a value which we adopt for the subsequent analysis of these H3O+
$3_{0}-2_{0}$ measurements. We have also extracted the continuum emission
associated with these measurements using the continuum subtraction and imaging
process described in Martín et al. (2021).

###  2.3 Spectral Line Overlap

Local thermodynamic equilibrium (LTE) modelling of the ALCHEMI data (Martín et
al., 2021) showed that the moderate line widths ($\Delta$v $\simeq 75$ km/s
FWHM) of the transitions detected towards NGC 253 cause many transitions to
overlap. As a result, measuring the emission above the noise in the vicinity
of each of our transitions will often result in total integrated intensities
that include flux from other transitions.

In order to estimate which of our SO and H3O+ transitions are expected to be
significantly blended with an interloper, we used the modeling results
described by Martín et al. (2021) as a guide, with inspection of the spectra
from each GMC. Where a target SO or H3O+ transition was found to be
potentially blended with the emission from another line, we estimate the
amount of overlap between the integrated intensities of these two transitions.
We then take the integrated intensity to be only the fraction of the emission
which comes from the target transition.

The overlap estimate is derived by performing multiple Gaussian fits to
spectra drawn through beam-sized areas centered on the position of each GMC.
The interloper overlap correction using these sample spectra is defined as

$\mathrm{Interloper~{}Correction}\equiv
1-\frac{\sum_{overlap}y_{interloper}}{\sum_{overlap}(y_{target}+y_{interloper})}$
(1)

where $y_{interloper}$ and $y_{target}$ are the intensities of the relative
frequency-constrained Gaussian fits to the target and interloper transitions,
respectively, within a spectral channel. The sum $\Sigma_{overlap}$ is taken
over all spectral channels with signal larger than the RMS noise in the image
cube under consideration. An example of the Gaussian fit analysis used to
derive the interloper correction factors defined by Equation 1 and listed in
Table 1 is shown in Fig. 2.

Figure 2: Observed spectrum from GMC 7 covering the H3O+ 11–21 transition in
grey with a multiple gaussian fit plotted in black. The colour traces show the
individual Gaussians used in the fit labelled by the transition which
corresponds to the central frequency. The derived interloper correction
factors, using Equation 1, are listed in Table 1.

As is apparent from the interloper correction factors listed in Table 1 and
the spectrum shown in Figure 2, the overlap between the H3O+ $1_{1}-2_{1}$ and
CH3OH $4_{1}-4_{0}$ transition is quite large, as the two transitions are
separated by just 26 MHz ($\sim 2.5$ spectral channels). As the H3O+
$1_{1}-2_{1}$ transition is one of only three H3O+ transitions available in
this study, it is important to properly assess the quality of our H3O+
spectral line extraction for this transition. The contribution from CH3OH to
the observed spectral profile was derived from LTE modelling using MADCUBA
(Martín et al., 2019) to both the CH3OH molecular emission across the entire
ALCHEMI frequency range as well as using only nearby surrounding (in frequency
and energies) CH3OH transitions. This LTE modeling procedure resulted in a
derived H3O+ to CH3OH spectral line peak intensity ratio of 0.2. Using the
procedure described above, with this derived Gaussian-peak relative intensity,
we derive the overlap corrections listed in Table 1.

###  2.4 Spectral Line Moment Extraction

Once all spectral interlopers were identified and their influence on the
spectrally-integrated emission from our target SO and H3O+ transitions
assessed, we proceeded to calculate the zeroth, first, and second moments of
all detected SO and H3O+ transitions using the
CubeLineMoment111https://github.com/keflavich/cube-line-extractor script
introduced for this same purpose by Mangum et al. (2019). CubeLineMoment uses
a series of spectral and spatial masks to extract integrated intensities for a
defined list of target spectral frequencies. As noted by Mangum et al. (2019),
the CubeLineMoment masking process uses a bright spectral line whose velocity
structure is representative of the emission over the galaxy as a ”tracer” of
the gas under study. In almost all of the SO and H3O+ transitions studied in
this analysis, we use the target transition itself as the tracer transition as
our target transitions are generally sufficiently intense. Final zeroth,
first, and second moment images were generated using a signal limit of three
times that of the spectral channel baseline rms for the respective transition
under study. Figure 1 shows representative samples of the H3O+ and SO
integrated intensity (zeroth moment) images resultant from this analysis. In
Appendix A we show the remaining H3O+ integrated intensity images and samples
of SO emission from the other ALMA receiver bands from our measurements (Bands
3, 5, 6, and 7).

As a result of this extraction, we obtain the integrated line intensity of
every transition in Table 1 for each GMC in the CMZ. This is done by taking
the average intensity over a beam sized region centred on the GMC centres
defined by Leroy et al. (2015). We then adopt uncertainties on these
intensities which combine the recommended 15% (for ALCHEMI measurements) or
20% (for H3O+ $3_{0}-2_{0}$ measurements) absolute calibration uncertainty and
the integrated noise added in quadrature. Finally, we adjust each integrated
intensity by the relevant interloper correction factor given in Table 1. We do
not adjust the uncertainty of a transition, choosing instead to retain the
total uncertainty on the blended emission. This effectively overestimates the
uncertainty on the intensity of the transition to account for our uncertainty
on the interloper correction factor. All integrated intensities, uncertainties
and other relevant transition properties are uploaded in a supplementary
table.

##  3 Modelling Approach

###  3.1 The Forward Model

To infer the cosmic ray ionization rate, we require a model which can take a
small number of free parameters, including the cosmic ray ionization rate, and
ultimately return line intensities for all the detected lines of H3O+ and SO.
To do this, we combine the gas-grain chemical model UCLCHEM (Holdship et al.,
2017) with RADEX (van der Tak et al., 2007) which is used via the Python
package SpectralRadex222spectralradex.readthedocs.io. We use recently
published collisional rates between H3O+ and p-H2 (Demes et al., 2021) as well
as SO collisional data (Lique et al., 2006), both taken from the
LAMDA333https://home.strw.leidenuniv.nl/~moldata/ database (Schöier et al.,
2005) .

In the chemical models, we assume the GMCs can be modelled as uniform clouds
with fixed physical conditions and a visual extinction that is sufficiently
high to make photo-processes negligible. This simple model is justified in the
next section. The simplified picture allows us to use a single point model as
the only depth dependent effect in UCLCHEM is the UV attenuation which becomes
unimportant if the majority of the gas considered is at high visual
extinction. This is important as even a 1D modelling approach would be
computationally unfeasible when combined with extensive sampling of the
necessary parameter space. One caveat of this model is that it is known that
the CRIR is attenuated by column density (Padovani et al., 2018). However,
that attenuation is weak for the column densities considered in this work and
so inferring an average value is reasonable.

The above chemical modelling requires three main inputs: the gas volume
density ($n_{H2}$), the gas kinetic temperature ($T_{kin}$) and the cosmic ray
ionizaton rate ($\zeta$) which will be given in units of
$\zeta_{0}$=$1.36\text{\times}{10}^{-17}\text{\,}{\mathrm{s}}^{-1}$ throughout
this work as that is the normalization factor used in the chemical network
(McElroy et al., 2013). In addition, we assume the initial elemental
abundances are the depleted values from Jenkins (2009) or that they are scaled
from those values by a constant metallicity factor ($Z$). The exception to
this is sulfur for which we treat the elemental abundance at $Z=1$ as a free
parameter (see Sect. 3.3). The model returns the equilibrium abundances of
H3O+ and SO. By adding the additional parameters of the total H2 column
density (NH2) of the GMC and the linewidth ($\Delta V$), we can generate line
intensities for all transitions of these two species by assuming a spherical
GMC and utilizing RADEX. RADEX requires the column density of each species
which is obtained by multiplying the fractional abundance from the chemical
model by the H2 column density. The collisional excitation rate file for
H3O+only contains rates for collisions with p-H2. For this reason we set the
p-H2 density in RADEX to $n_{H2}$ for H3O+calculations on the grounds that
assuming the o-H2 collisional rates are equal to the p-H2 rates must be a
better approximation than assuming no collisions with o-H2.

###  3.2 Model Justification

The GMCs in NGC253 are complex and many physical processes are at play. Thus,
two things must be demonstrated to justify the use of this modelling procedure
to infer the cosmic ray ionization rate. First, we must show the forward model
and, in particular, the H3O+ to SO ratio is sensitive to the cosmic ray
ionization rate as suggested by Bayet et al. (2011). Second, we must show that
other physical or chemical processes can be neglected. At a minimum, the outer
regions of any dense clouds that compose the single objects we observe at a
resolution of 28 pc will be UV irradiated. Further, on larger angular scales
($\gtrsim 15$″) low-velocity shocks have been shown to be a dominant heating
mechanism in the NGC 253 CMZ (Martín et al., 2006) and so shock chemistry may
be at work. Thus our simple, single point model with negligible influence from
an external UV field requires justification.

To address the model dependence on $\zeta$, a grid of UCLCHEM models were
generated in which the density, temperature and cosmic ray ionization rate of
a cloud were varied assuming a large visual extinction. Figure 3 shows the
steady state abundances of H3O+ and SO across the models of this grid as well
as their ratio. Both species show a clear and strong dependence on the cosmic
ray ionization rate with the H3O+ abundance increasing and the SO abundance
decreasing in response to an increasing $\zeta$.

In the case of H3O+ this is driven by the fact that H3O+ is formed through the
following chain of reactions,

$\displaystyle{\rm H_{2}+OH^{+}\longrightarrow H_{2}O^{+}+H}$ (2)
$\displaystyle{\rm H_{2}+H_{2}O^{+}\longrightarrow H_{3}O^{+}+H}$

where each reaction is the primary destruction route of the ionic reactant.
Therefore, the H3O+ abundance primarily depends on the OH+ abundance which is
primarily produced from both reactions between O+ and H2 and between H+ and
OH. In both cases, the ions are formed directly by cosmic rays and thus
increasing $\zeta$ generally increases the overall H3O+ production rate. This
ceases to be the case once $\zeta$ is such that the ionization fraction of the
gas is very high. At that point, it has been shown that increasing
dissociative reactions with electrons actually drive a decrease in the
H3O+abundance (Gerin et al., 2010). This is most noticable in the Tkin = 50 K
case in Figure 3 where the H3O+ abundance actually starts to decrease at the
highest $\zeta$.

SO is a much simpler case. It forms efficiently through neutral-neutral and
ion-neutral reactions in the gas-phase, obtaining relatively high abundances
($\sim$ ${10}^{-6}$) in most gas conditions. However, its primary destruction
routes are reactions with ions and cosmic rays, the majority being destroyed
in reactions with C+ and H+. Since the abundance of these reactants is
directly tied to the CRIR, SO is destroyed more efficiently as $\zeta$
increases. Fortunately, this appears to be even more efficient at low
temperatures, when the H3O+abundance is least sensitive to $\zeta$.

The opposite responses of these species to the CRIR results in their ratio
being highly sensitive to the value of $\zeta$. In fact, it varies by seven
orders of magnitude over the explored range of $\zeta$. Moreover, the
variation in the ratio due to $\zeta$ is much larger than the variation due to
temperature, so uncertainty in the gas temperature will not prevent us from
inferring the cosmic ray ionization rate. A final point to note about this
preliminary modelling is that the steady state abundances are reached quickly
- typically within ${10}^{5}$ years. This justifies the use of equilibrium
abundances for the GMCs.

Figure 3: Equilibrium H3O+ and SO abundances and their ratio as a function of
cosmic ray ionization rate given in units of
$\zeta_{0}$=$1.3\text{\times}{10}^{-17}\text{\,}{\mathrm{s}}^{-1}$. The
abundances are averaged over models with different densities in the range
${10}^{4}$ – ${10}^{7}\text{\,}{\mathrm{cm}}^{-3}$ and the colour of the lines
indicates the gas temperature.

To determine whether UV processes can be neglected, we use
UCL_PDR444https://github.io/UCL_PDR (Bell et al., 2006; Priestley et al.,
2017) to determine whether the assumption that H3O+ and SO primarily arise
from gas where the visual extinction is high is appropriate. This is a 1D
model which solves the equilibrium temperature and abundances for a semi-
infinite slab of gas. Fig. 4 demonstrates that for a broad range of conditions
which are reasonable for the GMCs under study, the vast majority of the H3O+
and SO columns in these PDR models arise from deeper within the model cloud.
Thus any H3O+ and SO emission from these objects should primarily trace the
higher column density regions where the visual extinction is greater than 5
magnitudes.

Figure 4: Abundance as a function of visual extinction for H3O+and SO from
UCL_PDR models. The majority of these species come from within the high AV
parts of the cloud. This is particularly true for SO, which is as much as
eight orders of magnitude higher in abundance once the $A_{V}$ is sufficiently
high to reduce photo-proccesses to zero than it is at the cloud edge.

Finally, the possibility that shocks throughout the gas affect the ratio of
H3O+and SO can also be addressed through preliminary modelling. To evaluate
this scenario, we generated a range of shock models using the C-shock
parameterization of Jiménez-Serra et al. (2008) which includes sputtering of
the ice mantles. We consider shock velocities between 5 and
$40\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ through gas with a preshock
density of ${10}^{5}$ and ${10}^{6}\text{\,}{\mathrm{cm}}^{-3}$. We find that
whilst the shock passage does tend to enhance the abundances of both species,
the ratio of the two species is relatively constant at H3O+/SO $\sim$
${10}^{-4}$ in shocked gas. Thus it will be important to test whether the
abundance ratios we obtain are close to this value but otherwise, we can
assume the ratio is not dominated by shocks.

In summary, our preliminary modelling shows that the ratio of H3O+ and SO is
strongly dependent on the ionization rate of the gas and comparatively weakly
dependent on the temperature. Furthermore, SO almost entirely arises from the
inner regions of clouds with high visual extinction and the H3O+ abundance is
also higher in these regions than close to the cloud edge so neglecting the UV
processes is justified. It is likely there are shocks present in the region
under study, but they tend to produce a constant H3O+ to SO ratio. If the
chemistry of these species is shock dominated, this can be checked a
posteriori by evaluating the H3O+ to SO ratio found by the inference.

###  3.3 Parameter Inference

With the model described in Section 3.1, we can use Bayesian inference to find
the probability distribution of the values of the model parameters given the
data that we observed (Section 2). From Bayes’ theorem, we know this
probability distribution is given by,

$p(\theta|d)=\frac{\mathcal{L}(d|\theta)p(\theta)}{p(d)}$ (3)

where $\theta$ represents the parameters and $d$ the data. The likelihood
($\mathcal{L}(d|\theta)$) is easy to formulate as our data have normally
distributed uncertainties, which allows us to use the standard Gaussian
likelihood,

$\mathcal{L}(d|\theta)=\exp\left(-\frac{1}{2}\sum_{i}{\frac{(d_{i}-M_{i})^{2}}{\sigma_{i}^{2}}}\right)$
(4)

where $d_{i}$ is the measured intensity of one transition, $\sigma_{i}$ is the
uncertainty of that measurement, and $M_{i}$ is the model spectral line
intensity for a given model parameter $\theta$.

The evidence ($p(d)$) is a constant for all $\theta$ for a given model and so
we neglect it in this work. To obtain the final probability distribution
$p(\theta|d)$, the probability distribution obtained from the numerator of
Equation 3 is simply normalized.

Finally, a prior distribution $p(\theta)$ must be chosen. The large amount of
previous work on NGC 253 provides a fantastic opportunity to set informed
priors on the density (Leroy et al., 2018), column density (Mangum et al.,
2019), and CRIR (Holdship et al., 2021; Harada et al., 2021). Despite this, we
chose in the first instance to use uninformative priors which are uniform
within limits based on those works. Finding that we adequately constrained our
parameters with these uninformative priors and the data at hand, we did not
progress to using informed priors. The ranges of the uniform priors are given
in Table 2.

A Bayesian inference approach gives us a simple way to deal with several
unknown parameters. In Table 2 we list the parameters of interest
($n_{H},T_{kin},N_{H2},\zeta$) as well as several nuisance parameters ($\Delta
V,S,Z,o/p$) which are unknown parameters that are not of interest to our work
but may affect our model fits. For example, it is well known that in Galactic
environments as little as 1% of sulfur is accounted for in dense gas
(Charnley, 1997). Thus it is unclear how much sulfur will be available for
gas-phase reactions in NGC 253 and so we leave the elemental abundance of
sulfur at a metallicity of 1 as a free parameter (S). Including nuisance
parameters in this way will account for the increased uncertainty in our
inferred values for the parameters of interest due to our lack of knowledge.

In addition to the sulfur depletion factor, we consider that the exact
metallicity (Z) of NGC 253 is unknown but is similar to the solar metallicity
(Marble et al., 2010) and so allow it to vary within a small range. We also
require an ortho:para ratio of H3O+ which UCLCHEM does not provide. Given that
this ratio has limiting values of 2 at temperatures below
$50\text{\,}\mathrm{K}$ and 1 at temperatures above $100\text{\,}\mathrm{K}$,
we allow the ratio to vary freely in this range.

Table 2: Prior distributions used for each parameteraaS $\equiv$ elemental sulfur abundance at Z=1 with $S_{\odot}$ its solar value. Symbol | Name | Prior Type | Range
---|---|---|---
$n_{H}$ | Gas density | log-uniform |  ${10}^{4}$–${10}^{7}\text{\,}{\mathrm{cm}}^{-3}$
$T_{kin}$ | Gas temperature | uniform | 50–300 K
$N_{H2}$ | H2 column density | log-uniform |  $3\text{\times}{10}^{22}$–${10}^{25}\text{\,}{\mathrm{cm}}^{-2}$
$\zeta$ | CRIR | log-uniform |  ${10}^{3}$–${10}^{7}$ $\zeta_{0}$
$\Delta V$ | Line width | uniform | 50–$150\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$
S | S abundance | log-uniform | 0.01 -1 S⊙
Z | Metallicity | uniform | 0.5–3 Z⊙
$o/p$ | H3O+ortho:para ratio | uniform | 1-2

Finally, in order to obtain the probability distribution $P(\theta|d)$, the
right hand side of Equation 3 was sampled using Ultranest (Buchner, 2021)
which is a nested sampling package. Given the relatively computationally
intensive nature of the forward model, the parameters utilized in UCLCHEM were
limited in precision by rounding. The temperature was rounded to the nearest
$1\text{\,}\mathrm{K}$ and other UCLCHEM input parameters were rounded to the
nearest 0.1 dex. This allowed model results to be saved and then if very
similar parameters were requested by the sampler, the abundances could be read
from file rather than repeatedly running UCLCHEM for very small changes in
parameter values. As a result the posterior distribution of $\theta$ is
limited in precision to the same degree.

##  4 Results

###  4.1 Model Spectral Line Intensities

Ultranest provides a set of samples from the posterior distribution weighted
so that points appear in proportion to how likely they are. We use this to
evaluate the goodness of fit of our model by generating the line fluxes using
all points in this sample and then comparing them to the measured fluxes. The
observed intensities with 1$\sigma$ errorbars are plotted in Figure 5
alongside the 16th to 83rd percentile range of the posterior sample fluxes.
This range was chosen because for a Gaussian distribution it is equivalent to
the 1$\sigma$ range.

Figure 5: Measured line intensities for the GMCs plotted as error bars which
cover the 1$\sigma$ uncertainty range. Shaded bars show the 16th to 83rd
percentile of fluxes predicted by our model when sampling from the posterior
distribution.

Ignoring clear outliers, the majority of transitions show an overlap between
the measurement uncertainties and the range of fluxes predicted by the model.
Some outlying transitions are to be expected considering that we are using a
highly simplified model of a homogeneous sphere to fit complex regions with
varying densities and temperatures. Points that do not overlap tend to be low
EU SO transitions which are underpredicted by the model. This is perhaps a
similar effect to that seen in C2H emission in the same GMCs (Holdship et al.,
2021) where low EU emission was clearly excited by an additional gas
component. However, since these deficiencies are small, we consider that the
majority of the gas emitting these species is well characterized by our model.

The H3O+ transition at Eu = $79.5\text{\,}\mathrm{K}$ ($1_{1}-2_{1}$ near 307
GHz) is also often underfit. However, this is the transition which suffers
from a large degree of overlap with a much stronger CH3OH line as shown in
Figure 2. It is possible that the overlap correction factors are too large and
we are assigning CH3OH flux to the H3O+ line. However, the H3O+ transition at
Eu = $129.8\text{\,}\mathrm{K}$ ($3_{2}-2_{2}$ near 365 GHz) is also much
stronger than our models predict. No potential interloper is strong enough to
explain the additional flux, but it is possible non-thermal excitation effects
are contributing to the excitation of this transition. Under certain
conditions, this transition can be infra-red pumped (Phillips et al., 1992;
Martín et al., 2021) which would explain the additional flux.

We checked that these outliers do not bias our results by fitting the GMCs
excluding the 307 and 365 GHz transitions. The H3Op abundance in these models
is then constrained by just the 396 GHz transition, which is not expected to
be radiatively pumped (Phillips et al., 1992). We find that we obtain similar
results, indicating these outliers are not strongly biasing our fit and it is
primarily the fit to SO that is driving our results.

It is possible that the chemistry creates a complicating factor which prevents
us finding adequate column densities to fit the data or to simultaneously fit
both species. We rule this out by fitting the data with RADEX only, treating
the SO and H3O+ column densities as free parameters alongside the gas
temperature and density. We find the same best fit fluxes are obtained,
indicating that no simple RADEX model can fit these data and it is not the
case that we simply cannot obtain the required abundances in our chemical
model.

Overall, whilst a better radiative transfer model, or perhaps simply a
multiple gas component fit to the data, would likely give better results, we
consider that our simple model fits are sufficient. Moreover, as discussed in
Sect. 4.3, the inference of the physical parameters allows us to validate our
results using previous measurements from each GMC modelled.

###  4.2 Posterior Distributions

The specific values of each posterior distribution differ from GMC to GMC but
the general trends are similar. We discuss here the corner plot for GMC 4
shown in Fig. 6 and include the corner plots for the other GMCs modelled in
Appendix B for completeness.

Figure 6: Corner plot showing the marginalized posterior probability
distribution of every parameter and the joint distributions of all parameter
pairs for GMC 4. Note the ortho:para ratio of H3O+is given as the fraction of
the H3O+in ortho form.

From the 1D marginalized posterior distributions, it is clear that we
constrain three of our four parameters of interest. The gas density, the H2
column density, and the cosmic ray ionization rate are well constrained whilst
the temperature is not. Each of the three constrained parameters has a strong
peak indicating a most likely value though the H2 column density peak becomes
flat below a certain limit. This is due to the fact that below a certain value
all column densities produce intensities well below the rms noise in our
measurements, and thus give identically poor fits.

The 2D joint posteriors show some interesting degeneracies. In particular, the
gas density and cosmic ray ionization rate shows a log-log relationship which
arises from a degeneracy in our chemical model. In UCLCHEM, the abundance can
remain unchanged for increasing CRIR if the density also increases. However,
since the density is also a parameter of RADEX, it is constrained by radiative
transfer considerations. Thus, the range of acceptable densities is greatly
limited allowing us to determine $\zeta$ despite the chemical degeneracy.
Thus, without an accurate measurement of the density, we would not be able to
estimate the cosmic ray ionization rate.

Finally, there are the four nuisance parameters we included in our inference.
The elemental sulfur abundance is constrained but the peak indicates that,
whatever the metallicity, the elemental sulfur abundance should have a similar
ratio to the other elements as found in the Sun. This means the sulfur cannot
be heavily depleted onto the grains in these regions. The other nuisance
parameters are unconstrained; the marginalized posteriors are very similar to
the priors. This is to the expected result for the line width in the case of
optically thin lines. Whilst constraining the metallicity or ortho:para ratio
of H3O+ would have been useful, the nuisance parameters were largely included
to account for their effect on the uncertainty of the cosmic ray ionization
rate and this is achieved regardless of the fact they are unconstrained.

###  4.3 The GMC Properties

In this section, we present the likely values of the three parameters of
interest which we have been able to constrain. Namely, the gas density, the H2
column density and the cosmic ray ionization rate. We present the most likely
value of each parameter in each region as well as a most likely interval that
contains 67% of the probability density, similar to a 1$\sigma$ uncertainty.
These are given in Table 3.

Table 3: Most likely values of each well constrained parameter with ranges containing 67% of the probability density. GMC |  $n_{H2}$ / ${10}^{5}\text{\,}{\mathrm{cm}}^{-3}$ |  $n_{H2}$ Range | NH2 / ${10}^{23}\text{\,}{\mathrm{cm}}^{-2}$ | NH2 Range |  $\zeta$ / ${10}^{4}$ $\zeta_{0}$**$\zeta_{0}$ = $1.36\text{\times}{10}^{-17}\text{\,}{\mathrm{s}}^{-1}$. |  $\zeta$ Range
---|---|---|---|---|---|---
3 | 1.3 | 0.4 - 4.9 | 0.8 | 0.4 - 2.1 | 0.2 | 0.1 - 1.2
4 | 2.0 | 0.9 - 5.9 | 3.7 | 1.0 - 9.3 | 0.8 | 0.2 - 3.2
5 | 2.2 | 1.0 - 7.2 | 3.0 | 1.1 - 6.5 | 1.1 | 0.3 - 4.5
6 | 2.6 | 1.1 - 9.8 | 2.5 | 0.7 - 7.5 | 0.5 | 0.2 - 2.6
7 | 0.3 | 0.1 - 1.2 | 5.3 | 1.1 - 21.7 | 0.8 | 0.2 - 5.9

The goal of this study is ultimately to infer the cosmic ray ionization rate
in the GMCs of NGC 253. It is therefore promising that we have constrained the
cosmic ray ionization rate in every case. Each of the five GMCs has a cosmic
ray ionization rate of $\zeta\sim$${10}^{4}$ $\zeta_{0}$. GMC 3 has a most
likely value that is a factor of 2-3 lower than the other GMCs and a
significant portion of the probability density is at low values. However, it
is consistent with the lower end of the range of likely values for the other
GMCs and so the difference cannot be said to be significant.

Figure 7: Violin plots showing the posterior probability distribution of
$\frac{\zeta}{n}$ in each GMC. The width of the blue area shows the relative
probability of a $\frac{\zeta}{n}$ value and the internal box plot marks the
quartiles.

The column densities we obtain are typically an order of magnitude lower than
those obtained from measurements of the dust (Mangum et al., 2019). However,
they are consistent with the H2 column densities derived from observations of
C18O (Harada et al., 2021). This agreement is important as RADEX relies on
species column densities. Thus the underlying abundances must be accurate if a
good fit is achieved with a reasonable H2 column density. This in turn gives
us confidence in our inferred CRIR as we know the SO and H3O+ abundances are
largely set by this parameter.

The gas densities are consistent with previous measurements. The range of
likely density values for each region overlap very well with those derived
from ALCHEMI observations of C2H (Holdship et al., 2021). They also fall well
within the range of ${10}^{5}$ \- ${10}^{6}\text{\,}{\mathrm{cm}}^{-3}$ found
in other ALCHEMI work (Harada et al., 2021) and elsewhere (Leroy et al.,
2018). Given how strongly correlated $\zeta$ is with this parameter, it is key
that these values are accurate. The chemical degeneracy is such that if the
density were a factor of 10 lower, our inferred $\zeta$ would be also be a
factor of ten smaller. Nevertheless, the good agreement with previously-
derived values of the gas number density from the literature should give
confidence in our derived CRIR.

However, we can also combine the density and CRIR to obtain a posterior
distribution on the value of $\frac{\zeta}{n}$. In Figure 7, we show these
distributions as violin plots. We actually constrain this value very well
indicating that a large part of the uncertainty we report on the CRIR is due
to the degeneracy with density. Future use of SO and H3O+as probes of the CRIR
could benefit from additional data to constrain the density or stronger priors
based on previous observations.

##  5 Discussion

###  5.1 The ALCHEMI view of the CRIR in NGC 253

We obtain a large cosmic ray ionization rate that is nevertheless consistent
with previous measurements. In previous ALCHEMI studies, Harada et al. (2021)
found that the cosmic ray ionization rate must be larger than ${10}^{3}$
$\zeta_{0}$ across the GMCs in the CMZ of NGC 253 for chemical models to
adequately explain the measured HCO+ and HOC+ emission. Holdship et al. (2021)
also find that a cosmic ray ionization rate of $\zeta$ = ${10}^{3}$-${10}^{6}$
$\zeta_{0}$ is required to explain high column densities of C2H. With H3O+and
SO, we are now able to constrain the CRIR to within an order of magnitude for
the first time, with all GMCs having a CRIR $\sim$${10}^{4}$ $\zeta_{0}$ or
${10}^{-13}\text{\,}{\mathrm{s}}^{-1}$.

It is useful to consider whether a consistent picture of the cosmic ray
ionization rate in NGC 253 is emerging from these studies. We therefore
validate our inferred parameter distributions using previously reported
abundances. To achieve this, we used the posterior samples provided by
Ultranest to run UCLCHEM and probability distributions on the abundances of
HCO+, HOC+, and C2H according to our parameter distributions. One could
include the line intensities of these species in the parameter inference but
validating in this way allows us to effectively check our fits on held out
data and evaluate the combination of H3O+and SO as a probe of the CRIR.

In analogy to a 1$\sigma$ interval, we define a likely range of model
abundances of HCO+, HOC+, and C2H by taking the 16th to 83rd percentile range
of their respective distributions. We then compare these ranges to the
observed abundances of those species as shown in Fig. 8. In almost every case,
the model and observed abundances overlap. The only exception is the C2H
abundance in GMC 4 but, statistically, one would not expect the 1$\sigma$
interval of every measurement to overlap with its true value so a small
difference between the model and one measurement is not surprising. The
overall agreement between the model and observations indicates that we have
found a consistent picture of cosmic ray driven chemistry in the GMCs of NGC
253.

As a further check, we also consider the abundance of H2, CO and e- and
include them in Fig. 8. This allows us to check how strongly ionized the gas
is. We find almost all H nuclei are in the form of H2 indicating that the gas
is still highly molecular under these conditions. Furthermore CO, we find that
the CO abundance is similar to the elemental C abundance. This is important
because Harada et al. (2021) reported a lower limit on the CRIR but limited it
from above by the observation that a sufficiently high CRIR will dissociate
CO, whilst observations of these regions indicate CO is abundant.

Figure 8: Abundance relative to total H nuclei of a selection of species for
each region. The model range is the range of abundances obtained when sampling
from the most likely portion of our parameter space using the constraints from
this work. We also show observed ranges for some species which are the
1$\sigma$ ranges from Holdship et al. (2021) (C2H) and Harada et al. 2021
(HCO+ and HOC+). The CO, H2 and e- abundances are shown to demonstrate that
despite the high CRIR, the gas is largely molecular.

###  5.2 Alternatives to a High CRIR

Holdship et al. (2021) found that either a high CRIR or a scenario where the
GMCs were sufficiently clumpy to allow PDR chemistry over a large column
density could reproduce the observed C2H abundance. Similarly, Harada et al.
(2021) found either PDR or cosmic ray dominated chemistry could reproduce
observed abundances of HOC+ and HCO+. This ambiguity is somewhat resolved by
the current study, as our H3O+and SO model fits greatly favour the high CRIR
scenario, largely due to the fact H3O+ is destroyed by UV but enhanced by
cosmic rays.

In Sect. 3.2, we showed that the SO and H3O+ abundances are very low in PDR
regions. However, when we use the posterior samples to generate model
abundances for the likely parameter range, we find H3O+ abundances between
${10}^{-10}$ and ${10}^{-9}$ which are too high for a PDR. Given that we have
a consistent picture of CR-driven chemistry which reproduces the abundances of
C2H, H3O+, SO, HOC+ and HCO+, a UV-dominant scenario that works for some
species but not others is disfavoured.

An alternative explanation for the high C2H abundance was the presence of
ubiquitous shocks in the GMCs. In Sect. 3.2, we found that in shock models,
the H3O+:SO abundance ratio was approximately ${10}^{-4}$ across a wide range
of physical parameters. Therefore, if our inferred ratio is much different to
this, it is unlikely shocks dominate the chemistry of H3O+ and SO. Using our
posterior samples, we find that the ratio of these species’ abundances varies
between 0.2 and 22.9 across all regions when we consider the most likely
parameter ranges given in Table 3. Since our inferred abundance ratio is
always at least three orders of magnitude larger than would be found in a
shock, it seems unlikely that it is in fact shock chemistry controlling the
abundance of these species.

One final possibility that has not yet been discussed in this work is that of
X-ray driven chemistry. Due to their similar ionizing effect and weak
attenuation with column density, X-rays can have very similar effects to
cosmic rays on the chemistry of a gas (Viti et al., 2014). Therefore, it is
possible some or all of the chemistry ascribed to cosmic rays in this work is
due to X-rays.

Whilst no X-ray flux was included in the PDR models discussed in Sect. 3.2,
UCL_PDR does treat X-rays if a flux is provided. We therefore run a model of
GMC 5 using our best fit gas density and an estimate of the largest X-ray flux
that can be motivated from observations of X-ray sources around the GMCs.
Lehmer et al. (2013) observed an X-ray source which they call Source B to be
in the brightest in NGC 253 with a luminosity of
$\sim$${10}^{39}\text{\,}\mathrm{erg}$. It is also very close to GMC 5, at a
distance of $\sim$$20\text{\,}\mathrm{pc}$. Combining these values, we obtain
a flux of
$0.05\text{\,}\mathrm{erg}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$.
By modelling a cloud that is subjected to this flux, we can observe the
maximum effect of X-rays on our GMC chemistry.

Figure 9: Lines show the modelled fractional abundance of SO and H3O+as a
function of visual extinction under an external X-ray flux of
$0.05\text{\,}\mathrm{erg}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$
and under no X-ray flux. The shaded regions show the measured 67% most likely
range of fractional abundances.

We plot the results in Fig. 9. We find that increasing the X-ray flux
increases the H3O+ abundance and decreases the SO abundance in a similar
manner to cosmic rays. However, the effect of an X-ray flux of this magnitude
is not sufficient to bring the model abundances in line with observations. One
could argue the H3O+ abundance is reasonable given the uncertainties in
chemical models but the SO abundance is two orders of magnitude higher than
that found from our RADEX fits. Since rectifying this with an XDR model would
require an input X-ray flux that is orders of magnitude higher than could be
expected from any observed X-ray source, we conclude that cosmic rays must be
driving most SO destruction.

##  6 Conclusions

Previous work has shown the cosmic ray ionization rate (CRIR) in the CMZ of
NGC 253 is large but has failed to constrain it to within an order of
magnitude. Therefore, in this work emission from SO and H3O+ toward several
positions in the CMZ of NGC 253 was analysed due to the fact that their
abundances are strongly tied to the CRIR.

We find that we can constrain the CRIR in these regions and that in every
location it is $\sim$${10}^{-13}\text{\,}{\mathrm{s}}^{-1}$
(${10}^{4}$$\zeta_{0}$), and is extremely unlikely to fall outside the range
${10}^{-14}\text{\,}{\mathrm{s}}^{-1}{10}^{-12}\text{\,}{\mathrm{s}}^{-1}$. We
validate this result by using our inferred parameter values to model the
abundances of previously detected species that have been used as probes of the
CRIR and find that we reproduce the observed abundances.

Previous ALCHEMI studies have been unable to rule out alternative processes
that could explain the observed abundances of their CRIR probes. However, our
inferred chemical abundances strongly disfavour shock dominated chemistry and
the high H3O+ abundance found in the GMCs of NGC 253 disfavours PDR chemistry.
Further, the low SO abundance is unattainable with purely X-ray driven
chemistry. Since all of the CRIR tracers presented by ALCHEMI can be
consistently modelled through cosmic ray dominated chemistry, we now consider
this to be the most likely scenario.

We find no evidence for variation in the CRIR between the GMCs we have
modelled. Whilst GMCs 3 and 7 have the lowest most likely CRIR values, the
uncertainties in those values are sufficiently large that they overlap
considerably with the CRIR values derived for other GMCs. Given that the model
already fails to capture the complexity of the data, it is likely a more
accurate model is required rather than more data in order to better constrain
the CRIR in each GMC.

We thank F. Priestley for helpful discussions on XDR modelling. We also thank
the anonymous reviewer for their insightful comments on this manuscript. This
work is part of a project that has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme MOPPEX 833460. This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2017.1.00161.L, ADS/JAO.ALMA#2018.1.00162.S., and 2016.1.01285.S.
ALMA is a partnership of ESO (representing its member states), NSF (USA) and
NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio
Astronomy Observatory is a facility of the National Science Foundation
operated under cooperative agreement by Associated Universities, Inc. The work
of YN was supported by NAOJ ALMA Scientific Research Grant No. 2017-06B and
JSPS KAKENHI grant No. JP18K13577. V.M.R. and L.C. have received funding from
the Comunidad de Madrid through the Atracción de Talento Investigador
(Doctores con experiencia) Grant (COOL: Cosmic Origins Of Life;
2019-T1/TIC-15379). L.C. has also received partial support from the Spanish
State Research Agency (AEI; project number PID2019-105552RB-C41).

## References

  * Bayet et al. (2011) Bayet, E., Williams, D. A., Hartquist, T. W., & Viti, S. 2011, Monthly Notices of the Royal Astronomical Society, 414, 1583, doi: 10.1111/j.1365-2966.2011.18500.x
  * Bell et al. (2006) Bell, T. A., Hartquist, T., Viti, S., & Williams, D. A. 2006, \Aap, 459, 805, doi: 10.1051/0004-6361:20065624
  * Buchner (2021) Buchner, J. 2021, Journal of Open Source Software, 6, 3001.  https://arxiv.org/abs/2101.09604v2
  * Charnley (1997) Charnley, S. B. 1997, The Astrophysical Journal, 481, 396, doi: 10.1086/304011
  * Demes et al. (2021) Demes, S., Lique, F., Faure, A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 509, 1252, doi: 10.1093/MNRAS/STAB3015
  * Gerin et al. (2010) Gerin, M., De Luca, M., Black, J., et al. 2010, Astronomy and Astrophysics, Volume 518, id.L110, ¡NUMPAGES¿5¡/NUMPAGES¿ pp., 518, L110, doi: 10.1051/0004-6361/201014576
  * Goldsmith (2001) Goldsmith, P. F. 2001, The Astrophysical Journal, 557, 736, doi: 10.1086/322255
  * González-Alfonso et al. (2013) González-Alfonso, E., Fischer, J., Bruderer, S., et al. 2013, Astronomy and Astrophysics, 550, A25, doi: 10.1051/0004-6361/201220466
  * Harada et al. (2021) Harada, N., Martín, S., Mangum, J. G., et al. 2021, Apj
  * Holdship et al. (2017) Holdship, J., Viti, S., Jiménez-Serra, I., Makrymallis, A., & Priestley, F. 2017, Astronomical Journal, 154, 38, doi: 10.3847/1538-3881/aa773f
  * Holdship et al. (2021) Holdship, J., Viti, S., Martín, S., et al. 2021, Astronomy & Astrophysics, 654, A55, doi: 10.1051/0004-6361/202141233
  * Indriolo et al. (2015) Indriolo, N., Neufeld, D. A., Gerin, M., et al. 2015, Astrophysical Journal, 800, 40, doi: 10.1088/0004-637X/800/1/40
  * Jenkins (2009) Jenkins, E. B. 2009, Astrophysical Journal, 700, 1299, doi: 10.1088/0004-637X/700/2/1299
  * Jiménez-Serra et al. (2008) Jiménez-Serra, I., Caselli, P., Martín-Pintado, J., & Hartquist, T. 2008, Astronomy and Astrophysics, 482, 549, doi: 10.1051/0004-6361:20078054
  * Lehmer et al. (2013) Lehmer, B. D., Wik, D. R., Hornschemeier, A. E., et al. 2013, The Astrophysical Journal, 771, 134, doi: 10.1088/0004-637X/771/2/134
  * Leroy et al. (2015) Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015, Astrophysical Journal, 801, 25, doi: 10.1088/0004-637X/801/1/25
  * Leroy et al. (2018) —. 2018, The Astrophysical Journal, 869, 126, doi: 10.3847/1538-4357/aaecd1
  * Lique et al. (2006) Lique, F., Dubernet, M.-L., Spielfiedel, A., & Feautrier, N. 2006, A&A, 450, 399, doi: 10.1051/0004-6361:20054520
  * Mangum et al. (2019) Mangum, J. G., Ginsburg, A. G., Henkel, C., et al. 2019, The Astrophysical Journal, 871, 170, doi: 10.3847/1538-4357/aafa15
  * Marble et al. (2010) Marble, A. R., Engelbracht, C. W., Van Zee, L., et al. 2010, Astrophysical Journal, 715, 506, doi: 10.1088/0004-637X/715/1/506
  * Martín et al. (2019) Martín, S., Martín-Pintado, J., Blanco-Sánchez, C., et al. 2019, Astronomy & Astrophysics, 631, A159, doi: 10.1051/0004-6361/201936144
  * Martín et al. (2006) Martín, S., Mauersberger, R., Martín-Pintado, J., Henkel, C., & García-Burillo, S. 2006, The Astrophysical Journal Supplement Series, 164, 450, doi: 10.1086/503297
  * Martín et al. (2021) Martín, S., Mangum, J. G., Harada, N., et al. 2021, Astronomy & Astrophysics, 656, A46, doi: 10.1051/0004-6361/202141567
  * McElroy et al. (2013) McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, Astronomy & Astrophysics, 550, A36, doi: 10.1051/0004-6361/201220465
  * Muller et al. (2016) Muller, S., Müller, H. S., Black, J. H., et al. 2016, Astronomy and Astrophysics, 595, A128, doi: 10.1051/0004-6361/201629073
  * Padovani et al. (2009) Padovani, M., Galli, D., & Glassgold, A. E. 2009, Astronomy and Astrophysics, 501, 619, doi: 10.1051/0004-6361/200911794
  * Padovani et al. (2018) Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, Astronomy & Astrophysics, 614, A111, doi: 10.1051/0004-6361/201732202
  * Phillips et al. (1992) Phillips, T. G., van Dishoeck, E. F., & Keene, J. 1992, The Astrophysical Journal, 399, 533, doi: 10.1086/171945
  * Priestley et al. (2017) Priestley, F., Barlow, M. J., & Viti, S. 2017, Monthly Notices of the Royal Astronomical Society, 472, 4444, doi: 10.1093/MNRAS/STX2327
  * Rekola et al. (2005) Rekola, R., Richer, M. G., McCall, M. L., et al. 2005, Monthly Notices of the Royal Astronomical Society, 361, 330, doi: 10.1111/j.1365-2966.2005.09166.x
  * Sakamoto et al. (2011) Sakamoto, K., Mao, R. Q., Matsushita, S., et al. 2011, Astrophysical Journal, 735, 19, doi: 10.1088/0004-637X/735/1/19
  * Schöier et al. (2005) Schöier, F. L. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, Astronomy & Astrophysics, 432, 369, doi: 10.1051/0004-6361:20041729
  * Turner et al. (1985) Turner, J. L., Ho, P. T. P., Turner, J. L., & Ho, P. T. P. 1985, ApJL, 299, L77, doi: 10.1086/184584
  * Ulvestad et al. (1997) Ulvestad, J. S., Antonucci, R. R. J., Ulvestad, J. S., & Antonucci, R. R. J. 1997, ApJ, 488, 621, doi: 10.1086/304739
  * van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H. H., Schöier, F. L. L., Jansen, D. J. J., & Van Dishoeck, E. F. 2007, Astronomy and Astrophysics, 468, 627, doi: 10.1051/0004-6361:20066820
  * Viti et al. (2014) Viti, S., García-Burillo, S., Fuente, A., et al. 2014, Astronomy & Astrophysics, 570, A28, doi: 10.1051/0004-6361/201424116
  * Williams & Viti (2013) Williams, D. A., & Viti, S. 2013, Observational Molecular Astronomy, Vol. 1 (Cambridge: Cambridge University Press), doi: 10.1017/CBO9781139087445

##  Appendix A The ALCHEMI H3O+and SO Data

In this section, we present the data analysed in this work. In Figures 10 and
11, we show additional moment 0 maps of NGC 253’s CMZ. We also provide an
example of the data included in the supplementary data table. Table 4 shows
the transitions and measured line intensities detected in GMC 4, the
equivalent values for the other GMCs are available in the online article.

Figure 10: Additional H3O+ integrated intensity (moment 0) images toward NGC
253. Markings, intensity scaling, and contours in each panel same as for
Figure 1. The continuum RMS values for the transitions shown are 10.0 (top)
and 1.0 (bottom) mJy/beam, respectively. Note that the field of view for the
H3O+ $3_{0}-2_{0}$ transition is 22 arcsec centered at
$\alpha=00^{h}47^{m}33.134^{s}$,
$\delta=-25^{\circ}17^{\prime}19.68^{\prime\prime}$ (ICRS, shown as a grey
dashed circle; see Section 2.2).

Figure 11: Additional SO integrated intensity (moment 0) images toward NGC 253. Markings, intensity scaling, and contours in each panel same as for Figure 1. The continuum RMS values for the transitions shown are 0.17 (top left), 0.33 (top right), 0.33 (lower left), and 0.33 (lower right) mJy/beam, respectively. Table 4: Detected transitions, their properties and measured intensities for GMC 4. Region | SpeciesaaFor machine readable format ”+” replaced by ”p” and transitions listed as J.K-1.K+1 and J.K for H3O+ and SO, respectively. | TransitionaaFor machine readable format ”+” replaced by ”p” and transitions listed as J.K-1.K+1 and J.K for H3O+ and SO, respectively. | Frequency | EU | RA | Dec | Corrected Integrated | Integrated Intensity | Interloper Correction
---|---|---|---|---|---|---|---|---|---
|  |  | (GHz) | (K) |  |  | Intensity (K km/s) | Uncertainty (K km/s) |
4 | H3Op | 1.1.1-2.1.0 | 307.192410 | 79.5 | 11.8874 | -25.28888 | 13.509450 | 8.817727 | 0.23
4 | H3Op | 3.2.0-2.2.1 | 364.797430 | 139.3 | 11.8874 | -25.28888 | 61.744377 | 10.782524 | 0.86
4 | H3Op | 3.0.1-2.0.0 | 396.272410 | 169.1 | 11.8874 | -25.28888 | 18.992481 | 3.812451 | 1.00
4 | SO | 2.3-1.2 | 99.299870 | 9.2 | 11.8874 | -25.28888 | 47.964352 | 7.209908 | 1.00
4 | SO | 1.1-0.1 | 286.340152 | 15.2 | 11.8874 | -25.28888 | 0.591833 | 0.354076 | 1.00
4 | SO | 3.4-2.3 | 138.178600 | 15.9 | 11.8874 | -25.28888 | 45.781056 | 6.875517 | 1.00
4 | SO | 2.2-1.1 | 86.093950 | 19.3 | 11.8874 | -25.28888 | 13.188165 | 2.028540 | 1.00
4 | SO | 2.2-1.2 | 309.502444 | 19.3 | 11.8874 | -25.28888 | 0.973434 | 0.366215 | 1.00
4 | SO | 3.2-2.1 | 109.252220 | 21.1 | 11.8874 | -25.28888 | 20.209262 | 3.069187 | 1.00
4 | SO | 4.5-3.4 | 178.605403 | 24.4 | 11.8874 | -25.28888 | 32.354178 | 6.111951 | 1.00
4 | SO | 3.3-2.2 | 129.138923 | 25.5 | 11.8874 | -25.28888 | 20.025210 | 3.028742 | 1.00
4 | SO | 4.3-3.2 | 158.971811 | 28.7 | 11.8874 | -25.28888 | 28.146298 | 4.251155 | 1.00
4 | SO | 4.4-3.3 | 172.181403 | 33.8 | 11.8874 | -25.28888 | 20.281011 | 3.075113 | 1.00
4 | SO | 5.6-4.5 | 219.949442 | 35.0 | 11.8874 | -25.28888 | 34.456452 | 5.173552 | 1.00
4 | SO | 5.4-4.3 | 206.176005 | 38.6 | 11.8874 | -25.28888 | 27.974328 | 4.237034 | 1.00
4 | SO | 5.4-4.4 | 100.029640 | 38.6 | 11.8874 | -25.28888 | 4.935441 | 0.872624 | 1.00
4 | SO | 5.5-4.4 | 215.220653 | 44.1 | 11.8874 | -25.28888 | 18.644894 | 2.810657 | 1.00
4 | SO | 6.7-5.6 | 261.843721 | 47.6 | 11.8874 | -25.28888 | 28.611328 | 5.660309 | 0.76
4 | SO | 6.5-5.4 | 251.825770 | 50.7 | 11.8874 | -25.28888 | 20.277435 | 5.747689 | 0.53
4 | SO | 6.5-5.5 | 136.634799 | 50.7 | 11.8874 | -25.28888 | 26.897403 | 4.049473 | 1.00
4 | SO | 6.6-5.5 | 258.255826 | 56.5 | 11.8874 | -25.28888 | 18.888161 | 2.867306 | 1.00
4 | SO | 7.8-6.7 | 304.077844 | 62.1 | 11.8874 | -25.28888 | 38.849257 | 5.843079 | 1.00
4 | SO | 7.6-6.5 | 296.550064 | 64.9 | 11.8874 | -25.28888 | 17.696586 | 2.673156 | 1.00
4 | SO | 7.7-6.6 | 301.286124 | 71.0 | 11.8874 | -25.28888 | 17.585570 | 2.662084 | 1.00
4 | SO | 8.9-7.8 | 346.528481 | 78.8 | 11.8874 | -25.28888 | 18.149714 | 2.795038 | 1.00
4 | SO | 8.7-7.6 | 340.714155 | 81.2 | 11.8874 | -25.28888 | 14.417925 | 2.531354 | 0.87
4 | SO | 8.8-7.7 | 344.310612 | 87.5 | 11.8874 | -25.28888 | 11.763875 | 1.843587 | 0.99
4 | SO | 9.8-8.8 | 254.573628 | 99.7 | 11.8874 | -25.28888 | 29.755258 | 4.477287 | 1.00

##  Appendix B Corner Plots

In this section, we include the corner plots for GMCs 3, 5, 6, and 7 which are
not in the main article.

Figure 12: Corner plot showing the marginalized posterior probability
distribution of every parameter and the joint distributions of all parameter
pairs for GMC 3. Figure 13: Similar to Fig. 12 for GMC 5 Figure 14: Similar to
Fig. 12 for GMC 6 Figure 15: Similar to Fig. 12 for GMC 7