# Shock-induced HCNH+ abundance enhancement in the heart of the starburst galaxy NGC 253 unveiled by ALCHEMI

Y. Gong, C. Henkel, C T Bop, J. G. Mangum, E. Behrens, F J Du, S.B.

Zhang, S. Martin, K. M. Menten, N Harada, et al.

► To cite this version:

Y. Gong, C. Henkel, C T Bop, J. G. Mangum, E. Behrens, et al. Shock-induced HCNH+ abundance enhancement in the heart of the starburst galaxy NGC 253 unveiled by ALCHEMI. Astronomy & Astrophysics - A&A, 2025, 696, pp.A31. 10.1051/0004-6361/202452835 hal-05080104v2

HAL Id: hal-05080104

https://hal.science/hal-05080104v2

Submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Shock-induced HCNH⁺ abundance enhancement in the heart of the starburst galaxy NGC 253 unveiled by ALCHEMI

Y. Gong<sup>1,2,*</sup> ◐, C. Henkel<sup>2,3</sup> ◐, C. T. Bop<sup>4,5</sup> ◐, J. G. Mangum<sup>6</sup> ◐, E. Behrens<sup>6,7</sup> ◐, F. J. Du<sup>1</sup>, S. B. Zhang<sup>1</sup>, S. Martin<sup>8,9</sup> ◐, K. M. Menten<sup>2</sup>, N. Harada<sup>10,11,12</sup> ◐, M. Bouvier<sup>13</sup> ◐, X. D. Tang<sup>3,14</sup> ◐, K. Tanaka<sup>15</sup> ◐, S. Viti<sup>13,16,17</sup>, Y. T. Yan<sup>2</sup> ◐, W. Yang<sup>18</sup> ◐, R. Q. Mao<sup>1</sup> ◐, and D. H. Quan<sup>19</sup>

Received 31 October 2024 / Accepted 26 February 2025

## ABSTRACT

**Context.** Understanding the chemistry of molecular clouds is pivotal to elucidate star formation and galaxy evolution. As one of the important molecular ions, HCNH⁺ plays an important role in this chemistry. Yet, its behavior and significance under extreme conditions, such as in the central molecular zones (CMZs) of external galaxies, are still largely unexplored.

**Aims.** We aim to reveal the physical and chemical properties of the CMZ in the starburst galaxy NGC 253 with multiple HCNH⁺ transitions to shed light on the molecule's behavior under the extreme physical conditions of a starburst.

**Methods.** We employed molecular line data including results for four rotational transitions of HCNH⁺ from the ALMA Comprehensive High-resolution Extragalactic Molecular Inventory (ALCHEMI) large program to investigate underlying physical and chemical processes.

**Results.** Despite weak intensities, HCNH⁺ emission is widespread throughout NGC 253’s CMZ, which suggests that this molecular ion can effectively trace large-scale structures within molecular clouds. Using the quantum mechanical coupled states’ approximation, we computed rate coefficients for collisions of HCNH⁺ with *para*-H₂ and *ortho*-H₂ at kinetic temperatures up to 500 K. Using these coefficients in a non-local-thermodynamic-equilibrium (non-LTE) modeling framework and employing a Monte Carlo Markov chain analysis, we find that HCNH⁺ emission originates from regions with H₂ number densities of ~ 10<sup>2.80</sup>–10<sup>3.55</sup> cm<sup>-3</sup>, establishing HCNH⁺ as a tracer of low-density environments. Our analysis reveals that most of the HCNH⁺ abundances in the CMZ of NGC 253 are higher than all values reported in the Milky Way. We perform static, photodissociation region, and shock modeling, and found that recurrent shocks could potentially account for the elevated HCNH⁺ abundances observed in this CMZ.

**Conclusions.** We propose that the unexpectedly high HCNH⁺ abundances may result from chemical enhancement, primarily driven by the elevated gas temperatures and cosmic ray ionization rates of shocked, low-density gas in the nuclear starburst regions of NGC 253.

**Key words.** ISM: clouds – ISM: molecules – radio lines: ISM – ISM: individual objects: NGC 253

<footnote n="1">¹ Purple Mountain Observatory, and Key Laboratory of Radio Astronomy, Chinese Academy of Sciences, 10 Yuanhua Road, Nanjing 210023, PR China</footnote>

<footnote n="2">² Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote n="3">³ Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 830011 Urumqi, PR China</footnote>

<footnote n="4">⁴ Univ Rennes, CNRS, IPR (Institut de Physique de Rennes) – UMR 6251, 35000 Rennes, France</footnote>

<footnote n="5">⁵ Nantes Université, CNRS, CEISAM, UMR 6230, 44000 Nantes, France</footnote>

<footnote n="6">⁶ National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475, USA</footnote>

<footnote n="7">⁷ Department of Astronomy, University of Virginia, PO Box 400325, 530 McCormick Road, Charlottesville, VA 22904-4325, USA</footnote>

<footnote n="8">⁸ European Southern Observatory, Alonso de Córdova, 3107, Vitacura, Santiago 763-0355, Chile</footnote>

<footnote n="9">⁹ Joint ALMA Observatory, Alonso de Córdova, 3107, Vitacura, Santiago 763-0355, Chile</footnote>

<footnote n="10">¹⁰ National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="11">¹¹ Institute of Astronomy and Astrophysics, Academia Sinica, 11F of AS/NTU Astronomy-Mathematics Building, No.1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan</footnote>

<footnote n="12">¹² Department of Astronomy, School of Science, The Graduate University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo 181-1855, Japan</footnote>

<footnote n="13">¹³ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="14">¹⁴ Key Laboratory of Radio Astronomy, Chinese Academy of Sciences, 830011 Urumqi, PR China</footnote>

<footnote n="15">¹⁵ Department of Physics, Faculty of Science and Technology, Keio University, 3-14-1 Hiyoshi, Yokohama, Kanagawa 223-8522, Japan</footnote>

<footnote n="16">¹⁶ Transdisciplinary Research Area (TRA) ‘Matter’/Argelander-Institut für Astronomie, University of Bonn, Bonn, Germany</footnote>

<footnote n="17">¹⁷ Physics and Astronomy, University College London, London, UK</footnote>

<footnote n="18">¹⁸ School of Astronomy & Space Science, Nanjing University, 163 Xianlin Avenue, Nanjing 210023, PR China</footnote>

<footnote n="19">¹⁹ Research Center for Intelligent Computing Platforms, Zhejiang Laboratory, Hangzhou 311100, PR China</footnote>

<footnote>* Corresponding author; ygong@pmo.ac.cn</footnote>
# 1. Introduction

Understanding the chemistry of molecular clouds in various environments affected by different physical and chemical boundary conditions is pivotal for gaining insights into star formation and galaxy evolution. Molecular ions are crucial in shaping the chemistry of molecular clouds through ion-neutral reactions (e.g., Agúndez & Wakelam 2013). As the simplest protonated nitrile, HCNH⁺ was first detected in the interstellar medium (ISM) almost four decades ago (Ziurys & Turner 1986). Previous studies suggested that the abundance of HCNH⁺ varied with different stages of star formation (Fontani et al. 2021; Gong et al. 2023), indicating its potential as an evolutionary signpost. Numerical simulations further proposed that HCNH⁺, along with its chemically related species HCN, could help probe ambipolar diffusion in star formation regions (Tritsis et al. 2023). This molecular ion was also suggested to play an important role in the formation of formamide (NH₂CHO; Kerkeni & Simmie 2023) which is a key prebiotic precursor (López-Sepulcre et al. 2019). Despite its importance in star formation and astrochemistry, HCNH⁺ received comparatively less attention than more commonly studied molecular ions like HCO⁺ and N₂H⁺. Because of its small dipole moment of 0.29 D (Botschwina 1986) and the resultant low H₂ critical densities (Gong et al. 2023), HCNH⁺ was expected to trace large-scale structures within molecular clouds across galaxies. However, previous observations of HCNH⁺ were predominantly confined to Galactic molecular clouds on linear scales of ≤10 pc (e.g., Ziurys & Turner 1986; Schilke et al. 1991; Quénard et al. 2017; Gong et al. 2023), leaving its behavior on larger scales and its role under extreme conditions largely unexplored.

At a distance of ~3.5 Mpc (Karachentsev et al. 2003; Rekola et al. 2005), the Sculptor galaxy NGC 253 stands as one of the nearest and most prominent starburst galaxies. This archetypal starburst galaxy boasts a total star formation rate of ~4.2 M<sub>☉</sub> yr<sup>-1</sup> (Sanders et al. 2003) with ~2 M<sub>☉</sub> yr<sup>-1</sup> concentrated in the nuclear starburst region, its central molecular zone (CMZ) (Bendo et al. 2015; Leroy et al. 2015). The energetic starburst activity drives large-scale galactic outflows (e.g., Turner 1985; Strickland et al. 2000; Bolatto et al. 2013; Bao et al. 2024). The starburst is fueled by inflowing gas from a bar (Iodice et al. 2014) and is characterized by the CMZ which harbors a total molecular gas mass of ~10<sup>8</sup> M<sub>☉</sub> (e.g., Sakamoto et al. 2011). The CMZ is significantly influenced by starburst feedback (e.g., Harada et al. 2021), and has proven to be a rich reservoir of molecular species (e.g., Martín et al. 2006; Aladro et al. 2015; Martín et al. 2021; Harada et al. 2024). Recent studies have highlighted elevated cosmic-ray ionization rates (e.g., Harada et al. 2021; Holdship et al. 2022; Behrens et al. 2022, 2024; Phan et al. 2024), large-scale shocks (e.g., Humire et al. 2022; Harada et al. 2022; Huang et al. 2023), and a large amount of warm or hot molecular gas (e.g., Mangum et al. 2019; Rico-Villas et al. 2020; Huang et al. 2023; Tanaka et al. 2024; Bouvier et al. 2024) within the CMZ. These extreme conditions make the CMZ of NGC 253 an exceptional target for investigating the behavior of HCNH⁺ under extreme physical environments.

Despite its potential, the weak intensities of its rotational transitions pose significant challenges for mapping its large-scale distribution of HCNH⁺ in galaxies. Thanks to the unprecedented sensitivity of the ALMA Comprehensive High-resolution Extra-galactic Molecular Inventory (ALCHEMI) survey (Martín et al. 2021), HCNH⁺ has been successfully discovered in the CMZ of NGC 253 (Harada et al. 2024) where only one single HCNH⁺ transition has been reported because of the selection criteria of

the study. Recent advancements in the calculation of excitation rate coefficients for HCNH⁺ collisions with H₂ have enabled the modeling of multiple HCNH⁺ transitions to determine physical conditions within molecular clouds (Bop & Lique 2023; Bop et al. 2024). However, such modeling has yet to be applied in an extragalactic environment. In this study, we aimed to investigate physical and chemical properties of the CMZ of NGC 253 with multiple transitions of HCNH⁺, shedding light on its behavior under the extreme starburst physical conditions.

