LETTER TO THE EDITOR

# HCNH⁺ abundance in cold dense clouds based on the first hyperfine resolved rate coefficients*

Cheikh T. Bop<fnref n="1" />, Marcelino Agúndez<fnref n="2" />, Jose Cernicharo<fnref n="2" />, Bertrand Lefloch<fnref n="3" />, and François Lique<fnref n="1" />

Received 13 December 2023 / Accepted 4 January 2024

## ABSTRACT

The protonated form of hydrogen cyanide, HCNH⁺, holds significant importance in astrochemistry, serving as an intermediate species in ion-neutral reactions occurring in the cold molecular clouds. Although it plays a crucial role in the chemistry of HCN and HNC, the excitation rate coefficients of this molecular cation by the dominant interstellar colliders have not been thoroughly investigated, leading to limitations in the radiative transfer models used to derive its abundance. In this work, we present the first hyperfine-resolved excitation rate coefficients for HCNH⁺ induced by collisions with both He and H₂ at low temperatures, addressing a crucial requirement for precise modeling of HCNH⁺ abundance in typical cold dense molecular clouds. Using non-local thermodynamic equilibrium (non-LTE) radiative transfer calculations, we reproduced the 1 → 0 and 2 → 1 observational spectra of HCNH⁺ fairly well and derived updated molecular column densities. For the TMC-1 molecular cloud, the new HCNH⁺ abundance is twice as large as suggested by previous LTE modeling, whereas the column density of this molecular cation is improved only by 10% in the case of the L483 proto-star. The factor of two in the case of TMC-1 most likely arises from an error in the early analysis of observational spectra rather than an effect of the LTE assumption, given that the HCNH⁺ lines are predominantly thermalized at densities higher than 2 × 10⁴ cm⁻<fnref n="3" />. For multiline studies of clouds of moderate densities, we strongly recommend using the collisional rate coefficients reported in this work.

**Key words.** molecular data – molecular processes – radiative transfer – scattering – ISM: abundances – ISM: molecules

## 1. Introduction

Protonated hydrogen cyanide, also known as iminomethylium (HCNH⁺), is the simplest protonated nitrile. It has a dipole moment of ~0.29 D (Botschwina 1986), large enough to allow its detection in radio astronomy. The first detection of this molecular cation in space took place toward Sgr B2, thanks to the observation of its three lowest rotational lines (Ziurys & Turner 1986). Since then, it has been observed in several cold star-forming regions such as the TMC-1 dark cloud (Schilke et al. 1991; Ziurys et al. 1992), the DR 21(OH) H II region (Hezareh et al. 2008), the L1544 pre-stellar core (Quénard et al. 2017), the L483 proto-star (Agúndez et al. 2022), and in 16 high-mass star-forming cores (Fontani et al. 2021). Recently, Gong et al. (2023) reported a comprehensive distribution analysis of HCNH⁺ within the Serpens filament and Serpens South, suggesting that this molecular cation is abundant in cold quiescent regions and deficient toward active star-forming regions. These observations present HCNH⁺ as a ubiquitous molecular cation in the cold interstellar medium (ISM) and reinforce the interest in understanding its chemistry.

HCNH⁺ is classified among the most important molecular ions in the ISM since it is the precursor of the widespread HCN and HNC. Ionic compounds play a crucial role in interstellar chemistry, serving as indispensable intermediates in ion-neutral reactions that govern gas-phase chemistry within cold cores (Agúndez & Wakelam 2013). Extensive research has been conducted to explore the chemistry of HCNH⁺ within dense, cold regions. This molecular cation is mostly formed in low temperature regions (Loison et al. 2014) through the following reactions:

$$ \mathrm{HCN}^{+}/\mathrm{HCN}^{+} + \mathrm{H}_{2} \rightarrow \mathrm{HCNH}^{+} + \mathrm{H}, \qquad (1) $$

and it undergoes destruction via a dissociative recombination with electrons (Loison et al. 2014; Semaniak et al. 2001):

