© ESO, 2013

1. Introduction

Characterizing the nitrogen chemistry is of great interest in studies of star-forming regions. Moreover, the possible existence of ¹⁵N-fractionation effects could have strong implications for our understanding of the current solar system. Indeed, large variations have been found in the ¹⁴N/¹⁵N ratio among different solar system bodies (see e.g. Jehin et al. 2009). While this elemental ratio is ~270 on Earth, it is found to be of the order of 450 in the solar wind and Jupiter atmosphere. The latter value is thought to be representative of the protosolar nebula (see e.g. Jehin et al. 2009). Finally, this ratio is in the range 100−150 in cometary ices, as deduced from the analysis of HCN or CN observations. In the coma of comets, most of the CN radicals come from the photodissociation of HCN molecules, and as a consequence, the ¹⁴N/¹⁵N ratios derived from both molecules are necessarily identical. To date, remote estimates of the nitrogen isotopic elemental ratio in comets only rely on these two molecules. In contrast, the analysis of the samples of the Stardust mission did not report such high differences. Indeed, the ¹⁴N/¹⁵N ratio was only found to differ by ~20% with respect to the Earth value (Stadermann et al. 2008). An explanation of such a discrepancy could come from a possible enrichment of the HCN molecule in ¹⁵N atoms, due to zero-point chemistry. This would lead to a lower value for the HC¹⁴N/HC¹⁵N ratio compared to the elemental ¹⁴N/¹⁵N ratio. Since the nitriles are not the main career of nitrogen atoms, such a process could reconcile the values obtained from the HCN or CN observations with the in-situ measurements.

Still today, the nitrogen fractionation chemistry is poorly constrained. Pioneering astrochemical models were computed by Terzieva & Herbst (2000), who predicted that the ¹⁵N-enrichment should remain low enough to prevent detection. However, it was subsequently shown by Charnley & Rodgers (2002) and Rodgers & Charnley (2003, 2004, 2008) that, in the regions where CO is highly depleted, the level of fractionation could be high for some species. These models predict that N₂H⁺ should show the highest degree of fractionation and, more generally, that the nitrogen hydride family should be highly fractioned. More recently, Wirström et al. (2012) has extended the previous studies and included the dependence on the ortho and para symmetries of H₂ in the chemical network. They show that, contrary to what was predicted earlier, the nitriles should show the highest degree of fractionation. This latest result is particularly interesting since it would explain the low values obtained for the HC¹⁴N/HC¹⁵N ratio in comets. Additionally, since the nitrogen hydrides, which are the main carriers of nitrogen atoms, are not highly fractioned, this would enable reconciling the in-situ and remote observations.

Observationally, evidence in favour or against possible fractionation effects are still rare. In part, this comes from few studies having been dedicated to this problem. Previous works concerned NH₃ (Lis et al. 2010), NH₂D (Gérin et al. 2009), HCN (Dahmen et al. 1995; Ikeda et al. 2002; Hily-Blant et al. 2013), HCN and HNC (Tennekes et al. 2006), HNC and CN (Adande & Ziurys 2012), or N₂H⁺ (Bizzocchi et al. 2010, 2013). Moreover, most of the studies performed to date focussed on a single molecular species (or a combination of two species from the same chemical family), which prevents distinguishing source-to-source variations in the elemental ¹⁴N/¹⁵N abundance ratio from fractionation effects related to the species under study. Another uncertainty concerns the methodology used. Indeed, many studies rely on the so-called double isotope method in order to retrieve the column density of the main isotopologue, which leads, to some extent, to an uncertainty linked to the value assumed for the fundamental ¹²C/¹³C ratio. Additionally, this method should only work as long as the ¹³C-fractionation effects remain small for the species considered.

Here, we focus on the interpreting rotational lines observed towards B1b for the ¹³C, ¹⁵N, and D-substituted isotopologues of N₂H⁺, NH₃, CN, HCN, and HNC. This source is located in the Perseus molecular cloud, which has been the subject of many studies. The whole molecular cloud has been mapped, leading to a cartography of its associated extinction (Cernicharo et al. 1985; Carpenter 2000; Ridge et al. 2006). Additionally, large-scale maps have been obtained of the CO isotopologues (Bachiller & Cernicharo 1986; Ungerechts & Thaddeus 1987; Hatchell et al. 2005; Ridge et al. 2006) and of the continuum emission at various wavelengths (Hatchell et al. 2005; Enoch et al. 2006; Jørgensen et al. 2006). From these studies, it appears that this region contains several active star-forming sites, with a few hundred protostellar objects identified. Moreover, it was found that the stars tend to form within clusters. Apart from these extended maps, a large number of dense cores were mapped in other molecular tracers, like NH₃, C₂S (Rosolowsky et al. 2008; Foster et al. 2009), or N₂H⁺ (Kirk et al. 2007).

The source B1b has been the subject of many observational studies because of its molecular richness. Indeed, many molecules were observed for the first time in this object, like HCNO (Marcelino et al. 2009) or CH₃O (Cernicharo et al. 2012). Additionally, B1b shows a high degree of deuterium fractionation and has been associated with first detections of multiply deuterated molecules, such as ND₃ (Lis et al. 2002) or D₂CS (Marcelino et al. 2005). The B1b region was mapped in many molecular tracers, e.g. CS, NH₃, ¹³CO (Bachiller et al. 1990), N₂H⁺ (Huang & Hirano 2013), H¹³CO⁺ (Hirano et al. 1999), or CH₃OH (Hiramatsu et al. 2010; Öberg et al. 2010). This source consists of two cores, B1-bS and B1-bN, which were first identified by Hirano et al. (1999) and initially classified as class 0 objects. However, since none of these cores seems to be the driving source of molecular outflows (Walawender et al. 2005; Hiramatsu et al. 2010) and because of the properties of their spectral energy distribution (Pezzuto et al. 2012), they may correspond to less evolved objects and have been proposed as candidates for the first hydrostatic core stage.

In the present study, we report observations of various isotopologues of N₂H⁺, NH₃, CN, HCN, and HNC toward B1b (Sect. 2). In Sect. 3, the analysis of these observations is described. The physical model of the source is based on the analysis of previously published continuum data (Sect. 3.1), which serves as input for the non-local radiative transfer analysis of the molecular data (Sect. 3.2). The current results are then discussed in Sect. 4 and the concluding remarks are given in Sect. 5.

2. Observations

The list of molecular transitions observed, as well as the parameters that describe the beam of the telescope at the corresponding line frequencies, are summarized in Table 1. This dataset contains previously published data of the NH₃ and NH₂D isotopologues, which are respectively taken from Lis et al. (2010) and Gérin et al. (2009). The central position of the observations is α = 3^h33^m20.90^s, δ = 31°07′34.0″ (J2000).

The observations described in this paper were performed with the IRAM 30 m telescope, except for the p-NH₃ and p-¹⁵NH₃ data, which were obtained with the Green Bank telescope and described by Lis et al. (2010). The IRAM 30 m data were obtained in various observing sessions in November 2001, August 2004, and December 2007 using the heterodyne receivers A100/B100, C150/D150, A230/B230, C270/D270 for spectral line transitions in the 80−115 GHz, 130−170 GHz, 210−260 GHz, and 250−280 GHz frequency ranges, respectively. The tuning of the receivers was optimized to allow simultaneous observations at two different frequencies simultaneously. These receivers were operated in lower sideband, with a rejection of the upper sideband, which is larger than 15dB. The data obtained at frequencies lower than 80 GHz required a specific calibration, as the receivers were operated in a double sideband mode at these low frequencies. The sideband ratio was calibrated using planets (Mars, Uranus, Neptune) as gain calibrators. The spectra were analysed with the Versatile SPectrometer Array, achieving a spectral resolution ranging from 20 kHz for the lower frequency lines (0.07 km s^-1) to 80 kHz (0.09 km s^-1) for the higher frequency transitions like N₂H⁺(3–2). For the weakest lines, the spectra were smoothed to a spectral resolution of ~0.12 km s^-1 to increase the signal-to-noise ratio. We used two observing modes: frequency switching with a frequency throw of ~7.2 MHz at 100 and 150 GHz and of ~14.4 MHz at 230 and 270 GHz, during the 2001 and 2004 observation campaigns (i.e. for the N₂H⁺ and NH₂D maps, as well as for the deuterated species). To limit the contamination of the line profiles by baseline ripples and ensure flat baselines, we also used the wobbler switching mode during the 2007 campaign, using a throw of 4′, since the line and continuum maps showed that most of the emission of the B1b core is concentrated in a region of less than 2′.

The pointing was checked on nearby quasars, while we used Uranus for checking the focus of the telescope. The spectra were calibrated using the applicable forward and main beam efficiencies of 0.95/0.75, 0.93/0.64, 0.91/0.58, and 0.88/0.5 for the A100/B100, C150/D150, A230/B230, and C270/D270 mixers, respectively. Since the radiative transfer model includes the convolution with the telescope beam, we present all spectra in the antenna temperature scale ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}$.

The spectra were processed with the CLASS (Continuum and Line Analysis Single-dish Software) software¹. We removed linear baselines from the wobbler-switched spectra, while higher order polynomials were necessary for the frequency switched data.

3. Continuum and molecular line modelling

The distance of the source is taken to be 235 pc (Hirota et al. 2008). At this distance, an angular size of 10″ corresponds to a physical size of 0.0114 pc.

3.1. Continuum

To derive the gas temperature and H₂ density profiles, we used continuum observations at λ = 350 μm and λ = 1.2 mm, respectively obtained at the CSO and IRAM telescopes. The corresponding continuum maps are reported in Figs. 1 and 2. In these figures, we indicate the positions of the B1-bS and B1-bN cores identified by Hirano et al. (1999), as well as the Spitzer source [EDJ2009] 295 (α = 03^h33^m20.34^s ; δ = 31°07′21.4″ (J2000)) referenced in Jørgensen et al. (2006) and Evans et al. (2009).

The peak position of the 1.2 mm observations is found at α = 03^h33^m21.492^s ; δ = 31°07′32.87″ (J2000), which is offset by (+8″, –1″) with respect to the position of the molecular line observations given in Sect. 2. The 350 μm observations peak at the B1-bS position and are thus offset by ~6″ to the south-south-west, with respect to the 1.2 mm peak position. In the modelling, we assumed that the central position of the density distribution is half-way between the 350 μm and 1.2 mm peaks, which corresponds to the coordinates α = 03^h33^m21.416^s; δ = 31°07′29.79″ (J2000). To characterize the continuum radial behaviour, we performed an annular average of the map around this position. In the B1 region, there are two other cores, labelled B1-d and B1-c by Matthews & Wilson (2002), which are respectively offset by ~80″ to the SW and ~125″ to the NW from the central position. While performing the average, we discarded the points distant by less than 30″ from these two cores. The fluxes defined in this way are represented in Fig. 3, and the error bars represent one standard deviation, σ.

3.1.1. Modelling

