ALCHEMI Finds a "Shocking" Carbon Footprint in the Starburst Galaxy NGC 253

The abundance of interstellar molecules depends on the balance between their formation and destruction processes in the gas phase and on grain surfaces. Exchange processes can transform molecules in one phase to another; gas-phase molecules can freeze onto dust grains as ice (adsorption), while molecules on grain surfaces can sublimate into the gas phase (desorption). Knowing both the gas- and ice-phase abundances is, therefore, necessary for a comprehensive understanding of their chemical composition and their related physical conditions. Carbon dioxide (CO₂) is one of the most dominant forms of ice on interstellar dust (Öberg et al. 2011) together with H₂O and CO. CO₂ is, in fact, one of the ice species detectable in extragalactic sources (e.g., Lahuis et al. 2007; Yamagishi et al. 2015). It is obvious from strong CO rotational emission that the gas-phase CO is abundant, but CO₂ is thought to reside more preferentially on dust due to inefficient gas-phase formation (Garrod & Pauly 2011). Because CO₂ has a higher desorption energy than that of CO, the presence of abundant gas-phase CO₂ requires stronger desorption mechanisms. While CO₂ can be detected via rotational–vibrational transitions in warm gas (≳several 100 K) (e.g., Boonman et al. 2003) or in ice with broader line features (e.g., Yamagishi et al. 2015), it cannot be observed in cold gas due to the lack of a permanent dipole moment.

Although the gas-phase abundances of such species without a dipole moment cannot be directly measured through commonly observed rotational transitions, it has been proposed that they can be estimated from their protonated forms (Herbst et al. 1977; Agúndez et al. 2015; Rivilla et al. 2019). The protonated form of CO₂, HOCO⁺, was first detected in Sgr B2 by Thaddeus et al. (1981), but its line identification required spectroscopic confirmation by Defrees et al. (1982). Since then, HOCO⁺ has been detected in translucent clouds (Turner et al. 1999), low-/high-mass star-forming regions (Sakai et al. 2008; Vastel et al. 2016; Fontani et al. 2018; Majumdar et al. 2018), the Galactic Center (Sgr A and B2 clouds; Minh et al. 1988, 1991), starburst galaxies (Martín et al. 2006; Aladro et al. 2015) and a z ∼ 0.9 molecular absorber (Muller et al. 2013). Among them, Minh et al. (1991) and Armijos-Abendaño et al. (2015) found about two orders of magnitude higher fractional abundances of HOCO⁺ in the Galactic Center than in spiral-arm molecular clouds.

The main formation paths of HOCO⁺ are gas-phase reactions, including the protonation of CO₂ such as (Vastel et al. 2016; Bizzocchi et al. 2017)

$$\begin{eqnarray*}&&{\mathrm{CO}}_{2}+{{\rm{H}}}_{3}^{+}\longrightarrow {\mathrm{HOCO}}^{+}+{{\rm{H}}}_{2},\end{eqnarray*}$$

or an ion-neutral reaction

$$\begin{eqnarray*}&&{{\rm{HCO}}}^{+}+{\rm{OH}}\longrightarrow {{\rm{HOCO}}}^{+}+{\rm{H}}.\end{eqnarray*}$$

In the former route, the HOCO⁺ abundance can increase due to the evaporation of CO₂ from grain surfaces because CO₂ is one of the most abundant forms of carbon on grain surfaces. Ice can sublimate thermally (e.g., in the vicinity of protostars), or nonthermally (e.g., photodesorption, cosmic-ray-induced evaporation, or shock sputtering). While the latter route from HCO⁺ is considered dominant in at least some parts of high- or low-mass star-forming regions (Fontani et al. 2018; Majumdar et al. 2018), the protonation of CO₂ can be the dominant route when there is a fast mechanism of CO₂ sublimation.

Shocks, one such driver of the ice sublimation process, are ubiquitous in galactic centers. In many barred-spiral galaxies, galactic centers host intersections of x₁ and x₂ orbits ¹⁹ (Athanassoula 1992; Seo & Kim 2013; Sormani & Li 2020). At these orbital intersections, shocks are naturally expected. The abundances of typical shock molecular tracers such as CH₃OH and SiO have been found to be enhanced at locations of orbital intersections in IC342 and M83 (Meier & Turner 2005; Harada et al. 2019).

