© E. Redaelli et al. 2020

1 Introduction

Nitrogen is the fifth element in the Universe for abundance. In molecular gas, it is believed that its main reservoir is N₂. However this molecule, similarly to H₂, is not directly observable in cold environments due to the lack of permanent dipole moment. Diazenylium (N₂H⁺) forms directly from molecular nitrogen through reaction with ${\text{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, N-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) show that N₂H⁺ and the deuterated N₂D⁺ also 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 ¹⁴N and the much less abundant ¹⁵N. In recent decades, its isotopic ratio ¹⁴N/¹⁵N has been extensively studied both in the Solar System and in the interstellar medium (ISM) since it allows us 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). Unlike for instance oxygen isotopic ratios, which are fairly constant with respect to the elemental value, nitrogen ratios range widely. In the primitive solar nebula, in situ measurements in the solar wind find ¹⁴N ∕¹⁵N = 440–450 (Marty et al. 2011), and a similar value is also reported in Jupiter’s atmosphere (Fouchet et al. 2004). On the contrary, molecular nitrogen in the Earth’s atmosphere is enriched in ¹⁵N (¹⁴N ∕¹⁵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 momentof the birth of our planetary system, multiple nitrogen reservoir were present (see e.g. Hily-Blant et al. 2017).

Measurements of ¹⁴N/¹⁵N in the ISM also yielded a variety of results, depending on the environments and on the tracer. Usually, nitriles (HCN, HNC, and CN) appear enriched in ¹⁵N with respect to the protosolar value: Wampfler et al. (2014) found ¹⁴N∕¹⁵N = 160–460 in a sample of low-mass protostars; Hily-Blant et al. (2013a,b) found values of 140–360 in HCN/HNC and 470–510 in CN in prestellar cores. In high-mass star forming regions, Colzi et al. (2018) found ¹⁴N ∕¹⁵N = 250–650. It is important to highlight, however, that all these measurements were derived using the double-isotope methods, i.e. observing the ¹³C -bearing species instead of the more abundant and optically thick ¹²C-bearing ones. As a consequence, these results depend on the assumed carbon isotopic ratio ¹²C ∕¹³C, and they can be off by up to a factor of 2, as recently shown by Colzi et al. (2020). ¹⁴N/¹⁵N observations with diazenylium, on the contrary, do not depend on assumptions of this kind, but are made more difficult by the intrinsic weakness of the N¹⁵NH⁺ and ¹⁵NNH⁺ lines, requiring long integration times, especially for low-mass star forming regions, and by the hyperfine anomalies shown by the N₂H⁺ (1-0) line, in particular in cold and dense environments. As a consequence, studies of the N-fractionation in N₂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 ¹⁴N ∕¹⁵N = 190–380. Bizzocchi et al. (2013) derived extremely high levels of ¹⁵N depletion in the prestellar core L1544. This result was later confirmed by Redaelli et al. (2018), who found ¹⁴N/¹⁵N = 580–1000 in a small sample of prestellar sources. In the high-mass regime, Fontani et al. (2015) found ¹⁴N ∕¹⁵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 ¹⁵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 ¹⁴N/¹⁵N values shown in N₂H⁺ is still puzzling.

In the last year, however, two possible solutions to give explanation to the N₂H⁺ fractionation observations have been proposed. Furuya & Aikawa (2018) suggested that the depletion in ¹⁵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 ¹⁵N is not abundant enough for self-shielding. This leads to a ¹⁵N enrichment in the atomic nitrogen (N) gas. When N freezes out onto dust grains, where it is rapidly transformed in NH₃ ices, the bulk gas results depleted in heavy nitrogen, while the NH₃ ices are enriched. The weak point of this theory is that the selective photodissociation works from low to moderate visual extinction values (A_V ≲ 1.5 mag), when ices cannot efficiently form, so that only a small fraction of ¹⁵N can be effectively trapped in ices.

Another possible explanation has been proposed by Loison et al. (2019), who used the three-phase chemical model Nautilus (Ruaud et al. 2016) to follow the nitrogen chemistry during the star formation process. They suggest that the ¹⁵N-antifractionation seen in N₂H⁺ can be due to a difference in the dissociative recombination (DR) rates for the different isotopologues. The main N₂H⁺ destruction pathways, according to their model, are the following: $$\begin{array}{l}{\text{N}}_{\text{2}}{\text{H}}^{\text{+}}\text{+CO}\to {\text{HCO}}^{\text{+}}{\text{+N}}_{\text{2}},\\ {\text{N}}_{\text{2}}{\text{H}}^{\text{+}}\text{+e}\to {\text{N}}_{\text{2}}\text{+H}.\end{array}$$

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 DR reaction is the dominant one. If the DR rate (κ_DR) of N₂H⁺ is lower than those of N¹⁵NH⁺ and ¹⁵NNH⁺, the N₂H⁺/N¹⁵NH⁺ and N₂H⁺/¹⁵NNH⁺ ratios can be significantly higher than the elemental value. On the other hand, when CO abundance is high, its reaction with N₂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₂H⁺ isotopologues exhibit DR rates that vary by up to 20% in value, a discrepancy which is of the same order of magnitude as that hypothesised by Loison et al. (2019), who assumed a DR rate of N₂H⁺ 50% lower than that of N¹⁵NH⁺. The laboratory results are in the direction opposite to the one required by Loison’s theory (i.e. the DR rate of N₂H⁺ is higher than that of N¹⁵NH⁺). However, nodata are available at the ISM low temperatures since the experiments were performed at room temperature. More recently, Hily-Blant et al. (2020) investigated this topic further. 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 to reproduce the ¹⁴N/¹⁵N values observed in N₂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₂H⁺) = 2−3 × κ_DR(N¹⁵NH⁺) is needed to reproduce the observations.

