Deuterium fractionation in cold dense cores in the low-mass star forming region L1688 ★

In this work, we study deuterium fractionation in four starless cores in the low-mass star-forming region L1688 in the Ophiuchus molecular cloud. We study how the deuterium fraction ( 𝑅 𝐷 ) changes with environment, compare deuteration of ions and neutrals, core centre and its envelope, and attempt to reproduce the observed results with a gas-grain chemical model. We chose high and low gas density tracers to study both core centre and the envelope. With the IRAM 30 m antenna, we mapped N 2 H + (1–0), N 2 D + (1–0), H 13 CO + (1–0) and (2–1), DCO + (2–1), and 𝑝 -NH 2 D(1 11 –1 01 ) towards the chosen cores. The missing 𝑝 -NH 3 and N 2 H + (1–0) data were taken from the literature. To measure the molecular hydrogen column density, dust and gas temperature within the cores, we used the Herschel/SPIRE dust continuum emission data, the GAS survey data (ammonia), and the COMPLETE survey data to estimate the upper limit on CO depletion. We present the deuterium fraction maps for three species towards four starless cores. Deuterium fraction of the core envelopes traced by DCO + /H 13 CO + is one order of magnitude lower ( ∼ 0.08) than that of the core central parts traced by the nitrogen-bearing species ( ∼ 0.5). Deuterium fraction increases with the gas density as indicated by high deuterium fraction of high gas density tracers and low deuterium fraction of lower gas density tracers and by the decrease of 𝑅 𝐷 with core radii, consistent with the predictions of the chemical model. Our model results show a good agreement with observations for 𝑅 𝐷 (N 2 D + /N 2 H + ) and R 𝐷 (DCO + /HCO + ) and underestimate the 𝑅 𝐷 (NH 2 D/NH 3 ).


INTRODUCTION
Pre-stellar cores are cold, dense, quiescent ( ≃ 10 K, (H 2 )=10 4 -10 7 cm −3 , e.g., Ward-Thompson et al. 1999;Pineda et al. 2023) structural elements of molecular clouds.Under these physical conditions, CO, the most abundant molecule after H 2 in cold gas, freezes onto dust grains.With low abundance of CO in the gas phase, deuteronproton exchange reaction, that starts deuteration, can proceed efficiently (Dalgarno & Lepp 1984), which leads to an increase in deuterium fraction (e.g., Crapsi et al. 2005).Deuterium fraction becomes an important instrument to study the pre-stellar phase; it is measured as a column density ratio of deuterated and hydrogenated isotopologues (e.g.Bacmann et al. 2003;Crapsi et al. 2005).Study of deuterium fractionation is important to understand the chemical processes in pre-stellar cores, the first stages of star formation, and to put constraints on chemical models (e.g., Caselli & Ceccarelli 2012;Ceccarelli et al. 2014).
In molecular clouds, at high gas densities and low temperatures, various species are adsorbed onto the dust grain surfaces depending on their desorption energies and formation routes.Carbon-bearing species form early in the molecular cloud and are known to start † E-mail: petrashkevich.igor@gmail.comdepleting catastrophically at gas densities above 10 4 cm −3 (Caselli et al. 1999).Some nitrogen-bearing species, such as NH 3 and N 2 H + , need molecular nitrogen to be formed, that is they need longer chemical evolution of the cloud (Hily-Blant et al. 2010).Hily-Blant et al. (2010) showed that at densities of 10 4 cm −3 , the typical time scale to convert N to N 2 is of the order of 10 6 yr, while the chemical age of a typical pre-stellar core is ∼ 10 5 yr (e.g., Tafalla & Santiago 2004;Jiménez-Serra et al. 2021;Punanova et al. 2022).This implies that there is a large gas-phase N 2 reservoir to form N-bearing species in pre-stellar cores, while N 2 is actively formed.Thus such C-bearing species as HCO + and their isotopologues mostly trace a less dense gas of core shell, where CO is abundant, while such N-bearing species as N 2 H + and NH 3 , the products of N 2 , mostly trace dense core gas.In centrally concentrated prestellar cores, N-bearing species are also expected to freeze out onto dust grains deep inside (within the central few thousand a.u.), leaving only species lighter than He in the gas phase (Caselli et al. 2022).Measuring the deuterium fraction (  ) in different tracers gives us an independent way to study deuterium fractionation under different physical conditions, in addition to the direct study of the correlations between the deuterium fraction and the physical conditions (gas density, gas and dust temperature) or linear scale (distance to the dust peak) within the cores.
The target of our study, L1688, is a site of low-mass star formation, a part of the -Ophiuchus molecular cloud, where active star formation takes place (Loren et al. 1990).L1688 is located at a distance of 119 pc (Zucker et al. 2019), which makes it a good target for an observational study.This region contains tens of starless and protostellar cores (e.g., Di Francesco et al. 2008).Friesen et al. (2017) and Choudhury et al. (2021) evaluated the gas temperature in all L1688 through ammonia observations and found it to be in a range of 10 K (towards cold dense cores) to 30 K (towards the clumps with protostars).Ladjelate et al. (2020) showed the dust temperature (12-24 K) of the entire L1688 region based on the Herschel/SPIRE data.L1688 consists of several clumps, A, B, C, D, E, F, H, I that have different size, mass, number of dense cores and embedded protostars, temperature, level of turbulence (Motte et al. 1998;André et al. 2007;Di Francesco et al. 2008;Friesen et al. 2017).
Deuterium fractionation in L1688 was studied in several works.Parise et al. (2011) and Bovino et al. (2021) studied o-H 2 D + /p-D 2 H + towards the pre-stellar core Oph-H-MM1 and other cores in Ophiuchus.Friesen et al. (2010) presented the deuterium fraction map in N 2 D + /N 2 H + for Oph-B2, a star-forming clump with cold dense cores and embedded protostars.Harju et al. (2017) presented a detailed study of deuteration in N 2 H + and NH 3 towards Oph-H-MM1 based on single-pointing observations of all their possible isotopologues including doubly and triply deuterated ammonia.Punanova et al. (2016) presented single-pointing observations of N 2 H + , N 2 D + , and C 17 O towards 40 dense cores.They revealed high level of deuteration (  =0.02-0.43)and low level of CO depletion in all sub-regions of L1688.Some of the dense cores showed very high deuteration with very low CO depletion which contradicts the current understanding of deuterium chemistry in cold cores (  =0.30-0.43 with CO depletion ≤ 4.4; Punanova et al. 2016).
We present for the first time a spatial distribution of deuterium fraction in four dense cores in three different pairs of species, N 2 H + and N 2 D + , NH 3 and NH 2 D, H 13 CO + and DCO + , and search for correlations between deuterium fractions and physical parameters of the gas.We consider the N-bearing species as tracers of dense gas and the C-bearing species as tracers of less dense core envelopes.For the mapping, we chose four of the highly deuterated cold dense cores that show low CO depletion: Oph-C-N, Oph-E-MM2, Oph-F, Oph-H-MM1, located in four sub-regions of L1688, three of them located far from protostars, to avoid the protostellar feedback.We also do a single-pointing probe towards the cores Oph-H-MM2 and Oph-I-MM2, newly identified via ammonia observations (Friesen et al. 2017).The mapped cold dense cores are stationary, bound by pressure (the virial parameters of the cores were estimated in Pattle et al. 2015;Kerr et al. 2019;Singh et al. 2021).).We model the chemical composition for the studied cold dense cores using the gasgrain chemical model pyRate3 presented in Sipilä et al. (2015a,b) and compare two different approaches to reproduce the physical profiles of the cores and two initial abundances of nitrogen.
Section 2 describes the observations presented in this work and other observations data taken from the literature.Section 3 describes the processing of observational data and analysis of the spectral lines.Section 4 describes the results of the line analysis, the obtained column densities of the species and the deuterium fraction.Section 4.7 describes the chemical modeling of the studied dense cores and the comparison of the model results with the observations.In Sect.5, we discuss our line analysis method, the obtained deuterium fraction values, the difference in deuterium fraction between the tracers, and the performance of the chemical model.In Section 6 we present our conclusions.

