© A. Jacob et al. 2022

1. Introduction

Over the duration of its operation, from 2009–2013, the Herschel Space Observatory (HSO) enabled observations of the fundamental rotational transitions of a variety of molecular hydrides and hydride ions, several of them being newly discovered. Many of these transitions cannot be observed from the ground at all because these high frequency lines lie in parts of the sub-millimetre (sub-mm) and far-infrared (FIR) wavelength range that are blocked by absorption in the Earth’s atmosphere. One of the many highlights of the Herschel mission, and a real surprise, has been the fortuitous detection of the J = 1−0 and J = 2−1 transitions of argonium, ArH⁺, in emission towards the Crab Nebula by Barlow et al. (2013). Following it, Schilke et al. (2014) were able to successfully assign, previously unidentified absorption features near 617 GHz to ArH⁺, along the lines-of-sight (LOS) towards five star-forming regions. It turned out that, fortunately, the J = 1−0 transition of ArH⁺ lies at a wavelength that is accessible with ground-based telescopes at high mountain sites under exceptional weather conditions, and Jacob et al. (2020) were recently able to detect ArH⁺ towards seven more sight lines in the inner Galaxy using the Atacama Pathfinder Experiment (APEX) 12 m sub-mm telescope. All of these observations confirm, as first discussed by Schilke et al. (2014), that the ArH⁺ molecular ion exclusively probes diffuse atomic material and that it is ubiquitously present in the Milky Way. This raises questions on the existence and nature of ArH⁺ in extragalactic sources.

Towards external galaxies, both the ArH⁺J = 1−0 and J = 2−1 lines have remained undetected in single side band observations covering their corresponding line frequencies, which were carried out using the Herschel Spectral and Photometric Imaging REceiver (SPIRE, Griffin et al. 2010) on board the HSO. The apparent non-detection of these lines is likely caused by a combination of effects including smearing by the spectrometer which results in unresolved spectral line profiles alongside effects of blending from nearby lines (such as HCN-v2 (7–6) at 623.3635 GHz and H₂O 2_2,0–2_1,1 at 1228.79 GHz). In addition, the ringing noise introduced by uncertainties in the fitting and subtraction of strong lines in the SPIRE FTS spectra using sinc functions can affect the line profiles of the underlying weak absorption. It is for these reasons that ArH⁺ has remained undetected in observations of extragalactic sources other than PKS 1830−211. Therefore, very little is known about the nature and abundance of ArH⁺ outside of the Milky Way. To date, there exists only a single detection of ArH⁺ in extragalactic space, which was carried out by Müller et al. (2015). Using the Atacama Large Millimetre/sub-millimetre Array (ALMA, Wootten & Thompson 2009), these authors were able to detect the J = 1−0 transitions of ³⁶ArH⁺ and ³⁸ArH⁺ through the intermediate redshift z = 0.8858 foreground galaxy absorbing the continuum of the gravitational lens-magnified blazar, PKS 1830−211 along two different sight lines.

Primarily residing in atomic gas with molecular hydrogen fractions, f_H2, between 10⁻² and 10⁻⁴ (Schilke et al. 2014; Neufeld & Wolfire 2016; Jacob et al. 2020), the abundance of ArH⁺ is sensitive to the X-ray and cosmic-ray fluxes that permeate the surrounding media as its formation is initiated by the reaction between H₂ and atomic argon ionised by either X-rays and/or cosmic-ray particles. Therefore, observations of the ground state transitions of ArH⁺ provide a unique tool for probing atomic gas and estimating ionisation rates. Regions permeated by a high flux of cosmic-rays can be heated by them to high gas temperatures, which in turn can strongly influence the initial conditions of star-formation and the initial mass function (IMF, Papadopoulos et al. 2011).

In this paper we present our search for ArH⁺ towards three luminous galaxies, Arp 220, NGC 253 and NGC 4945, using the SEPIA660 receiver on the APEX 12 m telescope. Both systems have been extensively studied over a wide range of wavelengths. In particular, a plethora of molecules have been observed towards these sources, including common hydrides and their cations, for example, OH, OH⁺, H₂O, H₂O⁺ (González-Alfonso et al. 2013, 2018; van der Tak et al. 2016). Arp 220 is the archetypical ultra-luminous infrared galaxy (ULIRG). A merging system, it hosts two compact nuclei (Baan & Haschick 1995; Rodríguez-Rico et al. 2005) that are surrounded by an immense amount of gas and dust (Scoville et al. 1997; Engel et al. 2011) with dust temperatures between 90 and 160 K (Sakamoto et al. 2008) and a luminosity of 0.2–1 × 10¹² L_⊙. Notably, the intense starburst activity within the dense interstellar medium (ISM) of its nuclear regions causes stars to form at a rate of up to 50–100 M_⊙ yr⁻¹ (Smith et al. 1998), which is > ​50 times that in the disk of the Milky Way galaxy today (Robitaille & Whitney 2010) and ∼​1000 times that in its central molecular zone (Immer et al. 2012). NGC 253 is a barred prototypical starburst galaxy part of the Sculptor group, with an infrared luminosity of 1.7 × 10¹⁰ L_⊙ (Radovich et al. 2001). Its strong nuclear starburst, drives a ∼​100 pc-scale molecular gas outflow/wind as seen, for example in observations of CO (Bolatto et al. 2013a). It has been suggested that at the centre of this barred spiral a weak AGN coexists with the strong starburst, an issue that is still under debate (for example see, Müller-Sánchez et al. 2010; Gutiérrez et al. 2020). NGC 4945 is an infrared-bright galaxy with a Seyfert nucleus, signifying an accreting supermassive black hole, in the Centaurus group with a luminosity of 2.4 × 10¹⁰ L_⊙ (Brock et al. 1988). It is the brightest Seyfert 2 galaxy and hosts a deeply enshrouded AGN at its centre which is revealed by X-ray emission in the 100-keV sky (Iwasawa et al. 1993). The AGN is surrounded by a strongly absorbing, inclined circumnuclear starburst ring with a radius of ∼​50 pc (Chou et al. 2007). With comparable star-formation rates of a few times M_⊙ yr⁻¹ (Bolatto et al. 2013a; Bendo et al. 2016) the similarities (or lack thereof) in the abundances and gas properties traced by ArH⁺ in these sources will shed light on their nuclear environments.

The observations are described in Sect. 2, followed by a qualitative and quantitative analysis of the data and a discussion of the results in Sects. 3 and 4. Finally, in Sect. 5 we discuss our main findings and summarise our results.

2. Observations

Using the Swedish-ESO PI (SEPIA660) receiver (Belitsky et al. 2018; Hesper et al. 2018) of the APEX 12 m sub-mm telescope, we were able to carry out observations of the J = 1−0 transition of ³⁶ArH⁺ (hereafter ArH⁺) between 2019 July and August (Project Id: M9519C_103). The SEPIA660 receiver is a two sideband (2SB), dual polarisation receiver that covers a bandwidth of 8 GHz, per sideband, with a sideband rejection level > 15 dB. The observations were carried out in wobbler switching mode, using a secondary wobbler throw of 180″ in azimuth at a rate of 1.5 Hz. Combined with the atmospheric stability at the high APEX site, this observing method allows reliable recovery of the sources’ continuum levels. Properties of our source sample are summarised in Table 1.

We tuned the upper sideband (USB) to a frequency of 618.5 GHz to cover the ArH⁺J = 1−0 transition at a rest frequency of 617.525 GHz. This allowed us to simultaneously observe the N_KaKc = 1₁₀ − 1₀₁, J = 3/2−1/2, and J = 3/2−3/2 transitions of p-H₂O⁺ at 604.678 and 607.227 GHz in the lower sideband (LSB), centred at a frequency of 606.5 GHz. The USB also covers an atmospheric absorption feature close to 620 GHz, however falls just short of covering the corresponding J = 1−0 transition of ³⁸ArH⁺ at 616.648 GHz. Our observations were carried out under excellent weather conditions, with precipitable water vapour (PWV) levels between 0.25 and 0.41 mm, corresponding to an atmospheric transmission better than or comparable to 0.5 in both sidebands and a mean system temperature of 1874 K, at 617 GHz. On average we spent a total (on+off) observing time of 4.6 h towards each source. The typical half power beam-width (HPBW) is $10\stackrel{″}{.}3$ at 617 GHz. The HPBW at this frequency corresponds 0.15 kpc and 0.19 kpc at the distances towards NGC 253 and NGC 4945, respectively. Therefore, the beam probes only the nucleus and foreground disk gas of these two galaxies. We converted the spectra into main-beam brightness temperature scales by using an antenna forward efficiency, F_eff, of 0.95 and a main-beam efficiency, B_eff, of 0.41 (determined from observations of Mars). As a backend we used an evolved version of the MPIfR built Fast Fourier Transform Spectrometer (XFFTS, Klein et al. 2012) which provided, over the entire 2 × 8 GHz bandwidth, a generic spectral resolution of 61 kHz (corresponding to 30 m s⁻¹ at 617 GHz). The calibrated spectra, smoothed to a velocity resolution appropriate for our sources (Δυ ∼ ​4.5 km s⁻¹), were subsequently processed using the GILDAS/CLASS software¹ and first order polynomial baselines were removed. Our search for ArH⁺ and p-H₂O⁺ towards Arp 220 was unsuccessful and we do not detect any lines in the sideband down to a noise level of 5 mK at a spectral resolution of 4.5 km s⁻¹. Therefore, we do not include Arp 220 in our analysis but quote a 3σ upper limit of 15 mK for the ArH⁺J = 1−0 and both p-H₂O⁺ transitions in this source. The non-detection of ArH⁺ towards Arp 220 might be because this system’s X-ray emission is not capable of ionising Ar or a combined effect of sensitivity and the low continuum level.

