Studying magnetic fields and dust in M17 using polarized thermal dust emission observed by SOFIA/HAWC+

We report the highest spatial resolution measurement of magnetic fields in M17 using thermal dust polarization taken by SOFIA/HAWC+ centered at 154 $\mu$m wavelength. Using the Davis-Chandrasekhar-Fermi method, we found the presence of strong magnetic fields of $980 \pm 230\;\mu$G and $1665 \pm 885\;\mu$G in lower-density (M17-N) and higher-density (M17-S) regions, respectively. The magnetic field morphology in M17-N possibly mimics the fields in gravitational collapse molecular cores while in M17-S the fields run perpendicular to the matter structure and display a pillar and an asymmetric hourglass shape. The mean values of the magnetic field strength are used to determine the Alfv\'enic Mach numbers ($\mathcal{M_A}$) of M17-N and M17-S which turn out to be sub-Alfv\'enic, or magnetic fields dominate turbulence. We calculate the mass-to-flux ratio, $\lambda$, and obtain $\lambda=0.07$ for M17-N and $0.28$ for M17-S. The sub-critical values of $\lambda$ are in agreement with the lack of massive stars formed in M17. To study dust physics, we analyze the relationship between the dust polarization fraction, $p$, and the thermal emission intensity, $I$, gas column density, $N({\rm H_2})$, and dust temperature, $T_{\rm d}$. The polarization fraction decreases with intensity as $I^{-\alpha}$ with $\alpha = 0.51$. The polarization fraction also decreases with increasing $N(\rm H_{2})$, which can be explained by the decrease of grain alignment by radiative torques (RATs) toward denser regions with a weaker radiation field and/or tangling of magnetic fields. The polarization fraction tends to increase with $T_{\rm d}$ first and then decreases when $T_ {\rm d}>50$ K. The latter feature seen in the M17-N, where the gas density changes slowly with $T_{d}$, is consistent with the RAT disruption effect.

