Understanding the links among magnetic fields, filament, the bipolar bubble, and star formation in RCW57A using NIR polarimetry

The influence of magnetic fields (B-fields) in the formation and evolution of bipolar bubbles, due to the expanding ionization fronts (I-fronts) driven by the Hii regions that are formed and embedded in filamentary molecular clouds, has not been well-studied yet. In addition to the anisotropic expansion of I-fronts into a filament, B-fields are expected to introduce an additional anisotropic pressure which might favor expansion and propagation of I-fronts to form a bipolar bubble. We present results based on near-infrared polarimetric observations towards the central $\sim$8'$\times$8' area of the star-forming region RCW57A which hosts an Hii region, a filament, and a bipolar bubble. Polarization measurements of 178 reddened background stars, out of the 919 detected sources in the JHKs-bands, reveal B-fields that thread perpendicular to the filament long axis. The B-fields exhibit an hour-glass morphology that closely follows the structure of the bipolar bubble. The mean B-field strength, estimated using the Chandrasekhar-Fermi method, is 91$\pm$8 {\mu}G. B-field pressure dominates over turbulent and thermal pressures. Thermal pressure might act in the same orientation as those of B-fields to accelerate the expansion of those I-fronts. The observed morphological correspondence among the B-fields, filament, and bipolar bubble demonstrate that the B-fields are important to the cloud contraction that formed the filament, gravitational collapse and star formation in it, and in feedback processes. The latter include the formation and evolution of mid-infrared bubbles by means of B-field supported propagation and expansion of I-fronts. These may shed light on preexisting conditions favoring the formation of the massive stellar cluster in RCW57A.

