Covariant polarized radiative transfer on cosmological scales for investigating large-scale magnetic field structures

TL;DR

This paper develops a covariant polarized radiative transfer framework for cosmological scales, enabling detailed modeling of polarization signals to study large-scale magnetic fields and aiding future observations like those from SKA.

Contribution

It introduces a fully covariant CPRT formulation including Faraday effects, and proposes an all-sky ray-tracing algorithm for modeling polarization maps in cosmological simulations.

Findings

Verified the CPRT code with multiple tests.

Demonstrated applications using galaxy cluster and magnetized universe models.

Discussed implications for large-scale magnetic field studies.

Abstract

Polarization of radiation is a powerful tool to study cosmic magnetism and analysis of polarization can be used as a diagnostic tool for large-scale structures. In this paper, we present a solid theoretical foundation for using polarized light to investigate large-scale magnetic field structures: the cosmological polarized radiative transfer (CPRT) formulation. The CPRT formulation is fully covariant. It accounts for cosmological and relativistic effects in a self-consistent manner and explicitly treats Faraday rotation, as well as Faraday conversion, emission, and absorption processes. The formulation is derived from the first principles of conservation of phase-space volume and photon number. Without loss of generality, we consider a flat Friedmann-Robertson-Walker (FRW) space-time metric and construct the corresponding polarized radiative transfer equations. We propose an all-sky CPRT calculation algorithm, based on a ray-tracing method, which allows cosmological simulation results to be incorporated and, thereby, model templates of polarization maps to be constructed. Such maps will be crucial in our interpretation of polarized data, such as those to be collected by the Square Kilometer Array (SKA). We describe several tests which are used for verifying the code and demonstrate applications in the study of the polarization signatures in different distributions of electron number density and magnetic fields. We present a pencil-beam CPRT calculation and an all-sky calculation, using a simulated galaxy cluster or a model magnetized universe obtained from GCMHD+ simulations as the respective input structures. The implications on large-scale magnetic field studies are discussed; remarks on the standard methods using rotation measure are highlighted.

Figures

Tables

Table 1. Table 1: Properties of different intervening plasma models used in this paper. To test the ability of our CPRT equation solver to handle the extreme limits, the total electron number density n e , tot subscript 𝑛 e tot n_{\rm e,tot} for models A is set equal to the mean electron number density of the Universe (see Appendix  E for details); temperature of the thermal electrons in the IGM-like and ICM-like plasma models are assumed to take the lower-end values typical to IGM and ICM. ℱ nt subscript ℱ nt {\mathcal{F}}_{\rm nt} denotes the non-thermal relativistic electron fraction, p 𝑝 p denotes the power-law index of the energy spectrum of the non-thermal relativistic electrons, which relates to the spectral index of the synchrotron radiation α = ( p − 1 ) / 2 𝛼 𝑝 1 2 \alpha=(p-1)/2 . γ i subscript 𝛾 𝑖 \gamma_{i} is the electrons’ low-energy cutoff Lorentz factor. | B | B |\textbf{{B}}| denotes the magnetic field strength. The magnetic field direction along the line-of-sight is described by cos ⁡ θ ∈ [ − 1.0 , 1.0 ] 𝜃 1.0 1.0 \cos{\theta}\in[-1.0,1.0] which is set random for Models A-II and B-II.

	IGM-like plasma		ICM-like plasma
	A-I	A-II	B-I	B-II
$n_{e,\mathrm{tot}}$ (${\mathrm{cm}}^{-3}$)	$2.1918\times {10}^{-7}$		$1.00\times {10}^{-3}$
${ℱ}_{\mathrm{nt}}$ ($\%$)	$1.00$		$1.00$
$T_{e,\mathrm{th}}$ (K)	$1.875\times {10}^3$		$5.00\times {10}^5$
$p$	$4.00$		$2.50$
corresponds to $\alpha$	$1.50$		$0.75$
${\gamma }_i$	$10.0$		$30.0$
$\left|\text{B}\right|$ (G)	$1.00\times {10}^{-9}$		$1.00\times {10}^{-6}$
$\cos \theta$	0.5	$[-1.0,1.0]$	0.5	$[-1.0,1.0]$

Table 2. Table 2: Results of the multiray code-verification test where two Gaussian profiles, originating at z ori = 5.94 subscript 𝑧 ori 5.94 z_{\rm ori}=5.94 and at z ori = 1.00 subscript 𝑧 ori 1.00 z_{\rm ori}=1.00 respectively, are cosmologically transported in a vacuum but an expanding flat space–time. Four parameters are compared against their theoretical values; the empirical results are found to be consistent with the expected values up to machine floating-point precision.

		Profile I	Profile II
$z_{\mathrm{ori}}$		$5.93623409775142$	$1.00082825012323$
$z_{\mathrm{central}}^{\mathrm{refine}}$		$5.90499449730497$	$0.9992434724235666$
Peak frequency	Input ${\nu }_{\mathrm{central}}^{\mathrm{in}}$ (GHz)	$1.41400848881021$	$1.41928070323115$
	Output ${\nu }_{\mathrm{peak},0}$ (GHz)	$1.41400848881021$	$1.41928070323115$
	Fractional difference $\left({\nu }_{\mathrm{peak},0}-{\nu }_{\mathrm{central}}^{\mathrm{in}}\right)/{\nu }_{\mathrm{central}}^{\mathrm{in}}$	$-8.99280649946616\times {10}^{-15}$	$-1.110223037256493\times {10}^{-16}$
Dispersion	Input ${\sigma }^{\mathrm{in}}$	$0.000866648853601$	$-0.000123250853536$
	Output ${\sigma }_0$	$0.000866648853601$	$-0.000123250853536$
	Fractional difference $\left({\sigma }_0-{\sigma }^{\mathrm{in}}\right)/{\sigma }^{\mathrm{in}}$	$1.59377719355\times {10}^{-17}$	$5.66495635124\times {10}^{-18}$
Peak intensity ratio	Analytical $r_I^{\mathrm{ana}}$	$0.00299659998884$	$0.12484483163785$
	Empirical $r_I^{\mathrm{emp}}$	$0.00299659998884$	$0.12484483163785$
	Fractional difference$\left(r_I^{\mathrm{emp}}-r_I^{\mathrm{ana}}\right)/r_I^{\mathrm{ana}}$	$4.34172933267\times {10}^{-16}$	$1.06713879307\times {10}^{-14}$
Power-law index of	Analytical $m^{\mathrm{ana}}$	$3.0000$
$(I_{\mathrm{peak},0}^{z_{\mathrm{ori}}=5.94}/$	Empirical $m^{\mathrm{emp}}$	$3.00000000000099$
$I_{\mathrm{peak},0}^{z_{\mathrm{ori}}=1.00})$	Fractional difference $\left(m^{\mathrm{emp}}-m^{\mathrm{ana}}\right)/m^{\mathrm{ana}}$	$-3.29218134236\times {10}^{-13}$

Table 3. Table 3: Numerical results of the CPRT calculations for the demonstrative cases where bright point source is (i) absent, (ii) located at z = 6.0 𝑧 6.0 z=6.0 or (iii) located at z = 0.206 𝑧 0.206 z=0.206 ; magnetic fields orientate along the line-of-sight at random angles (see Section 5.1.1 ). The Stokes parameters are in units of erg ​ s − 1 ​ cm − 2 ​ Hz − 1 ​ str − 1 erg superscript s 1 superscript cm 2 superscript Hz 1 superscript str 1 {\rm erg}\,{\rm s}^{-1}\,{\rm cm}^{-2}\,{\rm Hz}^{-1}\,{\rm str}^{-1} , φ 𝜑 \varphi is measured in radian, and Π l subscript Π l \Pi_{\rm l} , Π c subscript Π c \Pi_{\rm c} , and Π tot subscript Π tot \Pi_{\rm tot} are expressed in percentages. Note that for case (i) the resulting I 𝐼 I has an order of magnitude 10 − 23 superscript 10 23 10^{-23} , which is much smaller than the specific intensity of the cosmic microwave background of 10 − 18 superscript 10 18 10^{-18}   erg ​ s − 1 ​ cm − 2 ​ Hz − 1 ​ str − 1 erg superscript s 1 superscript cm 2 superscript Hz 1 superscript str 1 {\rm erg}\,{\rm s}^{-1}\,{\rm cm}^{-2}\,{\rm Hz}^{-1}\,{\rm str}^{-1} at the same observed frequency. This suggests that emission and polarization signals would be overwhelmed by the CMB background in real observations.

	Without point source	Point source at $z=6.0$	Point source at $z=0.206$
initial Stokes parameters ${\left[\begin{array}{c}I\\ Q\\ U\\ V\end{array}\right]}_{z=6.0}$	$\left[\begin{array}{c}0.0000\\ 0.0000\\ 0.0000\\ 0.0000\end{array}\right]$	$\left[\begin{array}{c}8.7096\times {10}^{-14}\\ -4.5372\times {10}^{-16}\\ 2.5731\times {10}^{-15}\\ 4.3548\times {10}^{-17}\end{array}\right]$	$\left[\begin{array}{c}0.0000\\ 0.0000\\ 0.0000\\ 0.0000\end{array}\right]$
final Stokes parameters ${\left[\begin{array}{c}I\\ Q\\ U\\ V\end{array}\right]}_{z=0.0}$	$\left[\begin{array}{c}1.0438\times {10}^{-23}\\ 6.7617\times {10}^{-25}\\ 5.1681\times {10}^{-26}\\ -1.7613\times {10}^{-32}\end{array}\right]$	$\left[\begin{array}{c}2.5392\times {10}^{-16}\\ -2.8095\times {10}^{-18}\\ 7.0807\times {10}^{-18}\\ 1.2696\times {10}^{-19}\end{array}\right]$	$\left[\begin{array}{c}2.5392\times {10}^{-16}\\ -1.3221\times {10}^{-18}\\ 7.5021\times {10}^{-18}\\ 1.2696\times {10}^{-19}\end{array}\right]$
initial $\phi \left(z=6.0\right)$	0.0000	0.8727	0.0000
final $\phi \left(z=0.0\right)$	$3.8142\times {10}^{-2}$	0.9731	0.8726
initial ${\Pi }_l\left(z=6.0\right)$	0.0000	$3.0000$	0.0000
final ${\Pi }_l\left(z=0.0\right)$	$6.4969$	$3.0000$	$3.0000$
initial ${\Pi }_c\left(z=6.0\right)$	0.0000	$5.0000\times {10}^{-2}$	0.0000
final ${\Pi }_c\left(z=0.0\right)$	$1.6874\times {10}^{-7}$	$5.0000\times {10}^{-2}$	$5.0000\times {10}^{-2}$
initial ${\Pi }_{\mathrm{tot}}\left(z=6.0\right)$	0.0000	$3.0004$	0.0000
final ${\Pi }_{\mathrm{tot}}\left(z=0.0\right)$	$6.4969$	$3.0004$	$3.0004$

Table 4. Table 4: Statistics of the input and output parameters at z = 0 𝑧 0 z=0 of the demonstrative pencil-beam CPRT calculation using the simulated galaxy cluster obtained from a GCMHD+ cosmological MHD simulation; see Section 5.2 . n e , tot subscript 𝑛 e tot n_{\rm e,\,tot} is in units of cm − 3 superscript cm 3 {\rm cm}^{-3} , while | B | B |{\textbf{{B}}}| is in G. The Stokes parameters are in units of  erg ​ s − 1 ​ cm − 2 ​ Hz − 1 ​ str − 1 erg superscript s 1 superscript cm 2 superscript Hz 1 superscript str 1 {\rm erg}\,{\rm s}^{-1}\,{\rm cm}^{-2}\,{\rm Hz}^{-1}\,{\rm str}^{-1} , Δ ​ φ Δ 𝜑 \Delta\varphi is in radian, and Π l subscript Π l \Pi_{\rm l} , Π c subscript Π c \Pi_{\rm c} , Π tot subscript Π tot \Pi_{\rm tot} are in per cent. All values are corrected to four decimal places for compactness.

	Mean	Standard Deviation	Minimum	Maximum
Input
$n_{e,\mathrm{tot}}$	$6.3561\times {10}^{-5}$	$2.7001\times {10}^{-4}$	$6.3577\times {10}^{-7}$	$2.6508\times {10}^{-2}$
$\left|\text{B}\right|$	$1.5585\times {10}^{-8}$	$5.0621\times {10}^{-8}$	$2.5227\times {10}^{-14}$	$2.5175\times {10}^{-6}$
$\cos \theta$	$-5.9235\times {10}^{-3}$	$5.7434\times {10}^{-8}$	-1.0000	1.0000
Output
$I$	$7.0217\times {10}^{-18}$	$4.4204\times {10}^{-18}$	$2.2564\times {10}^{-27}$	$1.3000\times {10}^{-16}$
$Q$	$1.6270\times {10}^{-19}$	$1.0358\times {10}^{-18}$	$-2.96286\times {10}^{-19}$	$3.3099\times {10}^{-17}$
$U$	$-2.3702\times {10}^{-21}$	$9.3950\times {10}^{-19}$	$-2.9899\times {10}^{-17}$	$3.2602\times {10}^{-17}$
$V$	$-2.6325\times {10}^{-25}$	$4.2232\times {10}^{-23}$	$-1.5752\times {10}^{-21}$	$1.9561\times {10}^{-21}$
$\Delta \phi$	$1.1957\times {10}^{-3}$	$0.2543$	$-3.1240$	$3.1401$
${\Pi }_l$	$68.5725$	$8.0302$	$1.1191$	$70.5876$
${\Pi }_c$	$2.3278\times {10}^{-4}$	$2.8580\times {10}^{-4}$	$6.4526\times {10}^{-10}$	$4.2846\times {10}^{-3}$
${\Pi }_{\mathrm{tot}}$	$68.5725$	$8.0302$	$1.1191$	$70.5876$

Table 5. Table 5: Statistics of the input and output parameters at z = 0 𝑧 0 z=0 of the demonstrative all-sky CPRT calculation using a model magnetized universe obtained from a cosmological GCMHD+ simulation; see Section 5.3 . n e , tot subscript 𝑛 e tot n_{\rm e,\,tot} is in units of cm − 3 superscript cm 3 {\rm cm}^{-3} , while | B | B |{\textbf{{B}}}| is in G. The Stokes parameters are in units of  erg ​ s − 1 ​ cm − 2 ​ Hz − 1 ​ str − 1 erg superscript s 1 superscript cm 2 superscript Hz 1 superscript str 1 {\rm erg}\,{\rm s}^{-1}\,{\rm cm}^{-2}\,{\rm Hz}^{-1}\,{\rm str}^{-1} , Δ ​ φ Δ 𝜑 \Delta\varphi is in radian, and Π l subscript Π l \Pi_{\rm l} , Π c subscript Π c \Pi_{\rm c} , Π tot subscript Π tot \Pi_{\rm tot} are in per cent. All values are corrected to four decimal places for compactness.

	Mean	Standard Deviation	Minimum	Maximum
