EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 687 (July 2024) A&A, 687 (2024) A70 Full HTML
Open Access
Issue
A&A
Volume 687, July 2024
Article Number A70
Number of page(s) 14
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202449960
Published online 28 June 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Cosmic rays (CRs) are a fundamental component of the interstellar medium (ISM), but they may be the most puzzling component as well. Low-energy CRs (E < 1 TeV) can penetrate the cold (T ≲ 20 K), dense gas (n ≳ 104 cm−3) of molecular clouds and interact with H2. The chain reactions following the H2 ionisation produce H3+${\rm{H}}_3^ + $ (see Indriolo & McCall 2012, for a review), a key molecule for the production of most ionic species within molecular clouds (Herbst & Klemperer 1973). The number of ions over the number of neutral species, namely the ionisation fraction, is of paramount importance for star formation because it regulates the coupling between magnetic fields and the gas (Padovani et al. 2020), which may slow down or even halt the gravitational collapse of the cloud (Bergin & Tafalla 2007). In addition, the number of ions, H3+${\rm{H}}_3^ + $ in particular, drives the chemical evolution of the coldest regions within a cloud (e.g. Caselli & Ceccarelli 2012). However, while observed in infrared absorption in diffuse clouds (Indriolo & McCall 2012; Oka et al. 2019), H3+${\rm{H}}_3^ + $ cannot be detected in the millimeter regime, preventing a direct measure of the degree of ionisation in dense molecular clouds. The cosmic-ray ionisation rate (CRIR) is directly connected to H3+${\rm{H}}_3^ + $ and has been therefore determined over the years through several proxies, such as radicals (Black & Dalgarno 1977), ions (Caselli et al. 1998; Ceccarelli et al. 2014a; Redaelli et al. 2021b), neutrals (Fontani et al. 2017; Favre et al. 2018), and a combination of the latter (e.g. Luo et al. 2023). The values of CRIR determined in molecular clouds span three orders of magnitude (~10−17−10−14 s−1) and show a dependence on the column density of H2, N(H2), as expected from theoretical models (Padovani et al. 2009, 2022). The variety of techniques (most of them model dependent), tracers, and resolutions, however, makes for an inhomogeneous sample of scattered estimates, especially in dense clouds (N(H2) > 1022 cm−2; e.g. Bovino & Grassi 2024).

Bovino et al. (2020) recently proposed o-H2D+ (hereafter o-H2D+) as a proxy for the CRIR relative to molecular hydrogen (ζH2ion )$\left( {\zeta _{{{\rm{H}}_2}}^{{\rm{ion }}}} \right)$. H2D+ is the first product of the deuteration of H3+${\rm{H}}_3^ + $ and its ortho form was successfully detected in the millimeter regime at both low (Caselli et al. 2003, 2008; Giannetti et al. 2019; Sabatini et al. 2020; Miettinen 2020; Bovino et al. 2021) and high resolution (Redaelli et al. 2021a, 2022; Sabatini et al. 2023). Because the deuterium enrichment is favoured in cold dense gas (Millar et al. 1989; Walmsley et al. 2004), o-H2D+ is a reliable proxy of the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ in young pre-stellar environments (t ≲ 0.2 Myr; Bovino et al. 2020). This analytic derivation was recently employed by Sabatini et al. (2020) and Sabatini et al. (2023), who determined ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ from o-H2D+ 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 (Padovani et al. 2024). Whilst the estimated ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ reproduces the theoretical expectations, the detection of o-H2D+ has been a challenge for more than a decade (e.g. Caselli et al. 2008). As a consequence, the method proposed by Bovino et al. (2020) has only been applied to a limited range of column densities and scales, in particular, for N(H2) ≳ 1023 cm−2 and R ≲ 0.5 pc (Sabatini et al. 2023).

We aim to derive the column density of o-H2D+ and then the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ for column densities within N(H2) ~ 1022–1023 cm−2 and for volume densities around n ~ 105 cm−3. In the cold gas of molecular clouds where these volume densities are achieved, CO depletes onto dust grains, ions become abundant, and deutera-tion processes are boosted in general (Caselli & Ceccarelli 2012). In this context, the abundance of H2D+ is expected to grow for increasing depletion factors, as is also observed in different samples of both low- (Caselli et al. 2008) and high-mass (Sabatini et al. 2020) star-forming regions. Following previous studies, we determine N(o-H2D+) from the C18O depletion factor to provide an independent estimate of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ on parsec scales. Despite its limitations, which are discussed throughout the paper, our approach produces results that are consistent with previous works while at the same time, extending the dynamic range of column densities in the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ estimation.

The paper is then structured as follows: We introduce the method for the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ determination and its main assumptions (Sect. 2). We then present the observations, both novel (Sect. 3.1) and archival (Sect. 3.2), which are used throughout the paper, and we desscribe the related source selection. We next describe the analysis and results obtained for the degree of depletion (Sect. 4.1), the expected abundance of o-H2D+ (Sect. 4.2), the deuteration fraction (Sect. 4.3), and the ionisation rate (Sect. 4.6), along with their error budget (Sect. 4.5) and the limits of applicability (Sect. 4.4). Finally, we discuss the connection between ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ and N(H2) based on our analysis and in comparison with the theoretical predictions, together with the main limitations on the ionisation rate estimate (Sect. 5). Section 6 presents the main findings of our study.

2 Method

The method introduced by Bovino et al. (2020) provides a model-independent analytical expression for estimating ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$, which reads as follows (see also Sabatini et al. 2023):

ζH2ion =kCOo−H3+N(o−H2D+)×N(CO)3RD×N(H2)×l.$\zeta _{{{\rm{H}}_2}}^{{\rm{ion }}} = k_{{\rm{CO}}}^{{\rm{o}} - {\rm{H}}_3^ + }{{N\left( {{\rm{o}} - {{\rm{H}}_2}{{\rm{D}}^ + }} \right) \times N({\rm{CO}})} \over {3{R_{\rm{D}}} \times N\left( {{{\rm{H}}_2}} \right) \times l}}.$(1)

In this equation: kCOo−H3+$k_{{\rm{CO}}}^{{\rm{o}} - {\rm{H}}_3^ + }$ is the destruction rate of o−H3+${\rm{o}} - {\rm{H}}_3^ + $ by CO, considered here as main destruction path for o−H3+${\rm{o}} - {\rm{H}}_3^ + $ (Redaelli et al. 2024); N(CO) is the total column density of CO; RD (=N(DCO+)/N(HCO+)) is the deuteration fraction; l is the path length over which CO is depleted, and over which we estimate the column densities of our tracers; and N( o-H2D+) and N(H2) are the column densities of o-H2D+ and H2, respectively.

The robustness and reliability of Eq. (1) have recently been confirmed by Redaelli et al. (2024) through synthetic observations of the molecular tracers simulated in a three-dimensional grid. The analysis provided ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ estimates that are accurate within a factor of 2–3 compared to the input value of the simulation. These findings confirm the confidence level suggested by Bovino et al. (2020) for the method.

To apply Eq. (1), we relied on two main assumptions. First, we assumed o-H2D+ to be efficiently traced by the C18O depletion factor (ƒD). This hypothesis is based on the correlation between the o-H2D+ abundance and the degree of CO depletion, a result that was consistently derived in earlier studies through observations (e.g. Crapsi et al. 2005; Caselli et al. 2008; Sabatini et al. 2020) and through theoretical predictions (e.g. Sipilä et al. 2015; Bovino et al. 2019). We therefore proxy the abundance of o-H2D+ via ƒD using the correlation reported by Sabatini et al. (2020).

Second, we made an assumption on the path length l in Eq. (1). This path length is the physical extent of the o-H2D+ emission along the line of sight. Largely unknown, l is usually taken as the physical size of the o-H2D+ emission in the plane of the sky (e.g. Sabatini et al. 2023) or in the box where it is simulated (e.g. Redaelli et al. 2024), and it constitutes one of the main limitations of our method. We used l = 0.050 ± 0.015 pc, the typical width of filamentary substructures seen in N2H+ in the Orion Molecular Clouds OMC-2 and OMC-3 (Socci et al. 2024). The molecular ion N2H+ becomes abundant when CO freezes-out on grains (Bergin & Tafalla 2007) and its emission has been proven to correlate with that of o-H2D+ (Redaelli et al. 2022). Thus, we take it as the representative scale for our study. Although unresolved in our observations (see Sect. 3.1), this scale is the most accurate estimate in the regions without relying on direct o-H2D+ observations. In Sect. 5, we discuss how different assumptions on l affect the determination of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$.

3 Observations

3.1 IRAM-3Om observations

Applying the method from Bovino et al. (2020) to estimate ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ outside high-mass cores and clumps implies exploring column densities N(H2) ≲ 1023 cm−2 at parsec scales. To fulfil this goal, the intregral shaped filament (ISF) in Orion (Bally et al. 1987) appeared as a natural choice. The ISF is a ~13 pc long filament in the Orion A molecular cloud. It is prominent in both dust continuum (e.g. Johnstone & Bally 1999) and molecular lines (e.g. Bally et al. 1987), and it is composed by several condensations, the OMCs (see Bally 2008, for a review). The ISF is also the closest site of high-mass star-formation in the Solar neighbourhood (414 pc; Menten et al. 2007) harbouring the Orion Nebula Cluster (ONC; Hillenbrand & Hartmann 1998) within OMC-1. The extended size, the wide dynamic range in H2 column densities, and the connection to high-mass star formation make the ISF the best candidate for our study.