We describe our observations in Sect. 2. Our findings are presented in Sect. 3, with discussions of the results in Sect. 4. The summary and conclusions are outlined in Sect. 5.

# 2. Observations and data reduction

The data used in this study were directly taken from the ALCHEMI survey (co-PIs.: S. Martín, N. Harada, and J. Mangum; Martín et al. 2021). A comprehensive description of the survey was already provided in Martín et al. (2021), so we only offer a brief summary here. The ALCHEMI survey (project codes: 2017.1.00161.L, 2018.1.00162.S) spans ALMA Bands 3–7, covering frequencies from 84 to 375 GHz, excluding regions affected by atmospheric absorption features. The spectroscopic information was taken from the entry in the Cologne Database for Molecular Spectroscopy (CDMS<fnref n="1" />, Endres et al. 2016), which is based on the laboratory results of Silva et al. (2024). HCNH⁺ (2–1), HCNH⁺ (3–2), HCNH⁺ (4–3), and HCNH⁺ (5–4) have rest frequencies at 148221.450(17) MHz, 222329.277(8) MHz, 296433.652(17) MHz, and 370533.375(23) MHz, respectively, where uncertainties in the last digits are given in parentheses. These four HCNH⁺ rotational transitions fall within the frequency range covered by the ALCHEMI survey. The phase center of these observations is (α<sub>J2000</sub>, δ<sub>J2000</sub>)=(00h47m33s28, -25°17′17″6). The surveyed area of 50″×20″ encompasses the entire CMZ of NGC 253. All data products were convolved to a common circular beam of 1″6 (i.e., 27 pc). In this work, we mainly use the transitions of HCNH⁺ for our analysis, and their maximum recoverable scales of ≥11″ was achieved by combining the 12 m Array and Atacama Compact Array data. The data cubes were extracted from the ALCHEMI continuum-subtracted FITS cubes using the SpectralCube<fnref n="2" /> python package. The observed flux densities were converted to the brightness temperature scale using the Rayleigh-Jeans formula. The absolute flux calibration uncertainties were found to be ≤10% (Martín et al. 2021). The typical rms noise levels are about 4–20 mK at a channel width of 10 km s<sup>-1</sup>. The velocity scale was calibrated with respect to the local standard of rest (LSR) in this work.

# 3. Results

## 3.1. Spatial distribution and spectra

Figure 1 illustrates the elongated distribution of HCNH⁺ (2–1) emission within the CMZ of NGC 253, with the nine most prominent giant molecular clouds (GMCs) labeled (see Table 1 in Harada et al. 2024). Figure 1a is adapted from Fig. 21 in Harada et al. (2024), while Fig. 1b shows the integrated intensity map convolved to a circular beam of 4″ FWHM, without sigma clipping. Although Fig. 1b does not recover the large-scale emission filtered out by the ALMA array, the smoothed image

<footnote n="1">¹ https://cdms.astro.uni-koeln.de/</footnote>

<footnote n="2">² https://spectral-cube.readthedocs.io/en/latest/</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Distribution of HCNH⁺ (2–1) in the CMZ of NGC 253. (a) HCNH⁺ (2–1) integrated intensity map (taken from Harada et al. 2024), clipped at the 3σ level. The contours start from 2.7 K km s⁻¹ (3σ) and increase by 2.7 K km s⁻¹. (b) HCNH⁺ (2–1) integrated intensity map that has been convolved to a circular beam of 4′′. The smoothed image manifests the large-scale structures traced by HCNH⁺ transitions. The contours start from 0.9 K km s⁻¹ (3σ) and increase by 0.9 K km s⁻¹. In each panel, the beam is shown in the lower-left corner, and the pluses indicate the positions of the nine most prominent GMCs.</figcaption></figure>

demonstrates that HCNH⁺ emission is more widespread than previously thought, with a radial extent of ~300 pc in the CMZ. This is because the sensitivity is not high enough in Fig. 1a to detect the weak signals seen in Fig. 1b. Especially, the emission peaks are largely offset with respect to the positions of GMC1 and GMC2 in the low-angular-resolution image (see Fig. 1b), which was caused by the smoothing process where low-density gas close to GMC1 and GMC2 substantially contributes to the HCNH⁺ distribution. This behavior is expected due to the rather low upper-state energy levels and critical densities (see Table 1 in Gong et al. 2023). Consequently, HCNH⁺ emission effectively traces large-scale structures within molecular clouds.

Figure 2 presents the observed HCNH⁺ spectra for the nine GMCs (see their positions in Fig. 1). HCNH⁺ is detected in eight of the GMCs, with at least two transitions in each GMC exceeding the 3σ level, except for GMC5, the nucleus and kinematic center of the galaxy, which is associated with the brightest radio continuum emission (Turner & Ho 1985). The peak velocities of HCNH⁺ transitions agree well with those of all other molecular species, supporting the assignment of these spectral features to HCNH⁺. Additionally, there are spectral features near the HCNH⁺ (3–2) and HCNH⁺ (4–3) lines. The features close to HCNH⁺ (3–2) correspond to multiple k-ladder transitions of CH₃CCH (12k − 11k), while those near HCNH⁺ (4–3) are likely due to the blend of SO (67 − 56) with a rest frequency of 296550.064(30) MHz, the highly excited SO₂ (244.20 − 243.21) with a rest frequency of 296535.4225(14) MHz, and ³⁰SiO (7–6) with a rest frequency of 296575.730(20) MHz (Endres et al. 2016). Thus, the classic integrated intensity maps over

a specific velocity range (i.e., from 100 to 400 km s⁻¹) inevitably introduce contamination to the HCNH⁺ (3–2) and (4–3) integrated intensity maps. However, the CH₃CCH features are redshifted by ~240 km s⁻¹ relative to HCNH⁺ (3–2), and the blended lines around 296550 MHz are blueshifted by ~100 km s⁻¹ with respect to HCNH⁺ (4–3), allowing us to isolate the observed properties of HCNH⁺ transitions. To minimize this contamination, we applied a masking method before producing the integrated-intensity maps.

Since the HCNH⁺ (2–1) line is brighter than the other HCNH⁺ transitions and less contaminated by emission from other species, we used emission in this line to create the mask both spatially and spectrally for the data cubes of all the HCNH⁺ transitions. Our masking method<fnref n="3" /> was a modified version of that previously adopted by Harada et al. (2024). We labeled the connected components above the 3σ threshold in the three-dimensional (α₂₀₀₀, δ₂₀₀₀, vₐLSR) data cube of HCNH⁺ (2–1) using a 26-connected neighborhood<fnref n="4" /> algorithm<fnref n="5" />. These labeled components were subsequently utilized as masks for all the HCNH⁺ transitions. The mask areas were carefully expanded in the spatial dimensions to include surrounding emission and in the spectral dimension to capture line wings while minimizing contamination from other spectral lines. Furthermore, we manually masked the regions near HCNH⁺ (4–3) emission toward GMC5 and GMC6 to mitigate potential contamination in the spectral axis. Data outside these masked regions were trimmed to focus the analysis on the relevant areas within the mask.

Based on the trimmed data cubes<fnref n="6" />, we derived the integrated intensity distributions of the four HCNH⁺ transitions over a range from 100 km s⁻¹ to 400 km s⁻¹. The results are presented in Fig. 3. The integrated intensities of Figs. 1a and 3a are consistent within uncertainties, with Fig. 3a showing a higher dynamical range and a more widespread distribution than Fig. 1a. This consistency suggests that our masking approach accurately reproduces the intrinsic integrated intensities and enhances the detection of low-intensity emissions. Figure 3 not only confirms the widespread distribution of HCNH⁺ emission but also shows that the four HCNH⁺ transitions exhibit nearly identical spatial distributions, with emission peaking at GMC6. The main difference is the different integrated intensities of these transitions, with HCNH⁺ (2–1) being the brightest. This indicates that HCNH⁺ transitions predominantly trace the same gas structures, albeit with varying intensities. While weak HCNH⁺ emission appears to be present in GMC5 in Fig. 3, its spectra lack sufficient signal-to-noise ratios in Fig. 2. Therefore, HCNH⁺ emissions in GMC5 were excluded from further analysis.

We fitted the observed HCNH⁺ transitions from the selected GMCs using a single-component Gaussian model, with the results presented in Table 1. The fitted velocity centroids range from ~180 to ~330 km s⁻¹ and full width at half-maximum (FWHM) line widths are in the range of 40–70 km s⁻¹. These values are consistent with previously reported velocities determined for other molecular tracers within 3σ uncertainties (e.g., Haasler et al. 2022; Humire et al. 2022; Bouvier et al. 2024). Based on the fitted results, we find that all HCNH⁺ transitions

<footnote n="3">³ The masking script is available from this link.</footnote>

<footnote n="4">⁴ In a single spectral plane, each pixel has 9 connected neighbors, including itself. In a data cube, each spectral plane has two adjacent spectral planes. Hence, each pixel has 27 connected neighbors in total. Excluding the pixel itself, 26-connected neighbors remain. If the central pixel has an intensity above 3σ threshold, the 26-connected neighbors are included in the mask.</footnote>

<footnote n="5">⁵ https://github.com/seung-lab/connected-components-3d</footnote>

<footnote n="6">⁶ The data cubes are available via this link.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Observed H<sup>13</sup>CN and HCNH<sup>+</sup> spectra (all with a common circular beam of 1.′′6) of the nine GMCs indicated in Fig. 1 and Table 1. H<sup>13</sup>CN (1–0), HCNH<sup>+</sup> (2–1), HCNH<sup>+</sup> (3–2), HCNH<sup>+</sup> (4–3), and HCNH<sup>+</sup> (5–4) are indicated by black, red, green, blue, and purple lines, respectively. The fits to the HCNH<sup>+</sup> spectra are indicated by black dotted lines. H<sup>13</sup>CN (1–0) spectra are scaled by a factor of 0.1 for better visualization. In each panel, the systemic LSR velocity is indicated by the vertical black dashed line. The potentially blended lines are marked by short vertical lines in the corresponding panels, with the associated species labeled in the central panel for GMC5.</figcaption></figure>

have peak intensities of ≤200 mK, which are much lower than those of H<sup>13</sup>CN (1–0) (see Fig. 2). Given the widespread distribution of HCNH<sup>+</sup> emission in the CMZ (see Fig. 1b), beam dilution effects might be negligible. Assuming an excitation temperature of 15 K (as discussed below) and neglecting beam dilution effects, we find that the low peak intensities correspond to optical depths of <0.1. Even accounting for a beam dilution factor as low as 0.1, the peak optical depths remain below 0.3. Therefore, we conclude that the HCNH<sup>+</sup> transitions are optically thin in the CMZ of NGC 253.

## 3.2. LTE analysis

