Analysis of Galactic molecular cloud polarization maps: a review of the methods

The Davis-Chandrasekhar-Fermi (DCF) method using the Angular Dispersion Function (ADF), the Histogram of Relative Orientations (HROs) and the Polarization-Intensity Gradient Relation (P-IGR) are the most common tools used to analyse maps of linearly polarized emission by thermal dust grains at submilliter wavelengths in molecular clouds and star-forming regions. A short review of these methods is given. The combination of these methods will provide valuable tools to shed light on the impact of the magnetic fields on the formation and evolution of subparsec scale hub-filaments that will be mapped with the NIKA2 camera and future experiments.


Introduction
This article gives a short review of the main methods used to analyse linear polarization maps. Given the more and more prolific literature on the subject, not all the articles using each technique are cited but the generic articles introducing the conceptual ideas behind each procedure are used as references. Additional methods employing wavelet transform analysis technics and Bayesian inference statistical tools currently in development are not discussed here, and we refer the interested reader to the method investigated by [1] and to the IMAGINE consortium project [2], respectively. The review focuses on the analysis of maps obtained from the observation of linearly polarized thermal dust emission at submillimeter (submm) wavelengths toward Galactic Molecular Clouds (MCs) and protostellar cores, but the same methods can be used on maps obtained from simulations. Independently of the dust grain alignment mechanism that is considered (see [3] for a review on this subject), the main accepted current picture is that of dust grains aligned perpendicular to the local magnetic field direction pervading the Interstellar Medium (ISM). Each polarization pseudo-vector displayed in one pixel of a polarization map is therefore an average measurement of the weighted contribution by all dust grains along a given Line-Of-Sight (LOS) in a direction perpendicular to the average magnetic field on the Plane-Of-Sky (POS). From the measurements of the Stokes parameters I, Q and U the total fraction of polarization (p = √ (Q 2 +U 2 ) I ) is often represented by the length of the pseudo-vector and the polarization angle (P.A.; θ = 1 2 × arctan( U Q )) by its orientation with respect to a given reference frame. Other representations of p exists in the literature and maps showing only drapery patterns of the POS magnetic field lines, or of the P.A.s, are more and more common. An example of submm polarization is shown in figure 1 (left panel).  Figure 1 from [5] showing the Submillimeter Array (SMA) 870 µm polarization map of the Collapsing core W51 e2 obtained by [4]. The thick red segments show the magnetic field orientation after rotation by 90 • of the polarization pseudo-vectors (see Introduction) where polarization was detected such that p/σ p > 3. The degree of polarization is not represented in this map. The blue pseudo-vectors show the gradient directions (see Section 3 and Section 4) of the dust emission continuum mapped at 1.3 mm by [6] with the Berkeley-Illinois-Maryland Association (BIMA). Right: Figure  2 from [14] showing the ADF of the SCUBA JCMT DR21 850 µm polarization maps (see Section 2).

