hammurabi X: Simulating Galactic Synchrotron Emission with Random Magnetic Fields

The hammurabi code (Waelkens et al. 2009) is an astrophysical simulator based on 3D models of the components of the magnetized ISM such as magnetic fields, thermal electrons, relativistic electrons, and dust grains. It performs an efficient LoS integral through the simulated Galaxy model using a HEALPix-based¹¹ (Górski et al. 2005) nested grid to produce observables such as the Faraday rotation measure and diffuse synchrotron and thermal dust emission¹² in full Stokes I, Q, and U, while taking into account beam and depth depolarization as well as Faraday effects. The updated version, hammurabi X (Wang et al. 2019),¹³ has been developed to achieve higher computing performance and precision. Version 2.3.0 is available at doi:10.5281/zenodo.3522599. Specific effort has been devoted to the parallel computation of the LoS integral and vector-field FFT.

hammurabi X currently uses the HEALPix library (Górski et al. 2005) for observable production, where the LoS integral accumulates through several layers of spherical shells with adaptable HEALPix resolutions. We provide two modes of integral shell arrangements. In the auto-shell mode, given R as the maximum simulation radius, the nth shell out of N total shells covers the radial distance from 2^{(n − N−1)}R to 2^(n − N)R, except for the first shell which starts at the observer. The nth shell is by default set up with the HEALPix resolution-controlling parameter N_side = 2^(n − 1)N_min,¹⁴ where N_min represents the lowest simulation resolution at the first shell. Alternatively in the manual-shell mode, shells are defined explicitly by a series of dividing radii and HEALPix N_side's. The radial resolution along the LoS integral is uniformly set by the minimal radial distance for each shell. The auto-shell mode follows the idea that the integral domain is discretized with elemental bins of the same volume, while the manual-shell mode allows users to refine specific regions to meet special realization requirements.

The LoS integral is carried out hierarchically: at the top level, the integral is divided into multiple shells with given spherical resolution settings, while at the bottom level inside each shell (where the spherical resolution is fixed), the radial integral is carried out with the midpoint rule for each radial bin. The accumulation of observable information from the inner to outer shells is applied at the top level. We emphasize that in hammurabi X, the simulation spherical resolution for each shell can be independent of that in the outputs, which means that we can simulate with an arbitrary number of shells and assign each shell a unique N_side value. During each step of the shell-accumulating process, we interpolate (with the linear interpolation provided by HEALPix library) the current result into the output resolution. Consequently, such interpolation between different angular resolutions will inevitably create a certain level of precision loss.

Previously in hammurabi, the generation of the anisotropic component of the random field as well as the modulation of the field strength following various parametric forms led to artificial magnetic field divergence. Now we propose two improved solutions for simulating the random magnetic field. On Galactic scales, a triple Fourier transform scheme is proposed to restore the divergence-free condition via a cleaning process. This imposes the divergence-free property in the random magnetic field (unlike in Planck Collaboration et al. 2016b), which will be discussed in detail in Section 3.3, with its observational implication in Section 4. Alternatively, in a given local region¹⁵ , a vector-field decomposition scheme is capable of simulating more detailed random-field power spectra.

FFTs are necessary for translating the power spectra of random fields into discrete magnetic field realizations on 3D spatial grids. Random-field generators in hammurabi X currently use the FFTW library.¹⁶ The detailed implementation will be discussed in Section 3. In cases where the field is input from an external or internal discrete grid, e.g., a random GMF, the LoS integral at a given position does a linear interpolation (in each phase-space dimension) from nearby grid points. The interpolation algorithm has been calibrated, so the high-resolution outputs are no longer contaminated by the numerical flaws in earlier versions of hammurabi. As illustrated in Figure 1, the interpolation process in the earlier version of hammurabi did not properly calculate the volume of elemental discretization, which resulted, for example, in negative values of the simulated dispersion measure and incorrect small-scale features in comparison to the corrected method in hammurabi X. In this new version, unit tests for linear interpolation can be found in the public repository.

Zoom In Zoom Out Reset image size

Generally speaking, the precision of the linear interpolation (and the corresponding discretization) can in principle be characterized by the goodness of the approximation. This is explicitly affected by the discretization resolution and the arrangement of the sampling/supporting points, and also by the smoothness (as measured by the inverse of the second-order derivative) of the approximation target. In hammurabi X, the interpolation affects the precision in realizing the power spectrum of the random magnetic field generation. This can be improved by increasing the sampling resolution. Furthermore, linear interpolation does not preserve the divergence, but the precision can be improved either by increasing the sampling resolution¹⁷ or by matching the elemental discretization volume in the LoS integral with that in the field generation (as discussed by Waelkens et al. 2009).

Profiling¹⁸ the numerical precision in producing observables is critical in guiding practical applications. A standard hammurabi X simulation routine consists of two major building blocks. The first part is the numerical implementation of specific physical processes like synchrotron emission and Faraday rotation, and the second part is the LoS integral that is universal to all observables. In the following integrated precision check, the correctness of both will be verified and profiled together.

A given magnetic field vector ${\boldsymbol{B}}$ can be decomposed into directions parallel (horizontal) and perpendicular (vertical/poloidal) to the Galactic disk, or to be specific, the $\{\hat{{\boldsymbol{x}}},\hat{{\boldsymbol{y}}}\}$ plane (with $\hat{{\boldsymbol{y}}}$ pointing toward Galactic longitude l = 90°) in the hammurabi X convention, i.e., ${\boldsymbol{B}}$_∥ and ${\boldsymbol{B}}$_⊥ at a given Galactic longitude–latitude position $\{l,b\}$. The LoS direction $\hat{{\boldsymbol{n}}}$ from the observer to the target field position reads

$$\begin{eqnarray}&&\hat{{\boldsymbol{n}}}=\cos (b)\cos (l)\hat{{\boldsymbol{x}}}+\cos (b)\sin (l)\hat{{\boldsymbol{y}}}+\sin (b)\hat{{\boldsymbol{z}}},\end{eqnarray} \tag{ 1 }$$

where $\hat{{\boldsymbol{x}}}$ is conventionally pointing from the observer to the Galactic center. In the same observer-centric Cartesian frame, we can explicitly write down two field components as

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{\parallel }={B}_{\parallel }(\cos ({l}_{0})\hat{{\boldsymbol{x}}}+\sin ({l}_{0})\hat{{\boldsymbol{y}}}),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{\perp }={B}_{\perp }\hat{{\boldsymbol{z}}},\end{eqnarray} \tag{ 3 }$$

where l₀ represents the projected direction of ${\boldsymbol{B}}$ in the $\{\hat{{\boldsymbol{x}}},\hat{{\boldsymbol{y}}}\}$ plane. Then, it is straightforward to calculate two key quantities for the calculation of synchrotron emissivity and Faraday rotation, respectively,

