© Y. Gong et al. 2021

1 Introduction

Extensive observational work has established that filamentary structures are ubiquitous in the interstellar medium (e.g., Schneider & Elmegreen 1979; Myers 2009; André et al. 2010; Molinari et al. 2010; Miville-Deschênes et al. 2010; Li et al. 2013, 2016; Schisano et al. 2020; Wang et al. 2020). The properties of these filaments are rather diverse, spanning ~0.1–100 pc in length and ~1–10⁵ M_⊙ in mass (e.g., Schneider et al. 2010; Hacar et al. 2013, 2018; Li et al. 2013, 2016; Kainulainen et al. 2016; Mattern et al. 2018; Zhang et al. 2019). A strong correlation between the spatial distribution of dense cores and dense filaments was observed in molecular clouds (e.g., Könyves et al. 2015; Marsh et al. 2016; Lin et al. 2019). Molecular line observations suggest that material is flowing along filaments to feed nascent stars and clusters as well as dense cores in filament hubs (e.g., Kirk et al. 2013; Peretto et al. 2014; Liu et al. 2015; Yuan et al. 2018; Lu et al. 2018; Gong et al. 2018). In nearby molecular clouds, the filament mass function and filament line mass (mass per unit length) function are consistent with a Salpeter-like mass function, indicating a possible link to the initial mass function (André et al. 2019). The core and star formation occurring in filaments may also cause its nearly constant efficiency in dense gas (~4.5 × 10⁻⁸ yr⁻¹, Shimajiri et al. 2017), which in turn regulates the so-called star-formation law in galaxies (Gao & Solomon 2004; Wu et al. 2005). These results strongly support the idea that filaments play a crucial role in the star-forming processes (e.g., André et al. 2010, 2014; Könyves et al. 2015; André 2017).

Molecular line observations have shown that filaments can consist of bundles of velocity-coherent transonic substructures (“fibers”) with a typical length of ~0.5 pc in both low-mass and high-mass filaments (e.g., Hacar et al. 2013, 2017, 2018; Shimajiri et al. 2019a; Chen et al. 2020a). Such transonic structures can be as long as a few parsecs (e.g., Hacar et al. 2016a; Gong et al. 2018). Most dense cores (i.e., star-formingseeds) are found in these transonic structures. The formation of transonic structures marks the initial conditions upon which gravitational collapse occurs. Many mechanisms have been proposed to explain filament formation and evolution, which can be divided into three classes of models in which gravity, turbulence, or magnetic fields play a dominant role, respectively (e.g., Padoan et al. 2001; Burkert & Hartmann 2004; Nakamura & Li 2008; Smith et al. 2014; Li et al. 2014; Gómez & Vázquez-Semadeni 2014; Inoue et al. 2018; Vázquez-Semadeni et al. 2019; Shimajiri et al. 2019b).

The Serpens filament is a velocity-coherent, transonic filament (Gong et al. 2018). At a distance of ~440 pc (Ortiz-León et al. 2017, 2018; Zucker et al. 2019), the Serpens filament, which is one of the nearest infrared dark clouds, is a part of the Serpens main molecular cloud (Burleigh et al. 2013; Gong et al. 2018; Su et al. 2019, 2020). This filament is found to be at the onset of slightly supercritical collapse with a line mass of 36–41 M_⊙ pc⁻¹ and H₂ column densities of >6 × 10²¹ cm⁻² (Gong et al. 2018). Its relative proximity allows for detailed analyses even with single-dish observations. Figure 1 presents a three-color composite image of the Serpens filament. Three embedded young stellar objects (YSOs) and seven dust cores have been revealed by the 1.1 mm Bolocam continuum and Spitzer legacy surveys (Enoch et al. 2007, 2009). Follow-up molecular line observations have demonstrated the filament’s overall extremely quiescent nature and also exceptionally slow radial and longitudinal infall motions (Gong et al. 2018). This is likely due to a scenario in which most of the filament’s material is still at an early evolutionary stage when collapse has just begun. Furthermore, the filament exhibits rather simple geometric and velocity structures, and its southeastern part is truly starless and thus not affected by stellar feedback. Given these characteristics, the Serpens filament lends itself for detailed studies of its physical/chemical properties and its underlying dynamical state, which, in turn, should shed light on the general mechanism of filament formation and early evolution. However, the angular resolution of previous molecular line observations (~1′, or 0.13 pc) was insufficient to resolve its radial distribution. We thus performed further observations with higher angular resolution (~30′′, ~0.07 pc), and the outcome is reported in the present paper.

In Sect. 2 we describe our observation and the archival data used in our study. In Sect. 3 we present the results, which are extensively discussed in Sect. 4. Our findings are summarized in Sect. 5.

Fig. 1 Three-color composite image of the Serpens filament (Herschel 70 μm: blue, Herschel 250 μm: green, Herschel 500 μm: red). The three green pluses indicate the positions of the three embedded YSOs, emb10 (also known as IRAS 18262+0050), emb16, and emb28 (Enoch et al. 2009), while the crosses mark the positions of 1.1 mm dust cores (Enoch et al. 2007). The beam size (36.′′3) for the 500 μm emission is shown in the lower right corner. The regions mapped with the IRAM-30 m and APEX telescopes are indicated by the red and white dashed boxes, respectively.

2 Observations and data reduction

2.1 IRAM-30 m and APEX-12 m observations

We carried out ¹³CO (1–0), C¹⁸O (1–0), and C¹⁷O (1–0) observations with the Institut de Radioastronomie Millimétrique 30-m (IRAM-30 m¹) telescope during 2018 May 9–11 and 2019 May 16–20. The EMIR dual-sideband and dual-polarization receiver (E090) was used as frontend (Carter et al. 2012), while the autocorrelator VErsatile SPectrometer Assembly (VESPA) spectrometers and Fast Fourier Transform Spectrometers (FFTSs) were simultaneously used as backend. VESPA was employed to achieve high spectral resolution, while the FFTSs were used to cover wider frequency ranges. Each VESPA spectrometer, with a usable bandwidth of ~35.5 MHz, provides 1817 channels, resulting in a channel width of 19.5 kHz, corresponding to 0.052 km s⁻¹ for C¹⁸O (1–0). Each FFTS provides spectral channels of 50 kHz, translating to a velocity spacing of 0.14 km s⁻¹ at 112 GHz.

Observations of ¹³CO (2–1), C¹⁸O (2–1), and C¹⁷O (2–1) were performed with the Atacama Pathfinder EXperiment 12 meter submillimeter telescope (APEX² ; Güsten et al. 2006) during 2019 April 29–May 1 (project code: M9511A_103). The PI230 two sideband (2SB) and dual-polarization receiver, built by the Max-Planck-Institute for Radio Astronomy, was used as frontend, while the facility FFTS was used as backend (Klein et al. 2012). Each of the 4 FFTS modules has a bandwidth of ~4 GHz and provides 65 536 channels, which results in a channel width of 61 kHz, corresponding to 0.083 km s⁻¹ for the C¹⁸O (2–1) line (219.6 GHz). During the IRAM-30 m and APEX observations, we performed on-the-fly (OTF) maps toward the Serpens filament, and the mapping was scanned alternatively along the long and short axes of the filament in order to reduce striping effects in the image restoration. The adopted off position (α_J2000 = 18^h27^m42.^s72, δ_J2000 = 00°44′54.′′49) was checked to be emission-free in ¹³CO (1–0) with a 1σ rms noise level of 96 mK at a channel width of 0.14 km s⁻¹.

The standard chopper-wheel method was used to calibrate the antenna temperature (Ulich & Haas 1976). The calibration was done about every 10 min. We converted the antenna temperature, $T_{\text{A}}^{*}$, to the main-beam brightness temperature, T_mb, with the relation $T_{\text{mb}}=T_{\text{A}}^{*}{\eta }_{\text{f}}/{\eta }_{\text{mb}}$, where η_f and η_mb are the forward and main beam efficiency, respectively. These efficiencies are based on the antenna efficiency reports of the IRAM-30 m³ and APEX ⁴ telescopes. Pointing was checked on nearby continuum sources every hour, and was found to be accurate to within ~3′′. The uncertainties of the absolute flux calibration are smaller than 10%. A dispersion of <5% between both polarizations is found between the peak intensities of the same molecular line.

Data reduction and analysis were performed with the GILDAS⁵ package (Pety 2005). Convolving with a Gaussian kernel of 1/3 half-power beam width (HPBW), we gridded the raw data into data cubes with regular separations of 5′′ between two adjacent pixels. A first-order baseline was subtracted from each spectrum. HPBWs after gridding, θ_B, typical noise levels, σ, and channel widths, δv, are given in Table 1. Velocities are given with respect to the local standard of rest (LSR) throughout this paper.

2.2 Archival data

We make use of the dust temperature and H₂ column density maps from the Herschel Gould Belt (GB) Survey⁶ (André et al. 2010; Könyves et al. 2015; Fiorellino et al. 2017). These maps were derived from spectral energy distribution (SED) fitting to data at four far-infrared bands, that is, 160, 250, 350, and 500 μm. The data were smoothed to a common angular resolution of 36.′′3. Based on the results of Könyves et al. (2015), typical uncertainties of 0.5 K in dust temperature and 10% in H₂ column density are assumed in our analysis.

3 Results

3.1 Overall distribution

Figure 2 presents maps of the integrated intensity of the emission in the ¹³CO (1–0), C¹⁸O (1–0), C¹⁷O (1–0), ¹³CO (2–1), C¹⁸O (2–1), and C¹⁷O (2–1) lines resulting from our observations. They provide a factor of two finer linear resolution than our previous C¹⁸O (1–0) observations (Gong et al. 2018). The red dashed line in Fig. 2d divides two distinct regions which we call northwest (NW) and southeast (SE) hereafter. The C¹⁸O and C¹⁷O line emission shows similar distributions, which trace the filamentary structure well. In contrast, the morphology of the emission of the two ¹³CO lines is markedly different. Much of the filamentary structure is not seen in both ¹³CO transitions, in line with previous observations (Burleigh et al. 2013). Also, ¹³CO emission is significantly weaker in SE than in NW. Inspecting their corresponding spectra (see Figs. 3a–c and also Appendix A), we find that C¹⁸O (2–1) is generally brighter than ¹³CO (2–1) near the systemic velocity of ~8 km s⁻¹. This is likely because of significant ¹³CO self-absorption, which is reinforced by the presence of the extremely narrow dip in ¹³CO spectra at velocities where C¹⁸O and C¹⁷O peak. Such narrow ¹³CO (2–1) self-absorption features were not reported in previous ¹³CO (2–1) observations (Burleigh et al. 2013), probably due to the coarser spectral resolution of these observations. The widespread ¹³CO (2–1) narrow self-absorption suggests high optical depths and lower excitation temperature in the outskirts of the filament or in the ambient gas.

3.2 A new molecular outflow driven by Ser-emb 28

Based on our APEX ¹³CO (2–1) observations which are more sensitive than our ¹³CO (1–0) observations, we uncovered a molecular outflow based on its line wing emission (see Fig. 4a). Figure 4b displays the distribution of the blueshifted and redshifted components which are bluer than 5.5 km s⁻¹ and redder than 9.4 kms⁻¹, respectively.Although both blueshifted and redshifted lobes are only marginally resolved, the morphology is suggestive of a bipolar molecular outflow which is nearly oriented in the north-south direction. The outflow direction is at ≲10° with respect to the filament. There is an overlap between the blueshifted and redshifted lobes, either caused by our viewing angle or by our limited angular resolution. The YSO Ser-emb 28 is located at this overlap, indicating that the molecular outflow is driven by Ser-emb 28 which was identified as a Class I YSO (Enoch et al. 2009). Figure 4c presents the position-velocity diagram along the outflow axis. High-velocity wing emission is evident. The brightest blueshifted and redshifted wing emission is located at about 20′′ (i.e., 8800 au) and 10′′ (i.e., 4400 au) away from Ser-emb 28, respectively. The blueshifted lobe is found to be stronger than the redshifted lobe. Furthermore, strong ¹³CO (2–1) self-absorption is observed around the systemic cloud velocity derived from C¹⁸O emission.

