S Feng

C Codella

C Ceccarelli

P Caselli

A Lopez-Sepulcre

R Neri

F Fontani

L Podio

B Lefloch

H B Liu

R Bachiller

S Viti

Typeset using

Keywords: Stars: formation, Stars: low-mass, ISM: lines and bands, Submillimeter: ISM

We present observations of SO and SO 2 lines toward the shocked regions along the L1157 chemically rich outflow, taken in the context of the Seeds Of Life In Space IRAM-NOrthern Extended Millimeter Array Large Program, and supported by data from Submillimeter Array and IRAM-30 m telescope at 1.1-3.6 mm wavelengths. We simultaneously analyze, for the first time, all of the brightest shocks in the blueshifted lobe, namely, B0, B1, and B2. We found the following. (1) SO and SO 2 may trace different gas, given that the large(-scale) velocity gradient analysis indicates for SO 2 a volume density (10 5 -10 6 cm -3 ) denser than that of the gas emitting in SO by a factor up to an order of magnitude. (2) Investigating the 0.1 pc scale field of view, we note a tentative gradient along the path of the precessing jet. More specifically, χ(SO/SO 2 ) decreases from the B0-B1 shocks to the older B2. (3) At a linear resolution of 500-1400 au, a tentative spatial displacement between the two emitting molecules is detected, with the SO peak closer (with respect to SO 2 ) to the position where the recent jet is impinging on the B1 cavity wall. Our astrochemical modeling shows that the SO and SO 2 abundances evolve on timescales less than about 1000 years. Furthermore, the modeling requires high abundances (2 × 10 -6 ) of both H 2 S/H and S/H injected in the gas phase due to the shock occurrence, so pre-frozen OCS only is not enough to reproduce our new observations.

INTRODUCTION

In the earliest stage of star formation, supersonic jets generated by protostars impact on the natal dense parent cloud and create shocks (e.g., [START_REF] Frank | Protostars and Planets VI[END_REF]. Shocks compress and heat the surrounding gas, driving complex chemistry in a short timescale, including endothermic chemical reactions, ice mantle sublimation, and sputtering. These chemical pro-the evolutionary stages of star-formation (i.e., constraining the age of the outflow and associated shocks; Charnley et al. 1997;Hatchell et al. 1998;Buckle & Fuller 2003;Wakelam et al. 2004a). For example, in agreement with theoretic predictions, observations revealed that the SO 2 /SO ratio is enhanced in young (< 10 3 yr) shocked region compared with the cloud, due to gaseous endothermic reactions; while this ratio decreases for the conversion of SO 2 to atomic sulfur in the old shocks (> 3 × 10 4 yr, e.g., Pineau des Forets et al. 1993;Hatchell et al. 1998;Wakelam et al. 2004a;Feng et al. 2015). However, the feasibility of taking the relative abundances of S-bearing species as a star formation dating tool depends on the physical condition of the source and a proper assumption of the S grain reservoir.

Moreover, understanding of the S-bearing chemistry is challenged by a collection of mysteries, the major one being the "S-depletion problem" (e.g, Tieftrunk et al. 1994); i.e., although S is known to be significantly depleted in molecular dark clouds, its depleted form (the main reservoir on grain mantles) is far from being established. Several species have been proposed to be the main carriers of S on dust grains, such as H 2 S (e.g., Li et al. 2015;Holdship et al. 2016) and OCS (e.g., Hatchell et al. 1998;Wakelam et al. 2004a;Liu et al. 2012;Podio et al. 2014;Minh et al. 2016;Holdship et al. 2019). Historically, H 2 S has never been definitively detected on ices, which calls for further and more comprehensive S-chemistry study using telescopes at higher sensitivity. Recent observations and models have shown that, at the very least, H 2 S is not the only S-bearing species on grains (e.g., Wakelam et al. 2004b;Codella et al. 2005). By far, OCS is detected (e.g., [START_REF] Geballe | [END_REF]Palumbo et al. 1995;Podio et al. 2014;Boogert et al. 2015), and SO 2 (e.g., Boogert et al. 1997) is tentatively detected in interstellar ices. Therefore, the main reservoir of S on dust is still unclear.

The blueshifted outflow lobe of L1157

Regions containing successive young ( 10 3 yr) shocks driven by mass losses from protostars are ideal laboratories for studying time-dependent S-chemistry, as the short timescales of the shocks provide the possibility to either directly observe the material coming off grains, or to trace the chemistry back to the original content of interstellar ices.

This well applies to the target of the present study, the shocks associated with the L1157 outflow region. The Class 0 protostar L1157-mm ( 3 L , Tobin et al. 2010), located at d ∼352 pc (Zucker et al. 2019), drives an episodic and precessing jet (e.g. Gueth et al. 1996Gueth et al. , 1998;;Podio et al. 2016), creating a bipolar outflow (e.g., Bachiller et al. 2001;Nisini et al. 2010). The impacts between the jet and the highdensity (10 4 -10 5 cm -3 , Lefloch et al. 2012;Benedettini et al. 2013;Gómez-Ruiz et al. 2015) cavity walls triggered several shocks in the last ∼ 2 × 10 3 yr (e.g., Lefloch et al. 2010).

The targets of the present work are the B0, B1, and B2 shocked regions, which are the brightest shocks associated with the southern blueshifted outflow lobe, located within 0.1 pc from the protostar in the plane of the sky. Specifically, the B1 region has been studied in detail by using either single-dish telescopes (e.g., Herschel-CHESS, Ceccarelli et al. 2010; 30 m-ASAI, Lefloch et al. 2018) or interferometers (SMA, Very Large Array, NOEMA, Tafalla & Bachiller 1995;Codella et al. 2009;Lefloch et al. 2012;Benedettini et al. 2013), revealing a clumpy bow-like structure associated with several chemical processes, such as the sputtering of dust mantles, the disruption of the dust refractory cores, and warm gas-chemistry reactions. The abundance of molecular species, from diatomic-such as CO, CS-to molecular ions-such as HOCO + , SO + , and HCS + Podio et al. (2014)and more complex molecules such as CH 3 OH, CH 3 CHO, and NH 2 CHO (e.g., Bachiller et al. 2001;Benedettini et al. 2007;Arce et al. 2008;Codella et al. 2010Codella et al. , 2015Codella et al. , 2017;;Busquet et al. 2014;Fontani et al. 2014a;Lefloch et al. 2017Lefloch et al. , 2018) ) are consequently enhanced. Specifically, S-bearing species (e.g., CS, SO, SO 2 , H 2 S, OCS) have been detected toward this region (see, e.g., Bachiller et al. 2001;Benedettini et al. 2013;Podio et al. 2014;Holdship et al. 2016).

Kinematically, B2 is older than B1 (e.g., Bachiller et al. 2001;Podio et al. 2016); B1 mainly consists of three shock layers produced by episodes of the precessing jet, with the oldest layer located the farthest from the protostar, i.e. in the southernmost edge of B1 (Codella et al. 2015). Observing the jet radial velocity in the inner knots, the precessing model (Podio et al. 2016) indicates that 1 B0 is ∼1340 yr old, B2 is ∼2530 yr old, and the kinematic age of B1 knots ranges from 1550 yr at the north edge to 1760 yr at the southern edge. Therefore, the entire B0-B1-B2 region provides us one of the best laboratories to study the time-dependent S-chemistry.