Star formation is a complex process that involves selfgravity, turbulence, magnetic fields, and stellar feedback. Understanding the exact role of magnetic fields in the evolution of molecular clouds (MCs) and star formation process is a challenge of modern astrophysics. In the past decades, there are emerging evidences suggesting the importance of magnetic fields (B-fields) in the evolution of MCs and star formation (McKee & Ostriker 2007;Derek Ward-Thompson & McKee 2020). There are two types of B-field models in which B-fields play contrasting roles. First, the strong B-field models support a paradigm of magnetic pressure acting against gravitational collapse of MCs. This magnetic-driven support is suppressed when the ratio of the core mass to the magnetic flux exceeds a critical value that turns the cloud into a state of gravitational collapse to form a new star (Nakano & Nakamura 1978). Second, in the weak-field models, B-fields are sufficiently weak in MCs and dominated by turbulence. Star formation takes place in filaments that are likely produced by intersection of turbulent supersonic flows (e.g., Padoan & Nordlund 1999;Elmegreen 2000;Mac Low & Klessen 2004;Crutcher 2012). It is, thus, necessary to carry out observations of B-fields in specific MCs (e.g., M17 in this work) to investigate their effects on star formation and the subsequent evolution of the entire clouds to testify theoretical predictions (Seifried & Walch 2015;Federrath 2016;Li & Klein 2019;Derek Ward-Thompson & McKee 2020).
Dust polarization induced by aligned grains is widely used to map B-fields (see e.g., Crutcher 2012). This method is based on the fact that the polarization direction of thermal dust emission is perpendicular to the B-fields. The strength of B-fields can then be estimated using the Davis-Chandrasekhar-Fermi (DCF) method (Davis 1951;Chandrasekhar & Fermi 1953). The DCF method measures indirectly the strength based on the fluctuations of B-fields that are encoded in the dispersion of dust polarization directions. Although we have a full-sky map of dust polarization provided by Planck at 353 GHz (see e.g., Planck Collaboration: et al. 2015) and many studies of B-fields at MC scales (see e.g., filament structures and star-forming cores ;Dotson 1996;Houde et al. 2002;Pellegrini et al. 2007;Chen et al. 2012;Pattle et al. 2017Pattle et al. , 2018Wurster & Li 2018;Chuss et al. 2019;Hennebelle & Inutsuka 2019;Sugitani et al. 2019), dust polarization data at higher resolutions are still lacking. Recently, we start to have highresolution polarisation observations made with the Atacama Large Millimeter/sub-Millimeter Array (Liu et al. 2020;Beuther et al. 2020;San 2021).
Dust polarization also allows us to get insights into fundamental properties of dust grains such as grain shapes, size distribution, and alignment. Grain alignment is a longstanding problem of astrophysics, and the leading theory for grain alignment is the Radiative Torque Alignment theory (RATA; Draine & Weingartner 1997;Lazarian & Hoang 2007;Hoang & Lazarian 2016) (see reviews e.g., Andersson et al. 2015;Lazarian et al. 2015). Note that radiative torques (RATs) arising from the interaction of an anisotropic radiation field with an irregular grain were first suggested by Dolginov & Mytrophanov (1976) and later numerically demonstrated by Draine & Weingartner (1996), and analytically modelled by Lazarian & Hoang (2007). The grain alignment efficiency by RATs is found to increase with increasing radiation field or dust temperature (see Hoang et al. 2021), which results in the increase of the polarization fraction with the dust temperature (Lee et al. 2020). Furthermore, Hoang et al. (2019) realized that large grains will be disrupted and depleted around a very strong radiation source via a mechanism, the so-called Radiative Torque Disruption (RATD; see Hoang 2020 for a review). The basic idea of the RATD mechanism is that an intense radiation field can spin up dust grains to extremely fast rotation, such that the centrifugal stress can exceed the tensile strength of the grain material and disrupts the dust grain into smaller fragments. The RATD effect is found to decrease the polarization fraction predicted by the RATA theory (Lee et al. 2020). The combination of RATA and RATD was demonstrated to successfully reproduce the observed dust polarization data in various regions of strong radiation fields such as Oph-A (Tram et al. 2021a), 30 Doradus . Therefore, dust polarization observations toward strong radiation sources like M17 are crucial to test grain alignment and disruption by RATs.
M17 is a well-known star-forming region (Povich et al. 2009;Lim et al. 2020) which is located in the Omega Nebula or Swan Nebula (also known as Horseshoe Nebula) in the constellation of Sagittarius at a distance of 1.98 kpc (Xu et al. 2011). Figure 1 is a RGB image of the region using Spitzer 1 Galactic Legacy Infrared Midplane Survey Extraordinaire (GLIMPSE) mosaic data. Besides the study of B-fields and dust physics, since M17 is the closest giant H II region to Earth it is thus an excellent laboratory for the investigation of stellar feedback from a nearby massive star cluster and a photodissociation region (PDR) as well as infrared sources in the regions such as UC 1, IRS 5, CEN 92, Anon 1, and Anon 3 (Lim et al. 2020). In the present work, we use M17 N M17 S N E Figure 1. An RGB image of M17 observed by Spitzer (R = 8 µm, G = 4.5 µm, and B = 3.6 µm). The white rectangle is the observed region of SOFIA/HAWC+ used in this work. data taken by the High-resolution Airborne Wideband Camera Plus (HAWC+; Harper et al. 2018) accommodated on Stratospheric Observatory for Infrared Astronomy (SOFIA; Temi et al. 2018) to serve our purposes.
The structure of the paper is organized as follows. The SOFIA/HAWC+ observations of M17 are presented in Section 2. In Section 3, we describe the use of the DCF method to estimate the strengths of B-fields. We then present the obtained results, including the B-field morphologies and strengths, and discuss the implications of dust polarization fraction for grain alignment and disruption in Section 4. A summary of our main findings is presented in Section 5.  The thermal dust polarization data of M17 were obtained by SOFIA/HAWC+ instrument centered at 154 µm with a designed beam size of 13.6 . The observed region is shown in Figure 1. The inferred maps have an original pixel-size of 6.9 (Harper et al. 2018). Then, Nyquist samplings were applied during the data reduction processes to have the final Stokes parameter maps with a re-sampled pixel-size of ∼ 3.4 . The Stokes I, Q, and U maps are shown in Figure 2.
The biased polarization fraction, p bias , is calculated as follows (Gordon et al. 2018): where I p = Q 2 + U 2 is the polarization intensity. The associated error on the biased polarization fraction is where σ Ip and σ I are the uncertainties on I p and I, respectively. The error propagation on the non-linear function I p (Q, U ) is expanded as Assuming Q and U are uncorrelated, namely the covariance term σ QU = 0, we obtain the error for the polarization intensity where σ Q and σ U are the errors on Stokes Q and U , respectively. The debiased polarization fraction, p, is calculated following the approach by Wardle & Kronberg (1974) The polarization angle, θ, is defined as and the error on the polarization angle is calculated as follows A detailed exploration of the raw data is presented in Appendix A. In the following, we apply two quality cuts to the data, requiring a high signal-to-noise ratio (SNR) on the measurements of (1) the intensities SNR I > 236 and (2) the polarization fraction SNR p > 3. In what follows, we call this 'master cut'. This quality cut automatically satisfies the third criterion recommended by the SOFIA collaboration for high-quality scientific data, namely p < 50% (Gordon et al. 2018). The choice of the cut on SNR I is to have the measurement uncertainties on the polarization fraction, σ p , better than 0.6%. The approximate correlation between these two quantities is (Gordon et al. 2018). After applying the master cut, the remaining pixels are 5139, corresponding to 23% of an original map of the size 138×162 (22356 pixels). The master cut excludes the low level emission and significantly improves the quality of the data. Its performance is illustrated in Figure A.1 and  Figure 3 shows the inferred B-field orientation map of M17. The lengths of the line segments are proportional to the polarization fraction, and their orientations are obtained by rotating 90 • from the polarization angles to trace the B-fields. These segments are called 'half-vectors' since the directions of the vectors are not known. The map gives the first impression of the general B-field morphology and the polarization properties of the 154 µm emission of the region.