The velocity interval of the outflow emission is determined by manually investigating the ¹³CO (2–1) channel map (see Appendix B). Assuming that emission in the high velocity wings of the ¹³CO (2–1) line is optically thin and is in local thermodynamic equilibrium (LTE) at an excitation temperature of 20 K, we can derive ¹³CO column densities in the outflow lobes following Mangum & Shirley (2015). Previous measurements suggest a range of 10–50 K for the excitation temperature in molecular outflows (e.g., Wu et al. 2004). For this range of excitation temperatures, our assumption of a constant 20 K for the excitation temperature would give a ≲30% uncertainty in the derived ¹³CO column densities. In order to derive the total molecular mass of the outflow lobes, the ¹³CO fractional abundance relative to H₂ is assumed to be 1.1 × 10⁻⁶ (e.g., Cabrit & Bertout 1992; Beuther et al. 2002) and the mean molecular weight per hydrogen molecule is taken to be 2.8 (e.g., Kauffmann et al. 2008). Arbitrarily, an inclination angle of 45° with respect to the line of sight is assumed to correct the outflow parameters. Following previous studies (e.g., Beuther et al. 2002), we also determine the momentum, the kinetic energy, the dynamical timescale, the mass entrainment rate, the mechanical force, and the mechanical luminosity. These physical properties are given in Table 2. They are comparable to those of outflows from other Class I YSOs (e.g., Mottram et al. 2017). Except for the dynamical timescale, these properties have much lower values than found for outflows in high-mass star-forming regions (e.g., Yang et al. 2018). Comparing the outflow’s kinetic energy with the turbulent energy of this filament (see Appendix C), we find that the kinetic energy contributes ~5% of the filament’s total turbulent energy and ~7% of NW’s turbulent energy. The outflow mechanical luminosity is comparable to the turbulent dissipation rate of this filament (see also Appendix C). Hence, the outflow will be able to maintain the observed level of turbulence. We also note that the outflow energy might be dissipated outside the filament, which would reduce the outflow’s contribution to the observed turbulence.

Fig. 2 Integrated intensity maps of 13CO (1–0) (a), C18O (1–0) (b), C17O (1–0) (c), 13CO (2–1) (d), C18O (2–1) (e), and C17O (2–1) (f). Their integrated velocity ranges are [6.5, 10] km s−1 for 13CO (1–0), C18O (1–0), 13CO (2–1), and C18O (2–1), [7.5, 9] km s−1 for C17O(1–0), and [6, 10] km s−1 for C17O (2–1). The contours start at 0.4 K km s−1 with a step of 0.4 K km s−1 for 13CO (1–0), C18O (1–0), 13CO (2–1), and C18O (2–1). For the two C17O maps, the contours start at 3σ with a step of 3σ, where 1σ = 0.08 K km s−1 for C17O (1–0) and 1σ = 0.06 K km s−1 for C17O (2–1). The color bar represents the C18O (1–0) and C18O (2–1) integrated intensities in units of K km s−1. The 13CO (1–0), 13CO (2–1), C17O (1–0), and C17O (2–1) maps have been scaled by 0.4, 0.5, 3, and 3 to match the color bar. Panel d: the red dashed line is used to divide the Serpens filament into two parts, NW and SE, and the blue line represents the crest of the Serpens filament extracted in the corresponding H2 column density map (Gong et al. 2018). In all panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses indicate the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (see Enoch et al. 2009, and Fig. 1), and the white crosses mark the positions of the seven dust cores (see Enoch et al. 2007, and Fig. 1).

3.3 Decomposition

Since ¹³CO spectra suffer from self-absorption (see Sect. 3.1), we only use the C¹⁷O and C¹⁸O data to study the kinematics of the Serpens filament. Following previous studies (e.g., Gong et al. 2018), we apply Gaussian decomposition to the C¹⁸O (1–0) and C¹⁸O (2–1) data to obtain their peak intensities, LSR velocities, and line widths for pixels that have at least three adjacent channels with intensities higher than 3σ. For C¹⁷O data, which have hyperfine structure (HFS), we apply the HFS fitting routine embedded in GILDAS to decompose them. Due to rather low signal-to-noise ratios, the opacities derived from the HFS fitting have large uncertainties, and are thus not used in the following analysis. The spatial distributions of these observed parameters are displayed in Figs. 5 and 6.

Figures 5a and d present the peak intensity maps of C¹⁸O (1–0) and C¹⁸O (2–1), whileFigs. 6a and d present the peak intensity maps of C¹⁷O (1–0) and C¹⁷O (2–1). Similar distributions are evident. The peak intensities are 0.42–3.49 K for C¹⁸O (1–0), 0.24–2.72 K for C¹⁸O (2–1), 0.28–0.97 K for C¹⁷O (1–0), and 0.21–0.87 Kfor C¹⁷O (2–1). Assuming an excitation temperature of 7 K (see Sect. 3.4) and neglecting beam dilution effects, we obtain optical depths of 0.07–2.51 for C¹⁸O (1–0), 0.09–3.85 for C¹⁸O (2–1), 0.08–0.30 for C¹⁷O (1–0), and 0.08–0.38 for C¹⁷O (2–1). Adopting a lower excitation temperature will result in even higher optical depths. Therefore, opacity effects cannot be neglected in the following analysis. We also find features extending to the west of the filament and directed nearly perpendicular to the filament, indicated by the dotted lines in Figs. 5a and d. These features may be indicative of substructures that providemass accretion (Kirk et al. 2013; Palmeirim et al. 2013). However, these features are not evident in the C¹⁸O integrated intensity and H₂ column density maps, probably due to their low contrast.

Figures 5b, e, and 6b, e show the velocity centroid maps derived from the C¹⁸O (1–0), C¹⁸O (2–1), C¹⁷O (1–0), and C¹⁷O (2–1) spectra. These plots show global velocity gradients similar to those reported by Gong et al. (2018). The inner part of SE with v_lsr = 8.1–8.3 km s⁻¹ is blue-shifted relative to its outskirts where values of v_lsr = 8.4–8.6 km s⁻¹ are found, whereas NW shows an opposite trend with v_lsr = 7.8–8.1 km s⁻¹ in its inner part and v_lsr = 7.5–7.9 km s⁻¹ in its outskirts. NW, showing more clumpy structures and more active star formation, is more blueshifted than SE with respect to the ambient gas derived from C¹⁸O (1–0) at a linear scale of 0.13 pc (v_LSR ~ 8.6 km s⁻¹, Gong et al. 2018).

Figures 5c, f, 6c, and f exhibit line width maps derived from the C¹⁸O (1–0), C¹⁸O (2–1), C¹⁷O (1–0), and C¹⁷O (2–1) data. Due to their lower signal-to-noise ratios, the fitting results for C¹⁷O are not as good as those for C¹⁸O data and for the former isotopologue line width fits with a significance lower than 3σ are discarded. We obtain median errors of 0.02 km s⁻¹, 0.02 km s⁻¹, 0.05 km s⁻¹, and 0.04 km s⁻¹ for the fitted line widths of C¹⁸O (1–0), C¹⁸O (2–1), C¹⁷O (1–0), and C¹⁷O (2–1), respectively. We see that the C¹⁸O (1–0) and C¹⁸O (2–1) lines are broader than the C¹⁷O lines in a significant number of pixels, despite the relatively large uncertainties in the C¹⁷O fitting results. This difference can be attributed to opacity broadening effects (e.g., Hacar et al. 2016b). Corrections for the opacity broadening effect are discussed in Appendix C. The opacity broadening corrected velocity dispersion analysis confirms the (tran-)sonic nature of the inner part of the whole filament unveiled by our previous observations (Gong et al. 2018). Such (tran-)sonic nature has already been widely found in both high-mass and low-mass star-forming regions (e.g., Pineda et al. 2010; Hacar et al. 2016a; Sokolov et al. 2017; Li et al. 2020). The opacity-corrected velocity dispersion, σ_nt,c, cannot be derived for low-brightness regions, where the C¹⁸O opacity is not valid due to low signal-to noise ratios (see Appendix C). For consistency, we use non-opacity-corrected velocity dispersions, σ_nt, for discussions unless specified otherwise. Corrections for the low-brightness regions should be small due to their low opacities anyway. Most of the molecular gas has line widths of <0.6 km s⁻¹ in the inner region of SE, corresponding to σ_nt <0.25 km s⁻¹ at an assumed kinetic temperature of 7 K. Higher line widths of >0.6 km s⁻¹ (σ_nt >0.25 km s⁻¹) are observed in NW and the outer region (N_H2 ≲ 10²² cm⁻²) of SE, suggestive of higher turbulent motions. Generally, NW is more turbulent than SE, probably because NW has more powerful accretion flows. This is supported by larger velocity gradients (see Sect. 4.1). Furthermore, the higher velocity dispersion around emb28 is also attributed to the outflow feedback because of their spatial coincidence and high turbulent energy contributions (see Sect. 3.2).

Fig. 3 Observed 13CO (2–1), C18O (2–1), and C17O (2–1) spectra toward the selected positions that are indicated in the upper left corner. For the origin of the coordinate system, see the caption to Fig. 2. 13CO self-absorption is evident in these spectra.

Fig. 4 Outflow driven by Ser-emb 28. Left: same as Fig. 2e. Right: (a) observed 13CO (2–1) spectra of the two positions indicated by the blue and red pluses in panel b overlaid with the C18O spectrum of the position indicated by the yellow star in panel b. (b) 13CO (2–1) outflow map of Ser-emb 28. The blueshifted emission integrated from 2.0 to 5.5 km s−1 is shown in blue dashed contours which start at 0.20 K km s−1 (5σ) and increase by 0.20 K km s−1. The redshifted emission integrated from 9.4 to 12 km s−1 is shown in red dashed contours which start at 0.20 K km s−1 and increase by 0.20 K km s−1. The yellow star marks the position of Ser-emb 28. The beam size is shown in the lower left corner. (c) Position-velocity diagram of the 13CO (2–1) emission along the cut indicated by the black line in panel b. The contours start at 0.18 K (3σ) and increase by 0.36 K. The horizontal dashed line indicates the position of Ser-emb 28, and the vertical dashed line marks its systemic velocity of 7.87 km s−1. The resolution element is shown in the lower left corner.

Fig. 5 Maps of peak intensities (a), LSR velocities (b), and line widths (c) derived from single-component Gaussian fits to our C18O (1–0) data. Panels d–f are similar to panels a–c, but for C18O (2–1). In panels a and d, the contours represent the peak intensities which start at 0.6 K and increase by 0.6 K, while in c and f the contours represent the H2 column densities that start at 1 × 1021 cm−2 and increase by 3 × 1021 cm−2. In panels a and d, the dotted lines indicate possible substructures. The beam size is shown in the lower right corner of each panel. In all panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses indicate the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (Enoch et al. 2009), and the white crosses mark the positions of the seven dust cores (Enoch et al. 2007).

Fig. 6 Similar to Fig. 5, but for C17O (1–0) and C17O (2–1). In panels a and d, the contours represent the peak intensities which start at 0.3 K and increase by 0.3 K. The pixels with large fitting uncertainties are masked out in panels b–c and e–f.

3.4 Molecular excitation

As mentioned above, opacity effects should be taken into account in the analysis. According to the radiative transfer equation, the main beam temperature, T_mb, can be expressed as $$T_{\text{mb}}=\eta (J(T_{\text{ex}})-J(T_{\text{bg}}))(1-{\text{e}}^{-\tau }),$$(1) $$J(T)=\frac{h\nu }{k}\frac{1}{{\text{e}}^{\frac{h\nu }{kT}}-1},$$(2)

where η is the beam dilution factor that is assumed to be unity, T_ex is the excitation temperature, T_bg is the background temperature that is taken to be 2.73 K, τ is the optical depth, h is the Planck constant, k is the Boltzmann constant, and ν is the transition’s rest frequency. Because both C¹⁷O transitions are affected by HFS splitting, we used the F = 7/2–5/2 component of C¹⁷O (1–0) and the combined component, Group I (see details in Appendix D) of C¹⁷O (2–1) in Eq. (1), and their optical depths are 0.4444 and 0.6933 times the total optical depths of their corresponding total C¹⁷O transitions. We also note that these components can be overestimated due to overlapping HFS components. This effect can result in up to 28% and 3% uncertainties (see details in Appendix D). We thus introduce additional 28% and 3% uncertainties for the peak intensities for the respective component of C¹⁷O (1–0) and C¹⁷O (2–1) in the following calculations.