In this paper, we present the data in Section 2, which are obtained from the "Seeds Of Life In Space"(SOLIS2 ; Ceccarelli et al. 2017) large program and supplementary Submillimeter Array (SMA) observations. Our results are presented in Section 3 and analyzed in Section 4. We study the map of relative abundance ratio between SO and SO 2 in Section 5.1, and discuss the possible origin of SO and SO 2 by using our chemical model to fit the observation measurement in Section 5.2. We conclude in Section 6.

Using the NOrthern Extended Millimeter Array (NOEMA) and the SMA, we carried out a line-imaging survey at 1.1-3.6 mm toward our targets. Given that the antennae of NOEMA and SMA have different diameters, and that the molecular lines we observed are at different wavelengths, the primary beam attenuation is different. Therefore, the B0, B1, and B2 shocked regions have been imaged by using one point or two mosaic points with NOEMA or SMA at different wavelengths.

NOEMA

Using NOEMA, we observed B0-B1 at 3.6 mm, 3.1 mm, and 1.5 mm in its A (7 antennae), C (6-8 antennae), and D (5-8 antennae) configurations from 2015 July to 2017 January. With the shortest projected baseline of 24 m, the observations are sensitive to structures up to 19 (at 3.6 mm), 16 (at 3.1 mm), and 8 (at 1.5 mm), respectively. For all observations, a common phase center at 20 h 39 m 10 s .20, 68 • 01 10 .50 (J2000) and the systemic velocity V sys ∼ 2.7 km s -1 was used. The precipitable water vapor (PWV) varied between 1 and 40 mm during the observations. Standard interferometric calibrations were performed during the observations. The bandpass, gain (phase/amplitude), and flux calibrators used during the observations are recorded in Table A1. The uncertainty of absolute flux scale is estimated to be correct to within ∼ 5%.

With the wide-band receiver (WIDEX), three spectral setups cover the frequency ranges of 80. 790-84.407, 95.794-99.499, and 203.962-207.632 GHz, at a frequency resolution of 1.95 MHz. We detected one SO 2 line and three SO ( 34 SO) lines (with E u /k B < 40 K). From our ASAI-30 m single-point observations at a frequency resolution of 0.195 MHz (Lefloch et al. 2018), these lines have full-width half-maximum (FWHM) linewidth of > 8 km s -1 toward the B1 shocked region (Mendoza et al. 2014), so the velocity resolution from WIDEX observations (∼ 7 km s -1 at 83.688 GHz, see Table 1 and Figure 1) is sufficient for our chemical analysis. Moreover, we use the narrow-band correlator and set four windows in each spectral setup (80 MHz backends with a spectral resolution of 156 kHz) to target on several complex organics lines (see the spectral setups in Ceccarelli et al. 2017;Codella et al. 2017, for details). Other lines detected toward this region will be presented in the follow-up papers (e.g., CS/CCS/OCS/NS/SiS/H 2 S/H 2 CS, cyanopolyynes; S. Feng et al. in prep.).

Furthermore, aiming at a robust analysis on the possible chemical differentiations between SO 2 and SO toward different time-dependent shock layers of B1, we also include one 34 SO and four SO 2 ( 34 SO 2 ) lines from our previous NOEMA-WIDEX observations (obsID: W068, Fontani et al. 2014a;obsID: X058, Podio et al. 2017; see observation details therein). The spectroscopic parameters of all these lines are taken from the Cologne Database for Molecular Spectroscopy (CDMS, Müller et al. 2002Müller et al. , 2005) ) and listed in Table 1.

Calibration and imaging were performed by using the standard procedures with the CLIC and MAPPING softwares of the GILDAS3 package. After testing different weighting algorithms, we used the natural weighting to achieve a better signal-to-noise ratio (S/N). However, our NOEMA data recover only 40-60% of the emission when compared with our ASAI-30 m pointing observations, indicating that these SO 2 and SO lines should have extended emissions toward our target.

Averaging the line-free channels of the WIDEX band, we did not detect significant (> 5σ rms) continuum emission at 3.6 mm with a synthesized beam of 5.08 × 4.69 (σ=0.016 mJy beam -1 ), at 3.1 mm with a synthesized beam of 3.12 × 2.30 (σ=0.012 mJy beam -1 ), or 1.5 mm with a synthesized beam of 1.20 × 0.85 (σ=0.036 mJy beam -1 ).

From our NOEMA observations, the detected SO isotopologue lines have low-J levels (E u /k B =9-18 K), while the SO 2 isotopologue lines have mid-J or high-J levels (E u /k B =16-102 K). They may selectively trace gas with different temperatures and therefore show different spatial distributions. With the aim of investigating any chemical difference from the spatial distribution maps of both species, it is necessary to observe more SO 2 and SO lines, to cover the same range of upper energy level (E u /k B ).

SMA

To minimize the excitation effect on the spatial distribution differentiation between SO and SO 2 , we carried out four-track observations toward our source with the SMA at 230/255 GHz (1.14-1.40 mm), by using the compact configuration. From 2017 June to July, we observed B1 and B2 with a two-point mosaic per track. The phase centers of B1 and B2 are at 20 h 39 m 09 s .500, 68 • 01 10 .500 and 20 h 39 m 11 s .358, 68 • 00 51 .00 (J2000), respectively. For all of the observations, the phase and amplitude calibrations were performed via frequent switch (every 20 minutes) on quasars 1849+670 and 1927+739.

The baselines range from 16 m to 76 m with 7 or 8 antennae at different dates. B1 and B2 share the bandpass (3C273 or 3C454.3) and flux (Callisto and MWC349a) calibrators. The zenith opacities, measured with water vapor monitors mounted on the James Clerk Maxwell Telescope (JCMT), were satisfactory during all tracks with τ(225 GHz) ∼ 0.05-0.3. The system temperature was 100-200 K at 230 GHz and 200-400 K at 255 GHz. Further technical descriptions of the SMA and its calibration schemes can be found in Ho et al. (2004).

We employed the dual-RX, single polarization and double-sideband observing mode. Our frequency coverages were sampled using the SWARM correlator with a 0.140 MHz frequency resolution (∼0.182 km s -1 at 230.538 GHz). We tuned 230.0 GHz and 255.8 GHz to be at the center of the spectral chunk 1 in the upper sideband of the receivers Rx230 and Rx240, respectively. To increase the S/N ratio per channel, we smooth the frequency resolution to 4.469 MHz (a velocity resolution of 5.812 km s -1 at 230.538 GHz).