In addition to the results of the ArH⁺ and p-H₂O⁺ observations newly presented in this work, we use complementary archival HIFI/Herschel data of other well known tracers of atomic gas, namely, OH⁺ and o-H₂O⁺, which were acquired as a part of the Herschel EXtraGALactic (HEXGAL) guaranteed time key project (PI: R. Güsten), published by van der Tak et al. (2016) and also of the diffuse molecular gas tracer HF observed by Monje et al. (2014). The HIFI HPBWs at the frequencies of the OH⁺, o-H₂O⁺, and HF transitions studied in the above works are comparable to one another at 22″, 20″, and 17″, respectively.

The spectrum of the N, J = 1, 2−0, 1 OH⁺ transitions near 971 GHz towards NGC 253, observed under the HEXGAL project (presented by van der Tak et al. 2016) is saturated over almost the entire range of heliocentric velocities associated with the central parts of the galaxy, between 185 and 235 km s⁻¹. This makes the subsequent determination of OH⁺ column densities from the observed spectrum extremely difficult. Therefore, in this work, towards NGC 253 we instead use observations of the N, J = 1, 1−0, 1 transitions of OH⁺ near 1033 GHz observed in 2010 July using the now decommissioned dual-channel 1.05 THz receiver on the APEX 12 m telescope (Leinz et al. 2010). The observations were carried out under excellent weather conditions with PWV between 0.15 and 0.25 mm. The fast Fourier Transform spectrometer (FFTS, Klein et al. 2006) provided a spectral resolution of ∼​0.053 km s⁻¹ (183 kHz) over a 2.4 GHz bandwidth for the 1.05 THz channel which was later smoothed to a velocity resolution of 4.5 km s⁻¹. The spectra were calibrated on main-beam temperature scale using a main-beam efficiency of ( ∼ ​21 ± 7)% as determined from observations of Uranus. The observed spectrum shows a narrow absorption feature at 320 km s⁻¹ which likely arises from blending with contaminant species along the sight line. In order to confirm the association of this feature to a species other than OH⁺ we cross-checked the 1033 GHz OH⁺ spectrum along the LOS towards the extensively studied Galactic source, Sgr B2 (M), observed under the HEXOS Herschel guaranteed time key project (Bergin et al. 2010). From this comparison we assign this absorption feature near 1032.783 GHz as being likely associated with the ¹³CH₂CHCN (21_10,11–20_9,12) transition. In addition we also observe a weaker absorption feature near 12 km s⁻¹ which potentially arises from the high-lying (37_9,29–37_6,32) transition of C₂H₅CN at 1033.868 GHz. Furthermore the high excitation lines of similar complex organic molecules (COMs) found in ALMA Band 6 and 7 observations towards NGC 253 by Mangum et al. (2019) makes it likely that the OH⁺ spectrum at 1033 GHz is contaminated by the high-lying transitions of these species.

We also compare the line profiles of the ArH⁺ line with those of the HI 21 cm line and subsequently determine the ArH⁺ abundances. For this we use archival interferometric data of HI absorption and emission for NGC 253 and NGC 4945, observed using the Australia Telescope Compact Array (ATCA) with beam sizes of $4\stackrel{″}{.}9\times 10\stackrel{″}{.}3$ and $7\stackrel{″}{.}9\times 9\stackrel{″}{.}4$, respectively. The HI column densities were determined as described in Winkel et al. (2017), by combining the absorption profiles with emission line data. The results of this HI analysis, resulting in determinations of the optical depth, spin temperatures, and HI column densities along with the corresponding HI emission and absorption spectra is given in Appendix A. In addition we also present CO J = (1−0) emission spectra towards NGC 253 and NGC 4945, previously published in Houghton et al. (1997) and Curran et al. (2001), respectively, for comparison². In both cases the spectra were obtained using the Swedish–ESO 15 m Sub-millimetre Telescope (SEST) at La Silla, Chile with a beam size of 43″ at 115 GHz. The spectroscopic parameters of all the lines described above are summarised in Table 2.

We analysed the stability of the quoted continuum levels at 617 GHz across scans (or time) and found the fluctuations to lie within 14%. This is illustrated in Fig. B.1 and is not surprising as the observations were carried out using a wobbling secondary with a fast switching rate which guaranteed the removal of any drifts that may arise due to atmospheric instabilities. Furthermore, we find the 617 GHz continuum flux to be well correlated with the continuum flux at 870 μm observed using the Large APEX Bolometer Camera (LABOCA) at the APEX telescope presented in Weiß et al. (2008) which leads us to conclude that the continuum levels used are fairly reliable.

3. Results

3.1. Spectral line profiles

The calibrated and baseline-subtracted spectra towards NGC 253 and NGC 4945 are presented in Fig. 1. In the following paragraphs we qualitatively discuss the observed ArH⁺ and p-H₂O⁺ features and compare them to spectra of OH⁺, o-H₂O⁺, HF, CO, and HI. As mentioned above, our tuning setup simultaneously covers both the J = 3/2−1/2 as well as the J = 3/2−3/2 fine structure transitions from the (1_1,0 − 1_0,1) level of p-H₂O⁺ at 604.684 and 607.225 GHz. However, towards both NGC 253 and NGC 4945, we do not detect the J = 3/2−1/2 fine structure transition near 604 GHz above a noise level of ∼3 and 7 mK at a spectral resolution of 4.5 km s⁻¹, respectively. Therefore this transition is not discussed any further in the text. For both NGC 253 and NGC 4945, we assume molecular source sizes of ∼20″, which were previously determined by Wang et al. (2004), Aladro et al. (2015), Pérez-Beaupuits et al. (2018) and references therein, using molecular emission maps of abundant species like CO. The beam sizes of all the species studied here with the exception of CO (see Table 2) are either smaller than or comparable to this source size.

Fig. 1. From top to bottom: normalised absorption spectra of ArH+, p-H2O+, o-H2O+, OH+ (at 971 and 1033 GHz), HI and HF as well as the CO (1–0) emission line spectrum for comparison, towards NGC 253 (left) and NGC 4945 (right), respectively. In the spectra for NGC 253, the dotted light blue curves display the Gaussian fit to the emission component. The dashed black line represents the individual absorption profiles after subtracting the Gaussian fit with their Wiener filter fits overlaid in red for all species except HI and CO. The relative intensities of the hyperfine structure (HFS) components of the o-H2O+ and OH+ transitions are shown in black above their respective spectra and the grey shaded regions display their HFS deconvolved spectra. The vertical dashed black lines mark the systemic velocity of NGC 253 and NGC 4945 at 240 km s−1 and 563 km s−1, respectively. The 971 GHz OH+ line profile towards NGC 253 is saturated at blueshifted velocities while the 1033 GHz OH+ spectrum is potentially contaminated by the 13CH2CHCN (2110,11–209,12) line at 1032.783 GHz (fit shown by the dark blue curve). Also marked here in blue is contamination from the C2H5CN (379,29–376,32) line at 1033.868 GHz. The 971 GHz OH+ spectrum towards NGC 4945 is contaminated by the CH3OH Jk = 94,6 − 8,53 E line near 959.900 GHz originating from the image sideband (marked in blue).

3.1.1. NGC 253

