Measurement of Interstellar Magnetization by Synchrotron Polarization Variance

Since synchrotron polarization fluctuations are related to the fundamental properties of the magnetic field, we propose the polarization intensity variance to measure the Galactic interstellar medium (ISM) magnetization. We confirm the method’s applicability by comparing it with the polarization angle dispersion and its reliability by measuring the underlying Alfvénic Mach number of magnetohydrodynamic turbulence. With the finding of the power-law relation of A∝MA2 between polarization intensity variance A and Alfvénic Mach number M A, we apply the new technique to the Canadian Galactic Plane Survey data, achieving the Alfvénic Mach number of the Galactic ISM. Our results show that the low-latitude Galactic ISM is dominated by sub-Alfénic turbulence, with M A approximately between 0.5 and 1.0.


INTRODUCTION
It is generally accepted that interstellar medium (ISM) is turbulent and magnetized, and is ubiquitous in the universe (Armstrong et al. 1995;Elmegreen & Scalo 2004).Magnetohydrodynamic (MHD) turbulence plays a crucial role in many key astrophysical processes, such as star formation, cosmic ray acceleration and propagation, as well as magnetic reconnection (e.g., Larson 1981;Lazarian & Vishniac 1999;Schlickeiser 2003;McKee & Ostriker 2007;Lazarian et al. 2015).Determining the properties of magnetized turbulence is a crucial step in exploring astrophysical activity.However, the measurement of ISM magnetic field strength is extremely challenging.
At present, several techniques measuring magnetic field strength have been proposed.As is well known, the Zeeman effect caused by spectral line splitting is a classic method in measuring the magnetic field component along the line of sight (LOS), requiring high instrument sensitivity and a longer observational integration time (Goodman et al. 1989;Crutcher 1999;Bourke et al. 2001;Crutcher et al. 2010;Mao et al. 2021;Inoue 2023).The traditional Faraday rotation synthesis can well present the magnetic structure of the Milky Way and nearby galaxies.However, this jfzhang@xtu.edu.cn(JFZ), hpxiao@xtu.edu.cn(HPX), jcho@cnu.ac.kr(JC), xjyang@xtu.edu.cn(XJY) method cannot determine the magnetic field direction due to the ambiguity of the phase angle (Burn 1966;Beck & Wielebinski 2013).The Davis-Chandrasekhar-Fermi (DCF, Davis 1951;Chandrasekhar & Fermi 1953) method utilizes fluctuations in observed polarization and gas Doppler shift velocity to estimate the projected magnetic field strength by B ∝ f 4πρδv/δϕ, where f is an adjustable factor (typically f = 1/2), δv is the LOS velocity dispersion, and δϕ is polarization angle dispersion.In general, this method results in the measured magnetic field strength deviating from the expected value.Note that this method is based on the assumption of incompressible isotropic turbulence.In fact, the theoretical and observational results imply that ISM turbulence is compressible and anisotropic (Goldreich & Sridhar 1995, henceforth GS95;Cho & Lazarian 2002;Spangler & Spitler 2004).
Recently, different improved versions of the DCF method have been enhanced.
For instance, Lazarian et al. (2022) proposed a differential measure approach where the adjustable factor in the DCF method is related to turbulence properties rather than a constant.In the framework of the DCF method, different effects have been included into the DCF method such as the average effect of independent eddies in the LOS direction (Cho & Yoo 2016), atomic arrangement effect (Pavaskar et al. 2023), and polarization angle difference as a function of the distance between a pair of measurement points (Falceta-Gonc ¸alves et al. 2008;Hildebrand et al. 2009;Houde et al. 2009).For the latter, due to the angle dispersion function affected by the integration effect along the LOS, the spatial range to measure should be limited within scales smaller than 0.1 pc (Cho & Yoo 2016;Liu et al. 2021).
With the modern understanding of the MHD turbulence theory, other techniques for measuring the strength of magnetic fields have been proposed.First, the anisotropy of the structure functions has been proposed considering the velocity channel (Esquivel & Lazarian 2011;Esquivel et al. 2015), and the velocity centroid and rotation measurement (Xu & Hu 2021;Hu et al. 2021).Second, gradient measurements involve velocity (Lazarian et al. 2018) and synchrotron polarization (Carmo et al. 2020).Third, the ratio of the E and B modes arises from dust polarization (Cho 2023).
In the earlier numerical work, Zhang et al. (2016) proposed the synchrotron polarization dispersion for measuring the scaling index of underlying MHD turbulence (see also Lazarian & Pogosyan 2016 for theoretical predictions).This method is named polarization frequency analysis (PFA) which considers the variance of polarization intensity as a function of the square of wavelength λ2 .They found that the turbulent spectral index can be successfully recovered from the polarization intensity variance in the region dominated by the mean field.In this work, we want to know whether polarization intensity variance can be used as a complementary method for measuring magnetic field strength.At the same time, to enhance the reliability of our techniques, we will also explore the use of polarization angle dispersion to measure magnetization levels. 1he structure of this paper is regularized as follows.In Section 2, we provide theoretical fundamentals including the properties of MHD turbulence, synchrotron polarization radiation, and statistical measure techniques.Section 3 describes simulation methods with basic parameter setup.Numerical results and application of observational data are presented in Sections 4 and 5, respectively.Sections 6 and 7 give discussion and summary, respectively.

