First sample of $\rm N_2H^+$ nitrogen isotopic ratio measurements in low-mass protostars

Context. The nitrogen isotopic ratio is considered an important diagnostic tool of the star formation process, and $N_2H^+$ is particularly important because it is directly linked to molecular nitrogen $N_2$. However, theoretical models still lack to provide an exhaustive explanation for the observed $^{14}N/^{15}N$ values. Aims. Recent theoretical works suggest that the $^{14}N/^{15}N$ behaviour is dominated by two competing reactions that destroy $ N_2H^+$: dissociative recombination and reaction with CO. When CO is depleted from the gas phase, if $N_2H^+$ recombination rate is lower with respect to the $N^{15}NH^+$ one, the rarer isotopologue is destroyed faster. This implies that the $N_2H^+$ isotopic ratio in protostars should be lower than the one in prestellar cores, and consistent with the elemental value of ~440. We aim to test this hypothesis, producing the first sample of $N_2H^+ / N^{15}NH^+$ measurements in low mass protostars. Methods. We observe the $N_2H^+$ and $N^{15}NH^+$ lowest rotational transition towards six young stellar objects in Perseus and Taurus molecular clouds. We model the spectra with a custom python code using a constant $T_{ex}$ approach to fit the observations. We discuss in appendix the validity of this hypothesis. The derived column densities are used to compute the nitrogen isotopic ratio. Results. Our analysis yields an average of $\rm ^{14}N/^{15}N|_{pro} = 420 \pm 15$ in the protostellar sample. This is consistent with the protosolar value of 440, and significantly lower than the average value previously obtained in a sample of prestellar objects. Conclusions. Our results are in agreement with the hypothesis that, when CO is depleted from the gas-phase, dissociative recombinations with free electrons destroy $N^{15}NH^+$ faster than $N_2H^+$, leading to high isotopic ratios in prestellar cores, where CO is frozen on dust grains.