The ArH⁺ line profile towards NGC 253 displays blueshifted absorption with respect to the systemic velocity of the source at 243 km s⁻¹, covering a velocity range from ∼​90 to 270 km s⁻¹ and centred at a velocity of 210 (±3.7) km s⁻¹. For comparison, other species like OH⁺, o-H₂O⁺ and HF all show spectra with P Cygni profiles with the absorption seen at comparable blueshifted velocities. The corresponding HI profile shows an asymmetric absorption component, centred at a blueshifted velocity of ∼​200 km s⁻¹. The velocity shift observed in the HI spectrum has been interpreted as evidence for a rotating nuclear disk of cold gas in this galaxy (Koribalski et al. 1995). The P Cygni profile observed towards the molecular ions studied here and HF, characterises the radial motion of the gas, which is indicative of outflows from the central region. The gas associated with the outflow component near 360 km s⁻¹ is unsurprisingly not traced by ArH⁺ and shows stronger emission in CO (particularly in its higher J-transitions) relative to the gas near 180 km s⁻¹. Previous studies have determined the outflow component to be kinematically distinct from the surrounding gas and is associated with the so-called western-superbubble located north-west of the central starburst region (Sakamoto et al. 2006) while the gas at 180 km s⁻¹ is co-located with the central disk (Krieger et al. 2019). Such an asymmetric emission is revealed by Hα observations whose emission is dominant on the approaching side of the outflow as reported by Bolatto et al. (2013a).

3.1.2. NGC 4945

Towards NGC 4945, the ArH⁺ line displays an asymmetric absorption profile between 455–715 km s⁻¹, similar to what is observed in the spectra of OH⁺, o-H₂O⁺ and HF lines. As discussed in Monje et al. (2014), van der Tak et al. (2016) and references therein, there are at least two velocity components, however unlike the other molecules and molecular ions, the broader component in ArH⁺ is centred at 510 km s⁻¹ and the narrower one at 605 km s⁻¹ while for the other species the broader component is the redshifted absorption feature.

In the p-H₂O⁺, o-H₂O⁺, OH⁺³, and HF spectra, the two absorption components are observed to be almost symmetric about the galaxy’s systemic velocity at 585 km s⁻¹ (Chou et al. 2007). This may indicate that the observed absorption arises from non-circular motions associated with the galaxy’s bar. A similar absorption dip is seen in the emission line profiles of lines from higher density gas tracers like HCN, HCO⁺, CN (Henkel et al. 1990), H₂CO (Gardner & Whiteoak 1974) and also from CO (Whiteoak et al. 1990) near 640 km s⁻¹. It likely traces foreground molecular gas that is moving towards the nucleus.

The HI absorption spectrum against the nuclear continuum of the source, also displays a similar profile with two absorption components with dips at ∼540 and ∼620 km s⁻¹. The observed shift in the centres of both the ArH⁺ absorption components in comparison to that of the other molecules, suggests that ArH⁺ does not trace the same layers of infalling molecular gas as the other molecules and molecular ions but rather traces mostly or exclusively atomic gas layers, as expected.

3.1.3. PKS 1830−211

As discussed in Sect. 1, the J = 1−0 transition of ArH⁺ was first detected in extragalactic environments by Müller et al. (2015) towards the intermediate redshift (z = 0.8858) lensing galaxy located in front of the blazar, PKS 1830−211. These authors study absorption spectra extracted from two separate lines-of-sight towards two magnified images located on the south-west (SW) and north-east (NE) sides of its nucleus. Owing to differences in the nature of the continuum between PKS 1830−211’s absorber and the sources studied here, a direct comparison of the spectral line profiles with those discussed in our study is not straightforward. While the background continuum for NGC 253 and NGC 4945 arise from the galaxies themselves, that for the sight lines studied towards foreground lensing galaxy in front of PKS 1830−211 arises from the distant quasar. Therefore, unlike the case for the nearby galaxies the spatial resolution for the latter is set by the small size of the background continuum emission (which is a few parsec in the plane of the z = 0.8858 galaxy) but nonetheless we briefly describe the properties of the ArH⁺ spectra observed towards both sight lines studied by Müller et al. (2015). The spectrum along the SW magnified image comprises of a single component with a line width of ∼57 km s⁻¹ and a weaker but broader blueshifted component, a feature that is also seen in the 607 GHz and 634 GHz p-H₂O⁺ spectra. However, the spectral line profile of ArH⁺ towards the NE image shows multiple narrow absorption features spanning 200 km s⁻¹ which is comparable to the spread over which absorption from other diffuse gas tracers is seen in this direction, namely CH⁺, HF, OH⁺ and H₂O⁺ (Muller et al. 2016, 2017).

3.2. Column densities

The line profiles can be expressed in terms of the optical depth, τ, using the radiative transfer equation, which for the particular case of absorption spectroscopy is given by T_b = T_ce^−τ, where T_b and T_c are the line, and the background continuum brightness temperatures, respectively. For NGC 253, we compute optical depths for the absorption profile obtained after subtracting the emission component. The emission component is modelled by fitting a Gaussian profile centred at the systemic velocity of the source. We fitted the optical depth profile on heliocentric velocity scales, that is to say τ versus. υ_helio, using the Wiener filter fitting technique as described in Jacob et al. (2019). This fitting procedure first fits the spectrum using the Wiener filter kernel by minimising the mean square error between the model (T_ce^−τ) and observations. For species like OH⁺ and o-H₂O⁺, whose rotational transitions further undergo hyperfine structure (HFS) splitting, this algorithm additionally deconvolves the HFS from the observed spectrum using the relative spectroscopic weights of the different HFS components. When fitting lines that do not exhibit HFS, the procedure simply assumes that there is only a single HFS component whose frequency corresponds to that of the fine-structure transition itself. Other than the observed spectrum and the spectroscopic parameters of the line to be fit, the only other input parameter required by the Wiener filter technique is the spectral noise, which is assumed to be independent of the observed signal. However, the Wiener filter faces singularities in portions of the spectrum in which the observed line profiles saturate or the line-to-continuum ratio tends to zero. This is the case for the OH⁺ spectrum towards NGC 253, which shows saturated absorption at blueshifted velocities between 185 and 235 km s⁻¹ (see Fig. 1). Therefore, as discussed in Sect. 2 we instead model the 1033 GHz transition of OH⁺. Prior to performing HFS deconvolution we remove contributions from the ¹³CH₂CHCN contamination using the WEEDS package in the GILDAS software.

The resulting optical depth profiles (τ_decon versus. υ_helio) are used to derive column densities as follows, assuming that the foreground absorption covers the background continuum source entirely:

$$\begin{array}{c}\hfill N_{\text{tot}}=\frac{8\pi {\nu }^3}{b_{\text{ff}}{c}^3}\frac{Q(T_{\text{ex}})}{g_{\text{u}}{A}_{\text{E}}}{\text{e}}^{E_{\text{u}}/T_{\text{ex}}}{[\text{exp}\left(\frac{h\nu }{k_{\text{B}}{T}_{\text{ex}}}\right)-1]}^{-1}{\displaystyle \int }{\tau }_{\mathrm{decon}}\text{d}\upsilon .\end{array}$$(1)

For a given species, all the spectroscopic terms in Eq. (1) namely, the upper level energy, E_u, the upper level degeneracy, g_u, and the Einstein A coefficient, A_E, all remain constant, except for the partition function, Q, which itself is a function of the rotation temperature, T_rot. Under conditions of local thermodynamic equilibrium (LTE), T_rot is equal to the excitation temperature, T_ex. The excitation of the molecules is straightforward as most of the particles are expected to occupy the ground state level, owing to the large Einstein A coefficients of all the lines involved, resulting in high critical densities of the order of a few 10⁷ cm⁻³ (see Table 2 for the Einstein A coefficients). Therefore, assuming a complete ground state occupation, we can approximate the excitation temperature to be less than the energy of the upper level above the ground and equal to the radiation temperature of the cosmic microwave background (T_CMB = 2.73 K). Similar assumptions for the excitation temperatures have been made by van der Tak et al. (2016) in their analysis. Since, T_CMB < T_ex < E_u/k_B, the column densities derived for the different species studied here strictly represent only lower limits. Unlike those of other species, the CO column densities are determined from the integrated intensities of their respective emission line profiles, for an excitation temperature of 20 K, as was assumed by Houghton et al. (1997) and Curran et al. (2001). Furthermore, the derived column densities are corrected to a first order for beam dilution effects, using the beam filling factor, b_ff, following $b_{\text{ff}}=\left[\left({\theta }_{\text{s}}^2+{\theta }_{\text{b}}^2\right)/{\theta }_{\text{s}}^2\right]$ where, θ_s, and θ_b represent the molecular source size and the beam size, respectively. The column density profiles hence derived per velocity channel are displayed in Fig. 2 and the total column densities derived by integrating between 55 and 295 km s⁻¹ for NGC 253, and 445 and 725 km s⁻¹ for NGC 4945, are summarised in Table 3. Using the column densities we determine from our spectra of the 607 GHz p-H₂O⁺ line, we predict, assuming conditions of LTE at 2.73 K, integrated intensities of 1.3 and 2.9 K km s⁻¹ for the 604 GHz transition of p-H₂O⁺, values that are lower than our 3σ upper limits of 2.1 and 5.2 K km s⁻¹ towards NGC 253 and NGC 4945, respectively.