The flagging and calibration was done with the MIR package (Scoville et al. 1993;[START_REF] Qi | SFChem 2002: Chemistry as a Diagnostic of Star Formation[END_REF] 4 . The mosaic imaging was conducted with MIRIAD package by using the natural weighting [START_REF] Sault | Astronomical Data Analysis Software and Systems IV[END_REF].

We detected four SO lines in the entire frequency coverage (E u /k B =35-56 K), while no SO 2 line is detected with > 3σ emission (σ=42.7 mJy beam -1 at an angular resolution of ∼ 3.4 ). Their spectroscopic parameters taken from CDMS are listed in Table 1 as well. In addition, we also detected lines from SiO, HNCO, CO, CS, H 2 CO in these observations.

Averaging the line-free channels, the sensitivity of continuum emission at 1.1-1.4 mm achieves σ=0.17 mJy beam -1 with a synthesized beam of 3.37 ×3.13 . At such sensitivity, we have marginally detected 3-4 σ of continuum emission toward B1-B2 region (see H. B. Liu in prep. for a detailed continuum study).

IRAM-30 m Telescope

To compensate the missing flux, we use IRAM-30 m telescope and performed the observations in the on-the-fly (OTF) mode from 2017 May to November. We mapped a 2.5 × 2.5 area (at least twice that of the primary beam of the interferometers for successful combination, according to Pety et al. 2013, covering both B1 and B2), centered at 20 h 39 m 09 s .50, 68 • 01 10 .5 (J2000). By superpositioning different spectral tunings, the broad bandpass of EMIR (16 GHz bandwidth simultaneously for each spectral tuning) covers the frequency range of 79. 302-87.084, 91.821-102.764, 106.902-114.683, 130.703-138.485, 161.102-169.882, 203.925-211.703, 213.422-227.385, 229.105-236.882, 239.420-247.204, and 255.103-262.880 GHz, with a frequency resolution of 0.195 MHz (by using the FTS200 backend). The FWHM beam of the 30 m telescope at 83.688 GHz is ∼ 30.86 .

To combine the emissions of SO 2 and SO isotopologue lines obtained from NOEMA or SMA at high spatial resolution with their large-scale emission data obtained from IRAM-30 m, we aligned the phase center and smoothed the velocity resolution of lines obtained from IRAM-30 m to those from NOEMA or SMA. Afterwards, we use the task "UV-SHORT" in the MAPPING software for NOEMA-30 m combination; and use the MIRIAD package for SMA-30 m combination. The synthesized beam and the 1σ rms value of the data cube of each line, after NOEMA-30 m or SMA-30 m combination, are listed in Table 1.

MOLECULAR SPATIAL DISTRIBUTION

For the first time, the interferometric data of L1157 B0-B1-B2 are complemented with the IRAM-30 m mapping, providing us images with full spatial scale coverage, as well as high angular resolutions. Moreover, the wide frequency coverage allows us to study five SO 2 lines (including one 34 SO 2 line) and eight SO lines (including two 34 SO line), with E u /k B in the range of 16-102 K. At a linear resolution of < 1000 au, we will be able to search for chemical differentiations from the systematic differentiations in the spatial distribution of both species.

All lines show a significant blueshifted emission toward B0-B1-B2, and the S/N > 3 part of each line is in the velocity range of -10 to +5 km s -1 (Figure 1 andA1). Specifically, the peaks of all SO 2 lines are close to the V sys , when observing at a velocity resolution of 3.6-7.0 km s -1 . Observed at higher velocity resolutions (2.8-5.9 km s -1 ), some SO lines seem to peak at ∼ -2 km s -1 toward positions such as B0e, B1f, and B1a. However, such a velocity offset (∼ 5 km s -1 ) only corresponds to 1-2 channel(s), and it smears out when these SO lines are binned to a coarser velocity resolution (7.0 km s -1 ). Therefore, whether these SO and SO 2 lines trace different gas calls for further observations at a highervelocity resolution.

We imaged the spatial distribution of each line by integrating their emissions over the velocity range5 of -10 to +5 km s -1 . Specifically, Figure 2 show the spatial distributions of SO 2 (8 1,7 -8 0,8 ) and four SO lines, i.e., 6 5 -5 4 , 5 5 -4 4 , 7 6 -6 5 , and 6 6 -5 5 , which are sensitive to a comparably larger field than the rest of the lines (see their primary beams in Table 1). Overlaid with CO (1-0) emission contours (Gueth et al. 1996), we found that all of these SO and SO 2 line emissions trace the extended structure from the eastern cavity wall B0 to two cavities associated with B1 and B2 (denoted as C2 and C1, respectively in Gueth et al. 1998). These lines show the same shock knots in B0 and B1 as those in SiO (5-4) from the same SMA-30 m observations6 (Figure 2) as well as SiO (2-1) observed by NOEMA-30 m (Gueth et al. 1998), which is consistent with the prediction of a precessing jet model (Podio et al. 2016).

Aside from the above five lines, the rest SO and SO 2 isotopologue lines in this work are only sensitive to the B0-B1 region. Nevertheless, they offer images observed at higher angular resolutions. At a spatial resolution of 500-1400 au, six clumpy substructures appear toward the B0-B1 region in Figure 3. These substructures are already identified from previous observations (e.g., by using line emissions from CH 3 CN, HCN, CH 3 OH, SiO, CS, H 2 CO, HCO + , and HNCO in Benedettini et al. 2007;Codella et al. 2009;Gómez-Ruiz et al. 2013;Burkhardt et al. 2016). In this work, we use the consistent names of B0e, B1a, B1b, B1c, and B1f with those defined by Benedettini et al. (2007); Codella et al. (2009); Burkhardt et al. (2016). Moreover, we denote a position B1i, which was noted previously but not labeled as a clumpy substructure from emission of other species. For any particular substructure traced by more than one molecule, the difference in the measurement of its absolute coordinate by different molecular tracers is about 1 -2 . Based on Figure 3, the absolute coordinates of these clumpy substructures are listed in Table 2.

Morphologically, B0e outlines the eastern wall of the northernmost blueshifted bow-shock cavity (e.g., named missed toward B2 in previous SMA-only observations (Gómez-Ruiz et al. 2013).

as "C2" in Gueth et al. 1998), showing higher abundance of species such as SiO and HCN than the western wall and the connecting ridge B1f (Gómez-Ruiz et al. 2013) . B0e and B1a are suggested as knots where episodic ejection impacting against the cavity wall (Podio et al. 2016), and high-velocity SiO and CS emissions (Benedettini et al. 2007;Gómez-Ruiz et al. 2013) are detected toward them. A U-shape structure connecting B1a-B1c-B1b pins down the shock front, and the apex B1i is coincident with the "finger" structure detected on the maps of CS and low-velocity SiO, which is suggested as the magnetic precursor of the shock (Gueth et al. 1998).

According to the spatial distributions of detected complex organic molecules (COMs), several chemical layers are defined in B1. One chemical layer is composed of B0e-B1a-B1f-B1b, where the detected COMs include HDCO (Fontani et al. 2014b), CH 3 CHO (Codella et al. 2017), and CH 3 OH (Codella et al. 2017), indicating the place where the current dust mantle and core are released by sputtering. To the southern shock front, B1c and B1i denote another chemical layer where NH 2 CHO is the main COM emitter (Codella et al. 2017), and this layer contains the gas older than B1a-B1f-B1b (Gueth et al. 1998;Podio et al. 2016).

To compare the line emissions along the path of the processing jet, we also denote a position B2a, where SO and SO 2 lines show the emission local maximum within B2.

All of the SO and SO 2 lines show the strongest emission toward B1a. More specifically, 34 SO, SO 2 , and 34 SO 2 lines show the U-shape structure toward the wall and southern rim 2). All lines are extracted from images smoothed to the same pixel size and the same angular resolution (5.1 ). The line profile of SO (6 5 -5 4 ) shown in black or dashed blue has a velocity resolution of 6.09 km s -1 (original spectral resolution) or 7.01 km s -1 (the same as that of SO 2 8 1,7 -8 0,8 ), respectively. In each panel, two red vertical lines at -10.0 km s -1 and +5.0 km s -1 indicate the velocity range for which we integrate the intensity; the gray dashed vertical line indicate the V sys = 2.7 km s -1 of the cloud. The horizontal cyan line indicates the baseline (T B =0 K). of B1, linking B0e-B1a-B1c-B1b. In contrast, SO lines suggest an anticorrelated spatial distributions, having a bright emission toward the eastern B1a. This finding could either reflect a chemical segregation or be due to optical depth effects, and it will be investigated in the next sections.