$$\begin{eqnarray}&&| {\boldsymbol{B}}\times \hat{{\boldsymbol{n}}}| =\sqrt{{B}_{\parallel }^{2}+{B}_{\perp }^{2}-| {\boldsymbol{B}}\cdot \hat{{\boldsymbol{n}}}{| }^{2}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\boldsymbol{B}}\cdot \hat{{\boldsymbol{n}}}={B}_{\parallel }\cos (b)\cos (l-{l}_{0})+{B}_{\perp }\sin (b).\end{eqnarray} \tag{ 5 }$$

It is obvious that Faraday rotation is more sensitive to ${\boldsymbol{B}}$_∥ at low Galactic latitudes, and to ${\boldsymbol{B}}$_⊥ at high latitudes. On the contrary, synchrotron emissivity, which is proportional to some power of $| {\boldsymbol{B}}\times \hat{{\boldsymbol{n}}}|$, is more sensitive to ${\boldsymbol{B}}$_⊥ at low Galactic latitudes and to ${\boldsymbol{B}}$_∥ at high latitudes.

Precision checks require a baseline model for each field, from which analytic descriptions of the observables can be explicitly derived. Here we assume spatially homogeneous distributions for the CR electrons (CREs), TEs, and GMF within a given radial distance to the observer. The spectral index of the CRE energy distribution is assumed to be a constant, and consequently, the CRE density N(γ) is described by

$$\begin{eqnarray}&&N(\gamma )={N}_{0}{\gamma }^{-\alpha },\end{eqnarray} \tag{ 6 }$$

where γ represents the CRE Lorentz factor and α represents the constant spectral index of CRE. With the assumed homogeneity in all fields, we can calculate the intrinsic synchrotron total intensity I₀ and polarization Stokes parameter Q₀ and U₀ (in the IAU convention¹⁹ ) before applying the Faraday rotation (Rybicki & Lightman 1979):

$$\begin{eqnarray}&&{I}_{0}={J}_{i}{R}_{0},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{Q}_{0}={J}_{\mathrm{pi}}{R}_{0}\cos (2{\chi }_{0}),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{U}_{0}={J}_{\mathrm{pi}}{R}_{0}\sin (2{\chi }_{0}),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}\begin{array}{rcl}{J}_{i} & = & \displaystyle \frac{\sqrt{3}{e}^{3}| {\boldsymbol{B}}\times \hat{{\boldsymbol{n}}}| {N}_{0}}{4\pi {m}_{{\rm{e}}}{c}^{2}(\alpha +1)}{\left(\displaystyle \frac{2\pi \nu {m}_{{\rm{e}}}c}{3e| {\boldsymbol{B}}\times \hat{{\boldsymbol{n}}}| }\right)}^{\tfrac{1-\alpha }{2}}\\ & & \times {\rm{\Gamma }}\left(\displaystyle \frac{\alpha }{4}+\displaystyle \frac{19}{12}\right){\rm{\Gamma }}\left(\displaystyle \frac{\alpha }{4}-\displaystyle \frac{1}{12}\right),\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\begin{array}{rcl}{J}_{\mathrm{pi}} & = & \displaystyle \frac{\sqrt{3}{e}^{3}| {\boldsymbol{B}}\times \hat{{\boldsymbol{n}}}| {N}_{0}}{16\pi {m}_{{\rm{e}}}{c}^{2}}{\left(\displaystyle \frac{2\pi \nu {m}_{{\rm{e}}}c}{3e| {\boldsymbol{B}}\times \hat{{\boldsymbol{n}}}| }\right)}^{\tfrac{1-\alpha }{2}}\\ & & \times {\rm{\Gamma }}\left(\displaystyle \frac{\alpha }{4}+\displaystyle \frac{7}{12}\right){\rm{\Gamma }}\left(\displaystyle \frac{\alpha }{4}-\displaystyle \frac{1}{12}\right),\end{array}\end{eqnarray} \tag{ 11 }$$

where e is the electron charge, and m_e is the electron mass, R₀ is the spherical LoS integral depth, and ν is the observational frequency. The intrinsic polarization angle χ₀ can be derived from

$$\begin{eqnarray}&&\tan ({\chi }_{0})=\displaystyle \frac{{B}_{\perp }\cos (b)-{B}_{\parallel }\sin (b)\cos (l-{l}_{0})}{{B}_{\parallel }\sin (l-{l}_{0})},\end{eqnarray} \tag{ 12 }$$

as illustrated in Figure 2. With the same modeling, the Faraday depth ϕ can be described by

$$\begin{eqnarray}&&\phi (l,b)={\phi }_{0}{R}_{0},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\phi }_{0}=-{N}_{{\rm{e}}}({\boldsymbol{B}}\cdot \hat{{\boldsymbol{n}}})\left(\displaystyle \frac{{e}^{3}}{2\pi {m}_{{\rm{e}}}^{2}{c}^{4}}\right),\end{eqnarray} \tag{ 14 }$$

where N_e represents the constant homogeneous TE density assumed within spherical radius R₀. In the end, the observed synchrotron polarization Stokes parameters Q and U reflect the Faraday rotation as

$$\begin{eqnarray}&&Q+{iU}=({Q}_{0}+{{iU}}_{0}){\int }_{0}^{{R}_{0}}\displaystyle \frac{{e}^{2i{\phi }_{0}{\lambda }^{2}r}}{{R}_{0}}{dr},\end{eqnarray} \tag{ 15 }$$

which indicates that the polarized intensity receives a correction factor $| \sin (\phi {\lambda }^{2})/(\phi {\lambda }^{2})|$ known as the Faraday depolarization. The formulae above analytically derive calculable results for reference in verifying the numerical outputs. In real applications, the magnetic field and CRE spectral index are not constant, and the methods used by hammurabi X for calculating synchrotron emissivity and Faraday rotation can be more generic, as presented in Appendix A.

Zoom In Zoom Out Reset image size