$$ \begin{aligned} \mathrm{HCNH}^{+} + \mathrm{e}^{-} & \rightarrow \mathrm{HCN} + \mathrm{H} \\ & \rightarrow \mathrm{HNC} + \mathrm{H} \\ & \rightarrow \mathrm{CN} + \mathrm{H} + \mathrm{H}. \end{aligned} \qquad (2) $$

Based on the reactions in Eqs. (1) and (2) for HCNH⁺ and similar reactions related to HC₃NH⁺, the chemical model of Quénard et al. (2017) was unable to simultaneously reproduce the observed abundance of these protonated molecules toward the cold L1544 pre-stellar core. The authors pointed out that their model potentially underproduces HCNH⁺. Fontani et al. (2021) included, in addition to the previous reactions, Eq. (3) in their

<footnote n="1">¹ Univ. Rennes, CNRS, IPR (Institut de Physique de Rennes) – UMR 6251, 35000 Rennes, France
e-mail: cheikhtidiane.bop@ucad.edu.sn; francois.lique@univ-rennes.fr</footnote>

<footnote n="2">² Instituto de Física Fundamental, CSIC, Calle Serrano 123, 28006 Madrid, Spain
e-mail: marcelino.agundez@csic.es</footnote>

<footnote n="3">³ Université de Bordeaux – CNRS Laboratoire d'Astrophysique de Bordeaux, 33600 Pessac, France</footnote>

<footnote>* Hyperfine resolved rate coefficients are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/681/L19</footnote>
model as part of the dominant HCNH⁺ formation paths for cold high-mass star-forming cores:

$$ \mathrm{NH_3} + \mathrm{C^+} \rightarrow \mathrm{HCNH^+} + \mathrm{H}. \quad (3) $$

They found that changing the initial conditions (the hydrogen column density) or key parameters (the cosmic ray ionization rate) of their chemical model does not lead to a better agreement with the observations. In both attempts, their prediction underestimates the observed HCNH⁺ abundance. In TMC-1, Agúndez et al. (2022) found that the [HCNH⁺]/[HCN]+[HNC] abundance ratio is underestimated by the chemical model. This finding was seen as a general trend for protonated-to-neutral abundance ratios. More recently, Gong et al. (2023) investigated the impact of the hydrogen number density and the two energy barriers available in the literature for reaction (4), which plays an important role in the HCNH⁺ chemistry:

$$ \mathrm{HNC} + \mathrm{O} \rightarrow \mathrm{NH} + \mathrm{CO}. \quad (4) $$

Although the authors successfully reproduced the observed abundance of HCNH⁺ in cold high-mass star-forming cores, it is worth noting the wide range of the reference data they used, [HCNH⁺]/[H₂] = 3×10⁻¹¹–10⁻⁹, which covers all the individual target sources. A selective comparison would suggest a chemical model that yields less HCNH⁺ than observed for five out of the nine studied sources.

The HCN and HNC molecules, which are two of the most important species found in star-forming regions, have been extensively studied both from a chemical perspective and in terms of rotational energy transfer. For example, the HCN and HNC abundance profiles derived by Daniel et al. (2013) through non-local thermodynamic equilibrium (non-LTE) analysis of observational spectra align well with the chemical predictions of Gérin et al. (2009). Part of this agreement between chemistry and observation can be attributed to the non-LTE modeling of observational spectra, which has been made possible by the availability of accurate HCN and HNC collisional rate coefficients (Ben Abdallah et al. 2012; Sarrasin et al. 2010). In the case of HCNH⁺, observed abundances have been determined using the LTE approximation, which can lead to discrepancies and consequently contribute to the disagreements with chemical models (Quéñard et al. 2017; Fontani et al. 2021; Agúndez et al. 2022; Gong et al. 2023) discussed above. Therefore, it is of high interest to reevaluate the observed HCNH⁺ abundances employing a non-LTE approach, complemented with high accurate collisional rate coefficients.

