© ESO, 2014

1. Introduction

Cosmic rays (CRs) are energetic charged particles that reach the Earth as an isotropic flux. They pervade the Galaxy and play a crucial role in the evolution of the interstellar medium, because they dominate the ionisation of molecular clouds where the gas is shielded from the UV radiation field. The ionisation degree in the molecular gas is a fundamental parameter throughout the star and planet forming process. First, ions couple the gas to the magnetic field and, therefore, they regulate the gravitational collapse of the cloud. In addition, ions sustain turbulence within protoplanetary discs and introduce non-ideal magnetohydrodynamics effects, which influence the accretion rate onto the protostar (Balbus & Hawley 1998; Lesur et al. 2014). Finally, the CR induced ions initiate efficient chemical reactions in the cold molecular clouds, which eventually lead to the formation of complex molecules, which enrich the gas even up to the first stages of planet formation.

Fig. 1 Left: the W28 complex on large scales. Grayscale (in σ) and thick contours show TeV emission as seen by HESS (levels are 4−6σ). Red contours show the CO(1−0) emission (Dame et al. 2001) integrated over 15−25 km s-1 and magenta contours trace the emission integrated over 5−15 km s-1 (levels are 40−70 K km s-1 by 5 K km s-1). Crosses show the positions observed with the IRAM 30m telescope and discussed in this paper. The blue contours show the 20 cm free-free emission in the M20 region (Yusef-Zadeh et al. 2000). The blue circle gives the approximate radio boundary of the SNR W28 (Brogan et al. 2006). Right: the northern cloud in the W28 complex (zoom on the black box). The red contours show the CO(3−2) emission in K km s-1, integrated over 15−25 km s-1 (levels are 15−130 K km s-1 by 5 K km s-1) (Lefloch et al. 2008). Diamonds show the locations of OH masers in the region (Claussen et al. 1997).

Open with DEXTER

However, in order to fully understand the influence of CRs on the above processes across the Galaxy, it is necessary to know where CRs are accelerated and how they propagate through the gas. Unfortunately, since CRs are scattered by magnetic fields all along their path through the Galaxy, the production sites of CRs cannot be observed directly, and the diffusion of CRs also makes the evolution of the energy spectrum during their propagation difficult to determine observationally.

However, we can detect indirect signatures of the interaction of hadronic CRs (essentially protons) with matter. Protons above a kinetic energy threshold of ≈280 MeV produce π⁰ pions when they collide with particles in the molecular cloud. Each pion then decays into two γ-ray photons (π⁰ → 2γ), each with a typical energy that thoseis ~10% that of the colliding proton. Bright γ-ray sources thus indicate regions with a large density of protons with energies above 0.28 GeV. In these regions, the observed γ-ray photon spectrum can, in addition, be used to derive the spectrum of the parental CR particles, before the scattering and propagation within the Galaxy (Ackermann et al. 2013).

Supernova remnants (SNR) are thought to be the sources of CRs. In this scenario, protons are accelerated in the expanding shell of the SNR, following the diffusive shock acceleration process (Bell 1978). Supporting this scenario, there is now clear evidence that SNR are spatially associated with GeV to TeV sources (Aharonian 2013). Moreover, several SNR are close to the molecular cloud that gave birth to the SN precursor. These molecular clouds now act as reservoirs of target material for the freshly accelerated protons, thus enhancing the production rate of γ rays. Although it is compatible with γ-ray observations, the hadronic scenario is challenged by the leptonic scenario involving electron CRs. In this alternative scenario, the γ-ray emission can be explained mainly by inverse Compton scattering of the cosmic microwave background (e.g. Morlino et al. 2009; Abdo et al. 2011). Yet, this scenario cannot explain the spatial correlation of TeV emission with molecular clouds. Moreover, recent observations of the IC443 and W44 SNR with the Fermi-LAT telescope (Ackermann et al. 2013) specifically support a hadronic origin of γ rays, consistent with the so-called SNR paradigm for the origin of primary CR (see e.g. Hillas 2005,for a review).

Cosmic ray protons with kinetic energy below the ≈280 MeV threshold of π⁰ production cannot be traced by the emission of γ-rays. Nevertheless, recent calculations suggest that the ionisation of UV-shielded gas is mostly due to keV-GeV protons (Padovani et al. 2009). Accordingly, low-energy CR protons can be traced indirectly by measuring the ionisation fraction of the dense gas. It has thus been proposed that an enhanced electron abundance in molecular clouds located in the vicinity of SNR could be the smoking gun for the presence of freshly accelerated CRs, with energies ≲1 GeV.

This idea was put forward by Ceccarelli et al. (2011, hereafter CC2011), who measured the ionisation fraction x_e = n(e⁻)/n_H in the W51C molecular cloud, located in the vicinity of the W51 SNR. The detection of TeV emission by both HESS and MAGIC telescopes close to the molecular cloud is evidence of a physical interaction with the SNR. This supports the idea of the pion-decay production of γ rays with W51C acting as a γ-ray emitter. In CC2011, an enhanced ionisation fraction was reported towards one position, W51C-E, which required a CR ionisation rate two orders of magnitude larger than the typical value of 1 × 10^-17 s^-1 in molecular clouds. This observational evidence strongly supports the hadronic scenario of γ-ray production, at least for W51.

Complementary studies of the CR ionisation rate in several diffuse clouds close to SNR have been carried out using different techniques, such as absorption (McCall et al. 2003). These studies also show an enhancement of a factor of 10−100 of the CR ionisation rate (Indriolo et al. 2010; Indriolo & McCall 2012) with respect to the canonical value. However, the interpretation is not straightforward, as Padovani et al. (2009) showed that the penetration into the cloud of high energy CRs results into an enhanced CR ionisation rate in low density molecular clouds even in absence of an increased CR flux.

The combined observations of two extreme energy ranges, namely TeV and millimetre, seems a powerful method to characterise an enhanced concentration of proton CRs. It also gives additional evidence supporting a physical interaction of the SNR shock with molecular clouds. From a theoretical point of view, it is expected that the most energetic CR protons diffuse at larger distances ahead of the SNR shock front, whilst the low-energy tail of the distribution remains closer. As a consequence, one expects that any ionisation enhancement by low energy CRs should be localised accordingly. In CC2011, however, only one location could be used to derive the ionisation fraction, and no constraint could be given regarding the spatial distribution of the ionisation and therefore the diffusion properties of CRs.