Figure 3 presents the absolute and relative numeric error distributions of the synchrotron total intensity from a single LoS integral shell. For an observable X, the absolute error is defined as the difference between the simulated output X_sim and the analytic reference X_ref as (X_sim − X_ref), while the relative error is defined by 2(X_sim − X_ref)/(X_sim + X_ref). The Faraday depth calculator shares a similar error distribution to the calculator of synchrotron total intensity. Meanwhile, Figure 4 presents the absolute and relative numeric error distributions of synchrotron Stokes Q also from a single LoS integral shell, which serves as an example for illustrating the numeric precision in calculating tensor fields. With constant field models in testing, the numeric errors are mainly induced by the integration and interpolation methods and therefore independent of the LoS resolution. Even with simple field settings, we can observe a few percent relative error appearing in Figure 4. Considering the future use of hammurabi X in inferring Galactic component structures with astrophysical measurements, if the magnitude of such numerical errors is larger than the observational uncertainties, a Bayesian analysis with hammurabi X will consequently suffer from higher uncertainties and bias in parameter estimation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In terms of the multishell arrangement in real application, the output precision is affected by the spherical surface interpolation provided by the HEALPix library. The motivation of allowing different resolution settings along with the divided LoS integral shells is to save computing resources as mentioned in Planck Collaboration et al. (2016b). It is worth noticing that in the simulation, the pixel values are calculated along their central spherical coordinates. This is different from the actual astrophysical measurements where each pixel value is estimated based on many observational hits. And thus, for quickly comparing low-resolution simulation results with high-resolution data, we recommend interpolating data on the sky, accounting for simulations' sample directions, instead of downgrading data by averaging over high-resolution pixels. In this way, we avoid comparing exactly predicted values of simulation to region-averaged values of measurements. Alternatively, a very stringent simulation should be designed to mimic the true observation beams, which is computationally heavy without hammurabi's method. But even with our method, no simulation can capture reality perfectly, and the user must always be careful to test that the simulation resolution is sufficient for probing the observational property in question.

The testing cases displayed above are prepared by assuming a constant magnetic field and TE field distributions. The numerical errors would inevitably grow larger when the input Galactic components have small-scale features near or below the discretization resolution. This issue can be handled efficiently in the future by an adaptively refined mesh/pixelization.

The computationally heavy processes in hammurabi X are the LoS integration for HEALPix map pixels, the random-field generation with FFTs, and the linear interpolation for fields prepared in grids (e.g., internal random fields and other external fields). Massive observable production, HEALPix map distribution, and recycling of physical fields require MPI²⁰ parallelization and therefore are beyond our scope in this report. In this work, multithreading is always essential at the bottom level of parallelism. Figure 5 presents the strong scaling²¹ in observable production with various GMF and TE field combinations. The strong scaling with either computationally heavy (with random-field generation) or light (without random-field generation) pipelines follows the Amdahl law (Amdahl 1967) with around 2% serial remnants. Note that the speedup properties are not very sensitive to the resolution setting in various simulation routines, because the workload of pure numerical operations is proportional to the discretization resolution.

Zoom In Zoom Out Reset image size

The realization of turbulent magnetic field is a major module in hammurabi X, as the correctness of most simulations relies on a physically motivated and accurate description of the turbulent fields in the multiphase ISM. In this section, we present two Gaussian random GMF generators that are by definition divergence-free and capable of realizing field alignment and/or strength modulation on Galactic scales or an anisotropic²² power spectrum on small scales.

There are several criteria that a random GMF generator should satisfy. That it be divergence-free (or solenoidal) is always the prime feature of any magnetic field. Absolute zero divergence is hard to define under discretization, but in principle either a vector-field decomposition or a Gram–Schmidt process in the frequency domain is capable of cleaning field divergence. In realistic cases when a large-scale spatial domain is expected to be filled with random magnetic fields, the field strength and alignment need to be correlated with the large-scale structures in the Galaxy. This requirement complicates the generating process, because the divergence-free property should also be satisfied simultaneously. It is straightforward to generate a divergence-free Gaussian random field, and equally simple to then rescale or stretch it, as done in Jaffe et al. (2010). But the reprofiling process destroys the divergence-free property if it is applied naively.

A triple Fourier transform scheme is thus proposed mainly to reconcile these two requirements. At Galactic scales, the new scheme allows modification of the Gaussian random realization by a given inhomogeneous spatial profile for the field strength. Note that aligning the magnetic field to a given direction is easy to implement in the spatial domain, but locally varying anisotropy in the energy power spectra is not feasible by a single FFT. In studies of Galactic emission from MHD plasma, the dependence of local structure on a varying direction profile breaks the symmetry required for using the FFT. To perform more detailed modeling of the turbulent GMF power spectrum, we provide a local generator ("local" in the sense that the mean field can be approximated in a uniform distribution) with explicit or implicit vector decomposition.

Consider a magnetic field distribution ${\boldsymbol{B}}({\boldsymbol{x}})$ = ${\boldsymbol{B}}$₀(${\boldsymbol{x}}$) + ${\boldsymbol{b}}({\boldsymbol{x}})$ and its counterpart $\tilde{{\boldsymbol{B}}}({\boldsymbol{k}})$ in the frequency domain, where ${\boldsymbol{B}}$₀ and ${\boldsymbol{b}}$ represent the regular and turbulent fields, respectively. The simplest turbulent power spectrum is represented by the trace of the isotropic magnetic field spectrum tensor in scalar form, $P(k)\propto \langle \tilde{{\boldsymbol{B}}}({\boldsymbol{k}})\cdot {\tilde{{\boldsymbol{B}}}}^{* }(k){\rangle }_{{\boldsymbol{B}}}$.²³ This kind of spectrum is widely used as a first approach to the turbulent field realization where the spectral shape is important. In general, we could parameterize the basic scalar spectrum as

$$\begin{eqnarray}\begin{array}{rcl}P(k) & = & \displaystyle \frac{{P}_{0}}{4\pi {k}^{2}}\left[{\left(\displaystyle \frac{{k}_{0}}{{k}_{1}}\right)}^{{\alpha }_{1}}{\left(\displaystyle \frac{k}{{k}_{1}}\right)}^{6}{ \mathcal H }({k}_{1}-k)\right.\\ & & +{\left(\displaystyle \frac{k}{{k}_{0}}\right)}^{-{\alpha }_{1}}{ \mathcal H }(k-{k}_{1}){ \mathcal H }({k}_{0}-k)\\ & & +\left.{\left(\displaystyle \frac{k}{{k}_{0}}\right)}^{-{\alpha }_{0}}{ \mathcal H }(k-{k}_{0})\right],\end{array}\end{eqnarray} \tag{ 16 }$$

where ${ \mathcal H }$ represents the Heaviside step function. The last term in Equation (16) represents the forward magnetic cascading of MHD turbulence from the injection scale k₀ to small scales (k > k₀), while the first two terms describe the inverse cascading (Pouquet et al. 1976) in MHD turbulence from k₀ to scale ${k}_{1}\simeq 1/L$, which corresponds to the physical size L of the MHD system. According to the simulation results from Brandenburg et al. (2019), we set ${k}_{1}=0.1\,{\mathrm{kpc}}^{-1}$ and α₁ = 0.0 by default in this work if not specified. Note that although not explicitly mentioned here, the Nyquist frequency cutoff k_max requires an extra Heaviside factor ${ \mathcal H }({k}_{\max }-k)$ in Equation (16).