In the frame of rotational energy transfer, the excitation of HCNH⁺ was first studied by Nkem et al. (2014). They used helium as a projectile and reported rate coefficients for temperatures up to 300 K, considering transitions among the 11 low-lying rotational energy levels. Recently, Bop & Lique (2023) extended the rotational basis considered in the previous work to the 16 first energy levels and also reported the HCNH⁺ scattering data due to collisions with H₂, the most abundant species in the ISM. They demonstrated that both ortho-H₂(j₂ = 1) and para-H₂(j₂ = 0), where j₂ represents the rotational quantum number of H₂, result in remarkably similar cross sections. This finding supports the idea that only rate coefficients induced by collisions with para-H₂(j₂ = 0) are necessary for modeling the abundance of molecular cations in cold star-forming regions.

We revisit the excitation of this molecular cation by para-H₂ taking into account the influence of the higher rotational energy level of the projectile (j₂ = 2), as previously done by Hernández Vera et al. (2017) for HCN and HNC. Furthermore,

we investigate the hyperfine splitting of the HCNH⁺ collisional rate coefficients resulting from the nonzero nuclear spin of nitrogen, as this effect is clearly resolved in the 1 → 0 observational line spectrum (Ziurys et al. 1992; Quéñard et al. 2017).

The structure of this paper is as follows: Sect. 2 provides a concise overview of the scattering calculations. Section 3 is dedicated to the astrophysical modeling, while Sect. 4 presents the concluding remarks.

## 2. Collisional rate coefficients

Accurate potential energy surfaces (PESs; 4D PES for HCNH⁺-H₂ and 2D PES for HCNH⁺-He), which are computed using the explicitly correlated coupled cluster method with single, double, and triple non-iterative excitation [CCSD(T)-F12] (Knowles et al. 1993, 2000) in conjunction with the augmented-correlation consistent-polarized valence triple zeta Gaussian best set (aug-cc-PVTZ) (Dunning 1989), are available in the literature (Bop & Lique 2023). These authors computed cross sections for the 16 low-lying rotational energy levels of HCNH⁺ due to collisions with both He and para-H₂ (hereafter denoted as H₂) using the “exact” close-coupling quantum mechanical approach (Arthurs & Dalgarno 1960), implemented in the MOLSCAT scattering code (Hutson & Green 1994).

The rotational energy levels were calculated using the spectroscopic constant of H₂ [B₀ = 59.322 cm⁻¹ and D₀ = 0.047 cm⁻¹] and HCNH⁺ [B₀ = 1.2360 cm⁻¹ and D₀ = 1.6075 × 10⁻⁶ cm⁻¹] (Huber 2013; Amano et al. 2006). The calculations were performed from a total energy of 2.5 cm⁻¹ to 800 cm⁻¹, using a fine step size. To better treat the couplings between open and closed channels for convergence reasons, it was necessary to take into account the 31 low-lying rotational energy levels of HCNH⁺ (j₁ = 0–30) in the calculations. The hybrid log derivative-airy propagator was used to solve the coupled equations (Alexander & Manolopoulos 1987). The integration limits were adjusted automatically for each total angular momentum (J), and the switching point from the log derivative to the airy integrator was set to 16 a₀. The integration step was maintained below 0.2 a₀ by adjusting the STEP parameter depending on the collision energy. Here, we address two important aspects of HCNH⁺-H₂ collisions that have not been considered before: (i) the effect of the H₂ rotational basis, specifically the inclusion of j₂ = 2, and (ii) the HCNH⁺ hyperfine structure arising from the nonzero nuclear spin of nitrogen.

### 2.1. The impact of the H₂ rotational basis

The excitation of HCNH⁺ by collisions with H₂ has been investigated using only j₂ = 0 (Bop & Lique 2023). The authors roughly estimated the mean deviation of the cross sections due to the inclusion of j₂ = 2 in the calculations to ~20%. Considering the huge anisotropy of the HCNH⁺-H₂ potential energy surface and its deep global minimum of ~1426.6 cm⁻¹, revisiting the collisional excitation of HCNH⁺ by H₂ with the most accurate level of precision becomes a necessary endeavor.

