Parsec-scale cosmic-ray ionisation rate in Orion
A. Socci G. Sabatini M. Padovani S. Bovino A. Hacar
 


   

   

   

   Institute for Astronomy (IfA), University of Vienna,
              Türkenschanzstrasse 17, A-1180 Vienna

              andrea.socci@univie.ac.at
         INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
          Chemistry Department, Sapienza University of Rome, P.le A. Moro, 00185 Rome, Italy
         INAF, Istituto di Radioastronomia, Via Gobetti 101, I-40129 Bologna, Italy
        Departamento de Astronomía, Facultad Ciencias Físicas y Matemáticas, Universidad de Concepción Av. Esteban Iturra s/n Barrio Universitario, Casilla 160, Concepción, Chile
            

   



 
  
  
   Cosmic rays are a key component of the interstellar medium since they regulate the dynamics and the chemical processes in the densest and coldest regions of molecular clouds. Still, the cosmic-ray ionisation rate of H_2 () is one of the most debated parameters characterising molecular clouds, due to the uncertainties in the adopted chemical networks and the analysis techniques.
  
   This work aims to homogeneously estimate the  at parsec scales, towards the Orion Molecular Clouds OMC-2 and OMC-3. We explore the change in  across a whole star-forming region by probing a range of column densities never explored before. The significant increase in statistics obtained by studying an entire region allows us to put stronger constraints on the range of  values and exploit its connection with the physical properties of the ISM.
  
   The most recent  estimates are based on o-H_2D^+, a direct product of the interaction between cosmic rays and H_2 in cold clouds. Since observations of o-H_2D^+ are challenging, in this work we proxy its abundance through CO depletion by employing C^18O (2-1) observations towards OMC-2 and OMC-3, taking advantage of the existing correlation between the two parameters. Using additional observations of HCO^+ (1-0) and DCO^+ (3-2), we determine the deuteration fraction and, we finally derive the map of  in these two regions.
  
   The C^18O depletion correlates with both the total column density of H_2 and the N_2H^+ emission across OMC-2 and OMC-3. The obtained depletion factors and deuteration fractions are consistent with previous values obtained in low- and high-mass star-forming regions. In addition, they show a functional dependence within the coldest fields in our sample. The cosmic-ray ionisation rate we derive shows values of ∼5×10^-18-10^-16 s^-1, in good agreement with previous estimates based on o-H_2D^+ observations. The  also shows a functional dependence with the column density of H_2 across a full order of magnitude (∼10^22-10^23 cm^-2). The estimated values of  overall decrease for increasing N(H_2), consistently with the predictions from theoretical models.
  
   The approach we follow delivers results comparable with theoretical predictions and previous independent studies, confirming both the robustness of the analytical framework as well as CO depletion being a viable proxy of o-H_2D^+. We also explore the major limitations of the method by varying the physical size of the gas crossed by the cosmic rays (i.e., the path length). By employing a path length obtained from low resolution observations, we recover values of the  well below any existing theoretical and observational prediction. This discrepancy suggests interferometric observations as mandatory in order to reliably constrain the  also at parsec scales.

   

   Parsec-scale cosmic-ray ionisation rate in Orion
    A. Socci
          1
          G. Sabatini2
          M. Padovani2
          S. Bovino3,4,5
          A. Hacar1
          
    —
==============================================================================================================================



labelfont=bf




§ INTRODUCTION



Cosmic rays (CRs) are a fundamental component of the interstellar medium (ISM), yet, possibly, the most puzzling. Low-energy CRs (E<1 TeV) can penetrate the cold (T≲20 K), dense gas (n≳10^4 cm^-3) of molecular clouds and interact with H_2. The chain reactions following the H_2 ionisation produce H_3^+ <cit.>, a key molecule responsible for the production of most ionic species within molecular clouds <cit.>. The amount of ions over the number of neutral species, namely the ionisation fraction, is of paramount importance for star formation, as it regulates the coupling between magnetic fields and the gas <cit.>, slowing or even halting the gravitational collapse of the cloud <cit.>. In addition, the amount of ions, H_3^+ in particular, drives the chemical evolution of the coldest regions within a cloud <cit.>. However, while observed in infrared absorption in diffuse clouds <cit.>, H_3^+ cannot be detected in the mm regime, preventing a direct measure of the degree of ionisation in dense molecular clouds. The cosmic-ray ionisation rate (CRIR), directly connected to H_3^+, has been therefore determined over the years through several proxies, such as radicals <cit.>, ions <cit.>, neutrals <cit.>, and a combination of the latter <cit.>. The values of CRIR determined in molecular clouds span three orders of magnitude (∼10^-17-10^-14 s^-1) and show a dependence with the column density of H_2, N(H_2), as expected from theoretical models <cit.>. Yet, the variety of techniques (most of them model-dependent), tracers and resolutions makes for a non-homogeneous sample of scattered estimates, especially in dense clouds <cit.>.


<cit.> has recently proposed ortho-H_2D^+ (hereafter o-H_2D^+) as a proxy for the CRIR relative to molecular hydrogen (). H_2D^+ is the first product of the deuteration of H_3^+ and its ortho form has been successfully detected in the millimetre regime at both low <cit.> and high resolution <cit.>. Since the deuterium enrichment is favoured in cold dense gas <cit.>, o-H_2D^+ is a reliable proxy of the  in young pre-stellar environments <cit.>. This analytic derivation has been recently employed by <cit.> and <cit.>, who determined  from o-H_2D^+ in a sample of high-mass clumps and cores. The ionisation rates obtained by the authors are consistent with the theoretical predictions for a high-density regime <cit.>. Whilst the estimated  reproduces the theoretical expectations, the detection of o-H_2D^+ has been a challenge for more than a decade <cit.>. As a consequence, the method proposed by <cit.> has been applied only to a limited range of column densities and scales, in particular for N(H_2)≳10^23 cm^-2 and R≲0.5 pc, respectively <cit.>. 