Introduction
Nitrogen is the fifth element in the universe for abundance. In molecular gas, it is believed that its main reservoir is N 2 . However this molecule, similarly to H 2 , is not directly observable in cold environments due to the lack of permanent dipole moment. Diazenylium (N 2 H + ) forms directly from molecular nitrogen through reaction with H + 3 , and it emits bright lines at millimetre wavelengths. It is therefore traditionally considered a good probe of the bulk of nitrogen gas in molecular clouds. Furthermore, nitrogen bearing species are less affected by depletion at very high densities with respect to C-and O-bearing ones, and hence nitrogen hydrides are usually used to trace the innermost parts of star forming regions. Recent works (Caselli et al. 2002;Bergin et al. 2002;Pagani et al. 2007;Redaelli et al. 2019) showed that in fact also N 2 H + and the deuterated N 2 D + experience freeze-out onto dust grains in very dense gas, but to a much lesser extent with respect to CO, HCO + and isotopologues.
Nitrogen is present in two stable isotopes, the main 14 N and the much less abundant 15 N. In the last decades, its isotopic ratio 14 N/ 15 N has been extensively studied both in the Solar System This work is based on observations carried out with the IRAM 30 m Telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). and in the interstellar medium (ISM), since it allows to link the early phases of star formation with more evolved stellar systems like our own, where pristine material is still preserved in icy bodies (e.g Altwegg et al. 2019). In fact, unlikely for instance oxygen isotopic ratios, the nitrogen one is spread over a wide range. In the primitive Solar Nebula, in situ measurements in the Solar wind find 14 N/ 15 N = 440 − 450 (Marty et al. 2011), and a similar value is reported also in the Jupiter atmosphere (Fouchet et al. 2004). Molecular nitrogen in the Earth's atmosphere is on the contrary enriched in 15 N ( 14 N/ 15 N = 272, Nier 1950). Hotspots in carbonaceous chondrites, which are believed to represent some of the most pristine materials in the Solar System, can present isotopic ratios as low as 50 (Bonal et al. 2010). These results led to the idea that at the moment of the birth of our planetary system, multiple nitrogen reservoir were present (see, for instance, Hily-Blant et al. 2017).
Measurements of 14 N/ 15 N in the ISM also yielded a variety of results, depending both on the environments and on the tracer. Usually nitriles (HCN, HNC, and CN) appear enriched in 15 N with respect to the protosolar value: 14 N/ 15 N = 160 − 460 were found in a sample of low-mass protostars by Wampfler et al. (2014); in prestellar cores, Hily-Blant et al. (2013a,b) found values of 140 − 360 in HCN/HNC and 470 − 510 in CN. In high-mass star forming regions, Colzi et al. (2018) found Article number, page 1 of 14 arXiv:2010.05954v1 [astro-ph.GA] 12 Oct 2020 A&A proofs: manuscript no. NFrac_protostars_ACC 14 N/ 15 N = 250 − 650. It is important to highlight, however, that all these measurements were derived using the double-isotope methods, i.e. observing the 13 C-bearing species instead of the more abundant and optically thick 12 C-bearing ones. As a consequence, these results depend on the assumed carbon isotopic ratio 12 C/ 13 C, and they can be up to a factor of 2 off, as recently shown by Colzi et al. (2020). 14 N/ 15 N observations with diazenylium, on the contrary, do not depend on assumptions of this kind, but are made more difficult by a) the intrinsic weakness of the N 15 NH + and 15 NNH + lines, requiring long integration times, especially for low-mass star forming regions; b) the hyperfine anomalies shown by the N 2 H + (1-0) line, in particular in cold and dense environments. As a consequence, studies of the N-fractionation in N 2 H + are rarer. Daniel et al. (2013) found 400 − 600 towards Barnard 1b, which hosts two very young protostellar objects (Gerin et al. 2015). In OMC-2, a protocluster containing several protostars, Kahane et al. (2018) derived 14 N/ 15 N = 190 − 380. Bizzocchi et al. (2013) derived extremely high levels of 15 N depletion in the prestellar core L1544. This result was later confirmed by Redaelli et al. (2018), who found 14 N/ 15 N = 580 − 1000 in a small sample of prestellar sources. In the high mass regime, Fontani et al. (2015) found 14 N/ 15 N = 180 − 1300.
From the theoretical point of view, it is currently difficult to interpret all these observational results. Chemical models are often able to reproduce the 15 N enrichment in nitriles with respect to the elemental value (assumed to be equal to the protosolar value of 440, see e.g. Roueff et al. 2015), even though more recent results seem to be in disagreement (Wirström & Charnley 2018). On the contrary, the case of the high 14 N/ 15 N values shown in N 2 H + is still puzzling.
In the last year, however, two possible solutions to give explanation to the N 2 H + fractionation observations have been proposed. Furuya & Aikawa (2018) suggested that the depletion in 15 N is inherited from the initial stages of the core evolution, when the gas density is low enough that UV photons can penetrate. Molecular nitrogen is selectively photodissociated, since the rare N 15 N is not abundant enough for self-shielding. This leads to a 15 N enrichment in the atomic nitrogen (N) gas. When N freezes out onto dust grains, where it is rapidly transformed in NH 3 ices, the bulk gas results depleted in heavy nitrogen, while the NH 3 ices are enriched. The weak point of this theory is that the selective photodissociation works at low-to-moderate visual extinction (A V 1.5 mag), when ices cannot efficiently form, so that only a small fraction of 15 N can be effectively trapped in ices.
Another possible explanation has been proposed by Loison et al. (2019), who used the 3-phase chemical model Nautilus (Ruaud et al. 2016) to follow the nitrogen chemistry during the star formation process. They suggest that the 15 Nantifractionation seen in N 2 H + can be due to a difference in the dissociative recombination (DR) rates for the different isotopologues. The main N 2 H + destruction pathways, according to their model, are the following: When gas-phase CO abundance is low, such as in cold and dense prestellar cores, where CO is mainly frozen onto dust grain surfaces, the dominant reaction is the DR one. If the DR rate (κ DR ) of N 2 H + is lower than those of N 15 NH + and 15 NNH + , the N 2 H + /N 15 NH + and N 2 H + / 15 NNH + ratios can be significantly larger than the elemental value. When on the other hand CO abundance is high, its reaction with N 2 H + becomes the dominant one and the molecular isotopic ratio decreases back to the elemental value. Recent laboratory work (Lawson et al. 2011) showed that N 2 H + isotopologues exhibit DR rates that vary up to 20% in value, a discrepancy which is of the same order of magnitude of the one hypothesised by Loison et al. (2019), who assumed a DR rate of N 2 H + 50% lower than that of N 15 NH + . The laboratory results are in the direction opposite to the one required by the Loison's theory (i.e. the DR rate of N 2 H + is higher than that of N 15 NH + ). However, no data are available at the ISM low temperatures, since the experiments were performed at room temperature. More recently, Hily-Blant et al. (2020) investigated further this topic. They model the nitrogen isotopic fractionation of several species during the collapse of a core. The novelty of this study is that the chemistry and the dynamics are run simultaneously, whilst in most other works the chemical evolution is simulated in a quasi-static fashion. Their results show that the isotopic dependence of the adsorption rates plays an important role on the evolution of the nitrogen fractionation during the collapse, but the model still fails in reproducing the 14 N/ 15 N values observed in N 2 H + . The paper concludes, in agreement with Loison et al. (2019), that different DR rates for the distinct diazenylium isotopologues could explain the observed values. In particular, they find that κ DR (N 2 H + ) = 2 − 3 × κ DR (N 15 NH + ) is needed to reproduce the observations. There are a few observational hints that point towards this direction: for instance, the 14 N/ 15 N values measured in OMC-2 and Barnard 1b, which hosts YSOs and are therefore warmer, are lower than the ones measured in the prestellar cores sample of Redaelli et al. (2018). In the high mass regime, Fontani et al. (2015) found an anti-correlation between N 2 H + /N 15 NH + and N 2 H + /N 2 D + . Since the deuteration process is highly favoured by the CO depletion (Dalgarno & Lepp 1984), these results also suggest that N 2 H + /N 15 NH + is larger where CO is mostly absent from the gas-phase.
In order to test the hypothesis of Loison et al. (2019), we performed observations of N 2 H + and N 15 NH + (1-0) lines towards six young stellar objects, hence obtaining the first sample of 14 N/ 15 N in diazenylium in low-mass prostostellar cores. If the theory is correct, we expect a lower 14 N/ 15 N with respect to prestellar cores.

