EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 700 (August 2025) A&A, 700 (2025) A141 Full HTML
Open Access
Issue
A&A
Volume 700, August 2025
Article Number A141
Number of page(s) 14
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202555336
Published online 13 August 2025

© The Authors 2025

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.

Open Access funding provided by Max Planck Society.

1 Introduction

Ammonia and its deuterated isotopologues are frequently used as probes of physical and chemical conditions in dense interstellar clouds through their rotation-inversion spectra at centimetre and millimetre wavelengths. These molecules contain two or three identical protons (fermions) or deuterons (bosons), and the symmetrisation postulate of quantum mechanics imposes certain restrictions on their complete internal wave functions, which describe the quantum state of an isolated molecule, including the electronic and rotational-vibrational states as well as the electronic and nuclear spins. For example, the complete internal wave function of NH2D must change sign under the permutation of the two protons, while that of NHD2 must remain unchanged by the interchange of the two deuterons. NH2D and NHD2 have both symmetric and anti-symmetric nuclear spin wave functions, and to comply with the symmetrisation postulate, these can only be associated with certain rotation-inversion states. The restrictions described above give rise to different nuclear spin modifications of NH3, NH2D, NHD2, and ND3, which can be treated as different chemical species. It is customary to use the Greek appellation ‘ortho’ for the modification with the highest nuclear spin statistical weight (number of symmetrised nuclear spin functions that have appropriate symmetry) and ‘para’ for the modification with the lowest weight. ND3 has a third symmetry species called ‘meta’ (e.g. Bunker & Jensen 2005; Hugo et al. 2009; Daniel et al. 2016).

Chemical models that distinguish between different spin species indicate that their relative abundances do not necessarily correspond to the nuclear spin statistical weights. Deviations from the statistical weights originate in chemical reactions that form ammonia if they proceed through an intermediate reaction complex where protons and deuterons can be completely mixed. The possible spin states of the reaction complex and the dissociation products are determined by selection rules that derive from the conservation of spin angular momentum and permutation symmetry (Hugo et al. 2009; Sipilä et al. 2015; Schmiedt et al. 2016; Hily-Blant et al. 2018). A practical reason to have a priori knowledge of the spin ratios is that it is not always easy to determine the abundances of all the spin modifications of a molecule and thus its total abundance directly from observations. For example, the fundamental transitions of ortho-NH3 and ortho-ND3 are at unfavourable frequencies (572 and 615 GHz, respectively) for ground-based observations. The former line has been observed in the pre-stellar core L1544 with Herschel but is difficult to interpret due to self-absorption (Caselli et al. 2017). Deviations from the statistical ortho/para (o/p) ratios of water and ammonia (the latter derived from NH2) have been reported in comets (Mumma & Charnley 2011; Shinnaka et al. 2016). These were previously thought to reflect the temperatures at which the ices in the comet nucleus were formed or condensed. Currently, the idea that possible deviations from statistical ratios are caused by gas-phase reactions in comae is perhaps more favoured (Faure et al. 2013; Hama et al. 2018; Faure et al. 2019). Experimental and theoretical studies suggest that water desorbing from cometary nuclei should have a statistical o/p ratio of 3 (Hama et al. 2018; Cheng et al. 2022 and references therein). A couple of recent observational estimates for cometary comae find no significant deviation from this value when optical thickness effects are carefully taken into account (Cheng et al. 2022; Woodward et al. 2025).

We have previously determined the o/p ratios of NH2D and NHD2 in two pre-stellar cores in Ophiuchus, H-MM1 and Oph D, using the Large APEX sub-Millimeter Array (LAsMA) multi-beam receiver of the Atacama Pathfinder EXperiment (APEX) telescope (Harju et al. 2024). These observations show, with a higher accuracy than previously achieved, that the mentioned spin ratios are within 20% (taking 3 σ limits into account) of their high-temperature statistical values, 3 and 2, respectively. The results were explained by a chemistry model that assumes the proton hop (PH) mechanism (see e.g. Oka 2004) in proton transfer and hydrogen abstraction reactions that contribute to the formation of ammonia (Sipilä et al. 2019a). The alternative to PH is the so-called full scrambling (FS) mechanism, in which chemical reactions proceed via the formation of a relatively long-lived intermediate complex, where the light H and D nuclei can be mixed. According to our tests with both the PH and FS scenarios, the efficiency of ammonia desorption from grains does not affect its spin ratios in the gas. The insensitivity to the desorption rate was explained by a rapid reprocessing of released material in ion-molecule reactions, and this was considered to rule out the possibility that statistical spin ratios are caused by formation on grains. In contrast to PH, the FS model predicts clear deviations from the statistical spin ratios. In FS simulations, the o/p-NH2D ratio generally decreases well below 3 or even below 2 in the gas, and o/p-NHD2 always stays clearly above 2 (Sipilä et al. 2015; Hily-Blant et al. 2018).

The Ophiuchus molecular cloud complex, where H-MM1 and Oph D are located, is exposed to strong radiation from the Sco OB2 association and local B stars (e.g. Ladjelate et al. 2020; Miret-Roig et al. 2022), and it is possible that the formation history and external conditions of the two cores have affected their composition in a way that current chemistry models cannot account for. Therefore, it seems helpful to repeat the measurement of the spin ratios in a quiescent environment. To this end, we present observations of the ortho and para modifications of NH2D and NHD2 towards the pre-stellar core L1544, which lies on the outskirts of the Taurus molecular cloud complex, far from strong sources of radiation and stellar winds. The core is characterised by a high central density (n(H2) > 106 cm−3) and a low central temperature (~7 K; e.g., Crapsi et al. 2007; Caselli et al. 2019; Caselli et al. 2022). An additional motivation for observing L1544 is that the column densities of ortho- and para-NH2D (oNH2D and pNH2D) derived in the spectroscopic survey of Melosso et al. (2021) indicate an o/p-NH2D ratio of 2.31 ± 0.11, which would agree with the predictions of FS models. This o/p ratio was derived from the ortho- and para-lines at 85.9 GHz and 110.2 GHz (λ ~ 3.5 and 2.7 mm) and includes uncertainty caused by the different beam sizes and the different Einstein coefficients and critical densities of the two transitions. To confirm the low o/p-NH2D ratio, and to also determine the o/p-NHD2 ratio in this object, we observed the ground-state lines of both species at λ ~ 0.9 mm with APEX/LAsMA, and combined these with lines of the same molecules around λ = 3 mm, obtained with the 30 m telescope of the Institut de radioastronomie millimétrique (IRAM).

thumbnail Fig. 1

Footprint of the APEX/LAsMA array on the 8 μm surface brightness images of L1544 observed by the Spitzer satellite. The position of the central pixel is RA 05h04m17⋅s1$\[05^{\mathrm{h}} 04^{\mathrm{m}} 17^\text{s}_\cdot1\]$, Dec +25°10′43″ (J2000.0). The circle diameter corresponds to the HPBW of APEX at 334 GHz. The core appears as a shadow against the weak mid-infrared background. The bright spots are foreground stars.

2 Observations

2.1 APEX