Modern understanding of MHD Turbulence
The modern understanding of the MHD turbulence theory can be traced back to the work of GS95, which focused on incompressible and trans-Alfvénic (with Alfvénic Mach number M A ∼ 1) turbulence and predicted the scale-dependent anisotropy expressed as l ∥ ∝ l 2/3 ⊥ , where l ∥ and l ⊥ are the scales parallel and perpendicular to the local magnetic fields, respectively.Essentially, the turbulent cascade in the perpendicular direction and wave-like motion in the parallel direction are mutually coupled, establishing a critical balance condition by between turbulent eddy turnover time and the Alfvén wave period.Here, V A = B/ ⟨4πρ⟩ is the Alfvénic velocity related to the magnetic field of B and gas density of ρ.
With this critical balance assumption, GS95 found that the cascading rate v 3 l /l remains constant, i.e., v l ∼ l 1/3 , which results in the classical Kolmogorov's spectrum of E(k) ∼ k −5/3 (Kolmogorov 1941).
In later studies, GS95's work, in relation to trans-Alfvénic turbulence, was generalized to two scenarios: (Lazarian & Vishniac 1999;Lazarian 2006), where V L is the injection velocity at the injection or outer scale L inj .For the former, M A < 1, Lazarian & Vishniac (1999) found a transition scale In the range of [L inj , l trans ], with l ∥ remaining unchanged while l ⊥ continuously decreasing until reaching the critical equilibrium condition at the scale l trans , the turbulent cascade presents the weak turbulence properties.In the range of [l trans , l diss ], where l diss is the turbulent dissipation scale 2 , turbulence cascade shows the strong turbulence properties, with the parallel and perpendicular scale relationship of and with velocity and scale relationship Note that when M A = 1, Eqs. (3) and (4) will return to the original GS95 relationship of l ∥ ∝ l 2/3 ⊥ and v ⊥ ∝ l 1/3 ⊥ , respectively.
For the latter, M A > 1, Lazarian (2006) predicted that this super-Alfvénic turbulence has the transition scale of where the anisotropy of turbulent eddies occurs.When M A >> 1, the magnetic field is not dynamically important at the maximum scale, and within the scale range of [L inj , l A ], turbulence follows an anisotropic Kolmogorov cascade (see Lazarian (2006) for more details).

Synchrotron Polarization
When relativistic electrons move in a ubiquitous magnetic field, the synchrotron radiation process inevitably arises.With an assumption of a uniformly isotropic power-law energy distribution N e ∝ E 2α−1 for relativistic electrons, where E is the electron energy and α is the photon spectral index related to the electron spectral index p by the relation of α = (p − 1)/2, we have the synchrotron intensity I per unit frequency interval dν (Ginzburg & Syrovatskii 1965) where L represents the integration length along the LOS.Based on Eq.( 6), we can get the linearly polarized intensity by The Stokes parameters Q and U are related to the polarization intensity by Q = P cos 2ϕ 0 , U = P sin 2ϕ 0 , where ϕ 0 is the polarization angle.Defining the complex polarization vector we can write the polarization intensity as Furthermore, when synchrotron radiation experiences a radiative transfer, the Faraday rotation effect will result in the depolarization of the radiation signal.The resulting polarization angle can be re-written as where the rotation measure (RM) is given by Here, n e in units of cm −3 is the thermal electron density, B ∥ in units of µG is the parallel component of the magnetic field along the LOS, and dz in units of pc is the path length.

