Shock and Cosmic-Ray Chemistry Associated with the Supernova Remnant W28

Supernova remnants (SNRs) exert a strong influence on the physics and chemistry of the nearby molecular clouds through shock waves and the cosmic rays (CRs) they accelerate. To investigate the SNR–cloud interaction in the prototype interacting SNR W28 (G6.4−0.1), we present new observations of the HCO+, HCN , and HNC J = 1–0 lines, supplemented by archival data of CO isotopes, N2H+ and H13CO+. We compare the spatial distribution and spectral line profiles of different molecular species. Using the local thermodynamic equilibrium assumption, we obtain an abundance ratio of N(HCO+)/N(CO) ∼10−4 in the northeastern shocked cloud, which is higher by an order of magnitude than the values in unshocked clouds. This can be accounted for by the chemistry jointly induced by shock and CRs, with the physical parameters previously obtained from observations: preshock density n H ∼ 2 × 105 cm−3, CR ionization rate ζ = 2.5 × 10−15 s−1, and shock velocity V s = 15–20 km s−1. Toward a point outside the northeastern boundary of W28 with a known high CR ionization rate, we estimate the abundance ratio N(HCO+)/N(N2H+) ≈ 0.6–3.3, which can be reproduced by a chemical simulation if a high density n H ∼ 2 × 105 cm−3 is adopted.


INTRODUCTION
When massive stars end their lives as energetic supernova (SN) explosions, many of them have not yet moved far away from their natal molecular clouds (MCs).The physical and chemical conditions in the MCs could be strongly modified by the energy feedback of the supernova explosions and the following supernova remnant (SNR) evolution.Many SNRs exhibit evidence of interaction with nearby MCs (e.g.Seta et al. 1998;Jiang et al. 2010;Kilpatrick et al. 2016;Zhou et al. 2023).These interactions are crucial to understanding the feedback of SNe and SNRs into the ISM.
The shock waves of SNRs propagating into dense MCs will compress and heat the gas (Draine & McKee 1993).The chemistry in MCs is sensitive to temperature and density, and thus sensitive to the passage of the shock wave.The shock chemistry has been extensively studied in shocked regions of protostellar outflow (e.g.Poygchen@nju.edu.cndio et al. 2014;Lefloch et al. 2017).There are also a few studies concerning the chemistry of SNRs.For example, van Dishoeck et al. (1993) presented a comprehensive discussion of the chemistry in the SNR IC443.Lazendic et al. (2010) found an enrichment of HCO + and SO in the SNR G349.7+0.2.Zhou et al. (2022) found high abundance ratio between HCO + and CO in the SNR W49B and attributed it to shock wave, cosmic ray (CR) ionization, and possibly the X-ray precursor of the shock.Mazumdar et al. (2022) revealed formation of H 2 CO and CH 3 OH in the shocked region of the SNR W28.However, observational study of shock chemistry is still scarce and more investigations are required.
Meanwhile, SNR shocks are believed to be the prime accelerators of Galactic CRs following the diffuse shock acceleration process (Bell 1978).Many SNRs are found to be located close to γ-ray sources (e.g., Aharonian 2013, and references therein).The high-energy emissions originate from either the decay of π 0 mesons produced by collision between high-energy (≳ 280 MeV) CR protons with H nuclei in MCs, which is the so-called hadronic scenario, or the inverse Compton scattering of background radiation by high-energy electrons, which is the so-called leptonic scenario.In the case that the SNR is associated with MCs, the hadronic scenario can naturally explain the γ-ray emission if the observed γ-ray spectra can be reasonably fitted (e.g.Aharonian et al. 2008;Fukui et al. 2021).
Low-energy CR protons will not contribute to γ-ray emission, but they are the dominating ionization sources in dense MCs into which ionizing UV emission is difficult to penetrate (Padovani et al. 2009).CR protons ionize molecular hydrogen gas mainly through: followed by The H + 3 ion starts various chemical reactions in MCs.Chemical simulations have shown that CR ionization exerts strong influence on the chemistry of MCs (e.g.Albertsson et al. 2018).The usability of chemical codes and networks requires observational tests.However, few observations have been done concerning the chemical effects of high CR ionization rate, especially in the MCs associated with SNRs.
SNR W28 is one of the prototypes for studying SNRcloud interaction, with an estimated age of ∼ 3.3-4.2× 10 4 yr (Velázquez et al. 2002;Rho & Borkowski 2002;Li & Chen 2010) at a distance of ∼ 1.9 kpc (Velázquez et al. 2002;Ranasinghe & Leahy 2022).Located in a complex of MCs, W28 is believed to be interacting with MCs.Evidence for the interaction includes molecular line broadening spatially coincidence with the radio continuum and TeV γ-ray emission (Wootten 1981;Arikawa et al. 1999;Reach et al. 2005;Nicholas et al. 2011Nicholas et al. , 2012;;Gusdorf et al. 2012;Maxted et al. 2016;Mazumdar et al. 2022), 1720 MHz OH masers (Claussen et al. 1997;Hewitt et al. 2008), infrared H 2 emissions (Reach & Rho 2000;Reach et al. 2005;Yuan & Neufeld 2011), etc.In the γ-ray band, GeV and TeV emissions are found associated with W28 (Aharonian et al. 2008;Cui et al. 2018), which are expected to originate from the hadronic scenario.The overall TeV emission shows two parts: HESS J1801−233 towards the northeastern MCs, and HESS J1800−240 towards the southern MCs outside the boundary of W28.The latter is further divided into three components: A, B and C, all of which are coincident with dense MCs.The γ-ray emission outside the radio boundary is ascribed to the illumination by the escaped CRs from W28 (e.g., Li & Chen 2010;Cui et al. 2018).Vaupré et al. (2014) estimated the CR ionization rate with the abundance ratio N (DCO + )/N (HCO+) in some regions to be ≳ 10 15 s −1 , which is higher than the standard value by two orders of magnitude (Glassgold & Langer 1974).This makes W28 the third SNR directly measured to exhibit high CR ionization rate after IC443 (Indriolo et al. 2010) and W51C (Ceccarelli et al. 2011).Besides, the existence of 1720 MHz OH masers and class I methanol CH 3 OH masers Pihlström et al. (2014) also requires high CR ionization rate (Nesterenok 2022).In addition, discovery of the Fe I Kα line at 6.4 keV (Nobukawa et al. 2018;Okon et al. 2018) further indicates the ionization induced by low energy CRs.All of the mentioned observations render W28 an ideal site to study the chemistry of SNR shock and CR ionization.
The J = 1-0 lines of HCO + , HCN and HNC molecules are typical tracers of dense gas (e.g.Shirley 2015).The abundance ratio between HCO + and other specific molecules, for example, CO (e.g., Zhou et al. 2018;Bisbas et al. 2023), has been proposed to be a tracer of CR ionization (e.g.Albertsson et al. 2018), while the line ratio I(HCN)/I(HNC) serves as a tracer of gas kinetic temperature (e.g.Hacar et al. 2020).Therefore, observation of these three lines towards W28 provides information of the physical and chemical effects of W28 on its adjacent MCs.
In this paper, we present new HCO + , HCN and HNC J = 1-0 observations towards W28.Using the new data together with the archival data and chemical models, we investigate the chemistry brought by shock wave and CRs.The paper is organized as follows.In Section 2, we describe the new observations and archival data.In Section 3, we present the results of the observations.We derive the abundance ratios between molecular species, and then present the results and relevant discussions of chemical simulations in Section 4. The conclusions are summarized in Section 5.