Assuming local thermodynamic equilibrium (LTE), we employed the rotational diagram method (Goldsmith & Langer 1999) to derive the rotational temperatures, $T_{\text{rot}}$, and molecular column densities, $N_{\text{HCNH}^+}$. Since all data had been convolved to the same angular resolution of 1.′′6, we did not correct for the beam dilution effects, which might be negligible (see above). While this ensured that rotational temperatures were unaffected, the derived column densities were inherently beam-averaged.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Distributions of HCNH⁺ (2–1) (a), HCNH⁺ (3–2) (b), HCNH⁺ (4–3) (c), and HCNH⁺ (5–4) (d) emission in the CMZ of NGC 253. All intensity maps are integrated from 100 to 400 km s⁻¹ and share a common integrated-intensity scale indicated by the color bar. In each panel, the overlaid contours represent the integrated intensities of respective HCNH⁺ transitions, starting at 1.2 K km s⁻¹ and increasing by 1.2 K km s⁻¹. The positions of the nine GMCs are marked with plus symbols. In each panel, the beam is shown in the lower-left corner.</figcaption></figure>

The rotational diagrams of the eight GMCs are shown in Fig. 4. Linear fits were applied to these diagrams, with uncertainties quantified using a Monte Carlo approach. Specifically, 10 000 realizations were generated by sampling Gaussian distributions of the integrated intensities, adopting the combined uncertainties from Table 1 and 10% of the absolute flux calibration uncertainties as the respective standard deviations. These were propagated through linear fits, and the resulting distributions were analyzed to derive uncertainties in both the rotational temperatures and column densities.

The fitting results of the eight GMCs were given in Table 2. The derived rotational temperatures were in the range of 12–21 K, aligning closely with those of C¹⁸O but exceeding those (<10 K) of H¹³CO⁺ (see Table 4 in Harada et al. 2021). This discrepancy was likely attributed to the similarly low critical densities of HCNH⁺ and C¹⁸O transitions, which are much lower than those of H¹³CO⁺. The rotational temperatures were also lower than the expected gas temperatures (e.g., Mangum et al. 2019; Rico-Villas et al. 2020; Huang et al. 2023; Tanaka et al. 2024), which implied that HCNH⁺ transitions were subthermally excited. The HCNH⁺ column densities derived from our analysis ranged from 8.3 × 10¹⁴ cm⁻² to 4.0 × 10¹⁵ cm⁻², which is almost comparable to those of HCO⁺ in this CMZ (Harada et al. 2021).

$H_2$ column densities, $N_{H_2}$, are essential to derive the fractional molecular abundances, $X_{HCNH^+}$, of HCNH⁺. In this work, we directly utilized values for $N_{H_2}$ from the low-density hierarchical Bayesian analysis of the ALCHEMI data at a circular beam of 1.′′6 (Tanaka et al. 2024). The resulting $X_{HCNH^+}$, presented in Table 2, were remarkably high, ranging from 6.3 × 10⁻⁹ to 2.0 × 10⁻⁸. In comparison, the observed $X_{HCNH^+}$ values in the Milky Way were ≲2 × 10⁻⁹ in all previous studies (e.g., Ziurys & Turner 1986; Schilke et al. 1991; Fontani et al. 2021; Gong et al. 2023). Particularly, $X_{HCNH^+}$ was only ~3 × 10⁻¹⁰ in Sgr B2, located in the CMZ of the Milky Way (Ziurys & Turner 1986), which is more than an order of magnitude lower than our observed values in the CMZ of NGC 253. This stark contrast highlights the significantly elevated abundance of HCNH⁺ in the CMZ of NGC 253, indicating that its physical conditions differ significantly from those in the Milky Way.

Based on Table 2, we observed a trend of decreasing HCNH⁺ abundances in GMCs as they approach the center of the CMZ, which is similar to the gradient of HOCO⁺ abundances (see Fig. 4d in Harada et al. 2022). To further explore this, we examined the relationship between HCNH⁺ abundance and $H_2$ column density, as shown in Fig. 5. Our analysis reveals that the anti-correlation between HCNH⁺ abundances and $H_2$ column densities, previously observed in molecular clouds of the
<table><caption>**Table 1.** Observational parameters of the selected GMCs.</caption><thead><tr><th>Target<br>(1)</th><th>(α<sub>J2000</sub>, δ<sub>J2000</sub>)<br>(h:m:s, °:′:′′)</th><th>Transition<br>(3)</th><th>v<sub>LSR</sub><br>(km s<sup>-1</sup>)<br>(4)</th><th>Δv<sub>LSR</sub><br>(km s<sup>-1</sup>)<br>(5)</th><th>T<sub>p</sub><br>(mK)<br>(6)</th><th>∫ T<sub>b</sub> dv<br>(K km s<sup>-1</sup>)<br>(7)</th></tr></thead><tbody><tr><td rowspan="4">GMC1</td><td rowspan="4">(00:47:31.93, -25:17:29.0)</td><td>HCNH<sup>+</sup> (2-1)</td><td>301 ± 4</td><td>66 ± 9</td><td>52 ± 12</td><td>3.63 ± 0.45</td></tr><tr><td>HCNH<sup>+</sup> (3-2)</td><td>302 ± 4</td><td>42 ± 9</td><td>59 ± 6</td><td>2.66 ± 0.47</td></tr><tr><td>HCNH<sup>+</sup> (4-3)</td><td>305 ± 3</td><td>45 ± 8</td><td>43 ± 4</td><td>2.03 ± 0.29</td></tr><tr><td>HCNH<sup>+</sup> (5-4)</td><td>315 ± 6</td><td>40 ± 13</td><td>43 ± 8</td><td>1.83 ± 0.54</td></tr><tr><td rowspan="4">GMC2</td><td rowspan="4">(00:47:32.36, -25:17:18.8)</td><td>HCNH<sup>+</sup> (2-1)</td><td>327 ± 3</td><td>43 ± 8</td><td>61 ± 15</td><td>2.75 ± 0.40</td></tr><tr><td>HCNH<sup>+</sup> (3-2)</td><td>332 ± 3</td><td>36 ± 10</td><td>47 ± 6</td><td>1.81 ± 0.37</td></tr><tr><td>HCNH<sup>+</sup> (4-3)</td><td>339 ± 3</td><td>45 ± 9</td><td>32 ± 7</td><td>1.54 ± 0.24</td></tr><tr><td>HCNH<sup>+</sup> (5-4)</td><td>...</td><td>...</td><td>&lt;36</td><td>&lt;1.62</td></tr><tr><td rowspan="4">GMC3</td><td rowspan="4">(00:47:32.81, -25:17:21.2)</td><td>HCNH<sup>+</sup> (2-1)</td><td>287 ± 1</td><td>52 ± 3</td><td>181 ± 15</td><td>10.06 ± 0.51</td></tr><tr><td>HCNH<sup>+</sup> (3-2)</td><td>286 ± 1</td><td>48 ± 2</td><td>143 ± 6</td><td>7.30 ± 0.30</td></tr><tr><td>HCNH<sup>+</sup> (4-3)</td><td>282 ± 1</td><td>53 ± 3</td><td>123 ± 17</td><td>6.92 ± 0.30</td></tr><tr><td>HCNH<sup>+</sup> (5-4)</td><td>284 ± 3</td><td>51 ± 8</td><td>92 ± 13</td><td>5.02 ± 0.64</td></tr><tr><td rowspan="4">GMC4</td><td rowspan="4">(00:47:32.95, -25:17:19.8)</td><td>HCNH<sup>+</sup> (2-1)</td><td>243 ± 3</td><td>49 ± 6</td><td>208 ± 10</td><td>10.76 ± 1.21</td></tr><tr><td>HCNH<sup>+</sup> (3-2)</td><td>250 ± 1</td><td>46 ± 2</td><td>133 ± 11</td><td>6.54 ± 0.25</td></tr><tr><td>HCNH<sup>+</sup> (4-3)</td><td>247 ± 1</td><td>46 ± 2</td><td>143 ± 9</td><td>6.96 ± 0.20</td></tr><tr><td>HCNH<sup>+</sup> (5-4)</td><td>244 ± 2</td><td>47 ± 4</td><td>120 ± 13</td><td>6.00 ± 0.40</td></tr><tr><td rowspan="4">GMC5</td><td rowspan="4">(00:47:33.16, -25:17:17.3)</td><td>HCNH<sup>+</sup> (2-1)</td><td>...</td><td>...</td><td>&lt;41</td><td>&lt;1.85</td></tr><tr><td>HCNH<sup>+</sup> (3-2)</td><td>...</td><td>...</td><td>&lt;38</td><td>&lt;1.71</td></tr><tr><td>HCNH<sup>+</sup> (4-3)</td><td>...</td><td>...</td><td>&lt;42</td><td>&lt;3.24</td></tr><tr><td>HCNH<sup>+</sup> (5-4)</td><td>...</td><td>...</td><td>&lt;54</td><td>&lt;2.43</td></tr><tr><td rowspan="4">GMC6</td><td rowspan="4">(00:47:33.33, -25:17:15.7)</td><td>HCNH<sup>+</sup> (2-1)</td><td>179 ± 1</td><td>51 ± 3</td><td>225 ± 11</td><td>12.14 ± 0.57</td></tr><tr><td>HCNH<sup>+</sup> (3-2)</td><td>180 ± 1</td><td>53 ± 2</td><td>177 ± 14</td><td>10.05 ± 0.31</td></tr><tr><td>HCNH<sup>+</sup> (4-3)</td><td>178 ± 1</td><td>57 ± 2</td><td>190 ± 6</td><td>11.66 ± 0.36</td></tr><tr><td>HCNH<sup>+</sup> (5-4)</td><td>182 ± 2</td><td>59 ± 6</td><td>137 ± 13</td><td>8.63 ± 0.71</td></tr><tr><td rowspan="4">GMC7</td><td rowspan="4">(00:47:33.65, -25:17:13.1)</td><td>HCNH<sup>+</sup> (2-1)</td><td>178 ± 4</td><td>74 ± 8</td><td>82 ± 9</td><td>6.53 ± 0.68</td></tr><tr><td>HCNH<sup>+</sup> (3-2)</td><td>179 ± 2</td><td>66 ± 5</td><td>83 ± 6</td><td>5.83 ± 0.38</td></tr><tr><td>HCNH<sup>+</sup> (4-3)</td><td>179 ± 2</td><td>70 ± 5</td><td>71 ± 9</td><td>5.27 ± 0.31</td></tr><tr><td>HCNH<sup>+</sup> (5-4)</td><td>180 ± 5</td><td>56 ± 11</td><td>69 ± 15</td><td>4.12 ± 0.70</td></tr><tr><td rowspan="4">GMC8</td><td rowspan="4">(00:47:33.94, -25:17:10.9)</td><td>HCNH<sup>+</sup> (2-1)</td><td>196 ± 3</td><td>50 ± 10</td><td>64 ± 13</td><td>3.40 ± 0.44</td></tr><tr><td>HCNH<sup>+</sup> (3-2)</td><td>197 ± 2</td><td>43 ± 5</td><td>39 ± 5</td><td>1.81 ± 0.18</td></tr><tr><td>HCNH<sup>+</sup> (4-3)</td><td>196 ± 3</td><td>49 ± 7</td><td>37 ± 10</td><td>1.93 ± 0.27</td></tr><tr><td>HCNH<sup>+</sup> (5-4)</td><td>...</td><td>...</td><td>&lt;43</td><td>&lt;2.13</td></tr><tr><td rowspan="4">GMC9</td><td rowspan="4">(00:47:34.14, -25:17:12.0)</td><td>HCNH<sup>+</sup> (2-1)</td><td>202 ± 6</td><td>66 ± 15</td><td>46 ± 11</td><td>3.21 ± 0.61</td></tr><tr><td>HCNH<sup>+</sup> (3-2)</td><td>207 ± 2</td><td>43 ± 4</td><td>51 ± 5</td><td>2.32 ± 0.18</td></tr><tr><td>HCNH<sup>+</sup> (4-3)</td><td>210 ± 3</td><td>50 ± 7</td><td>40 ± 7</td><td>2.08 ± 0.24</td></tr><tr><td>HCNH<sup>+</sup> (5-4)</td><td>...</td><td>...</td><td>&lt;45</td><td>&lt;2.27</td></tr></tbody></table>