The data from the literature
To measure deuterium fraction, we used the data on the hydrogenated isotopologues where they were available in the literature.The N 2 H + (1-0) data for Oph-F and Oph-E-MM2, observed with the IRAM 30 m antenna and SIS heterodyne receiver with an autocorrelation spectrometer as backend, were taken from André et al. (2007).These observations have spectral resolution 20-40 kHz and the beam size 26.4 ′′ , rms≃0.1-0.2K ( mb ).The p-NH 3 (1,1) and (2,2) observations for Oph-C-N, Oph-E-MM2, and Oph-F observed with the GBT telescope (32 ′′ beam) and the KFPA spectrograph, were taken from the GAS survey by Friesen et al. (2017).These observations have spectral resolution of 5.7 kHz and rms of 0.01 K ( mb ).To search for the correlations between deuterium fraction and physical conditions within the cores, and to construct the core profiles for modelling, we    Friesen et al. 2017).
used dust temperature and molecular hydrogen column density maps based on Herschel/SPIRE observations from Ladjelate et al. (2020), with effective resolution of 18.2 ′′ .The gas temperature was also taken from the GAS survey (Friesen et al. 2017), as mentioned above.
The 13 CO(1-0) observations of the COMPLETE survey (Ridge et al. 2006) obtained with the FCRAO antenna (46 ′′ beam and 0.07 km s −1 spectral resolution) were used to estimate the lower limit for the CO column density and the upper limit for the CO depletion.

Spectral maps
We convolved all maps to the same angular resolution of 33.6 ′′ , corresponding to the largest beam size in our data set.The pixel size was 12 ′′ , consistent with the Nyquist criterion.The map of H 13 CO + (2-1) for Oph-C-N was smoothed to the velocity resolution of 0.134 km s −1 to increase the signal-to-nose ratio.Up to the stage of spectral cubes, the spectra were processed with GILDAS 1 software.

