© P. Hily-Blant et al. 2022

1 Introduction

The partitioning of interstellar elements between the solid, gas, and ice phases is the result of the complex interplay of refractory grain formation, gas-phase reactive collisions, and adsorption and desorption processes, which all depend in a non-linear fashion on the environmental conditions (density, temperature, radiation field, etc.). Beyond the important successes of astrochemistry, the current understanding of these processes is still incomplete and large uncertainties on the abundances of important elements remain (Le Gal et al. 2014).

In the interstellar medium (ISM), some metals are consumed to build up the solid, refractory material of interstellar grains. The propensity of a given chemical element toward condensation into solids has been constrained primarily through absorption-line studies in the ultraviolet (Jenkins 2009), showing that it depends on the element. Atoms not locked into grain cores remain in the gas phase and constitute the volatile phase of the ISM (gas and ices). All gas-phase species are expected to freeze-out at sufficiently high density (Walmsley et al. 2004; Friesen et al. 2014), not with the same degree of efficiency, however (Tafalla et al. 2004; Hily-Blant et al. 2008).

The total abundance (solid, ice, and gas) of any element is equal to its cosmic abundance, or cosmic abundance standard (CAS), also called elemental abundance. The values for CAS are obtained through laboratory measurements in chondrites (Lodders 2003) and photospheric spectroscopy of the Sun (Asplund et al. 2009) or of young, B-type stars located within a ~ 1 kpc from the Sun (Nieva & Przybilla 2012). It may be emphasized that photospheric determinations are a difficult task. The compositionof the refractory phase is constrained indirectly by the wavelength-dependent reddening of background starlight (e.g., Draine & Li 2007). The composition of ices, which become an important reservoir at gas densities n_H typically above a few 10 cm⁻³ (Whittet 2010, Fig. 3), is mostly derived from infrared absorption spectra and has been so far limited to bright, star-forming regions (Boogert et al. 2015). At the lower densities of the diffuse ISM where ices are not present and chemical elements are essentially atomic, elemental gas-phase abundances can be measured directly through visible and UV spectroscopy (Sofia & Meyer 2001). In contrast, elemental abundances in high-density regions, such as dense cores, star-forming regions, or protoplanetary disks, have never been obtained directly because the main carriers in the gas and ice of the major metals are usually not observable. Thus, the partitioning between dust, gas-phase, and ice can only be inferred by combining observations of trace species and chemical models. The resulting uncertainties are usually larger than the statistical ones, although it is usually not possible to properly quantify by what factor. Part of the uncertainties of the models predictions come from the volatile abundance, that is CAS minus refractory, of the main elements: carbon, oxygen, nitrogen, and sulfur. These uncertainties then propagate in a complex way throughout the entire chemical network. Obtaining direct constraints on the volatile abundance of these elements in the dense ISM is thus of paramount importance to obtain quantitative and accurate predictions from astrochemical models.

The volatile abundances (gas plus ice) of carbon and oxygen are uncertain by typically a factor of two (Whittet 2010; Jenkins 2014), leading to large variations of the abundances predicted by astrochemical models (Bergin et al. 1997; Pratap et al. 1997). These uncertainties can be traced back to measurements of the refractory or gas-phase reservoirs. In particular, it was shown that the C/O ratio, namely, the relative amounts of gas-phase carbon and oxygen, can have a strong impact on the abundances of gas-phase N and N₂, presumably the main carriers of nitrogen in dense clouds (Hily-Blant et al. 2010b; Le Gal et al. 2014).

In contrast with carbon and oxygen, sulfur shows no tendency toward depletion into the solid phase (Jenkins 2009). Thus, it is expected that sulfur is entirely volatile that is, sulfur-bearing species are either in the gas or in the ice. The cosmic abundance of sulfur is well determined from spectroscopy of solar and B-type stars, with measurements converging toward parts per million (ppm) (Daflon et al. 2009). In this paper, we thus adopt 14 ppm as a reference value for the volatile abundance of sulfur. This value is in harmony with the abundance of gas-phase atomic sulfur in the Orion Molecular Cloud (Goicoechea & Cuadrado 2021).

In the dense ISM, only trace sulfur-bearing species are observed, which sum up to ~ 1% of the cosmic abundance, and chemical models are unable to consistently reproduce the observed abundances (Tieftrunk et al. 1994; Woods et al. 2015). Solutions to the problem generally involve species in ices, such as organo-sulfur molecules or sulfur chains or rings (Jiménez-Escobar & Muñoz Caro 2011; Jiménez-Escobar et al. 2014; Martín-Doménech et al. 2016; Laas & Caselli 2019). However, in ices, it is only OCS that has been detected so far, toward massive YSOs, with an abundance of at most 0.06 ppm, while only upper limits have been obtained for H₂S, which is not consistent with it being the main carrier of sulfur (Smith 1991; Boogert et al. 2015). Specifically, gas-grain models predict that S-bearing species in ices represent a minor reservoir (Woods et al. 2015). Recent chemical models using an updated gas-phase and gas-grain chemical network suggest that sulfur is, indeed, only weakly depleted (by a factor of three) from the gas-phase in dense clouds, but leaving open the question of the main sulfur reservoirs (Vidal et al. 2017). This is also consistent with the predictions from PDR models, which best reproduce the observations toward the Horsehead nebula when assuming that sulfur is not depleted and in the form of gas-phase S atoms (Goicoechea et al. 2006). Further supporting evidence comes from the finding that gas-phase atomic sulfur is an important, if not dominant, reservoir of sulfur in the dense ISM, derived from the detection of the fine-structure line emission of SI and leading to a lower limit on the abundance of S atoms of 5–10% with respect to the cosmic abundance (Anderson et al. 2013). However, constraints on the abundance of atomic sulfur in pre-stellar cores are scarce and highly indirect.

In this paper, we present observations of the NS molecule toward a sample of low-mass, starless, and pre-stellar cores that are combined with published N₂H⁺ column densities. Based on astrochemical model calculations over a wide parameter range, we show that the NS:N₂H⁺ abundance ratio is a direct tracer of the abundance of atomic sulfur in the gas-phase, a result that can be understood through a simple analysis of the chemistry of NS.

Fig. 1 Large-scale environment of the four cores studied in this work. The color scale is the dust emission at 500 μm measured with the Herschel/SPIRE instrument (André et al. 2010). At the distance of the Taurus molecular cloud (140 pc), 1° corresponds to 2.4 pc.

2 Observations

2.1 Source sample

Our source sample is made of low-mass dense cores located in the Taurus molecular cloud (see Fig. 1) and includes both starless cores (which do not show in-fall signatures) and one pre-stellar core (L1544), and none of the observed cores harbor a protostar. The reference position and some basic properties of each core are summarized in Table 1.

L1521E is a young, starless core located in the B218 filament in the Taurus molecular cloud (Tafalla & Santiago 2004) that is characterized by a very low depletion level, exemplified by its high CO:N₂H⁺ column density ratio. L1521B is another starless core located in the B216 filament in the Taurus molecular cloud and is classified as a chemically young core on the basis of its large abundance of carbon-chain molecules (e.g., CCS) over late-type species suchas ammonia (Hirota et al. 2004). Furthermore, the B216 region is believed to be the most recently condensed filament in the L1495 filamentary complex (Seo et al. 2015). Next is L1498, a well-studied starless core situated in the periphery of the Taurus molecular cloud. L1498 is believed to be in a more advanced stage of evolution than the young L1521E and L1521B, on the basis of clear CO-depletion patterns and high abundance of N₂ -daughter molecules (Tafalla et al. 2004; Magalhães et al. 2018). Finally, L1544, the most evolved of the cores in our sample, is a pre-stellar object also located at the outskirts of the Taurus molecular cloud. Contrary to the other cores of our sample, it shows signposts of infall in the very inner regions, as traced by inverse P-Cygni profile of water lines (Caselli et al. 2012). It is also the core with the highest central density.

In view of computing abundances, we have estimated the total molecular hydrogen column density, , in 26′′ -beam or 35′′-beam, from published radial density profiles or from multi-wavelength dust emission maps analyzed in a standard fashion assuming a black-body radiation combined with a frequency-dependent dust opacity. The minimization was performed using a Monte-Carlo Markov-chains (MCMC) approach, as described in detail in Appendix A.

2.2 IRAM/30m observations

The observed millimeter transitions of the NS radical (mass number 46) involve the rotational levels of the ²Π_1∕2 manifold of the fundamental electronic and vibration state and have a structure similar to the NO radical (Gerin et al. 1992). The lowest rotational level of the ²Π_3∕2 manifold lies 322 K above the J = 1∕2 lowest level of the ²Π_1∕2 manifold, and is not significantly populated at the low temperatures of dark clouds. The unpaired electron leads to Λ-type doubling of each rotation level leading to two rotational levels of opposite parity. Each Λ-level further splits into a set of hyperfine levels due to the coupling of the nuclear spin of the nitrogen atom (I = 1) with the orbital angular momentum J. The hyperfine levels are described by the quantum number F where F = I + J. Radiative transitions are subject to selection rules: ΔJ = ±1, ΔF = 0, ±1, and a change of parity. This results into a set of hyperfine transitions as shown in Fig. C.1.