Fig. 2. From top to bottom: column density distributions (black) of ArH+, p-H2O+, o-H2O+, OH+, HI, HF alongside the corresponding scaled H2 profile ([HF]/[H2] = 5×10−10; Emprechtinger et al. 2012) and CO towards NGC 253 (left) and NGC 4945 (right), respectively. The corresponding uncertainties are displayed by the blue shaded region.

3.3. Cosmic-ray ionisation rates

Cosmic-rays represent the dominant source of heating and ionisation in the inner parts of molecular clouds that are not penetrated by UV photons, making them an important driver for ion-molecular reactions, including those responsible for the formation of ArH⁺ (Roach & Kuntz 1970):

$$\begin{array}{c}\hfill \text{Ar}\stackrel{\text{CR}}{\to }{\text{Ar}}^{+}+{\mathrm{e}}^{-}\stackrel{{\text{H}}_2}{\to }{\text{ArH}}^{+}+\text{H}(\mathrm{\Delta }E=1.436\text{eV}).\end{array}$$

Hence, as a key ingredient in the ensuing chemistry, the abundance of ArH⁺ is sensitive to the cosmic-ray ionisation rate, ζ_p(H), and to the molecular fraction, f_H2, of the gas probed. The cosmic-ray ionisation rates in both NGC 253 and NGC 4945 have previously been determined by analysing the steady-state chemistry of oxygen-bearing ions like OH⁺ and H₂O⁺ by van der Tak et al. (2016). Following the steady-state analysis presented by these authors and Indriolo et al. (2015), we derive revised cosmic-ray ionisation rates by including contributions from p-H₂O⁺. The value of N_H used in our calculations is given by N(HI) + 2×N(H₂), where N(H₂) is obtained by using N(HF) as a surrogate for molecular hydrogen, where the abundance of HF in dense gas is assumed to be 5 × 10⁻¹⁰ (Emprechtinger et al. 2012). From the column density profiles presented in Fig. 2, it is clear that the total hydrogen content along the studied sight lines are dominated by dense molecular gas with pockets or bubbles of diffuse clouds as suggested by van der Tak et al. (2016). Including contributions from p-H₂O⁺ we derive mean cosmic-ray ionisation rates of 2.2 × 10⁻¹⁶ and 7.5 × 10⁻¹⁷ s⁻¹ across NGC 253 and NGC 4945, respectively, values that are not far from those derived by van der Tak et al. (2016). The derived cosmic-ray ionisation rates represent only lower limits because of uncertainties in the derived column densities, particularly in that of the HI and OH⁺ absorption profiles. In addition other assumptions made in this calculation, such as the values of the gas density (n_H = 35 cm⁻³), the electron fraction (x_e = 1.5 × 10⁻⁴) and the OH⁺ formation efficiency parameter (ϵ = 7%) are poorly constrained in these sources. The general impact of the uncertainties associated with the assumptions made in such an analysis, on the derived cosmic-ray ionisation rates of Galactic sources is discussed in more detail in Schilke et al. (2010) and Indriolo et al. (2015).

3.4. Gas properties

The cosmic-ray ionisation rates discussed in the previous section, ζ_p(H), describe only the total number of primary ionisations per hydrogen atom, per second. In contrast, the total cosmic-ray ionisation rate, ζ_t, includes contributions from both primary ionisation by cosmic-rays as well as ionisation by secondary electrons (resulting from the former). As discussed by Neufeld & Wolfire (2017), the exact relation between the two rates is dependent on the molecular fraction and fractional ionisation, but for the typical conditions of diffuse clouds the two are roughly related as ζ_p(H) = ζ_t/1.5. Therefore not only do cosmic-rays play an important role in heating the gas, they also increase the electron abundance.

In the following sections we evaluate the impact of cosmic-ray ionisation and heating through the physical properties traced by ArH⁺. We estimate the molecular fraction characteristic of the gas traced by ArH⁺ and explore the impact of collisional excitation by electrons on the observed ArH⁺ abundance.

3.4.1. Molecular fraction

Under the assumption that the cloud volumes containing ArH⁺ are exposed to the same cosmic-ray ionisation flux as those traced by both OH⁺ and H₂O⁺, we can derive the molecular fraction, f_H2, of the gas probed by ArH⁺, using the ArH⁺ abundances (with respect to HI) derived from observations as a constraint. We first present the relation between ζ_p(H), X(ArH⁺) and f_H2 by analysing the steady-state ion-molecular chemistry of ArH⁺ as discussed in Schilke et al. (2014). ArH⁺ is destroyed primarily via proton transfer reactions with H₂ or atomic oxygen, and photodissociation:

$$\begin{array}{cc}\hfill {\text{ArH}}^{+}& +{\text{H}}_2\to \text{Ar}+{\text{H}}_3^{+};k_5=8\times {10}^{-10}{\text{cm}}^3{\text{s}}^{-1}\hfill \\ \hfill {\text{ArH}}^{+}& +\text{O}\to \text{Ar}+{\text{OH}}^{+};k_6=8\times {10}^{-10}{\text{cm}}^3{\text{s}}^{-1}\hfill \\ \hfill {\text{ArH}}^{+}& +\text{h}\nu \to {\text{Ar}}^{+}+\text{H};k_7=1.0\times {10}^{-11}{\chi }_{\text{UV}}{f}_{\text{A}}{\text{s}}^{-1}.\hfill \end{array}$$

The photodissociation rate of ArH⁺ was estimated by Alekseyev et al. (2007) to be ∼ 1.0 × 10⁻¹¹f_A s⁻¹ for an unshielded cloud model that is uniformly surrounded by the standard Draine UV interstellar radiation field. The attenuation factor, f_A, is given by an exponential integral and is a function of visual extinction, A_v. For a cloud model with A_v = 0.3, Schilke et al. (2014) derived values for f_A between 0.30 and 0.56 that increase as you move outwards from the centre of the cloud. For our analysis, we use 0.43 for f_A, which lies mid-way through the computed range of values.

We further assume an atomic oxygen abundance (relative to H nuclei) of 3.9 × 10⁻⁴ (Cartledge et al. 2004) and an argon abundance close to this element’s solar abundance of 3.2 × 10⁻⁶ (Lodders 2008). Observations by van der Tak et al. (2016) suggest that the dense and diffuse gas phases are well mixed in the galaxies under consideration here, a notion based on the similarities between the observed profiles of OH⁺ and H₂O⁺ lines to those of H₂O and HI lines. Therefore, akin to these authors, who derived the cosmic-ray ionisation rates for NGC 253 and NGC 4945 by analysing the steady-state chemistry of OH⁺ and H₂O⁺, we assume a gas density, n(H), of 35 cm⁻³ (Indriolo et al. 2015) for the OH⁺–H₂O⁺ absorbing clouds in our analysis.

Using these values, the cosmic-ray ionisation rate can be approximated as

$$\begin{array}{cc}\hfill {\zeta }_{\text{p}}(\text{H})({\mathrm{s}}^{-1})& =\frac{N({\text{ArH}}^{+})}{N_{\text{H}}}\left(\frac{k_{5}n({\text{H}}_2)+k_{6}n(\text{O})+k_7}{11.42}\right)({\mathrm{s}}^{-1}),\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill & =\frac{N({\text{ArH}}^{+})}{N_{\text{H}}}\left(\frac{0.5005+448f_{{\text{H}}_2}}{1.2\times {10}^6}\right).\hfill \end{array}$$(3)

Equation (3) is re-arranged and expressed in terms of f_H2, as

$$\begin{array}{c}\hfill f_{{\text{H}}_2}=2.68\times {10}^3(\mathrm{s})[\frac{{\zeta }_{\text{p}}(\text{H})}{X({\text{ArH}}^{+})}-4.17\times {10}^{-7}]({\mathrm{s}}^{-1}).\end{array}$$(4)