There are a few observational hints that point in this direction; for instance, the ¹⁴N/¹⁵N values measured in OMC-2 and Barnard 1b, which hosts YSOs and are therefore warmer, are lower than those 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₂H⁺/N¹⁵NH⁺ and N₂H⁺/N₂D⁺. Since the deuteration process is highly favoured by the CO depletion (Dalgarno & Lepp 1984), these results also suggest that N₂H⁺/N¹⁵NH⁺ is higher 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₂H⁺ and N¹⁵NH⁺ (1-0) lines towards six young stellar objects, hence obtaining the first sample of ¹⁴N/¹⁵N in diazenylium in low-mass prostostellar cores. If the theory is correct, we expect a lower ¹⁴N/¹⁵N with respect to prestellar cores.

2 Observations

The observations of N₂H⁺ and N¹⁵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 the 3 mm observations during the summer (0.30 < τ_{225 GHz} < 0.70), and very good in winter (τ_{225 GHz} < 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 setups, the first centred on the N₂H⁺ (1-0) frequency (93 173.3991 MHz) and the second centred on the N¹⁵NH⁺ (1-0) frequency (91 205.6953 MHz). EMIR was combined with the VESPA backend, set to high spectral resolution (Δν = 20 kHz, corresponding to 0.06 km s⁻¹ 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 (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 N₂H⁺ in YSOs. The choice was first led by bright emission in the N₂H⁺ (1-0) transition, maximising the chance of detecting the much weaker N¹⁵NH⁺ line. Furthermore, we tried to select objects with different recorded dust temperatures and CO depletion factors 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₂H⁺ and between3 and 8 h for N¹⁵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).

The data were reduced using the GILDAS/CLASS package¹, and they were calibrated into main beam temperature (T_MB) using the tabulated beam efficiency (η_MB = 0.80) and forward efficiency (F_eff = 0.95).

Fig. 1 Observed spectra (black histograms) and best-fit models (red curves) in IRAS 03282. The top panels refer to N2H+, the lower ones to N15NH+. The right panels show the zoomed-in images of the grey-shaded area in the corresponding plots on the left.

3 Results

The observed spectra are shown in Figs. 1–6 as black histograms. The top panels present the main isotopologue, whilst the bottom panels show the N¹⁵NH⁺ transition. The N₂H⁺ (1-0) is clearly 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-to-noise ratio S∕N > 75. The N¹⁵NH⁺ (1-0) line is much weaker, which justifies our choice of selecting only YSOs with previously detected bright N₂H⁺ emission. The peak brightness T_MB ranges from 20 to 40 mK. The transition is however detected in all sources, with S∕N going from 3 (L1527) to S∕N = 8 (Barnard 5).

Fig. 2 Same as in Fig. 1, but for L1448 IRS2.

Fig. 3 Same as in Fig. 1, but for L1448 C.

4 Analysis

The goal ofthis work is to derive the nitrogen isotopic ratio of N₂H⁺ in each source. In the assumption that the two isotopologues are co-spatial, ¹⁴N/¹⁵N is computed as the ratio of column densities N_col(N₂H⁺) /N_col(N¹⁵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 fashion similar to 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 us to easily fit multiple velocity-components and to set a subset of the parameters as constraints. Accurate frequencies of the individual hyperfine components are taken from Dore et al. (2009) for N¹⁵NH⁺ and from our own calculations, based on the results from Cazzoli et al. (2012) for N₂H⁺. The values of the partition functions come from our own calculations.

Fig. 4 Same as in Fig. 1, but for L1455 A1.

Fig. 5 Same as in Fig. 1, but for Barnard 5.

Fig. 6 Same as in Fig. 1, but for L1527.

4.1 Single velocity-component analysis

The first implemented fitting strategy is based on an analysis one velocity-component. For each source, we first analyse the N₂H⁺ (1-0) line, which has a total optical depth τ_tot ≫ 1 and is therefore optically thick, which allows us to derive N_col(N₂H⁺) and T_ex at the same time, using one velocity-component. The weaker N¹⁵NH⁺ (1-0) transition is optically thinner (τ_tot ≪ 1), and therefore N_col(N¹⁵NH⁺) and T_ex are degenerate parameters that cannot be fitted simultaneously. In this case we fix T_ex to the value derived in the corresponding N₂H⁺ analysis. Table 2 summarises the best-fit parameters. The computed models are also shown in Figs. 1–6 as red curves overlaid on 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 can be seen especially from the zoomed-in image of the central triplet of the N₂H⁺ (1-0) transition. L1448 C (Fig. 3, top panels) presents a broad wing on one side of the spectrum, which is not reproduced by a one-component model. Furthermore, the model also fails to reproduce the hyperfine flux ratios. The situation is similar for L1455 A1 (Fig. 4, top panels), where the intensities of the components are not well matched; 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 Sect. 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₂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.

4.2 Multiple velocity-component analysis

As seen in Sect. 4.1, the fitting routine is able to reproduce reasonably well the line profiles of the N₂H⁺ (1-0) line for four sources (IRAS 03282, L1448 IRS2, Barnard 5, and L1527). For the remaining two the obtained fits are not as good, 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 the uncertainties with the simpler one-component analysis. In the discussion (Sect. 5) we therefore focus on the results from Sect. 4.1.

4.2.1 L1448 C

As visible in the top panels of Fig. 3, the N₂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 (FWHM = 0.8 km s⁻¹ 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.

In principle, the fitting code now has eight free parameters. However, since the two components are not strongly separated invelocity, 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₂H⁺ (1-0) in the source. We would like to note, however, that a change in T_ex of ≈ 2 K translates into 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 Fig. 7. The two components are separated by ≈ 0.1 km s⁻¹, 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 is the above-mentioned broad blue wing. In order to improve the fit further, it is necessary 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₂H⁺ (1-0) line is fitted, we can model the N¹⁵NH⁺ line with two components. The T_ex values must be fixed to the N₂H⁺ values, due to the low optical depth (see Sect. 4.1), so the free parameters are six, in principle. However, the S/N 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₂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¹⁵NH⁺) for each component is then the only free parameter.

The results are shown in the bottom panels of Fig. 7, and they are summarised in Table 3. We note that the column density of the weaker component is almost always unconstrained (relative uncertainty: 80%). From the obtained values of column densities, we derive ¹⁴N/¹⁵N as the ratio of the total column densities (summing together the two components). The total column densities are significantly higher than those derived with the single-component fit; in the N₂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₂H⁺) is hence higher. However, the derived isotopic ratio is ¹⁴N/¹⁵N|_2c = 470 ± 60 and it is consistent within the uncertainties with the one obtained in Sect. 4.1.

4.2.2 L1455 A1

The recorded N₂H⁺ (1-0) line towards L1455 A1 presents several spectral features, as visible in Fig. 4. The most evident one is found shifted by ≈ + 0.9 km s⁻¹ 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. However, its S/N is insufficient to model it independently, and it is undetected in the N¹⁵NH⁺ (1-0) transition. We therefore decided to focus on the main brighter feature alone. Similarly to the L1448 C case, this also presents a wing feature on the red side of the spectrum.

As in the L1448 C case, 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₂H⁺ (1-0) line is presented in the top panels of Fig. 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⁻¹. 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¹⁵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 Fig. 8. The best-fit values for both isotopologues are summarised in Table 4.

Unlike the L1448 C case, the total column densities of N₂H⁺ and N¹⁵NH⁺ derived with the two-component approach are consistent within the uncertainties with the values computed using the simpler one velocity-component method. As a consequence, the derived isotopic ratio (¹⁴N/¹⁵N|_2c = 560 ± 50) is consistent with the one presented in Sect. 4.1.

Fig. 7 Obtained best-fit models towards L1448 C using the two velocity-component method for N2H+ (top panels) and N15NH+ (bottom panels), overlaid on the observations (black histograms). The two components are shown as the dashed blue curve and the dash-dotted red curve. The total modelled spectrum is shown as the solid green curve. The right panels show the zoom-in on the central hyperfine group, highlighted with the grey shaded area in the corresponding panels on the left.

5 Discussion

Since the multiple-component analysis of L1448 C and L1455 A1 yields isotopic ratios consistent within the uncertainties with those coming from the one-component analysis, we focus on the results of the latter method (see Sect. 4.1 and Table 2). As already mentioned, the kinematic parameters (V_lsr and FWHM) of the N₂H⁺ and N¹⁵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 and/or 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 3 mm, in addition, makes our observations more sensitive to the lower density envelope.

The last column of Table 2 reports the derived ¹⁴N/¹⁵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¹⁵NH⁺), which are ≈5–15%, computed in the assumption that T_ex(N¹⁵NH⁺) = T_ex (N₂H⁺). Among the six observed objects, four of them present isotopic ratios consistent with or lower than the elemental value ¹⁴N/¹⁵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 ¹⁴N/¹⁵N|_pro = 420 ± 15.

In Fig. 9, we compare the just obtained measurements of ¹⁴N/¹⁵N with the observations towards prestellar cores from Redaelli et al. (2018). It is important to note 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 than the protostellar values (30% versus 10% on average). The errors reported in Redaelli et al. (2018) represent confidence ranges, 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 sample. The weighted average² is ¹⁴N/¹⁵N|_pre = 700 ± 90, which is significantly higher than ¹⁴N/¹⁵N|_pro.

In the theory of Loison et al. (2019), as mentioned in the introduction, the main parameter that influences the ¹⁴N/¹⁵N of N₂H⁺ is the CO abundance. When CO is heavily depleted in the gas phase due to freeze-out onto the dust grains, N¹⁵NH⁺ is selectively destroyed by reaction with free electrons and, as a consequence, ¹⁴N/¹⁵N increases. Protostellar cores, being warmer than prestellar ones, are expected to present lower ¹⁴N/¹⁵N, since CO starts to evaporate back into the gas phase at a temperature of ≈ 20 K. In Fig. 10, we show the relation of the nitrogen 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 error bars. The CO depletion factors are taken from Crapsi et al. (2005). In Fig. 10 (and in Fig. 11), for those sources where both N¹⁵NH⁺ and ¹⁵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 cores, and they also show the highest values of CO depletion (f_D > 10). Our data thus seem to confirm the hypothesis of Loison et al. (2019) and Hily-Blant et al. (2020).

Based on the data provided by Emprechtinger et al. (2009), we can test the correlation between T_dust (and f_D) and ¹⁴N/¹⁵N also within our protostar sample. No clear trend is visible. Furthermore, if we exclude from the analysis Barnard 5, which is the only Class I object in our sample, the correlation appears opposite to the expected one: the higher the value of T_dust (and hence lower f_D), the higher the isotopic ratio. We fit a linear relation to both pairs of datasets: $$\begin{array}{ll}{}^{\text{14}}\text{N/}{}^{\text{15}}{\text{N|}}_{\text{pro}}\hfill & =(10\pm 8)\times T_{\text{dust}}+(130\pm 250),\hfill \\ {}^{\text{14}}\text{N/}{}^{\text{15}}{\text{N|}}_{\text{pro}}\hfill & =(-110\pm 40)\times f_{\text{D}}+(770\pm 130).\hfill \end{array}$$

The Pearson correlation coefficients of p₁ = 0.60 and p₂ = − 0.75, respectively. Even excluding the outlier Barnard 5, the correlation between ¹⁴N/¹⁵N|_pro and T_dust is thus poor, whilst that between ¹⁴N/¹⁵N|_pro and f_D is more significant. However, the correlation coefficient decreases to p₂ = −0.13 if the whole sample is considered.

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₂H⁺ could trace the whole envelope, which ranges from warm to cold, surrounding the young stellar objects.

We also look for a correlation between the nitrogen isotopic ratio and the hydrogen ratio in N₂H⁺, using the results of Emprechtinger et al. (2009) and Crapsi et al. (2005). The data are shown in Fig. 11. It is clear that in prestellar cores the ¹⁴N/¹⁵N value, and also that of D/H, is higher than in protostellar cores. 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 only on the protostellar cores a tentative trend is seen, in the sense that the higher the deuteration level, the lower the nitrogenisotopic ratio. A linear fit to the data, excluding L1527 for which only an upper limit in D/H is known, yields $${}^{\text{14}}\text{N/}{}^{\text{15}}{\text{N|}}_{\text{pro}}\text{=}(\text{−455}\pm \text{480})\times \text{D/H+}(\text{450}\pm \text{60}),$$(5)

with a Pearson correlation coefficient of p₃ = −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). The deuteration level is sensitive to the CO depletion, since carbon monoxide can effectively destroy H₂ D⁺, the precursor of all deuterated gas species. We would therefore expect that when the D/H value is high (and hence CO is depleted from the gas phase) the value of ¹⁴N/¹⁵N also shows high values. However, we note that the correlation that we found is weak and only tentative, and we cannot draw conclusions based on it until further data are available.