The ground-state rotational lines of oNH2D, pNH2D, ortho-NHD2 (oNHD2), and para-NHD2 (pNHD2) near 333 and 335 GHz were observed towards L1544 with the LAsMA array (Güsten et al. 2008) on the APEX telescope (Güsten et al. 2006). The observations were made in August 2023, under project number M-0112.F-9506A-2023. These were single-point position switching observations with the central pixel of LAsMA pointed at the column density peak of L1544 (RA 05h04m17⋅s1$\[05^{\mathrm{h}} 04^{\mathrm{m}} 17^\text{s}_\cdot1\]$, Dec. +25°10′43″ (J2000.0)). The reference position (OFF) was selected 10′ to the equatorial east of the ON position. The array has 7 pixels separated by approximately 40″. Figure 1 shows the footprint of the instrument superposed on the 8 μm map of the core measured by the InfraRed Array Camera (IRAC) of the Spitzer Space Telescope. The beam size of APEX is 18⋅′′7$\[18^{\prime \prime}_\cdot 7\]$ (full width at half maximum) at 333 GHz. The 12 m APEX telescope is located at an altitude of 5100 m at Llano de Chajnantor in Chile.

The backend was composed of fast Fourier transform spectrometers (FFTSs; called the FFTS4G) that covered two 4 GHz intermediate frequency bands resolved into 65536 spectral channels of width 61.03 kHz (corresponding to ~55 m s−1). Lines of ortho- and para-NH2D and NHD2 were observed in the lower sideband of the receiver, which was centred at 334.145 GHz. The frequencies of the observed transitions are listed in Table 1. The energies of the upper transition level and the critical densities, ncrit, of the transitions are also listed. Spectroscopic data were taken from Melosso et al. (2021). The (de)excitation rate coefficients for NH2D and NHD2 in collisions with H2 were adopted from Daniel et al. (2014) and Daniel et al. (2016).

The critical densities in Table 1 were estimated at 10 K, using the optically thin formula presented by Shirley (2015, Eq. (4)). This formula takes into account all collisionally induced transitions starting from the upper transition level of the emission line. In addition to the downward collisions to the lower transition level, it is primarily the collisions 101 → 111 and 111 − 000 that affect the critical densities of the 101 → 000 and 111 → 101 emission lines of oNH2D and pNH2D. For oNHD2 and pNHD2, the collisions 111 → 101 compete with collisions 111 → 000, which correspond to the observed ground-state lines.

Table 1

Observed transitions.

2.2 IRAM 30 m

The L1544 core was mapped in the 85.9 GHz and 110.2 GHz ortho- and para-NH2D lines with the IRAM 30 m telescope in the on-the-fly mode. The mapping was performed in 2015. The scan directions alternated between the east-west equatorial and the north-south direction. The separation between rows and columns of the 85.9 GHz map was in general 5″; the half-power beam width (HPBW) is 28⋅′′6$\[28^{\prime \prime}_\cdot 6\]$ at this frequency, so the separation is ~0.17 × HPBW. The area covered was 120″ × 120″. In the 110.2 GHz map, the rows and columns were separated by 8″ (~0.36 × HPBW; the HPBW is 22⋅′′3$\[22^{\prime \prime}_\cdot 3\]$), and the map covered an area of ~160″ × 160″. Both maps were centred on the column density peak of the core, RA 05h04m17⋅′′2$\[05^{\mathrm{h}} 04^{\mathrm{m}} 17^{\prime\prime}_\cdot2\]$, Dec. +25°10′43″ (J2000.0).

For imaging, the maps were interpolated and averaged to regular grids using the cell sizes 5″ and 8″ for the 85.9 GHz and 110.2 GHz maps, respectively. As gridding function, we used a Bessel function of the first kind, J1, tapered with a Gaussian (Mangum et al. 2007; Sawada et al. 2008). In practice, the observation representing a certain grid position was obtained by accumulating all observations with a radius of one HPBW from this position, weighted by the function J1(πr/a)/(πr/a) × exp(−(r/b)2), where r is the distance from the grid position in units of 1/3 of the beam width, and the constants are a = 1.55, b = 2.52. This gridding function causes a small broadening of the effective beam. The effect is approximately 10% for the 110.2 GHz map, but negligible for the 85.9 GHz map, where the cell size is small compared to the beam (see Sawada et al. 2008, their Fig. 6). In addition to the gridding function, the accumulated spectra were weighted by the inverse square of the system noise temperature TSYS (the integration time was the same for all samples).

The observations were performed with the IRAM Eight Mixer Receiver (EMIR) operating at millimetre wavelengths. The two NH2D transitions required different spectral setups of the E90 channel of the receiver; the line at 85.9 GHz was observed in the lower-outer subband, whereas the 110.2 GHz line was placed in the upper-outer subband. The EMIR receiver was connected to FFTSs. The FFTS back-end was configured to give 7.28 GHz of instantaneous bandwidth and a channel width of 50 kHz. The FFTS broadband spectra also covered the 110 − 101 lines of ortho- and para-NHD2 near 111 GHZ. These lines were only marginally detected and were not useful in determining the o/p ratio.

2.3 ALMA

We included in the analysis the para-NH2D(111 − 101) map of L1544 that was observed in 2016 with the Atacama Large (Sub)millimeter Array (ALMA) and presented in Caselli et al. (2022). The observations consisted of a three-point mosaic using the 12 m array, which is the main array, and a single pointing with the 7 m array, the Atacama Compact Array (ACA). The project ID was 2016.1.00240.S (PI Caselli). Here, we used an image cube that was produced by using natural weighting of the combined (12 m array and ACA) visibilities, and the CLEAN algorithm in the deconvolution of the dirty image. The final image was obtained by convolving the CLEAN component image with a Gaussian beam of 2⋅′′3×1⋅′′6$\[2^{\prime \prime}_\cdot 3 \times 1^{\prime \prime}_\cdot 6\]$, which corresponds to the main lobe of the synthesised beam. The rms noise in the central parts of the image is 0.15 K per 42 m s−1 channel. The reduction procedure is described in Sect. 2.1 of Caselli et al. (2022). The extent of pNH2D emission was found to be approximately 1′×0⋅′5$\[1^{\prime} \times 0^{\prime}_\cdot5\]$ (see Fig. A.1 in Caselli et al. 2022, and Fig. B.1 below).

2.4 Observed spectra and maps

The NH2D and NHD2 spectra observed with APEX and IRAM towards the H2 column density maximum of L 1544 are shown in Figs. 2 and 3. The APEX spectra correspond to pixel ① of the LAsMA array. The integrated intensity maps of the 85.9 GHz and 110.2 GHz lines of NH2D observed with the 30 m telescope are shown in Fig. 4. Spectra and maps are presented on the antenna temperature, TA∗$\[T_{\mathrm{A}}^{*}\]$, scale.