The aim of this paper is to present measurements of the ionisation fraction within the molecular clouds in the vicinity of the W28 SNR. The paper is organized as follows. In Sect. 2, the W28 association is presented, with particular emphasis on the physical link between the SNR and the molecular clouds. In Sect. 3, the millimetre observations are described. The derivation of the physical conditions is presented in Sect. 4. The derivation of the ionisation fraction and the CR ionisation rates are described in Sects. 5 and 6, where we stress the strengths and limitations of the method. The results are discussed in Sect. 7.

2. The W28 association

The W28 SNR has an age greater than 10⁴ yr, and is likely in the Sedov or radiative phase (Westerhout 1958; Lozinskaya 1974). Its distance is estimated between 1.6 kpc and 4 kpc, based on kinematic determinations and Hα observations (Goudis 1976; Lozinskaya 1981). The LSR velocity, based on Hα and [NII] observations, is estimated to be 18 ± 5 km s^-1 (Lozinskaya 1974). In the remainder of the text, all velocities will refer to the local standard of rest (LSR) and projected distances will be given for distances of both 1.6 and 4 kpc. The boundary diameter of the SNR is 42 arcmin, corresponding to a linear radius of 9.6 to 24 pc (Fig. 1).

The large-scale region towards W28 contains a variety of objects such as HII regions (e.g. M8, M20, W28A2), molecular clouds, SNR, and new SNR candidates (Brogan et al. 2006). Molecular gas, as seen in CO(1−0) (Wootten 1981; Dame et al. 2001), coincides spatially with the W28 SNR, suggesting a physical association, and supporting a view in which the W28 SNR is interacting with its parental molecular cloud. Probably related is the fact that ongoing massive star formation has been observed in these molecular clouds, consistent with the triggered star formation scenario (Elmegreen 1998). The molecular gas located towards the north-east of the SNR boundary was mapped at high spatial resolution, in the CO(3−2) rotational line, by Lefloch et al. (2008) revealing a fragmented filamentary structure elongated north-south (Fig. 1).

In γ rays, the high spatial resolution H.E.S.S. imaging array of Cherenkov telescopes has revealed the presence of extended TeV emission (Aharonian et al. 2008), which splits into two components, separated by 14−34 pc: HESS J1801-233, in the north, and HESS J1800-240 in the south. The latter further splits into three well-separated components (Fig. 1). In projection, the entire TeV emission coincides with the molecular gas seen in CO(1−0) which also appears to bridge the northern and southern TeV components. The molecular gas coinciding with the northern TeV component J1801-233 shows velocities predominantly from 15 km s^-1 to 25 km s^-1, as does the southern J1801-240 A component observed in CS(1−0) by Nicholas et al. (2012). However, the southern components J1801-240 B and C coincide with molecular emission characterised by somewhat lower velocities, from 5 km s^-1 to 15 km s^-1.

Whether the TeV emission is physically associated with the molecular clouds is of utmost importance, for the question of pion-decay production. However, there are several indications, based on kinematic information, that this may well be the case.

First, when inspected in velocity space, the CO(1−0) emission covers the entire range from 5 to 25 km s^-1 continuously (Fukui et al. 2012a). This indicates that the molecular emission, traced either by CS or CO, is physically linked over the entire region, and not only in projection.

Second, OH masers have been reported towards the northern component (Claussen et al. 1997; Hewitt et al. 2008), with velocities ranging from 7.1 km s^-1 to 15.2 km s^-1. Such OH masers are thought to trace the interaction of the SNR shock with the molecular gas, an interpretation which is consistent with the velocity range of the CO(1−0) emission.

Finally, velocity differences observed between the northern and southern clouds are compatible with the differences up to ~6 km s^-1 observed in the M20 map of CO(3−2), encompassing the northern TeV component (Lefloch et al. 2008), indicating that velocity shifts of several km s^-1 are found which could be due to an interaction with the SNR.

Taken all together, these data strongly suggest a 3D picture in which the SNR is interacting with surrounding molecular clouds covering a wide and continuous velocity range, typically from 5 km s^-1 to 25 km s^-1, and characterised by large variations along the line of sight.

It is therefore most likely that the W28 association displays the interaction of the SNR with molecular clouds, and thus is an excellent target in which to study the ionisation by CRs.

3. Observations and data reduction

Observations were carried out over 40 h in December, 2011, and March, 2012, with the IRAM 30 m telescope. We used the EMIR bands with the Fast Fourier Transform Spectrometer as a backend in the position-switching mode, using OFF positions about 1000′′ to the east. Amplitude calibration was done typically every 15 min, and pointing and focus were checked every 1 and 3 h, respectively, ensuring ≈2′′ pointing accuracy. All spectra were reduced using the CLASS package of the GILDAS¹ software (Pety 2005). Residual bandpass effects were subtracted using low-order (≤3) polynomials. The weather was good and T_sys values were lower than 220 K. The observed molecular transitions used in the present work are listed in Table 1, along with the associated system temperature and sensitivity ranges T_sys and σ_rms obtained during the successive runs of observations. All spectra are presented on the main-beam temperature scale, , with F_eff and B_eff the forward and main-beam efficiencies of the telescope, respectively.

We observed towards 16 positions, of which 10 are located in the northern cloud and 6 in the southern cloud. The coordinates of these positions are listed in Table 2. The two lowest rotational transitions of ¹³CO and C¹⁸O are used to derive the physical conditions in the cloud, while H¹³CO⁺(1−0) and DCO⁺(2−1) are used to derive the cosmic ray ionisation rate.

4. Results

4.1. Observed spectra