**Notes.** (1) Selected GMC. (2) Equatorial coordinate. (3) Transition. (4) Velocity centroid. (5) Full width at half-maximum line width. (6) Peak main beam brightness temperature. (7) Integrated intensity.

solar neighborhood (Gong et al. 2023), also holds in the CMZ of NGC 253. However, the anti-correlation appeared to arise from the associated cosmic ray ionization rates (see Sect. 4), because GMCs close to the center of the CMZ tend to have higher cosmic ray ionization rates (Harada et al. 2021; Holdship et al. 2022; Behrens et al. 2022, 2024). In the solar neighborhood, cosmic ray ionization rates are likely governed by H₂ column densities. As suggested by previous studies (e.g., Neufeld & Wolfire 2017; Padovani et al. 2018), cosmic ray ionization rates decrease with increasing H₂ column density. If these rates were the primary factor regulating HCNH⁺ abundances, one would expect higher HCNH⁺ abundances in regions with higher H₂ column densities. However, this expectation was not consistent with recent

observations. Therefore, we propose that the physical mechanisms driving the observed anti-correlation between HCNH⁺ abundance and H₂ column density in the CMZ of NGC 253 differ from those in the solar neighborhood.

### 3.3. Non-LTE analysis

#### 3.3.1. Collisional rate coefficients

Rates for collisions of HCNH⁺ induced by encounters with *para*-H₂($j_2^7 = 0$), denoted hereafter as *p*-H₂, were first investigated by Bop & Lique (2023) using the time-independent quantum

<footnote n="7">⁷ $j_2$ stands for the rotational quantum number of H₂.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Rotational diagrams of HCNH⁺ for the selected GMCs. In each panel, the dashed line represents the linear least squares fit to the observed data points. The fit results for each GMC are shown in the lower-left corner of each panel.</figcaption></figure>

<table><caption>**Table 2.** Physical parameters of the selected GMCs.</caption><thead><tr><th rowspan="2">Target<br>(1)</th><th rowspan="2">N<sub>H2</sub><br>(cm<sup>-2</sup>)<br>(2)</th><th colspan="3">LTE</th><th colspan="4">RADEX</th></tr><tr><th>T<sub>rot</sub><br>(K)<br>(3)</th><th>N<sub>HCHN+</sub><br>(cm<sup>-2</sup>)<br>(4)</th><th>X<sub>HCHN+</sub><br>(5)</th><th>N<sub>HCHN+</sub><br>(cm<sup>-2</sup>)<br>(6)</th><th>X<sub>HCHN+</sub><br>(7)</th><th>T<sub>K</sub><br>(K)<br>(8)</th><th>log<sub>10</sub>[n<sub>H2</sub>]<br>(cm<sup>-3</sup>)<br>(9)</th></tr></thead><tbody><tr><td>GMC1</td><td>(5.0 ± 1.1) × 10<sup>22</sup></td><td>17.4±2.5</td><td>(1.0 ± 0.1) × 10<sup>15</sup></td><td>(2.0 ± 0.5) × 10<sup>-8</sup></td><td>4.0<sup>+0.9</sup><sub>-0.7</sub> × 10<sup>14</sup></td><td>8.6<sup>+2.4</sup><sub>-2.5</sub> × 10<sup>-9</sup></td><td>249.83<sup>+101.35</sup><sub>-117.63</sub></td><td>3.10<sup>+0.17</sup><sub>-0.16</sub></td></tr><tr><td>GMC2</td><td>(6.5 ± 1.2) × 10<sup>22</sup></td><td>12.9±1.8</td><td>(8.3 ± 1.3) × 10<sup>14</sup></td><td>(1.3 ± 0.3) × 10<sup>-8</sup></td><td>3.2<sup>+1.1</sup><sub>-0.7</sub> × 10<sup>14</sup></td><td>5.4<sup>+1.7</sup><sub>-1.7</sub> × 10<sup>-9</sup></td><td>240.34<sup>+110.20</sup><sub>-121.93</sub></td><td>3.01<sup>+0.23</sup><sub>-0.23</sub></td></tr><tr><td>GMC3</td><td>(2.2 ± 0.5) × 10<sup>23</sup></td><td>17.7±1.4</td><td>(3.0 ± 0.2) × 10<sup>15</sup></td><td>(1.4 ± 0.3) × 10<sup>-8</sup></td><td>9.8<sup>+1.2</sup><sub>-1.0</sub> × 10<sup>14</sup></td><td>4.8<sup>+1.1</sup><sub>-1.2</sub> × 10<sup>-9</sup></td><td>278.12<sup>+83.78</sup><sub>-106.70</sub></td><td>3.23<sup>+0.12</sup><sub>-0.11</sub></td></tr><tr><td>GMC4</td><td>(4.3 ± 0.8) × 10<sup>23</sup></td><td>19.0±1.5</td><td>(2.9 ± 0.2) × 10<sup>15</sup></td><td>(6.7 ± 1.3) × 10<sup>-9</sup></td><td>7.2<sup>+0.8</sup><sub>-0.6</sub> × 10<sup>14</sup></td><td>1.8<sup>+0.4</sup><sub>-0.4</sub> × 10<sup>-9</sup></td><td>252.49<sup>+102.07</sup><sub>-132.07</sub></td><td>3.55<sup>+0.22</sup><sub>-0.15</sub></td></tr><tr><td>GMC6</td><td>(6.3 ± 1.8) × 10<sup>23</sup></td><td>20.9±1.7</td><td>(4.0 ± 0.2) × 10<sup>15</sup></td><td>(6.3 ± 1.8) × 10<sup>-9</sup></td><td>1.2<sup>+0.1</sup><sub>-0.1</sub> × 10<sup>15</sup></td><td>2.1<sup>+0.6</sup><sub>-0.6</sub> × 10<sup>-9</sup></td><td>256.86<sup>+99.19</sup><sub>-123.33</sub></td><td>3.46<sup>+0.17</sup><sub>-0.11</sub></td></tr><tr><td>GMC7</td><td>(1.8 ± 0.3) × 10<sup>23</sup></td><td>19.2±2.0</td><td>(2.1 ± 0.2) × 10<sup>15</sup></td><td>(1.2 ± 0.2) × 10<sup>-8</sup></td><td>6.7<sup>+0.8</sup><sub>-0.7</sub> × 10<sup>14</sup></td><td>3.9<sup>+0.7</sup><sub>-0.8</sub> × 10<sup>-9</sup></td><td>250.09<sup>+102.97</sup><sub>-117.32</sub></td><td>3.35<sup>+0.17</sup><sub>-0.13</sub></td></tr><tr><td>GMC8</td><td>(5.9 ± 1.1) × 10<sup>22</sup></td><td>13.1±1.7</td><td>(9.5 ± 1.2) × 10<sup>14</sup></td><td>(1.6 ± 0.4) × 10<sup>-8</sup></td><td>4.5<sup>+3.6</sup><sub>-1.5</sub> × 10<sup>14</sup></td><td>9.4<sup>+4.9</sup><sub>-4.5</sub> × 10<sup>-9</sup></td><td>254.10<sup>+103.72</sup><sub>-118.16</sub></td><td>2.80<sup>+0.27</sup><sub>-0.33</sub></td></tr><tr><td>GMC9</td><td>(5.4 ± 1.1) × 10<sup>22</sup></td><td>13.9±2.1</td><td>(9.9 ± 1.5) × 10<sup>14</sup></td><td>(1.8 ± 0.5) × 10<sup>-8</sup></td><td>3.1<sup>+1.1</sup><sub>-0.7</sub> × 10<sup>14</sup></td><td>6.3<sup>+2.1</sup><sub>-2.1</sub> × 10<sup>-9</sup></td><td>245.05<sup>+107.46</sup><sub>-119.90</sub></td><td>3.16<sup>+0.25</sup><sub>-0.23</sub></td></tr></tbody></table>

**Notes.** (1) Selected GMC. (2) Beam-averaged H₂ column density for a circular beam of 1"6, taken from Tanaka et al. (2024). (3) Rotational temperature. (4) Beam-averaged HCNH⁺ column density, again for a circular beam of 1"6. (5) HCNH⁺ abundance. (6) Beam-averaged HCNH⁺ column density. (7) Non-LTE HCNH⁺ abundance. (8) Kinetic temperature. (9) H₂ number density.

mechanical close-coupling (CC) method (Arthur & Dalgarno 1960) and later revised to include the effect of the H₂ rotational basis (Gong et al. 2023; Bop et al. 2024). The former dataset, which was less precise, covered kinetic temperatures, T<sub>K</sub>, up to 100 K, while the latter more precise dataset was limited to 30 K due to the computational cost of the CC method. Given that the CMZ of NGC 253 is expected to be significantly warmer than 30 K (e.g., Mangum et al. 2019; Tanaka et al. 2024), a comprehensive set of collisional rate coefficients covering a wider temperature range is required for a robust interpretation of warm or hot HCNH⁺ emission lines in the CMZ. This would help in accurately modeling the excitation conditions of HCNH⁺ in such extreme environments.

To address the computational cost of the CC method, we employ the quantum mechanical coupled states' (CS) approximation (McGuire & Kouri 1974), which has long been tested

and validated for the computation of high-temperature collisional rate coefficients. We investigate the scattering of HCNH⁺ due to collisions with both *p*-H₂ and *ortho*-H₂ (*j*₂ = 1), denoted hereafter as *o*-H₂. The use of the latter collision partner is motivated by the warm/hot regime which we are interested in. Indeed, above 75 K, the population of molecular hydrogen is dominated by the *o*-H₂ spin-symmetry.

In practice, we implemented the four-dimensional potential energy surface calculated by Bop & Lique (2023) in the MOLSCAT molecular scattering code (Hutson & Green 1994) to compute state-to-state inelastic cross sections. We focused on transitions among the 30 lowest-lying energy levels of HCNH⁺, which, for convergence reasons, required an exhaustive rotational basis (*j*₁<fnref n="8" /> = 0–48). In addition, we included the two