ANALYSIS

To quantify the possible chemical differentiations between the SO and SO 2 molecules by excluding excitation effects, we need to constrain their column densities over the entire targeted region. Uncertainty in the molecular column density measurement is mainly from two parts, whether the optical depth of each line can or has to be corrected, and whether these lines can be assumed in the local thermal equilibrium (LTE) condition.

Optical depths

In the shocked region, we cannot exclude the possibility that some low-J SO and SO 2 lines may be optically thick. We note that 34 SO (3 2 -2 1 ) and 32 SO (3 2 -2 1 ) are a pair of lines that have similar upper-level energies (E u /k B ∼ 9 K) and Einstein coefficients (A i j ∼ 10 -6 s -1 ). Therefore, their brightness temperature ratio, T B, 34 SO(3 2 -2 1 ) /T B, 32 SO(3 2 -2 1 ) , can be used to measure the optical depth τ of 32 SO (3 2 -2 1 ) line (Myers et al. 1983) 7 .

Figure 4 shows that such ratio is in the range of 0.1-0.15 toward the pixels where both lines show > 5σ emissions. Considering the 10% systematic uncertainty per pixel, T B, 32 SO(3 2 -2 1 ) /T B, 34 SO(3 2 -2 1 ) , toward pixels within B0-B1 can be treated as uniform. Therefore, if we assume that 34 SO(3 2 -2 1 ) is optically thin, we can reasonably assume that 32 SO(3 2 -2 1 ) has uniform optical depths toward all valid pixels over the entire B0-B1 region.

Quantitatively, using the standard isotopic ratio R32 S/ 34 S ∼ 22 in the solar system (e.g., Wilson & Rood 1994;Chin et al. 1996;Lodders 2003), τ of the 32 SO (3 2 -2 1 ) line in B0-B1 can be estimated. Solving numerically the Myers et al. (1983) approximation (see footnote 7), we present the range of τ on Figure 4 as well. In general, τ of the 32 SO (3 2 -2 1 ) line is not negligible (≤ 1.5 toward B0e and B1a, as well as ≤ 3 toward B1f, B1c, and B1i, and ≤ 4 toward the western edge), but it has small dynamic range in B0-B1 region. As the absolute value of τ is an approximation and depends on the assumption of R32 S/ 34 S , τ toward individual pixels is not im-7 Assuming that fractionation (isotopic exchange reaction) of 34 S ⇔ 32 S is stable, according to Myers et al. (1983), we can estimate the optical depth at the line center of a given line τ α,0 by measuring the ratio between the main beam brightness temperature of the main line T mb, α,0 and its isotopologue T mb, β,0 as 1-exp(-τ α,0 /Rα)

1-exp(-τ α,0 ) ≈ T mb, β,0 T mb, α,0
. Here, R α is the intrinsic abundance of the main isotope (e.g., 32 S) compared to its rare isotope (e.g., 34 S) in the ISM. A numerical solver to this approximation indicates that τ α,0 is appropriate to ln( T mb, β,0 T mb, α,0

) at each pixel. All of the maps are obtained from the combination of NOEMA (or SMA) and IRAM-30 m, and line emissions with < 5σ rms are blanked. The red contours, starting from 8σ (σ = 0.27 Jy beam -1 ) and increasing with the step of 8σ, show the CO (1-0) emission, obtained from NOEMA and IRAM-30 m combination (Gueth et al. 1996). Gray circle(s) in each panel indicates the primary beam of the NOEMA/SMA at corresponding line rest frequency, and the angular resolution after IRAM-30 m combination is plotted at the bottom left in each panel. The black curve guidelines the path of the precessing jet from L1157-mm as modeled by Podio et al. (2016). The black squares indicate the driving source, L1157-mm, and the shock knots in B0, B1, and B2.

portant in this work. Instead, the relative gradient of τ from B0e to B1a and B1i is ∼1, indicating that the relative gradient of SO column density in B0-B1 does not change significantly before and after line optical depth correction.

We are not able to estimate the τ of any 32 SO 2 line, due to the lack of detections on their 34 SO 2 line counterparts (with the same transition level). Nevertheless, the relative abundance ratio between this isotopologue pair is comparable with the standard Solar R32 S/ 34 S (see Section 4.2), which indicates that the optical depths of the SO 2 lines may not lead to an underestimation of the SO 2 column density in this work.

Rotation temperature and molecular column density from Rotational Diagram (RD)

Although the line emissions are clumpy, their spatial distributions toward the positions list in Table 2 are more extended than their synthesized beams, so we do not take the beam filling factor into account in the following estimates. As the first-order approximation, we assume that all of the detected SO and SO 2 lines are optically thin and under LTE. Smoothing the line-intensity maps to the same angular resolution (i.e., 5.30 ×4.86 ) and the same pixel size (i.e., 0.5 ), we use the RD method, and fit the relative population distribution of all the detected energy levels by assuming a single component toward each clumpy substructure (see RD fit toward B1a in Figure 5 as an example). The least-square solutions of the molecular rotational temperature T rot , as well as total column density N rot from fittings are list Table 3.

When taking the optical depth of the SO (3 2 -2 1 ) lines into account, the SO column density in pixels within B0-B1 increases by a factor of 1.5-2.5 (Table 3). Taking into account the uncertainties, it seems that the difference between the SO column density throughout the B0-B1 structure is little, ranging from (3.9 ± 1.7) × 10 14 cm -2 (B1a) to (1.5 ± 1.0) × 10 14 cm -2 ( B1i). Also, the SO 2 column densities, as derived from RD, looks quite uniform, being around 1 × 10 14 cm -2 . SO seems to have uniform T rot toward the B0-B1 region (10-14 K), and the optical depth correction for the SO (3 2 -2 1 ) line does not bring in significant change in the T rot estimates. In contrast, precise T rot estimation of SO 2 depends on whether the high-J line (with E u /k B >100 K) is included, i.e., fitting to the mid-J SO 2 lines yields the same T rot as that from fitting the low-J and mid-J SO lines (with E u /k B <70 K). A possible reason could be that the high-J SO 2 may trace a warmer component.

