Evolution of Chemistry in the envelope of HOt corinoS (ECHOS): I. Extremely young sulphur chemistry in the isolated Class 0 object B 335

Astronomy & Astrophysics A&A, 678, A199 (2023) https://doi.org/10.1051/0004-6361/202346721 © The Authors 2023 Evolution of Chemistry in the envelope of HOt corinoS (ECHOS) I. Extremely young sulphur chemistry in the isolated Class 0 object B 335 G. Esplugues1, M. Rodríguez-Baras1, D. San Andrés2, D. Navarro-Almaida3, A. Fuente2, P. Rivière-Marichalar1, Á. Sánchez-Monge4,5, M. N. Drozdovskaya6, S. Spezzano7, and P. Caselli7 1Observatorio Astronómico Nacional (OAN), Alfonso XII 3, 28014 Madrid, Spain e-mail: [email protected] 2Centro de Astrobiología (CAB), INTA-CSIC, Carretera de Ajalvir Km. 4, Torrejón de Ardoz, 28850 Madrid, Spain 3Département d’Astrophysique (DAp), Commissariat à l’Énergie Atomique et aux Énergies Alternatives (CEA), Orme des Merisiers, Bât. 709, 91191 Gif-sur-Yvette, Paris-Saclay, France 4Institut de Ciències de l’Espai (ICE, CSIC), Can Magrans s/n, 08193, Bellaterra, Barcelona, Spain 5Institut d’Estudis Espacials de Catalunya (IEEC), Barcelona, Spain 6Center for Space and Habitability, Universität Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland 7Max-Planck-Institut für extraterrestrische Physik, 85748 Garching, Germany Received 21 April 2023 / Accepted 30 August 2023 ABSTRACT Context. Within the project Evolution of Chemistry in the envelope of HOt corinoS (ECHOS), we present a study of sulphur chemistry in the envelope of the Class 0 source B335 through observations in the spectral range λ= 7, 3, and 2 mm. Aims. Our goal is to characterise the sulphur chemistry in this isolated protostellar source and compare it with other Class 0 objects to determine the environmental and evolutionary effects on the sulphur chemistry in these young sources. Methods. We have modelled observations and computed column densities assuming local thermodynamic equilibrium and large velocity gradient approximation. We have also used the code Nautilus to study the time evolution of sulphur species, as well as of several sulphur molecular ratios. Results. We have detected 20 sulphur species in B 335 with a total gas-phase S abundance similar to that found in the envelopes of other Class 0 objects, but with significant differences in the abundances between sulphur carbon chains and sulphur molecules containing oxygen and nitrogen. Our results highlight the nature of B335 as a source especially rich in sulphur carbon chains unlike other Class 0 sources. The low presence or absence of some molecules, such as SO and SO+, suggests a chemistry not particularly influenced by shocks. We, however, detect a large presence of HCS+that, together with the low rotational temperatures obtained for all the S species (<15 K), reveals the moderate or low density of the envelope of B 335. Model results also show the large influence of the cosmic ray ionisation rate and density variations on the abundances of some S species (e.g. SO, SO2, CCS, and CCCS) with differences of up to ∼4 orders of magnitude. We also find that observations are better reproduced by models with a sulphur depletion factor of 10 with respect to the sulphur cosmic elemental abundance. Conclusions. The comparison between our model and observational results for B 335 reveals an age of 104<t<105yr, which highlights the particularly early evolutionary stage of this source. B 335 presents a different chemistry compared to other young protostars that have formed in dense molecular clouds, which could be the result of accretion of surrounding material from the diffuse cloud onto the protostellar envelope of B 335. In addition, the theoretical analysis and comparison with observations of the SO2/C2S, SO/CS, and HCS+/CS ratios within a sample of prestellar cores and Class 0 objects show that they could be used as good chemical evolutionary indicators of the prestellar to protostellar transition. Key words. astrochemistry – ISM: abundances – ISM: molecules – ISM: clouds – stars: formation 1. Introduction Hot corinos are regions that arise during the formation of solartype protostars. The progressive collapse of a prestellar core leads to the heating of the infalling envelope of dust and gas to temperatures of hundreds of Kelvin, while the densities increase up to ∼108–109cm−3in the inner ∼100 au around the central protostar. These are also the scales at which protoplanetary disks are expected to arise owing to the conservation of angular momentum. As the temperatures increase, the water-rich ice mantles sublimate, injecting molecules into the gas phase. Hot corinos are therefore chemically rich regions in Solar-like young protostars (Ceccarelli 2004;Ceccarelli et al. 2007b). In particular, these regions are relatively rich in interstellar complex organic molecules (Herbst & van Dishoeck 2009). Studying, therefore, the inventory and abundances of molecules in the earliest stages of star formation, that is in hot corinos, is essential to probe the initial conditions of planet and star formation processes. In spite of the detection of various hot corinos in the last decades (e.g. Cazaux et al. 2003;Okoda et al. 2023), just a few of them (such as IRAS16293-2422, NGC 1333, and IRAS 4A van Dishoeck et al. 1995;Ceccarelli et al. 2007a;Jaber et al. 2014;Caselli & Ceccarelli 2012;Jørgensen et al. 2016) have been investigated in detail until now. Many questions, therefore, still remain open about their nature and molecular composition: for instance, chemical differences between the earliest stages A199, page 1 of 21 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. This article is published in open access under the Subscribe to Open model.Subscribe to A&A to support open access publication. Esplugues, G., et al.: A&A, 678, A199 (2023) of star formation (Class 0 versus Class I objects), the role of the chemical characteristics inherited from the parent cloud, as well as the environmental effect on the chemical composition of hot corinos. One way to understand the dynamics and evolution of the earliest stages of star formation is through molecular depletion. At temperatures of roughly 10 K and densities above 104cm−3, several molecules, such as CO and CS, condense out onto dust grain surfaces (e.g. Caselli et al. 1999;Crapsi et al. 2005), and depletion increases with time (e.g. Bergin & Langer 1997;Aikawa et al. 2003). Therefore, depletion can be used as a clock, since evolved cores should be more depleted of certain species than younger cores. Regarding depletion, sulphur plays a fundamental role since, although sulphur is one of the most abundant species in the Universe and plays a crucial role in biological systems on Earth (e.g. Leustek 2002;Francioso et al. 2020), S-bearing molecules are not as abundant as expected in the interstellar medium (ISM). In particular, unlike the diffuse ISM and photon-dominated regions where the observed gaseous sulphur accounts (e.g. Goicoechea et al. 2006;Howk et al. 2006) for its total solar abundance (S/H ∼1.5 ×10−5;Asplund et al. 2009), sulphur is strongly depleted in more dense molecular gas, since the sum of the detectable S-bearing molecules only accounts for a very small fraction of the elemental S abundance (Tieftrunk et al. 1994). In fact, to reproduce observations in hot corinos and hot cores, one needs to assume a significant sulphur depletion of at least one order of magnitude lower than the solar elemental sulphur abundance (e.g. Wakelam et al. 2004; Esplugues et al. 2014,2022;Crockett et al. 2014;Vastel et al. 2018;Bulut et al. 2021;Navarro-Almaida et al. 2021;Hily-Blant et al. 2022;Fuente et al. 2023). Most of the sulphur seems to be locked on the icy grain mantles (e.g. Millar & Herbst 1990; Ruffle et al. 1999;Vidal et al. 2017;Laas & Caselli 2019) in the form of organo-sulphur molecules, also present in the comet 67P (Calmonte et al. 2016) and in meteoritic material (e.g. Naraoka et al. 2022;Ruf et al. 2021). However, OCS and SO2are the only S-bearing molecules unambiguously detected in ice mantles to date (Geballe et al. 1985;Palumbo et al. 1995;McClure et al. 2023), while H2S has not been detected (Yang et al. 2022; McClure et al. 2023). Theoretical models also suggest that sulphur could be locked into pure allotropes of sulphur, especially S8(Shingledecker et al. 2020). In any case, sulphur reservoirs in dense regions remain as an open question, as well as the time at which sulphur depletion occurs. Given the low number of discovered hot corinos and, therefore, the limited knowledge about them, the project Evolution of Chemistry in the envelope of HOt corinoS (ECHOS) aims to provide a complete inventory of the chemistry of hot corinos (with a sample including IRAS16293-2422, NGC 1333 IRAS 4A, HH 212, and B 335). ECHOS is carried out through homogeneous and systematic spectral surveys between 30 and 200 GHz using the Yebes-40 m and IRAM-30 m telescopes. In this paper, we focus on the analysis of B 335 and present a study of sulphurbearing molecules in this region. Observations and a description of the source are presented in Sects. 2and 3. In Sect. 4, we present the data and calculate column densities and abundances using different methods, depending on the number of detected transitions for each species, and also depending on the availability of collisional coefficients. A discussion of the results (including a comparison with other similar objects, as well as the use of a chemical code to study the fractional abundance evolution of sulphur species) is presented in Sect. 5. We finally summarise our conclusions in Sect. 6. 2. Observations The observations were carried out with the Yebes-40 m and IRAM-30 m telescopes pointing to the position with coordinates αJ2000.0= 19h37m0.93s,δJ2000.0= +07o34′09.9′′. For both telescopes, data were reduced and processed by using the CLASS and GREG packages provided within the GILDAS software1, developed by the IRAM institute. 2.1. Yebes-40 m telescope The observations from the Yebes-40 m radiotelescope located in Yebes (Guadalajara, Spain) were obtained using a receiver that consists of two cold high electron mobility transistor amplifiers covering the 31.0–50.3 GHz Q-band with horizontal and vertical polarisations. The backends are 2 ×8×2.5 GHz fast Fourier transform spectrometers with a spectral resolution of 38.15 kHz, providing the whole coverage of the Q-band in both polarisations. The observations were performed in the positionswitching mode with (-400′′, 0′′) as the reference position. The main beam efficiency varies from 0.6 at 32 GHz to 0.43 at 50 GHz. Pointing corrections were derived from nearby quasars and SiO masers, and errors remained within 2–3′′. 2.2. IRAM-30 m telescope The observations from the IRAM-30 m radio telescope located in Pico Veleta (Granada, Spain) cover the 73–118 GHz and 130–175 GHz spectral ranges. These observations were carried out in one session in July 2021 with the Wobbler-Switching mode (±120′′), using the broad-band EMIR (E090 and E150) receivers and the fast Fourier transform spectrometer in its 200 kHz of spectral resolution. The intensity scale in antenna temperature (T⋆ A, which is corrected for atmospheric absorption and for antenna ohmic and spillover losses) was calibrated using two absorbers at different temperatures and the atmospheric transmission model (ATM; Cernicharo et al. 1985;Pardo et al. 2001). Calibration uncertainties were adopted to be 10%. To convert to main beam brightness temperature (TMB), we used the expression TMB =(Feff/Beff)×T⋆ A=T⋆ A/ηMB,(1) where Feffis the telescope forward efficiency and Beffis the main beam efficiency2,3. See Table 1for a summary of the telescope beam sizes depending on frequency. 3. The source B 335 is an isolated dense globule associated with the embedded far-infrared source IRAS 19347+0727 (Keene et al. 1980, 1983), located at a distance of ∼100 pc (Olofsson & Olofsson 2009). The central source is a Class 0 source with a luminosity of 0.72 L⊙(Evans et al. 2015), whose mass is estimated to be 0.08 M⊙(Yen et al. 2015). Complex organic molecules (COMs), such as acetaldehyde (CH3CHO), methyl formate (HCOOCH3), and formamide (NH2CHO), have been detected (Imai et al. 2016; Okoda et al. 2022) in the vicinity of the protostar (r∼10 au). It reveals that B 335 harbors a hot corino (Cazaux et al. 2003; Herbst & van Dishoeck 2009), which warms up the surrounding regions leading to the sublimation of grain mantles and injecting 1http://www.iram.fr/IRAMFR/GILDAS 2https://rt40m.oan.es/rayo/index.php 3http://www.iram.es/IRAMES/mainWiki/Iram30mEfficiencies A199, page 2 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) Table 1. Yebes-40 m and IRAM-30 m telescope efficiency data along the covered frequency range. Telescope Frequency ηMB HPBW (GHz) (′′) Yebes-40 m 32.4 0.61 54.4 34.6 0.58 51.0 36.9 0.57 47.8 39.2 0.55 45.0 41.5 0.53 42.5 43.8 0.52 40.3 46.1 0.49 38.3 48.4 0.47 36.4 IRAM-30 m 86 0.82 29.0 100 0.79 22.0 145 0.74 17.0 170 0.70 14.5 complex organic molecules into the gas-phase (e.g. Ceccarelli 2004;Caselli & Ceccarelli 2012). Emission of CCH and c-C3H2 has also been detected around the protostar, which indicates the presence of carbon-chain molecules in this source (Sakai et al. 2008;Sakai & Yamamoto 2013). This carbon-chain molecules could suggest the presence of a so-called warm carbon-chain chemistry (WCCC; Sakai et al. 2008), characterised by various carbon-chain molecules being produced by reactions of CH4 (evaporated from grain mantles in a lukewarm region) and where the formation of HCnN is slower than that of CnHmmolecules (Sakai et al. 2008,2009;Oya et al. 2017). Particularly interesting are the single-dish observations of CS and H2CO in B335 (Zhou et al. 1993;Choi et al. 1995), since their optically thick lines show asymmetric profiles. This implies the presence of infalling motions in the envelope around B335. A more detailed analysis with C18O, SO, HCN, and HCO+ ALMA observations (Evans et al. 2015,2023;Yen et al. 2015) confirms an infalling gas motion of the inner part of the protostellar core, but rotation motion has not been detected at a spatial resolution of 0.3′′ (∼30 au). A rotation structure is only found at a scale of about 10 au from CH3OH and HCOOH emission (Imai et al. 2019). Single-dish and interferometric 12CO (1–0) and 13CO (1– 0) observations (e.g. Hirano et al. 1988,1992;Stutz et al. 2008;Bjerkeli et al. 2019) show also the presence of an east-west molecular outflow on ∼0.2 pc (Cabrit et al. 1988;Moriarty-Schieven & Snell 1989) with an inclination between 3o–10oon the plane of the sky and an outflow opening angle of 45o(Hirano et al. 1988,1992;Stutz et al. 2008;Yen et al. 2010). A dynamical analysis of the outflow molecular structure (Yıldız et al. 2015) implies that dynamical timescales for the CO emitting gas are of the order of 104yr. This result, along with the absence of a rotationally supported disk on scales greater than ∼10 au in size (Evans et al. 2015;Yen et al. 2015), suggest that the B 335 system is a particularly young source compared to other Class 0 objects, such as L1527 which contains an embedded IRAS source and has a dynamical age of 2 ×104–105yr (Tamura et al. 1996;André et al. 2000;Agúndez et al. 2019). The detection of abundant CCS (see Sect. 4.3) and other carbon-chain molecules in B335 also confirms the early stage of this source since it is generally recognised that CCS and other carbon-chain molecules are abundant in the early stage of chemical evolution and become deficient at the advanced stage (Sakai et al. 2008). So far, B 335 is one of the few hot corinos identified in a Bok globule, isolated from a large molecular cloud complex. This is therefore an ideal source for the study of the earliest stages of star formation with a chemical composition that could be regarded as a standard template for isolated protostellar cores. 4. Data analysis and results We have detected 81 emission lines from 21 different sulphurbearing species (including isotopologs) in B 335. In particular, we have detected three transitions of CS, three of 13CS, three of C34S, three of C33S, nineteen of CCS, three of CC34S, seven of CCCS, three of CCC34S, five of OCS, three of HCS+, two of HC34S+, seven of SO, one of SO2, one of 34SO, one of H2S, one of H234S, ten of H2CS, one of HSCN, one of HNCS, and six of NS (see Table A.1 and Figs. A.1–A.21). These transitions span an energy range of Eup =2.0–65.3 K. 4.1. Line profiles All the lines are observed in emission. We first fitted the observed lines with Gaussian profiles using the CLASS software to derive the radial velocity (vLSR), the line width, and the intensity for each line. Results are given in Table A.1. We observed that most of the detected molecules show narrow line profiles (∆v≲1.5 km s−1) that can be adjusted by a single velocity component. It is well known that the gravitational collapse of a molecular cloud when a new star is formed is accompanied by the development of highly supersonic outflows (e.g. Bachiller & Pérez Gutiérrez 1997;Holdship et al. 2019). The material ejected through these jets collides with the surrounding cloud, compressing and heating the gas, which leads to a drastic alteration of the chemistry as endothermic reactions become efficient and dust grains are partially destroyed. Many observations have shown that several sulphur-bearing species are enhanced in the shocked regions with respect to other areas of the cloud (e.g. Pineau des Forets et al. 1993;Bachiller & Pérez Gutiérrez 1997;Codella & Bachiller 1999;Wakelam et al. 2004,2005; Esplugues et al. 2013). However, the sulphur lines observed in B 335 present very narrow line widths (Figs. A.1–A.21), which is more consistent with the emission arising from the ambient quiescent cloud. Only for H2S, formed through the hydrogenation of sticked atomic sulphur on grains (Hatchell et al. 1998; Garrod et al. 2007;Esplugues et al. 2014;Oba et al. 2019) and thought to be the main sulphur reservoir in ices (Vidal et al. 2017;Navarro-Almaida et al. 2020), it is necessary to also consider a wider component (with ∆v = 2.5km s−1) to fit its only detected line, since the line profile shows the presence of faint wings. This could indicate that part of the H2S emission arises from the outflow, however, the line width of this component is much smaller (∼2.5 kms−1) than the ones observed for 12CO (∼10–20 km s−1) in the outflow (Yen et al. 2010). Therefore, the wider H2S emission may be likely related to infall motions. On the other hand, there are some species (CCS, CCCS, and H2CS) whose low transition lines present a double peak emission. The peak with the lowest intensity is located for all the cases at lower velocity (vLSR =7.39–8.1 km s−1, blue peak) than for the more intense peak (vLSR =8.80–8.97 km s−1, red peak). Line widths are slightly larger (∆v≲0.83–1.63 km s−1) in the blue peaks than in the red ones (∆v≲0.55–0.83 km s−1). The comparison of these profiles with those from the less abundant A199, page 3 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) isotopologues shows that these two-peak profiles are due to selfabsorption due to the large optical depths. All the detected lines present 0.005 <TMB <1.5 K. 4.2. Rotational diagrams For each molecule, we have computed a representative rotational temperature (Trot) and a column density (N) by constructing a rotational diagram, assuming a single rotational temperature for all energy levels (Goldsmith & Langer 1999). Rotational diagrams represent a technique to analyse cloud properties from molecular line emission assuming local thermodynamic equilibrium (LTE). The standard relation for the rotational diagram analysis is ln Nu gu!=lnN−lnZ−Eu kTrot ,(2) with Nu/gugiven by Nu gu =8kπ hc3× ν2 ul Aulgu ×f−1×ZTMBdv,(3) where Nuis the column density of the upper level in the optically thin limit, Nis the total column density, guis the statistical weight of the upper state of each level, Zis the partition function evaluated at a rotational temperature Trot,Eu/kis the energy of the upper level of the transition, νul is the frequency of the u  l transition, RTMB dv is the velocity-integrated line intensity corrected from beam efficiency, and fis the beam filling factor. Assuming that the emission source has a 2D Gaussian shape, f is equal to fbf =θ2 s/(θ2 s+θ2 b) being θbthe HPBW of the telescope in arcsec and θsthe diameter of the Gaussian source in arcsec. In order to derive the values of fin each case, we have used the HPBW values listed in Table 1, and we have assumed a source size of 10′′ according to interferometric maps of C17O emission in B 335 (Cabedo et al. 2021). The resulting rotational diagrams are shown in Fig. 1, while the column densities and rotational temperatures obtained are listed in Table 2. The uncertainties shown in Table 2indicate the uncertainty obtained in the least squares fit of the rotational diagrams. These uncertainties include the 10%error in the observed line intensities due to calibration. The uncertainty obtained in the determination of the line parameters with the Gaussian fitting programme is included within the error bars at each point of the rotational diagram. We did not applied the rotational diagram method to NS since its four observed lines are close to the detection limit and/or blended with other species, which would introduce a high uncertainty in the results. The analysis of NS and others molecules with only one emission line detected (e.g. SO2) has been done only with a LVG code (see next section) instead. The obtained rotational temperatures are ≲15 K. In the case of CCS, the molecule with the largest number of detected transitions, we observe a trend variation (Fig. 1) between the low and the high transition energies. In particular, the transitions with Eu< 40 K are well adjusted for a Trot ∼5 K, while transitions with Eu> 40 K are better adjusted for a Trot ∼15 K. This suggests that CCS arises from different regions characterised by distinct excitation temperatures. The rotational temperatures for the rest of molecules with several detected transitions (e.g. C3S, OCS, SO, H2CS) have 4.5 ≲Trot ≲7.4 K. These low temperatures together with their narrow line widths are consistent with emission arising from a cold and extended part of the dense core. Table 2. Column densities (N), rotational temperatures, and abundances (X) relative to H2of the sulphur-bearing molecules detected in B335. Species Nrot ×1013 Trot (a)NLVG ×1013 (b)X×10−9 (cm−2) (K) (cm−2) CS 8 ±2 3.6 ±0.8 12.00 (34.2) 11.0 13CS 0.8 ±0.1 2.8 ±0.9 0.73 0.24 C34S 2.1 ±0.8 2.7 ±0.7 1.40 0.45 C33S 0.18 ±0.09 2.9 ±0.5 0.13 0.04 CCS 5.4 ±0.8 5 ±0.5 4.40 1.42 0.19 ±0.06 15 ±1 CC34S 0.27 ±0.01 6.4 ±0.3 0.90 0.29 CCCS 1.8 ±0.4 7.4 ±0.8 0.65 0.02 CCC34S 0.14 ±0.07 2.8 ±0.9 0.12 0.04 OCS 3.8 ±0.5 7.0 ±0.8 2.00 0.65 HCS+0.5 ±0.2 3.8 ±0.4 0.55 (1.5) 0.05 HC34S+0.08 ±0.02 6 ±1 0.06 0.01 SO 1.8 ±0.6 7.2 ±0.1 6.10 1.97 34SO – – 0.25 0.08 o-H2CS 1.6 ±0.7 4.5 ±0.5 1.00 0.32 p-H2CS 0.6 ±0.2 5.4 ±0.2 0.70 0.22 NS – – 0.13 0.04 SO2– – 0.47 0.21 H2S – – 6.50 2.10 H234S – – 0.11 0.04 HSCN – – 0.10 0.03 HNCS – – 0.47 0.15 Total N6.2 ×1014 Total X2.0 ×10−8 Notes. Numbers in parentheses are obtained from its corresponding 34S isotopologue considering 32S/34S = 24.4 (Mauersberger et al. 2004). (a)LVG calculations assuming nH2= 5 ×104cm−3and TK= 15 K (Shirley et al. 2011;Cabedo et al. 2021). (b)Considering NH2= 3.1 × 1022 cm−2(Cabedo et al. 2021). Gas density at these scales is probably not too high as also indicated by the low Trot obtained for most of the S-molecules. All this suggests that the gas is sub-thermally excited with a moderate or low (<105cm−3) density, in agreement with Frerking et al. (1987). In addition, the low number of transitions detected for several species (CS, 13CS, C34S, C33S, HCS+, and HC34S+) and the optically thick emission of some transitions of CCS, CCCS, and H2CS, also indicate that the rotational diagram method not suitable to derive column densities. 4.3. Column densities and abundances We have therefore calculated molecular column densities for the detected S-species using a large velocity gradient (LVG) approximation. In particular, we have derived them by fitting the emission line profiles using the LVG code MADEX (Cernicharo 2012), which assumes that the radiative coupling between two relatively close points is negligible, and the excitation problem is local. The LVG models are based on the Goldreich & Kwan (1974) formalism. MADEX also includes corrections for the beam dilution of each line, depending on the different beam sizes at different frequencies. The final considered fit is the one that reproduces more line profiles better from all the observed A199, page 4 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) 0 2 4 6 8 10 12 14 E u (K) 25 26 27 28 29 30 31 32 Ln(N u /g u ) CS N =(8±2)1013 cm−2 Trot =(3.6±0.8) K 0 2 4 6 8 10 12 14 E u (K) 22 23 24 25 26 27 28 29 Ln(N u /g u ) 13CS N =(8±1)1012 cm−2 Trot =(2.8±0.9) K 0 2 4 6 8 10 12 14 E u (K) 24 25 26 27 28 29 30 Ln(N u /g u ) C34S N =(2.1±0.8)1013 cm−2 Trot =(2.7±0.7) K 0 2 4 6 8 10 12 14 E u (K) 21 22 23 24 25 26 27 28 Ln(N u /g u ) C33S N =(1.8±0.9)1012 cm−2 Trot =(2.9±0.5) K 0 10 20 30 40 50 60 70 E u (K) 16 18 20 22 24 26 28 30 Ln(N u /g u ) CCS N =(5.4±0.8)1013 cm−2 Trot =(5.0±0.5) K N =(1.9±0.6)1012 cm−2 Trot =(15±1) K 0 2 4 6 8 10 12 14 E u (K) 23.0 23.5 24.0 24.5 25.0 25.5 26.0 Ln(N u /g u ) CC34S N =(2.7±0.1)1012 cm−2 Trot =(6.4±0.3) K 0 5 10 15 20 25 30 35 40 E u (K) 19 20 21 22 23 24 25 26 27 Ln(N u /g u ) CCCS N =(1.8±0.4)1013 cm−2 Trot =(7.4±0.8) K 0 2 4 6 8 10 12 E u (K) 20.0 20.5 21.0 21.5 22.0 22.5 23.0 23.5 24.0 Ln(N u /g u ) CCC34S N =(1.4±0.7)1012 cm−2 Trot =(2.8±0.9) K 0 5 10 15 20 25 30 E u (K) 23 24 25 26 27 28 29 Ln(N u /g u ) OCS N =(3.8±0.5)1013 cm−2 Trot =(7.0±0.8) K 0 5 10 15 20 25 E u (K) 21 22 23 24 25 26 27 28 29 Ln(N u /g u ) HCS+ N =(5±2)1012 cm−2 Trot =(3.8±0.4) K 0 5 10 15 20 E u (K) 20 21 22 23 24 25 26 27 Ln(N u /g u ) HC34S+ N =(8±2)1011 cm−2 Trot =(6±1) K 0 5 10 15 20 25 30 35 40 E u (K) 23 24 25 26 27 28 Ln(N u /g u ) SO N =(1.8±0.6)1013 cm−2 Trot =(7.2±0.1) K 0 5 10 15 20 25 E u (K) 21 22 23 24 25 26 27 28 29 30 Ln(N u /g u ) o-H2CS N =(1.6±0.7)1013 cm −2 Trot =(4.5±0.5) K 0 5 10 15 20 25 E u (K) 20 21 22 23 24 25 26 27 Ln(N u /g u ) p-H2CS N =(6±2)1012 cm −2 Trot =(5.4±0.2) K Fig. 1. Rotational diagrams of the detected sulphur-bearing molecules in B335. Fitted values of the rotational temperature, Trot, column density, N, and their respective uncertainties are also indicated for each molecule. transitions. To carry out the fits, we assume uniform physical conditions (kinetic temperature, density, line width, radial velocity). In particular, we have assumed a gas density of 5×104cm−3 according to results from Sect. 4.2 and to results from Cabedo et al. (2023), who deduced a density of ∼104cm−3for the outer regions of B335 based on interferometric observations. We have also assumed a source size of 10′′ according to interferometric maps (Cabedo et al. 2021), and a kinetic temperature of the gas TK∼15 K estimated by Shirley et al. (2011) from the formula for dust temperature considering only central heating by the B 335 protostar, and that dust and gas are coupled. In the case of H2CS and H2S, we have used the code RADEX (van der Tak et al. 2007) since there are not available collisional rates in MADEX for these two molecules. The results are shown A199, page 5 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) in Table 2. For those species for which one transition is observed (34SO, SO2, H2S, H234S, HSCN, and HNCS), we have derived their column densities using only the LVG approximation and not the rotational diagram technique described in Sect. 4.2. For NS, since there are not available collisional rates, we have considered LTE approximation. Column densities derived using the rotational diagram technique and LVG calculations agree within a factor ≲3. Sources of uncertainty are the low angular resolution of the used telescopes (with beam sizes between ∼36′′ and ∼54′′ for the 40 m-Yebes telescope, and ∼14′′ and ∼29′′ for the 30m-IRAM telescope depending on frequency), which implies that the emission from the inner region of the B 335 is blended with the outer envelope, the possible overlap with the emission from other species (e.g. the case of H234S), the limited number of detected transitions (one transition in many cases, such as for H2S, H234S, SO2, HSCN, HNCS, and 34SO), the lack of collisional rates for some species (e.g. NS), and the assumed source size derived from C17O interferometric observations. This is one of the most important uncertainties due to the lack of interferometric observations of S-species in B335, since the emission from the different S-molecules may arise from distinct regions with slightly different sizes leading to a significant impact on the fits. In order to derive the uncertainty associated with the LVG results, we have run the LVG code varying the considered values for density, temperature, and source size. In particular, we have run models with density 5 ×104cm−3and 2.5 ×104cm−3(variation of a factor 2), temperatures of 15 K and 12 K, and source sizes of 10′′ and 8′′ (reduction of 20%). We found variations in the S-column densities ≲10%between these two TKwhich is comparable to the calibration errors, while the variations when changing the density and the source size were ≲20%and ∼30%, respectively. Taking all this into account, we estimate the uncertainties in the column densities derived with LVG approximation to be of the order of ∼50%. The total column density of the sulphur-bearing species observed in B 335 is 6.2 ×1014 cm−2, with the highest contribution coming from CS (calculated from C34S considering 32S/34S = 24.4 since CS is optically thick, Mauersberger et al. 2004), H2S, SO, and CCS. The lowest values are found for NS and HSCN. We have also observed several isotopologues relative to 13C, 34S, and 33S. In particular, from CS, we obtain 34S/33S = 10.7 (which is slightly higher than the solar abundance value 34S/33S = 5.6–6.3, Anders & Grevesse 1989;Mauersberger et al. 2004), and 12C/13C = 46.8, which is similar to that found in dark clouds by Cernicharo & Guelin (1987). Table 2also shows the fractional abundances relative to H2 (with NH2obtained from Cabedo et al. 2021) of the sulphur molecules detected in the range λ= 2, 3, and 7 mm. Considering the detected S-molecules, we have obtained a total molecular sulphur abundance of 2.0 ×10−8, which is alike to the one found in similar sources, such as L 1544 (1.1 ×10−8, Vastel et al. 2018), L 483 (9.1 ×10−9,Agúndez et al. 2019), and B1-b (1.2 ×10−8,Fuente et al. 2016). The highest S-abundances are >10−9(CS, H2S, CCS, and SO), with also several relatively abundant (>10−10) species, such as OCS and H2CS. In general, B 335 presents a great variety of S-bearing molecules, but less than other similar sources. In particular, although the total sulphur abundance in B 335 is similar to other sources, such as L483, we detect a significant lower number of sulphur molecules in B 335 than in L483. We are also missing some other S-species already reported towards L 483, including SO+, HCS, and CH3SH, for which we have derived upper limits for their column densities of N(SO+)<1.5 ×1011 cm−2, CS 13CS C34S C33S CCS CC34S CCCS CCC34S OCS SO 34SO SO2 HCS+ HC34S+ H2CS NS HSCN H2S H234S SO+ NS+ 10−12 10−11 10−10 10−9 10−8 10−7 N / N (H2) B335 L483 B1 -b Fig. 2. Comparison of fractional abundances of sulphur-bearing species between three Class 0 objects B335 (this work), L 483 (Agúndez et al. 2019), and B1-b (Fuente et al. 2016). N(HCS) <3×1012 cm−2, and N(CH3SH) <2×1012 cm−2. To derive these upper limits, once the physical parameters are fixed, we vary (for each species) the column density until the model fit reaches the observed intensity peak of any of the observed lines. We do not allow the model fit to be greater than any observed line. 5. Discussion 5.1. Chemical comparison Even though the total sulphur abundance in B 335 is similar to that from other Class 0 sources, the abundances of some specific molecules significantly vary among them. Figure 2shows a comparison between abundances of different sulphur species in three Class 0 sources, B335, L 483, and B1-b. The first thing that stands out is that sulphur carbon chains (CS, CCS, and CCCS) in B 335 are about one order of magnitude higher than the ones observed in L 483 and B1-b. Although high abundances of carbon-chain molecules have already been found in other cold clouds, such as TMC1 (Agúndez & Wakelam 2013), it is not a universal characteristic of this type of early regions, since many of them present different degrees of carbon chain richness (e.g. Suzuki et al. 1992;Hirota et al. 2009). This therefore highlights the nature of B 335 as a source especially rich in sulphur carbon chains. Unlike carbon chains, sulphur molecules containing oxygen in B 335 have a similar (OCS) or even lower (SO and SO2) abundance by up to one order of magnitude than those found in L 483 and B1-b (see also Fig. 3). These oxygen-sulphur molecules are usually detected in cold dense regions with abundances between 10−10–10−8, while in hot cores and hot corinos they are significantly enhanced (e.g. Tercero et al. 2010;Esplugues et al. 2013). In particular, SO2is a better tracer of warm gas than SO (Esplugues et al. 2013), while SO is a well-known outflow tracer (e.g. Codella & Scappini 2003;Lee et al. 2010;Tafalla et al. 2010). SO seems to be more enhanced than SO2in shocks for timescales ∼104yr (e.g. Codella & Bachiller 1999;Viti et al. 2004;Jiménez-Serra et al. 2005). The low SO abundance found in B 335 with respect to the other two sources L483 and B1-b suggests that its chemistry may not be significantly affected by shocks associated with the bipolar outflow (e.g. Hirano et al. 1988,1992;Stutz et al. 2008;Bjerkeli et al. 2019). Other sulphur molecules detected in B 335 with low abundances (Fig. 3) are those containing nitrogen, that is NS and A199, page 6 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) CS 13CS C34S C33S CCS CC34S CCCS OCS SO 34SO SO2 HSCN HCS+ HC34S+ H2CS NS H2S 10−2 10−1 100 101 102 XB 335 / XL 483 XB 335 / XB 1 − b Fractional ab ndance (X) ratio Fig. 3. Fractional abundance (X) ratios between B 335 and L483, and B 335 and B1-b for all the detected sulphur-bearing species in B 335. HSCN (4 ×10−11 and 3 ×10−11, respectively). HSCN, unlike HNCS which is the most stable isomer, is metastable and only detected on Earth under specific conditions (through UVphotolysis of HNCS Wierzejewska & Mielke 2001 or through a low-pressure discharge Brünken et al. 2009). In contrast, HSCN is more easily observed in the interstellar medium. For instance, in Sgr B2 (where HSCN was first time detected) the ratio HNCS/HSCN∼3 (Halfen et al. 2009) and in TMC1 this ratio is ∼1 (Adande et al. 2010), while in L 1544 Vastel et al. (2018) only detect HSCN, but not HNCS. This shows the important role of non-equilibrium chemistry in dense clouds. In B 335, we only detect one transition of HNCS with a signal to noise ratio σ∼3.1, so we have considered its column density (NHNCS ≤4.7 ×1012 cm−2) as an upper limit. This leads to a ratio HNCS/HSCN <4.7 consistent with the one found in Sgr B2. 5.2. Sulphur ions Apart from neutral sulphur species, we have also observed the presence of some ions in B335. Molecular ions are commonly observed in dense prestellar or protostellar regions (e.g. Caselli et al. 1998,2002;Hogerheijde et al. 1998;van der Tak et al. 2005). To a lesser extent, ions (such as HCO+and N2H+) have also been observed in protostellar shocks (e.g. Bachiller & Pérez Gutiérrez 1997;Hogerheijde et al. 1998;Caselli et al. 2002; van der Tak et al. 2005). The presence of ions in shocked regions is a balance between their enhancement by the sputtering of dust grain mantles, which ejects molecules into the gas phase that undergo reactions forming ions, and their destruction by electronic recombination (e.g. Neufeld & Dalgarno 1989;Viti et al. 2002). For this reason, those ions that usually decrease in shocks (such as N2H+) may indicate the presence of a pre-shock chemistry (Codella et al. 2013), while the ions that are enhanced by shock chemistry will be effective tracers of shocked gas. This is the case of HOCO+and SO+(Podio et al. 2014) since their abundances are predicted to be very low in the quiescent gas, but significantly enhanced in shocks following the release of CO2 and S-species from dust grain mantles (Minh et al. 1991;Turner 1992,1994;Deguchi et al. 2006). Among the sulphur ions, in B 335 we do not detect SO+ (strong tracer of shocked gas which forms from S+and OH Neufeld & Dalgarno 1989), but we detect HCS+and HC34S+ with abundances of 4.7 ×10−10 and 1.9 ×10−11, respectively. The value we have obtained for HCS+is similar to the one observed in the bow shock L1157-B1 (Podio et al. 2014) and one order of magnitude higher than the ones obtained by Agúndez et al. (2019) and Fuente et al. (2016) in the Class 0 objects L 483 and B1-b, respectively, as also found for CS. The abundances of HCS+and CS are strongly related, since this ion forms through the reaction of CS with HCO+, H+ 3, H3O+(Millar et al. 1985). In addition, large abundances of CS and HCS+may also indicate that OCS is one of the main sulphur carriers on dust grains upon sputtering from the grain mantles according to previous theoretical (Podio et al. 2014) and observational (Wakelam et al. 2005; Codella et al. 2005) results. Regarding the ratio R= HCS+/CS, we obtain from observations a value of 0.043. Thaddeus et al. (1981) determined an abundance ratio R∼0.01–0.03 in ten molecular sources located in Sgr B2 and Ori A, while Millar (1983) found that Rcould be as large as 0.1 in cold dark clouds. This result was also verified by Irvine et al. (1983). A similar Rvalue to that from B 335 is found as well in the pre-stellar core L1544 (0.03, Vastel et al. 2018), while in the protostellar sources L 483 and B1-b the HCS+/CS ratio is 0.017 and 0.008, respectively (Agúndez et al. 2019;Fuente et al. 2016). According to Clary et al. (1985), rate coefficients for ion-CS reactions rapidly increase at low temperatures suggesting that the sources B 335 and L 1544 are characterised by lower temperatures than L 483 and B1-b, which can be due to an earlier evolutionary stage of B 335 and L 1544 compared to L 483 and B1-b. 5.3. Time evolution of fractional abundances Figures 4–7show the evolution of several S-bearing species along 10 Myr for different physical conditions of the cosmicray (CR) ionisation rate (ζH2), gas temperature (Tg), hydrogen number density (nH), and initial sulphur abundance (S+ init). Given that the CR ionisation rate is still uncertain and spreads over a range of values (van der Tak & van Dishoeck 2000) determined as ζH2∼3×10−17 s−1in dense clouds, Indriolo & McCall (2012) obtained a range (1.7–10.6) ×−16 s−1in a sample of diffuse clouds, and Neufeld & Wolfire (2017) derived a ζH2of the order of a few 10−16 s−1in the Galactic disk), we have considered two values for the cosmic ionisation rate ζH2= 1.3 ×10−16 s−1 and 1.3 ×10−17 s−1. For the gas temperature and density, given the results obtained in Sects. 4.2 and 4.3, we have run models considering Tg= 7K and 15 K, and nH= 2 ×104cm−3 and 2 ×105cm−3. For the case of the initial sulphur abundance, we have considered the solar elemental sulphur fractional abundance (S+ init = 1.5 ×10−5) and also one factor of S depletion (1.5 ×10−6) since to reproduce observations in hot corinos one needs to assume a significant sulphur depletion of at least one order of magnitude lower than the solar elemental sulphur abundance as previously stated in Sect. 1. We have obtained these models using the Nautilus timedependent chemical code (Ruaud et al. 2016). Nautilus is a three-phase model in which gas, surface icy mantle, bulk icy mantle, and their interactions are considered. Nautilus solves the kinetic equations for both the gas phase species and the surface species of interstellar dust grains and computes the evolution with time of chemical abundances for a given physical structure. The chemical network is based on the KInetic Database for Astrochemistry (KIDA4). The used version here is detailed in Wakelam et al. (2021). In particular, this Nautilus version considers gas-phase processes including neutral-neutral and ion-neutral reactions, direct cosmic ray ionisation or dissociation, ionisation or dissociation by UV photons, ionisation or dissociation produced by photons induced by cosmic-ray interactions with the 4https://kida.astrochem-tools.org/ A199, page 7 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - ζ =1.3x10 −17 s−1 -- ζ =1.3x10 −16 s−1 CS CCS C3S OCS 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - ζ =1.3x10 −17 s−1 -- ζ =1.3x10 −16 s−1 SO SO2 HCS+ H2CS NS 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - T g=7 K -- T g=15 K CS CCS C3S OCS 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - T g=7 K -- T g=15 K SO SO2 HCS+ H2CS NS 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - n H=2x104 cm −3 -- n H=2x105 cm −3 CS CCS C3S OCS 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - n H=2x104 cm −3 -- n H=2x105 cm −3 SO SO2 HCS+ H2CS NS Fig. 4. Evolution of fractional abundances of CS, CCS, C3S, OCS, SO, SO2, HCS+, H2CS, and NS as a function of time for an initial sulphur abundance S+ init = 1.5 ×10−6, a hydrogen number density nH= 2 ×104cm−3,Tgas = 7 K, and two values of CR ionisation rate ζ= 1.3 ×10−17 s−1(solid line) and ζ= 1.3 ×10−16 s−1 (dashed line). 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - T g=7 K -- T g=15 K CS CCS C3S OCS 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - T g=7 K -- T g=15 K SO SO2 HCS+ H2CS NS Fig. 5. Evolution of fractional abundances of CS, CCS, C3S, OCS, SO, SO2, HCS+, H2CS, and NS as a function of time for an initial sulphur abundance S+ init = 1.5 ×10−6, a hydrogen number density nH= 2 ×104cm−3,Tgas = 7 K, and two values of CR ionisation rate ζ= 1.3 ×10−17 s−1(solid line) and ζ= 1.3 ×10−16 s−1 (dashed line). 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - n H=2x104 cm −3 -- n H=2x105 cm −3 CS CCS C3S OCS 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - n H=2x104 cm −3 -- n H=2x105 cm −3 SO SO2 HCS+ H2CS NS Fig. 6. Evolution of fractional abundances of CS, CCS, C3S, OCS, SO, SO2, HCS+, H2CS, and NS as a function of time for an initial sulphur abundance S+ init = 1.5 ×10−6, a hydrogen number density nH= 2 ×104cm−3,Tgas = 7 K, and two values of CR ionisation rate ζ= 1.3 ×10−17 s−1(solid line) and ζ= 1.3 ×10−16 s−1 (dashed line). 101102103104105106107 t (yr) 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3 n species/ n H - S + init=1.5x10 −6 -- S + init=1.5x10 −5 CS CCS C3S OCS 101102103104105106107 t (yr) 10 17 10 15 10 13 10 11 10 9 10 7 10 5 10 3 n species/ n H - S + init=1.5x10 −6 -- S + init=1.5x10 −5 SO SO2 HCS+ H2CS NS Fig. 7. Evolution of fractional abundances of CS, CCS, C3S, OCS, SO, SO2, HCS+, H2CS, and NS as a function of time for an initial sulphur abundance S+ init = 1.5 ×10−6, a hydrogen number density nH= 2 ×104cm−3,Tgas = 7 K, and two values of CR ionisation rate ζ= 1.3 ×10−17 s−1(solid line) and ζ= 1.3 ×10−16 s−1 (dashed line). medium (Prasad & Tarafdar 1983), and electronic recombinations. For species on the surfaces, there is a distinction between species in the most external layers (surface species) and species below these layers (mantle species). The species are adsorbed on the surface and become the mantle during the construction of the ices. Similarly, when the species desorb, only the species from the surface can desorb, but surface species are gradually replaced by the mantle species. Regarding surface reactions, the model includes thermal desorption, photo-desorption, chemical desorption, and cosmic-ray heating. See Wakelam et al. (2021) for extensive details about the surface parameters considered in each surface reaction. In all models, we have adopted the initial abundances shown in Table 3and a visual extinction AV= 15 mag5. In general, we observe that the abundances of all the species considered in Figs. 4–7increase with time up to t∼5×105yr, moment at which the abundances of various species decrease 5Derived from results in Table 2and the expression Ni=Xi×AV×1.6×1021,(4) where Niis the column density of the species i,Xiits abundance, AVthe visual extinction, and 1.6 ×1021 the hydrogen column density at 1 mag of extinction (Bohlin et al. 1978). A199, page 8 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) Table 3. Abundances with respect to total hydrogen nuclei considered in the chemical code Nautilus. Species Abundance Reference He 9.0 ×10−2(1) O 2.4 ×10−4(2) Si+8.0 ×10−9(3) Fe+3.0 ×10−9(3) S+1.5 ×10−5, 1.5 ×10−6(1) Na+2.0 ×10−9(3) Mg+7.0 ×10−9(3) P+2.0 ×10−10 (3) Cl+1.0 ×10−9(3) F+6.7 ×10−9(4) N 6.2 ×10−5(5) C+1.7 ×10−4(5) References. (1) Taken from Asplund et al. (2009) and Wakelam & Herbst (2008). (2) Taken from Hincelin et al. (2011). (3) Taken from Graedel et al. (1982). (4) Taken from Neufeld et al. (2005). (5) Taken from Jenkins (2009). abruptly. This is the case of CS, SO, and SO2, whose abundances decrease by about 2–4 orders of magnitude between t∼5×105yr and t∼5×106yr since they react with other species to form more complex molecules. In particular, for t>105yr, CS SO, and SO2are mainly destroyed by reacting with H+ 3and forming HCS+, HSO+, and HSO+ 2, respectively. For the abundances of the rest of sulphur species, this decrease after reaching the maximum value is lower and not greater than two orders of magnitude. From these figures, we also deduce the formation time scales of S-species. In particular, we have observed that CS is one of the molecules reaching first its maximum abundance value in all the models at t∼5×104yr, followed by HCS+with its maximum abundance reached at t∼105yr. Other molecules, such as H2CS, can be considered as late molecules, since their maximum abundance values are found when the chemical evolution exceeds one million years. Regarding the influence of the physical conditions (cosmicray ionisation rate, temperature, density, and initial sulphur abundance) on the S-abundances, the parameter with the highest impact is the hydrogen nuclei number density, nH, since varying nHby one order of magnitude leads to abundance variations of up to three orders of magnitude for CS, CCS, OCS, HCS+, SO, and SO2given a specific evolutionary time (see Fig. 6). In general terms, we observe in Fig. 6that the chemical evolution of each considered species is accelerated when the density increases since collisions are more frequent than in a low density regime. This results in the abundance peak of the species being reached at earlier times for larger densities, although the chemical behaviour along time is roughly maintained. In particular, for t∼104yr, we obtain that SO2is mainly formed by the reaction O+SO for nH= 2 ×104and 2 ×105cm−3, while for t∼105yr (when the SO2peak is reached) this molecule mainly forms by the reaction O+SO for nH= 2 ×104cm−3and by OH+SO for nH= 2 ×105cm−3. The most responsible species of destroying SO2for both densities is C for t∼104yr, and H+ 3 for t≥105yr. Regarding CS, for t∼104yr, it is mainly formed by the reaction HCS++e−for the low density value, while CS is mainly formed through the reaction between atomic S and C4 for nH= 2 ×105cm−3. For t∼105yr (when the CS abundance decreases for any model), CS is mainly destroyed by its reaction with H+ 3and with HCO+. Another parameter with a big impact on S-abundances is the cosmic-ray ionisation rate, ζ, since the variation from ζ= 1.3 ×10−17 to 1.3 ×10−16 s−1changes by about two orders of magnitude the abundances of all the sulphur species considered in the sample (Fig. 4). For the case of carbon-chains (CCS, C3S), we observe that the larger ζ, the larger the S-abundances for t≲105yr. Similar results are found for the species SO, SO2, NS, H2CS, and OCS. For t∼104yr, we find that CCS is mainly formed by C+HCS for the low ζmodel, while it is mainly formed through HC3S++e−in a more efficient way for a model with larger ζ(1.3 ×10−16 s−1). In the case of SO (one of the most affected molecules by ζbetween 104≲t≲106yr), it is mostly formed by the reaction O+HS in a model with low ζ, while the reaction between S+OH becomes more important to form SO2 when increasing ζby one order of magnitude. On the other hand, changing the temperature (from Tgas = 7K to 15 K) and especially the initial sulphur abundance (from S+ init = 1.5 ×10−5to 1.5 ×10−6) in the models leads to abundance variations smaller than two orders of magnitude for t>104yr. In particular, we have found that an increase of the S+ init value by one order of magnitude also enhances the S-abundances by one order of magnitude during the entire duration of core evolution. For the particular case of HCS+and CS, the main physical parameters significantly influencing their abundances are the density (the lower density, the higher the HCS+and CS abundances for t>104yr) and the cosmic-ray ionisation rate (the higher the ζ, the higher the HCS+and CS abundances for t>5×105yr). By contrast, other physical parameters, such as the gas temperature, barely affect the HCS+and CS abundances for t<106yr. 5.4. Comparison with observations We now compare the different models considered in the previous section with observational results to derive the current evolutionary stage of B335. Figures A.22–A.25 show the ratio between the abundances obtained from the models (Xmodel) and from the observations (Xobs) for several sulphur-bearing species. In particular, we have considered those S-species for which we observe more than one emission line and those that are not close to the limit of detection. The results are shown for five specific evolution times (t= 104, 5 ×104, 105, 5 ×105, and 106yr). In order to find the cases with the closest values between the theoretical and observational results, we have considered for each model the largest number of S-species for which −0.5 ≤log(Xmodel/Xobs)≤0.5 in order to only consider discrepancies between model and observational results smaller than a factor of ∼3. Results show that the evolutionary stages that better reproduce the observed abundances of the largest number of S-species in B335 are those between t= 104–105yr. This result obtained from a chemical approach agrees with the results obtained from Evans et al. (2015) who derived an age of ∼5×104yr for B 335 through a dynamical perspective, by applying an inside-out collapse model that is compared with ALMA observations. Altogether this reveals, therefore, the particularly early evolutionary stage of B 335, which is comparable to the ages of pre-stellar condensations (e.g. Caselli & Ceccarelli 2012). As B 335 is an isolated globule, this may be due to the fact that the condensation has formed from very little dense gas (∼103cm−3) without a dense pre-phase between the formation of the globule and the beginning of the collapse. This gives B 335 a different chemistry compared to other young protostars that have form in dense molecular clouds, such as L483 and B1-b. A199, page 9 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) Fig. A.1. Observed lines of CS in B 335 (black) with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.2. Observed lines of 13CS in B 335 (black) with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.3. Observed lines of C34S in B 335 (black) with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. G. Esplugues et al.: Extremely young sulphur chemistry in the isolated Class 0 object B 335 Fig. A.1. Observed lines of CS in B 335 (black) with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.2. Observed lines of 13CS in B 335 (black) with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.3. Observed lines of C34S in B 335 (black) with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.4. Observed lines of C33S in B 335 (black) with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Article number, page 17 of 22 Fig. A.4. Observed lines of C33S in B 335 (black) with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. A199, page 16 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) Fig. A.5. Observed lines of CCS in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.6. Observed lines of CC34S in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. A199, page 17 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) Fig. A.7. Observed lines of CCCS in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.8. Observed lines of CCC34S in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.9. Observed lines of OCS in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.10. Observed lines of HCS+in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. A199, page 18 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) Fig. A.11. Observed lines of HC34S+in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.12. Observed lines of SO in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.13. Observed lines of 34SO in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.14. Observed lines of o-H2CS in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. A199, page 19 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) Fig. A.15. Observed lines of p-H2CS in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.16. Observed lines of NS in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.17. Observed line of SO2in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.18. Observed line of H2S in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.19. Observed line of H234S in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. A199, page 20 of 21 Esplugues, G., et al.: A&A, 678, A199 (2023) Fig. A.20. Observed line of HSCN in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. Fig. A.21. Observed line of HNCS in B 335 with rms>3σ. The red dashed line indicates the systemic velocity of the source vLSR=8.3 km s−1. CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=104 yr ζ =1.3x10−17 s−1 ζ =1.3x10−16 s−1 CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=5x104 yr ζ =1.3x10−17 s−1 ζ =1.3x10−16 s−1 CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=105 yr ζ =1.3x10−17 s−1 ζ =1.3x10−16 s−1 CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=5x105 yr ζ =1.3x10−17 s−1 ζ =1.3x10−16 s−1 CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=106 yr ζ =1.3x10−17 s−1 ζ =1.3x10−16 s−1 Fig. A.22. Ratios between the abundances obtained from the models, Xmodel, and the observations, Xobs, for different sulphur species at specific evolution times (t=104, 5×104, 105, 5×105, and 106yr), when varying the cosmic-ray rate (ζ). CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=104 yr T =7 K T =15 K CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=5x104 yr T =7 K T =15 K CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=105 yr T =7 K T =15 K CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=5x105 yr T =7 K T =15 K CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=106 yr T =7 K T =15 K Fig. A.23. Ratios between the abundances obtained from the models, Xmodel, and the observations, Xobs, for different sulphur species at specific evolution times (t=104, 5×104, 105, 5×105, and 106yr), when varying the gas temperature. CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=104 yr n =2x104 cm−3 n =2x105 cm−3 CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=5x104 yr n =2x104 cm−3 n =2x105 cm−3 CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=105 yr n =2x104 cm−3 n =2x105 cm−3 CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=5x105 yr n =2x104 cm−3 n =2x105 cm−3 CS CCS CCCS H2CSOCS SO HCS+ −2 −1 0 1 2 log( X model/ X obs) t=106 yr n =2x104 cm−3 n =2x105 cm−3 Fig. A.24. Ratios between the abundances obtained from the models, Xmodel, and the observations, Xobs, for different sulphur species at specific evolution times (t=104, 5×104, 105, 5×105, and 106yr), when varying the hydrogen number density (nH). CS CCS CCCS H2CSOCS SO HCS+ −3 −2 −1 0 1 2 3 log( X model/ X obs) t=104 yr S +=1.5x10−6 S +=1.5x10−5 CS CCS CCCS H2CSOCS SO HCS+ −3 −2 −1 0 1 2 3 log( X model/ X obs) t=5x104 yr S +=1.5x10−6 S +=1.5x10−5 CS CCS CCCS H2CSOCS SO HCS+ −3 −2 −1 0 1 2 3 log( X model/ X obs) t=105 yr S +=1.5x10−6 S +=1.5x10−5 CS CCS CCCS H2CSOCS SO HCS+ −3 −2 −1 0 1 2 3 log( X model/ X obs) t=5x105 yr S +=1.5x10−6 S +=1.5x10−5 CS CCS CCCS H2CSOCS SO HCS+ −3 −2 −1 0 1 2 3 log( X model/ X obs) t=106 yr S +=1.5x10−6 S +=1.5x10−5 Fig. A.25. Ratios between the abundances obtained from the models, Xmodel, and the observations, Xobs, for different sulphur species at specific evolution times (t=104, 5×104, 105, 5×105, and 106yr), when varying the initial sulphur abundance (S+ init). A199, page 21 of 21