<footnote n="8">⁸ *j*₁ stands for the rotational quantum number of HCNH⁺.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. HCNH⁺ abundances as a function of H₂ column density. The colors represent the distances of GMCs to GMC5 which is located at the center of the CMZ of NGC 253. The circles represent HCNH⁺ abundances derived using the LTE approach (see Sect. 3), while the squares represent those obtained from the non-LTE approach (see Sect. 3.3.2). The cyan-shaded region indicates the observed HCNH⁺ abundances of $\le 2 \times 10^{-9}$ previously reported in the Milky Way, while the brown-shaded region marks the abundances reported toward Sgr B2.</figcaption></figure>

low-lying rotational transitions of H₂, namely, $j_2 = 0-2$ and $j_2 = 1-3$ for p- and o-H₂, respectively. All additional integration parameters were defined to ensure subpercent convergence of the cross sections. The calculations were performed for total energies up to 3500 and 3620 cm⁻¹ using p-H₂ and o-H₂ as colliders, respectively. Therefore, using the Maxwell-Boltzmann velocity distribution, we derived rate coefficients for HCNH⁺ induced by collisions with both p- and o-H₂ for kinetic temperatures up to 500 K, and the results are shown in Fig. 6. The collisional rate coefficients will be accessible through databases such as Basecol<fnref n="9" /> (Dubernet et al. 2024), EMAA<fnref n="10" />, and LAMDA<fnref n="11" /> (Schöier et al. 2005; van der Tak et al. 2020).

To assess the accuracy of the CS approximation, we compare in Fig. 7 the new HCNH⁺ collisional rate coefficients induced by p-H₂ with those previously computed by Bop et al. (2024) using the “exact” CC method. The collisional rate coefficients using the CS approximation turn out to underestimate the CC data by up to a factor of two, especially for the low-magnitude rate coefficients. To quantify the disagreement between the CC and the CS results, we calculate for each temperature the weighted mean error factor (Loreau et al. 2018; Pirlot Jankowiak et al. 2023) as follows:

$$WMEF = \frac{\sum_i k_i^{CC} \times r_i}{\sum_i k_i^{CC}}, \qquad (1)$$

where $r_i$ is defined as the maximum of $k_i^{CC}/k_i^{CS}$ or $k_i^{CS}/k_i^{CC}$, ensuring that $r_i \ge 1$. The WMEF decreases with increasing temperature, suggesting that the agreement is expected to be reasonable above 30 K, good for $T > 50$ K, and excellent above 80 K. For further validation tests, we refer readers to

<figure><img src="image_7.png" /><figcaption>Fig. 6. Temperature dependence of the HCNH⁺ rate coefficients induced by collisions with p-H₂ (dashed lines) and o-H₂ (solid lines).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Comparison of two datasets of rate coefficients for HCNH⁺ induced by collisions with p-H₂($j_2 = 0$). The x-axis represents the first dataset, computed using the “exact” CC method (Bop & Lique 2023), while the y-axis represents the second dataset, calculated in this work using the CS approximation. The dashed lines indicate an agreement region within a factor of 2.0. Different colors represent their weighted mean error factors (WMEF) which are defined in Eq. (1).</figcaption></figure>

Appendix A, where the impact of low-temperature discrepancies on radiative transfer calculations is shown.

### 3.3.2. Non-LTE modeling

The detection of multiple HCNH⁺ transitions in the CMZ of NGC 253 allows the first non-LTE analysis of HCNH⁺ in an extragalactic environment. Here, we use the non-LTE radiative transfer code RADEX<fnref n="12" /> for statistical equilibrium calculations (van der Tak et al. 2007). The collisional rates used in the analysis are based on the data discussed in Sect. 3.3.1. In our modeling, RADEX adopts an expanding spherical geometry, the so-called Sobolev or large velocity gradient (LVG) approximation (e.g., Sobolev 1960; Elitzur 1992), to calculate the escape probability. Given that HCNH⁺ transitions are optically thin in the CMZ of NGC 253 (see Sect. 3), the results are not affected by

<footnote n="9">⁹ https://basecol.vamdc.eu/</footnote>

<footnote n="10">¹⁰ https://emaa.osug.fr/</footnote>

<footnote n="11">¹¹ https://home.strw.leidenuniv.nl/~moldata/</footnote>

<footnote n="12">¹² https://home.strw.leidenuniv.nl/~moldata/radex.html</footnote>
the assumed geometry. Thus, an expanding spherical geometry providing localized source functions is sufficiently robust.

In our RADEX modeling, we adopted an H₂ ortho-to-para ratio of 0.25, following the values reported in previous studies (e.g., Neufeld et al. 2006). The ortho-to-para ratio also does not affect the modeling results significantly, because the rate coefficients for collisions with p- and o-H₂ are comparable (see Fig. B.1). The line widths are taken to be the averages of the fitted values for the four transitions (see Table 1). We assume that the emission from all HCNH⁺ transitions originates from the same regions. Given that the beam sizes are identical for all transitions, beam dilution effects primarily influence the modeling of HCNH⁺ column densities when the transitions are optically thin. For simplicity, we neglected the beam dilution effects (i.e., the beam dilution factor is set to be unity). As a result, the HCNH⁺ column densities derived from our non-LTE analysis are beam-averaged values.

We minimized our modeling results using the emcee<fnref n="13" /> code (Foreman-Mackey et al. 2013) and the affine-invariant ensemble sampler (Goodman & Weare 2010) for Monte Carlo Markov chain (MCMC) calculations, enabling us to obtain the posterior probability distributions of the modeled parameters. We assumed uniform priors for $T_K$, $n_{H_2}$, and $N_{HCNH^+}$. Based on previous studies on the CMZ of NGC 253 (e.g., Mangum et al. 2019; Holdship et al. 2022; Behrens et al. 2022; Tanaka et al. 2024), we set $T_K$, $n_{H_2}$, and $N_{HCNH^+}$ to be in the range of 50–500 K, $10^2–10^7$ cm⁻³, and $10^{13}–10^{17}$ cm⁻², respectively. The posterior distribution of these parameters was determined by the product of the prior distribution function and the likelihood function. The likelihood function was assumed to be $e^{-\chi^2/2}$ with

$$ \chi^2 = \frac{\sum_i (I_{\text{obs},i} - I_{\text{mod},i})}{\sigma_i}, \quad (2) $$

where $I_{\text{obs},i}$, $I_{\text{mod},i}$, and $\sigma_i$ are the observed integrated intensity ratio, the modeled integrated intensity, and the $1\sigma$ uncertainty in $I_{\text{obs},i}$, respectively. The MCMC simulations were conducted using 20 walkers for 4500 steps following a burn-in period to ensure that the walkers achieved convergence. The $1\sigma$ uncertainties of the physical parameters were derived from the 16th and 84th percentiles of the posterior distribution.

Figure C.1 shows an example of the fitting results for GMC6. The plot demonstrates that H₂ and HCNH⁺ column densities are well constrained, whereas the gas temperatures are less precisely determined. The large uncertainties associated with the gas temperatures suggest that HCNH⁺ transitions are relatively insensitive to variations in gas temperature within the parameter space considered in the models. The complete fitting results for the eight GMCs are summarized in Table 2.

The beam-averaged HCNH⁺ column densities of (0.3–1.2) × 10¹⁵ cm⁻², derived from the best-fit RADEX modeling, were about a factor of three lower than those obtained from the LTE approach. This discrepancy suggested that the HCNH⁺ gas deviated from LTE conditions, despite its low critical densities. Such behavior indicates that HCNH⁺ predominantly arises from low-density regions, where H₂ number densities are below $10^4$ cm⁻³ (Bop et al. 2024). Therefore, we considered the HCNH⁺ column densities and abundances derived from our RADEX modeling to be more reliable. Even with the low HCNH⁺ column densities derived from non-LTE modeling, the HCNH⁺ abundances were revised to (1.8–9.4)×10⁻⁹ in the CMZ of NGC 253, most of which remained higher than all reported Galactic HCNH⁺ abundances (≤2 × 10⁻⁹) (e.g., Ziurys & Turner 1986; Schilke et al. 1991; Fontani et al. 2021; Gong et al. 2023).

Given the low dipole moment of HCNH⁺ (0.29 D; Botschwina 1986), the critical densities of its transitions are relatively low (see Table 1 in Gong et al. 2023). Utilizing the new collisional rate coefficients, we calculated the H₂ critical densities of the four HCNH⁺ transitions at high kinetic temperatures using Eq. (5) in Shirley (2015). At $T_K$ = 100 K, the critical densities for HCNH⁺ (2–1), HCNH⁺ (3–2), HCNH⁺ (4–3), HCNH⁺ (5–4) in the optically thin case are $1.0 \times 10^3$ cm⁻³, $3.6 \times 10^3$ cm⁻³, $8.6 \times 10^3$ cm⁻³, and $1.7 \times 10^4$ cm⁻³, respectively, which are much lower than those of the transitions of other molecules like HCN, HNC, CN, HCO⁺, and N₂H⁺ with comparable J quantum numbers (see Table 1 in Shirley 2015). When $T_K$ > 100 K, the critical densities decrease slightly. These critical densities indicate that they trace low-density gas with $n_{H_2} \le 10^4$ cm⁻³. Our non-LTE modeling analysis indeed reveals that the derived H₂ number densities are ~$10^{2.80}–10^{3.55}$ cm⁻³ (see Table 2), aligning with the values of the low-density component in the CMZ of NGC 253 (e.g., Tanaka et al. 2024). Furthermore, the H₂ number densities of the inner GMCs (i.e., GMC4 and GMC6) exceed those in the outer regions, consistent with the trend of the low-density component characterized in Fig. 5 of Tanaka et al. (2024). This agreement confirms that HCNH⁺ primarily traces the low-density regions of GMCs, highlighting its potential as a reliable tracer of such environments.

# 4. Discussion

In this study, we derived surprisingly high HCNH⁺ abundances in the CMZ of NGC 253, exceeding all reported values in the Milky Way. This stark contrast suggests that HCNH⁺ can be more abundant than previously thought, which underscores the environmental influences on its abundance. Compared to previous studies in the Milky Way, the extreme physical conditions in the CMZ of NGC 253 might potentially enhance the formation of HCNH⁺. Therefore, we performed astrochemical calculations to study the environmental influences on HCNH⁺ abundances using models of static clouds, photodissociation regions (PDRs), and shocks.

## 4.1. Static clouds

We first used the gas-grain astrochemical code, Chemp1¹⁴ (Du 2021), for calculations of a static cloud, employing the UMIST RATE12 chemical network (McElroy et al. 2013) enhanced with updated chemical reactions as described by Gong et al. (2023). The initial elemental abundances are the same as in Table 3 of McElroy et al. (2013), who tested the network using a dark cloud model. Although their initial conditions may not precisely match those of the CMZ of NGC 253, the environmental influences on HCNH⁺ abundances are similar for different initial conditions. Therefore, our modeling results can provide an illustrative case for understanding the environmental influences on the trends in abundance variations.

