Unveiling the importance of magnetic fields in the evolution of dense clumps formed at the waist of bipolar H II regions: a case study on Sh2-201 with JCMT SCUBA-2/POL-2

We present the properties of magnetic fields (B-fields) in two clumps (clump 1 and clump 2), located at the waist of the bipolar H II region Sh2-201, based on JCMT SCUBA-2/POL-2 observations of 850 $\mu$m polarized dust emission. We find that B-fields in the direction of the clumps are bent and compressed, showing bow-like morphologies, which we attribute to the feedback effect of the H II region on the surface of the clumps. Using the modified Davis-Chandrasekhar-Fermi method we estimate B-fields strengths of 266 $\mu$G and 65 $\mu$G for clump 1 and clump 2, respectively. From virial analyses and critical mass ratio estimates, we argue that clump 1 is gravitationally bound and could be undergoing collapse, whereas clump 2 is unbound and stable. We hypothesize that the interplay between thermal pressure imparted by the H II region, B-field morphologies, and the various internal pressures of the clumps (such as magnetic, turbulent, and gas thermal pressure), has the following consequences: (a) formation of clumps at the waist of the H II region; (b) progressive compression and enhancement of the B-fields in the clumps; (c) stronger B-fields will shield the clumps from erosion by the H II region and cause pressure equilibrium between the clumps and the H II region, thereby allowing expanding I-fronts to blow away from the filament ridge, forming bipolar H II regions; and (d) stronger B-fields and turbulence will be able to stabilize the clumps. A study of a larger sample of bipolar H II regions would help to determine whether our hypotheses are widely applicable.