We surveyed OMC-1, OMC-2, OMC-3, and OMC-4 in C18O (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 3 mm. The observations were carried out at the 30-metre Insti-tut de Radioastronomie Millimétrique telescope (IRAM-30 m) in Granada (Spain) during November 2013 (Project: 032-13, PI: A. Hacar). The large-scale IRAM-30m mosaics centred on Orion BN/KL (Genzel & Stutzki 1989, α(J2000),δ(J2000) = 05h:35m:14s.20, −05º22′21″.5; see Fig. 1, panel a) were obtained with single 200″ × 200″ or 100″ × 100″ Nyquist-sampled on-the-fly (OTF) maps in position switching (PSw) mode (see Hacar et al. 2017, 2020, for additional information on the observations). Our suite of tracers have been observed with the same receiver, the Eight Mixer Receiver (EMIR; Carter et al. 2012), but with different backends: C18O(2–1) was observed with the VErsatile SPectrometer Array (VESPA) at a resolution of 20 kHz (~0.027 km s−1 at 219.560 GHz), while DCO+ (3–2) and HCO+ (1–0) were observed with the Fast Fourier Transform Spectrometer (FTS; Klein et al. 2012) at a resolution of 200 kHz (~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 (Ta*)$\left( {T_{\rm{a}}^*} \right)$ were converted into main beam temperature (Tmb) with the relation Tmb=Ta*Fef/Bef${T_{{\rm{mb}}}} = T_{\rm{a}}^*{F_{{\rm{ef}}}}/{B_{{\rm{ef}}}}$. Here, we used the standard forward (Fef) and backward (Bef) efficiencies for the IRAM-30m telescope1 (C18O (2–1): Bef = 0.61, Fef = 0.93; HCO+ (1–0): Bef = 0.81, Fef = 0.95; DCO+ (3–2): Bef = 0.62, Fef = 0.93). Finally, all the observations were convolved to a final resolution of 30″ (~0.06 pc at the Orion distance). Figure 1 shows the integrated-intensity maps of the three molecules obtained by integrating the spectral cubes in the same velocity range ΔVlsr = [6, 15] km s−1, around the typical velocity expected for the ISF (Vlsr ~ 10 km s−1; Bally 2008). Similarly, the rms noise was estimated across the whole map, integrated in spectral windows with same velocity range but devoid of emission.

C18O and HCO+ are readily detected and show widespread emission across the ISF. The high-sensitivity IRAM-30 m observations grant a peak-to-noise dynamic range of ≳100 for the intensities of both tracers. The more extended component of C18O and HCO+ emission across the ISF shows values of ~6– 10 K km s−1 on average. These values increase along the ISF crest up to peaks of ~35 K km s−1 for C18O 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 from the Trapezium stars (Hillenbrand 1997) and the ONC heats up the gas to temperatures above ~40 K (Friesen et al. 2017). While all tracers share OMC-1 as the region of peak intensity, the C18O 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 C18O emission is scattered and fragmented, while HCO+ shows a coherent filamentary shape (see Fig. 1, panels a, b); and second, in OMC-4, where the C18O emission is instead brighter than that of HCO+. The behaviour in OMC-2/OMC-3 is of interest as it indicates extended C18O depletion in these two regions.

thumbnail Fig. 1

Integral shaped filament as seen through our suite of tracers. Panel a: HCO+ (1–0) integrated-intensity map of the four OMCs composing the ISF (separated by white lines; Bally 2008). The contours correspond to [10, 30, 50, 70, 90] K km s−1. The red circle represents the position of the Orion BN/KL object. Panel b: C18O (2–1) integrated-intensity map of the OMCs composing the ISF. We include the O-type (yellow stars) and B-type (white stars) stars in the region, gathered from the SIMBAD Astronomical Database (Wenger et al. 2000), and the outflows detected by Tanabe et al. (2019) (red diamonds; see text for full discussion). The contours correspond to [10, 15, 20] K km s−1. Panel c: DCO+ (3–2) integrated-intensity map of the OMCs composing the ISF. The contours correspond to [0.5, 2, 4, 6] K km s−1, where the lowest contour corresponds to a signal-to-noise ratio of 3. The red box is the region selected for the analysis, comprising part of OMC-2 and the whole OMC-3 (see text for a discussion).

3.2 Archival observations and source selection

The ISF is a well-known region with a plethora of archival data. For our study, we took advantage of a few of these ancillary observations with a resolution comparable to that of our data. Among these, we are interested in a N(H2) map of the region to determine the candidate sites for depletion. Large-scale depletion is in fact expected for column densities N(H2) ≳ 1022 cm−2 (Bergin & Tafalla 2007; Tafalla et al. 2023). The column density map of H2 was obtained from the 850 µm dust opacity, observed by the Herschel and Planck observatories, following the prescriptions of Lombardi et al. (2014) (see also Appendix A). Figure A.1, panel a, shows the ISF as seen in column density of H2 (36.2″). The map shows material with column densities N(H2) > 1022 cm−2 extended across the whole filament with peaks of N(H2) ~ 1023 cm−2 along its crest, in correspondence to the OMCs. We focus on these regions, where the column density has a dynamic range of almost two orders of magnitude, and strong CO depletion is expected.

As mentioned in the previous section, N2H+ emission is strongly anti-correlated with CO emission because it is destroyed by the latter (e.g. Tobin et al. 2013). To further assess the presence of large-scale depletion in the ISF, we qualitatively compared the total column density of H2 with the IRAM-30 m observations of N2H+ (1–0) from Hacar et al. (2017) (Fig. A.1, panel b). Although not as extended as the material with N(H2) ~ 1022 cm−2, the N2H+ emission is widespread across the ISF and is closely correlated with N(H2) towards the OMCs (panel a).

All four regions are thus candidates for efficient depletion, and as a consequence, they are good targets for estimating ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$.

Our focus instead shifted towards OMC-2 and OMC-3 alone, as suggested above. There are two main reasons for considering these two of the four OMCs. First, the OMC-1 region hosts several O-B stars (stars in Fig. 1, panel b; Wenger et al. 2000) and a plethora of highly embedded outflows (e.g. Rivilla et al. 2013). The radiation and heating from these objects can modify the C18O abundance by its selective photodissociation or its release in the gas phase from the ice coatings of dust grains (see Draine 2011, for a more thorough discussion). For a reliable estimate of the depletion, we excluded OMC-1, where the above processes are the most efficient in the ISF. Second, DCO+ (3– 2) emits weakly across the ISF (see Sect. 3.1), and in particular, there is no emission above a signal-to-noise ratio of 3 towards OMC-4. We are therefore limited in our analysis to OMC-2 and OMC-3 (see Peterson & Megeath 2008, for a review).

OMC-2 and OMC-3 themselves show variation of more than an order of magnitude in N(H2) and temperatures between T ~ 10–40 K, according to different estimates (Friesen et al. 2017; Hacar et al. 2020). Because outflows are also detected in these two regions (diamonds in Fig. 1, panel b; Tanabe et al. 2019), we minimised their contamination to the C18O depletion estimate by imposing a temperature cut in the analysis of TK = 10–25 K (using TK(NH3) at 32″, Fig. A.1, panel c; Friesen et al. 2017). Given the close correlation between NH3 and N2H+ (e.g. Tafalla et al. 2002), we expect TK(NH3) to be the representative temperature for the gas in which CO freezes out onto the grains. The range TK = 10–25 K is consistent with previously explored ranges in studies of ƒD in high-mass regions (e.g. Sabatini et al. 2019). In addition, ~20–25 K corresponds to the temperature limit at which CO evaporates from dust grains (Bergin & Tafalla 2007) for the volume density regime expected here (i.e. n ≳ 105 cm−3).

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(H2) ~ 1023 cm−2, and temperatures mostly below ≲25 K. Throughout these ~1–2 pc north of the ISF, we expect C18O to be depleted and not photodissociated by stellar activity. Under these conditions, the depletion is expected to correlate to the o-H2D+ abundance (Caselli et al. 2008; Sabatini et al. 2020), and it may therefore be used to estimate the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$.

4 Analysis and results

4.1 Depletion map of OMC-2 and OMC-3

As we wish to determine the onset of depletion, is mandatory to properly estimate the column density of C18O (N(C18O)). First insights for this come from the inspection of I(C18O) and ƒ(N2H+) and their comparison with N(H2). We can qualitatively explore the N(H2) regime for which C18O depletes by computing the temperature-corrected line intensity (Tafalla et al. 2021),

IT=IJv(TK)−Jv(Tbg).${I_{\rm{T}}} = {I \over {{J_v}\left( {{T_{\rm{K}}}} \right) - {J_v}\left( {{T_{{\rm{bg}}}}} \right)}}.$(2)

Through Eq. (2), the column density of H2 alone now contributes to the line excitation, and therefore, to its intensity IT. Figure 2 shows the comparison between IT and N(H2) for our suite of tracers. We focus on C18O (panel a) and N2H+ (panel b). These two tracers show a continuous distribution with N(H2), but they behave differently for increasing column densities. Across the dynamic range probed by our analysis (N(H2) ~ 1022–1023 cm−2), the flattening of IT(C18O), on one hand, and the linear increase in IT(N2H+), on the other hand, represent a clear chemical differentiation of the two tracers, in particular, of the effect of depletion on C18O (e.g. Tafalla et al. 2023).

To quantitatively describe the degree of C18O depletion in OMC-2 and OMC-3, we determined its column density (N(C18O)). It can be computed as follows (Goldsmith & Langer 1999):

Ntot =8πkv2hc3AulQ(Tex)eEul/kTexguIτ1−e−τ,${N_{{\rm{tot }}}} = {{8\pi k{v^2}} \over {h{c^3}{A_{{\rm{ul}}}}}}{{Q\left( {{T_{{\rm{ex}}}}} \right){e^{{E_{{\rm{ul}}}}/k{T_{{\rm{ex}}}}}}} \over {{g_{\rm{u}}}}}{{I\tau } \over {1 - {e^{ - \tau }}}},$(3)

with the corresponding parameters h, the Planck constant, k, the Boltzmann constant, c, the light speed, v, the rest frequency of the line, Aul, the Einstein coefficient for spontaneous emission, Q, the partition function calculated at the excitation temperature of the line (Tex), gu, the degeneracy of the upper level, Eu, the excitation energy of the upper level, and τ, the optical depth of the line. We can simplify the above equation for N(C18O) by introducing two approximations. First, we assume C18O (2– 1) to be in local thermodynamic equilibrium (LTE) at the gas kinetic temperature TK. Given the presence of N2H+ in the gas phase, the volume densities are expected to be n ~ 105 cm−3, which is higher by at least an order of magnitude than the critical density of C18O (e.g. Tafalla et al. 2002). We therefore took Tex(C18O) = TK(NH3) for the rest of the analysis. Second, we work under the optically thin assumption (i.e. τ << 1), adducing the flattening in C18O intensity only to its freeze-out onto the dust grains (see Tafalla et al. 2023, for a discussion). With these two approximations, Eq. (3) can be rewritten as

Ntot(C18O)~8πkv2hc3AulQ(TK)eEul/kTKguI(C18O).${N_{{\rm{tot}}}}\left( {{{\rm{C}}^{18}}{\rm{O}}} \right)\~{{8\pi k{v^2}} \over {h{c^3}{A_{{\rm{ul}}}}}}{{Q\left( {{T_{\rm{K}}}} \right){e^{{E_{{\rm{ul}}}}/k{T_{\rm{K}}}}}} \over {{g_{\rm{u}}}}}I\left( {{{\rm{C}}^{18}}{\rm{O}}} \right).$(4)

The ratio of the expected abundance of C18O at the galacto-centric distance of Orion (i.e. X0(C18O) = 1.8 × 10−7; Giannetti et al. 2017) and the abundance measured from the C18O column density (Xobs = Ntot(C18O)/N(H2)) gives the depletion factor

fD=X0Xobs.${f_{\rm{D}}} = {{{X_0}} \over {{X_{{\rm{obs}}}}}}.$(5)

Figure 3, panel a, shows this depletion factor, which, despite our source selection (see Sect. 3.2), is extended for ~1 pc from OMC-2 to the tip of OMC-3. Depletion factors ƒD > 1 correlate extremely well with N2H+ intensities above ƒ(N2H+) > 4 K km s−1 (signal-to-noise ratio of 10, black contours) and with the column densities N(H2) > 1022 cm−2 in the two regions (see Fig. A.1). The depletion factors in our maps range from ƒD ~ 1–5, which is consistent with previous studies at comparable resolution (i.e. ~30″) in infrared dark clouds (IRDCs; e.g. Hernandez et al. 2011; Sabatini et al. 2019), single low- (e.g. Crapsi et al. 2005) and high-mass (e.g. Fontani et al. 2006) star-forming regions. Similarly to these studies, we also determine fields in our map with ƒD slightly below 1 (~0.85 on average), which would imply an over-abundance of C18O compared to the expected X0. Since these fields are restricted to the edges of the map, and given the small deviations from unity, these values ƒD < 1 could result from uncertainties in their derivation rather than from a physical over-abundance of C18O. We therefore masked these fields in Fig. 3, panel a, to consider only the regions in which depletion is effective (ƒD ≥ 1). This constrains our map to the filamentary structure that extends throughout OMC-2 and OMC-3.

thumbnail Fig. 2

Temperature-corrected integrated intensities (IT) for our suite of tracers: C18O (panel a), N2H+ (panel b), HCO+ (panel c), and DCO+ (panel d). All the IT values are computed following Tafalla et al. (2021) and are compared to N(H2) pixel by pixel (Lombardi et al. 2014).

4.2 Expected abundance of o-H2D+ on parsec scales

The challenging observations of o-H2D+ have limited its study to surveys of individual cores and clumps within low- (e.g. Caselli et al. 2008; Bovino et al. 2021) and high-mass (e.g. Sabatini et al. 2020; Redaelli et al. 2021a, 2022) star-forming regions. The opportunity to proxy its abundance with ƒD, under proper conditions, opens up the path for parsec-scale studies of the ionisation processes. We thus applied the empirical correlation between ƒD and X(o-H2D+) determined by Sabatini et al. (2020),

log10(X(o−H2D+))=0.05×fD−10.46.${\log _{10}}\left( {X\left( {{\rm{o}} - {{\rm{H}}_2}{{\rm{D}}^ + }} \right)} \right) = 0.05 \times {f_{\rm{D}}} - 10.46.$(6)

Figure 3, panel b shows the X(o-H2D+) estimated with Eq. (6). The dynamic range of X(o-H2D+) is constrained within a factor of 3 (~3–8 × 10−11) and agrees excellently with the same abundances from Sabatini et al. (2020), as expected, because the depletion factors also agree well. These abundances of o-H2D+ are broadly consistent with independent estimates for low-mass cores (Caselli et al. 2008), but they are lower by an order of magnitude on average than the values determined for high-mass regions (Redaelli et al. 2021a; Sabatini et al. 2023). The difference between our estimates and those in high-mass regions can be explained by two main factors. First, the column densities probed in these regions (N(H2) ≳ 5 × 1023 cm−2) are higher compared to our work. The second factor is the beam dilution effect when comparing our low-resolution IRAM-30 m observations (30″) to their high-resolution ALMA observations (~2″). A comparison between observations and 3D modelling of the CO freeze-out suggests that depletion factors of ƒD ~ 3–8 are typical at scales of ~0.1 pc (Bovino et al. 2019). These scales are close to the resolution of our maps (~0.06–0.07 pc), and therefore, X(o-H2D+) shows the same beam dilution effect as ƒD, from which is inferred.

The o-H2D+ abundances we derived agree with those determined in previous studies for the most part. This highlights the CO depletion factor as a viable proxy for X(o-H2D+). However, independent measurements of these two parameters are needed to ultimately break the degeneracy and evaluate the empirical correlation found by Sabatini et al. (2020). The increasing sensitivity and spatial resolution of modern astronomical facilities, such as ALMA, will allow us to achieve this goal in the near future.

thumbnail Fig. 3

Results for OMC-2 and OMC-3. Panel a: depletion map (ƒD) obtained from the ratio of the expected abundance of C18O (X0(C18O)) in Orion and the measured abundance (Ntot(C18O)/N(H2)). The black contours represent I(N2H+) > 4 K km s−1 (signal-to-noise ratio of 10; see Fig. A.1, panel b). Panel b: abundance of o-H2D+ (X(o-H2D+)) inferred from the C18O depletion factor using the correlation from Sabatini et al. (2020). Panel c: deuteration fraction, determined as RD = N(DCO+)/N(HCO+), in OMC-2 and OMC-3 (see text for its determination). Panel d: ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ determined with Eq. (1). The white star represents a B-type star within the selected region boundaries. The dotted red circle is a candidate expanding CO-shell (Feddersen et al. 2018) possibly connected to the star.

4.3 Deuteration fraction in OMC-2 and OMC-3

The deuteration fraction (RD) is directly connected to all D-bearing isotopologues of H3+${\rm{H}}_3^ + $. These molecules are the main deuterium donors in the cold and dense gas of molecular clouds (e.g. Ceccarelli et al. 2014a). We estimated the deuteration fraction as N(DCO+)/N(HCO+) through the DCO+ (3–2) and HCO+ (1–0) lines. Our first approach could again be to estimate the column densities of both tracers using Eq. (3) in the optically thin limit (τ << 1) and to assume LTE at the kinetic temperature TK. While reasonable for C18O, the LTE assumption is not ideal for DCO+ and HCO+ because of the higher critical densities of these two molecules (e.g. Sanhueza et al. 2012; Keown et al. 2016). In addition, the HCO+ (1–0) can easily become optically thick (e.g. Vasyunina et al. 2012), and therefore, a correction for τ is expected.

Similarly to C18O, we inspected the temperature-corrected intensities of HCO+ and DCO+ and compared them to N(H2) (see Fig. 2). IT(HCO+) (panel c) shows a clear dependence on N(H2) with a compact distribution. Only a handful of fields show significant scatter, which implies the almost complete absence of self-absorbed profiles in our sample (see Tafalla et al. 2023). However, while not self-absorbed, the HCO+ intensities likely become opaque as the distribution flattens towards N(H2) ~ 1023 cm−2. IT(DCO+) (panel d) also clearly depends on N(H2). Compared to those of HCO+, the intensities of DCO+ increase almost linearly with N(H2) for a distribution resembling that of N2H+. However, IT(DCO+) shows a higher degree of scatter compared to IT(N2H+), suggesting deviations from the LTE assumption used to apply Eq. (2).

To estimate the excitation effects on DCO+ and HCO+, we ran a grid of radiative transfer models for each pixel using RADEX (van der Tak et al. 2007) in two different steps2. In the first step, we input Ntot(DCO+) and Ntot(HCO+) calculated in the LTE optically thin limit and varied n, the volume density, in logarithmic intervals between 104 and 106 cm−3. In the second step, we fixed these same volume densities and let the two column densities vary by an order of magnitude around the LTE optically thin values. We then took from the output of RADEX the two Ntot that minimised the observed intensities. In both steps, we input TK as the kinetic temperature and fixed the line widths of DCO+ and HCO+ to ∆υ = 0.75 ± 0.10 km s−1 and ∆υ = 2.5 ± 0.4 km s−1, respectively. These line widths are the mean values obtained by fitting a few representative spectra in our fields, and their error was taken as the standard deviation seen in these fits.

The deuteration fraction (Fig. 3, panel c) shows variation by more than an order of magnitude (~0.005–0.1) in OMC-2 and OMC-3. These values are consistent with previous studies in high-mass star-forming regions (Sabatini et al. 2020, 2023; Pazukhin et al. 2023), that used the same molecular species, and with theoretical predictions from chemical models (Albertsson et al. 2013). The connection found between ƒD and RD in these and other studies (e.g. in low-mass regions; Crapsi et al. 2005) demonstrates the enhancement of the deuteration process in the cold and dense gas where CO is depleted (e.g. Caselli & Ceccarelli 2012, for a review). Simulations of turbulent and magnetised parsec-scale filaments, such as the ISF (e.g. Pattle et al. 2017; Hacar et al. 2018; Zielinski & Wolf 2022), show increasing levels of deuteration as the column density increases (Körtgen et al. 2018). The simulations from Körtgen et al. determined deuteration fractions within ~0.01–0.1 on timescales t ≲ 0.4 Myr, which is consistent with the expected survival time of the dense gas in Orion (~0.5 Myr; Hacar et al. 2024). Since o-H2D+ is the driver of deuteration in the early phases of the cloud evolution (Bovino et al. 2020), and given the dependence between ƒD and RD, a correlation between X(o-H2D+) and ƒD would be readily explained.

The variation in RD in our maps shows a progressive increase from OMC-2 in the south to OMC-3 towards the north. This gradient reflects the morphology of the HCO+ and DCO+ integrated intensities shown in Figs. 1a, c. The variation could be linked to a change in volume density which would cause the higher-J levels of the two molecules to be more populated. The non-LTE analysis returns a limited range in volume densities, however, which varies within ~0.3–2 × 105 cm−3 (with mean value and standard deviation of n = 1.0 ± 0.2 × 105 cm−3). Therefore, the significant variation in RD is likely the result of a more efficient deuteration process towards OMC-3.

It is well known (since Watson 1974) 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 (see Pazukhin et al. 2023, for a recent example). Thus, the gradient we observe in RD is possibly linked to the slight change in TK from OMC-3 to OMC-2 (see Fig. A.1, panel c). We explore this hypothesis in Fig. 4, where we display the comparison of RD with TK (left panel) and ƒD (right panel). The left panel shows an anti-correlation between the deuteration fraction and the kinetic temperature that was also reported in previous studies (e.g. Pazukhin et al. 2023). This anti-correlation becomes most prominent for fields with low depletion factors (ƒD ≲ 2). When the temperatures increase towards and above TK ≳ 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 (ƒD ≳ 2–3) also show a higher degree of deuteration (RD ≳ 0.03). These same points may be identified in the right panel and they show a positive correlation between RD and ƒD. This dependence between the two parameters was previously reported by Crapsi et al. (2005) for a similar range of values and volume density regime (n ~ 105 cm−3). We therefore expect the correlation between RD and ƒD to be robust below ≲20 K and to become weaker for higher temperatures with the corresponding ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ estimates possibly uncertain.

As a final test, we determined a scaling relation between RD and TK in Fig. B.1. First, we sampled the median value and interquartile range (IQR) of RD and TK in bins of 1 K (red circles and error bars, respectively). Then, we fitted a power-law dependence to these estimates, for which the best fit (black line) reads as RD~TK−4${R_{\rm{D}}}\~T_{\rm{K}}^{ - 4}$ (see Eq. (B.1)). This sharp decrease in RD(HCO+) with TK was previously observed in different clouds (e.g. Pazukhin et al. 2023), but it was never before condensed in a scaling relation. The direct comparison of this scaling relation with estimates of RD(HCO+) and TK in massive clumps (Miettinen et al. 2011) and high-mass regions (Gerner et al. 2015) confirms the goodness of our fit (see the full discussion in Appendix B). The systematic nature of the effect, which appears to be independent from the star-formation regime of the region, suggests the use of temperature-corrected deuteration fractions for future estimates of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ using the method described in Sect. 2.

thumbnail Fig. 4

Comparison between RD, ƒD and TK in our targets. Left panel: deuteration fraction (RD) compared to the kinetic temperature (TK). Each point is colour-coded with the corresponding depletion factor (ƒD). Right panel: deuteration fraction (RD) now compared to the depletion factor (ƒD) and colour-coded with the corresponding kinetic temperature (TK).

4.4 Destruction rate kCOo−H3+$k_{{\rm{CO}}}^{{\rm{o}} - {\rm{H}}_3^ + }$

Equation (1) assumes that o−H3+${\rm{o}} - {\rm{H}}_3^ + $ is primarily destroyed by CO (Bovino et al. 2020; Redaelli et al. 2024). While this is one of the main processes in dense cores, other destruction pathways can occur when the ionisation rate is studied at parsec scales. In particular, electrons can reduce the H3+${\rm{H}}_3^ + $ abundance through dissociative recombination. This process becomes important at extreme densities, when the CO is almost entirely depleted (e.g. Redaelli et al. 2024), and at lower densities, especially in the regions that are externally illuminated by stellar radiation, such as the ISF (e.g. Pabst et al. 2019).

We therefore considered estimates of the ionisation fraction, x(e−), in these two density regimes and compared them with the analytical limit ƒD < 9 × 10−7/x(e−) proposed by Redaelli et al. (2024) for Eq. (1). First, we considered the estimates in massive dense cores, including some in Orion (Bergin et al. 1999). Berginetal. estimated ionisation fractions x(e−) ≲ 10−7 at N(H2) ≳ 1022 cm−2 for these sources, leading to a depletion limit for the method of ƒD ≳ 9. Second, we considered the estimates in OMC-2 and OMC-3 using the C2H and HCN lines (Salas et al. 2021). Salas et al. instead determined ionisation fractions x(e−) ≤ 3 × 10−6 at n ~ 5 × 103 cm−3, leading to an unphysical ƒD > 0.3. Given our range in both column and volume density (N(H2) ~ 1022–1023 cm−2 and n(H2) ~ 0.3–2 × 105 cm−3, respectively), we can safely assume that the destruction through CO is the main process affecting the H3+${\rm{H}}_3^ + $ abundance for our degree of depletion.

We finally determined the destruction rate kCOo−H3+$k_{{\rm{CO}}}^{{\rm{o}} - {\rm{H}}_3^ + }$ in our regions. We computed the latter pixel by pixel using TK and the coefficients found in the KInetic Database for Astrochem-istry (KIDA3; Wakelam et al. 2012) for the ion-polar reaction H3++CO→HCO++H2${\rm{H}}_3^ + + {\rm{CO}} \to {\rm{HC}}{{\rm{O}}^ + } + {{\rm{H}}_2}$. The adopted kCOo−H3+$k_{{\rm{CO}}}^{{\rm{o}} - {\rm{H}}_3^ + }$ 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.

4.5 Error budget in the determination of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$

Our indirect determination of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ presents a number of error sources that are discussed in the following. Excluding X(o-H2D+), for which the CO depletion is assumed as an efficient proxy (see Sect. 2), all the other parameters in Eq. (1) bear a corresponding uncertainty, as we list below.

  • kCOo−H3+$k_{{\rm{CO}}}^{{\rm{o}} - {\rm{H}}_3^ + }$, N(C18O) both only depend on the temperature (see Eq. (4) and Sect. 4.4). We therefore propagated their error accordingly using the uncertainty on TK (Friesen et al. 2017).

  • [16O]/[18O] was used to convert N(C18O) into N(CO) in Eq. (1). This parameter is equal to 560 ± 26 for the local ISM (see Table 4 in Wilson & Rood 1994).

  • RD has, given its definition and its derivation with RADEX, the combined uncertainties on the DCO+ and HCO+ line widths as the total uncertainty (see Sect. 4.3);

  • l has the typical size of N2H+ structures in OMC-2 and OMC-3 (see Sect. 2), along with their corresponding error (see also Socci et al. 2024).

The error per field in our analysis reaches up to ~50% of the corresponding ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ in OMC-2 and OMC-3. However, this error contains a strong flat contribution from the conservative choices we made to estimate the uncertainty on some parameters (e.g. line widths of DCO+ and HCO+). A more thorough analysis would therefore significantly reduce the uncertainty in several of our fields. Since the determination of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ remains indirect and this error budget does not hinder our further discussion and conclusions, we continued with these choices for our analysis.

4.6 Cosmic-ray ionisation rate map of Orion

The map of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ (Fig. 3, panel d) shows a dynamic range of more than one order of magnitude in the OMC-2/OMC-3 region. The global south-north decrease in ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ is likely connected to the gradient seen in the deuteration fraction, which is further connected to the effect of temperature (see Fig. 4, left panel). However, the local variations in ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ are likely induced by the change in fD throughout the filament (see Fig. 4, right panel). As a result of these two effects, ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ varies within ~5 × 10−18–2 × 10−16 s−1 within the two regions.

The median value we derive of ζH2ion ~2.4×10−17 s−1$\zeta _{{{\rm{H}}_2}}^{{\rm{ion }}}\~2.4 \times {10^{ - 17}}{{\rm{s}}^{ - 1}}$ agrees well with the estimates in high-mass star-forming regions using the same method (Sabatini et al. 2020, 2023). The slight differences reflect the variation in X(o-H2D+) (see Sect. 4.2), which may be ascribed to a resolution effect. The highest values of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ we determine (i.e. ≳ 2 × 10−16 s−1) are seen south of the map, and in particular, in the region close to the OMC-2 FIR4 protocluster (e.g. Johnstone et al. 2003) and the B-type star HD 37060 (see Fig. 3, panel d). These two heating sources can indirectly affect our estimate through an increase in TK, with a corresponding decrease in RD. In addition, the outflows from the protoclus-ter (e.g. López-Sepulcre et al. 2013; Chahine et al. 2022) and the UV-irradiation from the star could promote the photodisso-ciation of C18O, inducing an increase in ƒD (see Fig. 4, right panel). However, higher values of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ could also result from a local enhancement in the CR flux (see Padovani et al. 2015,2016, 2019) caused by these objects, as suggested by previous findings: CRs ionisation rates as high as ζH2ion ~10−14 s−1$\zeta _{{{\rm{H}}_2}}^{{\rm{ion }}}\~{10^{ - 14}}{{\rm{s}}^{ - 1}}$ were previously reported in OMC-2 FIR4 (Ceccarelli et al. 2014b; Fontani et al. 2017). A candidate expanding CO-shell was detected in close correspondence to the B-type star (Feddersen et al. 2018), for which CRs could travel beyond its edge and promote the ionisation in the southern part of our filament.

Despite all the influences on the ionisation rate, ζH2ion ≈10−17 s−1$\zeta _{{{\rm{H}}_2}}^{{\rm{ion }}} \approx {10^{ - 17}}{{\rm{s}}^{ - 1}}$, originally found by Spitzer & Tomasko (1968), was instead considered for several decades as the reference value for the ISM. While this value is reproduced by our results, our parsec-scale map shows variation of at least an order of magnitude for ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ in 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. Of the possible physical influences on the CR ionisation rate, we consider the effect of increasing N(H2) as the prime driver of its variation.

5 Discussion

The different models for the energy loss of CR particles (protons, electrons and heavy nuclei) with increasing N(H2) were originally studied by Padovani et al. (2009). These models account for the ionisation by primary CR protons and electrons, plus secondary electrons from primary CRs. For the column density regime we probed, the models of CR propagation (see Fig. 5, coloured lines) are largely dominated by CR protons and their secondaries. As a general trend, all models predict a decrease in ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ with increasing N(H2). Different slopes come from the comparison between the CR propagation models and observational results: model ℒ, from the CR spectrum observed by the Voyager missions (Cummings et al. 2016; Stone et al. 2019); model ℋ, from the average ionisation rate estimated in diffuse clouds (Indriolo & McCall 2012; Neufeld & Wolfire 2017); and model 𝒰, from the upper limit of the ionisation rate in diffuse clouds. Observational results show a decreasing trend as well, but the large number of methods and tracers used to infer ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ produces a significant scatter in the data, especially for dense clouds (see Fig. B.1 in Padovani et al. 2024).

Our new results considerably increase the statistics for the available ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ estimates compared to previous studies (e.g. Caselli et al. 1998; van der Tak & van Dishoeck 2000). In particular, we probed the CR ionisation rate within a poorly explored column density regime. Figure 5 shows the comparison between ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ and N(H2), colour-coded by Tĸ, with their corresponding errors (see Lombardi et al. 2014, and Sect. 4.5, respectively). In addition, we plot the models from Padovani et al. (2022) for different slopes of the CR proton energy spectrum (coloured lines). Almost the entire sample of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ determined in our study is enclosed within models ℒ and ℋ, with only a fraction of fields outside of this range. Despite their departure from the bulk of our points, these ionisation rates above model ℋ and below model ℒ are still consistent with the theoretical predictions within the errors. None of the predicted propagation models clearly reproduces our results. Our estimates show the functional dependence between ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ and N(H2) expected from these models, however. While the temperature dependence of RD (see Sect. 4.3) contributes to the scatter, the fields with Tĸ ≲ 20 K show decreasing ionisation rates for increasing N(H2).

As we included the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ determined in two high-mass star-forming regions using the same method (diamonds and squares; Sabatini et al. 2023), the ionisation rates in the two studies agree well. Our estimates show values closer to those in AG531 overall, while only the highest ionisation rates we determined are comparable to those in AG354. The combined distribution of ionisation rates from the two studies is continuous with N(H2) and mildly depends on it within ~1022–5 × 1023 cm−2. The method proposed by Bovino et al. (2020) therefore allows homogeneous estimates of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ in a wide range of scales and column density regimes, when o-H2D+ is either successfully detected or reasonably inferred.

The two estimates from our paper and Sabatini et al. (2023) agree well, but spurious effects have to be considered in their determination. While resolution affects the degree of depletion on the one hand (see Sect. 4.2), it also influences the minimum physical size we can probe with our observations on the other hand. This physical size is condensed in the path length l from Eq. (1). As mentioned in Sect. 2, it is complex to measure either the line-of-sight length or the volume density of the tracers, and therefore, l is usually taken as the radial size of the structures identified in o-H2D+ (Sabatini et al. 2020, 2023). We explored the influence of l on the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ by measuring the typical radial size of the DCO+ emission, our least extended tracer. To do this, we performed 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 Appendix A for a full discussion). The FWHM of the profile (0.18 ± 0.15 pc) was taken as representative l to compute the estimate of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$. Since this FWHM is three to four times the l adopted in our analysis (0.05 pc; see Sect. 2), the corresponding ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ values are lower by the same factor on average (see Eq. (1), ζH2ion ∝1/l$\zeta _{{{\rm{H}}_2}}^{{\rm{ion }}} \propto 1/l$). The median ionisation rate with l ~ 0.18 pc is ζH2ion ~7×10−18 s−1$\zeta _{{{\rm{H}}_2}}^{{\rm{ion }}}\~7 \times {10^{ - 18}}{{\rm{s}}^{ - 1}}$, and the majority of the fields is below 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 (Pabst et al. 2019), and they strongly disagree with independent estimates, which indeed suggest higher values (e.g. Fontani et al. 2017; Salas et al. 2021). As ƒD and RD are positively correlated, their ratio is only mildly affected by resolution effects. The path length, l, on the other hand, is directly affected by a poor resolution, and thus represents the main limitation in the determination of ionisation rates with Eq. (1) (see also Bovino et al. 2020; Redaelli et al. 2024, for different ways to estimate l and a detailed analysis on its effect to the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$).

Our results highlight the CO depletion as available proxy for the o-H2D+ abundance at parsec scales. However, caution should be taken when interpreting these results: First, the determination of ƒD can be hard and the use of multiple CO lines is usually advised (e.g. Sabatini et al. 2022). Second, the connection between ƒD and the o-H2D+ abundance was derived at sub-parsec scales (Sabatini et al. 2020) and works within a constrained range of TK, outside which multiple processes may become temperature dependent (e.g. the degree of deuteration). In turn, estimates based on o-H2D+ are confirmed as an effective method for the homogeneous sampling of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ over more than an order of magnitude in column density. Ultimately, in the absence of local sources of CRs, N(H2) alone appears to drive the ionisa-tion degree in molecular clouds (Padovani et al. 2022). Finally, our results also highlight the importance of high-resolution observations to estimate the degree of depletion and to estimate l reliably. Both parameters can in fact significantly affect the determination of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$. Independent interferometric observations of both CO isotopologues and o-H2D+ are needed to validate the previous correlation by Sabatini et al. (2020) and to provide a more robust estimate of the cosmic-ray ionisation rate.

thumbnail Fig. 5

Comparison between ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ and the column density N(H2). The theoretical predictions for the CR propagation in molecular clouds are displayed as colour-coded lines (blue, green, and red lines; Padovani et al. 2022), while the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ estimates in high-mass regions are displayed in grey (squares and diamonds; Sabatini et al. 2023). The estimates of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ from our study are displayed as circles, colour-coded by Tĸ.

6 Conclusions

We presented the first parsec-scale estimate of the cosmic-ray ionisation rate determined by employing the most recent analytical framework proposed by Bovino et al. (2020). To this end, we presented novel IRAM-30 m observations of the ISF in Orion using the following suite of tracers: C18O (2– 1), DCO+ (3–2) and HCO+ (1–0). We first determined the C18O depletion factor under density and temperature regimes for which we expect significant freeze-out and no photodissociation effects. The C18O depletion factor was then used as a proxy of the o-H2D+ abundance, following the correlation determined by Sabatini et al. (2020) and to ultimately determine the ionisation rate ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$. The main results of the study are listed below:

  • The depletion factor map is extended from north of OMC-2 to the tip of OMC-3 with values ranging within ƒD ~ 1–5. This dynamic range is consistent with previous estimates for whole clouds (Hernandez et al. 2011; Sabatini et al. 2019; Feng et al. 2020) and within single star-forming regions (Crapsi et al. 2005; Fontani et al. 2006; Sabatini et al. 2022, 2023).

  • The o-H2D+ abundance inferred from ƒD has a constrained dynamic range between ~3–8 × 10−11, in agreement with previous independent observations (Caselli et al. 2008; Giannetti et al. 2019; Miettinen 2020; Redaelli et al. 2021a; Sabatini et al. 2023). Although estimated indirectly, the good agreement of our X(o-H2D+) compared to those measured in these studies suggests that ƒD is a reliable proxy, at least at the spatial scales we studied.

  • The deuteration fraction (RD) we determined shows values within ~0.005–0.1. These values are consistent with previous estimates using the same molecular species (Sabatini et al. 2023) and show a positive correlation with ƒD for the volume density regime we explored (n ~ 105 cm−3). A positive correlation between ƒD and RD was previously reported for similar volume densities (Crapsi et al. 2005). The deuter-ation fraction decreases with temperature, in particular for TK ≳ 20 K. Below this value, its positive correlation with ƒD is more robust.

  • The map of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ shows a south-north gradient with values that vary within ~5 × 10−18–2 × 10−16 s−1. This gradient is likely connected to the temperature dependence of RD, but local variations, especially for TK ~ 16–19 K are driven by the degree of depletion. The wide dynamic range in ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ overall is ascribed to the variation of local physical properties (e.g. N(H2)), as expected from theoretical models (Padovani et al. 2009).

  • The ionisation rates we determined for the column density regime we probed (~1022–2 × 1023 cm−2) are comparable to the theoretical predictions for the same column densities (Padovani et al. 2022). They are mildly anti-correlated with the column density, as expected from the models of CR propagation within clouds, extended to values of N(H2) ≳ 1023 cm−2 when complemented by the estimates of Sabatini et al. (2023).

  • Finally, we tested the influence of the estimated path length l on ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$. When the path length is directly determined from our low-resolution DCO+ observations, its value is ~0.18 pc. This l produces results no longer consistent with those in high-mass regions (Sabatini et al. 2023) and with independent studies in the ISF (e.g. Fontani et al. 2017) when applied to our analysis. In addition, these lower ionisation rates are only marginally reproduced by the theoretical predictions (Padovani et al. 2022). All these discrepancies suggest that high-resolution data are needed to provide reliable estimates of l. This path length becomes otherwise the main limiting factor in the determination of the CR ionisation rates.

Acknowledgements

The authors thank the anonymous referee for their suggestions to improve the manuscript. 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 proTo-stars to planEts”; PI: G. Sabatini). SB acknowledges the ANID BASAL project FB210003. This work is based on IRAM-30 m 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.