The upper energy level of the detected 34 SO lines are too close to derive the molecular rotational temperature map, and we are not able to obtain the rotational temperature of 34 SO 2 from only one detected line (7 1,7 -6 0,6 ). With the LTE and optically thin assumptions, we estimate the column densities of 34 SO and 34 SO 2 by using the T rot of their 32 S-isotopologues toward the clumpy substructures in B0-B1. The relative abundance ratio between each pair of 32 S-and 34 S-isotopologues is uniform (with a small standard deviation among pixels), especially toward B1. The ratio of χ( 32 SO/ 34 SO) is 9 ± 2 over the entire B0-B1 (∼ 7 toward B0e and 15 ± 4 toward B1 after SO (3 2 -2 1 ) optical depth correction), while the ratio of χ( 32 SO 2 / 34 SO 2 ) is ∼ 6 toward B0e and ∼ 15 toward B1a-B1c-B1b (assuming τ ≤ 1 of all of the SO 2 lines detected in this work). These ratios are consistent with the χ(C 32 S/C 34 S) ∼ 15 measured toward B1 from our previous single-dish observations (Bachiller & Pérez Gutiérrez 1997) 8 , and is comparable with the R32 S/ 34 S ∼ 22 in the solar system. Therefore, the optically thin assumption for the SO 2 lines can be used to measure the SO 2 column density in this work.

Kinetic temperature and molecular column density from Large Velocity Gradient (LVG)

In the temperature range of 10-100 K, the critical densitie of all of the detected SO and SO 2 isotopologue lines are in the range of 10 5 -10 6 cm -3 (estimated by using the Einstein coefficient and the collision rate given in Leiden Atomic and Molecular Database, Schöier et al. 2005, which are consistent with the effective critical density given by, e.g., [START_REF] Reiter | [END_REF]. According to previous studies, such as Benedettini et al. (2013);Gómez-Ruiz et al. (2015), the number density toward the entire B1 region is 10 5 -10 6 cm -3 . Therefore, excitation condition of SO and SO 2 gas in this shocked region need a LVG consideration.

Using the collision coefficients of SO-(para) H 2 and SO 2 -(para and ortho) H 2 obtained from the BASECOL database Map of the integrated intensity ratio T B, 34 SO (3 2 -2 1 ) /T B, 32 SO (3 2 -2 1 ) (top widget) as well as the optical depth of 32 SO (3 2 -2 1 ) (right widget) toward each pixel. The gray dashed contours, starting from 8σ (σ = 0.27 Jy beam -1 ) and increasing with the step of 8σ, show the CO (1-0) emission (Gueth et al. 1996). The pixels where 34 SO (3 2 -2 1 ) shows < 6σ emissions are blanked. Six clumpy structures toward B1 are labeled. (Dubernet et al. 2013), and applying the code of LVG approximation to non-LTE analysis described in Ceccarelli et al. (2002), we run grids of parameters and find the leastsquare (χ 2 ) best-fitting solution toward four positions along the precessing path, i.e., B0e, B1a (as well as B1b as a double-check), and B1i. B2a is not fit because less than four lines show > 5σ emission. In the fittings, we adopted a FWHM linewidth of 7 km s -1 toward each position. Then, we vary three parameters in a relatively large range: H 2 number density n (10 4 cm -3 -10 8 cm -3 ), gas kinetic temperature T kin (30 K-300 K), molecular column density N T (10 13 cm -2 -10 18 cm -2 ). Note that the FWHM linewidth we adopt in this work is a compromise made to obtain the best fit by including all SO and SO 2 lines observed at a coarse velocity resolution (6-7 km s -1 ). According to the previous analysis based on CO observations at a velocity resolution of 0.1-0.5 km s -1 (Lefloch et al. 2012), at least three components overlap in the velocity range from -30 to 5 km s -1 , with their intensity peaks close to the V sys . Before we resolve them with observations at higher-velocity resolution, the adopted FWHM linewidth in the range of 5-10 km s -1 provides results that represent the average behavior of the different kinematical components within this velocity range.

Comparing the best-fit results for the column density by using LVG and RD methods (Table 3), we have a good agree-ment. In addition, we find that the SO column density toward B1i is slightly lower than those measured toward the rest of the clumps, even after taking into account the uncertainties (see Table 3). The kinetic temperatures are in the range of 30-200 K, being substantially in agreement with the kinetic temperatures found by Lefloch et al. (2012) in the different gas components coexisting in the B1 shocked gas. On the other hand, the gas volume density n from LVG fit of SO lines is well constrained, from (3 -5) × 10 4 cm -3 (B0e) to 10 5 cm -3 ( B1i), while n as traced by SO 2 (10 5 -10 6 cm -3 ) is denser than that traced by SO by a factor up to an order of magnitude (see Table 3). The higher values of n are in agreement with those derived by Gómez-Ruiz et al. ( 2015) using a CS multiline analysis. Therefore, SO and SO 2 may be tracing different gas portions. Provided that the column densities as derived by the LVG analysis are in agreement with those derived from the LTE-RD analysis (after correction due to the optical depths), to investigate the possible chemical differentiation between SO and SO 2 along the path of the precessing jet in the entire B0-B1-B2 (in a time frame of thousand years), we fit the onecomponent RD toward pixels in B0-B1-B2 region.

When investigating the 0.1 pc scale field of view, we note that there is a tentative gradient along the path of the precessing jet in the χ(SO/SO 2 ) map before optical depth correction9 (Figure 6): χ(SO/SO 2 ) < 1 in B2, while 1.5 < χ(SO/SO 2 ) < 3 from the younger B0e. Unfortunately, only two SO 2 lines show > 5σ emission toward B2, so SO 2 column density toward B2a is only an upper limit, i.e., χ(SO/SO 2 ) is a lower limit there. Only further SO 2 observations toward B2 will allow us to verify the occurrence of such a spatial trend.

To test whether the tentative gradient of χ(SO/SO 2 ) is due to optical depth effects, we zoom into the region where SO column density can be corrected by considering the optical depth of SO (3 2 -2 1 ) line. Even after the correction for optical depth effects, Figure 7 allows us to speculate a possible trend for the χ(SO/SO 2 ) ratio decreasing moving from north (younger gas) to south (older gas). Again, the present results can only suggest the occurrence of such a spatial trend and only future multiline observations, possibly on smaller spatial scales, will verify this speculation.

--c --c --c SO 2 < 1.7 d ... ... < 64 d --c --c --c --c
Note. a. Column density estimated by taking the optical depth of SO (3 2 -2 1 ) into account. b. Rotation temperature derived by only taking the lines with E u /k B < 70 K into account. c. Less than four lines show emission with S/N > 5 in the primary beam. Upper limit is given by using two lines which show emission with S/N > 5 .

Chemical modeling

Our previous chemical study, especially the observations of the ions and COMs, unveiled at least three chemically different layers with gas-grain components from the north (B1a) to south (B1i) toward B1 (Codella et al. 2017(Codella et al. , 2020)). To trace back the possible forming path of SO and SO 2 from the χ(SO/SO 2 ) gradient through B0-B1-B2, we apply the same time-dependent gas-phase chemical model as that introduced by Codella et al. (2017). Briefly, we used the MyNahoon code, which is based on the Nahoon code originally developed by Wakelam et al. (2012). The code follows the gaseous abundances evolution in time. Note that it does not consider the reactions possibly occurring on the grain surfaces, except the formation of H 2 molecules, as our goal is to study a subsequent phase: namely, when the grain mantles are already formed, in conditions and timescales where the surface reactions are completely negligible (see later). In short, we simulate the physical effects caused by the shock passage by suddenly increasing the H 2 gas number density n from 10 4 to 10 5 cm -3 , and the temperature from 10 to 70 K.

In addition, because of the sputtering and/or grain-grain collisions caused by the shock, a number of species previously frozen on the grain mantles are injected into the gas phase. In our model, we assumed the injected species and their abundance as those in Codella et al. 2017 (their Table B.1), which is based on several previous studies of the same source (Tafalla & Bachiller 1995;Benedettini et al. 2013;Podio et al. 2014;Codella et al. 2015) or known properties of Figure 6. SO (left) and SO 2 (middle) column density maps and their relative abundance ratio map (right) derived by using the RD method, with the assumption that all lines are optically thin and under LTE. In the left panel, the pixels where less than three SO lines show > 5σ emission are blanked. In the middle panel, the pixels where less than two SO 2 lines show > 5σ emission are blanked, and the black contour denotes the region within which at least three lines show > 5σ emission. In the right panel, the pixels where less than two SO 2 lines show > 5σ emission are blanked, and the gray line indicates the path of the precessing jet (Podio et al. 2016). These maps are derived by smoothing all lines to the same angular resolution, i.e., 5.30 × 4.86 , shown as the synthesised beam in the bottom left of each panel. The gray dashed contours and the labeled clumpy structures are the same as shown in Figure 4.

the ice mantle compositions (Boogert et al. 2015;Busquet et al. 2014).

In particular, with respect to the S-bearing frozen species, Podio et al. (2014) could only roughly constrain them, based on the S-bearing ions SO + and HCS + . These authors found that OCS is injected in the gas phase with an abundance of 6 × 10 -6 (with respect to the H nuclei) but could not constrain the injected abundances of H 2 S and S, two species that are expected to be also present in the grain mantles (see e.g., Wakelam et al. 2004a;Vidal & Wakelam 2018). Therefore, we ran models with different injected abundances of H 2 S and S and tried to constrain them based on the new observations presented in this work.

The results of the modeling are shown in Figure 8. In the figure, we report the model predicted SO/SO 2 abundance ratio as a function of the time, namely the age of the shocked gas and for different values (2 × 10 -6 and 1 × 10 -8 ) of the injected H 2 S/H and S/H abundances. In the same figure, we also report the values measured in the previous sections (the uncertainty is given based on RD optically thin assumption, RD optical correction for the low-J SO line, and LVG) toward B1a, B1b and B1i and their ages estimated by Gueth et al. (1996) and Podio et al. (2016). Finally, we also show the model predictions for values slightly different of the density and cosmic ray ionization rate to have a rough "sensitivity" study and understand how much the predictions also depend on these, only roughly constrained, quantities.

The figure shows that the model predictions compare rather well with the observations and, remarkably, with the increasing ages of B1a ( B1b) to B1i. Based on the predictions, the injected H 2 S/H and S/H abundances are equal to 2 × 10 -6 : injecting pre-frozen OCS only is not enough to reproduce our new observations, an additional injection of previously frozen S-bearing species in form of H 2 S and S in approximately the same amount is also needed (a similar result was also found by Holdship et al. 2016). Second, the cosmic ray ionization rate is around 3 × 10 -16 s -1 , namely slightly (a factor 2) lower than that used in previous studies. This may mean that it slightly decreases going outward from the central object from which the outflow/jets emanate, which is consistent with the hypothesis that cosmic rays like particles are accelerated by the dense shocks at the jet base (Padovani et al. 2016;Ivlev et al. 2019).

Of course our model, although so far is successful in reproducing several species (see, for example, Codella et al. 2017Codella et al. , 2020)), is rather simplified; it does not take into account, for Figure 7. SO column density map (upper) and map of its relative abundance ratio with respect to SO 2 (lower) derived by using the RD method, with the LTE assumption and optical depth correction for the SO (3 2 -2 1 ) line. The pixels where 34 SO (3 2 -2 1 ) shows < 6σ emissions are blanked. In the right panel, the gray line indicates the path of the precessing jet (Podio et al. 2016). These maps are derived by smoothing all lines to the same angular resolution, i.e., 5.30 × 4.86 , shown as the synthesised beam in the bottom left of each panel. The gray dashed contours and the labelled clumpy structures are the same as shown in Figure 4. example, the short-lived hot/warm phase of the post-shocked gas, which, in the case of SO 2 could have an impact on the predictions. Having said that, it is encouraging to see that such a simplified model succeeds to grab the basic gaseous chemical composition changes caused by the shock.