DATA ANALYSIS
The Davis-Chandrasekhar-Fermi method is one of the most well-known techniques for the estimation of the B-field strength. There is some controversy around the method (e.g., Pattle & Fissel 2019;Liu et al. 2021), however the DCF method still provides a mean for estimating the magnetic field strength from dust polarization measurements. The method is based on the fact that turbulent motions induce a turbulent component on top of the mean, large-scale magnetic field and the assumption that the turbulent kinetic energy is balanced by the induced magnetic energy. The strength of B-fields in the plane-of-sky, B POS , is calculated using the following equation (Crutcher 2004) where Q c is a factor of the order of unity to correct for the light-of-sight and beam-integration effects, ρ is the gas density in g cm −3 , ∆V = σ ν √ 8 ln 2 = 2.355σ ν is the FWHM of the one dimensional non-thermal velocity dispersion, σ ν , in km s −1 , n(H 2 ) is the volume density of the molecular hydrogen in cm −3 , and σ θ is the polarization angle dispersion in degrees. In general, the magnetic field strength in MCs is in the range of microto milligauss. Figure 3 indicates two distinct regions encompassed by the two ellipses in terms of physical conditions, including density, temperature, and dynamics, as can be seen in Figures 5, 7, and 11 in the next sections. The ellipse in the north has a center at RA ∼ 18 h 20 m 31 s .6 and DEC ∼ −16 • 08 04 .3, major × minor axes of 137 × 69 , and a position angle of 14 • . The ellipse in the south has a center at RA ∼ 18 h 20 m 23 s .2 and DEC ∼ −16 • 12 24 , major × minor axes of 200 × 143 , and a position angle of 0 • . Therefore, we divide the map into two corresponding northern and southern regions coined M17-N and M17-S for further analyses. Subsequently, we are going to identify the input parameters for the DCF method to calculate the B-field strengths of the two regions.

Polarization angle dispersion: Structure function
Structure function method was initially proposed by Falceta-Gonçalves et al. (2008) and Hildebrand et al. (2009) for evaluating the polarization angle dispersion, σ θ , in the plane-of-sky. Fundamentally, this method calculates a two-point correlation function for pairs of the polarization angles as follows, where ... denotes an average; N ( ) is the number of pairs of pixels having a displacement between the two pixels | | = ; x and x + are the location vectors of the two considered pixels with corresponding polarization angles θ(x) and θ(x+ ), respectively. We note that only the pairs with ∆θ( ) = |θ(x) − θ(x + )| < 90 • are retained for further analyses because the directions of the half-vectors are unknown. Assuming that the B-fields have two independent components (a large-scale structure field, B 0 , and a turbulent one, δB), the structure function can be written as below for a small displacement (Hildebrand et al. 2009;Crutcher 2012) where b is the contribution of turbulence, m of the large-scale structure component; and σ M ( ) is the measurement uncertainties calculated for each pair of polarization angles. In practice, σ M ( ) is taken into account while calculating ∆θ ( ) 1/2 in Equation 9. The ratio of the turbulent to large-scale structure components is expressed as (Hildebrand et al. 2009) From the definition of polarization angle dispersion, σ θ , in Equation 9 combining with Equation 10 neglecting the contribution of the large-scale structure component, m , we have σ θ = b/ √ 2. We calculate the structure functions for M17-N and M17-S separately. In general, they display similar rising tendencies at small displacements and reach the value of random field of 52 • at large displacements (see Figure  4). The extent of M17-N is smaller than that of M17-S, hence the maximum displacement of M17-N is smaller. We then fit the calculated structure functions with the -dependence of the polarization angle dispersion using Equation 10. The fits were carried out for > 13.6 to avoid the beam effects. The upper -limits were chosen to achieve the best fits (black curves in Figure 4). Table  1  which means σ θ = 3.5 • ± 0.2 • and 6.2 • ± 0.5 • for M17-N and M17-S, respectively. The uncertainties on b are from the fits. We note these values are smaller than that found by Hildebrand et al. (2009), σ θ ∼ 10 • , using data from Caltech Submillimeter Observatory at 350 µm with a spatial resolution of ∼ 20 for M17. Table 1 also lists the ratios of the turbulent to large-scale structure B-field components calculated using Equation 11 which shows that the turbulent fields are much smaller than the large-scale one in M17.

Velocity dispersion
We estimate the velocity dispersion of M17 using the publicly available archival data 13 CO (J = 1 → 0) taken by Nobeyama 45-m telescope ). The measurements were made at the cen- tral frequency ν central ∼ 110.201 GHz with a spectral resolution of 0.1 km s −1 . The resulting map has a pixel size of 7.5 and a beam size of 14.9 which was then cropped to match with the M17 observed region by SOFIA/HAWC+. Figure 5 presents the velocityintegrated intensity (left) and the mean velocity (right) maps. The integration was done over the entire velocity range from −20 to 60 km s −1 in the local standard of rest (LSR) frame. Figure 6 displays the integrated spectra averaged over the number of pixels for M17-N (left) and M17-S (right). They can be welldescribed by two-Gaussian fits giving the mean positions and corresponding standard deviations, σ ν; 13 CO , of the peaks. The spectra of M17-N and M17-S show the main peaks at 22.7 ± 1.6 and 19.9 ± 2.3 km s −1 , respectively. These main peaks are much higher than the second ones. Therefore, we only use the velocity dispersion of these two peaks for the calculation of B-field strengths in M17-N and M17-S. Table 1 2020), although their results are for wider areas of M17. The statistical uncertainties of the velocity dispersion from the two-Gaussian fits are smaller than the spectral resolution of Nobeyama, therefore, we take the uncertainties equal to the spectral resolution of 0.1 km s −1 . We note that M17-N is affected by the outflows and shocks from G015.128 2 (Lim et al. 2020) which is expected to be more turbulent than M17-S.
The obtained standard deviations are then converted to the non-thermal velocity dispersion using mass equal to 29 amu, k B is the Boltzmann constant, and T is the gas temperature. It turns out that the thermal contribution to the velocity dispersion is negligible if we adopt the average gas temperature of 20 K according to Nguyen-Luong et al. (2020). Therefore, we take the non-thermal velocity dispersion for the M17-N and M17-S regions equal to the standard deviations of their corresponding main peaks listed in Table 1.