The line parameters and column densities derived from Gaussian fits to the hyperfine structure of the lines observed towards the column density peak of L1544 are presented in Table 2. The observed spectra were converted to the main-beam brightness temperature scale, TMB, using the beam efficiencies provided by the telescope teams. In column density estimates, we assumed uniform beam filling, line-of-sight homogeneity, and that the excitation temperature, Tex, is constant for all rotational transitions of the molecule. Furthermore, the calculation involved the assumption that populations of hyperfine states at a certain rotational level are proportional to their statistical weights. The strong quadrupole interaction of the N nucleus with the electric field of the electrons gives rise to five distinct components of the JKa, Kc = 111 − 101 lines of NH2D, and three components of the ground-state lines (101 − 000 and 111 − 000) of NH2D and NHD2. These are further split into several unresolved components owing to quadrupole interaction of the deuteron, and weaker spin-rotation and spin-spin couplings (Coudert & Roueff 2006; Melosso et al. 2021). The line width listed in Table 2 refers to the intrinsic width of an individual hyperfine component. The derived local standard of rest velocities of the 111 − 000 lines of oNHD2 and pNHD2 differ by approximately 210 m s−1, whereas velocity difference of the oNH2D and pNH2D lines in the same band is only 43 m s−1. This suggests that the rest frequencies of the NHD2 lines are not as accurate as those of the NH2D lines. The velocity difference 210 m s−1 corresponds to ~240 kHz.

The strong rotation-inversion lines of NH2D at λ = 3 mm with five hyperfine groups have very high signal-to-noise ratios, and the formal errors of the column densities of NH2D derived from these lines are clearly smaller than the estimates from the λ = 0.9 mm lines. Although the column densities of pNH2D from the 111 − 101 and 101 − 000 lines agree, the total oNH2D column densities derived from the two ortho-lines differ by approximately 15%. This is due in part to the different beam sizes for the two oNH2D lines (29″ vs 19″). For the two pNH2D lines the difference is smaller (22″ vs 19″). The o/p-NH2D ratio from the IRAM observations is 2.7 ± 0.2, whereas from the APEX spectra we obtain o/p-NH2D = 2.3 ± 0.5. The hyperfine structure fit to the pNHD2(111 − 000) line does not give reliable estimates for the Tex and τ, and we estimated the column density using the optically thin approximation, adopting the Tex derived for the ortho-line. The resulting o/p-NHD2 ratio is 2.3 ± 1.2. Depending on the distribution of the gas density and molecular abundances, the effective source size can be different for each observed line. This causes an uncertainty in the beam-source coupling efficiency, which propagates to Tex. In addition to the uncertainties involved in the determination of Tex and τ, especially when the line only has three resolved components, the assumption of a constant Tex along the line of sight in the presence of large density gradients is an additional source of uncertainty. In what follows, we estimate the column densities and the fractional abundances of the observed molecules applying radiative transfer calculations to a spheroidal model of the L1544 core.

thumbnail Fig. 2

NH2D spectra observed with IRAM and APEX towards the column density peak of L1544, which is the position of pixel ① in Fig. 1. The IRAM spectra are shifted up by 0.8 K for clarity.
thumbnail Fig. 3

NHD2 spectra observed with APEX towards the column density peak of L1544 (pixel ① in Fig. 1). The pNHD2 spectrum is shifted up by 0.3 K for display purposes.

3 Core model

We adopted the density and gas and dust temperature distributions from the model used for L1544 by Sipilä et al. (2019b) and Riedel et al. (2023), which, in turn, was developed from a model derived by Keto & Caselli (2010, see Fig. 2 of Riedel et al. 2023). However, because the core is clearly elongated, we deformed the originally spherical model to a prolate spheroid. A Gaussian fit to the Spitzer 8 μm image of the core gave an axis ratio of 2.0 and a long axis tilt angle of −36⋅∘1$\[-36^\circ_\cdot1\]$ from the declination axis. For the purpose of radiative transfer calculations, we constructed a cubical grid with an edge length of ~180″ or 30 600 au. The transformation to the spheroidal model was done so that the volume within an isosurface (a surface representing constant values of the density, temperature, etc.) is the same as that of the corresponding isosurface of the spherical model. Consequently, the gas properties at point (x, y, z) in the grid, measured from the origin of the principal axis system, correspond to those of the spherical model at radius r, which is obtained from r = 21/3{x2 + y2 + z2/γ2}1/2, where γ ≡ c/a is the axis ratio of the spheroid (here γ = 2).

We first tested how well the observed spectra can be reproduced by assuming that the fractional abundances of NH2D and NHD2 are constant along the line of sight. Line simulations were performed using the Line radiative transfer with the OpenCL (LOC) program (Juvela 2020). The hyperfine structure of the lines was modelled by assuming local thermodynamic equilibrium between the components, that is, the effects of line overlap were not taken into account. The fractional abundances were determined by comparing the simulated spectra with the observed ones; the method that applies minimum chi-square estimation is described in Sect. 3.2 of Harju et al. (2024). Only the statistical uncertainties based on the noise of the observed spectra are taken into account in this method. The best-fit abundances and their 1 σ uncertainties obtained in the centre of the map (0,0) and the second strongest position (−23″, 35″) are listed in Table 3. In the latter position, the NHD2 lines are weak, and only the NH2D abundances are given. The NH2D and NHD2 spectra observed in the seven LAsMA positions, including (0,0) and (−23″, 35″), are shown in Figs. A.1 and A.2. The derived fractional NH2D abundance is clearly higher towards the offset position than towards the centre of the core. We note that the N(H2) estimate used for the fractional abundance comes from the spheroidal core model, which may miss a local condensation near position (−23″, 35″). The other alternative is that the assumption of constant abundances throughout the core is not true, which would also render the assumption of line-of-sight constancy illogical. In what follows, we apply a chemistry model to predict the abundance distributions in the core and simulate spectra from the obtained models.

thumbnail Fig. 4

Integrated JKa, Kb = 111 − 101 line intensity maps of oNH2D and pNH2D towards L1544 observed with the IRAM 30 m telescope. The intensity is on the TA∗$\[T_{\mathrm{A}}^{*}\]$ scale. The beam sizes at the line frequencies (85.9 GHz and 110.2 GHz) are indicated in the bottom right.
Table 2

Line parameters and column densities from Gaussian fits to the hyperfine structure.

4 Chemistry model

For chemical simulations, we used the gas-grain chemical model pyRate described in Sipilä et al. (2019a), Riedel et al. (2023), and Harju et al. (2024), employing the PH and FS scenarios for proton donation and abstraction reactions. In the PH network, proton transfer reactions involving O, C, and N are assumed to proceed via the exchange of one proton (or deuteron) between the reactants, while in the FS network multiple atom exchanges are possible. However, the H3++H2$\[\mathrm{H}_{3}^{+}+\mathrm{H}_{2}\]$ reaction system and its deuterated variants form an exception in the PH network. For the reactions that involve only the H or D nuclei, the rate coefficients and branching ratios adopted from Hugo et al. (2009), follow the FS scenario in both networks. Some deuterated variants of the H3++H2$\[\mathrm{H}_{3}^{+}+\mathrm{H}_{2}\]$ system were recently studied by Jiménez-Redondo et al. (2024) through experiments and simulations. The results indicate that the growth of the abundances of D2H+, D2H+, and D3+$\[\mathrm{D}_{3}^{+}\]$ is suppressed more than previously thought in warm conditions by hydrogenation reactions. Since we are dealing with very cold gas in L1544, we did not make any changes to the original reaction set of Hugo et al. (2009).

The efficiency of chemical desorption was assumed to be constant, with 1% exothermic surface reactions leading to desorption. The dust grains were assumed to have uniform size characterised by the grain radius a and the density of the material ρ, assumed to be 2.5 g cm−3. The simulations were performed using the grain radii a = 0.1, 0.2, and 0.3 μm. The value a = 0.1 μm corresponds to a grain surface area of 5.0 × 10−22 cm−2 per H atom, with the adopted material density and the dust-togas mass ratio (0.01; Sipilä et al. 2020). The grain surface area is inversely proportional to the grain radius for particles of uniform size. We also varied the cosmic-ray ionisation rate, testing the values 0.5, 1, 2, 4, and 8 times 1.3 × 10−17 s−1.