Fig. 8 Same as Fig. 7, but for L1455 A1.

Fig. 9 Nitrogen isotopic ratios in N2H+ obtained in this work (red stars) compared with the values obtained towards prestellar cores (blue dots, Redaelli et al. 2018). L1544 and L694-2 present two values because in those sources both N15NH+ and 15NNH+ were detected. The blue dashed curve shows the weighted average of the prestellar sample, and the shaded area the associated uncertainty. The red dot-dashed curve is the average for the protostellar objects. The black solid line shows the elemental value of 14N/15N = 440.

Fig. 10 Scatterplots of the nitrogen isotopic ratio as a function of Tdust (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). Tdust and fD values for YSOs are taken from Emprechtinger et al. (2009). For prestellar objects, the Tdust values are taken from Chacón-Tanarro et al. (2019) and Redaelli et al. (2018); the values of fD are taken from Crapsi et al. (2005).

Fig. 11 Correlation between nitrogen and hydrogen isotopic ratios in the six protostellar cores investigated in this work (red stars) and in the prestellar sample of Redaelli et al. (2018) (blue dots). The D/H of prestellar cores are from Crapsi et al. (2005), whilst the protostellar values are from Emprechtinger et al. (2009).

6 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₂H⁺ and N¹⁵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 one-component fit. However, the ¹⁴N/¹⁵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 ¹⁴N/¹⁵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 seemto confirm the theory of Loison et al. (2019), which states that when CO is depleted the dominant destruction pathway of diazenylium isotopologues is through dissociative recombination, and that the DR rate for N₂H⁺ is lower than that of N¹⁵NH⁺. This would have profound implications for the use of diazenylium to trace molecular nitrogen, since it means that in cold gas the nitrogen isotopic ratio of N₂H⁺ and N₂ are not equal.