Using Eq. (4), the molecular fractions of the gas probed by ArH⁺ are found to be a few times 10⁻³ towards both sources, values that are comparable to what is derived along the NE line-of-sight towards PKS 1830−211 (which amongst the two sight lines studied by Müller et al. 2015, has been shown to probe more diffuse and atomic environments, Koopmans & de Bruyn 2005) with X(ArH⁺) = 2.8 × 10⁻⁹ and ζ_p(H) = 3 × 10⁻¹⁵ s⁻¹. The molecular gas fraction hence derived, per velocity channel over velocity intervals most relevant to the sources studied (as discussed in Table 3) is displayed in Fig. 3. Towards both sources we find the ArH⁺ bearing clouds to host abundances between 10⁻¹¹ and 10⁻¹⁰ and trace molecular gas fractions between 3 × 10⁻⁴ and a few times 10⁻². In comparison, OH⁺, whose abundance is almost two orders of magnitude higher than that of ArH⁺, mainly traces gas with molecular fractions of the order of a few times 10⁻¹ whereas Muller et al. (2016) derived f_H2 values for OH⁺-bearing gas that varies between 0.01 and 0.07 towards both sight lines in their study. Therefore, the different atomic gas tracers may not be spatially co-existent. Furthermore, the transition from atomic to molecular gas is illustrated in the top panel of Fig. 3 through the distribution of the molecular gas fractions traced by ArH⁺-, OH⁺-, and CO-bearing gas volumes. The molecular fractions in this analysis are computed using Eq. (4), Eq. (12) of Indriolo et al. (2015) and by assuming a CO-to-H₂ conversion factor of 2 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹ as recommended by Bolatto et al. (2013b), albeit with a factor of 2 uncertainty. By combining previous estimates of the cosmic-ray ionisation rates derived for Arp 220, of ζ_p(H) > 10⁻¹³ s⁻¹ (González-Alfonso et al. 2013; van der Tak et al. 2016), with our upper limit for the ArH⁺ abundance of 2.5 × 10⁻¹¹ (determined by integrating the 3σ detection limit quoted in Sect. 2 over a line width of 743 km s⁻¹, Mirabel 1982 and N(H₂) = 1 × 10²⁴ cm⁻², Downes & Eckart 2007), we are unable to derive reasonable values for the molecular fraction traced by ArH⁺.

Fig. 3. Top: distribution of molecular gas fraction (fH2) traced by ArH+ (blue), OH+ (red) and CO (orange). Bottom: contours of X(ArH+) (blue) and X(OH+) (red) abundances with respect to NH in the fH2–ζp(H) plane. Blue triangles and red diamonds represent the corresponding values derived from the LOS observations presented in this work computed per channel in the velocity interval (Δυ ∼ ​4.5 km s−1) between 55 and 295 km s−1, and 445 and 725 km s−1 for NGC 253 and NGC 4945, respectively.

3.4.2. Excitation by electrons

In this section we explore the effects of collisional excitation of ArH⁺ by electrons. We expect ArH⁺ to reside in cloud layers with electron fractions, x_e of ≥10⁻⁵–10⁻⁴, making them competitive collision partners to atomic hydrogen in such cloud environments. Moreover, the destruction of ArH⁺ via dissociative recombination reactions with electrons have, alongside photodisocciation pathways, been shown to have small to negligible impacts on the ArH⁺ abundances in astrophysical environments (Alekseyev et al. 2007; Mitchell et al. 2005; Abdoulanziz et al. 2018), making electron-impact excitation effects of significant importance for ArH⁺. In order to explore this and to evaluate the validity of our assumptions for the electron density and excitation temperature, we performed non-LTE radiative transfer calculations using the statistical equilibrium radiative transfer code RADEX (van der Tak et al. 2007). The models were run for a uniform sphere geometry (with the offline version of RADEX, which allows a choice of input geometries) under the large velocity gradient (LVG) approximation. Based on the escape probability formalism, the code computes level populations, line intensities, excitation temperatures, and optical depths as a function of the physical conditions (kinetic temperature and density) and radiative transfer parameters (column density and line width) specified as inputs. The models were run using rate coefficients recently computed by Hamilton et al. (2016) for collisions between ArH⁺ and electrons as well as those between ArH⁺ and atomic hydrogen computed by Dagdigian (2018).

By constraining the models using the ArH⁺ column densities as determined in Sect. 3.2, we model the excitation temperature of the observed 617 GHz ArH⁺ line as a function of the electron density and gas temperature, T_kin. The models are run across a 100 × 100 grid, with x_e values between 10⁻⁶, and 10⁻¹ and gas temperatures between 10 and 3000 K. Furthermore, the total gas density, n_H, approximated by the sum of the electron and atomic gas densities in the models, is fixed to values of 35, 100 and 1000 cm⁻³. The modelled results are visualised in Fig. 4. Similar to Dagdigian (2018), we find the excitation temperature for the models with total gas densities fixed at 100 and 1000 cm⁻³ to be greater than our assumed value of T_ex at 2.73 K, across the entire x_e parameter space, with values of T_ex only slightly higher than 2.73 for x_e ⪅ 10⁻³, beyond which the value of the excitation temperature increases quite rapidly. While the general trend for those models with a total gas density of 35 cm⁻³ (equal to the gas densities assumed in Sect. 3.3 to compute the cosmic-ray ionisation rates) is similar to those with higher gas densities, that is, 100 and 1000 cm⁻³, the ArH⁺ excitation temperature over the assumed electron fraction (≈1.5 × 10⁻⁴ see Sect. 3.3) remains very close to a value of 2.73 K. This implies electron densities of ∼5.3 × 10⁻³, consistent with the low electron abundances expected in gas volumes bearing OH⁺, a species formed in cloud environments similar to those that contain ArH⁺. Low electron densities are prerequisites for the formation of detectable amounts of OH⁺, whose formation pathway via ${\text{H}}_3^{+}$ and atomic oxygen competes with the dissociative recombination of ${\text{H}}_3^{+}$ with electrons. Observationally however, the range of electron densities is poorly constrained with a possible upper limit set by assuming the electron densities to be consistent with the photoionisation of neutral carbon.

Fig. 4. RADEX modelled excitation temperature as a function of electron density, for fixed values of the total gas density at 35 (red), 100 (orange) and 1000 cm−3 (dark blue) for the ArH+ (1–0) transition. Furthermore, each model is run for fixed gas temperatures at 30 (solid), 300 (dashed) and 3000 K (dotted), respectively. The horizontal dashed grey line marks an excitation temperature of 2.73 K. The inset zooms in on 10−4 < xe < 10−2 with the vertical dashed purple line indicating xe = 1.5 × 10−4 which corresponds to the value of xe used in our calculations.

Using the same physical conditions as used to model the ArH⁺ 1−0 transition, we model the 2−1 transition of ArH⁺ near 1234.602 GHz. The results are displayed in Fig. 5. For all the models the excitation temperature for the 2−1 line is higher, between ∼8 and 18 K, while the optical depths are extremely low, of the order of a few 10⁻⁶ for total gas densities of 35 cm⁻³ and only an order of magnitude higher for models with total gas densities of 1000 cm⁻³. Therefore, as discussed by Dagdigian (2018) given the low optical depths reproduced by the models, it is highly unlikely that the 2−1 transition of ArH⁺ can be observed in interstellar gas; so far it has only been detected (in emission) in the extreme circumnebular environment of the Crab (Barlow et al. 2013; Priestley et al. 2017).

Fig. 5. RADEX modelled excitation temperatures (left) and optical depths (right) as a function of electron density, for fixed values of the total gas density at 35 (red), 100 (orange) and 1000 cm−3 (dark blue) for the ArH+ (2–1) transition. Furthermore, each model is run for fixed gas temperatures at 30 (solid), 300 (dashed) and 3000 K (dotted), respectively. The horizontal dashed grey line marks an excitation temperature of 2.73 K (in the left-hand panel) while the vertical dashed purple line indicates, xe = 1.5 × 10−4.

3.5. H₂ O⁺ ortho-to-para ratio

The H₂O⁺ molecular ion exists in two symmetric states, o- and p-H₂O⁺, that have opposing parities due to the interaction between the magnetic moment of the unpaired electron and protons. Studying the ratio of molecules in these two states which is reflected by the spin temperature, can provide insight into the formation pathway and thermodynamic properties of the gas. In the following paragraphs we determine the ortho-to-para ratio (OPR) of H₂O⁺ and derive the nuclear spin temperature.