We assumed that the gas has been processed under low-density ‘dark cloud’ conditions (n(H2) = 104 cm−3, T = 10 K, AV = 5mag) for half a million years before the ‘dense core’ phase. The initial abundances used for the first simulation step are the same as those listed in Table 1 of Riedel et al. (2023). The chemical abundances attained at the end of the dark cloud phase were in turn used as initial values of the dense core phase. The abundances in the dense core were calculated for concentric spherical shells using the one-dimensional model described in Riedel et al. (2023). These were interpolated to the spheroidal model for spectral line simulations using the transformation described above.

The evolution of the average fractional abundances of NH2D and NHD2 in the FS and PH models is shown in Fig. 5, with the ortho- and para-species separated. Here, the grain radius is a = 0.2 μm, and the cosmic-ray ionisation rate is ζH2 = 2.6 × 10−17 s−1. The curves represent the column density ratios N(mol)/N(H2) towards the centre of the core. These column densities are derived from the spheroidal models and convolved with the telescope beam. The ratios are in this sense comparable to the estimates presented in Table 3 for the (0,0) position. In Fig. 6 we show the o/p ratios of NH3, NH2D, and NHD2, and the fractionation ratios NH2D/NH3, NHD2/NH2D, and ND3/NHD2 as functions of time predicted by the same PH and FS models as shown in Fig. 5. The PH model predicts constant o/p ratios corresponding to their high-temperature statistical values, whereas in the FS model the ratios show some time variation. One can see that the predicted values of the o/p-NH2D and o/p-NHD2 ratios are almost swapped between the PH and FS models. In agreement with the simulations of Faure et al. (2013) and Le Gal et al. (2014), the o/p ratio of NH3 remains below 1 in the FS model. The deuterium fractionation ratios will be discussed briefly in Sect. 7.2.

thumbnail Fig. 5

Evolution of the fractional abundances of oNH2D, pNH2D, oNHD2, and pNHD2 in the core model. The abundances correspond to the column density ratios N(mol)/N(H2) towards the centre of the core, smoothed to the resolution of the 18⋅′′7$\[18^{\prime \prime}_\cdot 7\]$ APEX beam at 333 GHz. Predictions from the PH scenario in the chemical network (a) and the FS scenario (b). In this simulation, a = 0.2 μm and ζH2 = 2.6 × 10−17 s−1. The horizontal dashed lines and thick semi-transparent lines represent the fractional abundances and their formal 1σ errors derived towards (0,0) in Sect. 3. These abundances are listed in Table 3.
Table 3

Fractional abundances, X, and abundance ratios in L1544 derived from the LAsMA spectra in pixels ① (0,0) and ④ (−23″, 35″) assuming that X is constant throughout the core.

thumbnail Fig. 6

Evolution of the o/p-NH3, o/p-NH2D, and o/p-NHD2 ratios and the fractionation ratios NH2D/NH3, NHD2/NH2D, and ND3/NHD2 according to the PH (a) and FS (b) models visualised in Fig. 5. The fractional abundances used for the ratios are calculated in the same way as those shown in Fig. 5. The horizontal straight dashed lines and semi-transparent thick lines represent the abundance ratios derived towards (0,0) from APEX observations (Sect. 3). Observational data for NH3 and ND3 are not available at the same resolution.

5 Simulated spectra

We first attempted to reproduce the spectra observed with APEX and IRAM using the abundance distributions predicted by the chemistry models described above. Almost all models tested here can reasonably well reproduce the two pNH2D lines and the oNH2D(101 − 000) line at 333 GHz at one simulation time. Examples of simulated NH2D spectra towards the centre of the core are shown in Fig. 7. These are the best-fit spectra from the FS and PH models with certain combinations of a and ζH2. For models assuming a = 0.2 μm and ζH2 = 2.6 × 10−17 s−1 (top panel), the best-fit simulation time is 3.2 × 105 yr. Like for most other models, the 86 GHz line of the oNH2D line is clearly brighter than observed whenever the three other lines agree. The discrepancy depends mainly on the relatively large telescope beam (29″ vs 19″−22″) in the 86 GHz observations. It suggests either that the predicted NH2D distribution is overly extended or that the model underpredicts the abundances near the centre of the core. In both cases, the lines simulated with smaller beams would be excessively weak with respect to the simulated 86 GHz line. Taking a closer look at the spectra one notices that the PH model can better reproduce both 333 GHz lines simultaneously. In the FS model, one of these lines is either too weak or too strong when the other agrees. The same is true for ortho- and para-NHD2. This difference is also seen in Figs. 5 and 6. For the PH model, the curves showing the fractional abundances of oNH2D and pNH2D cross the horizontal lines indicating the values derived from observations approximately at the same time, and the predicted o/p ratios (constantly 3 and 2) are close to the observationally derived values. For the FS model, the best-fit times are different, and the o/p-NH2D and o/p-NHD2 ratios are almost reversed with respect to those in the PH model. The best-fit simulation times in Fig. 7 do not correspond to the times when the curves shown in Fig. 5 cross the horizontal lines because the ‘observed levels’ shown in the latter figure were derived assuming constant abundances. Examples of simulated abundance distributions are shown below.

The agreement between the simulated and observed 86 GHz lines of oNH2D improves with an increase in the cosmic-ray ionisation rate. In the bottom panel of Fig. 7, we show the spectra obtained from the models with ζH2 = 1.04 × 10−16 s−1 for the grain radius a = 0.2 μm. The agreement is also good for the model with a = 0.3 μm. At high values of ζH2 chemical reactions proceed rapidly and the 333 GHz lines reach the observed intensities already after 3–4 × 104 yr of simulation time. Also in this case, the relative intensities of the ortho- and para-lines of NH2D at 333 GHz are better reproduced by the PH model. At the best-fit time, the o/p-NH2D and o/p-NHD2 ratios predicted by the FS model are 2.3 and 2.9, respectively, while in the PH model they are always 3.0 and 2.0. Figure 8 shows the abundance profiles of NH2D and NHD2 in the core models used for the left-hand side of Fig. 7. One can see that in the model with a high cosmic-ray ionisation rate, the abundances of NH2D and NHD2 are higher near the centre of the core than in the model with a lower ζH2. In fact, the abundances of ortho- and para-NH2D are almost constant in a wide range of densities and radial distances greater than ~10″ from the centre.

The ortho- and para-NHD2 lines at 335 GHz can be fitted with models that also reproduce the 333 GHz lines of NH2D. However, as suggested in Fig. 5, NHD2 reaches the best-fit abundance earlier than NH2D. We think the difference is attributable to the static core model, which neglects the contraction of the core at the same time as the chemistry is evolving. This makes it difficult to reach agreement with observations simultaneously for two different species. The simulations of NHD2 spectra from the models that approximately reproduce all four NH2D spectra (a = 0.2 μm, ζH2 = 1.04 × 10−16 s−1) are shown in Fig. 9.

6 Simulations with an ad hoc abundance profile