Column and volume densities
The column density, N (H 2 ), has been derived based on a graybody (i.e., modified blackbody) fit with the filter weighted opacity, κ ν , and blackbody radiation, B ν (T d ). As following the techniques described in Lim et al. (2016), Herschel 160, 250, and 350 µm data of M17 have been convolved to the beam size of Herschel 500 µm images (∼ 36 ) before executing the pixel-by-pixel graybody fitting via all four band data. The template T d of angular resolution ∼ 36 is re-gridded to match to the pixel geometry of SCUBA2 850 µm map (angular resolution ∼ 14 and pixel size=4 ). We then repeated the graybody fit by utilizing SCUBA2 850 µm map (Reid & Figure 11. The blue ellipses are the same as in Figure 3. In M17-S, the relation of T d and N (H2) is complex. T d is positively correlated to N (H2) up to ∼ 43 K, but turns to negatively otherwise. Table 1. Summary of the obtained results. b is the rms contribution of the turbulent component and m of the large-scale B-fields to the polarization angle dispersion, σ θ . δB/B0 is the ratio of the turbulent to the large-scale fields. Detailed description of the parameters can be found in the text.
(9.1 ± 4.0) × 10 21 (6.2 ± 6.5) × 10 22 Volume density, n(H2) [cm −3 ] (9.3 ± 4.1) × 10 3 (4.2 ± 4.4) × 10 4 Wilson 2006) of M17 with the re-gridded T d map. Figure 7 shows the resulting N (H 2 ) map superimposed by the dust temperature ( Figure 11). Using this map, we calculate the total column densities, N total,H2 , the average column densities, N (H 2 ) , as well as their associated uncertainties, σ N (H2) , for the M17-N and M17-S regions. The results are presented in Table 1. We then follow the approach described in Lee et al. (2012); Li et al. (2014); Ngoc et al. (2020) to estimate the volume density, n(H 2 ). The mass of a region is calculated by M = βm H2 N total,H2 (D∆) 2 , where β = 1.39 accounts for He of 10% abundance in addition to H 2 in the total mass, m H2 is the mass of hydrogen molecules, D = 1.98 kpc is the distance to M17 (Xu et al. 2011), and ∆ = 4 is the pixel size of the N (H 2 ) map. Then the volume density n(H 2 ) is expressed as where R is the radius of a considered clump equal to R = (ab)/2 (a and b are the major and minor axes of an ellipse covering the clump (Li et al. 2014)). Equation 12 is deduced assuming spherical geometry of the molecular clump. As shown in Figure 3 and 7, we define two ellipses that encompass the observed M17-N and M17-S regions with corresponding angular radii of 69 and 120 , respectively. The physical radii of the two clumps are obtained by multiplying R with the distance D to M17. Employing Equation 12, the volume densities, n(H 2 ), are calculated, and the results are shown in Table 1 for both M17-N and M17-S regions. The uncertainty on n(H 2 ) is calculated assuming that its relative uncertainty is the same as that of N total,H2 , namely σ n(H2) n(H 2 ) = σ N (H2) N (H 2 ) .

RESULTS AND DISCUSSION
In this section, we first describe the magnetic field morphology of the observed region and report the magnetic strengths in the plane of the sky, B POS , using the DCF method for the M17-N and M17-S regions. We then explore the relative contribution of B-fields, gravity, and turbulence by calculating the mass-to-flux ratios and Alfvénic Mach numbers. Finally, we investigate dust grain alignment and disruption based on the polarization fraction in the observed regions.