The present study aims to derive the column density of o-H_2D^+ and then the  for column densities within N(H_2)∼10^22-10^23 cm^-2 and for volume densities around n∼10^5 cm^-3. In the cold gas of molecular clouds where these volume densities are achieved, CO depletes onto dust grains, ions become abundant and deuteration processes are overall boosted <cit.>. In this context, the abundance of H_2D^+ is expected to grow for increasing depletion factors, as also observed in different samples of both low- <cit.> and high-mass <cit.> star-forming regions. Following the aforementioned studies, in this work we determine N(o-H_2D^+) from the C^18O depletion factor to provide an independent estimate of  on parsec scales. Despite its limitations, discussed throughout the paper, our approach produces results consistent with previous works, while extending the dynamic range of column densities explored so far in the  estimation.


The paper is then structured as follows: we introduce the method for the  determination and its main assumptions (Sect. <ref>); then, we present the observations, both novel (Sect. <ref>) and archival (Sect. <ref>), used throughout the paper and the related source selection; next, we describe the analysis and results obtained for the degree of depletion (Sect. <ref>), the expected abundance of o-H_2D^+ (Sect. <ref>), the deuteration fraction (Sect. <ref>), and the ionisation rate (Sect. <ref>), along with their error budget (Sect. <ref>) and limits of applicability (Sect. <ref>); finally, we discuss the connection between  and N(H_2) both from our analysis and in comparison with the theoretical predictions, together with the major limitations to the ionisation rate estimate (Sect. <ref>). Section <ref> wraps up the major findings of the present study.




§ METHOD