Unlike for H₂O, the lowest energy state of H₂O⁺ corresponds to the lowest rotational energy level of its ortho state. This is a result of the molecule’s C_2v symmetry and ²B₁ ground state configuration. The fine structure levels of the o-H₂O⁺ spin state, with a nuclear spin, I, of 1, further undergo HFS splitting while, p-H₂O⁺ (with I = 0), does not. The lowest rotational energy level of the ortho spin state, N_KaKc = 1_1,1 − 0_0,0, lies at an energy level which is 30.1 K lower energy than that of the lowest p-H₂O⁺ level, 1_1,0 − 1_0,1. Under the assumption that the rotational temperature is close to T_CMB = 2.73 K, one would expect that most of the ions occupy either one of the H₂O⁺ ground state levels. Similar to H₂ and H₂O, H₂O⁺ is expected to exhibit an OPR of at least 3:1 (Townes & Schawlow 1975). Typically, the conversion between the ortho and para states of H₂O⁺ occurs via gas phase reactions with atomic hydrogen, as follows:

$$\begin{array}{c}\hfill {\text{p-H}}_2{\text{O}}^{+}+\text{H}\rightleftharpoons {\text{o-H}}_2{\text{O}}^{+}+\text{H}\end{array}$$(5)

while H₂O⁺ energetically reacts with molecular hydrogen to yield H₃O⁺.

Moderately coupled via collisions, the OPR between the two spin states of H₂O⁺ is given by,

$$\begin{array}{c}\hfill \text{OPR}\equiv \frac{Q_{\text{ortho}}}{Q_{\text{para}}}\text{exp}(-\mathrm{\Delta }E/T_{\text{ns}}),\end{array}$$(6)

where Q_ortho and Q_para represent the partition functions of the respective spin isomers, ΔE is the energy difference between them, ΔE = −30.1 K and T_ns is the nuclear spin temperature. The value of ΔE is expressed as a negative quantity because, as discussed above, the lowest ortho state has a lower energy than the lowest para state. As discussed in Appendix A of Schilke et al. (2010), at low temperatures, the partition functions of the two states are governed by the degeneracy of their lowest fine-structure and HFS levels. Having the same quantum numbers and upper level degeneracies, the ratio of the partition functions approaches unity as the rotational temperature tends to 0 K.

Figures 6 and 7 display the distribution of the derived column densities per velocity intervals (corresponding to the spacing of the velocity channel bins) for o- and p-H₂O⁺, the OPR and the nuclear spin temperature. The OPRs determined for NGC 253 are quite large, with a mean value of roughly 7.0 ± 3.1, whereas, the mean value derived towards NGC 4945 is 3.5 ± 1.2. The derived OPRs are consistent with the equilibrium value of three within the quoted error bars and correspond to nuclear spin temperatures between 15 and 24 K.

Fig. 6. Top panel: column density per velocity interval distribution of o-H2O+ (solid purple) and p-H2O+ (dashed blue) over the entire velocity interval covering absorption towards NGC 253. Middle panel: OPR H2O+ distribution (left) and corresponding histogram (right). Bottom panel: nuclear spin temperature distribution (left) and corresponding histogram (right). The median and 1 σ levels of the OPR and nuclear spin temperature are marked by solid and dotted black lines, respectively.

Fig. 7. Same as Fig. 6, but towards NGC 4945. The median value for Tns derived from the distribution is the same as that obtained from deriving Tns using the median OPR for H2O+.

The very large values found for NGC 253 for LSR velocities > ​150 km s⁻¹ can be partly attributed to the uncertainties resulting from modelling the emission component observed in the P-Cygni profiles of its o- and p-H₂O⁺ spectra. Theoretical studies carried out by Tanaka et al. (2013) on the ortho-to-para exchange reactions of H₂O⁺ reveal quite low reaction rates, which implies that the OPR of H₂O⁺ does not significantly deviate from the equilibrium value of three, particularly at low temperatures, lower than the values we find for NGC 253. In diffuse regions affected by X-rays and/or cosmic-rays, the chemistry leading to H₂O⁺ is initiated by the charge transfer between ionised hydrogen and oxygen atoms to form O⁺, which is subsequently followed by reactions with H₂ to first form OH⁺ and then H₂O⁺. Alternatively, in regions with higher molecular fractions, both OH⁺ as well as H₂O⁺ can be formed via reactions between ${\text{H}}_3^{+}$ and oxygen atoms (Hollenbach et al. 2012, and references therein). Lastly, H₃O⁺ can also be formed from hydrogen abstraction reactions of H₂O⁺ in denser environments. If the exchange reaction is faster than other competing reactions then the observed upper limits of the OPR, point to gas kinetic temperatures that are greater than ∼20 K.

Since the diffuse and dense gas volumes along the sight lines towards the galaxies under study here are well mixed, the OPR of parent species such as ${\text{H}}_3^{+}$, H₂, and H₂O, from which H₂O⁺ originates, may impact the observed OPR of H₂O⁺ see for example the studies on the OPR of H₂ presented by Flower et al. 2006. Destruction reactions of both the ortho- and para-forms of ${\text{H}}_3^{+}$ via dissociative recombination have been experimentally investigated by Glosík et al. (2010), who observe a preferential recombination of this molecule’s para-state compared to that of its ortho-state, which can also lead to OPRs >3. Recently, Novotný et al. (2019) found a similar dependence on the rotational state, for the dissociative recombination of HeH⁺ ions. However, a better understanding of the OPR of H₂O⁺ requires a detailed understanding (not only) of what fraction of the observed H₂O⁺ is formed along each of the different chemical pathways towards its formation, alongside the efficiency with which either H₂O⁺ state (ortho or para) is formed in both quiescent and turbulent and/or shocked regions, but also the destruction of this ionic species which occurs primarily via dissociative recombination and hydrogen abstraction reactions with H₂. In Appendix C we briefly evaluate deviations in the observed OPR of H₂O⁺ from the original OPR at formation, by carrying out a steady-state analysis of H₂O⁺ chemistry.

Furthermore, investigating the influence of cosmic-ray ionisation on the OPR of H₂O⁺ for star-forming regions within the Milky Way, Jacob et al. (2020) find that the OPR of H₂O⁺ thermalises to the equilibrium value of three with increasing rates of cosmic-ray ionisation up to a value of ∼​2 × 10⁻¹⁶ s⁻¹, beyond which there exists no credible correlation. The correlation between cosmic-ray ionisation rates and OPRs suggests that the higher abundances of atomic hydrogen in regions exposed to higher cosmic-ray fluxes are able to efficiently drive the proton-exchange reaction causing a change in the OPR. Putting the current results into context, we find the data points for both studied galaxies to lie in the thermalised region of Fig. 8, with OPRs near three within the uncertainties. Given the uncertainties in the column densities derived from both the HI and the OH⁺ line profiles, which are used to compute the cosmic-ray ionisation rates, we represent the cosmic-ray ionisation rates as lower limits. In analogy to the Milky Way, we would expect the cosmic-ray ionisation rates to be higher towards the centres of these galaxies, however in contrast, we derive values closer to the canonical values derived towards the disk of the Milky Way (Indriolo et al. 2015; Jacob et al. 2020), outside of the central molecular zone.

Fig. 8. Observed OPR of H2O+ versus the cosmic-ray ionisation rates derived using the steady-state analysis of the OH+ and H2O+. The results from this work are marked by teal pentagons while the blue diamonds and black triangles display the results for sight lines studied in the Milky Way by Indriolo et al. (2015) and Jacob et al. (2020). For comparison we display the OPR of H2O+ also derived by Indriolo et al. (2015) towards Galactic centre (GC) sight lines using purple stars. The dashed red line marks the equilibrium OPR value of three.

4. Discussion

