First observations of the magnetic field inside the Pillars of Creation: Results from the BISTRO survey

We present the first high-resolution, submillimeter-wavelength polarimetric observations of -- and thus direct observations of the magnetic field morphology within -- the dense gas of the Pillars of Creation in M16. These 850$\,\mu$m observations, taken as part of the BISTRO (B-Fields in Star-forming Region Observations) Survey using the POL-2 polarimeter on the SCUBA-2 camera on the James Clerk Maxwell Telescope (JCMT), show that the magnetic field runs along the length of the pillars, perpendicular to, and decoupled from, the field in the surrounding photoionized cloud. Using the Chandrasekhar-Fermi method we estimate a plane-of-sky magnetic field strength of $170-320\,\mu$G in the Pillars, consistent with their having been formed through compression of gas with initially weak magnetization. The observed magnetic field strength and morphology suggests that the magnetic field may be slowing the pillars' evolution into cometary globules. We thus hypothesize that the evolution and lifetime of the Pillars may be strongly influenced by the strength of the coupling of their magnetic field to that of their parent photoionized cloud -- i.e. that the Pillars' longevity results from magnetic support


Introduction
One of the most iconic images taken by the Hubble Space Telescope (HST) was of the "Pillars of Creation" in M16 (Hester et al. 1996). These photoionized columns are typical of those found in high-mass star-forming regions throughout the interstellar medium. M16 is a relatively local (1.8 ± 0.1 kpc; Dufton et al. 2006), well-resolved site of active ongoing star formation (Oliveira 2008), typical of regions forming highmass (>8 M e ) stars (Zinnecker & Yorke 2007). We present the first detailed measurements of the magnetic field (hereafter B-field) in the densest parts of the Pillars, taken as part of the B-Fields in Star-forming Region Observations (BISTRO) survey (Ward- ) on the James Clerk Maxwell Telescope (JCMT) using the Submillimeter Common-User Bolometer Array 2 (SCUBA-2) camera and its polarimeter POL-2.
Young massive stars produce sufficient high-energy photons to ionize a volume of their parent molecular cloud, thereby driving a shock into the cloud (Strömgren 1939;Zinnecker & Yorke 2007). These photoionized regions indicate ongoing high-mass star formation. Complex structures can form in the dense gas at the shock interfaces (Spitzer 1954)-particularly, dense neutral columns are frequently seen protruding into photoionized regions, most famously in M16. The formation and evolution of these pillars remain disputed (White et al. 1999;Williams et al. 2001;Ryutov et al. 2005; hereafter Wh99; Wi01; R05 respectively), with the role of the B-field neither observationally nor theoretically well constrained (Williams 2007;hereafter Wi07). Near-infrared extinction observations of M16 suggest a difference in B-field direction between the Pillars and the surrounding photoionized cloud (Sugitani et al. 2007), but cannot probe the dense gas of the Pillars themselves.
The heads of the Pillars are dense star-forming molecular condensations (Wh99) interacting with the shock front associated with the young (∼1.3 Myr; Bonatto et al. 2006) high-mass cluster NGC 6611 (Hillenbrand et al. 1993). Whether these condensations predate, or were formed by, the shock interaction is uncertain (Wh99; Wi01). The heads are being destroyed by the interaction with NGC 6611, with a lifetime of 3×10 6 year (McLeod et al. 2015), and are thus likely to be considerably longer-lived than the lower-density pillars, which have an estimated lifetime of a few ×10 5 year (Wi01), suggesting that they will become disconnected cometary globules (Bertoldi & McKee 1990) unless another mechanism, such as a B-field, is at work.
We observed the Pillars of Creation in 850 μm polarized light with the POL-2 polarimeter (Friberg et al. 2016) on the SCUBA-2 camera (Holland et al. 2013), giving a map of the B-field in the dense gas of photoionized pillars unprecedented in sensitivity, area, and resolution. We observed Pillars I, II, and III (Hester et al. 1996) at high signal-to-noise ratio (S/N), and Pillar IV, the Spire, and SFO30 (not shown) at lower S/N.

