Measurement of Interstellar Magnetization by Synchrotron Polarization Variance
Ning-Ning Guo Jian-Fu Zhang Hua-Ping Xiao Jungyeon Cho Xue-Juan Yang





Department of Physics, Xiangtan University, Xiangtan, Hunan 411105, People’s Republic of China;



Department of Physics, Xiangtan University, Xiangtan, Hunan 411105, People’s Republic of China;

Key Laboratory of Stars and Interstellar Medium, Xiangtan University, Xiangtan 411105, People’s Republic of China

Department of Astronomy and Space Science, Chungnam National University, Daejeon, Republic of Korea



Department of Physics, Xiangtan University, Xiangtan, Hunan 411105, People’s Republic of China;

Key Laboratory of Stars and Interstellar Medium, Xiangtan University, Xiangtan 411105, People’s Republic of China



Department of Astronomy and Space Science, Chungnam National University, Daejeon, Republic of Korea



Department of Physics, Xiangtan University, Xiangtan, Hunan 411105, People’s Republic of China;

Key Laboratory of Stars and Interstellar Medium, Xiangtan University, Xiangtan 411105, People’s Republic of China


jfzhang@xtu.edu.cn (JFZ), hpxiao@xtu.edu.cn (HPX), jcho@cnu.ac.kr(JC), xjyang@xtu.edu.cn (XJY)




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 MHD turbulence. With the finding of the power-law relation of 𝒜∝ M_ A^2 between polarization intensity variance 𝒜 and Alfvénic Mach number M_ A, we apply the new technique to the Canadian Galactic Plane Survey (CGPS) data, achieving 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 (; ). 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., ). 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 (). The traditional Faraday rotation synthesis can well present the magnetic structure of the Milky Way and nearby galaxies. However, this method cannot determine the magnetic field direction due to the ambiguity of the phase angle (). The Davis-Chandrasekhar-Fermi (DCF, ) 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 (, henceforth GS95; ; ).

Recently, different improved versions of the DCF method have been enhanced. For instance,
<cit.> 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 (), atomic arrangement effect (), and polarization angle difference as a function of the distance between a pair of measurement points (). 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 (). 

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 (), and the velocity centroid and rotation measurement (). Second, gradient measurements involve velocity () and synchrotron polarization (). Third, the ratio of the E and B modes arises from dust polarization ().

In the earlier numerical work, <cit.> proposed the synchrotron polarization dispersion for measuring the scaling index of underlying MHD turbulence (see also  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.[When mentioning the magnetization or magnetization levels in this paper, we specifically refer to it as the Alfvénic Mach number M_ A.] 

The structure of this paper is regularized as follows. In Section <ref>, we provide theoretical fundamentals including the properties of MHD turbulence, synchrotron polarization radiation, and statistical measure techniques. Section <ref> describes simulation methods with basic parameter setup. Numerical results and application of observational data are presented in Sections <ref> and <ref>, respectively. Sections <ref> and <ref> give discussion and summary, respectively.



§ THEORETICAL CONSIDERATIONS
 



 §.§ 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

    v_l/l_⊥=V_ A/l_∥

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_l^3/l remains constant, i.e., v_l∼ l^1/3, which results in the classical Kolmogorov’s spectrum of E(k)∼ k^-5/3 ().  

In later studies, GS95's work, in relation to trans-Alfvénic turbulence, was generalized to two scenarios: M_ A=V_ L/ V_ A<1 and M_ A >1 (; ), where V_L is the injection velocity at the injection or outer scale L_ inj. For the former, M_ A <1, <cit.> found a transition scale 

    l_ trans=L_ injM_ A^2.

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[In general, at the scale smaller than l_ diss, turbulence dissipation is dominated by numerical dissipation (≲ 10 cells, ). If the transition scale is smaller than the dissipation scale, the numerical results connecting with the astrophysical environment may cause problems ().], turbulence cascade shows the strong turbulence properties, with the parallel and perpendicular scale relationship of

    l_∥≈ L_ inj(l_⊥/L_ inj)^2/3 M_ A^-4/3,

and with velocity and scale relationship 

    v_⊥=V_ A(l_⊥/L_ inj)^1/3M_ A^4/3=V_ L(l_⊥/L_ inj)^1/3M_ A^1/3.

Note that when M_ A=1, Eqs. (<ref>) and (<ref>) 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 () predicted that this super-Alfvénic turbulence has the transition scale of

    l_ A=L_ injM_ A^-3,

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 () 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ν ()

    I(ν) ∝∫_0^L B_⊥^1-αν^αd L^',

where L represents the integration length along the LOS. Based on Eq.(<ref>), we can get the linearly polarized intensity by 

    P(ν) ∝ I ×3-3α/5-3α=3-3α/5-3α∫_0^L B_⊥^1-αν^αd L^'.

The Stokes parameters Q and U are related to the polarization intensity by Q=Pcos2ϕ_0, U=Psin2ϕ_0, where ϕ_0 is the polarization angle. Defining the complex polarization vector

    P=Q+iU,

we can write the polarization intensity as 

    P= √(Q^2+U^2), ϕ_0=1/2arctanU/Q.


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

    ϕ=ϕ_0+ RM  λ^2,

where the rotation measure (RM) is given by 

    RM=0.81 ∫_0^Ln_e B_∥dz   rad  m^-2.

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 <cit.>, we define the correlation function of the polarization intensity as

    𝒜=⟨ P(R_1)P(R_2) ⟩,

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 𝒜 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., ), we define the polarization angle dispersion by

    ℬ=√((ϕ-⟨ϕ⟩)^2).