Massive stars with mass >8 M influence their surroundings via (a) energetic jets and outflows during their initial stages, (b) stellar winds, radiation pressure, and H ii regions (which drive shocks and ionization fronts (I-fronts)) during their intermediate stages, and (c) supernova explosions at the end of their lives (e.g., Tan et al. 2014;Motte et al. 2018). These factors impact the second generation of stars through the resultant injection of energy and momentum into the ambient medium. Stellar feedback has two potential consequences -first by injecting the turbulence into the cloud it stabilizes the cloud against its own gravity and maintains it in a state of quasi-static equilibrium (Krumholz et al. 2005;Krumholz & Tan 2007;Federrath 2013); and second, by triggering star formation it reduces the lifetime of a cloud to a few free-fall time scales (Elmegreen 2007;Dobbs et al. 2011;Vázquez-Semadeni et al. 2009). These effects, which are known as negative and positive feedback, respectively, result in reduced or enhanced levels of star formation in a cloud. The key agents involved in the processes described above include magnetic fields (hereafter B-fields), turbulence, gravity, and H ii region feedback. However the relative importance of B-fields in comparison to the other parameters, and the complex interplay between them is poorly understood. Deharveng et al. (2015) and Samal et al. (2018) identified several bipolar H ii regions in our Galaxy using Herschel and Spitzer data analyses. They have suggested that such regions form due to the anisotropic expansion of H ii region in flat-or sheet-like filamentary clouds, in accordance with recent 2D numerical simulations (Wareing et al. 2017(Wareing et al. , 2018. In addition, Samal et al. (2018) found that the most massive and compact clumps are always formed at the waist the bipolar H ii regions (see their Figure 3) showing the signatures of high-mass star formation. Since such clumps are possible sites of massive star and cluster formation, understanding the role of B-fields along with stellar feedback, turbulence, and gravity, as well as the interplay among them holds a key to understand the star formation in such environments. While all other parameters can be relatively well constrained, B-fields are difficult to probe, quantify, and constrain.
Dust grains are believed to be aligned with respect to the B-field orientation via the "radiative alignment torques (RAT)" mechanism (Lazarian 2007;Lazarian et al. 2015;Andersson et al. 2015). The RAT model predicts that asymmetric, nonspherical dust grains rotate as a result of radiative torques imparted by their local radiation field, and so align themselves with their long axis perpendicular to the ambient B-fields (Dolginov & Mitrofanov 1976;Draine & Weingartner 1997;Weingartner & Draine 2003;Lazarian & Hoang 2007). The polarized thermal dust emission yields two quantities -the polarization fraction and the polarization position angle, which reveal polarized dust characteristics and the plane-of-sky component of the B-field morphology, respectively.
Several studies have attempted to probe B-fields in regions with stellar feedback using optical, near-infrared, and sub-millimeter polarization observations Wisniewski et al. 2007;Kusune et al. 2015;Chen et al. 2017;Pattle et al. 2017;Liu et al. 2018). These studies have demonstrated that initially weak Bfields become stronger as a consequence of feedbackdriven compression. These stronger B-fields play a crucial role in the formation and evolution of a variety of structures around H ii regions. Eswaraiah et al. (2017) have carried out NIR polarimetry towards RCW57A, a bipolar H ii region hosting a filament and dense clumps at the waist of H ii region, and found that B-fields are not only important to the formation and evolution of the filamentary cloud, but also strong enough to constrain the flows of expanding I-fronts to form the bipolar H ii regions. However, they could not probe the B-fields in the deeply embedded clumps under the influence of early stellar feedback, due to their heavy extinction. In this study, we probe B-fields in the dense clumps located at the waist of a geometrically simple bipolar H ii region Sh2-201 (hereafter S201).

S201,
with central coordinates of RA (J2000)=03 h 03 m 17. s 9, Dec (J2000) = +60 • 27 52 , is located to the east of the W5-E star-forming complex, as shown in Figure 1. This region is located at a distance of 2 kpc in the Perseus arm (Megeath et al. 2008;Hachisuka et al. 2006). It is a part of an elongated (∼15 ) filamentary cloud of mass 3.3×10 4 M , as seen in 13 CO (Niwa et al. 2009, see their Figure 1). The local standard of rest velocities (V LSR ) for the clumps of the W5-E region, as well as of S201, lie between ∼ −38 km s −1 and ∼ −40 km s −1 (Niwa et al. 2009). Similarly, the velocities of radio recombination lines (RRLs) of the   Fich et al. 1990) are also in close agreement with those of the molecular gas of W5-E. Based on the distributions of (a) young stellar objects (Class 0, Class I, and Class II) from Spitzer (Koenig et al. 2008) and Herschel (Deharveng et al. 2012) observations, (b) the physical conditions of the cold dust, (c) H ii regions, and (d) exciting OB type stars, Deharveng et al. (2012) suggested that the entire W5-E complex and S201 are formed along a same parental dense, sheet-like, filamentary molecular cloud (see Figure 1). These results corroborate the hypothesis that S201 is a part of W5-E (see Figure 1). Figure 2 shows a zoomed-in view of S201. NIR observations (Ojha et al. 2004) reveal that S201 hosts a compact embedded star cluster containing more than one hundred stars and the most luminous member of the cluster is an O6-O8 zero age main sequence star (green plus symbol; Figures 2, 3, and 5(a)). As can be seen from Figure 2, S201 is made up of two lobes extending from the center of the H ii region and two dense clumps (namely, clump 1 and clump 2) at its waist. Radio, molecular hydrogen, and Brγ images have revealed arc-like or bow-like structured photo-dissociation regions at the interface between the H ii region and the clumps, highlighting the interaction between them (Ojha et al. 2004). Several candidate Class 0 and Class I sources have been found within the vicinity of the clumps (Koenig et al. 2008;Deharveng et al. 2012). Since the life-time of Class 0/I YSOs is the order of 10 5 yr (e.g., Evans et al. 2009), the clumps age is 10 5 yr. Therefore the clumps are likely to be in the early stages of their evolution, and so ideal candidates for investigating the interplay between B-fields, turbulence, gravity, and thermal pressure, and their implications for the formation and evolution of dense clumps and bipolar H ii regions.
This paper is organized as follows. Section 2 describes JCMT SCUBA-2/POL-2 observations, data reduction, and analyses. This section also presents molecular line ( 13 CO and C 18 O) data from JCMT/HARP. Results based on the detailed analysis of B-field morphology and correlations between B-fields and intensity gradients (based on the VLA 21cm data) are presented in Section 3. In this section, we also derive various parameters such as dust properties (gas column and number densities, and mass), gas kinematics (velocity dispersion and turbulent pressure), angular dispersion in B-fields (using structure function and auto-correlation function analyses), estimated B-field strength, and ionized gas properties (thermal and radiation pressure). Section 4 discusses the interplay between various parameters, stability analyses based on virial and critical mass estimates, and the consequences for the formation and evolution of clumps, and the formation of bipolar H ii regions. Finally, the conclusions of our current study are summarized in Section 5. Dust continuum polarization observations have been conducted using the POL-2 polarimeter installed on the SCUBA-2 camera (hereafter SCUBAPOL2) at the James Clerk Maxwell Telescope (JCMT; Holland et al. 2013), a 15 m single dish submillimeter observatory located on the summit of Mauna Kea in Hawaii, USA. POL-2 observations of the S201 region (project code: M17BP041; PI: Eswaraiah Chakali) were carried out on 2017 November 18 using the POL-2 DAISY mapping mode (Holland et al. 2013;Friberg et al. 2016). Three sets of observations were acquired under JCMT Band 1 weather conditions, during which the atmospheric optical depth at 225 GHz, τ 225 , was 0.03. Each set was observed for 30 minutes, resulting in a total integration time of ∼1.5 hr. . Column density map is overlaid using white contours at levels of [4,6,12,18,24,36,48,60,72, 84]% of the peak column density of 6.94×10 22 cm −2 . The Class 0 sources, based on Herschel 100/160 µm data, are marked with cyan circles (Deharveng et al. 2012). The Class I, II, and III sources, respectively, are marked with square, diamond, and cross symbols (Koenig et al. 2008). Positions of the two clumps, of swept-up matter around the two ionized lobes of the Hii region, and of the ionizing source (green plus) are shown. Both dust temperature and column density maps are provided by Deharveng et al. (2012). A scale length of 2 ∼ 1.2 pc is shown.
The POL-2 DAISY scanning mode observes a fully sampled circular region of 15 arcmin diameter. The rms noise is uniform within the central 3 -diameter region of the DAISY map, increasing towards the outer parts of the map. POL-2 data were taken at 450 and 850 µm simultaneously, with a resolution of 9. 6 arcsec and 14. 1, respectively. Here we present the results of 850 µm data only, due to the low sensitivity of the 450 µm data. A flux calibration factor (FCF) of 725 Jy pW −1 beam −1 was applied to the 850 µm Stokes I, Q, and U parameters. This FCF value was derived by multiplying the typical SCUBA-2 FCF of 537 Jy pW −1 beam −1 (Dempsey et al. 2013) by a transmission correction factor of 1.35 measured in the laboratory and confirmed empirically by the POL-2 commissioning team using observations of Uranus (Friberg et al. 2016).
The POL-2 data were reduced using the Starlink pol2map 1 routine (adapted from the SCUBA-2 data reduction procedure makemap), which was recently added to the Starlink Smurf mapmaking software suite (Berry et al. 2005;Chapin et al. 2013;Currie et al. 2014). To correct for the instrumental polarization at JCMT/850 µm, we employed the 2018 January IP model during the data reduction, which was extensively tested by the POL-2 commissioning team (Friberg et al. 2016(Friberg et al. , 2018. Instrumental polarization is typically ∼1.5% of the measured total intensity (Friberg et al. 2018). More details on the equations and procedures used to derive the polarization measurements: the debiased degree of polarization [P (%)], polarization angles [θ ( • )], Stokes parameters [Q (%) and U (%)], intensity [I (mJy/beam)], and polarized intensity [P I (mJy/beam)] along with their uncertainties, can be found in recentlypublished work (Wang et al. 2019;Coudé et al. 2019;Liu et al. 2019, and references therein).
The final Stokes I, Q, and U maps have a pixel size of 4 ; however, the polarization vector catalog was created by setting the bin size parameter in the third step of pol2map to 12 , in order to achieve better sensitivity. The mean rms noise in the Stokes I measurements, σ I , is ∼5 mJy/beam (note that the mean rms noise in the 4 pixel-size Stokes I map is ∼14 mJy/beam). In order to infer the B-field orientation in the clumps, we have excluded the data corresponding to fainter regions whose polarization measurements are generally noisedominated. Therefore, we adopted the following data selection criteria: ratio of intensity to its uncertainty, I/σ I , > 10; and ratio of polarization fraction to its uncertainty, P/σ P , > 2, yielding a total of 62 polarization measurements, which are listed, along with their coordinates, in Table 1. We also listed I and P I along with their uncertainties. It should be noted that the θ values listed have a correction of 90 • , and hence infer the B-field orientation 2 projected on the plane of sky.
In this work, we investigate the morphology and strength of B-fields. Results based on the polarization characteristics and alignment efficiency of the dust grains will be published elsewhere (Eswaraiah et al. in prep), and will consist of various analyses of the relationship between P and I using the POL2 data of S201.

Molecular lines data from JCMT HARP
The JCMT is also equipped with the Heterodyne Array Receiver Program (HARP)/Auto-Correlation Spectral Imaging System (ACSIS) high-resolution heterodyne spectrometer, capable of observing molecular lines between 325 and 375 GHz (or 0.922 mm and 0.799 mm). HARP is a 4 × 4 detector array that can be used in combination with ACSIS to rapidly produce large-scale velocity maps of astronomical sources (Buckle et al. 2009). In this paper, we use archival 13 CO (3-2) and C 18 O (3-2) molecular line data (∼14 resolution, project ID: M09BU04, PI: Mark Thompson, observed on 2009-08-25) to examine the distribution of gas and to extract the gas velocity dispersion values in clumps 1 and 2.

B-field morphology
The measured P values range from ∼2% to ∼25% with a mean and standard deviation of ∼7±5%, while the B-field orientations (θ) range from ∼4 • to ∼177 • with a mean and standard deviation of ∼99±50 • ; a higher standard deviation implies a widely distributed B-field morphology with multiple components. The mean measured uncertainties in P and θ are ∼2±1% and ∼7±4 • , respectively. The mean uncertainties in Stokes parameters, σ I , σ Q , and σ U , are found to be ∼5, ∼2, and ∼2 mJy beam −1 , respectively. Similarly, the mean uncertainties in polarization (σ P ) and B-field orientation (σ θ ) are found to be ∼2% and ∼7 • , respectively.
Our aim is to derive various parameters for the two clumps at the waist of Sh201. We thus separate the polarization data according to the areas covered by individual clumps, and by the ionized medium. Of the 62 total measurements, we find that 36 and 18 are in the direction of clumps 1 and 2, respectively, while the remaining 8 measurements are located between or away from the two clumps, and are excluded assuming that they are not representatives of either clump.
The B-field geometry, based on our 62 measurements, is superimposed on the color composite of POL-2 Stokes I, Herschel/SPIRE 250 µm, and Herschel/PACS 70 µm images shown in Figure 3. Red and gray contours correspond to the distributions of dust emission (based on POL-2 850 µm I map) and of the H ii region (based on VLA 1.45 GHz/21 cm continuum 3 ), respectively. Evidently, B-fields in clump 1 follow a bow-like morphology and are conspicuously compressed at the interface region between the dust emission and the ionized medium. This interaction can be witnessed from the closely spaced 21-cm contours. B-fields also seem to be compressed in clump 2 but with a lower degree of curvature.
Histograms of the B-field orientation, shown in Figure 4, reveal the existence of two major components in clump 1. One component, located near the interaction region of the H ii region and the dust emission (POL2 Stokes I), peaks at ∼150 • and is oriented northwestsoutheast. The other component, located on the eastern side of clump 1, is oriented at ∼115 • , approximately along the east-west direction, nearly parallel to the major axis (position angle ∼83 • ) of clump 1. Conversely, the B-fields in clump 2 exhibit a single component with a prominent peak at ∼50 • and is oriented northeast-  [3,6,12,24,48,96,192] × the rms noise of 14 mJy/beam. Gray contours, corresponding to the VLA/21 cm continuum emission representing the distribution of the ionized medium of H ii region, are drawn at [1,3,6,12,24,48,96,206] × the rms noise of 2.3×10 −4 mJy/beam (where beam size ∼ 17 × 13 ). In both panels, reference vectors with a B-field orientation of 90 • , along with mean uncertainty of 7 • are shown. southwest. This B-field component is neither parallel nor perpendicular to the major axis (position angle of ∼15 • ) of clump 2. Figures 3 and 4 imply the presence of multiple B-field components in clumps 1 and 2.

Intensity (ionized gas) gradients versus B-fields
In order to examine whether the multiple B-field components of S201 are shaped by the expanding I-front, we construct intensity gradients using the VLA 21 cm continuum intensity map. More details on making the intensity gradient map are given in Appendix A. To compare the orientation of the intensity gradients (θ IG ) with those of the B-fields (θ B ), we estimate mean θ IG over a ∼14 diameter region (corresponding to the beam size of JCMT) around each θ B vector. Figure 5(a) shows the pairs of θ IG (cyan vectors) and θ B (yellow vectors) overlaid on the VLA 21 cm continuum intensity map. Figure 5(b) shows the offset angle (∆θ) between θ B and θ IG as a function of radial distance (from clump 1 to clump 2) along the magenta line shown in Figure 5(a). The radial profiles of dust and H ii emission, extracted along the magenta line, are also shown to examine their correlation with the ∆θ values. At radial distances < 70 and > 160 , where the dust and H ii emission peaks of clump 1 and 2, respectively, are found, the majority of the ∆θ's lie between ∼50 • and ∼90 • . Notably, perpendicular components prevail around clump 1, and also at clump 2, albeit with less prominence due to the lack of ∆θ values in the range 70 • to 90 • . Here, ∆θ = 90 • refers to perpendicular alignment of B-fields with intensity gradients or parallel alignment of B-fields with intensity contours. Between clumps 1 and 2, i.e., from ∼ 80 to ∼ 150 , the ∆θ values are neither parallel nor perpendicular, lying between ∼ 30 • and ∼ 40 • . These values contribute towards the random component as evident from Figures 5(b) and (c). In addition, there also exist a few parallel components (∼ 0 • -∼ 30 • ) near clump 1 and clump 2.
Therefore, at clump 1, because of the prominent interaction between the dust and H ii emissions, B-fields tend to be perpendicular to the intensity gradients (θ IG ) or parallel to the intensity contours. Whereas at clump 2, as the interaction between the dust and the H ii region is less prominent (based on their peak brightness in comparison to those of clump 1), a lower degree of alignment between B-fields and intensity contours is evident. The cumulative distribution of ∆θ values also confirms the possibility of both perpendicular and random components, as shown in Figure 5(c). In order to study the projection effect from alignments between B-fields and intensity gradients in three-dimensional space to those in the plane of the sky, we also show cumulative distributions based on Monte Carlo simulations (Hull et al. 2014). These simulations randomly select pairs of orientations in three dimensions that have alignments in the ranges 0 • -20 • , 0 • -45 • , 70 • -90 • , or random alignment; then their ∆θ are measured on the plane of the sky. The resulting cumulative distribution functions of the simulations are shown in Figure 5(c). Kolmogorov-Smirnov statistics suggest that, with a probability of 82.5%, our data is consistent with the model distribution corresponding to a ∆θ of 70 -90 • , while there is a probability of only 26.5% that our data have a random component. Therefore, we conclude that the B-fields in the clumps are shaped by the H ii region.

Dust properties of the clumps: column density, number density, and mass
For an idealized cloud, the dust emission at frequencies where optical depth is small can be described by (Hildebrand 1983; see also Li et al. 1999) where S(ν) is the flux density (erg s −1 cm −2 Hz −1 or Jy) from a cloud at distance D, N is the number of spherical grains included in the cloud volume subtended by the beam, Q(ν) is the dimensionless absorption coefficient, σ is the geometric cross section of a dust grain, and B(ν, T d ) is the Planck function for a blackbody at temperature T d . The above equation can be written as follows to construct the column density map from the POL2 850 µm  Figure 3. Integrated fluxes within the green contours correspond to 26σ and 70σ flux levels of 0.006 Jy/beam and 0.016 Jy/beam around clumps 1 and 2, respectively, are used for estimating the thermal pressures exerted on the respective clump surfaces (c.f. Section 3.6 for more details). Locations of the clumps and reference scale length are also shown. Radial profiles of dust and H ii emission are extracted along the magenta line and are drawn in panel (b). (b) The offset between the position angles of B-fields and intensity gradients, i.e., ∆θ = |(θB − θIG)| as a function of radial distance (filled circles). The zero radial distance points to the left edge of the magenta line close to the clump 1. Also, in the right-hand of y-axis, we plotted the radial variation of radio emission (VLA 21 cm; blue dashed lines) as well as of the dust emission (SCUBAPOL2 Stokes I; red dashed line) along the same magenta line shown in panel (a). (c) Cumulative distribution of ∆θ = |(θB − θIG)| is shown with filled circles. Model cumulative distributions (lines with different colors) for ∆θ values of 20 • , 45 • , 70-90 • , and random angles, adopted from Monte Carlo simulations (Hull et al. 2014), are also shown.
We have performed CASA 2D Gaussian fits on the POL2 850 µm Stokes I map, specifically on the pixels around each clump having I > 140 mJy/beam (i.e., >10σ, where σ = 14 mJy beam −1 is the rms noise in the I map with pixel size = 4 ), and extracted the spatial extents of the clumps. The resultant dimensions of the clumps, along with their central coordinates, are given in Table 2. The effective radius of each clump where σ a and σ b are the extents of the semimajor and semiminor axes, respectively, and are found to be 13.3±0.3 (or 0.13±0.01 pc) for clump 1 and 15.2±0.4 (or 0.15±0.01 pc) for clump 2. The mean dust temperatures (T d ), within the dimensions of clump 1 and clump 2, are estimated to be 27 K and 29 K (Deharveng et al. 2012), and are used in the above Equation to estimate the respective column density maps.
We note that the masses determined above are likely to be lower limits on the true masses because they have been estimated using the average dust temperature over the clump areas. We thus measured masses within the areas corresponding to 10σ Stokes I but from the column density map constructed from the Herschel temperature map (for details see Deharveng et al. 2012). The resulting masses are found to be 191±13 M and 30±3 M for clump 1 and clump 2, respectively. Although these values agree within a factor of two with the masses derived from 850 µm Stokes I map, they are likely to better represent the true masses. We thus used these values for further analyses. The number densities and masses estimated above are listed in Table 2. 3.4. Gas properties: velocity dispersion Figure 6 shows the JCMT HARP/ACSIS 13 CO(3-2) and C 18 O(3-2) spectra of the two clumps, averaged over the clump dimensions, in units of brightness temperature (T b , K) versus velocity (V LSR , km s −1 ). We performed Gaussian fitting on each spectrum, resultant peak brightness temperature (T b,p ), central V LSR , and velocity dispersion (σ VLSR ) values are given in Table 3. The σ VLSR values for clumps 1 and 2 are 1.05±0.01 km s −1 and 1.06±0.01 km s −1 in 13 CO, and 0.69±0.02 km s −1 and 0.60±0.04 km s −1 in C 18 O. For clump 1, using C 18 O (J=1-0) measurements, Niwa et al. (2009, clump 9 in their work) derived the mean velocity dispersion over an area somewhat larger than that we consider as 0.71 km s −1 , which is in close agreement with our estimate derived from C 18 O(3-2).
The 13 CO gas traces the extended low-density gas around the clumps, whereas C 18 O traces the highly compact central dense regions of the clumps (see Figure 7). To demonstrate this, we estimate the optical depths, column densities of 13 CO and C 18 O , and the resultant H 2 column densities in Appendix B; these values are listed in Table 4. The optical depths suggest that C 18 O is optically thin and hence traces the densest parts of the clumps, thereby revealing the level of turbulence at the densities traced by 850 µm dust emission.
We further derive the one-dimensional thermal velocity dispersion (σ T ) implied by the kinetic temperature of the C 18 O gas using the relation where k is the Boltzmann constant, T kin is gas kinetic temperature, which is equivalent to the mean dust temperatures (T d = 27 K and 29 K for clump 1 and 2; Table 2) under the assumption that the gas and dust are in local thermal equilibrium. M C 18 O is the mass of the C 18 O molecule and is considered to be 30 amu. The estimated σ T values are found to be 0.087±0.024 km s −1 Figure 6. 13 CO (3-2) and C 18 O (3-2) brightness temperature versus VLSR spectra, averaged over the extents of clumps 1 (top) and 2 (bottom). To view clearly, the spectra of 13 CO and C 18 O are shifted by adding and multiplying by arbitrary numbers as shown in the figure. Best-fit Gaussians are shown with red dashed lines. and 0.090±0.017 km s −1 for clumps 1 and 2. Finally, the non-thermal velocity dispersion (σ NT ), which is due to the turbulence, is estimated by correcting for thermal velocity dispersion using the relation The derived σ NT values are 0.68±0.02 km s −1 and 0.59±0.04 km s −1 for clumps 1 and 2, respectively. Therefore, in the following analysis we used our derived C 18 O non-thermal velocity dispersions, σ NT , to estimate the B-field strength and pressure, turbulent pressure, Alfvénic velocity, Alfvén Mach number, virial mass, etc.

Magnetic field strength
In the subsections above, number densities and gas velocity dispersion values were extracted. Here, we derive the angular dispersion of the B-fields in order to estimate the B-field strength and other crucial parameters.
The plane-of-sky component of B-field strength (B pos ) can be estimated using the Davis-Chandrasekhar-Fermi method (Davis 1951;Chandrasekhar & Fermi 1953, hereafter DCF method), which is based on the assumption that turbulence-induced Alfvén waves will distort B-field orientations. This method assumes that the following two conditions hold: (a) the ratio of turbulent (δB) to large-scale ordered (B o ) B-field components is proportional to the ratio of one-dimensional non-thermal velocity dispersion (σ v ) to Alfvénic veloc- where σ θ is the dispersion in the measured B-field orientation. The DCF method can be applied to sets of Gaussian-distributed polarization angles with angular dispersions less than 25 • (Ostriker et al. 2001). However, B-fields in regions altered by H ii regions or dragged by gravity will generally exhibit multiple components, and so have a significantly higher angular dispersions (e.g., Arthur et al. 2011;Chapman et al. 2011;Santos et al. 2014;Eswaraiah et al. 2017;Chen et al. 2017). As shown in Figure 4, the histograms of B-fields being impacted by the H ii region feedback exhibit either a widely spread distribution of angles or multiple distributions with conspicuously separate peaks. In such regions alternative methods must be employed to extract the underlying dispersion in polarization angles caused by magnetized turbulence. Progress has recently been made towards the accurate estimation of δB/B o through statistical analyses of polarization angles. These analyses include the "structure function" ) and the "auto-correlation function" ) of polarization angles. These are referred to as modified DCF methods, and used to estimate the Bfields in such regions.