Observations
The Eagle Nebula was observed in 20 separate 40-minute exposures between 2017 June 6 and 2017 July 27, with a total integration time of 14 hr. The observations were taken in JCMT Band 2 weather, with atmospheric optical depth at 225 GHz, τ 225 , of 0.05<τ 225 <0.08. The BISTRO survey's observing strategy is described by Ward-Thompson et al. (2017).
The 850 μm POL-2 data were reduced using the pol map 2 routine, 12 recently added to SMURF (Berry et al. 2005;Chapin et al. 2013). The reduction process is described in detail by Kwon et al. (2018). The output Stokes Q, U, and I maps are gridded to 4″ pixels and are calibrated in mJy beam −1 . The output vectors are debiased using the mean of their Q and U variances to remove statistical biasing in regions of low signal-to-noise.
Our final map has a FWHM resolution of 14 1 (0.12 pc; ∼25,000 au), a diameter of 12′, and an rms noise level of 0.9 mJy beam −1 in Stokes Q and U intensity on 14 1 pixels. Figure 1 shows the observed B-field morphology in the Pillars. We detect Pillars I, II, and the material between their bases (the "Ridge") in polarized light, and marginally detect Pillar III. The B-field clearly runs along the length of the Pillars, apparently turning at the tips of the Pillars (best seen in the head of Pillar I). "Pillar I" has two separate components: Pillar Ia (northwest), located further along the line of sight than II and III, behind the source of ionizing photons; and Pillar Ib (southeast), approximately equidistant with II and III (Pound 1998;McLeod et al. 2015). The apparent change in field direction seen between Pillars Ia and Ib represents different field directions in the two Pillars.

Results
The B-field geometry in the Pillars is significantly different to that in the surrounding photoionized region, as measured using near-infrared extinction polarimetry (Sugitani et al. 2007), as shown in Figure 2. The near-infrared vectors vary smoothly across the photoionized region, producing a singly peaked distribution (at ∼90°east of north). The B-field in the dense gas shows more complex behavior, with field lines running roughly parallel to the length of the Pillars. The B-field distribution in the dense gas is bimodal, peaking at ∼70°(head of Ia, Ib, base of IV, Ridge) and ∼140°(length of Ia, II, IV), compared to mean pillar directions of 134°±17°in I, 132°±12°in II, 144°±16°in IV, and 48°±19°in the Ridge. The B-field vectors observed in Pillar II-upon which our subsequent analysis focuses-are shown in detail in Figure 3. The near-infrared polarization vectors observed by Sugitani et al. (2007) in the vicinity of Pillar II are shown alongside.

B-field Strength
We estimated the plane-of-sky B-field strength in Pillar II-the most well-defined Pillar, with velocity-coherent structure and a linear plane-of-sky morphology-using the Chandrasekhar-Fermi (CF) method (Chandrasekhar & Fermi 1953).
The CF method provides an estimate of the plane-of-sky B-field strength by assuming that the variation in B-field around the mean field direction represents distortion of the B-field lines by non-thermal motions in the gas. The plane-of-sky field strength (B pos ) is given by where ρ is the gas density, σ v is the non-thermal gas velocity dispersion, σ θ is the standard deviation in polarization angle about the mean field direction, and Q is a factor of order unity that accounts for variation in the field on scales smaller than the beam. We take Q=0.5 throughout (Ostriker et al. 2001;Crutcher et al. 2004). The second form of the expression takes number density of molecular hydrogen (n(H 2 )) to be in cm −3 , FWHM non-thermal gas velocity dispersion ( to be in km s −1 , and σ θ to be in degrees (Crutcher et al. 2004).
We gridded the data to 14 1 (statistically independent) pixels, and selected pixels with S/N in total intensity I of I/δI>10 associated with Pillar II. Of these, 16 have S/N in polarization fraction P of P/δP>2, and 11 have P/δP>3. In order to mitigate against small sample size effects potentially introduced by using only the P/δP>3 sample, we found the weighted standard deviations of both samples. The P/δP>3 sample has a weighted dispersion in angle of σ θ =14°. 4, while the P/δP>2 sample has a very similar σ θ =14°. 1. We thus adopt σ θ ∼14°. 4 as being a representative value. We assume that all dispersion in the position angles of the vectors associated with the Pillar represents dispersion about a uniform mean field direction. As the measured angular dispersion is greater than the uncertainty on angle in our vectors, it is not necessary to correct the angular dispersion for measurement uncertainty ). The P/δP values of our data for 14 1 pixels are shown in Figure 4.
We took the gas density in the Pillar to be n(H 2 )=5×10 4 cm −3 (R05), and the FWHM gas velocity dispersion to be in the range Δv=1.2-2.2 km s −1 , as measured by Wh99 in various dense gas tracers. These linewidths are highly supersonic (Wh99), and so the correction for the thermal component is negligible.
We thus estimated a plane-of-sky B-field strength of ∼170-320 μG in Pillar II. This value is intermediate between the B-field strengths of ∼10 μG observed in relatively unperturbed gas in low-mass star-forming regions (Crutcher 2012), and of ∼10 3 μG observed in massive, gravitationally unstable structures in high-mass star formation sites (e.g., Curran & Chrysostomou 2007;Hildebrand et al. 2009;Pattle et al. 2017).