In terms of more physical parameterization, we are interested in realizing theoretical descriptions of turbulence in the compressible plasma recently discussed by Cho & Lazarian (2002), Caldwell et al. (2016), and Kandel et al. (2017). In a compressible plasma, turbulence can be decomposed into Alfvén, fast, and slow modes. Two critical plasma status parameters are the ratio β and the Alfvén Mach number M_A. The plasma β is the ratio of gas pressure to magnetic pressure, which represents compressibility of the plasma, with $\beta \to \infty$ indicating the incompressible regime. The Alfvén Mach number is the ratio of the injection velocity to the Alfvén velocity, with M_A > 1.0 representing the super-Alfvénic regime while M_A < 1.0 means sub-Alfvénic turbulence. The general form of the compressible MHD magnetic field spectrum tensor trace reads

$$\begin{eqnarray}&&P(k,\alpha )=\sum _{i}{P}_{i}(k){F}_{i}({M}_{{\rm{A}}},\alpha ){h}_{i}(\beta ,\alpha ),\end{eqnarray} \tag{ 17 }$$

where i = {A, f, s} denotes one of the three MHD modes as described in detail in Section 4.1. In hammurabi X, compressible MHD is only realized by the local generator, thus $\cos (\alpha )=\hat{{\boldsymbol{k}}}\cdot {\hat{{\boldsymbol{B}}}}_{0}$ is adopted, with ${{\boldsymbol{B}}}_{0}$ taken as the regular field near the observer. A detailed application example of F_i and h_i is presented in Section 4. Some additional information can be found in Appendix B for readers who are interested in the technical shortcuts in random-field generation and the sampling precision.

One major task of hammurabi X is to generate a random GMF that can cover a specific scale in the spatial domain. However, an inhomogeneous correlation structure is not diagonal in the frequency domain. In this case, we try to impose an energy density and alignment profile in the spatial domain after the random realization is generated in the frequency domain with an isotropic spectrum. Then, the field divergence can be cleaned back in the frequency domain with the Gram–Schmidt process. The whole procedure of this scheme requires two backward and one forward FFTs.

After a Gaussian random magnetic field is realized in the frequency domain, each grid point holds a vector ${\boldsymbol{b}}$ drawn from an isotropic field dispersion. The key to the triple transform is the large-scale alignment and energy density modulation process. The alignment direction $\hat{{\boldsymbol{H}}}$ at different Galactic positions should be predefined, like the energy density profile. We introduce the alignment parameter ρ for imposing the alignment profile as

$$\begin{eqnarray}&&{\boldsymbol{b}}({\boldsymbol{x}})\to \displaystyle \frac{({{\boldsymbol{b}}}_{\parallel }\rho +{{\boldsymbol{b}}}_{\perp }/\rho )}{\sqrt{\tfrac{1}{3}{\rho }^{2}+\tfrac{2}{3}{\rho }^{-2}}},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{\boldsymbol{b}}}_{\parallel }=\displaystyle \frac{({\boldsymbol{b}}\cdot \hat{{\boldsymbol{H}}})}{| \hat{{\boldsymbol{H}}}{| }^{2}}\hat{{\boldsymbol{H}}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{\boldsymbol{b}}}_{\perp }=\displaystyle \frac{\hat{{\boldsymbol{H}}}\times ({\boldsymbol{b}}\times \hat{{\boldsymbol{H}}})}{| \hat{{\boldsymbol{H}}}{| }^{2}}.\end{eqnarray} \tag{ 20 }$$

ρ = 1.0 means no preferred alignment direction, while $\rho \to 0$ ($\rho \to \infty$) indicates extremely perpendicular (parallel) alignment with respect to $\hat{{\boldsymbol{H}}}$. (Previously, the alignment operation in hammurabi was carried out by regulating ${\boldsymbol{b}}$_∥ only (Jaffe et al. 2010), which is phenomenologically equivalent to our approach presented here.) Note that ρ and $\hat{{\boldsymbol{H}}}$ can either be defined as a global constant or as a function of other physical quantities such as the regular magnetic field and the Galactic ISM structure (a detailed description can be found in the hammurabi X wiki page).

For regulating the field energy density, a simple example with an exponential scaling profile (which can be customized in future studies) is proposed as

$$\begin{eqnarray}&&S({\boldsymbol{x}})=\exp \left(\displaystyle \frac{{R}_{\odot }-r}{{h}_{r}}\right)\exp \left(\displaystyle \frac{| {z}_{\odot }| -| z| }{{h}_{z}}\right),\end{eqnarray} \tag{ 21 }$$

where (r, z) is the coordinate in the Galactic cylindrical frame, and (R_⊙, z_⊙) represents the solar position in the Galactic cylindrical frame. The energy density modulation acts on the vector-field amplitude through

$$\begin{eqnarray}&&{\boldsymbol{b}}({\boldsymbol{x}})\to {\boldsymbol{b}}({\boldsymbol{x}})\sqrt{S({\boldsymbol{x}})}.\end{eqnarray} \tag{ 22 }$$

The above operations of reorienting, stretching, and squeezing magnetic field vectors in the spatial domain do not promise a divergence-free result. To clean the divergence, we transform the reprofiled field forward into the frequency domain and apply the Gram–Schmidt process:

$$\begin{eqnarray}&&\tilde{{\boldsymbol{b}}}\to \sqrt{3}\left(\tilde{{\boldsymbol{b}}}-\displaystyle \frac{({\boldsymbol{k}}\cdot \tilde{{\boldsymbol{b}}}){\boldsymbol{k}}}{| k{| }^{2}}\right),\end{eqnarray} \tag{ 23 }$$

where $\tilde{{\boldsymbol{b}}}$ indicates the frequency-domain complex vector. The coefficient $\sqrt{3}$ is for preserving the spectral power statistically. A second backward Fourier transform is then carried out to provide the final random GMF vector distribution in the spatial domain.

Note that separating the divergence-cleaning process from spatial reprofiling comes with a cost. Strong alignment with ρ ≪ 1 or ρ ≫ 1 are not realizable because the Gram–Schmidt process reestablishes some extra spatial isotropy according to Equation (23). Figure 6 presents typical results of the global random generator in the form of magnetic field probability density distributions, where we assume a Kolmogorov power spectrum. The distributions of b_y and b_z are expected to be identical, with the imposed alignment direction being $\hat{{\boldsymbol{H}}}=\hat{{\boldsymbol{x}}}$. Note that the global generator is designed to realize the inhomogeneity and anisotropy in both spatial and frequency domains, which we then have to process with divergence cleaning to provide conceptually acceptable realizations.

Zoom In Zoom Out Reset image size

A local generator is proposed to realize random GMFs in small-scale regions, like the solar neighborhood, where the regular field can be approximated as homogeneous with a uniform direction, or more precisely speaking, where the random magnetic field two-point correlation tensor can be approximated to be independent of the spatial position. With this assumption, random fields can be realized with a single FFT. Here we describe the vector decomposition method for realizing a Gaussian random magnetic field with a generic anisotropic power spectrum tensor P_ij(${\boldsymbol{k}}$, α), where α represents extra parameters in addition to the wave vector. By assuming Gaussianity, the power spectrum tensor reads