Though bipolar bubbles are the natural outcome of anisotropic expansion of ionizing fronts from Hii regions hosting massive O/B-type star(s) located in filaments (Minier et al. 2013;Deharveng et al. 2015;Zhang et al. 2016), details involved in their formation and evolution processes remain poorly understood. According to 2D (Bodenheimer et al. 1979) and 3D (Fukuda & Hanawa 2000) hydrodynamic simulations, the evolution of an Hii region in a filamentary cloud induces supersonic and subsonic I-fronts along the minor and major axes of the filament, respectively. This results the distribution of low and high density ma-terial along the minor and major axes. Thus, the I-fronts experience more hindered flows along the major axis than along the minor axis. As a consequence of this anisotropic expansion, a bipolar bubble will be formed in a filament (Figure 1 of Deharveng et al. 2015). Though a few studies (e.g., Minier et al. 2013;Deharveng et al. 2015) were devoted to the study of bipolar bubbles, there exists no study focused on exploring the connection between the B-fields anchored in the filaments and their role in the formation and evolution of bipolar bubbles. Fukuda & Hanawa (2000) did include B-fields (parallel to the filament long axis) in their simulations, but they studied only the influence of B-fields on shape, separation, and formation epoch of the cores that formed along the filament.
B-fields are believed to guide the contraction of cloud material to form filaments in molecular clouds, and thus play a crucial role in controlling cloud stability and collapse, fragmentation into cores, and internal star formation (Pereyra & Magalhães 2004;Alves et al. 2008;Li et al. 2013;Planck Collaboration et al. 2014;Franco & Alves 2015). A few studies (e.g., Pereyra & Magalhães 2007 towards the IRAS Vela Shell and Wisniewski et al. 2007 towards supershell NGC 2100 in the LMC) suggest the importance of B-fields in the dynamical evolution of the expanding shells driven by O/B-type stars. The interplay between the expansion of an ionized nebula and the influence of B-fields has been the subject of several studies (e.g., Pavel & Clemens 2012;Santos et al. 2012Santos et al. , 2014Planck Collaboration et al. 2015;Kwon et al. 2010Kwon et al. , 2011Kwon et al. , 2015. These studies hint that B-fields permeated through filamentary molecular clouds would also influence the expanding I-fronts or outflowing gas, from an Hii region created in the filament. RCW57A (also known as NGC 3576, G291. 27-0.70, or IRAS 11097-6102) is an Hii region associated with a filament and bipolar bubble and is located at a distance of 2.4 -2.8 kpc (Persi et al. 1994;de Pree et al. 1999). We adopt 2.4 kpc, which is within uncertainties of both the kinematic and spectroscopic determinations (see Persi et al. 1994). Figure 1 depicts the overall morphology of RCW57A. It contains optically bright nebulosity with several dark globules and luminous arcs (Persi et al. 1994). It is one of the massive star -Tricolor image of the RCW57A region made using WISE 4.6 µm (red), 2MASS K s -band (green), and DSS2 R-band (blue) images. Various contours, extracted from Purcell et al. (2009), are overplotted. Yellow contours enclosing white arrows represent the bipolar bubble-like structure observed in mid-infrared (MSX, Spitzer, and WISE) images. This feature appears as an extended widened loop in the R-band, as shown with an extended yellow contour in the Northern part. Thick cyan contour represents the 3.4 cm radio free-free emission from ionized gas, delimiting the extent of the Hii region. The white contours represent the 1.1 mm dust continuum emission (Hill et al. 2005), showing an elongated dusty filament from NE to SW which crosses the H ii region. The extent of the molecular cloud, from 13 CO(1 -0) line emission (Purcell et al. 2009), is depicted with magenta outline (the horizontal line at the bottom of the magenta outline is the limit of the mapped area). Infrared sources (IRS; Frogel & Persson 1974) and water/methanol maser sources (cf. Purcell et al. 2009) are shown with open black squares and green crosses, respectively. Possible directions of the outflowing gas (cf. de Pree et al. 1999) are shown with white arrows. Red square box denotes the area of 7. 7 × 7. 7 observed with SIRPOL.
forming regions in the southern sky, hosting an Hii region (cyan contour) embedded in a filament (white contour) from which a widely-extended bipolar bubble (yellow contours) is emerging. A deeply embedded near-IR cluster, consisting of more than 130 young stellar objects (YSOs) is associated with this region (Persi et al. 1994). The observed ratios of the infrared fine-structure ionic lines (Ne ii, Ar iii, and S iv) (Lacy et al. 1982) indicate that at least eight O7.5V stars are necessary to account for the ionization of the region. However, even these stars may not be sufficient to account for the Lyα ionizing photons inferred from radio data (Figuerêdo et al. 2002;Barbosa et al. 2003;Townsley 2009). Based on the newly discovered cluster of stars, using X-ray data, Townsley et al. (2014) suggested that an additional cluster of OB stars might be deeply embedded that were not known before because of heavy obscuration. This cluster is located slightly NW of the center of the near-IR cluster. The 10 µm map (cf., Frogel & Persson 1974), reveals the presence of five infrared sources (IRS; black squares in Figure 1) near the center of the Hii region. These, together with water and methanol maser sources (green crosses in Figure 1) distributed along the filament, are indicative of active ongoing star formation in RCW57A (Frogel & Persson 1974;Caswell 2004;Purcell et al. 2009). Therefore, RCW 57A is an ideal target to investigate the morphological links among filaments, bipolar bubbles, and B-fields so as to understand the star formation history.
Here, near-infrared (NIR) polarimetric observations towards the central region of RCW57A have been carried out to map the plane-of-sky B-field geometry in the star forming region. When unpolarized background starlight passes through an interstellar cloud, dust particles partially aligned by a B-field linearly polarize the light by a few percent. Virtually all possible mechanisms responsible for the dust grain alignment with respect to the B-fields (Davis & Greenstein 1951;Andersson et al. 2015;Lazarian et al. 2015) yield directions of the polarized light tracing the average B-field orientation projected on the plane-of-sky.
The aim of this paper is to understand whether (a) the B-field structure in the molecular cloud guides feedback processes, such as expansion and propagation of outflowing gas and ionization fronts, which lead to the formation of bipolar bubbles, or (b) feedback processes regulate the resulting B-field structure. The B-field is active in the former scenario, while it is passive in the latter.
The structure of this paper is as follows. Section 2 describes the observations and data reduction. Detailed analyses are presented in Section 3, the main purpose of which is to identify and remove the YSOs from the list of fore-ground and background stars by utilizing NIR and mid-infrared (MIR) colors, NIR polarization measurements, and polarization efficiencies. Results are presented in Section 4. This section presents estimates of the B-field strength, and detailed comparisons between B-field pressure to turbulent and thermal pressures. A possible evolutionary scenario, based on the morphological correspondences among the filament, bipolar bubble, and B-fields is discussed in Section 5. Furthermore, based on the scenario that best characterizes our observed B-fields in RCW57A, we discuss the possible preexisting conditions that might favored the formation of star cluster in RCW57A. Conclusions are summarized in Section 6.
One set of observations consisted of 10 s exposures at four HWP position angles (0 • , 22.5 • , 45 • and 67.5 • ) at 10 dithered sky pointings. Such sets of 4×10 images were repeated towards the same sky coordinates to increase signal-to-noise. Sky frames were also obtained between target observations. The total integration time was 400 s per HWP angle. The average seeing was 1. 54 (J), 1. 44 (H) and 1. 31 (K s ).
Master flats were created by utilizing the evening and twilight flat-field frames on the same night of the observations. We processed the data using the dedicated data reduction pipeline 'pyIRSF' 2 . This pipeline used the raw data from the target, sky, dark, and flat field frames as inputs. The data reduction tasks included dark subtraction, flat-field correction, median sky subtraction, frame registration, and averaging. The final products of this pipeline were the average combined four intensity images (I 0 , I 22.5 , I 45 , and I 67.5 ) corresponding to the four positions of the HWP.

Aperture photometry of point sources
We performed aperture photometry of pointlike sources on the intensity images for the HWP angles (I 0 , I 22.5 , I 45 , and I 67.5 ) of each of the JHK s -bands using iraf 3 and the idl Astronomy Library (Landsman 1993). Point sources with peak intensities above 5σ (where σ is the rms uncertainty of the particular image pixel values) of the local sky background were detected using DAOFIND. Aperture photometry was performed using the PHOT task of the DAOPHOT package in IRAF. An aperture radius was chosen to be nearly the FWHM of the point-like sources, i.e., 3.4, 3.2, and 2.9 pixels in the J-, H-, and K s bands, respectively. The sky annulus inner radius was set to be 10 pixels with a 5 pixel width. A source that matched in center-to-center positions to within 1 pixel radius (i.e., 2 pixel diameter equivalent to ∼ 1 ) in all four HWP images was judged to be same star.
The Stokes I intensity of each point-like source was calculated using

Aperture polarimetry
Polarimetry of point-like sources (aperture polarimetry) was performed on the combined intensity images at each HWP angle. Extracted intensities were used to estimate the Stokes parameters of each star using The aperture and sky radii were the same as those used in the aperture photometry on total intensity (I) images. To obtain the Stokes parameters (Q and U ) in the equatorial coordinate system, an offset rotation of 105 • (Kandori et al. 2006;Kusune et al. 2015) was applied. We calculated the degree of polarization P and the polarization angle θ as follows: The errors in P and θ were also estimated by propagating the errors in Stokes parameters according to equations 4 and 5. Since the polarization degree, P , is a positive definite quantity, the derived P value tends to be overestimated, especially for low S/N cases. To correct for this bias, we calculated the de-biased P using P db = √ P 2 − δP 2 (Wardle & Kronberg 1974), where δP is the error in P . It should be noted here that, hereafter, we consider P db as P . The absolute accuracy of the offset polarization angle for SIRPOL was estimated to be less than 3 • (Kandori et al. 2006). The polarization efficiencies of the wire grid polarizer are 95.5%, 96.3%, and 98.5% at the J-, H-, and K s -bands, respectively, and the instrumental polarization is less than 0.3% over the field of view at each band (Kandori et al. 2006). Due to these high polarization efficiencies and low instrumental polarization, no additional corrections were made to the data. To obtain the astrometric plate solutions, we matched the pixel coordinates of the detected sources and the celestial coordinates of their counterparts in the 2MASS Catalog (Skrutskie et al. 2006), and applied ccmap, ccsetwcs and cctran tasks of the iraf imcoords package to the matched lists. The rms error in the returned coordinate system was ∼0.03 . However, since the 2MASS astrometry system is limited about to 0.2 , the stellar positional uncertainties should take that into account. A total of 919 stars had polarization detection in at least one band with P/σ P ≥ 1.

Complete data sample
The polarimetric data for the 919 stars were merged with the photometry of the 702 stars, and the resultant data are presented in Table 1. Star IDs are given in column 1. The photometric data of the stars, those do not have them from IRSF, are extracted from 2MASS (Skrutskie et al. 2006) and are denoted with an asterisk along with star IDs in column 1. The coordinates, aperture polarimetric and photometric results of each source are presented in columns 2 -12. Column 13 reports our classifications of the stars, as described in section 3 below. Though we have presented all the stellar data satisfying the P/σ P ≥ 1 criterion in Table 1  Black dots are the total sample in the entire region. The blue and red filled circles represent Class I and Class II sources, respectively in the RCW57A region. The black circles surrounding some of the points are three Class I and thirty two Class II sources having 2σ polarization measurements. The slanted dotted line denotes the boundary between Class I and Class II YSOs. The average error in the colors is shown.

Identification of YSOs and foreground stars
The B-field inside the molecular cloud can be probed by using polarimetry of background stars whose light becomes differentially extincted and linearly polarized by the dust in the cloud. Hence, a sample of background stars needs to be identified among all the stars with measured polarizations. For this purpose, we identify and exclude the following sources: YSOs (Class I, Class II, HAe/Be stars and NIR excess sources) based on MIR and NIR color-color diagrams, stars with Lband excess using the data from Maercker et al. (2006), stars projected against the nebulous region where their polarizations might be contami-nated by scattering (especially at the cluster center), the stars with non-interstellar polarization origin due to reflected light from their circumstellar disks/envelopes (whose properties are not discerned based on MIR and NIR colors), and the stars that have suffered depolarization effects due to multiple B-field components (associated with multiple dust layers) along their lines of sight. The latter two are identified using polarization efficiency diagram.
Below, we present MIR and NIR two color diagrams (TCDs) and polarization efficiency (P versus H − K s ) diagram to identify and exclude the above mentioned stars with intrinsic polarization. Finally, the remaining sample, that are free from intrinsic polarization, is sub-categorized further into foreground (FG) and background (BG) stars, based on their NIR colors and polarization characteristics.

YSOs identified by their mid-infrared colors
Since the RCW57A region is still enshrouded in its natal molecular cloud, as evident from the 13 CO(1 -0) map (magenta contour in Figure 1 and red background in Figure 7), YSOs in the region can be deeply embedded. Therefore, the Spitzer MIR observations were used to probe deeper insight into the embedded YSOs. These occupy distinct regions in the Spitzer IRAC color plane, which makes MIR TCDs a useful tool for the classification of YSOs. Since 8.0 µm data are not available for the region, we used [[K s ] − [3.6]] 0 and [[3.6] − [4.5]] 0 (cf. Gutermuth et al. 2009) to identify deeply embedded YSOs, as shown in Figure  3. The zones of Class I and Class II YSOs, based on the color criteria of Gutermuth et al. (2009), are also depicted. Three sources were identified as Class I and thirty two sources as Class II and so were excluded from analyses. Identified Class I and Class II sources are mentioned in column 13 of Table 1.

YSOs identified by their near-infrared colors
NIR colors are also important tools to classify YSOs in star forming regions. Figure (Bessell & Brett 1988), respectively. The dotted blue line indicates the locus of unreddened CTTSs. The parallel blue dashed lines are the reddening vectors drawn from the tip (spectral type M4) of the giant branch (left reddening line), from the base (spectral type A0) of the MS branch (middle reddening line), and from the tip of the intrinsic CTTS line (right reddening line) with the extinction ratios A J /A V = 0.265, A H /A V = 0.155, and A K /A V = 0.090 (Cohen et al. 1981). The blue plus symbols on the reddening vectors show increments of A V by 5 mag each. The sources located in the "F" region could be either reddened field stars or Class II and Class III sources with small NIR excesses. The sources distributed in the "T" and "P" regions are considered to be mostly Classical T-Tauri stars (CTTSs; Meyer et al. 1997) or Class II sources with relatively large NIR excesses and likely Class I sources, respectively (for details see Pandey et al. 2008;Chauhan et al. 2011a,b). Typical errors in the plotted colors are comparable to the size of the symbols.  (Maercker et al. 2006). Two of the stars found to be Class I also have L-band excess, and similarly 3 stars found to be Class II also have L-band excess. Star identified as having NIR excess, Lband excess, and nebulous background are noted in column 13 of Table 1.
The remaining 342 stars that appear to be free from intrinsic polarization and distributed in the zone "F" in Figure 4 are (i) lightly reddened foreground field stars (FG), (ii) moderately reddened cluster members and highly reddened background stars (hereafter both cluster members and background stars referred to as BG), and (iii) weak line T-Tauri stars (WTTSs) or Class III sources. The locations of the WTTSs in the NIR TCD overlap with reddened background stars. Generally WTTSs exhibit almost negligible NIR-excess, as their disks will have been evaporated already or very optically thin (eg., Adams et al. 1987;Andrews & Williams 2005;Cieza et al. 2007). Therefore, it is reasonable to assume that these WTTSs are not intrinsically polarized by scattered light from the negligible amount of circumstellar material.
In order to estimate the level of FG contamination, a control field (α = 11 h 18 m 12. s 22, δ = −61 • 59 38. 92 [J2000]) of angular size 7. 7×7. 7 located ∼60 from the RCW57A region was chosen. Comparing the distributions of [H − K s ] colors of the control field with colors in the RCW57A region shows that the field stars to have [H − K s ] < 0.15. Therefore, the 108 RCW57A stars with [H − K s ] < 0.15 and lying in the 'F' region of the NIR TCD are judged to be FG star candidates and remaining the 234 stars with [H − K s ] ≥ 0.15 and lying in the 'F' region are judged to be BG star candidates.

Final confirmed FG and BG stars based on polarization efficiency diagrams
Polarization efficiency is defined as the ratio of polarization degree to interstellar reddening, such as P (λ)/E(λ 1 − λ 2 ). Observed polarization efficiencies of dust grains are found to vary from one line of sight to the another. Reddening (E(B−V )) and polarization in V -band (P V ) are correlated and can be represented by the observational upper limit relation: (Serkowski et al. 1975). Polarization measurements that exceed this upper limit relation are generally attributed to the intrinsic polarization (non-magnetic polarization) caused by scattered light from dust grains in circumstellar disks or en-  (Serkowski et al. 1975).
We also seek to identify and exclude the stars that suffer depolarization due to their starlight's passing through multiple B-field components along the line of sight. For optically thin clouds, P (λ) is expected to be proportional to the column density if the field structure is uniform. In contrast, presence of multiple B-fields with different orientations along a line of sight will result depolarization (Martin 1974) and hence weakens the correlation between polarization and extinction (eg., Eswaraiah et al. 2011Eswaraiah et al. , 2012. Therefore, polarization efficiency diagrams (eg., Kwon et al. 2011;Sugitani et al. 2011;Hatano et al. 2013) will help identify stars with either intrinsic polarization or depolarization. Figure 5 show P versus [H−K s ] diagrams of FG and BG candidate stars (Section 3.2) with H-band data. These diagrams contain no known YSOs as we excluded them in the above Sections 3.1 and 3.2. The BG candidate stars show an increasing trend in polarization in accordance with an increasing [H − K s ] and well distributed below the upper limit interstellar polarization efficiency relation (gray line), Serkowski et al. 1975, also see Hatano et al. 2013 for detailed description on deriving polarization efficiency relations in NIR wavelengths.). To find the stars with non-magnetic polarizations (non interstellar origin of polarizations) and depolarization affect, first, we fit a straight line to polarization versus extinction for the BG stars with σ P (H) < 0.3% (encircled filled circles). Second, we drew dashed lines with slope of twice (upper dashed line) and half (lower dashed line) of the fitted slopes. Third, we identified the stars falling above and below these lines as probable BG stars having influenced by non-magnetic origin of polarization and depolarization affect, respectively. The remaining stars falling between the two dashed lines are ascertained as confirmed BG stars. Therefore, the final number of confirmed BG stars was found to be 178 in H-band. Additional 64 candidate BG stars falling above and below the dashed lines (plus symbols) are mentioned as probable BG stars with intrinsic polarization or depolarization in Table 1. FG candidate stars (open circles) exhibit two distributions and those with P (H) > 3% (plus symbols on open circles) may be affected by additional polarization caused by scattering or simply with larger uncertainties. We consider them as probable FG stars. Remaining 97 stars with P (H) ≤ 3% are considered as confirmed FG stars in H-band. These, confirmed and probable FG stars with P (H) > 3% are mentioned in Table  1.
In the further analysis, we used only the confirmed 97 FG and 178 BG stars with Hband data. The fitted slope of P (H)/E(H − K s ) = 5.5±0.8% mag −1 (thick line) suggest that polarization efficiency of dust grains in the star forming cloud RCW57A are relatively higher than those of other Galactic line of sights (e.g., Hatano et al. 2013;Kusune et al. 2015).

Separation of foreground and background magnetic fields
FG stars towards any distant cluster are, generally, less polarized and extincted than BG stars (Eswaraiah et al. 2011(Eswaraiah et al. , 2012Pandey et al. 2013). This is because the dust layer(s) between the observer and the FG star contains fewer aligned dust grains, and because the degree of polarization and extinction are linearly correlated (Aannestad & Purcell 1973;Serkowski et al. 1975;Hatano et al. 2013). If the B-field orientation of the foreground dust layer is sufficiently different from that of the cloud, then the foreground stars would exhibit different polarizations (P and θ) from those of the background stars. Therefore, foreground and background stars can also be distinguished based on their polarization characteristics (P vs θ or Q versus U ), as well as with their NIR colors (Medhi et al. 2008(Medhi et al. , 2010Eswaraiah et al. 2011Eswaraiah et al. , 2012Medhi & Tamura 2013;Pandey et al. 2013;Santos et al. 2014).
Figures 6(a)-(d) depict that the 97 confirmed FG and 178 confirmed BG stars are well separated in their polarization characteristics and [H − K s ] colors. Foreground stars have [H − K s ] < 0.15, P (H) < 3% with a peak around 1.0%, and a single Gaussian distribution of θ(H) (Gaussian mean of 65 • with a standard deviation of 14 • ). Background stars have [H − K s ] ≥ 0.15, P (H) between 0.5 -12% with a peak ∼ 2% and an extended tail towards greater P values and widely distributed θ(H) (range from 0 • to 180 • ) indicating the presence of complex B-field structure in RCW57A with a prominent peak at θ(H) ∼163 • . About 55% of stars are distributed between 110 -180 • , whereas the remaining 45% of stars are distributed between 0 -110 • . This suggests that the plane-of-sky B-field structure of RCW 57A is dominated by a component with θ(H) ∼163 • , which is nearly orthogonal to the foreground B-field component (65 • ) as well as being nearly orthogonal to the position angle of the major axis of the filament (∼60 • ). Figure 6(d) also reveals that the distributions of θ(H) of the confirmed FG and BG stars In order to account for the influence of FG polarization on the polarization angles of the confirmed BG stars, below we performed the FG subtraction. First, the polarization measurements of confirmed 97 FG stars were converted into Stokes parameters. Second, weighted mean FG Stokes parameters were computed and found to be Q F G = −0.332±0.008% and U F G = 0.845±0.007% (with standard deviation of 0.64% in Q F G and 0.70% in U F G ) and were subtracted vectorially from those of each confirmed BG star. Third, FGsubtracted Stokes parameters of the BG stars were then converted back to θ(H) FG−corrected . The Gaussian mean and standard deviation, for the distribution of offset polarization angles of the BG stars (θ(H) − θ(H) FG−corrected ), are found to be 1 • and 7 • , respectively. Hence, there was no obvious change in the FG corrected B-field morphology inferred from the BG stars. Therefore, we ignored the FG contribution to the BG polarization measurements. Figure 7 shows the vector map of the H-band polarization measurements of the 178 confirmed BG stars. The polarization vector map reveals a systematically ordered B-field that is configured into an hour-glass morphology which closely resembles the structure of the bipolar bubble.

Magnetic field geometry
The direction of expanding I-fronts or outflowing gas is revealed based on a large NS velocity gradient in a radio recombination line observations (de Pree et al. 1999). These high velocity outflows signatures are further traced in 13 CO(1 -0) position-velocity diagrams ( Figure 11; shown below). Shocked gas is visible at the SE and NW ends of the white arrows in Figure 7 (also see Figure 1), highlighting the outflow-cloud interaction. Townsley et al. (2011, see their Figures 3 and 4) have also showed that outflows from the deeply embedded protostars are responsible for the soft X-ray emission observed in the SE part of the cloud. These outflows have created the cavities depicted with enclosed yellow contours ( Figure  1). Interestingly, the directions of the expanding I-fronts and/or outflowing gas, as depicted with white arrows (Figure 7), follow the B-field orientations.
B-field are compressed along the rim of BRC (located inside of white square box in Figure 7) and more details of which are addressed in Appendix B.

Extracting essential parameters for estimating magnetic field strength
To test whether B-fields play an active role in the formation of bipolar bubbles, it is essential to estimate the B-field strength by using the Chandrasekhar & Fermi (1953) relation (hereafter CF Method). For this purpose, below we determine the dispersion in θ(H) for the background stars, the gas velocity dispersion using CO data, and the dust and gas volume density using Herschel data.

Dispersion in polarization angles: σ θ(H)
The presence of multiple components of B-fields in RCW57A can be seen from the H-band polar-ization vector map overlaid on the 13 CO(1 -0) (Purcell et al. 2009) total integrated intensity image ( Figure 8) and the histogram of θ(H) (continuous histogram in Figure 6(d)). Each component follows its own spatial distribution and exhibits its own θ(H) distribution and dispersion. For this reason, we divided the observed area into five regions, named, A, B, C, D, and E, as shown in Figure 8 and Table 2. The spatial extent of each region is visually chosen in such a way that it should not include multiple components of θ(H). This criteria constrains the dispersion of θ(H) to be less than 25 • (Ostriker et al. 2001), so as to apply the CF method to estimate the B-field strength. Some portions are not included because they either have polarization data but not CO data (the area right to the C region), or have CO data but has multiple θ(H) distributions (area slightly above the E region), or contains few randomly oriented vectors (area distributed between D and E regions). We believe that the bias involved in the selection of five regions may not have significant impact on the final scientific results.
A Gaussian function was fit to the θ(H) distribution for each region, as shown in Figure 9. The fitted means and dispersions in polarization angles, along with their errors, are listed in columns 7 and 8 of Table 2, respectively. Mean θ(H) (=µ θ(H) ) values of regions A, B, C, D, and E regions found to be 26 • , 164 • , 153 • , 132 • , and 74 • respectively. Dispersions in θ(H) (=σ θ(H) ) of five regions lie between 10 -16 • .

Velocity dispersion: σ V LSR
To derive the velocity dispersion in regions A, B, C, D and E, we have used 13 CO(1 -0) emission line mapping (Figure 8) performed using the Mopra telescope by Purcell et al. (2009). The observed 13 CO(1 -0) data has a velocity resolution of 0.4 km s −1 and a spatial resolution of 40 . CASA software was used to extract the mean velocity (V LSR ) versus brightness temperature (T ) spectrum of each region, which are plotted in Figure 10 using small filled circles. All the regions exhibit at least two velocity components, with clear asymmetric spectra towards higher velocities. These asymmetric, spatially extended, spectra correspond to outflowing gas from the embedded YSOs, expanding I-fronts, or both. Each of these spectra (T vs. V LSR ) were then fitted by Fig. 7.-H-band polarizations (red vectors) of the 178 confirmed BG stars overlaid on the tricolor image constructed from the 13 CO(1 -0) total integrated emission map (red) (Purcell et al. 2009), Spitzer IRAC Ch2+Ch3 combined image (green), and SIRPOL J + H + K s -band combined image (blue). Reference vector with P = 6% and θ = 90 • is plotted. Thick black contour represents SIMBA 1.2-mm dust emission (Hill et al. 2005) with a flux level of ∼3 Jy/beam. Thin black contours denote ATLASGAL 0.87-mm thermal dust emission ranging from 2 Jy/beam to 22 Jy/beam (beam size of 18. 2) with an interval of 1 Jy/beam. Cyan contour represent the extent of the Hii region, as traced in 3.4-cm radio free-free emission (de Pree et al. 1999). Yellow contours delineate the morphology of the bipolar bubble. Magenta contour depicts the extent of the molecular cloud, as traced by 13 CO(1 -0) emission using Mopra telescope (Purcell et al. 2009). White square box shows the location of a Bright Rimmed Cloud. White arrows represent possible orientations of expanding I-fronts or outflowing gas. Green plus marks and black squares distributed along the filament are the water/methanol masers and IRS sources, respectively. Blue diamonds correspond to the seven massive cores (namely, S1-M1, S1-M2, S3-M4, S3-M5, S3-C3, S4-M6 and S5-M6) identified by André et al. (2008). All the information on various contours, and the location of masers and IRS sources is extracted from Purcell et al. (2009), by permission. Yellow thick lines represent 'Cut 1' and 'Cut 2' along the ionization front or outflowing gas and cloud main axis, respectively, as described in the text and used as the basis for Figure 11. A reference length scale with 2 corresponding to 2.2 pc is shown as the horizontal white, labeled line. a multi-Gaussian function using the custom IDL routine 'gatorplot.pro' 4 . Details regarding this multi-Gaussian fitting are described in Appendix C.
The resulting multi-Gaussian fitting parameters of each spectrum are presented in Table 3 and are plotted with black lines in Figure 10. The resulting combined spectrum of each region, the sum of all fitted multi-Gaussians (orange line), closely match the observed spectrum. The difference between observed and fitted spectra, the residuals, are shown with open squares and are closely distributed around zero T . Among the multiple velocity components of each region, we attribute the one related to the peak of the spectrum to the turbulent cloud component of that region (red Gaussian curve). These are used to estimate the B-field strength. The remaining Gaussian velocity components (black lines) are ascribed to the expanding 4 http://www.astro.ufl.edu/~warner/GatorPlot/ I-fronts and outflowing gas. These values are also given in column 9 of Table 2.
The dashed line drawn at V LSR = −26.5 km s −1 (Figure 10) corresponds to the line center velocity over the areas covered by ABCD regions of RCW57A. This V LSR = −26.5 km s −1 component closely matches center of the red Gaussian components of all regions except for region E. However, based on the dense gas tracers (NH 3 , N 2 H + , and CS), Purcell et al. (2009) witnessed that the peak V LSR of the filament remain constant around −24 km s −1 (see their Figure 13), which is slightly different from the V LSR = −26.5 km s −1 . This is mainly because of the following reasons: (a) CO is optically thick and hence traces the less dense outer parts of the cloud, and (b) the ABCD regions, located slightly away from the cloud center, contain less dense gas that is influenced by the Hii region.
Below, we try to resolve and quantify the contributions from turbulent cloud component, outflowing gas, and expanding I-front using positionvelocity (PV) diagrams. PV diagrams are useful diagnostics of gas kinematics in a star forming molecular clouds (e.g., Purcell et al. 2009;Zhang et al. 2016). Figures 11(a) and (b), show the P V diagrams corresponding to the two cuts shown and labeled as '1' and '2', respectively, in Figure 8. The cut '1' is along the orientation of expanding ionization front or outflowing gas (parallel to the minor axis of the filament) and the cut '2' is along the major axis of the filament.
The velocity range over the whole region, based on the extent of the contours shown in Figures  11(a) and (b), spans −28.5 to −16.5 km s −1 a width ∆V 12 km s −1 . This implies that the entire region is dynamically active under the influence of Hii region. However, among all parts of the cloud that we observed, one of the most quiescent clumps S4 (André et al. 2008;Purcell et al. 2009) is situated to SE part of the filament. This clump is situated at ∼250 and at ∼ −26.5 km s −1 in Figure 11(b). Dust temperature T d (∼19 -24 K; André et al. 2008) and kinetic temperature T k (∼15 -25 K; Purcell et al. 2009) of the cores in S4, suggest that this clump is dynamically quiescent as being far from the Hii region. Gas velocity around this clump spans over ∼ −25 to ∼ −28 km s −1 (centered at ∼-26.5 km s −1 ) and a width of 2 -3 km s −1 (Figure 11(b)). Therefore,  Table 3. According to the collect and collapse model of triggered star formation (see Torii et al. 2015, and references therein), the velocity separation due to expanding ionization front based on the P V diagrams to be ∼4 km s −1 . Therefore, Observed data points are shown with filled dots. The spectrum of each region was fit with multiple Gaussians shown with black lines. The combined spectrum is shown with a thick orange line. The mean velocity corresponding to the RCW57A region is drawn at −26.5 km s −1 using a dashed line.
the outflowing gas may have remaining contribution of ∼5 -6 km s −1 in addition to the cloud turbulent component of ∼2 -3 km s −1 .
Below we decompose the relative contribution of various velocity components in each region. In Figure 11(a), there exist two prominent redshifted, high-velocity components with peak intensities at −20 km s −1 and −23 km s −1 . These components are concentrated close to the cloud center. While the former component is attributed to the outflowing gas, the latter one is to the expanding I-front. The signatures of both the components prevail in A, B, C, and D regions. However, as shown in the bottom panel of Figure  10, the E region seems to be less influenced by the above two components, but has other two prominent components: (a) ∼ −25 km s −1 is related to the expanding I-front, and (b) ∼ −23 km s −1 is attributed to the turbulent cloud region. Moreover, spatially, the E region, in comparison to the others, is located far from the area influenced by the outflowing gas as depicted with white arrows (cf. Figure 7). This implies that the E region, containing a BRC, has a contribution from the expanding I-fronts alone, while the rest of the regions had the influences of both expanding I-fronts and outflowing gas. Therefore, the gas velocity dispersion due to cloud turbulent component (red lines in Figure  10) is well decomposed from other components using multi-Gaussian spectrum fitting.

Number density: n(H 2 )
Herschel PACS (Poglitsch et al. 2010) 160 µm and SPIRE (Griffin et al. 2010) 250 µm, 350 µm, and 500 µm images of RCW57A were used to construct a dust column density map by using the modified black body function (details are given in Appendix D). The resolution of the final column density map is 36.4 . To obtain the number densities, we further used Plummer-like profile fitting on column density profiles to obtain the volume density map as described below.
To construct the column density (N (H 2 )) profiles, first we manually identify the cloud spine using the pixels that have maximum column densities and distributed along the filament. Thus identified spine, depicted with green squares in Figure 12(c), divides the entire cloud region into two parts NW and SE. Second, for each spine element, we identify the N (H 2 ) profile elements distributed on a line that bisects the filament long axis and passes through that spine element. The dimensions of spine elements as well as N (H 2 ) profile elements are chosen as 36.4 × 36.4 . To achieve the Nyquist sampling, separation between adjacent elements was chosen as 18.2 . Third, we have constructed N (H 2 ) profiles as a function of projected distances from the filament spine elements towards both NW and SE regions, and are plotted using gray circles in Figures 12(a) and (b), respectively. Every point and its uncertainty, in this plots, corresponds to the mean and standard deviation of column density values in a particular N (H 2 ) profile element. The mean column density profiles, using distance wise averaged column densities of the N (H 2 ) profile elements distributed parallel to the spine, are shown with red filled circles. Similarly, the error bars are corresponding to the standard deviations in corresponding mean column densities of N (H 2 ) profile elements.
Assuming that the filament follows an idealized cylindrical model, we have extracted the volume density (n(r)) profile correspond to each N (H 2 ) profile by using Plummer-like profile (Nutter et al. 2008;Arzoumanian et al. 2011;Juvela et al. 2012;Palmeirim et al. 2013) The factor, A p = π, is estimated using the equation A p = +∞ −∞ [(1 + u 2 ) p/2 cos i] −1 du (Arzoumanian et al. 2011) and by assuming that the filament is on the plane of the sky implying an inclination angle of i = 0, and p = 2. In equation 7, N (r) is the column density profile as a function of offset distance (r) from the cloud spine. The fitting was performed using IDL mpfit nonlinear least-squares fitting programme (Markwardt 2009) and extracted three parameters: (i) n c the central mean volume density of the spine point, (ii) R f lat the radius within which the N (H 2 ) profile is flat, and (iii) p the power-law index that determines the slope of the power law fall beyond R f lat . Plummer-like profile fit was performed only to the data spanning to ∼5 pc (7.3 arcmin) from the spine, beyond which background column density dominates (horizontal dashed lines in Figures 12(a) & (b)). The fitted values of n c , R f lat , and p are (7.6±1.3)×10 4 cm −3 , 0.24±0.11 pc and 2.1±0.1 for NW part ( Figure  Fig Contours are plotted at CO intensity ranging from −3 K to 10 K with a step of 1 K. Blue dotted line connects the three brightest intensity points along the y-axis denoting a velocity gradient of 13 km s −1 pc −1 . (b) Same as (a) but for a cut along the cloud major axis (cut '2' with a position angle = 65 • , and central offset coordinates are α = 11 h 12 m 04 s , δ = −61 • 18 44 [J2000]). Contours are plotted at CO intensity ranging from −2 K to 26 K with a step of 2 K. Blue dotted line connects the brightest intensity points along the x-axis denoting a velocity gradient of 0.35 km s −1 pc −1 . Unit of the color scale in both panels is K. 12(a)), and (5.1±0.9)×10 4 cm −3 , 0.37±0.18 pc and 2.0±0.1 for SE part (Figure 12(b)). These fit results, as indicated in the plots, are similar.
Volume density map, created using n(r) profiles of all the filament spine elements, is shown in Figure 12(c). Since we considered spine point on the filament as the zeroth element while performing Plummer-like profile fitting for both the NW and SE regions, each spine point will have two volume density values and are averaged in the final map. The minimum and maximum n(H 2 ) values on the filament spine lie between 2.48×10 3 cm −3 and 1.52×10 5 cm −3 . The mean n(H 2 ) values for pixels at the extreme-end of the NW and SE regions are 111±26 cm −3 and 228±74 cm −3 , respectively. This implies that the RCW57A is still embedded in a diffuse background cloud.
We found that more than 95% of the polarization vectors are centered outside of the blue isodensity contour corresponds to volume density n(H 2 ) = 13,000 cm −3 . It implies that our NIR polarimetry fails to trace adequately the highly extincted parts as as shown in Figure 12(c). There-fore, we exclude the pixels with n ≥ 13,000 cm −3 while estimating the mean volume densities in A, B, C, D, and E regions. Means and Poisson errors (standard deviation divided by square root of the number of measurements) of the volume densities for the five regions (A, B, C, D, and E) are estimated and are listed in column 10 of Table  2. More details on producing volume density map from column density can also be found at Smith et al. (2014) and Hoq et al. (2017).
It should be noted here that there exist 12 stars within the high column density regions (within the contour of volume density n(H 2 ) = 13,000 cm −3 or column density N (H 2 ) 1×10 23 cm −2 ). Only one star (star ID = 290; Table 1)  to the cloud. Remaining 11 stars are consistent with either reddened cluster members embedded in the cloud or reddened field stars lying in the outskirts of the cloud but projected on the high density parts of the cloud region.

Magnetic field strength using Chandrasekhar-Fermi method
Using the dispersion in the polarization angles (σ θ(H) ; column 8 of Table 2), velocity dispersion values (σ V LSR ; column 9 of Table 2) and mean volume densities (n(H 2 ); column 10 of Table 2), for the five regions of RCW57A, we estimated the B-field strength using the Chandrasekhar & Fermi (1953) relation: The mass density ρ = n(H 2 ) m H µ H2 , where n(H 2 ) is the hydrogen volume density, m H is the mass of the hydrogen atom, and µ H2 ≈ 2.8 is the mean molecular weight per hydrogen molecule and includes the contribution from helium. The correction factor Q = 0.5 is included based on the studies using synthetic polarization maps generated from numerically simulated clouds (Ostriker et al. 2001;Heitsch et al. 2001) which suggest that for σ θ ≤ 25 • , B-field strength is uncertain by a factor of two. Uncertainties in B-field strength were estimated by propagating the errors in σ θ(H) , σ V LSR and n(H 2 ) values. The estimated B-field strengths and corresponding uncertainties are given in column 11 of Table 2. The uncertainties in B-fields, derived by propagating errors in n(H 2 ), σ V LSR , and σ θ H , are considerably large (see column 11 of Table 2). We also estimate the B-field strengths and corresponding uncertainties without considering errors in σ θ H values. The resultant B-field strengths, given in column 12 of Table 2, for A, B, C, D, and E regions are to be 123±21µG, 89±13µG, 63±15µG, 108±27µG, and 74±8µG, respectively. Mean B-fields strength over five regions is estimated to be 91±8µG.
To understand the importance of B-fields with respect to turbulence, we estimated the magnetic pressure and turbulent pressure using the relations P B = B 2 /8π and P turb = ρσ turb 2 (where σ turb = σ V LSR is given in Table 2), respectively, and the results are given in Table 4. The mean P B /P turb is estimated to be 2.4±0.6. These estimated parameters suggest that the magnetic pressure is more than the turbulent pressure at least by a factor of ∼2 in all five regions, signifying the dominant role of B-fields over turbulence. For the five regions together, the mean magnetic and turbulent pressures are estimated to be (35±7)×10 −11 (dyn cm −2 ) and (15±3)×10 −11 (dyn cm −2 ), respectively. Danziger (1974) estimated the mean electron density, n e , as ∼20 cm −3 , and mean electron temperature, T e , as ∼10000 K for the RCW57A region. For the entire region the thermal pressure P th using the relation P th 2n e kT e (where k is the Boltzmann constant), is estimated to be 6×10 −11 dyn cm −2 , which is smaller than the mean magnetic pressure. Ratio of mean magnetic to thermal pressure, P B /P th , is found to be ∼6. A lower value of P th again suggests dominant magnetic pressure over thermal pressure. Therefore, in comparison to thermal and turbulent pressures, magnetic pressure is playing a crucial role in guiding the expanding ionization fronts or outflowing gas in RCW57A. The implications on the relative importance of B-field pressure in comparison to turbulent and thermal pressures on the formation and evolution of bipolar bubble is discussed in the Section 5.3.

Anisotropic distribution of material in the cloud and the bipolar bubble
Based on the CO and Spitzer images, number of stars with polarization detection, and the column density profiles we infer on the distribution of gas and dust in RCW57A. The quiescent molecular cloud gas traced by CO (red background in Figure 7) is absent at the foot points of, as well as along, the bipolar bubbles. Moreover, the spatial distributions of the CO gas and the bubbles (green background by Spitzer images) appear to be anti coincident with each other. This implies that cloud material has been either eroded or blown out by the expanding I-fronts, outflowing gas, or strong stellar winds from the embedded early type protostars. Therefore, CO and Spitzer maps suggest that bipolar bubble likely blown by the massive stars emanated stellar winds, outflowing gas, and expanding ionization fronts (see Townsley 2009).
While the SE bubble appears to be closed towards the SE, the NW bubble is widely extended towards the NW (Figures 1 and 7) and exhibits large-scale loops, likely due to anisotropic distribution of material towards SE and NW regions. This is further corroborated by the following facts. Although the area together covered by regions A, B and C (located in NW) is similar to that of D (located in SE), the total number of background stars with polarization detection (51) in A, B, and C exceeds those in D (21) (see Table 2, and Figures  7 and 8). Furthermore, the mean N (H 2 ) profiles suggests that anisotropic distribution of gas and dust in RCW57A i.e., SE part (N (H 2 ) ∼1.3×10 22 cm −2 ) has relatively more background material than in the NW (N (H 2 ) ∼0.70×10 22 cm −2 ) as depicted with dashed lines in Figures 12(a) and (b). Similar asymmetric distribution of column density profiles was seen towards other region (e.g., Taurus B211/3 filament; Palmeirim et al. 2013).
An X-ray study (Townsley et al. 2011, see their Figures 3 and 4) revealed the presence of OB stellar association (NGC 3576OB) towards NW of the main embedded star forming region RCW57A (or NGC 3576GHIIR). Furthermore, presence of hard X-ray emission and a pulsar (+PSRJ1112-6103 & PWN) in NGC 3576OB region may suggests a fast occurrence of supernovae event. This supernova together with OB stars of NGC 3576OB stellar association, most probably, evacuated the material in the NW part of the RCW57A. This could be the plausible reason behind the low density material observed in the NW region. The observed large scale loops in the NE, possibly, formed due to the newly formed material as a result of expanding Ifronts driven by the embedded massive protostars of RCW57A.
A high velocity outflow signature discernible close to the eastward from the center, while the same feature is absent to the westward as shown in Figure 11(b). This could attribute to the confined ionized gas flow towards west from center possibly due to the presence of a dense clump (black thin contours in Figure 7). High velocity component is only seen towards the east which implies dearth of dense material there and hence ionized gas expanded with relatively more velocity. de Pree et al. (1999) and Purcell et al. (2009) have also shown that the Hii region expands more freely towards the east and is bounded to the west.

Column density profiles: filament width and B-field support
We have fitted a Gaussian to the mean column density profiles spread over ∼5 pc radii and resultant FWHM, corresponds to the filament characteristic width, is found to be 1.55±0.09 pc. However, a constant filament width ∼0.1 pc has been determined for the central dense sections of filaments of the nearby star forming regions in the Herschel Gould Belt survey by fitting a Gaussian to the inner sections (radii upto 0.3 -0.4 pc) of the filaments by (Arzoumanian et al. 2011). Therefore, we fit a Gaussian to the inner 1 pc width of the clubbed mean column density profiles of NW and SE regions and the resultant σ and FWHM are found to be 0.21±0.08 pc 0.50±0.18 pc, respectively. Within the error, the filament width (0.50 pc; green lines in Figures 12(a) and (b)) of RCW57A is higher than the constant filament width ∼0.1 pc (Arzoumanian et al. 2011). Filament widths larger than 0.1 pc (spanning over 0.1 pc -1 pc with the typical widths of ∼0.2 -0.3 pc) has also been witnessed in a number of studies (eg., Hennemann et al. 2012;Juvela et al. 2012). However, the filament width of RCW57A is not resolved because of the limited resolution (36. 4 corresponds to 0.42 pc as shown with cyan lines) of the column density map.
The mean power-law index p = 2.02±0.05, derived from Plummer-like profile fitting for RCW57A, is consistent with other studies (Arzoumanian et al. 2011;Juvela et al. 2012;Palmeirim et al. 2013;Hoq et al. 2017). This similar value of index (p = 2) is expected if the filament is supported by B-fields (Hennebelle 2003;Tilley & Pudritz 2003) which is in accordance with our observational results supporting active role of Bfields in the formation and evolution of filament as described in Section 5.3 below.

Evolutionary scenario of filament and bipolar bubble
Based on the morphological correlations among (i) dense filamentary cloud structure seen in sub-mm (at 0.87-mm ATLASGAL and 1.1-mm SIMBA dust continuum emission maps), (ii) bipolar bubble seen in mid infrared images, and (iii) the B-field morphology of RCW57A trace by NIR polarimetry, we postulate a scenario as illustrated with four schemas in Figure 13. Each schema describes our observational evidence along with our prediction at a particular evolutionary state of RCW57A.
Schemas A and B: As per the background star polarimetry, the filament (position angle ∼ 60 • ; Figure 7) is oriented perpendicular to the Bfields (the dominant component of θ(H) = ∼ 163 • ; Figure 6(d)). From this observational evidence, we speculate the formation history of filament and the role of B-fields in them. Initially a subcritical molecular cloud is supported by B-fields as shown in schema A, later compresses into a filamentary molecular cloud due to B-fields guided gravitational contraction as per schema B.
Although the role of B-fields in the formation and evolution of filamentary molecular clouds is still a matter of debate, its morphology with respect to the geometry of filaments is well established. The geometry of B-fields in a molecular cloud is mainly governed by the relative dynamical importance of magnetic forces to gravity and turbulence.
If the B-fields are dynamically unimportant compared to the turbulence, the random motions dominate the structural dynamics of the clouds, and the field lines would be dragged along the turbulent eddies (Ballesteros-Paredes et al. 1999). In that case, B-fields would exhibit chaotic or disturbed B-field structures. Based on the facts that the B-field structure is more regular with an hourglass morphology and B-field pressure dominates over turbulent pressure (see Section 4.4), we do not favor weak B-field scenario in RCW57A.
If B-fields are dynamically important then the support to the molecular cloud against gravity Fig. 13.-Schematic diagram illustrates the possible connection between filament, bipolar bubble and B-fields. Initially, subcritical molecular cloud permeated by uniform B-fields (schema A) could have gravitationally compressed along the B-fields into a filament (schema B). Later, as mass accumulated at the center of the filament, gravitationally contracting cloud could have pulled the B-fields along the cloud material towards the central region. Thus, B-fields followed a clear hour-glass morphology when the collapsing cloud reached supercritical state (schema C). Dense and high pressurized Hii region formed due to the formation of massive stars at the cloud center (schema D). Further, propagation of ionization fronts and stellar winds is rather controlled by the an isotropic distribution of gas and dust in the filament. Presence of B-fields in an hour-glass morphology is further introduced an additional anistropic pressure in favor of expanding ionization fronts or outflowing gas so as to form the bipolar bubble (schema D). is rendered predominantly by the B-fields. In this case, B-fields would be aligned, preferentially, perpendicular to the major axis of the cloud (Mouschovias 1978). This is expected as the cloud tends to contract more in the direction parallel to the B-fields than in the direction perpendicular to the it. Recent polarization studies have also demonstrated that B-fields are well ordered near dense filaments and perpendicular to their long axis (Moneti et al. 1984;Pereyra & Magalhães 2004;Alves et al. 2008;Chapman et al. 2011;Sugitani et al. 2011;Palmeirim et al. 2013;Li et al. 2013;Franco & Alves 2015;Cox et al. 2016;Planck Collaboration et al. 2016). Numerical simulations have also shown that filamentary molecular clouds form by gravitational compression guided by B-fields (e.g., Nakamura & Li 2008;Li et al. 2013). Therefore, we believe that the fil-ament in RCW57A could have formed due to the gravitational compression of cloud along the Bfields.
Schema C: Due to the accumulation of more mass at the center, the filament is achieved a supercritical state from its initial B-field supported subcritical state. B-fields are dragged along the gravity-driven material contraction towards the strong gravitational potential at the cloud center, and subsequently they configured into an hour-glass morphology (schema C). Similar hourglass morphology of B-fields were observed towards Orion molecular cloud 1 (Schleuning 1998;Vaillancourt et al. 2008;Ward-Thompson et al. 2017) and Serpens cloud core (Sugitani et al. 2010) signifying ongoing gravitational collapse.
In RCW57A, B-fields in some portions exhibit clear evidence for hour-glass morphology. In the NE part of the filament (covered by two regions A and D; see Figure 8), B-fields are bent by ∼ 90 • as the mean θ(H) for region A is ∼26 • and for D is ∼132 • (see column 7 of Table 2). Therefore, B-fields are distorted in the NE part of the filament, signifying a collapsing filamentary molecular cloud. However, this similar distortion is not clearly apparent in the SW part of the filament, partly because of the light from the background stars is obscured by the dense cloud material. Therefore, the NIR polarimetric observations were unable to detect sufficient numbers of background stars in the SW part to delineate the detailed structure of the B-field there.
Schema D: Gravitational instability caused the filament fragmentation into seven cores (blue diamonds in Figure 7) as evident from dust continuum emission maps (André et al. 2008) as well as dense gas tracers (Purcell et al. 2009). Star formation has began at the central region due to the onset of gravitational collapse during which B-fields were not strong enough to stop the cloud collapse. Ionization and shock fronts, driven by the Hii region, are propagated in to the ambient medium. Interaction between the Hii region and dense cloud material is evident from Figure 7 and as a consequence of this, several IRS sources are formed at the boundary of interaction signifying triggered star formation at the edges of Hii region of RCW 57A.
Our observational results suggest that B-field pressure dominates over thermal pressure (see Sec-tion 4.4), implying active role of B-fields in governing the feedback processes. Since the B-fields observed to be responsible for the formation of filament in RCW57A (Schema B), it is expected that B-fields anchored through the filament may also have significant influence on the expanding I-fronts. In addition to the anisotropic pressure that the I-fronts experience while they expand anisotropically in the filament, they undergo an additional anisotropic pressure in the presence of strong B-fields (Bisnovatyi-Kogan & Silich 1995). As a result I-fronts experience accelerated flows along the B-fields (i.e., low pressure) and hindered flows (i.e., high pressure) in the perpendicular direction to the B-fields (Tomisaka 1992;Gaensler 1998;Pavel & Clemens 2012;van Marle et al. 2015). Therefore, I-fronts expand into greater extent along the hour-glass shaped B-fields and eventually form the bipolar bubble as shown in schema D of Figure 13.
Similar applications for the influence of B-fields on I-fronts have found at other environments. Bubbles around young Hii regions and supernova remnants have been found to be elongated along the Galactic B-field orientation (Tomisaka 1992;Gaensler 1998;Pavel & Clemens 2012). Falceta-Gonçalves & Monteiro (2014) have studied that the strong B-fields can lead to form the bipolar planetary nebula. van Marle et al. (2015) have investigated the influence of B-fields on the expanding circumstellar bubble around a massive star using magneto-hydrodynamical simulations. They found that that the weak B-fields cause circumstellar bubbles become ovoid rather than spherical shape. On the other hand, strong B-fields lead to the formation of a tube-like bubble.
Therefore, according to the proposed evolutionary scenario (Figure 13), initially B-fields in molecular clouds may be important in guiding the gravitational compression of cloud material to form a filament. In the later stage, hour-glass B-fields might have guided the expansion and propagation of I-fronts to form the bipolar bubble.

B-fields favoring the formation of massive cluster at the center of the filament
Based on the presence of 29 NIR excess sources (Persi et al. 1994) and finding of 51 stars earlier than A0 (Maercker et al. 2006), it has been found that RCW57A hosts a massive infrared cluster at the center (within the cyan contour encloses Hii region in Figure 7) of the filament. Furthermore, the Chandra X-ray survey by Townsley et al. (2014, see also Townsley et al. 2011 has provided a first detailed census about the members in the embedded massive young stellar cluster (MYSC; see their Figures 9c and 9d). Their study unraveled many highly-obscured but luminous X-ray sources (one source exhibits photon pile-up) that are likely the cluster's massive stars ionizing its Hii region. This highly-obscured X-ray cluster is situated at the SW edge of the infrared cluster. Formation history of this cluster is poorly understood. Below, we propose that the large scale B-fields traced by NIR polarimetry may also govern crucial role in the formation of this massive cluster.
Initial turbulence in the molecular clouds may decay rapidly. To maintain turbulence in the cloud and hence to govern cloud stability against rapid gravitational collapse, supersonic turbulence should be replenished by means of prostellar outflows. In the model of outflow driven turbulence cluster formation (Li & Nakamura 2006;Nakamura & Li 2007, and references therein), Bfields are dynamically important in governing two processes. First, the outflow energy and momentum are allowed by the B-fields to escape from the central star forming region into the ambient cloud medium. Second, B-fields guide the gravitational infalling material (just outside of the outflow zone) to the central region. Therefore, the competition between protostellar outflow-driven turbulence and gravitational infall regulates the formation of stellar clusters in the presence of dynamically important B-fields. Observational signatures towards the Serpens cloud core (Sugitani et al. 2010) is in accordance with the scenario of B-field regulated cluster formation: (a) B-fields should align with the short axis of the filament, outflows and gravitational infall motions, and (b) outflow injected energy should be more than the dissipated turbulence energy in order to maintain the supersonic turbulence in the cloud. In the case of RCW57A, the presence of a deeply embedded cluster, and alignment of the I-fronts, outflowing gas, and bipolar bubbles with the NIR-traced Bfields, trace the pre-existing conditions in favor of cluster formation according to the protostellar outflow driven turbulence cluster formation in the presence of dynamically strong B-fields.

Summary and Conclusions
Though there exist few studies regarding the formation of bipolar bubbles, none of them has explored the importance of B-fields. Our aim was to understand the morphological correlations among the B-fields, filament, and bipolar bubbles, and their implications for the star formation history in RCW57A. We conducted NIR polarimetric observations in the JHK s -bands using SIRPOL to delineate the B-field structure in RCW57A. We employed various means (NIR and MIR color-color diagrams, polarization efficiency diagrams, and the polarization characteristics) to exclude YSOs with possible intrinsic polarization. Through our analyses, we separated 97 confirmed foreground stars from 178 confirmed background stars having H-band polarimetry. Below we summarize the results of our present work.
• The foreground dust dominated by a single component B-field having a mean θ(H) ∼ 65 • , which is different from the position angle (101 • ) corresponding to the Galactic plane.
• The polarization values of the background stars consistent with the dichroic origin of dust polarization and exhibit relatively higher polarization efficiencies, suggesting efficiently aligned dust grains in RCW57A.
• The polarization angles of the reddened background stars reveal that the B-field in RCW57A is configured into an hour-glass morphology which follows closely the structure of the bipolar bubbles. The dominant component of the B-field (θ(H) ∼ 163 • ) is perpendicular to the filament major axis. The orientations of both outflowing gas from the embedded protostars and the expanding ionization fronts from the Hii region are aligned with the B-fields.
• The B-field strength averaged over five regions across RCW57A based on the Chandrasekhar-Fermi method, is 91±8µG. The B-field pressure is estimated to be more than turbulent and thermal pressures.
• Morphological correlations among B-fields, filament, and bipolar bubbles as well as the dominance of B-field pressure over turbulent and thermal pressures suggest a scenario in which B-fields not only play an important role in formation the filamentary molecular cloud but also in guiding the expansion and propagation of I-fronts and/or outflowing gas to form bipolar bubbles. In this picture, B-fields impart additional anisotropic pressure to expanding I-fronts, from Hii regions, in order to be expanded and propagated into greater extent.
• Protostellar outflow driven turbulence and gravity in the presence of dynamically important B-fields might be responsible for the cluster formation. Therefore, our study, based on the link between B-fields, filament, and bipolar bubble, traces the preexisting conditions where B-fields might also be important in the formation of massive cluster in RCW57A.

Acknowledgments
We thank Prof D. P. Clemens for helpful discussions and providing constructive suggestions in the improvement of the draft. We thank the anonymous referee for his/her careful reading and insightful comments which have improved the contents of the paper. CE, SPL and JWW are thankful to the support from the Ministry of Science and Technology (MoST)   globule prior to the formation of BRC would be ∼163 • (which corresponds to the dominant component of B-fields in the cloud before the onset of star formation, see Figure 6(d) and Section 4.1) and is shown with a broken cyan line in Figure 15. The direction of propagating ionizing radiation from the Hii region with respect to the BRC makes a position angle of 95 • as shown with a green arrow. This implies that initial B-fields around BRC (∼ 163 • ) are oriented nearly orthogonal to the direction (∼ 95 • ) of I-fronts. Thus, the estimated B-field strength (74±8 µG) and its initial orientation around the BRC are similar to those of R5 models of Mackey & Lim (2011a) involving perpendicularly oriented B-fields with moderate strength (50µG). We believe that this BRC might be formed due to a portion of ionizing radiation that is leaked through the low density parts of the molecular cloud, otherwise it would have evaporated if all the radiation from Hii region is impinged on BRC. Therefore, this BRC might experienced the influence of thermal pressure equivalent to that emanated from a single O-type star. Mackey & Lim (2011a) also used only a single O-type star as the ionizing source in their models. According Mackey & Lim (2013), at 400 kyr, the models with moderate B-fields (∼ 50 µG), vertical motions along the B-fields are initiated (see left panel of Figure 1 of Mackey & Lim 2013). Due to this, gas flows in response to photo-evaporation flows along the B-fields result in the formation of dense ridge filament elongated parallel to the B-fields. At 500 kyr a clear cometary globule is formed in which B-fields are aligned with their bright rim similar to our results. Thus, our results are in accordance with RMHD models of Mackey & Lim (2011a with moderate B-field strength (initially oriented perpendicular to the direction of the radiation propagation). Observational evidences similar to our present results includes compressed B-fields along but behind the bright rim SFO 74 (Kusune et al. 2015).

C. Multi-Gaussian fitting on T vs. V LSR spectra
Initial guess values were provided by visual inspection of observed spectra. The number of Gaussian components are determined based on the number of peaks appeared in each spectrum. Since the 'gatorplot.pro' does not provide the errors in the fitted multi-Gaussians, we further used an output of 'gatorplot.pro' as an input to 'mpfitpeak.pro' of Markwardt IDL suite of functions (Markwardt 2009) to ascertain the uncertainties to the fitted measurements. Instead of varying all three parameters of one Gaussian (peak:T peak , center: V LSR and width: σ V LSR ) we fixed two of them and varied one (similarly, in case of three Gaussians, out of 9 parameters three were varied and six were kept constant). Likewise we constrained all the multi-Gaussian components of each spectrum along with their fitted errors using reduced chi-square minimization method. Equal weights are given to all points of an observed spectrum while fitting was performed.

D. Column density map using Herschel data
Column density map has been constructed using the following relation (Kauffmann et al. 2008) where S beam ν is the flux per beam (mJy beam −1 ), N H2 is the column density and K ν (=0.1 ν(GHz) where A pixel (Sr) is the pixel area. All the images were smoothed to match the resolution of 36. 4 corresponds to the SPIRE 500µm image. Equation D1 is fitted to the pixel-wise flux values of four images using 'mpfit.pro of Markwardt IDL suite of functions (Markwardt 2009). Equal weight is given to four input fluxes while fitting is performed.  (Skrutskie et al. 2006) The polarization values given in columns 4-6 are de-biased The uncertainties in the polarization angles given in columns 7-9 are not accounted for overall angular calibration uncertainty of 3 • , which is a systematic uncertainty which would affect all the measured polarization angles Column 13: 'Foreground' means confirmed foreground stars Column 13: 'Background' means confirmed background stars or cluster members Column 13: 'Foreground † ' means probable foreground stars with P (H) > 3% (see section 3.3) Column 13: 'Background † † ' means probable background stars or cluster member with excess polarization (see section 3.3) Column 13: 'Background ‡ ' means probable background star or cluster member with depolarization (see section 3.3) The entries those do not have 2MASS photometric uncertainties are left blank A portion of the table is given here, entire table will be available online Table 2: Central coordinates, widths, number of stars, mean (µ θ(H) ) and dispersion (σ θ(H) ) in polarization angles (using Gaussian fitting), velocity dispersion (σ V LSR ) using 13 CO(1 -0) data, volume density (n(H 2 )) using Herschel data and Plummer-like profile fitting on column density maps, and magnetic field strength estimated using Chandrasekhar-Fermi method for the regions A, B, C, D and E. Note: † Errors in B-fields strengths were estimated by propagating uncertainties in n(H2), σ θ(H) , and σ V LSR . ‡ Errors in B-fields strengths were estimated by propagating uncertainties only in n(H2) and σ V LSR . Table 3: Multiple Gaussian fitted components of 13 CO brightness temperature (T ) versus V LSR spectra of Table 4: Magnetic and turbulent pressures along with their uncertainties in the regions A, B, C, D and E. Region P B P turb P B /P turb 10 −11 (dyn cm −2 ) 10 −11 (dyn cm −2 ) (1) (2) (3) (4) A 60 ± 20 15 ± 5 4 ± 2 B 32 ± 9 20 ± 6 2 ± 1 C 16 ± 8 5 ± 2 3 ± 2 D 46 ± 23 22 ± 11 2 ± 1 E 22 ± 5 11 ± 3 2 ± 1 Note: In this table, B-field pressure values and their uncertainties are derived using B-field strength and corresponding uncertainties given in column 12 of Table 2.