To infer the radial structure of the source from these observations, we performed calculations using a ray-tracing code. As input, we fix the H₂ density and temperature profiles which gives, as an output, the intensity as a function of the impact parameter. The model predictions are then compared to the observations by convolving with the telescopes beams, which are approximated by Gaussians of 9″ and 15″ FWHM for the CSO and IRAM telescopes, respectively.

As a starting point, we performed a grid of models assuming that the source is isothermal. We fixed the dust temperature at T_d = 12 K, which is the typical temperature of the region, according to the analysis of NH₃ observations reported by Bachiller et al. (1990) or Rosolowsky et al. (2008). We parameterized the density profile as a series of power laws, i.e. $\begin{array}{ccc}{\mathit{n}}_{\mathrm{0}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}{\mathit{n}}_{\mathrm{0}}& \mathrm{if} \mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{0}}\\ {\mathit{n}}_{\mathit{i}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}{\mathit{n}}_{\mathit{i}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}}\mathrm{\left(}{\mathit{r}}_{\mathit{i}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}}\mathrm{\right)}\mathrm{\times }{\mathrm{(}{\frac{{\mathit{r}}_{\mathit{i}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}}}{\mathit{r}}}^{\mathrm{)}}}^{{\mathit{\alpha }}_{\mathit{i}}}& \mathrm{if} {\mathit{r}}_{\mathit{i}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}}\mathit{<}\mathit{r}\mathit{<}{\mathit{r}}_{\mathit{i}}\end{array}$and we computed a grid of models with the following parameters kept fixed: r₀ = 5″, r₁ = 35″, r₂ = 125″, r₃ = 450″, and α₃ = 3.0. The density profile is thus described by three free parameters, the central H₂ density n₀, and the slopes of the first and second regions, α₁ and α₂. This choice is based on some preliminary models, for which we found that these delimiting radii were accurate in order to reproduce the morphology of the SED. Moreover, the slope of region 3 was found to be poorly constrained, which motivated us to keep it fixed.