HCO + , HCN and HNC observations and data reduction
Observations of the J = 1-0 emission lines of HCO + , HCN and HNC were performed with the 13.7m millimeter-wavelength telescope of the Purple Mountain Observatory at Delingha (hereafter PMOD; PI: Yang Chen).The HCO + and HCN lines were simultaneously observed with the SIS receiver in two epochs during 2018 June and 2021 August-October.These observations used on-the-fly (OTF) mapping to cover the northeastern MCs of W28 with a mapping area of 28 ′ × 29 ′ and the southern MCs of W28 with two mappings, each with an area of 40 ′ ×45 ′ and 21 ′ ×24 ′ .HNC observations were conducted separately in 2021 June-August to cover the northeastern MCs of W28 with a mapping area of 28 ′ ×27 ′ .Positions and coverage of the regions are shown in Figure 1.The Fast Fourier Transform Spectrometers with 1 GHz bandwidth and 16,384 channels were used as the back ends, providing a velocity resolution of 0.21 km s −1 at 89 GHz.The half-power beamwidth (HPBW) of the telescope at 89 GHz is ≈ 60 ′′ .The main beam efficiency was around 0.60 in 2018 and 0.58 in 2021 1 .The typical system temperature T sys is between 150 and 310 K, depending on the weather condition and the elevation of the source.The pointing accuracy of the antenna better than 5 ′′ .The raw data were reduced with the GILDAS/CLASS package2 .The data cubes of HCO + , HCN and HNC were all resampled to have the same velocity channel width of 0.25 km s −1 and the same pixel size of 30 ′′ .The RMS noise is ∼ 0.1 K for HCO + and HCN and ∼ 0.05 K for HNC.

Other archival data
We used other archival data to support our analysis for a multiwavelength view.We obtained 12 CO, 13 CO and C 18 O J = 1-0 data from the Milky Way Image Scroll Painting (MWISP) line survey project.The HPBW is about 50 ′′ at 115 GHz and the pixel size is 30 ′′ .The velocity channel width is 0.16 km s −1 for 12 CO and 0.17 km s −1 for 13 CO and C 18 O, and the typical noise measured in T mb is 0.5 K for 12 CO and 0.3 K for 13 CO and C 18 O.Detailed description of the project can be found in Su et al. (2019).
We also retrieved the Very Large Array (VLA) 327 MHz continuum data from the MAGPIS (The Multi-Array Galactic Plane Imaging Survey) website3 .The 13 CO J = 2-1 data cube was taken from the SEDIGISM (Structure, Excitation and Dynamics of the Inner Galactic Interstellar Medium, Schuller et al. 2021) survey, which is observed with the APEX telescope with an angular resolution of 30 ′′ and a sensitivity of 0.8-1.0K at 0.25 km s −1 velocity resolution.Supplementary data cubes of H 13 CO + and N 2 H + were taken from the MALT90 survey (Millimetre Astronomy Legacy Team 90 GHz Survey, Jackson et al. 2013) with angular resolution of ≈ 38 ′′ and a pixel size of ≈ 9 ′′ .The typical sensitivity measured in antenna temperature (T * A ) is 0.2-0.25 K at a velocity resolution of 0.11 km s −1 .The antenna temperature is transferred to main beam temperature T mb with a main beam efficiency of 0.49.The Herschel column density map was obtained from Marsh et al. (2017) who fitted the data of Hi-GAL (Herschel infrared Galactic Plane) survey in 70-500 µm with the PPMAP procedure (point process mapping, Marsh et al. 2015).The angular resolution of the map is 12 ′′ .
All the processed data were further analyzed with Python packages Astropy (Astropy Collaboration et al. 2022) and Spectral-cube (Ginsburg et al. 2015).The data cubes were reprojected with Montage4 package when necessary.We visualized the data with Python package Matplotlib5 .