Davis-Chandrasekhar and Fermi (DCF) method and Angular Dispersion function (ADF)
The DCF method was introduced in the middle of the twentieth century by [7] and [8]. It was first designed to get estimates of the POS magnetic field strength, B pos , assuming a turbulent diffuse ISM. It was applied on polarimetry data obtained on star fields by [9] at visible wavelengths. In this regime, the polarization is produced by dichroic absorption of starlight by dust grain layers pervading the diffuse ISM, in a direction perpendicular to the one produced in emission at submm wavelengths by identical polarizing dust grains (see [10]). Inferred from the data is a large-scale uniform magnetic field along the Galactic Plane (GP). The fluctuations around the mean of the distribution of the polarization pseudo-vectors are assumed by [7] and [8] to be produced by Alfvén waves that are coupled to the gas such that there is equipartition between kinetic and perturbed magnetic energies. With this model B pos is a function of the ISM gas density ρ, of the gas velocity dispersion σ v , and of the polarization angle dispersion σ θ , such that: B pos = Q 4πρ σ v σ θ , where Q is a factor of proportionality. With this method, [8] calculated a diffuse ISM B pos estimate of a few µG consistent with those obtained with other independent methods. Later on, in the 1990s, when polarimetry detectors at submm wavelengths started to be sensitive enough, the DCF method was used on submm maps of molecular clouds, meaning transposed to spatial scales 1000 to 10000 times smaller than the GP scale in regions of density two orders or more magnitudes that of the diffuse ISM density. The first comprehensive study was done by [11] from 10 measurements obtained at 100 µm with the Kuiper Airborne Observatory (KAO) in the Orion Nebula, leading to B pos estimates of a few mG. The DCF method has been tested numerically by [12] and [13] and some refinements to the calculations of B pos have been proposed. A review on the values of Q has been given by [14].
Improvements to the Method: Further major improvements to the method were designed by [15], [16] and [17]. The Angular Dispersion Function (ADF) was introduced by [15] to avoid inaccurate estimates of magnetohydrodynamic or turbulent dispersion, as well as to avoid inaccurate estimates of B pos , due to the large-scale, non turbulent field structure.
is the angle asssociated to the projected POS magnetic field vector B(x) at position x in a map. The difference in angle between two points is obtained by ∆Φ(l) ≡ Φ(x) − Φ(x + l), and is calculated between the N(l) pairs of vectors separated by displacement, or lag, l. < ... > denotes an average and l = |l|. The square of the ADF, a second order structure function, is also often used (e.g. [13]). One example of ADF is shown in Figure 1 (right panel), where the angular dispersion b is the intercept of the fit at l = 0. Correlations in polarization angles at lags l smaller than the telescope beam (1.22λ/D) or than the turbulent correlation length (δ) have to be avoided. The fit is ideally applied on the set of points calculated by taking into account the measurement uncertainties and such that the lag distance l is smaller that the typical length scale d for variations in the large-scale magnetic fields. Once b is estimated from a map this method also provides an estimate of the turbulent to large-scale magnetic field strength ratio such that The method to take into account the effect of the signal integration through the thickness of the clouds as well as across the area sustended by the telescope is fully incorporated by [16]. The authors also show how to evaluate the turbulent magnetic field correlation length scale from polarization maps obtained with sufficiently high resolution and high enough spatial sampling rate. Further examples are given and discussed by [17] as well as the application of the technique to interferometry measurements.
Results: The DCF method has been applied to data from many Galactic MCs obtained on sky patches, including Gould belt MCs (e.g. [18]) and some of the closests low, intermediate or high star-forming regions. Estimates of B pos lie typically in the range of a few µG to a few mG. In OMC-1 the turbulent correlation length is estimated to δ ≈ 10 mpc (e.g. [16]). Independently of the uncertainties coming from the propagation of the errors, the estimates of B pos and δ will rely on the choice (or availabilty) of the gas tracers and of the value of Q used for the calculation of B pos (e.g. see discussion in [19]). In principle accurate and reliable results can be obtained with polarization data of sufficient spatial resolution and high enough spatial sampling rate ( [16]). The DCF method has been applied to a large fraction of the sky by [20] and on the full sky by [21]. Using the Planck Release 3 [21] have obtained the following relation between the ADF, S , and the fraction of polarization, p, as a function of the map resolution, w: < S p >= 0 • .31 p w 160 ′ . The results are displayed in Figure 2 (top-left panel), and show that down to a resolution of 10 ′ the systematic decrease in p with N H is determined mostly by the magnetic-field structure. At a lower resolution of 2.5 ′ , using the 500µm BLASTPol polarization map of Vela C, [22] discuss the dependence of p on the dust temperature and on N H and show that p ∝ N −0.45 H S −0.60 , suggesting that dust grain alignment properties may also contribute to the decrease of p in some conditions.