The dust absorption coefficient is parametrized according to $\begin{array}{ccc}{\mathit{\kappa }}_{\mathrm{abs}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{\mathit{\kappa }}_{\mathrm{1300}}{\left(\frac{\mathrm{1300}}{\mathit{\lambda }}\right)}^{\mathit{\beta }}& & \end{array}$(3)with λ expressed in μm and κ₁₃₀₀ the absorption coefficient at 1.3 mm. Values for the dust absorption coefficient were determined by Ossenkopf & Henning (1994), for the case of dust grains surrounded by ice mantles. The values were found to be in the range 0.5 < κ₁₃₀₀ < 1.1 cm²/g, where κ₁₃₀₀ is given as a function of the dust mass. We therefore fixed the gas-to-dust mass ratio to 100 and furthermore fixed the absorption coefficient to κ₁₃₀₀ = 0.005 cm²/g, where κ₁₃₀₀ is this time referred to the mass of gas. The results of the models only depend on the product κ₁₃₀₀ × n₀. The grid is then computed by considering β as an additional free parameter.

From the results of the grid, it appears that the best models are found for parameters in the range: 3 × 10⁶ < n₀ < 6 × 10⁶ cm^-3, 2.0 < α₁ < 2.4, 1.2 < α₂ < 1.5, and 1.7 < β < 2.4. However, we could not find a model that would satisfactorily fit the radial profiles at both wavelengths. The origin of the problem is illustrated in Fig. 3. In this figure, the blue curve corresponds to a model that satisfactorily reproduces the 1.2 mm radial profile (i.e. within 30%). It can be seen that this model is able to reproduce the flux at 350 μm, for radii r < 20″. However, for larger radii, the flux is overestimated by a factor 3. Obviously, such a radial dependence cannot be accounted for by a modification of the n₀, κ₁₃₀₀ or β parameters, since changes in these parameters would modify the relative fluxes at both wavelengths and independently of the radius.

A solution to this problem is to introduce a temperature gradient in the model. Indeed, an increase in the dust temperature will affect the fluxes at 350 μm and 1.2 mm differently, producing a larger increase at the shorter wavelength. We thus modified the best model obtained from the grid analysis by introducing a gradient in the dust temperature, from 17 K in the centre to 10 K in the outermost region. The resulting fit is shown by the red curve in Fig. 3. The parameters of this model are n₀ = 3 × 10⁶ cm^-3, α₁ = 2, α₂ = 1.2, and β = 1.9. The dust temperature and H₂ density are represented in Fig. 4. With the density profile defined this way, the peak value of the H₂ column density is N(H₂) = 2.1 × 10²³ cm^-2. The average column density, within a 30″ FWHM beam, is N(H₂) = 7.6 × 10²² cm^-2, in good agreement with the estimate of Johnstone et al. (2010), i.e. N(H₂) = 8.2 × 10²² cm^-2, obtained assuming the same beam. The FHWM of the column density distribution is ~16″. Finally, the current density profile corresponds to enclosed masses of 0.5, 1.8, 4.5, and 13.1 M_⊙ for delimiting radii of 10, 30, 60, and 120″, respectively (~0.01, 0.03, 0.07, and 0.14 pc, assuming a distance to B1 of 235 pc).

3.1.2. Comparison with previous studies

B1b has been the subject of many studies with both ground- and space-based telescopes. In particular, interferometric observations by Hirano et al. (1999) lead to identifying two cores, separated by ~20″ and designated as B1-bS and B1-bN. These objects are respectively offset by (+6″, –7″) and (+5″, +10″) with respect to the central position of our observations (see Sect. 2) and are thus both at a projected distance ~10″. By modelling the SED from 350 μm to 3.5 mm, these authors found that both objects can be characterized by a dust temperature of T_d = 18 K, with a density distribution consistent with n(r) ∝ r^-1.5 and without central flattening. They estimated the bolometric to submm luminosities to L_bol/L_submm ~ 10, well below the limit of 200 for class 0 objects (Andre et al. 1993). Thus, these characteristics would suggest that the two sources are very young class 0 objects.

Recently, Pezzuto et al. (2012) characterized B1-bS and B1-bN by combining Spitzer, Herschel, and ground-based observations. The modelling of the SEDs showed that the two objects are more evolved than prestellar cores, but have not yet formed class 0 objects. In particular, the non-detection of B1-bS at 24 μm and detection at 70 μm place this source as a candidate for the first hydrostatic core stage. However, the bolometric luminosity of this source, estimated as L_bol ~ 0.5 L_⊙, is above the limit of 0.1 L_⊙ that characterizes the maximum luminosity at this evolutionary stage (Omukai 2007). The non-detection of B1-bN at these two wavelengths gives this source the status of a less evolved object. Pezzuto et al. (2012) describe these objects by using a combination of blackbody plus modified blackbody components, which respectively characterize the central object and the surrounding dust envelope. For the two objects, the central temperature was found to be of the order of 30 K and the envelope has a typical temperature ~9 K. However, as noted by these authors, in a realistic model, the emission of the central object will be absorbed by the dusty envelope, which will lead to a stratification in the dust temperature.

Even more recently, Huang & Hirano (2013) have presented 1.1 mm continuum interferometric maps of the B1b region where the B1-bS and B1-bN objects are spatially resolved. For those two objects, they quote the parameters from the upcoming work by Hirano et al. (in prep.) that describe the fit of the SED. B1-bS and B1-bN are described by indexes β = 1.3 ± 0.2 and 1.8 ± 0.4 and T_d = 18.6 ± 1.6 K and 15.6 ± 2.2 K, respectively.

3.2. Molecular excitation modelling

The line radiative transfer (RT) is solved using the short-characteristics method and includes line overlap for the molecules with hyperfine structure. The underlying theory and numerical code are described in Daniel & Cernicharo (2008). For most of the molecules, the spectroscopy is retrieved from the splatalogue database² which centralizes the spectroscopic data available in other databases like JPL³ (Pickett et al. 1998) or CDMS⁴ (Müller et al. 2001, 2005). The rate coefficients are retrieved from the BASECOL⁵ database (Dubernet et al. 2013). The references for both the spectroscopy and collisional rate coefficients are further described in what follows when considering the individual molecules. In the case of linear molecules that have a hyperfine structure resolved by the observations, but for which the collisional rate coefficients just take the rotational structure into account (e.g. the HNC isotopologues), we emulate hyperfine rate coefficients using the IOS-scaled method described in Neufeld & Green (1994). It has been shown by Faure & Lique (2012) that this method of computing hyperfine rate coefficients leads to accurate results when compared to rate coefficients obtained with a more accurate recoupling scheme. In the case of the symmetric or asymmetric tops (i.e. NH₃ and NH₂D), the hyperfine rate coefficients are obtained using the M-random method. Moreover, in the case of collisional systems that consider He as collisional partner, we scale the rate coefficients in order to emulate rate coefficients with H₂, using a scaling factor that depends on the ratio of the reduced masses of the two colliding systems. The same approximation is used to determine the rate coefficients of the rarest isotopologues.

In the modelling, the centre of the density distribution corresponds to the one adopted when modelling the continuum radial profiles and is given is Sect. 3.1. The offset with the position of the molecular line observations is accounted for by computing the molecular spectra at a projected distance of 8″ from the centre. The abundance of a given molecule is obtained dividing the radial profile in a maximum of eight regions. These regions are chosen as a subdivision of the regions chosen to describe the temperature profile. Additionally, for each molecular species, the grid is defined to minimize the abundance variations between two consecutive shells. The outermost radius considered in the model is 450″ for all the molecules. However, in Figs. 5–9 where the abundances that result from the modelling are plotted, we truncate the radius at 150″ since the abundance is constant outside 120″. The abundance in each region is then determined using a Levenberg-Marquardt algorithm, which searches for the best value of the abundances on the basis of a χ² minimization. Typically, a solution is found in 20−30 iterations. In the case of the ¹³C and ¹⁵N isotopologues, the fit is performed simultaneously with the main isotopologue. We force the abundance profiles of the various isotopologues to be identical within a scaling factor. In the case of the D-substituted isotopologues, we require that the abundance in every region is in the range χ(HX)/1000 <χ(DX) < 10 × χ(HX).

For N₂H⁺ and o-NH₂D, the abundance profiles are better constrained because of the availability of emission maps, which is not the case for the other molecules. From the models computed for o-NH₂D and N₂H⁺, it appears, however, that the abundance profile derived using just the observations at the central position should give reliable abundance profiles. The main uncertainty would concern the balance of the molecular abundance between consecutive shells. However, the overall shape of the abundance profile would remain qualitatively unchanged. That the abundance profile is sensitive to observations at a single position comes from the steepness of the density/temperature profile, in conjunction with the availability of multiple lines for a single species with different critical densities.

The error bars of the abundances are obtained by considering each region separately. For each region, the upper and lower limits of the abundance are set so that the χ² does not exceed a given threshold. This way of computing the error bars does not account for the inter-dependence of the different regions, but gives indications of the influence of each region on the emergent intensities. For each molecule, the best abundances are represented in Figs. 5−9 as red lines. The grey regions correspond to the range of abundances delimited by the upper and lower limits. In the figures that compare the observed and modelled spectra, the grey regions correspond to line intensity uncertainties expected from the error bars set on the abundances. From the results of the non-local radiative transfer modelling, we calculate an averaged excitation temperature ${\overline{\mathit{T}}}_{\mathrm{ex}}$ for each molecular line, as well as the line-centre opacity. These quantities are given in Table 1. The averaged excitation temperature is calculated by taking the regions of the cloud that contribute to the emerging line profile into account, the average being given by Eq. (2) of Daniel et al. (2012).

The N₂H⁺ maps of Huang & Hirano (2013) show that the B1-bS and B1-BN cores have different V_LSR, with B1-bS at 6.3 km s^-1 and B1-bN at 7.2 km s^-1. In the present observations, the B1-bN component appears as a redshifted shoulder in the spectra of some species. However, with the current angular resolutions, it is impossible to disentangle the exact contributions of the two emitting cores, since they spatially overlap. Additionally, our model is 1D spherical, which does not allow us to describe the complexity of this region in detail. We thus assumed a single V_LSR ~ 6.7 km s^-1. For some molecules, such as N₂H⁺ or NH₃, the individual hyperfine component have symmetric shapes, and with this assumption, it is possible to obtain an overall good fit of the lines. However, for HCN or HNC, the emission at 7.2 km s^-1 is strong, and this assumption is no longer valid. In such cases, we have chosen to calculate the χ² by considering the blue part of the lines. The channels used to calculate the χ² are shown as black thick histograms in the figures that compare the models to the observations.

In the model, we introduce two regions in terms of infall velocity and turbulent motion. The infall velocity in the region with r < 120″ was set to 0.1 km s^-1 and the turbulent motion to 0.4 km s^-1. These values were chosen since they reproduce the widths of the molecular lines, as well as the asymmetry seen in some molecular lines, e.g. in the HNC or HCN spectra. Outside 120″, the gas is considered to be static, and the turbulent motion is increased to 1.5 km s^-1 in order to reduce the emerging intensities of the HCN and HNC lines. Indeed, with such high turbulence, the lines are absorbed, but no self-absorption features are seen in the synthesized spectra. Thus, this allows to increase the amount of absorbing material to a larger extent than would be possible by simply fitting the observed self-absorption features. Consequently, this allows increasing the column density of the main isotopologues with respect to the ¹³C isotopologues since the latter are less affected by self-absorption effects. This is motivated by the fact that we derive low HCN/H¹³CN and HNC/HN¹³C column density ratios by comparison to isotopic ratios in the local interstellar medium (see Sect. 3.4). Moreover, such an assumption is reasonable in view of the observations related to the Perseus molecular cloud. Indeed, Kirk et al. (2010) show that the turbulent motion in this region follows Larson’s law (Larson 1981), with a turbulent motion that increases with the distance to the core centre and typically reaches a value of 1 km s^-1 on a linear scale of 1 pc. Finally, we emphasize that this outermost region of the model is poorly constrained by the data, and its parameters should be treated with caution and seen as averaged values that would describe the overall foreground material.

Finally, in Table 2, we give the column densities of all the molecular species considered in this study, as well as the column density ratio of the rarest isotopologues with respect to the main ones. Since for a given species, the models are based on various lines observed with different spatial resolutions, we give column densities that simply correspond to an integration along the diameter of the sphere and which are not convolved with any telescope beam.

3.3. N₂H⁺

3.3.1. Spectroscopy and collisional rate coefficients

The spectroscopic coupling constants for the hyperfine structure of N₂H⁺ are taken from Caselli et al. (1995). We adopt the scaling of the rotational constants described in Pagani et al. (2009), where the way we calculate the line strengths is also described. In the case of the ¹⁵N isotopologues, the hyperfine structure due to the ¹⁵N nucleus is not resolved. In this case, we just introduce the spin of the ¹⁴N nucleus into the spectroscopy. The coupling constants are from Dore et al. (2009), and we use the same methodology as Pagani et al. (2009) to derive the line strengths of the two isotopologues.

The collisional rate coefficients are from Daniel et al. (2005), and they consider He as a collisional partner. In the case of the ¹⁵N isotopologues, these rate coefficients are also used, but we apply the recoupling method to one nuclear spin.

3.3.2. Results

In the case of N₂H⁺, in addition to the observations of the J = 1 − 0 and 3−2 lines at the central position, the abundance profile is constrained by taking the J = 1 − 0 map shown in Fig. 11 into account. To do so, we compute the integrated intensity over the hyperfine components for each position of the map. These integrated intensities are then averaged radially, and the χ² is calculated by taking both the central position and the integrated intensity radial profile into account.

The various abundance profiles derived from the modelling are shown in Fig. 5. The comparison of the observations with the modelled line profiles is shown in Fig. 10, and the map of the J = 1 − 0 observations is compared with the model in Fig. 11. The column densities for each isotopologue, as well as the column density ratios, are given in Table 2.

3.3.3. Comparison with previous studies

Kirk et al. (2007) have previously reported N₂H⁺ maps of the whole Perseus molecular cloud and various studies focussed on a sample of prestellar or protostellar cores in this region (Roberts & Millar 2007; Emprechtinger et al. 2009; Johnstone et al. 2010; Friesen et al. 2013). Additionally, a few studies were specifically devoted to the B1b region (Gérin et al. 2001; Huang & Hirano 2013). The column density estimates range from ~10¹³ to 4.5 × 10¹³ cm^-2, and the current estimate of ~3 × 10¹³ cm^-2 thus falls between the previous ones.

Moreover, the N₂D⁺/N₂H⁺ column density ratio was estimated in various studies. This ratio was found to be 0.15 by Gérin et al. (2001), 0.25 by Roberts & Millar (2007), 0.18 by Emprechtinger et al. (2009), and 0.1 by Friesen et al. (2013). More recently, from interferometric observations, Huang & Hirano (2013) estimate the ratio in the B1-bS and B1-bN cores and derive values of 0.13 and 0.27 towards those two positions. In the current study, we obtain a column density ratio ~0.34 at the core centre, a value higher than all the previous estimates.

Part of the differences comes from the assumptions made in the respective analysis. Prior to this work, the analysis was carried out using the assumption of local thermodynamic equilibrium (hereafter referred to as LTE approach, and the details related to this method can be found in e.g. Caselli et al. 2002), and in most cases, the same excitation temperature was assumed for the two isotopologues. In the current model, by considering the excitation temperatures given in Table 2, we can see that this assumption is true within 20% and 10% for the J = 1−0 and J = 3−2 lines, respectively. While this might not lead to substantial errors in the column density estimate, if based on the J = 1−0 line, the error can however be larger for the J = 3−2 line, since the J = 3 energy level is at ~26 K. We refer the reader to the discussion performed in the o-NH₂D case, given in Sect. 3.6, for further details. Finally, from the current model, we infer that the average excitation temperature of N₂D⁺ is higher than that of N₂H⁺, which is a consequence of a more centrally peaked distribution of χ(N₂D⁺) compared to χ(N₂H⁺). On the other hand, the excitation temperatures of the two ¹⁵N isotopologues of N₂H⁺ are lower than for N₂H⁺. This comes from different line-trapping effects. Thus, detailed radiative transfer models are required in order to accurately interpret the line intensities.

Additionally, regardless of the method of analysis, the models show that the respective abundances of N₂D⁺ and N₂H⁺ vary strongly with radius. Indeed, in the innermost region, i.e. r < 15″, we obtain χ(N₂D⁺)/χ(N₂H⁺) greater than unity, while for greater radii, this ratio falls below 0.1. This radial variation of the two isotopologue abundances explains, to some extent, the dispersion in the column density ratio estimates. The reason is that the various studies were not all based on the same observational position, and moreover, the telescope beams were not necessarily similar. The strong radial variations of χ(N₂D⁺)/χ(N₂H⁺) thus make beam dilution effects an important factor when estimating the column density ratio.

The abundance profiles we derive for N₂H⁺ and N₂D⁺ (given in Fig. 5) show that for r < 15′′ (the region that encloses both B1-bS and B1-bN), the abundance of these species are reduced by factors ~100 and 10, respectively, with respect to the maximum abundances that are found at radii 15″ < r < 50″. On the other hand, the interferometric maps of Huang & Hirano (2013) show that the peak intensities of the N₂H⁺(J = 3 − 2) and N₂D⁺(J = 3 − 2) maps are found at the positions of the B1-bN and B1-bS sources. Both findings may seem incompatible. It is thus important to emphasize that even if the N₂H⁺ and N₂D⁺ abundance profiles show signs of depletion in the current analysis, the abundance is still high enough in this depleted region for the lines to be sensitive to the gas density, especially in the excited rotational lines. In other words, despite the low abundances, the J = 3 − 2 lines will still trace the density enhancements. This can be seen when considering the upper boundaries of the abundance represented in Fig. 5. Indeed, the upper limits of the N₂H⁺ and N₂D⁺ abundances in the region r < 15″ are still factors >10 and >5 below the maximum abundances found at 15″ < r < 50″. This is further illustrated by the fact that the synthetic map of the N₂H⁺J = 3−2 emission shown in Fig. 12 has a maximum in the inner 15″ region.

Additionally, we recall that the description of the innermost part of the B1b region is just a crude approximation of the exact geometry of the source, and a quantitative estimate of the amount of depletion would require a more detailed analysis, preferably based on a 3D modelling of the region. Indeed, this point is illustrated in Fig. 12 where a synthetic map of the N₂H⁺J = 3 − 2 emission is shown. This map is obtained assuming an FWHM telescope beam of 4″, which corresponds to the angular resolution of the map shown in Fig. 2c of Huang & Hirano (2013). Compared to this interferometric map, we see that the current model can reproduce, at least qualitatively, the radial variation of the J = 3−2 intensity for radii r > 20″. However, the maximum integrated intensity predicted by the current model is 2 Jy km s^-1/beam, and the integrated intensity is roughly constant within the central 15″. On the other hand, the interferometric map shows that at the positions of B1-bS and B1-bN, the intensity can reach values of ~3 Jy km s^-1/beam, and is above 2 Jy km s^-1/beam close to these sources, in a filament that connects the two objects. The current model is thus reliable enough to describe the envelope on large scales, but fails to mimic the complexity of the B1b region in the vicinity of the B1-bS and B1-bN sources. To account for the geometry of the H₂ density distribution, it would be necessary to adopt a 3D model.

The N₂D⁺/N₂H⁺ column density ratio is mainly influenced by the abundance of the deuterated isotopologues of H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$, whose abundances will mainly be governed by two processes. The first one is linked to the amount of CO that remains in the gas phase. Indeed, if CO is abundant, it will tend to react with the H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ isotopologues, thus limiting their abundances. As a consequence, the N₂H⁺ deuteration will indirectly depend on the size of the dust grains, since smaller grains favour the depletion of CO (Flower et al. 2006a). A second process that governs the N₂H⁺ deuteration is linked to the ortho-to-para ratio of H₂ (Flower et al. 2006b; Pagani et al. 2013). Indeed, at low temperatures, the H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ deuteration is initiated by the reaction H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ + HD → H₂D⁺ + H₂, which has an exothermicity of ~230 K. The backward reaction is thus negligible at low temperatures (T_k ~ 10 K), which leads to a high deuterium fractionation for H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$. However, once created, an efficient destruction of H₂D⁺ is possible if it collides with an o-H₂ molecule, since the fundamental rotational level for this symmetry is ~170 K above the para-H₂ ground state. In other words, the o-H₂ internal energy helps overcome the barrier of the backward reaction. A key parameter for the N₂H⁺ deuteration is thus the H₂ ortho-to-para ratio. In that respect, as demonstrated by Pagani et al. (2013), the timescale at which the collapse proceeds plays an important role. Indeed, if the collapse is fast, the H₂ ortho-to-para conversion induced by the collisions with H⁺ or H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ does not have enough time to reach the steady-state value. Thus, both the initial H₂ ortho-to-para ratio and the age of the cloud (i.e. the time elapsed since the collapse started) will be key parameters when deriving the final N₂D⁺/N₂H⁺ ratio.

By considering the models of Pagani et al. (2013) for the N₂H⁺ deuteration, it appears that the N₂D⁺/N₂H⁺ abundance ratio will strongly vary in the innermost region of the cloud. In their Fig. 13, we can see that in the innermost 0.02 pc, this ratio will vary from a value in the range 1–10 in the centre to around 0.1 at 0.02 pc. The corresponding models have a central density of ~10⁶ cm^-3, which is close to the central density obtained in the current modelling. However, the comparison is only qualitative, since these models are aimed at describing prestellar cores where the central gas temperature is far below the temperature that characterizes the centre of the B1b region. At the distance of B1, the size of 0.02 pc corresponds to an angular size of ~17″. By examining Fig. 5, where the N₂H⁺ and N₂D⁺ abundance profiles are plotted, we indeed derive χ(N₂D⁺)/χ(N₂H⁺) > 1 for r < 15″ and of the order of 0.1 outside this radius. Thus the current results are consistent with what can be expected from the Pagani et al. (2013) model prediction, for a fast collapse model and for a cloud older than 5 × 10⁵ years. This relatively old age is consistent with the presence of a more evolved YSO in the vicinity of the B1b region, while the fast collapse may be related to the triggering effect of outflows from other YSOs in the B1 cloud.

3.4. HCN and HNC

3.4.1. Spectroscopy and collisional rate coefficients

For the HCN isotopologues, we used the HCN/H₂ collisional rate coefficients of Ben Abdallah et al. (2012) that take the hyperfine structure into account. In the particular case of HC¹⁵N, the ¹⁵N nucleus does not lead to a resolved hyperfine structure. Thus, we only considered the rotational energy structure, and the rate coefficients are derived from the hyperfine rate coefficients of Ben Abdallah et al. (2012) by summing over the hyperfine levels. For the HNC isotopologues, we used the rotational collisional rate coefficients for HNC / p-H₂ by Dumouchel et al. (2011). The hyperfine rate coefficients are determined using the IOS-scaled method.

It was shown by Sarrasin et al. (2010) and Dumouchel et al. (2010), when considering the rate coefficients for the HCN and HNC molecules with He, that accounting for the specificities of the two isotopomers leads to significant differences in the collisional rate coefficients. In particular, the HNC and HCN rate coefficients were found to have different propensity rules, since the HNC rate coefficients are higher for odd ΔJ values, while the HCN rate coefficients are higher for even ΔJ values. This conclusion still holds when considering the rate coefficients with H₂ by Ben Abdallah et al. (2012) and Dumouchel et al. (2011). The importance of the new HNC rate coefficients in astrophysical models was emphasized by Sarrasin et al. (2010), who show that an abundance ratio χ(HNC)/χ(HCN) > 1 that would be derived assuming the same rate coefficients for both isomers would instead be compatible with an abundance ratio χ(HNC)/χ(HCN) ~ 1, if one takes the specific rate coefficients for the two isomers into account.

3.4.2. Results

The various isotopologues abundance profiles derived from the modelling are shown in Fig. 6 for HCN and in Fig. 7 for HNC. The comparison of the observations with the modelled line profiles are shown in Fig. 13 for the HCN isotopologues and in Fig. 14 for the HNC ones. In the case of D¹³CN, we could only derive an upper limit for the molecular abundance. The column densities for each isotopologue, as well as the column density ratios, are given in Table 2.

The current results show that the HCN/HNC column density ratio ranges from 2 to 4, depending on the isotopologue considered (omitting the D-substituted isotopologues). If we consider the ¹³C and ¹⁵N isotopologues, we obtain a mean ratio HCN/HNC ~ 2.3. Additionally, we found that HCN shows the lowest degree of deuteration fractionation of all the molecules considered in this work and, in particular, the degree of deuteration is lower than for the HNC isomer. We thus obtained a ratio DCN/DNC ~ 0.6. All of these findings are in good agreement with the predictions of chemical models (see e.g. Fig. 2 of Gérin et al. 2009).

A main uncertainty in the modelling of HCN and HNC comes from the fact that we obtained ¹²C/¹³C ratios, which are respectively factors of 2 and 3 lower than expected for the local ISM (Lucas & Liszt 1998). A main issue while modelling HCN is that the lines have high opacities. Indeed, as we can see in Table 2, the opacities we obtained for the J = 1 − 0, 2−1 and 3−2 lines are ~56, 90, and 37, respectively. With such high opacities, we can expect that part of the HCN molecules will be hidden, since their emission will be subsequently absorbed by the outermost layers of the cloud.

In Fig. 15, we give the contribution of the various layers of the model to the emerging intensity, using the quantity: $\begin{array}{ccc}{\mathit{I}}_{\mathit{v}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{x}}{\mathit{S}}_{\mathit{v}}\mathrm{\left(}\mathit{y}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{v}}\mathrm{\left(}\mathit{y}\mathrm{\right)}}\hspace{0.17em}\mathrm{d}\mathit{y}& & \end{array}$(4)which gives the intensity at a velocity offset v, for an impact parameter that corresponds to the diameter of the cloud, D. The origin of the y-axis corresponds to the outer edge of the cloud where the photons escape. The integration from y = 0 to y = x thus indicates the number of photons created along that path and accounts for the subsequent absorption of the photons by the outermost layers and until the photons escape from the cloud at y = 0.