Northeastern molecular clouds/HESS J1801−233
In Figure 2, we show the integrated intensity map of HCO + from −30 km s −1 to +40 km s −1 in the northeastern (NE) region overlaid with VLA 327 MHz continuum and HESS TeV γ-ray emission.The range of localstandard-of-rest (LSR) velocity is chosen according to the velocity span of the shocked HCO + emission.This velocity interval is consistent with the broadened line profile observed with other molecular transitions (Reach et al. 2005;Nicholas et al. 2011Nicholas et al. , 2012)).The broad line emission of HCO + is generally aligned with the shock front traced by radio continuum, which is solid evidence of shock-cloud interaction.Notably, the HCO + emission also appears essentially coincident with the TeV emission, which further supports the hadronic scenario of the γ-ray emission and provides a hint for CR ionization.The main body of the HCO + emission is in the center of the figure, while small portions of emission are also found in the northwestern (surrounding the re- gion N) and northeastern (outside the radio emission) of the figure.The emission outside the radio continuum is spatially coincident with a massive star-forming region, G6.796−0.256(hereafter G6.796 for short), in the AT-LASGAL survey (Urquhart et al. 2018), and is observed by the MALT90 survey.We will discuss this region in Section 3.2.
The main body of the HCO + emission shows bright clumps coincident with the brightest part of the radio continuum spatially covering Region 1.The integrated intensity decreases towards the northeast.
Although the 1720 MHz OH masers are generally within the spatial extent of the HCO + emission, most of them are not coincident with the brightest part of the HCO + emission.The HCO + emission found towards W28F, the northeastern part of the radio continuum where a cluster of OH masers is located, is not so bright as the emission around Region 1.However, observations of Gusdorf et al. (2012) show bright high-J (up to J up = 11) CO emission towards this region, indicating strong shock disturbance and high temperature.The weak HCO + emission may be attributed to the higher energy level population due to the high temperature, or beam dilution since the beam of PMOD is larger than that of Gusdorf et al. (2012).
In Figure 3 we show the Herschel column density map as well as the integrated intensity maps of HCN, HNC, .Herschel column density map as well as integrated intensity of HCN J = 1-0, HNC J = 1-0, 12 CO J = 1-0, 13 CO J = 1-0, and 13 CO J = 2-1 in the "NE" region.The velocity range of integration is −30 km s −1 to +40 km s −1 for HCN, HNC and 12 CO while −15 km s −1 to +30 km s −1 for 13 CO J = 1-0 and J = 2-1 lines.Pixels with integrated intensity weaker than 5σ are masked out.The light brown contours show the VLA 327 MHz continuum emission (levels are 60 and 180 mJy beam −1 ).
12 CO and 13 CO J = 1-0 and 13 CO J = 2-1.Morphologies of the HCN and HNC emission are both similar to that of the HCO + emission.The HNC emission is also detected in the westernmost part of the figure (coincident with HII region W28A, Goudis 1976) and between the main body and G6.796.The non-detection of HCO + and HCN in these regions might be attributed to lower sensitivity of the detector.The spatial distribution of 12 CO is also generally similar to that of HCO + , although 12 CO is much stronger and more extended.Emission of 13 CO J = 1-0 is stronger outside the radio continuum than inside, which is different from the case of HCO + .The reason may be that 13 CO J = 1-0 is often optically thin and traces the unshocked cloud (e.g.Su et al. 2011).
The 13 CO J = 2-1 data has fine angular resolution and shows detailed structure of the molecular clouds.Several clumpy structures are seen towards the margin of the radio continuum.The distribution of 13 CO J = 2-1 is broadly coincident with the Herschel column density map.
To further investigate the properties of the shocked molecular clouds, we extract the spectra of the regions marked in Figure 2, namely Regions 1 to 5 and Region N. The Regions 1 to 5 are located towards the outside of the SNR to show the spatial variation of the shock interaction.The spectra are shown in Figure 4. Significant line broadening can be seen in the HCO + and HCN line profiles of Regions 1 to 4 and in Region N, which originate from the disturbance of the SNR blast wave.The linewidth of HCO + decreases from Region 1 to Region 5, indicating the disturbance effect of the shock is weakened towards the outer part of the SNR.In Regions 1 and 2, an absorption feature is in the center of the HCO + spectrum, where lines of HNC and 13 CO show pure emission features.This is caused by the cold and dense molecular gas in the line of sight (Reach et al. 2005).In Regions 3 and 4, the HCO + spectra exhibit double-peak features at ∼ +6 and ∼ +12 km s −1 (see the dotted vertical lines in Figure 4), with their blue sides extending to ≤ −3 km s −1 (see the dashed vertical lines in Figure 4).
The spectra of HCN are generally similar to the HCO + spectra in the shocked regions (1, 2 and N).The hyperfine structures of the HCN J = 1-0 line, at −7.1 and +4.9 km s −1 offset from the strongest line in the center, can be seen in Regions 3, 4 and 5, but are covered up by strong line broadening in Regions 1, 2 and N. .Spectra of HCO + J = 1-0 (orange), HCN J = 1-0 (black), HNC J = 1-0 (cyan), 12 CO J = 1-0 (grey), 13 CO J = 1-0 (magenta) and 13 CO J = 2-1 (green) lines averaged in the regions marked in Figure 2.Each region takes up two vertical figures to fully display the spectra.The main beam temperature T mb of HCO + , HCN and HNC is multiplied by a factor of 10 for better display.The dashed vertical lines in the figures of Regions 1 to 5 mark a velocity of −3 km s −1 , while the dotted vertical lines in the figures of Regions 3 to 5 mark a velocity of +6 and +12 km s −1 (see the text for details).
Similar to HCO + and HCN, HNC is also a dense gas tracer, but the line profiles of HNC are rather different in some of the regions.In Regions 1 and N, the line profiles of HNC are similar to those of HCO + and HCN, although they are relatively weak.In Region 2, HNC does not show absorption as HCO + and HCN do.In Regions 3 to 5, the HNC emission has only one peak at ∼ +6 km s −1 but still shows non-Gaussian line profile with a sharp peak.
Although the 12 CO lines also exhibit large linewidth, several overlapping components are present and prohibit us from identifying the shocked component.However, the 12 CO spectra in Regions 1 to 4 show strong emission at −3 km s −1 (see the dashed vertical lines in Figure 4), and their line profiles resemble that of HCO + around −3 km s −1 .These features of the 12 CO spectra are directly related to shock perturbation.In Region N, the line broadening of 12 CO is more prominent and is consistent with that of HCO + .Both 13 CO J = 1-0 and J = 2-1 are found in all of the selected regions.Strong unshocked 13 CO is found in Region 3, which is also consistent with what we present in Figure 3.The 13 CO molecules surrounding SNRs are commonly considered as a tracer of unshocked gas in SNRs because of their small optical depth compared with 12 CO, but the line profiles of 13 CO in Regions 1 and 2 are moderately broadened.Since Nicholas et al. (2012) have found moderately broadened C 34 S line close to the regions we select, it is possible that 13 CO can show broadened line profile originated from shock perturbation, since both C 34 S and 13 CO are rare isotopes and can trace dense gas.