Structure function (SF) analysis
In the structure function (SF) analysis ), the B-field is assumed to consist of a largescale structured field, B o , and a turbulent component, δB. The SF analysis demonstrates the variation of dispersion in position angles as a function of vector separation l. At some scale larger than the turbulent scale δ, δB should reach its maximum value. At scales smaller than a scale d, the higher-order terms of the Taylor expansion of B 0 can be canceled out. When δ < l d, the SF follows the form: In this equation, ∆Φ 2 (l) tot , the square of the total measured dispersion function, consists of b 2 , a constant turbulent contribution, m 2 l 2 , the contribution from the large-scale structured field, and σ 2 M (l), the contribution of the measured uncertainty. The ratio of the turbulent to the large-scale component of the magnetic field is given by:  And B o is estimated using the modified DCF relation: Then the estimated plane-of-sky magnetic field strength is corrected by a factor Q B pos = Q B 0 where Q is considered as 0.5 based on studies using synthetic polarization maps generated from numerically simulated clouds (Ostriker et al. 2001), which suggest that B-field strength is uncertain by a factor of two when the B-field dispersion is ≤ 25 • .
The blue filled circles plotted in Figure 8 represent the angular dispersions corrected for uncertainties ( ∆Φ 2 (l) tot − σ 2 M (l)) as a function of length scale, measured from the polarization data. The bin size of 12 , used in Figure 8, corresponds to the grid size of the polarization catalog yielded by pol2map. The maximum value in the current SF is lower than the value expected for a random field (52 • , Poidevin et al. 2010). Equation 6 is fitted to the data using the IDL MPFIT nonlinear least-squares fitting algorithm (Markwardt 2009). The resulting δB 2 1/2 Bo values are 0.40±0.02 and 0.40±0.07 for clumps 1 and 2, and the corresponding B pos strengths are derived using equations 7, 8, and 9 to be 147±15 µG and 65±6 µG.