The agreement between the simulated and observed spectra reached by increasing the cosmic-ray ionisation rate to ζH2 ~ 10−16 s−1 is not satisfactory because a clearly lower average value of ζH2 (~3 × 10−17 s−1) has been estimated in L1544 based on observations of molecular ions (Redaelli et al. 2021). On the other hand, the abundance profile of NH2D predicted by the high-ionisation-rate model is likely to have an element of truth. In the inner part of the core, it agrees with the ALMA observations and the modelling results of Caselli et al. (2022), which show that pNH2D becomes heavily depleted within 10″ (1700 au) from the centre of L1544 (see their Fig. 5). We ran another series of simulations assuming that the fractional abundances of ortho- and para-NH2D and NHD2 are constant in the density range nH2 = 5 × 104−5 × 105 cm−3 (corresponding to radial distances 10″ − 30″), but decrease both towards the centre and the outer parts. Here, the simulations included both the single-dish observations with APEX and IRAM and the pNH2D map observed with ALMA. We first assumed that fractional abundances are proportional to n−1 and n in these two regions. A similar distribution was used to model the observations of NH2D and NHD2 in the pre-stellar cores H-MM1 and Oph D in Ophiuchus (Harju et al. 2024), except that it was assumed that the density threshold of the region dominated by accretion is lower. The high threshold used in L1544 is based on the observations of Caselli et al. (2022). Although the APEX and IRAM spectra were not very sensitive to the slope of the density dependence in the inner part of the core, to reproduce the ALMA observations it was necessary to assume that the decrease in the NH2D abundance near the centre is steeper than n−1. We obtained rather a good agreement with all the observed spectra by assuming that the fractional NH2D and NHD2 abundances are proportional to n−2 in the inner parts of the core. The steepness of the slope does not affect the o/p ratios obtained from the APEX data.

The fractional abundances of oNH2D, pNH2D, oNHD2, and pNHD2 (in the flat region) were independently varied and the simulated spectra of the spheroidal core model were compared with observations. The best-fit abundances and their uncertainties were estimated by evaluating the χ2 statistic as explained in Harju et al. (2024). The results of these simulations are presented in Table 4. The simulated spectra from the best-fit model towards the seven positions observed with APEX/LAsMA are shown in Figs. A.1 and A.2, together with the observed spectra at these positions. The observed and simulated ALMA spectra of pNH2D along two orthogonal strips near the centre of the core are shown in Fig. B.1.

thumbnail Fig. 7

Observed (blue) and simulated (red) NH2D spectra towards (0,0). Left: simulations from the PH model with two different sets of model parameters. Right: results from the FS model with the same parameter values. The model parameters are given at the top of each panel. The upper plot in each of the four panels shows the 111 − 101 lines and the lower plot shows the 101 − 000 lines of oNH2D and pNH2D. The para-lines are shifted downwards for clarity.
thumbnail Fig. 8

Factional abundance profiles in the PH model core used to calculate the spectra shown in the left panels of Fig. 7.
thumbnail Fig. 9

Observed and simulated NHD2 spectra towards the column density peak of L1544. The red lines show simulations from the PH (top panel) and FS (bottom panel) reaction schemes. The same assumptions about the effective grain radius and the cosmic-ray ionisation rate have been used (a = 0.2 μm, ζH2 = 1.04 × 10−16 s−1) in the two models. The simulation time has been chosen to provide an approximate agreement with the pNHD2 spectrum.

7 Discussion

7.1 Distributions of NH2D and NHD2 in L1544

The efforts to simulate the observed NH2D and NHD2 lines towards L1544 have shown that it is fairly easy to construct a model that approximately reproduces the observed NH2D lines at 333 GHz and the pNH2D line at 110 GHz, and the lines of NHD2 at 335 GHz. Because the beam sizes used to observe the aforementioned lines are similar (~20″) their relative intensities do not strongly depend on the spatial extent of the molecular distribution as long as the ortho and para modifications of NH2D and NHD2 follow the same pattern. However, we encounter severe problems when including the 86 GHz line of oNH2D, which was observed with a larger beam (~29″). For most models, the 86 GHz line is too strong whenever the intensities of the other lines agree with the observations, indicating that the abundance of NH2D is excessively high in the outer parts of the core relative to that in the central parts. Because the chemistry model predicts heavy depletion of molecules in the inner parts of the core as a result of accretion onto the grains, it is natural to suspect that the depletion is too efficient in the model. One way to decrease accretion is to increase the effective grain radius, a, which reduces the total surface area of the grains. Simulations with constant grain radii of a = 0.2 and 0.3 μm gave better results than our standard assumption a = 0.1 μm. We also performed simulations with a varying grain radius according to the coagulation model of Ormel et al. (2009, with icy grains processed for 105 years), but these did not result in significant changes relative to our standard model. This is probably because in the varying grain size model the effective radius exceeds a = 0.3 μm only in a few central shells of the core model, within an angular distance of 5″ from the centre. The 86 GHz and 333 GHz lines of oNH2D could not be fitted simultaneously with any of the grain size distributions tested if the other model parameters were left unchanged.

By increasing the cosmic-ray ionisation rate, from the standard value ζH2 = 1.3 × 10−17 by a factor of 8, to ζH2 = 1.04 × 10−16 s−1, we achieved, at an early simulation time, spatial distributions of NH2D and NHD2 that reproduce all of the observed lines more or less simultaneously. An increased cosmic-ray ionisation rate speeds up gas-phase chemical reactions leading to high abundances of NH2D and NHD2 in the inner parts of the core before depletion takes effect. The resulting distribution of NH2D is fairly flat at the beginning, except for a drop near the very centre (Fig. 8 right). The shape of the distribution in the inner parts of the core agrees with the ALMA observations and modelling results of Caselli et al. (2022), which show that pNH2D becomes heavily depleted within 10″ (1700 au) from the centre of L1544 (see their Fig. 5). As discussed in this paper, the detection of such a small, almost complete freeze-out zone requires interferometric observations. That the present chemical model seems to overpredict molecular depletion and that single-dish observations can be modelled assuming constant fractional abundances have been previously noted by Spezzano et al. (2025). In the present study, increasing the cosmic-ray ionisation rate is a way to achieve high abundances of NH2D and NHD2 in the inner parts of the core, but we are not sure how seriously the value ζH2 ~ 10−16 s−1 should be taken. It would agree with the estimate of Ivlev et al. (2019) based on the gas temperature measurements in L1544, and a model for the equilibrium gas temperature and size-dependent dust temperature. This model assumes the canonical grain size distribution of Mathis et al. (1977) with dn(a)/da ∝ a−3.5 in the range amin ≤ a ≤ amax, where amin = 5 nm and amax = 250 nm. As noted in Redaelli et al. (2021), this size distribution is unlikely to be valid in dense molecular clouds where the smallest grains are eliminated by coagulation (Silsbee et al. 2020), and taking this into account, the estimate from the theory of Ivlev et al. (2019) is lowered by a factor of 3. The revised value agrees with the results of Redaelli et al. (2021) giving an average ζH2 of 3 × 10−17 s−1 in L1544. The abundances of the molecular ions observed by Redaelli et al. (2021) are highly dependent on the cosmic-ray ionisation rate, and the estimate of these authors is probably more trustworthy than the high value obtained here. The lower value would be in line with recent estimates of the cosmic-ray ionisation rate in nearby diffuse clouds (<10−16 s−1; Obolentseva et al. 2024; Neufeld et al. 2024). These are based on updated gas density estimates from three-dimensional extinction maps and the excitation of diatomic carbon. On the other hand, one could expect a local enhancement of the cosmic-ray ionisation rate in the Taurus molecular cloud complex, including L1544, because of their apparent location in the intersection of large-scale structures known as the Local Bubble and the Per-Tau Shell (Zucker et al. 2022; Bialy et al. 2021). As discussed in Zucker et al. (2022), these structures are likely to have been generated by multiple supernova events in the past. Shock acceleration in the blast waves of supernova remnants is a major source of cosmic rays (e.g., Diesing & Caprioli 2019 and references therein).