Magnetic field morphology
The magnetic field morphology is imprinted on star formation processes. We thus can use the observational data to test theoretical predictions of star formation.
Overall, over the whole region, B-fields in the outer, low-density regions are perpendicular to the density structure. Figure 3 and its close-up version ( Figure A.3) illustrate well this point with many half-vectors in the low-density regions orthogonal to the intensity contours. Another prominent feature is that when the field lines pass through the H II region, east of the observed region (see Figure 1 of Povich et al. (2009) for a more precise location of the H II region), they tend to run parallel threading M17-N and M17-S. However, in M17-S, the fields run perpendicular to these fields behind the highdensity areas at the center (see Figure 8).
In the M17-N region, the main B-fields line up along the north-south direction through the highest density area and curve in the eastern and western sides of the region (see Figures 8 and A.3). The morphology of the B-fields in M17-N mimics the fields of a gravitationally collapsed molecular core with the presence of strong magnetic fields where we see the classical hourglass geometry with straight fields in the middle and curved ones at the two wings. This is a common structure which we often find in observations and simulations (see e.g., Kandori et al. 2018;Wurster & Li 2018;). On the other hand, M17-N is only a small part of a much larger M17-N region (see Figure 9) which is not fully covered by the current SOFIA/HAWC+ observations. Therefore, the observed structure may be just a coincidence with the fields following the emission by Polycyclic Aromatic Hydrocarbons (PAHs) at 8 µm wavelength observed by Spitzer (the bright pink color structure in the northern region in Figure 9). Distinguishing these two scenarios requires observations of higher angular resolution and larger areas of M17-N.
In the M17-S region, generally, the fields run perpendicular to the elongation of the high-density structure. In the region, the fields also form an asymmetric largescale hourglass shape (see Figures 8,9,and 10). In the southern side of the hourglass, just below Anon 1 and Anon 3, the fields run from east toward west while gradually bending the to south. The fields are highly curved in the most southern part of the region and perturbed in the high density regions at the center. On the northern side of the hourglass, just above UC 1 and CEN 92, we see the curved fields (see Figures 8 and 10). This hourglass structure seems to be clearly caused by the gravitational contraction of the massive cores at the center of M17-S. Moreover, the asymmetry is due to the complex distribution of matter in the center and north-western part of the region. Another prominent feature of M17-S, which is commonly found in PDR regions (see e.g., Pattle et al. (2018)), is the pillar structure. Here, we find a 'triangular' pillar with the top-end coincident with the positions of UC 1 and IRS 5 and a 'base' in the western direction (see Figures 8 and 10). It is interesting to note  (Hanson et al. 1997). There are several infrared sources present in M17-S such as UC 1, IRS 5, CEN 92, Anon 1, and Anon 3. The KW (Kleinmann-Wright) object is a binary-star system Kleinmann & Wright (1973). The locations of the sources are adopted from Lim et al. (2020).  Figure 9. A RGB image of M17 using Spitzer data as described in Figure 1. The white box is the observed area of SOFIA/HAWC+ toward M17. The half-vectors are plotted with uniform lengths using SOFIA/HAWC+ data and rotated by 90 • to trace the magnetic fields. In this figure, we can see the larger scale structure mentioned in the text connected to the M17-N region. Figure 10. Same as Figure 8 with the 'drapery' patterns produced by the line integral convolution (LIC) tool (Cabral & Leedom 1993).
that the magnetic field morphology found here is similar to the ones reported by Pattle et al. (2018) where the fields run parallel to the pillar in the H II region side and perpendicular to the pillar on and behind the base. Dotson (1996) reported that magnetic field lines are elongated into the cloud core and bulged away from the H II region, heated by OB stars. This could be an evidence that the H II region is expanding into its surrounding medium (Zeng et al. 2013). As a conclusion, our magnetic field morphology is consistent with the scenario that the magnetic field is distorted by H II region (see more in Section 4.5), several IR sources, and stellar clusters.

Magnetic field strengths
Given the estimated values of the polarization angle dispersion, σ θ , the one dimensional non-thermal velocity dispersion, σ ν , and the volume densities, n(H 2 ), in Table  1, we calculate the strengths of B-fields in the plane of the sky using Equation 8, yielding B POS = 980±230 µG and 1665 ± 885 µG for M17-N and M17-S, respectively (see more details in Table 2). The magnetic field strength of M17-S is greater than that of M17-N because of its higher density.
We note that the DCF method is applicable here since all the polarization angle dispersion is smaller than 25 • (Crutcher 2004). The uncertainties are propagated from the uncertainties on σ θ , σ ν , and n(H 2 ) using the following equation where δn(H 2 ), δ∆V, and δσ θ are the uncertainties on σ ν , ∆V = 2.355σ ν , and σ θ , respectively. Table 2. Results of the magnetic field strengths, Alfvénic Mach numbers, mass-to-flux ratios for M17-N and M17-S. γ is the inclination angle of the B-fields with respect to the line of sight.

BPOS Mach number Mass-to-flux ratio [µG]
MA λ M17-N 980 ± 230 0.12 sin γ 0.07 ± 0.04 M17-S 1665 ± 885 0.22 sin γ 0.28 ± 0.33 Previous study by Pellegrini et al. (2007) found the magnetic field strength around the southwestern part of the M17 photodissociation region with the peak value of ∼ 600 µG. Brogan et al. (1999) measured directly the magnitude of the B-fields along the line-of-sight (LOS) using the Zeeman observations of H I absorption lines toward the H II region and obtained B LOS ∼ −450-550 µG. Chen et al. (2012) carried out a rough estimation of the magnetic field strength based on the polarization of point sources in NIR and FIR using a sampling rectangle in the M17-S region, and they found the total magnetic field strength of B = B 2 LOS + B 2 POS ∼ 230 µG and the inclination angle of the magnetic field vector with respect to the plane-of-sky of ∼ 40 • . Therefore, our estimated values of B POS = 980±230 µG and 1665±885 µG with rather large uncertainties are still comparable with previous results.
Several factors affect the estimation of the magnetic field strengths. One is that we assumed spherical geometries of the molecular clumps to estimate the radii of the observed regions when estimating the volume densities. Depending on the real 3D geometry of the clouds, this assumption could lead to a large uncertainty in B POS values. In addition, polarization observations integrate over all of the structure along the line of sight, which leads to the reduction of σ θ , therefore, an overestimate of B-field strengths. Also it is important to note that, statistically, margin of overestimation by the DCF method on the magnetic field strengths can be up to a factor of 2 for an individual cloud (Crutcher 2004). There are efforts to improve the DCF method such as of Cho & Yoo (2016) who modified the method to reduce the overestimate of B POS strength by a factor equal to the ratio of the average line-of-sight velocity dispersion and the standard deviation of centroid velocities. We calculated the standard deviation of centroid velocities for the M17-N and M17-S which are equal to 0.59 km s −1 and 1.46 km s −1 , respectively. This means that the conventional DCF method overestimates the strength of the mean magnetic field in the plane-of-the-sky by a factor of (1.6/1.46)=1.1 and (2.3/0.59)=3.9 in M17-N and M17-S, respectively (the line of sight velocity dispersion is taken from Table 1). Therefore, the strengths of Bfields now are 891±209 µG and 427±227 µG for M17-N and M17-S, respectively.
The calculated B-field strengths have quite large uncertainties which are mostly propagated from the uncertainties on column densities averaged over large regions (see the results in Table 1). We, therefore, examine the B-field strength limiting to the highest density regions with N (H 2 ) > 10 22 cm −2 . With this new cut, only the densest regions of M17-S are retained (see Figure A.4) and the new results are N (H 2 ) = (7.6 ± 6.5) × 10 22 cm −2 , n(H 2 ) = (6.9 ± 5.9) × 10 4 cm −3 , while ∆V = (2.3 ± 0.1) km s −1 and σ θ = 6.5 • ± 0.5 • are the same as before the cut is applied. Therefore, the magnetic field strength obtained is B POS = 2134 ± 738 µG and λ = 0.27 ± 0.25 with smaller uncertainties as ex-pected. This calculation serves as a reference providing a more precise measurement of the average magnetic strength in the dense region of M17 and a sense of possible changes of its value depending on the choices of the regions of interest. In reality, the magnetic field strength is expected to change with the local gas density.