The resulting spectra towards all positions are shown in Fig. 2. The ¹³CO and C¹⁸O spectra show multiple velocity components, most likely associated with several clouds along the line of sight. In some instances, negative features are apparent, which are due to emission from the reference position at different velocities. Isotopologues ¹³CO and C¹⁸O are detected towards 12 of the 16 positions, and the spectra show clearly two main components. The velocity of the dominant component varies between the northern (≲21 km s^-1) and southern (≳7 km s^-1) cloud, as presented in Sect. 2. At most positions, the rarer C¹⁷O isotopologue is also detected, although the hyperfine structure of the (1−0) transition is not always well resolved. The H¹³CO⁺(1−0) and DCO⁺(2−1) emission lines have only one velocity component, at 21 km s^-1. The C¹⁸O(1−0) and (2−1) lines show a clear distinction between the northern positions, where the line emission is the strongest, and southern positions. The compound H¹³CO⁺ is detected towards all northern positions but N8. However, DCO⁺ is detected only towards three positions, N2, N5, and N6. Southern positions display much weaker emission and SE1 is detected in both H¹³CO⁺ and DCO⁺.

The analysis was essentially driven by the H¹³CO⁺ and DCO⁺ lines which are the main focus of the present work. The velocity of the dominant CO transition always corresponds to the velocity of the H¹³CO⁺ line when detected. When more velocity components are detected in CO, we limited the Gaussian fits to the first two dominant components. Upper limits are given at the 1σ level for T_peak, and for the integrated intensity W by assuming Δv = 3 km s^-1. The results from the Gaussian fits are summarised in Table 3 (CO isotopologues) and Table 4 (H¹³CO⁺, DCO⁺).

In the following, we describe the two-step analysis of the measured line intensities. First, the physical conditions are derived according to the ¹³CO and C¹⁸O lines. Second, the HCO⁺/DCO⁺ abundance ratio is compared with model predictions computed using the derived physical conditions.

Fig. 2 Observations of millimetre emission lines towards all positions. The intensities are in units of main-beam temperature (K). For readibility, a multiplicative factor was applied to the spectra and is given under each transition line label. This factor was decreased for position N2 (protostar) where the signal is very strong. Vertical dashed lines indicate the two extreme velocity components of the complex at 7 and 21 km s-1.

Open with DEXTER

4.2. Determination of physical conditions

The physical conditions prevailing at the various locations were determined based on the ¹³CO (1−0) and (2−1) lines and the C¹⁸O (1−0) and (2−1) lines. To accomplish this, we performed non-local thermal equilibrium (non-LTE) calculations of the rotational level populations under the large velocity gradient (LVG) approximation (Ceccarelli et al. 2003). The H₂ density, gas kinetic temperature T_kin, and total column density of each species covered a large parameter space. For each set of input parameters, the expected line intensities and integrated intensities were computed adopting the linewidth Δv derived from the Gaussian fits (Table 3), and taking into account beam dilution effects by varying the size of the emitting regions. A simple χ² minimization was then used to constrain the physical conditions that best reproduce the observed intensities. The LVG analysis was performed towards positions where all four lines were detected. We considered collisions of ¹³CO and C¹⁸O with both para and ortho H₂, assuming an ortho-to-para ratio in Boltzmann equilibrium. Hence, in the temperature range considered here, H₂ is mainly in its para configuration. We used the collisional cross sections from Yang et al. (2010).

In this process, the relative abundances of C¹⁸O to ¹³CO need to be known, and we have assumed that the molecular isotopic ratios are equal to the elemental isotopic ratios, namely that (1)It is known that isotopic ratios may vary with the position within the Milky Way, resulting from the gradual depletion of ¹²C and enrichment of ¹³C with the cycling of gas through stars (e.g. Wilson & Rood 1994; Frerking et al. 1982). In addition, local variations are also possible, resulting from the competition of chemical fractionation and selective photodissociation (van Dishoeck & Black 1988; Federman et al. 2003). The dependence of the ¹²C/¹³C isotopic ratio on galactocentric distance was studied by Milam et al. (2005). Applying their results to W28, which is 4−6 kpc from the Galactic center, one gets ¹²C /¹³C = 50 ± 7. In practice, in our non-LTE analysis, we varied the ¹³C /¹²C ratio. The best χ² values were obtained using an isotopic ratio of 50, consistent with the above expectation, which we therefore adopted in what follows. Regarding the ¹⁶O /¹⁸O isotopic ratio, we adopted the value of 500 representative of the solar neighborhood.

The results of the LVG analysis are summarised in Table 5. At each position, the size of the emitting regions was found to be larger than the beam size. Position N2, which is considered a protostar in Lefloch et al. (2008), shows a peculiar behaviour with a visual extinction at least 5 times higher than for the other positions. Overall, the density and kinetic temperature we derived are typical of dense molecular clouds with visual extinctions larger than 10 mag. We also note that the H₂ densities we derived in the northern cloud are consistent with values published by Lefloch et al. (2008).

4.3. The DCO⁺/HCO⁺ abundance ratio

The species DCO⁺ was detected towards four positions, for which it was possible to derive values of the DCO⁺/ HCO⁺ abundance ratio. For all other positions where H¹³CO⁺ was detected, upper limits at the 1σ level on the abundance ratio we derived. We determined the column densities of H¹³CO⁺ and DCO⁺ from the same non-LTE LVG calculation, using the collisional cross sections from Flower (1999) and the physical conditions derived from the CO observations (Table 5). We derived the observed DCO⁺/ HCO⁺ abundance ratio for each set of physical conditions (n_H2,T), assuming ¹²C /¹³C = 50 (see Sect. 4.2). The uncertainty in the abundance ratio is dominated by the uncertainties in the physical conditions. Results are listed in Table 5 and will be used in the next section to constrain the CR ionisation rate.

5. Methods to measure the CR ionisation rate ζ

5.1. Analytical method

In their seminal paper, Guélin et al. (1977, hereafter G77) suggested that the abundance ratio of DCO⁺ to HCO⁺, which we denote R_D = DCO⁺/ HCO⁺, can be used to measure the ionisation fraction in molecular clouds, x_e = n(e⁻) /n_H. Subsequently, Caselli et al. (1998, hereafter C98) proposed using the R_H = HCO⁺/ CO abundance ratio in combination with R_D to derive both x_e and ζ in dark clouds. The basic idea is that DCO⁺ and HCO⁺ are chemically linked, and depend on a reduced number of chemical reactions where the CR ionisation rate plays a crucial role, through the ionisation of H₂ into H (Herbst & Klemperer 1973), leading to the formation of the pivotal H ion. The fast ion-neutral reaction of H with HD produces the deuterated ion H₂D⁺ which then initiates the formation of several deuterated species, including DCO⁺. In a similar fashion, HCO⁺ is formed by the reaction of H with CO. This forms the basis of the method of C98 which uses CO, DCO⁺, and HCO⁺ to derive x_e and ζ. The full set of reactions used in the C98 analysis is given in Table A.1 with updated reaction rates.