ATLASGAL clump G6.796
ATLASGAL clump G6.796 is a massive star-forming region (Urquhart et al. 2018) located to the northeast of the radio shell of W28 (see Figure 2).Lefloch et al. (2008) proposed that its star-forming activity might be triggered by W28.The MALT90 survey has detected significant line emission of HCO + , H 13 CO + , N 2 H + , HCN, HNC, C 2 H and HC 3 N from this clump.We plot the integrated intensity map of N 2 H + line emission overlaid with contours of the Herschel column density map in Figure 5.The clump is highly compact and barely resolved by the beam of Mopra (∼ 38 ′′ ).But it is clear that the peak of the molecular emission is offset from the peak of the dust emission represented by G6.796.2014) obtained high CR ionization rate (≳ 10 −15 s −1 ).Among the three points, only N6 shows prominent HCO + and N 2 H + line emissions.The spectra of the HCO + , H 13 CO + and N 2 H + lines averaged in a 27 ′′ × 27 ′′ region towards N6 and the results of Gaussian fitting are shown in Figure 6.We report only a marginal detection of H 13 CO + , which is inconsistent with the result of Vaupré et al. (2014) who obtained a peak temperature T peak of 0.53 K for H 13 CO + .This may result from the different beam sizes of the MALT90 data (38 ′′ ) and the IRAM 30m data of Vaupré et al. (2014) (29 ′′ ), as well as the relatively low sensitivity of the MALT90 data.The T peak of H 13 CO + in the MALT90 data is ≈ 0.33 K.If the inconsistency of is due to the beam dilution effect, we expect that the T peak detected by IRAM 30m should be approximately 0.33 × (38/29) 2 = 0.57 K, which is roughly consistent with the result of Vaupré et al. (2014).The spectrum of the N 2 H + emission shows clear hyperfine structures.

Southern molecular clouds/HESS J1800−240
In Figure 7, we show the integrated intensity map of the HCO + emission from 0 km s −1 to +30 km s −1 in the southern molecular clouds overlaid with the HESS TeV γ-ray emission.This velocity interval covers all the HCO + emission in this region.Clumpy structures are found widespread in the extended region and are marked with orange boxes.Except for clump G5.7, all other molecular clumps revealed by our HCO + line observation are coincident with the molecular cores found in Nicholas et al. (2011), and thus we follow their nomenclature.All of the marked clumps including G5.7 are related to infrared sources found by the ALTAS-GAL survey (Urquhart et al. 2018).Cores 4 and 4a are probably related to HII regions G6.225−0.569and G6.1−0.6, respectively (Nicholas et al. 2011).Triple cores SE and Cen are coincident with ultra-compact HII regions G5.89−0.39Aand B, in which G5.89−0.39Bharbours an explosive outflow confirmed by the ALMA high-resolution observation (Zapata et al. 2020).Clump G5.7 is classified as a young stellar object in Urquhart et al. (2018) using infrared observations.However, it is also spatially coincident with a supernova remnant, G5.7−0.1 (Brogan et al. 2006), and an 1720 MHz OH maser (Hewitt & Yusef-Zadeh 2009), which definitely shows shock-cloud interaction.
In Figure 8 we show the spectra of HCO + J = 1-0, HCN J = 1-0, 12 CO J = 1-0, 13 CO J = 1-0, and 13 CO J = 2-1 lines in core 4, triple core SE, triple core Cen, core 5 SW and clump G5.7.A significant contrast to the spectra in the northeastern molecular clouds is that  all of the spectra in the selected regions show relatively narrow line width.The hyperfine structures of the HCN line are clear in all of the regions.In the 12 CO spectra of some regions, a strong peak with T peak > 20 K stands out with weaker emission in its blue and red sides.We tend to attribute them to irrelevant unshocked components instead of line wings because most of them show counterparts in 13 CO J = 1-0 lines.In all of the regions, emission of C 18 O is found.
In triple core Cen, significant asymmetric line profiles of 12 CO and HCO + can be seen in the blue side of the main peak.It is likely to be associated with the explosive outflow of G5.89−0.39B.

Line ratios
Molecular line ratio diagnostic is often used to investigate the physical and chemical condition of the ISM.We herein show the line ratio maps of the line ratios Line ratio I( 13 CO J = 2-1)/I( 13 CO J = 1-0) is also a tracer of gas temperature since it involves two transition lines of one molecule.We present the I( 13 CO J = 2-1)/I( 13 CO J = 1-0) line ratio map in Figure 10.An enhancement of the line ratio is shown in the shocked cloud in the northeastern molecular clouds.However, the spatial distribution of I( 13 CO J = 2-1)/I( 13 CO J = 1-0) is different from that of I(HCN)/I(HNC).Relatively high value is also found in the UC-HII region G5.89−0.39Aand B, i.e.Triple core SE and Cen as la-beled in Figure 7, which might arise from the heating of radiation and the explosive outflow.
The line ratio I( 12 CO J = 2-1)/I( 12 CO J = 1-0) has long been used as a tracer of shock heating (e.g.Seta et al. 1998;Jiang et al. 2010).However, I( 13 CO J = 2-1)/I( 13 CO J = 1-0) can act as the same role.This diagnostic provides another usable way to identify SNRcloud interaction and other heating processes, which could take the advantage of the MWISP (Su et al. 2019) and SEDISISM (Schuller et al. 2021) surveys towards the inner Galaxy.But we caution that I( 13 CO J = 2-1)/I( 13 CO J = 1-0) also depends on density and other excitation effects.Morphological agreement between the shock front and the enhanced line ratio should be taken into consideration when using this evidence.