Histogram of Relative Orientations (HROs) between ISM tracer structures and magnetic field structures
This method has been designed by [23] and subsequently applied to many molecular cloud regions (e.g. [24], [25]). The concept of the method relies on the calculation of the gradient of the column density structure provided by a given material tracer (e.g. N H ) in the pixels of a polarization map. As an illustration the blue pseudo-vectors displayed in Figure 1 (left panel), show the gradient orientations obtained from a dust emission continuum map. Once the gradients are quantified their relative orientation with the POS magnetic field orientations inferred from the P.A.s can be estimated and the HROs can be analysed. Several methods exist to calculate the gradients in a map as a function of the morphology structure of the  Figure 11 from [21] showing the S × p relation as a function of column density N H where w is the map resolution parameter (see Section 2). Bottom-Left: Figure 6 from [27] showing the PRS Z x corrected for oversampling obtained with several molecular tracers (see Section 3). Right: Figure 3 from [5] showing the 90 • rotated P.A. α, and the angle associated to the gradient orientation ψ, used to retrieve information on the magnetic field strength and significance (see Section 4).
object one is interested to isolate (see e.g. [26] and references therein). HROs can be built by considering the relative orientation of P.A.s with respect to the average orientation of a large-scale structure identified in a map (e.g. [26]). More frequently, HROs are calculated by considering the relative orientation angle φ between the POS magnetic field < B pos > and a line tangent to the local iso-contour (see [23], [24], [27]) which is equivalent to the angle between the polarization direction E and the intensity gradient ▽I: φ = arctan(| ▽ I × E|, ▽I. E). Statistical measures of HROs have been quantified using the histogram shape parameter (e.g. [23], [24]), ζ = A 0 −A 90 A 0 +A 90 , where A 0 is a measure of the total number of points in the quartile 0 • < φ < 22 • .5, and A 90 a measure of the same quantity in the quartile 67 • .5 < φ < 90 • . Smoother and more accurate definitions of ζ have been investigated by [28]. For this Rayleigh statistics Z is used to test whether a given set θ i of n independent angles are uniformly distributed within the range[0, 2π]. Using the relation θ = 2φ, where φ is the relative orientation angle discussed above, [28] showed that the Projected Rayleigh Statistic (PRS) of Z: Z x = Σ n ind i cosθ i / √ n ind /2, where n ind is the number of independent data samples in the map, can be used to test for a preference for perpendicular or parallel alignment between the magnetic field orientations and the iso-contours. The variables ζ or Z x are often calculated on samples of data lying in different ranges of intensity of the material tracer map, and used to explore variations of the relative orientations between the magnetic field structures and the ISM morphology structures as a function of the column density or number density parameters. Figure 2 (bottom-left panel) shows the PRS obtained by [27] using several molecular tracers. For each tracer, Z x > 0(< 0) indicates the I structure preferentially aligns parallel (perpendicular) to < B pos >.
Results: The main trend found from the analysis of HROs derived with lines tangent to local iso-contours is that the POS magnetic field orientations change from mostly parallel (or not clearly defined) to perpendicular with respect to the molecular cloud structures probed with a given material tracer ( [24], [25], [27]). This translates by a change of orientation from lower to higher column densities N H . In Vela C the transition is estimated to occur at molecular hydrogen number density of approximately n H ≈ 10 3 cm −3 ( [27]). This local iso-contours approach has been tested by simulations (e.g. [23]). In their study [26] first produce a component separation analysis of the magnetic fields in the diffuse ISM and in higher column density regions, and use image analysis technics to extract filaments and clumps embedded in different background column densities, i.e. showing density contrast varying with their environment. Their analysis of the HROs obtained between filaments, embedded clumps and internal and background magnetic field orientations lead to a more complex picture than in other studies. Overall their results support the possibility of magnetic fields strong enough to influence the formation of molecular clouds and also of their embedded clumps.

Polarization-Intensity Gradient Relation (P-IGR)
This method has been designed by [5] to study star-forming regions. In the case of negligible viscosity and infinite conductivity (ideal MHD case) the force equation is given by the fol- where ρ and v are the dust density and velocity, respectively. B is the magnetic field, P is the hydrostatic dust pressure, φ is the gravitational potential resulting from the total mass contained in the region of interest and, ▽ denotes the gradient. The left-hand side term in the equation represents the resulting motion of the dust produced by the right-hand side terms which are the gradients of the hydrostatic pressure terms of the gas, the magnetic field, and the gravitational potential together as well as the magnetic field tension term (last term). The force equation can be transformed and under several assumptions the magnetic field strength can in principle be derived geometrically at each position of a polarization map and expressed by: B = sin(ψ) sin(α) (▽P + ρ ▽ φ)4πR, where α = P.A.-90 • , and ψ is the angle associated to the gradient orientation. Figure 2 (right panel), illustrates the terms displayed in this equation. The red and blue pseudo-vectors displayed in the figure can be compared to those displayed in figure 1 (left panel) where the gradients in intensity are estimated assuming central symmetry towards the brightest pixel in the map, but the method can be generalized to arbitrary cloud shapes. An important outcome of the method is given by the magnetic field significance: which gives a direct physical meaning to the factor sin(ψ) sin(α) derived geometrically from a map. Results: The method provides a quantification of the local significance of the magnetic field force compared to the other forces in a model-independent way. In W51 e2 it allows derivations of the azimuthally averaged radial profile B(r) ∼ r −1/2 ( [5]). The potential of the method is explored in additional works (e.g. [29] and references therein). In their study, [30] propose that in G34 the varying relative importance between magnetic field, gravity, and turbulence from large-to-small-scale drives and explains the different fragmentation types seen at subparsec scales (no fragmentation, aligned fragmentation and clustered fragmentation).

Conclusions and perspectives
The methods discussed above are complementary to each other. If submm polarization maps are available at different resolutions they are suited to explore the role of magnetic fields, turbulence and gravity on different spatial scales. Their combination starts to give insights on the interplay between magnetic fields, gravity and turbulence (e.g. [30]) and they are promising tools to shed light on the physics of hubs-filaments (e.g. [31]) detected at subparsec scales in molecular clouds. In addition to these methods, getting estimates of the mean inclination with respect to the LOS of the large-scale ordered magnetic field considered as the main factor regulating the mean level of polarization in a map can help guiding the overall interpretation. In this regard, a method first proposed by [14] has been further explored by other authors (e.g. [34]). Multi-wavelengths submm polarimetry of a region can give insights about polarization dust grain properties and also help to put constraints on the interpretation of the maps. On this matter we refer the reader to [3], [32], [33] and references therein.