Statistical Measure Techniques
Since magnetic field fluctuations cause a change in polarization signal, the properties of the magnetic field can be explored through analyzing the polarization observational information.According to Zhang et al. (2016), we define the correlation function of the polarization intensity as where R 1 and R 2 denote the spatial locations on the plane of the sky.Considering the correlation of polarization intensities at the same spatial separation ∆R = R 2 − R 1 =0 on the plane of the sky, we name it as the polarization variance, i.e., one-point statistics.Practically, we normalize the polarization intensity variance A in units of the mean of the polarization intensity.This procedure can eliminate the influence of dimensionality on the statistics for realistic observations in different astrophysical environments.Moreover, according to the analysis of the DCF method (e.g., Lazarian et al. 2022), we define the polarization angle dispersion by With the above two definitions (Eqs.( 12) and ( 13)), this paper will explore how A and B connect with the magnetization M A by the Bayesian analysis method.To implement Bayesian statistics and fitting algorithms, we use the Python package PyMC to solve general Bayesian statistical inference and prediction problems (Patil et al. 2010).With the PyMC, we can determine the best fitting level and fitting error from two variables.From a theoretical view, with fundamental parameter assumptions such as the parameter Y depending on the observable X, a linearly fitted parameter model β, and a normal distribution error ϵ, we have the linear regression relationship of Y = βX + ϵ.Therefore, this model can be represented as a probability distribution Y ∼ N(βX, σ 2 ), where Y denotes a normally distributed random variable, βX a mean of linear predictor and a square variance of σ 2 .
In practice, we consider the measured quantity M A as the independent variable for linear fitting, and direct statistical quantities (A and B) as the dependent variables.The assumed distributions and parameters of all priors in this linear model are listed as follows We use Markov Chain Monte Carlo methods in Bayesian analysis, after choosing 20, 000 posterior 3 , with each sample including 4 parallel chains for 2, 000 tune, the output with the highest density interval 94% can provide the mean, standard deviation, and intercept of the linear fitting slope.

GENERATION OF MHD TURBULENCE SIMULATION DATA
To simulate magnetized turbulent ISM, we consider a set of governing equations of MHD turbulence as follows: where ρ is gas density, t the evolution time of the fluid, p = c 2 s ρ gas pressure, J = ∇ × B the current density, and f a random driving force.At the same time, we use the isothermal state equation to close the above equations.
When numerically solving the above equations, we use a single-fluid, operator-split, and staggered-grid MHD Eulerian code ZEUS-MP/3D (Hayes et al. 2006).We set a non-zero uniform magnetic field B 0 along the x-axis, resulting in a total magnetic field B = B 0 + δB, where δB is a fluctuating magnetic field, and an average density ⟨ρ⟩ ≃ 1.With the periodic boundary conditions and solenoid turbulence driving, we run two sets of simulations with different resolutions: Group A with a numerical resolution of 792 at a driving wavenumber of k inj ≃ 2.5, and Group B with the resolution of 480 at k inj ≃ 3.5, by changing the uniform magnetic field B and injection energy.When reaching a quasi-steady state at t ≃ 15 t A (t A being Alfvénic timescale), we terminate the simulation and then output the information of the last snapshot, which includes the 3D velocities, magnetic fields, and density.Using the definitions of Alfvénic Mach number M A = ⟨V/V A ⟩, sonic Mach number of M s = ⟨V/c s ⟩, and the plasma parameter β = 2M 2 A /M 2 s , we characterize different data groups as listed in Table 1.The latter, β, describes the compressibility of turbulence, with a low β corresponding to the compressible regime dominated by magnetic pressure, and a high β to the nearly incompressible one dominated by gas pressure.