The steady-state abundance ratios R_H = HCO⁺/ CO and R_D = DCO⁺/ HCO⁺ can be analytically derived from this network, provided that HCO⁺ and DCO⁺ are predominantly formed and destroyed by reactions 1−4 and 7−9, respectively. One then finds that (2)and (3)where n(X) is the number density of species X and x(X) = n(X) /n_H its fractional abundance. In the following, we will assume that the gas is fully molecular such that n_H = 2n(H₂). The coefficients β, β′, and k are the reaction rates listed in Table A.1. Finally, is the total destruction rate of or H₂D⁺ by neutrals. Provided that R_H, R_D, n_H, and the kinetic temperature are known, x_e and ζ can then be derived as

where ΔE = 220 K, such that at sufficiently low temperatures the last term in brackets in Eq. (4) becomes negligible. Equation (4) demonstrates that x_e only depends on the abundance ratio R_D and the gas kinetic temperature, as originally proposed by G77.

Figure 3 shows x_e as a function of R_D, as predicted from Eq. (4), assuming a kinetic temperature of 20 K (blue dashed line). We note two regimes in the dependence of x_e on R_D. For low R_D values (≲10^-2), x_e is proportional to 1 /R_D with a factor that depends on the chemical reaction rates and the HD abundance. For higher R_D values, x_e drops sharply. In this regime, x_e varies by more than two orders of magnitude when R_D is changed by only a factor of two. The slope becomes steeper with increasing temperature, due to the predominance of the reverse reaction of the formation of H₂D⁺ (Table A.1). This indicates that accurate values of x_e in dark clouds through this method require extremely accurate values of R_D. It also suggests that in regions with higher ionisation, where x_e is proportional to 1 /R_D, DCO⁺ will be difficult to detect.

5.2. Numerical models

Fig. 3 xe = n(e−) /nH as a function of RD. Values from the analytical method (G77) in Sect. 5.1 are given at T = 20 K (dashed line). Results from our calculation (see Sect. 5.2) for temperatures between 7 K and 20 K are contained within the solid lines. ζ values are given for nH = 104 cm-3, although the model only depends on the ζ/nH ratio (see text). RD values or upper limits as derived from observations are indicated by the red symbols; see Sect. 6.1.

Open with DEXTER

To assess the validity of the analytical approach of G77 and C98, we have solved the OSU 2009² chemical network for the abundances of DCO⁺, HCO⁺, and e⁻, using the astrochem³ code. The time evolution of the gas-phase abundances was followed until a steady state was reached, after typically 10 Myr. The underlying hypothesis is that the cloud was already molecular when it was irradiated by the CRs emitted at the SN explosion. The chemical changes caused by the sudden CR irradiation are dominated by ion-neutral reactions, whose timescale is only ~10² yr, much shorter than the W28 SNR age (~10⁴ yr; see Introduction). Details of the models are given in Appendix A. In one calculation, the gas is shielded by 20 mag of visual extinction, such that UV photons can safely be ignored.

As anticipated from Eq. (2), we found that the abundances mostly depend on the ζ/n_H ratio, rather than on ζ and n_H separately. This behaviour is similar to photon-dominated regions in which a good parameter is the ratio of the UV radiation field to the total density. For each temperature from 7 K to 20 K, a series of calculations with ζ/n_H increasing from 10^-22 cm³ s^-1 to 10^-18 cm³ s^-1 were performed, and the steady-state values of x_e and R_D were recorded.

In these calculations, we assumed a standard CO abundance of ~7.3 × 10^-5 in the cloud. Indeed, since the density of the studied clouds is relatively low (Table 5), we do not expect CO depletion to play an important role here. However, in general, one has to take into account this uncertainty. In addition, our calculations do not consider separately the H2 ortho-to-para ratio, which is known to affect the DCO⁺/HCO⁺ abundance ratio when it is larger than about 0.1 (e.g. Pagani et al. 2011). Our choice here is based on the published observations that indicate that the H2 ortho-to-para ratio is smaller than about 0.01 in molecular clouds (e.g Troscompt et al. 2009; Dislaire et al. 2012).

5.3. A new view of the DCO⁺/HCO⁺ method

The results of the numerical models are shown in Fig. 3. As expected, the ionisation fraction x_e increases with ζ/n_H, with R_D decreasing in the process. There is good overall agreement between the analytical and numerical predictions for R_D ≳ 2 × 10^-2. In this high-R_D regime, the small differences between analytical and numerical values are due to the abundances of HD and CO not being constant as originally assumed by G77 and C98. Instead, as ζ/n_H and x_e increase, atomic deuterium becomes more abundant. Similarly, the CO abundance decreases because of the dissociating action of CRs. When R_D decreases and reaches ≈2 × 10^-2, the abundances predicted by the numerical model change dramatically to a regime characterised by large values of x_e and low values of R_D. This jump corresponds to the well-known transition from the so-called low ionisation phase (LIP) to the high ionisation phase (HIP) (Pineau des Forêts et al. 1992; Le Bourlot et al. 1993), and is due to the sensitivity of interstellar chemical networks to ionisation. The LIP is associated with R_D larger than 10^-2, whilst the HIP is characterised by R_D ≲ 10^-4. In our calculations, the LIP-HIP transition occurs at ζ/n_H ~ 3 × 10^-19 cm³ s^-1, and we note that this value depends only slightly on the temperature, although it is known to depend on other parameters such as the gas-phase abundance of metals (Wakelam et al. 2006a). A detailed analysis of the LIP-HIP transition is, however, not the aim of this study. Here, it is rather the existence of this instability which is of interest since it produces a sharp difference between the analytical and the numerical predictions from the ionisation point of view. Application to a practical case shows that what changes is not the jump itself but (slightly) the ζ at which it occurs (see e.g. Ceccarelli et al. 2011). In the former, the variations of x_e and R_D are continuous and, as already mentioned, predict x_e ~ 1 /R_D in the low-R_D regime. This scaling is, however, not observed in the numerical models, and is replaced by a discontinuous variation of both x_e and R_D. The present calculations show that the LIP is characterised by R_D = 10^-2− 10^-1, x_e ≲ 5 × 10^-7, and the HIP is characterised by R_D ≈ few 10^-5 and x_e ≈ few 10^-5.