CONCLUSIONS

We simultaneously analyze the chemical properties of SO and SO 2 toward all of the brightest shocked regions along the L1157 chemically rich blueshifted outflow, namely, B0, B1, and B2. Based on the data taken in the context of the SOLIS IRAM-NOEMA Large Program, and the supplementary observations from SMA interferometer as well as IRAM-30 m telescope at 1.1-3.6 mm wavelengths, our conclusions are as follows:

1. At an angular resolution of 1 -3 (500-1400 au at the source distance of 352 pc), SO and SO 2 may trace different gas toward B1. This is indicated from the different spatial distributions of five SO 2 lines (including one 34 SO 2 line) and eight SO lines (including two 34 SO lines) which we detected. The emission peak of the SO lines are in general toward northeast, where the gas is impacted by the recent jet impinging on the B1 cavity wall; while the emission peak of the SO 2 lines are all toward the southern edge of the older shock front. Moreover, LVG constraints on gas volume density in- The model predictions are also obtained for two values of the cosmic rays ionization rates, 6 × 10 -16 s -1 (solid curves) and 3 × 10 -16 s -1 (dashed curves), and a smaller H density, 2 × 10 4 cm -3 (dotteddashed curve). The results from the analysis of the observations are reported in the dashed areas toward three positions along the B1 region: B1a, B1b, and B1i, respectively.

dicates that, the gas traced by SO 2 has a volume density of 10 5 -10 6 cm -3 , which is denser than that traced by SO by a factor up to an order of magnitude.

2. Investigating the 0.1 pc scale field of view, there is a tentative gradient in the relative abundance ratio of χ(SO/SO 2 ) along the path of the associated precessing jet, showing smooth decrease from B0-B1a (> 2) to B1i-B2 (< 1). These findings call for further analysis at higher spatial resolutions to verify the tentative trends.