In addition of the cosmic-ray ionisation rate, there are several other parameters in the chemistry model that influence the accretion and desorption of molecules and thereby their radial distribution in a dense core. However, the objective of the present study is to determine the o/p ratios of NH2D and NHD2. In what follows, we leave the difficulties in modelling the spatial distribution of these molecules aside and concentrate on their spin ratios.

Table 4

Fractional abundances, X, and abundance ratios in L1544 derived using a model where the abundances are constant in the density range nH2 = 5 × 104−5 × 105 cm−3 but decrease towards higher and lower densities.

7.2 Ortho/para ratios

Judging from the difficulty in fitting the oNH2D line at 86 GHz together with the other lines of NH2D, we think that observations of the ground-state lines of NH2D and NHD2 at 333 and 335 GHz provide the most reliable estimates for the o/p ratios of these molecules. In addition to the fact that the ground-state lines are observed with the same beam, their relative intensities are free from calibration issues, because they lie close in frequency and are observable in the same spectrometer band. As can be seen in the spectra shown in Figs. 7 and 9, simulations from the FS reaction scheme underpredict the intensity of the oNH2D(101 − 000) when the pNH2D(101 − 000) line agrees, and conversely overpredict the oNHD2(111 − 000) intensity when pNHD2(111 − 000) fits. Simulations from the PH reaction scheme give reasonably good agreement with both ortho- and para-lines at the same time. The differences between the simulated line intensities from the PH and FS reaction schemes are not large, but the FS model produces systematically discrepant line ratios of the ortho- and para-species. The line simulations described in Sect. 6, which used a fixed fractional abundance profile but let the absolute values float freely, support the PH scenario. The best-fit fractional abundances imply the spin ratios o/p-NH2D = 2.85 ± 0.05 and o/p-NHD2 = 2.10 ± 0.06 (taking into account 1σ errors). This means that it is reasonable to assume that the o/p ratios of NH2D and NHD2 in L1544 are simply determined by the nuclear spin statistical weights. We arrived previously at the same conclusion in the Ophiuchus cores H-MM1 and Oph D. This also agrees with the earlier results of Daniel et al. (2016) towards the dense cores B1b and I16293E.

Based on the arguments presented in Harju et al. (2024), we assume that the o/p ratios of NH2D and NHD2 observed in the gas are not inherited from molecules released from grains, but result from ion-molecule reactions in the gas, which thus must be regarded as PH reactions rather than reactions where the H and D nuclei are completely scrambled. The viability of mixing hydrogen or deuterium atoms in reactions leading to deuterated ammonia was briefly discussed in Rist et al. (2013) and Harju et al. (2017), with reference to theoretical calculations of Herbst et al. (1991) and Ischtwan et al. (1992) of the potential energy surface of the system [NH3..H2/HD/D2]+, which is the dominant pathway to the ammonium ion NH4+$\[\mathrm{NH}_{4}^{+}\]$. As discussed in Ischtwan et al. (1992), interchange of H or D atoms is possible through certain transition structures of the reaction complex, some of which have relatively low energies while others constitute substantial barriers. However, in some cases, a barrier can be avoided due to the nearly free internal rotation of the NH4+$\[\mathrm{NH}_{4}^{+}\]$ moiety in the exit channel complex [NH4..H]+. We are not aware of theoretical calculations that would address the reactions NH3 + H2D+/D2H+/D3+$\[\mathrm{D}_{3}^{+}\]$, which probably produce most of the deuterated ammonia in dense cores. We note that when only one hydrogen or deuterium atom is transferred from one molecule to another, the spin ratios maintain their statistical values even if nuclei are mixed within the two parts of the entrance and exit channel complexes.

As discussed in Sipilä et al. (2019a), the choice of the reaction scheme also affects the deuterium fractionation in addition to the spin ratios. Reactions where two deuterium atoms are transferred, such as NH3 + D2H+ → NHD2 + H2, are possible in the FS scheme but not allowed in the PH scheme (see Fig. 1 in Sipilä et al. 2019a). In the present simulations, the PH model predicts higher NH2D/NH3 ratios and lower NHD2/NH2D and ND3/NHD2 ratios than the FS model, especially at late times (t ⪆ 106 yr; see Fig. 6). Because the abundance of NHD2 peaks earlier than that of NH2D (see Fig. 5), the predicted NHD2/NH2D ratio generally does not correspond to the observed value at the time when the NH2D spectra agree with the observations. In contrast to previous simulations of Sipilä et al. (2019a), the ratioNH2D/NH3 remains low (~0.1) also in the PH model. The present chemistry models include chemical desorption, which enhances the replenishment of the gas with molecules formed in the icy mantles of grains, where deuterium fractionation is weaker. These simulations also show that the molecules released from the grains are quickly reprocessed in the gas, and the spin ratios are determined by the ion-molecule reactions no matter how efficient the desorption is (Harju et al. 2024). As evident in Fig. 6, the predicted deuterium fractionation ratios of ammonia vary significantly during the probable lifetime of a dense core (⪅106 yr). The figure suggests that fractionation ratios can be useful for chemical dating, while spin ratios do not serve this purpose.

8 Conclusions

We have determined the o/p ratios of the NH2D and NHD2 molecules in one of the most sheltered and coldest places in the local Universe, the inner parts of the pre-stellar core L1544 in the Taurus molecular cloud complex. In accordance with similar measurements towards two starless cores in Ophiuchus, the o/p ratios of both molecules are found to correspond to the nuclear spin statistical weights, o:p-NH2D = 3 : 1 and o:p-NH2D = 2 : 1. Because the ratios are likely to be determined by gas-phase ion-molecule reactions, these reactions are best described as single-proton or single-deuteron hops that maintain the statistical spin ratios. Ammonia formation on grain surfaces is predicted to proceed through H or D atom additions that also produce ortho- and para-species according to their statistical weights, so it seems reasonable to assume that the spin ratios of ammonia and its deuterated isotopologues are always statistical. Another consequence of adopting the PH reaction scheme is that singly deuterated ammonia is favoured at the cost of NHD2 and ND3. This scheme therefore leads to fractionation ratios NH2D/NH3, NHD2/NH2D, and ND3/NHD2 that are different from those predicted by the FS scheme.

Acknowledgements

The authors thank APEX staff for performing the observations, the anonymous reviewer for helpful comments on the manuscript, and the Max Planck Society for financial support.

Appendix A Observed and simulated APEX and IRAM spectra