Observations
The observations of N 2 H + and N 15 NH + were performed with the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescope, located at Pico Veleta (Spain), during two different sessions (August 2019 and December 2019). The weather was good for 3mm observations during the summer (0.30 < τ 225GHz < 0.70), and very good in winter (τ 225GHz < 0.30). The pointing was frequently checked on bright nearby sources, and found to be usually accurate within 5 . Mercury and Uranus were used as focus calibrators. We used the EMIR E0 frontend, in two different frequency setup, the first centred on the N 2 H + (1-0) frequency (93173.3991 MHz) and the second centred on the N 15 NH + (1-0) frequency (91205.6953 MHz). EMIR was combined with the VESPA backend, set to high spectral resolution (∆ν = 20 kHz, corresponding to 0.06 km s −1 at 93 GHz). The beam size at 93 GHz is θ beam = 27 .
The source sample consists of six young stellar objects (YSOs), five of which belong to the Perseus molecular cloud. The sixth one (L1527) is located in the Taurus cloud. The protostellar cores were selected from the sample of Emprechtinger et al. (2009), who investigated the deuterium fractionation of Notes. (a) This is computed as the ratio between the canonical CO abundance and the observed one, the former being obtained from H 2 column density via X mol (C 18 O)/X mol (H 2 ) = 1.7 × 10 −7 (Frerking et al. 1982) . Notes. (a) The excitation temperature for N 15 NH + (1-0) is fixed to the value derived from the corresponding N 2 H + fit.
N 2 H + in YSOs. The choice was first led by bright emission in the N 2 H + (1-0) transition, maximising the chance of detecting the much weaker N 15 NH + line. Furthermore, we tried to select objects with different recorded dust temperature and CO depletion factor, since these two parameters may play a role in the nitrogen fractionation (see Introduction). Finally, in order to minimise environmental effects, the initial sample consisted of cores from a single molecular cloud. However, due to its high elevation from Pico Veleta, Perseus was not continuously observable, and hence L1527 was added as a filler source. At the cloud distances -295 pc for Perseus and 135 pc for Taurus (Zucker et al. 2018;Schlafly et al. 2014)-the telescope angular resolution corresponds to 0.04 and 0.02 pc, respectively. The integration times were 23 min for N 2 H + and between 3 and 8 h for N 15 NH + , resulting in rms = 20 − 25 mK (main isotopologue) and rms = 4 − 6 mK (rare isotopologue). Table 1 summarises the source sample, including coordinates, dust temperatures (T dust ), and CO depletion factors ( f D ).