Discussion
Simulations of photoionized regions suggest that B-field orientation is largely unchanged by the free passage of a planeparallel shock front (Henney et al. 2009). Hence, we assume that the B-field in the photoionized region is representative of the B-field direction in the unshocked gas-approximately parallel to the shock front. For a weak initial B-field, field lines are predicted to become aligned parallel to the Pillar's length in the Pillar itself, while remaining approximately perpendicular in the surrounding photoionized region (Wi07; Mackey & Lim 2011;hereafter ML11). This prediction results from otherwise quite different scenarios of magnetized pillar formation.
Wi07 finds that, in two dimensions, when a shock propagates into a dense medium (10 4 cm −3 ) in which a denser core (10 5 cm −3 ) is embedded, a pillar forms behind the core, and the weak, plane-parallel B-field in the dense medium is compressed. Thus, the B-field strength is enhanced by pillar formation, with the field "bowing" into the material behind the pillar. The pillar has a density of a few ×10 4 cm −3 , while the surrounding ionized material has a density ∼10 2 cm −3 . (The pillar head has higher density.) Arthur et al. (2011) found similar behavior in three-dimensional simulations of expanding H II regions, although with lower resolution.
ML11 find that when a shock impinges on a set of approximately co-linear dense globules embedded in a much lower-density medium (200 cm −3 ; c.f. Mackey & Lim 2010) threaded by a weak, plane-parallel B-field, a pillar-like feature forms behind the globules due to radiation-driven implosion and the rocket effect (Oort & Spitzer 1955). These effects orient the B-field along the length of the forming pillar on timescales ∼100 kyr.
Wi07 and M11 agree that a strong plane-parallel initial B-field should deviate significantly from its initial orientation only in the pillar head (see also Henney et al. 2009). Our results do not match this scenario, strongly suggesting that the B-field in M16 was dynamically unimportant in the formation of the Pillars.  (Hester et al. 1996). Vectors are gridded to 4″ (note oversampling), and have polarized intensity S/N PI/δPI>2. Polarization angles are rotated by 90°to show B-field direction. Vector length scale is arbitrary. Black lines delineate the Pillars. Beam size is shown in lower right-hand corner.
ML11 predict a B-field strength in the material around the Pillars of < 50 μG, but do not quantitatively predict the B-field strength inside the Pillars. Our plane-of-sky B-field is in the ML11 "strong-field" regime, which they exclude for M16. It is not clear how the gas compression necessary to increase the B-field could occur in this model. Henney et al. (2009) predicted volume-averaged B-field strengths to remain approximately constant with time within pillars formed behind individual globules. ML11 also show a B-field which, while broadly orientated parallel to the Pillar's length, shows considerable disorder, whereas our observations show an ordered (albeit not well-resolved) B-field along the length of the Pillars.
The Pillars are anchored to a larger cloud (Hester et al. 1996), similar to the Wi07 scenario. Moreover, the Wi07 simulations show the B-field compression necessary to significantly strengthen an initially dynamically unimportant Figure 2. Distribution of the B-field vectors in the dense gas (blue: this work; 850 μm dust polarization, 14″1 pixels, P/δP>3, I/δI>10, I>50 mJy beam −1 ) and in the photoionized region (red: H-band extinction polarization; Sugitani et al. 2007). Gray lines and shaded areas show the approximate orientations of Pillars I, II, and IV and the Ridge, with the range derived from the Pillars' plane-of-sky aspect ratios. Note how the red histogram peaks around ∼90°and the blue histogram peaks either side (roughly parallel to the Pillars and the Ridge, respectively).  Sugitani et al. (2007); excerpt from their Figure 6 (© PASJ, reproduced with permission). 850 μm vectors (this work) have P/δP>2 and I/δI>10. The HST composite is the same as in Figure 1. The B-field runs roughly parallel to the Pillar's axis. No polarization is detected at the Pillar's tip-this depolarization is consistent with a horseshoe-shaped B-field morphology on scales smaller than the beam. field, albeit qualitatively and two-dimensionally. We thus consider the Wi07 scenario to be broadly more consistent with our observations, and so illustrate it in Figure 5. However, both mechanisms could be involved in creating the observed B-field, and neither model quantitatively predicts B-field strength inside the pillars. Detailed, three-dimensional quantitative modeling of the B-field inside photoionized columns is needed to fully distinguish between these mechanisms.