Alfvénic Mach number M A
The Alfvénic Mach number, M A , represents the relative contribution of turbulence to magnetic fields. M A is an important parameter for describing the evolution of molecular clouds (Kritsuk et al. 2017). A sub-Alfvénic Mach (M A < 1) means the cloud has a strong magnetic field while a super-Alfvénic (M A > 1) implies a weak magnetic field compared to turbulence. In the super-Alfvénic case, the magnetic field morphology is significantly affected by turbulence due to it subdominance. In this scenario, the morphology of the magnetic field is therefore expected to be random. M A can be calculated from the velocity dispersion following Padoan et al. (2001); Nakamura & Li (2008); Wang et al. (2019) as where v A is the Alfvénic velocity. Combining with Equation 8, we obtain where γ is the inclination angle of B-fields with respect to the line-of-sight in the range of [0 • , 90 • ] with B POS = B sin γ, and Q c = 0.5 (Ostriker et al. 2001).
Using the inferred polarization angular dispersion in Table 1, we obtained M A = 0.12 sin γ and 0.22 sin γ for M17-N and M17-S, respectively. Thus, M17 is sub-Alfvénic, implying that the magnetic fields in the region dominate turbulence. The sub-Alfvénic Mach numbers are also in agreement with the generally well-ordered magnetic field morphology in the region (see Figure 8).

Mass-to-flux ratio
The mass-to-flux ratio, M/Φ, refers the ratio of the mass to the magnetic flux. This is a crucial parameter describing the role of B-fields relative to gravity in star formation (Crutcher 2012). Usually, it is determined by the critical value, defined as follows (Crutcher 2004) where the critical mass-to-flux ratio (M/Φ) critical = 1 2π √ G (Nakano & Nakamura 1978), N (H 2 ) is the average column density of the considered region in units of cm −2 , and B POS is given in units of µG. 3 The errors on the mass-to-flux ratio are calculated using the following equation A super-critical value of λ > 1 means that the gravity dominates the magnetic field pressure, and the cloud can undergo gravitational collapse to form a protostar. Conversely, a sub-critical value of λ < 1 indicates that the magnetic field is strong enough to counteract the gravitational collapse (Crutcher 2004;Pattle et al. 2017).
We obtained λ = 0.07 and 0.28 for the M17-N and M17-S regions, respectively (see Table 2). Overall, the M17 cloud is sub-critical. It means that M17 is magnetically supported and belongs to the strong magnetic field model of star formation theory (Nakano & Nakamura 1978). However, the λ values here are calculated by averaging over the whole considered regions. Therefore, the highest density cores can still become super-critical and then gravitationally collapse to form new stars. The inferred values of the mass-to-flux ratio are compatible with theoretical predictions of the am-bipolar diffusion model where the B-field morphology is dragged inward from the outer layers of the cloud to the massive cores in both M17-N and M17-S. In addition, these results are compatible with the deficiency of forming massive stars and the lack of gravitationally bound clumps in the regions (Nguyen-Luong et al. 2020).