$$\begin{eqnarray}&&{P}_{{ij}}({\boldsymbol{k}},\alpha ){\delta }^{3}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} )=\langle {\tilde{b}}_{i}({\boldsymbol{k}}){\tilde{b}}_{j}^{* }({\boldsymbol{k}}^{\prime} ){\rangle }_{\tilde{{\boldsymbol{b}}}},\end{eqnarray} \tag{ 24 }$$

where $\tilde{{\boldsymbol{b}}}$ represents the complex magnetic field vectors in the frequency domain. Depending on the specific form of the given power spectrum tensor, the vector-field decomposition can be either explicit or implicit.

The implicit vector decomposition sets up two modes (vector bases) for a complex Fourier vector $\tilde{{\boldsymbol{b}}}$, which means

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{b}}}}^{\pm }({\boldsymbol{k}})={\tilde{b}}^{\pm }({\boldsymbol{k}}){\hat{{\boldsymbol{e}}}}^{\pm },\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\hat{{\boldsymbol{e}}}}^{\pm }=\displaystyle \frac{{\hat{{\boldsymbol{e}}}}_{1}\pm i{\hat{{\boldsymbol{e}}}}_{2}}{\sqrt{2}},\end{eqnarray} \tag{ 26 }$$

where the two orthogonal basis vectors ${\hat{{\boldsymbol{e}}}}^{\pm }$ bind with the complex scalar ${\tilde{b}}^{\pm }$ respectively. The vectors $\{{\hat{{\boldsymbol{e}}}}_{1},{\hat{{\boldsymbol{e}}}}_{2},{\hat{{\boldsymbol{e}}}}_{3}\}$ form a Cartesian frame, and to ensure the divergence-free property of the resulting fields, we choose ${\hat{{\boldsymbol{e}}}}_{3}=\hat{{\boldsymbol{k}}}$. During the Fourier transform of $\tilde{{\boldsymbol{b}}}({\boldsymbol{k}})$ into the spatial domain, we have to consider an orthogonal base aligned with the Cartesian grid of ${\boldsymbol{b}}({\boldsymbol{x}})$, and here we adopt one convenient base representation as

$$\begin{eqnarray}&&\hat{{\boldsymbol{k}}}=\left(\displaystyle \frac{{k}_{x}}{k},\displaystyle \frac{{k}_{y}}{k},\displaystyle \frac{{k}_{z}}{k}\right),\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\hat{{\boldsymbol{e}}}}^{-}=\left(\displaystyle \frac{-{k}_{y}}{\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}},\displaystyle \frac{{k}_{x}}{\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}},0\right),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\hat{{\boldsymbol{e}}}}^{+}=\left(\displaystyle \frac{{k}_{x}{k}_{z}}{k\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}},\displaystyle \frac{{k}_{y}{k}_{z}}{k\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}},\displaystyle \frac{-({k}_{x}^{2}+{k}_{y}^{2})}{k}\right),\end{eqnarray} \tag{ 29 }$$

where $k=\sqrt{{k}_{x}^{2}+{k}_{y}^{2}+{k}_{z}^{2}}$. Then, we can proceed by projecting the complex field amplitude into this spatial frame:

$$\begin{eqnarray}&&\tilde{{\boldsymbol{b}}}\cdot \hat{{\boldsymbol{x}}}={\tilde{b}}^{+}({\hat{{\boldsymbol{e}}}}^{+}\cdot \hat{{\boldsymbol{x}}})+{\tilde{b}}^{-}({\hat{{\boldsymbol{e}}}}^{-}\cdot \hat{{\boldsymbol{x}}}),\end{eqnarray} \tag{ 30 }$$

where $\hat{{\boldsymbol{x}}}$ represents the spatial Cartesian coordinate. Implicit decomposition is irrelevant to the choice of $\{{{\boldsymbol{e}}}^{+},{{\boldsymbol{e}}}^{-}\}$ base and useful in the case where only the spectrum trace Tr[P_ij(${\boldsymbol{k}}$)] (over the i, j indices) is given. The amplitude of ${\tilde{b}}^{\pm }$ can be drawn from Gaussian distributions with zero mean and variances ${\sigma }_{\pm }^{2}$, which satisfy

$$\begin{eqnarray}&&{\sigma }_{+}^{2}+{\sigma }_{-}^{2}=\mathrm{Tr}[{P}_{{ij}}({\boldsymbol{k}})]{d}^{3}k,\end{eqnarray} \tag{ 31 }$$

with ${d}^{3}k$ representing the frequency-domain discretization resolution. Equation (31) indicates that the field amplitudes ${\tilde{b}}^{\pm }$ should have a joint power spectrum equal to the trace of the total power spectrum.

The explicit decomposition should be used when the power spectrum tensor is available along with the explicitly defined base $\{{{\boldsymbol{e}}}^{+},{{\boldsymbol{e}}}^{-}\}$, where

$$\begin{eqnarray}&&{\sigma }_{\pm }^{2}={P}^{\pm }({\boldsymbol{k}}){d}^{3}k.\end{eqnarray} \tag{ 32 }$$

A practical example is realizing the Alfvén, fast, and slow modes of an MHD turbulent magnetic field in a compressible plasma. Given a local regular GMF field ${\boldsymbol{B}}$₀, an Alfvén wave propagates along ${\hat{{\boldsymbol{B}}}}_{0}$ with magnetic turbulence in direction ${{\boldsymbol{e}}}^{+}=\hat{{\boldsymbol{k}}}\times {\hat{{\boldsymbol{B}}}}_{0}$, while slow and fast waves generate magnetic turbulence in the direction ${{\boldsymbol{e}}}^{-}={{\boldsymbol{e}}}^{+}\times \hat{{\boldsymbol{k}}}$. A detailed parameterization of compressible MHD turbulent power spectra will be introduced in Section 4 following the corresponding references therein. Note that when the wave vector ${\boldsymbol{k}}$ is aligned with ${\boldsymbol{B}}$₀, the amplitudes of the Alfvén and slow modes vanish, and the fast-mode realization requires an implicit decomposition as the base $\{{{\boldsymbol{e}}}^{+},{{\boldsymbol{e}}}^{-}\}$ is undefined.

Figure 7 presents typical examples of the distribution of the random GMF from the local generator. In comparison to the magnetic field distribution from the global generator where the spatial anisotropy is defined by the orientation alignment, the local generator is capable of realizing more subtle field properties, e.g., the spectral anisotropic MHD wave types described in Section 4. At the phenomenological level, the global generator can mimic the random magnetic field orientation alignment of the local realizations as illustrated by Figure 6 and Figure 7, but the spectral anisotropy is uniquely realizable by the local generator.