We tried to verify whether a correlation between T_dust and f_D with ¹⁴N/¹⁵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 reliably constrain 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 a significant trend of the nitrogen isotopic ratio with the kinetic temperature. We highlight the need of a complete and coherent sample of ammonia observations, also comprising the (3,3) transition (which is needed to correctly determine the temperature of warm or 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 note that further measurements of the DR rates of N₂H⁺ and N¹⁵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 presented in this work 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. To date, such a study has been performed in only 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 truly reveal the role of CO freeze out or desorption in driving the nitrogen isotopic ratio.

Appendix A Constant T_ex (C_ T_ex) assumption

Strictly speaking, a local thermodynamic equilibrium (LTE) analysis assumes that the excitation temperature of all the hyperfine components in all the rotational lines is the same, and that it is equal to the gas kinetic temperature T_K. In our approach, however, we assume that there is one T_ex value shared among all the hyperfine components of the (1-0) transition, which may not coincide with T_K. This assumption, known as C_ T_ex, is more relaxed than LTE, but it still needs justification. Daniel et al. (2006) extensively studied the problem of anomalies in the hyperfine intensity ratios of N₂H⁺, originally reported by Caselli et al. (1995), and found that these effects are more severe at very high volume densities, low temperatures, and high optical depths. These are the physical conditions found in evolved prestellar cores, where n >10⁵ cm⁻³ and T <10 K. This justified the choice made by the authors in Redaelli et al. (2018) to implement a fully non-LTE radiative transfer analysis to model N₂H⁺ and N¹⁵NH⁺ in the sample of prestellar cores.

Figure A.1 shows the results of the fit with our custom code on the N₂H⁺ (1-0) transition in L1544, a very evolved prestellar core. The routine is unable to reproduce the different hyperfine component 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₂H⁺), and also of N_col(N¹⁵NH⁺) (since this is estimated using the T_ex value of the main isotopologue).