Stellar feedback from the massive star cluster
As shown in Figures 8 and 10, the magnetic fields in the vicinity of the star cluster (yellow symbols in Figure 10) appear to be bent and aligned with the density structure (see more in Section 4.1). This suggests the effects of stellar feedback of the massive star cluster.
To understand the importance of feedback, we estimate the ratio of the ram pressure, P ram , by stellar winds (or expansion of HII region) to the magnetic pressure, P mag , which is given by where the wind speed v wd ∼ 20 km s −1 is adopted from Pellegrini et al. (2007) Figure 11. Dust temperature map of M17 with the contours indicating the corresponding temperatures. In M17-S, the dust temperature tends to decrease from the east to the west, in the direction away from the star cluster.
2.8m H n(H 2 ). The volume densities, n(H 2 ), and magnetic field strengths, B, close to those from Table 1 and 2 have been used. The dominance of the ram pressure over the magnetic pressure as shown in Equation 18 implies that the winds can bend the magnetic fields frozen in the gas. Stellar feedback due to the star cluster may also trigger star formation in the M17-S region, as revealed by the presence of numerous IR sources. 4.6. Dust temperature map As described in Section 3.3, the dust temperature, T d , map has been derived via the graybody fit with Herschel 160-500 µm data. Then, the template T d -map has been re-gridded to match to SCUBA2 850 µm map where the final map has the pixel size of 4 . Figure 11 shows the temperature map of the observed region by Herschel. The average temperatures are 74 ± 11 and 43 ± 10 K for the M17-N, and M17-S, respectively. Globally, the M17 dust temperature gradually decreases from east to west and from north to south. This general trend is consistent with the scenario that dust is heated by the massive star cluster (see e.g., Figure 10). The higher temperatures in M17-N and the eastern and northeastern parts of M17-S arise from the fact that these regions are located next to the star cluster and thus strongly heated by this intense radiation source as well as the H II region in the east. M17-N is additionally affected by shocks from G015.128 resulting in an obvious temperature excess. In the M17-S, the dust temperatures decrease along the direction away from the star cluster, which stems from dust absorption within dense structures of an order of magnitude higher than that of M17-N. 4.7. Polarization fraction vs. emission intensity, column density, and dust temperature We now study the dependence of the polarization fraction, p, on the total intensity, I, the column density, N (H 2 ), and the dust temperature, T d , which reveal basic properties of grain alignment and disruption in M17.
Firstly, we investigate the relationship between the polarization fraction, p, and the total intensity, I. It is clearly seen from Figure 3 that polarization fraction decreases when intensity increases. Generally, the dependence of p ∝ I −α describes the variation of the grain alignment efficiency and the magnetic field geometry across the cloud. For a uniform magnetic field, the slope α = 1 implies that the grain alignment is present only in the outer layer of the cloud and becomes completely lost in the inner region. Figure 12 shows a fitted power-law model p ∝ I −α with the power index α = 0.51 ± 0.01 for the whole M17 region. The values of α do not change much when we fit only with pixels from M17-N (α = 0.54 ± 0.02) or M17-S (α = 0.55 ± 0.01) separately. A best-fitted value of α = 0.51 reveals that grain alignment is significant toward the high emission intensity. As a conclusion, we observe that the polarization fraction, p, decreases with increasing the emission intensity, I.
To better understand the variation of p, we now study the relationship between the polarization fraction, p, and the column density, N (H 2 ), as well as the dust temperature, T d . Figure 13 (top-left) shows the p − N (H 2 ) relation. Here we derive the slopes by fitting piece-wise linear functions to the data. In the M17-S region, the polarization fraction gradually decreases with increasing the gas column density, before it experiences a steep drop with a slope of −0.71 for N H2 > 3 × 10 22 cm −2 or the visual extinction 4 A V > 37 mag for the typical total-to-selective extinction R V ≈ 3.1. In the M17-N region, a steep decrease with slope of −0.67 occurs at N (H 2 ) > 5 × 10 21 cm −2 or A V > 6 mag. The strong depolarization occurs in the region where the gas density is relatively low and the dust temperature is relatively high compared to that of M17-S (see Figure 11), which is unexpected from the RATA theory.
The top-right panel of Figure 13 shows the p − T d relation. It appears that the polarization fraction first decreases, then increases before decreasing again as the dust temperature increases. The decrease-increase feature originates from the M17-S region, while the depolarization is from M17-N.
The bottom panel of Figure 13 shows the gas column density as a function of the dust temperature. The column density decreases rapidly with increasing T d in M17-S, but it varies slowly in M17-N. The first decrease of p for T d < 40 K corresponds to the highest gas density (A V ∼ 100 mag). This depolarization is caused by decrease of grain alignment due to the attenuation of the radiation field and the enhancement of collisional damping. Towards higher dust temperatures, the gas density drops, and the grain alignment efficiency enhances toward the central luminous source, which results in the following increment of p. These features are expected from the context of the RATA theory. In M17-N, the dust temperature could be up to 150 K and the gas density becomes more diffuse at A V < 10 mag, the polarization fraction, however, appears to monotonically decrease to higher T d , which is completely contradictory to what expected from RATA theory.