In Fig. 1, we evaluate the influence of the H₂ rotational basis in the HCNH⁺ collisional cross sections for selected total energies. The incorporation of j₂ = 0–2 into the H₂ rotational basis results in deviations of up to a factor of three compared to the restriction to j₂ = 0, while using an extended basis (j₂ = 0–4) instead of j₂ = 0–2 leads to moderate improvements, that is, deviations less than a factor of 1.5. For a total energy of 250 cm⁻¹, the root mean square errors obtained when using j₂ = 0 and j₂ = 0–4 in comparison with the use of j₂ = 0–2
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Comparison of the HCNH⁺ cross sections for selected total energies. We note that $\sigma_{J=0-5}$ is the sum of partial cross sections over total angular momenta up to $J = 5$. The stars depict the impact of including $j_2 = 0-2$ in comparison to the constraint of the H₂ rotational basis to $j_2 = 0$, while the empty circles estimate the influence of a more exhaustive H₂ rotational manifold ($j_2 = 0-4$) with respect to $j_2 = 0-2$. The dashed diagonal lines delimit an agreement region of a factor of 1.5.</figcaption></figure>

are 65% and 9%, respectively. Therefore, restricting the H₂ rotational basis to the ground level is insufficient to derive accurate collisional rate coefficients and the inclusion of $j_2 = 0-4$ slightly improves the results obtained using $j_2 = 0-2$, but the computational cost increases by a factor of approximately ten. Therefore, all calculations are performed including $j_2 = 0-2$ in the H₂ rotational basis.

## 2.2. Hyperfine resolved rate coefficients

Fully exploiting the information embedded within hyperfine resolved observational spectra requires an explicit description of the hyperfine splitting in the collisional rate coefficients. In this work, we only consider the coupling between the HCNH⁺ rotation and the nitrogen nuclear spin ($I = 1$). This coupling results in a slight splitting of each HCNH⁺ rotational level into three hyperfine components, with the exception of the ground energy level, which remains unsplit. The hyperfine components are identified by a quantum number $F$ defined as $|I - j_1| \le F \le I + j_1$. We employed the nearly exact recoupling method (Alexander & Dagdigian 1985) in the scattering matrix, which produced the results presented in the previous section. In this manner, we computed hyperfine resolved rate coefficients for the HCNH⁺ 25 low-lying energy levels, ($j_1, F) \le (8, 9)$, at low temperatures ($T = 5-30$ K).

Figure 2 displays the HCNH⁺ hyperfine resolved rate coefficients obtained using both He and H₂ as collision partners at 10 K, the typical temperature of cold star-forming regions. The magnitude of the H₂-induced rates vary between ~$10^{-10}$ and ~$10^{-12}$ cm⁻³ s⁻¹, whereas the He-rates drop drastically down to ~$10^{-14}$ cm⁻³ s⁻¹. The hyperfine resolution does not alter the existing disparity between the He- and H₂-rate coefficients, as discussed by Bop & Lique (2023), for the rotational transitions. The new insight into this plot is the unveiled propensity rule, $\Delta j_1 = \Delta F$, which applies to both projectiles. The data presented in this section are available in electronic supplementary material via the CDS and they will be accessible through databases such as Basecol, LAMDA, and EMAA.

## 3. Modeling astronomical lines of HCNH⁺

The collisional rate coefficients calculated here can be applied to model the lines of HCNH⁺ in those astronomical sources where lines are narrow, so that the hyperfine structure can be resolved. In cold dense clouds, line widths are typically below 1 km s⁻¹ (e.g., Agúndez et al. 2023), and thus if observed with a good enough spectral resolution, the hyperfine structure of the low-$j_1$ lines of HCNH⁺ can be resolved.

