© P. Hily-Blant et al. 2020

1 Introduction

Protoplanetary disks provide the material out of which planets and primitive cosmic objects, such as meteorites and comets, form, and establishing their chemical composition has become a key objective (van Dishoeck et al. 2014; van Dishoeck 2018). More specifically, recent years have seen renewed interest in the isotopic ratio of interstellar nitrogen, with a view to establishing the degree to which planetary systems inherit their chemical composition from their parent interstellar clouds (Hily-Blant et al. 2013a). In this context, the unprecedented sensitivity of the ALMA interferometer plays an important role, in that it enables the nitrogen isotopic ratios to be measured in circumstellar disks. Similarly, the new generation of receivers at the IRAM 30 m telescope has allowed isotopic ratios to be measured in acceptable amounts of time. Isotopic ratios have long been used to trace the history of cosmic objects such as planets (e.g., Marty et al. 2017). In the case of nitrogen, however, isotopic ratios measured in a variety of solar-system objects present a rather confusing picture, with values ranging from 441 in Jupiter and the proto-Sun, down to approximately 50 in micron-sized inclusions in chondrites (Bonal et al. 2010). Interestingly, no ratio higher than the protosolar value of 441 has been reported so far in the Solar System. The current state of affairs is that the origin of nitrogen in the Solar System remains elusive (Hily-Blant et al. 2013a; Füri & Marty 2015). One of the current challenges is to determine the sources of isotopic ratio variations of nitrogen in star-forming regions and protoplanetary disks.

Progress in establishing the origin of nitrogen in the Solar System relies on a combination of theoretical and observational investigations of star-forming regions and protoplanetary disks, with the aim of obtaining an accurate picture of the processes that determine the isotopic ratios of nitrogen-bearing species. In so doing, care must be taken when comparing the 4.6 Gyr-old protosolar nebula to present-day star- and planet-forming regions since stellar nucleosynthesis, or the migration of the Solar System from a different galactocentric radius (Minchev et al. 2013; Martínez-Barbosa et al. 2015), could change the bulk ratio $R_{\text{tot}}$. Indeed, isotopic ratios well below 441, and as low as 270, have been measured in the diffuse (Lucas & Liszt 1998) and dense (Adande & Ziurys 2012; Hily-Blant et al. 2013a; Daniel et al. 2013) molecular local interstellar medium (ISM). When comparing the ratios in the present-day ISM to the value of 441 for the proto-Sun, it was realized that the above effects – stellar nucleosynthesis or the migration of the Solar System – must be included (Hily-Blant et al. 2017). Indeed, accurate measurements of nitrogen isotopic ratios in the local ISM performed in recent years (see Fig. 1 and Appendix A) have been found to be in very good agreement with predictions of the current bulk ratio from galactic chemical evolution models (Romano et al. 2017). The isotopic ratio in the local ISM is now proposed to be 330, thus significantly lower than the protosolar value of 441 (Hily-Blant et al. 2017, 2019). A low isotopic ratio in the solar neighborhood has received further observational support (Kahane et al. 2018; Colzi et al. 2018). Yet, the 330 value is still uncertain, and could be as low as ≈ 270 (Lucas & Liszt 1998; Adande & Ziurys 2012). This low value of the isotopic ratio is close to the terrestrial value (272), although this is likely a chance agreement rather than indicative of a true chemical link. In this paper, we will adopt $R_{\text{tot}}$ = 330 as the bulk isotopic ratio.

From a theoretical perspective, two processes are known that could lead to fractionation – that is, a deviation of the isotopic ratio in some species, such as HCN/HC¹⁵N, from the bulk value – in star-forming regions and protoplanetary disks. One is isotope-selective photodissociation of N₂ (Heays et al. 2014, 2017), which is expected to be effective in environments where UV radiation has a significant impact on chemistry, such as the upper layers of protoplanetary disks. A recent measurement of the ¹⁴N:¹⁵N ratio in the CN and HCN molecules in the circumstellar disk of the 8 Myr-old T Tauri object TW Hya demonstrated that the HCN:HC¹⁵N ratio increases with distance from the central star (Hily-Blant et al. 2019). This observational result lends support to the hypothesis that selective photodissociation of N₂ (Heays et al. 2014) could be the leading fractionation process in such disks. Disk models predict some enrichment of CN and HCN in ¹⁵N (Visser et al. 2018) but lead to lower fractionation levels than observed, which could result from several factors, including the UV flux or the dust properties (size distribution, mass fraction). In shielded gas, such as the collapsing clouds considered in the present study, isotope-selective photodissociation is relevant only through the differential dissociation of ¹⁴N₂ and ¹⁵N¹⁴N by the cosmic-ray induced ultraviolet radiation.

The second process is known as chemical mass fractionation (Watson 1976): heavier isotopic forms have lower zero-point energies, ℏω∕2, where ω ∝ μ^−1∕2 is the characteristic vibrational frequency, and μ is the reduced mass of the molecule. As a consequence of this energy difference, molecules can become enriched in a heavier isotope, through isotope-exchange reactions, when the kinetic temperature of the medium is comparable to, or lower than, the energy difference. Whilst this effect is most pronounced in the case of light molecules, such as H₂ and its deuterated forms, for which the differences in reduced mass, and hence in zero-point energies, are substantial, the rare isotopes of other elements, such as ¹³C and ¹⁵N, can also become enriched when the temperature is sufficiently low (Terzieva & Herbst 2000; Roueff et al. 2015). It has been proposed that the ¹⁵N-rich nitrogen reservoirs observed in the Solar System could be formed by chemical mass fractionation in the parent interstellar cloud (Charnley & Rodgers 2002; Hily-Blant et al. 2013a,b; Furuya & Aikawa 2018); but more sophisticated modeling of the observations has called this interpretation of nitrogen fractionation into question (Roueff et al. 2015). That chemical fractionation is not efficient, even in cold gas, is also supported by state-of-the-art determinations of the HCN:HC¹⁵N ratio in the L1498 pre-stellar core by Magalhães et al. (2018). These authors obtained HCN:HC¹⁵N = 338 ± 28, which is in harmony with the bulk nitrogen isotopic ratio $R_{\text{tot}}$ = 330 in the ISM (Hily-Blant et al. 2017) and indicates that chemical mass fractionation is inefficient, at least for HCN.

From the picture drawn above, one might conclude that models and observations agree qualitatively, in that ¹⁵N-rich reservoirs are formed in situ in protoplanetary disks by the action of isotope-selective photodissociation of N₂, inherited from a non-fractionated interstellar reservoir. However, measurements of the N₂ H⁺ isotopic ratio in dense clouds are a factor of ≈2−3 times higher than the bulk value $R_{\text{tot}}$ = 330, with no clear origin for the discrepancy; this is the most striking observational evidence that our models of nitrogen fractionation are still incomplete. Attempts to reproduce the observed N₂ H⁺:N¹⁵NH⁺ column density ratio with chemical models have failed so far (Roueff et al. 2015; Wirström & Charnley 2018; Loison et al. 2019). It may be noted that the fractionation reactions initially proposed by Terzieva & Herbst (2000) work no better in this respect. Accuracy in observational measurements of nitrogen isotopic ratios, in N₂ H⁺ and/or other species, has also improved significantly, suggesting that the discrepancy between observations and models could be due to an incomplete understanding of nitrogen fractionation, or even mistaken concepts of the nitrogen chemistry.

The main originality of the present study is to present time-dependent calculations of the isotopic ratio of nitrogen-bearing species during the early phases of the gravitational collapse of a pre-stellar core. Indeed, all previous studies (Terzieva & Herbst 2000; Charnley & Rodgers 2002; Rodgers & Charnley 2008; Hily-Blant et al. 2013b; Roueff et al. 2015; Wirström & Charnley 2018; Loison et al. 2019) have considered constant density models while focusing on the chemical aspects. Here, we consider a collapsing core starting from steady-state initial conditions, and explore various chemical processes with the hope of finding potential solutions to the problems outlined above. The following Sect. 2 summarizes the salient characteristics of the collapse model. In Sect. 3, we compare the computed degrees of fractionation of nitrogen with the results of previous calculations and with observations. In Sect. 4, the obtained results are discussed and compared to previous studies. Our conclusions are given in Sect. 5.