The observations have been performed in August 2016 (project 008-16) with the IRAM-30m telescope. The EMIR receivers were connected to the VESPA backends facility to measure simultaneously the + → − and − → + hyperfine transitions between the J = 5∕2 and J = 3∕2 rotational levels (see Table C.1 and Fig. C.1) with 20 kHz spectral resolution, or 0.05 km s⁻¹ at 115 GHz. Weather conditions were good and opacity at these frequencies is dominated by atmospheric oxygen. The median receiver temperature, system temperature, and zenith opacity were respectively 60 K, 294 K, and 0.4. In some cases, two J = 7∕2 → 5∕2 transitions could be observed. Assuming no beam dilution and negligible contribution from the error-beams, the antenna temperature was converted into main-beam temperatureas with F_eff and B_eff the forward and beam efficiencies. Values of F_eff and B_eff are 0.94 and 0.78 at 115 GHz,respectively, and 0.93 and 0.71 at 161 GHz. The half-power beam widths at the frequencies relevant to the present study—including the N₂H⁺ lines from the literature – are 26′′, 21′′, and 15′′, at 93, 115, and161 GHz, respectively.

3 Line properties

The hyperfine (hf) structure of the two sets of lines corresponding to each parity of the ² Π_1∕2 manifold of NS is detected in all sources, and toward most positions (see Figs. 2 and 3). The properties of each hyperfine component obtained through Gaussian fitting are summarized in Tables B.1 to B.5. At a few positions toward L1521E and L1498 (see Fig. B.3), transitions within the J = 7∕2–5/2 manifold at 161 GHz were observed and detected. Toward L1521B, the observed positions, based on the reference position of Hirota et al. (2004), are shifted by more than 80′′ from the dust emission peak.

The lineintensity varies among the sources, although by small factors, while the variation within a given source are up to a factor two. There is no apparent systematic dependence of the intensity with the dust emission, but in L1498, the line intensity increases when moving away from the central region.

The line properties have been obtained from Gaussian fits, the results of which are summarized in Tables B.1–B.5. The line width of the strongest hyperfine transition within J = 5∕2–3/2 toward L1521E and L1521B are typically 0.2–0.3 km s⁻¹. Toward L1498 and L1544, the lines are narrower, with a full width ≈ 0.15–0.25 km s⁻¹, which would correspond to a gas temperature above 20 K, well above the values in these sources (see references in Table 1). Thus, the 115 GHz lines are emitted in regions where non-thermal (turbulent) broadening is actually dominant.

It is instructive to compare the relative line intensities, normalized to that of the brightest component, to the theoretical line ratios expected for thermalized level populations in the optically thin case (see Fig. C.2). It is found that the intensity of the weakest line is compatible, in many cases, with its theoretical relative strength of 0.375. The same is true for the second brightest hyperfine line (with a relative intensity 0.63) although some departures are observed. Overall, this suggests that the NS linesare not optically thick and that the level populations are close to be thermalized at a single excitation temperature.

Fig. 2 Maps of the dust emission and positions of the NS spectra for each source in the sample. Left: dust emission maps at λ = 1.3 mm (IRAM/30 m) toward L1521E, L1498, and L1544 (Tafalla et al. 2002, 2004), and at λ = 0.5 mm (Herschel/SPIRE) toward L1521B. Contour levels (in MJy sr−1) and HPBW are indicated in each panel. The reference positions (0′′, 0 offsets) are given in Table 1. Right: map of the NS emission line (only the main hyperfine component of the Π+ manifold at 115 GHz is shown; see also Fig. 3). The intensity scale (, −0.1 to 0.5 K) is the same for all sources.

4 Methodology

4.1 Column density of N₂H⁺

When the total opacity of a N₂H⁺ transition is high, measuring the column density must take into account deviations from the single excitation temperature assumption because of hyperfine trapping effects (Tafalla et al. 2004; Daniel et al. 2007). To do so requires both a source structure (density, temperature, velocity, line broadening) and an abundance profile (e.g., Magalhães et al. 2018).

For L1498, we have used the value of Daniel et al. (2007), N(N₂H⁺) = 7.3(3) × 10¹² cm⁻², obtained by taking into account above effects. Indeed, this column density is a factor 1.5 larger than the one derived by Tafalla et al. (2004) who neglected hyperfine trapping. The use of updated, and more accurate, collision rate coefficients also explain the different column density. The N₂H⁺ column density obtained by Daniel et al. (2007) corresponds to offsets (+ 10′′, 0′′) relative to the dust peak, or (+ 20′′, +20′′) in our reference frame. Nevertheless, from the constant N₂H⁺ abundance derived by Tafalla et al. (2004) and the H₂ column density profile derived from the density profile from Magalhães et al. (2018), we infer the variation of N(N₂H⁺) over 20” to be 10%.

In L1521E, we derive the 26′′-beam averaged column density of N₂H⁺ toward offsets (− 10′′, −10′′) using the constant abundance, N₂H ⁺:H₂ = 2 × 10⁻¹¹ (or [N₂H ⁺]= 10⁻¹¹) from Tafalla & Santiago (2004), in combination with our H₂ column density, (see Appendix A), and thus N(N₂H⁺) = 4.8 × 10¹¹ cm⁻². The abundance in L1521E being 8.5 times lower than in L1498 using the same method as in Tafalla et al. (2004), the lines are probably not optically thick and we thus did not apply the correction factor of 1.5 found by Daniel et al. (2007) in L1498. Our estimate of N(N₂H⁺) is ≈ 50 times lower than the peak column density estimated by Nagy et al. (2019) from their local thermal equilibrium (LTE) study. However, for the sake of consistency with the H₂ column density and because the analysis from Tafalla et al. uses a 1D non-LTE radiative transfer, we use their results. For L1521B, we use the 1σ upper limit obtained by Hirota et al. (2004) of 1.9 × 10¹² cm⁻².

Finally, toward L1544, the 26′′-beam-averaged column density is N(N₂H⁺) = 4.1 × 10¹³ cm⁻² (Bizzocchi et al. 2013). In this source, we also estimated the N₂H⁺ column density at the offsets where NS was observed, by assuming a constant N₂H ⁺:H₂ ratio of 6.5 × 10⁻¹⁰, as in Bizzocchi et al. (2013), combined with their source model (from their Fig. 2), which corresponds very closely to a H₂ density profile of the form , with n₀ = 2 × 10⁷ cm⁻³ and r₀ = 316 au. The H₂ column density corresponding to each offset can be computed analytically, from which the column density N₂H⁺ follows.

4.2 Column density of NS

Two complementary, methods were applied to derive the total column density of NS. One is the usual HFS method in the CLASS software (e.g., Hily-Blant et al. 2013), which assumes a single excitation temperature within a hyperfine manifold. The second uses the escape probability approximation to relax the single excitation temperature assumption. However, as shall be discussed below, degenerate sets of solutions are obtained in some cases and the results from the HFS method are biased toward the high optical depth solutions. Therefore, only the results from the non-LTE approach are discussed and used in what follows.

Non-equilibrium level populations were computed with the public code RADEX (van der Tak et al. 2007) using new, dedicated NS-H₂ collision rate coefficients (see Appendix C.1) and spectroscopic properties from the Cologne Database for Molecular Spectroscopy (CDMS, Müller et al. 2005). In this procedure, the FWHM Δv of the lines (obtained from Gaussian fitting and taken to be equal for all hf lines) and the background temperature are held fixed, along with the H₂ density, the NS column density, N¹, and the kinetic temperature are the free parameters.

The three-dimensional parameter space exploration may be performed with a forward-fitting approach, through a regular grid at each point of which the reduced χ² is determined. Alternative methods include MCMC constrained by a prior probability on the meaningful parameter space. In this approach, we obtain the probability of the parameters’ given measured values. Both approaches lead to complementary views, but the MCMC is more optimized in that it naturally refines the parameter space exploration close to the various local minimum locations, while the forward-fitting requires progressively refined regular grids. In this study, we used the publicly-available implementation of MCMC of Foreman-Mackey et al. (2013).

It is widely known that line inversion is often degenerate in cold cores especially when only one rotational transition is available. In the present study, most sources were observed only through the J = 5∕2–3∕2 rotational manifold and a degenerate set of solutions, from high-/low-N to low-/high-N was effectively found for several lines of sight (see Fig. C.4). In contrast, when detected, the J =7∕2–5/2 usually breaks this degeneracy. Furthermore, the density is known to be larger than typically 10⁴ cm⁻³ in the sources of our sample, which also helps in alleviating the degeneracy. Eventually, the NS column density could be constrained in each source.