Protonated HCN has been observed in different types of molecular clouds (e.g., Schilke et al. 1991). Here we focus on two cold dense clouds, TMC-1 and L483, where low-$j_1$ lines of HCNH⁺ have been recently observed with a high spectral resolution. In the case of TMC-1, the $j_1 = 1 \to 0$ and $j_1 = 2 \to 1$ lines have been observed with the 30 m telescope of the Institut de Radio-Astronomie Millimétrique (IRAM) with a spectral resolution of 49 kHz. The observations of the $j_1 = 1 \to 0$ line at 74.1 GHz are part of a 3 mm line survey (Marcelino et al. 2007; Cernicharo et al. 2012), while those of the $j_1 = 2 \to 1$ line at 148.2 GHz are part of the Astrochemical Surveys At IRAM (ASAI) program (Lefloch et al. 2018). In the case of L483, the $j_1 = 1 \to 0$ line at 74.1 GHz was observed with the IRAM 30 m telescope with a spectral resolution of 49 kHz during a 3 mm line survey of this cloud (Agúndez et al. 2019, 2021). The observed lines are shown in Fig. 3.

To model the lines of HCNH⁺, we carried out excitation and radiative transfer calculations under the large velocity gradient (LVG) formalism (Goldreich & Kwan 1974). The code used is similar to MADEX (Cernicharo et al. 2012). We implemented the rate coefficients calculated here for inelastic collisions of HCNH⁺ with H₂ and He, where the hyperfine structure of HCNH⁺ is taken into account. The adopted abundance of He relative to H₂ is 0.17, based on the cosmic abundance of helium, which implies that collisional excitation is dominated by H₂, with He playing a minor role. We adopted the physical conditions of TMC-1 and L483 from the study of Agúndez et al. (2023). For TMC-1 we adopted a gas kinetic temperature of 9 K and a volume density of H₂ of $1.0 \times 10^4$ cm⁻³, while for L483 the adopted gas temperature is 12 K and the H₂ volume density $5.6 \times 10^4$ cm⁻³. The adopted line width, $0.46$ km s⁻¹ for TMC-1 and $0.39$ km s⁻¹ for L483, was taken directly from the arithmetic mean of the values measured on the spectrum of the $j_1 = 1 \to 0$ line. We then varied the column density of HCNH⁺ until matching the velocity-integrated intensity of the observed lines. The calculated line profiles are compared to the observed ones in Fig. 3. The column densities derived for HCNH⁺ in TMC-1 and L483 are $4.2 \times 10^{13}$ cm⁻² and $2.4 \times 10^{13}$ cm⁻², respectively. Previous determinations of the column density based on LTE are $1.9 \times 10^{13}$ cm⁻² in TMC-1 (Schilke et al. 1991) and $2.7 \times 10^{13}$ cm⁻² in L483 (Agúndez et al. 2019). The values determined here differ with respect to the previous values by a factor of two for TMC-1 and by 10% for L483.

To understand the different improvements resulting from non-LTE modeling in comparison to the previous LTE-based abundances of HCNH⁺ obtained for TMC-1 and L483, we investigated the deviation of the brightness temperature ($T_B$) with respect to LTE as a function of the gas density. Figure 4 shows that the use of He as a collision partner tends to delay the thermalization of the lines, whereas employing H₂ as a projectile, which is the major gas component in cold dense clouds, suggests that LTE assumption is valid for densities larger than $2 \times 10^4$ cm⁻³. Furthermore, the $j_1 = 2 \to 1$ line from Schilke et al. (1991) is reported to have an intensity (in main beam temperature) of 0.48 K, which aligns well with our $T_B$ of
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** HCNH⁺ hyperfine resolved rate coefficients for the (8, 9) → (j'₁, F') transitions at 10 K. The blue (red) line stands for collisional data obtained using H₂ (He) as a projectile.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Lines of HCNH⁺ in TMC-1 (j₁ = 1 → 0 and j₁ = 2 → 1 in the left panels) and L483 (j₁ = 1 → 0 in the right panel). Black histograms correspond to the observed line profiles, while the red lines correspond to the synthetic line profile calculated with the LVG model (see text).</figcaption></figure>

0.50 K obtained by correcting the 0.40 K antenna temperature with the beam efficiency of the IRAM 30 m telescope. We thus reinterpreted the observations of HCNH⁺ by employing the LTE assumption. As depicted in Table 1, using this approximation to determine the HCNH⁺ column density introduces an error of approximately 5% in both TMC-1 and L483. The factor of two observed in the case of TMC-1, when comparing the LVG-based column density with the data reported by Schilke et al. (1991), is likely to result from the assumptions incorporated in their analysis rather than an effect of the adopted LTE approximation.