Zoom In Zoom Out Reset image size

A realistic formulation of the local turbulent GMF is essential in this work, where simple random-field generators usually cannot take into account the anisotropy imprinted on the wave-vector phases of the power spectrum. The local generator we have designed in hammurabi X is capable of carrying out a theoretical parameterization of MHD turbulent modes, which have been discussed by Cho & Lazarian (2002), Caldwell et al. (2016), and Kandel et al. (2017, 2018). As described in these references, the turbulent field power spectra for Alfvén, fast, and slow modes can be formulated as

$$\begin{eqnarray}&&{P}_{i}(k,\alpha )={P}_{i}(k){F}_{i}({M}_{{\rm{A}}},\alpha ){h}_{i}(\beta ,\alpha ),\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}\begin{array}{rcl}{P}_{i}(k) & = & \displaystyle \frac{{p}_{i}}{4\pi {k}^{2}}\left[{\left(\displaystyle \frac{{k}_{0}}{{k}_{1}}\right)}^{{\alpha }_{1}}{\left(\displaystyle \frac{k}{{k}_{1}}\right)}^{6}{ \mathcal H }({k}_{1}-k)\right.\\ & & +{\left(\displaystyle \frac{k}{{k}_{0}}\right)}^{-{\alpha }_{1}}{ \mathcal H }(k-{k}_{1}){ \mathcal H }({k}_{0}-k)\\ & & +\left.{\left(\displaystyle \frac{k}{{k}_{0}}\right)}^{-{\delta }_{i}}{ \mathcal H }(k-{k}_{0})\right],\end{array}\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{h}_{{\rm{A}}}=1,\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{h}_{f}=\displaystyle \frac{2}{{D}_{++}(1+{\tan }^{2}\alpha {D}_{-+}^{2}/{D}_{+-}^{2})},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{h}_{s}=\displaystyle \frac{2}{{D}_{-+}(1+{\tan }^{2}\alpha {D}_{++}^{2}/{D}_{--}^{2})},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{D}_{\pm \pm }=1\pm \sqrt{D}\pm 0.5\beta ,\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&D={\left(1+0.5\beta \right)}^{2}-2\beta {\cos }^{2}\alpha ,\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{F}_{f}=1,\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{F}_{{\rm{A}},s}=\exp \left[-\displaystyle \frac{| \cos \alpha | }{{\left({M}_{{\rm{A}}}^{2}\sin \alpha \right)}^{2/3}}\right],\end{eqnarray} \tag{ 41 }$$

where $i\in \{{\rm{A}},f,s\}$ representing Alfvén, fast, and slow modes, respectively.²⁵ The two critical MHD parameters are the Alfvén Mach number M_A and the plasma β which is the ratio of gas pressure to magnetic pressure. In the sub-Alfvénic (M_A < 1) low-β (β < 1) regime, the spectral indices in Equation (34) can be approximated as δ_A = δ_s = 5/3 and δ_f = 3/2 (Cho & Lazarian 2002). The Alfvén speed v_A, which should appear in h_i(α), is absorbed by the normalization factor _pi for simplicity.

With the improved precision in hammurabi X, we present high-resolution Galactic synchrotron emission simulations with analytic models as described above. Presented in Figure 8 are the examples of synchrotron polarization at high Galactic latitudes predicted by a uniform regular GMF parallel to the Galactic disk and a random component from the global generator with a Kolmogorov power spectrum. Maps of synchrotron polarization from the same regular GMF but the local generator using a compressible MHD model are presented in Figure 9. Because we are presenting only illustrative models, the absolute strength of regular and random GMF is not essential here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The most prominent feature of the high-latitude synchrotron polarization is the quadrupolar structure that results from the GMF orientation at the solar neighborhood. As in the examples displayed in Figure 8, the quadrupole direction is largely determined by the regular field, but on top of which we can observe a flip in the polarization between the regimes when ρ > 1.0 versus ρ < 1.0. When the random GMF has no preferred alignment, i.e., the ρ = 1.0 case, the quadrupole pattern is undermined by the isotropic random-field contribution. This is visually clear because the random-field strength dominates. In Figure 9, the quadrupole pattern is well preserved with MHD turbulence injection scale k₀ = 10 kpc⁻¹, and also a flip in the polarization can be observed with the pure Alfvén mode when the random field dominates. When the spatial distribution or random GMF is close to spatially isotropic²⁶ with P_A/P_f,s = 3.0 (and Alfvén Mach number M_A = 0.5, plasma parameter β = 0.1) as displayed by the top panel in Figure 7, we observe a similar trend of weakening quadrupole pattern as demonstrated by Figure 9.

The synchrotron polarization fraction (or the degree of linear polarization) is mainly determined by the CRE spectral shape when a uniformly distributed regular GMF dominates. Assuming a constant CRE spectral index α = 3.0, the synchrotron polarization fraction Π = (3α + 3)/(3α + 7) is much higher than that observed from Planck data (Planck Collaboration et al. 2016b). Figures 10 and 11 demonstrate that the synchrotron polarization fraction can be suppressed by a Gaussian random field as long as the random field is not strongly anisotropic in the spatial domain. The suppression in polarization fraction grows with increasing random-field strength but depends on the specific field modeling. Recall that the addition of a random component to the magnetic field direction functions as a random walk in the polarization plane, which means that even for a purely turbulent field, the polarized intensity continues to increase with the number of turbulent cells added along the LoS. In principle, the increase goes as the square root of the number of cells, while the total intensity increases linearly, so the fraction should decrease accordingly. In practice, the precise trend is complicated by the effect of the observational beam and the locally varying anisotropy. The shape of the polarization fraction for the ρ = 0.5 model in Figure 10, for example, is due to the anisotropic random field canceling with the regular field before beginning to dominate. An inhomogeneous distribution (by field strength modulation) of the random field can change the efficiency of suppression differently depending on the field alignment, but the common features described above are preserved.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The above analyses imply that interpreting the synchrotron polarization toward the poles as due to the local field direction ignores the possible effects of anisotropic turbulence, which can mimic or flip the morphology. Though the physical process is different, the geometry of the field and its effect on the observables is the same for polarized dust emission. This work illustrates the opportunity for retrieving useful information on the local magnetic turbulence structure with high-latitude Galactic polarized emission and also shows the challenge from the degeneracy between random and regular magnetic field orientations when using emission data alone. It suggests that we need to be careful about realizing the local GMF structure to avoid misleading conclusions. For example, it has been proposed recently by Alves et al. (2018) that according to observations, the regular magnetic field structure may play a dominant role in Galactic dust emission near the solar neighborhood. We also emphasize that the Galactic synchrotron emission is also affected by the warm ISM in the Galactic thick disk and even the halo. The random-field generators in hammurabi X can be used to bridge the gap between simple large-scale field models and computationally intensive MHD simulations, and push toward more realistic analysis and modeling than previous methods.