Appendix A Additional maps

thumbnail Fig. A.1

ISF as seen in different types of emission. Panel (a): Column density map of H2 (N(H2)) obtained from dust opacity Herschel+Planck observations at 850 µm (Lombardi et al. 2014). The white contours correspond to [1022, 5 × 1022, 1023] cm−2, and the blue contour represents the footprint of our IRAM-30m observations. Panel (b): N2H+ (1–0) integrated-intensity map of observed with IRAM-30m at a resolution of 30″ (see Hacar et al. 2017, for additional information). Panel (c): Kinetic temperature of the gas (TK) obtained from observations of the two inversion transitions NH3 (1,1) and NH3 (2,2) at 32″ with the GBT (Friesen et al. 2017). The white and yellow stars represent the B- and O-type stars in the region, respectively, as in Fig. 1, panel (b).

Figure A.1 shows the OMC-1/-2/-3/-4 regions as seen through the total column density of H2 (panel a; Lombardi et al. 2014) through the N2H+ (1–0) integrated emission (panel b; Hacar et al. 2017) and through the kinetic temperature (panel c; Friesen et al. 2017). The latter was directly delivered by the authors alongside with the corresponding error4, while the other two were computed as follows:

  • N(H2). We first determined the visual extinction in the K band (AK) from the dust opacity at 850 µm (τ850),

    AK=γ×τ850+δ, ${A_{\rm{K}}} = \gamma \times {\tau _{850}} + \delta {\rm{, }}$(A.1)

    where γ = 2460 mag and δ = 0.012 mag for Orion A5. Then we converted the extinction in the K band into visual extinction (AV) using AK/AV = 0.112 (Rieke & Lebofsky 1985). Finally, we obtained N(H2) using the conversion (Bohlin et al. 1978)

    N(H2)=0.94×1021×AVcm−2.$N\left( {{{\rm{H}}_2}} \right) = 0.94 \times {10^{21}} \times {A_{\rm{V}}}{\rm{c}}{{\rm{m}}^{ - 2}}.$(A.2)

    The corresponding error on the column density of H2 per pixel was computed following the same procedure and using the error on τ850.

  • ƒ(N2H+). The intensity map of N2H+ (1–0) was obtained by integrating all the hyperfine components of the line (e.g. Caselli et al. 1995) in the velocity range ΔVlsr = [−5, 20] km s−1.