Fig. 3 NS(J = 5∕2–3/2) spectra, in the temperature scale (in K), toward L1521E, L1521B, L1498, and L1544 (from left to right). In each panel, the hf set of lines (identified by their quantum number F) is split into the Π+ (top) and Π− (bottom) (see Table C.1). For the sake of clarity, hf transitions are shown individually. The position number (see Fig. 2) and the offsets (in arcsec), with respect to the reference position given in Table 1, are indicated on the left of each spectrum. The vertical dotted lines indicate the adopted systemic velocity. The horizontal and vertical scales are indicated in bold face.

5 Results

The MCMC minimization procedure is described in Appendix C and the results are shown in Figs. C.3–C.6. The resulting column densities and ratios are given in Table 2.

5.1 L1521E

The physical properties and column density of NS were constrained toward four positions for which the signal-to-noise ratio (S/N) was sufficient. The results of the minimization process are shown in Fig. C.3. At offsets (− 10′′, −10′′), the J =7∕2–5∕2 line was observed enabling the H₂ density and NS column density to be determined. The most likely density is dex² (or ) to which corresponds³ N(NS) = 12.8(2) dex (or ).

Toward the other offsets, only the J =5∕2–3∕2 manifold was observed and an ensemble of density and column density are found that can reproduce the line properties, while the kinetic temperature is unconstrained. Nevertheless, the density at these positions is most likely lower than at the peak position at (− 10′′, −10′′). At offsets (− 30′′, −10′′), although a continuumof solutions cannot be excluded, the priors distributions for the density and column density favor two solutions, at low (≈ 10⁴ cm⁻³) and high (> few 10⁵ cm⁻³) density. Thus, focusing on the low-density one gives a column density of 13.3(4) dex. A conservative approach consists in adopting 4.75 dex as a lower limit on the H₂ density, leading to lower limits on the NS column density of 12.8, 12.9, and 12.6 dex, toward offsets (− 30′′, −10′′), (+ 10′′, −10′′), and (− 10′′, −30′′) respectively. Hence, a trend may be apparent of an outward increase of the NS column density.

5.2 L1521B

In the case of L1521B, where only J = 5∕2–3∕2 was observed, a continuous set of density and column density is found at all positions (see Fig. C.4). Nevertheless, at offsets (60′′, 0′′), (40′′, 0′′), and (50′′, −20′′), which are also the closest to the dust emission peak (see Fig. 2), solutions with H₂ density above 10⁵ cm⁻³ are favored, in agreement with the previous study of Tafalla et al. (2004) who found . In this density regime, the column density N(NS), is well constrained. The results are summarized, at each position, in Table 2, showing that the column density varies little among the various positions, with N(NS) = 12.3(2) dex.

Toward the other positions, high H₂ density (above ~10⁵ cm⁻³) solutions are also more likely with associated values of N(NS) again very close to 12.3 dex. However, it is not possible to firmly exclude lower densities (), for which N(NS) ≈13.2 dex.

5.3 L1498

In contrast to L1521B, the H₂ density and NS column density in L1498 are well constrained toward several positions in L1498, in particular, at the four locations where the J =7∕2–5∕2 Π⁻ and Π⁺ manifolds are well detected (see Fig. B.3). Toward the dust peak position at offsets (+10′′, +20′′), we obtain N(NS) = 12.7(2) dex, with a weak dependence on the temperature. We find significant spatial variations of the NS column density, from 13.3(5) at (+20′′, −10′′) and (+40′′, −20′′) to 12.0(1) at (−20′′, +10′′) and (−40′′, +20′′). However, nosystematic radial variation is apparent. At other locations, the solutions are degenerate.

5.4 L1544

L1544 is the densest of the four cores considered in this study, with a central density of (Bizzocchi et al. 2013). The NS lines from L1544 are weaker than in the other three sources and although NS was detected at the five observed positions, the derivation of the column density was performed toward only three, (0′′, 0′′), (20′′, 20′′), and (20′′, −20′′). Toward (0′′, 0′′), a two-component Gaussian fit was performed and the minimization was done using the strongest of the two. Toward the other two positions, the line intensity was obtained with a single-component Gaussian fit. The results, shown in Fig. C.6, confirm that the density is higher than 10⁶ cm⁻³, and most likely closer to 3 × 10⁷ cm⁻³. The corresponding column density are low, N(NS) ≈ 12.0(1) dex.

5.5 Summary

The main observational result is that the NS:N₂H⁺ column density ratio varies by more than two orders of magnitude, from ≈ 15 in L1521E to 0.02 in L1544. This variation is due to a combination of a decrease of the NS abundance and an increase of that of N₂H⁺. Although typical uncertainties of 50–100% may apply on the H₂ column density (see Fig. C.3), not only these are much smaller than the source-to-source variations reported here, but the consistent method used in determining the H₂ column density mitigates the impact of such uncertainties when comparing the NS:N₂H⁺ ratio in our sample. There are trends (at 1σ level) of spatial variations of the NS:N₂H⁺ ratio within L1544, and of the abundance of NS in L1521B and L1498. The abundance of NS vary between ≈ 10⁻¹¹ and 1.5 × 10⁻¹⁰. These values are typically a factor of 2 or more lower than to those obtained in TMC-1 and L134N (McGonagle et al. 1994). This may be due to the LTE assumption which, as already noted, is biased toward larger column densities.

6 Abundance of gas-phase sulfur

6.1 Chemical models

We now turn to chemical model calculations to obtain constraints on the abundance of gas-phase sulfur, S_gas – namely, the total abundance of gas-phase sulfur-bearing species – from the measured NS:N₂H⁺ column density ratio. We consider here only steady-state calculations, allowing us to explore a parameter space in which the C/O ratio, the density, the cosmic-ray ionization of molecular hydrogen , and the volatile sulfur abundance S_gas, are varied. The abundance of gas-phase oxygen O_gas is varied while that of carbon, C_gas, is kept constant and equal to 83 ppm, in such a manner that the C/O ratio (equal to C_gas/O_gas) covers the range from 0.4 to 1.2 by steps of 0.2. Three representative cosmic-ray ionization rates are considered to account for the lack of accurate determinations in the dense gas: , 3, and 10. Finally, the value of S_gas covers two orders of magnitude, namely 0.14, 0.443, 1.4, 4.43, and 14 ppm. In what follows (except otherwise stated), the gas:grain processes (depletion, hydrogen saturation, and desorption) were not included, although varying the value of O_gas may be seen as a simple mean of simulating such processes (Le Gal et al. 2014). The UGAN network (Hily-Blant et al. 2018) has been updated and upgraded, using the KIDA database (Wakelam et al. 2015) for those reactions involved in the chemistry of NS. The updated and new reactions are listed in Table D.1.

In the present calculations, the pre-stellar core is modeled as an isothermal gaseous sphere with purely atomic initial conditions (see Table 3) and the abundances are computed until a steady-state is reached. As detailed above, a grid of models has three free parameters (C/O, ζ₁₇, S_gas), while the gas temperature and density are held fixed with T_kin = 10 K, and n₄ = n_H∕10⁴ cm⁻³ = 1, 10, and 100.

6.2 Results

Figure 4 shows the NS:N₂H⁺ abundance ratio as a function of S_gas for two gridsof models, with n₄ = 1 (top) and n₄ = 10 (bottom). A close-to linear power-law is evident, which is found to depend only weakly on C/O and ζ₁₇, while the intercept does depend on these environmental parameters. The origin of the correlation can be understood based on simple chemical considerations (see Sect. 7).

The primary result of this study is that the tight NS:N₂H⁺–S_gas correlation shows that the sources in our sample have distinct gas-phase sulfur abundances: L1521E has the highest value of S_gas, L1544 the lowest, L1498 is intermediate, and the lower limit in L1521B indicates a high value of S_gas.

To explore the results in more details, let us first consider the model with n₄ = 1 and ζ₁₇ = 3 (top and middle panels). For each value of S_gas, the predicted NS:N₂H⁺ ratio depends only weakly on C/O and ζ₁₇, and only marginally modifies the trend of a decreasing S_gas from L1521E to L1544. Rough estimates of S_gas are summarized in Table 4.

Toward L1521E, our observations and models indicate that volatile sulfur abundance is compatible with the cosmic abundance of 14 ppm, provided that C/O = 0.6–0.8, while C/O = 0.4 would require an unrealistic value of S_gas (>30 ppm). Toward L1521B, there is a lower limit on S_gas of 1.4 ppm for C/O = 0.6–0.8, and S_gas > 3 ppm for lower or higher C/O ratios. In L1498, the observed NS:N₂H⁺ ratio constrains S_gas to 1.4–3 ppm over the entire range of C/O. L1544 is the most sulfur-depleted source in this model, with S_gas < 0.14 ppm toward the innermost position over the entire n₄ = 1 grid, two orders of magnitude lower than in L1521E.