Fig. 1 Overview of measurements of the nitrogen isotopic ratio in low- to intermediate-mass cores and Class0. The horizontal line indicates the value of the bulk isotopic ratio of nitrogen adopted in this study ($R_{\text{tot}}$ = 330). Black dots indicate indirect measurements based on the double-isotope method. Colored points refer to direct measurements, with N2H+:N15NH+ outlined in blue. The N2H+:15NNH+ ratios are indicated by an asterisk. See Appendix A.1 and Table A.1 for details.

2 The model

2.1 Physical and chemical model

The model used in this study was first presented in Hily-Blant et al. (2018a; hereafter Paper I). Detailed information on the properties and validation of the model are given in Paper I and we here summarize its central assumptions and features. Our model follows the evolution of the chemical composition of a nonrotating, spherical cloud of gas and dust undergoing isothermal collapse. The chemistry is followed self-consistently with the collapse, whose treatment was inspired by the results of Larson (1969) and Penston (1969): The cloud has a free-falling core of radius R_c and uniform, but increasing, density ϱ_core, which is surrounded by an envelope with a density profile given by ${\varrho }_{\text{env}}(R)={\varrho }_{\text{core}}{(R/R_{\text{c}})}^{-2}$. As the collapse proceeds, the mass of the inner core decreases whilst that of the envelope increases by advection of gas and dust across the interface with the core, such that the total mass of the cloud $M=4\pi R_0^3{\varrho }_0/3$ remains constant. In this expression, the initial external radius is denoted by R₀ and the initial mass density by ϱ₀.

The simple model above captures quite well the Larson-Penston (L-P) solution for the time-dependent density during the early stages of collapse, before the formation of the first hydrostatic (Larson) core at a typical density n_H = n(H) + 2n(H₂) ~ 10¹⁰ cm⁻³ (ϱ_core ~ 10⁻¹⁴ g cm⁻³) (Larson 1969; Young et al. 2019) beyond which our model loses validity. In what follows, we study the initial phase of collapse and focus on the evolution of chemical abundances up to the stage corresponding to a central density no higher than 10⁸ cm ⁻³. The central density, which is a monotonically increasing function of time, parametrizes the evolution of abundances. We note also that the model allows us to compute separately the contributions of the envelope and of the inner core to the total column density of any species.

The initial chemical conditions are derived assuming that the cloud had time to reach a steady state prior to collapse. The subsequent chemical evolution of the 3-fluid medium, which comprises neutral, positively and negatively charged species, is followed by means of conservation equations of the form $$\frac{1}{n}\frac{\text{d}n}{\text{d}R}=\frac{N}{nv}-\frac{2}{R},$$(1)

where the number density of the species is denoted n(R) [ cm ⁻³], and N(R) [ cm ⁻³ s⁻¹] is its rate of formation per unit volume through chemical reactions; the flow speed is given by v(R) = dR∕dt. Grain coagulation is followed self-consistently during the collapse (Flower et al. 2005). The mass density of each of the three fluids is determined by the corresponding conservation equation $$\frac{1}{{\varrho }_{\text{env}}}\frac{\text{d}{\varrho }_{\text{env}}}{\text{d}R}=\frac{S}{{\varrho }_{\text{env}}v}-\frac{2}{R},$$(2)

where ϱ_env(R) [g cm ⁻³] denotes the mass density, and S(R) [g cm ⁻³ s⁻¹] is the rate of production of mass per unit volume. Throughout the paper, abundances are expressed relative to hydrogen nuclei, n(X)∕n_H, and the isotopic ratio of any species X, noted $R(X)$, may refer to abundance or column density.

2.2 Chemical network

For the purpose of the present study, the condensed¹ UGAN network (see Paper I) has been extended to include the ¹⁵N isotopologs (see Appendix B). A “sanity” check of the duplicated network was undertaken to verify that the ¹⁴N and ¹⁵N variants are independently conserved when no fractionation reactions are included.

We have considered two sets of fractionation reactions and rate coefficients, compiled by Terzieva & Herbst (2000) and Roueff et al. (2015) (hereafter TH00 and R15, respectively; see Table B.3), with a two-fold objective: first, to enable a direct comparison with our previous study (Paper I); and second, to compare the dependence of our results on the two sets available in the literature. Our fiducial network does not take into account the possibility that CN + N could lead to isotopic exchange (see also Wirström & Charnley 2018). The modifications to R15 made by Loison et al. (2019) have also been tested but will not be considered in detail since they were found to have a negligible effect on the isotopic ratios.

The bi-molecular reaction rate coefficients in our network are expressed as $$k(T)=\alpha {(T/300)}^{\beta }\exp (-\gamma /k_{\text{B}}T),$$(3)

and hence the form of the rate coefficient for the N⁺ + N₂ fractionation reaction, adopted by R15, could not be used. Nevertheless, because the kinetic temperature in our models is assumed fixed and equal to 10 K, our simplified expression for the reverse reaction, k(T) = 2.4 × 10⁻¹⁰ exp(−28.3∕T), deviates from theirs by only 3% at T = 10 K, and by less than 20% for T ≤ 30 K. The rate of the critical reaction N⁺ + H₂ → NH⁺ + H was approximated to the rate with o-H₂ (Dislaire et al. 2012), which is appropriate for the ortho-to-para ratio of H₂ (~ 10⁻³) and the low kinetic temperature considered here. The rate for the deuterated variant is taken from Marquette et al. (1988).

Since Hily-Blant et al. (2013b), new calculations of the rates of dissociation and ionization by cosmic-ray induced secondary photons have become available (Heays et al. 2017), and the results of these calculations have been incorporated in the present model (see also Paper I). Furthermore, compared to Paper I, cosmic-ray induced desorption of water and ammonia has been revised by Faure et al. (2019). As a result, desorption of ammonia (and water) via direct impact of cosmic-rays is found to be non-negligible even at 10 K owing to the prompt desorption mechanism (Bringa & Johnson 2004). In this context, we emphasize that grain-surface chemistry in the UGAN network is handled in a simplified way, where sticking atoms are assumed to react with hydrogen on the surface, forming saturated molecules such as CH₄, H₂O, and NH₃ (see Table B.1). The binding energies of the ¹⁴N and ¹⁵N variants of adsorbed species have equal values.

The initial (deuterated) UGAN network consists of 1196 reactions linking 156 species². After including the ¹⁵N variants of nitrogen-bearing species and the fractionation reactions, the ¹⁵N version of UGAN contains ≈1720 reactions and 205 species.

2.3 Model parameters and initial composition

As in Paper I, our model calculations of the evolution of atomic and molecular abundances proceed in two steps: first, a calculation at constant density (hereafter step 1) is followed until a steady state is attained. In this first step, gas-grain processes (adsorption and desorption) and grain-surface chemistry are ignored. The radius of the grains includes the contributions of a refractory core and of an initial ice mantle; it remains constant and intervenes only in establishing the free-electron concentration, through electron attachment. Second, the abundances are computed following our L-P prescription of the collapse (hereafter step 2). The initial abundances are the steady-state values, and adsorption and desorption processes and grain-surface chemsitry are included.