Column densities and Abundance ratio between
12 CO and HCO + The abundance ratio between HCO + and 12 CO (hereafter CO if not specifed) is a tracer of the CR ionization rate (e.g., Zhou et al. 2022;Bisbas et al. 2023).To further test this notion in the SNR W28, we calculate the column density of HCO + and CO and their abundance ratio.In the shocked region of the northeastern MCs, we assume the emission of HCO + and CO is optically thin because of their large linewidths (which in turn mean large velocity gradients).We assume local thermodynamic equilibrium (LTE) in our calculation.For an optically thin transition, the column density of a molecular species is (Mangum & Shirley 2015): where Q rot ≈ kT /hb + 1/3 is the rotational partition function, J ν (T ) = (hν/k)/[exp (hν/kT ) − 1] is the Rayleigh-Jeans equivalent temperature, and f , the beam filling factor, is assumed here to be unity.Other parameters were explained in detail by Mangum & Shirley (2015).Under the LTE assumption, the excitation tem-T ex is equal to the kinetic temperature T k .We adopt the kinetic temperature estimated by Nicholas et al. (2011) and Maxted et al. (2016) with NH 3 observations.NH 3 is a good thermometer of dense molecular gas (Ho & Townes 1983).The molecular constants are taken from the Cologne Database for Molecular Spectroscopy (CDMS)6 .
To distinguish the emission of the shocked cloud from the intricate overlapping spectra, especially for the case of 12 CO, we carry out multi-Gaussian fitting of the spectra.We present the fitting results in Table 1 and Appendix A. Our fitting is based on the assumptions that (1) each 13 CO J = 1-0 component should have a corresponding component of 12 CO, (2) the shocked 12 CO component does not have 13 CO counterpart, and (3) the spectrum of the shocked 12 CO component can be approximated as a Gaussian.Although assumptions (1) and ( 2) hold true in most cases, we note that assumption (3) is just an expedient, which may bring extra uncertainty.
Among the shocked regions in the northeastern MCs, we choose Regions 3 and N for typifying the fitting because both of them show clear shocked components in the HCO + spectra with FWHM ≳ 20 km s −1 .The region sample is small but representative, while it is hard to decompose the 12 CO emissions in other regions.We include 7-8 components to account for the 12 CO spectra, which is reasonable because all of the assumptions mentioned in the previous paragraph are satisfied.In Region 3, we also notice an unshocked component of HCO + at +6 km s −1 (see Figure 14).We regard it as preshock gas and calculate the column density (shown in Table 1).Since H 13 CO + data is not available, we assume that the HCO + emission is optically thin despite its small linewidth, and hence the derived column density would be a lower limit.We use C 18 O emission to determine the column density of the unshocked CO molecules other than using 12 CO emission that is optically thick.
For the unshocked molecular clump in the southern MCs, the optically thin assumption may not hold true for the 12 CO and HCO + lines.In these cases, we calculate the column densities of their isotopic counterparts, i.e.C 18 O and H 13 CO + , which are expected to be optically thin.We assume the isotope ratios 12 C/ 13 C = 50 (Milam et al. 2005) and 16 O/ 18 O = 500 (Vaupré et al. 2014).We conduct multi-Gaussian fitting to the C 18 O and H 13 CO + spectra of the unshocked clouds.The results are presented in Table 1 and Appendix A. Among all the unshocked molecular clump in the southern MCs, we do not fit the spectra in the triple core SE and G5.7.In the triple core SE, the peak of the C 18 O line com- Note-a The kinetic temperature of the broad components of Region 3 is taken from Maxted et al. (2016), while the temperature of the narrow component of Region 3 is just an assumption.The temperature of Region N is set equal to that of the Region 3 broad component because of their similar physical conditions.The kinetic temperature used for the southern MCs are taken from Nicholas et al. (2011).Under the LTE assumption, T k = Tex is used to derive the column density.b Column density ratio between HCO + and 12 CO.
prises two components but they are difficult to decompose (see Figure 8).In G5.7, no H 13 CO + data is available.In the other regions, the components are clearly distinguishable (see Figure 15).The resulting column densities and abundance ratios are listed in Table 1.We note that it is almost impossible to determine the uncertainty brought by the LTE assumption due to some factors such as unknown gas density, so we do not give an uncertainty estimation here.Generally, the abundance ratio N (HCO + )/N (CO) is of the order of 10 −5 or lower in the unshocked MCs, which is consistent with the value found in cold MCs (10 −6 -10 −4 ) (e.g.Agúndez & Wakelam 2013;Miettinen 2014;Fuente et al. 2019).In the shocked clouds where the HCO + and 12 CO lines show large linewidths, N (HCO + )/N (CO) of order 10 −4 .Although this value is not significantly higher than the normal values, it is an order of magnitude higher than that in unshocked clouds, which indicates that there is an enhancement of HCO + /CO in the shocked molecular clouds of W28.
An important formation route of HCO + is: in which H + 3 is the key product of CR ionization.Therefore, the abundance of HCO + is sensitive to the CR ionization rate.In comparison, the abundance of CO does not vary significantly with CR ionization rate (Albertsson et al. 2018;Zhou et al. 2022).In the northeastern MC of W28, a high CR ionization rate is detected by Vaupré et al. (2014).There may be a link between the high CR ionization rate and the enhancement of N (HCO + )/N (CO).Also, since the enhancement of N (HCO + )/N (CO) is found in shocked clouds, we expect that shock chemistry may play an important role.
To further investigate the joint chemical effect of shock wave and CRs, we apply the Paris-Durham magnetohydrodynamic (MHD) shock model (Flower & Pineau des Forêts 2003, 2015;Godard et al. 2019) to simulate the chemical processes in MCs subjected to an MHD shock wave.The model is a 1D plane-parallel MHD model coupled with more than 1000 chemical reactions of over 100 species.This model has been used by Gusdorf et al. (2012) to study the multi-level CO emissions in W28F.The entire procedure of calculation consists of two steps.First, a steady state condition is obtained with a fixed density to self-consistently mimic the preshock condition.Then, the shock condition is exerted onto the preshock gas and the code begins to calculate both the dynamic and chemical changes of the cloud.
The parameters we choose to vary are: the preshock density of hydrogen n H , the CR ionization rate per H 2 molecule ζ, and the shock velocity V s .We fix the magnetic parameter, β = 2, which is defined as B 0 /µG = β n H /cm −3 0.5 , where B 0 is the strength of the initial traverse magnetic field.We ignore the radiation field because of the high column density found by Herschel in the northeastern MCs (see Figure 3).The the propagation time of the cloud shock is fixed at 10 kyr, which is roughly consistent with the age of W28 (Velázquez et al. 2002).When V s ≥ 40 km s −1 , the shock wave is J-type, which would strongly dissociate the molecules on its passage, while the shock wave is Cor CJ-type when V s < 40 km s −1 , which is consistent with the type given by Gusdorf et al. ( 2012) towards W28F.We therefore limit our discussions to V s ≤ 40 km s −1 .The column densities of HCO + and CO are obtained from integrating the number densities along the shock layer.
We present the results of the MHD simulation in Figure 11.
Generally, the abundance ratio N (HCO + )/N (CO) is higher for lower preshock density, larger CR ionization rate and smaller shock velocity.The model with n H = 2 × 10 5 cm −3 , ζ = 2.5 × 10 −15 s −1 and V s = 15-20 km s −1 can explain the observed abundance ratio N (HCO + )/N (CO), which is consistent with the n H derived by Nicholas et al. (2012), ζ derived by Vaupré et al. (2014) and the shock velocity derived by Gusdorf et al. (2012).However, other combinations of parameters can also fit the observed abundance ratio N (HCO + )/N (CO) with low CR ionization rate, but require lower density than that observed.
There are a few caveats in our comparison between the MHD simulation and the observations.First, the 1D simulation cannot reflect the 3D structure of the SNR.Second, the shock does not evolve to steady state except for high velocity shock with high density, meaning that the results of the simulation are sensitive to the propagation time of the cloud shock which is hard to determine.Third, the LTE assumption that the excitation temperature of both of the CO and HCO + emission is equal to the kinetic temperature is not necessarily satisfied in the shocked MCs, and the optically thin assumption may lead to an underestimation of the molecular column densities.Despite these caveats, our simulation shows that the observed N (HCO + ) can be explained by the parameters estimated in previous studies, and that the chemistry induced by the shock and CRs can explain the enhancement of N (HCO + )/N (CO).