Pressure Balance
Magnetic pressure is given by P B =B 2 /8π. Our measured plane-of-sky B-field, 170-320 μG, implies P B /k B ∼(0.9-3.0)× 10 7 K cm −3 . R05 gave an ablation pressure on the heads of the pillars of 1.6×10 8 K cm −3 , an order of magnitude higher than our inferred P B . This suggests that the B-field cannot support the Pillars against longitudinal erosion by the shock front, unless the field is compressed in the Pillar heads (which are not resolved by our observations).
The effective gas pressure within the Pillars is = P nk T g,int B eff , where T eff is the effective gas temperature and n is number density of particles. Taking n≈n(H 2 )=5×10 4 cm −3 and T=20 K (Wh99), P g,int /k B =1.0×10 6 K cm −3 , an order of magnitude lower than our P B . However, Wh99 and Wi01 argued that nonthermal gas motions create an effectively hydrostatic pressure within the Pillars. The Wh99 FWHM gas velocity dispersion range 10 7 K cm −3 , very comparable to our inferred P B . (This result follows naturally from the assumptions of the CF analysis.) Hester et al. (1996) argue that atomic hydrogen number density n(H)∼29 cm −3 in the M16 photoionized region. Simulations take n(H)∼10 2 cm −3 (assumed by Wi01 and ML11; predicted by Wi07). Using » = ( ) n n 2 H 58 cm −3 (assuming » ( ) n n H e and that the number fraction of helium atoms is small), and taking T=8000 K (Hester et al. 1996;García-Rojas et al. 2006), this implies an external gas pressure~Ṕ k 4.6 10 g,ext B 5 K cm −3 on the Pillars. Using » = ( ) n n 2 H 400 cm −3 ,P k g,ext B 3.2 10 6 K cm −3 , still an order of magnitude lower than our P B and P g,int values. Higgs et al. (1979) found a non-thermal velocity dispersion in the photoionized gas of M16 of σ v =11.5 km s −1 . M e pc −1 , assuming cylindrical symmetry (taking μ=2.8 and n=n(H 2 )=5×10 4 cm −3 ). If non-thermal gas motions within the Pillars create hydrostatic pressure, the critical line mass (Stodólkiewicz 1963;Ostriker 1964 For c eff ≈0.51-0.93 km s −1 , (M/L) crit ≈120-400 M e pc −1 , comparable to the observed M/L. Thus, there may be some concentration of mass toward its axis, somewhat lowering P g,int at the H II region boundary. However, the B-field will provide significant support against radial gravitational collapse, with the observed B-field geometry resisting radial motion of material. We note that these estimates are accurate only to order of magnitude. Our results broadly suggest that the Pillar walls are in approximate pressure equilibrium, with P B and P g,int supporting against P g,ext , and also that, contrary to common assumptions, the Pillar's self-gravity is non-negligible. Both P g,int and P g,ext require a non-thermal component in order to be comparable to our inferred P B . Other sources of external pressure could include ram pressure due to flow of material across the ionization front into the Pillar (e.g., Henney et al. 2009).

The Alfvén Velocity
Our favored scenario requires (a) the flux-frozen (infinite conductivity) approximation (Alfvén 1942; Crutcher 2012) to The ionization front is slowed by the over-density. The flux-frozen B-field "bows" into the forming pillar. (c) The compressed B-field supports the pillar against radial collapse, but cannot support against longitudinal erosion by the shock interaction. Dark blue represents molecular gas; light blue represents ionized material; black line indicates the shock front. Gray dashed lines indicate local B-field direction. Red arrows represent photon flux/ablation pressure, black arrows represent magnetic and internal gas pressure, and green arrows represent confining gas pressure, possibly supplemented by ram pressure.
hold (neutral and ionized material are collisionally coupled; flow across field lines is forbidden), and (b) the pillars to form faster than the compressed B-field can relax to a lower-energy configuration, i.e., the photoionized region must expand faster than the Alfvén velocity (v A ; Alfvén 1942 1.5 km s −1 . Because µ B n in flux-frozen plasma, this value should apply throughout the Pillars' lifetimes, suggesting that they could have formed too quickly for the B-field to react, allowing the observed highly pinched geometry to form (see Figure 5).
The (flux-frozen) B-field geometry should allow longitudinal motion of material along the Pillars, but strongly resist motion across the Pillars that would lead to radial collapse. This suggests that the predicted evolution of the Pillars into disconnected cometary globules (Bertoldi & McKee 1990) may be considerably slowed by the effects of the B-field geometry.

Summary
We have observed the dense gas of the Pillars of Creation in M16 in 850 μm polarized light using the POL-2 polarimeter on the JCMT. We find that the B-field in the Pillars is ordered, running along the length of the Pillars, with a plane-of-sky field strength of ∼170-320 μG, estimated using the CF method. The observed morphology is consistent with the field being dynamically negligible in the Pillars' formation. However, the current B-field strength suggests that magnetic pressure provides significant support against both gravitational and pressure-driven radial collapse of the Pillars, and may be slowing the Pillars' evolution into cometary globules. We hypothesize that the persistence of such photoionized columns as objects connected to their parent molecular cloud may be related to the geometry of their B-fields, and specifically to the relative orientation of the Bfields in the Pillars and their surrounding photoionized regions. The BISTRO project is currently surveying B-fields in the dense gas of many nearby high-mass star-forming regions, thus allowing further testing of this hypothesis in the immediate future.