In order to constrain molecular excitation, we use the rotational diagram method under the assumption of LTE. By using the optical depth correction factor, C_τ = $\frac{\tau }{1-e^{-\tau }}$, the level populations result in the following relation (e.g., Goldsmith & Langer 1999): $$N_{\text{tot}}=\frac{3kW}{8{\pi }^3\nu {\mu }^{2}S}Q(T_{\text{rot}})e^{\frac{E_{\text{u}}}{k{T}_{\text{rot}}}}{C}_{\tau },$$(3)

where μ is the permanent dipole moment of the corresponding molecule, S is the transition’s intrinsic strength, N_tot is the total molecular column density, W is the integrated intensity of the corresponding transition, T_rot is the rotational temperature, Q is the partition function for the corresponding molecule, and E_u is the upper level energy of the transition. We took the values of Q and μ² S from the CDMS (Müller et al. 2005). We also assume a homogeneous [¹⁸O/¹⁷O] isotopic ratio, r_O, in this filament. This gives us the following relations: $$r_{\text{O}}=N_{18}/N_{17},$$(4) $$r_{\text{O}}={\tau }_{18,i}/{\tau }_{17,i},$$(5)

where N₁₈ and N₁₇ are the molecular column densities of C¹⁸O and C¹⁷O, and τ_18,i and τ_17,i are the peak optical depths of the corresponding C¹⁸O and HFS-corrected C¹⁷O transitions. In Eq. (5), we neglect the small difference (~3%) between the column density ratios and opacity ratios. In order to make use of Eqs. (4) and (5), we first determined the [¹⁸O/¹⁷O] isotopic ratio in low brightness regions. For this, we created a mask where the C¹⁸O (1–0) and C¹⁸O (2–1) peak intensities are within 0.27–1.0 K and 0.2–0.7 K. This enables the pixels in the mask to have at least 3σ while their optical depths are lower than 0.3 at an assumed excitation temperature of 7 K. Therefore, we can assume all C¹⁸O and C¹⁷O transitions to be almost optically thin in the mask. In order to improve the signal-to-noise ratios (S/N), we averaged the observed spectra in the low-brightnesstemperature mask and the average spectra are presented in Fig. 7. We determined that the integrated intensity ratios $\frac{{\displaystyle \int T_{18}}\text{d}v}{{\displaystyle \int T_{17}}\text{d}v}$ are 3.94 ± 0.23 and 4.41 ± 0.24 for J = 1–0 and J = 2–1, respectively.Using an opacity cut of 0.3 to define the mask means that we may underestimate the true [¹⁸O/¹⁷O] isotopic ratio by a factor of 1.15. Using a cut at 0.2 instead would reduce this potential bias. A test shows that the resulting values of the ratios do not change significantly, but the uncertainties increase by a factor of 4. The [¹⁸O/¹⁷O] isotopic ratio is thus 3.94–4.41, which is in agreement with previous measurements in other nearby molecular clouds (Wouterloot et al. 2005, 2008; Zhang et al. 2007, 2020; Giannetti et al. 2014). We therefore use a constant r_O value of 4.1 in Eqs. (4) and (5).

In order to solve Eqs. (1)–(5) for both molecules simultaneously, we assume that all C¹⁸O and C¹⁷O transitions have the same excitation (or rotational) temperature. As already mentioned above, we fixed the r_O value to be 4.1. For the purpose of improving the signal-to-noise ratios and the following comparison with the dust temperature and dust-based H₂ column density maps, we first convolved C¹⁸O and C¹⁷O peak and integrated intensity maps to the angular resolution of the Herschel maps (i.e., 36.′′3). For pixels with signal-to-noise ratios higher than 3, we solved Eqs. (1)–(5) with the emcee⁷ code (Foreman-Mackey et al. 2013) to perform the Monte Carlo Markov chain (MCMC) calculations with the affine-invariant ensemble sampler (Goodman & Weare 2010) in order to obtain the physical parameters and their uncertainties simultaneously. We assumed uniform priors for $N_{{\text{C}}^{18}\text{O}}$, T_rot, τ_18,1, and τ_18,2. The posterior distribution of these parameters is given by the product of the prior distribution function and the likelihood function. The likelihood function is assumed to be $\propto {\text{e}}^{-{\chi }^2/2}$ with $${\chi }^2={\Sigma }_i{(I_{\text{obs},\text{i}}-I_{\text{mod},\text{i}})}^2/{\sigma }_{\text{obs},\text{i}}^2,$$(6)

where I_obs,i and I_mod,i are observed and modeled parameters including peak and integrated intensity, and σ_obs,i is the standard deviation of I_obs,i. The MCMC simulations were run with 20 walkers and for 4500 steps after the burn-in period. The 1σ uncertainties are given by the 16th and 84th percentiles of the posterior distribution. An example of the fitting result is shown in Fig. 8.

The optical depth maps of C¹⁸O (1–0) and C¹⁸O (2–1) are shown in Fig. 9. The C¹⁸O (1–0) opacities lie in the range of 0.3–2.9, while the C¹⁸O (2–1) opacities are 0.4–2.5. It is surprising that we found more than 80% of pixels with the C¹⁸O (1–0) and C¹⁸O (2–1) optical depths higher than unity (see Fig. 10), because these C¹⁸O lines are usually assumed to be optically thin in the literature (e.g., Ikeda & Kitamura 2009; Gong et al. 2016, 2018). The highest opacities are found around Bolo6. Such high C¹⁸O opacities of >1 were also reported in other infrared dark clouds (e.g., Sabatini et al. 2019).

Figure 11a shows the distribution of the rotational temperature that spans 4.6–7.6 K with a median of 6.2 K. All rotational temperatures are lower than 8 K. A pixel-by-pixel comparison between the dust temperature and the rotational temperature is shown in Fig. 12a. The dust temperature is generally higher than the rotational temperature, and decreases with increasing rotational temperatures. A linear fit results in the empirical relation, T_d (K) = (−0.8 ± 0.1)T_rot (K) + (18.5 ± 0.2). Theoretical studies suggest that the gas temperature is comparable to or slightly higher than the dust temperature at an H₂ density of > 10^4.5 cm⁻³ due to the gas-dust coupling, while the gas temperature becomes much higher than the dust temperature at lower H₂ densities due to the relatively inefficient gas cooling (Goldsmith 2001). The observed empirical relation appears to be inconsistent with the theoretical prediction. The reasons can be: (1) the background and foreground contributions to the dust emission can lead to a higher dust temperature than gas temperature from the SED fitting approach. This is supported by previous radiative transfer simulations (e.g., Nielbock et al. 2012); (2) C¹⁸O and C¹⁷O transitions become sub-thermally excited in low-density regions (n_H2 ≲ 10³ cm⁻³), leading to the result that the rotational temperature is lower than the corresponding kinetic temperature. The latter effect becomes more prominent in the lower-density regime, leading to the observed anticorrelation in Fig. 12a. This is also supported by the result that rotational temperatures correlate with H₂ column densities and high rotational temperatures are found to be associated with high H₂ column densities (see Fig. 12b). We performed a linear fit to Fig. 12b, and obtain an empirical relation of $T_{\text{rot}}(\text{K})=(1.04\pm 0.03)\frac{N({\text{H}}_2)}{1\times {10}^{22}{\text{cm}}^{\text{-2}}}$ + (5.17 ± 0.03). The fact that the rotational temperature is lower than the dust temperature in high H₂ column density regions can also be partially caused by CO depletion (see discussions below). Because the densest regions with higher rotational temperatures cannot be probed by CO due to its depletion, the measured rotational temperatures come from less dense regions that have lower rotational temperatures.

Figure 11b presents the distribution of C¹⁸O column densities⁸. C¹⁸O column densities range from 9.1 × 10¹⁴ cm⁻² to 3.0 × 10¹⁵ cm⁻² with a median of 1.7 × 10¹⁵ cm⁻², which almost doubles the values reported in Gong et al. (2018). This is because Gong et al. (2018) did not take C¹⁸O opacity effects into account and assumed a slightly higher excitation temperature of 10K. In addition,the improved angular resolution may also contribute to the higher beam-averaged C¹⁸O column densities. The C¹⁸O fractional abundances relative to H₂ are directly calculated by the ratio between C¹⁸O and H₂ column densities. The latter is derived from the SED fitting with Herschel data (see Sect. 2.2). The distribution of the C¹⁸O fractional abundances is shown in Fig. 11c. The C¹⁸O fractional abundances range from 5 × 10⁻⁸ to 3.5 × 10⁻⁷ with a median of 1.6 × 10⁻⁷. The map shows that the C¹⁸O abundance drops toward the 1.1 mm dense cores. The region toward bolo3 and Ser-emb10 has the lowest abundance of ~5 × 10⁻⁸. This confirms the presence of C¹⁸O depletion in this filament(Gong et al. 2018). Setting a C¹⁸O fractional abundance of 2 × 10⁻⁷ as a threshold for C¹⁸O depletion in Fig. 11c, the areas (projected onto the plane of the sky) showing depletion are estimated to be about 0.5 pc × 0.12 pc in NW and 0.3 × 0.1 pc in SE. NW hasa larger depletion size than SE, which is likely because NW has more widespread dense gas. This reason also explains the fact that the C¹⁸O depletion size in NW is larger than the typical depletion size (≲0.1 pc) found ina single prestellar core (e.g., Bergin & Tafalla 2007).

A pixel-by-pixel comparison between H₂ column densities and C¹⁸O fractional abundances is displayed in Fig. 13. We find that the C¹⁸O fractional abundance decreases with increasing H₂ column densities. This suggests that CO depletion onto dust grains becomes more efficient in denser regions, in agreement with theoretical predictions (e.g., Bergin & Tafalla 2007). We performed a linear fit to obtain the empirical relation characterizing the decrease in the C¹⁸O fractional abundances with increasing H₂ column densities, which gives $\frac{{\chi }_{{\text{C}}^{18}\text{O}}}{1\times {10}^{-7}}=(-0.93\pm 0.02)\frac{N({\text{H}}_2)}{1\times {10}^{22}{\text{cm}}^{-2}}+(2.67\pm 0.03)$. Furthermore, the lowest C¹⁸O abundances of ~5 × 10⁻⁸ are about a factor of 6 lower than those (~3.0 × 10⁻⁷) in the ambient gas. We also find that different regions show distinct distributions in Fig. 13. Intriguingly, there are two branches at H₂ column densities of >1.5 × 10²² cm⁻². The upper branch arises from the bolo2 region, while the lower branch originates in the bolo3 region. Despite similar H₂ column densities, dust temperature, and C¹⁸O rotational temperature, the bolo2 region shows a factor of ≳2 higher C¹⁸O fractional abundance than the bolo3 region. This might be caused by their different evolutionary stages with the bolo3 region being more evolved than the bolo2 region, (see discussion in Sect. 4.4). Alternatively, the H₂ number density of the bolo3 region might be enhanced by the outflow shock, resulting in more efficient CO depletion in the bolo3 region.

Fig. 7 C18O (1–0), C17O (1–0), C18O (2–1), and C17O (2–1) spectra averaged over the low-brightness-temperature mask (details are discussed in Sect. 3.4). The two C17O spectra have been scaled by a factor of 4. The hyperfine structure components of both C17O transitions are indicated by vertical dashed lines.

Fig. 8 Posterior probability distributions of C18O column density, $N_{{\text{C}}^{18}\text{O}}$, rotational temperature, Trot, the C18O (1–0) opacity, τ18,1, and the C18O (2–1) opacity, τ18,2, toward the dense core, Bolo12, with the maximum posterior possibility point in the parameterspace shown in orange lines and points. Contours denote the 0.5, 1.0, 1.5, and 2.0σ confidence intervals.

Fig. 9 Optical depth maps of C18O (1–0) (a) and C18O (2–1) (b) overlaid with contours of signal-to-noise ratio (τ∕στ) of 2 and 3 indicated by the red and black lines, respectively. The beam size is shown in the lower right corner of each panel. In all panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses show the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (Enoch et al. 2009), and the white crosses mark the positions of the seven dust cores (Enoch et al. 2007).

Fig. 10 Histogram and cumulative distribution of the C18O (1–0) (a) and C18O (2–1) (b) opacities.

Fig. 11 (a) C18O rotational temperature map. (b) C18O column density map. (c) C18O fractional abundance map overlaid with H2 column density contours. The contours start at 6 × 1021 cm−2, and increase by 3 × 1021 cm−2. The beam size is shown in the lower right corner of each panel. In all panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses indicate the positions of the three embedded YSOs, and the white crosses mark the positions of the seven dust cores (Enoch et al. 2007).