Implications for grain alignment and rotational disruption by RATs
We now discuss the implications of the observed polarization fraction toward M17 for physics of grain alignment and disruption based on radiative torques.
According to the RATA theory (see Lazarian & Hoang 2007), the polarization fraction of thermal emission increases with decreasing alignment size, a align , which is the minimum size of aligned grains. The alignment size is determined by the balance between spin-up by RATs and spin-down by gas collisional damping, which is a function of the local dust temperature (or radiation intensity) and gas density as a align ∼ n 2/7 H T −12/7 d . As a result, a denser gas and lower dust temperature increase a align (see details in Hoang et al. 2021), which results in the decrease of the polarization fraction. Such a prediction by the RATA theory can explain the decrease of the polarization fraction with increasing the gas column density observed in the M17-S region (see Figure 13; left panel). However, it can not explain the decrease of p with increasing T d in the M17-N region where the gas density changes slowly (see Figure  13; right panel), which reveals evidence of the RATD.
The polarization fraction at far-IR/submm is very sensitive to the maximum size of the grain size distribution because large grains dominate long-wavelength emission. Following the RATD mechanism (Hoang et al. 2019), the maximum grain size above which large grains are disrupted by RATs is determined by the local dust temperature (radiation field), gas density, and the grain tensile strength, S max , which follows by a disr ∼ n 1/2 H T −3 d S 1/4 max . For the average estimated volume density n H2 3 × 10 4 cm −3 (see Table 1), and the temperature at the breaking point in Figure 13 of T d 53 K, one obtains a disr = 0.16, 0.28, and 0.49 µm for S max = 10 7 , 10 8 , and 10 9 erg cm −3 , assuming the mean wavelength of radiation field of 1 µm. The decrease of the disruption size by RATD leads to the decrease of the dust polarization as the dust temperature increases (Lee et al. 2020;Tram et al. 2021b), which successfully reproduces the observed trend toward the M17-N (see Figure 13).
To quantify the magnetic field tangling at small scales on the depolarization, we calculate the polarization angle dispersion function, S, (see Section 3.3 in Planck Collaboration: et al. 2015). For each pixel at location x, S is calculated as the standard deviation of the polarization angle difference, S xi , of pixel x and pixel i which lies on a circle having x as the center and a radius of δ where S xi = θ(x) − θ(x + δ) is the polarization angle difference and N is the number of pixels lying on the circle. For the current SOFIA/HAWC+ data set, we calculate S for δ = 27.2 (∼ 2 beam sizes). The left panel of Figure 14 shows the S -T d relation where it can be seen that the angular dispersion function S strongly correlates with the dust temperatures but not the highest ones, T d > 90 K, in M17-N. Thus, the decline of p at T d > 90 K (see the upper-right panel of Figure 13) cannot be explained by the magnetic field tangling and presents an evidence for the RATD effect. The right panel of Figure 14 shows a general anti-correlation of the angular dispersion function S and polarization fraction p, except, again, for the highest temperatures pixels (red dots) whose temperatures are greater than 90 K. For T d < 35 K located in M17-S (blue dots in Figure 13 and 14), we can see that the rapid decrease of p is due to B-field tangling in the high density region (see bottom panel of Figure 13). For the temperature range from 35 K to 65 K, in both M17-N and M17-S, N (H 2 ) and S decrease, therefore, p increases. This clearly shows the important contribution of field tangling to the depolarization at low dust temperatures. A detailed modeling of dust polarization using the grain alignment and disruption by RATs for M17 is beyond the scope of this paper, and will be addressed in a follow-up study.

SUMMARY
In this study, we use the dust polarization data taken by SOFIA/HAWC+ at 154 µm to study magnetic fields and dust physics in the M17 nebula. Our main results are summarized as follows: 1. Using the DCF method, we estimated the magnetic field strengths to be 980 ± 230 µG in the lower-density region (M17-N) and 1665 ± 885 µG in the higher-density region (M17-S). The strong magnetic field could be a result of the pressure exerted by the H II region in the eastern part of the observed region by SOFIA/HAWC+ toward M17. In the M17-N, the B-field morphology can be that of a gravitational collapse molecular core. The morphology of B-fields in M17-S displays a well-organized elongated structure along with the gravitational collapse directions from the outer regions to the denser central regions. The fields are dragged inward to the gravitational center. The whole B-field morphology represents an asymmetric large-scale hourglass structure. We also found a pillar structure which is one of the common features of PDR regions. A rough estimation of the contribution of ram over the magnetic pressure suggests that the wind from the PDR region is strong enough to impact on the B-field morphology of M17.
2. The Mach numbers are determined to be sub-Alfvénic (M A < 1) which indicate that the magnetic field dominates turbulence. In addition, the inferred sub-critical values of mass-to-flux ratios, λ < 1, imply that the magnetic fields in the regions are strong enough to resist gravitational collapse. These results are consistent with the deficiency of the formation of massive stars in the region from previous studies.
3. To study dust physics, we analyzed the relation of the polarization fraction, p, with the emission intensity, I, gas column density, N (H2), and dust temperature, T d . The power index of the p vs I relation, α = 0.51, implies that dust grains could still be aligned by radiation in the region. The decrease of the dust polarization with the column density in the M17-S can be explained by the RATA theory as well as the tangling of the magnetic field.
4. To study the effect of magnetic field tangling on the dust polarization, we also analyzed the variation of the polarization angle dispersion function, S, with dust temperature and gas column density. In the M17-N, the decrease of p with T d at high temperatures when both N (H 2 ) and S decrease is most consistent with the theoretical predictions of dust polarization by both RATA and RATD effects.
5. There are large statistical biases on the estimation of the polarization angle dispersion, σ θ , as well as the volume density, n(H 2 ), due to the fact that the magnetic field strengths are estimated over a large areas. These biases lead to large uncertainties on the measurement of the magnetic field strengths. It is also known that for a random sample of magnetic field orientations, using the DCF method the mass-to-flux ratio will on average be overestimated by a factor ∼ 3. In fact, using a method proposed by Cho & Yoo (2016) we found the overestimate factors equal to 1.1 and 3.9 for the M17-N and M17-S, respectively. Therefore, new methods of estimating magnetic field strength which are able to identify the fields and accompanied useful parameters such as Alfvénic Mach numbers and mass-to-flux ratios in the basis of pixel by pixel would provide detailed and more precise conditions of star formation in the M17 molecular cloud.

APPENDIX
A. APPENDIX This appendix shows the characteristics of the observational data and the behavior of the master cut applied to the data. Figure A.1 shows the distributions of the SNRs of I, p, and I p before and after applying the master cut defined in the text. Their corresponding mean and RMS values are listed in Table A.1. The mean values of the SNRs increase significantly after the cut by 45%, 18%, and 17% for I, p, and I p , respectively.  Figure A.2 presents one dimensional histograms of the raw data of I, σ I , p, σ p , I p , σ Ip , θ, σ θ . The associated mean and RMS values of these quantities presented in Table A.2.