Our investigation focused on the effects of various parameters, including gas temperatures, gas densities, visual extinctions, interstellar radiation, and cosmic ray ionization rates on HCNH⁺ abundances. In these models, dust temperatures were assumed to be equal to gas temperatures. Based on previous studies, we constrained gas temperatures to $T_K \le 500$ K (Mangum et al. 2019; Rico-Villas et al. 2020; Huang et al. 2023; Tanaka et al. 2024), gas densities to $n_{\text{gas}} = n_H + n_{H_2} \le 10^5$ cm⁻³ (see Sect. 3.3.2 and Tanaka et al. 2024), interstellar radiation fields to $G_0 \le 10^9$ in

<footnote n="13">¹³ https://emcee.readthedocs.io/en/stable/</footnote>

<footnote n="14">¹⁴ https://github.com/fjdu/chempl</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Environmental influence on HCNH⁺ abundances relative to H₂ as a function of time, calculated from the static models of Chempl (Du 2021). (a) Modeling results for different gas temperatures with nH₂ = 10³ cm⁻³, Aᵥ = 10, G₀=1, and ζ = ζ₀. (b) Modeling results for different H₂ number densities with Tₖ = 300 K, Aᵥ = 10, G₀=1, and ζ = ζ₀. (c) Modeling results for different visual extinctions with Tₖ = 300 K, nH₂ = 10³ cm⁻³, G₀=1, and ζ = ζ₀. (d) Modeling results for different interstellar radiation fields with Tₖ = 300 K, nH₂ = 10³ cm⁻³, Aᵥ = 20, and ζ = ζ₀. (e) Modeling results for different cosmic ray ionization rates with Tₖ = 300 K, nH₂ = 10³ cm⁻³, Aᵥ = 20, and G₀ = 10⁵. In all panels, the cyan-shaded regions indicate the observed HCNH⁺ abundances derived from our non-LTE analysis (see Table 2).</figcaption></figure>

Habing units (Pérez-Beaupuits et al. 2018), and cosmic-ray ionization rates to $\zeta \lesssim 10^5 \zeta_0$ (e.g., Harada et al. 2021; Holdship et al. 2022; Behrens et al. 2022, 2024) with $\zeta_0 = 1.36 \times 10^{-17} \text{ s}^{-1}$ (e.g., van der Tak & van Dishoeck 2000).

Figure 8 illustrates the modeling results under different physical conditions. In Fig. 8a, HCNH⁺ abundances at Tₖ = 100 K can be even lower than those at Tₖ = 10 K at $10^6$ yr, but HCNH⁺ abundances increase with increasing gas temperatures above Tₖ = 100 K at $\ge 10^6$ yr. This transition is related to the evaporation temperatures, Tₑvap, of HCNH⁺'s precursor molecules such as HCN, above which HCN can be effectively released from dust grains. The evaporation temperature, at which the concerned molecule is half in gas phase and half on grain surfaces, can be estimated with the following formula (e.g., Hollenbach et al. 2009; Du et al. 2012):

$$T_{\text{evap}} \approx E_D / \left\{ 60 + \ln \left[ \left( \frac{10^5 \text{ cm}^{-3}}{n_H} \right) \left( \frac{20 \text{ K}}{T_{\text{gas}}} \right)^{1/2} \left( \frac{m}{28 \text{ amu}} \right)^{1/2} \right] \right\} \quad (3)$$

where E_D is the binding energy barrier of the species on grain surfaces, n_H is the hydrogen number density, T_gas is the gas temperature, and m is the molecular mass of the species. HCN binding energies were estimated to vary within the range of 2057–7421 K, depending on the different surface models (Ferrero et al. 2020). Based on the results from experiments

(Bergner et al. 2022), we adopted HCN–HCN and HCN-H₂O binding energies of E_D = 3207 ± 197 K and E_D = 4192 ± 68 K for our calculations. For HCN, m = 27 amu. Therefore, we estimate the evaporation temperatures of about 50–70 K for HCN molecules. In the CMZ of NGC 253, the gas temperature ranges from $\le 50$K on larger scales to $\ge 300$ K on 10 pc scales (Mangum et al. 2019). On a scale of 27 pc (i.e., our beam size of 1''6), the gas temperatures are likely higher than the HCN evaporation temperatures (Mangum et al. 2019; Huang et al. 2023; Tanaka et al. 2024; Bouvier et al. 2024). These elevated gas temperatures can significantly enhance the formation of the gas-phase precursors of HCNH⁺ (i.e., HCN and HCN⁺), thereby increasing its abundance.

When varying nH₂ in Fig. 8b, the peak HCNH⁺ abundances remain relatively constant, with differences within an order of magnitude. HCNH⁺ reaches its peak abundance earlier at higher nH₂, as chemical timescales are inversely proportional to nH₂ for two-body reactions. In Fig. 8c, we find that HCNH⁺ abundances decrease with decreasing Aᵥ at Aᵥ $\le 5$ due to photodissociation. At Aᵥ $\ge 5$, the HCNH⁺ abundances are nearly independent of Aᵥ. In Fig. 8d, the model results also show no dependence on the interstellar radiation field at Aᵥ = 20 (see Table 2). The independence can be explained by the fact that regions with such high Aᵥ are well shielded from ultraviolet radiation.

Figure 8e shows that higher cosmic-ray ionization rates result in lower HCNH⁺ abundances, consistent with Fontani et al. (2021). This is because high cosmic-ray ionization rates can
efficiently ionize H₂, leading to an enhancement of free electrons. These electrons can subsequently destroy molecular ions, including HCNH⁺ and key precursors to its formation, via dissociative recombination reactions, thereby reducing the HCNH⁺ abundance.

As shown above, the increased H₂ number densities, visual extinctions, interstellar radiations, and cosmic-ray ionization rates do not significantly increase the HCNH⁺ abundances. Given the high gas temperatures in the CMZ of NGC 253 (Mangum et al. 2019; Rico-Villas et al. 2020; Huang et al. 2023; Tanaka et al. 2024), we suggest that the unexpectedly high HCNH⁺ abundances could result from chemical enhancement triggered by high gas temperatures. However, recent ALCHEMI studies show that the cosmic-ray ionization rates are very high, ≥10³ζ₀, in these GMCs (e.g., Harada et al. 2021; Holdship et al. 2022; Behrens et al. 2022, 2024). At such high cosmic-ray ionization rates, none of the static models can reproduce the high HCNH⁺ abundances observed at their gas temperatures (e.g., Tₖ ≤ 300 K; Mangum et al. 2019). In contrast, chemically related species such as HCN and HNC can be successfully modeled with these static models (Behrens et al. 2022, 2024). We also note that gas and dust may be not well coupled in low-density regions traced by HCNH⁺. In these regions, dust temperatures are expected to be lower than gas temperatures, as confirmed in the CMZ of our Galaxy (e.g., Ao et al. 2013; Ginsburg et al. 2016). Since gas evaporation from dust grain mantles is regulated by dust temperatures, lower dust temperatures would lead to even lower HCNH⁺ abundances. Nevertheless, our main conclusion that the static models cannot reproduce the observed high HCNH⁺ abundances remains unchanged.

## 4.2. PDRs

HCNH⁺ transitions are likely to trace low-density regions (see Sect. 3.3.2), indicating that they may arise from PDRs. To explore this scenario, we used the Chempl code for PDR model calculations (Du 2021). As shown in Sect. 4.1, high gas temperatures can effectively increase the HCNH⁺ abundances. In order to reproduce the observed HCNH⁺ abundances, we used Tₖ = 500 K and a gas density of 10³ cm⁻³ as our fiducial parameters in the PDR models where dust temperatures are approximated following Eq. (9.18) in Tielens (2005),

$$T_d^5 = 8.9 \times 10^{-11} v_0 G_0 \exp[-1.8 A_V] + 2.78^5 + 3.4 \times 10^{-2} (0.42 - \ln[3.5 \times 10^{-2} \tau_{100\mu m} T_0]) \tau_{100\mu m} T_0^6 \quad (4)$$

where v₀ = 3 × 10¹⁵, τ₁₀₀µm = 10⁻³, and T₀ = 12.2G₀¹/⁵. The interstellar radiation field strength at the PDR edge, G₀, and the cosmic-ray ionization rate, ζ, were the only parameters varied in our calculations.

The modeling results are shown in Figure 9. As the visual extinction, A_V, increases, HCNH⁺ abundances generally increase for a given G₀ and ζ. For different G₀, HCNH⁺ abundances are nearly constant at high A_V (i.e., A_V ≥ 15). This is because high interstellar radiation fields are likely to destroy HCNH⁺ at low A_V, making HCNH⁺ most abundant in well-shielded regions. However, it should be noted that high values of ζ significantly reduce HCNH⁺ abundances. While the observed HCNH⁺ abundances can be reproduced for ζ ≤ 10²ζ₀, PDR models fail to match the observed values at the expected high cosmic-ray ionization rates, ζ ≥ 10³ζ₀ (e.g., Harada et al. 2021; Holdship et al. 2022; Behrens et al. 2022, 2024).

<figure><img src="image_10.png" /><figcaption>Fig. 9. HCNH⁺ abundance as a function of visual extinctions at different interstellar radiation fields and cosmic ray ionization rates, calculated from the PDR models of Chempl (Du 2021). The corresponding values are shown in the legend. The cyan-shaded region indicates the observed HCNH⁺ abundances in this study.</figcaption></figure>

## 4.3. Shocks

In Sect. 4.1, high kinetic temperatures of Tₖ ≥ 300 K are required to produce the observed high HCNH⁺ abundances in environments associated with high cosmic-ray ionization rates (see Fig. 8e). Such high kinetic temperatures can be caused by shock heating (Huang et al. 2023). Furthermore, Huang et al. (2023) found that shocked gas, traced by SiO transitions, had low H₂ number densities of ≤10⁴ cm⁻³, similar to the density regime probed by HCNH⁺ (see Sect. 3.3.2). This led us to speculate that shocks in the CMZ may enhance HCNH⁺ abundances. To investigate the impact of shocks on HCNH⁺ abundances, we employed the UCLCHEM<fnref n="15" /> code (Holdship et al. 2017) to run shock models. In these models, we also augmented the UMIST RATE12 chemical network with C₂N⁺ + H₂ → HCNH⁺ + C, as this reaction may contribute to the formation of HCNH⁺ (Gong et al. 2023). This code involves two phases. In phase 1, we modeled a static cloud of isothermal gas with Tₖ = 10 K and a gas density of 10³ cm⁻³ in atomic/ionic form with no molecules, leading to a model of a molecular cloud used in phase 2. In phase 2, we adopted C-type shocks for our calculations, parameterized according to Jiménez-Serra et al. (2008), similar to the approach in Huang et al. (2023). In this approach, the preshock magnetic field strength was approximated using Eq. (63) of Draine et al. (1983).