Fig. 12 (a) Dust temperature as a function of the C18O rotational temperature. The gray error bars represent the measured uncertainties in the measured dust temperature and rotational temperature. The markers’ colors correspond to their H2 column densities indicated by the color bar. The Pearson correlation coefficient is shown in the lower left. (b) C18O rotational temperature as a function of the H2 column density. The gray error bars represent the measured uncertainties in the measured rotational temperature and H2 column density. The markers’ colors correspond to their dust temperatures indicated by the color bar. The Pearson correlation coefficient is shown in the lower right. In both panels, the black dashed line represents the linear fit to the observed trend.

4 Discussion

4.1 Local velocity gradients

Thanks to the improved angular resolution, we are able to investigate the local velocity gradients at a linear resolution of ~0.07 pc. Because the velocity centroid maps agree well with each other (see Figs. 5–6), we only use the highest signal-to-noise C¹⁸O (2–1) data to derive the local velocity gradients. Following previous studies of dense cores (e.g., Goodman et al. 1993), we determine the local velocity gradients in the Serpens filament by fitting the linear function: $$v_{\text{LSR}}=v_0+a\Delta \alpha +b\Delta \delta ,$$(7)

where v_LSR is the observed LSR centroid velocity, v₀ is the systemic LSR velocity, Δα and Δδ are the offsets in right ascension and declination, and a and b are the components of the velocity gradient along the directions of right ascension and declination. The Levenberg-Marquardt algorithm was employed to fit this function toward each block with adjacent 3 × 3 pixels where the LSR velocity has been successfully derived in order to calculate the local velocity gradient of the central pixel. The local velocity gradient ∇V has a magnitude of $|\nabla V|=\sqrt{a^2+b^2}$ and a position angle of θ_pa = arctan(a∕b). The position angle, θ_pa, increases counter-clockwise with respect to the north. The uncertainties in the local velocity gradients and their position angles are calculated through error propagation.

The derived ∇V distributions are shown in Fig. 14. The magnitude of |∇V | varies from ~0.1 km s⁻¹ pc⁻¹ around Bolo11 to ~9 km s⁻¹ pc⁻¹ in the east of emb10. Most |∇V | are higher than the values estimated from a large-scale analysis (≲1.5 km s⁻¹ pc⁻¹, Gong et al. 2018). Generally speaking, SE has a lower mean |∇V | than NW. Because the velocity gradient can be regarded as the inverse of a timescale that it takes for the edge of the filament to fall into the center, the high |∇V | may suggest that NW has higher accretion rates, leading to more active star-forming activities and higher turbulent motions as observed. This is also supported by the fact that higher |∇V | are observed in active star-forming regions (Lee et al. 2014; Chen et al. 2019, 2020b). Furthermore, it is evident that |∇V | is generally higher in the filament’s outskirts than its crest. NW shows higher |∇V | and more chaotic θ_pa than SE, which can be attributed to the YSO feedback and/or a greater impact of accretion onto the filamentfrom ambient gas. Toward the more quiescent region SE, the local velocity gradient, |∇V |, decreases from the ambient cloud (|∇V | >2 km s⁻¹ pc⁻¹) to the crest (|∇V | < 0.5 km s⁻¹ pc⁻¹). Figure 15a shows the statistical distribution of the angle between the filament’s long axis (θ_f = 147°, Roccatagliata et al. 2015)and the local velocity gradient vectors. We obtained 47% for angles ranging from 60° to 90°. If we only consider local velocity gradients with |∇V | >1.5 km s⁻¹ pc⁻¹, the percentage increases to 62%.

In order to test the observed trend, we carried out 3D Monte Carlo simulations which were then projected onto 2D to simulate the cumulative distribution function (CDF) of the expected angles that are nearly parallel (0–20°), nearly perpendicular (70–90°), or completely random (0–90°), following the method of unimodal and bimodal simulations introduced by previous studies (Hull et al. 2013; Stephens et al. 2017). The bimodal distribution includes the expected angles that are both nearly parallel and nearly perpendicular. We considered 99 bimodal cases in steps of 1% starting with 1% parallel and 99% perpendicular angles. Figure 15b compares the CDF of the observed angles and the simulated models. It is evident that the CDFs of the observed angles deviate from the unimodal model that is nearly parallel. The two CDFs of the observed angles are inconsistent with a purely parallel alignment at a >99% confidence level with a corresponding p-value of <0.001 given by the Kolmogorov-Smirnov (KS) test. However, we cannot reject the null hypothesis for the other unimodal and bimodal simulations. For the observed CDFs with ∇V > 0 and ∇V >1.5 km s⁻¹, the KS test gives a p-value of 0.82 and 0.93 against the unimodal model that is completely random, and 0.12 and 0.99 against the unimodal model that is nearly perpendicular. For the observed CDF with ∇V > 1.5 km s⁻¹, we also found high p-values of 0.90 against the bimodal simulations with >95% perpendicular angles. This demonstrates that the velocity gradients with ∇V >1.5 km s⁻¹ are more likely perpendicular to the filament’s long axis. As observed in other filaments (e.g., Beuther et al. 2015; Dhabal et al. 2018), these vectors may be indicative of accreting molecular gas from its ambient cloud.

The velocity gradients are indicative of motions of mass flows. Because material can be accreted from either blueshifted or redshifted velocities, the accretion flow is either in the same direction as the arrows representing the local velocity gradients in Fig. 14b (in the case of accretion from redshifted velocities, i.e., gas coming from the foreground), or in the opposite direction (in the case of accretion from blueshifted velocities, i.e., gas coming from the background). Inspecting Fig. 14b, we find that the vectors are clearly converging toward at least two dense cores, Bolo2 and Bolo12, and their magnitudes are changing significantly, ruling out the possibility of bulk motions being dominated by the filament’s solid-body rotation. Furthermore, ∇V is also converging toward ~30′′ south of Bolo6 where the H₂ column density is not enhanced. This indicates core accretion at a very early evolutionary stage. The ∇V vectors are also wrapping the filament’s southern end, providing additional observational support for the edge effect of a finite filament (Burkert & Hartmann 2004; Heitsch 2013; Yuan et al. 2020).

Because of the coarse angular resolution (~5′), the polarization of the Planck 353 GHz thermal dust emission is only employed to measure the mean magnetic field direction (Planck Collaboration Int. XX 2015; Planck Collaboration Int. XIX 2015). The result is shown in Fig. 14b. It is evident that the mean magnetic field direction is perpendicular to the filament’s long axis, and parallel to the velocity gradients of the mass flows accreted onto the filament’s crest. Given the coarse angular resolution of the Planck polarization measurements, further higher angular resolution polarization observations are needed to better assess the role of magnetic fields in the Serpens filament.

Fig. 13 Left: C18O fractional abundance map with different labeled masks which correspond to the legend of the right panel. Right: C18O fractional abundance as a function of H2 column density. The gray error bars indicate the measured uncertainties in C18O fractional abundance and H2 column density. The colors of the markers correspond to different regions indicated in the legend. The Pearson correlation coefficient isshown in the lower left. The black dashed line represents the linear fit to the observed relation.

Fig. 14 (a) Local velocity gradient magnitude map overlaid with the H2 column density contours. (b) Herschel H2 column density map overlaid with the normalized velocity vector maps. The arrows represent the estimated local velocity gradients which are rotated by 180° in order to better visualize the accretion directions in SE. The polarized angle of the Planck 353 GHz thermal dust emission has been rotated by 90° to trace the magnetic field direction which is indicated by the red line. In both panels, the contours correspond to H2 column densities from 1 × 1021 cm−2 to 1.6 × 1022 cm−2 with a step of 3 × 1021 cm−2, and the beam sizes for the H2 column density and velocity gradient maps are shown as the gray and blue circles in the lower right corner. The three green plusesgive the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (Enoch et al. 2009), and the white crosses mark the positions of the seven dust cores (Enoch et al. 2007).

Fig. 15 (a) Histograms of the angle of the vector of the observed local velocity gradient relative to the long axis of the filament. The blue histograms show all the angles of the observed velocity gradients in SE, while the yellow histograms represent the velocity gradients with |∇V | >1.5 km s−1 pc−1. (b) Cumulative distribution function of the angle between the velocity gradients and the major axis of the filament. The solid blue and yellow lines represent the observed cumulative distribution function of relative angles with corresponding and |∇V | >1.5 km s−1 pc−1 and |∇V | >1.5 km s−1 pc−1, respectively. The three blue dashed lines present the unimodal simulations of the expected angles which are three-dimensionally parallel (0–20°), perpendicular (70–90°), or random (0–90°), while the red dashed lines show the bimodal simulations of the angles that are parallel and perpendicular with increments of 10% (i.e., the top red line is 90% parallel and 10% perpendicular, the next one is 80% parallel and 20% perpendicular).

4.2 Longitudinal profiles

Figure 16 presents the C¹⁸O rotational temperature, C¹⁸O column density, C¹⁸O fractional abundance, dust temperature, H₂ column density, nonthermal velocity dispersion, and velocity centroid⁹ profile along the crest of the Serpens filament. We find that the rotational and dust temperatures have a small variation at offsets of <0.9 pc and their median values are 6.8 and 12.9 K, respectively. This indicates that the crest is nearly isothermal with a kinetic temperature of ~7 K, because the C¹⁸O level populations are likely thermalized in the crest due to its high H₂ number density of ~4 × 10⁴ cm⁻³ (Gong et al. 2018). The velocity centroids and nonthermal velocity dispersions are also nearly constant at offsets of <0.9 pc where the velocity gradients are ≲2 km s⁻¹ pc⁻¹. On the other hand, we neither see apparent velocity oscillation along the crest at a linear resolution of ~0.07 pc, which were reported in other filaments (e.g., Hacar & Tafalla 2011; Liu et al. 2019) nor velocity discontinuities expected in shock-induced structures (e.g., Clarke et al. 2018). The observed nonthermal velocity dispersions become sonic at offsets of <0.9 pc when the sonic speed is assumed to be 0.17 km s⁻¹ at a kinetic temperature of 7 K. This implies that most parts of the crest have a very low level of turbulence. These properties indicate that this filament can be approximately regarded as a finite isothermal and quiescent filament.

Around the YSOs and 1.1 mm dust cores except Bolo4 in Fig. 16, we find a bimodal distribution of the position angle, θ_pa, of ∇V. Around these targets, we find either a nearly flat θ_pa distribution or a jump function. Depending on different geometries, accretion flows can come from either blueshifted or redshifted velocities from both sides, so such a bimodal distribution can be readily interpreted as a converging mass flow accreted by the YSOs or cores in an inhomogeous medium. The flat distribution may also be caused by cores’ rotation. However, a recent study suggests that the velocity gradient at such a linear scale of 0.07 pc (14 000 au) may not be dominated by rotation (Gaudel et al. 2020). On the other side, |∇V | is low with a typical value of <2 km s⁻¹ pc⁻¹ around these YSOs and 1.1 mm dust cores except Bolo4. Bolo4 has the highest |∇V | among these 1.1 mm dust cores, suggesting that Bolo4 has the highest accretion rate and will form the next protostar in this filament. Together with the fact that YSOs are located in its north and SE is more pristine, we speculate that star formation proceeds from north to south. Around Ser-emb28 (i.e., at offsets of ~1.19 pc), we find elevated rotational temperature and dust temperature, but with lower H₂ column density. This is readily explained by the feedback of the outflow from Ser-emb28 (see Sect. 3.2). The outflow can sweep up the ambient molecular gas, resulting in a $N_{{\text{H}}_2}$ drop at an offset of ~1.2 pc and two enhanced $N_{{\text{H}}_2}$ peaks at offsets of ~1.1 pc and ~1.35 pc. Furthermore, outflow shocks can heat the surrounding molecular gas, leading to elevated C¹⁸O rotational and dust temperatures. One should note that the elevated dust temperature and reduced H₂ column density might alternatively arise from their degeneracy during the SED fitting. As mentioned above, C¹⁸O molecules have been frozen out onto dust grains in most parts of the crest, but these C¹⁸O molecules are sublimated or enhanced by sputtering of related species in shock regions, producing the observed higher C¹⁸O column density and fractional abundance than in Ser-emb28’s nearby siblings, emb10 and Bolo3. The high nonthermal velocity dispersions can also be explained by outflow shocks which drive additional turbulence (see Sect. 3.2 and Appendix C). Furthermore, enhanced infall velocities may also contribute to the high nonthermal velocity dispersions.