While models with n₄ = 1 are representative of the outer parts of the cores, examining models with a higher density (n₄ = 10 and 100) enable a study of regions closer to the inner plateau. Compared to models with n₄ = 1, the predicted NS:N₂H⁺ ratios for n₄ = 10 are more scattered across the C/O range. Nevertheless, the above trends remain valid. Thus, in L1521E, C/O = 0.4 requires a very large S_gas and is thus excluded by these models, while C/O = 0.6 to 1.2 are all consistent with S_gas = 14 ppm. In L1521B, the lower limit on the NS:N₂H⁺ ratio gives S_gas > 8 ppm for C/O = 0.4, while for C/O = 0.6–0.8, a lower limit of 1.4 ppm on S_gas is found. In L1498, a range of S_gas = 0.8–6 ppm is derived for C/O = 0.4–1.2, while S_gas = 1–2 ppm for C/O between 0.6 and 0.8. In L1544, low values of S_gas, as in modelswith n₄ = 1, are derived withS_gas < 0.6 ppm. L1544 is thus the most sulfur-depleted core, followed by L1498, while L1521E and L1521B are both compatible with a non-depleted sulfur gas phase.

The spatial distribution of NS and N₂H⁺ in L1544 is also of interest. As can be seen from Table 2, the abundance of NS and the NS:N₂H⁺ ratio increase outward by a factor of two over 30′′ or 4200 au. We note that at these offsets, where the density in L1544 is = 10⁵ cm⁻³, similar to that in L1498, S_gas becomes ≈ 0.14 ppm for C/O = 0.6 (and n₄ = 10), still significantly lower than in L1498 (2 ppm). The spatial variation of the NS:N₂H⁺ abundance ratio in L1544 provides clear indication that depletion of sulfur increases with density. In L1521B, the abundance of NS also decreases when moving toward the peak position of the core, but the lack of N₂H⁺ information does not allow us to constrain the NS:N₂H⁺ ratio and thereby the spatial variation of S_gas in this source.

Fig. 4 [NS]:[N2H+] ratio as a function of Sgas at steady-state from chemical models with with n4 = 1 (top) and 10 (bottom) and ζ17 = 1, 3, and 10 (left to right). In each panel, five values of the C/O ratio (0.4–1.2) are shown. The observed NS:N2H+ column density ratios (with ± 1σ uncertainties) toward L1521E, L1498, and L1544 (central position), are shown as hashed areas, and the lower limit in L1521B is indicated with an arrow. The black, straight, line shows our analytical prediction given in Eq. (13).

7 Chemical considerations

7.1 Chemistry of NS

The reason why the NS:N₂H⁺ abundance ratio provides a direct constraint on S_gas may be understood by analyzing the chemistry of NS under typical dark cloud physical conditions. The main chemical routes involved in the formation and destruction of NS, along with its chemical link with atomic sulfur, are summarized in Fig. 5. Nitrogen sulfide is formed primarily by the neutral-neutral reactions: (1) (2)

both with temperature-independent rates of 10⁻¹⁰ cm³ s⁻¹. The dissociative recombination of HCNS⁺ also contributes to the formation ofNS, but to a lesser extent. In parallel, NS is destroyed via reactions with atomic oxygen and carbon: (3) (4)

We note that rates for reactions (1)–(4) have been estimated and evaluated by the KIDA committee, but measurements at a low temperature are required.

The relative weight of reactions (1) and (2) depends primarily on S_gas: reaction (1) is always the main formation route (more than 50%) while the importance of reaction (2) increases with S_gas, becoming typically half that of reaction (1) for S_gas=14 ppm.

As already emphasized by McGonagle et al. (1994) and Cernicharo et al. (2018), the abundance of NS is sensitive to the carbon to oxygen gas-phase abundance ratio. Indeed, the relative importance of reactions (3) and (4) depend on the relative abundance of atomic oxygen and carbon in the gas phase: the oxygen route will dominate if n(O):n(C) is larger than the ratio of the two rates, > 2 × 10⁻¹⁰∕3 × 10⁻¹¹ = 6.7. Reaction (3) is found to dominate destruction (>55%) when C/O ≤ 0.6; at higher C/O ratios, reaction (4) becomes the main destruction route (≈ 35%) but other routes, involving collisions with H⁺ or atomic nitrogen, N, also become non-negligible.

Fig. 5 Main chemical routes involved in the formation and destruction of NS under dark cloud conditions and its chemical link with atomic sulfur. We note that SECPHO stands for secondary photons.

7.2 Steady-state abundance at low C/O

From the previous remarks, it transpires that a single analytical expression predicting the abundance of NS over the entire parameter space explored in this study is highly unlikely. Nevertheless, focusing on the low C/O-regime (≤ 0.6), a simple expression can be obtained. Values of C/O ratios lower than 1 were found in conditions typical of pre-stellar cores but C/O ratio larger than 1 cannot be ruled out (Pratap et al. 1997; Bergin et al. 1997; Le Gal et al. 2014).

From reactions (1)–(4), the steady-state abundance of NS is given by (5)

where we note [X] = n(X)∕n_H the abundance relative to hydrogen nuclei. The abundance of atomic sulfur may then be written as: (6)

Over most of the parameter space, reaction (1) is actually the primary formation route – although its importance reduces to ≈60% when S_gas=14 ppm – such that the right-hand side in Eq. (6) can reasonably be reduced to its first term. Regarding the destruction, ours limiting to C/O ≤ 0.6 ensures that reaction (3) is the main route such that Eq. (6) can be simplified into: (7)

The steady-state abundance of NH may be obtained considering that the formation is primarily through the dissociative recombination of N₂H⁺ (Hily-Blant et al. 2010a; Vigren et al. 2012) and its destruction proceeds by the neutral-neutral reaction with atomic oxygen NH + O → NO + H, such that: (8)

where is the dissociative recombination rate toward NH. It thus follows that in the low-C/O regime: (9)

At 10 K, k_DR = 2.6 × 10⁻⁷ cm³ s⁻¹, and the (dimensionless) numerical coefficient is 7.6 × 10⁻⁵. Since atomic sulfur is, in all our models, the main sulfur carrier, this expression shows how the NS:N₂H⁺ abundance ratio can be used to measure S_gas. We now discuss the term.

7.3 Electronic abundance

The electronic density, n_e, results from the competition between ionization and recombination processes and is thus driven by the cosmic-ray ionization rate and the total density. The abundance of electrons also depends on the total grain area per unit volume although to a lesser extent (Walmsley et al. 2004; Flower et al. 2007). The formation rate of electrons is due to the ionization of H₂ and He by cosmic rays, with rates and 0.5, respectively, while the destruction proceeds primarily through dissociative recombination with molecular ions with rate k_R. This leads to (assuming n_H = 2n(H₂) and [He]=0.1): (10)

where we have noted [m⁺] the fractional abundance of positively charged molecular ions. Provided that molecular ions dominate over atomic ions, charge neutrality translates into [m⁺] = [e⁻], and the usual dependence is recovered (e.g., Walmsley et al. 2004). Numerically, with k_R ≈ 2 × 10⁻⁶ cm³ s⁻¹ at 10 K, we obtain (11)

with n₄ the proton density in units of 10⁴ cm⁻³ and ζ₁₇ the H₂ cosmic-ray ionization rate in units of 10⁻¹⁷ s⁻¹. These values compare well with the results of our numerical models (see Tables D.2 and D.3) as long as S_gas and the C/O ratio remain lower than a few ppm and ≈0.6, respectively.

However, the high-S_gas regime deserves special attention because S⁺ becomes the main charge carrier so that the approximation [e⁻] = [m⁺] is no longer valid. In such instances, we obtain the following: (12)

instead of Eq. (11). Moreover, the chemical state of the gas-phase – and more specifically its ionization state – may become increasingly sensitive to the amount of molecular oxygen and charged species (Pineau des Forêts et al. 1992; Le Bourlot et al. 1995; Boger & Sternberg 2006).

7.4 Abundance of atomic oxygen

The abundance of atomic oxygen in the gas phase depends on the partitioning of oxygen among its main carriers in solid, ice, and gas forms, and there are indications that this partitioning is not well understood (Whittet 2010; Jenkins 2014). In cores exposed to the ambient interstellar UV field, the abundance of oxygen in the gas-phase is essentially driven by the competition between adsorption and photo-desorption of water ice (Hollenbach et al. 2009). On observational grounds, the abundance of oxygen-bearing species in ices in low-mass cores are reviewed in Boogert et al. (2015). The abundance of water, CO, and CO₂ in ices are 39, 12, and 13 ppm, respectively, which, together with other trace species, account for ≈ 80 ppm. On the other hand, 240 ppm of oxygen is captured into refractory material (Jenkins 2009), so that the amount of oxygen locked in solid (ice and refractory) is ≈320 ppm (see also Hily-Blant et al. 2018).

The abundance of elemental oxygen in the local ISM is not precisely known, but it is most likely in the range 442 to 575 ppm (Sofia & Meyer 2001; Nieva & Przybilla 2012), from which the corresponding total gas-phase oxygen abundance is O_gas = 442 − 320 = 120 ppm and O_gas = 575 − 320 = 255 ppm, respectively.