The N (HCO
The abundance ratio N (HCO + )/N (N 2 H + ) has also been used to estimate the CR ionization rate (Morales Ortiz et al. 2014;Ceccarelli et al. 2014).Both molecular species are expected to be sensitive to CR ionization rate.The main formation route of N 2 H + is: which is similar to that of HCO + (see Reaction 5).The abundance ratio N (HCO + )/N (N 2 H + ) can be used in chemical simulation to estimate the CR ionization (e.g.Ceccarelli et al. 2014).
Enhanced CR ionization rate is detected in points N5, N6 and N7 close to the clump G6.796 (see Figure 5).These points are located outside the boundary of the W28 SNR, so they have not been disturbed by the shock wave and make it possible to investigate the pure chemical effect of CRs.Towards N6, we derive the column density of HCO + and N 2 H + using the LTE assumption and the temperature derived by Vaupré et al. (2014) who used non-LTE analysis.Since H 13 CO + shows only marginal detection and hence the estimated N (H 13 CO + ) may be unreliable, we derive two values of N (HCO + ): (1) N (HCO + ) = 2.7 × 10 12 cm −2 for the optically thin case, and (2) N (HCO + ) = 50 × N (H 13 CO + ) = 1.5 × 10 13 cm −2 for the optically thick case.For N 2 H + , we assume that it is optically thin because the integrated intensity ratio of the blue and middle hyperfine structures (centered at 13 and 21 km s −1 respectively) is 0.16:1, which is close to the optical-thin value 0.2:1 (e.g.Purcell et al. 2009), and get N (N 2 H + ) = 4.5 × 10 12 cm −2 .The abundance ratio N (HCO + )/N (N 2 H + ) is then 0.6-3.3.This value is within the range found in massive molecular clumps (Hoq et al. 2013) but lower than that found in typical dark clouds like TMC 1 and L134N (Agúndez & Wakelam 2013;Fuente et al. 2019).
To further investigate the relation between the abundance ratio N (HCO + )/N (N 2 H + ) and the CR ionization rate, we use the public version of Nautilus-1.1 chemical code (Ruaud et al. 2016) to model the chemical evolution of a MC.Nautilus is a three-phase gas-grain model which can compute the gas and ice abundances of molecules under various ISM conditions.We adopt the gas-phase reactions from kida.uva.2014(Wakelam et al. 2015).The grain reactions are taken from the original Nautilus code package.
Such discrepancy is expected to result from either the simulation or the observation.The chemical models is not always perfect because of the incomplete chemical network, uncertain reaction coefficients, undiscovered chemical processes, and unknown initial conditions, all of which are hard to deal with.From the observational point of view, we note that the density is derived using 13 CO J = 1-0, J = 2-1 and C 18 O J = 1-0 and J = 2-1 lines (Vaupré et al. 2014).Using the coefficients taken from Leiden Atomic and Molecular Database 7 (van der Tak et al. 2020), we find that the critical density of the used 13 CO and C 18 O lines are below ∼ 10 4 cm −3 , while the critical density of HCO + and N 2 H + 1-0 are both ∼ 10 5 cm −3 .This means that HCO + and N 2 H + trace the gas with higher density than the gas traced by CO isotopes, which will result in the inconsistency of the observed gas density and the density in the simulation.Future observations of higher-J CO transitions, whose critical densities are higher, may shed light on this problem.

Implications on the origin of HESS J1800−240
The origin of the γ-ray emission towards HESS J1800−240 is still under debate.Some argue that all the three sources (A, B, and C) are originated from the escaped CRs from the W28 SNR (Li & Chen 2010;Cui et al. 2018), while others emphasize contribution of local CR accelerators including the massive stellar cluster toward region B (Gusdorf et al. 2015;Hampton et al. 2016) and SNR candidate G5.7−0.1 towards region C (Joubert et al. 2016).If local CR contributors do exist, they are expected to drive fast shock waves that can accelerate particles to > TeV energies.From the view of MC, we do not find either line broadening induced by shock disturbance (see Figure 4) nor significant shock heating (see Figure 10) except towards UC-HII region G5.89−0.39.However, we cannot rule out the possibility that the local accelerator is too young to have exerted significant and observable influence (line broadening and heating effect) on the MCs (Sano & Fukui 2021).Future high-resolution and high-sensitivity observations of multiple molecular transitions are needed to search for local CR-accelerating shock waves interacting with MCs.

X-ray emission towards W28F
The morphology of X-ray emission associated with W28 features an ear-like structure towards the northeastern MC (Rho & Borkowski 2002).In Figure 13, we display a zoom-in view of the HCO + emission towards W28F overlaid with green contours of the X-ray ear.The X-ray emission is located on the western side of the HCO + emission.The close spatial coincidence between the molecular and X-ray emission suggests that the SNR blast wave which has heated the gas and generated an X-ray emitting plasma may be propagating into a dense molecular clump, inducing the broad HCO + line (see the subfigure in Figure 13) and the 1720 MHz OH masers.Assuming a crude pressure balance between the cloud shock and the X-ray emitting gas, we have (McKee & Cowie 1975) 1.4n m m H v 2 m ∼ 2.3n H kT , where n m is the number density of the hydrogen atoms ahead of the cloud shock, n H and T are the hydrogen density and temperature of the X-ray emitting gas respectively (n H ≈ 2.7 cm −3 and kT ≈ 0.3 keV according to Zhou et al. (2014) where k is the Boltzmann  2012)).Then we obtain n m ∼ 340 n H /2.7 cm −3 (kT /0.3 keV) v m /20 km s −1 −2 cm −3 .This typical density is lower than the proposed preshock density (∼ 10 4 cm −3 ) given by Gusdorf et al. (2012) and the critical density of the HCO + 1-0 line (∼ 10 5 cm −3 ).
The production of 1720 MHz OH masers in SNR-MC interaction requires high CR ionization rate (∼ 10 −15 s −1 (Nesterenok 2022)).However, X-ray is also possible to induce enhanced ionization rate in MCs, which facilitates the formation of OH masers (Wardle 1999).The luminosity of the X-ray ear is L ≈ 3.6 × 10 34 erg s −1 (Zhou et al. 2014).The X-ray ionization rate in the molecular clump adjacent to the X-ray ear is approximately (Maloney et al. 1996): Assuming the angular distance between the X-ray ear and the OH masers is 1× the beam size of PMOD, i.e. 1 ′ , and the attenuating hydrogen column density is ∼ 10 22 cm −2 according to the Herschel column density map (see Figure 3), we obtain an X-ray ionization rate of ∼ 10 −16 s −1 , which is lower than the CR ionization rate by an order of magnitude.We note that this is just a rough estimation in order of magnitude, because ζ X depends strongly on the X-ray spectrum (Zhou et al. 2018), while Equation 7 applies to AGNs with harder spectra than SNRs.But indeed, our estimation shows that X-ray ionization does not take a dominant role in the formation of the 1720 MHz OH masers in W28F.