The large angular scale Galactic synchrotron polarization pattern driven mainly by the GMF orientation at the solar neighborhood is quite evident as illustrated in Figures 8 and 9. However, the small angular structures can be analyzed with the angular power spectrum, which can be decomposed by rotation-invariant components, i.e., the T, E, and B modes (Hu & White 1997a). With the two random-field generators proposed in this work, we intend to figure out which properties of the random GMF are imprinted on the synchrotron B/E ratio. Specifically, we are interested in verifying whether MHD turbulence modes are capable of producing B/E < 1.0 in both the perturbative and the nonperturbative regimes. Because we are focusing on high-latitude polarization, pixels at Galactic latitude within ±60° are masked out. We also set a lower limit to the radius in the LoS integral according to the random-field grid resolution and the spherical mode range. Technical details of the precision checks for the pseudo-C_ℓ estimation are discussed in Appendix C.

We present in Figure 12 the B/E ratio distribution (by collecting results from an ensemble of realizations with each given parameter set) for varying random-field strengths and alignments of the global random GMF. Figure 12 implies that to reproduce B/E < 1.0, we either need random GMF in the nonperturbative regime (b/B₀ > 1.0) or parallel alignment (ρ > 1.0). We also note that the divergence-cleaning step is what leads to $B/E\ne 1.0$. As illustrated in the same figure, all realizations end up with B/E = 1.0 regardless of random-field alignment, when the Gram–Schmidt process is switched off. This is expected given that a simple Gaussian random field should have E = B on average, whereas a magnetic field must be divergence-free and therefore the difference between the naive random vector field and the magnetic field, which has been ignored in many previous analyses, is crucial in studying Galactic emissions. Now we conclude that the divergence-free random magnetic field can provide synchrotron ${\text{}}B/{\text{}}E\ne 1.0$. The Gram–Schmidt cleaning method is computationally useful and correct for reproducing the divergence-free random magnetic field (which in the simplest case can alternatively be obtained from a Gaussian random vector potential as shown in Appendix D, where synchrotron B/E < 1 arises naturally out of either method in the nonperturbative regime) and has the added benefit that we can spatially modulate its strength and orientation.

Zoom In Zoom Out Reset image size

By contrast, the C_ℓ estimated from the local MHD realizations have a clear analytic representation, to which simulations can be directly compared. To look for the low B/E ratio according to Kandel et al. (2018), we keep the random GMF strength at the perturbative level and tune the MHD Mach number M_A = 0.2 and plasma parameter β = 0.1. As illustrated in Figure 13, we find clear evidence that a Gaussian realization of MHD turbulence can provide a synchrotron B/E ratio smaller than 1.0, in both perturbative and nonperturbative regimes. The fast mode in a sub-Alfvénic low-β plasma has a unique power spectrum shape and is less affected by the anisotropy function h(α) than the slow mode. By assuming equal power in the turbulence modes at the injection scale, the observed angular power spectra are mainly influenced by the fast mode and so the B/E ratio has a different behavior for the case where slow and Alfvén modes dominate. With the given MHD Mach number and plasma parameter, slow-mode turbulence results in a much lower B/E ratio than that from the Alfvén mode, while fast mode prefers B/E ≃ 0.8 in the perturbative regime. These features are conceptually consistent with analytic predictions by Kandel et al. (2018) as demonstrated in the top panel of Figure 13, where the differences between two estimations are likely because of the simplification in analytic derivation, e.g., the Limber and flat-sky approximations. Beyond the perturbative regime, we observe that the B/E ratio evolves with the growth of the random-field strength and suggests an upper limit for the random-field strength to achieve the observed B/E ratio with solely MHD turbulence.

Zoom In Zoom Out Reset image size

The observational implications of the Galactic synchrotron emission from the above two types of random-field realizations are that both the divergence-free and MHD turbulent nature of the field are important for producing synchrotron B/E < 1.0 (aside from the fact that the divergence-free condition is physically required). It is possible to use directly the angular power spectra estimated in the way presented here for studying Galactic components like the work by Vansyngel et al. (2018), but we should be aware of the numeric uncertainty if the simulation resolution is lower than that of astrophysical measurements, in addition to the fundamental difference between simulation and observation mentioned in Section 2.

With the CRE differential flux distribution Φ(E, ${\boldsymbol{r}}$), synchrotron total and polarized emissivities at a given observational frequency ν and spatial position ${\boldsymbol{r}}$ read

$$\begin{eqnarray}&&{j}_{\mathrm{tot}/\mathrm{pol}}(\nu ,{\boldsymbol{r}})=\displaystyle \frac{1}{4\pi }{\int }_{{E}_{1}}^{{E}_{2}}{dE}\displaystyle \frac{4\pi }{\beta c}{\rm{\Phi }}(E,{\boldsymbol{r}})2\pi {P}_{\mathrm{tot}/\mathrm{pol}}(\omega ),\end{eqnarray} \tag{ 42 }$$

where ${P}_{\mathrm{tot}/\mathrm{pol}}(\omega )$, which represents the emission power from one electron at frequency ν = ω/2π, is calculated (Rybicki & Lightman 1979) through synchrotron functions $F(x)\,=x{\int }_{x}^{\infty }{K}_{\tfrac{5}{3}}(\xi )d\xi$ and $G(x)={{xK}}_{\tfrac{2}{3}}(x)$ (with ${K}_{\tfrac{5}{3}}(x)$ and ${K}_{\tfrac{2}{3}}(x)$ known as two of the modified Bessel functions of the second kind) as

$$\begin{eqnarray}&&{P}_{\mathrm{tot}}(\omega )=\displaystyle \frac{\sqrt{3}{e}^{3}{B}_{\mathrm{per}}}{2\pi {m}_{{\rm{e}}}{c}^{2}}F(x),\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{pol}}(\omega )=\displaystyle \frac{\sqrt{3}{e}^{3}{B}_{\mathrm{per}}}{2\pi {m}_{{\rm{e}}}{c}^{2}}G(x),\end{eqnarray} \tag{ 44 }$$

where e is the electron charge, m_e the electron mass, and B_per (defined as $| {\boldsymbol{B}}\times \hat{{\boldsymbol{n}}}|$ in Section 2) represents the strength of the magnetic field projected in the direction perpendicular to the LoS direction. Statistically, we assume that the synchrotron emission at a given position is isotropic, and so an observer only receives 1/4π of the emission power, which explains the 1/4π coefficient in front of the right-hand-side in Equation (42). In addition, we place an extra 2π before ${P}_{\mathrm{tot}/\mathrm{pol}}(\omega )$, due to the relation P(ν) = 2πP(ω). The term $\tfrac{4\pi }{\beta c}{\rm{\Phi }}(E,{\boldsymbol{r}})$, with β representing the relativistic speed, is actually N(E, ${\boldsymbol{r}}$), the CRE differential density.