In this paper, we report an enhancement of HOCO⁺ at the orbital intersections near the center of the starburst galaxy NGC 253. NGC 253 is one of the nearest (d = 3.5 Mpc; Rekola et al. 2005) and most studied starburst galaxies. Its nuclear ring forms a central molecular zone (CMZ) within a few hundred parsec scale, with a mass of ∼2 × 10⁸ M_⊙ within a radius of r = 150 pc (Leroy et al. 2015). This large reservoir of molecular gas enables active star formation (Sakamoto et al. 2006, 2011; Bolatto et al. 2013; Krieger et al. 2019; Rico-Villas et al. 2020), which affects the properties of the molecular gas (e.g., through heating; Mangum et al. 2019). Gas in the bar orbit (x₁ orbits) is being fed to the center of NGC 253 (x₂ orbits), and shocks are expected when this gas flow collides with the nuclear ring, as discussed above for other galaxies. Figure 1 (left) shows IRAC 8 μm (Dale et al. 2009; LVL team 2009) and CO(2-1) integrated-intensity (Atacama Large Millimeter/submillimeter Array, ALMA, data 2018.1.01321.S; PI: Faesi) images covering most of NGC 253 to indicate the locations of these orbits. The association between shocks at orbital intersections and chemistry has been proposed by García-Burillo et al. (2000) and Meier et al. (2015).

Zoom In Zoom Out Reset image size

At the center of NGC 253, a large number of molecular species are detectable (Martín et al. 2006; Meier et al. 2015). To fully explore the chemical complexity in this galaxy, we conducted the ALMA large program ALMA Comprehensive High-resolution Extragalactic Molecular Inventory (ALCHEMI; Martín et al. 2021). ALCHEMI is a wide-frequency, unbiased spectral scan mosaic toward the CMZ of NGC 253 at a common 16 resolution. This survey has discovered high cosmic-ray ionization rates (Holdship et al. 2021; Harada et al. 2021; Holdship et al. 2022), detected a phosphorus-bearing species for the first time in an extragalactic source (Haasler et al. 2022), and detected methanol masers (Humire et al. 2022) in the CMZ of NGC 253. We utilize the ALCHEMI data to study the multiple transitions of HOCO⁺.

This paper is organized as follows. Section 2 describes our observational parameters and data analysis, and images of the integrated intensities are presented in Section 3. The derived column densities of HOCO⁺ are shown in Section 4, while the HOCO⁺/CO₂ ratios are discussed using chemical models in Section 5. In Section 6, we discuss our results including a comparison with ice abundances. Our results are summarized in Section 7.

The ALCHEMI spectral survey mosaic of the NGC 253 CMZ was performed between 2017 and 2019, with ALMA's 12 m antenna and 7 m antenna arrays. It covered a broad frequency range between 84 and 373 GHz (Bands 3 to 7, avoiding deep atmospheric lines), down to sensitivities of ∼10 mK per 10 km s⁻¹ channel. The ALCHEMI data products have a uniform angular resolution of 16 and cover a field of view of 50'' × 20'' centered on the CMZ of NGC 253 (phase center: α = 00^h47^m3328, $\delta =-25^\circ 17^{\prime} 17\buildrel{\prime\prime}\over{.} 7$, in the International Celestial Reference System, ICRS). The extent of the largest recoverable angular scale is greater than or equal to 15''. A detailed description of the ALCHEMI survey products can be found in Martín et al. (2021).

We extracted the line cubes around the transitions of HOCO⁺ with a velocity resolution binned to 10 km s⁻¹. Within the ∼290 GHz coverage of the ALCHEMI survey, we find 14 detectable HOCO⁺ transitions, occurring every ∼21.4 GHz. However, some transitions are severely blended with transitions from other species and are not used in this analysis. The spectroscopic parameters and spectral channel rms values for the HOCO⁺ transitions used in this paper are listed in Table 1.

Table 1. HOCO⁺ Spectroscopic Properties and rms Noise Values

Transition	${\nu }_{\mathrm{rest}}^{(a)}$	Eup (b)	$\mathrm{log}({A}_{\mathrm{ul}})$ (c)	rms (d)	rms (e)	Blending (f)
	(GHz)	(K)	(s−1)	(mJy beam−1)	(mK)
40,4-30,3	85.531	10.3	−4.63	0.19	12.	Potential blending with CH3CCH. Minor blending with U-lines
50,5-40,4	106.914	15.4	−4.33	0.20	8.4	N/A
60,6-50,5	128.295	21.6	−4.08	0.33	9.5	Minor blending with U-lines
70,7-60,6	149.676	28.7	−3.88	0.42	9.0	Potential blending with U-lines
80,8-70,7	171.056	36.9	−3.70	0.75	12.	Potential blending with CH3CCH.
120,12-110,11	256.566	80.0	−3.16	0.99	7.2	Potential blending with CH3CCH, HC3N v7 = 2
Transitions below were not used for analysis
90,9-80,8	192.435	46.2	−3.54			Blended with U-line
100,10-90,9	213.813	56.4	−3.40			Blended with C2H5OH
110,11-100,10	235.190	67.7	−3.28			Blended with SO2
130,13-120,12	277.941	93.4	−3.06			Blended with U-line
140,14-130,13	299.314	107.7	−2.96			Non detection
150,15-140,14	320.686	123.1	−2.87			Non detection
160,16-150,15	342.056	139.6	−2.78			Non detection
170,17-160,16	363.424	157.0	−2.70			Non detection