In the present calculations, the rate of cosmic-ray ionization of H₂ was ${\zeta }_{{\text{H}}_2}=3\times {10}^{-17}$ s⁻¹. The refractory core has a radius a_c = 0.1 μm and a mass density ϱ_c = 2 g cm⁻³, while that of the mantles is ϱ_m = 1 g cm⁻³. The total grain radius, including the mantles, is thus initially a_gr = 0.128 μm (see Eq. (C.2). The constant density models have n_H = 10⁴ cm⁻³, and the kinetic temperature is 10 K in both steps. As in Paper I, the total mass M is equal to the Jeans mass of a cloud with density n_H = 10⁴ cm⁻³ and temperatureT_kin = 10 K, viz. M = M_J ≈ 7M_⊙, which is also close to the critical mass of a pressure-truncated, isothermal, Bonnor-Ebert sphere. The free-fall timescale of such a sphere of gas is ${\tau }_{\text{ff}}={(3\pi /32G{\varrho }_0)}^{1/2}=4.4\times {10}^5$ yr.

3 Results

3.1 CO depletion and the timescale of collapse

The time required for nitrogen-bearing species to reach steady state (often called chemical equilibrium), through gas-phase reactions can exceed 1 Myr for a cloud of density n_H = 10⁴ cm⁻³. Since the quasi-static evolution of cores may last 0.5–1 Myr (Brünken et al. 2014), it is not guaranteed that the initial abundances of these species in the gas phase, prior to collapse, are equal to their values in steady state. Nevertheless, with the possible exception of the role played by the N₂:N abundance ratio in the nitrogen chemistry, the abundances when collapse sets in are not crucial to the calculation of the composition of the gas during its subsequent evolution providing that its composition adapts rapidly to the prevailing physical conditions. As the cloud collapses, and its density increases, this condition will tend to be satisfied. However, freeze-out on to the grains occurs simultaneously, reducing the abundances of molecules in the gas phase (and hence the likelihood of their being observed). It follows that the assumption of an initial chemical equilibrium may influence the predicted composition of the gas during the early stages of its gravitational collapse; this caveat should be borne in mind when considering the results in Sect. 3.3.

A key question in star-formation studies is the dynamical timescale, which may be constrained by time-dependent models of collapse that include the gas-grain chemistry, since both gas-phase and adsorption processes depend on the density. In our nonrotating, magnetic-pressure-free model, the dynamical timescale is that of the gravitational collapse, τ_ff. Nevertheless, the collapse may be slowed by magnetic fields and proceed on the ambipolar diffusion timescale, which can be up to 10 times larger than τ_ff. Such a delay can be introduced in our model calculations by increasing artificially the dynamical timescale (Flower et al. 2005). However, delaying the collapse enhances freeze-out which, if not counterbalanced by desorption, will lead to lower gas-phase abundances.

In Fig. 2, we compare the depletion of CO with increasing density³, as predicted by our L-P simulation, with the observations of L183 by Pagani et al. (2007, 2012). Results are shown for collapse timescales corresponding to 1, 2, and 3 times τ_ff and an initial grain-core radius a_c = 0.1 μm. It may be seen that, when the dynamical timescale is equal to the free-fall time, there is striking agreement with the observational results of Pagani et al. (2007, 2012). On the other hand, when the timescale for collapse is increased, the agreement with the observations deteriorates significantly. Short collapse timescales, ≤ 0.7 Myr, are also supported by an independent observational constraint, based on the D/H ratio in N₂ H⁺ (Pagani et al. 2013). Also displayed in this figure is the result of calculations with a_c = 0.05 μm for a dynamical timescale equal to τ_ff, showing that the agreement with the observed abundances is preserved.

In what follows. we thus consider cores collapsing on a free-fall timescale.

Fig. 2 Depletion of CO with increasing density, as predicted by our L-P collapse models for three values of the dynamical timescale, tdyn, and ac = 0.10 μm, compared with the observational results of Pagani et al. (2007, 2012) for L183. The L-P calculations corresponding to tdyn = tff, where tff is the free-fall time, are shown also for ac = 0.05 μm.

Fig. 3 Results from chemical models. In each row, left panel: temporal evolution during the quasi-static evolution of the core (step 1), and right panel: evolution of the inner core as a function of the central density during the collapse (step 2). The bulk 14N:15N isotopic ratio is $R_{\text{tot}}$ = 330. (a) Gas-phase abundances, relative to nH, of the main isotopologs of relevant species. (b) 14N/15N abundance ratios when adopting the TH00 set of fractionation reactions. (c) Same as (b) but with the fractionation reactions from R15.

3.2 Steady-state abundances and isotopic ratios

In the left panels of Fig. 3 are plotted the evolution toward steady state of the gas-phase fractional abundances of nitrogen-containing species during step 1 and the predicted isotopic ratios using the reactions and rate coefficients of either TH00 or R15. The initial gas- and ice-phase abundances of the elements are listed in Table 1, with hydrogen in molecular form and metals in atomic form, either neutral (¹⁴N, ¹⁵N, O) or singly ionized (C⁺, S⁺, Fe⁺). The results are summarized in Table 2.

Our calculations agree with previous studies (Wirström & Charnley 2018) in that the degrees of fractionation are globally low; this is in agreement with the directly measured ratios in starless cores and star-forming regions summarized in Table A.1. The degrees of fractionation are lower when the reaction rates of R15 are adopted, rather than those of TH00. We note significant differences between the present calculations and those of our previous study (Hily-Blant et al. 2013b), which are due primarily to the update of the rates of dissociation and ionization by the cosmic-ray induced ultraviolet radiation field, following the results of Heays et al. (2017). For example, the rate of dissociation of N₂ is smaller by afactor of approximately 25 than the rate used previously, which explains the higher ratio of the steady-state abundances of N₂ and N, n(N₂)∕n(N), in the present case.

We have investigated the effects of varying ill-defined parameters such as the grain-core radius, a_c, ${\zeta }_{{\text{H}}_2}$, the elemental S abundance, and the C:O elemental abundance ratio (in practice, the carbon abundance). The effects on the degrees of fractionation of ¹⁵N were generally modest, being somewhat more pronounced when the rate coefficients of TH00 were adopted, rather than those of R15.

3.3 Fractionation during gravitational collapse

The principal question that we wish to answer is whether the degrees of fractionation of the nitrogen-bearing species, computed in steady state, vary significantly in a subsequent gravitational collapse. Our setup is described in Sect. 2.3.

As the density of the medium increases, molecules ultimately freeze on to the grains, building up layers of ice that contribute, in addition to grain coagulation, to the grains’ overall radius. Cosmic-ray induced processes restore molecules to the gas phase from the grain mantles, but the rates of evaporation increase less rapidly than the rates of freeze-out, and hence molecular species become depleted from the gas phase.

The evolution of the abundances and ¹⁴N:¹⁵N abundance ratios with the central density (or, equivalently, time) during step 2, for selected species, is shown in the right panels of Fig. 3. As can be seen, the fractional abundances decrease with increasing density, for all species. However, not all species behave identically as the density increases, and the differences among species are more pronounced with the TH00 reaction rates. In this case, the isotopic ratio in atomic nitrogen decreases from its steady-state value of 440 to 246 when the density has reached n_H = 10⁷ cm⁻³. Other species, directly related to atomic nitrogen, namely CN and NO, also see their isotopic ratios decrease by significant factors. To a lesser extent, this is also true for species, such as NH₃ and HCN, that are more indirectly related to N. In contrast, the isotopic ratios of N₂ and N₂ H⁺ are approximately constant during the collapse. When adopting the R15 rates, the variations are obviously smaller because the steady-state ratios are closer to the isotopic ratio of bulk nitrogen. Nevertheless, there are variations among species.

We find that, with both the TH00 and R15 fractionation reactions, the gas becomes gradually enriched in the heavier isotope of ammonia in the early stages of the collapse, but this tendency ultimately reverses as desorption of the initial reservoirs⁴ of the ¹⁴NH₃ and ¹⁵NH₃ ices begins to be significant (see Fig. 3). We recall that ¹⁴NH₃ and ¹⁵NH₃ ices are assumed to be initially present in the grain mantles in an isotopic ratio of 330. At high density, ammonia experiences an increase in its isotopic ratio, which is also observed for both isotopologs of N₂ H⁺.

It may be noted that, in the course of the collapse, the gas-phase fractionation reactions have little influence, relative to adsorption to and desorption from the grains. Consequently, the differences in the degrees of fractionation, as calculated with the fractionation reactions of TH00 and R15, are attributable essentially to the differences in their initial, steady-state values.

4 Discussion

4.1 Depletion-driven fractionation

The decrease in the isotopic ratios with increasing density during collapse seen in Fig. 3 (step 2, right-hand panels) is due to adsorption: owing to their higher mass, ¹⁵N-containing species have slightly slower thermal speeds than their ¹⁴N variants, and hence their collision frequencies with the grains are lower (see also Loison et al. 2019). The later increase in the ratios at densities above ~ 5× 10⁶ cm⁻³ is instead a result of desorption by cosmic-rays and secondary photons.

One key parameter in the adsorption-driven fractionation is the freeze-out timescale, which can be written ${\tau }_{\text{fo}}=x_0{\tilde{m}}^{1/2}$ (see Eq. (C.7) with $\tilde{m}$ the atomic mass number and x₀ a factor which depends on the density, temperature and on the grain properties (size, mass density). For the grain parameters adopted in this study (see Sect. 2.3), one obtains x₀ = 7.9 × 10⁴ yr. In a constant density, isothermal, model, and ignoring adsorption/desorption processes, x₀ is constant: The difference in adsorption rates for the ¹⁴N and ¹⁵N variants is thus $\delta {\tau }_{\text{fo}}/{\tau }_{\text{fo}}=1/2(\delta \tilde{m}/\tilde{m})$ or 1.8% for HCN, and twice larger for N atoms. It follows that the isotopic ratio $R$ would decrease with time as $$R=R_0{e}^{-t/{\tau }_{\text{fr}}},$$(4)

where $R_0$ is the initial isotopic ratio and the e-folding timescale is given by (see Eq. (C.14) $${\tau }_{\text{fr}}\approx 2x_0{\tilde{m}}^{3/2}\approx 2\tilde{m}{\tau }_{\text{fo}}.$$(5)

The depletion-driven fractionation timescale, τ_fr, is thus significantly longer than the freeze-out one.

Our simple analytical estimate agrees well with the result of the constant-density study of nitrogen fractionation by Loison et al. (2019) for HCN, but not for HNC. Furthermore, these authors report variations by several tens of percent for the isotopic ratio of atomic nitrogen after typically 0.4 Myr, while the above derivation predicts a decrease of 13% only (see Appendix C). The above description of depletion-driven fractionation is also in good agreement with our constant-density calculations, shown in Fig. C.1. For our model parameters, and with a_gr = 0.128 μm, the timescale for depletion-driven fractionation is τ_fr = 8.3 Myr for N atoms, and 22 Myr for HCN molecules.

From Eq. (5) one would expect to see a dependence of τ_fr upon the mass of the depleting species, which is not observed in the models. Instead, all species follow a single exponential decay which coincides with that of atomic nitrogen. This result indicates that the temporal variation of $R$(t) is determined primarily by the differential depletion of ¹⁴N and ¹⁵N, which propagates to other species through gas-phase chemistry. It may be verified that the timescale for the formation of N-bearing gas-phase species (~ 4 Myr at a time t = 0.1 Myr and 1 Myr at t = 0.3 Myr) is significantly longer than the freeze-out timescale of N atoms but shorter than τ_fr for the N atoms, and much shorter than τ_fr for the heavy species (m ≈ 28 amu). Therefore, the freeze-out of N-bearing molecules will reflect $R$(t) of atomic N. Eventually, desorption brings the gas-phase isotopic ratios back toward the value of 330 (see Fig. C.1) that is adopted for the initial composition of the ices.

When collapse is taken into account, the density increases and hence the freeze-out timescale τ_fo ~ 1∕n_H decreases faster than the free-fall one, $t_{\text{ff}}~n_{\text{H}}{}^{-0.5}$. The effect of differential depletion is thus more pronounced than in a constant density model. Indeed, as Fig. 3 shows, $R$ decreases by more than 10%, when using the R15 rates, for most N-bearing molecules as the density increases from 10⁴ to 10⁶ cm⁻³. As the density increases above ~5 × 10⁶ cm ⁻³, desorption by direct cosmic ray impact liberates ¹⁵N-poor ammonia ices in the gas-phase thus raising the value of $R$(NH₃). The same mechanism also operates for other species but at higher densities (not shown in Fig. 3) The decrease in $R$ due to differential depletion is significantly larger – reaching ~30% for atomic nitrogen, NO, and CN – when the TH00 rates are used. Calculations at higher temperatures show that, in a warmer gas, desorption would occur earlier during collapse, at densities typically ~ 10 times smaller at a temperature T_kin = 30 K.

4.2 Comparison with observations

In Fig. 4, we compare the ¹⁴N:¹⁵N column density ratios of several nitrogen-containing species to the values obtained through direct measurements in pre-stellar cores (see Appendix A). The case of N₂ H⁺ is shown in Fig. 5 and is considered further in Sect. 4.2.2.

Fig. 4 Column density ratios of 14N- and 15N-containing species in the course of L-P collapse (step 2), compared with the values observed in pre-stellar cores through direct (filled symbols) and indirect (open symbols) measurements (see Fig. 4 and Appendix A.1).

4.2.1 Nitriles

Because direct measurements in CN and HCN yield essentially non-fractionated values of $R$ (see Table A.1), it is not surprising that our model predictions are in relatively good agreement with the observed ratios (see Wirström & Charnley 2018). However, the new result here is that this remains true during the collapse, which brings the gas-phase isotopic ratios back toward the value of 330.

Figure 4 appears to indicate that the direct measurement of $R$(CN) toward L1498 – which has been re-analyzed using a non-LTE model of the radiative transfer in spherical geometry (see Appendix A.2) – favors the reaction set of R15 to that of TH00. The ratio toward Barnard 1b is compatible with both reaction sets, although it suggests that $R$(CN) could be lower than 330. Regarding HCN, the only direct measurement was obtained toward L1498, and our model predictions, using either TH00 or R15, are in agreement with the observed value. The indirectly determined ratio toward L183 is significantly lower than the model predictions and the isotopic ratio for the bulk nitrogen. We note that it is impossible to anticipate the sense in which this ratio would change if a detailed radiative-transfer model were to be used in the analysis of the observations; the same holds true for the indirect measurement of $R$(CN) toward L1544.

4.2.2 N₂H⁺

The resultsof our fiducial model for N₂ H⁺ are shown as black lines in Fig. 5. In practice, our results fit the observations no better than previous attempts, by other groups. Thus, this species provides evidence of an unresolved issue with the nitrogen chemistry.

¹⁵NNH⁺ and N¹⁵NH⁺ form mainly in the proton-transfer reaction of ¹⁵NN with H${}_3^{+}$, for which the rate coefficient has been taken identical to that for the reaction N₂ (H${}_3^{+}$, N₂ H⁺)H₂. Were the rate coefficients of the reactions ¹⁵NN(H${}_3^{+}$, ¹⁵NNH⁺)H₂ and ¹⁵NN(H${}_3^{+}$, N¹⁵NH⁺)H₂ to be significantly smaller, the depletion of the ¹⁵N-bearing isotopic forms could be understood. However, the rate coefficients for these barrierless exothermic proton-transfer reactions are presumably independent of the isotopic form of nitrogen.

Alternatively, enhanced rates of dissociative recombination (DR) of ¹⁵NNH⁺ and N¹⁵NH⁺ with electrons at T = 10 K – the main destruction mechanism – might account for their depletion, as proposed also by Loison et al. (2019). While the electronic potential curves are identical for all isotopologs, the vibrational and rotational levels are shifted due to the differing masses, changing their position relative to the dissociation curve. The measurements of Lawson et al. (2011) have shown that, at 300 K, the DR of N₂ H⁺ proceeds faster than that of ¹⁵N₂H⁺, by about 20%. A larger effect, in the opposite sense, is possible at very low temperature, where an indirect mechanism dominates in the recombination of N₂ H⁺ (dos Santos et al. 2016). This tentative explanation requires experimental or theoretical confirmation, and we note in this context that calculations are in progress, using Multichannel Quantum Defect Theory (D. Talbi, priv. comm.).

Our model confirms the increase in the isotopic ratio of N₂ H⁺ when the DR rate of ¹⁵NNH⁺ and N¹⁵NH⁺ is multiplied by a factor of ɛ compared to the DR rate of N₂H⁺. Figure 5 shows that the isotopic ratio scales approximately linearly with ɛ; this is particularly true at high densities: the ratio rises from 330 to 560 for ɛ = 2, and reaches ≈825 when ɛ = 3. Thus, it appears that a factor of two to three enhancement in the DR rate coefficients yields values of $R$ for both N¹⁵NH⁺ and ¹⁵NNH⁺ which are compatible with those observed by Redaelli et al. (2018). We also note that the enhancement preserves the slightly lower value of the ratio for N₂ H⁺/N¹⁵NH⁺ compared to N₂ H⁺/¹⁵NNH⁺ (Kahane et al. 2018; Redaelli et al. 2018).

In Fig. 1, the isotopic ratios in N₂ H⁺ are shown with the sources sorted by increasing temperature from left to right, with L1544 and L183 being the coldest and OMC2-FIR4 the warmest. Although there are uncertainties in the temperature determination, the I16293E pre-stellar core is known to beinfluenced by the nearby I16293-2422 protostar and the kinetic temperature is significantly higher – increasing from 10 to 15 K from center to edge (Crapsi et al. 2007; Pagani et al. 2013; Bacmann et al. 2016) – than in the cold cores sampled by Redaelli et al. (2018). Similarly, the temperature in Barnard 1b and OMC-FIR4 is higher than in the cold cores of their sample (Daniel et al. 2013; Crimier et al. 2009). A trend is thus apparent whereby the isotopic ratio in N₂ H⁺ decreases as the temperature increases, reaching the bulk ratio of 330. We note that the value of $R$(N₂H⁺) measured toward the three warm sources is compatible with a value of ɛ close to 1, which might indicate a temperature dependence of ɛ, from ≈2 below ~ 10 K toward 1 above ~15 K. Alternatively, it could be that, in these sources, the destruction of N₂ H⁺ is not dominated by DR with electrons but by other reactions, such as with CO.

Finally,from Figs. 4 and 5, we see that, when ɛ > 1, both column density ratios of $R$(N₂H⁺) increase with density, whereas for other species, the ratios decrease from their steady-state value. This could produce some source-to-source variation of the isotopic ratios of different species.

Fig. 5 Isotopic ratio (column density) of N2H+ as a functionof the core density (nH) in the course of L-P collapse (step 2), compared with the values observed in pre-stellar cores (see also Table A.1). The abundance ratios N2 H+:N15NH+ and N2 H+:15NNH+ are shown for three values of the enhancement factor, ɛ, of the rate coefficient for dissociative recombination of N15NH+ (see Sect. 4.2.2), viz. ɛ = 1, 2, and 3. Results with rates from both R15 (full curves) and TH00 (dashed) are shown.

4.2.3 Ammonia

Ammonia is,after N and N₂, an important reservoir of nitrogen in pre-stellar cores and also in cometary ices. Establishing the value of $R$ in ammonia at early stages of star formation can therefore bring valuable clues as to the origin of the isotopic ratios measured in comets. The fundamental rotational transition of ammonia cannot be observed in cold gas from the ground. However, direct measurements of $R$ in the closely related deuterated variant NH₂D have been made by Gerin et al. (2009), using the 1_1,1−1_0,1 rotational transition of ortho-NH₂D at 86 GHz. Despite its very weak intensity, the main hyperfine component of the manifold of ¹⁵NH₂D was detected at the 3σ level on the integrated intensity toward four sources. The values of $R$ extend from 360 up to 810. Despite the large error bars, the depletion in ¹⁵N is manifest, as can be seen in Fig. 4 (see also Table A.1). We note that the excitation temperature of ¹⁵NH2D from the LTE analysis of Gerin et al. (2009) is likely to be overestimated (Daniel et al. 2013), and so would be the isotopic ratios, by up to a factor of two, bringing down their values toward the bulk isotopic ratio $R_{\text{tot}}$ = 330. Assuming a 15% decrease in the excitation temperature of the ¹⁵N variant, relative to the main isotopolog, as suggested by Daniel et al. (2013), we find $R$ = 193, 260, and 590, and >350, for NGC1333-DCO⁺, L134N(S), L1689N, and L1544, respectively. Given the large error bars (typically 30% downward and 30–100% upward; see Gerin et al. 2009), the ratios in the first two sources would then be marginally compatible with our prediction of no fractionation, as shown in Fig. 4. However, the ratio measured toward L1689N and the lower limit in L1544 indicate depletion of ¹⁵N which is not compatible with our results.

In the above context, we note that the charge transfer reaction between He⁺ and N₂ is another process with a known isotopic effect. This reaction has two output channels, $$\begin{array}{l}{\text{N}}_2+{\text{He}}^{+}\to {\text{N}}_2^{+}+\text{He}\\ {\text{N}}_2+{\text{He}}^{+}\to \text{N}+{\text{N}}^{+}+\text{He},\end{array}$$

and has attracted much interest because the energy needed to effect the transition ${\text{N}}_2({\text{X}}^1{\Sigma }_{\text{g}}^{+})\to N_2^{+}(C^2{\Sigma }_u^{+})$ nearly equals the helium ionization potential (see Hrodmarsson et al. 2019, and references therein). In the UGAN network, the total rate constant is 1.2 × 10⁻⁹ cm³ s⁻¹ and the branching ratio between the dissociative and non-dissociative channels is 1.94, for both N₂ isotopologs. However, this branching ratio is known to be sensitive to which isotopolog of N₂ is involved (Govers et al. 1975). Room temperature measurements suggest that the dissociative:non-dissociative branching ratios could be 1.5 for ¹⁴N¹⁴N and 0.69 for ¹⁴N¹⁵N, with a total rate coefficient of 1.3 × 10⁻⁹  cm ³ s⁻¹ (Govers, priv. comm.). Once again, theoretical or experimental investigations are required in order to establish the values of these branching ratios at low temperatures.

We have tested the impact of these branching ratios in our fiducial model (see Sect. 2.3 for its parameters). Our results (see Appendix D) indicate that the largest effect is on ammonia and NH₂D, which see their steady-state isotopic ratios increase to ≈450 and then decrease during the collapse. Thus, it appears that an explanation for a 30% depletion of ammonia in ¹⁵N can be found in the difference in the branching ratios for the reactions of ¹⁴N¹⁴N and ¹⁴N¹⁵N with He⁺. We note that, with the new branching ratio, $R$ increases also for HCN and CN, while that of N₂ H⁺ decreases, which conflicts with the observations. As to what could lead to the branching ratios varying among sources: temperature differences might have been considered, as for the factor ɛ in the case of N₂ H⁺; but, in the present context, L1544 and L1689N are respectively cold (6 K in the innermost regions) and relatively warm (10–15 K) cores.

5 Concluding remarks

This study of the fractionation of nitrogen-containing species is the first to consider the effects of the collapse on the $R$ abundance ratios. It was found that observations of CO tightly constrain the characteristic timescale of the collapse to be on the order of the free-fall timescale. Delaying further the collapse would imply more efficient desorption processes than those considered in our study (Vasyunin & Herbst 2013), as our model predicts that CO would be one order of magnitude smaller than is derived from the observations if the timescale was only a factor of two larger than in pure free-fall.

During the quasi-static evolution of the pre-stellar core, the degrees of fractionation of nitrogen-bearing species are modest, as found in earlier studies, and more so when the fractionation reactions and rate coefficients of R15, rather than those of TH00, are adopted. We have shown that the differential rates of photodissociation of N₂ and ¹⁵NN by the cosmic-ray induced ultraviolet radiation field are crucial to the degree of fractionation of molecular nitrogen, particularly under conditions of chemical equilibrium.

Regarding nitrogen fractionation, the main new result of this study is that, during the gravitational collapse of a pre-stellar core, differential enrichment of the gas in ¹⁵N-containing species occurs owing to their higher mass and hence lower rate of thermal collisions with and adsorption to grains. The increase in the density enhances the effect of depletion-driven fractionation, which can be up to 30% for species suchas atomic nitrogen, CN, and NO. Gas-phase fractionation reactions are relatively insignificant during gravitational collapse, as their associated timescales exceed the free-fall time. For CN, the predicted ratio is in good agreement with the revised observational value toward L1498. We also note that, similar to HCN, the CN/¹³CN ratio is ≈ 40, which is substantially lower than the value of ¹²C/¹³C = 70 for carbon in the local ISM; this underlines the danger of using the double-isotope-ratio method, discussed in Appendix A.1.

We have explored the influence of different rates for the dissociative recombination of N₂ H⁺ and its ¹⁵N-variants, showing that the observed depletion of N₂ H⁺ in ¹⁵N could be explained if N¹⁵NH⁺ and ¹⁵NNH⁺ dissociatively recombined more rapidly with electrons than ¹⁴N₂H⁺. Furthermore, our results suggest that the enhancement of the dissociative recombination rate of the ¹⁵N-variants increases as the temperature decreases. However, high-precision and -accuracy theoretical calculations are required to put this possibility on a firmer footing.

It is probably naive to expect to be able to reproduce the nitrogen isotopic ratios observed in different sources by means of a single model. In fact, there is evidence suggesting that variations of environmental parameters need to be taken into account. Figure 1 shows that, for N₂ H⁺, most measurement in pre-stellar cores yield ratios that are ~2–3 times larger than the bulk value $R_{\text{tot}}$ = 330 for nitrogen; only one measurement is consistent with 330. Further evidence is provided by the interferometric study of Colzi et al. (2019), showing spatial variation of the abundance ratio in a single protostar. In the context of low-mass pre-stellar cores, source-to-source variations – even when located in the same large-scale environment – were also proposed by Magalhães et al. (2018) to explain the HCN/HC¹⁵N ratio. Hily-Blant et al. (2017) have noted the large scatter of the directly measured CN/C¹⁵N ratio along various lines of sight through the diffuse ISM (Ritchey et al. 2015), with the values being mutually inconsistent by more than 2σ. Thus, data show that several species – CN, HCN, N₂ H⁺ – exhibit variations that suggest some sensitivity of the isotopic ratio to local physical conditions. In fact, that the environment could induce variations of the $R$ abundance ratios is not surprising, since we know that mass fractionation is a temperature-dependent process, and differential photodissociation depends on the cosmic-ray ionization rate and dust size distribution. Therefore, one might contemplate using nitrogen isotopic ratios to constrain the physical parameters, once our understanding of the interstellar nitrogen chemistry and fractionation processes is adequate for this purpose.

Acknowledgments

We thank the anonymous referee for useful comments and suggestions which helped improving the clarity of this article. We also thank Laurent Pagani for providing us with the density profile in L183, and Tomas Govers for recommending us the branching ratio between the dissociative and non-dissociative channels in ¹⁵N¹⁴N. This work was supported by the LABEX OSUG@2020 and by the Programme National Chimie Interstellaire (PCMI), and made use of the GRICAD infrastructure, which is supported by Grenoble research communities. D.R.F. acknowledges support from STFC (ST/L00075X/1), including provision of local computing resources.

Appendix A Observations

A.1 Review of existing measurements

We present, in Table A.1, an overview of measurements of the nitrogen isotopic ratio in pre-stellar cores, derived from observations of different molecular species. Direct and indirect measurements are shown separately. We recall that indirect measurements are based on the so-called double-isotope method, used to circumvent the problem of allowing for the optical depth in the line of the main isotopolog. For instance, one could use a ¹³C substitute to the ¹²C species $$\text{H{^}12\text{C}N/H{^}12\text{C}{^}15\text{N}}=(\text{H{^}12{\text{C}N/H}}^{13}\text{CN})\times {\text{H}}^{13}\text{CN}/\text{H{^}12{\text{C}}}^{15}\text{N}$$

or a deuterated variant $$\text{H{^}12\text{C}N/H{^}12{\text{C}}}^{15}\text{N}=({\text{H}}^{12}{\text{CN/D}}^{12}\text{CN})\times {\text{D}}^{12}{\text{CN/H}}^{12}{\text{C}}^{15}\text{N}$$

but one has to assume a value for the term in brackets, which is not measured. The usual assumption consists of adopting the isotopic ratio of the element, unless there is clear evidence for systematic deviation. Thus, the (H¹²C N/H¹³C N) term is usually taken equal to ¹²C/¹³C = 70 in the local ISM. However, such an approximation is not always justified and may lead to error.

Direct measurements can be made by several means. Gerin et al. (2009) proposed using the deuterated form of ¹⁵N-ammonia to measure the nitrogen isotopic ratio directly through the ratio NH₂D/¹⁵NH₂D. A more general way of obtaining direct ratios was put forward by Daniel et al. (2013), based on radiative transfer calculations with the ALICO code, which takes into account the contribution of overlap of the hyperfine transitions in the excitation of the hyperfine levels. This approach was further developed by Magalhães et al. (2018), who combined ALICO with a Markov-chain Monte Carlo (MCMC) exploration of the parameter space. This method enabled the first direct measurement of $R$ for HCN in a pre-stellar core; the result was found to differ from estimates based on the double-isotopic-ratio method. Direct measurements in pre-stellar cores have been made by various groups, as reported in Table A.1.

An illustration of the importance of making direct rather than indirect measurements of the degree of fractionation is shown in Fig. 1 (see also Table A.1). The indirect ratios span a broad range of values, from 150 to 510, with an average 283 ± 93, while direct measurements for the same species (nitriles) are less scattered, with an average 295 ± 50. Nevertheless, the 1∕σ²-weighted mean values of indirect ($\langle R\rangle =314\pm$10) and direct measurements ($\langle R\rangle =329\pm$14) are consistent within 1σ⁵. We note that the uncertainty of the Araki et al. (2016) indirect measurement is most likely underestimated. Yet, the weighted averages are in remarkable agreement with the CN/C¹⁵N ratio measured directly in the disk orbiting TW Hya, and also with the prediction of the galactic chemical evolution model of Romano et al. (2017), lending support to the value of 330 adopted in this study for the present-day ¹⁴N:¹⁵N ratio (Hily-Blant et al. 2017, 2019). We also note that the direct and indirect averages are consistent within 1σ.

From Fig. 1, one also notices that the N₂ H⁺ ratios, all obtainedwith the direct method in a sample of starless cores (Bizzocchi et al. 2013; Redaelli et al. 2018), are significantly higher than in other species. Yet, toward the relatively warmer OMC2-FIR4 and I16293E sources, the observed ratios agree well with 330. Indeed, the weighted average of all N₂ H⁺ ratios is 350 ± 34 (we note that unweighted average value is 610 ± 250). Furthermore, when excluding N₂H⁺, the weighted average of direct measurements becomes 325 ± 15.

Fig. A.1 Observed hyperfine manifolds of CN(1–0) at 113.5 GHz (top) and of CN(2–1) at (bottom) toward L1498 as measured with the IRAM-30m telescope (Hily-Blant et al. 2010). For CN(1–0), the five hyperfine lines are shown separately, with their relative intensity indicated at the top. The best-fit model is shown in red, and the residuals are shown below each spectrum. The offsets are indicated on the right of each panel. The x-axis is the LSR velocity, and the intensity is on a main-beam temperature scale.

A.2 Re-evaluation of the CN/C¹⁵N ratio toward L1498

As was emphasized above, direct measurements of isotopic ratios are essential because of the unquantifiable uncertainties introduced when using the double-isotope method (Gerin et al. 2009; Hily-Blant et al. 2017). In the L1498 pre-stellar core, the CN/C¹⁵N ratio was determined using the latter method (Hily-Blant et al. 2013b), in order to circumvent the optical-depth of the main isotopolog. The derived ratio was $R$ = 476 ± 70, where the error was only statistical and did not include the uncertainty intrinsic to the double-isotope method. To obtain a direct measurement, we used the source model of Magalhães et al. (2018) and followed the same methodology and assumptions. The minimization, performed using Markov chains, was done simultaneously on the CN, ¹³CN, and C¹⁵N, spectra obtained along cuts across the core (see Figs. A.1–]A.2). The inner plateau of constant density has a radius of (77${}_{-10}^{+5}$)′′, which is larger than the value of 47 ± 6′′ found for HCN; this appears to indicate that CN is more extended than HCN, but a comprehensive study of L1498, similar to that performed by Daniel et al. (2013) for Barnard1b, is beyond the scope of the present study. The resulting relative abundances are log₁₀(CN/H₂) = $-{8.51}_{-0.07}^{0.08}$, log₁₀(CN/¹³CN) = ${1.59}_{-0.06}^{0.08}$, and log₁₀(C¹⁵N/CN) = ${2.49}_{-0.06}^{0.07}$ (see Fig. A.3), corresponding to CN/H₂ = 3.2 ± 0.5 × 10⁻⁹, CN/C¹⁵N = 310${}_{-40}^{53}$, and CN/¹³CN = 39${}_{-4}^{+8}$.

Fig. A.2 Same as Fig. A.1 for 13CN(1–0) (left) and C15N(1–0) (right).

Fig. A.3 Correlations between the four parameters (from left to right): radius (expressed as a fraction of 8′) of the constant density plateau in the core density profile, log10(CN/H2), log10 (CN/13CN), and log10(CN/C15N). The first 50 MCMC steps were not included.

Appendix B Building the ¹⁵N network

Most reactions can be duplicated in an automatic way, and routine was developed to construct the ¹⁵N- from the ¹⁴N-network. But some reactions must be treated individually. We recall that our network is limited to singly substituted species.

Reactions such as $$\begin{array}{l}\text{N}+\text{crp}\to {\text{N}}^{+}+{\text{e}}^{-}\\ \text{NH}+\text{O}\to \text{OH}+\text{N}\\ \text{NO}+\text{N}\to {\text{N}}_2+\text{O}\end{array}$$

(with “crp” standing for cosmic ray particles) are simply duplicated. This is also the case for adsorption and desorption reactions (see Tables B.1 and B.2) where equal binding energies are adopted for the ¹⁴N and ¹⁵N species.

In other instances, such as $${\text{N}}_2+{\text{He}}^{+}\to {\text{N}}^{+}+\text{N}+\text{He},$$(B.4)

the corresponding ¹⁵N reaction rate coefficient is assumed to be divided equally among the output channels (in this example, each channel is attributed half of the N₂ reaction rate coefficient). Finally, for the following three reactions, which we assume to proceed through proton- and deuteron-hop (Hemsworth et al. 1974), the rate coefficients of each of the ¹⁵N-substituted reaction are taken equal to the corresponding ¹⁴N rate coefficient: $${\text{N}}_2{\text{H}}^{+}+{\text{NH}}_3\to {\text{NH}}_4^{+}+{\text{N}}_2$$(B.5) $$\begin{array}{l}{\text{N}}_2{\text{D}}^{+}+{\text{NH}}_3\to {\text{NH}}_3{\text{D}}^{+}+{\text{N}}_2\\ {\text{N}}_2{\text{H}}^{+}+\text{NO}\to {\text{HNO}}^{+}+{\text{N}}_2.\end{array}$$

The fractionation reactions of TH00 and R15 are listed in Table B.3.

Appendix C Depletion-driven fractionation

The freeze-out timescale for a species of mass $m=\tilde{m}{m}_{\text{H}}$⁶ is $${\tau }_{\text{fo}}={(n_{\text{d}}\pi a_{\text{gr}}{}^2{v}_{\text{th}}S)}^{-1},$$(C.1)

where n_d is the dust density[ cm⁻³] and S is the sticking coefficient, which is taken equal to 1 for all species in our isothermal models at 10 K. For a Maxwell-Boltzmann velocity distribution, the average thermal speed, v_th, is given by ${(8kT/\pi m)}^{1/2}$.

In our model, grains are decomposed into a core, made of refractory material and of constant composition, and mantleswhich result from the competition between adsorption and desorption processes. The grain size, a_gr, is thus the radius of the core, a_c, plus mantles. Noting X_m (X_c) the fractional abundance of a mantle (core) species, and A_m (A_c) its atomic mass number, the grain radius (core and mantles) is given by $$a_{\text{gr}}=a_{\text{c}}{\left(1+\frac{{\sum }_{\text{m}}{X}_{\text{m}}{A}_{\text{m}}}{{\sum }_{\text{c}}{X}_{\text{c}}{A}_{\text{c}}}\frac{{\varrho }_{\text{c}}}{{\varrho }_{\text{m}}}\right)}^{1/3},$$(C.2)

where the core and mantle mass density are ϱ_c = 2 g cm⁻³ and ϱ_m = 1 g cm⁻³, respectively. With the core and mantles composition listed in Table 1, the initial grain radius is a_gr = 1.28a_c or 0.128 μm. We note that the refractory:gas mass ratio (usually called dust:gas mass ratio) is $$Q={\sum }_c{X}_{\text{c}}{A}_{\text{c}}/\mu =0.0061,$$(C.3)

while the grain:gas mass ratio is $$Q_{\text{m}}=\left({\sum }_c{X}_{\text{c}}{A}_{\text{c}}+{\sum }_m{X}_{\text{m}}{A}_{\text{m}}\right)/\mu =0.0094,$$(C.4)

where μ = 1.4 is the mean molecular weight per hydrogen nucleus. The fractional dust abundance is then: $$X_{\text{d}}=\frac{n_{\text{d}}}{n_{\text{H}}}=Q\frac{\mu m_{\text{H}}}{\frac{4}{3}\pi a_{\text{c}}{}^3{\varrho }_{\text{c}}}.$$(C.5)

Numerically, this translates into $$X_{\text{d}}=2.77\times {10}^{-12}{\left(\frac{a_{\text{c}}}{0.1\mu \text{m}}\right)}^{-3}{\left(\frac{{\varrho }_{\text{c}}}{2\text{g}{\text{cm}}^{-3}}\right)}^{-1}\left(\frac{Q}{0.01}\right).$$(C.6)

In our model setup, X_d = 1.69 × 10⁻¹². Given the above, the freeze-out timescale can be written as $${\tau }_{\text{fo}}=x_0{\tilde{m}}^{1/2},$$(C.7)

where $$x_0^{-1}=\frac{3}{4}Q\frac{\mu n_{\text{H}}{m}_{\text{H}}}{{\varrho }_{\text{c}}}\frac{a_{\text{gr}}{}^2}{a_{\text{c}}{}^3}S\sqrt{\frac{8k_{\text{B}}T}{\pi m_{\text{H}}}}.$$(C.8)

Numerically, one obtains: $$\begin{array}{lll}{x}_0[\text{yr}]=\hfill & 7.91\times {10}^4{\left(\frac{n_{\text{H}}}{{10}^4{\text{cm}}^{-3}}\right)}^{-1}{\left(\frac{Q}{0.01}\right)}^{-1}\hfill & \hfill \\ \hfill & \times {\left(\frac{T}{10\text{K}}\right)}^{-1/2}{\left(\frac{a_{\text{gr}}}{0.1\mu \text{m}}\right)}^{-2}{\left(\frac{a_{\text{c}}}{0.1\mu \text{m}}\right)}^3\left(\frac{{\varrho }_{\text{c}}}{2\text{g}{\text{cm}}^{-3}}\right).\hfill & \hfill \end{array}$$(C.9)

Incidentally, with our model setup parameters, the value of x₀ happens to be also 7.91 × 10⁴ yr.

Fig. C.1 Left panel: temporal evolution of the abundances of N, N2, HCN, and NH3 in a constant density model with our fiducial parameters (see Sect. 2.3) and using the R15 rates. Middle panel: same for the isotopic ratios. The full curves are for models with adsorption only, while dashed curves are for models with both adsorption and desorption. Right panel: temporal evolution of the freeze-out timescales τfo of N and N2 and of the dust grain radius agr.

The freeze-out timescale for species such as N₂ or HCN is thus τ_fo ≈ 0.4 Myr, and for atomic nitrogen, this is 0.3 Myr, in agreement with Fig. C.1 showing the results of a constant-density model (n_H = 10⁴ cm ⁻³), but including adsorption. During the isothermal evolution of a starless core, the freeze-out timescale depends on the various quantities as: $$\delta {\tau }_{\text{fo}}/{\tau }_{\text{fo}}=1/2(\delta \tilde{m}/\tilde{m})-2\delta a_{\text{gr}}/a_{\text{gr}}-\delta n_{\text{d}}/n_{\text{d}}.$$(C.10)

In the constant-density model of Fig. C.1, the adsorption timescale depends only on m and a_gr as $$\delta {\tau }_{\text{fo}}/{\tau }_{\text{fo}}=1/2(\delta \tilde{m}/\tilde{m})-2\delta a_{\text{gr}}/a_{\text{gr}}.$$(C.11)

We note that grain coagulation – which would decrease the grain surface area, and thus increase the freeze-out timescale (see Fig. 2 of Flower et al. 2005) – is not taken into account in this simple analysis. Hence, for a given species, the timescale evolves due to the increase in a_gr, which is ≈ 10−15%. Now, for the ¹⁴N and ¹⁵N variants of a given species, $\delta \tilde{m}=1$, and at any instant, the relative difference in the freeze-out timescale reduces to $$\delta {\tau }_{\text{fo}}/{\tau }_{\text{fo}}=1/2(\delta \tilde{m}/\tilde{m})=1/(2\tilde{m}),$$(C.12)

or 1.8% for HCN and N₂, and twice larger for N atoms.

At constant density, the number of a gas-phase species decreases due to adsorption as n = n₀ exp(−t∕τ_fo) and the isotopic ratio of a ¹⁴N-species of mass ${\tilde{m}}_{14}$ decreases with time as $$R(t)=\frac{n_{14}}{n_{15}}=R_0\exp [-t/(2x_0{\tilde{m}}_{14}\sqrt{{\tilde{m}}_{15}})].$$(C.13)

Therefore, differential depletion of ¹⁴N- and ¹⁵N- variants will produce significant fractionation after a typical timescale (see Eq. (5) $${\tau }_{\text{fr}}=2x_0{\tilde{m}}_{14}\sqrt{{\tilde{m}}_{15}}\approx 2x_0{\tilde{m}}^{3/2}=2\tilde{m}{\tau }_{\text{fo}},$$(C.14)

where we approximated $\tilde{m}\approx {\tilde{m}}_{14}\approx {\tilde{m}}_{15}$. For our model parameters, the timescale for depletion-driven fractionation are τ_fr = 8.3 Myr for N atoms, and 23.6 Myr for N₂ or HCN molecules.

Those numbers are a factor of 2.7 smaller in the constant-density models of Loison et al. (2019), who adopted ϱ_c = 3 g cm⁻³, a_gr = 0.1 μm, n_H = 4 × 10⁴ cm ⁻³, and Q = 0.01, for which x₀ = 2.96 × 10⁴ yr. With these parameters, the fractionation timescales become τ_fr = 3.1 and 8.8 Myr for N and HCN, respectively. One thus expects that, after a time 0.4 Myr in the models of Loison et al. (2019) who considered an isotopic ratio of $R_0=$ $R_{\text{tot}}$ = 440, the relative decrease in $R$ should be $\delta R/R_0=0.4/5.1\approx 5\%$ for HCN, or an isotopic ratio HCN/HC¹⁵N ≈ 390. This is in good agreement for HCN (their Fig. 2), while for HNC their model leads to a ratio of ≈ 370. We also note that, for atomic N, the relative variation is ≈13% while these authors report variations by several tens of percent which cannot be explained by differential depletion alone.

From Eq. (C.9, we notice that x₀ depends on a_gr through the ratio a_gr∕a_c which is itself determined by the ratio of the core and mantles properties (compositions, mass density) to the power 1/3 (see Eq. (C.2). Hence this ratio is not expected to vary significantly (as already noticed from Fig. C.1). Also, the mass density of the refractory core is not expected to take values significantly different from our adopted value of 2 g cm⁻³. Thus one finds that x₀ ~ a_c∕(Qn_H). Grain growth would thus increase the freeze-out and the fractionation timescales, as expected since the surface of grains per unit volume would decrease. However, higher values of the density and dust:gas mass ratio, such as toward the midplane of protoplanetary disks (Williams & Best 2014), both shorten the fractionation timescale.

Appendix D Dissociation of N₂ and isotopologs by He⁺

In Fig. D.1, we illustrate, for the R15 rates, the effects on isotopic ratios of differing branching ratios for the reactions of ¹⁴N¹⁴N and ¹⁴N¹⁵N with He⁺, as discussedin Sect. 4.2.3. The modified rates for the two output channels of the N${}_2^{+}$ + He⁺ reaction are listed in Table D.1.

Fig. D.1 Effects on isotopic ratios of differing branching ratios, 1.5 and 0.69, respectively, for the reactions of 14N14N and 14N15N with He+ (Sect. 4.2.3). Results with different branching ratio (full curves) are compared to those obtained with both branching ratios equal to 1.5 (dashed curves). The R15 set of fractionation reactions was used in these calculations.

Impact of cosmic-ray induced ultraviolet radiation

We consider the influence of differential dissociation of ¹⁴N₂ and ¹⁵N¹⁴N by the cosmic-ray induced ultraviolet radiation on the relative abundances of ¹⁴N- and ¹⁵N-containing species, during their evolution to steady-state and in the course of gravitational collapse. As has been demonstrated (Heays et al. 2017, their Table 20), the rates of photodissociation of N₂ and ¹⁵NN, by the cosmic-ray induced radiation, could differ by approximately a factor of three when self-shielding of the principal isotopic form, N₂, is significant. Following these authors, we adopt a rate of photodissociation of N₂ of $39{\zeta }_{{\text{H}}_2}$, while the unshielded rate is 120${\zeta }_{{\text{H}}_2}$ (see their footnote c). In Fig. D.2, we thus compare results obtained with either the same rate or a three-times higher rate for ¹⁵NN. When the latter value is used, molecular and atomic nitrogen are respectively depleted and enriched in ¹⁵N, with consequences for the other nitrogen-containing species that depend on their routes of formation.

Fig. D.2 Influence of the higher rate of photodissociation of 15NN on the relative abundances of N- and 15N-containing species in the course of L-P collapse: using the fractionation reactions of TH00 (panel a), and R15 (panel b). Results with the enhanced rate are shown with full curves, while those using the default rate are in dashed.

References

All Tables

All Figures

Fig. 1 Overview of measurements of the nitrogen isotopic ratio in low- to intermediate-mass cores and Class0. The horizontal line indicates the value of the bulk isotopic ratio of nitrogen adopted in this study ($R_{\text{tot}}$ = 330). Black dots indicate indirect measurements based on the double-isotope method. Colored points refer to direct measurements, with N2H+:N15NH+ outlined in blue. The N2H+:15NNH+ ratios are indicated by an asterisk. See Appendix A.1 and Table A.1 for details.

In the text

	Fig. 2 Depletion of CO with increasing density, as predicted by our L-P collapse models for three values of the dynamical timescale, tdyn, and ac = 0.10 μm, compared with the observational results of Pagani et al. (2007, 2012) for L183. The L-P calculations corresponding to tdyn = tff, where tff is the free-fall time, are shown also for ac = 0.05 μm.
In the text

Fig. 3 Results from chemical models. In each row, left panel: temporal evolution during the quasi-static evolution of the core (step 1), and right panel: evolution of the inner core as a function of the central density during the collapse (step 2). The bulk 14N:15N isotopic ratio is $R_{\text{tot}}$ = 330. (a) Gas-phase abundances, relative to nH, of the main isotopologs of relevant species. (b) 14N/15N abundance ratios when adopting the TH00 set of fractionation reactions. (c) Same as (b) but with the fractionation reactions from R15.

In the text

	Fig. 4 Column density ratios of 14N- and 15N-containing species in the course of L-P collapse (step 2), compared with the values observed in pre-stellar cores through direct (filled symbols) and indirect (open symbols) measurements (see Fig. 4 and Appendix A.1).
In the text

Fig. 5 Isotopic ratio (column density) of N2H+ as a functionof the core density (nH) in the course of L-P collapse (step 2), compared with the values observed in pre-stellar cores (see also Table A.1). The abundance ratios N2 H+:N15NH+ and N2 H+:15NNH+ are shown for three values of the enhancement factor, ɛ, of the rate coefficient for dissociative recombination of N15NH+ (see Sect. 4.2.2), viz. ɛ = 1, 2, and 3. Results with rates from both R15 (full curves) and TH00 (dashed) are shown.

In the text

Fig. A.1 Observed hyperfine manifolds of CN(1–0) at 113.5 GHz (top) and of CN(2–1) at (bottom) toward L1498 as measured with the IRAM-30m telescope (Hily-Blant et al. 2010). For CN(1–0), the five hyperfine lines are shown separately, with their relative intensity indicated at the top. The best-fit model is shown in red, and the residuals are shown below each spectrum. The offsets are indicated on the right of each panel. The x-axis is the LSR velocity, and the intensity is on a main-beam temperature scale.

In the text

	Fig. A.2 Same as Fig. A.1 for 13CN(1–0) (left) and C15N(1–0) (right).
In the text

	Fig. A.3 Correlations between the four parameters (from left to right): radius (expressed as a fraction of 8′) of the constant density plateau in the core density profile, log10(CN/H2), log10 (CN/13CN), and log10(CN/C15N). The first 50 MCMC steps were not included.
In the text

Fig. C.1 Left panel: temporal evolution of the abundances of N, N2, HCN, and NH3 in a constant density model with our fiducial parameters (see Sect. 2.3) and using the R15 rates. Middle panel: same for the isotopic ratios. The full curves are for models with adsorption only, while dashed curves are for models with both adsorption and desorption. Right panel: temporal evolution of the freeze-out timescales τfo of N and N2 and of the dust grain radius agr.

In the text

	Fig. D.1 Effects on isotopic ratios of differing branching ratios, 1.5 and 0.69, respectively, for the reactions of 14N14N and 14N15N with He+ (Sect. 4.2.3). Results with different branching ratio (full curves) are compared to those obtained with both branching ratios equal to 1.5 (dashed curves). The R15 set of fractionation reactions was used in these calculations.
In the text

	Fig. D.2 Influence of the higher rate of photodissociation of 15NN on the relative abundances of N- and 15N-containing species in the course of L-P collapse: using the fractionation reactions of TH00 (panel a), and R15 (panel b). Results with the enhanced rate are shown with full curves, while those using the default rate are in dashed.
In the text