The emerging intensity is thus I_v(D), and in Fig. 15, we show the quantity (I_v(D) − I_v(x))/I_v(D), which indicates the percentage contribution of each radial point to the emerging intensity. In this figure, we show the contribution of the various regions to the emergent intensity, for the J = 1 − 0 and J = 3 − 2 lines and at various velocity offsets. When considering the abundance of HCN given in Fig. 6, we see that regions depleted in HCN, i.e. 0′′ < r < 15′′ and to a lesser extent 50′′ < r < 120′′, make no significant contribution to the emergent intensity. However, the whole region 15′′ < r < 50′′ contributes to the line profile, and despite the high opacities, the innermost part of this region is still visible. The contribution at a given radius, however, depends on the frequency considered. For example, if we consider the velocity offset v = 0.07 km s^-1 for the J,F = 1, 2–0,1 line, we see that 40% of the emergent flux is created in a thin layer, which is 1″ or 2″ wide, at a distance r ~ 50″.

When the velocity offset with respect to the hyperfine transition increases, the layer that contributes to the emergent flux becomes wider. This can be seen by considering, for example, the velocity offsets v = 0.6 km s^-1 for the J,F = 1,2−0,1 line or v = 5.4 km s^-1 for the J,F = 1,1−0,1 line. Additionally, we can see that the contribution of the backward hemisphere is more important for the blue part of the spectra, while the red part is more significantly influenced by the front hemisphere. This effect, which can be seen when considering the velocity offsets v = −7.4 km s^-1 and −6.6 km s^-1 for the J,F = 1,0−0,1 line, is due to the infall motion we have introduced in the model. If we consider the relative contribution of various layers to the J = 1 − 0 and J = 3 − 2 lines, we can see that the region 120′′ < r < 450′′ affects the J = 1 − 0, line while the J = 3 − 2 line remains mostly unaffected.

The above conclusions are reached by considering a single impact parameter that crosses the centre of the sphere. When comparing the models to the observations, various impact parameters have to be considered in order to perform the convolution with the telescope beam. However, the conclusions reached above remain true for other impact parameters. This explains why in Fig. 6 the abundance in every region has an upper boundary when we fully treat the coupling between the telescope beam and the emission that emerges from the cloud. Finally, given that all the regions where HCN is abundant contribute to some extent to the emerging profile, it is not possible, with the current model, to increase the HCN abundance in a given region without modifying the emerging intensity. With the current assumptions, it is therefore impossible to increase the ¹²C/¹³C ratio to obtain a value that would be consistent with the local ISM value. To perform such a calculation, a multi-dimensional model would be needed.

3.5. NH₃

3.5.1. Spectroscopy and collisional rate coefficients

We used the collisional rate coefficients of Maret et al. (2009), which are calculated by considering p-H₂ (J = 0) as a collisional partner. For both NH₃ and ¹⁵NH₃, accurate rest frequencies for the hyperfine transitions have been measured by Kukolich (1967, 1968). The Hamiltonian was subsequently refined by Hougen (1972). More recently, Coudert & Roueff (2009) have determined a more extensive set of frequencies and line strengths using ab-initio calculations. These values are reported in the splatalogue database. In the NH₃ modelling, we used the values from Coudert & Roueff (2009), except for the (1,1) and (2,2) lines. For these lines, the spectroscopy was retrieved from the CLASS software since the rest frequencies are adapted from the experimental measurements and have a greater accuracy.

When computing the excitation of ¹⁵NH₃, the RT is done by including the unsplit rotational levels of the molecule. The spectrum for the hyperfine structure is then obtained using the formula: $\begin{array}{ccc}{\mathit{T}}_{\mathrm{A}}^{\mathrm{hyp}}\mathrm{\left(}\mathit{v}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}}{\mathit{T}}_{\mathrm{A}}^{\mathrm{rot}}\mathrm{(}\mathit{v}\mathrm{-}{\mathit{v}}_{\mathit{i}}\mathrm{)}\mathrm{\times }\frac{\mathrm{1}\mathrm{-}\exp \left[\mathrm{-}{\mathit{s}}_{\mathit{i}}\mathrm{\times }{\mathit{\tau }}^{\mathrm{rot}}\mathrm{(}\mathit{v}\mathrm{-}{\mathit{v}}_{\mathit{i}}\mathrm{)}\right]}{\mathrm{1}\mathrm{-}\exp \left[\mathrm{-}{\mathit{\tau }}^{\mathrm{rot}}\mathrm{(}\mathit{v}\mathrm{-}{\mathit{v}}_{\mathit{i}}\mathrm{)}\right]}& & \end{array}$(5)where ${\mathit{T}}_{\mathrm{A}}^{\mathrm{hyp}}\mathrm{\left(}\mathit{v}\mathrm{\right)}$ is the antenna temperature of the hyperfine lines at a velocity offset v, ${\mathit{T}}_{\mathrm{A}}^{\mathrm{rot}}\mathrm{\left(}\mathit{v}\mathrm{\right)}$ the one of the rotational line, v_i is the velocity offset of the ith hyperfine component, s_i is its associated line strength, and τ^rot(v) is the opacity of the rotational line. The rest frequencies and line strengths of the (1,1) line correspond to the values given by Kukolich (1967, 1968).