thumbnail Fig. A.2

Radial profile sampling of the DCO+ integrated intensity in OMC-2 and OMC-3. Left panel: Integrated-intensity map of DCO+ (3–2) (see also Fig. 1), masked based on the criteria discussed in Sects. 3.2,4.1. The dashed red line represents the axis used for the radial sampling. Right panel: Average intensity (Ibs) and temperature (TK) radial profiles computed across the red line in the left panel. The Gaussian best fit to the profile is also displayed (solid black line), along with the derived FWHM (solid arrow).

Figure A.2 shows the radial sampling results onto the DCO+ (3–2) integrated-intensity map to determine the path length l (see Sect. 5). To perform the radial profile fitting, we employed the following approach: We first masked the integrated-intensity map using the same criteria as applied in Sects. 3.2,4.1. We then subtracted a baseline of 3 × σ(DCO+) from the profile (obtaining Ibs). Third, we drew the axis manually onto the map following the emission peaks. Fourth, we sampled all four segments with perpendicular cuts spaced by θbeam, and we sampled the map along each cut every pixel (i.e. every θbeam/2). We finally averaged the four profiles and fit the resulting radial profile with a Gaussian function. The average profile width is FWHM = 0.18 ± 0.15 pc.

Appendix B Scaling relation of the deuteration fraction