In the shock models, we mainly tested with different cosmic ray ionization rates and shock velocities with a fixed preshock gas density of 10³ cm⁻³. The modeling results are presented in Fig. 10. In these shock models, the temperatures reached their maximum around the nearly identical timescale at 6 × 10³ yr. In this time evolution plot, an enhancement of HCNH⁺ abundances was evident around the timescales close to the time of the maximum temperature and the HCNH⁺ abundances were significantly lower than the observed values before and after the shock phase. Both cosmic ray ionization rates and shock velocities affected HCNH⁺ abundances. In these models, elevated cosmic ray ionization rates are necessary to reproduce the observed HCNH⁺ abundance enhancements. With ζ = ζ₀, shock models cannot account for the observed HCNH⁺ abundance. Only by increasing ζ can the observed HCNH⁺ levels be achieved. This is because cosmic ray ionization rates enhanced the abundances

<footnote n="15">¹⁵ https://uclchem.github.io/</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Physical structures and the evolution of HCNH⁺ abundances in C-type shocks, calculated from the UCLCHEM models (Holdship et al. 2017). Top: C-shock physical structures of our models for the three different shock velocities. The black and orange lines represent the density and temperature profiles of shock models. Bottom: HCNH⁺ abundances as a function of time for C-type shocks. The preshock gas density is fixed to 10³ cm⁻³ in the models. Different lines correspond to different cosmic ray ionization rates and shock velocities which are shown in the legend. The cyan-shaded region indicates the observed HCNH⁺ abundances in this study. GIF animations of the shock models at a given ζ = 10³ζ₀ and ζ = 10⁵ζ₀ are available via link1 and link2, respectively.</figcaption></figure>

of molecular ions, which in turn led to the HCNH⁺ enhancement. However, when ζ ≥ 10³ζ₀, the modeled HCNH⁺ abundances decrease. This occurs because high cosmic ray ionization rates efficiently ionize H₂ and produce the elevated electron abundances, which in turn efficiently destroy molecular ions. In contrast, higher shock velocities led to significantly elevated HCNH⁺ abundances and sustained these enhanced levels over longer timescales.

We also compared the HCN and HNC abundances derived from Behrens et al. (2024) with our shock models. To reproduce the observed HCN and HNC abundances, higher cosmic ray ionization rates of ≥10³ζ₀ were required, consistent with the findings of previous studies (Behrens et al. 2022, 2024). However, it is important to note that the densities probed by HCN and HNC (Behrens et al. 2024) appeared to be higher than those traced by HCNH⁺, suggesting that HCNH⁺ might trace a distinct component of the molecular gas in the CMZ of NGC 253 compared to HCN and HNC.

While our static (Section 4.1) and PDR (Section 4.2) models failed to reproduce the observed HCNH⁺ abundances at high cosmic ray ionization rates of ζ ~ 10⁵ζ₀, our C-type

shock models could reproduce the high HCNH⁺ abundances at such extreme physical conditions with shock velocities of v_s ≥10 km s⁻¹. However, we note that a single shock event cannot sustain high HCNH⁺ abundances due to the rapid destruction of HCNH⁺ molecules in post-shock gas. Previous studies suggested that the CMZ might undergo multiple shock episodes, potentially triggered by outflows, star formation, and cloud-cloud collisions (Huang et al. 2023). These recurrent shocks could, in principle, sustain the elevated HCNH⁺ abundances observed. We therefore propose that HCNH⁺ primarily originates from shocked gas in the CMZ.

### 4.4. Formation and destruction paths

In the adopted chemical network, HCNH⁺ can be produced through ion-neutral reactions as follows:

$$H_3^+ + \text{HCN/HNC} \rightarrow \text{HCNH}^+ + \text{H}_2 \quad (5)$$

$$\text{HCO}^+ + \text{HCN/HNC} \rightarrow \text{HCNH}^+ + \text{CO} \quad (6)$$

$$\text{C}^+ + \text{NH}_3 \rightarrow \text{HCNH}^+ + \text{H} \quad (7)$$

$$\text{H}_2 + \text{HCN}^+/\text{HNC}^+ \rightarrow \text{HCNH}^+ + \text{H} \quad (8)$$

$$\text{C}_2\text{N}^+ + \text{H}_2 \rightarrow \text{HCNH}^+ + \text{C}, \quad (9)$$

while its destruction occurs via the following reactions,

$$\begin{align}
\text{HCNH}^+ + \text{e}^- &\rightarrow \text{HCN} + \text{H} \tag{10} \\
&\rightarrow \text{HNC} + \text{H} \nonumber \\
&\rightarrow \text{CN} + \text{H} + \text{H}. \nonumber
\end{align}$$

$$\text{CH} + \text{HCNH}^+ \rightarrow \text{HCN/HNC} + \text{CH}_2^+ \quad (11)$$

$$\text{CH}_2 + \text{HCNH}^+ \rightarrow \text{HCN/HNC} + \text{CH}_3^+ \quad (12)$$

$$\text{NH}_3 + \text{HCNH}^+ \rightarrow \text{HCN/HNC} + \text{NH}_4^+ \quad (13)$$

In cold regions ($T_K \sim 10$ K), reactions (5)–(7) are thought to dominate the formation of HCNH⁺, while dissociative recombination through reaction (10) serves as the main destruction pathway (e.g., Schilke et al. 1991). However, both HCN and HNC acted as reactants and products, presenting a “chicken-and-egg problem”. Chemical modeling by Fontani et al. (2021) suggested that reactions (7) and (8) were the primary formation pathways, with HCN⁺ mainly produced via CN + H₃⁺, effectively resolving the chicken-and-egg issue. Reaction (9) was previously overlooked because it was not included in the UMIST RATE12 chemical network. However, using the rate coefficient of reaction (9) (Knight et al. 1988; Loison et al. 2014), Gong et al. (2023) compared the chemical reaction rates of reactions (5)–(9) and suggested that reaction (9) could make non-negligible contributions. Fontani et al. (2021) also explored the chemistry of warm regions ($T_K \sim 27$ K), finding that reaction (6) became the primary formation mechanism, while the dominant destruction pathway remained unchanged.

We also explored the chemistry based on the shock models (see Sect. 4.3). The HCNH⁺ chemistry in the shock models is presented in Fig. 11. For instance, we employed a shock model
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Schematic diagram of the simplified chemical network that dominates the HCNH⁺ chemistry in the shock model.</figcaption></figure>

with a cosmic-ray ionization rate $\zeta = 10^5\zeta_0$ and a shock velocity of $v_s=30 \text{ km s}^{-1}$. In this model, the primary formation pathway for HCNH⁺ during the early phase of C-type shocks (at timescales of $\lesssim 1.4 \times 10^4$ years in Fig. 10) is reaction (8) involving $H_2+HCN^+$, similar to the cold model in Fontani et al. (2021). As the shock evolves into the second phase (at timescales of $(1.4-2) \times 10^4$ years in Fig. 10), reaction (9) became the dominant formation pathway, contributing to $\ge 80\%$ of the HCNH⁺ formation rate. Eventually, at later stages, reaction (8) once again took over as the primary formation mechanism. Additionally, reaction (5) via $H_3^+$ and HCN also significantly contributed to the formation of HCNH⁺ molecules during the early phase of shocks (at timescales of $5 \times 10^3$ to $1.4 \times 10^4$ years; Fig. 10), accounting for 16–27% of the HCNH⁺ formation rate during this period.

Another factor contributing to the elevated HCNH⁺ abundances is the temperature dependence of the destruction reaction (10). The rate coefficient for this reaction decreases with increasing gas temperature, following a power-law dependence of $-0.65$ (McElroy et al. 2013). At a gas temperature of 300 K, the reaction rate coefficient is about an order of magnitude lower than at 10 K. For demonstration, we once again employed the shock model with a cosmic-ray ionization rate $\zeta = 10^5\zeta_0$ and a shock velocity of $v_s=30 \text{ km s}^{-1}$. During the shock phase (on timescales of $2 \times 10^3$ to $3 \times 10^4$ years in Fig. 10), reactions (11)–(13) became increasingly efficient, but the overall destruction of HCNH⁺ was significantly reduced. Over the entire process, reaction (10) remains the dominant mechanism for the destruction of HCNH⁺. Thus, high temperatures induced by shock heating substantially mitigate the destruction of HCNH⁺, leading to an enhancement in its abundance.

**5. Summary and conclusions**

To explore the behavior and significance of HCNH⁺ under extreme conditions, we analyzed four rotational transitions of HCNH⁺ from the ALCHEMI large program, focusing on the physical and chemical properties of giant molecular clouds in the heart of the starburst galaxy NGC 253. Our main results are summarized as follows:

1. HCNH⁺ emission is widespread with a radial extent of $\sim 300$ pc in the central molecular zone, effectively tracing large-scale gas structures within the CMZ of NGC 253.

2. Our non-LTE analysis yields HCNH⁺ column densities of $(0.3–1.2) \times 10^{15} \text{ cm}^{-2}$, with corresponding HCNH⁺ abundances relative to $H_2$ of $(1.8–9.4) \times 10^{-9}$. Most of these HCNH⁺ abundances are higher than all reported values in the Milky Way. This highlights the elevated HCNH⁺ abundances in the CMZ of NGC 253. Furthermore, we observed a trend of decreasing HCNH⁺ abundances in GMCs as they approach the center of the CMZ.

3. Using the quantum mechanical coupled states' approximation, we calculated collisional rate coefficients for HCNH⁺ induced by *para*-H₂ and *ortho*-H₂ at kinetic temperatures up to 500 K. Incorporating these new coefficients into a non-LTE modeling framework and applying Monte Carlo Markov chain analysis, we suggest that HCNH⁺ emissions originate from [1]regions with H₂ number densities of $\sim 10^{2.80–10^{3.55}} \text{ cm}^{-3}$ with typical errors of $\sim 0.2$ dex, confirming HCNH⁺ as a tracer of low-density regions, even in nuclear starburst environments.

4. Astrochemical calculations assuming static, PDR, and shock environments have been performed. These modeling results indicate that the unexpectedly high HCNH⁺ abundances observed in the CMZ of NGC 253 can be attributed to chemical enhancement, driven by elevated gas temperatures and cosmic ray ionization rates in shocked, low-density gas.

The detection of HCNH⁺ in the CMZ of NGC 253 reveals significant differences compared to its presence in quiescent molecular clouds in the solar neighborhood. This study opens new avenues for understanding the complex chemistry and physical conditions in extreme environments of starburst galaxies. Given its role as a tracer of low-density regions, future investigations could benefit from mapping HCNH⁺ across various galactic environments to explore its potential as a diagnostic tool for studying various stages of molecular cloud evolution, the early stages of star formation, and the impacts of star formation feedback.

**Data availability**

Animated figures associated to Fig. 10 are available at https://gongyan2444.github.io/img/c-shock-zeta1e3.gif and https://gongyan2444.github.io/img/c-shock-zeta1e5.gif.

The HCNH⁺ data cubes and corresponding masks are available via https://zenodo.org/records/14959763.

The masking script is available at https://github.com/shbzhang/MaskingNGC253