3.5.2. Results

Figure 9 shows the isotopologue abundance profiles derived from the modelling. The comparison of the observations with the modelled line profiles is shown in Fig. 16. The column density for each isotopologue, as well as the column density ratios, is given in Table 2.

3.5.3. Comparison with previous studies

The p-NH₃ emission toward B1b was studied by Bachiller et al. (1990), Rosolowsky et al. (2008), Johnstone et al. (2010), and Lis et al. (2010). The current column density estimate, N(p-NH₃) = 5.5 × 10¹⁴ cm^-2, is in reasonable agreement with the previous works. The p-NH₃ column density was estimated to 1.25 × 10¹⁵ cm^-2 by Bachiller et al. (1990), 5.6 × 10¹⁴ cm^-2 by Rosolowsky et al. (2008), 1.4 × 10¹⁵ cm^-2 by Lis et al. (2010), and 6.5 × 10¹⁴ cm^-2 by Johnstone et al. (2010). Except in the latter study, the authors report the ortho+para NH₃ column density, with an ortho-to-para ratio of 1. The values reported here thus correspond to half of the column densities reported in these works. The column density estimates from the various studies agree within a factor ~2, with the current estimate at the low end of the reported values.

The NH₃ and ¹⁵NH₃ data used in the current study were previously analysed by Lis et al. (2010). The column density ratio we obtained, N(NH₃)/N(¹⁵NH₃) = 300$\begin{array}{c}\mathrm{+}\mathrm{55}\\ -40\end{array}$, agree, within the error bars, with the previous estimate of 334$\begin{array}{c}\mathrm{+}\mathrm{50}\\ -50\end{array}$.

3.6. NH₂D

3.6.1. Collisional rate coefficients

For o-NH₂D, we used the collisional rate coefficients with He of Machin & Roueff (2006). Recently, rate coefficients for ND₂H with H₂ have been determined by Wiesenfeld et al. (2011), and it was shown that for this molecule, the rate coefficients with H₂ are, on average, a factor 10 higher than the rate coefficients with He. We thus expect similar differences for NH₂D and, as a first guess, applied a scaling factor of 10 to the rate coefficients of Machin & Roueff (2006). However, with such an approach, the models fail to reproduce the spectral line shape of the 1_1,1s − 1_0,1a line⁶. In particular, in order to reproduce the overall intensities of the hyperfine lines, it is necessary to adopt an abundance that would produce high opacity for this line (i.e. τ > 10). Subsequently, since the excitation temperature decreases with radius, this would entail the main hyperfine component (i.e. F = 2 − 2) showing a self-absorption feature that is not seen in the observed spectra. Consequently, to reproduce the observed profiles, it would be necessary to increase the density to decrease the opacity needed to fit the hyperfine multiplets.