Line analysis
For line analysis, we took the parts of the spectral cubes with the uniform noise (shown in Table 2) and signal-to-noise ratio > 4.
All observed transitions have hyperfine structure (hfs).We fitted the spectra to find the line parameters using the python package pyspeckit2 (Ginsburg & Mirocha 2011a;Ginsburg et al. 2022).The package fits the hyperfine structures using the frequencies and the relative intensities of the hfs components under the assumption of local thermodynamic equilibrium.The fitted spectral parameters are the transition excitation temperature  ex , optical depth , radial velocity  LSR , and the velocity dispersion  (connected to the full width at half maximum, FWHM, as The package uses the gradient descent method and the radiation transfer equation to determine the transition excitation temperature and optical depth.
For each parameter map, pyspeckit creates a parameter uncertainty map.Estimations show that the radial velocity and line width have uncertainties between 0.001-0.010km s −1 .We corrected velocity dispersion for the channel width: where  obs -observed velocity dispersion, Δ res -channel width.
In the majority of the pixels in all maps except for NH 3 ,  and  ex had large uncertainties, 30-100% of the values.To avoid such large uncertainties, we fix the parameter that affects the resulting column density the least.We explored the parameter space to find any correlation between the parameters.We found strong correlation between optical depth and excitation temperature (see Fig. A1).The spread of excitation temperature values has a smaller range (about 4 K) than that of the optical depth (up to 10).The variation of excitation temperature affects the column density less than the variation of optical depth, which affects the column density linearly: at the temperatures of 4-8 K, the column density increases with temperature by 25-30% per 1 K for N 2 H + and N 2 D + , by 12-15% per 1 K for DCO + , and by 2-5% per 1 K for NH 2 D. Thus we decided to use a single excitation temperature value for each line towards each core to estimate the optical depth and get a better estimate of the column density.To find the most suitable excitation temperature value, we explored the parameter space of the optical depth and the excitation temperature and produced the probability field with the Monte-Carlo method (see details in Appendix A).
We show the  ex found via parameter space exploration in the last column of Table 2.The uncertainty of the excitation temperature was estimated as the standard deviation of the average excitation temperature in the most probable area and amounts to 10-20%.We could not find the most probable excitation temperature for the H 13 CO + (1-0) line.We considered the H 13 CO + (1-0) line optically thin, fixed  = 0.1, and used  ex =10 K to determine the rest of the parameters as was done in Punanova et al. (2016).Since the H 13 CO + (2-1) line is very narrow and unsaturated, we also considered it optically thin and fitted with a Gaussian.The optical depth  is > 0.1 in the maps for all other transitions.
To estimate the goodness of the fit, we subtracted the resulting spectrum model from the observed spectrum.The rms of the subtraction does not exceed the rms noise level of the observed spectrum (the difference is within 10%).Therefore, the model spectra with the fixed excitation temperature were sufficient.The spectra with the highest signal-to-noise ratio and their fits for each spectral cube (for each core and tracer) are shown in Fig. A2.
For NH 3 hfs line analysis, we used the method described in Rosolowsky et al. (2008); Ginsburg & Mirocha (2011b).We used the observational maps of para-NH 3 (1,1) and para-NH 3 (2,2) from the Green Bank ammonia survey (GAS, Friesen et al. 2017) to obtain the NH 3 column density.The spectral parameters were fitted with good accuracy (median relative discrepancy between the model and the spectra = 13%), since two spectral transitions are used in the fitting.The method 'fiteach' takes into account the / ratio for NH 3 as a parameter, which we set to / = 1:1 following Harju et al. (2017), and gives the total column density of all (ortho + para) NH 3 .

The distribution of emission
Figure 2 shows the 60% contours of the peak integrated intensity of the observed transitions that represent the distribution of the observed tracers.The colored dots show the integrated intensity peaks.In Oph-C-N, where we have all six tracers observed, we see that the emission of the carbon-bearing DCO + shows the most extended spatial distribution.Among the nitrogen-bearing species, deuterated isotopologues (N 2 D + and NH 2 D) show more compact emission compared to those with hydrogen.NH 2 D shows the most compact emission and appears to be the best tracer of dense gas in our selection of lines.H 13 CO + (2-1) shows a very small emission area because of the low signal-to-noise level of the spectra.The integrated intensity maps are presented in Fig. B1.

Spectrum parameters
Figures B2, B3 and B4 show the resulting fit parameters: the maps of the optical depth  (Fig. B2, except for NH 3 where the fit gives directly  tot ), the radial velocities  LSR (Fig. B3), and the velocity dispersion  (Fig. B4) are presented for each transition line in each core.The used excitation temperature is given in the last column of Table 2.The parameter maps are masked to show only the values with the relative uncertainty smaller than 1/3.

Velocity dispersions
The velocity dispersion for all lines in the studied dense cores ranges from 0.05 km s −1 to 1.00 km s −1 (see Fig. B4).The maximum velocity dispersions of 0.85 and 1.00 km s −1 are observed in N 2 H + and NH 3 towards Oph-F core (row 3, columns 3 and 6 in Fig. B4).
The rest of the species show moderate and low velocity dispersions, up to 0.34 km s −1 , with median values of 0.07-0.17km s −1 .The lowest velocity dispersion are observed in NH 2 D, 0.05-0.14km s −1 with a median value of 0.10 km s −1 (see the third column in Fig. B4).
The general picture of narrow lines (=0.10-0.15km s −1 ) is sometimes disturbed with spots of larger velocity dispersion (=0.20-0.30km s −1 ).Visual inspection of the spectra towards the northeastern side of Oph-C-N (see Fig. B4, row 1), the strip of increased velocity dispersion as seen in all lines, reveals the presence of a second velocity component.The velocities of the two peaks differ by 0.3-0.4km s −1 .The north-western part of Oph-F (see Fig. B4, row 3) also shows an increased velocity dispersions in N 2 H + and N 2 D + .The N 2 H + (1-0) lines show a second velocity component there (also reported in André et al. 2007;Punanova et al. 2016).Additional components that affect the overall velocity field were also found towards the Oph-H-MM1 core in Pineda et al. (2022).However, the study of gas kinematics is out of scope of this paper, thus we do not analyse the two velocity components.The cores Oph-E-MM2, Oph-H-MM1, Oph-H-MM2, Oph-I-MM2 do not show the presence of the second velocity components in our data.Ammonia maps (the far right column in Fig. B4) show up to 0.07 km s −1 larger velocity dispersion compared to the other Nbearing species within the cores (the three left columns in Fig. B4), and, outside of the cores, up to 0.7 km s −1 larger velocity dispersion compared to all other species, if we define the cores as the area with subsonic velocity dispersion (as in, e.g., Pineda et al. 2010;Chen et al. 2019).The large dispersions are due to the contribution of the surrounding cloud velocity component (Choudhury et al. 2020) which has not been taken into account in our analysis.However, towards the cores, the narrow lines emitted by the core dense gas dominate the emission, and the NH 3 velocity dispersion is similar to that of the other species (see Fig. B4) and shows stronger decrease towards the core centres than, e.g., N 2 H + (see the ratio of the described below non-thermal components in Fig. B5), similar to the result of Pineda et al. (2021), who found that, in dense gas, NH 3 shows smaller velocity dispersions than N 2 H + .
We estimated the non-thermal component of velocity dispersion as follows (Myers et al. 1991): where  is the Boltzmann's constant,   is the gas kinetic temperature taken from Friesen et al. (2017), and  obs is the mass of the observed molecule in a.m.u..The thermal component is 0.053-0.070km s −1 for the selected species.The velocity dispersion in all cores is low (median dispersions are 0.07-0.17km s −1 ), which indicates mostly thermal or subthermal motions, that is, low turbulence.Figure 3 compares the non-thermal components of velocity dispersions of all species (right panel) with the thermal velocity dispersion (black line) at 12 K, which is the median gas temperature in the studied cores.We focus on the emission of the deuterated species, N 2 D + , DCO + and NH 2 D as the least contaminated by the emission of the surrounding cloud gas (left panel).The velocity dispersion in the majority of the data points is subthermal.The velocity dispersion in DCO + , ∼0.16 km s −1 (the shells of the cores), is greater than that in N 2 D + and NH 2 D ∼0.11 km s −1 and ∼0.08 km s −1 (the centers of the cores).NH 2 D shows the narrowest lines which may indicate that they originate in even colder gas in the core centers.

Column density
The formula to calculate the column density can be obtained from the radiation transfer equation.For this it is necessary to write down the intensities of the rotational transition and the absorption coefficient in terms of the number of molecules.We assume that cold dense cores are in local thermodynamic equilibrium (e.g., Shirley 2015b).Column density of NH 3 was obtained from the fit; for the rest of the species (except for H 13 CO + ) we calculate the column density as it was done in Caselli et al. (2002) for optically thick transitions: where  is the velocity dispersion,   and   are the statistical weights of the lower and upper energy levels,  is the transition wavelength,  ul is the Einstein coefficient,  is the optical depth, ℎ is the Planck constant,  is the transition frequency,  is the Boltzmann constant,  ex is the excitation temperature,  rot is the partition function,   is the energy of the lower level.The partition function of a linear molecule (N 2 H + , N 2 D + , H 13 CO + , DCO + ) can be calculated as follows: where  is the rotational quantum number,   -rotational transition energy ( ( + 1)ℎ),  is the gas temperature.The partition function for deuterated ammonia was taken from the CDMS data base (Endres et al. 2016).All quantum coefficients used to calculate the column densities are listed in Table 3.Since we assumed that the H 13 CO + (1-0) and (2-1) lines are optically thin, the column density of H 13 CO + was calculated as in Caselli et al. (2002) for optically thin transitions: where  is the integrated intensity,   () is the equivalent Rayleigh-Jeans temperature,  bg is the background temperature of 2.7 K.
The column density of ammonia is calculated using the pyspeckit fitter (see Sect. 3).The / ratio of (NH 3 ) is introduced in the fitter, we set it to 1:1 (the last parameter of the ammonia fitter is fixed to 0.5) following Harju et al. (2017).
Figures 4, 5, and 6 show the column density maps of all species and Table 4 shows the column density for the single-pointing observations towards Oph-H-MM2 and Oph-I-MM2.In Oph-C-N and Oph-E-MM2, the column density of all species increases towards the dust peaks in cores.In Oph-F, the molecular emission peaks are slightly shifted to the south-east from the dust peak by one to half beamsize.In Oph-H-MM1, there are two molecular peaks (the two peaks apper prominently in the interferometric view of NH 3 ; Pineda et al. 2022), and N 2 D + shows strong emission across all the core.The column density is ∼10 12 cm −2 in N 2 H + , N 2 D + , ∼10 13 cm −2 in -NH 2 D, and ∼10 14 cm −2 in NH 3 for all cores.In H 13 CO + and DCO + , column densities are ∼10 11−12 cm −2 for Oph-E-MM2 and Oph-C-N.
In Table 5, we compare our column densities of N 2 H + and N 2 D + to those measured by Punanova et al. (2016) towards the dust emission peaks.The positions of their observations are marked with white crosses in our maps.Our column densities are lower then those in Punanova et al. (2016), although not systematically and some agree within the uncertainties, with the only exception of N 2 D + towards Oph-F where our column density is higher.This is likely due to

Note:
The partition functions  rot , lower level energies   , and Einstein coefficients  ul of the linear molecules were calculated as in Caselli et al. (2002).For that, the rotational constants and dipole moments are taken from the works, listed in the right column* and are also available through the CDMS data base (Endres et al. 2016).The Einstein coefficient ( ul ) and the lower level energy (  ) for -NH 2 D(1 11 -1 01 ) were taken from Harju et al. (2017). rot for -NH 2 D(1 11 -1 01 ) was taken from the CDMS data base (Endres et al. 2016).The critical densities ( 10K  ) of hydrogenated species were taken from Shirley (2015a) and those of deuterated species were calculated in the same manner as in Shirley (2015a) using the collisional rates from the works listed in the right column** and available through the LAMBDA data base (Schöier et al. 2005).
Table 3.Quantum coefficients used to calculate the column densities.Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.the fact that our adopted  ex or measured optical depths are lower than those measured in Punanova et al. (2016).The reason why their  tot (N 2 D + ) in Oph-F is lower is that they adopted  = 0.1 for the N 2 D + line while we measured =1.6 towards that position.Our maps also show that some of the pointings in Punanova et al. (2016) missed the line emission peaks.Our NH 3 column density maps agree with the map in Friesen et al. (2017), which is expected since we used their data and the same line analysis method.
Since some lines are subthermally excited, with the LTE ap-proximation, we might underestimate the column densities (e.g., Shirley 2015b).The higher transitions have higher critical densities (see Table 3) and also have more tendency to subthermal excitation, that is either way underestimated column density.The combination of a ground transition of hydrogenated isotopologue and higher transition of deuterated isotopologue would then lead to an underestimated deuterium fraction.However, the example of cores H-MM2 and I-MM2, where we have both (1-0) and (2-1) N 2 D + transitions observed, shows the opposite result (see where and agree within the errors.With RADEX (van der Tak et al. 2007), we estimated the column densities of the species towards the dust peaks of the cores, where volume density is measured, and found the difference within 10-20%, both smaller and larger values.Non-LTE estimates of the column densities require the non-LTE radiative transfer modelling and thus the data on volume density distribution, that is the farther from the dust peak, the more model assumptions on the source geometry are involved.Thus, we consider LTE approximation the most optimal for column density mapping of dense cores.

Deuterium fraction
To estimate the deuterium fraction, we divided pixel by pixel the resulting column density maps of deuterated isotopologues by those of hydrogenated isotopologues.To calculate the column density of HCO + , we applied the 12 C/ 13 C=77 isotopic ratio for local ISM (Wilson & Rood 1994) to our (H 13 CO + ) maps and suggested that fractionation of 13 C is negligible.We estimated the total deuterated ammonia column density based on our column density of para-NH 2 D using o/p(NH 2 D) = 3:1 obtained by Harju et al. (2017).Figure 7 shows the maps of the deuterium fractions in all tracers of the studied cold dense cores.Table 4 shows the deuterium fraction   (N 2 D + /N 2 H + ) derived from the single-pointing observations towards Oph-H-MM2 and Oph-I-MM2.
Deuterium fractions appear high over entire maps and show a significant difference in the values of nitrogen-bearing   (N 2 H + ) and   (NH 3 ) (peak value 0.6-0.8)and carbon-bearing   (HCO + ) (peak value 0.06-0.10)species.Despite the small extent, the obtained maps show that the deuterium fraction increases towards the dust peaks of the cores as it was observed in other dense cores (e.g. in L1544, L1521F, L183; Caselli et al. 2002;Crapsi et al. 2004Crapsi et al. , 2007;;Pagani et al. 2007Pagani et al. , 2009a;;Chacón-Tanarro et al. 2019;Redaelli et al. 2019), most prominently seen in the largest map of   (HCO + ) to-wards Oph-E-MM2.Only in Oph-F, where protostars are present, the maximum   (N 2 D + /N 2 H + ) is away from the dust peak.However,   (N 2 H + ) is still high (≃ 0.4) towards the protostars, and   (NH 3 ) is the highest towards one of the protostars.This means that the early stage protostars (Class 0-I) have not yet heated the dense gas of the core, or the heating of gas and dust has a very limited extent, diluted by our beam, or that they are not embedded in this core, being in the background.Towards the cores where the map size allows to see it (Oph-F and Oph-H-MM1), the deuterium fraction of N 2 H + is high (  > 0.4) over a large area (several beam sizes) which excludes the idea of very compact depletion zones with high deuteration proposed in Punanova et al. (2016).
We consider carbon-bearing H 13 CO + and DCO + core shell tracers since they are depleted at high gas densities due to CO freeze-out, and nitrogen-bearing NH 3 , NH 2 D, N 2 H + , and N 2 D + core centre tracers since they stay abundant at high densities of cold cores (e.g., Caselli et al. 2002;Tafalla et al. 2006) and only get depleted at densities ≥ 10 6 cm −3 (Caselli et al. 2022).However, the large-scale Green Bank Ammonia Survey showed that NH 3 traces also the lowdensity gas of the clouds (Friesen et al. 2017).We still consider NH 3 as a dense gas tracer since it stays undepleted within the cores (as indicated, e.g., by the lines, narrowing towards the core centre, see Fig. B5) and recall that it also traces the low-density gas on the line of sight.Ammonia begins to deplete only at high densities (> 2 × 10 5 cm −3 ) and smaller spatial scales (< 1000-3000 au Pineda et al. 2022;Caselli et al. 2022).
In Fig. 8 we assess the difference in the deuterium fractions of nitrogen-bearing and carbon-bearing species.The lines show the difference in deuterium fraction between N 2 H + and NH 3 (factors of 0.5, 1, and 2, left panel) and between the N-and C-bearing species (factors of 2, 5, and 10, central and right panels).The deuterium fraction in nitrogen-bearing species is 2-10 times greater (the median difference is a factor of 4.2) than the deuterium fraction in carbonbearing species.The deuterium fraction in ammonia is lower than  Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.that in N 2 H + (the median difference is a factor of 1.5).The difference in deuterium fraction between the tracers indicates that it is greater in the core center than in the envelope.

Deuterium fraction and the physical parameters of the cores
To study how the physical conditions within the cores affect deuterium fraction, we analysed the correlations between the deuterium fraction of each species and the molecular hydrogen column density ( tot (H 2 )), the dust temperature ( dust ), the gas temperature ( gas ) and gas turbulence, represented by the non-thermal component of the velocity dispersion ( NT ).Kinetic temperature rules the H 2 / ratio and thus the deuterium chemistry.The correlation between deuterium fraction and the temperature is predicted by theoretical works (e.g., Pagani et al. 2009a;Parise et al. 2011;Kong et al. 2015).We used the dust temperature and molecular hydrogen column density maps from Ladjelate et al. ( 2020) based on the Herschel/SPIRE dust observa-tions at 250 m, 350 m, and 500 m, and the gas temperature map of L1688 from Friesen et al. (2017).The non-thermal velocity dispersion was measured in our fits using the gas kinetic temperature from Friesen et al. (2017) as described in Section 4.2.The non-thermal component is associated with turbulence, which increases the collisions of particles increasing the chances of species destruction and thus the number of deuterated molecules is decreased.The physical parameter maps were convolved to the common resolution (33.6 ′′ ) of our data, for consistency.
We fitted the linear approximations to estimate the correlation between the deuterium fractions and the physical parameters.We found no significant correlations between the deuterium fraction, the gas and dust temperature, molecular hydrogen column density ( tot (H 2 )) and the non-thermal component.The gas temperature  gas dynamical range of only 2 K is probably too small to show any correlation.The lines are mainly subthermal, and the maps of deuterium fraction do not reach the 'transition to coherence' (Pineda et al. 2010; Chen et al. 2019), which means also small variation in turbulence.The correlations between deuterium fraction, molecular hydrogen column density ( tot (H 2 )) and dust temperature have large uncertainties in the parameters of approximating linear functions due to the small number of points.Maps of deuterium fraction spanning beyond the 'transition to coherence' are needed to show the correlation between deuterium fraction and physical parameters.

CO depletion
Deuterium fraction is expected to increase with CO freeze out (Bacmann et al. 2003;Caselli et al. 2008;Kong et al. 2015), since CO, the second most abundant molecule in molecular cloud gas after H 2 , reacts and destroys the H + 3 ion, which is needed to support the deuterium chemistry via the H + 3 +HD ⇌ H 2 D + +H 2 reaction.To test this prediction with our data set, we took the 13 CO(1-0) observations from the COMPLETE survey Ridge et al. (2006) and produced the CO column density maps.The spectral lines of 13 CO(1-0) were partly optically thick and had several unresolved velocity components.We were unable to estimate the optical depth of the line, thus we calculated column density as if the lines were optically thin, as described in Sect.4.3 and since the lines are partly optically thick the column density therefore represents its lower limit.We converted the column density of 13 CO to the column density of CO using the 12 C/ 13 C fractional abundance of 77 (Wilson & Rood 1994).We compared the obtained values of the column density towards the cores dust peaks to those measured via the C 17 O observations from Punanova et al. (2016).Our estimates of the CO column density done with 13 CO(1-0) are 3-8 times lower than the values from Punanova et al. (2016), depending on the core.Even with the upper limit, the distribution of   is consistent with theory, it is higher towards the dust peak where the deuterium fraction should be the highest.However we note that while 13 CO(1-0) is obviously optically thick towards the core dust peaks, its intensity decreases at lower densities, thus the estimated column density and   are more likely to be accurate at the edges of the core maps.We took the reference CO abundance ((CO) ref = 2.69×10 −4 ) for Ophiuchus from Lacy et al. (1994) and plotted the maps of the upper limit of the CO depletion: where (CO) obs =(CO)/(H 2 ). Figure B6 shows the maps of the Table 6.The upper limits of the CO (from 13 CO) depletion towards the observed starless cores.
CO depletion upper limit for our cores based on 13 CO.The maximum values of CO depletion are given in Table 6.Since we have only lower limits on the CO column density, the distribution of the CO depletion is rather affected by the column density of H 2 , and the correlation between the deuterium fraction and CO depletion factor is very similar to that between the deuterium fraction and molecular hydrogen column density.The analysis of C 17 O(1-0) presented in Punanova et al. (2016) suggested that there is no CO depletion towards the dust peaks of Oph-C-N, Oph-E-MM2, and Oph-F.Despite the upper limits of   towards the dense cores, the advantage of these maps is that they show moderate depletion (   =1-2) towards the cores outskirts where the 13 CO(1-0) line have low optical depth and where the column density of CO might be measured accurately.

Ionisation fraction
The ionization degree indicates the strength of the coupling between the magnetic field and the dense core.A strong magnetic field can affect the stability of the dense core, thereby preventing contraction or fragmentation.The simplest way to estimate the ionization degree is to count all available observed ions: as it is described in Caselli et al. (2002).HCO + and N 2 H + , and their deuterated isotopologues are the most abundant ions in dense gas after H + 3 , the latter we did not observe, so this method can only estimate the lower limit of the ionization degree.HCO + contributes the most to our estimate of the ionisation degree.In our data set, we observed all four ions towards two cores, Oph-E-MM2 and Oph-C-N.However, for Oph-C-N, HCO + , the most abundant of the four ions, is represented by the H 13 CO + (2-1) line at 173.5 GHz, where atmospheric absorption is strong, thus the noise temperature is high, and the high uncertainty in the integrated intensity propagates to the column density and the estimate of the ionization degree, so the uncertainty of the latter is very high.Thus we performed the ionization analysis only for Oph-E-MM2.The left panel of Fig. B7 shows the ionization degree in Oph-E-MM2.The maximum ionization degree (() = 3.05 ± 0.9 × 10 −9 ) is similar to that of the prototypical pre-stellar core L1544 (() ∼ 10 −9 ; Caselli et al. 2002).
We also estimated the ionization fraction and the cosmic ray ionization rate for Oph-E-MM2 using the analytical definition of [DCO + ]/[HCO + ]≡   and [HCO + ]/[CO]≡   based on a simple steady-state model, where   and   are defined via formation, destruction and dissociation reaction rates, ionization rate, ionization degree and HD abundance (for the details, see Sect. 3 in Caselli et al. 1998).Following this approach, we estimate the ionization degree at the temperature of 10 K, characteristic for a cold core and already implemented in the model: and cosmic rays ionization rate : where   is the depletion factor, the numerical coefficients 2.7×10 −8 and 1.2 × 10 −6 , 7.5 × 10 −4 and 4.6 × 10 −10 are obtained by substituting reaction constants following Caselli et al. (1998).The value of volume density (H 2 )=4×10 5 cm −3 for Oph-E-MM2 was taken from Ward-Thompson et al. (1999); for (CO) we used (C 17 O) from Punanova et al. (2016) and the isotopic ratios 16 O/ 18 O=560 (Wilson & Rood 1994) and 18 O/ 17 O=4.11(Wouterloot et al. 2005).
Figure B7 shows the ionization degree (central panels) and ionization rate (right panels) towards Oph-E-MM2 (the upper panels show the values, the lower panels show the uncertainties).The minimum values for the ionization degree are (2.0±0.2)×10−8 for the steady-state model from Caselli et al. (1998) and the lower limit of (1.00±0.03)×10−9 for observations.The ionization rate is the lowest close to the core centre: (2±8)×10 −17 s −1 .However, towards the central pixels, () and thus  become negative, since the   increases and the   remain constant (we use one value).In principle, the dust peak is the only position where we could correctly estimate () and  since it is the only position with reliable   measurement.The value is so small (   =1.25) that it leads to a negative result in Eq. 8.
To obtain here a positive (),   needs to be at least 3.83.This is another indication of the fact that the total undepleted abundance of CO in the Ophiuchus molecular cloud may be larger than a 'canonical' value, common for the local ISM (between 1 and 2×10 −4 , e.g., Frerking et al. 1982;Lacy et al. 1994) and possible variation of metallicity from one cloud to another.However, towards the core outskirts the value   =1.25 might be adequate, and the estimates of () and  may be correct.Unlike the first method that gives the lower limit of the ionization degree (3.8×10 −9 ), the model estimates the total value, which is two-three orders of magnitude higher, ∼ 10 −6 (similar to the one, ∼ 10 −6.5 , recently found towards NGC 1333 by Pineda et al., submt.), and shows that the most abundant ions are not counted, possibly H + 3 and its isotopologues.Although such high ionization degree (∼ 10 −6 ) is rather characteristic for a surrounding translucent cloud, not dense cores (e.g.Bron et al. 2021), there's some evidence that L1688 has a higher local ISRF than, e.g., towards L1544, given its nearby YSO population and warmer gas temperatures overall.To apply correctly the simple model from Caselli et al. (1998Caselli et al. ( , 2002) ) one needs to correctly estimate CO depletion that remains challenging for L1688 and requires mapping of the rare CO isotopologues down to low densities where CO is undepleted.

Chemical modelling
We use the gas-grain chemical model described in Sipilä et al. (2015a) to reproduce the observed column densities and deuterium fractions of the observed species.The model distinguishes between spin states of hydrogen and deuterium bearing species, with branching ratios derived by applying selection rules under the principle of the conservation of nuclear spin.Ortho-hydrogen affects the deuterium fractionation because the deuteron-proton exchange reaction with it is reversible even at low temperatures (Pagani et al. 1992).The model considers two phases, gas phase and dust surface, and connects the two via adsorption and several mechanisms of desorption.Chemical reactions on the dust surfaces are also included.
The elemental abundances are taken from Sipilä et al. (2015a), and correspond to the low-metal composition (e.g., Wakelam & Herbst 2008), which is more consistent with molecular cloud evolution and the low electron fractions observed there (Graedel et al. 1982).Harju et al. (2017) found that they can not reproduce the observed ammonia abundances with the same chemical model and suggested using a higher nitrogen abundance (N/H=5.3×10−5 instead of 2.14×10 −5 ), that improved the results but did not produce a perfect match.We explored how the initial nitrogen abundance affects the resulting abundances of the studied species and, in addition to the low-metal N/H=2.14×10−5 , used even higher than that in Harju et al. (2017) nitrogen abundance N/H=7.60×10−5 that corresponds to the highmetal elemental abundances (Cardelli et al. 1993;Wakelam & Herbst 2008), to find the best match with the observations.

Physical model
We model the dense cores as static spheres with gas and dust temperature and density profiles.The dust emission peaks were chosen to be the core centres.We used column density (H 2 ) and dust temperature  dust , based on the Herschel observations, from Ladjelate et al. (2020).To obtain discrete profiles of (H 2 ) and  dust , we averaged their values within concentric rings with a width of 1285 au which is equal to the pixel size (12 ′′ at the distance of 119 pc; Zucker et al. 2019) in our maps.We used the method described, e.g., in equation 1 of Hasenberger & Alves (2020) to obtain the gas density, (H 2 ), from the molecular hydrogen column density.We convert (H 2 ) to interstellar extinction,   , applying a factor of 9.4×10 20 (Frerking et al. 1982).For the density towards the centre, we took the volume densities obtained by Pattle et al. (2015), since they are based on the combination of the Herschel (250 m, 350 m, 500 m) and SCUBA2 (450 m and 850 m) observations, that is on dust emission at longer wavelength that better probe cold dust in core centre.We added one extra point at a distance of 642 au from the centre to have the central 1285 au-wide circle.This gives 9 points from the center to the edge of the core with a radius of 9000 au, so the model cores would cover entirely the extent of our maps.The dust temperature slightly decreases towards the core centers, by 2-5 K compared to the core edge (where   =15.5-17.5 K).The molecular hydrogen volume density increases towards the core centers by an order of magnitude.
We also used an alternative analytical profile based on the central density from Pattle et al. (2015) and theoretical radial density distribution from Tafalla et al. (2002).The dust temperature is constant in the alternative model.Figure 9 shows the gas density and the dust temperature profiles for our cores.Profiles 1 are based on the data from Ladjelate et al. (2020) and profiles 2 are based on the method from Tafalla et al. (2002).

Conversion of abundances to column densities
The chemical model produces the abundances of species in each position independently.To compare the model results to our observed column densities, we convert the modelled abundances to the line of sight column densities, considering the spherical geometry and convolution with the telescope beam.The total modelled column density on the line of sight is calculated as the sum of the column densities multiplied by the layer thickness in point  on the line of sight (the profile has  = 8 points): where   is the abundance in point ,   (H) is the concentration of atomic hydrogen in point , Δ  is the distance between points  and -1.For the beam towards the centre: where  0 is the abundance in the central point,  0 (H) is the concentration of atomic hydrogen in the central point, and Δ is the 1285 au step of our density profile which is here the distance between any two neighboring points on the line of sight towards the centre.

Convolution with a beam
For a correct comparison with observational results, we convolved the modelled column densities with the beam size of our observations, 33.6 ′′ .We have constructed Gaussian distribution normalized by one for each point with FWHM of 33.6 ′′ .The smoothed column density for each point along the profile is then: where   tot -column density in point ,  = 17 -points from −8 through 0 to 8 is needed for the sum over the entire length of the Gaussian,   beam -the intensity of the Gaussian in point .As a result, we have nine points of the column density convolved with our beam in the plane of the sky along the profile from the center to the edge of the core, for each species.

Model results: temporal evolution and chemical age
We ended up with four models with different input parameters: profile 1 with (N)=7.60×10−5 and (N)=2.14×10−5 and profile 2 with (N)=7.60×10−5 and (N)=2.14×10−5 .Figures B8, B9, B10, and B11 compare the modelled and observed column densities and deuterium fractions for all cores.The first and second (upper and middle upper) rows of panels show how the column densities and deuterium fractions change with the evolution time.The third and fourth (middle lower and lower) rows of panels show how the column densities and deuterium fractions change with the distance from the core centres at the chosen age.
In all four variants (physical profiles and initial abundance of N), the model showed an increase with time in the column densities of N 2 H + and of the deuterated species -N 2 D + , DCO + , and NH 2 D. N 2 D + showed the steepest increase with no decreasing at the early ages (< 10 5 yr).Both NH 3 and NH 2 D showed a plateau or even a decrease around a few 10 3 -10 4 yr.The peak of gas-phase NH 3 column density is around ∼ 10 3 yr.The HCO + column density slightly increases over the time.
To compare the model results with observations, we need to set the age of the cores.The CO depletion factor allows to estimate the core chemical age, since CO chemistry is fairly simple and well established (as was done in, e.g., Vasyunin et al. 2017).However, we have only upper limits of the CO depletion factor, so we did not use it to estimate the age.We still plot the modelled CO depletion factor (in the second rows of Fig. B8, B9, B10, and B11) over the chemical evolution to discuss later its correlation with deuterium fraction.To estimate the model CO depletion, we need to know the undepleted abundance of CO.For that, we run the chemical model with  gas =20 K,   =2, (H 2 )=200 cm −3 until 10 6 yr, and use the maximum gas-phase CO abundance reached in the diffuse cloud,  (CO) cloud = 2.5 × 10 −5 .We estimate the model CO depletion factor as where  (CO) gas is the simulated CO column density of CO in the gas phase,  (H 2 ) and  (H) are the column densities of molecular and atomic hydrogen, all produced by the chemical model of the core.
We then suggested that the observed high the deuterium fractions might be the highest possible in our cores; and that deuterium fraction of N 2 H + is the most sensitive to the change of the physical conditions.Thus, to compare the modelled and observed column density profiles, we chose the ages of the first   (N 2 D + /N 2 H + ) deuterium fraction peaks.Based on the peaks of the deuterium fraction   (N 2 D + /N 2 H + ), we obtained the chemical ages of cold dense cores.The ages vary from one the model to the other (see Fig. B8, B9, B10, B11), but not significantly.The median ages are 2.65×10 5 yr, 2.65×10 5 yr, 0.82×10 5 yr, and 1.80×10 5 yr for Oph-C-N, Oph-E-MM2, Oph-F and Oph-H-MM1, respectively.The deuterium fractions of N 2 H + , NH 3 , and HCO + show peaks at different chemical ages, all of them are around (0.8-4.1)×10 5 yr.The age of ∼10 5 yr is considered to be a typical chemical age for cold cores (e.g., Caselli et al. 1998;Vasyunin et al. 2017;Harju et al. 2017).We nevertheless analyse each variant of the model (with two nitrogen abundances and two physical profiles) and plot the radial distribution of the column densities at the ages found within each model.The ages are shown on top of each column in Fig. B8, B9, B10, and B11.Chemical ages are not related to dynamical ages because our model has a static physical structure.

Model results: column densities and deuterium fractions
While all the modelled column densities agree with the observed ones within one order of magnitude, which falls into the model accuracy (Sipilä et al. 2015a), we still see the following trends in the model results.We found that modelled column densities of HCO + and DCO + agree with the observed ones for both physical profiles and both nitrogen elemental abundances, with an underestimate for HCO + and an overestimate for DCO + at the edges of the cores (outer 6000-9000 au, see Fig. B8, B9, B10, and B11).Regardless of the density and temperature profiles applied, with the low nitrogen abundance, the modelled column densities of N 2 H + and N 2 D + are underestimated compared to the observations.The high nitrogen abundance gives an overestimate in the N 2 H + column density and either an overestimate or good agreement in the N 2 D + column density.The modelled column density of NH 3 is overestimated by all variants of the model except for the models with the low nitrogen abundance for Oph-F, where the modelled column densities agree with the observed ones.The modelled column densities of NH 2 D agree with the observed ones for the cores Oph-E-MM2 and Oph-C-N with small overand underestimates, with all four the models.All four models underestimate the column density of NH 2 D for Oph-F, models with the high nitrogen abundance give a better agreement with observations.
Although the model deuterium fractions at the chosen chemical ages agree with the observed deuterium fractions within the uncertainty of the model, we were able to identify the most suitable models.The modelled deuterium fraction   (N 2 D + /N 2 H + ) either agrees or underestimates the observed one for all cores, but profile 2 gives a better result for the entire core while the Herschel-based profile 1 shows agreement only towards the core centres.The modelled deuterium fractions   (N 2 D + /N 2 H + ) show a similar distribution along the core radii as the observed deuterium fractions.
The modelled CO depletion factor reaches 4-7 at the chosen chemical ages.Figure B6 shows the upper limit on the CO depletion factor (described in Sect.4.5) and these values are only slightly higher (   ≤8-10) than those in the model.
The analytical profile 2 causes an increased abundance of deuterated species.The high nitrogen abundance makes the biggest difference in the modelled column densities: all nitrogen-bearing species show higher column densities, the carbon-bearing species show lower column densities or remain unchanged.In cores Oph-F and Oph-H-MM1 the initial nitrogen abundance affects only the column densities of the N-bearing species, the column densities of the C-bearing species are almost the same in the four variants of the model.The results in these two cores do not depend on the chosen physical profile.

Our line fitting approach
The problem with fitting hfs structure to the spectra with moderate signal-to-noise ratio is the large uncertainty of the optical depth.One approach to treat such spectra of lines with low optical depths is to consider the lines optically thin (e.g., Crapsi et al. 2005).However, a reasonable value for optical depth is important to adequately estimate column densities.We used the Monte-Carlo method to explore the mutual correlation between  ex and , found the most probable  ex for each transition in each core, and applied a  ex -constrained fit to estimate .An incorrect  ex also affects the estimate of column density (see Sect. 3.2).As a result, towards the dust peaks, we obtained lower column densities of both N 2 H + and N 2 D + than those in Punanova et al. (2016) since either  or  ex in our fits were lower (see Sect. 4.3).
We compared our results of the deuterium fraction   (N 2 H + ) to the values from Punanova et al. (2016) (see Table 5).Despite the fact that we got lower values for the column densities, the deuterium fractions towards Oph-C-N, Oph-E-MM2, and Oph-H-MM1 agreed within the uncertainties (our values are higher) while towards Oph-F our value was higher since we measured the optical depth of the N 2 D + (2-1) line, and Punanova et al. (2016) assumed the line to be optically thin, that lead to an underestimate of  tot (N 2 D + ) and the deuterium fraction.
The map of   (N 2 H + ) towards Oph-H-MM1 was presented in Petrashkevich et al. (2020), and they also assumed the N 2 D + (1-0) emission was optically thin.Our results show the optical depth of the N 2 D + (1-0) line  ≤ 5.It means that Petrashkevich et al. (2020) underestimated the column density of N 2 D + due to the assumption of the optically thin lines.With our fitting method, we obtained a factor of ∼2 higher column densities of N 2 D + , 16% higher column densities in N 2 H + , and as a result, a factor of 1.7 higher deuterium fraction.
To conclude, a free-parameter fit gives a better result but requires high quality spectra.Our approach is a reasonable compromise to fit noisy spectra with moderate or low signal-to-noise ratio, which are always present in mapping, and results in more realistic estimates of column densities than a mere assumption of optically thin emission.
At the same time, in L1688, there are clumps with low or moderate deuterium fractions: Oph-A with   < 0.1 and Oph-B2 with   < 0.2 (Friesen et al. 2010;Punanova et al. 2016); these clumps contain many protostars (Di Francesco et al. 2008), which explains the lower values.Protostars slow deuterium fractionation and start the reverse process by heating gas and dust.Friesen et al. (2010) presented the map of deuterium fraction,   (N 2 D + /N 2 H + ), towards Oph-B2, a clump in L1688.They obtained the maximum deuterium fraction of 0.16.Several protostars embedded in Oph-B2, some of them associated with outflows (Kamazaki et al. 2003), could increase the temperature and ionization degree within the clump and cause the decrease in deuterium fraction (Friesen et al. 2010).Unlike Oph-B2, Oph-F and Oph-I-MM2, also containing embedded protostellar cores (Kirk et al. 2017), show high deuterium fractions, up to 0.3 and 0.7 in NH 3 and N 2 H + , respectively.Either the protostars are too young and have not yet affected the deuterium chemistry in the cores (or have affected small areas diluted by our 30 ′′ ≃ 4000 au beam), or they are not embedded but projected onto the cores.
The high deuterium fraction observed towards the cores supports the previously proposed idea that the cores in L1688 are close to the pre-stellar stage (Punanova et al. 2016) and consistent with the age of several free-fall times.This might be a result of the dominating pressure of the cloud where cores are fully or mostly pressure-confined (Pattle et al. 2015;Chen et al. 2019).The gas in the cores evolve chemically while the mass is not enough to start the contraction (the cores are gravitationally unbound; Pattle et al. 2015;Kerr et al. 2019), the relatively strong magnetic field ( tot =75-135 G in L1688; Hu & Lazarian 2023) also resists the contraction, so deuteration increases.The low CO depletion traced by C 17 O observed towards the core dust peak contradicts this idea as well as the fact of high deuterium fraction observed across the cores.Possibly, our estimate of CO depletion in L1688 is incorrect, and one should apply a larger reference CO abundance (we apply the largest published in literature).A recent work of Punanova et al. (2022) also showed that a larger reference CO abundance is needed to correctly estimate CO depletion in dense cores.Then the picture will be consistent: CO is depleted, deuteration is enhanced, that is the cores exist in this quasistatic phase for time much longer than free-fall time.

The three different gas density tracers
We studied the tracers of different gas density -carbon-bearing HCO + and nitrogen-bearing NH 3 and N 2 H + , and their deuterated isotopologues.
Carbon-bearing species form relatively fast and are abundant in the molecular cloud.The formation of nitrogen-bearing NH 3 and N 2 H + occurs in denser regions at later stages of molecular cloud evolution, because N 2 formation requires longer time (e.g., Hily-Blant et al. 2010).Figure 8 shows the difference between the three density tracers.Deuterium fraction of HCO + is smaller than that of NH 3 and N 2 H + by a factor of 2-12.This is due to the higher abundance of CO in the gas phase and higher temperature in the core envelope.
NH 3 is observed in molecular clouds at lower densities than N 2 H + since it is excited and probably formed at lower densities than N 2 H + .
Because of this, we trace more material with NH 3 than that with N 2 H + on the line of sight.At the same time, NH 2 D and N 2 D + are present and excited in almost the same compact and dense area, so using NH 3 to measure   , we divide by a larger column due to the excitation conditions (see Fig. 2).It is also possible that at high densities a decreasing ionization rate affects the difference in the deuteration of N 2 H + and NH 3 : N 2 D + dissociates in a reaction with a free electron, while NH 2 D needs a free electron to be formed out of the NH 3 D + ion (this reaction accounts for over 80% of the formed NH 2 D, according to the model; Sipilä et al. 2015a).
The deuterium fraction map of Oph-F is not significantly affected by the embedded protostar and   (N 2 D + /N 2 H + ) towards the protostars remain high (≃ 0.4).  (NH 2 D/NH 3 ) remains unaffected.While we can assume that reverse deuterium fractionation (e.g., Friesen et al. 2010) has started in N 2 H + , it has not started in NH 3 .This also can be caused by locally increased ionization degree (due to the presence of a protostar) that affects N 2 D + faster than NH 2 D.

Performance of the chemical model
We model the chemical composition of the studied cold dense cores using the gas-grain chemical model pyRate3 presented in Sipilä et al. (2015a,b).The model distinguishes the hydrogen isotopes, protium and deuterium, and takes into account hydrogen and deuterium spin states and thus the energy difference between the ortho-and parastates in the reactions containing up to five atoms of hydrogen or deuterium based on the spin selection rules.To reproduce the physical conditions, we used a model with two different molecular hydrogen density profiles and gas and dust temperatures profiles.The first type of profile fits the Herschel-based H 2 column densities and dust temperature from Ladjelate et al. (2020) and the second type is an analytical profile, described by Tafalla et al. (2002).We also used two different initial nitrogen abundances.
The chemical model reproduced the observed column densities within its accuracy (one order of magnitude; Sipilä et al. 2015a).However, there are trends in the model results which could be improved: the modelled deuterium fraction of HCO + is too high, while that of NH 3 is too low.The model overproduces the column density of DCO + , and as a result, the deuterium fraction   (DCO + /HCO + ).
The model overproduces NH 3 , resulting in a smaller deuterium fraction   (NH 2 D/NH 3 ) compared to observations.Previous works also could not produce large values of singly deuterated ammonia in cold cores (e.g.Roueff et al. 2005;Harju et al. 2017).The mismatch between the model and the observations may be due to several reasons.First of all, the observations account for the column of the emitting region, which is different for the observed NH 2 D and NH 3 lines due to different excitation conditions, while the model accounts for the column across the whole core.Second, there is an open question of how the deuteration reactions actually proceed; if they proceed via full scrambling (implemented in our work), one expects lower NH 2 D/NH 3 ratios than in the case of proton hop (see Sipilä et al. 2019, for the details).Third, in the model we obtain lower / ratio both in NH 3 and NH 2 D compared to those adopted to process the observing data: /(NH 3 )=0.6-0.8 and /(NH 2 D)=1.5-1.8 in the model, while we adopt /(NH 3 )=1 and /(NH 2 D)=3 (following Harju et al. 2017) to obtain the observed total column densities of NH 3 and NH 2 D. The simulated spin ratios are expected to change to statistical ones if proton hop is used (Sipilä et al. 2019).
We have the problem of comparing observational and modelling values of CO depletion.We can only estimate the upper limit of CO depletion from observations (Fig. B6).The model shows CO depletion   =4-7 (see the second row of the plots in Fig. B8, B9, B10, B11) around the deuteration peak, which is consistent with the observational values.The question of high deuteration along with the low CO depletion is still open.
We compared two density and temperature profiles: profile 1 based on Herschel observations and profile 2 based on analytical radial density distribution.The difference between profile 1 and profile 2 in the models is evident in the Oph-E-MM2, Oph-C-N and Oph-F cores, with profile 2 leading to a significant increase in the column densities of deuterated species.We assume that profile 1, constructed from Herschel observations, does not probe the dense gas on the line of sight where deuteration increases.If the physical model would include collapse, the density profile would evolve from the shape of profile 2 to the shape of profile 1, depending on the rate of infall, and would reduce the timescale of chemical evolution (e.g., Kong et al. 2015).
We compared two different initial nitrogen abundances.Changing the initial nitrogen abundance significantly affected the abundances of nitrogen-bearing species and, as a result, the deuterium fractions.The density and temperature profiles better probing the dense gas and based on the observations with higher spatial resolution could also improve the model results for the N-bearing and D-bearing species.Modeling with (N)=7.60×10−5 provided the column densities of N 2 D + , N 2 H + , and NH 2 D closer to the observed column densities than (N)=2.14×10−5 .The assumed static physical conditions limit the model performance.
(i) We analyse spectra with hyperfine structure and moderateto-low signal-to-noise ratio.We find the most probable excitation temperature and constrain it to find a realistic estimate of the optical depth.This may result in lower column density estimates compared to an unconstrained fit, but is better than treating the noisy spectra as optically thin.As auxiliary data, we present the parameter maps from our fits: optical depth, central velocity, and velocity dispersion of the lines.
(ii) We presented the column density maps of all six species towards Oph-C-N and Oph-E-MM2, the maps of N 2 H + , N 2 D + , NH 3 , and NH 2 D towards Oph-F, and the maps of N 2 H + and N 2 D + towards Oph-H-MM1.The column densities are ∼ 10 12 cm −2 in N 2 H + , N 2 D + , ∼ 10 13 cm −2 in -NH 2 D, and ∼ 10 14 cm −2 in NH 3 for all cores.In H 13 CO + and DCO + , the column densities are ∼ 10 11−12 cm −2 for Oph-E-MM2 and Oph-C-N.
(iii) We presented the deuterium fraction maps for three species:   (N 2 D + /N 2 H + ) for four cores,   (NH 2 D/NH 3 ) for three cores, and   (DCO + /HCO + ) for two cores.While the observed   (DCO + /HCO + ) is similar to deuterium fraction previously observed in other cold cores in C-bearing species (0.014-0.090), the observed   (N 2 D + /N 2 H + ) and   (NH 2 D/NH 3 ) are very high (0.26-0.81).The high deuterium fraction may be partly due to our line analysis method (in case the excitation temperature of the hydrogenated species was underestimated), but still the studied cores show some of the highest deuterium fractions among cold low-mass cores in the literature.
(iv) Deuterium fraction maps   (DCO + /HCO + ),   (N 2 D + /N 2 H + ) and   (NH 2 D/NH 3 ) show a significant, a factor of 2-10 (median difference is factor of 4.2) difference in the deuterium fraction values for nitrogen-bearing and carbon-bearing species, which confirms the fact that deuterium fraction increases with gas density.
(v) We analysed the correlations between the deuterium fractions and the physical parameters of the cores -H 2 column density, gas and dust temperature and the level of turbulence, represented by the non-thermal component of the velocity dispersion.We found no significant correlation with any of the parameters, probably because of the small temperature range covered (≃2 K) and small area of the   maps which do not reach the point of transition to coherence/turbulence.
(vi) We estimated the upper limit on the CO depletion factor (in all cores) and the lower limit on the ionization degree and cosmic ray ionization rate (in Oph-E-MM2).Both estimates reminded us about a well known paradox of undepleted CO in the dense cores of L1688.The observations of rare CO isotopologues, such as C 17 O or C 18 O both towards the dense cores and undepleted outskirts of the L1688 cloud are needed to measure the CO abundance in the core and study the correlation between deuterium fraction and CO depletion.
(vii) We run a model of chemical evolution of the cold cores (pyRate3; Sipilä et al. 2015a) trying to reproduce our observation data and constrain the model parameters.We analysed four models with two different gas density profiles and two initial nitrogen elemental abundances, 2.14×10 −5 and 7.60×10 −5 , comparing the observed and modelled column densities.We find that the variation in the density profile does not affect the model result significantly.The initial N-abundance, on the contrary, significantly affects the column densities of the N-bearing species and slightly affects the column densities of the C-bearing species.Taking into account the accuracy of the model, in all four variants the modeled column densities agree with the observed ones (within an order of magnitude).However, when it comes to the deuterium fraction, the modelled   (DCO + /HCO + ) is much higher than the observed one; the modelled   (N 2 D + /N 2 H + ) is underestimated but close to the observed ones; the modelled   (NH 2 D/NH 3 ) is significantly underestimated compared to the observed ones.The best modelling result for   (N 2 D + /N 2 H + ) is obtained with (N) = 2.14×10 −5 and the density profiles based on the Herschel data (profile 1).The best result for   (DCO + /HCO + ) was obtained with (N) = 7.6×10 −5 and the analytical profiles (profile 2).With all four models, NH 3 column density is overproduced.
(viii) We estimated the chemical ages of the dense cores from the deuterium fraction peaks, the median ages are 2.65×10 5 , 2.65×10 5 , 0.82×10 5 and 1.80×10 5 years for Oph-C-N, Oph-E-MM2, Oph-F and Oph-H-MM1 respectively.The obtained ages are very similar to each other.The chemical ages are not related to dynamical ages, due to the assumed static physical structure.
The spatial distribution of deuterium fraction measured for the three pairs of species in dense cores is presented for the first time and can serve as an observational benchmark for chemical models.Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.

Figure B4.
Velocity dispersion maps for all studied lines.The contours show the molecular hydrogen column density.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .The beam size is shown in the bottom left corner of each map.The stars show the positions of the YSOs (YLW15, Class 0+I, CRBR 2422.8-3423,Class I; Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.

Figure 3 .
Figure 3. Histogram of the non-thermal component of the velocity dispersion  NT of deuterated species (left) and for all species (right) for all cores.The step is 0.011 km s −1 .The black line shows the thermal velocity dispersion   at 12 K.

Figure 4 .
Figure 4.The column density maps of N 2 H + , N 2 D + , -NH 2 D, NH 3 , DCO + , and H 13 CO + towards Oph-C-N.The molecular hydrogen column density is shown by black contours.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .The beam size is shown in the bottom left corner of each map.The white cross shows the position observed in Punanova et al. (2016).

Figure 5 .
Figure 5.The column density maps of N 2 H + , N 2 D + , -NH 2 D, NH 3 , DCO + , and H 13 CO + towards Oph-E-MM2.The molecular hydrogen column density is shown by black contours.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .The beam size is shown in the bottom left corner of each map.The white cross shows the position observed in Punanova et al. (2016).

Figure 6 .
Figure 6.The column density maps of N 2 H + , N 2 D + towards Oph-H-MM1 and the column density maps of N 2 H + , N 2 D + , -NH 2 D, and NH 3 towards Oph-F.The molecular hydrogen column density is shown by black contours.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .The beam size is shown in the bottom left corner of each map.The white crosses show the positions observed in Punanova et al. (2016).The stars show the positions of the YSOs (YLW15, Class 0+I, CRBR 2422.8-3423,Class I;Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.

Figure 7 .
Figure 7.The deuterium fraction maps towards Oph-C-N, Oph-E-MM2, Oph-F and Oph-H-MM1.The molecular hydrogen column density is shown by black contours.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .Top: Oph-C-N; centre: Oph-E-MM2; bottom: Oph-F and Oph-H-MM1.The beam size is shown in the bottom left corner of each map.The crosses show the positions of the single pointing observations in Punanova et al. (2016).The stars show the positions of the YSOs (YLW15, Class 0+I, CRBR 2422.8-3423,Class I;Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.

Figure 9 .
Figure 9. Model profiles of gas density and dust temperature: profile 1 is based on the  (H 2 ) and  dust from Ladjelate et al. (2020) with a central density and temperature from Pattle et al. (2015); profile 2 is based on the analytical profile from Tafalla et al. (2002) and (H 2 ) from Pattle et al. (2015),  dust is constant, from Pattle et al. (2015).

Figure A1 .
Figure A1.Parameter space exploration of the optical depth and the excitation temperature of the observed lines performed with the Monte-Carlo method for Oph-C-N, Oph-E-MM2, Oph-F, and Oph-H-MM1.Abscissa and ordinate show the excitation temperature and the optical depth, respectively.The colour scale shows the uncertainty zones of the normal distribution.

Figure A2 .
Figure A2.The example spectra of the observed lines towards Oph-C-N, Oph-E-MM2, Oph-F, Oph-H-MM1, Oph-I-MM2 and Oph-H-MM2.The black lines show the observed spectrum, the red lines show the synthetic spectra obtained from the hyperfine fitting.

Figure B1 .
Figure B1.Integrated intensities of all studied lines.The contours show the molecular hydrogen column density.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .The beam size is shown in the bottom left corner of each map.The stars show the positions of the YSOs (YLW15, Class 0+I, CRBR 2422.8-3423,Class I;Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.

Figure B2 .
Figure B2.Optical depth maps for all studied lines.The contours show the molecular hydrogen column density.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .The beam size is shown in the bottom left corner of each map.The stars show the positions of the YSOs (YLW15, Class 0+I, CRBR 2422.8-3423,Class I;Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.

Figure B3 .
Figure B3.Radial velocity maps for all studied lines.The contours show the molecular hydrogen column density.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .The beam size is shown in the bottom left corner of each map.The stars show the positions of the YSOs (YLW15, Class 0+I, CRBR 2422.8-3423,Class I;Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.

Figure B6 .
Figure B6.The upper limits of the CO depletion factor for the observed cold dense cores.The contours show the molecular hydrogen column density.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .The beam size is shown in the bottom left corner of each map.The white crosses show the positions observed inPunanova et al. (2016).The stars show the positions of the YSOs (YLW15, Class 0+I, CRBR 2422.8-3423,Class I;Young et al. 1986;Kirk et al. 2017;Comeron et al. 1993;Bontemps et al. 2001) in Oph-F.

Figure B7 .
Figure B7.Ionization degree  () found as abundance of all observed ions (left), ionization degree and ionization rate  based of the observations and analytical model (centre and right).The contours show the molecular hydrogen column density.The first contour starts at 1.5×10 22 cm −2 with a contour step of 0.5×10 22 cm −2 .The beam size is shown in the bottom left corner of each map.The white crosses show the positions observed in Punanova et al. (2016).

Figure B8 .
Figure B8.The model column densities and deuterium fractions for Oph-C-N.The simulation was performed for two different nitrogen abundances (N)=7.60×10−5 and (N)=2.14×10−5 and two gas density and dust temperature profiles.Profile 1 is based on observations (Ladjelate et al. 2020), profile 2 is based on the analytical model from Tafalla et al. (2002).The column densities and deuterium fractions are convolved with the 33.6 ′′ beam.First (top) row: temporal evolution of the column densities towards the core centre.Second row: temporal evolution of the deuterium fractions.Third row: the column density profiles at the age indicated on top of each column.Fourth (bottom) row: the deuterium fractions at the age indicated on top of each column.The ages for the profiles are shown by gray vertical lines.

Figure B9 .
Figure B9.The model column densities and deuterium fractions for Oph-E-MM2.The simulation was performed for two different nitrogen abundances (N)=7.60×10−5 and (N)=2.14×10−5 and two gas density and dust temperature profiles.Profile 1 is based on observations (Ladjelate et al. 2020), profile 2 is based on the analytical model from Tafalla et al. (2002).The column densities and deuterium fractions are convolved with the 33.6 ′′ beam.First (top) row: temporal evolution of the column densities towards the core centre.Second row: temporal evolution of the deuterium fractions.Third row: the column density profiles at the age indicated on top of each column.Fourth (bottom) row: the deuterium fractions at the age indicated on top of each column.The ages for the profiles are shown by gray vertical lines.

Figure B10 .
Figure B10.The model column densities and deuterium fractions for Oph-F.The simulation was performed for two different nitrogen abundances (N)=7.60×10−5 and (N)=2.14×10−5 and two gas density and dust temperature profiles.Profile 1 is based on observations (Ladjelate et al. 2020), profile 2 is based on the analytical model from Tafalla et al. (2002).The column densities and deuterium fractions are convolved with the 33.6 ′′ beam.First (top) row: temporal evolution of the column densities towards the core centre.Second row: temporal evolution of the deuterium fractions.Third row: the column density profiles at the age indicated on top of each column.Fourth (bottom) row: the deuterium fractions at the age indicated on top of each column.The ages for the profiles are shown by gray vertical lines.

Figure B11 .
Figure B11.The model column densities and deuterium fractions for Oph-H-MM1.The simulation was performed for two different nitrogen abundances (N)=7.60×10−5 and (N)=2.14×10−5 and two gas density and dust temperature profiles.Profile 1 is based on observations (Ladjelate et al. 2020), profile 2 is based on the analytical model from Tafalla et al. (2002).The column densities and deuterium fractions are convolved with the 33.6 ′′ beam.First (top) row: temporal evolution of the column densities towards the core centre.Second row: temporal evolution of the deuterium fractions.Third row: the column density profiles at the age indicated on top of each column.Fourth (bottom) row: the deuterium fractions at the age indicated on top of each column.The ages for the profiles are shown by gray vertical lines.

Table 1 .
The (Di Francesco et al. 2008mission peaks of the studied cold dense cores(Di Francesco et al. 2008) and of the ammonia peaks (with asterisk*,

Table 2 .
Note:   obs is the native half-power beam width at the given frequency.The transition frequencies, the frequencies and relative intensities of the hyperfine structure components are taken from the following works:  from Pagani et al. (2009b),  from Daniel et al. (2016),  from Schmid-Burgk et al. (2004),  from Lattanzi et al. (2007a).The frequencies and relative intensities of the hyperfine structure components are also available via CDMS (Endres et al.The observed lines, observational parameters and the excitation temperature used for the fits.

Table 5 .
The comparison of our column densities and deuterium fractions (left part of the table) with those from Punanova et al. (2016, right part of the table).
Punanova et al. (2016)values of F-MM2 inPunanova et al. (2016).*The  ex and  are separated with a slash for N 2 H + and N 2 D + .