Images of polarization intensity and polarization angle with Faraday rotation
To generate synchrotron polarization observations, we use typical parameters from the Galactic ISM, such as thermal electron number density of n e ∼ 0.01 cm −3 , spectral index of relativistic electron of p = 2.5 (note that change in the index p cannot affect our results; see Zhang et al. (2018) for more details), magnetic field strength of B ∼ 1.23 µG, Faraday integral depth of L box = 1 kpc.Furthermore, we consider the z-axis along the LOS direction, which results in the plane of the sky parallel to the x-y plane.Based on Section 2.2, we obtain the images of the Stokes parameters I, Q, and U using the above parameters.
As is well known, when a linearly polarized electromagnetic wave propagates along a magnetic field, the polarization direction will change due to the influence of Faraday rotation.Equation (10) indicates that the change in polarization angle is wavelength-dependent.To understand the effect of Faraday rotation on polarization radiation, we present the polarization intensity (upper row) and polarization angle (lower row) maps at different frequencies in Figure 1.From left to right, images correspond to the frequencies of 10, 5, 1.42, and 1.0 GHz, respectively.As shown in the upper panels of Figure 1, with decreasing frequency, the maximum value of polarization intensity is slightly lower, while the minimum value is significantly lower.At the low-frequency regime (see panel (d)), we can see significantly elongated filament structures due to the enhancement of Faraday depolarization.
The lower panels of Figure 1 present the distribution of polarization angle.At the high frequency of 10 GHz (see panel (e)), it is −0.18 < ϕ < 0.18 radian, where the Faraday rotation effect is almost negligible.As the frequency decreases, we see that Faraday rotation increases the polarization angle.

Finding of the power-law relationship
Based on Sections 2.3 and 4.1, we analyze the polarization intensity variance A and polarization angle dispersion B using each MHD turbulence data listed in Table 1.Then, we use Bayesian analysis to search for the power-law relationship of A and B as a function of M A .Given the differences in the turbulence nature of MHD (see Section 2.1), we perform a separate fitting for the sub-and super-Alfvénic regimes.
The results are plotted in Figure 2 for A vs. M A (left column) and B vs. M A (right column).From top to bottom, the frequencies used are 10, 5, 1.42, and 1.0 GHz, respectively.The blue and green shaded areas represent the highest density interval 94%, i.e., the 94% fitted confidence level, in the sub-and super-Alfvénic regimes, respectively.In the sub-Alfvénic regime, we find the power-law relations of A = 0.05M 1.91 A and B = 0.15M 0.94 A , while in the super-Alfvénic regime, we obtain the relations of A = 0.06M 1.88  A and B = 0.15M 1.18 A .These relationships have been transformed from a logarithmic space to a linear space.In the range from 10 to 1.42 GHz, we find that the power-law relation remains the same, except for the change in fitting error as summarized in Table 2 due to the Faraday rotation effect.Although the power-law relationship between A and M A changes slightly when the frequency decreases to 1.0 GHz, we focus more on the results at a frequency of 1.42 GHz, so this has no impact on our observational application in Section 5.  Table 1.Simulation data arising from MHD turbulence.The plasma parameter β = 2M 2 A /M 2 s is related to M A and M s .B 0 is the uniform magnetic field set initially, and δB rms is the root mean square of the fluctuation field δB as a consequence of MHD turbulence evolution.Models A1 to A10 are simulations with a numerical resolution of 792, while models B1 to B9 are simulations of 480.

Figure 1. Maps of polarization intensities (upper row) and polarization angle distributions (lower row)
. From left to right, the images correspond to the frequencies of 10, 5, 1.42, and 1.0 GHz, respectively.Polarization intensity is normalized in units of the mean synchrotron intensity while polarization angle is in units of a radian.The calculation is based on A3 listed in Table 1.
As a result, since the polarization intensity variance and polarization angle dispersion exhibit a good power-law relationship with M A , we can employ these relations to determine the magnetization M A of the ISM.