In the gas-phase of dense cores, the main carriers are CO and atomic oxygen, with molecular oxygen being a minor species (Goldsmith et al. 2011). In low-C/O environments, approximately 83 ppm of the oxygen (according to our assumed abundance of gas-phasecarbon) will be locked into CO; namely, what remains in the gas phase, presumably in the form of atomic oxygen, ranges from 120 − 83 = 37 to 255 − 83 = 172 ppm, or [O ] = 3.7 × 10⁻⁵ to 1.7 × 10⁻⁴. As noticed by Whittet (2010), and highlighted by Jenkins (2014), the higher limit would imply that a large fraction of oxygen is locked into an unidentified form. In taking a conservative approach, we use the lower limit in what follows, namely, [O ] = 3.7 × 10⁻⁵. It is worth stressing that this reasoning applies only if C/O < 1. In addition, we expect that the abundance of gas-phase oxygen will decrease as the C/O ratio increases.

7.5 The abundance of atomic sulfur for C/O < 1

Inserting Eq. (11) into Eq. (9) together with [O ] = 37 ppm, and assuming a kinetic temperature of 10 K, the abundance of atomic sulfur in the gas phase is expressed as: (13)

In all our models, atomic sulfur is the main carrier of sulfur (see Tables D.2 and D.3) such that one may adopt [S] = S_gas. Therefore, the previous expression may be compared to the results of our model calculations (black line in Fig. 4).

Despite the various assumptions applied to derive Eq. (13), the overall agreement between the analytical and numerical predictions may appear surprisingly good. This is especially the case at low C/O, such as 0.4, for both n₄ = 1 and 10, and regardless of ζ₁₇. The scatter of the model predictions is less than 1 dex over most of the parameter space; however, the models with n₄ = 10 and ζ₁₇ = 1 has the largest dispersion of predicted ratios, and the gradient of S_gas among the sources in the sample is less clear unless a consistent value of C/O is used in the four sources. Incidentally, since the NS:N₂H⁺ ratio is directly related to [S], the high value of this ratio in L1521E is a direct proof that atomic sulfur is the main carrier of sulfur.

8 Discussion

The observable NS:N₂H⁺ abundance ratio has been shown to provide a direct constraint to the total gas-phase abundance of sulfur, S_gas, and, hence of the abundance of gas-phase atomic sulfur [S]. Furthermore, the case of L1521E provides a clear confirmation that atomic sulfur is the main carrier of sulfur, which is in agreement with our and published models (Boger & Sternberg 2006).

The differences of S_gas in our source sample may be interpreted in terms of evolutionary stage. In early-type–chemistry cores, such as L1521E and L1521B, our observations indicate that sulfur is essentially not depleted, with 14 ppm in the gas-phase in the form of atomic sulfur. As cores evolve chemically, the amount of atomic sulfur decreases, which may be interpreted as gas-phase depletion in ices. Thus, in the chemically more evolved core L1498, the gas-phase abundance of sulfur is significantly lower than in L1521E (see Table 4), while L1544 has the highest level of sulfur depletion.

Our results are consistent with the non-depletion of sulfur in refractory material as derived from observations of ionized sulfur in the diffuse medium (Jenkins 2009, 2014). Moreover, these results are in agreement with the high abundance of atomic sulfur, [S ] = 3.5 ppm, needed to reproduce the observations of CS toward the Horsehead nebula (Goicoechea et al. 2006). We note that this abundance is very close to that in L1498 in our model with n₄ = 10. Indeed, the density in their photo-chemical models () is also comparable to that toward L1498. The present study is also in agreement with the finding by Nagy et al. (2019) of a higher abundance of sulfur species in L1521E compared to L1544. However, our method provides a direct quantitative estimate of S_gas. Our observations also provide observational support to recent chemical models showing that the abundances of S-bearing species toward the TMC-1(CP) dense cloud can be reasonably reproduced assuming a S_gas = 5–14 ppm (Vidal et al. 2017), in agreement with our measurements toward L1521B, and L1521E. In the translucent part of TMC-1, observations of CS, SO, and isotopologues, and HCS⁺, indicate that sulfur is significantly depleted, with S_gas = 0.8 ppm (Fuente et al. 2019), which is interpreted as a signature of strong sulfur depletion during the atomic-to-molecular transition: this behavior is thus different to what was observed in the Horsehead nebula, although the ratio of the radiation field to the H₂ density is similar ( in the Horsehead and 2–10 × 10⁻⁴ in TMC-1). It has also been proposed that sulfur depletion increases with density as S/H (Rodríguez-Baras et al. 2021), based on the behavior of sulfur-bearing species. The NS:N₂H⁺ method could be used to confirm this is true also for S_gas.

Spatial variations of S_gas among a given source could only be established in L1544. In the outer parts of L1544, a core for which signatures of collapse were found in the innermost regions from water line observations (Caselli et al. 2012), we find lower levels of sulfur depletion than in its central region. The increase in NO synthesis, located at 4200 au, where complex organic molecules (COMs) have been detected (Vastel et al. 2014; Jiménez-Serra et al. 2016, 2021), may indicate that the production of NS could be enhanced by the same mechanism as COMs (Vasyunin et al. 2017). Alternatively, this could be a demonstration of the progressive depletion of sulfur toward the inner regions. Indeed, in L1521B, there is a trend of decrease of the NS abundance toward the dust peak. However, since the current spectra were measured by taking the N₂H⁺ observations of Hirota et al. (2004) as a reference, it is necessary to obtain new NS and N₂H⁺ observations covering thedust peak.

The cores included in this study being located within the same large-scale environment, the differences in S_gas among the various sources are not due to variations in the elemental abundances or in the cosmic-ray ionization rate. Furthermore, the C/O ratio may also be comparable in the four cores since this ratio is dictated by the values of C_gas and O_gas, which are primarily induced by the composition of the refractory material. However, the fact that the C/O ratio is effectively comparable in all sources should be established from the observations.

From Table 4, a C/O ratio of 0.4, to which would correspond an unrealistic amount of sulfur in L1521E, may be excluded. Assuming that all the sources have a similar carbon and oxygen volatile budgets, our results suggest C/O ratios of 0.6–0.8. For these values, the sulfur reservoir decreases by more than two orders of magnitude in L1521E, L1498, and L1544, respectively, reflecting their evolutionary stage. The young L1521B core is also compatible with an non-depleted sulfur reservoir although the non-detection of N₂H⁺ prevents any conclusion. A value of C/O ≈ 0.8 was also derived, but in a different environment, toward the envelope of the low-mass protostar (Le Gal et al. 2014). Yet, in the TMC-1 ridge, a lower ratio (0.4–0.5) was proposed (Pratap et al. 1997), while Fuente et al. (2019) suggest higher values up to unity. The present results are unlikely to be used to measure the C/O ratio; a more promising method may be through the CN:NO abundance ratio (Le Gal et al. 2014). In any instance, independent observational constraints on the C/O ratio are needed in combination with the new NS:N₂H⁺ method, to measure S_gas.

One important assumption in our model is the abundance of atomic oxygen in Eq. (9), for which we adopted the lower limit of 37 ppm from the studies of Jenkins (2009) and Whittet (2010). Nevertheless, our grid of models, while exploring a range of C/O ratios by varying the total oxygen abundance O_gas, which also influences [O], show that our analytical derivation remains overall a good approximation. We have also tested the influence of oxygen depletion by performing calculations with gas-grain chemistry producing saturated species which are returned into the gas-phase by photo-desorption (see Hily-Blant et al. 2018). The impact of oxygen depletion was, however, not as direct as anticipated from our chemical analysis for C/O < 1 because carbon-destruction routes then become more efficient, leading to essentially unchanged NS:N₂H⁺ ratios.

In Eq. (7), the abundance of sulfur was expressed in terms of the NS:NH ratio, and the abundance of NH was derived assuming pure gas-phase formation. It has been proposed that a fraction of the intermediate radicals NH and NH₂ could be released in the gas-phase during the hydrogenation of nitrogen in ices (Hickson et al. 2015). However, laboratory experiments failed to detect these radicals indicating that hydrogenation proceeds faster than desorption (Hidaka et al. 2011). Furthermore, gas-phase processes alone can reproduce the observed amounts of nitrogen hydrides, as well as their ortho:para ratios, and their deuterated forms in low-mass star-forming environments (Le Gal et al. 2014; Hily-Blant et al. 2018).