CONCLUSION
In this paper, we conduct HCO + , HCN, and HNC J = 1-0 observations towards the MCs related to SNR W28 using the PMOD 13.7m telescope.With a combination of the archival data of MWISP, SEDIGISM and MALT90, we investigate the spatial distribution of the molecule emissions, their spectra, line ratios, and the chemical processes concerning their abundance ratios.Our main findings are summarized below: 1.In the northeastern MCs, strong emissions of HCO + , HCN and HNC are found spatially coincident with the radio continuum.We find broadened molecular lines with FWHM > 20 km s −1 , which are induced by the shock-cloud interaction.In the MCs south to the SNR, we find clumpy emissions of HCO + and HCN.Spectra of Triple Core Cen (see Figures 1 and 8) show line wing structures indicating disturbance from star formation processes.
2. Enhancement of line ratios I(HCN)/I(HNC) and I( 13 CO J = 2-1)/I( 13 CO J = 1-0) in the shocked regions suggests shock heating.The kinetic temperature estimated with I(HCN)/I(HNC) is roughly consistent with the temperature estimated with NH 3 observations (Maxted et al. 2016), supporting it as a thermometer of molecular gas.
3. We find that the abundance ratio N (HCO + )/N (CO) is typically ∼ 10 −4 in the shocked clouds and ∼ 10 −5 in unshocked clouds.The Paris-Durham MHD shock code is applied to investigate the chemical effects of shock and CR ionization.We find that the abundance ratio N (HCO + )/N (CO) can be reproduced with the physical parameters previously obtained from observations: preshock density n H = 2 × 10 5 cm −3 , CR ionization rate ζ = 2.5 × 10 −15 s −1 and shock velocity V s = 15-20 km s −1 .Therefore, the enhancement of N (HCO + )/N (CO) is a joint effect of shock and CR chemistry.
4. Towards point N6 (see Figure 5) with known high CR ionization rate outside the radio boundary of the remnant, we estimate the abundance ratio N (HCO + )/N (N 2 H + ) ≈ 0.6-3.3.We use the public version of Nautilus chemical code to investigate the relation between N (HCO + )/N (N 2 H + ) and CR ionization rate.The results of the simulation show that the observed value of N (HCO + )/N (N 2 H + ) can be explained if higher density (n H ∼ 2 × 10 5 cm −3 ) is adopted compared to the density derived with multiple transitions of CO isotopes, probably because HCO + and N 2 H + trace denser gas with lower CR ionization rate compared with CO.
Further high-sensitivity and high-resolution observations of multiple molecular species towards regions which exhibit shock perturbation and/or CR bombardment can help us reveal the chemical effects of shocks and CRs, and test the usability and accuracy of current chemical codes.Generally, we get satisfactory results.Although we use four Gaussian components to account for the unshocked emissions in −10 to +30 km s −1 , each of the components has a corresponding 13 CO component, indicating that our estimation of the numbers of components is reasonable.We find a shocked component in 12 CO and HCO + lines without 13 CO counterpart in each region.
In Figure 15, we show the results of the spectral fitting of the C 18 O and H 13 CO + 1-0 lines in core 4, core 4a, triple core Cen, triple core NW, core 5 SW and core 5 NE in the southern MCs of W28.We only fit the spectra of H 13 CO + and C 18 O because the corresponding HCO + and 12 CO emissions are highly optically thick.In core 4a, we find two components in both of the spectra, but we use only the stronger one to calculate column density because the H 13 CO + emission is rather weak.In core 5 NE, we find two components in the C 18 O spectrum, but only one in the H 13 CO + spectrum.