As discussed in the Introduction, the observed HCNH⁺ column densities are underestimated by predictions from chemical models. In the case of TMC-1, for example, Agúndez et al. (2022) found that their chemical model underestimates the [HCNH⁺]/([HCN]+[HNC]) abundance ratio by a factor of ten. Since we revealed that the observed HCNH⁺ column density in this region is actually twice as large, the protonated-to-neutral abundance ratio turns out to be a factor of 20 higher than predicted by the chemical model. Although this matter remains

unresolved, we dispel any doubts that could implicate the observations, clearly identifying the chemical models as the sole factor. We suggest for future modelings using the new experimental rate constants of reactions in Eq. (1) measured for temperatures down to 17 K (Dohnal et al. 2023). Employing the latter results in the chemical models for cold dense clouds is more reasonable than the early reaction rates which were measured at room temperature (Wakelam et al. 2012).

## 4. Conclusion

We computed the first hyperfine resolved rate coefficients of HCNH⁺ induced by collisions with He and H₂. We used the most accurate recoupling method based on nuclear spin-independent scattering matrices calculated by the mean of the close-coupling quantum mechanical approach. When employing H₂ as the collision partner, we considered the coupling with the first excited rotational energy level of para-H₂, thereby improving the previously available nuclear spin-free rotational rate coefficients.
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Density dependence of the HCNH⁺ brightness temperature ratio for the $j_1 = 1 \rightarrow 0, 2 \rightarrow 1$, and $3 \rightarrow 2$ lines. Solid and dashed lines were obtained using the collisional rate coefficients of H₂ and He, respectively. The dash-dotted line refers to thermalization. The line width was set to 1.0 km s⁻¹. Thanks to the optically thin regime, these ratios are valid for column densities lower than $5 \times 10^{13} \text{ cm}^{-2}$.</figcaption></figure>

<table><caption>**Table 1.** Column density of HCNH⁺ (in $10^{13} \text{ cm}^{-2}$) under the LTE assumption and the LVG formalism for TMC-1 and L483.</caption><thead><tr><th>Sources</th><th>LTE</th><th>LVG</th></tr></thead><tbody><tr><td>TMC-1</td><td>1.93 <sup>(a)</sup></td><td>4.40</td><td>4.20</td></tr><tr><td>L483</td><td>2.70 <sup>(b)</sup></td><td>2.50</td><td>2.40</td></tr></tbody></table>

**Notes.** (a) and (b) stand for Schilke et al. (1991) and Agúndez et al. (2019), respectively.

Based on the new rate coefficients, we modeled HCNH⁺ emission lines observed toward TMC-1 and L483 using non-LTE radiative transfer calculations under the LVG formalism. With column densities of $4.2 \times 10^{13} \text{ cm}^{-2}$ and $2.4 \times 10^{13} \text{ cm}^{-2}$ for TMC-1 and L483, respectively, the synthetic spectra reproduced the observed ones quite well. The updated HCNH⁺ abundances differ by a factor of two and by 10% compared to the data previously available in the literature for TMC-1 and L483, respectively. It is worth noting that the large discrepancy observed in the case of TMC-1 is more likely due to an error in the early analysis of the observational spectra rather than an effect of the LTE assumption. The actual difference in the HCNH⁺ column density derived using LTE and LVG is approximately 5% for both TMC-1 and L483. Therefore, we confirm that the use of LTE to model the abundance of HCNH⁺ in cold, dense regions is reasonable. However, we strongly recommend employing the

rate coefficients reported in this work for multiline analysis and for observations toward regions of moderate densities.