Input
$n_{e,\mathrm{tot}}$	$2.1918\times {10}^{-7}$	0.0000	$2.1918\times {10}^{-7}$	$2.1918\times {10}^{-7}$
$\left|\text{B}\right|$	$5.2855\times {10}^{-9}$	$6.8763\times {10}^{-9}$	$4.5734\times {10}^{-11}$	$1.9587\times {10}^{-7}$
$\cos \theta$	$-1.8850\times {10}^{-3}$	$0.5769$	$-1.0000$	$0.9999$
Output
$I$	$9.6274\times {10}^{-24}$	$3.3006\times {10}^{-26}$	$9.5894\times {10}^{-24}$	$1.1695\times {10}^{-23}$
$Q$	$4.7385\times {10}^{-26}$	$2.6055\times {10}^{-26}$	$1.7357\times {10}^{-26}$	$1.6790\times {10}^{-24}$
$U$	$2.9300\times {10}^{-31}$	$3.5445\times {10}^{-28}$	$-2.7505\times {10}^{-26}$	$1.0346\times {10}^{-26}$
$V$	$-1.0209\times {10}^{-33}$	$2.0177\times {10}^{-30}$	$-2.1578\times {10}^{-28}$	$8.4035\times {10}^{-29}$
$\Delta \phi$	$1.1597\times {10}^{-5}$	$4.7925\times {10}^{-1}$	$-3.0689\times {10}^{-1}$	$3.4868\times {10}^{-1}$
${\Pi }_l$	$4.9131\times {10}^{-3}$	$2.5886\times {10}^{-3}$	$1.8100\times {10}^{-3}$	$1.4359\times {10}^{-1}$
${\Pi }_c$	$5.4946\times {10}^{-8}$	$1.8629\times {10}^{-7}$	$1.5424\times {10}^{-12}$	$1.8452\times {10}^{-5}$
${\Pi }_{\mathrm{tot}}$	$4.9131\times {10}^{-3}$	$2.5886\times {10}^{-3}$	$1.8100\times {10}^{-3}$	$1.4359\times {10}^{-1}$

Table 6. Table 6: Values of the transfer coefficients and optical depths computed using parameters of models A-I and B-I at radiation frequency ν = 1.4 𝜈 1.4 \nu=1.4  GHz. The transfer coefficients obtained have a very small order of magnitude, suggesting a scale length of a few Mpc.

	IGM-like model A-I	ICM-like model B-I
${ϵ}_{I,\mathrm{tot}}$	${ϵ}_{I,\mathrm{th}}+{ϵ}_{I,\mathrm{nt}}$	${ϵ}_{I,\mathrm{th}}+{ϵ}_{I,\mathrm{nt}}$
($\mathrm{erg}{s}^{-1}{\mathrm{cm}}^{-3}{\mathrm{Hz}}^{-1}{\mathrm{str}}^{-1}$)	$=2.59\times {10}^{-53}+2.62\times {10}^{-55}$	$=6.91\times {10}^{-47}+1.25\times {10}^{-38}$
	$=2.62\times {10}^{-53}$	$=1.25\times {10}^{-38}$
${\kappa }_{\mathrm{tot}}$	${\kappa }_{\mathrm{th}}+{\kappa }_{\mathrm{nt}}$	${\kappa }_{\mathrm{th}}+{\kappa }_{\mathrm{nt}}$
(${\mathrm{cm}}^{-1}$)	$=2.23\times {10}^{-38}+8.64\times {10}^{-52}$	$=2.23\times {10}^{-34}+7.11\times {10}^{-34}$
	$=2.23\times {10}^{-38}$	$=9.34\times {10}^{-34}$
$\tau ={\int }_{z_{\mathrm{init}}}^{0.0}{\kappa }_{\mathrm{tot}}\left(z\right)\cdot ds$	$2.72\times {10}^{-13}$	$1.14\times {10}^{-8}$
$f_{\mathrm{tot}}$	$f_{\mathrm{th}}+f_{\mathrm{nt}}$	$f_{\mathrm{th}}+f_{\mathrm{nt}}$
(${\mathrm{cm}}^{-1}$)	$=2.54\times {10}^{-30}+1.52\times {10}^{-33}$	$=1.16\times {10}^{-23}+8.54\times {10}^{-28}$
	$=2.54\times {10}^{-30}$	$=1.16\times {10}^{-23}$
$h_{\mathrm{tot}}$	$h_{\mathrm{th}}+h_{\mathrm{nt}}$	$h_{\mathrm{th}}+h_{\mathrm{nt}}$
(${\mathrm{cm}}^{-1}$)	$=3.76\times {10}^{-42}+8.12\times {10}^{-43}$	$=1.72\times {10}^{-32}+3.01\times {10}^{-32}$
	$=4.57\times {10}^{-42}$	$=4.72\times {10}^{-32}$

Equations

$$\frac{d{I}_{i,\nu }}{ds}=-{\kappa }_{ij,\nu }{I}_{j,\nu }+{ϵ}_{i,\nu },$$

$$\displaystyle\frac{\rm d}{\operatorname{d}\!{s}}\left[\begin{array}[]{c}I_{\nu}\\ Q_{\nu}\\ U_{\nu}\\ V_{\nu}\end{array}\right]=-\left[\begin{array}[]{cccc}\kappa_{\nu}&q_{\nu}&u_{\nu}&v_{\nu}\\ q_{\nu}&\kappa_{\nu}&f_{\nu}&-g_{\nu}\\ u_{\nu}&-f_{\nu}&\kappa_{\nu}&h_{\nu}\\ v_{\nu}&g_{\nu}&-h_{\nu}&\kappa_{\nu}\end{array}\right]\left[\begin{array}[]{c}I_{\nu}\\ Q_{\nu}\\ U_{\nu}\\ V_{\nu}\end{array}\right]+\left[\begin{array}[]{c}\epsilon_{I,\nu}\\ \epsilon_{Q,\nu}\\ \epsilon_{U,\nu}\\ \epsilon_{V,\nu}\end{array}\right]$$

$$\frac{d{\mathcal{I}}_{\nu }}{d{\tau }_{\nu }}=-{\mathcal{I}}_{\nu }+\frac{{\xi }_{\nu }}{{\zeta }_{\nu }}=-{\mathcal{I}}_{\nu }+{\mathcal{S}}_{\nu }$$

$$v^{\beta }=P^{\alpha \beta }{k}_{\alpha }=k^{\beta }+(k_{\alpha }{u}^{\alpha })u^{\beta }$$

$$\frac{ds}{d{\lambda }_{\mathrm{a}}}$$

$$\displaystyle\frac{k_{\alpha}u^{\alpha}\Big{|}_{\lambda_{\rm a,co}}}{k_{\beta}u^{\beta}\Big{|}_{\lambda_{\rm a,obs}}}=\frac{\nu_{\rm co}}{\nu_{\rm obs}}\ ,$$

$$\frac{\mathrm{d}{\mathcal{I}}_{\nu }}{d{\lambda }_{\mathrm{a}}}$$

$$\displaystyle\frac{\operatorname{d}\!{(}\mathcal{I}_{i,\rm co})}{\operatorname{d}\!{\lambda}_{\rm a}}=\frac{\operatorname{d}\!{(}I_{i,\rm co}/\nu_{\rm co}^{3})}{\operatorname{d}\!{\lambda}_{\rm a}}=-k_{\alpha}u^{\alpha}\Big{|}_{\lambda_{\rm a,co}}\left\{-\kappa_{ij,{\rm co}}\,\left(\frac{I_{j}}{\nu_{\rm co}^{3}}\right)+\frac{\epsilon_{i,{\rm co}}}{\nu_{\rm co}^{3}}\right\}\ .$$

$$\displaystyle k^{\alpha}=\left[\begin{array}[]{c}E\\ p_{r}\\ p_{\theta}\\ p_{\phi}\end{array}\right]=\nu\left[\begin{array}[]{c}1\\ 1\\ 0\\ 0\end{array}\right]\ ;\hskip 22.76228ptu^{\beta}=\gamma\left[\begin{array}[]{c}1\\ \beta_{r}\\ \beta_{\theta}\\ \beta_{\phi}\end{array}\right]\ ,$$

$$\displaystyle{k_{\alpha}u^{\alpha}\Big{|}_{z}}$$

$$\displaystyle\frac{{k_{\alpha}u^{\alpha}\Big{|}_{z}}}{k_{\beta}u^{\beta}\Big{|}_{z_{\rm obs}}}$$

$$\displaystyle\frac{{k_{\alpha}u^{\alpha}\Big{|}_{z}}}{k_{\beta}u^{\beta}\Big{|}_{z_{\rm obs}}}=\frac{\nu_{\rm z}}{\nu_{z_{\rm obs}}}\ ,$$

$$\frac{\mathrm{d}}{\mathrm{d}{\lambda }_{\mathrm{a}}}=\frac{\mathrm{d}{x}^0}{\mathrm{d}{\lambda }_{\mathrm{a}}}\frac{\mathrm{d}}{\mathrm{d}{x}^0}=E\frac{\mathrm{d}}{\mathrm{d}s}=E\frac{\mathrm{d}z}{\mathrm{d}s}\frac{\mathrm{d}}{\mathrm{d}z},$$

$$\frac{{\nu }_z}{{\nu }_{z_{\mathrm{obs}}}}=\frac{a_{\mathrm{obs}}}{a}=\frac{1+z}{1+z_{\mathrm{obs}}},$$

$$\displaystyle\frac{\rm d}{\operatorname{d}\!{z}}\left[\begin{array}[]{c}{\mathcal{I}}\\ {\mathcal{Q}}\\ {\mathcal{U}}\\ {\mathcal{V}}\end{array}\right]=({1+z})\left\{-\left[\begin{array}[]{cccc}\kappa&q&u&v\\ q&\kappa&f&-g\\ u&-f&\kappa&h\\ v&g&-h&\kappa\end{array}\right]\left[\begin{array}[]{c}{\mathcal{I}}\\ {\mathcal{Q}}\\ {\mathcal{U}}\\ {\mathcal{V}}\end{array}\right]\right.+\left.\left[\begin{array}[]{c}\epsilon_{I}\\ \epsilon_{Q}\\ \epsilon_{U}\\ \epsilon_{V}\end{array}\right]\frac{1}{\nu^{3}}\right\}\frac{\operatorname{d}\!{s}}{\operatorname{d}\!{z}}\ ,$$

$$\frac{ds}{dz}=\frac{c}{H_0}(1+z{)}^{-1}{\left[{\Omega }_{\mathrm{r},0}(1+z{)}^4+{\Omega }_{\mathrm{m},0}(1+z{)}^3+{\Omega }_{\Lambda ,0}\right]}^{-\frac{1}{2}},$$

$$\displaystyle\log_{10}\left(\frac{B^{2}_{\parallel}(z)}{8\pi}\right)=\left\{\begin{array}[]{ll}8.1737\,x^{4}-40.352\,x^{3}+73.647\,x^{2}-55.264\,x-12.16&:0.64<x<1.70\\ 0.67\tanh(-x/0.18+2.72)-26.14&:0.15\leq x\leq 0.64\\ -\tanh(x/0.52+0.28)-24.91&:-2.00\leq x<0.15\ ,\end{array}\right.$$

$$\displaystyle{M}_{ij}=\left(\begin{array}[]{ccc}1&0&0\\ 0&\cos{\theta}&\sin{\theta}\\ 0&-\sin{\theta}&\cos{\theta}\end{array}\right)\ ,$$

$$\displaystyle\left(\begin{array}[]{c}{x}\\ {y}\\ {z}\end{array}\right)=\left(\begin{array}[]{ccc}1&0&0\\ 0&\cos{\theta}&-\sin{\theta}\\ 0&\sin{\theta}&\cos{\theta}\end{array}\right)\left(\begin{array}[]{c}{\tilde{x}}\\ {\tilde{y}}\\ {\tilde{z}}\end{array}\right)\ .$$

$$\displaystyle\textbf{{E}}(z,t)=\textbf{{E}}_{0}\,e^{\pm i(\omega t-kz)}=\left(\begin{array}[]{c}E_{x}(z,t)\\ E_{y}(z,t)\end{array}\right)\ =\left(\begin{array}[]{c}E_{x0}\,e^{\pm i(\omega t-kz+\phi_{x})}\\ E_{y0}\,e^{\pm i(\omega t-kz+\phi_{y})}\end{array}\right)\ ,$$

$$\displaystyle\textbf{{E}}^{\rm RCP}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\ \mp i\end{array}\right)\ ,$$

$$V_{\nu }=$$

$$=$$

$$f_{\mathrm{th}}$$

$$h_{\mathrm{th}}$$

$${\kappa }_{\mathrm{th}}$$

$$q_{\mathrm{th}}$$

$$v_{\mathrm{th}}$$

$${\nu }_{\mathrm{c}}$$