3. Our astrochemical modeling shows that the SO and SO 2 abundances evolve on timescales less than about 1000 yr. The computed values of timescales toward individual pixels are consistent with the ages estimated in L1157 B0-B1 by previous studies. Furthermore, the modeling requires high abundances of both H 2 S/H and S/H (2×10 -6 ) injected in the gas phase due to the shock occurrence, so pre-frozen OCS only is not enough to reproduce our new observations.

We thank the NOEMA, SMA, and IRAM-30 m staff for their helpful support during the observations and data reduction of the SOLIS project. We also thank the referee for constructive suggestions.

S.F. acknowledges the support of National Natural Science Foundation of China No. 11988101, the support of the CAS International Partnership Program No.114A11KYSB20160008, and the support of the EACOA fellowship from the East Asia Core Observatories Association (EACOA). EACOA consists of the National Astronomical Observatory of China, the National Astronomical Observatory of Japan, the Academia Sinica Institute of Astronomy and Astrophysics, and the Korea Astronomy and Space Science Institute.

Figure A1. Profiles of identified SO and SO 2 isotopologue lines, averaged from a beam-sized region with the center toward each clumpy substructure in the plane of the sky (Table 2). All lines are extracted from images smoothed to the same pixel size and angular resolution (5.1 ). Their velocity resolutions are listed in Table 1. In each panel, two red vertical lines at -10 km s -1 and +5 km s -1 indicate the velocity range for which we integrate the intensity; the gray dashed vertical line indicate the V sys = 2.7 km s -1 of the cloud. The horizontal cyan line indicates the baseline (T B =0 K).

Figure 1 .

 1 Figure1. Profiles of two identified SO and SO 2 lines with the same E u /k B and similar beam attenuation. The profiles of all identified SO and SO 2 lines are shown in Appendix FigureA1. Lines are averaged from a beam-sized region with the center toward each clumpy substructure in the plane of the sky (Table2). All lines are extracted from images smoothed to the same pixel size and the same angular resolution (5.1 ). The line profile of SO (6 5 -5 4 ) shown in black or dashed blue has a velocity resolution of 6.09 km s -1 (original spectral resolution) or 7.01 km s -1 (the same as that of SO 2 8 1,7 -8 0,8 ), respectively. In each panel, two red vertical lines at -10.0 km s -1 and +5.0 km s -1 indicate the velocity range for which we integrate the intensity; the gray dashed vertical line indicate the V sys = 2.7 km s -1 of the cloud. The horizontal cyan line indicates the baseline (T B =0 K).

Figure 2 .

 2 Figure2. Outline of the southern (blueshifted) outflow lobe from the low-mass protostar L1157-mm. Color maps show the intensity maps of several shock tracer lines integrated over the velocity range of -10 to +5 km s -1 , which are sensitive to the extended B0-B1-B2 region. All of the maps are obtained from the combination of NOEMA (or SMA) and IRAM-30 m, and line emissions with < 5σ rms are blanked. The red contours, starting from 8σ (σ = 0.27 Jy beam -1 ) and increasing with the step of 8σ, show the CO (1-0) emission, obtained from NOEMA and IRAM-30 m combination(Gueth et al. 1996). Gray circle(s) in each panel indicates the primary beam of the NOEMA/SMA at corresponding line rest frequency, and the angular resolution after IRAM-30 m combination is plotted at the bottom left in each panel. The black curve guidelines the path of the precessing jet from L1157-mm as modeled byPodio et al. (2016). The black squares indicate the driving source, L1157-mm, and the shock knots in B0, B1, and B2.

Figure 3 .

 3 Figure3. Color maps show the intensity maps of all of the detected SO 2 and SO lines integrated over the velocity range of -10 to +5 km s -1 , zoomed in toward the central 1 × 1 region. All the maps are obtained from the combination of NOEMA (or SMA) and IRAM-30 m, and the emission < 5σ rms are blanked. Seven clumpy substructures labeled as B1a, B1b, B1c, B0e, B1f, B1i, and B2a show strong line emissions. Gray circle(s) in each panel indicates the primary beam of the NOEMA/SMA at corresponding line rest frequency, and the angular resolution after IRAM-30 m combination is plotted at the bottom left in each panel.

  Figure 4.Map of the integrated intensity ratio T B, 34 SO (3 2 -2 1 ) /T B, 32 SO (3 2 -2 1 ) (top widget) as well as the optical depth of 32 SO (3 2 -2 1 ) (right widget) toward each pixel. The gray dashed contours, starting from 8σ (σ = 0.27 Jy beam -1 ) and increasing with the step of 8σ, show the CO (1-0) emission(Gueth et al. 1996). The pixels where 34 SO (3 2 -2 1 ) shows < 6σ emissions are blanked. Six clumpy structures toward B1 are labeled.

  Map of relative abundance ratio between SO and SO 2

T)Figure 5 .

 5 Figure 5. Rotation diagrams of SO (le f t) and SO 2 (right) toward B1a as an example. The black dots correspond to the line emission integrated in the velocity range of -10 to +5 km s -1 . The error bar of each data point comes from the uncertainty in observations and data reduction (calibration, combination, and intensity integration). The black line is fitted with the assumption that lines from the same species are optically thin, under LTE, and trace a single component toward B1a. The gray dashed line and gray dot in the le f t panel is fitted by taking the optical depth of SO (3 2 -2 1 ) into account. The gray dashed line in the right panel is fitted by only taking the low-J and mid-J lines with E u /k B < 70 K into account.

Figure 8 .

 8 Figure8. SO/SO 2 abundance ratio predicted by our astrochemical model (see text) as a function of the time for different H 2 S/H and S/H abundances, 2 × 10 -6 (red curves) and 1 × 10 -8 (green curves). The model predictions are also obtained for two values of the cosmic rays ionization rates, 6 × 10 -16 s -1 (solid curves) and 3 × 10 -16 s -1 (dashed curves), and a smaller H density, 2 × 10 4 cm -3 (dotteddashed curve). The results from the analysis of the observations are reported in the dashed areas toward three positions along the B1 region: B1a, B1b, and B1i, respectively.