The method introduced by <cit.> provides a model-independent analytical expression to estimate , which reads as follows <cit.>:

    = k_CO^o-H_3^+N(o-H_2D^+× N(CO)/3 R_D× N(H_2)× l.

In the above equation: k_CO^o-H_3^+ is the destruction rate of o-H_3^+ by CO, considered here as main destruction path for o-H_3^+ <cit.>; N(CO) is the total column density of CO; R_D (= N(DCO^+)/N(HCO^+)) is the deuteration fraction; l is the path length over which CO suffers depletion and we estimate the column densities of our tracers; N(o-H_2D^+) and N(H_2) are the column densities of o-H_2D^+ and H_2, respectively. 


The robustness and reliability of Eq. (<ref>) have been recently confirmed by <cit.> through synthetic observations of the molecular tracers simulated in a three-dimensional grid. The analysis provided  estimates accurate within a factor of 2-3 compared to the input value of the simulation. These findings confirm the confidence level suggested by <cit.> for the method. 


To apply Eq. (<ref>), our study relies on two main assumptions. First, we assume o-H_2D^+ to be efficiently traced by the C^18O depletion factor (f_D). This hypothesis is based on the correlation between the o-H_2D^+ abundance and the degree of CO depletion, a result consistently derived in earlier studies both through observations <cit.> and through theoretical predictions <cit.>. We will therefore proxy the abundance of o-H_2D^+ via f_D using the correlation reported by <cit.>. 


Second, we make an assumption on the path length l in Eq. (<ref>). This path length is the physical extent of the o-H_2D^+ emission along the line of sight. Largely unknown, l is usually taken as the physical size of the o-H_2D^+ emission in the plane of the sky <cit.> or in the box where it is simulated <cit.> and constitutes one of the main limitations of the present method. We make use of l=0.050±0.015 pc, the typical width of filamentary substructures seen in N_2H^+ across OMC-2 and OMC-3 (Socci et al. submitted). N_2H^+ is a molecular ion which becomes abundant when CO freezes-out on grains <cit.> and whose emission has been proven to correlate with the one of o-H_2D^+ <cit.>. Thus, we take it as the representative scale for our study. Although unresolved in our observations (see Sect. <ref>), this scale is the most accurate estimate in the regions without relying on direct o-H_2D^+ observations. In Sect. <ref>, we will discuss how different assumptions on l affect the determination of . 



§ OBSERVATIONS




 §.§ IRAM-30m observations





Applying the method from <cit.> to estimate  outside high-mass cores and clumps implies exploring column densities N(H_2)≲10^23 cm^-2 at parsec scales. To fulfil this goal, the Intregral Shaped Filament (ISF) in Orion <cit.> appeared as a natural choice. The ISF is a ∼13 pc long filament in the Orion A molecular cloud, prominent in both dust continuum <cit.> and molecular lines <cit.>, and composed by several condensations, the Orion Molecular Clouds <cit.>. The ISF is also the closest site of high-mass star-formation in the Solar neighbourhood <cit.> harbouring the Orion Nebula Cluster <cit.> within OMC-1. The extended size, the large dynamic range in  H_2 column densities and the connection to high-mass star formation make the ISF the best candidate for the present study.


We surveyed OMC-1, OMC-2, OMC-3 and OMC-4 in C^18O (2-1) (219.560 GHz), DCO^+ (3-2) (216.113 GHz) and HCO^+ (1-0) (89.189 GHz) using single-dish observations at 1mm and 3mm. The observations were carried out at the 30-metre Institut de Radioastronomie Millimétrique telescope (IRAM-30m) in Granada (Spain) during November 2013 (Project: 032-13, PI: A. Hacar). The large-scale IRAM-30m mosaics centred on Orion BN/KL <cit.> are obtained with single 200”×200” or 100”×100” Nyquist-sampled On-The-Fly (OTF) maps in Position Switching (PSw) mode <cit.>. Our suite of tracers have been observed with the same receiver, the Eight Mixer Receiver <cit.>, but with different backends: C^18O (2-1) was observed with the VErsatile SPectrometer Array (VESPA) at 20 kHz resolution (∼0.027 km s^-1 at 219.560 GHz), while DCO^+ (3-2) and HCO^+ (1-0) with the Fast Fourier Transform Spectrometer <cit.> at 200 kHz resolution (∼0.27 km s^-1 at 216.113 GHz and ∼0.66 km s^-1 at 89.189 GHz, respectively). The spectra obtained in antenna temperature (T_a^*) have been converted in main beam temperature (T_mb) with the relation T_mb = T_a^* F_ef/B_ef. Here we used the standard forward (F_ef) and backward (B_ef) efficiencies for the IRAM-30m telescope[<https://publicwiki.iram.es/Iram30mEfficiencies>] (C^18O (2-1): B_ef = 0.61, F_ef = 0.93; HCO^+ (1-0): B_ef = 0.81, F_ef = 0.95; DCO^+ (3-2): B_ef = 0.62, F_ef = 0.93). Finally, all the observations have been convolved to a final resolution of 30” (∼0.06 pc at the Orion distance). Figure <ref> shows the integrated intensity maps of the three molecules obtained by integrating the spectral cubes in the same velocity range Δ V_lsr = [6, 15] km s^-1, around the typical velocity expected for the ISF <cit.>. Similarly, the rms noise has been estimated across the whole map, integrated in spectral windows with same velocity range but devoid of emission.


C^18O and HCO^+ are readily detected and show widespread emission across the ISF. The high-sensitivity IRAM-30m observations grant a peak-to-noise dynamic range of ≳100 for the intensities of both tracers. The more extended component of C^18O and HCO^+ emission across the ISF shows values on average of ∼6-10 K km s^-1. These values increase along the ISF crest up to peaks of ∼35 K km s^-1, for C^18O and ∼100 K km s^-1, for HCO^+, seen towards OMC-1. DCO^+, on the other hand, shows a limited peak-to-noise dynamic range of ∼50, along with a patchy and uneven emission across the ISF. The few peaks with intensities of ≳4 K km s^-1 are localised towards OMC-1 and OMC-2/OMC-3. All three tracers show their peak intensities towards OMC-1, where the radiation coming from the Trapezium stars <cit.> and the ONC heats up the gas to temperatures above ∼40 K <cit.>. While all tracers share OMC-1 as the region of peak intensity, the C^18O emission is often anti-correlated across the ISF compared to the emission of HCO^+. This change in morphology is prominent towards two regions: first, along OMC-2/OMC-3, where the C^18O emission is scattered and fragmented while HCO^+ shows a coherent filamentary shape (see Fig. <ref>, panels a, b); second, in OMC-4, where C^18O emission is instead brighter compared to the one of HCO^+. The behaviour in OMC-2/OMC-3 is of interest as it hints towards the presence of extended C^18O depletion in these two regions.



 §.§ Archival observations and source selection



The ISF is a well-known region with a plethora of archival data available. For our study, we will take advantage of a few of these ancillary observations with comparable resolution to our data. Among these, we are interested in a N(H_2) map of the region to determine the candidate sites for depletion. Large scale depletion is in fact expected for column densities N(H_2)≳10^22 cm^-2 <cit.>. The column density map of H_2 has been obtained from the 850 μm dust opacity, observed by the Herschel and Planck observatories, following the prescriptions of <cit.> (see also Appendix <ref>). Figure <ref>, panel (a), shows the ISF as seen in column density of H_2 (36.2”). The map shows material with column densities N(H_2)≥10^22 cm^-2 extended across the whole filament with peaks of N(H_2)∼10^23 cm^-2 along its crest, in correspondence of the OMCs. We will focus on these regions, where the column density has a dynamic range of almost two orders of magnitude and large CO depletion is expected.


As mentioned in the previous section, N_2H^+ emission strongly anti-correlates with the one of CO, as it gets destroyed by the latter <cit.>. To further assess the presence of large-scale depletion in the ISF, we qualitatively compare the total column density of H_2 with the IRAM-30m observations of N_2H^+ (1-0) from <cit.> (Fig. <ref>, panel b). Although not as extended as the material with N(H_2)∼10^22 cm^-2, N_2H^+ has emission widespread across the ISF and closely correlated with N(H_2) towards the OMCs (panel a). All four regions are thus candidates for efficient depletion and, as a consequence, good targets to estimate .


Our focus is instead shifted towards OMC-2 and OMC-3 only, as suggested previously. There are two main reasons to consider only these two among the four OMCs. First, the OMC-1 region hosts several O-B stars <cit.> and a plethora of highly embedded outflows <cit.>. The radiation and heating from these objects can modify the C^18O abundance by its selective photodissociation or its release in the gas phase from the ice coatings of dust grains <cit.>. For a reliable estimate of the depletion, we have to exclude OMC-1 where the above processes are the most efficient in the ISF. Second, DCO^+ (3-2) has weak emission across the ISF (see Sect. <ref>) and, in particular, no emission above a signal-to-noise of 3 towards OMC-4. We are therefore limited in our analysis to OMC-2 and OMC-3 <cit.>.


OMC-2 and OMC-3 themselves show more than an order of magnitude variation in N(H_2) and temperatures between T∼10-40 K, according to different estimates <cit.>. Since outflows are also detected in these two regions <cit.>, we minimise their contamination to the C^18O depletion estimate by imposing a temperature cut in the analysis of T_K = 10-25 K <cit.>. Given the close correlation between NH_3 and N_2H^+ <cit.>, we expect T_K(NH_3) to be the representative temperature for the gas in which CO freezes-out onto the grains. The range T_K = 10-25 K is consistent with those previously explored in studies of f_D in high-mass regions <cit.>. In addition, ∼20-25 K corresponds to the temperature limit for which CO evaporates from dust grains <cit.>, given n≳10^5 cm^-3, volume densities we expect to explore with the present work.


OMC-2 and OMC-3 match all the previous characteristics with a clear detection of DCO^+, a low number of O-B stars nearby, a few scattered outflows, column densities up to N(H_2)∼10^23 cm^-2 and temperatures mostly below ≲25 K. Across these ∼1-2 pc towards the North of the ISF, we expect C^18O to be depleted and not photodissociated by stellar activity. Under such conditions, the depletion is expected to correlate to the o-H_2D^+ abundance <cit.>, and, therefore, it may be used to estimate the .



§ ANALYSIS AND RESULTS




 §.§ The depletion map of OMC-2 and OMC-3





As we want to determine the onset of depletion, properly estimating the column density of C^18O (N(C^18O)) is mandatory. First insights in this direction come from the inspection of both I(C^18O), I(N_2H^+) and their comparison with N(H_2). We can qualitatively explore the N(H_2) regime for which C^18O depletes by computing the temperature corrected line intensity <cit.>:

    I_T = I/J_ν(T_K) - J_ν(T_bg).


Through Eq. (<ref>), the column density of H_2 is now the only factor contributing to the line excitation, therefore to its intensity I_T. Figure <ref> shows the comparison between I_T and N(H_2) for our suite of tracers. We will now focus on C^18O (panel a) and N_2H^+ (panel b): these two tracers show a continuous distribution with N(H_2), however they behave differently for increasing column densities. Across the dynamic range probed by our analysis (N(H_2)∼10^22-10^23 cm^-2), the flattening of I_T(C^18O), on one hand, and the linear increase of I_T(N_2H^+), on the other, represent a clear chemical differentiation of the two tracers and, in particular, the effect of depletion on C^18O <cit.>.


To quantitatively describe the degree of C^18O depletion in OMC-2 and OMC-3, we need to determine its column density (N(C^18O)). The latter can be computed as follows <cit.>:

    N_tot = 8π kν^2/hc^3A_ulQ(T_ex) e^E_ul/kT_ex/g_uI τ/1-e^-τ,

with the corresponding parameters: h, the Planck constant; k, the Boltzmann constant; c, the light speed; ν, the rest frequency of the line; A_ul, the Einstein coefficient for spontaneous emission; Q, the partition function calculated at the excitation temperature of the line (T_ex); g_u, the degeneracy of the upper level; E_u, the excitation energy of the upper level; τ, the optical depth of the line. We can simplify the above equation for N(C^18O) by introducing two approximations. First, we can assume C^18O (2-1) to be in Local Thermodynamic Equilibrium (LTE) at the gas kinetic temperature T_K. Given the presence of N_2H^+ in the gas phase, the volume densities are expected to be n∼10^5 cm^-3, which is at least an order of magnitude larger than the critical density of C^18O <cit.>. We therefore take T_ex(C^18O) = T_K(NH_3) for the rest of the analysis. Second, we work under the optically thin assumption (i.e., τ<<1) adducing the flattening in C^18O intensity only to its freeze-out onto the dust grains <cit.>. With these two approximations, Eq. (<ref>) can be rewritten as:

    N_tot(C^18O) ∼8π kν^2/hc^3A_ulQ(T_K) e^E_ul/kT_K/g_uI(C^18O).



The ratio between the expected abundance of C^18O at the galactocentric distance of Orion <cit.> and the one measured (X_obs = N_tot(C^18O)/N(H_2)) gives the depletion factor:

    f_D = X_0/X_obs.

Figure <ref>, panel (a), shows this depletion factor, which, despite our source selection (see Sect. <ref>), is extended for ∼1 pc from OMC-2 to the tip of OMC-3. Depletion factors f_D>1 correlate extremely well with N_2H^+ intensities above I(N_2H^+)>4 K km s^-1 (signal-to-noise of 10, black contours) and with the column densities N(H_2)≳10^22 cm^-2 across the two regions (see Fig. <ref>). Our map shows depletion factors ranging within f_D∼1-5, consistently with previous studies at comparable resolution (i.e., ∼30”) in infrared dark clouds <cit.>, single low- <cit.> and high-mass <cit.> star-forming regions. Similarly to these studies, we also determine f_D slightly below 1 (∼0.85 on average), which would imply an over-abundance of C^18O compared to the expected X_0. Since these fields are restricted to the edges of the map and given the small deviations from unity, such f_D<1 could result from uncertainties in their derivation rather than a physical over-abundance of C^18O.
We have therefore masked their corresponding fields in Fig. <ref>, panel (a), to consider only the regions where depletion is effective (f_D≥1), which constrains our map to the filamentary structure extended throughout OMC-2 and OMC-3.





 §.§ Expected abundance of o-H_2D^+ on parsec scales



The challenging observations of o-H_2D^+ have limited its study to surveys of individual cores and clumps within low- <cit.> and high-mass <cit.> star-forming regions. The opportunity to proxy its abundance with f_D, under proper conditions, opens up the path for parsec-scale studies of the ionisation processes. We thus apply the empirical correlation between f_D and X(o-H_2D^+) determined by <cit.>:

    log_10(X(o-H_2D^+)) = 0.05× f_D - 10.46

Figure <ref>, panel (b) shows the X(o-H_2D^+) estimated with Eq. (<ref>). The dynamic range of X(o-H_2D^+) is constrained within a factor of 3 (∼3-8×10^-11) and in excellent agreement with the same abundances from <cit.>, as expected since the depletion factors are also in close agreement. These abundances of o-H_2D^+ are broadly consistent with independent estimates for low-mass cores <cit.>, however, they are an order of magnitude, on average, lower than the values determined for high-mass regions <cit.>. The difference between our estimates and those in high-mass regions can be explained based on two major factors. First, the higher column densities probed in these regions (N(H_2)≳5×10^23 cm^-2) compared to our work. Second, the beam dilution effect when comparing our low-resolution IRAM-30m observations (30”) to their high-resolution ALMA observations (∼2”). A comparison between observations and 3D modelling of the CO freeze-out, in fact, suggests depletion factors of f_D∼3-8 to be typical at scales of ∼0.1 pc <cit.>. These scales are close to the resolution of our maps (∼0.06-0.07 pc), thus, since X(o-H_2D^+) is inferred from f_D, it suffers from the same beam dilution effect. 


The o-H_2D^+ abundances we derive are mostly in agreement with those determined in previous studies, highlighting the CO depletion factor as a viable proxy for X(o-H_2D^+). However, independent measurements of these two parameters are needed to ultimately break the degeneracy and evaluate the empirical correlation found by <cit.>. The increasing sensitivity and spatial resolution of modern astronomical facilities, e.g., ALMA, will allow us to achieve this goal in the near future.



 §.§ Deuteration fraction across OMC-2 and OMC-3



The deuteration fraction (R_D) is directly connected to all D-bearing isotopologues of H_3^+. These molecules are in fact the main deuterium donors in the cold and dense gas of molecular clouds <cit.>. In the present study, we estimate the deuteration fraction as N(DCO^+)/N(HCO^+) through the DCO^+ (3-2) and HCO^+ (1-0) lines. Our first approach could be again estimating the column densities of both tracers using Eq. (<ref>) in the optically thin limit (τ≪1) and assuming LTE at the kinetic temperature T_K. While reasonable for C^18O, assuming LTE for DCO^+ and HCO^+ is not ideal given the larger critical densities of these two molecules <cit.>. In addition, the HCO^+ (1-0) can easily become optically thick <cit.>, therefore a correction for τ is expected.


Similarly to C^18O, we inspect the temperature-corrected intensities of both HCO^+ and DCO^+, and compare them to N(H_2) (see Fig. <ref>). I_T(HCO^+) (panel c) shows a clear dependence with N(H_2) with a compact distribution. Only a handful of fields show a significant scatter, suggesting self-absorbed profiles to be almost absent in our sample <cit.>. However, while not self-absorbed, the HCO^+ intensities likely suffer from opacity as the distribution flattens towards N(H_2)∼10^23 cm^-2. I_T(DCO^+) (panel d) also shows a clear dependence with N(H_2). Compared to those of HCO^+, the intensities of DCO^+ increase almost linearly with N(H_2), for a distribution resembling the one of N_2H^+. However, I_T(DCO^+) shows a higher degree of scatter compared to I_T(N_2H^+), suggesting deviations from the LTE assumption used to apply Eq. (<ref>).


To estimate the excitation effects on both DCO^+ and HCO^+, we run a grid of radiative transfer models per each pixel using RADEX <cit.> in two different steps[The radiative transfer calculations have not been applied to C^18O because, without proper modelling of the cloud which includes a freeze-out threshold, depletion may be mistaken for opacity <cit.>]. In the first step, we input N_tot(DCO^+) and N_tot(HCO^+) calculated in the LTE, optically thin limit and vary n, the volume density, in logarithmic intervals between 10^4 and 10^6 cm^-3. In the second step, we fix these same volume densities and let the two column densities vary by an order of magnitude around the LTE, optically thin values. We then take from the output of RADEX the two N_tot which minimise the observed intensities. In both steps, we input T_K as the kinetic temperature and fix the linewidths of DCO^+ and HCO^+ to Δ v=0.75±0.10 km s^-1 and Δ v=2.5±0.4 km s^-1, respectively. These linewidths are the mean values obtained by fitting a few representative spectra in our fields and their error is taken as the standard deviation seen across these fits.


The deuteration fraction (Fig. <ref>, panel c) shows more than an order of magnitude variation (∼0.005-0.1) across OMC-2 and OMC-3. Such values are consistent with previous studies in high-mass star-forming regions <cit.>, using the same molecular species, and theoretical predictions from chemical models <cit.>. The connection found between f_D and R_D in these and other studies <cit.> demonstrates how the deuteration process is enhanced in the cold and dense gas where CO is depleted <cit.>. Simulations of turbulent and magnetised pc-scale filaments, such as the ISF <cit.>, show increasing levels of deuteration as the column density increases <cit.>. These authors determine deuteration fractions within ∼0.01-0.1 in t≲0.4 Myr, a timescale consistent with the expected survival time of the dense gas in Orion <cit.>. Since o-H_2D^+ is the driver of deuteration in the early phases of the cloud evolution <cit.>, and given the dependence between f_D and R_D, a correlation also between X(o-H_2D^+) and f_D would be readily explained.




The variation in R_D across our maps shows a progressive increase from OMC-2, in the South, to OMC-3 towards the North. Such gradient reflects the morphology of the HCO^+ and DCO^+ integrated intensities seen in Fig. <ref> (a,c). The variation could be linked to a change in volume density causing the higher-J levels of the two molecules to be more populated. Yet, the non-LTE analysis returns a limited range in volume densities, varying within ∼0.3-2×10^5 cm^-3 (with mean value and standard deviation of n=1.0±0.2×10^5 cm^-3). Therefore, the significant variation in R_D is likely the result of a more efficient deuteration process towards OMC-3. 


It is well known <cit.> that deuteration is affected by temperature variations, which cause a change in the reaction rates. The deuteration of HCO^+ into DCO^+ is particularly sensitive to this temperature effect <cit.>. Thus, the gradient we observe in R_D is possibly linked to the slight change in T_K from OMC-3 towards OMC-2 (see Fig. <ref>, panel c). We explore this hypothesis in Fig. <ref>, where we display the comparison of R_D with T_K (left panel) and f_d (right panel). The left panel shows an anti-correlation between the deuteration fraction and the kinetic temperature, already seen in previous studies <cit.>, most prominent for fields with low depletion factors (f_D≲2). When the temperatures increase towards and above T_K≳20 K, the deuteration fraction drops to values of ≲0.01. For temperatures ∼16-19 K however, there is a significant scatter, for which the fields with higher depletion (f_D≳2-3) also show a higher degree of deuteration (R_D≳0.03). These same points can be seen in the right panel as they show a positive correlation between R_D and f_D. This dependence between the two parameters has been previously reported by <cit.> for a similar range of values and volume density regime (n∼10^5 cm^-3). We therefore expect the correlation between R_D and f_D to be robust below ≲20 K, and become more loose for higher temperatures making the corresponding  estimates possibly uncertain.


As a final test, we determine a scaling relation between R_D and T_K in Fig. <ref>. First, we sample the median value, and Inter-Quartile Range (IQR), of R_D and T_K in bins of one Kelvin (red circles and error bars, respectively). Then, we fit a power law dependence to these estimates, for which the best fit (black line) reads as R_D∼ T_K^-4 (see Eq. <ref>). Such sharp decrease of R_D(HCO^+) with T_K has been previously observed across different clouds <cit.>, although never before condensed in a scaling relation. The direct comparison of this scaling relation with estimates of R_D(HCO^+) and T_K in massive clumps <cit.> and high-mass regions <cit.> confirms the goodness of our fit (see full discussion in Appendix <ref>). The systematic nature of the effect, which appears as independent from the star-formation regime of the region, suggests the use of temperature-corrected deuteration fractions for future estimates of  using the method described in Sect. <ref>.



 §.§ The destruction rate k^o-H_3^+_CO



Equation <ref> assumes o-H_3^+ to be primarily destroyed by CO <cit.>. While this is one of the main processes in dense cores, other destruction pathways can occur when studying the ionisation rate at parsec scales. In particular, electrons can reduce the H_3^+ abundance through dissociative recombination. This process becomes important both at extreme densities, when the CO is almost entirely depleted <cit.>, and at lower densities, especially in a region externally illuminated by stellar radiation, such as the ISF <cit.>.


We thus consider estimates of the ionisation fraction, x(e^-), in these two density regimes and compare them with the analytical limit f_D < 9×10^-7/x(e^-) proposed by <cit.> for Eq. (<ref>). First, we consider the estimates in massive dense cores, including some in Orion <cit.>. These authors estimate ionisation fractions x(e^-)≲10^-7 at N(H_2)≳10^22 cm^-2 for these sources, leading to a depletion limit for the method of f_D≳9. Second, we consider the estimates in OMC-2 and OMC-3 using the C_2H and HCN lines <cit.>. These authors determine ionisation fractions x(e^-)≤3×10^-6 at n∼5×10^3 cm^-3, leading to a non-physical f_D>0.3. Given our range in both column and volume density (N(H_2)∼10^22-10^23 cm^-2 and n(H_2)∼0.3-2×10^5 cm^-3, respectively), we can safely assume the destruction through CO as main process affecting the H_3^+ abundance for our degree of depletion.


We can finally determine the destruction rate k^o-H_3^+_CO in our regions. We computed the latter pixel-by-pixel using T_K and the coefficients found in the KInetic Database for Astrochemistry <cit.> for the ion-polar reaction . The adopted k^o-H_3^+_CO rates, in the temperature range of 10-25 K, are relatively constant and lie within ∼2.15-2.35×10^-9 cm^-3 s^-1.



 §.§ Error budget in the determination of 



Our indirect method to determine  presents a number of sources of error, which will be discussed in the following. Aside from X(o-H_2D^+), for which we assume the CO depletion is an efficient proxy (see Sect. <ref>), all the other parameters in Eq. (<ref>) bear a corresponding uncertainty.

    
  * k_CO^o-H_3^+, N(C^18O): both only depend on the temperature (see Eq. <ref> and Sect. <ref>), thus we propagated their error accordingly using the uncertainty on T_K <cit.>;
    
  * [^16O]/[^18O]: used to convert N(C^18O) to N(CO) in Eq. (<ref>), it is equal to 560±26 for the Local ISM <cit.>;
    
  * R_D: given its definition and its derivation with RADEX, we consider as total uncertainty the combined uncertainties on the DCO^+ and HCO^+ linewidths (see Sect. <ref>);
    
  * l: we assume the typical size of N_2H^+ structures in OMC-2 and OMC-3 (see Sect. <ref>), along with its corresponding error (see also Socci et al. submitted).



The error per field in our analysis reaches up to ∼50% of the corresponding  across OMC-2 and OMC-3. However, this error contains a major flat contribution coming from the conservative choices made to estimate the uncertainty on some parameters (e.g., linewidths of DCO^+ and HCO^+). A more thorough analysis would therefore significantly reduce the uncertainty in several of our fields. Since the determination of  remains indirect and this error budget does not hinder our further discussion and conclusions, we stick to these choices for the present analysis.



 §.§ The Orion cosmic-ray ionisation rate map



The map of  (Fig. <ref>, panel d) shows a dynamic range of more than one order of magnitude across the OMC-2/OMC-3 region. The global South-North decrease in  is likely connected to the gradient seen in the deuteration fraction, therefore to the effect of temperature (see Fig. <ref>, left panel). However, the local variations are likely induced by the change in f_ D along the filament (see Fig. <ref>, right panel). As a result of both effects,  varies within ∼5×10^-18-2×10^-16 s^-1 within the two regions.


The median value we derive of ∼2.4×10^-17 s^-1 is in good agreement with the estimates in high-mass star-forming regions using the same method <cit.>. The slight differences reflect the variation in X(o-H_2D^+) (see Sect. <ref>), thus we can ascribe these differences to a resolution effect. Among the  we determine, its highest values (i.e., ≳2×10^-16 s^-1) are seen towards the South of the map and, in particular, in the region close to the OMC-2 FIR4 protocluster <cit.> and the B-type star HD 37060 (see Fig. <ref>, panel d). The presence of these two heating sources can indirectly affect our estimate through an increase in T_K, with a corresponding decrease of R_D. In addition, the outflows from the protocluster <cit.> and the UV-irradiation from the star could promote the photodissociation of C^18O, inducing an increase in f_D (see Fig. <ref>, right panel). However, higher values of  could also result from a local enhancement in the CR flux <cit.> caused by these objects as suggested by previous findings: values of CRs ionisation rates as high as ∼10^-14 s^-1 has been previously reported in OMC-2 FIR4 <cit.>; a candidate expanding CO-shell has been detected in close correspondence of the B-type star <cit.>, for which CRs could travel beyond its edge and promote the ionisation in the South of our filament.


Despite all the influences to the ionisation rate, ≈10^-17 s^-1, originally found by <cit.>, has been instead considered for several decades as the reference value for the ISM. While such value is reproduced by our results, our pc-scale map shows at least an order of magnitude variation for  across OMC-2 and OMC-3. No constant value for the ionisation rate is therefore expected, but instead, local changes in the physical properties of the gas play a crucial role in its variation. Among the possible physical influences on the CR ionisation rate, we consider here the effect of increasing N(H_2) as the prime driver of its variation.



§ DISCUSSION





The different models for energy loss of CR particles (protons, electrons and heavy nuclei) with increasing N(H_2) have been originally studied by <cit.>. These models account for the ionisation by primary CR protons and electrons, plus secondary electrons from primary CRs. For the column density regime we probe, the models of CR propagation (see Fig. <ref>, coloured lines) are largely dominated by CR protons and their secondaries. As a general trend, all models predict a decrease of  with increasing N(H_2). Different slopes come from the comparison between the CR propagation models and observational results: model ℒ, from the CR spectrum observed by the Voyager missions <cit.>; model ℋ, from the average ionisation rate estimated in diffuse clouds <cit.>; model 𝒰, from the upper limit of the ionisation rate in diffuse clouds. Observational results show a decreasing trend as well, however, the large number of methods and tracers used to infer  produces a significant scatter in the data, especially for dense clouds <cit.>.


Our new results considerably increase the statistics for the available  estimates compared to previous studies <cit.>, probing a column density regime barely explored before. Figure <ref>, shows the comparison between  and N(H_2), colour-coded by T_K, with their corresponding errors <cit.>. 
In addition, we plot the models from <cit.> for different slopes of the CR proton energy spectrum (coloured lines). Almost the entire sample of  determined in our study is enclosed within the model ℒ and the model ℋ, with only a fraction of fields outside of this range. Despite their departure from the bulk of our points, these ionisation rates above the model ℋ and below the model ℒ are still consistent with the theoretical predictions within the errors. None of the predicted propagation models clearly reproduces our results. Our estimates show, however, the functional dependence between  and N(H_2) expected from these models. While the temperature dependence of R_D (see Sect. <ref>) contributes to the scatter, the fields with T_K≲20 K show ionisation rates decreasing with N(H_2).


As we include the  determined in two high-mass star-forming regions using the same method <cit.>, the ionisation rates across the two studies are in close agreement. Overall our estimates show values closer to those in AG531, while only the highest ionisation rates we determine are comparable to the ones in AG354. The combined distribution of ionisation rates from the two studies is continuous with N(H_2) and show a mild functional dependence across column densities ∼10^22-5×10^23 cm^-2. The method proposed by <cit.> therefore allows for homogeneous estimates of  in a large range of scales and column density regimes, when o-H_2D^+ is either successfully detected or reasonably inferred.


The two estimates from the present paper and <cit.> are in close agreement, however spurious effects have to be considered in their determination. While resolution on the one hand affects the degree of depletion (see Sect. <ref>), on the other hand it also influences the minimum physical size we can probe with our observations. Such physical size is condensed in the path length l from Eq. (<ref>). As mentioned in Sect. <ref>, measuring either the line-of-sight length or the volume density of the tracers is complex, therefore l is usually taken equal to the radial size of the structures identified in o-H_2D^+ <cit.>. We explore the influence of l on the  by measuring the typical radial size of the DCO^+ emission, our least extended tracer. We therefore perform a radial sampling on the DCO^+ integrated intensity map following the peaks of its emission and fit the average profile with a Gaussian function (see Sect. <ref> for a full discussion). The FWHM of the profile (0.18±0.15 pc) is taken as representative l to compute the estimate of . Being 3-4 times the l adopted in our analysis (0.05 pc; see Sect. <ref>), this FWHM results in  values on average lower by the same factor (see Eq. <ref>, ∝ 1/l). The median ionisation rate with l∼0.18 pc is ∼7×10^-18 s^-1 and the majority of the fields is below the model ℒ. Ionisation rates as low as ∼10^-18 s^-1, produced by correspondingly low-energy CRs, are unrealistic when considering the degree of feedback present in the ISF <cit.> and strongly disagree with independent estimates, which indeed suggest higher values <cit.>. As f_D and R_D are positively correlated, their ratio is only mildly affected by resolution effects. The path length, l, on the other hand, suffers directly from a poor resolution, thus representing the main limitation in the determination of ionisation rates with Eq. (<ref>) <cit.>.


Our results highlight the CO depletion as available proxy for the o-H_2D^+ abundance at parsec scales. However, caution should be taken when interpreting these results: first, determining f_D can be tricky and the use of multiple CO lines is usually advised <cit.>; second, the connection between f_D and the o-H_2D^+ abundance has been derived at sub-parsec scales <cit.> and works within a constrained range of T_K, outside which multiple processes may become temperature dependent (e.g., degree of deuteration). In turn, o-H_2D^+-based estimates are confirmed as an effective method to homogeneously sample  across more than an order of magnitude in column density. Ultimately, in the absence of local sources of CRs, N(H_2) alone appears to drive the ionisation degree in molecular clouds <cit.>. Finally, our results also highlight the importance of high resolution observations when estimating the degree of depletion and to reliably estimate l, which can significantly affect the  determined. Independent interferometric observations of both CO isotopologs and o-H_2D^+ are needed to validate the previous correlation by <cit.> and provide a more robust estimate of the cosmic ray ionization rate.



§ CONCLUSIONS


In this work, we present the first parsec-scale estimate of the cosmic-ray ionisation rate determined by employing the most recent analytical framework proposed by <cit.>. To this end, we presented novel IRAM-30m observations of the Integral Shape Filament in Orion using the following suite of tracers: C^18O (2-1), DCO^+ (3-2) and HCO^+ (1-0). We first determined the C^18O depletion factor under density and temperature regimes for which we expect significant freeze-out and no photodissociation effects. The C^18O depletion factor was then used to proxy the o-H_2D^+ abundance, following the correlation determined by <cit.>, and ultimately determine the ionisation rate . The main results of the study are the following:

    
  * The depletion factor map is extended from the North of OMC-2 to the tip of OMC-3 with values ranging within f_D∼1-5. Such dynamic range is consistent with previous estimates across whole clouds <cit.> and within single star-forming regions <cit.>.
    
  * The o-H_2D^+ abundance inferred from f_D has a constrained dynamic range between ∼3-8×10^-11, in agreement with previous independent observations <cit.>. Although estimated indirectly, the good agreement of X(o-H_2D^+) across studies suggests f_ D to be a reliable proxy, at least at the studied spatial scales.
    
  * The deuteration fraction (R_D) we determine shows values within ∼0.005-0.1. These values are consistent with previous estimates using the same molecular species <cit.> and show a positive correlation with f_D for the volume density regime we explore (n∼10^5 cm^-3), consistently with previous works <cit.>. The deuteration fraction shows a temperature dependence, most prominent for T_K≳20 K. Below this value, its positive correlation with f_D is more robust. 
    
  * The map of  shows a South-North gradient with values varying within ∼5×10^-18-2×10^-16 s^-1. Such gradient is likely connected to the temperature dependence of R_D, however, local variations, especially for T_K∼16-19 K are driven by the degree of depletion. Overall the large dynamic range in  suggests local physical properties (e.g., N(H_2)) to influence its value, as expected from theoretical models <cit.>.
    
  * The ionisation rates determined for the column density regime probed (∼10^22-2×10^23 cm^-2) are comparable to the theoretical predictions for the same column densities <cit.>. They show a mild anti-correlation with the column density, as expected from the models of CRs propagation within clouds, extended to values of N(H_2)≳10^23 cm^-2 when complemented by the estimates of <cit.>. 
    
  * Finally, we tested the influence of the estimated path length l on . When directly determined from our low-resolution DCO^+ observations, the path length l is ∼0.18 pc. Such value, when applied to our analysis, produces results no longer consistent with those in high-mass regions <cit.> and with independent studies in the ISF <cit.>. In addition, these lower ionisation rates are only marginally reproduced by the theoretical predictions <cit.>. All these discrepancies suggest the need of high resolution data, to provide reliable estimates of l, which represents the major limiting factor in the ionisation rates determination.




This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 851435).
GS acknowledges the projects PRIN-MUR 2020 MUR BEYOND-2p (“Astrochemistry beyond the second period elements”, Prot. 2020AFB3FX) and INAF-Minigrant 2023 TRIESTE (“TRacing the chemIcal hEritage of our originS: from proTostars to planEts”; PI: G. Sabatini). SB acknowledges the ANID BASAL project FB210003.
This work is based on IRAM-30m telescope observations carried out under project numbers 032-13. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). 
This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.
This research has made use of NASA’s Astrophysics Data System.