Results
The observed spectra are shown in Figures 1 to 6 as black histograms. The top panels present the main isotopologue, whilst the bottom panels show the N 15 NH + transition. The N 2 H + (1-0) is well detected in all sources, with the brightest hyperfine component presenting T MB values higher than 2.5 K (with the exception of L1527, where the central component reaches 1.75 K). The rms of these observations hence translates in a signal-tonoise ratio S NR > 75. The N 15 NH + (1-0) line is much weaker, which justifies our choice of selecting only YSOs with previously detected bright N 2 H + emission. The peak brightness T MB ranges from 20 to 40 mK. The transition is however detected in all sources, with S NR going from 3 (L1527) to S NR = 8 (Barnard 5).

Analysis
The goal of this work is to derive the nitrogen isotopic ratio of N 2 H + in each source. In the assumption that the two isotopologues are co-spatial, 14 N/ 15 N is computed as the ratio of column densities N col (N 2 H + )/N col (N 15 NH + ). Reliable estimates of the latter are therefore needed. In order to estimate N col , we use a custom python code that implements a constant T ex (C_T ex ) fit of the hyperfine structure (see also Melosso et al. 2020). The code works in a similar fashion as the HFS routine of the CLASS package, but the fit parameters are the molecular column density N col and the excitation temperature T ex (instead of the optical depth τ), together with the kinematic parameters (centroid velocity V lsr and line full width at half maximum FWHM). Furthermore, it allows easily to fit multiple velocity components and to set a subset of the parameters as constrains. Accurate frequencies of the individual hyperfine components are taken from Dore et al. (2009) for N 15 NH + and from our own calculations, based on the results from Cazzoli et al. (2012), for N 2 H + . The values of the partition functions come from our own calculations.

Single velocity component analysis
The first implemented fitting strategy is based on a single velocity component analysis. For each source, we first analyse the N 2 H + (1-0) line, which has a total optical depth τ tot >> 1 and is therefore optically thick, hence allowing to derive N col (N 2 H + ) and T ex at the same time, using a single velocity component. The weaker N 15 NH + (1-0) transition is optically thinner (τ tot << 1), and therefore N col (N 15 NH + ) and T ex are degenerate parameters that cannot be fitted simultaneously. In this case, then, we fix T ex to the value derived in the corresponding N 2 H + analysis. Table 2 summarises the best-fit parameters. The computed models are also shown in Figure 1-6 as red curves overlaid to the observations. For most sources, the fit is reasonably good, and the model is able to correctly reproduce the linewidths and the intensity ratios of the different hyperfine components. There are two protostars, however, for which the fit is not optimal, as one can see especially from the zoom-in on the central triplet of the N 2 H + (1-0) transition. L1448 C (figure 3, top panels) presents a broad wing on one side of the spectrum, which is not reproduced by a single velocity component model. Furthermore, the model fails also to reproduce the hyperfine flux ratios. The situation is similar for L1455 A1 (figure 4, top panels), where not only the intensities of the components are not well matched, but there is an evident feature on the red side of the line, which appears as a second, isolated velocity component. This was not detected by Emprechtinger et al. (2009) due to the limited sensitivity of the observations. We discuss in Sec. 4.2 the results of a multiple velocity component fit for these two objects.
A strong assumption made in our approach is that the different hyperfine components share the same excitation temperature (constant T ex assumption), which is often debatable in the case of N 2 H + . In Appendix A we discuss this hypothesis in greater detail. We demonstrate that due to a combination of broader linewidths and lower optical depths, the hyperfine anomalies are expected to be weak, and definitely less important than in prestellar cores.

Multiple velocity component analysis
As seen in Sec. 4.1, the fitting routine is able to reproduce reasonably well the line profiles of the N 2 H + (1-0) line for four sources (IRAS 03282, L1448 IRS2, Barnard 5, and L1527). For the remaining two the obtained fits are worse, since they cannot reproduce the linewidths and the hyperfine intensities correctly. In this Section, we report the modelling of L1448 C and L1455 A1 using a multiple velocity component fit. We will show that the derived isotopic ratios are consistent within uncertainties with the simpler, one-component analysis. In the discussion (Sec. 5) we will therefore focus on the results from Sec. 4.1.