**Acknowledgements.** The authors acknowledge the European Research Council (ERC) for funding the COLLEXISM project No 811363, the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of Centre National de la Recherche Scientifique (CNRS)/Institut National des Sciences de l'Univers (INSU) with Institut de Chimie (INC)/Institut de Physique (INP) co-funded by Commissariat à l'Energie Atomique (CEA) and Centre National d'Etudes Spatiales (CNES). F.L. acknowledges the Institut Universitaire de France. M.A. and J.C. acknowledge funding support from Spanish Ministerio de Ciencia e Innovación through grants PID2019-107115GB-C21 and PID2019-106110GB-I00. This work made use of ASAI "Astrochemical Surveys At IRAM".

## References

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

Agúndez, M., Marcelino, N., Cernicharo, J., Roueff, E., & Tafalla, M. 2019, *A&A*, **625**, A147

Agúndez, M., Roueff, E., Cabezas, C., Cernicharo, J., & Marcelino, N. 2021, *A&A*, **649**, A171

Agúndez, M., Cabezas, C., Marcelino, N., et al. 2022, *A&A*, **659**, L9

Agúndez, M., Marcelino, N., Tercero, B., Jiménez-Serra, I., & Cernicharo, J. 2023, *A&A*, **677**, A106

Alexander, M. H., & Dagdigian, P. J. 1985, *J. Chem. Phys.*, **83**, 2191

Alexander, M. H., & Manolopoulos, D. E. 1987, *J. Chem. Phys.*, **86**, 2044

Amano, T., Hashimoto, K., & Hirao, T. 2006, *J. Mol. Struct.*, **795**, 190

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

Ben Abdallah, D., Najar, F., Jaidane, N., Dumouchel, F., & Lique, F. 2012, *MNRAS*, **419**, 2441

Bop, C., & Lique, F. 2023, *J. Chem. Phys.*, **158**, 7

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

Cernicharo, J., Marcelino, N., Roueff, E., et al. 2012, *ApJ*, **759**, L43

Daniel, F., Gérin, M., Roueff, E., et al. 2013, *A&A*, **560**, A3

Dohnal, P., Jusko, P., Jiménez-Redondo, M., & Caselli, P. 2023, *J. Chem. Phys.*, **158**, 24

Dunning, T. H., Jr 1989, *J. Chem. Phys.*, **90**, 1007

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

Goldreich, P., & Kwan, J. 1974, *ApJ*, **189**, 441

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

Gérin, M., Marcelino, N., Biver, N., et al. 2009, *A&A*, **498**, L9

Hernández Vera, M., Lique, F., Dumouchel, F., Hily-Blant, P., & Faure, A. 2017, *MNRAS*, **468**, 1084

Hezareh, T., Houde, M., McCoey, C., Vastel, C., & Peng, R. 2008, *ApJ*, **684**, 1221

Huber, K. 2013, *Molecular Spectra and Molecular Structure: IV. Constants of Diatomic Molecules* (US: Springer)

Hutson, J., & Green, S. 1994, *Collaborative Computational Project*

Knowles, P. J., Hampel, C., & Werner, H.-J. 1993, *J. Chem. Phys.*, **99**, 5219

Knowles, P. J., Hampel, C., & Werner, H.-J. 2000, *J. Chem. Phys.*, **112**, 3106

Lefloch, B., Bachiller, R., Ceccarelli, C., et al. 2018, *MNRAS*, **477**, 4792

Loison, J.-C., Wakelam, V., & Hickson, K. M. 2014, *MNRAS*, **443**, 398

Marcelino, N., Cernicharo, J., Agúndez, M., et al. 2007, *ApJ*, **665**, L127

Nkem, C., Hammami, K., Halalaw, I. Y., Owono Owono, L. C., & Jaidane, N.-E. 2014, *Astrophys. Space Sci.*, **349**, 171

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

Sarrasin, E., Abdallah, D. B., Wernli, M., et al. 2010, *MNRAS*, **404**, 518

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

Semaniak, J., Minaev, B., Derkatch, A., et al. 2001, *ApJS*, **135**, 275

Wakelam, V., Herbst, E., Loison, J.-C., et al. 2012, *ApJS*, **199**, 21

Ziurys, L., & Turner, B. 1986, *ApJ*, **302**, L31

Ziurys, L., Apponi, A., & Yoder, J. 1992, *ApJ*, **397**, L123