Physical and chemical structure of the Serpens filament: Fast formation and gravity-driven accretion (Corrigendum)

Context. The Serpens filament, a prominent elongated structure in a relatively nearby molecular cloud, is believed to be at an early evolutionary stage, so studying its physical and chemical properties can shed light on filament formation and early evolution.
Aims. The main goal is to address the physical and chemical properties as well as the dynamical state of the Serpens filament at a spatial resolution of ~0.07 pc and a spectral resolution of ≲0.1 km s−1.
Methods. We performed 13CO (1–0), C18O (1–0), C17O (1–0), 13CO (2–1), C18O (2–1), and C17O (2–1) imaging observations toward the Serpens filament with the Institut de Radioastronomie Millimétrique 30-m and Atacama Pathfinder EXperiment telescopes.
Results. Widespread narrow 13CO (2–1) self-absorption is observed in this filament, causing the 13CO morphology to be different from the filamentary structure traced by C18O and C17O. Our excitation analysis suggests that the opacities of C18O transitions become higher than unity in most regions, and this analysis confirms the presence of widespread CO depletion. Further we show that the local velocity gradients have a tendency to be perpendicular to the filament’s long axis in the outskirts and parallel to the large-scale magnetic field direction. The magnitudes of the local velocity gradients decrease toward the filament’s crest. The observed velocity structure can be a result of gravity-driven accretion flows. The isochronic evolutionary track of the C18O freeze-out process indicates the filament is young with an age of ≲2 Myr.
Conclusions. We propose that the Serpens filament is a newly-formed slightly-supercritical structure which appears to be actively accreting material from its ambient gas.


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. 2013Li et al. , 2016Schisano 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 5 M in mass (e.g., Schneider et al. 2010;Hacar et al. 2013Hacar et al. , 2018Li et al. 2013Li et al. , 2016Kainulainen 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) The reduced datacubes are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http: //cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/646/A170 function are consistent with a Salpeter-like mass function, indicating a possible link to the initial mass function . The core and star formation occurring in filaments may also cause its nearly constant efficiency in dense gas (∼4.5 × 10 −8 yr −1 , 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. 2010André et al. , 2014Kö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. 2013Hacar et al. , 2017Hacar et al. , 2018Shimajiri 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-forming seeds) 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;A&A 646, A170 (2021) 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. 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. 2017Zucker 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. 2019Su et al. , 2020. This filament is found to be at the onset of slightly supercritical collapse with a line mass of 36-41 M pc −1 and H 2 column densities of >6 × 10 21 cm −2 (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(Enoch et al. , 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.

IRAM-30 m and APEX-12 m observations
We carried out 13 CO (1-0), C 18 O (1-0), and C 17 O (1-0) observations with the Institut de Radioastronomie Millimétrique 30-m (IRAM-30 m 1 ) 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 −1 for C 18 O (1-0). Each FFTS provides spectral channels of 50 kHz, translating to a velocity spacing of 0.14 km s −1 at 112 GHz.
Observations of 13 CO (2-1), C 18 O (2-1), and C 17 O (2-1) were performed with the Atacama Pathfinder EXperiment 12 meter submillimeter telescope (APEX 2 ; 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 −1 for the C 18 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 h 27 m 42. s 72, δ J2000 = 00 • 44 54. 49) was checked to be emission-free in 13 CO (1-0) with a 1σ rms noise level of 96 mK at a channel width of 0.14 km s −1 .
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 * A , to the main-beam brightness temperature, T mb , with the relation T mb = T * A η f /η 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 3 and APEX 4 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 5 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, δ , are given in Table 1. Velocities are given with respect to the local standard of rest (LSR) throughout this paper.

Archival data
We make use of the dust temperature and H 2 column density maps from the Herschel Gould Belt (GB) Survey 6 (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 2 column density are assumed in our analysis.

Overall distribution
Figure 2 presents maps of the integrated intensity of the emission in the 13 CO (1-0), C 18 O (1-0), C 17 O (1-0), 13 CO (2-1), C 18 O (2-1), and C 17 O (2-1) lines resulting from our observations. They provide a factor of two finer linear resolution than our previous C 18 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 18 O and C 17 O line emission shows similar distributions, which trace the filamentary structure well. In contrast, the morphology of the emission of the two 13 CO lines is markedly different. Much of the filamentary structure is not seen in both 13 CO transitions, in line with previous observations (Burleigh et al. 2013). Also, 13 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 18 O (2-1) is generally brighter than 13 CO (2-1) near the systemic velocity of ∼8 km s −1 . This is likely because of significant 13 CO self-absorption, which is reinforced by the presence of the extremely narrow dip in 13 CO spectra at velocities where C 18 O and C 17 O peak. Such narrow 5 http://www.iram.fr/IRAMFR/GILDAS/ 6 http://www.herschel.fr/cea/gouldbelt/en/index.php 13 CO (2-1) self-absorption features were not reported in previous 13 CO (2-1) observations (Burleigh et al. 2013), probably due to the coarser spectral resolution of these observations. The widespread 13 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.

A new molecular outflow driven by Ser-emb 28
Based on our APEX 13 CO (2-1) observations which are more sensitive than our 13 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 −1 and redder than 9.4 km s −1 , 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 13 CO (2-1) self-absorption is observed around the systemic cloud velocity derived from C 18 O emission.
The velocity interval of the outflow emission is determined by manually investigating the 13 CO (2-1) channel map (see Appendix B). Assuming that emission in the high velocity wings of the 13 CO (2-1) line is optically thin and is in local thermodynamic equilibrium (LTE) at an excitation temperature of 20 K, we can derive 13 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 13 CO column densities. In order to derive the total molecular mass of the outflow lobes, the 13 CO fractional abundance relative to H 2 is assumed to be 1.1 × 10 −6 (e.g., Cabrit & Bertout 1992;Beuther et al. 2002)  Their integrated velocity ranges are [6.5, 10] km s −1 for 13 CO (1-0), C 18 O (1-0), 13 CO (2-1), and C 18 O (2-1), [7.5, 9] km s −1 for C 17 O(1-0), and [6, 10] km s −1 for C 17 O (2-1). The contours start at 0.4 K km s −1 with a step of 0.4 K km s −1 for 13 CO (1-0), C 18 O (1-0), 13 CO (2-1), and C 18 O (2-1). For the two C 17 O maps, the contours start at 3σ with a step of 3σ, where 1σ = 0.08 K km s −1 for C 17 O (1-0) and 1σ = 0.06 K km s −1 for C 17 O (2-1). The color bar represents the C 18 O (1-0) and C 18 O (2-1) integrated intensities in units of K km s −1 . The 13 CO (1-0), 13 CO (2-1), C 17 O (1-0), and C 17 O (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 H 2 column density map (Gong et al. 2018). In all panels, the (0, 0) offset corresponds to α J2000 = 18 h 28 m 50. s 4, δ 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). 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.

Decomposition
Since 13 CO spectra suffer from self-absorption (see Sect. 3.1), we only use the C 17 O and C 18 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 18 O (1-0) and C 18 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 17 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   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 18 O (1-0), 0.09-3.85 for C 18 O (2-1), 0.08-0.30 for C 17 O (1-0), and 0.08-0.38 for C 17 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 A170, page 5 of 25 A&A 646, A170 (2021) 4.0 × 10 −6 1.0 × 10 −3 Redshifted-2 ( * ) [8.8, 9.4] 0.09 0.009 0.01 2.5 × 10 41 3.8 × 10 −7 6.6 × 10 −7 9.0 × 10 −5 Total 0.049 0.16 5.4 × 10 42 1.9 × 10 −6 5.5 × 10 −6 1.7 × 10 −3 Notes.  Panels d-f are similar to panels a-c, but for C 18 O (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 H 2 column densities that start at 1 × 10 21 cm −2 and increase by 3 × 10 21 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 = 18 h 28 m 50. s 4, δ 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).
perpendicular to the filament, indicated by the dotted lines in Figs. 5a and d. These features may be indicative of substructures that provide mass accretion Palmeirim et al. 2013). However, these features are not evident in the C 18 O integrated intensity and H 2 column density maps, probably due to their low contrast.  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. an opposite trend with lsr = 7.8-8.1 km s −1 in its inner part and lsr = 7.5-7.9 km s −1 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 18 O (1-0) at a linear scale of 0.13 pc ( LSR ∼ 8.6 km s −1 , Gong et al. 2018).
Figures 5c, f, 6c, and f exhibit line width maps derived from the C 18 O (1-0), C 18 O (2-1), C 17 O (1-0), and C 17 O (2-1) data. Due to their lower signal-to-noise ratios, the fitting results for C 17 O are not as good as those for C 18 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 −1 , 0.02 km s −1 , 0.05 km s −1 , and 0.04 km s −1 for the fitted line widths of C 18 O (1-0), C 18 O (2-1), C 17 O (1-0), and C 17 O (2-1), respectively. We see that the C 18 O (1-0) and C 18 O (2-1) lines are broader than the C 17 O lines in a significant number of pixels, despite the relatively large uncertainties in the C 17 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 18 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 −1 in the inner region of SE, corresponding to σ nt <0.25 km s −1 at an assumed kinetic temperature of 7 K. Higher line widths of >0.6 km s −1 (σ nt >0.25 km s −1 ) are observed in NW and the outer region (N H2 10 22 cm −2 ) 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).