Table 1 .

 1 Target lines in this work with their spectroscopic and observation parameters Unpublished data was taken from NOEMA observations W068 reported byFontani et al. (2014a). d. Unpublished data was taken from NOEMA observations X058 reported byPodio et al. (2017).

	Mol.	Freq. a	Transition a	Sµ 2	E u /k B	a	∆V	∆θ (P.A.) b	rms b	Primary beam Interferometric
		(GHz)		(D 2 )	(K)		(km s -1 )	( × , • )	(mJy beam -1 ch -1 )	( )	Observations
	SO 2	83.688	8 1,7 -8 0,8	17.01	37		7.01	5.30 × 4.86 (-94 • )	0.7	61	SOLIS-NOEMA
	SO 2	134.004	8 2,6 -8 1,7	15.18	43		4.37	2.64 × 2.43 (-89 • )	1.1	38	NOEMA c
	SO 2	135.696	5 1,5 -4 0,4	8.35	16		4.32	2.65 × 2.51 (75 • )	1.3	38	NOEMA c
	SO 2	163.606 14 1,13 -14 0,14 17.10	102		3.58	2.01 × 1.87 (43 • )	1.4	31	NOEMA d
	34 SO 2 162.776	7 1,7 -6 0,6	11.86	26		3.60	1.98 × 1.86 (47 • )	0.9	31	NOEMA d
	SO	99.230	3 2 -2 1	6.91	9		5.90	2.95 × 2.20 (25 • )	0.4	52	SOLIS-NOEMA
	SO	206.176	4 5 -3 4	8.91	39		2.84	1.15 × 0.80 (-171 • )	1.3	25	SOLIS-NOEMA
	SO	215.221	5 5 -4 4	11.31	44		6.23	3.91 × 3.66 (-2 • )	21.9	55	SMA
	SO	219.949	6 5 -5 4	14.01	35		6.09	3.82 × 3.58 (-3 • )	32.1	53	SMA
	SO	258.256	6 6 -5 5	13.74	56		5.19	3.23 × 3.04 (-4 • )	22.2	45	SMA
	SO	261.844	7 6 -6 5	16.38	48		5.12	3.18 × 3.00 (-4 • )	73.1	45	SMA
	34 SO	97.715	3 2 -2 1	6.92	9		5.99	2.98 × 2.33 (25 • )	0.4	52	SOLIS-NOEMA
	34 SO	135.776	4 3 -3 2	9.28	16		4.31	2.64 × 2.50 (75 • )	1.2	38	NOEMA c
	SiO	217.105	5 -4	48.0	32		6.17	3.99 × 3.75 (-3 • )	22.7	54	SMA

Note. a. Spectroscopic parameters are taken from the CDMS.

b. Data from NOEMA-30 m or SMA-30 m combination, without angular resolution smoothing, measured with the referred beam size ("beam") and channel width ("ch"). c.

Table 2 .

 2 Positions corresponding to the clumpy structures shown on the maps in Figure3.

	R.A. (J2000)	Dec (J2000)
	B0e 20 h 39 m 10 s .365	68 • 01 23.80
	B1a 20 h 39 m 10 s .301	68 • 01 13.51
	B1b 20 h 39 m 08 s .734	68 • 01 10.54
	B1c 20 h 39 m 09 s .677	68 • 01 06.45
	B1f 20 h 39 m 09 s .444	68 • 01 16.70
	B1i 20 h 39 m 10 s .000	68 • 01 01.11
	B2a 20 h 39 m 12 s .322	68 • 00 43.85

Table 3 .

 3 Best-fit results by using RD and LVG methods

				RD				LVG
	Location Species	N rot (10 14 cm -2 )	T rot (K)	N rot,τ (10 14 cm -2 ) a	T rot,c (K)	b	N T (10 14 cm -2 )	T kin (K)	n (10 5 cm -3 )	χ 2
	B0e	SO SO 2	2.1 ± 0.8 1.0 ± 0.8	13 ± 1 28 ± 9	1.8 ± 1.3 ...	... 12 ± 3	2.0 ± 0.2 --c	70-200 --c	0.3-0.5 --c	0.94 --c
	B1a	SO SO 2	2.4 ± 0.7 1.2 ± 1.0	14 ± 1 31 ± 8	3.9 ± 1.7 ...	... 15 ± 4	2.5 ± 0.5 1.4 ± 0.4	40-70 50 ± 10	0.8-0.9 3 ± 2	0.6 0.5
	B1b	SO SO 2	2.3 ± 1.4 1.2 ± 1.0	13 ± 1 25 ± 5	3.5 ± 2.0 ...	... 15 ± 4	2.0 ± 0.2 1.0 ± 0.2	30-100 30 ± 20	0.8-1 10 ± 5	1.2 0.7
	B1i	SO SO 2	1.4 ± 0.8 1.4 ± 0.7	13 ± 1 28 ± 9	1.5 ± 1.0 ...	... 12 ± 3	1.4 ± 0.2 1.2 ± 0.2	30-100 50 ± 5	0.8-1 3 ± 1	0.5 0.5
	B2a	SO	0.8 ± 0.5	19 ± 3	...	...		--c

OBSERVATIONS1Podio et al. (2016) estimated the kinematic ages of each knot by using d ∼250 pc. The kinematic ages of these knots in this work are updated by using d ∼352 pc(Zucker et al. 2019).2 The "Seeds Of Life In Space" large programme is aiming at investigating molecular formation and processes during the early stages of the star formation process, by using IRAM NOrthern Extended Millimeter Array (NOEMA).

http://www.iram.fr/IRAMFR/GILDAS

The MIR package was originally developed for the Owens Valley Radio Observatory, and is now adapted for the SMA, http://cfawww.harvard.edu/∼cqi/mircook.html.

We compared the integrated intensity of each line over the velocity range of -10 to +5 km s -1 (S/N > 3) with -20 to +15 km s -1 (line emission down to zero), and found the difference between these integrations is < 5%, which is less than the systematic uncertainty from observations and data reduction.

The combination of two-point mosaic of SMA observations and IRAM-30 m observations recover extended emission of SiO (5-4), which that was

The column density of C 34 S given inBachiller & Pérez Gutiérrez (1997) was a typo, and it is corrected as "1.8(13)".

Because of the line image sensitivity, the optical depth correction only works for the SO (3 2 -2 1 ) line toward the B0-B1 region. As a compromise for inspecting the map of column density ratio between SO and SO 2 in a larger field, we assume that all lines are optically thin and under LTE in Figure6.

Figure A1. (continued)

This work was also supported by (i) the French program "Physique et Chimie du Milieu Interstellaire" (PCMI) funded by the Conseil National de la Recherche Scientifique (CNRS) and Centre National d'Etudes Spatiales (CNES), (ii) by the Italian project PRIN-INAF 2016 The Cradle of Life -GENESIS-SKA (General Conditions in Early Planetary Systems for the rise of life with SKA), (iii) by the program PRIN-MIUR 2015 STARS in the CAOS -Simulation Tools for Astrochemical Reactivity and Spectroscopy in the Cyberinfrastructure for Astrochemical Organic Species (2015F59J3R, MIUR Ministero dell'Istruzione, dell'Università della Ricerca e della Scuola Normale Superiore), (iv) by the French Agence Nationale de la Recherche (ANR), under reference ANR-12-JS05-0005, and (v) by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, for the Project "The Dawn of Organic Chemistry" (DOC), grant agreement No. 741002, (vi) the Ministry of Science and Technology (MoST) of Taiwan (grant No. 108-2112-M-001-002-MY3). This research made use of NASA's Astrophysics Data System.

Software: GILDAS/CLASS (Pety 2005), MIR (Scoville et al. 1993), MIRIAD (Sault et al. 1995)

APPENDIX

Figure A1 shows the profile of the SO and SO 2 isotopologue lines we use to measure the molecular column densities toward seven positions.

Table A1 lists the NOEMA observation dates, spectral setup, configurations, calibrators, and weather conditions.

Feng et al.