Fig. A.1 Modelling of N2H+ (1-0) transitions in the prestellar core L1544. The data, from Redaelli et al. (2018), were observed with the same instrumental setup presented in Sect. 2. The right panel shows the zoomed-in image of the grey-shaded area in the left panel, highlighting the problem with the hyperfine intensity estimation in the central triplet.

In protostellar cores, however, the conditions are different, and they make the excitation anomaly effects less critical. Due to higher temperatureand more turbulent motions, lines are in general broader. The average linewidth of the sample presented in Redaelli et al. (2018) is ${\langle FWHM\rangle }_{\text{pre}}=0.3$km s⁻¹, whilst it is ${\langle FWHM\rangle }_{\text{pro}}=0.5$km s⁻¹ for the YSOs analysed in this paper. 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 also derives the line total optical depths, which have a mean value for the YSOs of ${\langle {\tau }_{\text{tot}}\rangle }_{\text{pro}}\approx 7.0$. In comparison, the optical depths in the prestellar sample are always >13.0, and their average is ${\langle {\tau }_{\text{tot}}\rangle }_{\text{pre}}\approx 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 1D, and can be inferred from spectroscopic data (see e.g. Keto et al. 2015). The structure of protostellar cores is on the contrary more complex. Due to the presence of warm or hot dust at the centre and a colder envelope surrounding it, many wavelengths are needed to constrain the dust thermal emission and thus its temperature and density distribution. Furthermore, the presence of molecular outflows and accretion motions make the structure deviate strongly for a 1D assumption.