Figure 1 .
Figure 1.Integrated intensity of 12 CO J = 1-0 from −50 km s −1 to +70 km s −1 .The brown contours show the VLA 327 MHz continuum emission (levels are 60 and 180 mJy beam −1 ).White contours trace the TeV emission as seen by HESS (levels are 4-6σ), with HESS J1801−233 in the northeast and HESS J1800−240 in the south (A, B and C of HESS J1800−240 are arranged from east to west).The blue boxes delineate the regions where we observed with PMOD the HCO + and HCN J = 1-0 lines, while the dashed cyan box delineates the region where we observed the HNC J = 1-0 line.The black boxes mark the regions from which we extract the spectra (see Figures 4 and 8)

Figure 2 .
Figure2.Integrated intensity map of HCO + J = 1-0 from −30 km s −1 to +40 km s −1 in the NE region (as delineated in Figure1).The brown contours show the VLA 327 MHz continuum emission (levels are 60 and 180 mJy beam −1 ).The dashed blue contours trace the TeV emission of HESS J1801−233 (levels are 4σ and 5σ).The red crosses are 1720 MHz OH masers detected byClaussen et al. (1997).The black boxes are regions where the HCO + spectra were extracted.
Figure4.Spectra of HCO + J = 1-0 (orange), HCN J = 1-0 (black), HNC J = 1-0 (cyan), 12 CO J = 1-0 (grey), 13 CO J = 1-0 (magenta) and 13 CO J = 2-1 (green) lines averaged in the regions marked in Figure2.Each region takes up two vertical figures to fully display the spectra.The main beam temperature T mb of HCO + , HCN and HNC is multiplied by a factor of 10 for better display.The dashed vertical lines in the figures of Regions 1 to 5 mark a velocity of −3 km s −1 , while the dotted vertical lines in the figures of Regions 3 to 5 mark a velocity of +6 and +12 km s −1 (see the text for details).

Figure 5
Figure 5 also shows the positions of three points, N5, N6 and N7, from whereVaupré et al. (2014) obtained high CR ionization rate (≳ 10 −15 s −1 ).Among the three points, only N6 shows prominent HCO + and N 2 H + line emissions.The spectra of the HCO + , H 13 CO + and N 2 H + lines averaged in a 27 ′′ × 27 ′′ region towards N6 and the results of Gaussian fitting are shown in Figure6.We report only a marginal detection of H 13 CO + , which is inconsistent with the result ofVaupré et al. (2014) who obtained a peak temperature T peak of 0.53 K for H 13 CO + .This may result from the different beam sizes of the MALT90 data (38 ′′ ) and the IRAM 30m data ofVaupré et al. (2014) (29 ′′ ), as well as the relatively low sensitivity of the MALT90 data.The T peak of H 13 CO + in the MALT90 data is ≈ 0.33 K.If the inconsistency of is due to the beam dilution effect, we expect that the T peak detected by IRAM 30m should be approximately 0.33 × (38/29) 2 = 0.57 K, which is roughly consistent with the result ofVaupré et al. (2014).The spectrum of the N 2 H + emission shows clear hyperfine structures.

Figure 13
Figure13.A zoom-in view of the integrated intensity of HCO + from −30 km s −1 to +40 km s −1 overlaid with green contours of XMM-Newton 0.3-7.0keV X-ray flux map taken fromZhou et al. (2014).The red crosses are 1720 MHz OH masers detected byClaussen et al. (1997).An example of the HCO + spectrum towards the maser points is shown in the upper left subfigure.
The authors thank the staff of Qinghai Radio Observing Station at Delingha for their help during the observation.T.-Y.T. thanks Qian-Cheng Liu for instructions on data reduction and analysis, Gao-Yuan Zhang and Benjamin Godard for helpful discussions, and Gavin Rowell for the HESS image of W28.Y.C. acknowledges the support from NSFC grants Nos.12173018 and 12121003.P.Z.acknowledges the support from NSFC grant No. 12273010.S. S.-H.acknowleges support from NSERC through the Canada Research Chairs & the Discovery Grants program.
use of the data from the Milky Way Imaging Scroll Painting (MWISP) project, which is a multi-line survey in 12CO/13CO/C18O along the northern galactic plane with PMO-13.7mtelescope.We are grateful to all the members of the MWISP working group for their support.This research also made use of Montage, funded by the National Science Foundation under Grant NumberACI-1440620, previously  funded by the National Aeronautics and Space Administration's Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology.The SEDIGISM data products are available from the SEDIGISM survey database located at https://sedigism.mpifr-bonn.mpg.de/index.html,which was constructed by James Urquhart and hosted by the Max Planck Institute for Radio Astronomy.The data was acquired with the Atacama Pathfinder Experiment (APEX) under programmes 092.F-9315 and 193.C-0584.APEX is a collaboration among the Max-Planck-Institut fur Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory.

Figure 14 .Figure 15 .
Figure14.Result of the spectral fitting for Region 3 (left panel) and Region N (right panel) in the northeastern molecular clouds of W28.In the upper panel, we show the observed spectra of 12 CO, 13 CO, C 18 O and HCO + 1-0 lines in grey, violet, skyblue and orange respectively.The spectra of C 18 O and HCO + are multiplied by a factor of 5 for better comparison.The result of multi-Gaussian fitting is presented in darker colors and dashed lines for each molecules.In the lower panel, we show the Gaussian components of each spectrum.The thicker lines of 12 CO and HCO + the broad lines induced by the SNR shock wave.

Table 1 .
Results from multi-Gaussian fitting of the emission lines of 12 CO, C 18 O, HCO + and H 13 CO + towards two regions in the northeastern molecular clouds and six regions in the southern molecular clouds.Molecule v0(km s −1 ) T peak (K) FWHM(km s −1 ) N (cm −2 ) N (HCO + )/N (CO) b