The NH2D spectra observed towards positions corresponding to the seven pixels of the LAsMA array are shown in Fig. A.1. These positions were observed in the 101 − 000 lines with APEX. The 111 − 101 spectra are extracted from the maps observed with the IRAM 30 m telescope. The 111 − 000 spectra of NHD2 observed with APEX are shown in Fig. A.2. In both figures, the predictions using the abundance profile described in Sect. 6 are shown in red. In this profile, the maximum fractional abundance is assumed to be found at intermediate densities, with decreasing tendencies both towards the centre and towards the outer boundary of the core. The maximum values are given at the top of the panels. The offset from the core centre is indicated on the top right of each panel.

thumbnail Fig. A.1

NH2D spectra observed with the APEX and IRAM 30 m telescopes (blue lines). The positions correspond to the seven pixels of the LAsMA array. Upper panels: 111 − 101 spectra from IRAM 30 m. Lower panels: 101 − 000 from APEX. Simulated spectra from the model discussed in Sect. 6 are shown in red.
thumbnail Fig. A.2

111 − 000 spectra of oNHD2 and pNHD2 observed with the LAsMA array (blue lines). The simulations from the model discussed in Sect. 6 are shown in red.

Appendix B Observed and simulated ALMA spectra

The integrated pNH2D(111 − 101) line intensity map of L1544 and a selection of spectra observed with ALMA by Caselli et al. (2022) are shown in Fig. B.1. The two line segments where the spectra are extracted start from the column density peak (our (0,0)), and are aligned with the minor and major axis of the elliptical core as determined in Sect. 3. The spectra calculated from the ad hoc model discussed in Sect. 6 are shown in red. The position angle of the core as delineated by NH2D emission (~ − 25°) differs from the angle determined by 8 μ absorption (~ − 36°). Similar figures were shown in Caselli et al. (2022, their Figs. 3 and A.1), except that the strip along the main axis was directed south-east. There, the observed line emission weakens more quickly towards the core edge than in the north-west, where the spheroidal model agrees well with the observations.

thumbnail Fig. B.1

111 − 101 spectra of pNH2D observed with ALMA along two orthogonal strips across L1544. The data are from Caselli et al. (2022). (a): Integrated line intensity map. The positions where the spectra are extracted are indicated with numbered circles. The synthesised beam is shown in the bottom left. (b): Six spectra along the minor axis of the core (from the centre to the south-west). The number in the top left of each panel refers to the circled numbers on the map. The positions are separated by 2″. (c): Same as panel (b) but along the major axis of the core (from the centre to the north-west). The spectra predicted from the model described in Sect. 6 are shown in red. Unlike the integrated intensity map of panel (a), the observed spectra in panels (b) and (c) have been corrected for the primary beam response.

References

  1. Bialy, S., Zucker, C., Goodman, A., et al. 2021, ApJ, 919, L5 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bunker, P. R., & Jensen, P. 2005, Fundamentals of Molecular Symmetry (London: IOP Publishing Ltd) [Google Scholar]
  3. Caselli, P., Bizzocchi, L., Keto, E., et al. 2017, A&A, 603, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Caselli, P., Pineda, J. E., Zhao, B., et al. 2019, ApJ, 874, 89 [NASA ADS] [CrossRef] [Google Scholar]
  5. Caselli, P., Pineda, J. E., Sipilä, O., et al. 2022, ApJ, 929, 13 [NASA ADS] [CrossRef] [Google Scholar]
  6. Cheng, Y. C., Bockelée-Morvan, D., Roos-Serote, M., et al. 2022, A&A, 663, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Coudert, L. H., & Roueff, E. 2006, A&A, 449, 855 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Crapsi, A., Caselli, P., Walmsley, M. C., & Tafalla, M. 2007, A&A, 470, 221 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Daniel, F., Faure, A., Wiesenfeld, L., et al. 2014, MNRAS, 444, 2544 [NASA ADS] [CrossRef] [Google Scholar]
  10. Daniel, F., Rist, C., Faure, A., et al. 2016, MNRAS, 457, 1535 [NASA ADS] [CrossRef] [Google Scholar]
  11. Diesing, R., & Caprioli, D. 2019, Phys. Rev. Lett., 123, 071101 [Google Scholar]
  12. Faure, A., Hily-Blant, P., Le Gal, R., Rist, C., & Pineau des Forêts, G. 2013, ApJ, 770, L2 [NASA ADS] [CrossRef] [Google Scholar]
  13. Faure, A., Hily-Blant, P., Rist, C., et al. 2019, MNRAS, 487, 3392 [Google Scholar]
  14. Güsten, R., Nyman, L. Å., Schilke, P., et al. 2006, A&A, 454, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Güsten, R., Baryshev, A., Bell, A., et al. 2008, SPIE Conf. Ser., 7020, 702010 [Google Scholar]
  16. Hama, T., Kouchi, A., & Watanabe, N. 2018, ApJ, 857, L13 [NASA ADS] [CrossRef] [Google Scholar]
  17. Harju, J., Daniel, F., Sipilä, O., et al. 2017, A&A, 600, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Harju, J., Pineda, J. E., Sipilä, O., et al. 2024, A&A, 682, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Herbst, E., Defrees, D. J., Talbi, D., Pauzat, F., & Koch, W. 1991, J. Chem. Phys., 94, 7842 [Google Scholar]
  20. Hily-Blant, P., Faure, A., Rist, C., Pineau des Forêts, G., & Flower, D. R. 2018, MNRAS, 477, 4454 [NASA ADS] [CrossRef] [Google Scholar]
  21. Hugo, E., Asvany, O., & Schlemmer, S. 2009, J. Chem. Phys., 130, 164302 [Google Scholar]
  22. Ischtwan, J., Smith, B. J., Collins, M. A., & Radom, L. 1992, J. Chem. Phys., 97, 1191 [Google Scholar]
  23. Ivlev, A. V., Silsbee, K., Sipilä, O., & Caselli, P. 2019, ApJ, 884, 176 [Google Scholar]
  24. Jiménez-Redondo, M., Sipilä, O., Jusko, P., & Caselli, P. 2024, A&A, 692, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Juvela, M. 2020, A&A, 644, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Keto, E., & Caselli, P. 2010, MNRAS, 402, 1625 [Google Scholar]
  27. Ladjelate, B., André, P., Könyves, V., et al. 2020, A&A, 638, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Le Gal, R., Hily-Blant, P., Faure, A., et al. 2014, A&A, 562, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Mangum, J. G., Emerson, D. T., & Greisen, E. W. 2007, A&A, 474, 679 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
  31. Melosso, M., Bizzocchi, L., Dore, L., et al. 2021, J. Mol. Spectrosc., 377, 111431 [NASA ADS] [CrossRef] [Google Scholar]
  32. Miret-Roig, N., Galli, P. A. B., Olivares, J., et al. 2022, A&A, 667, A163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471 [NASA ADS] [CrossRef] [Google Scholar]
  34. Neufeld, D. A., Welty, D. E., Ivlev, A. V., et al. 2024, ApJ, 973, 143 [Google Scholar]
  35. Obolentseva, M., Ivlev, A. V., Silsbee, K., et al. 2024, ApJ, 973, 142 [NASA ADS] [CrossRef] [Google Scholar]
  36. Oka, T. 2004, J. Mol. Spectrosc., 228, 635 [NASA ADS] [CrossRef] [Google Scholar]
  37. Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Riedel, W., Sipilä, O., Redaelli, E., et al. 2023, A&A, 680, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Rist, C., Faure, A., Hily-Blant, P., & Le Gal, R. 2013, J. Phys. Chem. A, 117, 9800 [NASA ADS] [CrossRef] [Google Scholar]
  41. Sawada, T., Ikeda, N., Sunada, K., et al. 2008, PASJ, 60, 445 [NASA ADS] [Google Scholar]
  42. Schmiedt, H., Jensen, P., & Schlemmer, S. 2016, J. Chem. Phys., 145, 074301 [Google Scholar]
  43. Shinnaka, Y., Kawakita, H., Jehin, E., et al. 2016, MNRAS, 462, S124 [Google Scholar]
  44. Shirley, Y. L. 2015, PASP, 127, 299 [Google Scholar]
  45. Silsbee, K., Ivlev, A. V., Sipilä, O., Caselli, P., & Zhao, B. 2020, A&A, 641, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Sipilä, O., Harju, J., Caselli, P., & Schlemmer, S. 2015, A&A, 581, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Sipilä, O., Caselli, P., & Harju, J. 2019a, A&A, 631, A63 [Google Scholar]
  48. Sipilä, O., Caselli, P., Redaelli, E., Juvela, M., & Bizzocchi, L. 2019b, MNRAS, 487, 1269 [Google Scholar]
  49. Sipilä, O., Zhao, B., & Caselli, P. 2020, A&A, 640, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Spezzano, S., Redaelli, E., Caselli, P., et al. 2025, A&A, 694, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Woodward, C. E., Bockélee-Morvan, D., Harker, D. E., et al. 2025, psj, 6, 139 [Google Scholar]
  52. Zucker, C., Goodman, A. A., Alves, J., et al. 2022, Nature, 601, 334 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Observed transitions.