L1448 C
As visible in the top panels of Figure 3, the N 2 H + (1-0) line in L1448 C presents a broad wing on the blue side of the spectrum. This is most probably due to the internal kinematics of the core, which is known to power an extended bipolar outflow (see e.g. Bachiller et al. 1990). The one-component fit routine then models a broad line (the FWHM = 0.8 km s −1 is the highest in the sample, see Table 2), but is not able to reproduce the narrower peaks of the hyperfine components, in particular in the central triplet. We therefore tried to fit the observed spectrum with two components, a narrow one and a broad one.
The fitting code has now in principle eight free parameters. However, since the two components are not strongly separated in velocity, it is not possible to derive the hyperfine intensity ratios for each group, independently. This information is crucial to derive simultaneously the optical depth (and hence the column density) and T ex . As a consequence, for one velocity component these two parameters result degenerate. We therefore had to fix the excitation temperature T ex,1 of one of the two components.
The choice of which value to set for T ex,1 is quite arbitrary, since we do not have other observations that can constrain this parameter. We therefore decided to fix T ex,1 to the value derived with the single-component modelling, following the idea that this should be indicative of at least the average T ex of N 2 H + (1-0) in the source. We would like to highlight, however, that a change in T ex of ≈ 2 K translates in a change of ≈ 15 in the isotopic ratio, since the excitation temperature is the same for both isotopologues.
The code is thus run with seven free parameters. The obtained best-fit values are presented in Table 3, and the resulting fit is shown in Figure 7. The two components are separated by ≈ 0.1 km s −1 , and the linewidth of the broad one is twice as large as the narrow one. The fit is still not ideal, but the hyperfine intensity ratios are better reproduced, as well as the aforementioned broad blue wing. In order to improve further the fit, one would need to model the kinematic structure of the core, which is beyond the scope of this work (see also comments in Appendix A).
Once the N 2 H + (1-0) line is fitted, we can model the N 15 NH + one with two components. The T ex values must be fixed to the N 2 H + values, due to low optical depth reasons (see Sec. 4.1), so the free parameters are in principle six. However, the signal to noise ratio of the data is too low to constrain all of them, and the uncertainties on the best-fit values are 50 − 100%. Hence, in the hypothesis that the two transitions arise from the same medium, we also fixed V lsr and FWHM for each component to the values found from the N 2 H + analysis. This approach is justified by the results from the single-component analysis, which shows that these two parameters are usually consistent within 3σ uncertainties between the two isotopologues. N col (N 15 NH + ) for each component is then the only free parameter.
The results are visually displayed in the bottom panels of Figure 7, and they are summarised in Table 3 column densities are significantly higher than those derived with the single-component fit. This is due to the fact that in the N 2 H + (1-0) line, the T ex,2 of the unconstrained component (which is five times denser than the other) is low (the lowest in the sample), and lower than T ex,1 . The resulting N col (N 2 H + ) is hence higher. However, the derived isotopic ratio is 14 N/ 15 N| 2c = 470 ± 60 and A&A proofs: manuscript no. NFrac_protostars_ACC N col /10 10 cm −2 0.5 ± 0.4 6.4 ± 0.7 6.9 ± 0.8 Notes. (a) Parameter fixed to the best-fit value found with the one-component analysis. (b) Parameter fixed to the corresponding one in N 2 H + analysis.
it is consistent within the uncertainties with the one obtained in Sec. 4.1.

L1455 A1
The recorded N 2 H + (1-0) line towards L1455 A1 presents several spectral features, as visible in Figure 4. The most evident one is found shifted by ≈ +0.9 km s −1 with respect to the main component and it is approximately six times weaker. Since it is present in all three hyperfine groups, it is most likely a second velocity component along the line of sight. Its signal to noise ratio is however insufficient to model it independently, and it is undetected in the N 15 NH + (1-0) transition. We therefore decided to focus on the main, brighter feature only. Similarly to the L1448 C case, this also presents a wing feature on the red side of the spectrum.
As for L1448 C, the fit routine is not able to converge if all parameters of the two components are unconstrained, and we therefore fix one T ex,1 value to the best-fit obtained with a single velocity component. The best fit of the N 2 H + (1-0) line is presented in the top panels of Figure 8. The model is now able to reproduce the hyperfine main beam temperatures within 15 − 20%. The two velocity components are separated in velocity by ≈ 0.5 km s −1 . The weaker one, with a column density one order of magnitude lower than the stronger component, is functional to reproduce the broad wing on the red side of the line. The fit of the N 15 NH + (1-0) line is done with the same approach illustrated for L1448 C. The only free parameters are the column densities of the two components. The results are shown in Figure  8. The best-fit values for both isotopologues are summarised in Table 4.
Unlike for L1448 C, the total column densities of N 2 H + and N 15 NH + derived with the two-component approach are consistent within uncertainties with the ones computed using the simpler, one velocity component method. As a consequence, also the derived isotopic ratio ( 14 N/ 15 N| 2c = 560 ± 50) is consistent with the one presented in Sec. 4.1.