In practice, the CRE spectral integral can be achieved in two technically different approaches with the same theoretical origin. If given a numerical CRE flux information Φ(E) prepared on a discrete grid, the integral Equation (42) can be directly evaluated by the numerical integral. Alternatively, we can start with an analytic differential density distribution $N(\gamma ,{\boldsymbol{r}})=4\pi {\rm{\Phi }}(E,{\boldsymbol{r}}){m}_{{\rm{e}}}c/\beta$, and by doing so, Equation (42) reads

$$\begin{eqnarray}&&{j}_{\mathrm{tot}/\mathrm{pol}}(\nu ,{\boldsymbol{r}})=\displaystyle \frac{1}{2}{\int }_{{\gamma }_{1}}^{{\gamma }_{2}}d\gamma N(\gamma ,{\boldsymbol{r}}){P}_{\mathrm{tot}/\mathrm{pol}}(\omega ).\end{eqnarray} \tag{ 45 }$$

The reason for keeping Equation (45) as an alternative method is to calculate the integral analytically once the CRE spectral index is constant at any given position as illustrated in Section 2. The detailed derivation follows the auxiliary definition of

$$\begin{eqnarray}&&{\omega }_{c}=\displaystyle \frac{3}{2}{\gamma }^{2}\displaystyle \frac{{{eB}}_{\mathrm{per}}}{{m}_{{\rm{e}}}c},\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&x=\displaystyle \frac{\omega }{{\omega }_{c}}.\end{eqnarray} \tag{ 47 }$$

Then, by assuming $N(\gamma )={N}_{0}{\gamma }^{-\alpha }$, Equation (45) ends up in the form of

$$\begin{eqnarray}\begin{array}{rcl}{j}_{\mathrm{tot}}(\nu ,{\boldsymbol{r}}) & = & \displaystyle \frac{\sqrt{3}{e}^{3}{B}_{\mathrm{per}}{N}_{0}}{8\pi {m}_{{\rm{e}}}{c}^{2}}{\left(\displaystyle \frac{4\pi \nu {m}_{{\rm{e}}}c}{3{{eB}}_{\mathrm{per}}}\right)}^{(1-\alpha )/2}\\ & & \times \int {dxF}(x){x}^{(\alpha -3)/2},\end{array}\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}\begin{array}{rcl}{j}_{\mathrm{pol}}(\nu ,{\boldsymbol{r}}) & = & \displaystyle \frac{\sqrt{3}{e}^{3}{B}_{\mathrm{per}}{N}_{0}}{8\pi {m}_{{\rm{e}}}{c}^{2}}{\left(\displaystyle \frac{4\pi \nu {m}_{{\rm{e}}}c}{3{{eB}}_{\mathrm{per}}}\right)}^{(1-\alpha )/2}\\ & & \times \int {dxG}(x){x}^{(\alpha -3)/2},\end{array}\end{eqnarray} \tag{ 49 }$$

where the spectral integrals can be analytically calculated by using

$$\begin{eqnarray}&&\int {dxF}(x){x}^{\mu }=\displaystyle \frac{{2}^{\mu +1}}{\mu +2}{\rm{\Gamma }}\left(\displaystyle \frac{\mu }{2}+\displaystyle \frac{7}{3}\right){\rm{\Gamma }}\left(\displaystyle \frac{\mu }{2}+\displaystyle \frac{2}{3}\right),\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&\int {dxG}(x){x}^{\mu }={2}^{\mu }{\rm{\Gamma }}\left(\displaystyle \frac{\mu }{2}+\displaystyle \frac{4}{3}\right){\rm{\Gamma }}\left(\displaystyle \frac{\mu }{2}+\displaystyle \frac{2}{3}\right).\end{eqnarray} \tag{ 51 }$$

Figure 14 illustrates the dependence of the synchrotron total emissivity T_tot and polarized emissivity T_pol on the CRE energy, with varying magnetic field strengths, observational frequencies, and CRE spectral shapes. The peaks in emissivities are inherited from F(x) and G(x), where the dimensionless parameter x is the ratio of the observational frequency to the CRE gyrofrequency.

Zoom In Zoom Out Reset image size

In this work, we focus on simulating synchrotron emission at the GHz level, for which the Galactic environment is optically thin (Rybicki & Lightman 1979; Schlickeiser 2002), and so we ignore both synchrotron self-absorption and free–free absorption. For readers who might be confused with the synchrotron emissivity calculation formulae presented above, please turn to the hammurabi X wiki page for more technical details.

Faraday rotation describes the phenomenological manifestation of the refractive index difference in the polarization directions for photons that propagate through a plasma with an external magnetic field. For a linearly polarized photon emitted with wavelength λ and intrinsic polarization angle χ₀, the observed polarization angle after traversing distance s₀ is

$$\begin{eqnarray}&&\chi ={\chi }_{0}+\phi ({s}_{0}){\lambda }^{2},\end{eqnarray} \tag{ 52 }$$

where ϕ, the Faraday depth, reads

$$\begin{eqnarray}&&\phi ({s}_{0})=\displaystyle \frac{{e}^{3}}{2\pi {m}_{{\rm{e}}}^{2}{c}^{4}}{\int }_{0}^{{s}_{0}}{{dsN}}_{{\rm{e}}}(s\hat{{\boldsymbol{p}}}){\boldsymbol{B}}(s\hat{{\boldsymbol{p}}})\cdot \hat{{\boldsymbol{p}}},\end{eqnarray} \tag{ 53 }$$

where $\hat{{\boldsymbol{p}}}$ represents the photon propagation direction and N_e represents the distribution of TE density. Note that the IAU convention²⁸ for polarization is adopted in hammurabi X, which means that the intrinsic synchrotron polarization angle is determined by the polarization ellipse semimajor axis perpendicular to magnetic field orientation. Under Faraday rotation at a given observational frequency ν, the observed emission accumulates Stokes parameters dQ and dU over a distance s₀ by

$$\begin{eqnarray}&&{dQ}+{idU}={{dI}}_{\nu }^{p}\exp \{2i\chi \},\end{eqnarray} \tag{ 54 }$$

where ${{dI}}_{\nu }^{p}$ represents the polarized intensity in the radial bin $[{s}_{0},{s}_{0}+{ds}]$. Though Faraday rotation brings in extra information about the TE distribution, a relatively high observational frequency is sometimes preferred for studying synchrotron emission, e.g., 30 GHz in this report, to suppress the complicated effects of TE turbulence, which will be addressed in our future studies with hammurabi X.