In the text
Table 2

Line parameters and column densities from Gaussian fits to the hyperfine structure.

In the text
Table 3

Fractional abundances, X, and abundance ratios in L1544 derived from the LAsMA spectra in pixels ① (0,0) and ④ (−23″, 35″) assuming that X is constant throughout the core.

In the text
Table 4

Fractional abundances, X, and abundance ratios in L1544 derived using a model where the abundances are constant in the density range nH2 = 5 × 104−5 × 105 cm−3 but decrease towards higher and lower densities.

In the text

All Figures

thumbnail Fig. 1

Footprint of the APEX/LAsMA array on the 8 μm surface brightness images of L1544 observed by the Spitzer satellite. The position of the central pixel is RA 05h04m17⋅s1$\[05^{\mathrm{h}} 04^{\mathrm{m}} 17^\text{s}_\cdot1\]$, Dec +25°10′43″ (J2000.0). The circle diameter corresponds to the HPBW of APEX at 334 GHz. The core appears as a shadow against the weak mid-infrared background. The bright spots are foreground stars.
In the text
thumbnail Fig. 2

NH2D spectra observed with IRAM and APEX towards the column density peak of L1544, which is the position of pixel ① in Fig. 1. The IRAM spectra are shifted up by 0.8 K for clarity.
In the text
thumbnail Fig. 3

NHD2 spectra observed with APEX towards the column density peak of L1544 (pixel ① in Fig. 1). The pNHD2 spectrum is shifted up by 0.3 K for display purposes.
In the text
thumbnail Fig. 4

Integrated JKa, Kb = 111 − 101 line intensity maps of oNH2D and pNH2D towards L1544 observed with the IRAM 30 m telescope. The intensity is on the TA∗$\[T_{\mathrm{A}}^{*}\]$ scale. The beam sizes at the line frequencies (85.9 GHz and 110.2 GHz) are indicated in the bottom right.
In the text
thumbnail Fig. 5

Evolution of the fractional abundances of oNH2D, pNH2D, oNHD2, and pNHD2 in the core model. The abundances correspond to the column density ratios N(mol)/N(H2) towards the centre of the core, smoothed to the resolution of the 18⋅′′7$\[18^{\prime \prime}_\cdot 7\]$ APEX beam at 333 GHz. Predictions from the PH scenario in the chemical network (a) and the FS scenario (b). In this simulation, a = 0.2 μm and ζH2 = 2.6 × 10−17 s−1. The horizontal dashed lines and thick semi-transparent lines represent the fractional abundances and their formal 1σ errors derived towards (0,0) in Sect. 3. These abundances are listed in Table 3.
In the text
thumbnail Fig. 6

Evolution of the o/p-NH3, o/p-NH2D, and o/p-NHD2 ratios and the fractionation ratios NH2D/NH3, NHD2/NH2D, and ND3/NHD2 according to the PH (a) and FS (b) models visualised in Fig. 5. The fractional abundances used for the ratios are calculated in the same way as those shown in Fig. 5. The horizontal straight dashed lines and semi-transparent thick lines represent the abundance ratios derived towards (0,0) from APEX observations (Sect. 3). Observational data for NH3 and ND3 are not available at the same resolution.
In the text
thumbnail Fig. 7

Observed (blue) and simulated (red) NH2D spectra towards (0,0). Left: simulations from the PH model with two different sets of model parameters. Right: results from the FS model with the same parameter values. The model parameters are given at the top of each panel. The upper plot in each of the four panels shows the 111 − 101 lines and the lower plot shows the 101 − 000 lines of oNH2D and pNH2D. The para-lines are shifted downwards for clarity.
In the text
thumbnail Fig. 8

Factional abundance profiles in the PH model core used to calculate the spectra shown in the left panels of Fig. 7.
In the text
thumbnail Fig. 9

Observed and simulated NHD2 spectra towards the column density peak of L1544. The red lines show simulations from the PH (top panel) and FS (bottom panel) reaction schemes. The same assumptions about the effective grain radius and the cosmic-ray ionisation rate have been used (a = 0.2 μm, ζH2 = 1.04 × 10−16 s−1) in the two models. The simulation time has been chosen to provide an approximate agreement with the pNHD2 spectrum.
In the text
thumbnail Fig. A.1

NH2D spectra observed with the APEX and IRAM 30 m telescopes (blue lines). The positions correspond to the seven pixels of the LAsMA array. Upper panels: 111 − 101 spectra from IRAM 30 m. Lower panels: 101 − 000 from APEX. Simulated spectra from the model discussed in Sect. 6 are shown in red.
In the text
thumbnail Fig. A.2

111 − 000 spectra of oNHD2 and pNHD2 observed with the LAsMA array (blue lines). The simulations from the model discussed in Sect. 6 are shown in red.
In the text
thumbnail Fig. B.1

111 − 101 spectra of pNH2D observed with ALMA along two orthogonal strips across L1544. The data are from Caselli et al. (2022). (a): Integrated line intensity map. The positions where the spectra are extracted are indicated with numbered circles. The synthesised beam is shown in the bottom left. (b): Six spectra along the minor axis of the core (from the centre to the south-west). The number in the top left of each panel refers to the circled numbers on the map. The positions are separated by 2″. (c): Same as panel (b) but along the major axis of the core (from the centre to the north-west). The spectra predicted from the model described in Sect. 6 are shown in red. Unlike the integrated intensity map of panel (a), the observed spectra in panels (b) and (c) have been corrected for the primary beam response.
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.