In conclusion, modelling the physical structure of our sample of protostellar cores, a fundamental step for non-C_T_ex analysis, requires extensive datasets, available for all the targeted sources, which is beyond the scope of this work. At the same time, for the above reasons, we do not expect significant hyperfine anomalies. We conclude that the assumption of constant T_ex holds reasonably well.

References

All Tables

All Figures

	Fig. 1 Observed spectra (black histograms) and best-fit models (red curves) in IRAS 03282. The top panels refer to N2H+, the lower ones to N15NH+. The right panels show the zoomed-in images of the grey-shaded area in the corresponding plots on the left.
In the text

	Fig. 2 Same as in Fig. 1, but for L1448 IRS2.
In the text

	Fig. 3 Same as in Fig. 1, but for L1448 C.
In the text

	Fig. 4 Same as in Fig. 1, but for L1455 A1.
In the text

	Fig. 5 Same as in Fig. 1, but for Barnard 5.
In the text

	Fig. 6 Same as in Fig. 1, but for L1527.
In the text

Fig. 7 Obtained best-fit models towards L1448 C using the two velocity-component method for N2H+ (top panels) and N15NH+ (bottom panels), overlaid on the observations (black histograms). The two components are shown as the dashed blue curve and the dash-dotted red curve. The total modelled spectrum is shown as the solid green curve. The right panels show the zoom-in on the central hyperfine group, highlighted with the grey shaded area in the corresponding panels on the left.