$$\displaystyle\Lambda=\Biggl{\{}\begin{array}[]{ll}\ \Big{(}\frac{2}{1.781}\Big{)}^{5/2}\Big{(}\frac{k_{\rm B}T_{\rm e}}{m_{\rm e}}\Big{)}^{1/2}\Big{(}\frac{k_{\rm B}T_{\rm e}}{e^{2}\omega}\Big{)}&,\>{\rm for}\ T_{\rm e}\leq 3.16\times 10^{5}~{}{\rm K}\\ \\ \ \frac{8\pi k_{\rm B}T_{\rm e}}{1.781h\omega}&,\>{\rm for}\ T_{\rm e}>3.16\times 10^{5}~{}{\rm K}\end{array}$$

$${ϵ}_{I,\mathrm{th}}={\kappa }_{\mathrm{th}}{B}_{\omega },{ϵ}_{Q,\mathrm{th}}=q_{\mathrm{th}}{B}_{\omega },\mathrm{and}{ϵ}_{V,\mathrm{th}}=v_{\mathrm{th}}{B}_{\omega },$$

$$\displaystyle\operatorname{d}\!{n}=[n_{\gamma}\gamma^{p}]\gamma^{-p}\Theta(\gamma-\gamma_{i})g(\Psi)\operatorname{d}\!{\gamma}\operatorname{d}\!{\Omega}_{\Psi}\$$

Full text

††pubyear: 2018††pagerange: Covariant polarized radiative transfer on cosmological scales for investigating large-scale magnetic field structures–F

1 Introduction

Signatures of magnetic fields are seen on all scales, from planets (see e.g.  Stevenson, 2003; Schubert & Soderlund, 2011) and stars (see e.g.  Parker, 1970; Brandenburg & Subramanian, 2005; Beck, 2008; Schrijver & Zwaan, 2008; Vallée, 1998, 2011b), to galaxies (Ferrière, 2009; Vallée, 2011a; Beck & Wielebinski, 2013; Planck Collaboration et al., 2015, 2016b) and galaxy clusters (see e.g.  Govoni et al., 2006; Guidetti et al., 2008; Vacca, V. et al., 2010; Pratley et al., 2013; Kronberg, 2016). Magnetic fields should also permeate large-scale structures such as superclusters (see e.g.  Xu et al., 2006), filaments (see e.g.  Ryu, Kang & Biermann, 1998; Brüggen et al., 2005; Ryu et al., 2008), walls, and voids (see e.g.  Beck et al., 2013), as the early-time magnetic seeds get amplified during the structure formation and evolution processes in the Universe (see e.g.  Widrow, 2002; Durrer & Neronov, 2013; Kronberg, 2016, for comprehensive reviews). However, observational evidence of these weak large-scale magnetic fields is scarce. Magnetic fields must have played a pivotal role in: (i) star formation by transporting angular momentum out from accretion discs and so allowing materials to accrete onto proto-stars (see e.g.  Balbus & Hawley, 1991), (ii) jet production by affecting the central accretion, as well as accelerating and collimating the materials that form jets (see e.g. Pudritz, Hardcastle & Gabuzda, 2012), (iii) cosmic-ray production through the acceleration of charged particles (Fermi, 1949), and (iv) cosmic-ray propagation through deflecting the ray or confining the charged particles (see e.g. Jokipii, 1966, 1967). However, the origins of large-scale magnetic fields, their co-evolution with astrophysical structures, and their properties at present are as of yet to be determined. How magnetic fields impact the formation of the first structures and their subsequent evolution remains a pressing problem in contemporary astrophysics and cosmology. The current understanding, however, will be revolutionized as high-quality all-sky polarization data set will become available from upcoming radio telescopes, such as the Square Kilometer Array (SKA)111https://www.skatelescope.org (see e.g. Beck & Gaensler, 2004; Feretti & Johnston-Hollitt, 2004; Gaensler, Beck & Feretti, 2004; Johnston-Hollitt et al., 2015).

Polarization surveys by emerging generation of radio telescopes, such as the GaLactic and Extragalactic All-sky MWA (GLEAM) Survey (Wayth et al., 2015) on the Murchison Widefield Array (MWA222http://www.mwatelescope.org/) (Tingay et al., 2013), the polarization Sky Survey of the Universe’s Magnetism (POSSUM) (Gaensler et al., 2010) on the Australian SKA Pathfinder (ASKAP333http://www.atnf.csiro.au/projects/askap/index.html) (Hotan et al., 2014), as well as the Multifrequency Snapshot Sky Survey (MSSS) (Heald et al., 2015b) on the Low-Frequency Array (LOFAR444http://www.lofar.org/) (van Haarlem et al., 2013), due to their improved sensitivities and resolutions, can already access a domain of weak magnetic field strengths that was unexplored before. They enable investigations of magnetism in a variety of astrophysical sources. These experiments pave the way for broad-band spectro-polarimetric surveys to be performed by the SKA, which will be a game-changer. The SKA is an interferometric radio telescope designated to have a total collecting area of about a square kilometer in its complete configuration. Its sensitivity, bandwidths, and field-of-view will provide a transformational polarization data set with which the detection of the very weak magnetic field of the cosmic web555Current upper limits on the intergalactic field strength are all model-dependent but generally fall within the range of $|\textbf{{B}}^{\rm IGM}|\leq 10^{-8}$ to $10^{-9}$ G (see e.g. Kronberg, 1994; Blasi, Burles & Olinto, 1999; Brown et al., 2017). may become possible (Giovannini et al., 2015; Vazza et al., 2015) – with the SKA the evolution of magnetism in galaxies and galaxy clusters may be traced (Gaensler et al., 2015), and the detailed internal structure of the magnetized cosmic plasmas (both across the medium and along the line-of-sight) may be mapped or imaged (Han et al., 2015; Heald et al., 2015a). An all-sky polarization survey performed with the SKA will significantly increase the density of known polarized background sources on the sky (Beck & Gaensler, 2004; Feretti & Johnston-Hollitt, 2004). These sources serve as distant radio backlights, illuminating the magnetized Universe via the effect of Faraday rotation, i.e. rotation of the polarization plane of radiation as it travels and interacts with the magnetic fields threaded in an ionized medium. Such a data set will be immensely rich, containing information of the polarized sources themselves, as well as the foreground sources lying along the line-of-sight. These sources can be the Milky Way, nearby or distant galaxies, galaxy clusters, and even the cosmic filaments connecting clusters of galaxies. With all the exciting opportunities opened up by observational advances, the pressing questions to be addressed now are: how do we uncover and characterize the polarization signals from data, and ultimately, use them to infer and quantify magnetic field properties? How do we confront our theoretical models of cosmic magnetism against observations? More specifically, how do we compare simulation results which encode physical model predictions to the results obtained by observational experiments?

This paper aims to address the second and third questions by providing a solid theoretical foundation and a polarized radiative transfer tool to investigate cosmic magnetism on large scales (i.e. Mpc scales and beyond). We present a new formulation of cosmological polarized radiative transfer (CPRT), which is fully covariant and is valid for polarization transfer in flat space–time. Our derivation is based on a covariant general relativistic radiative formulation stemmed from the first principles of conservation of phase–space volume and photon number (Fuerst & Wu, 2004; Younsi, Wu & Fuerst, 2012). The covariant CPRT equation allows the properties of the magnetic fields to be captured as they co-evolve with the structures in the expanding Universe. Furthermore, since our formulation accounts for the relativistic and cosmological effects in a self-consistent manner, polarization evolution in various cosmic media as a function of redshift can be investigated. The formulation preserves the basic structure of the conventional polarized radiative transfer (see e.g. Sazonov & Tsytovich, 1968; Sazonov, 1969; Jones & Odell, 1977a, b; Pacholczyk, 1977; Degl’innocenti & Degl’innocenti, 1985), making it easy to implement for practical calculations, as we demonstrate in example problems and applications. Moreover, the formulation is general: it can be reduced to the form from which the conventional rotation measure (RM) quantity (see e.g. Rybicki & Lightman, 1986) is derived, assuming the absence of emission and absorption, insignificant Faraday conversion, and negligible effects of non-thermal electrons in the medium (see On et al., 2019, for details and the generalization of the standard RM expression to account for an isotropic distribution of non-thermal relativistic electrons with a power-law energy spectrum). At the same time, since the CPRT formulation explicitly accounts for absorption, emission, and Faraday processes, its application is not restricted to any special cases. To our knowledge, our formulation of CPRT, which is applicable to study large-scale cosmic magnetism, is the first of its kind666Formulations and codes capable of computing general relativistic polarized radiative transfer (GRPRT) in the (curved) Kerr space–time metric have been extensively studied and presented Broderick & Blandford (2003, 2004); Shcherbakov & Huang (2011); Gammie & Leung (2012); Dexter (2016); Mościbrodzka & Gammie (2018). Their applications primary concern polarized emissions from magnetized accretion flows and jets around a spinning black hole. .

The CPRT formulation serves as a solid platform whereby, given some input distributions of electron number densities and magnetic fields, one can trace the rays and compute their intensities and polarization over redshifts. These inputs can either be generated by simple modeling or cosmological simulations. We devise and construct a ray-tracing algorithm that solves the CPRT equation, thereby constructing model templates and making theoretical all-sky intensity and polarization maps. These data outputs, when combined with advanced statistical methods for data analysis and characterization, will help us achieve a reliable interpretation of observational data, crucial for scientific extraction. Results obtained from such a forward approach also provide an experimental test-bed for assessing line-of-sight component separation methods and methods used for characterizing signals themselves and the underlying physical processes.

This paper focuses on laying the foundation of the cosmological polarized radiative transfer approach to study the structure of large-scale magnetic fields. In Section 2, we present the CPRT formulation and its derivation. Our ray-tracing algorithms for solving the CPRT equation are given in Section 3. Verification tests for the code implementation and their results are described in Section 4. Demonstrations of applying the CPRT calculations to practical astrophysical applications are discussed in Section 5. We perform a set of single-ray CPRT calculations, showcasing the ability of the tool to study the cosmological evolution of polarization with or without bright radio sources along the line-of-sight. We also demonstrate how to compute polarization maps of an astrophysical object and an entire polarized sky, interfacing cosmological MHD simulation results with the CPRT calculations. We highlight the implications of these calculations on large-scale magnetic field studies. In Section 6, we summarize the whole paper.

Unless otherwise specified, c.g.s. units and a $[\,-,+,+,+\,]$ signature are used throughout this work.

2 Cosmological polarized radiative transfer

The CPRT formulation is derived based on a covariant general relativistic radiative transfer (GRRT) formulation Fuerst & Wu (2004); Younsi, Wu & Fuerst (2012), stemming from the first principles of conservation of phase–space volume and photon number. We start off by reviewing the polarized radiative transfer equation and the GRRT formulation to derive the covariant CPRT formulation. Then, we construct the corresponding CPRT equations assuming a flat Friedmann–Robertson–Walker (FRW) space–time metric, without loss of generality.

2.1 Conventional polarized radiative transfer

We first set out the polarized radiative transfer (PRT) equation and show how the covariant formulation of radiative transfer can be directly generalized to that of the PRT.

In the absence of scattering, the transfer equation of polarized radiation, in tensor representation, can be written as

[TABLE]

or in the matrix form,

[TABLE]

(Sazonov, 1969; Pacholczyk, 1970, 1977; Jones & Odell, 1977a; Degl’innocenti & Degl’innocenti, 1985; Huang & Shcherbakov, 2011a; Janett et al., 2017; Janett, Steiner & Belluzzi, 2017; Janett & Paganini, 2018), where $s$ is the path length of the radiation; the tensor index $i$ or $j$ in equation (1) runs from 1 to 4, denoting the Stokes parameters $I_{\nu}$, $Q_{\nu}$, $U_{\nu}$ and $V_{\nu}$, respectively. The coefficient tensor $\kappa_{ij,\nu}$ accounts for the amount of absorption (through $\kappa_{\nu}$, $q_{\nu}$, $u_{\nu}$ and $v_{\nu}$), rotation (through $f_{\nu}$) and conversion (through $h_{\nu}$ and $g_{\nu}$) of the radiation along its direction of propagation, and $\epsilon_{i,\nu}$ accounts for the amount of emission. Essentially, absorption acts as a sink; emission serves as a source. Propagation effects of Faraday rotation and Faraday conversion are non-dispersive. Faraday rotation, due to circular birefringence (i.e. the slightly different speeds at which the left and right circularly waves travel in a magneto-ionic medium), results in the change of polarization angles as radiation propagates (i.e. $Q_{\nu}\leftrightarrow U_{\nu}$). Faraday conversion, due to linear birefringence, concerns with the interconversion between the linear and circular polarization modes of the radiation (i.e. $Q_{\nu}\leftrightarrow V_{\nu}$; $U_{\nu}\leftrightarrow V_{\nu}$). The equation for the transfer of polarized radiation presented at above, with variable transfer coefficients, is suitable for transport in a homogeneous or weakly anisotropic medium (Sazonov & Tsytovich, 1968; Sazonov, 1969; Pacholczyk, 1977; Jones & Odell, 1977a). Note that all quantities in equations (1) and (18) depend on the frequency of the radiation $\nu$.

It is useful to note that the Stokes parameters are observables fully describing the properties of light but are coordinate-system dependent quantities. They can be combined in the complex forms, i.e. $(Q_{\nu}\pm iU_{\nu})$, and be linearly transformed to so-called $E$- and $B$- modes, which describe, respectively, parity-odd polarization and parity-even polarization, and so are invariant under transform of coordinate systems. Stokes parameters alone are not rotationally invariant. Therefore, coordinate systems adopted, as well as the definitions and conventions of polarization, must be explicitly stated to remove any ambiguities in the interpretation of the Stokes results. We note that different handedness of coordinate systems (right-handed or left-handed), as well as the geometry of the problem, have been used in the literature that derived the polarized radiative transfer equations and the (thermal and non-thermal) transfer coefficients (see e.g. Sazonov, 1969; Pacholczyk, 1970; Melrose & McPhedran, 1991; Huang & Shcherbakov, 2011a). The sign of Stokes $V_{\nu}$ that describes the sense of the circular polarization also varies from paper to paper. Furthermore, different conventions have been used in the literature regarding the definition of the polarization angle777 Investigations of the polarization of the cosmic microwave background adopt the opposite convention to the International Astronomical Union (IAU) standard, for which polarization angle increases clockwise (counterclockwise) when looking at the source for the former (latter). To rectify the discrepancy requires an opposite sign applied to Stokes $U_{\nu}$ (see https://aas.org/posts/news/2015/12/iau-calls-consistency-use-polarization-angle)., the definition of the handedness of circular polarization, and the definition of $V_{\nu}$ (see Robishaw, 2008, for a compilation of the conventions used in radio polarization work). We thus define in Appendix A the coordinate systems and the geometry of the problem considered in this work, and discuss in Appendix B the intricacies of keeping a consistent polarization convention.

Here, we note that the $U_{\nu}$ components $u_{\nu}$, $g_{\nu}$ and $\epsilon_{U,\nu}$ can vanish (so $V_{\nu}$ couples only to $U_{\nu}$, i.e. $U_{\nu}\leftrightarrow V_{\nu}$ but $Q_{\nu}\nleftrightarrow V_{\nu}$) by a choice of a local coordinate system (see e.g. Sazonov, 1969; Pacholczyk, 1977). With the geometry defined in Fig. 13 in Appendix A, $u_{\nu}$, $g_{\nu}$ and $\epsilon_{U,\nu}$ become zero in the basis $(x,y)$ since the projection of the magnetic field onto the $(x,y)$-plane is parallel to $y$.

Another useful remark concerns the features of equations (1) and (18). They reduce to the usual scalar radiative transfer equation only when a specific intensity $I_{\nu}$ is considered, i.e. ${\rm d}I_{\nu}/{\operatorname{d}\!{s}}=-\kappa_{\nu}\,I_{\nu}+\epsilon_{\nu}$. Conversely, one can utilize the fact that all the Stokes parameters have the same physical units to easily include polarization in the covariant formulation of radiative transfer, as is outlined in the subsequent subsection.

2.2 Covariant general relativistic radiative transfer

From the first principles of conservation of photon number and phase space volume, it can be shown that the Lorentz-invariant intensity is given by $\mathcal{I}_{\nu}\equiv\ {I_{\nu}}/{\nu^{3}}$, and that the covariant formulation of the radiative transfer takes the form

[TABLE]

(see Appendix D), where $\tau_{\nu}=\int\kappa_{\nu}\,{\rm d}s$ is the optical depth, $\zeta_{\nu}=\nu\,\kappa_{\nu}$ and $\xi_{\nu}={\epsilon_{\nu}}/\nu^{2}$ are the Lorentz-invariant coefficients of absorption and emission respectively, and the Lorentz-invariant source function is defined by ${\mathcal{S}_{\nu}}\equiv{\xi_{\nu}}/{\zeta_{\nu}}$.

In relativistic settings, we want the covariant radiative transfer equation to be evaluated in space–time intervals instead of optical depth or path length. This can be achieved by introducing the mathematical affine parameter $\lambda_{\rm a}$. The problem is then translated into an evaluation of ${\operatorname{d}\!{s}}/{\operatorname{d}\!{\lambda}_{\rm a}}$ (i.e. the variation in the path length $s$ with respect to $\lambda_{\rm a}$), and asking the question of what is the co-moving 4-velocity $v^{\beta}$ of a photon traveling in a fluid that has 4-velocity $u^{\beta}$.

Assuming the photon has a 4-momentum $k^{\alpha}$, then the co-moving 4-velocity $v^{\beta}$ can be obtained by the projection of $k^{\alpha}$ on to the fluid frame, i.e.

[TABLE]

(Fuerst & Wu, 2004), where we have used the projection tensor $P^{\alpha\beta}=g^{\alpha\beta}+u^{\alpha}u^{\beta}$, with $g^{\alpha\beta}$ as the space–time metric tensor. The variation in $s$ with respect to $\lambda_{\rm a}$ is therefore

[TABLE]

(Younsi, Wu & Fuerst, 2012). Note that for a stationary observer positioned at infinity $k_{\beta}u^{\beta}=-E_{\rm obs}$. It follows that the ratio

[TABLE]

which corresponds to the relative energy shift of the photon between the observer’s frame and the comoving frame. Using the Lorentz-invariant properties of $\mathcal{I}_{\nu}$, $\zeta_{\nu}$ and $\xi_{\nu}$ yields the covariant relativistic radiative transfer equation

[TABLE]

(Younsi, Wu & Fuerst, 2012), where all the quantities are frequency dependent and are evaluated along the path of a photon, i.e. comoving as denoted by the subscript “co".

2.3 Cosmological polarized radiative transfer formulation

The CPRT formulation is constructed by making two generalizations to the GRRT: (i) by accounting for the polarization of the radiation and (ii) by incorporating a cosmological model to describe the space–time geometry of the Universe in which the radiation propagates. The former generalization is straightforward in the sense that the PRT equation takes the general form of radiative transfer (see Section 2.1) and that all the Stokes parameters have the same physical units. Therefore, similar to how one can obtain the Lorentz-invariant intensity by taking $\mathcal{I}_{\nu}\equiv{I_{\nu}}/{\nu^{3}}$, the invariant Stokes parameters are obtained by $\mathcal{I}_{\nu,i}=[{\mathcal{I}}_{\nu},{\mathcal{Q}}_{\nu},{\mathcal{U}}_{\nu},{\mathcal{V}}_{\nu}]^{\rm T}=[I_{\nu},Q_{\nu},U_{\nu},V_{\nu}]^{\rm T}/\nu^{3}$ where the tensor index $i$ runs from 1 to 4, and the superscript “T" denotes the transpose (for notational simplicity we drop the subscript $\nu$ of the Stokes parameters and in the coefficients of absorption and emission hereafter). It follows that the covariant polarized radiative transfer equation, in tensor notation, takes the form

[TABLE]

Next, to make the formulation appropriate in cosmological settings and, therefore, suitable for (but not limited to) the investigation of cosmological magnetic fields, the factor $k_{\alpha}u^{\alpha}$ is to be evaluated using the space–time metric of a chosen cosmological model such that equation (24) is evaluated in terms of a cosmological variable, e.g.  the redshift $z$, instead of the mathematical affine parameter $\lambda_{\rm a}$.

Without loss of generality, we consider a flat FRW universe whose space–time metric has the diagonal elements  ($-1,a^{2},a^{2},a^{2}$), where $a=1/(1+z)$ is the cosmological scale factor describing the expansion of the universe. For simplicity, we consider a photon propagating radially in a cosmological medium with 4-velocity $u^{\beta}$, i.e.

[TABLE]

where ${\textbf{{p}}}=(p_{r},p_{\theta},p_{\phi})$ denotes the 3-velocity of the photon, ${\beta}$ $=(\beta_{r},\beta_{\theta},\beta_{\phi})$ denotes the 3-velocity of the medium, and $\gamma=1/\sqrt{(1+\beta^{2})}$ is the corresponding Lorentz factor (here, we use $c=h=1$). Evaluating $k_{\alpha}u^{\alpha}$ then yields

[TABLE]

and the ratio

[TABLE]

If the motion of the medium can be neglected (i.e. $\beta$=0, $\gamma$=1), the ratio is simplified to

[TABLE]

which is the relative shift of energy (or frequency) of the photon, as one expects from equation (22). By defining $k^{\alpha}=(E,{\textbf{{p}}})={\operatorname{d}\!{x^{\alpha}}}/{\operatorname{d}\!{\lambda}_{\rm a}}$, one may also obtain

[TABLE]

and use this to also show that the photon’s energy $E\propto{a}^{-1}$ and thus

[TABLE]

in a flat FRW universe (see e.g.  Dodelson, 2003). In other words, the ratio in equation (40) corresponds to the relative energy shift of the photon due to the cosmic expansion.

Finally, by applying the chain rule given in equation (41) to equation (24), we obtain the CPRT equation defined in redshift space:

[TABLE]

where all the quantities are Lorentz invariant and ${\operatorname{d}\!{s}}/{\operatorname{d}\!{z}}$ for a flat FRW universe is given by

[TABLE]

(see e.g.  Peacock, 1999), where $H_{0}$ is the standard Hubble parameter, $\Omega_{{\rm r},0}$, $\Omega_{{\rm m},0}$ and $\Omega_{\Lambda,0}$ are the dimensionless energy densities of relativistic matter and radiation, non-relativistic matter, and a cosmological constant (dark energy with an equation of state of $w\equiv-1$), respectively. The subscript “0" denotes that the quantities are measured at the present epoch (i.e. $z=0$).

Note that the CPRT formulation is general and can adopt different cosmological models with flat space–time geometry through the $k_{\alpha}u^{\alpha}$ factor. Ray-tracing calculation for equation (59) can then be performed for arbitrary photon geodesics. For clarity, we reiterate that a flat space–time is considered in our derivation such that straightforward parallel transport of the polarization Stokes vector $S_{\nu,i}=[I_{\nu},Q_{\nu},U_{\nu},V_{\nu}]^{\rm T}$ of the radiation along the photon geodesics is enabled888 For completeness, we note that polarized radiative transfer in Kerr space–time has been extensively studied (Broderick & Blandford, 2003, 2004; Shcherbakov & Huang, 2011; Gammie & Leung, 2012; Dexter, 2016; Mościbrodzka & Gammie, 2018), for which the standard approach involves solving the photon geodesic, keeping track of the local coordinate system such that polarized emission is being added appropriately in the presence of a rotation of the coordinate system propagated along the ray, and finally, connecting these frames to the polarization frame at the point of observation. Difficulties stem from that the Stokes parameters are not rotationally invariant quantities. Working with rotationally invariant quantities, e.g. the spin-2 signals of $E$- and $B$-modes, might therefore be more favourable. . For radiation propagating in a curved space–time, the rotation of its polarization vector measured by the observer has a contribution caused not only by the Faraday rotation but also by the curvature of the embedded manifold, i.e. angle is not preserved transporting along the line-of-sight. Taking advantage of the flat geometry of the Universe (Planck Collaboration et al., 2016a), we therefore limit our evaluation and discussion to a cosmological model describing a flat universe only. The flatness of space–time ensures that the angles measured in the local comoving frame would be the same everywhere along the geodesic.

We highlight that the covariant nature of the CPRT formulation allows a straightforward transform of an observable between the comoving frame and the observer’s frame. Computation from the invariant Stokes parameters to the observable Stokes parameters in the comoving frame requires only a scalar multiplication of the cube of the radiation frequency, i.e. $[I_{\nu}(z),Q_{\nu}(z),U_{\nu}(z),V_{\nu}(z)]^{\rm T}=[{\mathcal{I}}_{\nu}(z),{\mathcal{Q}}_{\nu}(z),{\mathcal{U}}_{\nu}(z),{\mathcal{V}}_{\nu}(z)]^{\rm T}\times\nu(z)^{3}$. The results at $z=0$ are then what would be measured in the observer’s frame at the present time, provided that the transform of the local polarization frame to the instrument’s polarization frame are properly handled (as is noted in Appendix A), along with the corrections of instrumental effects and foregrounds, such as ionospheric effects.

2.4 Polarized transfer coefficients

Following the derivation of the CPRT equation, i.e. equation (59), in this subsection we discuss the corresponding transfer coefficients appropriate for the context of cosmic plasmas. The expressions of the coefficients considered in this paper are explicitly specified in Appendix C.

An astrophysical plasma generally consists of a population of thermal electrons and a population of non-thermal electrons, which could be relativistic electrons that gyrate around magnetic field lines, electrons accelerated by shocks, or electrons injected by cosmic rays. Given that dielectric suppression999Dielectric suppression, or known as the Razin effect or Razin-Tsytovich effect (Razin, 1960; Ramaty, 1968), is a plasma effect on synchrotron emission. Synchrotron radiation is suppressed exponentially below the Razin frequency $\omega_{\rm R}=\omega_{\rm p}^{2}/\omega_{\rm B}$, where $\omega_{\rm p}$ is the plasma angular frequency and $\omega_{\rm B}$ is the electron angular gyrofrequency, since the electrons can no longer maintain the phase with the emitted radiation as the wave phase velocity would increase to above the speed of light (see e.g. Melrose, 1980). (see e.g. Bekefi, 1966; Rybicki & Lightman, 1986) is insignificant, which is generally true for cosmic plasmas (see e.g. Melrose & McPhedran, 1991), here the transfer coefficients are expressed as the sum of their respective thermal and non-thermal components, i.e. $\kappa_{ij}=(\kappa_{ij,{\rm th}}+\kappa_{ij,{\rm nt}})$ and $\epsilon_{i}=(\epsilon_{i,{\rm th}}+\epsilon_{i,{\rm nt}})$, where “th" and “nt" denote the thermal and non-thermal components of the absorption and emission coefficients respectively.

In this work, we consider thermal bremsstrahlung and non-thermal synchrotron radiation process101010 In addition to thermal bremsstrahlung and non-thermal synchrotron radiation process considered in this work, we note that transfer coefficients appropriate for different astrophysical environments have been extensively studied in the literature. Accurate expressions for the coefficients of Faraday rotation and Faraday conversion in uniformly magnetized relativistic plasmas, such as those in jets and hot accretion flows around black-holes, are reported in Huang & Shcherbakov (2011b). Expressions of the transfer coefficients in the case of ultra-relativistic plasma that is permeated by a static uniform magnetic fields, for frequencies of high harmonic number limits, and for a number of distribution functions (isotropic, thermal, or power law) are presented in Heyvaerts et al. (2013). Emission and absorption coefficients for cyclotron process, that is important in accretion discs of compact objects, have also been studied by Chanmugam et al. (1989); Vaeth & Chanmugam (1995). Careful incorporation of the above would be a useful improvement to the current CPRT implementation, expanding the range of its applications and enabling a realistic modeling of the magnetized Universe. . For thermal bremsstrahlung, expressions of the Faraday rotation coefficient $f_{\rm th}$ and Faraday conversion coefficient $h_{\rm th}$, as well as the expressions of the absorption coefficients $\kappa_{\rm th}$, $q_{\rm th}$ and $v_{\rm th}$ follow Pacholczyk (1977)111111The same expressions of $\kappa_{\rm th}$, $q_{\rm th}$ and $v_{\rm th}$ are provided in Wickramasinghe & Meggitt (1985) but a typo of an extra factor of the square of angular frequency is found in the denominator of $v_{\rm th}$ via dimensional analysis. We also note that the sign of $q_{\rm th}$ in Wickramasinghe & Meggitt (1985) is also different to Pacholczyk (1977), which might be due to different polarization sign conventions or a sign error.. The emission coefficients are computed via Kirchoff’s law accordingly. For non-thermal synchrotron emission, we consider relativistic electrons that have a power-law energy distribution. We use the expressions of the transfer coefficients that follow Jones & Odell (1977a) and consider an isotropic distribution of relativistic electrons’ momentum direction.

As detailed in Appendix B, the sign of Stokes $V$ depends on its definition, polarization conventions, handedness of the coordinate systems used, as well as the time dependence of the electromagnetic wave (i.e. whether the exponent has $+i\omega t$ or $-i\omega t$) and the definition of the relative phase between the $x$ and $y$-components of the electric field of the radiation. However, some of these information were not explicitly stated in Jones & Odell (1977a), and inconsistent definitions of the time dependence of the electromagnetic wave were used in Pacholczyk (1977) in deriving the radiative transfer coefficients for bremsstrahlung and synchrotron radiation process (see their equations 3.33 and 3.93). We therefore expound on the strategy to eliminate ambiguity in Appendix B and present a consistent set of expressions of all the transfer coefficients in Appendix C, given the geometry explicitly defined in Appendix A and the polarization convention conforming to the IEEE/IAU standard.

3 Algorithms and Numerical Calculation

The CPRT equations given in equation (59) can, in principle, be either solved by direct integration via numerical methods, or by diagonalizing and determining the inverse of the transfer matrix operator. We adopt the former approach and employ a ray-tracing method in this paper. In this section, we present algorithms to solve the CPRT equation numerically. We first present the algorithm for computing the CPRT for a single ray, followed by the algorithm for an all-sky setting wherein cosmological MHD simulation results may be incorporated to generate a set of theoretical all-sky intensity and polarization maps.

3.1 Ray-tracing

The CPRT algorithm consists of three basic components concerning (i) the interaction of radiation with the line-of-sight plasmas, (ii) the cosmological effects on radiation and the co-evolution of plasmas with the Universe’s history, and (iii) numerical computation of the CPRT equation, which is a set of four coupled differential equations evaluated in the redshift $z$-space. In the following, we discuss each of these components, starting with the numerical method. We describe the implementation of the algorithm and highlight its specific designs to accommodate the inclusion of line-of-sight astrophysical sources and intervening plasmas of different properties.

3.1.1 Numerical method

The radiation propagation is parameterized by redshift $z$ and is sampled discretely into $N_{\rm cell}$ number of cells. We adopt a sampling scheme such that each $z$-interval corresponds to an approximately equal light travel distance. That is, between the initial redshift $z_{\rm init}$ at which we start evaluating the CPRT equation and the final redshift $z=0$ at which observation is made, the total light travel distance $s_{\rm tot}$ is first computed by solving equation (60) followed by finding the corresponding lower and upper boundary values of $z$; where in each $z$-interval, light travels a distance as close to $s_{\rm eq}=s_{\rm tot}/N_{\rm cell}$ as possible. Note that the light travel distance acts as a scaling factor in the context of numerical evaluation. For efficient numerical computation, its multiplications with the transfer coefficients would ideally be close to unity.

Each $z$ interval can be further refined to incorporate astrophysical structure(s) and their sub-structure(s). Our code implementation allows an option to switch on/off such a refinement scheme, as well as to incorporate multiple structures at different redshifts. Within the refinement zone, we employ a uniform sampling in the $\log_{10}{(1+z)}$ space which has the advantage of preserving the profile shape when multi-frequency calculations are to be carried out. In algorithmic terms, at the cell of index $ind^{\rm refine}$, the increment over each refined cell is given by $[\log_{10}{(1+z^{\prime})}-\log_{10}{(1+z)}]/N^{\rm refine}_{\rm cell}$, where $z^{\prime}$ and $z$ are, respectively, the upper and lower boundaries of the $z$-interval to be refined, and $N^{\rm refine}_{\rm cell}$ is the total number of refined cells.

A fourth-order Runge?Kutta (RK) differential equation solver is used to integrate equation (60), and to solve equation (59), which ultimately gives us the Stokes parameters $\{I,Q,U,V\}$ at $z=0$ in the observer’s frame. Parameters to be set for the solver include the total number of (coupled) differential equations to be solved $N_{\rm eqn}$, the number of steps for the RK solver $N_{\rm step}$; and the error tolerance level eps. The error estimation of the solver is carried out by comparing the solution obtained with a fourth-order RK formula to that obtained with a fifth-order RK formula. If the computed error is less than eps then the calculation proceeds; otherwise the algorithm halts, reports errors of non-convergence, and returns without further computation.

The upper and lower limits of the $z$-variables are updated along the ray. The outputs are passed into the next computation as the inputs (i.e. as updated initial conditions). Since the evaluation of the CPRT starts from a higher $z$ to a lower z value until the present $z_{0}=0$ is reached, a substitution of $z\rightarrow-z$ is made in equation (59) as we set that as the function to be evaluated by the RK solver.

3.1.2 Interaction of radiation with plasmas

Radiation is parameterized by frequency $\nu(z)$, which has a redshift dependence of $\nu(z)=\nu_{\rm obs}(1+z)$, where $\nu_{\rm obs}$ is the observed frequency at the present epoch $z=0$. Its intensity and polarization properties change when passing through the magnetized intervening plasmas. The strength of the radiative processes, captured through transfer coefficients in the CPRT equation, depends on the physical properties of the plasmas, in addition to the frequency of the radiation. In general, both thermal and non-thermal electrons are present in astrophysical plasmas. Parameters describing them include: $n_{\rm e,tot}$, fraction of non-thermal electrons ${\mathcal{F}}_{\rm nt}$ (and thus $n_{\rm e,th}$ and $n_{\rm e,nt}$), temperatures $T_{\rm e}$ for thermal electrons, power-law index of the non-thermal electrons’ energy spectrum $p$ and the electrons’ low energy cutoff described by the Lorentz factor $\gamma_{i}$. Added to this list are parameters describing the strength and direction of magnetic fields, B, which can be decomposed into two components. One component is decomposed along the line-of-sight direction $\textbf{{B}}_{\parallel}=|\textbf{{B}}|\cos{\theta}$, and another component is decomposed in the plane normal to the line-of-sight $\textbf{{B}}_{\perp}=|\textbf{{B}}|\sin{\theta}$, where $\theta$ is the angle between the direction of the magnetic field and the line-of-sight.

By specifying the observed frequency of radiation $\nu_{\rm obs}$ at $z=0$ and the radiation background at an initial redshift $z_{\rm init}$, and given some input distributions of electron number density $n_{\rm e}(z)$ and magnetic field strength $|\textbf{{B}}(z)|$ through which light travels, solving the CPRT equations yields the evolution of the intensity and polarization of the radiation as a function of $z$.

3.1.3 Cosmological effects

In this paper, we adopt the maximum likelihood cosmological parameters obtained by the Planck Collaboration et al. (2016a) with the present Hubble constant $H_{0}=100\,h_{0}=67.74~{}\rm{kms^{-1}Mpc^{-1}}$, the matter density today $\Omega_{\rm m,0}=0.3089$, and the cosmological constant or vacuum density today $\Omega_{\Lambda,0}=0.6911$ (Planck Collaboration et al., 2016a). The radiation density today is given by $\Omega_{\rm r,0}=4.1650\times 10^{-5}(h_{0})^{-2}$ (Wright, 2006).

We have already noted the frequency shift of the radiation due to the expansion of the Universe, i.e. $\nu(z)=\nu_{\rm obs}(1+z)$. The cosmological (expansion) effects on the temperatures, electron number densities, as well as the strengths of magnetic fields are given by, respectively, $T_{\rm e}(z)=T_{\rm e,0}(1+z)^{2}$, $n_{\rm e}(z)=n_{\rm e,0}(1+z)^{3}$, and $|\textbf{{B}}(z)|=|\textbf{{B}}_{0}|(1+z)^{2}$, assuming frozen-in flux condition. These properties, as well as the structures of magnetic fields, are also subjected to local structure formation, evolution and outflows, as well as to influences by external injections, such as cosmic rays. Consequently, the inter-stellar medium (ISM), intra-cluster medium (ICM), and intergalactic media (IGM) all exhibit different characteristic properties. The CPRT formulation is covariant and accounts for cosmological and relativistic effects self-consistently. Because of these advantages, theoretical predictions of the intensity and polarization of the radiation can be computed straightforwardly by incorporating simulation results describing the cosmic plasmas into the computation of the transfer coefficients, and then solving the CPRT equation.

3.2 All-sky polarization calculation

We construct an all-sky CPRT algorithm that can interface with cosmological simulation results, numerically solve the CPRT equation, and thereby generate theoretical all-sky polarization maps that serve as model templates. A schematic of the algorithm is shown in Fig. 1.

Fig. 2 illustrates the concept of the all-sky algorithm, in which the CPRT equation is solved in a spherical polar coordinate system $(r,\theta,\phi)$, where $(\theta,\phi)$ corresponds to the celestial sky coordinates and the radial axis $r$ corresponds to the redshift axis $z$. Note that outputs of the cosmological evolutions of plasma properties, e.g. $n_{\rm e}(z)$ and $|\textbf{{B}}(z)|$, obtained from a cosmological MHD simulation can be inputted to the CPRT calculations through the transfer coefficients. Spatial fluctuations of the plasma properties in a finite simulation volume, usually in Cartesian coordinate system $(i,j,k)$, can also be mapped to the spherical polar coordinate system $(r,\theta,\phi)$ at each sampled redshift $z$. A more rigorous treatment that guarantees the magnetic field is divergence-free is also possible within our all-sky framework. We will present these details in a forthcoming paper.

We add a remark here on the sampling scheme over a sphere for efficient follow-on data analysis. There is an option to compute rays that are randomly positioned over the entire celestial sphere. Alternatively, one may utilize the advantages of efficient spherical sampling schemes, such as the HEALPix sampling (Górski et al., 2005) and the sampling scheme devised by McEwen & Wiaux (2011) which affords exact numerical quadrature. In such a case, ray-tracing CPRT calculation is performed at each grid point on the sphere. Map data constructed this way allows efficient power spectrum analyses and spherical wavelet analyses (e.g. McEwen, Hobson & Lasenby, 2006; Sanz et al., 2006; Starck et al., 2006; Geller et al., 2008; Marinucci et al., 2008; Wiaux et al., 2008; Leistedt et al., 2013; McEwen, Vandergheynst & Wiaux, 2013; McEwen et al., 2015; McEwen, Durastanti & Wiaux, 2018; Chan et al., 2017) to characterize the spatial fluctuations of polarization, crucial for searching polarization signatures imprinted by large-scale magnetic fields in observational data.

4 Code verification

In this section, we present the single-ray and multiple-ray experiments performed for code verification121212Consistency test is also performed by comparing the results of light-travel time obtained by integrating equation (60) using our code (then dividing by the speed of light) to those that are obtained using the publicly available cosmological calculator by Wright (2006), http://www.astro.ucla.edu/~wright/CosmoCalc.html. The results agree with each other, up to the maximum digits displayed in Wright (2006), i.e. three decimal places.. The $z$-sampling scheme follows the recipe described in Section 3.1.1 (or see the related red boxes in Fig. 1). We consider polarized radiative transfer at frequencies $\nu_{\rm obs}=1.42$ GHz and $5.00$ GHz for illustrative purposes131313$\nu_{\rm obs}=1.4$ GHz is chosen since it lies within the operating range of many current and upcoming radio telescopes, such as the Arecibo radio telescope (http://www.naic.edu/), the Five hundred meter Aperture Spherical Telescope (FAST, http://fast.bao.ac.cn), the Australia Telescope Compact Array (ATCA, https://www.narrabri.atnf.csiro.au/), LOFAR, MWA, ASKAP, SKA, etc.. Properties of the intervening plasma considered are listed in Table 1, which can be IGM-like (model A) or ICM-like (model B) with magnetic field directions along the line-of-sight set at a fixed angle (models A-I and B-I) or set as randomly oriented (models A-II and B-II). Thermal bremsstrahlung and non-thermal synchrotron radiation process are accounted for.

4.1 Single-ray tests

To test the accuracy and precision of our CPRT integrator in handling polarized radiative transfer in scenarios investigated in this paper, we repeat the two integration tests presented in Section 3.2 in Dexter (2016). We solve the standard PRT equation (equation (18)) that is reduced from the general CPRT equation ((59), and compare the numerical solutions that we obtained with the analytic solutions, which are explicitly given in Appendix C of Dexter (2016) for the idealized situations with constant transfer coefficients along a ray.

In the first test, we consider pure emission and absorption in Stokes $I$ and $Q$. The light ray travels through the Faraday-thin IGM-like plasma or the Faraday-thick ICM-like plasma (models A-I or B-I) over a cosmological distance from $z=6$ to $z=0$. Detailed values of both the (thermal and non-thermal) emission and absorption transfer coefficients, as well as the optical depths used in the calculations are given in Table 6 in Appendix F. As is seen in Fig. 3, the numerical solution obtained by our CPRT integrator agrees with the analytical solution up to the machine floating-point precision throughout the entire light path.

In the second test, we consider radiation of observed frequencies $\nu_{\rm obs}=1.4$ GHz and $\nu_{\rm obs}=5.0$ GHz. The radiation travels through the Faraday-thick ICM-like plasma (B-I) of a few Mpc in length scale. Only pure Faraday rotation and Faraday conversion and polarized emission in $Q$ and $V$ are considered (note that $\epsilon_{U}$ is set to zero due to the choice of coordinate systems (see Appendix 2.1)). To ease checking the oscillatory behavior of the resulting $V$, we boost the Faraday conversion effect artificially by setting its transfer coefficient to the same order magnitude as the Faraday rotation coefficient. The results of the second test is presented in Fig. 4. An excellent agreement between the numerical and analytic solutions is obtained in both cases of different radiation frequencies. Machine floating-point precision is maintained over the ray despite that the residuals in $Q$, $U$ and $V$ increase with each oscillation. Similar trend is also found in Fig. 4 in Dexter (2016) and that in Mościbrodzka & Gammie (2018).

4.2 Multiple-ray test

To verify the redshift-refinement scheme and the whole code, we performed multiray calculations evaluating the cosmological radiative transfer of two Gaussian profiles that centered at two different redshifts. The two input profiles, originating at $z_{\rm ori}=5.94$ and $1.00$, have frequency samples assigned through the redshift-refinement scheme at that specific $z_{\rm ori}$. The central frequency of the profiles is then given by $\nu^{\rm in}_{\rm central}=\nu_{\rm obs}(1+z^{\rm refine}_{\rm central})$, where $z^{\rm refine}_{\rm central}$ is the redshift value of the $N^{\rm refine}/2$ cell, for $N^{\rm refine}=500$. The ray then freely propagates in vacuum afterwards, i.e. all the transfer coefficients are set to zero when computing the CPRT equation. As such, frequency shift of the radiation is the only cosmological effect which modifies the radiation properties in its transport. We compare the values of four quantities obtained from the CPRT calculations against the theoretical expected values. These quantities are (i) the frequency at which the resulting profile peaks, $\nu_{\rm peak,0}$, (ii) the standard deviation of the resulting profile, $\sigma_{0}$, (iii) the empirical ratio of the output to the input peak intensity $r^{\rm emp}_{I}=I^{\rm in}_{\rm peak}/I_{\rm peak,0}$, for each Gaussian profile, and (iv) the power-law index of the ratio of the output peak intensities of the two profiles, $m^{\rm emp}$. Analytically, the resulting profile obtained from the CPRT of each case (i.e. emission at $z_{\rm ori}=5.94$ or at $z_{\rm ori}=1.00$) should remain Gaussian and peak at the frequency of $\nu_{\rm obs}\times(1+z^{\rm refine}_{\rm central})/(1+z_{\rm ori})$ with $\nu_{\rm obs}=1.42$ GHz. The standard deviation of the normalized input and the output Gaussian profiles should remain the same. The ratio of the peak intensity of the output emission profile to that of the input profile follows $r^{\rm ana}_{I}=1/(1+z_{\rm ori})^{3}$. Furthermore, comparing the outputs of the two cases, the ratio of the peak intensity at zero redshift should follow a power law of $[(1+z_{\rm ori}^{\prime\prime})/(1+z_{\rm ori}^{\prime})]^{3}$, where $z^{\prime\prime}$ denotes the higher redshift. That is, the power-law index $m^{\rm ana}=3.0$.

We summarize the results in Table 2, from which one can see that the empirical results are consistent with the theoretical expectation up to machine floating-point precision. Furthermore, consistent results are obtained using the parallelized code (i.e.  with multiple threading using OpenMP) as those obtained by the serial execution.

5 Applications

Here, we present a set of CPRT calculations to demonstrate the ability of the algorithm in tracking the change of polarization on astrophysical and cosmological scales. Changes in polarization features caused by the frequency shift of the radiation, or those caused by the evolution of intervening cosmic plasmas can be separately investigated; direct studies of their combined effects can also be directly carried out.

We start with a set of single-ray calculations, showing in our case studies how polarization changes over cosmological distances with and without a bright line-of-sight point source. Then we demonstrate how to incorporate cosmological MHD simulation results into CPRT calculations to make polarization maps. We compute the polarization of a simulated galaxy cluster. We also compute the entire polarized sky using a model magnetized universe. Polarization maps generated in such a way, i.e. by CPRT calculations with an interface of simulation results, encapsulate theoretical predictions. They are crucial to aid our interpretation of observational data. Model templates of the entire sky are particularly important for comparison with future observational data, such as those from all-sky surveys of polarized emission with the SKA.

5.1 Cosmological evolution of polarization

We perform a set of ray-tracing calculations for radiation with observed frequency $\nu_{\rm obs}=\nu_{0}=1.42$ GHz propagating from $z=6.0$ through some distributions of $n_{\rm e}(z)$ and $|\textbf{{B}}(z)|$ as described below. The $z$-sampling scheme follows the recipe described in Section 3.1.1.

5.1.1 Point-source emissions

Bright polarized emitters such as quasars and radio galaxies may lie along the line-of-sight acting as back-light illuminating the foreground. Here, we calculate how the polarization and intensity of a fiducial quasar-like point source changes over a cosmological distance. Emissions of such a point source at $z$ observed at $1.42$ GHz is given by $[{I},{Q},{U},{V}]|_{z}=[{I},{Q},{U},{V}]|_{z=0}(1+z)^{3}$, where $[{I},{Q},{U},{V}]|_{z=0}=[\,2.54\times 10^{-16},-1.32\times 10^{-18},7.50\times 10^{-18},1.27\times 10^{-19}]$  ${\rm erg}\,{\rm s}^{-1}\,{\rm cm}^{-2}\,{\rm Hz}^{-1}\,{\rm str}^{-1}$, where we have assumed the degree of linear polarization to be $3.00\%$ (Jagers et al., 1982), the degree of circular polarization to be $0.05\%$ (Conway et al., 1971), and the polarization angle $\varphi=0.87$ rad. For demonstrative purposes, we adopt such a simple interpolation of $[{I},{Q},{U},{V}]|_{z}$ from $[{I},{Q},{U},{V}]|_{z=0}$, focusing on polarization effects caused by our input plasma of known properties.

Three cases are investigated, including (i) the control experiment where there is no bright point source lying along the line-of-sight, no radiation background, but the intervening medium is a self-emitting, absorbing, Faraday-rotating and Faraday-converting medium, (ii) the fiducial point source is placed at $z=z_{\rm init}=6.0$, serving as a bright distant radio back-light, and (iii) the fiducial point source is located much nearer, at $z=0.206$ (cf. Jagers et al. 1982). The prescription of the intervening plasma at $z=0$ follows model A-II described in Table 1; simple cosmological evolutions of $n_{\rm e}(z)$, $T_{\rm e}(z)$ and $|\textbf{{B}}(z)|$ described in Section 3.1.3 are now accounted for while the fraction of non-thermal relativistic electrons ${\mathcal{F}}_{\rm nt}$, their energy spectral index $p$ and the Lorentz factor of low-energy electron cut-off $\gamma_{i}$ are assumed to be constant over all redshifts. The results of the three different scenarios are displayed in parallel in Figs. 5 – 7 for comparison purposes. Numerical results are summarized in Table 3.

Differences in the results of the three cases indicate that on cosmological scale, polarized radiative transfer of light traveling through a foreground cosmologically-evolving IGM-like plasma, with or without a bright point source, can impart unique polarization features. Also, it can be readily seen from Figs. 5 and 6 that both the total emission and the polarized emission from the fiducial point source dominate over the contributions from the foreground plasma, as expected. The invariant intensity ${\mathcal{I}}$ of the radiation stays by and large constant from where the bright point source is positioned with a very small increase over increasing $z$ due to the emission of the line-of-sight plasma, which is calculated in case (i). Fluctuations in Stokes parameters are induced by random field orientations along the line-of-sight.

The observed change of polarization angle $\Delta\varphi$, which is a measure of the amount of Faraday rotation and is sensitive to the magnetic field directions along the line-of-sight, depends on the $z$-position of the point-source, as is seen in Fig. 7 and Table 3. In all three cases, we obtained $\Delta\varphi<\pi$. This indicates that the effect of Faraday rotation is weak, as is expected for a line-of-sight plasma that is threaded with a weak magnetic field of nG and has a low electron number density. Insignificant Faraday conversion is also observed in case (i), for which there is only the plasma but no bright sources lying along the line-of-sight. Note that ${\Pi}_{\rm c}$ is much weaker than ${\Pi}_{\rm l}$ by an order of magnitude of $10^{5}$. For case (ii), ${\Pi}_{\rm l}$ and ${\Pi}_{\rm c}$ are dominated by the contributions of the bright point source over the foreground plasmas. For case (iii), the sudden drops in ${\Pi}_{\rm l}$ and in ${\Pi}_{\rm tot}$, and the large rise in $\Delta\varphi$, and ${\Pi}_{\rm c}$ show the effects of having a foreground (nearby) source. Understanding the foregrounds, particularly any bright line-of-sight sources and their locations, is crucial for scientific inference of magnetic fields and their evolution.

In addition, depolarization effect is observed: there is a net drop in $\Pi_{\rm tot}$ as $z$ decreases (i.e. as path-length increases). By the experimental set up, this is mainly due to differential Faraday rotation (i.e. emission at different $z$ is rotated by different amount due to their magneto-ionized foreground, thus reducing the net polarization). Random magnetic fields has also been identified in the literature as another cause of depolarization (see e.g. Burn, 1966). Investigation of the effects of random fields is beyond the scope of this demonstration, but our results here illustrate how the effects on polarization can be quantified by performing a full cosmological polarized radiative transfer.

5.2 Single galaxy cluster

Here, we illustrate the making of intensity and polarization maps of an astrophysical object by carrying out pencil-beam (post-processing) CPRT calculations, where results obtained from a cosmological MHD simulation are incorporated. Each pixel of the maps corresponds to a solution obtained by the radiative transfer calculation.

We use the data of a simulated galaxy cluster obtained from the “cleaned" implementation of a higher resolution GCMHD+ simulation (see Section 4 in Barnes et al., 2018). The GCMHD+ simulations, designed to focus on the evolution of the magnetic field due to structure formation without the additional complications, are adiabatic, i.e. no radiative cooling, reionization, star formation and feedback from supernovae and Active Galactic Nuclei. The cluster obtained at $z=0$ from the simulation has a virial radius of $R_{\rm vir}=1.4439$ Mpc, and a gas mass of $m_{\rm gas}\sim 10^{13}$ ${\rm M}_{\odot}$. Non-thermal electrons has energy density that amounts to 1% of the thermal energy density (see Barnes et al., 2018). Simple statistics of the properties of the cluster are summarized in Table 4. In Fig. 8 we plot the central slices of the data cube viewing along the $z$-direction, illustrating the input structures of electron number density, magnetic field strength and orientation for the CPRT calculation.

Radiative transfer of a total number of $256^{2}=65536$ rays is computed from $z=6.0$ to $z=0.0$ through the galaxy cluster centered at $z_{\rm cluster}$. Without loss of generality, we choose $z_{\rm cluster}=0.5$ (i.e. placed between $z=0.500645$ and $z=0.499355$, corresponding to a length scale of $2.89$ Mpc $\approx 2R_{\rm vir}$). In order to study the intrinsic polarization emission of the cluster, no materials fill the line-of-sight outside the cluster and zero initial radiation background are assumed. Emission, absorption, Faraday rotation, and Faraday conversion by thermal bremsstrahlung and non-thermal synchrotron radiation process are taken into account.

Fig. 9 shows the resulting intensity and polarization maps obtained at $z=0$; simple statistics of those maps are summarized in Table 4. The simulated cluster is intrinsically polarized at the $\nu_{\rm obs}=1.42$ GHz with the mean value of degree of total polarization $\sim 68.57$ %, dominated by linear polarization. Emission is the highest in the cluster’s central region, where both magnetic field and electron number density are the highest (see Fig. 8). Faraday rotation is also strong in the central region, leading to a bigger change of polarization angle, as is seen in the map of $\Delta\varphi$ shown in Fig. 9. At the same time, depolarization in that region is also the most significant, where the degree of polarization is $\lesssim 30\%$ and the minimum reaches $\sim 1\%$. Strong differential Faraday rotation and the effect of random field orientations along the line-of-sight are the causes of depolarization in this demonstration. These results agrees with the observational trends of smaller degree of polarization for sources close to the cluster center (see e.g. Bonafede et al., 2011; Feretti et al., 2012).

CPRT calculation provides a rich set of data products, enables quantitative measures of polarization and intensity, and its algorithm allows interfacing with simulation results. While here we demonstrate the calculation of a simulated cluster at a fixed redshift and show only the intensity and polarization maps at $z=0$, the CPRT algorithm can generate maps at any sampled redshifts. Comparisons of the statistics of maps generated at different redshifts may provide a useful means to study the cosmological evolution of magnetic fields, as well as giving insights to tomographic studies of large-scale magnetic fields in real data. Mock data set obtained from CPRT calculations can also be used to test analysis tools used for magnetic field structure inference.

5.3 All-sky calculation

With the advent of the SKA surveys over a very large fraction of the celestial sky will be enabled. Here, we describe how applying the all-sky CPRT algorithm allows us to compute theoretical polarization maps of the radio sky, with a model magnetized universe obtained from a cosmological MHD simulation with GCMHD+ code (Barnes, Kawata & Wu, 2012; Barnes et al., 2018) as an input structure. Mock data of such a kind can be statistically characterized for comparison with observations, as well as serving as testbeds for validating analysis methods used for scientific inference.

Ray-tracing CPRT calculations are carried out for a total number of $N_{\rm ray}=12\times 64^{2}=49152$ rays distributed on $z$-spheres according to the HEALPix sampling scheme (Górski et al., 2005). We consider radiation frequency of $\nu_{\rm obs}=1.42$ GHz. Contribution from the redshifted CMB photons to the radiation background is neglected, and radio polarization is arisen from sources consisting both thermal and non-thermal electrons distributed across the entire universe in the post-reionization epoch141414We expect that non-linear growth in magnitudes and structures of electron number density during the reionization epoch would have imparted observational signatures to the traveling radiation, varying the statistics such as the polarization power spectrum. However, for demonstrative purpose we do not consider such an effect in this paper. (i.e. $\,z\leq 6.0$). Both thermal bremsstrahlung and non-thermal synchrotron radiation process are taken into account. To isolate the polarization signatures imparted by magnetic structures, electron number density $n_{\rm e,tot}(z,\theta,\phi)$ is assumed to be uniform across the entire sky at each $z$; its cosmological evolution over $z$ underwent a dilution in an expanding universe, i.e. $n_{\rm e,tot}(z)=n_{\rm e,tot,0}(1+z)^{3}$, where $n_{\rm e,tot,0}=2.1918\times 10^{-7}$ cm*-3* (see Appendix E for details). We assume that non-thermal relativistic electrons amounts to 1% of the total electron number density. We also assume that their energy spectrum follows a power law with a spectral index of $p=4.0$ (i.e. the non-thermal electrons have aged, steepening the spectrum), corresponding to a radiation power-law spectrum with index $\alpha=(p-1)/2=1.5$. The low cutoff of the electron energy is set to $\gamma_{i}=10.0$, and the high cut-off is set to infinity.

We use the GCMHD$+$ cosmological MHD simulation (Barnes, Kawata & Wu, 2012) to determine the evolution of the large-scale magnetic field as structures in the universe assemble. A cubic region of comoving volume ($40$ Mpc)3 was taken from a comoving ($100$ Mpc)3 volume in the simulation, which started at $z=47.4$ as determined by the initial condition generator grafic$++$. The magnetic field was assumed to be generated at some early epoch via a method that filled the volume of the simulation. It has a configuration of $B=(10^{-11},0,0)$ ${\rm G}$. We fit analytically the output of $B_{\parallel}(z)$ obtained from the GCMHD$+$ simulation by the piecewise function:

[TABLE]

with $x=\log_{10}{z}$. This fit151515Note that the anomalous bump in Fig. 10 at $\log{z}\approx-0.5$ is caused by the instantaneous infall and outflow of the simulation box. This structure does not appear in the other four simulations that ran with different initial conditions, and is therefore neglected., plotted in Fig. 10, is smoothed by interpolation using twenty-one-points averages to model the input of $B_{\parallel}(z)$ for the CPRT calculation. We assumed a log-normal spatial distribution of $B_{\parallel}(z,\theta,\phi)$ over each $z$-sphere, where the mean value is deduced from equation (64) multiplied by a factor of $10^{3}$ to match the expected observed field strength of $1.0$ nG typical to filaments (see e.g. Araya-Melo et al., 2012). The log-normal distribution ensures the magnetic field strength to be all positive. Directions of the magnetic fields, which are defined by the $\cos{\theta}$, are assumed to have random orientations along the line-of-sight.

5.3.1 Results and discussion

**(I) Along a randomly selected ray

**Fig. 11 shows the resulting cosmological evolution of both the invariant and co-moving Stokes parameters, as well as the cosmological evolution of $\Delta\varphi$, $\Pi_{\rm l}$, $\Pi_{\rm c}$ and $\Pi_{\rm tot}$ of a randomly selected ray. Notably, one can see that the fluctuations in $Q,U$ and $V$ increase significantly during the late time, i.e. when the structure formation and evolution processes (such as the assembly of galaxy clusters) in the cosmological simulation become prominent and that magnetic fields become significantly amplified along with these processes, hence imposing a Faraday screen (i.e. strong Faraday-rotating component). In addition, highly volatile behavior is observed in the change of polarization angle over $z$, i.e. throughout the entire radiation path. Volatility in the evolution of polarization angle increases the difficulty to distinguish between different Faraday depth components, limiting the usage of the standard approach to infer magnetic field properties using RM synthesis (see e.g. Brentjens & de Bruyn, 2005) in some cases. Results showing similar trends in polarization evolution are observed commonly in all the other randomly selected rays.

(II) All-sky maps

Theoretical all-sky polarization maps of $I$, $Q$, $U$, and $V$ can be generated at any chosen redshifts. In Fig. 12 we show the Stokes maps obtained at $z=0$. Their statistics are summarized in Table 5. As pointed out from the previous discussion on the single-ray results, the evolution of the change of polarization angle, which serves as a probe to Faraday rotation effect, is highly volatile and complex, demanding advanced statistical analyses at different redshifts to be performed for science extraction from map data.

(III) Remarks on the CPRT method and the existing RM techniques

Our demonstrative calculation results have two major implications in the study of large-scale magnetic fields, firstly on the future power spectrum analysis, and secondly on the validity of the current methodologies to investigate large-scale magnetic fields.

Our results show that a Faraday screen can be introduced when structure formation and evolution processes in the universe becomes prominent, as is seen in Fig. 11 where significant polarization fluctuations happened during the late time when galaxy clusters started to assemble in the simulation, boosting the mean magnetic field strength. This finding means that cosmological contributions from line-of-sight IGM-like media will likely be screened (or shielded) by fluctuations sourced from astrophysical structures like a galaxy cluster (i.e. ionized systems with relatively high magnetic field strengths and electron number density). This further implies that the polarization power spectrum of an all-sky map will be dominated by high frequency (small scale) signals. At the same time, it is worth noting that the morphology of ionized bubbles during the Epoch of Reionization, which has not been investigated in this paper, may imprint observable signatures onto the polarization maps, contributing to the power in low frequency (large scale) in polarization power spectrum as those ionized regions overlapped.

The highly volatile cosmological evolution of Stokes parameters suggests that analysis methods using RM is likely to be deemed inappropriate to study inter-galactic magnetic fields, particularly those that permeate emitting cosmic filaments. This is because Faraday rotation will no longer be the single important process that imprints the polarization signals, but also the emissions from filaments themselves, as well as the absorption processes along the line-of-sight. However, the quantity of RM is derived from a restrictive case of polarized radiative transfer, as is described in Section 1 and reviewed in details in our related paper (On et al., 2019). Interpretation of polarization signals in those cases, therefore, requires full CPRT consideration, so to correctly and accurately determine how large-scale magnetic fields have evolved and where they came from.

6 Discussion and Summary

In this paper, a covariant formulation of cosmological polarized radiative transfer, which provides a solid theoretical foundation to use polarized light as a probe of large-scale magnetic fields, is presented. Such a formulation naturally accounts for the space–time metric of an arbitrary cosmological model with a flat geometry. It is derived based on a covariant general relativistic radiative transfer formulation, which is derived from the first principles of conservation of phase–space volume and photon number. Without loss of generality, the corresponding polarized radiative transfer equations derived using a flat FRW space–time metric are constructed. In addition, we developed the (all-sky) CPRT algorithm that allows incorporation of the results from cosmological MHD simulation to the CPRT calculations, henceforth, a straightforward generation of theoretical polarization maps. Those maps serve as model templates, crucial for interpreting all-sky polarized data which will be measured by the next generation radio telescopes such as the SKA.

Sets of CPRT calculations are performed to validate the code implementation of the ray-tracing algorithm and to demonstrate its applications for practical astrophysical studies. We summarize below the richness of the polarization data product the CPRT algorithm offers, as well as the findings from our demonstrative sets of calculations.

Solving the CPRT equation yields the evolution of the Stokes parameters of radiation as a function of $z$, allowing tracking of how the intensity and polarization of radiation are modified on its way by local radiation processes (thermal bremsstrahlung and non-thermal synchrotron radiation process in this paper) in a cosmologically evolving universe. From the set of single-ray calculations presented in Section 5.1, we showed the resulting evolution of the polarization for cases where a bright radio point source is present or absent, and magnetic field orientations are random along the line-of-sight. It is seen that line-of-sight bright radio point sources dominate the intensity and polarization, and their locations lead to different signatures in the polarization evolution of the radiation. CPRT calculations provide quantitative studies of the intensity and polarization of radiation in their transport. They allow direct tracking of the change of polarization angle, which will help to resolve $n\pi$-ambiguity problem, aiding our interpretation of observational data. Evolution of the degree of linear, circular and total polarization can also be computed, allowing further investigation of Faraday rotation and depolarization.

Carrying out multiple-ray CPRT calculation yields data maps of intensity and polarization, where spatial fluctuations across the sky plane can be statistically studied and characterized for comparison with observational data for magnetic field structure inference. We performed such a calculation using simulated cluster data obtained from a GCMHD+ simulation as the input. Contributions from the galaxy cluster dominated over those from the inter-galactic space. This is as expected due to their much higher electron number density and magnetic field strength. Faraday screening effect may be dominated when performing common analysis methods that use RM as a quantitative measure to study intra-cluster magnetic fields. We highlight that carrying out a full CPRT calculation, which does not assume the relative strength of radiative transfer effect in emission, absorption, Faraday rotation, and Faraday conversion, allows a reliable assessment of the validity of the standard RM methods in different astrophysical scenarios in an expanding Universe.

In full cosmological settings, we performed an all-sky CPRT calculation using a model magnetized universe obtained from a cosmological GCMHD+ simulation as the input. Our results show that the cosmological evolution of the polarization components of propagating radiation is highly volatile, suggesting that full CPRT consideration is needed for accurate large-scale inter-galactic magnetic field studies, particularly for the fields that permeate emitting cosmic filaments. Another implication is that polarization power spectra obtained from all-sky measurements are likely to be dominated by the high frequency (small scale) signals caused by strong Faraday-rotating components, such as galaxy clusters. Impacts on the polarization signals due to the morphology of the cosmic reionization, which are not addressed in this paper but are important research problems, will be considered in our future work.

All in all, the CPRT formulation provides a reliable platform to compute polarized sky. Furthermore, with known input distributions of $n_{\rm e}(z)$ and $\textbf{{B}}(z)$, and full radiative transfer processes taken into account, results obtained from the forward computation of the CPRT algorithm will provide valuable data sets that may also serve as a testbed for assessing analysis tools used for large-scale magnetic field studies. Also, since the cosmological terms in the CPRT equation can be easily switched off in our algorithm and its code implementation, calculations in astrophysical contexts can be easily carried out; calculations for foreground contributions in cosmology studies can also be performed.

With the current version of the implementation, our next step is to characterize the polarization fluctuations for different input magnetic field and electron number density distributions, developing statistical methods for reliable scientific inference from data. Alongside, implementation of more accurate transfer coefficient expressions for a broader class of synchrotron distributions, as well as the inclusion of cyclotron process, is to be carried out and tested. Solving a stiff set of CPRT equations is foreseen to be one of the biggest numerical challenges. Nonetheless, the CPRT formulation, and its algorithm provide a solid theoretical foundation and a reliable platform to study large-scale magnetic fields. The CPRT formulation derived and the (all-sky) algorithm that has been developed enable more straightforward comparisons between theories and observations, ultimately guiding us to answers about the origins and the co-evolution of magnetic fields with structures in the Universe.

Acknowledgments

We thank Ziri Younsi for discussions of covariant radiative transfer. Jennifer Y. H. Chan is supported by UCL Graduate Research Scholarship and UCL Overseas Research Scholarship. Alvina Y. L. O is supported by the Government of Brunei under the Ministry of Education Scholarship. Thomas D. Kitching acknowledges support from a Royal Society University Research Fellowship.

Appendix A Adopted coordinate system

In this work, we adopt right-handed coordinate systems as depicted in Fig. 13 following Huang & Shcherbakov (2011a). In our notations, magnetic field B is directed along the $\tilde{z}$-axis, making an angle $\theta$ clockwise to the propagation direction of the radiation k. An orthonormal $(x,y,z)$ basis is defined such that $\textbf{{z}}\parallel\textbf{{k}}$, $\textbf{{x}}=C(\textbf{{B}}\times\textbf{{k}})$, where $C$ is a scalar that can be positive or negative, and $\textbf{{x}}\parallel\tilde{\textbf{{x}}}$, and $\textbf{{y}}=(\textbf{{k}}\times\textbf{{x}})$. Here, x is perpendicular to the plane of $(\textbf{{B}},\textbf{{k}})$, and $(\textbf{{B}},\textbf{{k}},\textbf{{y}})$ are coplanar. Electric field of an electromagnetic wave traveling along $\textbf{{k}}\parallel\textbf{{z}}$ oscillates in the $(x,y)$-plane. By such a choice of configuration (or by the choice of $\textbf{{y}}\parallel\tilde{\textbf{{y}}}$ in the systems defined in Sazonov (1969); Pacholczyk (1977)), absorption coefficient $u_{\nu}$, conversion coefficient $g_{\nu}$ and emission coefficient $\epsilon_{U,\nu}$ are zeros.

Note that the transfer coefficient matrices are commonly derived in the “magnetic-field" system, i.e. first in the $(\tilde{x},\tilde{y},\tilde{z})$ basis, and then projecting them onto $(x,y)$ for $\textbf{{k}}\parallel\textbf{{z}}$ and $\cos{\theta}=(\textbf{{k}}\cdot\textbf{{B}})/(|\textbf{{k}}||\textbf{{B}}|)$ (see e.g. Sazonov, 1969; Pacholczyk, 1970, 1977; Jones & Odell, 1977a; Huang & Shcherbakov, 2011a). Transformation between the coordinate systems $\tilde{e}_{i}=(\tilde{x},\tilde{y},\tilde{z})$ and ${e}_{j}=(x,y,z)$ is given by ${e}_{j}=\tilde{e}_{i}{M}_{ij}$, where

[TABLE]

i.e.

[TABLE]

It follows that the rotation of vectors is given by $A_{i}=({M}^{\rm T})_{ij}\tilde{A}_{j}$, and the rotation of tensors is given by $\sigma_{ij}=({M}^{\rm T})_{ik}\tilde{\sigma}_{km}({M})_{mj}$ (Huang & Shcherbakov, 2011a). In future studies where observational data are confronted with theoretical predictions obtained by CPRT calculations, it is also useful to introduce the “observer’s" (or polarimeter’s) system $(a,b)$, which is defined by rotating the $(x,y)$-plane about the $k$-direction. Such a transformation, i.e. between the local system (given by the local projection of the magnetic field) in the comoving frame and the frame in which polarimetric data are measured, invokes the use of rotational matrix ${\bf{\mathcal{R}}}(\chi)$, which follows the definition given in equations (50) and (51) in Huang & Shcherbakov (2011a), where the angle $\chi$ relates a and b to the magnetic field components perpendicular to k, i.e. $\textbf{{B}}_{\perp}=\textbf{{B}}-\textbf{{k}}(\textbf{{k}}\cdot\textbf{{B}})/k^{2}$, by $\sin{\chi}=(\textbf{{a}}\cdot\textbf{{B}}_{\perp})/|\textbf{{B}}_{\perp}|$ and $\cos{\chi}=-(\textbf{{b}}\cdot\textbf{{B}}_{\perp})/|\textbf{{B}}_{\perp}|$ respectively.

Appendix B Convention of polarization

Stokes parameters $I_{\nu}$, $Q_{\nu}$, $U_{\nu}$ are defined unambiguously once the $(x,y)$ coordinate system is specified. The different definitions of polarization angle adopted in the cosmic microwave background community and the International Astronomical Union (IAU) can be reconciled by a sign flip of $U_{\nu}$. However, interpretation of the sign of $V_{\nu}$ (and consequently the signs for the corresponding transfer coefficients $\epsilon_{V,\nu}$, $v_{\nu}$ and $h_{\nu}$) in the literature is often ambiguous. This is because the sign of $V_{\nu}$ depends not only on the definition of the senses of circular polarization (which also depends on the handedness of the coordinate systems used) and the definition of $V_{\nu}$, but also on the choice of sign in the time-dependent description of the electromagnetic wave, as well as the definition of the relative phase between the $x$ and $y$-components of the electric vector of the wave. Much variation in these dependences exist in the literature, or sometimes this information is inexplicitly assumed or left unstated. Another source of variation comes from the choice of the attachment of the sense of circular polarization to the helicity of the photon. Any confusion and ambiguity can easily cause a slip in the interpretation of $V_{\nu}$.

Here, we first describe the circular polarization sense defined by the Institute of Electrical and Electronics Engineers (IEEE) (IEEE, 1998), which is commonly adopted by radio astronomers (but opposite to classical physicists and optical astronomers’ common practice161616The right-handed circular polarization convention by the IEEE corresponds to the left-handed circular polarization convention in the classical sense, i.e. IEEE-RCP $=$ classical-LCP.), and the International Astronomical Union (IAU) convention of Stokes $V_{\nu}$ (Reid, 2007). Then we discuss the intricacies to test the conformity to the IEEE/IAU polarization convention. Finally we remark on the magnetic field direction of the system and state explicitly the Stokes V convention used in this paper.

IEEE/IAU polarization convention:

The exact quote of the IEEE (1998)’s definition171717The same definition was first introduced in 1942 when the IEEE was still known as the Institute of Radio Engineers (IRE). of a right-handed polarized wave reads “a circularly or an elliptically polarized electromagnetic wave for which the electric field vector, when viewed with the wave approaching the observer, rotates counter-clockwise in space". As pointed out by Hamaker & Bregman (1996), such a definition stipulates that the position angle $\varphi$ of the electric vector of the wave at any point increases with time, implying that the $y$-component of the filed, $\textbf{{E}}_{y}$, to lag the $x$- component, $\textbf{{E}}_{x}$. In other words, the electric field traces out a counter-clockwise helix (right-hand screw) in time at fixed position, whereas in space at any instant in time it forms a clock-wise helix (left-hand screw) (see e.g. Rochford, 2001). The IAU endorses the sense of circular polarization defined by IEEE and defines $V_{\nu}={\rm(RCP-LCP)}$, i.e. $V_{\nu}$ is positive for RCP (Reid, 2007). The $x$- and $y$- axes of a right-hand triad align with North and astronomical East, and the $z$- axis points towards the observer for standard IAU convention.

Conformity to IEEE/IAU convention:

It is important to note that even when the Stokes parameters are uniquely defined when combined the IEEE/IAU definition with the standard formulae for $I_{\nu}=\langle|E_{x}|^{2}+|E_{y}|^{2}\rangle$, $Q_{\nu}=\langle|E_{x}|^{2}-|E_{y}|^{2}\rangle$, and $U_{\nu}=2\langle|E_{x}|\,|E_{y}|\,\cos{\delta}\rangle$, two similar but distinct mathematical representations are allowed for the same physics of the problem, as is shown by Hamaker & Bregman (1996). One has the choice181818 Another choice is related to the attachment of the RCP and LCP to positive and negative helicity (see also Appendix III in Simmons & Guttmann, 1970, for a complete table of different conventions of RCP, including those that do not comply to IEEE/IAU convention).

of the sign of the time dependence of the electromagnetic wave, i.e. $e^{+i\omega t}$ or $e^{-i\omega t}$, for $\omega>0$. Both choices are equally valid, but once the sign is chosen for

[TABLE]

the following quantities must have the sign adjustments such that $\textbf{{E}}_{y}$ lags $\textbf{{E}}_{x}$ for a (unit amplitude) RCP wave:

[TABLE]

Hamaker & Bregman (1996), and that

[TABLE]

Hamaker & Bregman (1996), so that $V$ is positive for RCP, i.e. IEEE/IAU compliant. Note that the sign adjustment in equation (88) is equivalent to defining the sign of $\delta=\pm(\phi_{y}-\phi_{x})$ in equation (87) for $\delta\in(0,\pi)$, where time delays correspond to negative (positive) values of the phases $\phi_{x}$ and $\phi_{y}$ for $e^{\pm i(\omega t-kz)}$ according to Equation (83). It is apparent that one differing convention of the above would lead to a sign reversal. In brief, an unambiguous interpretation of the circular polarization from $V_{\nu}$ requires a clear specification of the adopted handedness of the coordinate systems, the convention of circular polarization, the definition of Stokes $V_{\nu}$, as well as the chosen mathematical representation of the traveling plane wave.

Remark on the B-field convention:

Given the coordinate systems and the geometry of the problem presented in Fig. 13, let’s consider the simple case where a uniform magnetic field B aligns with k, so ${\theta}=0$. An electron would then precess about B in the $(\tilde{x}-\tilde{y})$-plane, moving counter-clockwise as viewed along $\textbf{{k}}\parallel\textbf{{B}}$. The electric vector of the electromagnetic wave follows the electron motion, thus also rotating couter-clockwise as viewed by the observer. This results in IEEE-RCP, and according to the IAU convention, $V_{\nu}>0$.

In this paper, we adopt the conventions conforming to the IEEE/IAU standard and stick to the magnetic field convention where the magnetic field is positive when pointing towards the observer191919This is opposite to the astronomical convention that traditionally defines magnetic field direction as positive when pointing away from the observer (i.e. $\theta=0$ corresponds to a negative field while $\theta=\pi$ corresponds to a positive field).. We follow the same coordinate systems as Huang & Shcherbakov (2011a) and use it as the main reference paper to check against the signs of the Stokes parameters and their corresponding transfer coefficients. The transfer coefficients therefore all have positive signs in their expressions.

Appendix C Transfer coefficients

In this Appendix we present the transfer coefficients for both thermal bremsstrahlung and non-thermal synchrotron radiation process. The non-thermal relativistic electrons gyrating around magnetics field lines has a power-law energy spectrum. We adopt the expressions given in Pacholczyk (1977) and Jones & Odell (1977a) respectively, but the sign of the circular polarization described by Stokes $V$ are made to be consistent and complied to the IEEE/IAU convention, given the coordinate system explicitly shown in Appendix A. The emission coefficients have units of ${\rm erg}\,{\rm s}^{-1}\,{\rm cm}^{-3}\,{\rm Hz}^{-1}\,{\rm str}^{-1}$ and the absorption and Faraday coefficients have units of ${\rm cm}^{-1}$.

C.1 Thermal bremsstrahlung

Transfer coefficients of thermal bremsstrahlung have been presented in Pacholczyk (1977); Meggitt & Wickramasinghe (1982); Wickramasinghe & Meggitt (1985); Rybicki & Lightman (1986). In this paper, we adopt the expressions given in Pacholczyk (1977) and make certain changes such that the set of coefficients would follow the same conventions of polarization we have specified.

For a magnetized thermal plasma, the coefficients of Faraday rotation and Faraday conversion are respectively,

[TABLE]

(Pacholczyk, 1977), where $\omega=2\pi\nu$ is the radiation angular frequency, $\omega_{\rm p}=(4\pi n_{\rm e}e^{2}/m_{\rm e})^{{1}/{2}}$ is the plasma frequency, $\omega_{\rm B}=(eB/m_{\rm e}c)$ is the electron gyrofrequency, and $\theta$ is the angle between the radiation propagation and the magnetic field. The thermal bremsstrahlung components of the absorption coefficients are given by

[TABLE]

(Pacholczyk, 1977), where the collisional frequency is

[TABLE]

with the Coulomb logarithm factor, for $\omega\gg\omega_{\rm p}$,

[TABLE]

(Lang, 1974). $k_{\rm B}$ is the Boltzmann constant and $T_{\rm e}$ is the temperature of the electrons in thermal equilibrium. The emission coefficients in $I$, $Q$ and $V$ can be computed via the Kirchoff’s law:

[TABLE]

where the Planck function $B_{\omega}=k_{\rm B}T_{\rm e}\omega^{2}/(2{\pi}^{2}c^{2})$ by the Rayleigh-Jeans law.

It is interesting to note that both the frequency dependence and the dependence on the magnetic field are different for Faraday rotation and Faraday conversion. The strength of the Faraday rotation effect is proportional to $\nu^{-2}n_{\rm e,th}|\textbf{{B}}_{\parallel}|\delta s$, and the strength of Faraday conversion is proportional to $\nu^{-3}n_{\rm e,th}|\textbf{{B}}_{\perp}|^{2}\delta s$, where $|\textbf{{B}}_{\parallel}|=|\textbf{{B}}|\cos{\theta}$, $|\textbf{{B}}_{\perp}|=|\textbf{{B}}|\sin{\theta}$, and $\delta s$ is the photon propagation length.

Another useful remark concerns the use of rotation measure (RM) in the literature for quantifying the strength of Faraday rotation. RM is defined as ${\mathcal{R}}\equiv\Delta\varphi\,c^{2}/\nu^{2}$, where $\varphi=0.5\arctan({U/Q})$. A widely-used formula in RM analysis is ${\mathcal{R}}(s)=0.812\int^{s}_{s_{0}}\frac{{\rm d}s^{\prime}}{\rm pc}\,\frac{n_{e}(s^{\prime})}{\rm cm^{-3}}\,\frac{|\textbf{{B}}_{\parallel}(s^{\prime})|}{\mu{\rm G}}$ ${\rm rad}\,{\rm m}^{-2}$, which it can be shown that this is derived from the polarized radiative transfer equation (equation 18), under the assumptions that the effects of emission, absorption, Faraday conversion and contribution from non-thermal electrons are negligible (see On et al., 2019, for details). In a realistic situation, however, these assumptions do not hold. The intensity of $Q$ and $U$ of the observed polarized light is not solely dictated by Faraday rotation process. An accurate inference of magnetic field properties from the polarization signatures of observed light, therefore, demands a full polarized radiative transfer treatment.

C.2 Non-thermal synchrotron radiation

We adopt the expressions of the transfer coefficients for cosmic synchrotron sources from Jones & Odell (1977a), and make appropriate sign changes for the transfer coefficients at $V_{\nu}$ to keep a self-consistent polarization convention defined explicitly in this paper. For relativistic electrons following a power-law energy distribution with an index $p$ ,

[TABLE]

(Jones & Odell, 1977a), where $\Theta(\gamma-\gamma_{i})$ is the step function, $\gamma_{i}$ is the low-energy cutoff of electrons, and $g(\Psi)$ is the pitch-angle distribution, normalized to $\int\operatorname{d}\!{\Omega}_{\Psi}g(\Psi)=1$. The corresponding number density of electron is

[TABLE]

The normalization factor $[n_{\gamma}\gamma^{p}]$ and the index $p$ are related to the spectral index of the radiation by $\alpha=(p-1)/2$. The transfer coefficients for non-thermal synchrotron radiation are

[TABLE]

(Jones & Odell, 1977a, and references therein), where $\omega_{B_{\perp}}=\omega_{\rm B}\,\sin\theta$, $\kappa_{\perp}=(2\pi r_{\rm e}c)\,{\omega^{\,\,\,-1}_{B_{\perp}}}[4\pi g(\theta)][n_{\gamma}\gamma^{p}]$, $\epsilon_{\perp}=(m_{\rm e}c^{2})(r_{\rm e}/2\pi c)\,{\omega_{B_{\perp}}}[4\pi g(\theta)][n_{\gamma}\gamma^{p}]$ with the classical electron radius $r_{\rm e}=e^{2}/{m_{\rm e}c^{2}}$, and the fiducial frequency $\omega_{i}=\gamma_{i}^{2}\omega_{B_{\perp}}$. The dimensionless functions in the transfer coefficients are

[TABLE]

The transfer coefficients are derived from a nearly isotropic dielectric tensor, appropriate for cosmic plasmas with low electron densities and weak magnetic fields, such that $\omega>\omega_{i}$ and both $\omega$ and $\omega_{i}$ are above the gyro-frequency $\omega_{\rm B}$. The condition $\gamma_{i}^{2}>\cot^{2}\theta$ also has to be satisfied. In addition, dielectric suppression is assumed to be negligible, which generally holds valid for cosmic media (see Jones, O’dell & Stein, 1974; Melrose & McPhedran, 1991, for details). In this paper isotropic electron distribution is assumed so $g(\theta)=1/4\pi$. Comparing to the thermal bremsstrahlung expression in the high-frequency limit $(\omega\gg\omega_{\rm B})$, the non-thermal synchrotron Faraday rotation coefficient has an extra function factor \zeta(p,\gamma_{i})=\frac{(p-1)(p+2)}{(p+1)}\Big{(}\frac{\ln{\gamma_{\rm i}}}{\gamma_{\rm i}^{2}}\Big{)}, implying that Faraday rotation weakens with increasing electron energy (see also Melrose, 1997; Huang & Shcherbakov, 2011a).

Appendix D Derivation of the covariant radiative transfer formulation

Derivation of the covariant radiative transfer formulation has been presented in Rybicki & Lightman (1986); Fuerst & Wu (2004); Younsi, Wu & Fuerst (2012). Here, we repeat the derivation for clarity and completeness.

Consider a bundle of particles filling a phase–space volume element ${\rm d}\mathcal{V}\equiv\operatorname{d}\!{\bf x}^{3}\,\operatorname{d}\!{\textbf{{p}}}^{3}$, with 3-spatial volume element $\operatorname{d}\!{\bf x}^{3}=\operatorname{d}\!{x}\,\operatorname{d}\!{y}\,\operatorname{d}\!{z}$ and the 3-momentum volume element $\operatorname{d}\!{\textbf{{p}}}^{3}=\operatorname{d}\!{p}_{x}\,\operatorname{d}\!{p}_{y}\,\operatorname{d}\!{p}_{z}$ at given time $t$. According to the Liouville’s theorem, ${{\rm d}\mathcal{V}}/{{\rm d}\lambda_{\rm a}}=0$. Since ${\rm d}\mathcal{V}$ is conserved along the affine parameter $\lambda_{\rm a}$, it is Lorentz invariant.

The distribution function (or phase space density) of the particles in the bundle is represented by $f(x^{i},\,p^{i})=\operatorname{d}\!{N}/{\rm d}\mathcal{V}$, where $\operatorname{d}\!{N}$ is the number of particles in $\operatorname{d}\!{V}$. Since $\operatorname{d}\!{N}/\operatorname{d}\!{\mathcal{V}}$ is Lorentz invariant, $f(x^{i},\,p^{i})$ is also Lorentz invariant.

For photons, $v=c$ and $cp=E$, where $E$ is the photon energy. The spatial and momentum volume elements are $\operatorname{d}\!{\bf x}^{3}=\operatorname{d}\!{A}\,c\,\operatorname{d}\!{t}$ and $\operatorname{d}\!{\textbf{{p}}}^{3}=E^{2}\operatorname{d}\!{E}\,\operatorname{d}\!{\Omega}$, where $\operatorname{d}\!{A}$ is the area element through which the photons travel in the time interval $\operatorname{d}\!{t}$ and $\operatorname{d}\!{\Omega}$ corresponds to the direction of photon propagation. It follows that

[TABLE]

(see Rybicki & Lightman, 1986). The specific intensity of the radiation is

[TABLE]

Comparing the two expressions yields

[TABLE]

where $\mathcal{I}_{\nu}$ is the Lorentz-invariant intensity.

The Lorentz-invariant absorption and emission coefficients are $\zeta_{\nu}=\nu\,\kappa_{\nu}$ and $\xi_{\nu}={\epsilon_{\nu}}/\nu^{2}$, respectively (Rybicki & Lightman, 1986). It follows that the covariant radiative transfer equation takes the form

[TABLE]

where the source function ${\mathcal{S}}_{\nu}\equiv\xi_{\nu}/\zeta_{\nu}={\epsilon}_{\nu}/({\kappa_{\nu}\,\nu^{3}})$. Since $\zeta_{\nu}$ and $\xi_{\nu}$ are invariants under the Lorentz transformation, the transfer coefficients measured in the observer’s frame relates to those in the co-moving frame (i.e. the local rest frame of the medium) via $\nu\,\kappa_{\nu}=\nu_{{\rm co}}\,\kappa_{\nu,{\rm co}}$ and ${\epsilon_{\nu}}/\nu^{2}={\epsilon_{\nu,{\rm co}}}/\nu_{{\rm co}}^{2}$. Hence, the radiative transfer equation becomes

[TABLE]

(Fuerst & Wu, 2004; Younsi, Wu & Fuerst, 2012).

Appendix E Calculation of the total electron number density at the present epoch

The Universe is neutral as a whole and the most common atoms in it are Hydrogen and Helium. We can approximate $n_{\rm e,tot}=n_{\rm p,tot}=n_{\rm p,He}+n_{\rm p,H}$, where “p" stands for proton, “H" for Hydrogen and “He" for Helium; $n_{\rm p,He}\approx\rho_{\rm He}/m_{\rm He}$, and $n_{\rm p,H}\approx\rho_{\rm H}/m_{\rm H}$. By approximating the density of Hydrogen taking up 75 % of the density of baryons (i.e.  $\rho_{\rm H}=3\rho_{\rm b}/4$), and the density of Helium taking up the remainder, it gives $n_{\rm e}=7\rho_{\rm b}/8m_{\rm p}$. The value of $\rho_{\rm b,0}$ can be calculated from $\Omega_{\rm b,0}=\rho_{\rm b,0}/\rho_{\rm crit}$, with $\Omega_{\rm b,0}h^{2}=0.02230$ (Planck Collaboration et al., 2016a), and $\rho_{\rm crit}=3H_{0}/(8\pi G)=1.87882\times 10^{-29}h^{2}$. This gives $n_{\rm e,0}=2.1918\times 10^{-7}$ cm*-3*.

Appendix F Remarks on finding an appropriate scale length

Here, in Table 6, we present the numerical values of the absorption, emission and Faraday rotation coefficients used in the calculations presented in Section 4.1. In general, the very different properties of cosmic media lead to a wide range of orders of magnitude spanned by transfer coefficients in the CPRT equation, resulting in a stiff set of coupled differential equations to solve. It is therefore essential and important to test the capability of the equation solver and the stability of the numerical solution (see Section 4.1). We emphasize that finding an appropriate scale length is crucial to overcoming the stiffness issue. In this work, the very small order of magnitude of the transfer coefficients computed using parameters typical to an IGM and an ICM at $\nu_{\rm obs}=1.42$ GHz suggests a scale length of a few Mpc when determining the $z$-sampling scheme.

In addition, note that all the CPRT calculations for the situations discussed in this paper are optically thin (i.e. $\tau\ll 1$). While the media are optically thin, they can be Faraday thick at the same time, such as in the cases of ICM-like environments. Numerical values of the optical depths and Faraday conversion coefficients obtained using the IGM-like model A-I and the ICM-like model B-I are included in Table 6. Note also that the effect of Faraday conversion is usually much weaker than that of Faraday rotation. Hence, $V$ is nearly always zero in the cases of our interests.