Testing for the measurement techniques
To test the reliability of the measurements, we first apply the two techniques to the imaging data we obtained by MHD turbulence simulations.By knowing the values of M A in advance, we can compare their values with the M A measured by techniques.Figure 3 shows the results we obtained by polarization intensity variance (upper panel) and the polarization angle dispersion (lower panel).The error bars are calculated using the upper and lower error limits summarized in Table 2, which also indicates the range of fitting errors presented in Figure 2. From Figure 3, we can see that filled data points for circles and triangles are distributed near the dashed line with a slope equal to 1, which demonstrates that the M A measured by both methods is consistent with the realistic M A determined from 3D MHD turbulence simulations.
Consequently, the power-law relationships we obtained are reliable for measuring M A .
Next, we apply the power-law relationships summarized in Table 2 to measure the magnetization of the Galactic ISM.

APPLICATION TO OBSERVATIONS
Here, we apply our techniques to the Canadian Galactic Plane Survey (CGPS) data with arcminute-scale images of all major components of the ISM over a large portion of the Galactic disk.Note that the Dominion Radio Astrophysical Observatory (DRAO) Synthesis Telescope surveys have imaged a 73 section of the Galactic plane.The DRAO fields were combined into 1024 × 1024 pixel mosaic images (5.• 12 × 5. • 12 for each), with a resolution close to 1 ′ on a Galactic Cartesian grid at the 1.42 GHz.By randomly selecting 9 mosaics named MA1, MF2, MST2, MU1, MU2, MW1, MEV1, MEV2, and MEX1 (see Taylor et al. 2003 for detailed coordinate locations), we extract 800 × 800 pixels, corresponding to the region of 4. • 0 × 4. • 0, from individual mosaic images to avoid the margin effect of the images (see also Zhang et al. 2019).We consider a statistical average of different sub-regions from the whole mosaic image to measure expected M A from the whole mosaic image.In practice, we first decompose the selected 4. • 0 × 4. • 0 mosaic image into 400 and 256 sub-block regions, which correspond to 0. • 20 × 0. • 20 and 0. • 25 × 0. • 25 pixels, respectively.Then, we employ the power-law relationships summarized in Table 2 to obtain the M A values of each sub-block region.Finally, we obtain statistical measurements of M A for all sub-regions to characterize the magnetization of the whole mosaic image.
The M A measurement results using 0. The upper panels of Figure 5 show the comparison of the average M A values measured by methods A and B when selecting different sub-region sizes.As seen in the upper panels, the results measured by the two methods are relatively similar in most regions, except that the results of method B are about 0.1 larger for the MEV1 and MEX1 mosaics.This result shows that on the one hand, the measurement results of method A are more reliable, and on the other hand, the size of the selected sub-region has no effect on measurements featured by average values.The lower panels of Figure 5 compare the M A distribution peak values of the same method when changing the size of the sub-region sizes, with the error bar plotted by a Gaussian function fitting error.It can be seen that the measurement results of methods A and B are the same in different sub-region sizes, indicating that the choice of sub-region size has little effect on the measurement results of M A .As a result, we propose that the low-latitude Galactic ISM is mainly dominated by sub-Alfvénic turbulence, with a magnetization level approximately between 0.5 and 1.0.