Auto-correlation function (ACF) analysis
Auto-correlation function (ACF) analysis ) is an extension of SF analysis, which includes the effect of signal integration along the line of sight as well as within the beam. According to Houde et al. (2009), the ACF can be written as: where ∆Φ(l) is the difference in position angle of two vectors separated by a distance l, W is the beam radius (6.0 for the JCMT, i.e., the FWHM beam size of 14 divided by √ 8 ln 2), a 2 is the slope of the second-order term of the Taylor expansion, and δ is the turbulent correlation length. N is the number of turbulent cells probed by the telescope beam and is given by: where ∆ is the effective thickness of the cloud. The ordered magnetic field strength is given by The top panels of Figure 9 show the angular dispersion functions of the polarization vectors in clumps 1 and 2, while the bottom panels show the respective correlated components of the dispersion functions. In our fitting, ∆ is set to 31 ± 1 for clump 1 and 36 ± 1 for clump 2 (these correspond to the effective thickness of the clumps, derived using the relation √ δ a × δ b ; where δ a and δ b are the FWHMs of the major and minor axes, respectively (see Table 2).
Equation 10 is fitted to the ACF data shown in Figure  9 for clumps 1 (top panel) and 2 (bottom panel). The bin width for constructing the ACF (1 − cos[∆Φ(l)] ) function was chosen to be 9 . Various bin widths were investigated; we found that best fit was achieved with a bin width of 9 . The IDL MPFIT non-linear leastsquare fitting algorithm (Markwardt 2009) was used and simultaneously constrained the three fitting parameters (i) δ, (ii) δB 2 B 2 o , and (iii) a 2 ; the best-fit values are listed in Table 2. Using the fitted parameter δB 2 B 2 o , along with derived parameters such as number densities and velocity dispersions, we have estimated the B-field strength using the modified DCF relation (equation 12). The estimated B-field strengths are 266±32 µG and 61±31 µG for clumps 1 and 2 respectively. The turbulent correlation length δ is 13 ± 3 (126±29 mpc) and 7 ± 4 (68±36 mpc) for clumps 1 and 2. The number of turbulent cells (N) is derived to be 1.4±0.4 and 5.0±4.8 for clumps 1 and 2, respectively. The derived parameters are listed in Table 2.
In summary, the SF and ACF analyses yielded two B-field strengths: for clump 1 these differ by a factor of ∼2 (147±15 µG and 266±32 µG), while they are nearly identical for clump 2 (65±6 µG and 61±31 µG), although the ACF-derived B-field strength has ∼50% uncertainty. It should be noted here that the measured uncertainties in the B-fields strengths are 15 µG and 32 µG for clump 1 and 6 µG and 31 µG for clump 2. These are random errors, resulting from propagation of the uncertainties of various parameters (such as gas number density, gas velocity dispersion, and B-field angular dispersion) through the modified DCF formulae. We caution here that, in practice these measurement errors can be dominated by systematic errors associated with various parameters such as dust opacity, dust temperature, flux calibration, clump geometry (and hence depth), and inclination angle, etc. (e.g., Pattle et al. 2017;Wang et al. 2019;Liu et al. 2019). These factors could be more severe in clump 1 because of the existence of more extended emission around the peak of dust emission. Therefore the real uncertainties could be much larger than the quoted errors in B-field strengths.
It is clear from Sections 3.5.2 and 3.5.1 above (see also Table 2) that for clump 2 the ACF-yielded parameters have higher uncertainties in comparison to those yielded by SF. This could be attributed to the relatively small number of vectors in clump 2 causing ACF to fail constraining all the fitting parameters simultaneously (see column 3 of Table 2 for ACF). In addition, since the column density is relatively lower, the beam dilution and signal integration effects may not be important in clump 2. Conversely, for clump 1 these factors seem to be crucial and taken care of ACF. Therefore, we have used B-field strengths derived from ACF (266±32 µG) for clump 1 and from SF (65±6 µG) for clump 2 in the further analyses.
3.6. Ionized gas properties: thermal and radiation pressure in S201 Figure 10 shows the radio continuum view of the S201 H ii region at VLA 21 cm/1.4 GHz. The flux density (S ν ) of S201, estimated by integrating over 3σ contours, is ∼1.0±0.1 Jy, where σ is the rms noise of the 21 cm map. Our 21 cm flux density is, within uncertainty, consistent with flux density at 6 cm (1.2±0.2 Jy; Felli et al. 1987, and references therein). The flux densities at 21 cm and 6 cm reflect a flat spectrum, indicating that the nebula is optically thin at 21 cm. Considering 8302 where S ν is the radio continuum integrated flux at frequency ν, θ D is the angular diameter of the source, D (2 kpc) is the distance from the Sun, and T e is the electron temperature in the ionized plasma. Using this formalism, we estimated n e as 226±11 cm −3 . We then estimated the corresponding thermal pressure, due to ionized gas distributed over the area of diameter ∼120 , to be (5.2±0.3) × 10 −10 dyn cm −2 using the relation P te = 2n e k b T e . The mean thermal pressure given above my be valid for the entire region surrounded by the H ii emission. The radio emission is observed to be uneven in the region of S201 (cf. Figure 10 and also see Figure 5) in the sense that more H ii emission is concentrated close to clump 1 than is concentrated near clump 2. Because of this, the relative impact of H ii emissions, in terms of the thermal pressures acting on the clump surfaces, will be different. Therefore, we estimate average thermal pressures close to the clumps. The integrated fluxes, S ν , within the green contours (see Figure 5) extended over circular diameters of ∼66 and ∼48 , are found to be 0.4±0.2 Jy and 0.05±0.01 Jy for clumps 1 and 2, respectively. Using these parameters, along with the value of T e quoted above, we estimated the electron densities, n e , as 360±80 cm −3 and 207±14 cm −3 , and the resultant thermal pressures, P te , to be (8±2) × 10 −10 dyn cm −2 and (4.7±0.3) × 10 −10 dyn cm −2 for clumps 1 and 2, respectively.
We estimated the mean radiation pressure (P rad ) driven by an ionizing star of spectral type O6V (Ojha et al. 2004;Deharveng et al. 2012). The ionizing flux emitted by one O6V star is q 0 = 4.15×10 10 photons cm −2 s −1 (Sternberg et al. 2003), and each UV photon carries an energy hν = 20 eV, so the estimated P rad , using the relation P rad = hν q 0 /c, is 0.44 × 10 −10 dyn cm −2 . We note here that the same P rad value is used for both clumps and in both cases is negligible in comparison to the thermal pressures. The derived thermal and radiation pressure values are listed in Table 2.

DISCUSSION
Here we discuss the interplay amongst various key parameters, clump stability based on virial and critical mass estimations, and their relevance to the formation and evolution of clumps as well as to the feedback process.
The turbulent Alfvén Mach number (M A ) describes the relative importance of B-fields compared to that of turbulence, and hence it is a key parameter in models of cloud formation and evolution (e.g., Ostriker et al. 2001;Padoan et al. 2001;Nakamura & Li 2008). In the sub-Alfvénic case (M A 1), B-fields are strong enough to regulate turbulence, causing an organized B-field orientation. In the super-Alfvénic case (M A > 1), the turbulence is capable of perturbing the morphology of Bfields.
The Alfvénic velocity, , is estimated to be 1.6±0.2 km s −1 for clump 1 and 0.7±0.1 km s −1 for clump 2. The Alfvén Mach number, M A = √ 3( σNT VA ), is estimated to be 0.8±0.1 for clump 1 and 1.4±0.2 for clump 2. These estimates agree with the findings above based on the ratios of P B /P turb , and imply that turbulent motions are sub-Alfvénic and super-Alfvénic in clumps 1 and 2, respectively.

B-field versus Thermal pressure by H ii region
The ratio of magnetic (c.f. Section 4.1) to thermal pressures (c.f. Section 3.6), P B /P te , is found to be 3±1 and 0.4±0.1 in clumps 1 and 2, respectively. These results imply two contrasting scenarios, in that B-fields dominate over thermal pressure acting on clump 1, whereas thermal pressure dominates over B-fields in clump 2. Evidently, B-fields are strong enough to control the expanding I-front in clump 1 and, conversely, the I-front is strong enough to dictate the B-fields in clump 2. The consequences of this interplay between Bfields and thermal pressures will be discussed in Section 4.6.

Are the clumps gravitationally bound?
Considering the case in which clumps are not supported by B-fields and also not confined by external pressure (in the form of envelope material around the clumps), we estimate the viral masses and virial mass ratios.
In order for self-gravitating clumps to be in virial equilibrium, the following relation between gravitational potential energy (|G|) and internal kinetic energy (E) should hold (McKee & Zweibel 1992) where E = 3/2 Mσ 2 . The gravitational potential energy can be written as where r = R eff and G is the gravitational constant. α corresponds to a geometric factor, a function of eccentricity, and β is a function of the power-law index of the density profile (ρ ∝ r −a , where a = 1.6 for an isothermal cloud in equilibrium; Bonnor 1956). More details on deriving these factors, assuming that the clumps are prolate ellipsoids, can be found in Li et al. (2013a, and references therein). For our given values of σ and r, using the above equations, the virial mass (M vir ) can be estimated using the relation Taking σ = σ VLSR of C 18 O and r = R eff , M vir is estimated to be 50 M and 40 M for clump 1 and clump 2, respectively. For our estimated masses (M ; 191 M for clump 1 and 30 M for clump 2; c.f. Section 3.3), the derived virial mass ratios, R vir = M /M vir are 3.9 and 0.7 for clumps 1 and 2. Therefore, clump 1 is bound by gravity and may collapse once it becomes unstable, whereas clump 2 is gravitationally unbound.

Temperature, and B-fields
Critical mass M C is the maximum mass that can be supported by the combined contributions of internal velocity dispersion (contribution from turbulence and neutral gas temperature, i.e., non-thermal and thermal contributions) and B-field in the clump. The two effects can be represented as which is accurate within 5% to results from more rigorous calculations (McKee 1989). The Jean mass for a non-magnetic isothermal cloud (Bonnor 1956;McKee & Zweibel 1992) is where C eff = C s 2 + σ NT 2 . The thermal sound speed C s = k T kin µH mH was estimated to be 0.28±0.01 km s −1 for clump 1 and 0.29±0.01 km s −1 for clump 2. In this equation, T kin = T dust and is 27 K for clump 1 and 29 K for clump 2. The C eff values are estimated as 0.74±0.02 km s −1 and 0.66±0.04 km s −1 for clumps 1 and 2.
In Equation 19, the envelope pressure (or external pressure) caused by the low-density 13 CO gas can be estimated as where σ env is the velocity dispersion in the low-density envelope, which we treated as σ env = σ VLSR of the 13 CO gas (c.f. Table 3), and the corresponding mean gas number densities n env were estimated over larger extents to be ∼0.9×10 4 cm −3 and ∼0.5×10 4 cm −3 for clumps 1 and 2. The maximum mass that can be supported by B-fields will be where c φ ∼ 0.12 according to numerical simulations (Tomisaka et al. 1988). The estimated critical masses M C are 78 M and 48 M for clumps 1 and 2. The critical mass ratios R C = M/M C are found to be 2.5 and 0.6 for clumps 1 and 2. These results suggest that for clump 1 the support rendered by the combined contributions from thermal gas energy, turbulence, and B-fields is not sufficient to counteract gravity, whereas the opposite situation prevails in clump 2, such that its stability is determined by thermal energy, turbulence, and B-fields rather than by self-gravity. This picture is further corroborated by the presence of a larger number of Class 0 and I source in and around clump 1. In contrast, clump 2 is inactive as there exist no YSOs at its center, except for a few Class I sources formed at its boundary (see Figure 2).

Compressed B-fields and enhanced B-field strength in the clumps
The B-field strength in clump 1 (266±32 µG) is larger by a factor of ∼4 than that in clump 2 (65±6 µG). Additionally, the spread in ADF values (∼30 -∼ 50 • ; Figure  8) as well as ACF values (∼0.05 -∼0.37; top panels of Figure 9) for clump 1 are relatively higher than those of clump 2 (ADF range ∼30 -∼40 • ; while ACF range ∼0.03 -∼0.16). Furthermore, the spread in the offset angles (∆(θ)) between B-fields (θ B ) and intensity gradients (θ IG ) in clump 1 is relatively larger than that of clump 2. These signatures imply that B-fields are more curved and draped around clump 1, thereby following a bow-like structure (see Figure 3). Similar features have been seen in clump 2, but with a lesser degree of curvature. Alūzas et al. (2014), based on 2D MHD simulations, show that when an oblique shock interacts with an isolated cylindrical cloud, the B-fields wrap around the cloud, attaining a roughly circular shape (see their Fig. 1b) similar to the B-field morphologies observed in clump 1 (Figure 3). Based on 3D radiationmagnetohydrodynamic (RMHD) simulations of pillars and globules, Mackey & Lim (2011); Mackey (2012); Mackey & Lim (2013) have shown that in the case of initially strong B-fields (160 µG) oriented perpendicular to an I-front, field lines at the head of a cometary globule get compressed into a curved morphology, by closely following the bright rim in a similar manner to our present observations towards clumps 1 and 2 (see Figure 3). These findings are consistent with other observations and MHD simulations (e.g., Lyutikov 2006;Dursi & Pfrommer 2008;Pfrommer & Jonathan Dursi 2010;Arthur et al. 2011;Santos et al. 2014;Kusune et al. 2015;Klassen et al. 2017).
As can be seen from Figure 5(a), contours of dust and ionized emission are closely spaced in their zone of interaction, suggesting possible compression in the interface between the hot and cold media. This interaction has also compressed B-fields and enhanced their strength. Consequently, the stronger B-fields may shield the clumps from the H ii region. In the perpendicular field case, the B-fields get amplified directly upstream of the cloud where the flow stagnates against it, where Bfield pressure and field tension continue to build (Alūzas et al. 2014, see their Fig. 1c). In MHD simulations, when B-fields get compressed their strengths become enhanced in the shells or clumps by a factor of about 5 to 6 compared to those inside the expanding H ii regions (e.g., Mac Low et al. 1994;Gregori et al. 1999;Klassen et al. 2017, see also Wareing et al. 2017 for a similar enhancement of B-field strength in environments with mechanical stellar feedback). Therefore, our results provide evidence for enhanced B-field strengths in clumps, which we attribute to the effect of thermal pressure. This is because the more the Hii region interacts with the cloud, the more the field lines are compressed, resulting in a considerable enhancement in B-field strengths.
Higher values of B-field strength around ∼50 -∼400 µG, have been measured at the edges of H ii regions using HI/OH Zeeman measurements (Troland et al. 1986(Troland et al. , 2016Mayo & Troland 2012) as well as dust extinction (e.g., Kusune et al. 2015;Eswaraiah et al. 2017;Chen et al. 2017), and emission polarization (e.g., Vallée & Fiege 2005;Pattle et al. 2018). For example, Mayo & Troland (2012) have measured B-field strength of ∼80 µG using the HI Zeeman effect in the photodissociation region (PDR 4 ) DR 22, which they interpret as resulting from B-fields amplified by compression of the PDR due to absorption of the momentum of stellar radiation (also see Pellegrini et al. 2007, for a similar explanation). Based on NIR polarimetry towards the bright-rimmed cloud SFO 74, Kusune et al. (2015) have measured an enhanced B-field strength of ∼90 µG inside the tip rim due to UV-radiation-induced shock. Similarly, enhanced B-field strengths of 100 -300 µG have been inferred in PDRs around ionized regions using radio recombination lines (RRLs; Balser et al. 2016). SCUBAPOL2 observations towards M16 show that the derived B-field strength lies between 170 µG and 320 µG (Pattle et al. 2018). Based on SCUBAPOL observations towards S106, Vallée & Fiege (2005) estimated that its B-field strengths lie between 240 µG and 1040 µG. Similarly, Roberts et al. (1995) conducted OH Zeeman measurements towards S106; their derived B-field strengths range from 100 µG to 400 µG. Evidently, our derived B-field strengths (∼50 -∼200 µG) towards the clumps in S201 are in close agreement with the values quoted in the literature.

The role of H ii region feedback and its relation to the observed turbulence in the clumps
We compare the magnetic (P B ) and turbulent (P turb ) pressures of the clumps with the thermal pressure (P te ) exerted on them by the H ii region, and examine whether turbulence is being injected into the clumps by the H ii region in the presence of B-fields of varying strengths. The comparison among these pressures suggests two relations: (i) P B > P turb > P te in clump 1 and (ii) P te > P turb > P B in clump 2. This implies the dominance of B-fields in clump 1 and of thermal pressure on clump 2.
The stronger B-fields in clump 1 could be able to guide the I-front to escape from the filament-ridge and also shield the I-front from entering clump 1; as a consequence the shock strength will be reduced (e.g., ). Two-dimensional numerical simulations of the interactions between magnetized shocks and radiative clouds show that B-fields external to, but concentrated near, the surface of the cloud suppress the growth of destructive hydrodynamic instabilities (Chandrasekhar 1961;Mac Low et al. 1994;Jones et al. 1996;Fragile et al. 2005), thereby shielding the cloud from erosion or destruction. Conversely, non-magnetized, nonradiative clouds are destroyed on a few dynamical timescales through hydrodynamic Kelvin-Helmholtz, Richtmyer-Meshkov, and Rayleigh-Taylor instabilities (e.g., Klein et al. 1994;Nakamura et al. 2006). Eventually, the B-field dominates over turbulence as is seen in clump 1 (c.f. Section 4.1). In contrast, due to the limited impact of the H ii region on clump 2, the B-field strength has not been enhanced to higher values in this region. Therefore, we hypothesize that due to these relatively weak B-fields the expanding I-front might have driven a shock front into clump 2, resulting in a higher turbulence pressure compared to magnetic pressure, as is seen (c.f. Section 4.1). 4.6. Pressure balance between clumps and stellar feedback, and the consequences Assuming that the primordial filament within which S201 and W5E complexes formed (Deharveng et al. 2012), followed a Plummer-like column density profile (Arzoumanian et al. 2011;Juvela et al. 2012;Palmeirim et al. 2013;André et al. 2019) and also that the primordial B-fields were threaded perpendicular to the filament's long axis (e.g., Chapman et al. 2011;Sugitani et al. 2011;Li et al. 2013b;Wang et al. 2017;Planck Collaboration et al. 2016), we discuss below the formation of clumps, enhanced gas and magnetic pressures, the pressure balance between clumps and feedback, and their consequences for the evolution of clumps and star formation in them, and the formation of bipolar H ii regions.
The expanding I-front from a deeply embedded H ii region in a filament becomes anisotropic, such that the flows along the dense filament-ridge will become sonic as they are obstructed, while they are supersonic in the low-density regions both below and above the ridge. As a result of the natural anisotropic distribution of material in the filament, the H ii region is lead to form bipolar bubbles (Bodenheimer et al. 1979;Fukuda & Hanawa 2000). Moreover, inclusion of B-fields in the filament would introduce an additional anisotropic pressure (Tomisaka 1992;Gaensler 1998;Pavel & Clemens 2012;van Marle et al. 2015), because ionized material flowing along the B-fields will be accelerated, while those flows in the direction perpendicular to B-fields will be hindered due to the Lorenz force. , based on MHD simulations, have shown that Bfields suppress the sweeping up of gas perpendicular to field lines. As the H ii region expands further into the cloud, gas and dust in the filament-ridge will be swept up, and as a result the accumulated material is led to form the dense clumps at the waist of H ii region. It should be noted here that B-field strength will also be continuously enhanced due to flux freezing, as is seen in the clumps (c.f. Section 4.4).
Gas and magnetic pressures, and so total clump pressure, will increase with time. As a result, at certain point in time, the enhanced clump pressure stops further expansion of the I-front into the clump (e.g., Ferland 2008, 2009, andreferences therein). For this to occur, a near pressure-equilibrium must be achieved between clumps and feedback, i.e., the total pressure within the clumps should be equal to or higher than the pressure imparted by the feedback from the H ii region. Below we test this hypothesis.
The pressure balance equation can be written as where the left-hand side (LHS) corresponds to the clump internal pressure, P clump = P B + P turb + P Tg (e.g., equation 14 of Miao et al. 2006), which is the combination of magnetic (P B ), turbulent (or non-thermal; P turb ), gas thermal (or kinetic; P Tg ) pressures. The combination of the last two components can be treated as the total molecular gas pressure (P mol ) in the clumps, i.e., P turb + P Tg = P mol . Molecular gas pressure (P mol ) is estimated using the following relations (e.g., Liu et al. 2017) and Using the C eff values measured in Section 4.3.2), the estimated T eff values are 188±10 K and 148±11 K for clumps 1 and 2. Finally, P mol is derived to be (13 ± 2) × 10 −10 dyn cm −2 for clump 1 and (2.7 ± 0.5) × 10 −10 dyn cm −2 for clump 2.
The right-hand side (RHS) in Equation 22 corresponds to feedback pressure, P fb = P Te + P rad , due to the combination of the thermally ionized medium (from electron temperature; P T e ) and radiation (P rad ) components. P clump1 and P fb1 are estimated to be 41×10 −10 and 8×10 −10 dyn cm −2 for clump 1. Similarly, P clump2 and P fb2 are estimated to be 4.4×10 −10 and 5.1×10 −10 dyn cm −2 for clump 2. These parameters, holding to the relations P clump1 > P fb1 and P clump2 P fb2 , suggest that clump 1 stops further expansion of the ionized region, whereas feedback pressure is nearly in equilibrium with the pressure in clump 2.
Our analyses show that magnetic pressure dominates over thermal pressure (at least in clump 1; Section 4.2), and that B-fields within the clumps situated at the waist of the H ii region can constrain the paths of Ifronts to blow away from the filament ridge. In RMHD simulations (Mackey & Lim 2011) initially stronger Bfields are shown to confine the photoevaporation flow into a bar-shaped, dense, ionized ribbon which shields the I-front. These features are observed in both the clumps of S201, as is clear from the PACS/70µm image shown in Figure 3. Therefore, combined contributions from the anisotropic expansion of the I-front, additional anisotropic pressure introduced by the B-fields in the primordial filament, and the enhanced B-fields in the clumps would result in the formation of bipolar H ii regions (e.g., Deharveng et al. 2015;Samal et al. 2018).
Furthermore, enhanced B-fields not only guide the ionized gas and aid the formation of bipolar H ii regions but also shield the clumps from erosion. These signatures imply that the enhanced B-fields and underlying turbulence at the clump centers counteract gravitational collapse and hence delay the evolution of the clumps. Based on NIR polarimetry, Chen et al. (2017) have found that the B-field strength in a shell, N4, has been enhanced due to magnetically frozen-in gas being swept up into the expanding shell. As a result, the fragmented clumps in N4 remain in a magnetically subcritical state, indicating that the B-field is the dominant force, stronger than gravity. Similarly, based on SCUBAPOL2 observations, Pattle et al. (2018) have probed B-fields in the denser parts of the 'Pillars of creation' in M16 and found that initially B-fields are swept up, becoming aligned parallel to the pillars and later due to gas compression the B-field become stronger; and hence govern the evolution and longevity of the pillars.

Limitations of the current study
Due to limited sensitivity (mean rms noise in 850µm Stokes I is ∼5 mJy/beam with a bin size of 12 ) achieved in our study and the relatively small field-of-view (3 diameter with uniform sensitivity) of SCUBAPOL2 observations, we could not probe B-fields in the clumps on the far side of the H ii region. We also note that because of the small number of measurements made in clump 2, a reliable B-field strength has not been derived from ACF analysis.
According to Mackey & Lim (2011), initial weak (18 µG) and intermediate-strength (53 µG) B-fields that are initially oriented perpendicular to the I-front are swept into alignment with the pillar during its dynamical evolution, consistent with B-field observations in M16 (Sugitani et al. 2007;Pattle et al. 2018). Based on SCUBAPOL observations towards S106, Vallée & Fiege (2005) suggested that at large scales B-fields are roughly oriented along the north-south direction around the bipolar H ii region, but close to central region near the IR star the B-fields are twisted into a toroidal morphology. Due to a lack of observations over an extended area, we could not examine these scenarios in S201.
Based on MHD simulations of the evolution of sheetlike clouds due to mechanical stellar feedback from a single massive star (due both to stellar winds and to a supernova explosion), Wareing et al. (2017) showed that B-fields tend to follow a bipolar bubble-like structure similar to B-fields observed in RCW57A using NIR polarimetry . In order to examine whether B-fields follow and connect the structures of the photodissociation regions (PDRs) of the clumps (this work) as well as bipolar cavity walls , further observations probing B-fields towards similar targets with SCUBAPOL2 are desirable.
MHD simulations focusing on the time evolution of B-fields, turbulence, gravity, and thermal energies, and their impact on the formation and evolution of clumps and bipolar H ii regions would also be valuable.

SUMMARY AND CONCLUSIONS
We have performed the dust continuum polarization observations at 850 µm and probed the B-fields in the deeply embedded massive clumps (clump 1 and clump 2) located at the waist of the bipolar H ii region S201 using the JCMT SCUBA-2/POL-2. In addition, we have utilized JCMT/HARP molecular lines ( 13 CO (3-2) and C 18 O (3-2)) and VLA 21 cm radio data, respectively, to quantify the turbulent and thermal energies in the region. In this work, we have derived various parameters such as B-fields, turbulence, gravity, and thermal pressures, and studied their interplay in the context of H ii region-influenced star formation. The following are the main findings of our study: 1. The morphological correlation between the orientations of B-fields and intensity gradients of ionized gas (based on the VLA 21 cm continuum) suggests that B-fields are compressed and bent by the expanding ionization fronts from the H ii region.
2. Compressed B-fields in clump 1 follow a bow-like morphology, while the degree of compression as well as of bending in B-fields is lower prominent in clump 2. The observed degree of compression, and of enhancement in the B-field strengths in clumps are in accordance with the amount of ionized medium interacting with the clump surfaces.
3. We have employed structure function (SF) and auto-correlation function (ACF) analyses to derive the ratio of turbulent to ordered B-field in the clumps. Using the velocity dispersion from C 18 O data and column density derived from the POL2 Stokes I map, B-field strengths have been estimated by using the modified Davis-Chandrasekhar-Fermi relations. B-field strengths are found to be 266±32 µG (from ACF) for clump 1 and 65±6 µG (from SF) for clump 2.
4. We find that magnetic pressure dominates in clump 1 and thermal pressure (due to ionized gas) dominates in clump 2. The amount of turbulent pressure lies between those of B-fields and those of ionized gas thermal energies in both clumps.
5. The comparison between clump internal pressure (magnetic, gas thermal, and non-thermal or turbulence) and feedback pressure (ionized gas thermal and radiation) imply that clump 1 has stopped further expansion of the H ii region, while clump 2 maintains a near equilibrium with the feedback pressure.
6. Virial analyses suggest that clump 1 is bound by its gravity and may collapse once it becomes unstable, whereas clump 2 is gravitationally unbound. We suggest that clump 2 may become bound in future if it progressively accumulates sufficient mass.
7. Critical mass ratios reveal that the combined contribution from gas thermal energy, turbulence, and B-fields is not sufficient to counteract the gravity and hence is under collapse, whereas an opposite situation prevails in clump 2, in that its stability against collapse is maintained by these three factors. These results are consistent with the observed distribution of young stellar objects (YSOs), which suggests that star formation is ongoing in clump 1, while there is no star formation in clump 2.
8. Feedback from H ii region has the following consequences -(a) causes the formation of clumps in the filament ridge, i.e., at the waist of the H ii region, (b) enhances the B-field strength in the clumps and injects turbulence into the clumps, and (c) eventually the enhanced B-fields will be able to shield the clumps from erosion and so govern their stability, guiding the expanding I-fronts to be blown away from the filament ridge, and aiding in the formation of bipolar H ii regions.

ACKNOWLEDGMENTS
We thank the anonymous referee for useful suggestions. This work is supported by the National Natural Science Foundation of China (NSFC) grant Nos. 11988101 and 11725313, and the International Partnership Program of Chinese Academy of Sciences grant No. 114A11KYSB20160008. This work also supported by Special Funding for Advanced Users, budgeted and administrated by Center for Astronomical Mega-Science, Chinese Academy of Sciences (CAMS). C.E. and S.P.L. acknowledge the support from the Ministry of Science and Technology (MOST) of Taiwan with grant MOST 106-2119-M-007-021-MY3. C.E. thanks David Berry for the help with JCMT SCUBA2-POL2 data reduction and analyses. A.Z. thanks the support of the Institut Universitaire de France. D.K.O. acknowledges the support of the Department of Atomic Energy, Government of India, under project NO. 12-R&D-TFR-5.02-0200. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; Center for Astronomical Mega-Science. Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada.