Probing accretion of ambient cloud material into the Taurus B211/B213 filament

Herschel observations have emphasized the role of molecular filaments in star formation. However, the origin and evolution of these filaments are not yet well understood, partly because of the lack of kinematic information. To examine whether the B211/B213 filament is accreting background gas due to its gravitational potential, we produced a toy accretion model and compared its predictions to the 12CO(1--0) and 13CO(1--0) velocity patterns. We also examined the spatial distributions of Halpha, 857 GHz continuum, and HI emission to search for evidence of large-scale external effects. We estimated the depth of the cloud around the B211/B213 filament to be 0.3--0.7 pc under the assumption that the density of the gas is the same as the 13CO critical density. Compared to a linear extent of>10 pc in the plane of the sky, this suggests that the 3D morphology of the cloud is sheet-like. 12CO and 13CO PV diagrams perpendicular to the filament axis show that the emission from the gas surrounding B211/B213 is redshifted to the northeast of the filament and blueshifted to the southwest, respectively, and that the velocities of both components approach the filament velocity as the line of sight approaches the filament crest. The PV diagrams predicted by our accretion model are in good agreement with the observed 12CO and 13CO PV diagrams, supporting the scenario of mass accretion into the filament proposed by Palmeirim et al. Moreover, inspection of the distribution of the Halpha and 857 GHz emission in the Taurus-California-Perseus region suggests that the B211/B213 filament may have formed as a result of an expanding supershell generated by the Per OB2 association. Based on these results, we propose a scenario in which the B211/B213 filament was initially formed by large-scale compression of HI gas and then is now growing in mass due to the gravitational accretion of ambient cloud molecular gas.


Introduction
The observations of the Herschel Gould Belt survey (HGBS) have revealed an omnipresence of parsec-scale filaments in molecular clouds and emphasized their importance for solar-type star formation (e.g. André et al. 2010;Men'shchikov et al. 2010;Arzoumanian et al. 2011;Palmeirim et al. 2013). In particular, most Herschel prestellar cores are found to lie in dense (thermally supercritical) filaments, suggesting that cores generally form by filament fragmentation (eg. Könyves et al. 2015;Marsh et al. 2016;Benedettini et al. 2018). Molecular line observations of the velocity field around cores and filaments further support this view (Tafalla & Hacar 2015). Based on the HGBS results, André et al. (2014) proposed a filament paradigm for star formation, whereby large-scale compression of interstellar material in supersonic flows generates a quasi-universal web of ∼0.1-pc wide filaments in the cold interstellar medium (ISM) and then the denser filaments fragment into prestellar cores by gravitational instability. Recently, Shimajiri et al. (2017) found that the star formation efficiency in dense molecular gas (A v > 8), where filamentary structures dominate the mass budget, is remarkably uniform over a wide range of scales from 1-10 pc to >10 kpc (see also, Gao & Solomon 2004;Lada et al. 2010Lada et al. , 2012Chen et al. 2015). Furthermore, Shimajiri et al. (2017) proposed that this common star formation efficiency in dense gas results from the microphysics of star formation in filaments (see also André et al. 2014). This result suggests the existence of a universal "star formation law" converting dense molecular gas into stars along filaments. Therefore, unveiling how molecular filaments grow in mass and fragment is crucial to understanding star formation in filaments.
The B211/B213 filament system is located in the Taurus molecular cloud, which is one of the nearest star-forming regions (d∼140 pc, Elias 1978). Wide-field mapping observations in 12 CO, 13 CO, C 18 O, N 2 H + , and SO emission revealed a whole network of filamentary structures in the B211/B213 area (Goldsmith et al. 2008;Hacar et al. 2013;Panopoulou et al. 2014;Tafalla & Hacar 2015). Goldsmith et al. (2008) and Palmeirim et al. (2013) found that many low-density striations are elongated parallel to the magnetic field, and that blueshifted and red-A&A proofs: manuscript no. B211_accretion Fig. 1. (le f t) 12 CO (1-0) and 13 CO (1-0) emission observed toward the B211/B213 filament and (right) schematic picture of the velocity components. The 12 CO and 13 CO data are from Goldsmith et al. (2008). The panel (le f t) is adopted from Palmeirim et al. (2013). In panel (le f t), the red color shows the distribution of the 12 CO (1-0) emission with a velocity range of 6.6-7.4 km s −1 , the green color shows 13 CO (1-0) emission with a velocity range of 5.6-6.4 km s −1 , and the blue color shows 12 CO (1-0) emission with a velocity range of 4.2-5.5 km s −1 . The white box perpendicular to the filament axis shows the cut line for the position velocity diagrams shown in Fig. 6. shifted components in both 12 CO (1-0) and 13 CO (1-0) emission are distributed to the southwest and the northeast of the B211/B213 filament, respectively, as shown in Fig. 1. This morphology was suggestive of mass accretion along magnetic field lines into the B211/B213 filament. To quantify mass accretion, Palmeirim et al. (2013) assumed cylindrical geometry and used the observed mass per unit length M line to estimate the gravitational acceleration φ(R) = 2GM line /R of a piece of gas in freefall toward the filament (where R and G denote radius from filament center and the gravitational constant, respectively). The free-fall velocity of gas initially at rest at a cylindrical radius R init ∼2 pc was estimated to reach 1.1 km s −1 when the material reached the outer radius R out ∼0.4 pc of the B211/B213 filament. This estimation was consistent with the velocity observed in CO, suggesting that the background gas accretes into the B211/B213 filament owing to the gravitational potential of the B211/B213 filament. However, the velocity structure was not investigated in detail. Investigation of the velocity structure is crucial to confirm this suggested scenario from the kinematic viewpoint. This is the topic of the present paper.
The paper is organized as follows: in Sect. 2, we describe the 12 CO (1-0) and 13 CO (1-0) data, as well as complementary Hα, 857 GHz, and HI data. In Sect. 3, we estimate the optical depth of the 12 CO (1-0) line and present the 12 CO (1-0) and 13 CO (1-0) velocity structures observed in the B211/B213 cloud. In Sect. 4, we discuss the cloud structure, whether the surrounding material accretes onto the B211/B213 filament from the kinematic viewpoint, and whether the filament is formed by large-scale compression. In Sect. 5, we summarize our results.

Observational data
In this paper, we used the 12 CO (1-0) and 13 CO (1-0) data obtained by Goldsmith et al. (2008); Narayanan et al. (2008) with the 14 m diameter millimeter-wave telescope of the Five College Radio Astronomy Observatory (FCRAO). The half-power beam width of the telescope was 45 for 12 CO (1-0) and 47 for 13 CO (1-0). We applied Gaussian spatial smoothing to improve the signal to noise ratio, resulting in an effective beam resolution of ∼76 , corresponding to ∼0.05 pc at a distance of 140 pc. The velocity resolution of the data is 0.26 km s −1 for 12 CO (1-0) and 0.27 km s −1 for 13 CO (1-0). The rms noise level is 0.1 K (T * A ) for 12 CO (1-0) and 0.05 K (T * A ) for 13 CO (1-0), respectively. As complementary observations of the Taurus cloud region and its large-scale environment, we also used the Hα data 1 of Finkbeiner (2003), as well as Planck 857 GHz 2 (Planck Collaboration Int. I 2014) and HI data 3 (Kalberla et al. 2017) from the archive.
3.2. 12 CO (1-0) and 13 CO (1-0) velocity channel maps Figure 3 shows the velocity channel maps observed in 12 CO (1-0) and 13 CO (1-0). In the maps for V LSR < 3.7 km s −1 , both 12 CO (1-0) and 13 CO (1-0) emission is seen in the northeastern part of the maps (RA, DEC = 4:24:00, 28:15:00). In the channel maps for 4.0 < V LSR < 7.3 km s −1 , enhanced emission is seen toward the B211/B213 filament in both 12 CO (1-0) and 13 CO (1-0). The emission at these velocities is likely to be directly associated with the B211/B213 filament. Furthermore, while the emission at 4 km s −1 < V LSR < 6 km s −1 is distributed to the southwest of the B211/B213 filament, the emission at 6 km s −1 < V LSR < 7 km s −1 is distributed to the northeast of the filament. In the channel maps for V LSR > 7.3 km s −1 , the distribution of the 12 CO (1-0) and 13 CO (1-0) emission is suggestive of an arc-like structure around L1495. Figure 1 (right) is a sketch showing the location of each velocity component.

Modeling of the data and discussion
4.1. 3D morphology of the B211/B213 ambient cloud Here, we discuss the 3D morphology of the material surrounding the B211/B213 filament by comparing the extent of the gas in the plane of the sky and its depth along the line of sight. Hereafter, we refer to the system consisting of the B211/B213 filament and its surrounding gas as the B211/B213 cloud (i.e. red, green, and dark blue areas in Fig. 1 (right)). The projected extent of the B211/B213 cloud in the plane of the sky is more than ∼10 pc. Taking the viewing angle into account, the real extent of the cloud may be larger. At the same time, we can estimate the depth of the cloud along the line of sight under the assumption that the surrounding material is filled by gas with density exceeding the critical density of the 13 CO (1-0) line, since 13 CO (1-0) emission is observed over the entire mapped area. The critical density of 13 CO (1-0), n 13 CO critical , may be estimated as follows: where A ul , σ cross , ν, and T ex are the Einstein spontaneous emission coefficient, collision cross section, collision velocity, and line excitation temperature. The values of A 10 and σ cross in the LAMDA database 4 are 6.294×10 −8 s −1 and 10 −15 cm −2 . The collision velocity can be calculated as v= √ 3k B T ex /m = 10 4 √ T ex cm s −1 , where k B is the Boltzmann constant and m is hydrogen molecular mass. This leads to a value of 1.7 × 10 3 cm −3 for the critical density of 13 CO (1-0) assuming T ex 14 K. Here, we assumed that the excitation temperature T ex of the 13 CO (1-0) line is the same as the dust temperature T dust ∼ 14K derived by Palmeirim et al. (2013) from Herschel data in the ambient cloud around B211/B213 (red and dark blue area in Fig. 1 (right)). Palmeirim et al. (2013) also estimated the mean Herschel column density in the material surrounding the B211/B213 filament to be N H 2 1.4 × 10 21 cm −2 . Thus, the depth of the cloud (=N H 2 /n 13 CO critical ) is estimated to be 0.3 pc. Recently, Qian et al. (2015) independently estimated the depth of the whole Taurus molecular cloud and found a value of ∼0.7 pc using the core velocity dispersion (CVD) method. With a projected extent of more than 10 pc and a depth of ∼0.3 -0.7 pc, we conclude that the 3D morphology of the cloud resembles a sheet-like structure (see Fig. 4). From HC 3 N (2-1) and (10-9) observations, Li & Goldsmith (2012) found that the depth of the dense (∼10 4 cm −3 ) portion of the B213 region (i.e. the dense filament) was ∼0.12 pc. This is smaller than our estimate for the depth of the ambient molecular gas, but consistent with the view that the dense inner part of the B213 filament is a cylinder-like structure of ∼0.1 pc diameter (Palmeirim et al. 2013), embedded in a lower-density sheet-like cloud.

Accretion of background gas into the B211/B213 filament
Here, we compare the velocity pattern seen in 12 CO (1-0) and 13 CO (1-0) emission with the prediction of an accretion gas model, in order to investigate whether the B211/B213 filament accretes ambient cloud gas from a kinematic viewpoint.

Observed position-velocity diagrams
As mentioned in Sect. 3.2, the highly blueshifted and redshifted components at V LSR < 3.7 km s −1 and V LSR > 7.3 km s −1 do not seem to be directly connected to the B211/B213 cloud/filament. To focus on the velocity field of the gas associated with the B211/B213 filament, we subtracted these two components as follows. We applied Gaussian fitting with N Gaussian components to each pixel, where N=1, 2, 3, 4, or 5. Wherever the signal to noise (S/N) ratio of the residual peak intensity was less than 5, the fit was deemed to be acceptable and the corresponding spectrum was assumed to consist of N Gaussian components. Then, if the peak LSR velocity of a Gaussian component was lower than 4.0 km s −1 or higher than 7.0 km s −1 , the component was not considered to be associated with the B211/B213 filament or cloud and was subtracted from the data cube (see also  (1-0) velocity channel maps after subtracting these components. Hereafter, we used these subtracted data cubes. Figure 6 shows the resulting position-velocity (PV) diagrams in 12 CO (1-0) and 13 CO (1-0) along a direction perpendicular to the B211/B213 filament as indicated in Fig. 1. On these PV diagrams, distinct velocity pattern can be recognized in 12 CO (1-0) and 13 CO (1-0) toward the filament (|offset| < 10 ∼ 0.4 pc). This is probably due to differing optical depths in the two lines. As described in Sect. 3.1, the 12 CO (1-0) optical depth toward the filament is > 50 and much larger than the optical depth toward the outskirts of the filament, suggesting that the 12 CO (1-0) emission only traces the surface of the filament. In the outskirts of the B211/B213 filament (|offset| > 10 ), the blueshifted emission is distributed to the southwest (offset > 0 ) and the redshifted emission is distributed to the northeast (offset < 0 ) of the filament. It can be seen that the velocities of the blueshifted and redshifted components approach the velocity of the B211/B213 filament as the offset approaches 0 (i.e. the crest of the filament). Transverse velocity gradients perpendicular to the major axis of filaments have been also observed toward several dense filaments in the Serpens cloud (Dhabal et al. 2018) as well as the main filament in the northwestern part of the L1495 subregion (Arzoumanian et al. 2018b).

Gas accretion model
The PV diagrams in Fig. 6 show an asymmetric velocity distribution on either side of the 0 position (filament crest), suggesting that the sheet-like ambient cloud surrounding the B211/B213 filament has a different inclination to the plane of the sky to the northeast and the southwest of the filament. To investigate whether the B211/B213 filament accretes gas from the ambient cloud, we thus constructed a 3-component toy model (one filament component and two components for the northeastern and southwestern sheets) under the assumption that the sheet components to the northeast (red-shifted) and the southwest (blues-  shifted) lie on the near and far sides of the B211/B213 filament, respectively, as shown in Fig. 4. Our modeling procedure is summarized in the schematic picture shown in Fig. A.3.
• Model for the central filament component First, we produced a model for the filament. Herschel observations of nearby clouds have shown that the radial column density profiles of molecular filaments in the radial direction R (i.e. perpendicular to the filament crest) can be well described by the following "Plummer-like" function (Arzoumanian et al. 2011;Palmeirim et al. 2013): where ρ c , Σ p , µ, m H , N 0 H 2 , p, and R flat are the central density of the filament, the column density as a function of radius R , the mean molecular mass, the hydrogen atom mass, the central column density, the index of the power-law density profile at large radii (R R flat ), and the radius of the flat inner region, respectively. For the B211/B213 filament, we adopted N 0 H 2 =1.4×10 21 cm −2 , p=2.0, and R flat =0.03 pc from the fitting results of Palmeirim et al. (2013). We assumed that the filament itself lies in the plane of the sky and that the shape of the intensity profile of the B211/B213 filament as traced in 12 CO (1-0) and 13 CO (1-0) emission is the same as that found in the Herschel column density map. Then, we rescaled the peak integrated intensity to be 2 K km s −1 as observed in 13 CO (1-0).
Approximating the Plummer density profile of the filament by a broken power-law, the gravitational potential in the radial direction R can be expressed as follows 5 (cf. Hennebelle & André 2013): where ρ flat and R out are the density of the filament at R ≤ R flat and outer radius of the filament, respectively. We adopted • Models for the northeastern and southwestern sheet components Second, we produced models for the northeastern and southwestern sheet components assuming that the B211/B213 filament accretes the gas of the sheets as a result of its gravitational potential. Taking into account the pressure gradient force, conservation of energy for a parcel of unit mass of the ambient cloud falling onto the central filament may be expressed as follows (cf. Smith 1994Smith , 2012:, 1 2 (V 0 init,N/S ) 2 +φ(R init,N/S )+C 2 s,eff ln(ρ init ) where V(R) is the projected velocity and C s,eff is the effective sound speed. The projected velocity V(R) of the gas flow can thus be expressed as follows: V(R)=V filament ± 2 1 2 (V 0 init,N/S ) 2 +φ(R init,N/S )−φ(R )+C 2 s,eff ln ρ init ρ(R ) × cos(θ N/S ), (8) 5 Here, R denotes the radius corrected for inclination to the line of sight, where the relation between the corrected radius R and radius in the sky plane R is R = R/ sin(θ N/S ) and θ N/S is the inclination angle of the northeastern/southwestern sheet component to the line of sight. where V filament , V 0 init,N/S , R init,N/S , and ρ init are the systemic velocity of the filament, the velocity of the accreting gas at the initial point corrected for inclination, the initial radius of the accreting gas corrected for inclination, and the volume density at the initial point, respectively. Here, we define V 0 init,N/S as (V init,N/S − V filament )/ cos(θ N/S ), where V init,N/S is the projected velocity of the northeastern/southwestern sheet component at R init,N/S . We adopted C s,eff = 0.9 km s −1 from C s,eff δV FWHM ( 12 CO)/ √ 8 ln 2, where δV FWHM ( 12 CO) is the 12 CO (1-0) line width (=2.1 km s −1 ) observed toward the B211/B213 filament. Wherever the value of C 2 s,eff ln ρ init ρ(R ) was larger than 1 2 (V 0 init,N/S ) 2 + φ(R init,N/S ) − φ(R ), we adopted V(R) = V filament . The Herschel observations show that the density profile of the B211/B213 filament is proportional to R −2 at R ≤ R out and has a shallower slope at R ≥ R out (Palmeirim et al. 2013). Furthermore, the slope for the southwestern sheet component is slightly steeper than the slope for the northeastern sheet component. At R > R init , the gas density in the model was assumed to be constant.
To summarize, we assumed the following density distribution as a function of radial direction R (see Fig. 7): For the northeastern sheet component, For the southwestern sheet component, We also assumed that both sheet components have integrated intensities of ∼1 K km s −1 as observed in 13 CO (1-0). To get a good agreement between the models and the observations (see Appendix B), we adopted θ N = 70 • , V init,N = 6.8 km s −1 , and R init,N = 10 pc for the northeastern sheet component and θ S = 20 • , V init,S = 4.4 km s −1 , and R init,S = 10 pc for the southwestern sheet component. The parameters of our model are summarized in Table 1. For simplification, we assumed constant inclinations to the line of sight, of 70 • for the northeastern sheet component and 20 • for the southwestern sheet component, in our model. In reality, however, the inclinations of the two components may vary smoothly with radius from the B211/B213 filament and match on the filament crest.

• Combined 3-component model
We first generated an integrated intensity distribution and a peak velocity field for each of the three components with IDL (Interactive Data Language). Using the MIRIAD task velimage 6 , we then produced individual data cube components for the filament and the two sheet components assuming uniform velocity dispersions of 1.3 km s −1 for the filament and 0.9 km s −1 for the sheet components. The velocity dispersions were obtained from fitting the observed 13 CO (1-0) spectra. One of the reasons for the larger velocity dispersion observed in 13 CO toward the central filament may be that the B211/B213 filament contains several velocity subcomponents (or "fiber-like" structures, Hacar et al. 2013), possibly as a result of accretion-driven turbulence (cf. Hennebelle & André 2013;Heitsch 2013;André et al. 2014) 7 .
Finally, we used IDL to co-add the three individual data-cube components and produce a combined model data cube.
• Large-scale kinematic model We adopted initial velocities (corrected for inclination) of V 0 init,N = 1.8 km s −1 (= [V init,N − V filament ]/ cos(θ N )) for the northeastern sheet component and V 0 init,S = −1.9 km s −1 (= [V init,S − V filament ]/ cos(θ S )) for the southwestern sheet component. This almost symmetric initial velocity pattern after correction for inclination is suggestive of gravitational accretion. If the accreting gas comes from far away positions R far,N/S ( R init ) and is accelerated by the gravitational potential of the B211/B213 filament/cloud, the line of sight velocity at R far,N/S is likely to be similar to V filament . The positions R far,N/S for the northeastern and southwestern sheet components can be estimated from the equation of V init,N/S = 2[φ(R far,N/S ) − φ(R init,N/S )] × cos(θ N/S ) since the pressure density gradient is probably small and can be neglected. We adopted R init,N = 10 pc, and R init,S = 10 pc, respectively. Thus, assuming that the initial velocities are entirely generated by gravitational acceleration, the surrounding gas for the northeastern and southwestern sheet components would have to come from R far,N (R far,N ) = 270 (260) pc and R far,S (R far,S ) = 6 velimage makes a data cube output_cube(x, y, v centroid ) from an input integrated intensity image input_intensity(x, y), input centroid velocity image input_velocity(x, y), and dispersion σ. The v centroid -values are the centroid velocity and are input as an (x, y) image. The output cube image is produced as output_cube(x, y, v centroid ) = input_intensity(x, y) × exp(−(v centroid − input_velocity(x, y)) 2 /(2σ 2 ))). 7 In recent numerical simulations of this process, Seifried & Walch (2015) did find that the accretion flow increases the velocity dispersion of the central filament, and Clarke et al. (2017) suggested that fiber-like structures could be produced as a result of the vorticity generated by an inhomogeneous accretion flow. θ N =70 • and θ S =20 • are assumed. Fig. 6 corresponds to -2.4 pc < offset < 1.6 pc in this figure. The vertical dashed lines mark R init,N = -9.4 pc (R init,N = 10 pc) and R init,S = 3.4 pc (R init,S = 10 pc). 520 (180) pc (see Fig. 8). Here, for simplification, we did not include the mass of the sheets when estimating the gravitational potential. Thus, these R far,N/S values should be considered upper limits. The HI emission observed at V LSR ∼ 6 km s −1 , which corresponds to the systemic velocity of the B211/B213 filament, has an extended emission with an extent of several × 100 pc which is consistent with the above value of R far,N (see also Sect. 4.3 and Fig. A.4). Thus, one of the reasons why the initial velocities at R init,N/S in the northeastern and southwestern sheet components differ from the velocity of the filament may be the large-scale effect of the gravitational potential of the B211/B213 cloud/filament. We will discuss another possible explanation in Sect. 4.3.

Comparing the combined model with the observations
The synthetic PV diagram predicted by the model is shown in Fig. 6 (c) for comparison with the PV diagrams observed in 12 CO (1-0) and 13 CO(1-0) (Figs. 6(a) and 6(b)). A good quantitative agreement especially with the 13 CO (1-0) diagram can be seen. In particular, in both the model and the observed PV diagrams, the velocity of the gas surrounding the B211/B213 filament (red and dark blue area in Fig. 1(right)) approaches the systemic velocity of the B211/B213 filament as the positional offset approaches 0 (i.e. the filament crest). While the gas is accelerated by the gravitational potential of the filament/cloud at large scales (several×10 pc), it is decelerated by the pressure gradient force of the dense filament at small scales (several pc) (see Fig. 8 and Fig. 9). The good agreement between the model and the data indicates that observational kinematic constraints are consistent with the B211/B213 filament accreting background cloud material as a result of its gravitational potential. This provides strong support to the scenario of mass accretion along magnetic field lines into the filament proposed by Palmeirim et al. (2013). The mass accretion rate onto the B211/B213 filament was estimated to be 27-50 M pc −1 Myr −1 by Palmeirim et al. (2013), suggesting that it took ∼1-2 Myr to form the filament. Thus, accretion of gas from the ambient cloud in B211/B213 likely plays a key role in the evolution of the filament.

Formation of the B211/B213 filament by large-scale compression
As described in Sect. 4.2, we adopted different inclinations for the northeastern sheet component (i=70 • ) and for the southwestern sheet component (i=20 • ) in our model to get a good agreement with the observations. This suggests that the B211/B213 cloud is actually shaped like a shell (see Fig. 9). One possibility is that this shell-like structure was produced by large-scale compression.
In this section, we try to investigate whether the cloud surrounding the B211/B213 filament is affected by large-scale flow phenomena using wide-field Hα maps tracing gas ionized by massive stars (Finkbeiner 2003), the Planck 857 GHz dust continuum map tracing cold dust (Planck Collaboration Int. I 2014), and HI map tracing lower density atomic gas (HI4PI Collaboration et al. 2016). Figure 10 (see also Figs A.4 and A.5) compare the spatial distributions of the Hα and 857 GHz emission in the Taurus-California-Perseus region (e.g. Taurus, Auriga, California, and Perseus). The 857 GHz dust emission traces each molecular cloud and exhibits a hole-like structure. This hole-like structure can also be seen in HI emission as shown in Fig. A.4 and Fig. A.5. The Hα emission fills the hole-like structure seen in the 857 GHz dust emission near the center of the field. The Taurus, California, and Perseus molecular complexes traced by the 857 GHz dust emission are distributed at the edge of the hole-like structure. Lim et al. (2013) also found evidence of a shell-like structure using dust extinction and 12 CO (1-0) maps. The holelike structure may result from the expansion of a large-scale supershell produced by a supernova in the Per OB2 association that compresses the Taurus cloud from the far side (Olano & Poeppel 1987;Bally et al. 2008). An Hα absorption feature is detected toward the Taurus cloud (see Fig. 10 and Fig. A.6), suggesting that the Taurus cloud lies at the front surface of the large-scale supershell produced by the Per OB2 association. The distance to the Per OB2 association is estimated to be 340 pc from the Sun (Cernis 1993), while the distance to the Taurus cloud is ∼140 pc (Elias 1978). These distances are consistent with the Taurus cloud lying in front of the Per OB2 association. The B211/B213 filament also appears to be in front of the HI shell (see Fig. 10 in Chapman et al. 2011). This morphology suggests that the B211/B213 filament may have formed as a result of an expanding supershell. This may provide another reason for the different initial gas velocities differed for the northeastern and southwestern sheet components besides large-scale acceleration by the gravitational potential of the B211/B213 cloud (see Sect. 4.2.3). The Local Bubble surrounding the Sun might also compress the Taurus cloud from the opposite direction. The Local Bubble surrounding the Sun was produced by supernovae (Snowden et al. 1998;Sfeir et al. 1999) and the wall of the Local Bubble is located close to the Taurus cloud (Könyves et al. 2007;Lallement et al. 2014).
Interestingly, the Planck 353 GHz data show variations in the polarization fraction (i.e., polarized intensity/total intensity) across the B211/B213 filament, with lower and higher polarization fractions in the southwestern and northeastern parts of the filament, respectively (see Fig. 10 in Planck Collaboration Int. XXXIII 2016). If the gas surrounding the filament is shaped as a shell-like structure with an ordered magnetic field in the plane of each sheet component and if the southwestern sheet component is oriented closer to the line of the sight compared to the northeastern sheet component (cf. Fig. 4(right)), the polarization fraction is expected to be lower in the southwestern area (dark blue in Fig. 1(right)) than in the northeastern area (red in Fig. 1(right)) assuming uniform dust grain properties. Moreover, both the polarization fraction and the polarization angle show smooth variations across the filament, which is consistent with the northeastern and southwestern sheets being curved (i.e., shell-like). The Planck polarization results are therefore support the present model.
Using magnetic magnetohydrodynamic (MHD) numerical simulations, Inutsuka et al. (2015) and Inoue et al. (2018) have argued that multiple compressions associated with expanding bubbles can create star-forming filamentary structures within sheet-like molecular cloud. A similar model of anisotropic filament formation in shock compressed layers has been proposed by Chen & Ostriker (2014), also based on MHD simulations. Such anisotropic filament formation model naturally account for transverse velocity gradients across the B211/B213 filament (see Fig. 1) and other dense molecular filaments (Dhabal et al. 2018), and are good agreement with the observational picture presented.
Based on these considerations, we propose the following scenario for the formation and evolution of the B211/B213 filamentary system: 1. A large-scale flow associated with the Per OB2 supershell compressed and deformed the cloud centered on the B211/B213 filament and created a bent shell-like structure. 2. Owing to its strong gravitational potential, the B211/B213 filament is growing in mass due to accretion of background gas from the surrounding shell-like structure.

Conclusions
To examine whether the B211/B213 filament is accreting gas from the surrounding cloud, we investigated the velocity patterns observed in the 12 CO (1-0) and 13 CO (1-0) lines. Our main findings may be summarized as follows: 1. The optical depth of the 12 CO (1-0) line was estimated to be ∼3-300. The 12 CO optical depth toward the B211/B213 filament is much larger than that toward the outskirts of the filament. The position-velocity diagrams observed in 12 CO (1-0) and 13 CO (1-0) exhibit different velocity patterns close to the filament, which is likely due to different optical depths. 2. The 12 CO (1-0) and 13 CO (1-0) emission from the B211/B213 filament are seen at an LSR velocity of ∼6 km s −1 . In the northeastern and southwestern parts of the B211/B213 filament, the 12 CO (1-0) and 13 CO (1-0) emission are redshifted and blueshifted, respectively. The line of sight velocities are gradually approaching the systematic velocity of the filament as one gets closer to the filament. 3. The linear extent of the cloud around the B211/B213 filament is more than 10 pc in the plane of the sky. In contrast, the depth of the cloud along the line of sight is estimated to be ∼0.3-0.7 pc (=N H 2 /n 13 CO critical ) under the assumption that the density of the surrounding material is the same as the critical density of 13 CO (1-0). These results suggest that the 3D morphology of the gas cloud surrounding the B211/B213 filament is sheet-like. 4. To investigate whether the B211/B213 filament is in the process of accreting the surrounding gas material, we compared the velocity patterns observed in 12 CO (1-0) and 13 CO (1-0) with our 3-component model. The predictions of the model were found to be in good agreement with the distribution of 12 CO (1-0) and 13 CO (1-0) emission in the observed position-velocity diagrams, supporting the scenario of mass accretion along magnetic field lines into the B211/B213 filament proposed by Palmeirim et al. (2013).
5. From an inspection of the wide-field spatial distributions of Hα and 857 GHz dust emission in the Taurus-California-Perseus region, we concluded that the B211/B213 filament was probably formed as a result of the expansion of a largescale supershell originated in the Per OB2 association. This scenario provides a simple explanation for the different inclinations of the northeastern and southwestern sheet components inferred from our modeling analysis. 6. Based on these results, we propose that a) large-scale compression(s) generated by the Per OB2 association initially formed the B211/B213 filament system, and b) accretion of ambient gas material due to the gravitational potential of the filament is now responsible for the growth of the filament.