Although this work is devoted to the NS:N₂H⁺ ratio, as ratios are less sensitive to uncertainties in the derivation of the H₂ column density, we compare in Fig. 6 the estimated abundances of NS and N₂H⁺ listed in Table 2 to those predicted by our models. Our model predictions cover the entire range of observed values of both NS and N₂H⁺, while previous attempts failed to reproduce the NS abundance (Millar & Herbst 1990; Hasegawa & Herbst 1993; Lee et al. 1996). This is, at least in part, due to our NS abundances being smaller by factors of two or more, compared to the earlier values in low-mass dense cores (McGonagle et al. 1994). This figure shows that lowering the density in our model both increases and broadens the range of the predicted N₂H⁺ abundances, while only slightly increasing the NS abundance. On the other hand, it is evident that models with the n₄ = 1, ζ₁₇ = 1 and with n₄ = 10, ζ₁₇ = 10 are similar, suggesting that ζ₁₇∕n₄ is a good control parameter. It is then apparent that the spread in predicted abundances increases with ζ₁₇ ∕n₄. Furthermore, from the previous discussion, C/O ratios between 0.6 and 0.8 may be favored, suggesting that a model with (n₄, ζ₁₇) = (1, 10) and C/O = 0.8 provides a reasonable agreement for three of the four sources, while over-predicting N₂H⁺ and NS in L1544. Nevertheless, from the models with ζ₁₇∕n₄ = 0.1 and 0.3, it appears that the NS:N₂H⁺ ratio L1544 can be explained with a C/O ratio within 0.6 and 0.8. More generally, it seems possible to find, for a given C/O similar in all sources, a set of models with different ζ₁₇ ∕n₄ ratios accounting for the NS:N₂H⁺ ratio in our sample. Going further in reproducing the abundances is, however, beyond the present study, as this would also require us to follow the dynamical collapse (e.g., Hily-Blant et al. 2018), and, more importantly, to simultaneously reproduce more nitrogen- and sulfur-bearing species.

New observations of NS and N₂H⁺ should be performed toward the dust peak of L1521B in order to obtain a measure of the ratio and confirm the absence of depletion in this young object. It would also be interesting to obtain maps of the N₂H⁺ column density toward these cores, similarly to what was done for L1544, in order to confirm the depletion of sulfur with density. Detailed one-dimensional radiative transfer calculations of NS, similar to what Daniel et al. (2007) did for N₂H⁺, may be useful in assessing the level of deviation from single excitation temperature within a rotational manifold at the hyperfine level.

From a theoretical perspective, chemical models are needed to explore the high sulfur-abundance regime in more detail, especially the electronic fraction. Also, the abundance of nitrogen-bearing species was found to decrease with increasing S_gas in Le Gal et al. (2014) and the NH₂:NH and NH₃:NH abundance ratios, measured in the envelope of the I16293-2422 protostar, put an upper limit of 8 ppm on S_gas. According to the present results, this would indicate that the envelope is in an advanced stage of chemical evolution. The dynamics of the collapse may likely be impacted by the present findings as we would expect a strong coupling of the gas with magnetic fields in the high-S_gas regime.

9 Conclusion

We provide a new and direct method to measure the total abundance of sulfur, that is, of atomic sulfur, in the gas-phase of dense starless cores. We also provide a simplified analysis of the chemistry of NS, leading to a simple relation between the NS:N₂H⁺ ratio and the abundance of gas-phase atomic sulfur, which is in good agreement with the numerical results. This is strong suggestion that the new method is robust over a broad parameter space. Our results provide a clear demonstration of the progressive depletion of atomic sulfur as cores evolve and, for a given core, with density. Our observationsalso show that sulfur is entirely in the gas-phase in atomic form in chemically young cores. Comparisons with spatially resolved observations of sulfur-bearing species in ices at densities below and above cm⁻³ with the JWST may reveal the progressive uptake of sulfur in ices that we have evidence of in L1544, and to elucidate the main carriers of sulfur in ices.

Fig. 6 Steady-state abundance of NS and N2H+ from our grids of models (top: n4 = 1; bottom: n4 = 10) are compared to the observational estimates listed in Table 2. The value of n4, ζ17, and the ratioζ17∕n4 are indicated in each panel. Five values of Sgas (see Table 3) are associated with each C/O ratio.

Acknowledgements

We are indebted to the anonymous referee for very careful reading and for questions and comments which led to a significant improvement of the paper. The following public softwares were used: GILDAS, Gnuplot, emcee. We also thank Mélisse Bonfand who contributed to this work during an internship at IPAG in 2015. PHB and FL acknowledge the Institut Universitaire de France (IUF) for financial support. This work was supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the Netherlands Organization for Scientific Research, and the National Research Council of Canada.

Appendix A H₂ Column density

We present, in greater detail, the derivation of the H₂ column density for L1521E, L1521B, and L1498.

A.1 L1521E

The H₂ column density was estimated following two approaches. One is based on the radial density profile from Tafalla & Santiago (2004), , which is integrated to give the H₂ column density at a distance b from the peak, . With n₀ = 2.7 × 10⁵ cm⁻³ and r₀ = 30′′, the 11"-beam averaged N_{H_2} toward the dust peak (at offsets − 15”, −5”) is thus 5.3 × 10²² cm⁻², which compares well with the value of 3.9 × 10²² cm⁻² that can be derived directly from the dust emission, 24 mJy/beam (Tafalla et al. 2002). Averaged in the 26"-beam of the N₂H⁺ spectra, the column density at offsets (−10”, −10”) is reduced by only 6% compared to its peak value of 5.1 × 10²² cm⁻² (17% for b = 20′′), or 4.8 × 10²² cm⁻².

As a sanity check, we also performed a second calculation based on dust emission maps at several wavelengths, thus providing a measurement of independent of the radialdensity profile. Thus the dust emission spectral energy density (SED) was fitted by a modified black-body: (A.1)

with B_ν(T) the black-body spectrum (in MJy sr⁻¹) at a dust temperature T, μ the mean molecular weight, equal to 2.8, m_H is the proton mass (taken equal to the atomic mass unit), and is the H₂ column density. The mass absorption coefficient is and we adopted κ₀ = 0.09 cm² g⁻¹ at λ = 250 μm (Magalhães et al. 2018, and references therein). The T, , and β parameter space was explored using the emcee MCMC sampler.

Observations at wavelengths from 160 to 1300 μm from IRAM, JCMT, and Herschel facilities have been used (Tafalla & Santiago 2004; Kirk et al. 2005; André et al. 2010; Marsh et al. 2016). At 70 μm, the source is not detected. The maps, convoluted by a 35"-HPBW Gaussian kernel and re-sampled on the 500 μm grid, are shown in the left panel of Fig. A.1. The resulting specific intensity are listed in Table A.1. Prior to minimization, a background was subtracted to each intensity, which was estimated in a nearby region outside the core (indicated by a hatched histogram in Fig. A.1.) At some wavelengths for which the maps are too small, the background was estimated from the histogram of the intensity.

The results of parameter space exploration for the dust peak position at offsets (−15”, −5”) are shown in the middle panel of Fig. A.1 and the corresponding SED in the right panel. We note that the peak position is not strictly the same at all wavelengths, with ≈ 5” differences, much smaller than the 35" beam at 500 μm. Results toward the offsets observed in NS are summarized in Table A.1. The derived dust temperature and spectral index are similar at all positions and typical of starless cores. The H₂ column density toward the peak position is 22.37(7)dex, or 2.4 ± 0.3 × 10²² cm⁻², a factor of two smaller than the value derived above from the density profile. Given the uncertainties induced by the mass absorption coefficient (with β = 2.2, our value at 1.3 mm is 0.0024 cm² g⁻¹, half that adopted by Tafalla et al 2004 in their study), and the fact that the density profile results from fitting an azimuth average of the map, an agreement within a factor of two can be considered a good one. Our results agree well with those obtained by Makiwa et al. (2016), although our study includes longer wavelength maps, thus providing tight constraints to the spectral index.

A.2 L1521B

For L1521B, Herschel maps at wavelengths 160, 250, 350, and 500 microns have been used. As for L1521E, the source is not detected at 70 μm and only a lower limit was obtained at 850 μm (Kirk et al. 2005). The maps, convolved by a 35"-HPBW Gaussian kernel and re-sampled on the 500 μm grid, are shown in Fig. A.2. The specific intensity are listed in Table A.1. At 160 μm, we note that our intensity compares well with that reported by Kirk et al. (2007) based on Spitzer data.

The resulting physical and dust parameters obtained from our likelihood calculations are shown in the middle panel of Fig. A.2. In addition to positions observed in this study, with the (0”, 0”) position corresponding to the peak of the carbon-chain molecule CCS emission Hirota et al. (2004), we also determined the H₂ column density toward the dust emission peak which is located at (115”, −15”) or RA, DEC = +04:24:23,+26:36:43 (J2000.0) (see Figs. 2 and A.2).

The most likely dust temperatures are ≈ 10 − 11K, typical of starless cores, and in good agreement with the results from ammonia line studies indicating that L1521B has a lower gas kinetic temperature than other starless cores in the L1495 region (Seo et al. 2015). The spectral index, β =1.8 − 2.0, is also typical of the expected value in starless cores. The H₂ column density is well constrained, 22.10(7)dex at (0”, 0”) offsets and 22.24(5) toward the dust peak at (115.-15). These values are in very good agreement with that derived by Seo et al. (2015) toward their cores 35 and 37 (see their Table 1).