As shown in Fig. 4, the low values of R_D in the HIP are due to a very low abundance of DCO⁺, whilst the abundance of HCO⁺ decreases by a smaller amount. This has important consequences when using the DCO⁺/HCO⁺ method to derive the ionisation fraction and CR ionisation rate. First, it must be recognised that this method may provide a value of x_e only for LIP-dominated gas conditions. In other words, where DCO⁺ is detected, the line of sight is dominated by low-x_e gas. For lines of sight dominated by HIP gas, the abundance of DCO⁺ is expected to be well below detectable thresholds, such that only upper limits on R_D can be derived. Yet, an upper limit on R_D still provides essential information, since it is associated with a lower limit on x_e, which in turn corresponds to a lower limit of ζ/n_H. On the contrary, for LIP-dominated lines of sight, the value of x_e and ζ/n_H may be derived directly from R_D, although x_e is extremely sensitive to uncertainties on R_D in this regime.

6. The CR ionisation rate in W28

6.1. Constant density and temperature cloud analysis

A new view of the DCO⁺/HCO⁺ method thus emerges, which stresses its strengths and limitations. The method allows the determination of the ionisation fraction x_e and the ζ/n_H ratio for gas in the LIP, and provides lower limits of x_e and ζ/n_H for gas in the HIP. In the following, we apply this method to the sample of observed points, using the constraints on the gas temperature and density, and the R_D value in each point (Table 5). We emphasise that the model calculations summarised in Fig. 3 assumed constant density and gas temperature. In the next section, we will discuss how the DCO⁺/HCO⁺ method can be used to constrain x_e and ζ/n_H, taking into account the thermal structure of the cloud.

Of the 16 lines of sight initially observed in CO, 12 were also detected in H¹³CO⁺, of which 4 led to R_D determinations and 5 to upper limits (Table 5). The four points with measured R_D are N5, N6, SE1, and N2. In the following analysis, we exclude N2 as it coincides with a protostar, which means that a more accurate analysis taking into account the structure of the protostar and the inner ionisation is necessary. The values obtained towards N5, N6, and SE1 are shown in Fig. 3.

The SE1 point lies on the LIP branch, enabling a determination of the ionisation fraction x_e = (0.15−4) × 10^-7 and of the CR ionisation rate ζ = (0.2−20) × 10^-17 s^-1. On the contrary, the values of R_D towards N5 and N6 lie in the gap between the LIP and HIP branches, even when considering a kinetic temperature as high as 20 K, a temperature which is larger than the values derived for these positions. In these cases, adopting n_H = 2n(H₂) ≳ 4 × 10³ cm^-3 (Table 5), Fig. 3 provides the following lower limits: x_e ≳ 4 × 10^-7 and ζ ≳ 1.3 × 10^-15 s^-1  for both points. We note that the detection of DCO⁺ indicates that the line of sight includes a non-negligible amount of LIP, which can serve to further constrain the value of ζ/n_H. This can be seen in Fig. 4, which shows R_D and DCO⁺ as a function of ζ, for a density and temperature appropriate to position N5 (Table 5). The measured R_D intersects the model predictions at the edge of the LIP/HIP jump, at ζ ≈ 2.5 × 10^-15 s^-1. More importantly, the figure shows that the gas is neither entirely in the LIP nor HIP state as expected from the detection of DCO⁺. A similar plot has also been obtained for N6, leading to the same conclusion.

	Fig. 4 RD = DCO+/ HCO+ (thick line, left axis) and x(DCO+) = n(DCO+) /nH (dashed line, right axis) as a function of ζ, for T = 10 K and nH = 8 103 cm-3, i.e. the physical conditions characterising position N5. The HIP and LIP are marked. The hatched area shows the range of observed RD at that position.
Open with DEXTER

Finally, the non-detection of DCO⁺ in the other lines of sight leads to upper limits on R_D that are well outside the LIP branch. At these positions, the gas is very likely to be almost entirely in the HIP state, which means that ζ/n_H ≳ 3 × 10^-19 cm³ s^-1. An exception is the point SW4, where an extremely energetic outflow has been detected (Harvey & Forveille 1988), and where the HCO⁺ is therefore likely contaminated by the outflowing material.

6.2. Constant density cloud analysis

As discussed in the previous section, the points N5 and N6 are likely composed of a mixture of gas in the LIP and HIP state. This is similar to the situation observed in W51C-E (CC2011). In that case, ζ was estimated with a model that takes into account the thermal and chemical structure of a constant density cloud, where a fraction (the deepest) is in the LIP and the rest in the HIP (Fig. 2 in CC2011). Here, we do a similar analysis, using basic arguments instead of a sophisticated model, and we show that it leads to similar results, namely the determination of ζ to within a factor of 2. The advantage of this analysis is that it shows in a straightforward way the uncertainty due to the model parameters.

The crucial point is understanding what causes the gas to flip from the HIP to the LIP state going deeper into the cloud. Since the column density is too low to appreciably reduce ζ across the cloud, the only macroscopic quantity that changes is the gas temperature. Specifically, the temperature increases by a few K (in the UV-shielded region) going deeper into the cloud because the CO line opacity increases and, consequently, the line cooling becomes less efficient. The effect is larger for larger ζ as the heating, dominated by the CR ionisation, is less compensated by the line cooling.

It is instructive to see how the R_D ratio changes as a function of the gas temperature for different ζ. This is shown in Fig. 5, for a range of temperatures (5−80 K) and ζ/n_H (2–5 × 10^-19 cm³  s^-1, appropriate to the N5 and N6 points). In these calculations, we consider a cell of gas of constant density, shielded by 20 mag of visual extinctions as before, such that the ionisation is driven by the CRs. The figure shows important features:

• i)

  for ζ/n_H ≲ 2 × 10^-19 cm³ s^-1, the cloud is always in the LIP, regardless of the temperature;

• ii)

  for ζ/n_H ≳ 5 × 10^-19 cm³ s^-1, the cloud is always in the HIP for temperatures lower than 50 K;

• iii)

  for intermediate values of ζ/n_H, the gas flips from HIP to LIP with increasing temperature, and the larger ζ is, the larger the temperature where the flip occurs.

These calculations show that there is a range of ionisation rates in which the gas is extremely sensitive to temperature variations. A tiny increase in temperature is sufficient to make the gas flip from the HIP to the LIP. In particular, for position N5, there is such a combination of values of R_D, T_kin, and ζ/n_H (see Fig. 5) that places it precisely in a region where the transition from HIP to LIP can be triggered by an increase in T_kin as small as a few K. A similar argument applies to N6. In addition, in regions exposed to an enhanced CR ionisation rate, the outer part characterised by large ionisation fractions will be extended farther into the cloud, thus decreasing the relative amount of LIP with respect to the HIP.

Fig. 5 RD as a function of the gas temperature Tkin for different values of ζ/nH: from 2 to 5 × 10-19 s-1, as marked. Note that for ζ/nH ≤ 2 × 10-19 cm3 s-1 (thick solid line), the cloud is always in the LIP, regardless of the temperature. For ζ/nH> 5 × 10-19 cm3 s-1 (thin dashed curve), the cloud is always in the HIP for temperatures ≤50 K. Hatched areas show observations of N5 and N6. We assume AV = 20 mag.

Open with DEXTER

Based on the derived kinetic temperatures and values of R_D (Table 5) and using Fig. 5, we can further constrain the value of ζ/n_H. Towards N5, the temperature was found to be 10 ± 2 K, while R_D = 0.014−0.020. When inserted into Fig. 5, these delineate a region that is compatible with a narrow range of ζ/n_H = (2.8−3.0) × 10^-19 cm³ s^-1. For the N6 line of sight, we find similar values, (2.9−3.2) × 10^-19 cm³ s^-1. The densities derived from the analysis and their uncertainties then lead to cosmic ray ionisation rates of (1.3−3.3) × 10^-15 s^-1 and (1.3−4.0) × 10^-15 s^-1 for N5 and N6 respectively. Results are summarised in Table 5 and are included in Figs. 6 and 7.

7. Discussion

Fig. 6 Compilation of measured ζ in different objects (open squares), as reported by Padovani & Galli (2013). The black filled square denotes W51 (Ceccarelli et al. 2011). Red points and lower limits report the values derived in this work. The dashed lines show the range of column densities (0.5−10) × 1022 cm-2, typical of dense molecular clouds, corresponding to visual extinctions of 5 and 100 mag, respectively. On the left lie the diffuse clouds and on the right highly obscured environments such as infrared dark clouds or protoplanetary discs.

Open with DEXTER

Table 5 lists the observed positions and the corresponding CR ionisation rates derived using the method described in the previous section. With the exception of the SE1 point, in all other points ζ is at least 10 to 260 times larger than the standard value (1 × 10^-17 s^-1) in Galactic clouds. This is shown in Fig. 6, where we present a compilation of the ζ measured in various objects (from Padovani & Galli 2013), plus our measurements. In the range of column densities (0.5−10) × 10²² cm^-2, typical of dense molecular clouds, the points in which we derived ζ are those with the highest values, together with the CC2011 point (filled square). The first conclusion of this work is, therefore, that clouds next to SNR are indeed irradiated by an enhanced flux of CRs of relatively low energy (see below for a more quantitative statement on the CR particle energies).

	Fig. 7 CR ionisation rate ζ as a function of the approximate projected distance from the SNR radio boundary (blue circle in Fig. 1), assuming a W28 distance of 2 pkc. We note that the ζ error bars are dominated by the uncertainties on the H2 densities (see text).
Open with DEXTER

The dependence of ζ on the projected distance from the SNR radio boundary (assuming a W28 distance of 2 kpc) is shown in Fig. 7. Remarkably, the point farthest (~10 pc) from the SNR edge is the one with the lowest ζ. Actually, it is the only point where the gas is predominantly in the LIP state. All other points, at distances ≲3 pc, have at least a fraction of the gas in the HIP, namely they have a larger x_e and ζ. Of course, this analysis does not take into account the 3D structure of the SNR complex. Yet, this can still provide us with constraints on the propagation properties of CRs, as will be discussed in the following.

Valuable additional information is provided by observations in the γ-ray domain. Both the northern and southern clouds coincide with sources of TeV emission, as seen by HESS. This means that the clouds are illuminated by very high energy (≳10 TeV) CRs, which already escaped the SNR expanding shell and travelled the ≳10 pc (or more, if projection effects play a role) to the southern cloud. Conversely, the low CR ionisation rate measured in SE1 tells us that the ionising lower energy CRs remain confined closer to the SNR. In the same vein, GeV emission has been detected towards the northern region but only towards a part of the southern one. This difference between the GeV and TeV γ-ray morphology has been interpreted as a projection effect: the portion of the southern region that exhibits a lack of GeV emission is probably located at a distance from the shock significantly larger than the projected one, >10 pc, and thus can be reached by ≳TeV CRs but not by ≳GeV ones (Gabici et al. 2010; Li & Chen 2010; Nava & Gabici 2013). Remarkably, the SE1 point is located in the region where the lack of GeV emission is observed.

The picture that emerges is that of a stratified structure with CRs of larger and larger energies occupying larger and larger volumes ahead of the shock. Within this framework, it is possible to estimate the CR diffusion coefficient in the region. This can be done by recalling that in a given time t, CRs diffuse over a distance , where D is the energy dependent diffusion coefficient. For the situation under examination, one gets (6)where D_{(≈ 10 GeV)} is the diffusion coefficient of ≈10 GeV CRs, which are those responsible for the ≈GeV γ-ray emission, and t is the time elapsed since the escape of CRs from the SNR. The value obtained in Eq. (6) is in substantial agreement with more accurate studies (see e.g. Nava & Gabici 2013).