Since all the other molecules are satisfactorily fitted with the current density profile, we examined in detail the impact of the assumption made on the collisional rate coefficients on the resulting line intensities. To do so, we used the Pearson’s correlation coefficient r defined as $\begin{array}{ccc}\mathit{r}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathit{xy}}}{{\mathit{\sigma }}_{\mathit{x}}\hspace{0.17em}{\mathit{\sigma }}_{\mathit{y}}}& & \end{array}$(6)where σ_xy is the covariance between variables x and y, defined as $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{xy}}\mathrm{=}\frac{\mathrm{1}}{\mathit{N}}\sum _{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{x}}_{\mathit{i}}\mathrm{-}\mathit{x̅}\mathrm{)}\left({\mathit{y}}_{\mathit{i}}\mathrm{-}\mathit{y̅}\right)& & \end{array}$(7)and σ_x and σ_y are the standard deviations with, for example, $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{x}}\mathrm{=}\sqrt{\frac{\mathrm{1}}{\mathit{N}}\sum _{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathit{N}}{\mathrm{(}{\mathit{x}}_{\mathit{i}}\mathrm{-}\mathit{x̅}\mathrm{)}}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(8)In the above expressions, and stand for the mean values of variables x and y. Pearson’s coefficient defines the correlation between two variables. The closest | r | is to 1, the higher the correlation. When r = 0, the two variables are independent.

As previously stated, the ND₂H/He and ND₂H/H₂ rate coefficients differ, on average, by a factor 10. But, when looking closely to the differences, it appears that there is large scatter in the ratio, with values typically in the range 3–30. We thus performed RT calculations where we randomly multiplied the NH₂D/He rate coefficients by a factor in the range 5−20. The statistics are performed on a sample of 100 calculations. To assess the influence of a particular collisional rate coefficient on the resulting line intensity, we used Pearson’s coefficient previously defined, and we identified x to be the integrated line intensity and y the rate coefficient of a particular transition. The result of such a calculation is given in Fig. 17. In this figure, it can be seen that most of the transitions have a correlation coefficient such that | r | < 0.25. These transitions are thus only weakly correlated with the 1_1,1 − 1_0,1 line intensity. Additionally, no transition appears to be strongly correlated with the line intensity, i.e. with | r | → 1. This is because since NH₂D is an asymmetric rotor, the pumping scheme is complex, which implies that the line intensity variations are due to various collisional transitions. The transitions that show the strongest correlation with the line intensity are the 1_1,1 − 1_0,1, 2_0,2 − 1_0,1, 2_1,1 − 1_1,1, and 2_2,1 − 1_0,1 transitions. Except for the 2_1,1 − 1_1,1 transition, an increase in the rate coefficient is correlated with an increase in the line intensities. If we define a correlation coefficient that considers the sum of the rate coefficients for these transitions, weighted by r/ | r |, we obtain a correlation coefficient r = 0.88. This shows that most of the pumping scheme is related to these four transitions.

By comparing the rate coefficients for ND₂H with He and H₂, we can see that for the 1_1,1 − 1_0,1, 2_0,2 − 1_0,1, 2_1,1 − 1_1,1, and 2_2,1 − 1_0,1 transitions, the ratios are ~52, 12, 10, and 13 at 10 K and ~35, 10, 10, and 13 at 25 K. We thus refine the scaling of the NH₂D / He rate coefficients, by still applying an overall factor of 10 to the rates except for the 1_1,1 − 1_0,1, 2_0,2 − 1_0,1, and 2_2,1 − 1_0,1 for which we apply factors of 45, 12, and 13, independently of the temperature. With rate coefficients defined in this way, we are able to reproduce more satisfactorily the o-NH₂D observations shown in Fig. 18.

To test the influence of the above assumptions on the estimate of the NH₂D isotopologue column densities, we performed a few test calculations. More precisely, we ran models where the collisional rate coefficients of the various transitions were multiplied by a random factor between 1/1.5 and 1.5. From one set to the other, we found that the best models were obtained introducing variations in the χ(NH₂D) radial distribution. These variations introduce changes in the column density estimates at most of the order of 50%. On the other hand, the changes in the column density ratios were found to be modest, always less than 10%, with respect to the value of N(NH₂D)/N(¹⁵NH₂D) = 230.

3.6.2. Results

The isotopologue abundance profiles derived from the modelling are shown in Fig. 9. The comparison of the observations with the modelled line profiles are shown in Fig. 18 for both isotopologues, and in Fig. 19 we compare the 1₁₁ − 1₀₁ NH₂D map with the model.

3.6.3. Comparison with previous studies

The NH₂D and ¹⁵NH₂D observations considered in this work were previously analysed by Gérin et al. (2009). In that study, the column densities for both isotopologues were retrieved using LTE, which assumes a uniform distribution for the energy levels throughout the source and the respective populations of the upper and lower energy levels given by the excitation temperature T_ex. Since NH₂D has hyperfine components resolved in frequency, it is possible to obtain an estimate of both the excitation temperature and the total opacity of the rotational line. These parameters were retrieved using the hyperfine method of CLASS. The values obtained were T_ex = 6 K and τ = 5.24, which agree within ~20% with the averaged values T_ex = 6.5 K and τ = 6.4 derived in the present work. Using, say, Eq. (1) of Lis et al. (2002), both sets of values lead to column density estimates that agree within 15%. On the other hand, the hyperfine components of ¹⁵NH₂D are not resolved, and using this methodology, it is not possible to estimate both T_ex and τ from the observations. It was thus assumed that the excitation temperature of ¹⁵NH₂D is equal to that of NH₂D, which then enables the opacity of the line to be estimated from the observed spectra. However, within the current methodology, we obtain an excitation temperature for ¹⁵NH₂D, which is 30% lower than that of NH₂D, which is a consequence of greater line trapping effects for the more opaque isotopologue. This difference in the value of the excitation temperature for ¹⁵NH₂D is important when estimating the line opacity. As an example, using T_ex = 5 K rather than the value of 6 K used in Gérin et al. (2009), leads to an increase in the line opacity by a factor of 1.5. This subsequently leads to re-estimating the ¹⁵NH₂D column density, which is increased by a factor of 2. By examining Eq. (1) of Lis et al. (2002), we can see that the strong dependence of the column density on the excitation temperature comes from the term exp(E_u/k_b T_ex), and is due to the excitation temperature of the 1₁₁ −   − 1₀₁ transition being low by comparison with the energy of the upper level (E_u ~ 20.7 K). This strong dependence of the ¹⁵NH₂D column density on the excitation temperature explains why the ¹⁴N/¹⁵N column density ratio was estimated to be $\mathrm{47}{\mathrm{0}}_{-100}^{\mathrm{+}\mathrm{170}}$ in Gérin et al. (2009), while the current estimate is $\mathrm{23}{\mathrm{0}}_{-55}^{\mathrm{+}\mathrm{105}}$.

The o-NH₂D map in Fig. 19 shows the peak intensity close to the B1-bN core. The emission is elongated towards the position of B1-bS, but the line intensities are lower towards this position. It was found by Huang & Hirano (2013), while discussing the deuterium fractionation of N₂H⁺, that the deuterium enrichment is a factor of two higher towards the B1-bN position compared to B1-bS. The current NH₂D observations suggest a similar trend for this species. The fit of the SED performed by Hirano et al. (in prep.) suggests that B1-bS has a temperature higher by a few Kelvins compared to B1-bN, with a temperature close to 20 K for the former source. The lower level of deuterium enrichment towards B1-bS could thus be related to the fact that the temperature of the dust is close to the temperature at which CO desorbs. A higher amount of CO in the gas phase would drastically reduce the amount deuterated isotopologues of H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ and would thus limit the deuteration fraction.

3.7. CN

3.7.1. Spectroscopy and collisional rate coefficients

Since CN is an open shell molecule, its rotational energy levels are split into fine structure levels due to the interaction with the electron spin. These levels are further split because of the interaction with the nuclear spin of the nitrogen atom. In the case of the ¹³CN and C¹⁵N isotopologues, the half spins of the ¹³C and ¹⁵N nuclei lead to magnetic coupling, which have associated coupling constants large enough so that the degeneracy break can be spectroscopically resolved.

For CN, we have used the CN/He collisional rate coefficients of Lique & Kłos (2011). For either ¹³CN or C¹⁵N, the energy pattern is modified compared to CN, because of the large magnetic coupling constants associated with the ¹³C or ¹⁵N nuclei. Since all these molecules have specific energy structures, it is not possible to directly use the CN/He collisional rate coefficients to obtain the ¹³CN or C¹⁵N rates. In the present study, we therefore derived the rate coefficients for the rarest isotopologues from the diffusion matrix elements obtained for the CN/He system, but accounting for their specificities, while performing the recoupling of the matrix elements. The description is given in Appendix A, and these rates will be made available through the BASECOL database (Dubernet et al. 2013).

In the case of CN, rate coefficients for the CN / p-H₂ system have recently been calculated by Kalugina et al. (2012). They show that individual state-to-state rate coefficients with He or H₂ can differ by up to a factor 3. However, these differences mainly affect the rates of lowest magnitude, while the highest ones scale accordingly to the factor expected from the differences in the reduced masses of the two collisional systems. For consistency between the various isotopologues, we used the CN / He collisional rate coefficients of Lique & Kłos (2011) for the main isotopologue.

3.7.2. Results

The various isotopologue abundance profiles derived from the modelling are shown in Fig. 8. In Fig. 20 we compare the model with the observations of the N = 1 − 0 and N = 2 − 1 spectra of CN and ¹³CN and of the N = 1 − 0 transition of C¹⁵N. In the N = 1 − 0 CN spectrum, a transition due to the doubly deuterated thioformaldehyde isotopologue is seen at a velocity offset v ~ 22 km s^-1 (Marcelino et al. 2005). The column densities for each isotopologue, as well as the column density ratios, are given in Table 2. In this study, we report the first observations made for the various CN isotopologues toward B1b.

4. Discussion

In the upper panel of Fig. 21, we report the abundance profiles of the main isotopologues as a function of the H₂ density. These abundances are obtained from a polynomial fit to the discrete values reported in Figs. 5−9, where all the regions of the model, except the region with r > 120′′ are considered. In Figs. 5 to 9, it can be seen that the molecular abundances have large error bars when the abundance is much below the maximum abundance found for a given molecule. This means that in Fig. 21, we expect the abundances to be accurate within a factor of a few around 10⁵ cm^-3, while around 10⁴ cm^-3 and 10⁶ cm^-3, the error will be typically within an order of magnitude. This can be seen in the middle and bottom panels of Fig. 21 where the abundances used in the fit, as well as their estimated error bars, are reported for the hydrogenated and deuterated isotopologues of NH₃ and N₂H⁺. The abundances reported in this figure can be directly compared to the prediction of the gas-phase chemical model in Fig. 2 of Gérin et al. (2009). While the molecular abundances inferred from the observations agree with the predictions of this chemical model at low densities, the current modelling shows that all the molecules are expected to suffer depletion at densities above 10⁵ cm^-3 and are then considerably lower than predicted by the gas-phase chemistry.

4.1. Molecular depletion and consequences on the N₂H⁺ and NH₃ excitation

Johnstone et al. (2010) analysed the N₂H⁺ and 850 μm continuum observations of Kirk et al. (2007) again by combining them with the p-NH₃ observations of Rosolowsky et al. (2008) for a number of prestellar and protostellar objects in the Perseus molecular cloud. They find that the column density ratio of these two species is relatively constant regardless of the evolutionary stage of the object, with N(p-NH₃)/N(N${}_{\mathrm{2}}\mathrm{H}^{\mathrm{+}}\mathrm{)}\mathrm{=}\mathrm{22}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{10}$. Our results for the B1b region, i.e. N(p-NH₃)/N(N₂H⁺) = 18, are thus in good agreement with the earlier results. Additionally, Johnstone et al. (2010) find that the p-NH₃ abundance is fairly constant in all the prestellar cores, i.e. χ(p-NH₃) ~ 10^-8, while there is larger scatter for the protostellar objects, where it ranges from 10^-8 at low H₂ column densities and decreases by an order of magnitude at higher column densities.

Johnstone et al. (2010) also derived a tight correlation between the excitation temperatures of the N₂H⁺ (J = 1 − 0) and p-NH₃ (1,1) transitions, with a tendency for the p-NH₃ excitation temperature to be higher by 1K. Their sample is consistent with constant excitation temperatures, independently of the central density of the objects. Again, our current results for the B1b region, reported in Table 2; i.e., T_ex(N₂H⁺) ~ 6.8 K and T_ex(p-NH3) ~ 8.1 K, confirm the tendency outlined by Johnstone et al. (2010). As pointed out by these authors, the similarity of the excitation temperatures for these two species is puzzling when considering the critical densities of the respective transitions. More precisely, the excitation temperatures are below the estimated kinetic temperatures of the studied regions which indicates subthermal excitation. Indeed, Johnstone et al. (2010) found that the dust emission at 850 μm can be reproduced assuming T_d ~ 11 K. For densities higher than 10⁴ cm^-3, the kinetic temperature is believed to be coupled well to the dust temperature (Goldsmith 2001). However, considering the differences in the critical densities, the p-NH₃ excitation temperature should be higher than for N₂H⁺, if the molecules trace the same volume of gas. As pointed out by these authors, similar excitation temperatures would require densities that are high enough for the lines to be thermalized. They thus stress that the excitation temperatures for both species may have been underestimated. Our current results shed some light on the puzzling behaviour reported by Johnstone et al. (2010). Indeed, considering the spatial distribution of χ(p-NH₃) and χ(N₂H⁺) shown in Figs. 5 and 9, we see that N₂H⁺ has an abundance distribution that is more centrally peaked than p-NH₃. Since the innermost region of the cloud is warmer and denser, this enhances the mean excitation temperature of N₂H⁺ with respect to that of p-NH₃, thus reducing the impact of the differing critical densities. Moreover, our models show that both molecules are removed from the gas phase in the innermost part of the cloud. If we consider that both molecules suffer depletion above a given density, it means that the properties derived from both molecules will characterize the gas that is below this density threshold. On the other hand, the dust will still trace the total volume of gas. Thus, if depletion affects both p-NH₃ and N₂H⁺, it can explain why Johnstone et al. (2010) conclude that the p-NH₃ and N₂H⁺ excitation temperatures are independent of the central H₂ density.

Tobin et al. (2011) present interferometric observations of NH₃ and N₂H⁺ in a sample of class 0 and class I protostars. In particular, they observe that the position of the maximum intensity in these two tracers is offset from the position of the continuum peak, hence from the position of the protostar. Apart from depletion, an alternative way to explain the disappearance of these molecules from the gas phase can be chemical destruction. Above 20 K, CO will be released from the ice mantles and will react with N₂H⁺ to form HCO⁺. This reaction will thus efficiently destroy N₂H⁺. Such a process has been invoked by Tobin et al. (2013) to explain the central hole observed for both the N₂H⁺ and N₂D⁺ abundances in the protostellar object L1157. Additionally, in this object, N₂H⁺ is found to be more centrally peaked than N₂D⁺, which might be a consequence of variations in the H₂ ortho-to-para ratio. In the case of NH₃, an efficient destruction can come from the reaction with HCO⁺ (Woon & Herbst 2009). Thus, chemical destruction could explain the shape of the abundance profiles we obtain for NH₃ and N₂H⁺ in the B1b region. Two arguments can be invoked against such an possibility. First, the observations of H¹³CO⁺ performed by Hirano et al. (1999) have shown that at the positions of B1-bS and B1-bN, the abundance of this molecule was reduced by a factor ~3−5 with respect to the H¹³CO⁺ abundance of the surrounding envelope. Second, the dust temperature obtained from the SED fitting is found to be below 20 K. The hypothesis of depletion therefore seems to be favoured in order to explain the central decrease in the N₂H⁺ and NH₃ abundances.

That NH₃ is strongly depleted onto dust grains in the innermost region of B1b is strengthened by the analysis of the NH₃ ice features. At the position of the Spitzer source that is offset by ~15″ from the central position of our model, Bottinelli et al. (2010) infers a total column density of NH₃ locked in the ice of ~7 × 10¹⁷ cm^-2, while we derive a column density N(p-NH₃) ~ 5.5 × 10¹⁴ cm^-2 in the gas phase, at the central position of our observations. Converting this column density to the total amount of NH₃ requires knowledge of the respective abundances of ortho- and para-NH₃. It is expected that if NH₃ is formed in the gas phase, the ratio will be close to unity. On the other hand, a formation that occurred on dust grains would raise the ortho-to-para ratio above unity (Umemoto et al. 1999). If we assume a similar amount of ortho- and para-NH₃, this leads to a total gas phase column density of N(NH₃) ~ 1 × 10¹⁵ cm^-2. Since the observations of Bottinelli et al. (2010) are offset with respect to our core centre, we can infer that the gas phase NH₃ abundance is a factor >700 lower than inferred from the ice feature observations.

In principle, comparison of the respective abundances in the gas and solid phases should shed light on the evolutionary stage of the source. Indeed, it is predicted theoretically that the relative abundance of NH₃ in both phases will vary with time during the collapse (e.g. Rodgers & Charnley 2003; Lee et al. 2004; Aikawa et al. 2008, 2012). In particular, after the formation of a protostar, the NH₃ molecules locked in the ices will be released into the gas phase, because of the warm-up, and the ratio of the NH₃ gas-to-dust abundances will increase with time. However, in the present case, this point can only be discussed qualitatively, because of the simplicity of the geometry we assumed to describe the B1b region. We consider a 1D-spherical model while two objects, B1-bS and B1-bN are present and separated by a projected distance of only ~20″. Moreover, if we consider, say, Fig. 2 of Aikawa et al. (2012), we see that ~500 years after the formation of a protostar, the increase in the NH₃ abundance in the gas phase occurs on the 10 AU scale (i.e. ~0.04″ at 235 pc) and on the 100 AU scale at t = 10⁴ years. This makes the FHSC stage (t = −560 yr in Aikawa et al. 2012 where t = 0 corresponds to the formation of the second hydrostatic core, i.e. the protostar) hardly distinguishable from the stage that harbours a newly born protostar, at least in single-dish observations. Thus, the abundance profile we derived for the gas-phase p-NH₃ abundance would be consistent with the predictions of Aikawa et al. (2012) from t = −560 yr to t = 10⁴ yr, and it is thus impossible to infer the evolutionary stage of the two sources in the B1b region from these observations, i.e. to disentangle the FHSC stage from the class 0 stage.

4.2. Deuteration

The abundance profiles reported in Fig. 21 show that at the innermost radii, the behaviour of the p-NH₃ and o-NH₂D abundances changes, with a ratio χ(o-NH₂D) / χ(p-NH₃) below unity for densities lower than ~6 × 10⁴ cm^-3, and above unity for higher densities. We found that this enhancement is higher than an order of magnitude at the core centre. This behaviour will still be true if we consider the total column densities of NH₃ and NH₂D and if we consider standard ratios to relate the abundances of the ortho and para species. Figure 3 of Aikawa et al. (2012) shows that a ratio χ(NH₂D)/χ(NH₃) > 1 is expected at small radii, but only at an evolutionary stage close to the FHSC stage. At times greater than 500 yrs after the formation of the first hydrostatic core, the NH₃ abundance always exceeds the NH₂D abundance throughout the envelope. This is due to the fact that when the gas temperature becomes high enough, the backward reaction that leads to the deuterated isotopologues of H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ becomes efficient enough to inhibit the deuteration of this ion. Finally, the NH₂D map in Fig. 19 shows that the NH₂D peak is close to the position of the B1-bN core. This would suggest that this core is closer to the FHSC stage than B1-bS.

4.3. Nitrogen fractionation

To date, there have been few studies estimating the abundance ratio of ¹⁴N and ¹⁵N isotopologues in dark clouds. In the source L1544, the ¹⁴N/¹⁵N ratio was estimated to be 1000 ± 200 by Bizzocchi et al. (2010, 2013) for the N₂H⁺ isotopologues. Toward the same source, Hily-Blant et al. (2013) estimated 140 <¹⁴N/¹⁵N < 360 for the HCN molecule, using a double-isotope analysis of H¹³CN and HC¹⁵N. In the same study, a similar analysis was also performed towards L183 leading to 140 <¹⁴N/¹⁵N < 250. Ikeda et al. (2002) observed HCN and CO isotopologues in a sample of three dark clouds in the Taurus molecular cloud. They derived values of 151 ± 16 and >813 for L1521E and L1498, respectively, using the ¹²C/¹³C ratio derived from the analysis of the CO isotopologues. In the case of L1498, this ratio was determined to be ¹²C/¹³C > 117 from the analysis of ¹³C¹⁸O and ¹²C¹⁸O lines, a value significantly higher than the standard value in the local ISM. Assuming a ratio of 60 rather than 117 would lead to HCN / HC¹⁵N = 472 ± 55. Moreover, both Ikeda et al. (2002) and Hily-Blant et al. (2013) find that the H¹³CN/HC¹⁵N abundance ratio varies as a function position. As pointed out by Ikeda et al. (2002), this result can be interpreted either as due to an incomplete mixing of the parental gas or as variations in the chemistry of nitrogen fractionation within the clouds.

The NH₃ and NH₂D data considered in the present study have already been analysed by Lis et al. (2010) and Gérin et al. (2009), using an LTE approach. In the case of NH₃, the isotopologue abundance ratio was found to be 334 ± 50, in good agreement with the present determination of $\mathrm{30}{\mathrm{0}}_{-40}^{\mathrm{+}\mathrm{55}}$. However, in the case of NH₂D, we obtain a value of $\mathrm{23}{\mathrm{0}}_{-55}^{\mathrm{+}\mathrm{105}}$ compared to the earlier estimate of $\mathrm{47}{\mathrm{0}}_{-100}^{\mathrm{+}\mathrm{170}}$. As explained in Sect. 3.6, this occurs because the assumption of equal excitation temperatures for NH₂D and ¹⁵NH₂D made in Gérin et al. (2009) is not totally valid. Moreover, we found that variations as small as ~20% in the excitation temperatures lead to a factor 2 variations in the column density ratio. Apart from B1b, Lis et al. (2010) and Gérin et al. (2009) report column density ratio estimates in other dark clouds. A value of 344 ± 173 was obtained for NH₃ in NGC1333, and for NH₂D, values in the range 350 < ¹⁴N/¹⁵N < 850 were obtained in a sample of four clouds. In the latter case, as for B1b, the column density ratio may have been overestimated by a factor ~2.

Tennekes et al. (2006) report observations of various HCN and HNC isotopologues towards Cha-MMS1. They analysed the emission from these molecules both with a non-LTE non-local radiative transfer modelling and with the LTE approach. In their non-LTE analysis, they fixed the isotopologue abundances by assuming ¹²C/¹³C = 20 and ¹⁴N/¹⁵N = 280, on the basis of their LTE analysis. At that time, it has not yet been recognized that the HCN and HNC collisional rate coefficients could differ by a significant factor. Indeed, in their LTE analysis, they assumed a single excitation temperature for all the HCN and HNC isotopologues. However, due to differing HCN and HNC rate coefficients and because of line trapping effects that affect the main isotopologues, this assumption is not valid, as can be seen in Table 2.

An estimate of the ¹⁴N/¹⁵N galactic gradient was reported by Dahmen et al. (1995) from HCN observations and Adande & Ziurys (2012) from CN and HNC observations. The first study used a double isotope method with a ¹²C/¹³C ratio based on H₂CO observations. However, H₂CO might undergo carbon fractionation effects, and the results from Dahmen et al. (1995) were then reanalysed by Adande & Ziurys (2012) assuming a different ¹²C/¹³C ratio as derived from the CN observations of Milam et al. (2005). From both HCN, HNC, and CN, they obtain a value for the local ISM of ¹⁴N/¹⁵N = 290 ± 40. The sample of clouds used in these studies have typical kinetic temperatures in the range 25 K < T_K < 90 K (Adande & Ziurys 2012). At such high kinetic temperatures, fractionation effects for the nitrogen are unlikely to occur, and the quoted value should thus be representative of the elemental atomic abundance ratio. This estimate is in good agreement with the estimate of ¹⁴N/¹⁵N = 237$\begin{array}{c}\mathrm{+}\mathrm{27}\\ -21\end{array}$ obtained toward diffuse clouds of the local ISM observed in absorption toward compact extragalactic sources (Lucas & Liszt 1998).

The values obtained in the current study for the column density ratios are summarized in Fig. 22. In this figure, we also report the mean ratio (~296) computed from all the values except that of ¹⁵NNH⁺, for which we only obtained a lower limit. In the figure, for this molecule, we indicate this lower limit, as well as approximate error bars that are obtained by scaling that of N¹⁵NH⁺. The values reported for the nitriles come from the double isotope method rather than from a direct estimate (see Sect. 3.4). By examining Fig. 22, it can be seen that given the error bars, the column density ratios are consistent with a constant value. The only exception is N₂H⁺, and particularly ¹⁵NNH⁺, for which a higher ¹⁴N/¹⁵N ratio is obtained. The current result supports the conclusion that no fractionation affects the nitrogen chemistry in the B1b region. The mean value we obtain for B1b is in good agreement with the estimate of 290 ± 40 (Adande & Ziurys 2012) obtained toward warmer sources in the local ISM. The latter value could be representative of the gas associated with the Perseus molecular cloud. But, it was also recognized (Ikeda et al. 2002; Adande & Ziurys 2012) that the elemental ¹⁴N/¹⁵N atomic ratio can vary from source to source. The absence of large differences between the nitriles and the nitrogen hydrides obtained in the current study suggest that no fractionation effects occur in B1b, and a ratio ¹⁴N/¹⁵N ~ 300 should thus be representative of this region.

However, two things must be kept in mind while considering this conclusion. The first point is that the N₂H⁺ column density ratio show a tendency towards higher values, in particular for ¹⁵NNH⁺, for which we derived a lower limit N(N₂H⁺)/N(¹⁵NNH⁺) > 600. The second point concerns the nitriles. For this family, we derived the main isotopologue column densities from the ¹³C isotopologues assuming a ratio ¹²C/¹³C = 60. This value is representative of the local ISM, but as for the nitrogen case, can be affected by source-to-source variations.

The astrochemical modelling by Wirström et al. (2012) suggests that the nitrogen hydride and nitrile families should show different degrees of fractionation. Such a finding is supported by the results of Ikeda et al. (2002) and Hily-Blant et al. (2013) where low ¹⁴N/¹⁵N ratios were derived for HCN. Such differential fractionation effects are appealing since they would enable reconciliation of some of the observations made in the solar system. However, in the current study, no differential fractionation is found between the nitrogen hydride and nitrile families. Regarding the absence of fractionation in B1b, a possible explanation may be that this region contains two cores in a late evolutionary stage. The gas temperature in the core centre is estimated to be ~17 K, somewhat higher than the lowest temperatures ~7 K that can be reached during the prestellar phase. As for our definition of the mean excitation temperature given in Sect. 3.2, we can define a mean gas temperature, averaged over all the molecules considered in this work, which describes the mean temperature of the gas in which the emission originates. We find that ${\overline{\mathit{T}}}_{\mathrm{K}}\mathrm{=}\mathrm{13.}{\mathrm{7}}_{-2.6}^{\mathrm{+}\mathrm{1.9}}$ K. At such high temperatures, the nitrogen fractionation might already be ineffective. Indeed, as shown by Rodgers & Charnley (2008), decreasing the gas temperature from 10 K to 7 K can lead to an increase by a factor ~3 in the ¹⁵N-enrichment for the NH₃ and N₂ molecules. This effect is partly due to the low exothermicity of the exchange reactions that involve the ¹⁵N⁺ ions, which are typically in the range ΔE = 20 − 30 K. As a consequence, small variations in the temperature around the canonical value of 10 K, commonly used to describe prestellar cores, could result in substantial variations in ¹⁵N-fractionation effects, because of the variations in the backward reaction rates. However, to date, a quantitative estimate of the amount of ¹⁵N-fractionation for temperatures higher than 10 K has never been addressed by astrochemical models. The strong dependence of ¹⁵N-fractionation effects with temperature would imply that the ¹⁴N/¹⁵N ratio should vary within the sources with temperature gradients. However, this issue cannot be addressed in the present study, since to test this hypothesis, it would be necessary to constrain the models on the basis of emission maps of the ¹⁵N-substituted species. Finally, the current model shows that at the core centre, all the molecules are affected by depletion onto dust grains. Thus, it cannot be discarded that fractionation effects could have occurred earlier in the life of the B1b cloud. However, the ¹⁵N-enrichment would nowadays be unobservable, since the molecules formed in the previous colder evolutionary stages would subsequently have been incorporated into the ice mantles that surround the dust grains.

As discussed by Wirström et al. (2012), in order to link the prestellar phase chemistry with the observations of the solar system, the astrochemical models have to be able to produce high degrees of fractionation for both ¹⁵N- and D-substituted isotopologues. Another prerequisite is that they also have to account for the fact that the meteoritic hot-spots for these two isotopes are not necessarily coexistent. In the present case, while no ¹⁵N-fractionation is seen for both the nitriles and nitrogen hydrides in B1b, we conclude that the nitrogen hydrides are highly fractioned in deuterium. The D-fractionation for the nitriles is also found to be high, even if lower than for the nitrogen hydride family. The gas depleted at this evolutionary stage will create hot spots in D and will not be associated with high ¹⁵N enrichments.

5. Conclusion

We have presented a model based on continuum emission and molecular line observations for the B1b region, located in the Perseus molecular cloud. First, we modelled the 350 μm and 1.2 mm continuum radial profiles to obtain an estimate of the temperature and H₂ density distributions in the source. These density and temperature profiles were subsequently used as input for the non-local radiative transfer modelling of the molecular lines. The molecular data include D-, ¹³C-, and ¹⁵N-substituted isotopologues of the two molecules that belong to the nitrogen hydride family, i.e. N₂H⁺ and NH₃, as well as those from the nitrile family, HCN, HNC, and CN. For each molecule, the radial abundance profile was divided into several zones, in which the molecular abundance was considered as a free parameter. The best estimate for the abundance profile was then obtained from a Levenberg-Marquardt algorithm, which minimizes the χ² between models and observations. Additionally, since our ¹³C- and ¹⁵N-observations are only towards a single position, we assumed that the various isotopologues have abundance profiles that only differ by a scaling factor.

A common feature obtained for all the molecules is that in the innermost part of B1b, all the molecules have abundances, which are at least one order of magnitude lower than the peak abundance found in the envelope. Since all the molecules seem affected, this feature is interpreted as depletion onto dust grains rather than due to chemical destruction, a process that might be invoked for N₂H⁺ since the innermost region of B1b has a temperature that is close to the CO thermal desorption temperature. We derived high degrees of deuteration with the D/H ratio increasing inwards. On the other hand, we derived similar ¹⁴N/¹⁵N abundance ratios for all the molecules and did so independently of the chemical family, with an average ratio of ~300. This value is consistent with a recent estimate of the ¹⁴N/¹⁵N atomic elemental abundance ratio obtained for warmer sources in local ISM. We thus conclude that the ratio ¹⁴N/¹⁵N~300 inferred from the molecules is representative of the atomic elemental abundances of the parental molecular cloud. We propose that the absence of ¹⁵N-fractionation in B1b is a consequence of the relatively high gas temperature of this region, as compared to the low temperatures that can be reached during the prestellar phase. The current state of the chemistry in B1b would not produce any ¹⁵N-fractionation. On the other hand, it cannot be excluded by the current observations that conditions favourable to producing a ¹⁵N-enrichment were met during earlier evolutionary stages of the cloud. Since the molecules are found to be affected by depletion onto dust grains, the products of such a chemistry could have been incorporated into the ice mantles that surround the dust grains.

Given that ¹⁵N-fractionation effects can have strong implications for our understanding of the diversity observed in the various bodies of the solar system, it would be important to extend such a study focussing on earlier evolutionary stages with lower gas temperature, which is more favourable for producing nitrogen fractionation. Additionally, it was found in the current study that a direct estimate of the fractionation in the nitriles, i.e. HCN and HNC, is difficult because of the high opacities of the lines. This difficulty is not found for CN since this molecule is an open-shell molecule, which leads to a redistribution of the total opacity of the rotational line over a larger number of radiative transitions. Observations of less abundant species from this family, like HC₃N, would thus be particularly helpful. Moreover, we had to use the double isotope method to estimate the HCN and HNC column densities from the H¹³CN and HN¹³C isotopologues. We assumed the standard ¹²C/¹³C elemental ratio, which introduces uncertainties due to possible source-to-source variations. It would thus be helpful to complement the set of molecular observations with other carbon-bearing molecules, such as C¹⁸O.

Finally, we have shown that the isotopic ratios are sensitive to the assumptions made in deriving the column densities. Even relatively small differences in excitation temperatures of different isotopologues may lead to large variations in the derived isotopic ratios. The use of a sophisticated radiative transfer method is therefore recommended.

---

Acknowledgments

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). This work is based upon observations carried out at the Caltech Submillimeter Observatory, which is operated by the California Institute of Technology under cooperative agreement with the National Science Foundation (AST-0838261). This paper was partially supported within the programme CONSOLIDER INGENIO 2010, under grant Molecular Astrophysics: The Herschel and ALMA Era.- ASTROMOL (Ref.: CSD2009- 00038). We also thank the Spanish MICINN for funding support through grants AYA2006-14876 and AYA2009-07304.