A.3 L1498

Similarly to L1521E and L1521B, continuum maps made with IRAM/30m and Herschel/SPIRE were used to constrain the H₂ column density toward several offsets in L1498. In Fig. A.3, results for three offsets are shown, the dust peak, and two directions for which the fluxes differ by more than 10% compared to the peak. The resulting dust temperature is 11(1) K at all positions, in agreement with previous studies (Tafalla et al. 2004; Magalhães et al. 2018). The spectral index is ≈ 2. The total 35"-averaged H₂ column density are constrained with ≈ 30 − 40% accuracy (1σ), decreasing from 1.8(+8,-6) × 10²² cm⁻² toward the peak to 1.4(+8,-5) × 10²² cm⁻², toward (+80”, −40”). It is worth mentioning that background subtraction changes the minimized values only within 1σ.

We note that the 35"-averaged peak column density is thus 1.5 times smaller than the 26"-beam average value derived from the H₂ density profile of Magalhães et al. (2018), 2.8 × 10²² cm⁻², and a factor of two below the value from Tafalla et al. (2004) based on 1.3mm data only. Nevertheless, we consider this as a good agreement, and will adopt the here derived column density, taking into account that the profile obtained by Magalhaes et al is based on azimuth averaging.

A.4 L1544

Computing the H₂ column density toward L1544 could in principle be done from the source model used by Bizzocchi et al. (2013), itself adopted from Keto & Caselli (2010). This density profile has a similar form as for L1521E and L1498, with n₀ = 2 × 10⁷ cm⁻³, r₀ = 270 au (or 1.9", based on Fig. 1 of Keto et al, although no value is provided in these references), and α =2. The 26"-beam averaged H₂ column density is thus 8.9 × 10²² cm⁻², which is larger than the value of 6.6 × 10²² cm⁻² reported by Bizzocchi et al. However, a more observation-based determination of can be obtained as for the other three sources. In contrast with the other sources, the spectral index was kept fixed, β =2. If not, the most likely index is ≈ 1.6 and the dust temperatureis ≈ 11 K: since at long wavelength, I_ν ~ ν^2+β, this is likelya result of the 1.3mm intensity being overestimated (no background estimate). Instead, keeping β equal to 2 leads to dust temperatures ≈ 10 K. The 35"-beam–averaged column densities are summarized in Table A.1, with a peak value of 6.6 × 10²² cm⁻², in remarkable agreement with the value from Bizzocchi et al. Varying the spectral index by ± 0.2, the dust temperature varies by 0.5 K, between 10.1 and 9.1 for β =1.8 and 2.2 respectively, with corresponding column densities 22.68(1) and 22.95(1) dex. We thus adopt a peak column density .

Fig. A.1 Determination of in L1521E through SED-fitting by a modified black-body spectrum. Left Dust emission, measured in L1521E with the IRAM-30m/MAMBO, Herschel/PACS, and SPIRE facilities, used to compute the H2 column density (see Section A). The white cross-hairs locate the dust emission peak, at offsets (−15”, −5”). The contours levels (start, end, increment, in MJy sr−1) are indicated in each panel. The hatched area indicates the region used to compute the background level and rms. Middle Results toward offsets (−15”, −5”). Shown is a so-called "corner plot" showing the distribution of each free parameter (T, , and spectral index β) in the diagonal panels, and the correlation between any pair of them (off-diagonal panels). The most likely parameters based on the χ2 likelihood are shown by the central dashed, vertical line, in each diagonal panel (the other vertical lines show the 10% and 86% quantiles). Right: Observed intensities and fitted SED (see Table A.1).

Fig. A.2 Determination of in L1521B (see Fig. A.1). The ammonia cores (35 and 36) mapped by Seo et al 2015 are shown in the 250 μm map of the left panel.

Fig. A.3 Determination of in L1498 (see Fig. A.1).

Fig. A.4 Determination of in L1544 (see Fig. A.1). The spectral index was held fixed with β = 2.

Appendix B Analysis of NS lines

The hyperfine lines of each manifold were fitted independently, with no constraints on their relative velocity or line width (FWHM). In all cases, single-component Gaussian fits were used, except toward offsets (0”, 0”) and (−20”, 20”) in L1544 for which two-component Gaussian fits have been applied. The results, for each source, are shown in Fig. B.1–B.4 and the line properties are provided in Tables B.1–B.5.

Fig. B.1 Results of the Gaussian fitting to the J = 5∕2 − 3∕2 Π− (top left) and Π+ (top right) and J = 7∕2 − 5∕2 Π− (bottom left) and Π+ (bottom right) emission lines of NS toward L1521E ( scale). For each multi-fold, the three hyperfine lines (see Fig. C.1) are shown separately. In each panel, the vertical dashed line indicates the systemic velocity of 6.9 km s−1, and the Gaussian fit is shown in red.

Fig. B.2 Results of the Gaussian fitting for the J = 5∕2 − 3∕2 Π− (left) and Π+ (right) emission lines of NS toward L1521B (see Fig. B.1). The systemic velocity is 6.7 km s−1.

Fig. B.3 Results of the Gaussian fitting for the NS lines toward L1498. Top: J = 5∕2 − 3∕2 (left) and 7∕2 − 5∕2 (right) Π−. Bottom: same as top row for the 5∕2 − 3∕2 and 7∕2− 5∕2 Π+ manifolds. The systemic velocity is 7.9 km s−1.

Fig. B.4 Results of the Gaussian fitting for the J = 5∕2 − 3∕2 Π− (left) and Π+ (right) emission lines of NS in L1544 (see Fig. B.1). A 2-components Gaussian profile was adopted. The systemic velocity is 7.25 km s−1.

Appendix C NS column density

C.1 Collisional excitation of NS

The energy level diagram of the NS radical is shown in Fig. C.1 and the details on the hyperfine transitions of the J =5∕2 − 3∕2 and J =7∕2 − 5∕2 rotational transitions are summarized in Table C.1.

In the radiative transfer calculations used in this work, we used new collision rate coefficients for the NS molecule that will be presented ina forthcoming paper. Nevertheless, we provide here a brief description of the calculations of these new NS–He data. NS–He rate coefficients were computed from new rigid-rotor NS–He interaction potentials determined using the spin-unrestricted single and double excitation coupled cluster approach with non-iterative perturbation treatment of triples [UCCSD(T)] (Knowles et al. 1993, 2000). The NS bond distance was kept fixed at the equilibrium distance r_e = 2.82 bohr (Amano et al. 1969). The basis set for these electronic structure calculations was composed from augmented correlation-consistent polarized valence triple-zeta (aug-cc-pVTZ) basis set supplemented with additional mid-bond functions (Dunning 1989). Close-coupling calculations of the fine structure resolved rate coefficients (Alexander 1985a) were performed between NS energy levels up to J = 13∕2 in the ground spin orbit manifold. Hyperfine-structure–resolved rate coefficients were then derived using a re-coupling technique (Alexander 1985b). Finally, the NS–He rate coefficients were scaled by a factor of 1.4 (Roueff & Lique 2013) to obtain data for NS in collision with H₂.

C.2 Minimization

The results of the MCMC minimization using our NS–He collision rate coefficients are shown in Figs. C.3-C.6.

Fig. C.1 Energy level diagram (in K) for the first four rotational levels J of the 2 Π1∕2 spin sub-level of NS. The energy shift of the Λ-doubling states and the hyperfine levels with respect to the rotational level has been multiplied by 100 for the sake of visibility. The hyperfine quantum number F is indicated on the right of each hf level, where F = I + J with I the nuclear spin of the nitrogen atom (I = 1) and J the orbital angular momentum. Allowed transitions are restricted by the selection rules ΔJ = ±1, ΔF = 0, ±1, and + ↔−. The transitions shown in Figs. 3 and B.1–B.4, are outlined with arrows, where full lines indicate the 115 GHz lines, and the dashed lines those at 161 GHz.

Fig. C.2 Peak temperature (in mK) of the J = 5∕2 − 3∕2 Π− and Π+ manifolds (left and right columns resp.) toward L1521E, L1521B, L1498, and L1544 (from top to bottom) as obtained from Gaussian fits to the spectra. In each panel, the x-axis gives the position offsets (in ′′). Open symbols show the observations with 1σ error bars while the dashed lines and filled squares locate the theoretical peak intensity if optically thin lines and Boltzmann equilibrium assumptions hold.

Fig. C.3 Results of the RADEX-MCMC minimization of NS lines toward L1521E. Left: Results at offsets (−10”, −10”). Corner plot showing the distribution of Tkin, , and log10N within the prior intervals, [7:13], [4:7], and [11.5:14], respectively. Other panels: Cross-histograms toward position offsets (−30”, −10”), (+10”, −10”), and (−10”, −30”), from left to right. The inset shows the histograms of the column density for by step of 0.5. The red contours show 50% and 80% of the maximum.

Fig. C.4 Cross-histograms toward L1521B (see Fig. C.3).

Fig. C.5 Results of RADEX-MCMC minimization for L1498 (see also Fig. C.3). For offsets (0",0"): the first panel has three free parameters and the next two panels correspond to Tkin held fixed at 8 K and 10 K. Priors on Tkin, , and log10N taken in [7:13], [4:7], and [11:14], respectively. Toward (20",-10"), (10",20"), and (60",-30"), the kinetic temperature was fixed and equal to 10 K. For offsets (+10",+20"), the results are shown in a similar way as for L1521E.

Fig. C.6 Results of MCMC parameter space exploration (see Fig. C.3) for L1544 with priors on Tkin, , and log10N taken in [5:12], [5:8], and [11:15] respectively. The bottom row shows the results in a similar way as for L1521E.

Appendix D Chemical models

The list of the updated and new reactions included to the UGAN network are provided in Table D.1. The predicted steady-state abundances at density n_H = 10⁴, 10⁵, and 10⁶ cm⁻³ are listed in Tables D.2 to D.4, respectively.

References

All Tables

All Figures

	Fig. 1 Large-scale environment of the four cores studied in this work. The color scale is the dust emission at 500 μm measured with the Herschel/SPIRE instrument (André et al. 2010). At the distance of the Taurus molecular cloud (140 pc), 1° corresponds to 2.4 pc.
In the text

Fig. 2 Maps of the dust emission and positions of the NS spectra for each source in the sample. Left: dust emission maps at λ = 1.3 mm (IRAM/30 m) toward L1521E, L1498, and L1544 (Tafalla et al. 2002, 2004), and at λ = 0.5 mm (Herschel/SPIRE) toward L1521B. Contour levels (in MJy sr−1) and HPBW are indicated in each panel. The reference positions (0′′, 0 offsets) are given in Table 1. Right: map of the NS emission line (only the main hyperfine component of the Π+ manifold at 115 GHz is shown; see also Fig. 3). The intensity scale (, −0.1 to 0.5 K) is the same for all sources.

In the text

Fig. 3 NS(J = 5∕2–3/2) spectra, in the temperature scale (in K), toward L1521E, L1521B, L1498, and L1544 (from left to right). In each panel, the hf set of lines (identified by their quantum number F) is split into the Π+ (top) and Π− (bottom) (see Table C.1). For the sake of clarity, hf transitions are shown individually. The position number (see Fig. 2) and the offsets (in arcsec), with respect to the reference position given in Table 1, are indicated on the left of each spectrum. The vertical dotted lines indicate the adopted systemic velocity. The horizontal and vertical scales are indicated in bold face.

In the text

Fig. 4 [NS]:[N2H+] ratio as a function of Sgas at steady-state from chemical models with with n4 = 1 (top) and 10 (bottom) and ζ17 = 1, 3, and 10 (left to right). In each panel, five values of the C/O ratio (0.4–1.2) are shown. The observed NS:N2H+ column density ratios (with ± 1σ uncertainties) toward L1521E, L1498, and L1544 (central position), are shown as hashed areas, and the lower limit in L1521B is indicated with an arrow. The black, straight, line shows our analytical prediction given in Eq. (13).

In the text

	Fig. 5 Main chemical routes involved in the formation and destruction of NS under dark cloud conditions and its chemical link with atomic sulfur. We note that SECPHO stands for secondary photons.
In the text

	Fig. 6 Steady-state abundance of NS and N2H+ from our grids of models (top: n4 = 1; bottom: n4 = 10) are compared to the observational estimates listed in Table 2. The value of n4, ζ17, and the ratioζ17∕n4 are indicated in each panel. Five values of Sgas (see Table 3) are associated with each C/O ratio.
In the text

Fig. A.1 Determination of in L1521E through SED-fitting by a modified black-body spectrum. Left Dust emission, measured in L1521E with the IRAM-30m/MAMBO, Herschel/PACS, and SPIRE facilities, used to compute the H2 column density (see Section A). The white cross-hairs locate the dust emission peak, at offsets (−15”, −5”). The contours levels (start, end, increment, in MJy sr−1) are indicated in each panel. The hatched area indicates the region used to compute the background level and rms. Middle Results toward offsets (−15”, −5”). Shown is a so-called "corner plot" showing the distribution of each free parameter (T, , and spectral index β) in the diagonal panels, and the correlation between any pair of them (off-diagonal panels). The most likely parameters based on the χ2 likelihood are shown by the central dashed, vertical line, in each diagonal panel (the other vertical lines show the 10% and 86% quantiles). Right: Observed intensities and fitted SED (see Table A.1).

In the text

	Fig. A.2 Determination of in L1521B (see Fig. A.1). The ammonia cores (35 and 36) mapped by Seo et al 2015 are shown in the 250 μm map of the left panel.
In the text

	Fig. A.3 Determination of in L1498 (see Fig. A.1).
In the text

	Fig. A.4 Determination of in L1544 (see Fig. A.1). The spectral index was held fixed with β = 2.
In the text

Fig. B.1 Results of the Gaussian fitting to the J = 5∕2 − 3∕2 Π− (top left) and Π+ (top right) and J = 7∕2 − 5∕2 Π− (bottom left) and Π+ (bottom right) emission lines of NS toward L1521E ( scale). For each multi-fold, the three hyperfine lines (see Fig. C.1) are shown separately. In each panel, the vertical dashed line indicates the systemic velocity of 6.9 km s−1, and the Gaussian fit is shown in red.

In the text

	Fig. B.2 Results of the Gaussian fitting for the J = 5∕2 − 3∕2 Π− (left) and Π+ (right) emission lines of NS toward L1521B (see Fig. B.1). The systemic velocity is 6.7 km s−1.
In the text

	Fig. B.3 Results of the Gaussian fitting for the NS lines toward L1498. Top: J = 5∕2 − 3∕2 (left) and 7∕2 − 5∕2 (right) Π−. Bottom: same as top row for the 5∕2 − 3∕2 and 7∕2− 5∕2 Π+ manifolds. The systemic velocity is 7.9 km s−1.
In the text

	Fig. B.4 Results of the Gaussian fitting for the J = 5∕2 − 3∕2 Π− (left) and Π+ (right) emission lines of NS in L1544 (see Fig. B.1). A 2-components Gaussian profile was adopted. The systemic velocity is 7.25 km s−1.
In the text

Fig. C.1 Energy level diagram (in K) for the first four rotational levels J of the 2 Π1∕2 spin sub-level of NS. The energy shift of the Λ-doubling states and the hyperfine levels with respect to the rotational level has been multiplied by 100 for the sake of visibility. The hyperfine quantum number F is indicated on the right of each hf level, where F = I + J with I the nuclear spin of the nitrogen atom (I = 1) and J the orbital angular momentum. Allowed transitions are restricted by the selection rules ΔJ = ±1, ΔF = 0, ±1, and + ↔−. The transitions shown in Figs. 3 and B.1–B.4, are outlined with arrows, where full lines indicate the 115 GHz lines, and the dashed lines those at 161 GHz.

In the text

Fig. C.2 Peak temperature (in mK) of the J = 5∕2 − 3∕2 Π− and Π+ manifolds (left and right columns resp.) toward L1521E, L1521B, L1498, and L1544 (from top to bottom) as obtained from Gaussian fits to the spectra. In each panel, the x-axis gives the position offsets (in ′′). Open symbols show the observations with 1σ error bars while the dashed lines and filled squares locate the theoretical peak intensity if optically thin lines and Boltzmann equilibrium assumptions hold.

In the text

Fig. C.3 Results of the RADEX-MCMC minimization of NS lines toward L1521E. Left: Results at offsets (−10”, −10”). Corner plot showing the distribution of Tkin, , and log10N within the prior intervals, [7:13], [4:7], and [11.5:14], respectively. Other panels: Cross-histograms toward position offsets (−30”, −10”), (+10”, −10”), and (−10”, −30”), from left to right. The inset shows the histograms of the column density for by step of 0.5. The red contours show 50% and 80% of the maximum.

In the text

	Fig. C.4 Cross-histograms toward L1521B (see Fig. C.3).
In the text

Fig. C.5 Results of RADEX-MCMC minimization for L1498 (see also Fig. C.3). For offsets (0",0"): the first panel has three free parameters and the next two panels correspond to Tkin held fixed at 8 K and 10 K. Priors on Tkin, , and log10N taken in [7:13], [4:7], and [11:14], respectively. Toward (20",-10"), (10",20"), and (60",-30"), the kinetic temperature was fixed and equal to 10 K. For offsets (+10",+20"), the results are shown in a similar way as for L1521E.

In the text

	Fig. C.6 Results of MCMC parameter space exploration (see Fig. C.3) for L1544 with priors on Tkin, , and log10N taken in [5:12], [5:8], and [11:15] respectively. The bottom row shows the results in a similar way as for L1521E.
In the text