Figure B.1 shows the correlation between RD and TK in our study and in comparison with previous estimates. To reduce the scatter of our fields, we sampled the median value, with corresponding IQR, of both RD and TK 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 RD=RD0×(TK/10 K)a${R_{\rm{D}}} = R_{\rm{D}}^0 \times {\left( {{T_{\rm{K}}}/10{\rm{K}}} \right)^a}$, where 10 K was chosen as arbitrary normalisation temperature. The best fit to this dependence reads as follows (black line):

RD=0.24×(TK10 K)−3.8,${R_{\rm{D}}} = 0.24 \times {\left( {{{{T_{\rm{K}}}} \over {10{\rm{K}}}}} \right)^{ - 3.8}},$(B.1)

suggesting variations of two orders of magnitude already within ~ 10–40 K, which is the temperature range seen across the ISF (see Fig. A.1, panel c).

thumbnail Fig. B.1

Estimates of RD = N(DCO+)/N(HCO+) and TK across different studies: estimates in Orion, both single fields and median per bin (green and red circles, respectively; this study), in IRDCs (black diamonds; Miettinen et al. 2011), in high-mass star-forming regions (colour-coded squares; Gerner et al. 2014, 2015, see text for discussion on TK). The black line represents the best fit to the median values of our estimates in Orion and it is exemplified by Eq. B.1.