References

Appendix A: Inelastic rate coefficients for the CN isotopologues

We used the two dimensional CN-He PES of Lique et al. (2010) to determine hyperfine resolved excitation and de-excitation cross sections and rate coefficients of C¹⁵N and ¹³CN molecules by He. CN-He rate coefficients obtained from this PES were found to be in excellent agreement with the experimental rotational (with unresolved fine and hyperfine structure) rate coefficients of Fei et al. (1994). This agreement demonstrate the accuracy of the PES and suggests that the rate coefficients obtained from this potential will be accurate enough to allow for improved astrophysical modelling of the ISM.

The Alexander’s description (Alexander 1982) of the inelastic scattering between an atom and a diatomic molecule in a ${}^{\mathrm{2}}\mathrm{\Sigma }^{\mathrm{+}}$ electronic state and a fully-quantum close-coupling calculations was used in order to obtain the S^J(jl;j′l′) diffusion matrix between fine structure levels of CN. J and l denote the total angular momentum (J = j + l) and the orbital angular momentum quantum numbers, respectively. The nuclear spin-free S^J(jl;j′l′) matrix were computed recently in Lique & Kłos (2011) and we used these matrices in the following.

For C¹⁵N, The coupling between the nuclear spin (I = 1/2) of the ¹⁵N atom and the molecular rotation results in a weak splitting (Alexander & Dagdigian 1985) of each rotational level j. Each hyperfine level is designated by a quantum number F (F = I + j) varying between | I − j | and I + j. The integral cross sections corresponding to transitions between hyperfine levels of the C¹⁵N molecule can be obtained from scattering S-matrix between fine structure levels using the re-coupling method of Alexander & Dagdigian (1985). Inelastic cross sections associated with a transition from the initial hyperfine level (NjF) to a final hyperfine level (N′j′F′) were thus obtained as follow: $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{NjF}\mathrm{\to }{\mathit{N}}^{\mathrm{\prime }}{\mathit{j}}^{\mathrm{\prime }}{\mathit{F}}^{\mathrm{\prime }}}& \mathrm{=}& \frac{\mathit{\pi }}{{\mathit{k}}_{\mathit{NjF}}^{\mathrm{2}}}\mathrm{(}\mathrm{2}{\mathit{F}}^{\mathrm{\prime }}\mathrm{+}\mathrm{1}\mathrm{)}\\ & & \mathrm{\times }\sum _{\mathit{K}}{\left\{\begin{array}{c}\\ & & \\ & & \end{array}\right\}}^{\mathrm{2}}{\mathit{P}}^{\mathit{K}}\mathrm{(}\mathit{j}\mathrm{\to }{\mathit{j}}^{\mathrm{\prime }}\mathrm{)}\mathit{.}\end{array}$(A.1)The P^K(j → j′) are the tensor opacities defined by: ${\mathit{P}}^{\mathit{K}}\mathrm{(}\mathit{j}\mathrm{\to }{\mathit{j}}^{\mathrm{\prime }}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{K}\mathrm{+}\mathrm{1}}\sum _{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left|}{\mathit{T}}^{\mathit{K}}\mathrm{\right(}\mathit{jl}\mathrm{;}{\mathit{j}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}\mathit{.}$(A.2)The reduced T-matrix elements (where T = 1 − S) are defined by Alexander & Davis (1983): $\begin{array}{ccc}{\mathit{T}}^{\mathit{K}}\mathrm{\left(}\mathit{jl}\mathrm{;}{\mathit{j}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{-}\mathit{j}\mathrm{-}{\mathit{l}}^{\mathrm{\prime }}}\mathrm{(}\mathrm{2}\mathit{K}\mathrm{+}\mathrm{1}\mathrm{\left)}\sum _{\mathit{J}}\mathrm{\right(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{J}}\mathrm{(}\mathrm{2}\mathit{J}\mathrm{+}\mathrm{1}\mathrm{)}\\ & & \mathrm{\times }\left\{\begin{array}{c}\\ & & \\ & & \end{array}\right\}{\mathit{T}}^{\mathit{J}}\mathrm{\left(}\mathit{jl}\mathrm{;}{\mathit{j}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{.}\end{array}$(A.3)For ¹³CN, the situation is more complex since both the ¹³C and N atoms have non zero nuclear spin. The coupling between the nuclear spin (I₁ = 1/2) of the ¹³C atom and the molecular rotation is first taken into account. Each hyperfine level is designated by a quantum number F₁ (F − 1 = I₁ + j) varying between | I₁ − j | and I₁ + j. The coupling between the nuclear spin (I₂ = 1) of the N atom and the first hyperfine state is then taken into account. Each hyperfine level is designated by a quantum number F (F = I₂ + F₁) varying between | I₂ − F₁ | and I₂ + F₁. The integral cross sections corresponding to transitions between hyperfine levels of the ¹³CN molecule can be obtained from scattering S-matrix between fine structure levels also using the re-coupling method of Alexander & Dagdigian (1985). Inelastic cross sections associated with a transition from the initial hyperfine level (NjF₁F) to a final hyperfine level (${\mathit{N}}^{\mathrm{\prime }}{\mathit{j}}^{\mathrm{\prime }}{\mathit{F}}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{F}}^{\mathrm{\prime }}$) were thus obtained as follow (Daniel et al. 2004): $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{Nj}{\mathit{F}}_{\mathrm{1}}\mathit{F}\mathrm{\to }{\mathit{N}}^{\mathrm{\prime }}{\mathit{j}}^{\mathrm{\prime }}{\mathit{F}}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{F}}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathit{\pi }}{{\mathit{k}}_{\mathit{NjF}}^{\mathrm{2}}}\mathrm{(}\mathrm{2}{\mathit{F}}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{F}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{F}}^{\mathrm{\prime }}\mathrm{+}\mathrm{1}\mathrm{)}& & \\ \mathrm{\times }\sum _{\mathit{K}}{\left\{\begin{array}{c}\\ & & \\ & & \end{array}\right\}}^{\mathrm{2}}{\left\{\begin{array}{c}\\ & & \\ & & \end{array}\right\}}^{\mathrm{2}}{\mathit{P}}^{\mathit{K}}\mathrm{(}\mathit{j}\mathrm{\to }{\mathit{j}}^{\mathrm{\prime }}\mathrm{)}\mathit{.}& & \end{array}$(A.4)From the rotationally inelastic cross sections σ_α → β(E_c), one can obtain the corresponding thermal rate coefficients at temperature T by an average over the collision energy (E_c) (Smith 1980): $\begin{array}{ccc}{\mathit{k}}_{\mathit{\alpha }\mathrm{\to }\mathit{\beta }}\mathrm{\left(}\mathit{T}\mathrm{\right)}& \mathrm{=}& {\left(\frac{\mathrm{8}}{\mathit{\pi \mu }{\mathit{k}}_{\mathrm{B}}^{\mathrm{3}}{\mathit{T}}^{\mathrm{3}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}\\ & & \mathrm{\times }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{\sigma }}_{\mathit{\alpha }\mathrm{\to }\mathit{\beta }}\hspace{0.17em}{\mathit{E}}_{\mathrm{c}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\frac{{\mathit{E}}_{\mathrm{c}}}{{\mathit{k}}_{\mathrm{B}}\mathit{T}}}\hspace{0.17em}\mathrm{d}{\mathit{E}}_{\mathrm{c}}\end{array}$(A.5)where k_B is Boltzmann’s constant and μ is the reduced mass of the CN-He complex. α and β designate the initial and final states of the molecule, respectively.

The complete set of (de)excitation rate coefficients for rotational transitions considered in this studies is available on-line from the BASECOL website⁷.

All Tables

All Figures