The gas properties of the ArH⁺ bearing cloud volumes that we have derived in Sect. 3.4 towards the nearby galaxies studied in this work are comparable to those found by Schilke et al. (2014) and Jacob et al. (2020) for clouds along the LOS towards star forming regions in the disk of the Milky Way. Furthermore, the cosmic-ray ionisation rates derived towards both NGC 253 and NGC 4945 in Sect. 3.3 are comparable to the average cosmic-ray ionisation rate derived towards the disk of the Milky Way, of (2.3 ± 0.3)×10⁻¹⁶ s⁻¹ (Indriolo et al. 2015; Jacob et al. 2020). The molecular gas traced by HF has an almost two orders of magnitude higher column density than HI (as illustrated in Fig. 2), suggesting that most of the gas, by volume, along the LOS towards these galaxies is dense. Ostensibly, this brings into question whether ArH⁺ is a probe of diffuse atomic gas with a very small molecular fraction (see Sect. 3.4). However, as previously discussed by van der Tak et al. (2016), who observed OH⁺ and H₂O⁺, other well-known tracers of (predominantly) atomic gas, towards both NGC 253 and NGC 4945, the presence of gas with low f_H2 may suggest an inhomogeneous ISM with local pockets of diffuse clouds mixed with denser material. Perhaps, the marginally lower cosmic-ray ionisation rates that we find for the galaxies studied here, in comparison to the Milky Way, may be a result of the larger volumes of dense gas present along the LOS. We note, however, that values for the cosmic-ray ionisation rate estimated for the nuclear regions of ULIRGs like NGC 4418 and Mrk 231, which harbour much more extreme star bursts than NGC 253 and NGC 4945, are almost three orders of magnitude higher with a lower-limit of ζ_p(H)​ > ​10⁻¹³ s⁻¹ (González-Alfonso et al. 2013, 2018), similar to values derived towards the Milky Way’s Galactic centre region, indicative of higher levels of energetic phenomena present at the centres of these galaxies in comparison to their disks. This is also the case towards the diffuse atomic gas rich lines-of-sight studied towards PKS 1830−211, who derived column densities for atomic gas tracers such as ArH⁺ (Müller et al. 2015), p-H₂O⁺, and OH⁺ (Muller et al. 2016) that are roughly an order of magnitude greater than those that we derived towards NGC 253 and NGC 4945 while the column density values for molecular gas tracers such as HF are comparable between both studies (see Table 2 of Muller et al. 2016). This results in cosmic-ray ionisation rates of ∼2 × 10⁻¹⁴ s⁻¹ and s ∼ 3 × 10⁻¹⁵ s⁻¹ for the SW and NE lines-of-sight, respectively.

In addition, while we quote LOS averaged values for the cosmic-ray ionisation rates towards both NGC 253 and NGC 4945 using the steady-state analysis of OH⁺ and H₂O⁺, their interpretation is not trivial as the cosmic-ray ionisation rate is likely to vary spatially, across these galaxies. This variation is best demonstrated by observations towards different molecular cloud sight lines within the Milky Way by Indriolo et al. (2015) and Jacob et al. (2020) and reflects the proximity of the studied molecular clouds to the nearest supernova remnants, which are likely sites for particle acceleration amongst other propagation effects. Therefore we expect the cosmic-ray ionisation rates determined from compact regions of star-formation like the central molecular zones of the galaxies studied here, to be orders of magnitude greater. In general with star-formation rates higher than that of the Milky Way by almost two orders of magnitude, one would anticipate higher cosmic-ray ionisation rates towards both NGC 253 and NGC 4945 since the injection of cosmic-rays is governed by star-formation and accelerated by associated shocks like supernova remnants and stellar wind bubbles. Observationally, higher rates of cosmic-ray ionisation have been estimated towards the centre of NGC 253 through the detection of gamma-rays at TeV energies using ground based Cherenkov telescopes (Acero et al. 2009). Exhibiting both starburst and AGN activities the nuclear environment of NGC 4945 is more complex, as a significant portion of the injected cosmic-ray energy is used in the production of pions and γ-rays. Wojaczyński & Niedźwiecki (2017) have shown that Seyfert galaxies like NGC 4945 have larger γ-ray luminosity’s than the calorimetric limit which sets the total energy of cosmic-rays used for the production of pions assuming that each supernova explosion injects 10⁵⁰ ergs of cosmic-ray energy while that of starbursts like NGC 253 is 50 % lower.

Moreover, while ionisation by cosmic-rays plays an important role in heating the gas in the nuclear environment of galaxies (Bradford et al. 2003), heating by X-rays (Usero et al. 2004) and by stellar UV radiation (in widespread photodissociation regions (PDRs) Hollenbach & Tielens 1997), dynamical shock heating (García-Burillo et al. 2001) to a name a few, can all be prominent sources of heating in these regions and it is not possible to discriminate between their contributions merely by an analysis such as ours.

5. Conclusions

Well established as a tracer for purely atomic gas, outside of the Milky Way the noble gas species ArH⁺ had only been detected towards the intermediate redshift z = 0.89 gravitational lens system PKS 1830−211. In this work we extend the notion that ArH⁺ is an ubiquitous tracer of diffuse atomic gas towards external galaxies by presenting observations of its J​ = ​1−0 transition.

We present the detection of ArH⁺ in absorption towards two nearby galaxies, NGC 253 and NGC 4945 and the non-detection towards Arp 220. In addition to ArH⁺ we also report the detection of the J = 3/2−1/2 transition of p-H₂O⁺ at 607 GHz towards the former. We compare the observed profiles of ArH⁺ and p-H₂O⁺ lines with those of the HI 21 cm line and lines from the well known atomic gas tracers OH⁺ and o-H₂O⁺ and the molecular gas tracer, HF. Using the cosmic-ray ionisation rates derived by analysing the steady-state chemistry of OH⁺ and H₂O⁺ and by assuming that the cloud populations bearing ArH⁺ are exposed to the same cosmic-ray ionisation rates, we derive the molecular fraction of the gas traced by ArH⁺ towards both galaxies to be between ​10⁻³ and ∼​10⁻². This is consistent with estimates made within the Milky Way by Schilke et al. (2014), Neufeld & Wolfire (2016), Jacob et al. (2020). From the detection of p-H₂O⁺ we estimate the H₂O⁺ OPR of 7.0 and 3.5 towards NGC 253 and NGC 4945, and subsequently derive values for the nuclear-spin temperature that are greater than 15 and 24 K, respectively. The variation in the H₂O⁺ OPR with cosmic-ray ionisation rates follows the trend displayed by sight line components towards star-forming regions within the Milky Way other than those towards the Galactic centre. This is likely because our observations are not sensitive to spatial variations along the sight lines studied and only probe the average properties.

In order to comment on the ubiquity and chemistry of ArH⁺ and its possible role in probing the energetics of extragalactic sources, the searches for ArH⁺ need to be extended to a wider source sample. Since its chemistry and properties are widely studied, alongside that of OH⁺ and o-H₂O⁺, it may be promising to search for ArH⁺ towards a range of extragalactic sources hosting various levels of nuclear activity, which was a selection criterion for the Herschel EXtraGALactic (HEXGAL) survey which was aimed at studying the physical and chemical composition of the ISM in galactic nuclei using HIFI spectroscopy.

Acknowledgments

This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX) under the project id M9519C_103. APEX is a collaboration between the Max-Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. We would like to express our gratitude to the APEX staff and science team for their continued assistance in carrying out the observations presented in this work. We would like to thank Dr. Paul Dagdigian for providing us with the ArH⁺ collisional rate coefficients and the anonymous referee for their careful review of the article and valuable input. We are thankful to the developers of the C++ and Python libraries and for making them available as open-source software. In particular, this research has made use of the NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001) and matplotlib (Hunter 2007) packages.

References

Appendix A: HI analysis

In this appendix we discuss the HI 21 cm line analysis. The HI absorption (T_on) and emission spectra (T_off), along with derived quantities such as the optical depth, spin temperature and column density are presented in Figs. A.1 and A.2 for NGC 253 and NGC 4945, respectively. The solid black curve in the top panel of these figures, represents the continuum normalised on-source absorption profiles ($T^{\text{on,obs}}/T_{\text{cont}}^{\text{sou,obs}}$) smoothed to the resolution of the emission line data. The off-source spectrum is taken from a position next to the source of interest (well outside of the beam width). Similar to the Milky Way disk, the HI data towards these galaxies are subject to strong spatial fluctuations and using a single off position would be error prone. Therefore, following Winkel et al. (2017), a spatial filtering technique using a ring (or doughnut-shaped) kernel was applied to obtain an interpolated brightness temperature, as shown in the second panel (black curve). For comparison, in the panel displaying spin temperatures (computed for regions with τ > 0.01), we also plot the brightness temperature in blue and lastly, the column density panel also displays the uncorrected column density profile, $N_{\text{HI}}^{\ast }$. The blue and grey shaded regions indicate the 1σ (68% percentile) and 3σ (99.7% percentile) confidence intervals, respectively. A detailed description is presented in Winkel et al. (2017). As a caveat in this analysis, we mention here that this filtering technique may not be what is best suited to describe the often observed, ‘butterfly’-shaped distribution of gas in these external galaxies and a proper treatment of the same, would require the use of a tilted ring model. In the on-source HI absorption profile towards NGC 253, we note the appearance of an emission feature between 5 and 40 km s⁻¹, which is merely an artefact arising from the post-processing (cleaning) of the data or is caused by HI emission from the disk of the galaxy.

Fig. A.1. From top to bottom: HI absorption and emission spectra, optical depth, spin temperature, and HI column density towards NGC 253.

Fig. A.2. Same as Fig. A.1, but towards NGC 4945.

Appendix B: Continuum level uncertainties