To further asses the reliability of this scaling relation, we directly compared the scaling relation determined with deuter-ation fractions from previous studies. To perform an homogeneous comparison, we selected deuteration fractions determined as N(DCO+)/N(HCO+) and with an independent measurement of TK. Thus, the first comparison is against the estimates of RD and TK in massive clumps embedded within IRDCs (black diamonds in Fig. B.1; Miettinen et al. 2011, Tables 5-9). The deuteration fractions in these clumps shows a similar scaling with the temperature as our estimates in Orion, up to ~22 K. While the scaling is similar, the absolute values are lower by a factor of ~5 on average, which suggests an enhancement of HCO+ deuteration in our fields. Given the similar resolution in the observations (~20″ – 30″), column density (~1022 – 1023 cm−2) and volume density (~0.3 – 2 × 105 cm−3) in the two studies, these higher deuteration fractions might be connected to overall warmer temperatures in our fields, and maybe to unresolved depletion and density gradients within the massive clumps (Miettinen et al. 2011).

We also explored whether the variation in deuteration fraction is connected to the star-formation regime of the region. Thus, the second comparison was made with the estimates of RD and TK in high-mass regions in different evolutionary stages (i.e. High-Mass Protostellar Objects (HMPO), Cores (HMC) and Ultra-Compact HII regions (UCHII), colour-coded squares in Fig. B.1; Gerner et al. 2015). We gathered the column densities of DCO+ and HCO+ from CDS6, computed the median and IQR per evolutionary stage, and associated their average kinetic temperature determined from the chemical models of Gerner et al. (2014)7. Despite their different evolutionary stages, the high-mass regions all agree very well with the scaling relation determined in Eq. B.1. The prime factor driving the decrease in RD(HCO+) in molecular clouds therefore appears to be TK, while other factors, such as the star-formation history of the regions, affect its variation only marginally.