The diffusion coefficient obtained in Eq. (6) can then be rescaled to lower energies, according to D ∝ p^sβ, where p is the particle momentum, β = v/c its velocity in units of the speed of light, and s depends on the spectrum of the ambient magnetic turbulence. The typical value of s in the interstellar medium is poorly constrained to be in the range 0.3 to 0.7 (Castellina & Donato 2011). In the following, we adopt s = 0.5. To estimate the diffusion length of low energy CRs, one has to keep in mind that, while CRs with energies above ≈GeV are virtually free from energy losses (the energy loss time for proton–proton interactions in a density n_H ≈ 10³ cm^-3 is comparable to the age of the SNR), lower energy CRs suffer severe ionisation losses over a short timescale (Berezinskii et al. 1990): (7)This approximate expression is sufficiently accurate in the range of energies spanning 1−100 MeV. The diffusion length of low energy CRs can then be estimated by equating the diffusion time to the energy loss time τ_ion, which gives R_d ≈ 0.02, 0.3, and 3 pc for CRs of energy 1, 10, and 100 MeV, respectively. This implies that only CRs with energies ≳100 MeV can escape the shock and spread over a distance of 3 pc or more, and thus these are the CRs that play a major role in ionising the gas. Whether the ionisation of the gas is due directly to these CRs or to the products of their interaction with the gas (namely slowed down lower energy CRs) remains an open question. It is remarkable that the particle energies of ionising CRs (≈0.1−1 GeV) also make them capable of producing sub–GeV γ rays, given that the kinetic energy threshold for π⁰ production is ≈280 MeV.

Of course, the order of magnitude estimates discussed in this section cannot substitute in any way more sophisticated calculations, yet they clearly indicate an intriguing possible link between low and high energy observations of SNR environments. In fact, in the scenario described above, the very same CRs are responsible for both ionisation of the gas and production of low energy γ rays. If confirmed, such a link would constitute robust evidence for the presence of accelerated protons in the environment of the SNR W28, a thing that would bring further support to the idea that SNR are the sources of Galactic CRs. Additional theoretical investigations are needed in order to examine and possibly rule out alternative scenarios which may include other contributions to the ionisation rate (e.g. CR electrons, X-ray photons) or different means of propagation (e.g. straight-line or advective propagation of low energy CRs).

8. Conclusion

In this work, we presented new observations to measure the CR ionisation rate in molecular clouds close to supernova remnants (SNR). In doing so, the DCO⁺/HCO⁺ method was also revisited. The major results may be summarised as follows.

• 1)

  We observed the two lowest rotational transitions of¹³CO and C¹⁸O towards 16 positions in the northern and southern clouds close to theSNR W28. The four lines were detected in emission towards 12 of these positions, where we could, therefore, derive the physical conditions using a non-LTE LVG analysis. With the exception of one position (N2) coinciding with a proto-star in the region, we derived H₂ densities and temperatures typical of molecular clouds, namely n_H2 = (0.2−10) × 10³ cm^-3 and T = 6−24 K. We searched for H¹³CO⁺ and DCO⁺ line emission in the above 12 positions, and detected it in 9 and 4, respectively. From these data, we could derive the R_D = DCO⁺/ HCO⁺ in 4 positions, one of which coincides with the protostar, and give upper limits for the remaining 5 positions.

• 2)

  We reinvestigated the DCO⁺/HCO⁺ method used to derive the ionisation fraction x_e = n(e⁻) /n_H and the relevant CR ionisation rate ζ causing it. To this aim, we compared the steady-state abundances of HCO⁺, DCO⁺, and e⁻ as predicted by the analytical model of G77, to numerical calculations, assuming constant density and gas temperature. The numerical model leads to two well separated regimes of ionisation, also known as the low- and high-ionisation phases (LIP and HIP; Pineau des Forêts et al. 1992). In the context of this work, these two phases lead to two separated regimes in terms of ζ and R_D values:

  • i)

    for ζ/n_H ≲ 3 × 10^-19 cm³ s^-1, the gas is in the LIP, where R_D ≳ 2 × 10^-2 and x_e ≲ 6 × 10^-7. In this regime, the dependence of x_e on R_D is very steep leading to large uncertainties on x_e;

  • ii)

    for ζ/n_H ≳ 3 × 10^-19 cm³ s^-1, the gas is in the HIP, where R_D ≲ 10^-4 and x_e ≳ 2 × 10^-5. In this regime, DCO⁺ is not detectable and the numerical prediction for x_e differs significantly from the analytical one.

  Therefore, the DCO⁺/HCO⁺ abundance ratiocan provide a measure of x_e and ζ in the LIP, and only lower limits if the gas is in the HIP.

• 3)

  We found only one position, SE1, in the LIP, where R_D = 0.032−0.05, x_e = (0.3−4) × 10^-7 and ζ = (0.2−20) × 10^-17 s^-1. Two positions, N5 and N6, lie in the gap between the LIP and HIP, namely the gas is neither entirely in the LIP nor in the HIP, although it certainly contains a fraction of gas in the LIP, where DCO⁺ is detectable (and detected). The jump from the HIP to the LIP when penetrating farther into the cloud is associated with an increase in the temperature, and we showed that model calculations at several temperatures further constrain the value of ζ. The uncertainty in ζ towards these positions is dominated by the uncertainty in the H₂ density and the derived values are ζ = (1.3−3.3) and (1.3−4.0)× 10^-15 s^-1 for N5 and N6, respectively. Towards the remaining 5 positions with upper limits on R_D, the derived ζ values are at least 10 to 260 times higher than the standard value of 1 × 10^-17 s^-1.

• 4)

  The points of the northern cloud have the largest CR ionisation rates measured so far in the Galaxy. The point towards the southern cloud is, on the contrary, consistent with the average galactic CR ionisation rate of molecular clouds not interacting with a SNR. Since the northern and southern clouds have projected distances from the SNR shock of ≤3 and ~10 pc, respectively, this can be explained by the fact that the low energy ionising CRs have not reached the southern cloud yet. On the other hand, the observations show that both the northern and southern clouds coincide with TeV emission sources, suggesting that high ≳10 TeV CRs have reached both. This is also consistent with γ-ray emission sources inciding with the northern cloud but only partially with the southern cloud, indicating that the former is irradiated by ≈0.1−1 GeV CRs, while only the nearest portion of the southern cloud is so affected.