We now briefly assess the reliability of the absolute calibration of the continuum brightness temperatures used in our analysis that were derived from the sub-mm spectra of the ArH⁺ and p-H₂O⁺ lines taken with the APEX telescope. In Fig. B.1, for both NGC 253 and NGC 4945, we plot the continuum brightness temperatures measured for all individual scans versus the time over velocity intervals used to determine the continuum level and baseline. The scatter across the different scans displayed is found to lie between 11 and 17% for the two sources with an average value of 14%, thereby yielding a fairly stable continuum level.

Fig. B.1. Continuum fluctuations across scans for NGC 253 (purple triangles) and NGC 4945 (yellow pentagons) with the dashed lines indicating the median continuum levels for each source.

Appendix C: Steady-state chemistry of H₂O⁺

As discussed in Gerin et al. (2013), the relationship between the initial OPR of the H₂O⁺ ions when formed, OPR₀, with that which is observed, depends upon the relative rates of these processes. We adapt the steady-state analysis of H₂Cl⁺ presented in Neufeld et al. (2015) but for H₂O⁺, to estimate the relative importance of the different reactions in thermalising the OPR. Hence, the OPR of H₂O⁺ is given by,

$$\begin{array}{c}\hfill \text{OPR}=\frac{k_{\mathrm{po}}n(\mathrm{H})+k_{\mathrm{d}r}{n}_e{\text{OPR}}_0+k_{\mathrm{ha}}n({\mathrm{H}}_2)OP{R}_0}{k_{\mathrm{op}}n(\mathrm{H})+k_{\mathrm{dr}}{n}_e+k_{\mathrm{ha}}n({\mathrm{H}}_2)},\end{array}$$(C.1)

where k_po and k_op are the rate coefficients for the forward and backward reactions in (5) and k_dr and k_ha are the reaction rates for the dissociative recombination and hydrogen abstraction reactions, respectively. Assuming a detailed balance, such that k_po = k_opOPR_LTE(H₂O⁺), the above equation can be re-written as

$$\begin{array}{c}\hfill \text{OPR}=x{\text{OPR}}_{\mathrm{LTE}}+(1-x){\text{OPR}}_0,\end{array}$$(C.2)

where x = k_opn(H)/(k_opn(H)+k_drn_e+k_han(H₂)). The rate coefficient for the proton exchange reaction (5), k_op, has not been measured but is assumed to be of the order of 10⁻¹⁰ cm³ s⁻¹ in view of the rate coefficients determined for exchange reactions of other related species such as H₂, p-H₂ + H⁺ ⇌ o-H₂ + H⁺, for which k = 4.15 × 10⁻¹⁰ cm³ s⁻¹ (Honvault et al. 2011). While the rate coefficient for the dissociative recombination reaction is k_dr ∼ 7.44 × 10⁻⁷ cm³ s⁻¹, that for the abstraction reaction is k_ha = 6.1 × 10⁻¹⁰ cm³ s⁻¹ (computed at a gas temperature of 100 K, typical for diffuse gas environments), competes with k_op. Using values for n_H, n_H2, and n_e as discussed in Sects. 3.3 and 3.4 and k_op = 4.15 × 10⁻¹⁰ cm s⁻³, we estimate x to be 0.74. This implies that, for values of k_op > 4.15 × 10⁻¹⁰ cm s⁻³, x → 1 and that the OPR is close to the values at LTE. On the contrary, if k_op < 4.15 × 10⁻¹⁰ cm s⁻³ then x → 0 implying that the OPR is close to OPR₀. Given the many uncertainties involved in this calculation we can only conclude that the OPR derived using Eq. 6 merely represents a limit.

All Tables

All Figures

Fig. 1. From top to bottom: normalised absorption spectra of ArH+, p-H2O+, o-H2O+, OH+ (at 971 and 1033 GHz), HI and HF as well as the CO (1–0) emission line spectrum for comparison, towards NGC 253 (left) and NGC 4945 (right), respectively. In the spectra for NGC 253, the dotted light blue curves display the Gaussian fit to the emission component. The dashed black line represents the individual absorption profiles after subtracting the Gaussian fit with their Wiener filter fits overlaid in red for all species except HI and CO. The relative intensities of the hyperfine structure (HFS) components of the o-H2O+ and OH+ transitions are shown in black above their respective spectra and the grey shaded regions display their HFS deconvolved spectra. The vertical dashed black lines mark the systemic velocity of NGC 253 and NGC 4945 at 240 km s−1 and 563 km s−1, respectively. The 971 GHz OH+ line profile towards NGC 253 is saturated at blueshifted velocities while the 1033 GHz OH+ spectrum is potentially contaminated by the 13CH2CHCN (2110,11–209,12) line at 1032.783 GHz (fit shown by the dark blue curve). Also marked here in blue is contamination from the C2H5CN (379,29–376,32) line at 1033.868 GHz. The 971 GHz OH+ spectrum towards NGC 4945 is contaminated by the CH3OH Jk = 94,6 − 8,53 E line near 959.900 GHz originating from the image sideband (marked in blue).

In the text

	Fig. 2. From top to bottom: column density distributions (black) of ArH+, p-H2O+, o-H2O+, OH+, HI, HF alongside the corresponding scaled H2 profile ([HF]/[H2] = 5×10−10; Emprechtinger et al. 2012) and CO towards NGC 253 (left) and NGC 4945 (right), respectively. The corresponding uncertainties are displayed by the blue shaded region.
In the text

Fig. 3. Top: distribution of molecular gas fraction (fH2) traced by ArH+ (blue), OH+ (red) and CO (orange). Bottom: contours of X(ArH+) (blue) and X(OH+) (red) abundances with respect to NH in the fH2–ζp(H) plane. Blue triangles and red diamonds represent the corresponding values derived from the LOS observations presented in this work computed per channel in the velocity interval (Δυ ∼ ​4.5 km s−1) between 55 and 295 km s−1, and 445 and 725 km s−1 for NGC 253 and NGC 4945, respectively.

In the text

Fig. 4. RADEX modelled excitation temperature as a function of electron density, for fixed values of the total gas density at 35 (red), 100 (orange) and 1000 cm−3 (dark blue) for the ArH+ (1–0) transition. Furthermore, each model is run for fixed gas temperatures at 30 (solid), 300 (dashed) and 3000 K (dotted), respectively. The horizontal dashed grey line marks an excitation temperature of 2.73 K. The inset zooms in on 10−4 < xe < 10−2 with the vertical dashed purple line indicating xe = 1.5 × 10−4 which corresponds to the value of xe used in our calculations.

In the text

Fig. 5. RADEX modelled excitation temperatures (left) and optical depths (right) as a function of electron density, for fixed values of the total gas density at 35 (red), 100 (orange) and 1000 cm−3 (dark blue) for the ArH+ (2–1) transition. Furthermore, each model is run for fixed gas temperatures at 30 (solid), 300 (dashed) and 3000 K (dotted), respectively. The horizontal dashed grey line marks an excitation temperature of 2.73 K (in the left-hand panel) while the vertical dashed purple line indicates, xe = 1.5 × 10−4.

In the text

Fig. 6. Top panel: column density per velocity interval distribution of o-H2O+ (solid purple) and p-H2O+ (dashed blue) over the entire velocity interval covering absorption towards NGC 253. Middle panel: OPR H2O+ distribution (left) and corresponding histogram (right). Bottom panel: nuclear spin temperature distribution (left) and corresponding histogram (right). The median and 1 σ levels of the OPR and nuclear spin temperature are marked by solid and dotted black lines, respectively.

In the text

	Fig. 7. Same as Fig. 6, but towards NGC 4945. The median value for Tns derived from the distribution is the same as that obtained from deriving Tns using the median OPR for H2O+.
In the text

Fig. 8. Observed OPR of H2O+ versus the cosmic-ray ionisation rates derived using the steady-state analysis of the OH+ and H2O+. The results from this work are marked by teal pentagons while the blue diamonds and black triangles display the results for sight lines studied in the Milky Way by Indriolo et al. (2015) and Jacob et al. (2020). For comparison we display the OPR of H2O+ also derived by Indriolo et al. (2015) towards Galactic centre (GC) sight lines using purple stars. The dashed red line marks the equilibrium OPR value of three.

In the text

	Fig. A.1. From top to bottom: HI absorption and emission spectra, optical depth, spin temperature, and HI column density towards NGC 253.
In the text

	Fig. A.2. Same as Fig. A.1, but towards NGC 4945.
In the text

	Fig. B.1. Continuum fluctuations across scans for NGC 253 (purple triangles) and NGC 4945 (yellow pentagons) with the dashed lines indicating the median continuum levels for each source.
In the text