*Acknowledgement.* We acknowledge the ALMA staff for their assistance with our observations. Y.G. is supported by the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDB0800301. C.H. acknowledges support by the Chinese Academy of Sciences President's International Fellowship Initiative under grant No. 2025PVA0048. We also acknowledge financial support from the European Research Council (Consolidator Grant COLLEXISM, Grant Agreement No. 811363) and the support from CEA/GENCI for granting us access to the TGCC/IRENE supercomputer under the A0110413001 project. F.J.D. is supported by National Key R&D Program of China grant 2023YFA1608000 and NSFC grant 12041305. S.V and M.B acknowledge the support from the European Research Council (ERC) Advanced Grant MOPPEX 833460. X.D.T. acknowledges the support of the National Key R&D Program of China under grant No. 2023YFA1608002, the Chinese Academy of Sciences (CAS) "Light of West China" Program under grant No. xbzg-zdsys-202212, the Tianshan Talent Program of Xinjiang Uygur Autonomous Region under grant No. 2022TSY-CLJ0005, and the Natural Science Foundation of Xinjiang Uygur Autonomous Region under grant No. 2022D01E06. W.Y. acknowledges the support from the National Natural Science Foundation of China (12403027), China Post-doctoral Science Foundation (2024M751376), and Jiangsu Funding Programme for Excellent Postdoctoral Talent (2024ZB347). This research has made use of spectroscopic and collisional data from the EMAA database (https://emaa.osug.frandhttps://dx.doi.org/10.17178/EMAA). EMAA is supported by the Observatoire des Sciences de l'Univers de Grenoble (OSUG). This research has made use of NASA's Astrophysics Data System. This work also made use of the Cube Analysis and Rendering Tool for Astronomy (CARTA) software
(Comrie et al. 2021) and Python libraries including Astropy (https://www.astropy.org/, Astropy Collaboration 2013), NumPy (https://www.numpy.org/, van der Walt et al. 2011), SciPy (https://www.scipy.org/, Jones et al. 2001), Matplotlib (https:// matplotlib.org/, Hunter 2007). This article makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00161.L and 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), NSTC 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. During the review process, we were deeply saddened by the passing of our esteemed colleague, Prof. Karl M. Menten, whose contributions to this work and to the field will always be remembered. We thank the anonymous referee for the helpful comments that improved this manuscript.

## References

Agúndez, M., & Wakelam, V. 2013, Chem. Rev., 113, 8710

Aladro, R., Martín, S., Riquelme, D., et al. 2015, A&A, 579, A101

Ao, Y., Henkel, C., Menten, K. M., et al. 2013, A&A, 550, A135

Arthurs, A., & Dalgarno, A. 1960, Proc. R. Soc. London. Ser. A. Math. Phys. Sci., 256, 540

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Bao, M., Harada, N., Kohno, K., et al. 2024, A&A, 687, A43

Behrens, E., Mangum, J. G., Holdship, J., et al. 2022, ApJ, 939, 119

Behrens, E., Mangum, J. G., Viti, S., et al. 2024, ApJ, 977, 38

Bendo, G. J., Beswick, R. J., D'Cruze, M. J., et al. 2015, MNRAS, 450, L80

Bergner, J. B., Rajappan, M., & Öberg, K. I. 2022, ApJ, 933, 206

Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013, Nature, 499, 450

Bop, C. T., & Lique, F. 2023, J. Chem. Phys., 158, 074304

Bop, C. T., Agúndez, M., Cernicharo, J., Lefloch, B., & Lique, F. 2024, A&A, 681, L19

Botschwina, P. 1986, Chem. Phys. Lett., 124, 382

Bouvier, M., Viti, S., Behrens, E., et al. 2024, A&A, 689, A64

Comrie, A., Wang, K.-S., Hsu, S.-C., et al. 2021, Astrophysics Source Code Library [record ascl:2109.031]

Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485

Du, F. 2021, Res. Astron. Astrophys., 21, 077

Du, F., Parise, B., & Bergman, P. 2012, A&A, 538, A91

Dubernet, M. L., Boursier, C., Denis-Alpizar, O., et al. 2024, A&A, 683, A40

Elitzur, M. 1992, ARA&A, 30, 75

Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. P. 2016, J. Mol. Spectr., 327, 95

Ferrero, S., Zamirri, L., Ceccarelli, C., et al. 2020, ApJ, 904, 11

Fontani, F., Colzi, L., Redaelli, E., Sipilä, O., & Caselli, P. 2021, A&A, 651, A94

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Ginsburg, A., Henkel, C., Ao, Y., et al. 2016, A&A, 586, A50

Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209

Gong, Y., Du, F. J., Henkel, C., et al. 2023, A&A, 679, A39

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Haasler, D., Rivilla, V. M., Martín, S., et al. 2022, A&A, 659, A158

Harada, N., Martín, S., Mangum, J. G., et al. 2021, ApJ, 923, 24

Harada, N., Martín, S., Mangum, J. G., et al. 2022, ApJ, 938, 80

Harada, N., Meier, D. S., Martín, S., et al. 2024, ApJS, 271, 38

Holdship, J., Viti, S., Jiménez-Serra, I., Makrymallis, A., & Priestley, F. 2017, AJ, 154, 38

Holdship, J., Mangum, J. G., Viti, S., et al. 2022, ApJ, 931, 89

Hollenbach, D., Kaufman, M. J., Bergin, E. A., & Melnick, G. J. 2009, ApJ, 690, 1497

Huang, K. Y., Viti, S., Holdship, J., et al. 2023, A&A, 675, A151

Humire, P. K., Henkel, C., Hernández-Gómez, A., et al. 2022, A&A, 663, A33

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Hutson, J., & Green, S. 1994, MOLSCAT computer code, version 14, in Collaborative computational project

Iodice, E., Arnaboldi, M., Rejkuba, M., et al. 2014, A&A, 567, A86

Jiménez-Serra, I., Caselli, P., Martín-Pintado, J., & Hartquist, T. W. 2008, A&A, 482, 549

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python

Karachentsev, I. D., Grebel, E. K., Sharina, M. E., et al. 2003, A&A, 404, 93

Kerkeni, B., & Simmie, J. M. 2023, J. Phys. Chem. A, 127, 5382

Knight, J. S., Petrie, S. A. H., Freeman, C. G., et al. 1988, J. Am. Chem. Soc., 110, 5286

Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2015, ApJ, 801, 25

Loison, J.-C., Wakelam, V., Hickson, K. M., Bergeat, A., & Mereau, R. 2014, MNRAS, 437, 930

López-Sepulcre, A., Balucani, N., Ceccarelli, C., et al. 2019, ACS Earth Space Chem., 3, 2122

Loreau, J., Lique, F., & Faure, A. 2018, ApJ, 853, L5

Mangum, J. G., Ginsburg, A. G., Henkel, C., et al. 2019, ApJ, 871, 170

Martín, S., Mauersberger, R., Martín-Pintado, J., Henkel, C., & García-Burillo, S. 2006, ApJS, 164, 450

Martín, S., Mangum, J. G., Harada, N., et al. 2021, A&A, 656, A46

McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36

McGuire, P., & Kouri, D. J. 1974, J. Chem. Phys., 60, 2488

Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163

Neufeld, D. A., Melnick, G. J., Sonnentrucker, P., et al. 2006, ApJ, 649, 816

Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111

Pérez-Beaupuits, J. P., Güsten, R., Harris, A., et al. 2018, ApJ, 860, 23

Phan, V. H. M., Peretti, E., Cristofari, P., Gusdorf, A., & Mertsch, P. 2024, MNRAS, 531, 2930

Pirlot Jankowiak, P., Lique, F., & Dagdigian, P. 2023, MNRAS, 523, 3732

Quénard, D., Vastel, C., Ceccarelli, C., et al. 2017, MNRAS, 470, 3194

Rekola, R., Richer, M. G., McCall, M. L., et al. 2005, MNRAS, 361, 330

Rico-Villas, F., Martín-Pintado, J., González-Alfonso, E., Martín, S., & Rivilla, V. M. 2020, MNRAS, 491, 4573

Sakamoto, K., Mao, R.-Q., Matsushita, S., et al. 2011, ApJ, 735, 19

Sanders, D. B., Mazzarella, J. M., Kim, D. C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 1607

Schilke, P., Walmsley, C. M., Millar, T. J., & Henkel, C. 1991, A&A, 247, 487

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369

Shirley, Y. L. 2015, PASP, 127, 299

Silva, W. G. D. P., Bonah, L., Schmid, P. C., Schlemmer, S., & Asvany, O. 2024, J. Chem. Phys., 160, 071101

Sobolev, V. V. 1960, Moving Envelopes of Stars (Cambridge: Harvard University Press)

Strickland, D. K., Heckman, T. M., Weaver, K. A., & Dahlem, M. 2000, AJ, 120, 2965

Tanaka, K., Mangum, J. G., Viti, S., et al. 2024, ApJ, 961, 18

Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge: Cambridge University Press)

Tritsis, A., Basu, S., & Federrath, C. 2023, MNRAS, 521, 5087

Turner, B. E. 1985, ApJ, 299, 312

Turner, J. L., & Ho, P. T. P. 1985, ApJ, 299, L77

van der Tak, F. F. S., & van Dishoeck, E. F. 2000, A&A, 358, L79

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22

van der Tak, F. F. S., Lique, F., Faure, A., Black, J. H., & van Dishoeck, E. F. 2020, Atoms, 8, 15

Ziurys, L. M., & Turner, B. E. 1986, ApJ, 302, L31
# Appendix A: Assessing errors from coupled state rate coefficients in radiative transfer calculations

<figure><img src="image_13.png" /></figure>

**Fig. A.1.** Dependence of radiation temperature ratios on gas density for the 1 → 0, 2 → 1, 3 → 2, and 4 → 3 emission lines. Superscripts CC and CS denote that the rate coefficients used in the radiative transfer calculations were obtained using the close-coupling (CC) method and the coupled states' (CS) approximation, respectively. The dashed line indicates perfect agreement between the methods. All calculations are performed under the LVG formalism as implemented in the RADEX code, assuming a FWHM line width of 1.0 km s⁻¹ and an HCNH⁺ column density of 2.3 × 10¹³ cm⁻² (van der Tak et al. 2007).

# Appendix B: Systematic comparison of HCNH⁺ rate coefficients induced by collisions with *o*- and *p*-H₂

<figure><img src="image_14.png" /></figure>

**Fig. B.1.** Comparison of HCNH⁺ rate coefficients induced by collisions with *o*- and *p*-H₂ for selected temperatures. The x- and y-axes represent the *o*- and *p*-H₂ scattering data, respectively. The dashed lines indicate an agreement within a factor of 2.0.
# Appendix C: Modeling the physical properties of GMC6 using the Monte Carlo Markov Chain calculations

<figure><img src="image_15.png" /><figcaption>**Fig. C.1.** Posterior probability distributions of gas temperatures, H₂ number densities, and HCNH⁺ column densities for GMC6, with the maximum posterior possibility point in the parameter space highlighted by orange lines and points. Contours represent the 0.5, 1.0, 1.5, and 2.0σ confidence intervals. The vertical dashed lines represent the 1σ spread.</figcaption></figure>