DISCUSSION
Our current work has explored how polarization intensity variance can be used to measure the ISM magnetization, confirming the feasibility of this method.Considering the influence of Faraday depolarization on the measurement method, we find that in the explored parameter space, even if the polarization angle rotation reaches 2π radians (see Figure 1), the power-law relation we established is less affected by the Faraday rotation effect, that is, only the fitting error variation of Bayesian statistics (see Table 2).
To further confirm the complementarity of the polarization intensity variance measurements, we have also tested other methods such as the top-base ratio from polarization gradient statistics, inverted variance method, and polarization degree dispersion.These methods have been explored in Lazarian et al. (2018) and Carmo et al. (2020).Our testing results are similar to Figure 7 of Lazarian et al. (2018) using the Galactic Arecibo L-band Feed Array HI data.Note that the measurement of the polarization degree dispersion often provides a smaller measurement value, because the Faraday rotation effect will result in a relatively low observed polarization degree.We would like to emphasize that our current work is based on the framework of the modern understanding of MHD turbulence theory.The established power-law relationship is confined to the strong turbulence regime of the MHD turbulence cascade (see Section 2.1), that is, we ensure that the MHD simulation dataset listed in Table 1 has a sufficiently wide inertial range, with l trans and l A greater than l diss (see Equations ( 2) and ( 5 A for compressible turbulence.What we would like to clarify is that their polarization dispersion is defined as δϕ = δB ⊥ /B, which is different from our definition in Equation ( 13).Note that Falceta-Gonc ¸alves et al. ( 2008) mentioned that the polarization angle dispersion may be related to the angle between the LOS direction and the mean magnetic field.In our work, we do not consider the influence of the angle between the LOS and the mean magnetic field on our new method; this will be explored in the future.
In general, when the brightness temperature approaches the thermal temperature of relativistic electrons, the synchrotron self-absorption effect is expected to be important at low frequencies (Zhang & Wang 2022).However, our testing shows that this radiative transfer process is negligible for the parameters used in this paper.In addition, it is not important for the effects of nonadiabatic propagation for Faraday rotation in the Galactic ISM (Lee et al. 2019).

SUMMARY
In this work, we proposed a new method, polarization intensity variance, for measuring ISM magnetization.After confirming the reliability of the new technique using numerical simulation data, we applied the technique to measure the Galactic ISM magnetization with real observations.The main results are summarized as follows.
• We proposed that polarization intensity variance can achieve accurate measurement of the ISM magnetization.Compared with the polarization angle dispersion method, we verify the reliability of the polarization intensity variance method in measuring the ISM magnetization.
• We found a power-law relationship of A ∝ M 2 A between polarization intensity variance A and magnetization M A .
For the polarization angle dispersion, the power-law relationships are B ∝ M 0.9 A in the sub-Aflvénic turbulence and B ∝ M 1.2 A in the super-Aflvénic turbulence.
• Using simulation polarization data with the known M A values, we demonstrated that polarization intensity variance can successfully recover the underlying magnetization M A , supporting the applicability of this new method.
• With the application of polarization intensity variance, we propose that the low-latitude Galactic ISM is dominated by sub-Alfénic turbulence, with M A approximately between 0.5 and 1.0.
We would like to thank the anonymous referee for the valuable comments that improved our manuscript.

Figure 2 .
Figure 2. The polarization intensity variance A (left column) and polarization angle dispersion B (right column) as a function of Alfvénic Mach number M A .From upper to lower, the panels correspond to 10, 5.0, 1.42, and 1.0 GHz, respectively.
• 20 × 0. • 20 pixels are plotted in Figures 4, from which we can see that histogram distributions from methods A and B are approximately consistent.To obtain quantitative M A values, we characterize the measurement for each mosaic image by the average and peak values of histogram distributions.Fitting the M A distribution histogram with the Gaussian function at a 95% confidence level, we obtained the average and peak values of the M A distribution for each mosaic image.

Figure 3 .
Figure 3.Comparison between the measured M A and the realistic M A .The upper and lower panels are from the polarization intensity variance and polarization angle dispersion, respectively.The errors are calculated by formulae listed in Table2.
)).For polarization angle dispersion, our results arising from compressible turbulence show the power-law relation of B ∝ M 0.9 A and B ∝ M 1.2 A in the sub-and super-Alfvénic regimes, respectively.Differently, Skalidis et al. (2021) and Lazarian et al. (2023) claimed that the polarization angle dispersion maintains the relation of B ∝ M A for incompressible turbulence, and B ∝ M 2

Figure 4 .
Figure 4.The measured M A distribution of CGPS mosaic images with the size of 0. • 20 × 0. • 20 pixels.The red and green histograms correspond to the distribution of M A calculated by the polarization intensity variance and polarization angle dispersion, respectively.
J.F.Z.thanks for the support from the National Natural Science Foundation of China (grant No. 11973035), the Hunan Natural Science Foundation for Distinguished Young Scholars (No. 2023JJ10039), and the China Scholarship Council for the overseas research fund.H.P.X.thanks the support from the Scientific Research Foundation of the Education Bureau of Hunan Province (No. 23A0132).X.J.Y is supported in part by NSFC 12333005 and 12122302 and CMS-CSST-2021-A09.

Table 2 .
The power-law relationship of polarization intensity variance A and polarization angle dispersion B vs. M A at three different frequencies.