Depletion and fractionation of nitrogen in collapsing cores

Measurements of the nitrogen isotopic ratio in Solar System comets show a constant value, ~140, which is three times lower than the protosolar ratio, a highly significant difference that remains unexplained. Observations of static starless cores at early stages of collapse confirm the theoretical expectation that nitrogen fractionation in interstellar conditions is marginal for most species. Yet, observed isotopic ratios in N2H+ are at variance with model predictions. These gaps in our understanding of how the isotopic reservoirs of nitrogen evolve, from interstellar clouds to comets, and, more generally, to protosolar nebulae, may have their origin in missing processes or misconceptions in the chemistry of interstellar nitrogen. So far, theoretical studies of nitrogen fractionation in starless cores have addressed the quasi-static phase of their evolution such that the effect of dynamical collapse on the isotopic ratio is not known. In this paper, we investigate the fractionation of 14N and 15N during the gravitational collapse of a pre-stellar core through gas-phase and grain adsorption and desorption reactions. The initial chemical conditions, which are obtained in steady state after typically a few Myr, show low degrees of fractionation in the gas phase, in agreement with earlier studies. However, during collapse, the differential rate of adsorption of 14N- and 15N-containing species onto grains results in enhanced 15N:14N ratios, in better agreement with the observations. Furthermore, we find differences in the behavior, with increasing density, of the isotopic ratio in different species. We find that the collapse must take place on approximately one free-fall timescale, based on the CO abundance profile in L183 [see the end in the PDF file]


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 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 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. 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 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 15 N, from the bulk valuein star-forming regions and protoplanetary disks. One is isotopeselective photodissociation of N 2 (Heays et al. 2014(Heays et al. , 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 14 N: 15 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 15 N ratio increases with distance from the central star . This observational result lends support to the hypothesis that selective photodissociation of N 2 (Heays et al. 2014) could be the leading fractionation process in such disks. Disk models predict some enrichment of CN and HCN in 15 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 14 N 2 and 15 N 14 N by the cosmic-ray induced ultraviolet radiation.
The second process is known as chemical mass fractionation (Watson 1976): Heavier isotopic forms have lower zeropoint 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 2 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 13 C and 15 N, can also become enriched when the temperature is sufficiently low (Terzieva & Herbst 2000;Roueff et al. 2015). It has been proposed that the 15 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 15 N ratio in the L1498 pre-stellar core by Magalhães et al. (2018). These authors obtained HCN:HC 15 N=338±28, which is in harmony with the bulk nitrogen isotopic ratio R 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 15 N-rich reservoirs are formed in situ in protoplanetary disks by the action of isotope-selective photodissociation of N 2 , inherited from a nonfractionated interstellar reservoir. However, measurements of the N 2 H + isotopic ratio in dense clouds are a factor of ≈ 2 − 3 times higher than the bulk value R 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 2 H + :N 15 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 2 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 timedependent 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 Section 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.

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 env (R) = core (R/R 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πR 3 0 0 /3 remains constant. In this expression, the initial external radius is denoted by R 0 and the initial mass density by 0 .
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 2 ) ∼ 10 10 cm −3 ( core ∼ 10 −14 g cm −3 ) (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 8 cm −3 . 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 fol-lowed by means of conservation equations of the form where the number density of the species is denoted n(R) [ cm −3 ], and N(R) [ cm −3 s −1 ] 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 where env (R) [g cm −3 ] denotes the mass density, and S (R) [g cm −3 s −1 ] 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.

Chemical network
For the purpose of the present study, the condensed 1 UGAN network (see Paper I) has been extended to include the 15 N isotopologs (see Appendix B). A "sanity" check of the duplicated network was undertaken to verify that the 14 N and 15 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 and hence the form of the rate coefficient for the N + + N 2 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 −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 2 − −− → NH + + H was approximated to the rate with o-H 2 (Dislaire et al. 2012), which is appropriate for the ortho-to-para ratio of H 2 (∼ 10 −3 ) 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 Table 1. Partition of the elemental abundances † in the gas phase and in the grain mantles, as adopted when specifying the initial composition of the gas, prior to its evolution to chemical steady state.

Element
Total Volatile Refractory Gas Mantles Notes.-Numbers in parentheses are powers of 10. † Abundances are expressed relative to the hydrogen nuclei density, n H . § The initial composition of the mantles can be found in Table 3 of Paper I. ‡ The isotopic ratio of bulk nitrogen is the adopted present-day value R tot =330.
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 4 , H 2 O, and NH 3 (see Table B.1). The binding energies of the 14 N and 15 N variants of adsorbed species have equal values.
The initial (deuterated) UGAN network consists of 1196 reactions linking 156 species 2 . After including the 15 N variants of nitrogen-bearing species and the fractionation reactions, the 15 N version of UGAN contains ≈ 1720 reactions and 205 species.

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 2 was ζ H 2 = 3 ×10 −17 s −1 . The refractory core has a radius a c =0.1 µm and a mass density c = 2 g cm −3 , while that of the mantles is m = 1 g cm −3 . The total grain radius, including the mantles, is thus initially a gr = 0.128 µm (see Eq. (C.2)). The 2 Available upon request to the corresponding author.  Fig. 2. The depletion of CO with increasing density, as predicted by our L-P collapse models for three values of the dynamical timescale, t dyn , and a c = 0.10 µm, compared with the observational results of Pagani et al. (2007Pagani et al. ( , 2012 for L183. The L-P calculations corresponding to t dyn = t ff , where t ff is the free-fall time, are shown also for a c = 0.05 µm. constant density models have n H = 10 4 cm −3 , 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 4 cm −3 and temperature T 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 τ ff = (3π/32G 0 ) 1/2 = 4.4 ×10 5 yr.

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 4 cm −3 . 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 2 :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 Section 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 3 , as predicted by our L-P simulation, with the observations of L183 by Pagani et al. (2007Pagani et al. ( , 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. (2007Pagani et al. ( , 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 2 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 freefall timescale.

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 nitrogencontaining 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 ( 14 N, 15 N, O) or singly ionized (C + , S + , Fe + ).
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 pri- 3 The density to which reference is made in the models of gravitational collapse is that of the core (cf. Sect. 2.1). marily 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 2 is smaller by a factor of approximately 25 than the rate used previously, which explains the higher ratio of the steady-state abundances of N 2 and N, n(N 2 )/n(N), in the present case.
We have investigated the effects of varying ill-defined parameters such as the grain-core radius, a c , ζ 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 15 N were generally modest, being somewhat more pronounced when the rate coefficients of TH00 were adopted, rather than those of R15.

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 14 N: 15 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 7 cm −3 . 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 3 and HCN, that are more indirectly related to N. In contrast, the isotopic ratios of N 2 and N 2 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 4 of the 14 NH 3 and 15 NH 3 ices begins to be significant (see Fig. 3). We recall that 14 NH 3 and 15 NH 3 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 2 H + .
It may be noted that, in the course of the collapse, the gasphase 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.

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, 15 N-containing species have slightly slower thermal speeds than their 14 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 6 cm −3 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 τ fo = x 0m 1/2 (see Eq. (C.7)) withm the atomic mass number and x 0 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 0 = 7.9 ×10 4 yr. In a constant density, isothermal, model, and ignoring adsorption/desorption processes x 0 is constant: The difference in adsorption rates for the 14 N and 15 N variants is thus δτ fo /τ fo = 1/2(δm/m) or 1.8% for HCN, and twice larger for N atoms. It follows that the isotopic ratio R would decrease with time as where R 0 is the initial isotopic ratio and the e-folding timescale is given by (see Eq. (C.14)) τ fr ≈ 2x 0m 3/2 ≈ 2mτ fo .
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 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 14 N and 15 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 ff ∼ n −0.5 H . 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 4 to 10 6 cm −3 . As the density increases above ∼ 5 ×10 6 cm −3 , desorption by direct cosmic ray impact liberates 15 N-poor ammonia ices in the gas-phase thus raising the value of R(NH 3 ). 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 largerreaching ∼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.

Comparison with observations
In Fig. 4, we compare the 14 N: 15 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 2 H + is shown in Fig. 5 and is considered further in Sect. 4.2.2.

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. Fig. 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.

N 2 H +
The results of our fiducial model for N 2 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. 15 NNH + and N 15 NH + form mainly in the proton-transfer reaction of 15 NN with H + 3 , for which the rate coefficient has been taken identical to that for the reaction N 2 (H + 3 , N 2 H + )H 2 . Were the rate coefficients of the reactions 15 NN(H + 3 , 15 NNH + )H 2 and 15 NN(H + 3 , N 15 NH + )H 2 to be significantly smaller, the depletion of the 15 N-bearing isotopic forms could be understood. However, the rate coefficients for these barrierless exothermic protontransfer reactions are presumably independent of the isotopic form of nitrogen.
Alternatively, enhanced rates of dissociative recombination (DR) of 15 NNH + and N 15 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, chang-ing their position relative to the dissociation curve. The measurements of Lawson et al. (2011) have shown that, at 300 K, the DR of N 2 H + proceeds faster than that of 15 N 2 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 2 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 (Dahbia Talbi, private communication).
Our model confirms the increase in the isotopic ratio of N 2 H + when the DR rate of 15 NNH + and N 15 NH + is multiplied by a factor of compared to the DR rate of N 2 H + . Fig. 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 15 NH + and 15 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 2 H + /N 15 NH + compared to N 2 H + / 15 NNH + (Kahane et al. 2018;Redaelli et al. 2018).
In Fig. 1, the isotopic ratios in N 2 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, I16293E pre-stellar core is known to be influenced by the nearby I16293-2422 protostar and the kinetic temperature is significantly higher-increasing from 10 K 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 2 H + decreases as the temperature increases, reaching the bulk ratio of 330. We note that the value of R(N 2 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 ∼ 10K toward 1 above ∼ 15 K. Alternatively, it could be that, in these sources, the destruction of N 2 H + is not dominated by DR with electrons but by other reactions, such as with CO.
Finally, from Fig. 4 and Fig. 5, we see that, when > 1, both column density ratios of R(N 2 H + ) increase with density, whereas for other species, the ratios decrease from their steadystate value. This could produce some source-to-source variation of the isotopic ratios of different species.

Ammonia
Ammonia is, after N and N 2 , 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 2 D have been made by Gerin et al. (2009), using the 1 1,1 − 1 0,1 rotational transition of ortho-NH 2 D at 86 GHz. Despite its very weak intensity, the main hyperfine component of the manifold of 15 NH 2 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 15 N is manifest, as can be seen in Fig. 4 (see also We note that the excitation temperature of 15 NH 2 D 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 tot =330. Assuming a 15% decrease in the excitation temperature of the 15 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% to 100% upward; see Gerin et al.), 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 15 N which is not compatible with our results.
In the above context, we note that the charge transfer reaction between He + and N 2 is another process with a known isotopic effect. This reaction has two output channels, N 2 + He + − −− → N + N + + He (7) and has attracted much interest because the energy needed to effect the transition N 2 (X 1 Σ + g ) → N + 2 (C 2 Σ + 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 −9 cm 3 s −1 and the branching ratio between the dissociative and non-dissociative channels is 1.94, for both N 2 isotopologs. However, this branching ratio is known to be sensitive to which isotopolog of N 2 is involved (Govers et al. 1975). Room temperature measurements suggest that the dissociative:non-dissociative branching ratios could be 1.5 for 14 N 14 N and 0.69 for 14 N 15 N, with a total rate coefficient of 1.3 ×10 −9 cm 3 s −1 (Govers, private communication). 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 2 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 15 N can be found in the difference in the branching ratios for the reactions of 14 N 14 N and 14 N 15 N with He + . We note that, with the new branching ratio, R increases also for HCN and CN, while that of N 2 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 2 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.

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 2 and 15 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 15 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 such as 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/ 13 CN ratio is ≈ 40, which is substantially lower than the value of 12 C/ 13 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 2 H + and its 15 N-variants, showing that the observed depletion of N 2 H + in 15 N could be explained if N 15 NH + and 15 NNH + dissociatively recombined more rapidly with electrons than 14 N 2 H + . Furthermore, our results suggest that the enhancement of the dissociative recombination rate of the 15 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 2 H + , most measurement in pre-stellar cores yield ratios that are ∼2-3 times larger than the bulk value R 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 prestellar 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 15 N ratio. Hily-Blant et al. (2017) have noted the large scatter of the directly measured CN/C 15 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 2 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 15 N 14 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. DRF acknowledges support from STFC (ST/L00075X/1), including provision of local computing resources.
Most reactions can be duplicated in an automatic way, and routine was developed to construct the 15 N-from the 14 N-network. But some reactions must be treated individually. We recall that our network is limited to singly substituted species.
Reactions such as (with "crp" standing for cosmic ray particles) are simply duplicated. This is also the case for adsorption and desorption reactions (see Table B.1 and Table B.2) where equal binding energies are adopted for the 14 N and 15 N species.
In other instances, such as the corresponding 15 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 2 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 15 N-substituted reaction are taken equal to the corresponding 14 N rate coefficient: The fractionation reactions of TH00 and R15 are listed in Table B.3. models at 10 K. For a Maxwell-Boltzmann velocity distribution, the average thermal speed, v th , is given by (8kT/πm) 1/2 . In our model, grains are decomposed into a core, made of refractory material and of constant composition, and mantles which 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 where the core and mantle mass density are c = 2 g cm −3 and m = 1 g cm −3 , 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 while the grain:gas mass ratio is where µ = 1.4 is the mean molecular weight per hydrogen nucleus. The fractional dust abundance is then: Numerically, this translates into In our model setup, X d = 1.69 ×10 −12 . Given the above, the freeze-out timescale can be written as Numerically, one obtains: Incidentally, with our model setup parameters, the value of x 0 happens to be also 7.91 ×10 4 yr. The freeze-out timescale for species such as N 2 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 constantdensity model (n H = 10 4 cm −3 ), but including adsorption. During the isothermal evolution of a starless core, the freeze-out timescale depends on the various quantities as: δτ fo /τ fo = 1/2(δm/m) − 2δa gr /a gr − δn d /n d .
(C.10) A&A proofs: manuscript no. fractionation  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 14 N and 15 N variants of a given species, δm = 1, and at any instant, the relative difference in the freeze-out timescale reduces to δτ fo /τ fo = 1/2(δm/m) = 1/(2m), (C.12) or 1.8% for HCN and N 2 , and twice larger for N atoms.
At constant density, the number of a gas-phase species decreases due to adsorption as n = n 0 exp(−t/τ fo ) and the isotopic ratio of a 14 N-species of massm 14 decreases with time as (C.13) Therefore, differential depletion of 14 N-and 15 N-variants will produce significant fractionation after a typical timescale (see Eq. (5)) τ fr = 2x 0m14 m 15 ≈ 2x 0m 3/2 = 2mτ fo , (C.14) where we approximatedm ≈m 14 ≈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 2 or HCN molecules. Those numbers are a factor of 2.7 smaller in the constantdensity models of Loison et al. (2019), who adopted c = 3 g cm −3 , a gr =0.1µm, n H =4 ×10 4 cm −3 , and Q = 0.01, for which x 0 = 2.96 ×10 4 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 tot = 440, the relative decrease in R should be δR/R 0 = 0.4/5.1 ≈ 5% for HCN, or an isotopic ratio HCN/HC 15 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 0 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 −3 . Thus one finds that x 0 ∼ 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. 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.

Appendix D: Dissociation of N 2 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 14 N 14 N and 14 N 15 N with He + , as discussed in Sect. 4.2.3. The modified rates for the two output channels of the N + 2 + He + reaction are listed in Table D.1.

Appendix D.1: Impact of cosmic-ray induced ultraviolet radiation
We consider the influence of differential dissociation of 14 N 2 and 15 N 14 N by the cosmic-ray induced ultraviolet radiation on the relative abundances of 14 N-and 15 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 2 and 15 NN, by the cosmic-ray induced radiation, could differ by approximately a factor of three when self-shielding of the principal isotopic form, N 2 , is significant. Following these authors, we adopt a rate of photodissociation of N 2 of 39ζ H 2 , while the unshielded rate is 120ζ H 2 (see their footnote c). In Fig. D.2, we thus com- pare results obtained with either the same rate or a three-times higher rate for 15 NN. When the latter value is used, molecular and atomic nitrogen are respectively depleted and enriched in 15 N, with consequences for the other nitrogen-containing species that depend on their routes of formation. 2.40E-10 0.00E+00 2.83E+01 15 N + CN − −− → N + C 15 N 0.00E-10 1.667E-01 0.00E+00 (2) N + C 15 N − −− → 15 N + CN 0.00E-10 1.667E-01 2.29E+01 (2) 15 N + C 2 N + − −− → N + C 15 2 N + 3.80E-12 -1.00E+00 0.00E+00 N + C 15 2 N + − −− → 15 N + C 2 N + 3.80E-12 -1.00E+00 3.81E+01 Notes.-(1) Added for consistency with the UGAN network.
(2) Our fiducial network does not take into account the possibility that CN + N could lead to isotopic exchange (see also Wirström & Charnley 2018).