With the above two definitions (Eqs. (<ref>) and (<ref>)), this paper will explore how 𝒜 and ℬ 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 (). 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∼𝒩 (β 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 (𝒜 and ℬ) as the dependent variables. The assumed distributions and parameters of all priors in this linear model are listed as follows

    Y=𝒩(β M_ A,σ^2),


    β=𝒩(-1, 1),


    σ^2=| Cauchy(0.1,2) |.

We use Markov Chain Monte Carlo methods in Bayesian analysis, after choosing 20, 000 posterior [We have tested the influence of sample numbers on the results. When choosing posterior samples are greater than 10, 000, the numerical results will remain the same.], 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:

    ∂ρ/∂ t + ∇·(ρν)=0,


    ρ[∂ν/∂ t + (ν·∇)ν] + ∇ p - J×B/4π = f,


    ∂B/∂ t -∇×(ν×B)=0,


    ∇·B=0,

where ρ is gas density, t the evolution time of the fluid, p=c_s^2 ρ 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 (). 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_ A^2/M_ s^2, we characterize different data groups as listed in Table <ref>. 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.





§ RESULTS





 §.§ 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. () 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 <ref>, 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 (<ref>) 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 <ref>. 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 <ref>, 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 <ref> 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. When reduced to 1.42 GHz, polarization angle rotates about ϕ≃± 3π  radians (see panel (g)).






 §.§ Finding of the power-law relationship


Based on Sections <ref> and <ref>, we analyze the polarization intensity variance 𝒜 and polarization angle dispersion ℬ using each MHD turbulence data listed in Table <ref>. Then, we use Bayesian analysis to search for the power-law relationship of 𝒜 and ℬ as a function of M_ A. Given the differences in the turbulence nature of MHD (see Section <ref>), we perform a separate fitting for the sub- and super-Alfvénic regimes.

The results are plotted in Figure <ref> for 𝒜 vs. M_ A (left column) and ℬ 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 𝒜=0.05M_ A^1.91 and ℬ=0.15M_ A^0.94, while in the super-Alfvénic regime, we obtain the relations of 𝒜=0.06M_ A^1.88 and ℬ=0.15M_ A^1.18. 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 <ref> due to the Faraday rotation effect. Although the power-law relationship between 𝒜 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 <ref>.

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 <ref> 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 <ref>, which also indicates the range of fitting errors presented in Figure <ref>. From Figure <ref>, 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 <ref> 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×512 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  for detailed coordinate locations), we extract 800× 800 pixels, corresponding to the region of 40 × 40, from individual mosaic images to avoid the margin effect of the images (see also ).

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×40 mosaic image into 400 and 256 sub-block regions, which correspond to 0 20×020 and 0 25×025 pixels, respectively. Then, we employ the power-law relationships summarized in Table <ref> 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 20×020 pixels are plotted in Figures <ref>, from which we can see that histogram distributions from methods 𝒜 and ℬ 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. 

The upper panels of Figure <ref> show the comparison of the average M_ A values measured by methods 𝒜 and ℬ 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 ℬ are about 0.1 larger for the MEV1 and MEX1 mosaics. This result shows that on the one hand, the measurement results of method 𝒜 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 <ref> 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 𝒜 and ℬ 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 <ref>), 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 <ref>).  

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 <cit.> and <cit.>. Our testing results are similar to Figure 7 of <cit.> 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 <ref>), that is, we ensure that the MHD simulation dataset listed in Table <ref> has a sufficiently wide inertial range, with l_ trans and l_ A greater than l_ diss (see Equations (<ref>) and (<ref>)). For polarization angle dispersion, our results arising from compressible turbulence show the power-law relation of ℬ∝ M_ A^0.9 and ℬ∝ M_ A^1.2 in the sub- and super-Alfvénic regimes, respectively. Differently, <cit.> and <cit.> claimed that the polarization angle dispersion maintains the relation of ℬ∝ M_ A for incompressible turbulence, and ℬ∝ M_ A^2 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 (<ref>). Note that <cit.> 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 <cit.>. 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 <cit.>.



§ 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 𝒜∝ M_ A^2 between polarization intensity variance 𝒜 and magnetization M_ A. For the polarization angle dispersion, the power-law relationships are ℬ∝ M_ A^0.9 in the sub-Aflvénic turbulence and ℬ∝ M_ A^1.2 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. 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.




aasjournal