Discussion
Since the multiple-component analysis of L1448 C and L1455 A1 yields isotopic ratios consistent within uncertainties with the ones coming from the single-component analysis, we focus on the results of the latter method (see Sec. 4.1 and Table 2). As already mentioned, the kinematic parameters (V lsr and FWHM) of the N 2 H + and N 15 NH + lines are always consistent within 3σ for each object, supporting the assumption that the emission from the two isotopologues arises from the same spatial region.
The derived excitation temperatures are in the range 5 − 7 K. These values are significantly lower than the observed dust temperatures (23 − 50 K, see Table 1). These transitions are in fact subthermally excited, so that their T ex is lower than the local gas kinetic temperature. Furthermore, the T dust values were derived fitting continuum data at far infrared wavelengths (≈ 60 − 1000 µm), which are more sensitive to the warm/hot component of the dust envelope. Diazenylium, on the contrary, is expected to be destroyed by CO in the innermost part of the protostellar cores, and therefore traces preferentially a colder gas component. The large IRAM beam at 3mm, in addition, makes our observations more sensitive to the lower density envelope.
The last column of Table 2 reports the derived 14 N/ 15 N values. The associated uncertainties are computed via standard error propagation from the column densities values. They are dominated by the uncertainties on N col (N 15 NH + ), which are ≈ 5 − 15%, computed in the assumption that T ex (N 15 NH + ) = T ex (N 2 H + ). Among the six observed objects, four of them present isotopic ratios consistent with or lower than the elemental value 14 N/ 15 N = 440. Two protostars show instead fractionation ratios higher than the elemental value, namely L1448 IRS2 (3.5σ discrepancy) and L1455 (2.1σ). The weighted average across the whole sample is 14 N/ 15 N| pro = 420 ± 15.
In Figure 9 we compare the just obtained measurements of 14 N/ 15 N with the observations towards prestellar cores from Redaelli et al. (2018). It is important to notice that the prestellar sample was analysed using a different approach, i.e. a fully non-LTE radiative transfer analysis. This explains why the uncertainties on the prestellar isotopic ratios are significantly larger then the protostellar ones (30% versus 10% on average). The errors reported in Redaelli et al. (2018) Notes. (a) Parameter fixed to the best-fit value found with the one-component analysis. (b) Parameter fixed to the respective one in N 2 H + analysis.
whilst the uncertainties evaluated in this work are 1σ statistical errors.
The nitrogen isotopic ratios measured in the prestellar sample are systematically higher than in the protostellar one. The weighted average 2 is 14 N/ 15 N| pre = 700 ± 90, which is significantly higher than 14 N/ 15 N| pro .
In the theory of Loison et al. (2019), as mentioned in the introduction, the main parameter that influences the 14 N/ 15 N of N 2 H + is the CO abundance. When CO is heavily depleted in the gas phase due to freeze-out onto the dust grains, N 15 NH + is selectively destroyed by reaction with free electrons and, as a consequence, 14 N/ 15 N increases. Protostellar cores, being warmer than prestellar ones, are expected to present then lower 14 N/ 15 N, since CO starts to evaporate back into the gas phase at a temperature of ≈ 20 K. In Figure 10 we show the relation of the ni-trogen isotopic ratios with the dust temperature (left panel) and the CO depletion factor (right panel). The values of T dust and f D for the protostellar sample are taken from Emprechtinger et al. (2009). Concerning the prestellar cores, we report the central dust temperatures from Redaelli et al. (2018) for L183, L694-2, and L429, whilst for L1544 we use the value indicated by Chacón-Tanarro et al. (2019). Since these T dust values come from modelled profiles, they are shown without errorbars. The CO depletion factors are taken from Crapsi et al. (2005). In Figure 10 (and in Figure 11), for those sources where both N 15 NH + and 15 NNH + were observed, we report the weighted average of the two values. As expected, prestellar cores are significantly colder (T dust < 10 K) than protostellar ones, and they also show the highest values of CO depletion ( f D > 10). Our data seem hence to confirm the hypothesis of Loison et al. (2019) and Hily-Blant et al. (2020). Basing on the data provided by Emprechtinger et al. (2009), we can test the correlation between T dust (and f D ) and 14 N/ 15 N also within our protostar sample. No clear trend is visible. Furthermore, if one excludes from the analysis Barnard 5, which is also the only Class I object in our sample, the correlation appears opposite to the expected one: the higher T dust (and hence the lower f D ), the higher the isotopic ratio. We fit a linear relation to both pairs of datasets, finding:  Fig. 10. Scatterplots of the nitrogen isotopic ratio as a function of T dust (left panel) and CO depletion factor (right panel). Red stars represent protostellar cores from this work, whilst the blue dots represent the prestellar source analysed by Redaelli et al. (2018). T dust and f D values for YSOs are taken from Emprechtinger et al. (2009). For prestellar objects, the T dust values are taken from Chacón-Tanarro et al. (2019) and Redaelli et al. (2018); the values of f D are taken from Crapsi et al. (2005).
In order to assess if there is a correlation among these parameters, and which one holds, better data are needed. In particular, T dust and f D values are derived using observations from different telescopes and hence distinct spatial resolutions. Furthermore, the far-infrared data used to estimate T dust are known to be sensitive to the warmer component of the interstellar medium, whereas N 2 H + could trace the whole warm-to-cold envelope surrounding the young stellar objects.
We also look for a correlation between the nitrogen isotopic ratio and the hydrogen one in N 2 H + , using the results of Emprechtinger et al. (2009) andCrapsi et al. (2005). The data are shown in Figure 11. It is clear that in prestellar cores not only the 14 N/ 15 N ratio, but also the D/H one is higher than in protostellar ones. This is a well known chemical effect, due to the fact that deuteration processes are very effective at high densities and low temperatures (see Ceccarelli et al. 2014, and references therein).
Focusing on the protostellar cores only, a tentative trend is seen, in the sense that the higher the deuteration level, the lower the nitrogen isotopic ratio. A linear fit to the data, excluding L1527, for which only an upper limit in D/H is known, yields: with a Pearson correlation coefficient of p 3 = −0.56. This correlation is opposite to the one found by Fontani et al. (2015) in high mass cores, and also to what is expected from Loison et al. (2019). In fact, the deuteration level is sensitive to the CO depletion, since carbon monoxide can effectively destroy H 2 D + , the precursor of all deuterated gas species. One would therefore expect that when the D/H ratio is high (and hence CO is depleted from the gas phase), 14 N/ 15 N also shows high values. However, we stress that the correlation that we found is weak and only tentative, and we do not draw conclusions based on it, until further data are available.