In the text

	Fig. 8 Same as Fig. 7, but for L1455 A1.
In the text

Fig. 9 Nitrogen isotopic ratios in N2H+ obtained in this work (red stars) compared with the values obtained towards prestellar cores (blue dots, Redaelli et al. 2018). L1544 and L694-2 present two values because in those sources both N15NH+ and 15NNH+ were detected. The blue dashed curve shows the weighted average of the prestellar sample, and the shaded area the associated uncertainty. The red dot-dashed curve is the average for the protostellar objects. The black solid line shows the elemental value of 14N/15N = 440.

In the text

Fig. 10 Scatterplots of the nitrogen isotopic ratio as a function of Tdust (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). Tdust and fD values for YSOs are taken from Emprechtinger et al. (2009). For prestellar objects, the Tdust values are taken from Chacón-Tanarro et al. (2019) and Redaelli et al. (2018); the values of fD are taken from Crapsi et al. (2005).

In the text

	Fig. 11 Correlation between nitrogen and hydrogen isotopic ratios in the six protostellar cores investigated in this work (red stars) and in the prestellar sample of Redaelli et al. (2018) (blue dots). The D/H of prestellar cores are from Crapsi et al. (2005), whilst the protostellar values are from Emprechtinger et al. (2009).
In the text

	Fig. A.1 Modelling of N2H+ (1-0) transitions in the prestellar core L1544. The data, from Redaelli et al. (2018), were observed with the same instrumental setup presented in Sect. 2. The right panel shows the zoomed-in image of the grey-shaded area in the left panel, highlighting the problem with the hyperfine intensity estimation in the central triplet.
In the text