• 5)

  The emerging picture is that of energy-dependent diffusion properties of hadronic CRs. The high-energy CRs responsible for TeV γ-ray emission through π⁰-decay can diffuse far ahead of the SNR shock, while the low-energy CRs (0.1−1 GeV), responsible for both the low γ-ray emission and the ionisation of the gas, remain closer to the SNR shock. The present work thus gives first observational evidence to the theoretical predictions that hadrons of energy 0.1−1 GeV contribute most to the ionisation in dense gas (Padovani et al. 2009).

Acknowledgments

We warmly thank Marco Padovani for providing us with his compilation and for useful discussions. This work has been financially supported by the Programme National Hautes Énergies (PNHE). Based on observations carried out with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). S. Gabici acknowledges the financial support of the UnivEarthS Labex Program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02).

References

Appendix A: Chemical models

We solved the OSU 2009⁴ chemical network using the astrochem⁵ code.

astrochem is a numerical code that computes the time-dependent chemical abundances in a cell of gas shielded by a given visual extinction A_V and with given physical parameters: the total H density n_H, the gas kinetic temperature T_kin, and the dust temperature T_d. It also takes as an input the initial chemical abundances and the CR ionisation rate. We followed the abundance until a steady state was reached, for a grid of models covering a large range of physical conditions (n_H,T_kin) and CR ionisation rates ζ, at a given A_V = 20 mag, far inside the cloud, where the gas is shielded from the UV radiation field and the ionisation is dominated by CRs. The results only depend on the ζ/n_H ratio (see Sect. 5.2). The role of the dust in astrochem is limited to the absorption and desorption processes, namely no grain surface chemistry is considered. As discussed in Sect. 5.2, neither process is relevant to the present discussion. Initial conditions were taken from the low-metal abundances as in Graedel et al. (1982) and Wakelam et al. (2006b) using an updated He/H relative abundance of 0.09 (Asplund et al. 2009). In Table A.2, we list the range of physical parameters used in this study.

The OSU network contains 6046 reactions involving 468 species. We appended 12 reactions labeled 7−14 in Table A.1 involving the deuterated species: D, HD, H₂D⁺, and DCO⁺. Chemical rates are taken from Roberts & Millar (2000).

All Tables

All Figures

	Fig. 1 Left: the W28 complex on large scales. Grayscale (in σ) and thick contours show TeV emission as seen by HESS (levels are 4−6σ). Red contours show the CO(1−0) emission (Dame et al. 2001) integrated over 15−25 km s-1 and magenta contours trace the emission integrated over 5−15 km s-1 (levels are 40−70 K km s-1 by 5 K km s-1). Crosses show the positions observed with the IRAM 30m telescope and discussed in this paper. The blue contours show the 20 cm free-free emission in the M20 region (Yusef-Zadeh et al. 2000). The blue circle gives the approximate radio boundary of the SNR W28 (Brogan et al. 2006). Right: the northern cloud in the W28 complex (zoom on the black box). The red contours show the CO(3−2) emission in K km s-1, integrated over 15−25 km s-1 (levels are 15−130 K km s-1 by 5 K km s-1) (Lefloch et al. 2008). Diamonds show the locations of OH masers in the region (Claussen et al. 1997).
Open with DEXTER
In the text

	Fig. 2 Observations of millimetre emission lines towards all positions. The intensities are in units of main-beam temperature (K). For readibility, a multiplicative factor was applied to the spectra and is given under each transition line label. This factor was decreased for position N2 (protostar) where the signal is very strong. Vertical dashed lines indicate the two extreme velocity components of the complex at 7 and 21 km s-1.
Open with DEXTER
In the text

	Fig. 3 xe = n(e−) /nH as a function of RD. Values from the analytical method (G77) in Sect. 5.1 are given at T = 20 K (dashed line). Results from our calculation (see Sect. 5.2) for temperatures between 7 K and 20 K are contained within the solid lines. ζ values are given for nH = 104 cm-3, although the model only depends on the ζ/nH ratio (see text). RD values or upper limits as derived from observations are indicated by the red symbols; see Sect. 6.1.
Open with DEXTER
In the text

	Fig. 4 RD = DCO+/ HCO+ (thick line, left axis) and x(DCO+) = n(DCO+) /nH (dashed line, right axis) as a function of ζ, for T = 10 K and nH = 8 103 cm-3, i.e. the physical conditions characterising position N5. The HIP and LIP are marked. The hatched area shows the range of observed RD at that position.
Open with DEXTER
In the text

	Fig. 5 RD as a function of the gas temperature Tkin for different values of ζ/nH: from 2 to 5 × 10-19 s-1, as marked. Note that for ζ/nH ≤ 2 × 10-19 cm3 s-1 (thick solid line), the cloud is always in the LIP, regardless of the temperature. For ζ/nH> 5 × 10-19 cm3 s-1 (thin dashed curve), the cloud is always in the HIP for temperatures ≤50 K. Hatched areas show observations of N5 and N6. We assume AV = 20 mag.
Open with DEXTER
In the text

	Fig. 6 Compilation of measured ζ in different objects (open squares), as reported by Padovani & Galli (2013). The black filled square denotes W51 (Ceccarelli et al. 2011). Red points and lower limits report the values derived in this work. The dashed lines show the range of column densities (0.5−10) × 1022 cm-2, typical of dense molecular clouds, corresponding to visual extinctions of 5 and 100 mag, respectively. On the left lie the diffuse clouds and on the right highly obscured environments such as infrared dark clouds or protoplanetary discs.
Open with DEXTER
In the text

	Fig. 7 CR ionisation rate ζ as a function of the approximate projected distance from the SNR radio boundary (blue circle in Fig. 1), assuming a W28 distance of 2 pkc. We note that the ζ error bars are dominated by the uncertainties on the H2 densities (see text).
Open with DEXTER
In the text