Conclusions and summary
In this work we have observed the nitrogen isotopic ratio of diazenylium in the first sample of low-mass protostellar cores. We detected N 2 H + and N 15 NH + (1-0) lines above the 3σ level in all six targeted sources. We analysed the observations using a custom code that implements a constant T ex fit of the spectra. As illustrated in Appendix A we do not expect excitation anomalies to be important in protostellar cores, as opposed to prestellar ones. For two sources, a two-component analysis yields better results than the single-component fit. However, the 14 N/ 15 N values obtained with the two approaches are consistent within the errors, and we therefore discuss the results obtained with the one-component method.
The weighted average of the isotopic ratio is 14 N/ 15 N| pro = 420 ± 15, which is consistent within 2σ with the protosolar value of 440. On the contrary, this result is significantly lower than the prestellar values that are in the range 580 − 1000. These results seem to confirm the theory of Loison et al. (2019), according to which when CO is depleted the dominant destruction pathway of diazenylium isotopologues is through dissociative recombination, and that the DR rate for N 2 H + is lower than that of N 15 NH + . This would have profound implications concerning using diazenylium to trace molecular nitrogen, since it means that in cold gas the nitrogen isotopic ratio of N 2 H + and N 2 are not equal.
We tried to verify if a correlation between T dust and f D with 14 N/ 15 N is present within the YSOs sample. We do not find significant correlations, and in fact a weak trend opposite to the expected one is seen. We speculate that better data, tracing exactly the gas emitting the diazenylium lines, are needed to constrain reliably these correlations. For most sources (all but L1448 IRS2), estimations of the kinetic temperature from ammonia (1,1) and (2,2) transitions can be found in the literature (Jijina et al. 1999;Hatchell 2003). These values are in the range 10 − 15 K, hence always higher than the derived T ex , supporting the hypothesis of subthermal excitation for the analysed transitions. Similarly to the T dust case, we do not find significant trend of the nitrogen isotopic ratio with the kinetic temperature. We highlight the need of a complete and coherent sample of ammonia observations, comprising also the (3,3) transition (which is needed to correctly determine the temperature of warm/hot gas), in order to further investigate this point.
Our results provide new, valuable inputs for chemical modellers and also for laboratory studies. In particular, we stress that further measurements of the DR rates of N 2 H + and N 15 NH + in interstellar conditions are needed. These laboratory results would in fact provide definitive proof for the theory of Loison et al. (2019) and of Hily-Blant et al. (2020). We also plan to use the data here presented as a starting point for further investigation at higher spatial resolution using for instance the NOEMA interferometer. This would give us the chance to test how the isotopic ratio varies with temperature and density within the protostellar envelope. So far, such study has been performed only in one high-mass star forming region by Colzi et al. (2019). Since low-mass star forming regions are on average closer, we could have the resolution to really unveil the role of CO freeze out/desorption in driving the nitrogen isotopic ratio. N 2 H + and N 15 NH + in the sample of prestellar cores. Figure A.1 shows the results of the fit with our custom code on the N 2 H + (1-0) transition in L1544, a very evolved prestellar core. The routine is unable to reproduce the different hyperfine components intensities. In particular, with the exception of two components, the intensities are underestimated by 20 − 25%. As a result, the column density and/or the excitation temperature is underestimated, since the total flux is not reproduced. This translates directly in unreliable estimations of N col (N 2 H + ), and also of N col (N 15 NH + ) (since this is estimated using the T ex value of the main isotopologue).
In protostellar cores, however, conditions are different, and they make the excitation anomaly effects less critical. Due to higher temperature and more turbulent motions, lines are in general broader. The average linewidth of the sample presented in Redaelli et al. (2018) is < FWHM > pre = 0.3 km s −1 , whilst it is < FWHM > pro = 0.5 km s −1 for the YSOs here analysed. Since the optical depth is inversely proportional to the line width (see e.g. Eq. 3 in Redaelli et al. 2019), this means that for protostellar cores the total optical depth τ tot , summed over all the hyperfine components, is lower than for prestellar cores. Our code derives also the line total optical depths, which have a mean value for the YSOs of < τ tot > pro ≈ 7.0. In comparison, the optical depths in the prestellar sample are always > 13.0, and their average is < τ tot > pre ≈ 16. Furthermore, broader lines means that selective trapping effects, which contribute to hyperfine anomalies, are less severe. We therefore expect that the C_T ex assumption holds better for protostellar cores than for prestellar ones. This is supported by the fact that our fitting routine is overall able to reproduce the observed spectra, especially for those sources without evidence of multiple velocity components.
We also want to highlight the difficulties that performing a full non-LTE, non-C_T ex analysis carries. In order to implement it, the physical structure of the source in terms of temperature, density, and kinematics is needed. Prestellar cores can easily be modelled as spherically symmetric, especially around their centre, and far-infrared data (such as from Herschel) can be used to characterise the T dust and volume density profiles. The kinematics, usually due to infall or expansion motions, is also often one-dimensional, and can be inferred from spectroscopic data (see for instance Keto et al. 2015). The structure of protostellar cores is on the contrary more complex. Due to the presence of warm/hot dust at the centre and a colder envelope surrounding it, many wavelenghts are needed to constrain the dust thermal emission and thus its temperature and density distribution. Fur-thermore, the presence of molecular outflows and accretion motions make the structure deviate strongly for a 1-D assumption.
In conclusion, modelling the physical structure of our sample of protostellar cores, a fundamental step for non-C_T ex analysis, require extensive datasets, available for all the targeted sources, which is beyond the scopes of this work. At the same time, for the aforementioned reasons, we do not expect significant hyperfine anomalies. We conclude that the assumption of constant T ex holds reasonably.