References

  1. Albertsson, T., Semenov, D. A., Vasyunin, A. I., Henning, T., & Herbst, E. 2013, ApJS, 207, 27 [CrossRef] [Google Scholar]
  2. Bally, J. 2008, in Handbook of Star Forming Regions, Volume I, 4, ed. B. Reipurth, 459 [NASA ADS] [Google Scholar]
  3. Bally, J., Langer, W. D., Stark, A. A., & Wilson, R. W. 1987, ApJ, 312, L45 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339 [Google Scholar]
  5. Bergin, E. A., Plume, R., Williams, J. P., & Myers, P. C. 1999, ApJ, 512, 724 [NASA ADS] [CrossRef] [Google Scholar]
  6. Black, J. H., & Dalgarno, A. 1977, ApJS, 34, 405 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132 [Google Scholar]
  8. Bovino, S., & Grassi, T. 2024, Astrochemical Modeling: Practical Aspects of Microphysics in Numerical Simulations (Elsevier) [Google Scholar]
  9. Bovino, S., Ferrada-Chamorro, S., Lupi, A., et al. 2019, ApJ, 887, 224 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bovino, S., Ferrada-Chamorro, S., Lupi, A., Schleicher, D. R. G., & Caselli, P. 2020, MNRAS, 495, L7 [Google Scholar]
  11. Bovino, S., Lupi, A., Giannetti, A., et al. 2021, A&A, 654, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Carter, M., Lazareff, B., Maier, D., et al. 2012, A&A, 538, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Caselli, P., & Ceccarelli, C. 2012, A&A Rev., 20, 56 [NASA ADS] [CrossRef] [Google Scholar]
  14. Caselli, P., Myers, P. C., & Thaddeus, P. 1995, ApJ, 455, L77 [Google Scholar]
  15. Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234 [Google Scholar]
  16. Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, L37 [CrossRef] [EDP Sciences] [Google Scholar]
  17. Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014a, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 859 [Google Scholar]
  19. Ceccarelli, C., Dominik, C., López-Sepulcre, A., et al. 2014b, ApJ, 790, L1 [Google Scholar]
  20. Chahine, L., López-Sepulcre, A., Podio, L., et al. 2022, A&A, 667, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Crapsi, A., Caselli, P., Walmsley, C. M., et al. 2005, ApJ, 619, 379 [Google Scholar]
  22. Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18 [CrossRef] [Google Scholar]
  23. Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium [Google Scholar]
  24. Favre, C., Ceccarelli, C., López-Sepulcre, A., et al. 2018, ApJ, 859, 136 [Google Scholar]
  25. Feddersen, J. R., Arce, H. G., Kong, S., et al. 2018, ApJ, 862, 121 [NASA ADS] [CrossRef] [Google Scholar]
  26. Feng, S., Li, D., Caselli, P., et al. 2020, ApJ, 901, 145 [Google Scholar]
  27. Fontani, F., Caselli, P., Crapsi, A., et al. 2006, A&A, 460, 709 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Fontani, F., Ceccarelli, C., Favre, C., et al. 2017, A&A, 605, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Friesen, R. K., Pineda, J. E., co-PIs, et al. 2017, ApJ, 843, 63 [NASA ADS] [CrossRef] [Google Scholar]
  30. Genzel, R., & Stutzki, J. 1989, ARA&A, 27, 41 [NASA ADS] [CrossRef] [Google Scholar]
  31. Gerner, T., Beuther, H., Semenov, D., et al. 2014, A&A, 563, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Gerner, T., Shirley, Y. L., Beuther, H., et al. 2015, A&A, 579, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Giannetti, A., Leurini, S., König, C., et al. 2017, A&A, 606, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Giannetti, A., Bovino, S., Caselli, P., et al. 2019, A&A, 621, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209 [Google Scholar]
  36. Hacar, A., Alves, J., Tafalla, M., & Goicoechea, J. R. 2017, A&A, 602, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Hacar, A., Tafalla, M., Forbrich, J., et al. 2018, A&A, 610, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Hacar, A., Bosman, A. D., & van Dishoeck, E. F. 2020, A&A, 635, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Hacar, A., Socci, A., Bonanomi, F., et al. 2024, A&A, in press https://doi.org/10.1051/0004-6361/202348565 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Herbst, E., & Klemperer, W. 1973, ApJ, 185, 505 [Google Scholar]
  41. Hernandez, A. K., Tan, J. C., Caselli, P., et al. 2011, ApJ, 738, 11 [NASA ADS] [CrossRef] [Google Scholar]
  42. Hillenbrand, L. A. 1997, AJ, 113, 1733 [Google Scholar]
  43. Hillenbrand, L. A., & Hartmann, L. W. 1998, ApJ, 492, 540 [NASA ADS] [CrossRef] [Google Scholar]
  44. Indriolo, N., & McCall, B. J. 2012, ApJ, 745, 91 [NASA ADS] [CrossRef] [Google Scholar]
  45. Johnstone, D., & Bally, J. 1999, ApJ, 510, L49 [Google Scholar]
  46. Johnstone, D., Boonman, A. M. S., & van Dishoeck, E. F. 2003, A&A, 412, 157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Keown, J., Schnee, S., Bourke, T. L., et al. 2016, ApJ, 833, 97 [Google Scholar]
  48. Klein, B., Hochgürtel, S., Krämer, I., et al. 2012, A&A, 542, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Körtgen, B., Bovino, S., Schleicher, D. R. G., et al. 2018, MNRAS, 478, 95 [Google Scholar]
  50. Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, A&A, 566, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. López-Sepulcre, A., Taquet, V., Sánchez-Monge, Á., et al. 2013, A&A, 556, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Luo, G., Zhang, Z.-Y., Bisbas, T. G., et al. 2023, ApJ, 946, 91 [NASA ADS] [CrossRef] [Google Scholar]
  53. Menten, K. M., Reid, M. J., Forbrich, J., & Brunthaler, A. 2007, A&A, 474, 515 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Miettinen, O. 2020, A&A, 634, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Miettinen, O., Hennemann, M., & Linz, H. 2011, A&A, 534, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Millar, T. J., Bennett, A., & Herbst, E. 1989, ApJ, 340, 906 [NASA ADS] [CrossRef] [Google Scholar]
  57. Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163 [Google Scholar]
  58. Oka, T., Geballe, T. R., Goto, M., et al. 2019, ApJ, 883, 54 [Google Scholar]
  59. Pabst, C., Higgins, R., Goicoechea, J. R., et al. 2019, Nature, 565, 618 [NASA ADS] [CrossRef] [Google Scholar]
  60. Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Padovani, M., Hennebelle, P., Marcowith, A., & Ferrière, K. 2015, A&A, 582, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2016, A&A, 590, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Padovani, M., Marcowith, A., Sánchez-Monge, Á., Meng, F., & Schilke, P. 2019, A&A, 630, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Padovani, M., Ivlev, A. V., Galli, D., et al. 2020, Space Sci. Rev., 216, 29 [CrossRef] [Google Scholar]
  65. Padovani, M., Bialy, S., Galli, D., et al. 2022, A&A, 658, A189 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Padovani, M., Galli, D., Scarlett, L. H., et al. 2024, A&A, 682, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Pattle, K., Ward-Thompson, D., Berry, D., et al. 2017, ApJ, 846, 122 [Google Scholar]
  68. Pazukhin, A. G., Zinchenko, I. I., Trofimova, E. A., Henkel, C., & Semenov, D. A. 2023, MNRAS, 526, 3673 [NASA ADS] [CrossRef] [Google Scholar]
  69. Peterson, D. E., & Megeath, S. T. 2008, in Handbook of Star Forming Regions, Volume I, 4, ed. B. Reipurth, 590 [NASA ADS] [Google Scholar]
  70. Redaelli, E., Bovino, S., Giannetti, A., et al. 2021a, A&A, 650, A202 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Redaelli, E., Sipilä, O., Padovani, M., et al. 2021b, A&A, 656, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Redaelli, E., Bovino, S., Sanhueza, P., et al. 2022, ApJ, 936, 169 [NASA ADS] [CrossRef] [Google Scholar]
  73. Redaelli, E., Bovino, S., Lupi, A., et al. 2024, A&A, 685, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Rieke, G. H., & Lebofsky, M. J. 1985, ApJ, 288, 618 [Google Scholar]
  75. Rivilla, V. M., Martín-Pintado, J., Sanz-Forcada, J., Jiménez-Serra, I., & Rodríguez-Franco, A. 2013, MNRAS, 434, 2313 [NASA ADS] [CrossRef] [Google Scholar]
  76. Sabatini, G., Giannetti, A., Bovino, S., et al. 2019, MNRAS, 490, 4489 [Google Scholar]
  77. Sabatini, G., Bovino, S., Giannetti, A., et al. 2020, A&A, 644, A34 [EDP Sciences] [Google Scholar]
  78. Sabatini, G., Bovino, S., Sanhueza, P., et al. 2022, ApJ, 936, 80 [NASA ADS] [CrossRef] [Google Scholar]
  79. Sabatini, G., Bovino, S., & Redaelli, E. 2023, ApJ, 947, L18 [NASA ADS] [CrossRef] [Google Scholar]
  80. Salas, P., Rugel, M. R., Emig, K. L., et al. 2021, A&A, 653, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Sanhueza, P., Jackson, J. M., Foster, J. B., et al. 2012, ApJ, 756, 60 [Google Scholar]
  82. Sipilä, O., Caselli, P., & Harju, J. 2015, A&A, 578, A55 [Google Scholar]
  83. Socci, A., Hacar, A., & Bonanomi, F., et al. 2024, A&A, submitted [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Spitzer, Lyman, J., & Tomasko, M. G. 1968, ApJ, 152, 971 [Google Scholar]
  85. Stone, E. C., Cummings, A. C., Heikkila, B. C., & Lal, N. 2019, Nat. Astron., 3, 1013 [NASA ADS] [CrossRef] [Google Scholar]
  86. Tafalla, M., Myers, P. C., Caselli, P., Walmsley, C. M., & Comito, C. 2002, ApJ, 569, 815 [Google Scholar]
  87. Tafalla, M., Usero, A., & Hacar, A. 2021, A&A, 646, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Tafalla, M., Usero, A., & Hacar, A. 2023, A&A, 679, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Tanabe, Y., Nakamura, F., Tsukagoshi, T., et al. 2019, PASJ, 71, S8 [NASA ADS] [Google Scholar]
  90. Tobin, J. J., Bergin, E. A., Hartmann, L., et al. 2013, ApJ, 765, 18 [NASA ADS] [CrossRef] [Google Scholar]
  91. van der Tak, F. F. S., & van Dishoeck, E. F. 2000, A&A, 358, L79 [NASA ADS] [Google Scholar]
  92. van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Vasyunina, T., Vasyunin, A. I., Herbst, E., & Linz, H. 2012, ApJ, 751, 105 [NASA ADS] [CrossRef] [Google Scholar]
  94. Wakelam, V., Herbst, E., Loison, J. C., et al. 2012, ApJS, 199, 21 [Google Scholar]
  95. Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, 418, 1035 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Watson, W. D. 1974, ApJ, 188, 35 [NASA ADS] [CrossRef] [Google Scholar]
  97. Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191 [Google Scholar]
  99. Zielinski, N., & Wolf, S. 2022, A&A, 659, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
2

The radiative transfer calculations were not applied to C18O because without proper modelling of the cloud, which includes a freeze-out threshold, depletion may be mistaken for opacity (e.g. see Tafalla et al. 2021).

5

The linear correlation between the two parameters holds for most of Orion A, except for the region close to the Trapezium stars. Since we excluded OMC-1 from the analysis, we can safely apply Eq. (A.1) for our study (see Lombardi et al. 2014, for a full discussion).

7

Gerner et al. (2015) also provided updated kinetic temperatures for a chemical network including deuterated species. However, their models did not significantly improve the results of Gerner et al. (2014) and still did not converge for some of the molecular species in the chemical network. Thus, we opted for the TK from Gerner et al. (2014) given the more general set of reactions considered.

All Figures

thumbnail Fig. 1

Integral shaped filament as seen through our suite of tracers. Panel a: HCO+ (1–0) integrated-intensity map of the four OMCs composing the ISF (separated by white lines; Bally 2008). The contours correspond to [10, 30, 50, 70, 90] K km s−1. The red circle represents the position of the Orion BN/KL object. Panel b: C18O (2–1) integrated-intensity map of the OMCs composing the ISF. We include the O-type (yellow stars) and B-type (white stars) stars in the region, gathered from the SIMBAD Astronomical Database (Wenger et al. 2000), and the outflows detected by Tanabe et al. (2019) (red diamonds; see text for full discussion). The contours correspond to [10, 15, 20] K km s−1. Panel c: DCO+ (3–2) integrated-intensity map of the OMCs composing the ISF. The contours correspond to [0.5, 2, 4, 6] K km s−1, where the lowest contour corresponds to a signal-to-noise ratio of 3. The red box is the region selected for the analysis, comprising part of OMC-2 and the whole OMC-3 (see text for a discussion).
In the text
thumbnail Fig. 2

Temperature-corrected integrated intensities (IT) for our suite of tracers: C18O (panel a), N2H+ (panel b), HCO+ (panel c), and DCO+ (panel d). All the IT values are computed following Tafalla et al. (2021) and are compared to N(H2) pixel by pixel (Lombardi et al. 2014).
In the text
thumbnail Fig. 3

Results for OMC-2 and OMC-3. Panel a: depletion map (ƒD) obtained from the ratio of the expected abundance of C18O (X0(C18O)) in Orion and the measured abundance (Ntot(C18O)/N(H2)). The black contours represent I(N2H+) > 4 K km s−1 (signal-to-noise ratio of 10; see Fig. A.1, panel b). Panel b: abundance of o-H2D+ (X(o-H2D+)) inferred from the C18O depletion factor using the correlation from Sabatini et al. (2020). Panel c: deuteration fraction, determined as RD = N(DCO+)/N(HCO+), in OMC-2 and OMC-3 (see text for its determination). Panel d: ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ determined with Eq. (1). The white star represents a B-type star within the selected region boundaries. The dotted red circle is a candidate expanding CO-shell (Feddersen et al. 2018) possibly connected to the star.
In the text
thumbnail Fig. 4

Comparison between RD, ƒD and TK in our targets. Left panel: deuteration fraction (RD) compared to the kinetic temperature (TK). Each point is colour-coded with the corresponding depletion factor (ƒD). Right panel: deuteration fraction (RD) now compared to the depletion factor (ƒD) and colour-coded with the corresponding kinetic temperature (TK).
In the text
thumbnail Fig. 5

Comparison between ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ and the column density N(H2). The theoretical predictions for the CR propagation in molecular clouds are displayed as colour-coded lines (blue, green, and red lines; Padovani et al. 2022), while the ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ estimates in high-mass regions are displayed in grey (squares and diamonds; Sabatini et al. 2023). The estimates of ζH2ion$\zeta _{{{\rm{H}}_2}}^{{\rm{ion}}}$ from our study are displayed as circles, colour-coded by Tĸ.
In the text
thumbnail Fig. A.1

ISF as seen in different types of emission. Panel (a): Column density map of H2 (N(H2)) obtained from dust opacity Herschel+Planck observations at 850 µm (Lombardi et al. 2014). The white contours correspond to [1022, 5 × 1022, 1023] cm−2, and the blue contour represents the footprint of our IRAM-30m observations. Panel (b): N2H+ (1–0) integrated-intensity map of observed with IRAM-30m at a resolution of 30″ (see Hacar et al. 2017, for additional information). Panel (c): Kinetic temperature of the gas (TK) obtained from observations of the two inversion transitions NH3 (1,1) and NH3 (2,2) at 32″ with the GBT (Friesen et al. 2017). The white and yellow stars represent the B- and O-type stars in the region, respectively, as in Fig. 1, panel (b).
In the text
thumbnail Fig. A.2

Radial profile sampling of the DCO+ integrated intensity in OMC-2 and OMC-3. Left panel: Integrated-intensity map of DCO+ (3–2) (see also Fig. 1), masked based on the criteria discussed in Sects. 3.2,4.1. The dashed red line represents the axis used for the radial sampling. Right panel: Average intensity (Ibs) and temperature (TK) radial profiles computed across the red line in the left panel. The Gaussian best fit to the profile is also displayed (solid black line), along with the derived FWHM (solid arrow).
In the text
thumbnail Fig. B.1

Estimates of RD = N(DCO+)/N(HCO+) and TK across different studies: estimates in Orion, both single fields and median per bin (green and red circles, respectively; this study), in IRDCs (black diamonds; Miettinen et al. 2011), in high-mass star-forming regions (colour-coded squares; Gerner et al. 2014, 2015, see text for discussion on TK). The black line represents the best fit to the median values of our estimates in Orion and it is exemplified by Eq. B.1.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.