Molecular excitation
As mentioned above, opacity effects should be taken into account in the analysis. According to the radiative transfer A170, page 7 of 25 A&A 646, A170 (2021) equation, the main beam temperature, T mb , can be expressed as 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 17 O transitions are affected by HFS splitting, we used the F = 7/2-5/2 component of C 17 O (1-0) and the combined component, Group I (see details in Appendix D) of C 17 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 17 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 17 O (1-0) and C 17 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 τ = τ 1−e −τ , the level populations result in the following relation (e.g., Goldsmith & Langer 1999): 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 µ 2 S from the CDMS (Müller et al. 2005). We also assume a homogeneous [ 18 O/ 17 O] isotopic ratio, r O , in this filament. This gives us the following relations: where N 18 and N 17 are the molecular column densities of C 18 O and C 17 O, and τ 18,i and τ 17,i are the peak optical depths of the corresponding C 18 O and HFS-corrected C 17 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 [ 18 O/ 17 O] isotopic ratio in low brightness regions. For this, we created a mask where the C 18 O (1-0) and C 18 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 18 O and C 17 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  (Wouterloot et al. 2005(Wouterloot et al. , 2008Zhang et al. 2007Zhang et al. , 2020Giannetti 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 18 O and C 17 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 2 column density maps, we first convolved C 18 O and C 17 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 7 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 C 18 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 ∝ e −χ 2 /2 with 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 18 O (1-0) and C 18 O (2-1) are shown in Fig. 9. The C 18 O (1-0) opacities lie in the range of 0.3-2.9, while the C 18 O (2-1) opacities are 0.4-2.5. It is surprising that we found more than 80% of pixels with the C 18 O (1-0) and C 18 O (2-1) optical depths higher than unity (see Fig. 10), because these C 18 O lines are usually assumed to be optically thin in the literature (e.g., Ikeda & Kitamura 2009;Gong et al. 2016Gong et al. , 2018. The highest opacities are found around Bolo6. Such high C 18 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 2 density of >10 4.5 cm −3 due to the gas-dust coupling, while the gas temperature becomes much higher than the dust temperature at lower H 2 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 18 O and C 17 O transitions become sub-thermally excited in low-density regions (n H2 10 3 cm −3 ), 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 2 column densities and high rotational temperatures are found to be associated with high H 2 column densities (see Fig. 12b). We performed a linear fit to Fig. 12b, and obtain an empirical relation of T rot (K) = (1.04 ± 0.03) N(H 2 ) 1 × 10 22 cm −2 + (5.17 ± 0.03). The fact that the rotational temperature is lower than the dust temperature in high H 2 column density regions can A170, page 9 of 25 A&A 646, A170 (2021) Fig. 9. Optical depth maps of C 18 O (1-0) (a) and C 18 O (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 = 18 h 28 m 50. s 4, δ 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). 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 18 O column densities 8 . C 18 O column densities range from 9.1 × 10 14 cm −2 to 3.0 × 10 15 cm −2 with a median of 1.7 × 10 15 cm −2 , which almost doubles the values reported in Gong et al. (2018). This is because Gong et al. (2018) did not take C 18 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 18 O column densities. The C 18 O fractional abundances relative to H 2 are directly calculated by the ratio between C 18 O and H 2 column densities. The latter is derived from the SED fitting with Herschel data (see Sect. 2.2). The distribution of the C 18 O fractional abundances is shown in Fig. 11c. The C 18 O fractional abundances range from 5 × 10 −8 to 3.5 × 10 −7 with a median of 1.6 × 10 −7 . The map shows that the C 18 O abundance drops toward the 1.1 mm dense cores. The region toward bolo3 and Ser-emb10 has the lowest abundance of ∼5 × 10 −8 . This confirms the presence of C 18 O depletion in this filament (Gong et al. 2018). Setting a C 18 O fractional abundance of 2 × 10 −7 as a threshold for C 18 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 has a larger depletion size than SE, which is likely because NW has more widespread dense gas. This reason also explains the fact that the C 18 O depletion size in NW is larger than the typical depletion size ( 0.1 pc) found in a single prestellar core (e.g., Bergin & Tafalla 2007).
A pixel-by-pixel comparison between H 2 column densities and C 18 O fractional abundances is displayed in Fig. 13  increasing H 2 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 18 O fractional abundances with increasing H 2 column densities, which gives 1 × 10 22 cm −2 + (2.67 ± 0.03). Furthermore, the lowest C 18 O abundances of ∼5 × 10 −8 are about a factor of 6 lower than those (∼3.0 × 10 −7 ) in the ambient gas. We also find that different regions show distinct distributions in Fig. 13. Intriguingly, there are two branches at H 2 column densities of >1.5 × 10 22 cm −2 . The upper branch arises from the bolo2 region, while the lower branch originates in the bolo3 region. Despite similar H 2 column densities, dust temperature, and C 18 O rotational temperature, the bolo2 region shows a factor of 2 higher C 18 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 2 number density of the   bolo3 region might be enhanced by the outflow shock, resulting in more efficient CO depletion in the bolo3 region.

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-tonoise C 18 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: where LSR is the observed LSR centroid velocity, 0 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 |∇V| = √ 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 −1 pc −1 around Bolo11 to ∼9 km s −1 pc −1 in the east of emb10. Most |∇V| are higher than the values estimated from a large-scale analysis ( 1.5 km s −1 pc −1 , 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. 2019Chen et al. , 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 filament from ambient gas. Toward the more quiescent region SE, the local velocity gradient, |∇V|, decreases from the ambient cloud (|∇V| >2 km s −1 pc −1 ) to the crest (|∇V| < 0.5 km s −1 pc −1 ). 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 −1 pc −1 , 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 −1 , 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 −1 , we also found high p-values of 0.90 against the bimodal simulations with >95% perpendicular angles. This demonstrates that Herschel H 2 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 H 2 column densities from 1 × 10 21 cm −2 to 1.6 × 10 22 cm −2 with a step of 3 × 10 21 cm −2 , and the beam sizes for the H 2 column density and velocity gradient maps are shown as the gray and blue circles in the lower right corner. The three green pluses give 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).  the velocity gradients with ∇V >1.5 km s −1 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 2 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. Figure 16 presents the C 18 O rotational temperature, C 18 O column density, C 18 O fractional abundance, dust temperature, H 2 column density, nonthermal velocity dispersion, and velocity centroid 9 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 18 O level populations are likely thermalized in the crest due to its high H 2 number density of ∼4 × 10 4 cm −3 (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 −1 pc −1 . 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., 9 We find that generally the C 18 O (1-0) line's LSR velocities are systematically ≈0.04 km s −1 higher than the velocities of the C 18 O (2-1) line. The uncertainties in the rest frequencies of the C 18 O (1-0) and C 18 O (2-1) lines are 6.3 × 10 −6 GHz and 1.5 × 10 −6 GHz, respectively; see Table 1, whose values are taken from the CDMS (Müller et al. 2005). This results in 1σ uncertainties of 0.017 km s −1 and 0.002 km s −1 in velocity units. Hence, the systematic velocity difference may be caused (largely) by the fact the rest frequency we adopted for the 1-0 line is slightly too low. Given that the two lines were observed with different telescopes, there might also be a slight inaccuracy in the frequency calibration and calculation of one of the two telescopes. 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 −1 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.

Longitudinal profiles
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 −1 pc −1 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 2 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 H 2 drop at an offset of ∼1.2 pc and two enhanced N 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 18 O rotational and dust temperatures. One should note that the elevated dust temperature and reduced H 2 column density might alternatively arise from their degeneracy during the SED fitting. As mentioned above, C 18 O molecules have been frozen out onto dust grains in most parts of the crest, but these C 18 O molecules are sublimated or enhanced by sputtering of related species in shock regions, producing the observed higher C 18 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.

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 temperature profile increases monotonically toward the crest, which is caused by the variation of the H 2 number density. The observed rotational temperature profile validates the assumption used to interpret the observed blueskewed 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.  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 C 18 O (1-0) and C 18 O (2-1) are indicated by blue and red symbols, respectively.
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) 10 . Since the velocity dispersions drop to nearly sonic levels (see Fig. 17g) at this radius, this supports the scenario that the transition to the sonic turbulence can also be the reason for the origin of the star formation threshold. Figure 17d shows an increasing N C 18 O toward the crest while Fig. 17e shows a decreasing χ C 18 O toward the crest. The former one is expected with increasing H 2 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 In c, the yellow shaded region marks the star formation threshold of A v = 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 = c s and σ t = 2c s where c s is the sound speed at a kinetic temperature of 7 K. 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: where ∆ 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 1 while the velocity structure is likely induced by self-gravity when C 1. Based on Fig. 17f, we find ∆ h ∼0.17 km s −1 for r ∼ 0.17 pc where M(r)/L was found to be 45 M pc −1 . We thus obtain C ∼ 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. 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 A170, page 16 of 25 Y. Gong et al.: Witnessing filament formation 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).

Chemical model
As mentioned in Sect. 3.4, our observations confirm the widespread C 18 O depletion, and also reveal a trend that the C 18 O fractional abundances decrease with increasing H 2 column densities. This can be attributed to the H 2 density dependence of the CO freeze-out timescale. In order to simulate the observed trend, we run a set of models using the pseudo-time dependent 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 2 (see Table 3 of 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 −17 s −1 (e.g., van der Tak & van Dishoeck 2000). The self-shielding of H 2 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 18 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 [ 16 O/ 18 O] of 530 (Wilson & Rood 1994;Bockelée-Morvan et al. 2012) to convert the 12 CO fractional abundance to the C 18 O fractional abundance. Only two parameters, visual extinction (A v ) and H 2 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 2 column density profile (Gong et al. 2018). The H 2 number density can be determined by dividing the H 2 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 2 ) = 9.4 × 10 20 (A v /mag) (Bohlin et al. 1978).
The pseudo-time dependent chemical modeling results are shown in Fig. 18a. This demonstrates that higher H 2 column density regions tend to be depleted earlier, as expected. Based on the modeling results, we utilize the isochronic tracks in the N(H 2 )-χ c 18 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 2 number density used in the models. If we adopt a lower value for thickness, the H 2 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 18 O fractional abundances of ∼3 × 10 −7 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 2 )-χ c 18 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 2 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 2 ) but CO is not abundant enough for detectable emission.

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 −1 , the characteristic accretion timescale to reach the observed line mass of 36-41 M pc −1 is on the order of 0.5-1 Myr based on the timedependent 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 −1 (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.

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(Smith et al. , 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   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. (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 created by 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.  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 by green dashed lines.
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 sheetlike (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 −1 , 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 2 column density of ∼4 × 10 21 cm −2 . 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 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 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.

Summary and conclusion
We have performed 13 CO (1-0), C 18 O (1-0), C 17 O (1-0), 13 CO (2-1), C 18 O (2-1), and C 17 O (2-1) line imaging toward the Serpens filament with the IRAM-30 m and APEX-12 m telescopes A170, page 19 of 25 A&A 646, A170 (2021) 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 −1 . Our observations have led to the following main results: 1. We reveal widespread 13 CO (2-1) self-absorption, resulting in a markedly different morphology relative to the filamentary structure traced by C 18 O and C 17 O. We discover a new molecular outflow driven by the Class I protostar Ser-emb28 in 13 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 18 O and C 17 O spectra in the low-brightness temperature regions, we derive the isotopic ratio [ 18 O/ 17 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 [ 18 O/ 17 O] to be 4.1, we constrain the opacity, rotational temperature, and molecular column density for C 18 O and C 17 O simultaneously with Monte Carlo Markov chain calculations. We find that the opacities of C 18 O (1-0) and C 18 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 18 O column densities range from 9.1 × 10 14 cm −2 to 3.0 × 10 15 cm −2 . The C 18 O fractional abundances relative to H 2 range from 5.8 × 10 −8 to 3.5 × 10 −7 , 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 18 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 18 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 gravitydriven 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.  C 17 O transitions show hyperfine structure (HFS) splitting because of the nonzero nuclear spin in 17 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 17 O (1-0) and C 17 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: where τ 0 and 0 are the opacity and LSR velocity at the line center, and σ 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. σ should be the same for different components. Taking C 17 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 −1 and −2.7348705 km s −1 with respect to F = 7/2-5/2. The observed line widths are determined to be 0.25-0.89 km s −1 with a median of 0.52 km s −1 (see Fig. 5), suggesting that σ should be lower than 0.38 km s −1 in the Serpens filament. We therefore use a σ range of 0.10-0.38 km s −1 for our calculations. Since both C 17 O lines are likely to have low opacities, we use a peak total opacity of τ 0 ≤1 for our test. We are thus able to construct the synthetic C 17 O (1-0) and C 17 O (2-1) spectra with Eq. (1) where 0 is assumed to be 0 and T ex is set to be 7 K as the fiducial case. The synthetic C 17 O (1-0) and C 17 O (2-1) spectra are shown in Fig. D.1. Because C 17 O (1-0) and C 17 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 τ 0 and σ for C 17 O (1-0). The highest ratio becomes 0.78 when we use the highest values (τ 0 = 1, σ = 0.38 km s −1 ). This means that the peak intensity can be overestimated by 28% in this extreme case. For C 17 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 17 O (1-0) and C 17 O (2-1) in studying molecular excitation in Sect. 3.4. A170, page 24 of 25