*
aa






§ ADDITIONAL MAPS







Figure <ref> shows the OMC-1/-2/-3/-4 regions as seen through the total column density of H_2 <cit.>, through the N_2H^+ (1-0) integrated emission <cit.> and through the kinetic temperature <cit.>. The latter was directly delivered by the authors alongside with the corresponding error[<https://greenbankobservatory.org/science/gbt-surveys/gas-survey-2/>], while the other two have been computed as follows:

    
  * N(H_2): first, we determine the visual extinction in the K-band (A_K) from the dust opacity at 850 μm (τ_850):
    
    A_K = γ×τ_850 + δ,

    where γ = 2460 mag and δ = 0.012 mag for Orion A[The linear correlation between the two parameters holds across the majority of Orion A, aside for the region close to the Trapezium stars. Since we have excluded OMC-1 from the analysis, we can safely apply Eq. (<ref>) for our study <cit.>.]. Then we convert the extinction in the K-band to visual extinction (A_V) using A_K/A_V = 0.112 <cit.>. Finally, we obtain N(H_2) using the conversion <cit.>:
    
    N(H_2) = 0.94×10^21× A_V cm^-2.

    The corresponding error on the column density of H_2 per pixel is computed following the same procedure and using the error on τ_850.
    
  * I(N_2H^+): the intensity map of N_2H^+ (1-0) is obtained by integrating all the hyperfine components of the line <cit.> in the velocity range Δ V_lsr = [-5, 20] km s^-1.


Figure <ref> shows the radial sampling results onto the DCO^+ (3-2) integrated intensity map to determine the path length l (see Sect. <ref>). To perform the radial profile fitting, we employed the following approach: first, we mask the integrated intensity map using the same criteria applied in Sects. <ref>,<ref>; second, we subtract a baseline of 3×σ(DCO^+) to the profile (obtaining I_bs); third, we draw the axis manually onto the map following the peaks of emission; fourth, we sample all 4 segments with perpendicular cuts spaced by θ_beam and we sample the map along the cuts each pixel (i.e., each θ_beam); finally, we average the 4 profiles and fit the resulting radial profile with a Gaussian function. The average profile width is FWHM=0.18±0.15 pc.



§ THE DEUTERATION FRACTION SCALING RELATION





Figure <ref> shows the correlation between R_D and T_K in our study and in comparison with previous estimates. To reduce the scatter of our fields, we have sampled the median value, with corresponding IQR, of both R_D and T_K in bins of one Kelvin within 14-25 K (red circles and error bars, respectively). Then, we fitted these 11 median values with a power law dependence of the form R_D = R_D^0×(T_K/10 K)^a, where 10 K has been chosen as arbitrary normalisation temperature. The best fit to this dependency reads as follows (black line):

    R_D=0.24×(T_K/10 K)^-3.8,

suggesting variations of two orders of magnitude already within ∼10-40 K, temperature range seen across the ISF (see Fig. <ref>, panel c).


To further asses the reliability of this scaling relation, we directly compare the scaling relation determined with deuteration fractions from previous studies. To perform an homogeneous comparison, we have selected deuteration fractions determined as N(DCO^+)/N(HCO^+) and with an independent measurement of T_K. Thus, the first comparison is against the estimates of R_D and T_K in massive clumps embedded within IR Dark Clouds <cit.>. The deuteration fractions in these clumps show a similar scaling with the temperature as our estimates in Orion, up to ∼22 K. While the scaling is similar, the absolute values are, on average, a factor of ∼5 lower, which suggests an enhancement of HCO^+ deuteration in our fields. Given the similar resolution in the observations (∼20”-30”), column density (∼10^22-10^23 cm^-2) and volume density (∼0.3-2×10^5 cm^-3) across the two studies, these higher deuteration fractions are possibly connected to, overall, warmer temperatures in our fields, plus the possible presence of unresolved depletion and density gradients within the massive clumps <cit.>.


We now explore whether or not the deuteration fraction variation is connected to the star-formation regime of the region. Thus, the second comparison is against the estimates of R_D and T_K in high-mass regions in different evolutionary stages <cit.>. We have gathered the column densities of DCO^+ and HCO^+ from CDS[<https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/579/A80>], computed the median and IQR per evolutionary stage, and associated to them their average kinetic temperature determined from the chemical models of <cit.>[<cit.> also provides updated kinetic temperatures for a chemical network including deuterated species. However, their models do not improve significantly the results of <cit.> and still do not converge for some of the molecular species in the chemical network. Thus, we opted for the T_K from <cit.> given the more general set of reactions considered.]. Across the different evolutionary stages, the high-mass regions show an excellent agreement with the scaling relation determined in Eq. <ref>. The prime factor driving the decrease of R_D(HCO^+) in molecular clouds appears therefore to be T_K, while other factors, such as the star-formation history of the regions, only marginally affect its variation.