Note. (a) Rest frequency. (b) Upper level energy of the transition. (c) A_ul : Einstein coefficient of spontaneous emission. All values were taken from the Cologne Database for Molecular Spectroscopy (CDMS; https://cdms.astro.uni-koeln.de; Müller et al. 2001, 2005; Bizzocchi et al. 2017). (d) and (e) rms values of a single channel with Δv = 10 km s⁻¹ in mJy beam⁻¹ and mK units. (f) Presence of blending; transitions are shown with quantum numbers ${J}_{{K}_{a},{K}_{c}}$. Only K_a = 0 transitions are shown because transitions with K_a ≠ 0 are not detected due to their higher energy state and lower Einstein coefficients. The upper part of this table shows transitions used for the analyses of this paper, whose line shapes are separable from neighboring lines. The lower part lists transitions with severe blending or without reliable detections. "Potential blending" means the case where the line centers are separated by more than 200 km s⁻¹, but the line wings can contaminate the moment maps of the HOCO⁺ transitions.

Download table as:  ASCIITypeset image

Figures 2(a)–(d) show the velocity-integrated-intensity (moment 0) images of HOCO⁺ in multiple transitions (see Table 1 for their properties). Many transitions used in this work have neighboring lines, and it is not possible to make the moment 0 images simply by collapsing neighboring channels. To exclude this contamination, we applied a 3D mask made from position–position–velocity space by only including pixels with CO J = 1 − 0 detections above the 20σ level. The 20σ cutoff may sound unnecessarily high, but the signal-to-noise ratios of HOCO⁺ transitions are more than 100 times lower than that of CO(1-0). Therefore, this mask does not exclude any notable HOCO⁺ emission, but helps to exclude contamination from nitrogen sulfide transitions neighboring with CO(1-0) to be included in the integrated-intensity map. Despite the elimination of contamination with this mask, the only transition that is free from contamination is 5_0,5−4_0,4 (Figure 2(a)). Other images that are relatively less affected by blending are also shown in Figure 2 (panels b–d). The level of contamination is usually very low (<10%) except for giant molecular cloud (GMC) 5, where there is little HOCO⁺ emission and stronger emission from neighboring lines. Transitions 4_0,4−3_0,3 and 8_0,8−7_0,7 can still be used to obtain column densities as we use spectral fitting, but their images are not shown. A low-excitation line of HOCO⁺ (5_0,5−4_0,4; E_up = 15.4 K) shows peaks near the outer CMZ, in GMCs 1, 7, and 9 (GMC numbering is shown in panel (f); Leroy et al. 2015; see also Appendix C for coordinates). On the other hand, the higher-excitation transitions (E_up ≳ 30 K) peak closer to the center (GMCs 3 and 6).

Zoom In Zoom Out Reset image size

Panels (e) and (f) of Figure 2 show integrated-intensity images of CH₃OH(2_k −1_k ) (group of transitions at ν_rest ∼ 96.74 GHz with the strongest transition ${2}_{0}^{+}-{1}_{0}^{+}$) ²⁰ and H¹³CO⁺(1-0) (ν_rest = 86.75 GHz) for comparison. The emission distribution of CH₃OH(2_k −1_k ) is similar to that of the low-J transitions of HOCO⁺, as well as the low-J transitions of HNCO, and SiO (Meier et al. 2015; K.-Y. Huang et al., 2022, in preparation). On the other hand, the distribution of H¹³CO⁺(1-0), which is rather similar to that of molecules with strong emission (e.g., CO, HCO⁺, HCN, CS, etc.; Meier et al. 2015; Martín et al. 2021), is clearly different from that of HOCO⁺. The H¹³CO⁺(1-0) emission is concentrated near the center of NGC 253 (GMCs 3–7) instead of the outer CMZ.

The similarity between the integrated intensities of HOCO⁺ (5_0,5−4_0,4) and CH₃OH (2_k −1_k ) and the difference between those of HOCO⁺ (5_0,5−4_0,4) and H¹³CO⁺(1-0) are highlighted in Figure 3. While the HOCO⁺(5_0,5−4_0,4)/CH₃OH (2_k −1_k ) ratios are relatively constant, the HOCO⁺(5_0,5−4_0,4)/H¹³CO⁺(1-0) ratio varies significantly. All these transitions have relatively low upper state energies (HOCO⁺ (5_0,5−4_0,4): E_u = 15.4 K; CH₃OH (${2}_{0}^{+}-{1}_{0}^{+}$): E_u = 7.0 K; H¹³CO⁺(1-0): E_u = 4.2 K), and these ratio maps should be good proxies for the variations in column densities of these molecules.

Zoom In Zoom Out Reset image size

Figure 4(a) shows a column-density map of HOCO⁺. These column densities were derived using the public software CASSIS ²¹ (Vastel et al. 2015) supplied with spectroscopic constants from the spectroscopic database from the CDMS (Müller et al. 2001, 2005). CASSIS calculates molecular column densities based on input spectral line brightness temperatures with consideration of optical depths, either with an LTE assumption or, if collisional rates are available, non-LTE assumption. We used a Markov Chain Monte Carlo (MCMC) algorithm assuming LTE to fit the column densities, excitation temperatures, line velocities, and line widths (see Appendix A). The column densities were calculated only for pixels with a >3σ detection of both HOCO⁺(4_0,4−3_0,3) and HOCO⁺(5_0,5−4_0,4) at velocities within 30 km s⁻¹ from the line center. Line-center velocities used for this 3σ detection criterion are determined from the image cube of CO(1-0) from the ALCHEMI data. These line-center velocities from CO may be different from the HOCO⁺ velocities fitted from CASSIS. Instead of deriving the column densities and excitation temperatures on a pixel-by-pixel basis, we bin the intensities within hexagonal pixels with a horizontal length of 08, half of the image spatial resolution of 16, to reduce the computational time running CASSIS. Examples of the spectral fitting are shown for hexagonal pixels located at (x_offset, y_offset) = (−173, −104) and (20, 21) in Figure 5, where the offset is taken from the phase center. In general, the observed spectra fit well with the LTE spectra. The resulting column-density distribution appears similar to that of the moment 0 images of low-excitation transitions (e.g., Figure 2(a)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The excitation temperatures derived from the above spectral fitting are shown in Figure 4(b). While regions far from the kinematic center ²² of NGC 253 show low excitation temperatures of ≲10 K (GMCs 1, 2, 8, and 9), the excitation temperatures become higher, up to 40 K, near the center of NGC 253 (GMCs 3 and 6).

We also obtain the total hydrogen column densities (Figure 4(c)) to calculate the fractional abundances of HOCO⁺ (i.e., the HOCO⁺ column density divided by the total hydrogen column density N(H₂)). The total H₂ column densities were derived from the dust continuum image at 361.5 GHz shown by Harada et al. (2021) with the derivation method based on Hildebrand (1983) for pixels above 3σ detection. A simplified formula is given as Equation (3) in Mangum et al. (2019)

$$\begin{eqnarray}&&N({H}_{2})({\mathrm{cm}}^{-2})\sim 7.0\times {10}^{22}{R}_{{dg}}{\left(\displaystyle \frac{\lambda (\mathrm{mm})}{0.4}\right)}^{\beta }\displaystyle \frac{{T}_{R}({\rm{K}})}{{T}_{d}({\rm{K}})},\end{eqnarray} \tag{ 1 }$$

where N(H₂) is the molecular hydrogen column density, R_dg is the dust-to-gas mass ratio, λ is the wavelength, T_R is the radiation temperature, and T_d is the dust temperature. This formula is valid for h ν ≪ kT_d (the Rayleigh–Jeans approximation). We use an emissivity of β = 1.5, a dust temperature of T_dust = 30 K, and a dust-to-gas mass ratio of 150 following Mangum et al. (2019). This estimate of the dust temperature is close to the observed value, but some dust components may be warmer. Pérez-Beaupuits et al. (2018) derived dust temperature components of 37, 70, and 188 K in the central region of NGC 253 using their assumed source size of 173 × 92 from their Herschel and SOFIA observations. These components contain mass fractions of 65, 26, and 9%, respectively. If the dust is warmer, the actual column density should be smaller by a similar factor; for instance, a factor of 5 smaller if T_dust = 150 K. The column-density dependence on the dust temperature becomes larger than $\propto \tfrac{1}{{T}_{d}}$ when the dust is cold (T_d ≲ 20 K) and one cannot use the Rayleigh–Jeans approximation, but we expect that there is a very small amount of cold dust in the center of NGC 253.

The fractional abundance of HOCO⁺ is higher at larger distances from the center of NGC 253, and it decreases by more than an order of magnitude at the center (Figure 4(d)). At the peaks of HOCO⁺(GMCs 1, 8, and 9), the fractional abundance is ∼(1–2) × 10⁻⁹, similar to those observed in Galactic center clouds: (2–8) × 10⁻⁹ (Minh et al. 1991; Armijos-Abendaño et al. 2015). On the other hand, it is orders of magnitude higher than those observed in Galactic disk clouds, which range from 10⁻¹³ to 5 × 10⁻¹¹ (Vastel et al. 2016; Fontani et al. 2018; Majumdar et al. 2018).

To estimate the gas-phase abundances of CO₂ from HOCO⁺, we ran chemical abundance models based on Nautilus (Ruaud et al. 2016), accounting for gas, ice surface, and ice mantle phases. In addition to the thermal evaporation, desorption from dust heating due to cosmic rays (Hasegawa & Herbst 1993) is included in the model. Desorption through the cosmic-ray heating of dust is where the dust grain is temporarily heated to a certain maximum temperature for a very short timescale (∼10⁻⁵ s), which then cools down. The model also includes photodesorption (Öberg et al. 2009a, 2009b) both from direct UV photons and cosmic-ray-induced UV photons with a default yield of 10⁻⁴ for all the grain species. We also ran a model with a desorption yield of 10⁻³ for CO₂. Our models do not include shocks. We calculated grid models with varying densities (n = 10³–10⁶ cm⁻³) and cosmic-ray ionization rates (ζ = 10⁻¹⁷−10⁻¹² s⁻¹) following a similar approach as Harada et al. (2021). Temperatures were calculated with the Meudon photodissociation region (PDR) code (ver. 1.5.4) (Le Petit et al. 2006) (Figure 6(a)), and were fed to Nautilus to run chemical abundance models with a larger chemical network. Despite the high gas temperature with the high cosmic-ray ionization rate (T > 1000 K when ζ = 10⁻¹² s⁻¹ and n = 10³ cm⁻¹), the dust temperature calculated from the Meudon code remains cold, around 11 K. We adopted a maximum visual extinction of A_V = 20 mag with a turbulent velocity of 1 km s⁻¹, and used the temperature in the model at A_V = 10 mag (Figure 6(d)), where the effects of the PDRs are negligible. We note that, unlike in the description of the observational results, the fractional abundances are expressed as abundances of certain species over the total hydrogen abundance (N_Htotal = N_Hatom + 2N_H2), instead of as molecular hydrogen abundances.

Zoom In Zoom Out Reset image size

Figure 6(b) shows the HOCO⁺/CO₂ abundance ratios in the gas phase, that vary between 10⁻⁵ and 10⁻² for the most part. Previous ALCHEMI studies have suggested that the cosmic-ray ionization rates in NGC 253 are at least a few orders of magnitude higher than that in Galactic spiral-arm clouds (Holdship et al. 2021; Harada et al. 2021; Holdship et al. 2022), which increases the HOCO⁺/CO₂ abundance ratios due to an increased ${{\rm{H}}}_{3}^{+}$ abundance. Although the cosmic-ray ionization rates are expected to be high, extremely high rates (ζ ≳ 10⁻¹³ s⁻¹ for n = 10⁵ cm⁻³) would destroy HOCO⁺ (Figure 6(c)), especially in lower-density regions. Therefore, we consider HOCO⁺/CO₂ ratios of ∼ 10⁻³–10⁻² in the latter sections, which are taken from the parameter space where the HOCO⁺ abundance is moderately high. ²³

It is worth noting that cosmic-ray-induced desorption mechanisms do not change the abundances significantly. This is because cosmic-ray-induced photodissociation is a more efficient form of destruction than desorption if the desorption yield is on the order of 10⁻³ or lower. According to Öberg et al. (2009b), the photodesorption yield of CO₂ cannot go higher than a few times 10⁻³, which means photodesorption is not significant. The desorption of CO₂ due to cosmic-ray heating of dust has an even lower effect than the photodesorption for models run with the commonly used maximum dust temperature of 70 K.

In our model, the dominant formation reactions of HOCO⁺ vary with time. In general, the protonation of CO₂ is more dominant at early times (<10⁵ yr), and the gas-phase production with HCO⁺ and OH becomes more efficient in later times (see Section 1). However, the dominant formation routes also vary with physical conditions and it is difficult to conclude which one is more dominant.

Although our models do not include shocks, but we argue that this approach should be sufficient to estimate the HOCO⁺/CO₂ ratios, especially their upper limits. The HOCO⁺/CO₂ ratios are determined by the balance among the protonation of CO₂, electron recombination of HOCO⁺, proton transfer from HOCO⁺ to species with a higher proton affinity than CO₂, and the ion-neutral production of HOCO⁺(HCO⁺ + OH). These reactions occur regardless of shocks. If shocks evaporate CO₂ significantly, there should be less contribution to HOCO⁺ formation through HCO⁺ and OH compared with the protonation of CO₂, and the HOCO⁺/CO₂ ratios should be lower, while shocks should increase the fractional abundances of both HOCO⁺ and CO₂. Therefore, our models without shocks are likely sufficient to obtain upper limits of the HOCO⁺/CO₂ ratios, but more realistic modeling with shocks will be conducted as future work. It should also be noted that the models without shocks severely underproduce the HOCO⁺ fractional abundances compared with the observed peak values, which implies the need for shocks to explain the observed abundances (see Section 6.1). We include the gas-neutral reaction of HCO⁺ + OH and many other related reactions in the model. The fact that our model could not reproduce the observed fractional abundances suggests that the formation route through this reaction is not enough to explain our observations.

6.1. Origins of the HOCO⁺ Emission

6.2. Inferred Gas-phase CO₂ Fractional Abundances

6.3. Comparison with Ice Observations of CO₂

Here we compare the fractional abundance of CO₂ gas estimated above with that of the CO₂ ice observations. We utilize the CO₂ column densities of selected regions of AKARI observations made by Yamagishi et al. (2015). Although this reference describes these observations in detail, we have included a summary of them in Appendix B. These observations were made with a rectangle slit with a size of 5'' × 58 (Figure 7, left). Subsequently, we extracted the values of the continuum flux from the same regions to estimate the total H₂ column density using the method described in Section 4. We then derived the fractional abundances of CO₂ in the ice phase in these regions, as shown in Figure 7 (right). We note that the components traced by CO₂ ice likely come from relatively lower-column-density regions than the ones traced by the dust continuum. Therefore, we have to be aware of the caveat that our estimations of the CO₂ fractional abundances are rather crude, only accurate for an order-of-magnitude approximation.

Zoom In Zoom Out Reset image size

The CO₂ ice fractional abundances near the center of NGC 253 (rectangular regions 2 and 3) are lower than in regions 4 and 5 by about an order of magnitude or more. Away from the center, the CO₂ ice fractional abundances are similar to those of the Milky Way ISM of ∼10⁻⁶–10⁻⁵ (Boonman et al. 2003), which are also consistent with the fractional abundances of CO₂ ice in our chemical model. Therefore, the CO₂ ice abundances decrease toward the center of NGC 253, deviating from the Milky Way value. We note that, if we use a higher dust temperature, we would derive lower H₂ column densities than those shown in Figure 4. Subsequently, the fractional abundance estimate would increase. Although warm dust does not fully explain the large difference in the derived fractional abundances between the outer part (rectangular regions 4 and 5) and center (rectangular regions 2 and 3) of the CMZ, this uncertainty should be considered in the interpretation of the data.

If the CO₂ ice fractional abundance is indeed lower in the center than at the outer CMZ, one of the possible factors that may contribute to this suppression is photodissociation, either directly by UV photons or cosmic-ray-induced UV photons. Because of the high star formation rate in the rectangular regions 2 and 3, strong photodissociation is expected (Meier et al. 2015).

Warm dust can also lower the ice abundance of CO₂ through desorbing CO₂ ice into the gas phase, or desorbing the precursors of CO₂ ice. Desorption of many CO₂ precursors can take place with lower dust temperatures than the desorption of CO₂ itself. CO₂ ice is thought to be formed through (Minissale et al. 2013)

$$\begin{eqnarray}&&{\rm{CO}}+{\rm{OH}}\longrightarrow {{\rm{CO}}}_{2}+{\rm{H}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{CO}}+{\rm{O}}\longrightarrow {{\rm{CO}}}_{2},\end{eqnarray} \tag{ 3 }$$

and possibly (Minissale et al. 2015)

$$\begin{eqnarray}&&{{\rm{H}}}_{2}{\rm{CO}}+{\rm{O}}\longrightarrow {{\rm{CO}}}_{2}+{{\rm{H}}}_{2}.\end{eqnarray} \tag{ 4 }$$

Some of these reactants have lower binding energies than CO₂(E_b (CO₂) ∼ 3000 K, E_b (CO) ∼ 1300 K, and E_b (O) ∼ 1600 K (Wakelam et al. 2017; Minissale et al. 2022). ²⁴ Higher dust temperatures increase the desorption rates of ice species, and species with lower binding energies can desorb with lower dust temperatures (35–50 K for CO, ∼50 K for O, and ∼80 K for CO₂; Minissale et al. 2022). If the desorption rate of any of these reactants is faster than the reaction rate to form CO₂ there would be less abundances of CO₂ ice. Although the mass of warm or hot dust is smaller than that of cold dust, these warm/hot components are likely concentrated in rectangular regions 2 and 3 in Figure 7 where the star formation is active. This can explain the suppression of ice in these regions as well as the low HOCO⁺ fractional abundances in GMCs 4 and 5 (see Figure 2(f)).

Cosmic rays could also contribute to ice desorption. We argued in Section 5 that cosmic-ray-induced desorption is likely negligible, but its effect could possibly be enhanced in NGC 253. Hasegawa & Herbst (1993) estimated that the maximum dust temperature due to the cosmic-ray heating of the dust is 70 K, but this temperature can be different for the case of NGC 253, where some dust is already warm. With a higher maximum dust temperature, the evaporation rate of ice species can be enhanced (e.g., Kalvāns & Kalnin 2020). This dependence of desorption due to cosmic-ray heating on the dust temperature should be further explored with theoretical studies.

Variations of the CO₂ ice abundance can be also caused by the initial ice composition. For example, the ice may be rich in atomic or molecular hydrogen. If the ice is abundant in atomic H, frequent hydrogenation reactions can occur, and species such as water, CH₄, NH₃, and CH₃OH may be abundant. On the other hand, if atomic H is deficient, CO₂ formation may be a more dominant route of destroying CO ice than CH₃OH formation. Although it is difficult to assess how much the ice compositions differ between Galactic star-forming regions and the NGC 253 CMZ, we note that this is another factor that could affect the overall chemistry.

6.4. Ice and Gas-phase Chemistry in NGC 253

In this paper, we analysed the abundances of HOCO⁺, the protonated form of CO₂, in the CMZ of the starburst galaxy NGC 253, and discussed its relationship with the gas- and ice-phase CO₂. Below is a summary of our findings.

• 1.  
  The distribution of HOCO⁺ shows clear enhancements at the locations of the x₁ and x₂ orbital intersections where shocks are expected. This distribution is similar to that of methanol but is different from that of H¹³CO⁺. There are two formation routes of HOCO⁺; one is the ion-neutral reaction HCO⁺ + OH and the other is the protonation of CO₂. If the former route is dominant, the HOCO⁺ distribution should be similar to that of HCO⁺, while the latter route should cause similarity with the CH₃OH distribution. Therefore, HOCO⁺ is likely produced through the protonation of CO₂.
• 2.  
  We derive HOCO⁺ column densities across the CMZ using CASSIS, from which we also obtain its fractional abundances using the total H₂ column densities estimated from the dust emission. We find HOCO⁺ fractional abundances as high as ∼2 × 10⁻⁹, which are similar to those observed in the Galactic center, but orders of magnitude higher than those reported in Galactic spiral-arm molecular clouds.
• 3.  
  From the results of chemical modeling and the values of the cosmic-ray ionization rates derived from previous ALCHEMI works, we estimate that the gaseous HOCO⁺/CO₂ ratio is likely 10⁻³–10⁻². This ratio suggests that the gas-phase CO₂ fractional abundances are (1–20) × 10⁻⁷ at the peaks of HOCO⁺ emission.
• 4.  
  We also estimate fractional abundances of CO₂ ice from their column densities in the literature. The ice fractional abundance at the HOCO⁺ peak is similar to the value in the Galactic interstellar medium (10⁻⁶−10⁻⁵), but is lower ( ∼(1–3) × 10⁻⁷) near the NGC 253 galactic center.
• 5.  
  The increased gaseous and ice fractional abundances of CO₂ at the outer CMZ of NGC 253 imply that a large fraction of ice is sublimated. Because of the association of these locations with the evidence of shocks, we propose that this efficient sublimation is attributed to shock-induced sputtering.

High-spatial-resolution observations of molecular emission in external galaxies, such as those performed by the ALCHEMI survey toward the central regions of the starburst galaxy NGC 253 with ALMA, have greatly improved our understanding of gas-phase abundances and will continue to do so. Now, complementary observations of ice at high angular resolutions with the JWST are required to obtain a complete picture of the chemical processes in starburst galaxies.

We thank the anonymous referee for constructive comments. N.H. thanks Hideko Nomura for the helpful discussion on CO₂ transitions at infrared wavelengths, and Kotomi Taniguchi for the initial help using CASSIS. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00161.L, ADS/JAO.ALMA#2018.1.00162.S, ADS/JAO.ALMA#2018.1.01321.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Data analysis was in part carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan. N.H. acknowledges support from JSPS KAKENHI grant No. JP21K03634. V.M.R. has received support from the Comunidad de Madrid through the Atracción de Talento Investigador Modalidad 1 (Doctores con experiencia) grant (COOL:Cosmic Origins of Life; 2019-T1/TIC-5379), and the Ayuda RYC2020-029387-I funded by MCIN/AEI/10.13039/501100011033. L.C. has received partial support from the Spanish State Research Agency (AEI; project number PID2019-105552RB-C41). K. N. acknowledges support from JSPS KAKENHI grant No. 19K03937. P.H. is a member of and received financial support for this research from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. K.S. is supported by the grant MOST 111-2112-M-001-039 from the Ministry of Science and Technology, Taiwan. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.

Facility: ALMA - Atacama Large Millimeter Array.

Software: Astropy (Astropy Collaboration et al. 2013, 2018), CASA (McMullin et al. 2007), CASSIS (Vastel et al. 2015), Nautilus (Ruaud et al. 2016).

CASSIS constrains parameters such as the column density, line velocity, line width (FWHM), excitation temperatures, and source sizes by fitting the spectra. We used the MCMC algorithm provided by CASSIS to derive these parameters. We used the transitions in the spectra, as shown in the upper half of Table 1. When using the MCMC method, users provide acceptable ranges of these parameters as well as initial guesses. We used the column-density range of [10¹², 10¹⁶] cm⁻² with an initial guess of 10¹⁴ cm⁻², the excitation temperature range [5,50] K with an initial guess of 6 K, and the velocity range [v_CO − 30,v_CO + 30. ] km s⁻¹, where v_CO is the velocity obtained from the moment 1 image of CO(1-0), which is also used for the initial guess. The range of FWHM we used is [σ_V -20.,90.] km s⁻¹ with an initial guess of σ_V , where σ_V is the velocity dispersion of the CO(1-0) image. ²⁵ We note that the line width of HOCO⁺ is significantly different from that of CO(1-0), and this range of FWHM is simply determined by running CASSIS multiple times and checking the range to produce a reasonable fit.

Uncertainties from the CASSIS fit are reasonably small for most cases. Figure 8 (left) shows that the errors of the column densities are around 10% for most cases, and ∼20%–30% for a small fraction of pixels with low signal-to-noise ratios. The uncertainties in the excitation temperatures are within 0.5 K for most pixels, with a maximum of 3 K (Figure 8, right).

Zoom In Zoom Out Reset image size

There are also other sources of uncertainties, in addition to the spectral fitting. For example, the HOCO⁺ column densities derived with the spectroscopic constants from the Jet Propulsion Laboratory (JPL ²⁶ ) would yield up to a factor of 2 larger values compared with ones from CDMS. Although this is a large factor, the use of different spectroscopic constants changes the results uniformly within the field of view. There are also observational uncertainties of up to 15% (Martín et al. 2021). These uncertainties will not change our main conclusions.

The HOCO⁺ transitions are optically thin even where the column densities are high. CASSIS does not provide the optical depths, but we also ran MADCUBA (Martín et al. 2019), a similar spectral fitting software, to obtain optical depths. The optical depths are <0.1 for positions that we checked, which have high intensities of HOCO⁺.

The modeled intensities with CASSIS for pixels that are closest to each GMC position are shown in Table 2 if a HOCO⁺ transition for a GMC is detected.

Here we summarize the analyses of Yamagishi et al. (2015), and present spectra for the regions used in our analysis. The observed wavelength range is about 2.5–5.0 μm. Within this range, CO₂ ice, H₂O ice, Brα, and PAH 3.3 μm features were detected in addition to the continuum. The CO₂ ice features at 4.27 μm are fit using the data range 4.1–4.4 μm as there is only this narrow ice feature in this wavelength range. Because the ice composition changes the spectral shape, multiple ice compositions were tested to best fit the spectra. Consequently, an ice composition of H₂O:CH₃OH:CO₂= 9: 1: 2 was used in the final analysis.

Figure 9 shows the spectra used to derive the CO₂ ice column densities shown in Figure 7. Note that our Regions 1–5 correspond to IDs 47–51 in Table 3 of Yamagishi et al. (2015).

As already noted by Humire et al. (2022) and Behrens et al. (2022), the GMC nomenclature was adopted from Leroy et al. (2015), but with modified positions. How these modifications are made is explained in Behrens et al. (2022). These positions are shown in Table 3.

The large-scale CO(2-1) image of NGC 253 shown in Figure 1 was taken from the ALMA archive (project code #2018.1.01321.S). These data were obtained from a configuration consisting of the 7 m array complemented by the total power antenna. Pipeline-reduced image cubes (QA2 products) for the 7 m array and total power single-dish data were combined with the CASA command feather. This image was shown to indicate the rough positions of the x₁ and x₂ orbits only. We expect that the PI team will present the data with better imaging quality and scientific analysis.

Zoom In Zoom Out Reset image size