4.3 Radial profiles

In order to avoid the influence of stellar feedback from YSOs in NW, we only investigate the radial distribution toward SE that is more pristine than NW. Figure 17 presents the radial dependencies of the different physical and chemical properties. Figure 17a suggests that the rotational temperatureprofile increases monotonically toward the crest, which is caused by the variation of the H₂ number density. The observed rotational temperature profile validates the assumption used to interpret the observed blue-skewed profiles in filaments as the result of infall motions (e.g., Kirk et al. 2013; Gong et al. 2018). Unlike the rotational temperature profile, the dust temperature profile decreases monotonically toward the spine (see Fig. 17b), similar to the results reported in other studies (e.g., Arzoumanian et al. 2011). This can be because the contribution of the foreground and background warm dust emission decreases toward the spine.

In Fig. 17c, we also find that the radius of r ~ 0.05 pc corresponds to the star formation threshold of A_v = 7–10 mag, above which star formation can only occur (e.g., Johnstone et al. 2004; Lada et al. 2010; Heiderman et al. 2010; André et al. 2010; Könyves et al. 2015)¹⁰. Since the velocity dispersions drop to nearly sonic levels (see Fig. 17g) at this radius, this supports the scenario that the transition tothe sonic turbulence can also be the reason for the origin of the star formation threshold. Figure 17d shows an increasing $N_{{\text{C}}^{18}\text{O}}$ toward the crest while Fig. 17e shows a decreasing ${\chi }_{{\text{C}}^{18}\text{O}}$ toward the crest. The former one is expected with increasing H₂ column densities, while the latter one can be explained by CO depletion that becomes higher in denser regions.

Figure 17f presents the radial dependence of velocity centroids. The outer part is redshifted with respect to the inner part, pointing out that mass flows accrete material predominantly from the foreground gas. We can also infer that the velocity gradients are higher in the outer part. Furthermore, the rather smooth distribution also implies that there are no strong shocks. The radial velocity difference can be used to assess the motions in the filament. As introduced by Chen et al. (2020c), the ratio between the kinetic energy of the flow transverse to the filament and the gravitational potential energy of the filament gas is: $$C_v=\frac{\Delta v_{\text{h}}^2}{GM(r)/L},$$(8)

where Δv_h is half the velocity difference across the filament out to a radius r, and G is the gravitational constant, and M(r)∕L is the mass per unit length of the filament. They pointed out that the local turbulence is more important than self-gravity when C_v ≫ 1 while the velocity structure is likely induced by self-gravity when C_v ≲ 1. Based on Fig. 17f, we find Δv_h ~0.17 km s⁻¹ for r ~ 0.17 pc where M(r)∕L was found to be 45 M_⊙ pc⁻¹. We thus obtain C_v ~ 0.15, much lower than unity. Hence, this supports gravity-driven motions in the filament.

Figure 17g presents the radial dependence of the nonthermal velocity dispersions. We find a characteristic radius of ~0.1 pc. The nonthermal velocity dispersion is nearly sonic at r < 0.1 pc, while the nonthermal velocity dispersion can have higher values and becomes more scattered at larger radii. The observed low turbulence in the inner region surrounded by more turbulent outer layers can be readily explained by the fact that turbulence has higher energy dissipation rates in denser regions (e.g., Elmegreen & Scalo 2004; Li 2017). This also supports the sonic scale of turbulence as the origin of the universal width of ~0.1 pc found in nearby molecular filaments (Arzoumanian et al. 2011, 2019). On the other hand, the nonthermal velocity dispersion is nearly constant at <0.1 pc, which can be sustained by an internally driven turbulence, for instance, accretion-driven turbulence (Heigl et al. 2020; Xu & Lazarian 2020). On the other hand, our result is opposite to the result in the massive filament DR21 which shows increasing velocity dispersion toward the crest (Schneider et al. 2010). This may be due to the different dynamical state of the DR21 filament or limited angular resolutions which are unable to resolve (sub)sonic structures.

Figure 17h presents the radial profile of |∇V | that is decreasing toward the spine. Such a radial dependence of |∇V | was also reported in the velocity-coherent filaments in NGC 1333 (Chen et al. 2020a). They proposed that the radial dependence of |∇V | is caused by the ongoing accretion damped by the high density material as the accreting gas moves closer to the filament spine. Alternatively, this can be a geometric effect (more discussions are presented in Sect. 4.5). Figure 17i shows that the position angles of ∇V become ordered at a distance of r ~ 0.1 pc away from the spine while they are more randomly distributed in the inner part at r < 0.1 pc. This is also evident in Fig. 14b. Since the uncertainties of the position angles are even smaller in the inner part than in the outskirts, the position angle transition is thus an intrinsic result of accretion flows. The ordered ∇V can be channeled by magnetic fields at r ~ 0.1 pc, while magnetic fields may become less important in channeling accretion flows at r < 0.1 pc (for a more detailed discussion, see Sect. 4.5).

Fig. 16 Physical and chemical parameters (from top to bottom: position angle of local velocity gradient, magnitude of local velocity gradient, C18O rotational temperature, C18O column density, C18O fractional abundance, dust temperature, H2 column density, nonthermal velocity dispersion, and velocity centroid) along the crest of the Serpens filament indicated by the blue line in Fig. 2d. The origin of the offsets corresponds to the southern end of the blue line in Fig. 2d. The offset increases from south to north. The positions which are close to the seven 1.1 mm dust cores and the three YSOs are indicated by the cyan and orange vertical lines, respectively. In the two lowest panels, the data derived from C18O (1–0) and C18O (2–1) are indicated by blue and red symbols, respectively.

Fig. 17 C18O rotational temperature (a), dust temperature (b), H2 column density (c), C18O column density (d), C18O fractional abundance (e), velocity centroid (f), nonthermal velocity dispersion (g), magnitude of velocity gradient (h), and position angle of velocity gradient (i) plotted as a function of their distance to the crest of SE. In c, the yellow shaded region marks the star formation threshold of Av = 7–10 mag (e.g., Johnstone et al. 2004; Lada et al. 2010; Heiderman et al. 2010; André et al. 2010; Könyves et al. 2015). In all panels, the radius of 0.1 pc is marked with a red dashed line, which indicates the transition discussed in Sect. 4.3. We note that the data points become incomplete at larger radial distances in panels a, d, e, and f due to the 3σ clipping mentioned in Sect. 3.4. In panel g, the two red horizontal dashed lines mark σt = cs and σt = 2cs where cs is the soundspeed at a kinetic temperature of 7 K.

4.4 Widespread C¹⁸O depletion and its implications

4.4.1 Chemical model

As mentioned in Sect. 3.4, our observations confirm the widespread C¹⁸O depletion, and also reveal a trend that the C¹⁸O fractional abundances decrease with increasing H₂ column densities. This can be attributed to the H₂ density dependence of the CO freeze-out timescale. In order to simulate the observed trend, we run a set of models using the pseudo-timedependent chemical model, Chempl (Du 2020). The 2012 version of the UMIST network (McElroy et al. 2013) is used, augmented by adsorption, desorption, and grain surface reactions. The initial conditions are assumed to be the standard oxygen-rich initial elemental abundances with every element in atomic form except that hydrogen is in H₂ (see Table 3of McElroy et al. 2013). The interstellar ultra-violet radiation field is set to be the standard value in the solar neighborhood (Draine 1978) and the cosmic-ray ionization rate is assumed to be the canonical value of 1.36 × 10⁻¹⁷ s⁻¹ (e.g., van der Tak & van Dishoeck 2000). The self-shielding of H₂ is implemented using the analytical formula of Draine & Bertoldi (1996), and the self-shielding of CO is treated with the “one-band” approximation (Morris & Jura 1983). Based on the rotational temperature of the thermalized C¹⁸O level populations (see Fig. 16), the kinetic temperature is simply assumed to be 7 K as the fiducial case. We adopted an isotopic ratio [¹⁶O/¹⁸O] of 530 (Wilson & Rood 1994; Bockelée-Morvan et al. 2012) to convert the ¹²CO fractional abundance to the C¹⁸O fractional abundance. Only two parameters, visual extinction (A_v) and H₂ number density, are free parameters. Since the filament is likely formed from a sheet-like molecular cloud (see Sect. 4.5), we thus assumed a constant thickness, and the value is set to be 0.68 pc, four times the filament’s FWHM width derived from the dust-based H₂ column density profile (Gong et al. 2018). The H₂ number density can be determined by dividing the H₂ column density by the thickness. We ran 12 models with A_v ranging from 2.12 to 25.4 mag. We adopt the extinction relationship as N(H₂) = 9.4 × 10²⁰ (A_v /mag) (Bohlin et al. 1978).

The pseudo-time dependent chemical modeling results are shown in Fig. 18a. This demonstrates that higher H₂ column density regions tend to be depleted earlier, as expected. Based on the modeling results, we utilize the isochronic tracks in the N(H₂)–${\chi }_{{\text{c}}^{18}\text{o}}$ plot (see Fig. 18b) in order to constrain the chemical timescale of the Serpens filament. The two branches from the bolo2 and bolo3 regions (see Sect. 3.4) are characterized by ~0.4 Myr and ~0.8 Myr, respectively,suggesting that the bolo3 region is more evolved than the bolo2 region. Overall, the observed values lie within the chemical timescale of ≲1.5 Myr. Intriguingly, the trend of the observed distribution is roughly reproduced by the isochronic track of ~0.8 Myr. However, we also note that the assumed thickness is weakly constrained and can lead to a large uncertainty in the chemical timescale. The assumed thickness directly implies the H₂ number density used in the models. If we adopt a lower value for thickness, the H₂ number density will become higher, leading to a shorter timescale. On the contrary, a greater thickness will give rise to a longer timescale. In order to investigate this impact, we run the models with a smaller thickness of 0.34 pc, twice the filament’s FWHM width, and a greater thickness of 1 pc. The results are shown in Figs. 18c–f. This demonstrates that the CO depletion is dependent on the assumed thicknesses, that is, a greater thickness will result in a longer timescale. On the other hand, the high C¹⁸O fractional abundances of ~3 × 10⁻⁷ cannot be reproduced by the model using a thickness of 0.34 pc, and a greater thickness is more plausible for these regions. For the thickness of 1 pc, the chemical timescale increases to ≲2 Myr. Because the freeze-out timescale is weakly dependent on gas temperature (e.g., McElroy et al. 2013) and the reaction rates would not change significantly with a small change in gas temperature, effects of temperature variation are neglected in our case. We also neglect the inhomogeneity of the molecular cloud which can give rise to the scatter in the N(H₂)–${\chi }_{{\text{c}}^{18}\text{o}}$ plot. We have to note that the density evolution is not taken into account in our models, which could prolong the chemical timescale. H₂ formation is a prerequisite to CO formation in chemical models (e.g., Bergin & Tafalla 2007), and the chemical timescales do not include any earlier phase when gas is predominantly molecular (e.g., H₂) but CO is not abundant enough for detectable emission.

4.4.2 Evolutionary stage

Different accretion models can also allow us to estimate and constrain evolutionary timescales. If we assume that this filament starts accreting mass from its surrounding material with an initial thermal critical line mass of 16 M_⊙ pc⁻¹, the characteristic accretion timescale to reach the observed line mass of 36–41 M_⊙ pc⁻¹ is on the order of 0.5–1 Myr based on the time-dependent accretion model (see Appendix C in André et al. 2019). For a constant accretion model, the timescale can be as low as about 0.3 Myr if we assume a constant infall rate of 72 M_⊙ Myr⁻¹ (Gong et al. 2018). Because the initial gas in the chemical models is nearly atomic, the chemical models include processes prior to the accretion phase mentioned above. Therefore, these timescales are roughly consistent with our derived chemical timescale of ≲2 Myr. On the other hand, the presence of Class I YSOs indicates an age of about 1–1.5 Myr, because statistical studies suggest that the lifetime of prestellar and Class 0+I phases are about 1 Myr and 0.5 Myr (Dunham et al. 2014; Könyves et al. 2015), respectively. This age is comparable to our chemical timescale.

The fact that the accretion timescale coincides with the chemical timescale provides hints concerning the origin of the Serpens filament. First, we note that these timescales do not necessarily need to agree with each other. For example, for a non-accreting filament, its accretion timescale would reach infinity and would far exceed the chemical timescale. In our case, both timescales seem to be comparable, which means that the accretion is occurring and contributing to its mass growth significantly. This further suggests that the Serpens filament, although only slightly supercritical, is still growing and might end up as a much more massive filament before being dispersed.

4.5 Filament formation and evolution

The Serpens filament is a long velocity-coherent (tran)sonic filament, showing similar properties to fibers. Many models have been proposed to explain such (tran)sonic filaments in molecular clouds. In the “fray and gather” model (Smith et al. 2014, 2016), filaments form as a consequence of the turbulent cascade in molecular clouds. The simulations of Smith et al. (2016) suggest that fibers form first. As the large-scale turbulence decays, gravity becomes more dominant, gathering the adjacent fibers to form a main filament. In the “fray and fragment„ model (Tafalla & Hacar 2015; Clarke et al. 2017), a filament forms by the supersonic collision of two flows, and “frayed” to produce intertwined subsonic fibers. Gravitationally unstable fibers further fragment into dense cores that collapse to form nascent stars. Filaments can be also created by a magneto-hydrodynamic (MHD) shock wave in an inhomogeneous cloud (Inoue & Fukui 2013; Inoue et al. 2018; Arzoumanian et al. 2018; Bonne et al. 2020). We note that the MHD shock wave is usually assumed to be created by colliding flows, similar to the “fray and fragment” model. Furthermore, a dense filament can form in a shell createdby a large-scale compression (e.g., Shimajiri et al. 2019b). Although there is no indication for substructures parallel to the main filament in the Serpens filament, this could result from the limited angular resolutions of our observations. Therefore, we are not able to exclude or discriminate between the “fray and gather” and “fray and fragment” models. On the other hand, we did not find any indication of strong shocks which are expected in the MHD shock scenario. However, this could be because the shock-induced energies have been already dissipated at an earlier stage, leading to the observed (tran)sonic motions. We also note that MHD shocks tend to create filaments perpendicular to the magnetic field (see Fig. 3 in Inoue et al. 2018), as indicated by the observed large-scale polarization direction (see Fig. 14b). In our observations, we do not find shell structures at the current linear resolution (~0.07 pc), so the large-scale compression is not a likely scenario for the formation of the Serpens filament. Therefore, we propose that turbulence or colliding flows may best explain the filament formation.

Sheet-like structures are very common in numerical simulations, and they can be easily created by shocks (e.g., Vázquez-Semadeni et al. 2006; Seifried et al. 2020). Previous observations also indicate that molecular clouds can be sheet-like (e.g., Qian et al. 2015; Shimajiri et al. 2019b). Thus, we assume that the filament is formed in a sheet-like molecular cloud. In order to explain the observed local velocity gradients, we hypothesize that the filament is formed not in the midplane of the sheet, but inclined with SE at the far side and NW at the near side with respect to the observer, as shown in Fig. 19. The proposed turbulence or colliding flows in an inhomogeneous cloud could create a filament in such a geometry. Once the initial filament is formed, accretion flows can be driven onto the filament due to its greater gravitational potential. Thus, our geometry hypothesis could readily explain the observed redshifted ambient gas in SE and the blueshifted ambient gas in NW. The momentum of these flows would make the velocity difference between the filament and ambient gas greater with evolution. The current velocity difference is only about 0.5 km s⁻¹, implying an early evolutionary stage.

As shown in Fig. 14b, the observed ∇V in the outskirts is parallel to the mean magnetic field that is perpendicular to the filament’s long axis. In nearly magnetically critical regions, molecular gas tends to move along magnetic field lines, because the Lorentz force tends to resist gas flows in the perpendicular direction (Nakamura & Li 2008; Li et al. 2014). In case that the filament is close to magnetically critical conditions, we can provide a rough estimate of the magnetic strength. Following Eqs. (I)–(II) in Li et al. (2014), we obtain a magnetic strength of about 8–16 μG in the corresponding regions with an H₂ column density of ~4 × 10²¹ cm⁻². In the inner region closer to the crest, molecular gas becomes magnetically supercritical, that is, gravity dominates over the magnetic fields. This leads to the observed converging ∇V morphology around Bolo2 and Bolo12. This case can also explain the observed transition that the position angle of ∇V becomes more random at r ≲ 0.1 pc (see Fig. 17i). Alternatively, simulations show that the accretion flows are able to drag the magnetic field lines along the flow in the magnetically supercritical regime, causing it to be oriented perpendicular to the filaments (Gómez et al. 2018). In this case, the magnetic strength should be lower than 8–16 μG.

The radial dependence of v_lsr and |∇V | observed in Figs. 17f and h can be at least partially a result of geometric effects. In a symmetric cylinder, one would expect a nearly constant v_lsr which is the density-weighted 3D velocity averaged over the line of sight. At r < 0.1 pc, the geometry is dominated by such a symmetric cylinder, leading to the nearly constant v_lsr and low |∇V |. At r > 0.1 pc, accretion from ambient gas becomes more important. Hence, more redshifted v_lsr and higher ∇V are observed at a greater radial distance.

The Serpens filament is apparently truncated. Different from infinitely-long filaments, material tends to pile up at finite filaments’ ends due to the edge effect, also known as gravitational edge focusing (Burkert & Hartmann 2004; Heitsch 2013). Because of the accumulated mass, the southern end becomes magnetically supercritical. Hence, the observed ∇V is wrapping the southern end of this filament (see Fig. 14b), rather than parallel to the observed mean magnetic field. Furthermore, the edge effect also gives rise to the steep component in Figs. 17f and h because of greater gravitational potentials.

Fig. 18 (a) C18O fractional abundances, ${\chi }_{{\text{C}}^{18}\text{O}}$, as a function of time, as calculated with different H2 column densities and an assumed thickness of 0.68 pc. The shaded region corresponds to the observed C18O fractional abundance range. (b) Similar to Fig. 13, but overlaid with the isochronic lines which are based on the results in a. The markers represent the observed relationship between C18O fractional abundance and H2 column density with different regions labeled in different colors. c and e are similar to a, but for an assumed thickness of 0.34 pc and 1 pc, respectively. Panels d and f are similar to panel b, but for an assumed thickness of 0.34 pc and 1 pc, respectively.

Fig. 19 Schematic view of the geometry of the Serpens filament. The left panel gives a perspective view of its 3D geometry, and the right panels represent the projection on to the x-y, x-z, and y-z planes. The midplane and its parallel surfaces have the same distances to observers, and the line of sight is along the z direction. The local velocity gradients that are associated with blueshifted and redshifted flows are indicated by the blue and red arrows, while the magnetic field is indicated bygreen dashed lines.

5 Summary and conclusion

We have performed ¹³CO (1–0), C¹⁸O (1–0), C¹⁷O (1–0), ¹³CO (2–1), C¹⁸O (2–1), and C¹⁷O (2–1) line imaging toward the Serpens filament with the IRAM-30 m and APEX-12 m telescopes in order to characterize its physical and chemical structure with a linear resolution of ~0.07 pc and a spectral resolution of ≲0.1 km s⁻¹. Our observations have led to the following main results:

• 1.

  We reveal widespread ¹³CO (2–1) self-absorption, resulting in a markedly different morphology relative to the filamentary structure traced by C¹⁸O and C¹⁷O. We discover a new molecular outflow driven by the Class I protostar Ser-emb28 in ¹³CO (2–1). This outflow is likely to exert feedback on the filament, driving a small amount (~5%) of turbulent energy in this filament. The energy injection rate from the outflow is comparable to the filament’s energy dissipation rate, indicating that the outflow can sustain the observed turbulence.

• 2.

  Averaging C¹⁸O and C¹⁷O spectra in the low-brightness temperature regions, we derive the isotopic ratio [¹⁸O/¹⁷O] to be 3.94 ± 0.23 and 4.41 ± 0.24 for J = 1–0 and J = 2–1, respectively. Assuming local thermodynamic equilibrium and fixing the isotopic ratio [¹⁸O/¹⁷O] to be 4.1,we constrain the opacity, rotational temperature, and molecular column density for C¹⁸O and C¹⁷O simultaneously with Monte Carlo Markov chain calculations. We find that the opacities of C¹⁸O (1–0) and C¹⁸O (2–1) become higher than unity in most regions. The rotational temperature is determined to be 4.6–7.6 K, while the C¹⁸O column densities range from 9.1 × 10¹⁴ cm⁻² to 3.0 × 10¹⁵ cm⁻². The C¹⁸O fractional abundances relative to H₂ range from 5.8 × 10⁻⁸ to 3.5 × 10⁻⁷, confirming CO depletion in this filament.

• 3.

  Based on the velocity centroid map, we derive local velocity gradients toward this filament. Velocity gradients decrease toward the filament’s crest while the position angles are randomly distributed at a radius of <0.1 pc and become ordered at a radius of ~0.1 pc. We performed Monte Carlo simulations to confirm that the local velocity gradients have a tendency to be perpendicular to the filament’s long axis in the outskirts. Velocity gradient vectors with relatively high magnitudes are wrapping the filament’s southern end, providing additional observational support for the edge effect of a finite filament.

• 4.

  We investigate the longitudinal and radial profile of physical and chemical parameters of the Serpens filament. The southeastern part of the filament can be regarded approximately as a finite isothermal and quiescent filament, while its northwestern part is likely affected by the outflow feedback of Ser-emb28. Unlike the dust temperature, the C¹⁸O rotational temperature increases toward the filament’s crest. At a radius of ≲0.1 pc, the nonthermal velocity dispersion is nearly sonic, less turbulent than the region farther out. The local velocity gradients are lower and more randomly-oriented at r ≲ 0.1 pc than at r > 0.1 pc.

• 5.

  We utilize a pseudo-time dependent chemical model in order to model the CO depletion and constrain the chemical timescale of the Serpens filament. We obtain a timescale of ≲2 Myr, which is roughly comparable to the accretion timescale. This suggests that the filament is young and still accreting. The isochronic evolutionary track of the C¹⁸O freeze-out process may serve as an independent tool to evaluate the evolution of filaments.

On the basis of these results, we propose that the Serpens filament is a finite filament formed in a sheet-like molecular cloud, with a filamentary density enhancement first generated by turbulence or colliding flows, and then followed by gravity-driven accretion of ambient gas. Higher angular resolution observations of magnetic fields will be needed to further assess the role of the magnetic field which may affect the observed gravity-driven accretion flows.

Acknowledgments

sWe wish to thank the IRAM and APEX staff for their assistance with our observations. This work was supported by the National Key R&D Program of China under grant 2017YFA0402702. The research leading to these results has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No 730562 [RadioNet]. This work is based on observations carried out under project numbers 102-17 and 150-18 with the IRAM 30m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX). APEX is a collaboration between the Max-Planck-Institut fur Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. Y.G. thanks Henrik Beuther, Philippe André, and Alvaro Hacar for helpful discussions during the APEX2020 meeting. Y.G. thanks Michael Rugel and Aiyuan Yang for organizing the treasure hunt in Valle de la Luna during the APEX observations. Y.G. also thanks his wife, Wenjin Yang, for her suggestion to use the color style of Piet Mondrian in this paper. F.Du is financially supported by the Hundred Talents Program of Chinese Academy of Sciences and the National Natural Science Foundation of China (Project number 11873094). G.X.L. is supported by Yunnan University Grant No. C176220100028. This research has made use of NASA’s Astrophysics Data System. This work also made use of Python libraries including Astropy (https://www.astropy.org/; Astropy Collaboration 2013), NumPy (https://www.numpy.org/; van der Walt et al. 2011), SciPy (https://www.scipy.org/; Jones et al. 2001), Matplotlib (https://matplotlib.org/; Hunter 2007), APLpy (Robitaille & Bressert 2012), seaborn (https://seaborn.pydata.org/), Glue (Beaumont et al. 2015; Robitaille et al. 2017), emcee (Foreman-Mackey et al. 2013), and corner.py (Foreman-Mackey 2016). The Overleaf online editor (www.overleaf.com) was used to prepare this manuscript. A 3D view of the C¹⁸O (2–1) position-position-velocity cube can be found in https://gongyan2444.github.io/gallery.html. We thank the anonymous referee for the insightful comments that improve this manuscript.

Appendix A Zoom-in self-absorption spectra

Here, we present a closer view onto the ¹³CO and C¹⁸O spectra toward two positions in Fig. A.1. Absorption dips are visible in the two lowest-J ¹³CO transitions. The dips appear at 7.92 km s⁻¹ at the offset (−80′′, 130′′) and 8.20 km s⁻¹ at the offset (0′′, −120′′). The velocity difference corresponds to 8 channels in the ¹³CO (1–0) spectra and 4 channels in the ¹³CO (2–1) spectra. This demonstrates that the absorption features are not due to contamination in the off position because such a contamination would cause the dip to occur at the same LSR velocity toward all observed positions.

Fig. A.1 Detailed view on the spectra toward two selected positions, the offsets of which are indicated in the upper left corner. For the origin of the coordinate system, see the caption to Fig. 2. In each panel, the dashed line marks the velocity of the absorption dip seen in the 13CO spectra, highlighting that the velocity of the dip varies with position.

Appendix B Additional material for the identified molecular outflow

Here, we present the channel map of the molecular outflow driven by Ser-emb 28 which is shown in Fig. B.1. Based on this map, we are able to distinguish the outflow and ambient cloud emission at velocities of v_lsr ≤5.5 km s⁻¹ and v_lsr ≥9.4 km s⁻¹. We notice that the compact emission in the range of 8.8–9.3 km s⁻¹ can also arise from the molecular outflow. We investigated the integrated intensity profile (see Fig. B.2a) in the integrated intensity map (see Fig. B.2b). The intensity profiles suggest the presence of a compact component superimposed over extended ambient emission. We found that the ambient emission is characterized by a nearly constant integrated intensity of 0.24 K km s⁻¹ for each pixel. Subtracting the constantvalue of 0.24 K km s⁻¹ from Fig. B.2b, we obtained the residual map shown in Fig. B.2c whichwas then used to estimate the outflow contribution in the velocity range of 8.8–9.3 km s⁻¹. The results are shown as the redshifted-2 component in Table 2. We found that these values are not negligible in terms of mass, momentum, and kinetic energy, so its contribution is also taken into account when calculating the overall outflow properties.

Fig. B.1 13CO (2–1) channel map of the molecular outflow driven by Ser-emb 28. The LSR velocity is indicated in the top left corner of each panel. The contours start at 0.2 K (3σ), and increase by factors of 2. The offset (0,0) corresponds to the position of Ser-emb 28. The displayed area is the same as Fig. 4b.

Fig. B.2 (a) Integrated intensity profiles across the two cuts indicated by the two dashed lines in Fig. B.2b. The offsets increase from south to north for the red line, and from east to west for the blue line. (b) Intensity map integrated from 8.8 km s−1 to 9.4 km s−1. The contours start at 0.05 K km s−1 (3σ), and increase by 0.05 K km s−1. (c) Ambient-cloud-emission-subtracted integrated intensity map.

Appendix C Nonthermal velocity dispersions with opacity correction

As mentioned in Sect. 3.3, the observed C¹⁸O line widths are affected by opacity broadening, while the signal-to-noise ratios (S/N) for C¹⁷O lines are relatively low. In order to derive the intrinsic velocity dispersion, we need to correct the C¹⁸O line width for the opacity broadening effects. Following previous studies (e.g., Phillips et al. 1979; Hacar et al. 2016b), the opacity broadening effect can be expressed as: $$\frac{\Delta v}{\Delta v_{\text{int}}}=\frac{1}{\sqrt{\text{ln2}}}{\left[\text{ln}\left(\frac{{\tau }_0}{\text{ln}\left(\frac{2}{\text{exp}(-{\tau }_0)+1}\right)}\right)\right]}^{1/2},$$(C.1)

where Δv and Δv_int are the observed and intrinsic line widths, respectively. Based on our results shown in Figs. 5 and 9, we find that the ratio $\frac{\Delta v}{\Delta v_{\text{int}}}$ can become as high as 1.4. Subtracting the thermal contributions from the intrinsic line widths, we can derive the opacity-corrected nonthermal velocity dispersion with the relationship: $${\sigma }_{\text{nt},\text{c}}=\sqrt{\frac{\Delta v_{\text{int}}^2}{8\text{ln}2}-\frac{\text{k}{T}_{\text{kin}}}{m_{\text{i}}}},$$(C.2)

where $\sqrt{\frac{\text{k}{T}_{\text{kin}}}{m_{\text{i}}}}$ is the thermal velocity dispersion, σ_th; m_i is the C¹⁸O mass weight; k is the Boltzmann constant; and T_kin is the kinetic temperature. We assumed the kinetic temperature to be 7 K (see Sect. 4.2), corresponding to an H₂ sound speed of c_s = 0.17 km s⁻¹ and a C¹⁸O thermal dispersion of σ_th = 0.04 km s⁻¹. Applying Eqs. (C.1)–(C.2) to our observed line widths (see Fig. 5), we derived the opacity-corrected nonthermal velocity dispersion map which is shown in Fig. C.1. The velocity dispersions outside the 2σ boundary in Fig. 9 are also calculated by neglecting the opacity broadening effects. In the innerregion, the opacity-corrected nonthermal velocity dispersions range from 0.07 to 0.40 km s⁻¹, corresponding to Mach numbers, M = σ_nt,c∕c_s, of 0.4–2.4. This suggests that the Serpens filament is more quiescent than previously thought (Gong et al. 2018). On the other hand, the region around emb28 becomes supersonic (i.e., M > 2) likely due to the outflow feedback (see Fig. 4).

Fig. C.1 Opacity-corrected nonthermal velocity dispersion maps derived from C18O (1–0) and C18O (2–1) data. The beam size is shown in the lower right corner of each panel. In both panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses show the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (Enoch et al. 2009), and the white crosses mark the positions of the seven 1.1 mm dust cores (Enoch et al. 2007).

The turbulent energy, E_t, can be estimated with $E_{\text{t}}={\Sigma }_i\frac{1}{2}{m}_i{\sigma }_{i,\text{3d}}^2$, where m_i and ${\sigma }_{i,\text{3d}}=\sqrt{3}{\sigma }_{\text{nt},\text{c}}$ are the molecular mass and three-dimenional velocity dispersion of the corresponding pixel, i. Based on Fig. C.1, we obtain a total turbulent energy of 1.15 × 10⁴⁴ erg, where the turbulent energies of NW and SE are 7.3 × 10⁴³ erg and 4.2 × 10⁴³ erg, respectively. The difference in the turbulent energies of NW and SE is 3.1 × 10⁴³ erg, which is about six times the outflow kinetic energy from Ser-emb 28 (see Table 2). Therefore, the different levels of turbulence in NW and SE cannot be explained only by the outflow feedback.

The turbulent dissipation rate, Ė_t, is defined as, $${\dot{E}}_{\text{t}}=E_{\text{t}}/t_{\text{diss}}$$(C.3)

where t_diss is the turbulent dissipation time. Based on McKee & Ostriker (2007), the turbulent dissipation time is approximately given by $$t_{\text{diss}}\approx 0.5\frac{d}{{\sigma }_{\text{los}}}$$(C.4)

where d is the cloud diameter and σ_los is the turbulent velocity dispersion along the line of sight. We adopt the filament’s FWHM length as the cloud diameter, that is, d = 0.17 pc, and use the mean σ_nt,c of 0.2 km s⁻¹ as σ_los. We thus obtain the dissipation time of 0.42 Myr. Based on Eq. (C.3), we find that Ė_t is about 8.8 × 10³⁰ erg s⁻¹ which is comparable to the outflow mechanical luminosity of 5.6 × 10³⁰ erg s⁻¹. This suggests that the outflow feedback could sustain the observed turbulent motions.

Appendix D Overlapping hyperfine structure lines of C¹⁷O

C¹⁷O transitions show hyperfine structure (HFS) splitting because of the nonzero nuclear spin in ¹⁷O. The broadening introduced by velocity dispersion and eventually also by opacity can cause the blending of their HFS components. This will introduce additional uncertainties when we estimate the peak optical depth with Eq. (1) by assuming that a single component or a combined component make the dominant contributions to the peak.

In order to study this impact, we compute the synthetic C¹⁷O (1–0) and C¹⁷O (2–1) spectra with Eq. (1) (see Sect. 3.4) to compare the intrinsic peak intensity of a single component or a combined component with the peak intensity of the blended spectra. Assuming a Maxwellian velocity distribution, the opacity can be expressed by a Gaussian distribution: $$\tau (v)={\tau }_0{\text{e}}^{-{(v-v_0)}^2/2{\sigma }_v^2},$$(D.1)

where τ₀ and v₀ are the opacity and LSR velocity at the line center, and σ_v is the intrinsic velocity dispersion. For the different HFS lines, their opacities are proportional to their relative line strengths assuming LTE (Müller et al. 2005). Based on the rest frequency of HFS components (Müller et al. 2005), we are able to obtain their displaced LSR velocities. σ_v should be the same for different components. Taking C¹⁷O (1–0) for example, the F = 3/2–5/2, F = 7/2–5/2, and F = 5/2–5/2 lines have opacities of 2/9, 4/9, and 1/3 of the total opacity, while the F = 3/2–5/2 and F = 5/2–5/2 lines are displaced by 0.5469741 km s⁻¹ and −2.7348705 km s⁻¹ with respect to F = 7/2–5/2. The observed line widths are determined to be 0.25–0.89 km s⁻¹ with a median of 0.52 km s⁻¹ (see Fig. 5), suggesting that σ_v should be lower than 0.38 km s⁻¹ in the Serpens filament. We therefore use a σ_v range of 0.10–0.38 km s⁻¹ for our calculations. Since both C¹⁷O lines are likely to have low opacities, we use a peak total opacity of τ₀ ≤1 for our test. We are thus able to construct the synthetic C¹⁷O (1–0) and C¹⁷O (2–1) spectra with Eq. (1) where v₀ is assumed to be 0 and T_ex is set to be 7 K as the fiducial case.

The synthetic C¹⁷O (1–0) and C¹⁷O (2–1) spectra are shown in Fig. D.1. Because C¹⁷O (1–0) and C¹⁷O (2–1) are dominated by the F = 7/2–5/2 feature (see Figs. D.1a–b) and the Group I line (see the Group I line in Figs. D.1c–d), we use the ratio of the peak intensities of these lines and their corresponding blended lines to assess the difference due to the overlapping HFS lines. We find the ratio increases with increasing τ₀ and σ_v for C¹⁷O (1–0). The highest ratio becomes 0.78 when we use the highest values (τ₀ = 1, σ_v = 0.38 km s⁻¹). This means that the peak intensity can be overestimated by 28% in this extreme case. For C¹⁷O (2–1), the ratio does not change much, and the peak intensity is only overestimated by 3% in the extreme case. We also note that the ratio is independent of T_ex, suggesting that the ratios are robust within our modeling ranges. Therefore, we introduce additional 28 and 3% uncertainties for the peak intensities of C¹⁷O (1–0) and C¹⁷O (2–1) in studying molecular excitation in Sect. 3.4.

Fig. D.1 Synthetic C17O (1–0) and C17O (2–1) spectra with the modeled τ0 and σv shown in the upper left corners. For C17O (2–1), Group I contains four HFS lines, F = 5/2–5/2, F = 7/2–5/2, F = 9/2–7/2, and F = 1/2–3/2, while Group II contains three HFS lines, F = 5/2–7/7, F = 7/2–7/2, and F = 5/2–3/2.

References

All Tables

All Figures

Fig. 1 Three-color composite image of the Serpens filament (Herschel 70 μm: blue, Herschel 250 μm: green, Herschel 500 μm: red). The three green pluses indicate the positions of the three embedded YSOs, emb10 (also known as IRAS 18262+0050), emb16, and emb28 (Enoch et al. 2009), while the crosses mark the positions of 1.1 mm dust cores (Enoch et al. 2007). The beam size (36.′′3) for the 500 μm emission is shown in the lower right corner. The regions mapped with the IRAM-30 m and APEX telescopes are indicated by the red and white dashed boxes, respectively.

In the text

Fig. 2 Integrated intensity maps of 13CO (1–0) (a), C18O (1–0) (b), C17O (1–0) (c), 13CO (2–1) (d), C18O (2–1) (e), and C17O (2–1) (f). Their integrated velocity ranges are [6.5, 10] km s−1 for 13CO (1–0), C18O (1–0), 13CO (2–1), and C18O (2–1), [7.5, 9] km s−1 for C17O(1–0), and [6, 10] km s−1 for C17O (2–1). The contours start at 0.4 K km s−1 with a step of 0.4 K km s−1 for 13CO (1–0), C18O (1–0), 13CO (2–1), and C18O (2–1). For the two C17O maps, the contours start at 3σ with a step of 3σ, where 1σ = 0.08 K km s−1 for C17O (1–0) and 1σ = 0.06 K km s−1 for C17O (2–1). The color bar represents the C18O (1–0) and C18O (2–1) integrated intensities in units of K km s−1. The 13CO (1–0), 13CO (2–1), C17O (1–0), and C17O (2–1) maps have been scaled by 0.4, 0.5, 3, and 3 to match the color bar. Panel d: the red dashed line is used to divide the Serpens filament into two parts, NW and SE, and the blue line represents the crest of the Serpens filament extracted in the corresponding H2 column density map (Gong et al. 2018). In all panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses indicate the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (see Enoch et al. 2009, and Fig. 1), and the white crosses mark the positions of the seven dust cores (see Enoch et al. 2007, and Fig. 1).

In the text

	Fig. 3 Observed 13CO (2–1), C18O (2–1), and C17O (2–1) spectra toward the selected positions that are indicated in the upper left corner. For the origin of the coordinate system, see the caption to Fig. 2. 13CO self-absorption is evident in these spectra.
In the text

Fig. 4 Outflow driven by Ser-emb 28. Left: same as Fig. 2e. Right: (a) observed 13CO (2–1) spectra of the two positions indicated by the blue and red pluses in panel b overlaid with the C18O spectrum of the position indicated by the yellow star in panel b. (b) 13CO (2–1) outflow map of Ser-emb 28. The blueshifted emission integrated from 2.0 to 5.5 km s−1 is shown in blue dashed contours which start at 0.20 K km s−1 (5σ) and increase by 0.20 K km s−1. The redshifted emission integrated from 9.4 to 12 km s−1 is shown in red dashed contours which start at 0.20 K km s−1 and increase by 0.20 K km s−1. The yellow star marks the position of Ser-emb 28. The beam size is shown in the lower left corner. (c) Position-velocity diagram of the 13CO (2–1) emission along the cut indicated by the black line in panel b. The contours start at 0.18 K (3σ) and increase by 0.36 K. The horizontal dashed line indicates the position of Ser-emb 28, and the vertical dashed line marks its systemic velocity of 7.87 km s−1. The resolution element is shown in the lower left corner.

In the text

Fig. 5 Maps of peak intensities (a), LSR velocities (b), and line widths (c) derived from single-component Gaussian fits to our C18O (1–0) data. Panels d–f are similar to panels a–c, but for C18O (2–1). In panels a and d, the contours represent the peak intensities which start at 0.6 K and increase by 0.6 K, while in c and f the contours represent the H2 column densities that start at 1 × 1021 cm−2 and increase by 3 × 1021 cm−2. In panels a and d, the dotted lines indicate possible substructures. The beam size is shown in the lower right corner of each panel. In all panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses indicate the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (Enoch et al. 2009), and the white crosses mark the positions of the seven dust cores (Enoch et al. 2007).

In the text

	Fig. 6 Similar to Fig. 5, but for C17O (1–0) and C17O (2–1). In panels a and d, the contours represent the peak intensities which start at 0.3 K and increase by 0.3 K. The pixels with large fitting uncertainties are masked out in panels b–c and e–f.
In the text

	Fig. 7 C18O (1–0), C17O (1–0), C18O (2–1), and C17O (2–1) spectra averaged over the low-brightness-temperature mask (details are discussed in Sect. 3.4). The two C17O spectra have been scaled by a factor of 4. The hyperfine structure components of both C17O transitions are indicated by vertical dashed lines.
In the text

	Fig. 8 Posterior probability distributions of C18O column density, $N_{{\text{C}}^{18}\text{O}}$, rotational temperature, Trot, the C18O (1–0) opacity, τ18,1, and the C18O (2–1) opacity, τ18,2, toward the dense core, Bolo12, with the maximum posterior possibility point in the parameterspace shown in orange lines and points. Contours denote the 0.5, 1.0, 1.5, and 2.0σ confidence intervals.
In the text

Fig. 9 Optical depth maps of C18O (1–0) (a) and C18O (2–1) (b) overlaid with contours of signal-to-noise ratio (τ∕στ) of 2 and 3 indicated by the red and black lines, respectively. The beam size is shown in the lower right corner of each panel. In all panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses show the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (Enoch et al. 2009), and the white crosses mark the positions of the seven dust cores (Enoch et al. 2007).

In the text

	Fig. 10 Histogram and cumulative distribution of the C18O (1–0) (a) and C18O (2–1) (b) opacities.
In the text

Fig. 11 (a) C18O rotational temperature map. (b) C18O column density map. (c) C18O fractional abundance map overlaid with H2 column density contours. The contours start at 6 × 1021 cm−2, and increase by 3 × 1021 cm−2. The beam size is shown in the lower right corner of each panel. In all panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses indicate the positions of the three embedded YSOs, and the white crosses mark the positions of the seven dust cores (Enoch et al. 2007).

In the text

Fig. 12 (a) Dust temperature as a function of the C18O rotational temperature. The gray error bars represent the measured uncertainties in the measured dust temperature and rotational temperature. The markers’ colors correspond to their H2 column densities indicated by the color bar. The Pearson correlation coefficient is shown in the lower left. (b) C18O rotational temperature as a function of the H2 column density. The gray error bars represent the measured uncertainties in the measured rotational temperature and H2 column density. The markers’ colors correspond to their dust temperatures indicated by the color bar. The Pearson correlation coefficient is shown in the lower right. In both panels, the black dashed line represents the linear fit to the observed trend.

In the text

Fig. 13 Left: C18O fractional abundance map with different labeled masks which correspond to the legend of the right panel. Right: C18O fractional abundance as a function of H2 column density. The gray error bars indicate the measured uncertainties in C18O fractional abundance and H2 column density. The colors of the markers correspond to different regions indicated in the legend. The Pearson correlation coefficient isshown in the lower left. The black dashed line represents the linear fit to the observed relation.

In the text

Fig. 14 (a) Local velocity gradient magnitude map overlaid with the H2 column density contours. (b) Herschel H2 column density map overlaid with the normalized velocity vector maps. The arrows represent the estimated local velocity gradients which are rotated by 180° in order to better visualize the accretion directions in SE. The polarized angle of the Planck 353 GHz thermal dust emission has been rotated by 90° to trace the magnetic field direction which is indicated by the red line. In both panels, the contours correspond to H2 column densities from 1 × 1021 cm−2 to 1.6 × 1022 cm−2 with a step of 3 × 1021 cm−2, and the beam sizes for the H2 column density and velocity gradient maps are shown as the gray and blue circles in the lower right corner. The three green plusesgive the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (Enoch et al. 2009), and the white crosses mark the positions of the seven dust cores (Enoch et al. 2007).

In the text

Fig. 15 (a) Histograms of the angle of the vector of the observed local velocity gradient relative to the long axis of the filament. The blue histograms show all the angles of the observed velocity gradients in SE, while the yellow histograms represent the velocity gradients with |∇V | >1.5 km s−1 pc−1. (b) Cumulative distribution function of the angle between the velocity gradients and the major axis of the filament. The solid blue and yellow lines represent the observed cumulative distribution function of relative angles with corresponding and |∇V | >1.5 km s−1 pc−1 and |∇V | >1.5 km s−1 pc−1, respectively. The three blue dashed lines present the unimodal simulations of the expected angles which are three-dimensionally parallel (0–20°), perpendicular (70–90°), or random (0–90°), while the red dashed lines show the bimodal simulations of the angles that are parallel and perpendicular with increments of 10% (i.e., the top red line is 90% parallel and 10% perpendicular, the next one is 80% parallel and 20% perpendicular).

In the text

Fig. 16 Physical and chemical parameters (from top to bottom: position angle of local velocity gradient, magnitude of local velocity gradient, C18O rotational temperature, C18O column density, C18O fractional abundance, dust temperature, H2 column density, nonthermal velocity dispersion, and velocity centroid) along the crest of the Serpens filament indicated by the blue line in Fig. 2d. The origin of the offsets corresponds to the southern end of the blue line in Fig. 2d. The offset increases from south to north. The positions which are close to the seven 1.1 mm dust cores and the three YSOs are indicated by the cyan and orange vertical lines, respectively. In the two lowest panels, the data derived from C18O (1–0) and C18O (2–1) are indicated by blue and red symbols, respectively.

In the text

Fig. 17 C18O rotational temperature (a), dust temperature (b), H2 column density (c), C18O column density (d), C18O fractional abundance (e), velocity centroid (f), nonthermal velocity dispersion (g), magnitude of velocity gradient (h), and position angle of velocity gradient (i) plotted as a function of their distance to the crest of SE. In c, the yellow shaded region marks the star formation threshold of Av = 7–10 mag (e.g., Johnstone et al. 2004; Lada et al. 2010; Heiderman et al. 2010; André et al. 2010; Könyves et al. 2015). In all panels, the radius of 0.1 pc is marked with a red dashed line, which indicates the transition discussed in Sect. 4.3. We note that the data points become incomplete at larger radial distances in panels a, d, e, and f due to the 3σ clipping mentioned in Sect. 3.4. In panel g, the two red horizontal dashed lines mark σt = cs and σt = 2cs where cs is the soundspeed at a kinetic temperature of 7 K.

In the text

Fig. 18 (a) C18O fractional abundances, ${\chi }_{{\text{C}}^{18}\text{O}}$, as a function of time, as calculated with different H2 column densities and an assumed thickness of 0.68 pc. The shaded region corresponds to the observed C18O fractional abundance range. (b) Similar to Fig. 13, but overlaid with the isochronic lines which are based on the results in a. The markers represent the observed relationship between C18O fractional abundance and H2 column density with different regions labeled in different colors. c and e are similar to a, but for an assumed thickness of 0.34 pc and 1 pc, respectively. Panels d and f are similar to panel b, but for an assumed thickness of 0.34 pc and 1 pc, respectively.

In the text

Fig. 19 Schematic view of the geometry of the Serpens filament. The left panel gives a perspective view of its 3D geometry, and the right panels represent the projection on to the x-y, x-z, and y-z planes. The midplane and its parallel surfaces have the same distances to observers, and the line of sight is along the z direction. The local velocity gradients that are associated with blueshifted and redshifted flows are indicated by the blue and red arrows, while the magnetic field is indicated bygreen dashed lines.

In the text

	Fig. A.1 Detailed view on the spectra toward two selected positions, the offsets of which are indicated in the upper left corner. For the origin of the coordinate system, see the caption to Fig. 2. In each panel, the dashed line marks the velocity of the absorption dip seen in the 13CO spectra, highlighting that the velocity of the dip varies with position.
In the text

	Fig. B.1 13CO (2–1) channel map of the molecular outflow driven by Ser-emb 28. The LSR velocity is indicated in the top left corner of each panel. The contours start at 0.2 K (3σ), and increase by factors of 2. The offset (0,0) corresponds to the position of Ser-emb 28. The displayed area is the same as Fig. 4b.
In the text

Fig. B.2 (a) Integrated intensity profiles across the two cuts indicated by the two dashed lines in Fig. B.2b. The offsets increase from south to north for the red line, and from east to west for the blue line. (b) Intensity map integrated from 8.8 km s−1 to 9.4 km s−1. The contours start at 0.05 K km s−1 (3σ), and increase by 0.05 K km s−1. (c) Ambient-cloud-emission-subtracted integrated intensity map.

In the text

Fig. C.1 Opacity-corrected nonthermal velocity dispersion maps derived from C18O (1–0) and C18O (2–1) data. The beam size is shown in the lower right corner of each panel. In both panels, the (0, 0) offset corresponds to αJ2000 = 18h28m50.s4, δJ2000 = 00°49′58.′′72, the three green pluses show the positions of the three embedded YSOs, emb10, emb16, and Ser-emb 28 (Enoch et al. 2009), and the white crosses mark the positions of the seven 1.1 mm dust cores (Enoch et al. 2007).

In the text

	Fig. D.1 Synthetic C17O (1–0) and C17O (2–1) spectra with the modeled τ0 and σv shown in the upper left corners. For C17O (2–1), Group I contains four HFS lines, F = 5/2–5/2, F = 7/2–5/2, F = 9/2–7/2, and F = 1/2–3/2, while Group II contains three HFS lines, F = 5/2–7/7, F = 7/2–7/2, and F = 5/2–3/2.
In the text