U_{0}(\mu)=\left\\{\begin{array}[]{lrl}\displaystyle\frac{6}{3+2b}(\mu+b\mu^{2}),&0&<\mu<1;\\\
\displaystyle 0,&-1&<\mu<0,\\\ \end{array}\right.\\\ \end{split}$ (146)

where $b$ is a constant and set to be $2$ in the practical calculation.
Obviously $U_{0}(\mu)$ gives rise to a limb-darkening law for the incident
radiations.

Our results is shown in Figs. 17-20, corresponding to Figs. 2-6 of Haardt
(1993). Figs. 17-19 demonstrate the total Comptonized spectra together with
the spectra of each single scattering. The spectra vary with different optical
depth, plasma temperature and observational angle. Fig. 20 shows how the total
spectrum changes with the viewing angles, especially for $\mu<0$,
corresponding to the spectrum reflected back from the atmosphere. One can see
that our results are consistent with that of Haardt (1993).

Figure 17.— The Comptonized spectra of low-energy thermal radiations by a hot
plasma atmosphere. The total spectrum (solid black line) and each single
scattering contributions are shown together. The dimensionless temperature and
optical depth for the hot plasma are $\Theta=1.7$ and $\tau=0.15$. The cosines
of the viewing angles are $\mu=0.1$, $0.5$ and $0.9$ for panels from left to
right, respectively. The observed photon energies are in unit of
$kT_{\text{bb}}$, which is the temperature of the incident thermal radiation.
Compare to Fig. 2 of Haardt (1993). Figure 18.— Same as in Figure 17, but for
$\Theta=1.59$ and $\tau=0.01$. Compare to Fig. 3 of Haardt (1993). Figure 19.—
Same as in Figure 17, but for $\Theta=0.188$ and $\tau=0.06$. Compare to Fig.
4 of Haardt (1993). Figure 20.— The inverse Comptonized spectra of low-energy
thermal radiations by a hot plasma atmosphere for various viewing angles. The
parameters for the hot plasma are $\Theta=0.5$, $\tau=0.1$ and $\Theta=1.0$,
$\tau=0.1$ for left and right panels, respectively. Compare to Figs. 5 and 6
of Haardt (1993).

###  4.8. The First Order Green’s Function For Relativistic Compton Reflection

Figure 21.— The first-order Green’s function for monoenergetic incident
photons with 511 keV reflected from a semi-infinite atmosphere of cold
electron. The curves with different colors correspond to different values of
the cosines of viewing angle $\mu_{2}$. The solid and dotted lines represent
the results calculated by the analytical formulae of Hua & Lingenfelter (1992)
and our MC scheme, respectively. Compare to Fig. 1 in Hua & Lingenfelter
(1992).

With the above discussed strategy, we will reproduce the calculations of Hua &
Lingenfelter (1992) for the first order Green’s functions of relativistic
Compton reflection by a semi-infinite plane-parallel atmosphere of cold
electrons. Hua & Lingenfelter (1992) have shown that the first order Green’s
function for such a system has an analytical expression and is strongly
dependent on the viewing angles of the atmosphere. Since the cold electrons
can be approximately regarded as static, we do not have to consider any
electron sampling process and Lorentz transformation at all. The unpolarized
and monoenergetic (511 keV) photons are isotropically injected from all
directions to the surface of the atmosphere. The Green’s function can be
expressed as an expansion of functions characterized by the number of
scattering times (Hua & Lingenfelter, 1992)

$\displaystyle\displaystyle G(E,\mu_{2})=\sum_{n=1}^{\infty}g_{n}(E,\mu_{2}),$
(147)

where $E$ is the escaping energy of the photon, $\mu_{2}$ is the cosine of the
viewing polar angle $\theta_{2}$. Using the the MC method, one can easily
obtain the numerical result for any function $g_{n}$ in this series. In our
scheme, we simply have

$\displaystyle\begin{split}&g_{n}(E,\mu_{2})=\displaystyle\int\cdots\int
S(P_{0})K(P_{0}\rightarrow P_{1})\cdots\\\ &\times K(P_{n-1}\rightarrow
P_{n})f(P_{n})\delta(\mu-\mu_{2})dP_{0}dP_{1}\cdots dP_{n}.\end{split}$ (148)

Here we mainly care about the first order $g_{1}(E,\mu_{2})$, i.e.,

$\displaystyle\begin{split}&g_{1}(E,\mu_{2})=\displaystyle\iint
S(P_{0})K(P_{0}\rightarrow P_{1})\\\ &\times
f(P_{1})\delta(\mu-\mu_{2})dP_{0}dP_{1},\end{split}$ (149)

which has been analytically integrated by Hua & Lingenfelter (1992).

For a given viewing angle, there is a cut-off energy (Hua & Lingenfelter,
1992)

$\displaystyle\displaystyle
E_{\text{cutoff}}=\frac{511}{2-\sqrt{1-\mu_{2}^{2}}}\text{ keV},$ (150)

above which no reflected photon can be seen. The results of $g_{1}$ for
various $\mu_{2}$ are plotted in Fig. (21), where the analytical results of
Hua & Lingenfelter (1992) are also shown for comparison. One can see that the
agreement between these results is quite well.

###  4.9. Comptonization of Polarized Radiation from a Reflecting Disk with a
Hot Corona

In the last two applications, we mainly focus on the unpolarized Compton
scattering dominated RT. In this section we proceed to discuss the
applications with polarization and calculate the analogous spectra emerging
from the hot corona of an optically thick disk. The effect that radiations can
be reflected from the disk when they irradiate the disk surface is also
considered. The reflected radiation $\mathbf{I}_{r}$ in the direction of
$\mathbf{\Omega}$ equals the sum of contributions coming from all other
incident directions $\mathbf{\Omega}^{\prime}$, i.e., (Poutanen et al.,
1996b),

$\displaystyle\begin{split}\mathbf{I}_{r}(\nu,\mathbf{\Omega})&=\int\mathbf{G}(\nu,\mathbf{\Omega};\nu^{\prime},\mathbf{\Omega}^{\prime})\mathbf{I}(\nu^{\prime},\mathbf{\Omega}^{\prime})d\nu^{\prime}d\mathbf{\Omega}^{\prime}.\end{split}$
(151)

where $\mathbf{G}(\mathbf{\Omega},\nu;\mathbf{\Omega}^{\prime},\nu^{\prime})$
is the Green function obtained by solving the RT process in the disk
(Magdziarz & Zdziarski, 1995). The matrix function $\mathbf{QS}/4\mu$ in Eq.
(143) is actually such a Green function. We will use it to get the reflected
radiations. In Appendix E, we discussed how to incorporate the reflection term
$\mathbf{I}_{r}$ into the integral equation and sample the modified transfer
kernel.

As aforementioned the transformation between the incident and Compton
scattered Stokes parameters are given by Eqs. (69) and (70). Considering the
distribution of electron gas and the unpolarized situation, this
transformation formula should be modified as

$\displaystyle\begin{split}&\displaystyle\mathbf{\Psi}^{\prime}=\frac{1}{4\pi}\int_{1}^{\infty}d\gamma
N_{e}(\gamma)\int_{-1}^{1}d\mu_{e}(1-\mu_{e}v)\int_{0}^{2\pi}d\varphi_{e}\\\
&\times\frac{r^{2}_{e}}{\gamma^{2}(1-\mu_{e}v)^{2}}\left(\frac{\nu^{\prime}}{\nu}\right)^{2}\mathbf{FM}(\varphi)\mathbf{\Psi}.\end{split}$
(152)

Similarly the integrals in the above equation can be evaluated by sampling a
set of $\gamma,\mu_{e},\varphi_{e}$. With them and after the multiplication of
matrices $\mathbf{FM}(\varphi)$ and $\mathbf{\Psi}$, we have

$\displaystyle\begin{split}&\displaystyle\Psi^{\prime(X)}=A_{0}f_{X}(\mu^{\prime},\varphi^{\prime};\mu,\varphi)\Psi^{(I)},\\\
&\displaystyle\Psi^{\prime(V)}=A_{0}f_{V}(\mu^{\prime},\varphi^{\prime};\mu,\varphi)\Psi^{(V)},\end{split}$
(153)

where
$A_{0}=(r^{2}_{e}/2)(\epsilon^{\prime}/\epsilon)^{2}/\sigma_{\text{KN}}(\epsilon)$
and $X=(I,Q,U)$. Transforming $f_{X}$ and $f_{V}$ into the rest frame of the
electron, we have

$\displaystyle\left\\{\begin{array}[]{ll}&\displaystyle\widetilde{f_{I}}=F_{0}+F_{3}(q\cos
2\Phi^{\prime}+u\sin 2\Phi^{\prime}),\\\
&\displaystyle\widetilde{f_{Q}}=F_{3}+F_{33}(q\cos 2\Phi^{\prime}+u\sin
2\Phi^{\prime}),\\\ &\displaystyle\widetilde{f_{U}}=F_{11}(-q\sin
2\Phi^{\prime}+u\cos 2\Phi^{\prime}),\\\
&\displaystyle\widetilde{f_{V}}=F_{22},\end{array}\right.$ (158)

where $q=\Psi^{(Q)}/\Psi^{(I)},u=\Psi^{(U)}/\Psi^{(I)}$, and
$(\Psi^{\prime},\Phi^{\prime})$ are the scattered polar and azimuth angles in
the rest frame of the electron. We also use $A_{0}\widetilde{f_{I}}$ to
construct a PDF for $(\Psi^{\prime},\Phi^{\prime})$ and sample it to get the
scattered direction. After transforming it back into the static frame, we get
$(\mu^{\prime},\varphi^{\prime})$ and the scattered $\mathbf{\Psi}^{\prime}$
as

$\displaystyle\left\\{\begin{array}[]{ll}&\displaystyle\Psi^{\prime(X)}=\frac{f_{X}(\mu^{\prime},\varphi^{\prime};\mu,\varphi)}{f_{I}(\mu^{\prime},\varphi^{\prime};\mu,\varphi)}\Psi^{(I)},\\\
&\displaystyle\Psi^{\prime(V)}=\frac{f_{V}(\mu^{\prime},\varphi^{\prime};\mu,\varphi)}{f_{I}(\mu^{\prime},\varphi^{\prime};\mu,\varphi)}\Psi^{(V)}.\end{array}\right.$
(161)

The estimation functions for the Stokes parameters are given by

$\displaystyle\left\\{\begin{array}[]{ll}&\displaystyle
F_{X}=A_{0}\widetilde{f_{X}}(\mathbf{\Omega}_{\text{obs}})\Psi^{(I)}_{m}\exp\left[-\tau(\nu,\mathbf{\Omega}_{\text{obs}})\right],\\\
&\displaystyle
F_{V}=A_{0}\widetilde{f_{V}}(\mathbf{\Omega}_{\text{obs}})\Psi^{(V)}_{m}\exp\left[-\tau(\nu,\mathbf{\Omega}_{\text{obs}})\right].\end{array}\right.$
(164)

With these preparations, we can now proceed to complete the calculations
mentioned at the beginning of this section. The source thermal radiation has a
temperature $kT_{bb}=10$ eV and is located isotropically and uniformly on the
disk surface. The radiations scattered back into the disk will re-emerge with
an altered SPs determined by Eq. (143). The results are plotted in Figs. 22
and 23. The corona temperature and Thomson optical depth are $T_{e}$ and
$\tau=h\sigma_{T}n_{e}$, where $h$ is the height of corona. The values of
these parameters correspond to those of Figs. 5 and 6 of Poutanen & Svensson
(1996a) (also refer to Figs 12 and 13 of Schnittman & Krolik (2013)). One can
see that our results are similar but different with them, especially for the
curves corresponding to each single scattering. In our scheme, since each
scattering (including the reflection) can make contributions to the spectra,
we do not reset the variable of scattering number $N_{\text{scat}}$ to zero
after the photons are reflected from the disk. It is contrast to the treatment
of Schnittman & Krolik (2013). The reflection processes are described by Eq.
(143), where the photoelectric absorption and energy losses due to Compton
recoil are not considered. Therefore the reflection processes will not change
the energies of photons. We believe that these different treatments give rise
to the discrepancies in the spectra profiles and the polarization degrees.

One can see that the high energy parts of the spectra, almost greater than
$\sim 100$ keV, are quite noise and oscillate heavily, because the photons
received by the observer at that energy band have experienced a large averaged
scattering number. The weights of these photons have been attenuated heavily,
hence the polarization signals comprised by these weights are difficult to be
resolved (Schnittman & Krolik, 2013).

Figure 22.— The spectra (top panels) and polarization degrees (bottom panels)
of radiations emerging from a plane-parallel corona above an accretion disk.
The solid blue curves correspond to the total spectra and the solid red,
dotted, dashed, dot-dashed, triple-dot-dashed curves correspond to scattering
number of 0, 1, 2, 3, 4, respectively. Compare to Figure 5 of Poutanen &
Svensson (1996a). Figure 23.— The same as in Fig. 22, but for $T_{e}=353$ keV
and $\tau=0.05$. Compare to Figure 6 of Poutanen & Svensson (1996a).

###  4.10. The Scattering of Polarized Photon Beam by Relativistic Electrons

In this section we will consider the polarization states of a polarized photon
beam scattered off only once by a relativistic electron gas that can have an
arbitrary distribution function $N_{e}(\gamma)$. We denote the energies,
polarization and wave vectors for the incident and scattered photons as
$\epsilon,\boldsymbol{\varepsilon},\mathbf{k}$ and
$\epsilon^{\prime},\boldsymbol{\varepsilon}^{\prime},\mathbf{k}^{\prime}$
respectively, and the scattering angle between $\mathbf{k}$ and
$\mathbf{k}^{\prime}$ is $\theta^{\prime}$. Then the scattered intensity can
be expressed semi-analytically as follows (Bonometto et al., 1970; Bonometto &
Saggion, 1973a, b)

$\displaystyle\begin{split}J_{\boldsymbol{\varepsilon}^{\prime}}&=\pi\left(\frac{e^{2}}{4\pi}\right)\frac{\epsilon^{\prime}}{\epsilon}E_{\text{min}}\left\\{\left|\boldsymbol{\varepsilon}\cdot\boldsymbol{\varepsilon}^{\prime*}+\frac{(\mathbf{k}\cdot\boldsymbol{\varepsilon}^{\prime*})(\mathbf{k}^{\prime}\cdot\boldsymbol{\varepsilon})}{1-\cos\theta^{\prime}}\right|^{2}\Sigma_{1}\right.\\\
&\left.+\left|\boldsymbol{\varepsilon}\cdot\boldsymbol{\varepsilon}^{\prime}+\frac{(\mathbf{k}\cdot\boldsymbol{\varepsilon}^{\prime})(\mathbf{k}^{\prime}\cdot\boldsymbol{\varepsilon})}{1-\cos\theta^{\prime}}\right|^{2}\Sigma_{2}+\Sigma_{2}\right\\},\end{split}$
(165)

where

$\displaystyle\begin{split}&\displaystyle\Sigma_{1}=\int_{0}^{1}m(\gamma)\left(x^{2}-\frac{1}{x^{2}}+2\right)dx,\\\
&\displaystyle\Sigma_{2}=\int_{0}^{1}m(\gamma)\frac{(1-x^{2})^{2}}{x^{2}}dx,\\\
&\displaystyle x=\frac{E_{\text{min}}}{\gamma},\quad\quad
E_{\text{min}}=\sqrt{\frac{1}{2}\frac{\epsilon^{\prime}}{\epsilon}\frac{1}{(1-\cos\theta^{\prime})}},\end{split}$
(166)

and the quantity $m(\gamma)$ is related to $N_{e}(\gamma)$, and

$\displaystyle\begin{split}m(\gamma)=\frac{N_{e}(\gamma)}{\gamma^{2}}.\end{split}$
(167)

For simplicity we denote $\displaystyle
J_{0}=\pi\left(\frac{e^{2}}{4\pi}\right)\frac{\epsilon^{\prime}}{\epsilon}E_{\text{min}}$.
Using the base vectors $\Ket{\mathbf{e}_{x}},\Ket{\mathbf{e}_{y}}$, the
polarization vectors $\boldsymbol{\varepsilon}$ and
$\boldsymbol{\varepsilon}^{\prime}$ can be decomposed as

$\displaystyle\begin{split}\boldsymbol{\varepsilon}=\alpha\Ket{\mathbf{e}_{x}}+\beta\Ket{\mathbf{e}_{y}},\quad\boldsymbol{\varepsilon}^{\prime}=\alpha^{\prime}\Ket{\mathbf{e}^{\prime}_{x}}+\beta^{\prime}\Ket{\mathbf{e}^{\prime}_{y}}.\end{split}$
(168)

Then with $\alpha,\beta$ and $\alpha^{\prime},\beta^{\prime}$, Eq. (165) can
be recast as (Bonometto et al., 1970)

$\displaystyle\begin{split}J_{\boldsymbol{\varepsilon}^{\prime}}&=J_{0}\left\\{c_{1}\Sigma_{1}+c_{2}\Sigma_{2}+\Sigma_{2}\right\\}.\end{split}$
(169)

where $c_{1}=\left|\alpha\alpha^{\prime*}-\beta\beta^{\prime*}\right|^{2}$ and
$c_{2}=\left|\alpha\alpha^{\prime}-\beta\beta^{\prime}\right|^{2}$. Without
loss of generality, one can choose $\alpha=x,\,\beta=ye^{i\theta}$ (obviously
$x^{2}+y^{2}=1$), the coefficients $c_{1}$ and $c_{2}$ can be written as

$\displaystyle\begin{split}c_{1}&=x^{2}x^{\prime
2}-2xx^{\prime}yy^{\prime}\cos(\theta-\theta^{\prime})+y^{2}y^{\prime 2},\\\
c_{2}&=x^{2}x^{\prime
2}-2xx^{\prime}yy^{\prime}\cos(\theta+\theta^{\prime})+y^{2}y^{\prime
2}.\end{split}$ (170)

Usually we describe the polarization state by SPs: $I,Q,U,V$ or the normalized
SPs: $\xi_{1}=U/I,\,\xi_{2}=V/I,\,\xi_{3}=Q/I$. For simplicity we define a
vector $\mathbf{S}=(1,\xi_{3},\xi_{1},\xi_{2})$. By using Eq. (56), we can
immediately obtain $\xi_{i}$ from $x,y$ and $\theta$:

$\displaystyle\begin{split}&\xi_{1}=2xy\cos\theta,\quad\xi_{2}=2xy\sin\theta,\quad\xi_{3}=x^{2}-y^{2},\end{split}$
(171)

and vice versa

$\displaystyle\begin{split}&x=\sqrt{\frac{1+\xi_{3}}{2}},\quad
y=\sqrt{\frac{1-\xi_{3}}{2}},\\\
&\cos\theta=\frac{\xi_{1}}{\sqrt{\xi_{1}^{2}+\xi_{2}^{2}}},\quad\sin\theta=\frac{\xi_{2}}{\sqrt{\xi_{1}^{2}+\xi_{2}^{2}}}.\end{split}$
(172)

For any two orthogonal polarization states of
$\boldsymbol{\varepsilon}=(\alpha,\beta)$ and
$\boldsymbol{\varepsilon}_{\perp}=(\alpha_{\perp},\beta_{\perp})$, since
$\boldsymbol{\varepsilon}\cdot\boldsymbol{\varepsilon}_{\perp}^{*}=\alpha\alpha_{\perp}^{*}+\beta\beta_{\perp}^{*}=0$,
without loss of generality we can obtain
$\boldsymbol{\varepsilon}_{\perp}=(-\beta^{*},\alpha^{*})$. Using Eq. (56),
one can immediately obtain that the normalized SPs for
$\boldsymbol{\varepsilon}_{\perp}$ is
$\mathbf{S}_{\boldsymbol{\varepsilon}_{\perp}}=(1,-\xi_{3},-\xi_{1},-\xi_{2})$
if the corresponding ones for $\boldsymbol{\varepsilon}$ is
$\mathbf{S}_{\boldsymbol{\varepsilon}}=(1,\xi_{3},\xi_{1},\xi_{2})$. The total
intensity for the scattered radiation is the sum of
$J_{\boldsymbol{\varepsilon}^{\prime}}$ over any two orthonormal polarization
states, and one can easily obtain (Bonometto et al., 1970)

$\displaystyle\begin{split}J_{I^{\prime}}=J_{\boldsymbol{\varepsilon}^{\prime}}+J_{\boldsymbol{\varepsilon}^{\prime}_{\perp}}=J_{0}\left\\{\Sigma_{1}+3\Sigma_{2}\right\\}.\end{split}$
(173)

Due to the additive property of SPs (Chandrasekhar, 1960), the normalized SPs
for $J_{I^{\prime}}$ is
$\mathbf{S}_{I^{\prime}}=\mathbf{S}_{\boldsymbol{\varepsilon}^{\prime}}+\mathbf{S}_{\boldsymbol{\varepsilon}^{\prime}_{\perp}}=(1,0,0,0)$.

Similarly, the scattered intensity for the polarization state of $Q^{\prime}$
is given by

$\displaystyle\begin{split}J_{Q^{\prime}}&=J_{S^{\prime}_{+}}-J_{S^{\prime}_{-}}\\\
&=J_{0}\left[(c_{1}^{+}-c_{1}^{-})\Sigma_{1}+(c_{2}^{+}-c_{2}^{-})\Sigma_{2}\right],\end{split}$
(174)

where $\mathbf{S}^{\prime}_{\pm}=(1,\pm 1,0,0)$. The corresponding
coefficients $c_{1}^{\pm}$ and $c_{2}^{\pm}$ are obtained from Eq. (170).
Similarly, the scattered intensities for polarization states $U^{\prime}$ and
$V^{\prime}$ are given by

$\displaystyle\begin{split}J_{U^{\prime}}&=J_{S^{\prime}_{+}}-J_{S^{\prime}_{-}}\\\
&=J_{0}\left[(c_{1}^{+}-c_{1}^{-})\Sigma_{1}+(c_{2}^{+}-c_{2}^{-})\Sigma_{2}\right],\\\
J_{V^{\prime}}&=J_{S^{\prime}_{+}}-J_{S^{\prime}_{-}}\\\
&=J_{0}\left[(c_{1}^{+}-c_{1}^{-})\Sigma_{1}+(c_{2}^{+}-c_{2}^{-})\Sigma_{2}\right],\end{split}$
(175)

where the coefficients $c_{1}^{\pm}$ and $c_{2}^{\pm}$ are obtained from the
normalized SPs $\mathbf{S}^{\prime}_{\pm}=(1,0,\pm 1,0)$ for $U^{\prime}$ and
$\mathbf{S}^{\prime}_{\pm}=(1,0,0,\pm 1)$ for $V^{\prime}$, respectively.
Since any polarization state is completely determined by SPs, therefore once
$J_{I^{\prime}},J_{Q^{\prime}},J_{U^{\prime}},J_{V^{\prime}}$ are obtained,
the scattered polarization state is also determined.

With these preparations, we can readily obtain the scattered polarization
states semi-analytically, which can be used to verify our MC scheme. We first
reproduce the results in Fig. 2 of Mościbrodzka (2020), where four polarized
incident beams are scattered off by a group of hot electron gas distributed
according to the relativistic Maxwell function (see Eq. (27)). The
dimensionless temperature of the gas is $\Theta_{e}=kT_{e}/m_{e}c^{2}=100$.
The incident beams with fixed energy $\epsilon=2.5\times 10^{-11}$ propagate
along the positive $z$-axis and the viewing direction is
$(\theta^{\prime},\varphi^{\prime})=(85^{\circ},0)$. The normalized SPs of the
four incident beams are $\mathbf{S}_{\text{in},1}=(1,1,0,0)$,
$\mathbf{S}_{\text{in},2}=(1,-1,0,0)$, $\mathbf{S}_{\text{in},3}=(1,0,1,0)$,
and $\mathbf{S}_{\text{in},4}=(1,0,0,1)$, respectively. The scattered
intensities for all beams are the same and independent on the polarizations of
the incident beams (see Eq. (173)). The intensity (left panel) and fractional
polarizations (right panel) of the scattered photons distributed with respect
to the photon energy are shown in Fig. 24, where only the nonvanishing
normalized SPs are plotted. One can see that the numerical results produced by
our MC scheme agree with the semi-analytical results of Bonometto et al.
(1970) very well, especially in the low energy band. The noise of the high
energy part is relatively large due to the sparse productions of energetic
photons comparing to that of the low energy ones.

Next we consider the case where the electron gas is distributed in a power law
profile, i.e.,

$\displaystyle N_{e}(\gamma)=K\gamma^{-\alpha},\quad\gamma\geq\gamma_{1},$
(176)

where $K=(\alpha-1)/\gamma_{1}^{1-\alpha}$ is the normalization factor and
$\alpha$ is the power law index. Then $m(\gamma)=K\gamma^{-(\alpha+2)}$, one
can immediately obtain the expressions for quantities $\Sigma_{1}$ and
$\Sigma_{2}$ as (Bonometto et al. (1970))

$\displaystyle\begin{split}\Sigma_{1}=\frac{K}{E^{\alpha+1}_{\text{min}}}x_{1}^{\alpha+1}\left[\frac{x_{1}^{4}}{\alpha+5}-\frac{1}{\alpha+1}+\frac{2x_{1}^{2}}{\alpha+3}\right],\\\
\Sigma_{2}=\frac{K}{E^{\alpha+1}_{\text{min}}}x_{1}^{\alpha+1}\left[\frac{x_{1}^{4}}{\alpha+5}+\frac{1}{\alpha+1}-\frac{2x_{1}^{2}}{\alpha+3}\right],\end{split}$
(177)

where (notice that $E_{\text{min}}$ is dependent on the energies of scattered
photons)

$\displaystyle x_{1}=\left\\{\begin{array}[]{lr}\displaystyle
1,&\quad\quad\text{ if }E_{\text{min}}\geq\gamma_{1},\\\
\displaystyle\frac{E_{\text{min}}}{\gamma_{1}},&\quad\quad\text{ if
}E_{\text{min}}<\gamma_{1}.\end{array}\right.$ (180)

Our results in comparison with the analytic ones of Bonometto et al. (1970)
are shown in Figs. 25 and 26. In these calculations we set $\gamma_{1}=60.0$,
the power law index $\alpha=3$ and the incident monoenergy $\epsilon=2.5\times
10^{-11}$. Using the inverse CDF method, we can easily sample $\gamma$ from
$N_{e}(\gamma)$, i.e., $\gamma_{\xi}=\gamma_{1}(1-\xi)^{1/(1-\alpha)}$. From
Eq. (177) one can see that there exist a turning point
$\epsilon^{\prime}_{\text{tr}}$ for the scattered photon energy
$\epsilon^{\prime}$, above which all of the fractional polarizations:
$\xi_{3},\xi_{2},\xi_{1}$ are constants (but notice that the SPs $I,Q,U,V$ are
not). From the condition $E_{\text{min}}=\gamma_{1}$ one can easily obtain
that
$\epsilon^{\prime}_{\text{tr}}/\epsilon=2(1-\cos\theta^{\prime})\gamma_{1}^{2}$.
In Figs. 25 and 26, where $\theta^{\prime}=85^{\circ}$ and $\gamma_{1}=60.0$,
we have $\log_{10}(\epsilon^{\prime}_{\text{tr}}/\epsilon)\approx 3.818$. The
normalized SPs of the four incident beams in Fig. 25 are the same with those
in Fig. 24. While in Fig. 26, the incident SPs are
$\mathbf{S}_{\text{in}}=(1,0.25,0.433,0.866)$ and all of the scattered SPs are
nonvanishing.

Figure 24.— The intensity (left) and normalized SPs (right) of photons
scattered off a group of thermal relativistic electron gas with various
incident polarizations $\mathbf{S}_{\text{in}}$. The dimensionless temperature
of the electron gas is $\Theta_{e}=100$. The incident energy is monoenergetic
and $\epsilon=h\omega/m_{e}c^{2}=2.5\times 10^{-11}$. The solid lines and
histograms represent the results obtained by the semi-analytical formulae of
Bonometto et al. (1970) and our MC scheme, respectively. Compare to Fig. 2 of
Mościbrodzka (2020).

Figure 25.— Same as in Fig. 24, but the electron gas is distributed in a power
law profile given be Eq. (176), where $\gamma_{1}=60.0$, $\alpha=3$ and
$\epsilon=2.5\times 10^{-11}$.

Figure 26.— Same as in Fig. 25, but for $S_{\text{in}}=(1,0.25,0.433,0.866)$.

###  4.11. Polarized RT with Faraday Rotation and Conversion

In section 3.4 we discussed a scheme to solve the polarized RT with Faraday
rotation and conversion incorporated. Now we will test its validation by
applying it to solve two RT problems that have analytic solutions presented in
Appendix D. Similar tests have been carried out by Dexter (2016) and
Mościbrodzka & Gammie (2018) using alternative numerical methods. We emphasize
that if there has no scattering in these RT problems, the MC method just
provides another choice to solve them and has no any advantages comparing to
these numerical methods. The advantage of our method lies in solving the RT
with scattering. The results of the two tests are demonstrated in Fig. 27,
where the solid and dotted lines represent the analytical and our numerical
results, respectively. One can see that they are consistent with each other
very well.

Figure 27.— The tests of MC scheme for RT with Faraday rotation and conversion
incorporated against the analytical solutions. Left panel: only the SPs $I$
and $Q$ are involved, and the parameters are:
$\alpha_{I}=2.0,\alpha_{Q}=1.5,j_{I}=2.1,j_{Q}=1.2$ and $\alpha=3.0$. Right
panel: the Figure illustrates the evolution of SPs $Q$, $U$ and $V$. The
parameters are:
$\alpha_{IQUV}=0.0,\rho_{Q}=7.5,\rho_{U}=3.4,\rho_{V}=7.2,j_{Q}=-10.717,j_{U}=9.033,\rho_{V}=9.0587$
and $\alpha=200.0$. In the two panels, the solid and dotted lines represent
the analytical and MC numerical results, respectively.

##  5\. DISCUSSION

In this paper we have proposed a new MC scheme based on Neumann series
solutions of Fredholm integral equations of second kind to solve the RT
problems with scattering incorporated. From this perspective, the problem of
RT is equivalent to evaluate infinite terms of multiple integrals
simultaneously and the MC method is naturally introduced since it has a
remarkable efficiency in evaluating multiple integrals. Due to the special
structure of Neumann solution, the procedure of integral evaluation amounts to
generate random sequences in the phase space. The combination of these
sequences with their corresponding weight functions can immediately give the
estimations for the observational quantities.

We have reviewed and redescribed the procedure to solve an RTE by the MC
method completely and systematically from the perspective of integral equation
and its Neumann solution. We particularly emphasize that all complicated and
abstract MC sampling algorithms for the transport and scattering of a photon
can be understood and explained mathematically. Meanwhile we try to ignore the
corresponding physical explanations. Such a treatment not only makes the MC
method to be more understandable but also can avert mistakes caused by
physical intuitions. Most importantly it enables one to have the flexibility
of choosing various sampling PDFs and corresponding weight functions
conveniently. These choices are equivalent but with different computational
efficiency and accuracy. Particularly the widely adopted conventional photon
tracing scheme can be regarded as a special choice, where the PDF for position
transport has definitions in the whole space.

We emphasize the importance of recording function which connects the Neumann
solution with the observational quantities directly. Especially, any delta
function (which is exactly the reason why the photon tracing scheme has a low
computational efficiency, since tremendous calculations are inevitably
abandoned) contained in the recording function can be eliminated in advance.
After that we can choose anther recording function, which enables that any
sample in a random sequence can make contribution to the quantity that is
under calculation. This estimation strategy can significantly improve the
calculation accuracy, especially for systems with axial or spherical symmetry.
Usually under the same condition our results have a relatively higher
precision comparing to that of the photon tracing scheme. In some ideal
situations, the precision is even comparable with the semi-analytical results.
But for each estimation, our scheme requires to calculate the recording
function and total optical depth from the scattering site to the boundary of
the region. It will inevitably occupy additional computational resources.

This new scheme can be directly extended to deal with polarized radiative
transfer processes, especially for that the Faraday rotation and conversion
effects are incorporated. The polarizations are described by the SPs and the
RTEs become a set of differentio-integral equations. The key point is that one
should choose common transport PDFs (including position transport and
scattering) that is shared by all components of the SPs, meanwhile their
updating ways from the incident direction to the scattered direction and
estimation functions are different. One of the advantages of our scheme in
dealing with the polarized Compton scattering is that we do not have to carry
out the Lorentz transformations in terms of SPs between different references.
This procedure that involves complicated SPs transformations and rotations
(Krawczynski, 2012) now can be avoided.

According to the scheme, a public available MCRT code, named Lemon, has been
developed with Fortran. In the development, we find that the Object-Oriented
programming is very suitable and convenient to implement our scheme. We first
define a basic particle class, which contains the primary members and methods,
such as coordinates, momentum, energy, weights, etc.. Then we extend this
class gradually by class inheritance and adding new members and methods to the
subclass that can implement various special functions. Each extension
corresponds one function and the relevant variables and subroutines of which
are saved in a single module. These functions include initiating all
parameters, starting a sequence, position transport, photon scattering and
quantity estimation and so on. As one of the functions is modified, any other
will not be affected. Finally we can instantiate these subclasses and obtain
different objects. Calling the object’s methods we can generate a scattering
sequences and estimate the observational quantity as the generation goes on.

In order to verify the correctness of the scheme and the code, we have
reproduced a lot of classic results and made comparisons with them, which have
demonstrated the validation and advantages of our scheme. Since the algorithms
for sampling various PDFs are crucial for MC method, they have been discussed
in detail and the proofs for some of them are also provided in the appendix.

MCRT has an intrinsic parallelization property due to the calculations of
random sequences are independent from each other. In order to speed up the
calculations, the MPI parallel scheme has been incorporated into the code. In
the implementation of parallelization, all of the MPI processes are treated
equally, i.e., they are initiated with the same parameters, load the same data
files (that contain the tabulated total scattering cross section, the
temperature and the electron number density distributions, etc.) and share the
same emissivity, scattering and absorption coefficients. Which means that each
MPI process requires the access to the full radiating domain. And this
requirement can be easily fulfilled for all of our test examples, since they
are relatively simple (or the radiating domain can be described semi-
analytically) and do not involve any applications to complicated numerical
simulations. The sizes of these data files are small and can be easily loaded
by each process.

Without any applications to numerical simulations is really the shortage of
our code testing. Considering the extreme significance and widely applications
of magnetohydrodynamic simulations in nowadays, we really hope that we can
apply our code to these simulations in the near future. Numerical simulations
with large simulation domains can no longer be load by each process. To
overcome these difficulties, new parallelization algorithms and techniques are
in desperate needed (e.g., down sampling of the simulation grid and domain
decomposition). And the structure of MPI parallelization implemented in the
present code should be redesigned completely.

Anther limit of our code testing is that all of the test problems are confined
to static configurations. However the fluid flows in numerical simulations
usually have large nonuniform velocities that even can be comparable to the
speed of light. So can our code handle the RT processes in fluid flow with
irregular velocity field consistently? The answer is yes or in principle it
can. The procedure can be stated as follows.

First we must choose a global and static reference, since the Neumann solution
consists of integrals defined over the entire radiative domain. Such reference
obviously exists, in which the velocity field of the fluid flow is also
appropriately defined on a grid. Next our scheme requires the mathematical
expressions for all relevant quantities of the fluid flow (such as the
emissivity, absorption and scattering coefficients, the mass and particle
number density, temperature distribution) to be explicitly given in this
frame. These quantities are well defined in the comoving (or rest) frame of
the fluid flow. So simply through a Lorentz transformation, we can obtain
these expressions in the global static frame (one should also notice that in
order to get the values of these quantities at an arbitrary position from the
simulation data that defined on a grid, a suitable interpolation scheme should
be introduced). The Lorentz transformations for these quantities have been
extensively studied and can be found in many textbooks (e.g., in the appendix
of Pomraning (1973)). With these quantities are specified, we can write down
the RTE immediately and obtain the Neumann solution as well. Then the
following steps are just routines that have been discussed in the previous
sections.

Another problem needs be carefully treated when dealing with RT in optically
thick media with our scheme. In this special case, the aforementioned
criterion (as the weight of these photons falling below a threshold) to
truncate a scattering sequence will become extremely inefficient and therefore
should be rejected. The reason can be stated as follow. Since the optically
thick region tends to trap any photon that travels into it, the samples of a
scattering sequence could be disproportionately distributed over this region.
Due to the factor $\exp(-\tau)$ in the recording function, the large optical
depth will heavily attenuate the contributions made by these samples to the
observed flux. Thus the whole contribution of the sequence is quite low.
Meanwhile the evaluation of $\tau$ for each sample is time consuming in
optical thick region, especially for post-processing of finite volume
simulations with large grids. The situation can be even worse if the weight
decreases slowly, such as in cases where the true absorption can be neglected
or no absorption exists at all. If we still require the sequence to be
truncated as the weight is smaller than the threshold, the scattering times
will be a quite large number, which makes the computational cost for
evaluating $\tau$ becomes unacceptable. The cost almost offsets the limited
gain in the lower variance.

We think the best way to alleviate this difficulty is to truncate the
scattering sequence as soon as possible. Once the photon is scattered into
high optical depth region, the frequent scattering will trap it and prevent it
to escape from the region. Considering the high computational costs and low
contributions made by the subsequent samples, we can safely truncate the
sequence once the optical depth of a sample is larger than a threshold.

Another useful truncation standard is related to the scattering numbers. For
systems with low scattering optical depth, we approximately have
$\exp(-\tau_{s})\approx 1+\tau_{s}$ and $A=1-\exp(-\tau_{s})\approx\tau_{s}$,
then $I_{m}\propto A_{1}A_{2}\cdots A_{m}\exp(-\tau_{s})\approx\tau_{s}^{m}$,
which will decrease rapidly as $m$ increases. So we can truncate the sequence
as $m$ is greater than a particular number. For example, from Fig. 18, one can
see that the contributions made by photons that have been scattered for more
than 4 times can be safely rejected.

In present paper, the discussions are mainly restricted to the classic RT
problems in flat spacetime. Our scheme is dependent heavily on the specific
forms of RTEs, which should be given explicitly to obtain the Neumann solution
and the transport kernel. In curved spacetime, the RTEs are usually given in a
covariant form (Younsi et al., 2012; Gammie & Leung, 2012; Bronzwaer et al.,
2020), but the scattering incorporated RTEs are rarely provided completely and
explicitly, especially for the polarized situations. One exception is in
cosmology regime, where in order to study the polarized radiative transfer and
perturbations of CMB in the cosmic metrics, the RTEs with polarized Compton
scattering have been extensively discussed (Kosowsky, 1996; Portsmouth &
Bertschinger, 2004; Weinberg, 2008). For the MCRTs with scattering surround
Black holes, instead of writing down the RTEs and starting from the Neumann
solutions, one prefers to generate photons and trace them until to the endings
(Dolence et al., 2009; Schnittman & Krolik, 2013; Zhang et al., 2019), and the
polarizations are usually described by the so-called polarization vectors,
which are parallel transported along geodesics. In the future, we hope to
extend our scheme from flat to curved spacetime and incorporate it with the
results of GRMHD simulations to mimic more real RT processes, which are
crucial to interpret the polarized observations generated by accreting black
holes.

Finally we want to point out a difference between flat and curved spacetime
which may cause difficulties for our scheme. In flat spacetime, once the line
of sight of the observer, $\mathbf{\Omega}_{\text{obs}}$, is provided, the
scattered direction $\mathbf{\Omega}^{\prime}$ in the scattering kernel
$C(\nu,\mathbf{\Omega}\rightarrow\nu^{\prime},\mathbf{\Omega}^{\prime})$ can
be replaced by $\mathbf{\Omega}_{\text{obs}}$ directly and
$C(\nu,\mathbf{\Omega}\rightarrow\nu^{\prime},\mathbf{\Omega}_{\text{obs}})$
will be taken as a weight for estimations. While in curved spacetime, due to
the light bending, we can not carry out such a replacement anymore, unless we
determine the geodesic connecting the observer and the scattering point and
then $\mathbf{\Omega}^{\prime}$ should be replaced by the unit tangent vector
of the geodesic at the scattering point rather than
$\mathbf{\Omega}_{\text{obs}}$. It implies that for each estimation of our
scheme in curved spacetime, a geodesic with two fixed points (i.e., observer
at infinity and the scattering point) needs be determined additionally.
Fortunately, some excellent codes already exist for computing the photon
geodesics in curved spacetime. Especially the public available code geokerr
(Dexter & Agol, 2009) can compute the geodesics between two points in Kerr
metric, even without full integrations.

## Acknowledgments

We thank a lot for the anonymous referee’s valuable comments and suggestions
on our manuscript, which have improved the quality of the manuscript and the
results greatly. YXL thanks Yuan Zun-li, Zhang Guo-bao, Hou Xian and Zhang
Wen-da for helpful discussions. Especially, we thank Zhang Guo-bao and Hou
Xian for their generosities of allowing us to use their personal workstations
freely, without which the calculations in this manuscript are impossible to be
finished. We also thank the Yunnan Observatories Supercomputing Platform, on
which our code is partly tested. We acknowledge the financial supports from
the National Natural Science Foundation of China (grant Nos. 11573060,
11661161010, U1838116, 11673060, U2031111, U1931204, 12073069) and Yunnan
Natural Science Foundation (Nos. 2019FB008).

##  Appendix A the proof of the algorithm for position transport

Here we present a proof for Algorithm 3 that randomly generates a scattering
distance $s$. What we need to prove is that the probability $p(s)$ for a
sample $s$ is selected by this algorithm is proportional to
$\sigma_{s}(s)\exp\left(-\int_{0}^{s}\sigma_{s}(l)dl\right)$, where
$\sigma_{s}(s)=\sigma_{s}(\mathbf{r}+s\mathbf{\Omega},\nu,\mathbf{\Omega})$ is
the scattering coefficient. According to the algorithm, we know that an $s$
can be expressed as the sum of $s_{i}$: $s=s_{1}+s_{2}+\cdots+s_{n}$, where
$n=1,2,\cdots$, and $s_{i}$ are samples of a random variable $l$ with the PDF
given by
$\displaystyle\widetilde{p}_{1}(l)=\sigma_{\text{max}}\exp(-\sigma_{\text{max}}l)/N,\,0\leq
l\leq s_{\text{max}}$, and $N=1-\exp(-\sigma_{\text{max}}s_{\text{max}})$ is
the normalization factor. Suppose $p_{n}(s)$ is the probability that the
$s_{i}$ ($i=1,\cdots,n$) are generated and their sum is accepted as a sample
of $s$. Then the total probability corresponding an $s$ is selected can be
expressed as

$\displaystyle p(s)=\sum^{\infty}_{n=1}p_{n}(s).$ (A1)

Obviously for $n=1$, we have

$\displaystyle\displaystyle
p_{1}(s)=\widetilde{p}_{1}(s)\frac{\sigma_{s}(s)}{\sigma_{\text{max}}}=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s).$
(A2)

As $n=2$, we have $s=s_{1}+s_{2}$, and the probability of generating
$s_{1},s_{2}$ by sampling $\widetilde{p}_{1}(l)$ and accepting them can be
written as

$\displaystyle\begin{split}&p(s_{1})\displaystyle=\widetilde{p}_{1}(s_{1})\left(1-\frac{\sigma_{s}(s_{1})}{\sigma_{\text{max}}}\right)\widetilde{p}_{1}(s_{2})\frac{\sigma_{s}(s_{1}+s_{2})}{\sigma_{\text{max}}}N\displaystyle=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s)[\sigma_{\text{max}}-\sigma_{s}(s_{1})].\end{split}$
(A3)

Notice that $s_{1}$ can actually take any value on interval $(0,s)$ and each
value corresponds a probability $p(s_{1})$, thus the total probability should
be

$\displaystyle\displaystyle
p_{2}(s)=\int_{0}^{s}p(s_{1})ds_{1}=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s)\int_{0}^{s}[\sigma_{\text{max}}-\sigma_{s}(s_{1})]ds_{1}.$
(A4)

Similarly, for $n=3$, we have

$\displaystyle\begin{split}p(s_{1},s_{2})&\displaystyle=\widetilde{p}_{1}(s_{1})\left(1-\frac{\sigma_{s}(s_{1})}{\sigma_{\text{max}}}\right)\widetilde{p}_{1}(s_{2})\left(1-\frac{\sigma_{s}(s_{1}+s_{2})}{\sigma_{\text{max}}}\right)N\widetilde{p}_{1}(s_{3})\frac{\sigma_{s}(s)}{\sigma_{\text{max}}}N\\\
&\displaystyle=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s)[\sigma_{\text{max}}-\sigma_{s}(s_{1})][\sigma_{\text{max}}-\sigma_{s}(s_{1}+s_{2})].\end{split}$
(A5)

Notice that $s_{1}$ can also take any value on interval $(0,s)$ and $s_{2}$ on
$(0,s-s_{1})$, thus we have

$\displaystyle\begin{split}&p_{3}(s)\displaystyle=\int_{0}^{s}ds_{1}\int_{0}^{s-s_{1}}ds_{2}p(s_{1},s_{2})=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s)\int_{0}^{s}ds_{1}\int_{0}^{s-s_{1}}ds_{2}[\sigma_{\text{max}}-\sigma_{s}(s_{1})][\sigma_{\text{max}}-\sigma_{s}(s_{1}+s_{2})].\end{split}$
(A6)

Then the expression for $p_{n}(s)$ can be obtained directly as

$\displaystyle\begin{split}&p_{n}(s)\displaystyle=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s)\int_{0}^{s}ds_{1}\int_{0}^{s-s_{1}}ds_{2}\cdots\int_{0}^{s-(s_{1}+\cdots+s_{n-2})}ds_{n-1}f(s_{1})f(s_{1}+s_{2})\cdots
f(s_{1}+\cdots+s_{n-1}),\end{split}$ (A7)

where $f(s)=\sigma_{\text{max}}-\sigma(s)$. We introduce a set of new
variables $y_{1},\cdots,y_{n-1}$ defined by

$\displaystyle
y_{1}=s_{1},\,\,\,\,y_{2}=s_{1}+s_{2},\,\,\,\,\cdots,\,\,\,\,y_{n-1}=s_{1}+\cdots+s_{n-1},$
(A8)

we have
$s_{1}=y_{1},\,\,s_{2}=y_{2}-y_{1},\,\,\cdots,\,\,s_{n-1}=y_{n-1}-y_{n-2}$.
One can easily show that the Jacobian of this transformation is
$J=\left|\partial(s_{1},\cdots,s_{n-1})/\partial(y_{1},\cdots,y_{n-1})\right|=1$.
So the volume element is unchanged, i.e., $dy_{1}dy_{2}\cdots
dy_{n-1}=ds_{1}ds_{2}\cdots ds_{n-1}$. Then $p_{n}(s)$ can be rewritten as

$\displaystyle
p_{n}(s)\displaystyle=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s)\int_{0}^{s}dy_{1}\int_{y_{1}}^{s}dy_{2}\cdots\int_{y_{n-2}}^{s}dy_{n-1}f(y_{1})f(y_{2})\cdots
f(y_{n-1}).$ (A9)

Notice that

$\displaystyle\int_{0}^{s}dy_{1}\int_{y_{1}}^{s}dy_{2}\cdots\int_{y_{n-2}}^{s}dy_{n-1}f(y_{1})\cdots
f(y_{n-1})=\frac{1}{(n-1)!}\left[\int_{0}^{s}f(y)dy\right]^{n-1},$ (A10)

we have

$\displaystyle
p_{n}(s)\displaystyle=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s)\frac{1}{(n-1)!}\left[\int_{0}^{s}[\sigma_{\text{max}}-\sigma_{s}(y)]dy\right]^{n-1}.$
(A11)

Eventually we arrive at

$\displaystyle\begin{split}p(s)&=\sum_{n=1}^{\infty}p_{n}(s)=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s)\sum_{n=1}^{\infty}\frac{1}{(n-1)!}\left[\int_{0}^{s}[\sigma_{\text{max}}-\sigma_{s}(y)]dy\right]^{n-1}\\\
&=\frac{\sigma_{s}(s)}{N}\exp(-\sigma_{\text{max}}s)\exp\left[\sigma_{\text{max}}s-\int_{0}^{s}\sigma_{s}(y)dy\right]=\frac{\sigma_{s}(s)}{N}\exp\left[-\int_{0}^{s}\sigma_{s}(y)dy\right],\end{split}$
(A12)

i.e., $p(s)\propto\sigma_{s}(s)\exp\left(-\int_{0}^{s}\sigma_{s}(l)dl\right)$.

##  Appendix B The Sampling Algorithm For PDFs with Triangular functions

When dealing with scattering processes, one often encounters PDFs that consist
of triangular functions, for example
$p(\varphi)=a_{0}+a_{1}\cos\varphi+a_{2}\cos
2\varphi+b_{1}\sin\varphi+b_{2}\sin 2\varphi$. If one uses the inverse CDF
method to sample such an PDF, a transcendental equation which is analogous to
Kepler equation but more complicated needs to be solved. However, we can
construct a new algorithm that utilizes the properties of piecewise central
symmetry and the odd and even of triangular functions to sample a more general
triangular function PDF $p(\varphi)$ given by

$\displaystyle p(\varphi)=a_{0}+\sum_{m=1}^{M}a_{m}\cos
m\varphi+\sum_{n=1}^{N}b_{n}\sin n\varphi,$ (B1)

where $M$ and $N$ can be any positive integers, the coefficients $a_{i}$ and
$b_{i}$ are given in a way that $p(\varphi)$ is always positive. We introduce
some new constants defined as

$\displaystyle p_{m}=\frac{|a_{m}|}{a_{0}},\quad p_{0}=\displaystyle
1-\sum_{m=1}^{M}p_{m},\quad\displaystyle P_{m}=\sum_{i=0}^{m}p_{i},$ (B2)

which satisfy the conditions: $0\leq p_{i}<1$ and
$p_{0}+p_{1}+\cdots+p_{M}=1$. Then $p(\varphi)$ can be recast as

$\displaystyle p(\varphi)=F(\varphi)+H(\varphi),$ (B3)

where

$\displaystyle\displaystyle
F(\varphi)=a_{0}\left[p_{0}+\sum_{m=1}^{M}p_{m}(1+s_{m}\cos
m\varphi)\right],\quad\displaystyle H(\varphi)=\sum_{n=1}^{N}b_{n}\sin
n\varphi,$ (B4)

and $s_{m}=\text{sign}(a_{m})$. The sampling procedure of $p(\varphi)$ is
provided in Algorithm. 4, where $\xi_{i}$ are random numbers,
$n=[\varphi_{t}/(\pi/m)]$ and $[x]$ represents the function that gets the
maximum integer that is smaller than $x$. A test example of this algorithm is
given in Fig. 28.

① Let $\varphi_{t}=2\pi\cdot\xi_{1}$

if  _$\left(\displaystyle 0\leq\xi_{2} <P_{0}\right)$_ then

$\varphi_{\xi}=\varphi_{t}$

else if  _$\left(\displaystyle P_{m-1}\leq\xi_{2} <P_{m}\right)$_ then

if  _$2\xi_{3}\leq 1+s_{m}\cos m\varphi$_ then

$\varphi_{\xi}=\varphi_{t}$

else

$\varphi_{\xi}=(2n+1)\pi/m-\varphi_{t}$

if
_$\left(\displaystyle\xi_{4}-1\leq\frac{H(\varphi_{\xi})}{F(\varphi_{\xi})}\right)$_
then

$\quad\varphi=\varphi_{\xi}$

else

$\varphi=2\pi-\varphi_{\xi}$

Algorithm 4 Triangular PDF Sampling Algorithm Figure 28.— A testing example of
Algorithm. 4 for sampling a PDF given by Eq. (B1). The coefficients $a_{i}$
and $b_{i}$ are all random numbers and $M=40,\,N=30$. The function
$p(\varphi)$ is plotted in red solid line, while the normalized sample numbers
generated by the algorithm is plotted in black dotted line.

##  Appendix C The Sampling Procedures For Rayleigh’s Scattering matrix in a
static frame

In section 4.5, we discussed the angular dependent SPs of radiations that are
diffusely reflected from a plane-parallel semi-infinite atmosphere. The
equations describing the RT in such a system are given by Eqs. (129), in which
the Rayleigh’s scattering matrix $\mathbf{P}$ is given by (Eqs. (220)-(224) of
(Chandrasekhar, 1960) Chap. I, page 42)

$\displaystyle\displaystyle\mathbf{P}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})=\mathbf{Q}[\mathbf{P}^{(0)}(\mu,\mu^{\prime})+\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime
2}}\mathbf{P}^{(1)}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})+\mathbf{P}^{(2)}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})],$
(C1)

where

$\displaystyle\displaystyle\mathbf{Q}=\left(\begin{array}[]{cccc}\displaystyle
1&\quad 0&\quad 0&\quad 0\\\ \displaystyle 0&\quad 1&\quad 0&\quad 0\\\
\displaystyle 0&\quad 0&\quad 2&\quad 0\\\ \displaystyle 0&\quad 0&\quad
0&\quad 2\\\
\end{array}\right),\quad\quad\displaystyle\mathbf{P}^{(0)}(\mu,\mu^{\prime})=\frac{3}{4}\left(\begin{array}[]{cccc}\displaystyle
2(1-\mu^{2})(1-\mu^{\prime 2})+\mu^{2}\mu^{\prime 2}&\quad\mu^{2}&\quad
0&\quad 0\\\ \displaystyle\mu^{\prime 2}&\quad 1&\quad 0&\quad 0\\\
\displaystyle 0&\quad 0&\quad 0&\quad 0\\\ \displaystyle 0&\quad 0&\quad
0&\quad\mu\mu^{\prime}\\\ \end{array}\right),$ (C10)

$\displaystyle\displaystyle\mathbf{P}^{(1)}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})=\frac{3}{4}\left(\begin{array}[]{cccc}\displaystyle
4\mu\mu^{\prime}\cos(\varphi^{\prime}-\varphi)&\quad 0&\quad
2\mu\sin(\varphi^{\prime}-\varphi)&\quad 0\\\ \displaystyle 0&\quad 0&\quad
0&\quad 0\\\ \displaystyle-2\mu^{\prime}\sin(\varphi^{\prime}-\varphi)&\quad
0&\quad\cos(\varphi^{\prime}-\varphi)&\quad 0\\\ \displaystyle 0&\quad 0&\quad
0&\quad\cos(\varphi^{\prime}-\varphi)\\\ \end{array}\right),$ (C15)

$\displaystyle\displaystyle\mathbf{P}^{(2)}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})=\frac{3}{4}\left(\begin{array}[]{cccc}\displaystyle\mu^{2}\mu^{\prime
2}\cos 2(\varphi^{\prime}-\varphi)&\quad-\mu^{2}\cos
2(\varphi^{\prime}-\varphi)&\quad\mu^{2}\mu^{\prime}\sin
2(\varphi^{\prime}-\varphi)&\quad 0\\\ \displaystyle-\mu^{\prime 2}\cos
2(\varphi^{\prime}-\varphi)&\quad\cos
2(\varphi^{\prime}-\varphi)&\quad-\mu^{\prime}\sin
2(\varphi^{\prime}-\varphi)&\quad 0\\\ \displaystyle-\mu\mu^{\prime 2}\sin
2(\varphi^{\prime}-\varphi)&\quad\mu\sin
2(\varphi^{\prime}-\varphi)&\quad\mu\mu^{\prime}\cos
2(\varphi^{\prime}-\varphi)&\quad 0\\\ \displaystyle 0&\quad 0&\quad 0&\quad
0\\\ \end{array}\right).$ (C20)

Notice that the above matrices are defined with respect to SPs
$I_{l},I_{r},U,V$ rather than $I,Q,U,V$. In MCRT, we prefer to use $I,Q,U,V$
to describe the polarizations instead, so we need to make a transformation to
get a new matrix $\widetilde{\mathbf{P}}$ that is defined in terms of
$I,Q,U,V$. One can show that the elements of the new matrix are given by

$\displaystyle\left\\{\begin{array}[]{llll}\displaystyle\widetilde{P}_{11}=(P_{11}+P_{12}+P_{21}+P_{22})/2,&\quad\widetilde{P}_{12}=(P_{11}-P_{12}+P_{21}-P_{22})/2,&\quad\widetilde{P}_{13}=P_{13}+P_{23},&\quad\widetilde{P}_{14}=0,\\\
\displaystyle\widetilde{P}_{21}=(P_{11}+P_{12}-P_{21}-P_{22})/2,&\quad\widetilde{P}_{22}=(P_{11}-P_{12}-P_{21}+P_{22})/2,&\quad\widetilde{P}_{23}=P_{13}-P_{23},&\quad\widetilde{P}_{24}=0,\\\
\displaystyle\widetilde{P}_{31}=(P_{31}+P_{32})/2,&\quad\widetilde{P}_{32}=(P_{31}-P_{32})/2,&\quad\widetilde{P}_{33}=P_{33},&\quad\widetilde{P}_{34}=0,\\\
\widetilde{P}_{41}=0,&\quad\widetilde{P}_{42}=0,&\quad\widetilde{P}_{43}=0,&\quad\widetilde{P}_{44}=P_{44}.\end{array}\right.$
(C25)

The updating relationship for the quantity
$\mathbf{\Psi}=(\Psi^{(I)},\Psi^{(Q)},\Psi^{(U)},\Psi^{(V)})^{T}$ from $m$ to
$m+1$ can be obtained directly by multiplying $\mathbf{\Psi}$ with
$\widetilde{\mathbf{P}}$ from the left side, i.e.,

$\displaystyle\left.\begin{array}[]{ll}\displaystyle\Psi^{(I)}_{m+1}=\widetilde{P}_{11}\Psi^{(I)}_{m}+\widetilde{P}_{12}\Psi^{(Q)}_{m}+\widetilde{P}_{13}\Psi^{(U)}_{m},&\displaystyle\Psi^{(Q)}_{m+1}=\widetilde{P}_{21}\Psi^{(I)}_{m}+\widetilde{P}_{22}\Psi^{(Q)}_{m}+\widetilde{P}_{23}\Psi^{(U)}_{m},\\\
\displaystyle\Psi^{(U)}_{m+1}=\widetilde{P}_{31}\Psi^{(I)}_{m}+\widetilde{P}_{32}\Psi^{(Q)}_{m}+\widetilde{P}_{33}\Psi^{(U)}_{m},&\displaystyle\Psi^{(V)}_{m+1}=\widetilde{P}_{44}\Psi^{(V)}_{m}.\end{array}\right.$
(C28)

Or in a more compact form

$\displaystyle\begin{array}[]{ll}\displaystyle\Psi^{(I)}_{m+1}=f_{I}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})\Psi^{(I)}_{m},&\displaystyle\Psi^{(Q)}_{m+1}=f_{Q}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})\Psi^{(I)}_{m},\\\
\displaystyle\Psi^{(U)}_{m+1}=f_{U}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})\Psi^{(I)}_{m},&\displaystyle\Psi^{(V)}_{m+1}=f_{V}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})\Psi^{(V)}_{m},\end{array}$
(C31)

where

$\displaystyle\left\\{\begin{array}[]{ll}f_{I}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})=&\displaystyle\frac{3}{32\pi}\left[2(1-\mu^{2})(1-\mu^{\prime
2})(1+Q_{m})+(1+\mu^{2})(1+\mu^{\prime 2}-(1-\mu^{\prime 2})Q_{m})\right.\\\
&\displaystyle+4\mu\mu^{\prime}\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime
2}}(1+Q_{m})\cos(\varphi^{\prime}-\varphi)+4\mu\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime
2}}U_{m}\sin(\varphi^{\prime}-\varphi)\\\ &\left.+(1-\mu^{2})(1-\mu^{\prime
2})-(1-\mu^{2})(1+\mu^{\prime 2})Q_{m}\cos
2(\varphi^{\prime}-\varphi)-2(1-\mu^{2})\mu^{\prime}U_{m}\sin
2(\varphi^{\prime}-\varphi)\right],\\\
f_{Q}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})=&\displaystyle\frac{3}{32\pi}[(1-\mu^{2})(1-\mu^{\prime
2})(2+3Q_{m})-(1-\mu^{2})(1+\mu^{\prime
2})+4\mu\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime
2}}U_{m}\sin(\varphi^{\prime}-\varphi)\\\
&+4\mu\mu^{\prime}\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime
2}}(1+Q_{m})\cos(\varphi^{\prime}-\varphi)\\\ &+(1+\mu^{2})(\mu^{\prime
2}-1+(1+\mu^{\prime 2})Q_{m})\cos
2(\varphi^{\prime}-\varphi)+2(1+\mu^{2})\mu^{\prime}\sin
2(\varphi^{\prime}-\varphi)],\\\
f_{U}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})=&\displaystyle\frac{3}{32\pi}[4\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime
2}}U_{m}\cos(\varphi^{\prime}-\varphi)-4\mu^{\prime}(1+Q_{m}\mu)\sin(\varphi^{\prime}-\varphi)\\\
&+2\mu\mu^{\prime}U_{m}\cos 2(\varphi^{\prime}-\varphi)-2\mu(1+\mu^{\prime
2}(1+Q_{m}))\sin 2(\varphi^{\prime}-\varphi)],\\\
f_{V}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})=&\displaystyle\frac{3}{8\pi}[\mu\mu^{\prime}+\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime
2}}\cos 2(\varphi^{\prime}-\varphi)],\end{array}\right.$ (C41)

and $Q_{m}=\Psi^{(Q)}_{m}/\Psi^{(I)}_{m}$ and
$U_{m}=\Psi^{(U)}_{m}/\Psi^{(I)}_{m}$. As aforementioned in the Section. 3.3,
we will take $f_{I}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})$ as the PDF to
sample the scattered direction $(\mu,\varphi)$ when the incident direction
$(\mu^{\prime},\varphi^{\prime})$ are provided. We can recast it into a more
compact form given as follows

$\displaystyle\displaystyle
f_{I}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})=\frac{3}{32\pi}(f_{1}+g_{2}\cos\varphi+g_{3}\sin\varphi+g_{4}\cos
2\varphi+g_{5}\sin 2\varphi),$ (C42)

where

$\displaystyle\left\\{\begin{array}[]{l }\displaystyle
f_{1}(\mu,\mu^{\prime})=2(1-\mu^{2})(1-\mu^{\prime
2})(1+Q_{m})+(1+\mu^{2})[1+\mu^{\prime 2}-(1-\mu^{\prime 2})Q_{m}],\\\
\displaystyle
f_{2}(\mu,\mu^{\prime})=4\mu\mu^{\prime}\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime
2}}(1+Q_{m}),\\\ \displaystyle
f_{3}(\mu,\mu^{\prime})=4\mu\sqrt{1-\mu^{2}}\sqrt{1-\mu^{\prime 2}}U_{m},\\\
\displaystyle f_{4}(\mu,\mu^{\prime})=(1-\mu^{2})(1-\mu^{\prime
2})-(1-\mu^{2})(1+\mu^{\prime 2})Q_{m},\\\ \displaystyle
f_{5}(\mu,\mu^{\prime})=-2(1-\mu^{2})\mu^{\prime}U_{m},\\\ \displaystyle
g_{2}=f_{2}\cos\varphi^{\prime}+f_{3}\sin\varphi^{\prime},\quad\quad\quad
g_{3}=f_{2}\sin\varphi^{\prime}-f_{3}\cos\varphi^{\prime},\\\ \displaystyle
g_{4}=f_{4}\cos 2\varphi^{\prime}+f_{5}\sin 2\varphi^{\prime},\quad\quad
g_{5}=f_{4}\sin 2\varphi^{\prime}-f_{5}\cos
2\varphi^{\prime}.\end{array}\right.$ (C50)

Then one can readily show that PDF
$p(\mu,\varphi)=f_{I}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})$ is
automatically normalized. Integrating out $p(\mu,\varphi)$ over $\varphi$, we
get the marginal PDF of $\mu$ as

$\displaystyle\begin{array}[]{l}\displaystyle
p(\mu)=\int_{0}^{2\pi}p(\mu,\varphi)d\varphi=\frac{3}{16}f_{1}(\mu,\mu^{\prime})=\frac{3}{16}(c_{1}\mu^{2}+c_{2}),\end{array}$
(C52)

where $c_{1}=1+\mu^{\prime 2}-(1-\mu^{\prime 2})Q_{m}-2(1-\mu^{\prime
2})(1+Q_{m})$ and $c_{2}=1+\mu^{\prime 2}-(1-\mu^{\prime
2})Q_{m}+2(1-\mu^{\prime 2})(1+Q_{m})$. Hence we can use the inverse CDF
method to get the scattered $\mu$, which is the real root of the equation
$c_{1}\mu^{3}+3c_{2}\mu+c_{1}+3c_{2}-16\xi=0$ on interval $[-1,1]$, where
$\xi$ is a random number. When $\mu$ is obtained, we can immediately get the
PDF for $\varphi$ by using the Bayes formula
$p(\varphi|\mu)=p(\mu,\varphi)/p(\mu)$ as

$\displaystyle\displaystyle p(\varphi|\mu)=\frac{1}{2\pi
f_{1}}(f_{1}+g_{2}\cos\varphi+g_{3}\sin\varphi+g_{4}\cos 2\varphi+g_{5}\sin
2\varphi).$ (C53)

Now the inverse CDF method seems unfeasible to sample $p(\varphi|\mu)$, since
it will give rise to a transcendental equation in terms of $\varphi$ and more
complicated than the Kepler equation. But utilizing the symmetrical and
periodic properties of triangular functions, we have constructed an Algorithm.
4 that can circumvent this difficulty. In order to sample $p(\varphi|\mu)$
through Algorithm. 4, we introduce some relevant quantities given by

$\displaystyle\begin{array}[]{l}\displaystyle F(\varphi)=\frac{1}{2\pi
f_{1}}(f_{1}+g_{2}\cos\varphi+g_{4}\cos 2\varphi),\quad
H(\varphi)=\frac{1}{2\pi f_{1}}(g_{3}\sin\varphi+g_{5}\sin 2\varphi),\\\
\displaystyle p_{0}=1-p_{1}-p_{2},\quad p_{1}=\frac{|g_{2}|}{f_{1}},\quad
p_{2}=\frac{|g_{4}|}{f_{1}},\quad s_{1}=\text{sign}(g_{2}),\quad
s_{2}=\text{sign}(g_{4}).\end{array}$ (C56)

With $\mu$ and $\varphi$ are obtained, we can simply update the components of
scattered $\mathbf{\Psi}$ as

$\displaystyle\begin{array}[]{ll}\displaystyle\Psi^{(I)}_{m+1}=\Psi^{(I)}_{m},&\displaystyle\Psi^{(Q)}_{m+1}=\frac{f_{Q}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})}{f_{I}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})}\Psi^{(I)}_{m},\\\
\displaystyle\Psi^{(U)}_{m+1}=\frac{f_{U}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})}{f_{I}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})}\Psi^{(I)}_{m},&\displaystyle\Psi^{(V)}_{m+1}=\frac{f_{V}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})}{f_{I}(\mu,\varphi;\mu^{\prime},\varphi^{\prime})}\Psi^{(V)}_{m}.\end{array}$
(C59)

##  Appendix D The Analytic Solutions for RT equations with Faraday Rotation
and Conversion

Now we give the analytical solutions for two RT problems including Faraday
effects. In the first case, only the $j_{IQ}$ and $\alpha_{IQ}$ are presented
and the RTE reads as

$\displaystyle\begin{split}\frac{d}{ds}\begin{pmatrix}I\\\
Q\end{pmatrix}=\begin{pmatrix}j_{I}\\\
j_{Q}\end{pmatrix}-\begin{pmatrix}\alpha_{I}&\alpha_{Q}\\\
\alpha_{Q}&\alpha_{I}\end{pmatrix}\begin{pmatrix}I\\\
Q\end{pmatrix},\end{split}$ (D1)

the analytic solution of which can be written as

$\displaystyle\begin{split}I(s)=\frac{j_{+}}{2\alpha_{+}}(1-\text{e}^{-\alpha_{+}s})+\frac{j_{-}}{2\alpha_{-}}(1-\text{e}^{-\alpha_{-}s})+\frac{1}{2}\text{e}^{-\alpha_{+}s}(I_{0}+Q_{0})+\frac{1}{2}\text{e}^{-\alpha_{-}s}(I_{0}-Q_{0}),\\\
Q(s)=\frac{j_{+}}{2\alpha_{+}}(1-\text{e}^{-\alpha_{+}s})-\frac{j_{-}}{2\alpha_{-}}(1-\text{e}^{-\alpha_{-}s})+\frac{1}{2}\text{e}^{-\alpha_{+}s}(I_{0}+Q_{0})-\frac{1}{2}\text{e}^{-\alpha_{-}s}(I_{0}-Q_{0}),\end{split}$
(D2)

where $j_{\pm}=j_{I}\pm j_{Q}$, $\alpha_{\pm}=\alpha_{I}\pm\alpha_{Q}$, and
$I_{0},Q_{0}$ are the initial values of $I$ and $Q$ at $s=0$. In the second
case, all of the absorption coefficients $\alpha_{IQUV}$ and emissivity
$j_{I}$ vanish, the RTE can be written as

$\displaystyle\frac{d}{ds}\begin{pmatrix}Q\\\ U\\\
V\end{pmatrix}=\begin{pmatrix}j_{Q}\\\ j_{U}\\\
j_{V}\end{pmatrix}-\begin{pmatrix}0&\rho_{V}&-\rho_{U}\\\
-\rho_{V}&0&\rho_{Q}\\\ \rho_{U}&-\rho_{Q}&0\\\
\end{pmatrix}\begin{pmatrix}Q\\\ U\\\ V\end{pmatrix}.$ (D3)

Integrating the above equations directly, we have (Mościbrodzka & Gammie,
2018)

$\displaystyle\left\\{\begin{split}Q(s)=&\frac{\rho_{Q}}{\rho^{2}}\rho\cdot\boldsymbol{J}s+(j_{Q}\rho^{2}-\rho_{Q}\rho\cdot\boldsymbol{J})\frac{\sin(\rho
s)}{\rho^{3}}+(\rho_{U}j_{V}-\rho_{V}j_{U})\frac{1-\cos(\rho s)}{\rho^{2}}\\\
&+Q_{0}\cos(\rho
s)+2\frac{\rho_{Q}(\rho\cdot\boldsymbol{S})}{\rho^{2}}\sin^{2}(\rho
s/2)+\frac{\rho_{U}V_{0}-\rho_{V}U_{0}}{\rho}\sin(\rho s),\\\
U(s)=&\frac{\rho_{U}}{\rho^{2}}\rho\cdot\boldsymbol{J}s+(j_{U}\rho^{2}-\rho_{U}\rho\cdot\boldsymbol{J})\frac{\sin(\rho
s)}{\rho^{3}}+(\rho_{V}j_{Q}-\rho_{Q}j_{V})\frac{1-\cos(\rho s)}{\rho^{2}}\\\
&+U_{0}\cos(\rho
s)+2\frac{\rho_{U}(\rho\cdot\boldsymbol{S})}{\rho^{2}}\sin^{2}(\rho
s/2)+\frac{\rho_{V}Q_{0}-\rho_{Q}V_{0}}{\rho}\sin(\rho s),\\\
V(s)=&\frac{\rho_{V}}{\rho^{2}}\rho\cdot\boldsymbol{J}s+(j_{V}\rho^{2}-\rho_{V}\rho\cdot\boldsymbol{J})\frac{\sin(\rho
s)}{\rho^{3}}+(\rho_{Q}j_{U}-\rho_{U}j_{Q})\frac{1-\cos(\rho s)}{\rho^{2}}\\\
&+V_{0}\cos(\rho
s)+2\frac{\rho_{V}(\rho\cdot\boldsymbol{S})}{\rho^{2}}\sin^{2}(\rho
s/2)+\frac{\rho_{Q}U_{0}-\rho_{U}Q_{0}}{\rho}\sin(\rho s),\end{split}\right.$
(D4)

where $\rho^{2}=\rho_{Q}^{2}+\rho_{U}^{2}+\rho_{V}^{2}$,
$\rho\cdot\boldsymbol{S}=\rho_{Q}Q_{0}+\rho_{U}U_{0}+\rho_{V}V_{0}$,
$\mathbf{\rho}\cdot\boldsymbol{J}=\rho_{Q}j_{Q}+\rho_{U}j_{U}+\rho_{V}j_{V}$,
and $Q_{0},U_{0},V_{0}$ are the initial values of $Q,U,V$ at $s=0$.

##  Appendix E The Neumann solution as the disk reflection effect considered

Now we discuss what should the sampling procedure be modified when the
reflection term $\mathbf{I}_{r}$ given by Eq. (151) is incorporated. Although
the discussion is presented in a plane-parallel geometrical system, it is
possible to extend it to any other systems. Since the reflected radiations can
reenter the corona, they can be regarded as the seed photons that will
initiate another round of RT process. It means that the RT problem now becomes
an iterative one. Assume that $\mathbf{J}(\tau_{0},\mathbf{\Omega},\nu)$ is
the primary radiations emitted by the disk, the corresponding radiative flux
generated by $\mathbf{J}$ is $\mathbf{\Psi}^{(0)}$, which obviously satisfies
the integral equation

$\displaystyle\begin{split}&\mathbf{\Psi}^{(0)}(P)=\mathbf{J}(P)+\int\mathbf{K}(P^{\prime}\rightarrow
P)\mathbf{\Psi}^{(0)}(P^{\prime})dP^{\prime}.\end{split}$ (E1)

Their Neumann solutions are

$\displaystyle\begin{split}&\mathbf{\Psi}^{(0)}(P)=\sum_{m=0}^{\infty}\mathbf{\Psi}^{(0)}_{m}(P),\quad\mathbf{\Psi}^{(0)}_{0}(P)=\mathbf{J}(P),\quad\mathbf{\Psi}^{(0)}_{1}(P)=\int\mathbf{K}(P^{\prime}\rightarrow
P)\mathbf{J}(P^{\prime})dP^{\prime}=\mathbf{J}\mathbf{K},\cdots,\\\
&\mathbf{\Psi}^{(0)}_{m}(P)=\int\mathbf{K}(P_{m-1}\rightarrow
P)\mathbf{K}(P_{m-2}\rightarrow P_{m-1})\cdots\mathbf{K}(P_{0}\rightarrow
P_{1})\mathbf{J}(P_{1})dP_{0}\cdots
dP_{m-1}=\mathbf{J}\mathbf{K}\cdots\mathbf{K}=\mathbf{J}\mathbf{K}^{m}.\end{split}$
(E2)

For simplicity, we denote the multiple integrals in terms of any $P$ in a
compact form. Using $\mathbf{\Psi}^{(0)}$, we can immediately obtain the
radiations escaped from the up and down surfaces of the corona respectively as

$\displaystyle\begin{split}\begin{array}[]{lr}\displaystyle\mathbf{I}^{(0)}_{\text{up}}(0,\mathbf{\Omega}^{\prime},\nu^{\prime})=\sum_{m=0}^{\infty}\mathbf{I}^{(0)}_{\text{up},m}=\sum_{m=0}^{\infty}\int_{0}^{\tau_{0}}\mathbf{J}\mathbf{K}^{m}\exp\left(-\frac{\tau^{\prime}}{\mu^{\prime}}\right)\frac{d\tau^{\prime}}{\mu^{\prime}},&\quad\mu^{\prime}>0,\\\
\displaystyle\mathbf{I}^{(0)}_{\text{dw}}(\tau_{0},\mathbf{\Omega}^{\prime},\nu^{\prime})=\sum_{m=0}^{\infty}\mathbf{I}^{(0)}_{\text{dw},m}=\sum_{m=0}^{\infty}\int_{0}^{\tau_{0}}\mathbf{J}\mathbf{K}^{m}\exp\left(-\frac{\tau_{0}-\tau^{\prime}}{|\mu^{\prime}|}\right)\eta(-\mu^{\prime})\frac{d\tau^{\prime}}{|\mu^{\prime}|},&\quad\mu^{\prime}<0,\end{array}\end{split}$
(E3)

where $\tau_{0}$ is the optical depth of the corona. Substituting
$\mathbf{I}^{(0)}_{\text{dw}}$ into Eq. (151) gives the radiation intensity
reflected from the disk as

$\displaystyle\begin{split}\mathbf{I}^{(1)}_{r}(\tau_{0},\mathbf{\Omega},\nu)=\sum_{m=0}^{\infty}\mathbf{I}^{(1)}_{r,m}(\tau_{0},\mathbf{\Omega},\nu),\end{split}$
(E4)

where

$\displaystyle\begin{split}\mathbf{I}^{(1)}_{r,m}(\tau_{0},\mathbf{\Omega},\nu)=\int\int_{0}^{\tau_{0}}\mathbf{G}(\nu^{\prime}\rightarrow\nu,\mathbf{\Omega}^{\prime}\rightarrow\mathbf{\Omega})\mathbf{J}\mathbf{K}^{m}\displaystyle
p(\tau_{0},\tau^{\prime})\eta(-\mu^{\prime})\frac{d\tau^{\prime}}{|\mu^{\prime}|}d\mathbf{\Omega}^{\prime}d\nu^{\prime}=\mathbf{J}\mathbf{K}^{m}p\eta\mathbf{G},\end{split}$
(E5)

and $\displaystyle
p(\tau_{0},\tau^{\prime})=\exp\left(-\frac{\tau_{0}-\tau^{\prime}}{|\mu^{\prime}|}\right)$.
Now each $\mathbf{I}^{(1)}_{r,m}(\tau_{0},\mathbf{\Omega},\nu)$ can be taken
as the incident radiation from the bottom surface of the corona, and the
corresponding $\mathbf{\Psi}^{(1)}_{m}$ satisfies:

$\displaystyle\begin{split}&\mathbf{\Psi}^{(1)}_{m}(P)=\mathbf{I}^{(1)}_{r,m}\mathbf{K}(P)+\int\mathbf{K}(P^{\prime}\rightarrow
P)\mathbf{\Psi}^{(1)}_{m}(P^{\prime})dP^{\prime}.\end{split}$ (E6)

The Neuman solution of $\mathbf{\Psi}^{(1)}_{m}$ is
$\mathbf{\Psi}^{(1)}_{m,n}=\mathbf{I}^{(1)}_{r,m}\mathbf{K}\cdot\mathbf{K}^{n}=\mathbf{J}\mathbf{K}^{m}p\eta\mathbf{G}\mathbf{K}\cdot\mathbf{K}^{n}$.
Similarly we can get that the reflected radiations generated by
$\mathbf{\Psi}^{(1)}_{m,n}$ is
$\mathbf{I}^{(2)}_{r}(\tau_{0},\mathbf{\Omega},\nu)=\mathbf{J}\mathbf{K}^{m}p\eta\mathbf{G}\mathbf{K}\cdot\mathbf{K}^{n}p\eta\mathbf{G}$
and its Neuman solution is given by
$\mathbf{\Psi}^{(2)}_{m,n,l}=\mathbf{J}\mathbf{K}^{m}p\eta\mathbf{G}\mathbf{K}\cdot\mathbf{K}^{n}p\eta\mathbf{G}\mathbf{K}\cdot\mathbf{K}^{l}$.
We can repeat this procedure to obtain the Neumann solution
$\mathbf{\Psi}^{(i)}_{m,n,l,\cdots}$ for radiations reflected for arbitrary
times. Then the total Neumann solution for this reflection RT process can be
written as

$\displaystyle\left\\{\begin{split}\mathbf{\Psi}^{(0)}&=\sum_{m=0}^{\infty}\mathbf{JK}\cdot\mathbf{K}^{m}=\mathbf{JK}+\mathbf{JK}\cdot\mathbf{K}+\cdots\mathbf{JK}\cdot\mathbf{K}^{m}+\cdots\\\
\mathbf{\Psi}^{(1)}&=\sum_{m=0,n=0}^{\infty}\mathbf{JK}\cdot\mathbf{K}^{m}\cdot(p\eta\mathbf{GK})\cdot\mathbf{K}^{n},\\\
&=\mathbf{JK}(p\eta\mathbf{GK})+\mathbf{JK}(p\eta\mathbf{GK})\mathbf{K}+\mathbf{JK}(p\eta\mathbf{GK})\mathbf{K}^{2}+\cdots\mathbf{JK}(p\eta\mathbf{GK})\mathbf{K}^{m}+\cdots\\\
&+\mathbf{JK}^{2}(p\eta\mathbf{GK})+\mathbf{JK}^{2}(p\eta\mathbf{GK})\mathbf{K}+\mathbf{JK}^{2}(p\eta\mathbf{GK})\mathbf{K}^{2}+\cdots\mathbf{JK}^{2}(p\eta\mathbf{GK})\mathbf{K}^{m}+\cdots\\\
&\cdots\\\
&+\mathbf{JK}^{m}(p\eta\mathbf{GK})+\mathbf{JK}^{m}(p\eta\mathbf{GK})\mathbf{K}+\mathbf{JK}^{m}(p\eta\mathbf{GK})\mathbf{K}^{2}+\cdots+\mathbf{JK}^{m}(p\eta\mathbf{GK})\mathbf{K}^{m}+\cdots\\\
&\cdots,\\\
\mathbf{\Psi}^{(2)}&=\sum_{m=0,n=0,l=0}^{\infty}\mathbf{JK}\cdot\mathbf{K}^{m}\cdot(p\eta\mathbf{GK})\cdot\mathbf{K}^{n}\cdot(p\eta\mathbf{GK})\cdot\mathbf{K}^{l},\\\
\mathbf{\Psi}^{(3)}&=\sum_{m=0,n=0,l=0,k=0}^{\infty}\mathbf{JK}\cdot\mathbf{K}^{m}\cdot(p\eta\mathbf{GK})\cdot\mathbf{K}^{n}\cdot(p\eta\mathbf{GK})\cdot\mathbf{K}^{l}\cdot(p\eta\mathbf{GK})\cdot\mathbf{K}^{k},\cdots,\end{split}\right.$
(E7)

and the total Neumann solution is
$\mathbf{\Psi}=\mathbf{\Psi}^{(0)}+\mathbf{\Psi}^{(1)}+\mathbf{\Psi}^{(3)}+\cdots+\mathbf{\Psi}^{(m)}+\cdots$.
Combining those terms with equal power index of $\mathbf{K}$,
$\mathbf{\Psi}(P)$ can be reexpressed as:

$\displaystyle\begin{split}\mathbf{\Psi}(P)&=\sum_{k=0}^{\infty}\mathbf{J}\mathbf{K}\cdot(\mathbf{K}+p\eta\mathbf{GK})^{k},\end{split}$
(E8)

where

$\displaystyle\begin{split}(\mathbf{K}+p\eta\mathbf{GK})^{2}=&\mathbf{K}^{2}+\mathbf{K}(p\eta\mathbf{GK})+(p\eta\mathbf{GK})\mathbf{K}+(p\eta\mathbf{GK})^{2},\\\
(\mathbf{K}+p\eta\mathbf{GK})^{3}=&\mathbf{K}^{3}+\mathbf{K}^{2}(p\eta\mathbf{GK})+\mathbf{K}(p\eta\mathbf{GK})\mathbf{K}+(p\eta\mathbf{GK})\mathbf{K}^{2}+\mathbf{K}(p\eta\mathbf{GK})^{2}\\\
&+(p\eta\mathbf{GK})\mathbf{K}(p\eta\mathbf{GK})+(p\eta\mathbf{GK})^{2}\mathbf{K}+(p\eta\mathbf{GK})^{3},\cdots.\\\
\end{split}$ (E9)

Of course we have

$\displaystyle\begin{split}&\mathbf{J}\mathbf{K}(P)=\int\mathbf{K}(P_{0}\rightarrow
P)\mathbf{J}(P_{0})dP_{0},\\\
&\mathbf{J}\mathbf{K}\cdot(\mathbf{K}+p\eta\mathbf{GK})(P)=\mathbf{J}\mathbf{K}\mathbf{K}(P)+\mathbf{J}\mathbf{K}p\eta\mathbf{GK}(P)=\int\mathbf{K}(P_{1}\rightarrow
P)\int\mathbf{K}(P_{0}\rightarrow P_{1})\mathbf{J}(P_{0})dP_{0}dP_{1}\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\int\mathbf{K}(P^{\prime}_{1}\rightarrow
P)\mathbf{G}(P_{1}\rightarrow
P^{\prime}_{1})p(\tau_{0},\tau_{1})\eta(-\mu_{1})\int\mathbf{K}(P_{0}\rightarrow
P_{1})\mathbf{J}(P_{0})dP_{0}dP^{\prime}_{1}dP_{1},\cdots.&\end{split}$ (E10)

Comparing to the case without reflection, the Neumann solution is almost the
same except that the transport kernel now becomes
$\mathbf{K}+p\eta\mathbf{GK}$. For upward propagating radiations, we always
have $\eta(-\mu)\equiv 0$, then the transport kernel is $\mathbf{K}$. Thus
only the cases with $\mu<0$ need to be considered and the transport kernel
becomes $\mathbf{K}+p\mathbf{GK}$. Suppose that the initial position,
direction and frequency of the radiation are
$\tau^{\prime},\,\mathbf{\Omega}^{\prime}=(\mu^{\prime},\varphi^{\prime}),\,\nu^{\prime}$
(notice $\mu^{\prime}<0$), the transport kernel is given by
$\mathbf{K}=\mathbf{C}(\nu^{\prime}\rightarrow\nu,\mathbf{\Omega}^{\prime}\rightarrow\mathbf{\Omega})\cdot
T_{s}(\tau^{\prime}\rightarrow\tau)\cdot w_{s}$, where $\mathbf{C}$ and
$T_{s}=\exp[-(\tau-\tau^{\prime})/|\mu^{\prime}|]/A$ are the scattering and
position transport kernels respectively, and
$w_{s}=A=1-p(\tau_{0},\tau^{\prime})$ is the weight. Thus we have

$\displaystyle\begin{split}\mathbf{K}+p\mathbf{GK}=&T_{s}(\tau^{\prime}\rightarrow\tau)\mathbf{C}(\nu^{\prime}\rightarrow\nu,\mathbf{\Omega}^{\prime}\rightarrow\mathbf{\Omega}|\tau)\cdot(1-p)\\\
+&\mathbf{G}(\nu^{\prime}\rightarrow\nu^{\prime\prime},\mathbf{\Omega}^{\prime}\rightarrow\mathbf{\Omega}^{\prime\prime}|\tau_{0})T_{s}(\tau_{0}\rightarrow\tau)w^{\prime\prime}_{s}\mathbf{C}(\nu^{\prime\prime}\rightarrow\nu,\mathbf{\Omega}^{\prime\prime}\rightarrow\mathbf{\Omega}|\tau)\cdot
p,\end{split}$ (E11)

where $w^{\prime\prime}_{s}=1-\exp(-\tau_{0}/\mu^{\prime\prime})$. From them,
we can obtain the sampling procedure for the transport with reflection. We
first generate a random variable $\xi$, if $\xi\leq p$, then the radiation
will be reflected and we sample $\mathbf{G},T_{s}(\tau_{0}\rightarrow\tau)$
and $\mathbf{C}$ in turn to get the new point: $(\tau,\mathbf{\Omega},\nu)$;
otherwise if $\xi>p$ we sample $T_{s}(\tau^{\prime}\rightarrow\tau)$ and
$\mathbf{C}$ instead. This is exactly the revised sampling procedure as
reflection effects are considered in our MC scheme.

From the above discussions one can see that the sampling procedure can be
obtained correctly from mathematical deductions rather than physical
intuitions. It can help us to avoid some unnecessary mistakes.

## References

  * Akhiezer & Berestetskii (1969) Akhiezer, A. I., Berestetskii, V. B. 1969, Quantum Electrodynamics (Nauka, Moscow)
  * Bonometto et al. (1970) Bonometto, S., Cazzola, P., & Saggion, A. 1970, A&A, 7, 292
  * Bonometto & Saggion (1973a) Bonometto, S. A. & Saggion, A. 1973, A&A, 23, 9
  * Bonometto & Saggion (1973b) Bonometto, S. A. & Saggion, A. 1973, Astrophys. Lett., 13, 193
  * Böttcher & Liang (2001) Böttcher, M. & Liang, E. P. 2001, ApJ, 552, 248
  * Broderick & Blandford (2003) Broderick, A. & Blandford, R. 2003, MNRAS, 342, 1280
  * Broderick & Blandford (2004) Broderick, A. & Blandford, R. 2004, MNRAS, 349, 994
  * Bronzwaer et al. (2018) Bronzwaer, T., Davelaar, J., Younsi, Z., et al. 2018, A&A, 613, A2
  * Bronzwaer et al. (2020) Bronzwaer, T., Younsi, Z., Davelaar, J., et al. 2020, arXiv:2007.03045
  * Canfield et al. (1987) Canfield, E., Howard, W. M., Liang, E. P. 1987, ApJ, 323, 565
  * Chen et al. (2015) Chen, B., Kantowski, R., Dai, X., et al. 2015, ApJS, 218, 4
  * Chan et al. (2019) Chan, J. Y. H., Wu, K., On, A. Y. L., et al. 2019, MNRAS, 484, 1427
  * Chandrasekhar (1960) Chandrasekhar, S. 1960, Radiative Transfer (New York: Dover)
  * Chandrasekhar (1983) Chandrasekhar, S. 1983, The Mathematical Theory of Black Holes (New York: Oxford Univ. Press)
  * Connors & Stark (1977) Connors, P. A. & Stark, R. F. 1977, Nature, 269, 128
  * Connors et al. (1980) Connors, P. A., Piran, T., Stark, R. F. 1980, ApJ, 235, 224
  * Cunningham (1975) Cunningham, C. T. 1975, ApJ, 202, 788
  * Cunningham (1976) Cunningham, C. 1976, ApJ, 208, 534
  * Dauser et al. (2010) Dauser, T., Wilms, J., Reynolds, C. S., et al. 2010, MNRAS, 409, 1534. doi:10.1111/j.1365-2966.2010.17393.x
  * Davison (1957) Davison, B. 1957, Neutron Transport Theory (OXFORD: Clarendon Press)
  * Dexter & Agol (2009) Dexter, J. & Agol, E. 2009, ApJ, 696, 1616
  * Dexter et al. (2010) Dexter, J., Agol, E., Fragile, P. C., et al. 2010, ApJ, 717, 1092
  * Dexter (2016) Dexter, J. 2016, MNRAS, 462, 115
  * Dexter et al. (2020) Dexter, J., Jiménez-Rosales, A., Ressler, S. M., et al. 2020, MNRAS, 494, 4168
  * Dolence et al. (2009) Dolence, J. C., Gammie, C. F., Mościbrodzka, M., et al. 2009, ApJS, 184, 387
  * Dovčiak et al. (2004) Dovčiak, M., Karas, V., & Yaqoob, T. 2004, ApJS, 153, 205
  * Dovčiak et al. (2008) Dovčiak, M., Muleri, F., Goosmann, R. W., et al. 2008, MNRAS, 391, 32
  * Dovčiak et al. (2011) Dovčiak, M., Muleri, F., Goosmann, R. W., et al. 2011, ApJ, 731, 75
  * Fano (1949) Fano, U. 1949, Journal of the Optical Society of America (1917-1983), 39, 859
  * Fano (1957) Fano, U. 1957, Reviews of Modern Physics, 29, 74
  * Gammie & Leung (2012) Gammie, C. F. & Leung, P. K. 2012, ApJ, 752, 123
  * Goertzel (1958) Goertzel, G., Kalos, M. H. 1958, Progress in Nuclear Energy, Series I, p. 315-369, Pergamon Press, New York
  * Goosmann & Gaskell (2007) Goosmann, R. W., Gaskell, C. M. 2007, A&A, 465, 129, doi: 10.1051/0004-6361:20053555
  * Gorecki & Wilczewski (1984) Gorecki, A. & Wilczewski, W. 1984, Acta Astron., 34, 141
  * Haardt (1993) Haardt, F. 1993, ApJ, 413, 680
  * Hauschildt & Wehrse (1991) Hauschildt, P. H. & Wehrse, R. 1991, J. Quant. Spec. Radiat. Transf., 46, 81
  * Hua & Lingenfelter (1992) Hua, X.-M. & Lingenfelter, R. E. 1992, ApJ, 397, 591
  * Hua & Titarchuk (1995) Hua, X.-M. & Titarchuk, L. 1995, ApJ, 449, 188
  * Hua (1997) Hua, X. M. 1997, Comput. Phys., 11, 660
  * Huang & Shcherbakov (2011) Huang, L., Shcherbakov, R. V., 2011, MNRAS, 416, 2574
  * Janett et al. (2017a) Janett, G., Carlin, E. S., Steiner, O., et al. 2017, ApJ, 840, 107. doi:10.3847/1538-4357/aa671d
  * Janett et al. (2017b) Janett, G., Steiner, O., & Belluzzi, L. 2017, ApJ, 845, 104. doi:10.3847/1538-4357/aa7aa3
  * Janett (2019) Janett, G. 2019, A&A, 622, A162. doi:10.1051/0004-6361/201833984
  * Jiménez-Rosales & Dexter (2018) Jiménez-Rosales, A. & Dexter, J. 2018, MNRAS, 478, 1875
  * Kahn (1950) Kahn, H. 1950, Nucleonics, 6, 27
  * Kahn (1954) Kahn, H. Applications of Monte Carlo, AECU-3259 1954
  * Kojima (1991) Kojima, Y. 1991, MNRAS, 250, 629
  * Kosowsky (1996) Kosowsky, A. 1996, Annals of Physics, 246, 49
  * Krawczynski (2012) Krawczynski, H. 2012, ApJ, 754, 133
  * Laor et al. (1990) Laor, A., Netzer, H., & Piran, T. 1990, MNRAS, 242, 560
  * Lewis & Bridle (2002) Lewis, A. & Bridle, S. 2002, Phys. Rev. D, 66, 103511. doi:10.1103/PhysRevD.66.103511
  * Li et al. (2009) Li, L.-X., Narayan, R., & McClintock, J. E. 2009, ApJ, 691, 847
  * Lindquist (1966) Lindquist, R. W. 1966, Annals of Physics, 37, 487
  * López Ariste & Semel (1999) López Ariste, A. & Semel, M. 1999, A&A, 350, 1089
  * Magdziarz & Zdziarski (1995) Magdziarz, P. & Zdziarski, A. A. 1995, MNRAS, 273, 837
  * Meliani et al. (2017) Meliani, Z., Mizuno, Y., Olivares, H., et al. 2017, A&A, 598, A38
  * Morrison & McCammon (1983) Morrison, R., & McCammon, D. 1983, ApJ, 270, 119
  * Mościbrodzka & Gammie (2018) Mościbrodzka, M. & Gammie, C. F. 2018, MNRAS, 475, 43
  * Mościbrodzka (2020) Mościbrodzka, M. 2020, MNRAS, 491, 4807
  * Nagirner & Poutanen (1993) Nagirner, D. I. & Poutanen, J. 1993, A&A, 275, 325
  * Noebauer & Sim (2019) Noebauer, U. M. & Sim, S. A. 2019, Living Reviews in Computational Astrophysics, 5, 1. doi:10.1007/s41115-019-0004-9
  * Pihajoki et al. (2017) Pihajoki, P., Rantala, A., & Johansson, P. H. 2017, New Frontiers in Black Hole Astrophysics, 324, 347
  * Pihajoki et al. (2018) Pihajoki, P., Mannerkoski, M., Nättilä, J., et al. 2018, ApJ, 863, 8
  * Pomraning (1973) Pomraning, G. C. 1973, The Equations of Radiation Hydrodynamics (OXFORD. NEW YORK. TORONTO: PERGAMON Press)
  * Portsmouth & Bertschinger (2004) Portsmouth, J., Bertschinger, E., 2004, arXiv:astro-ph/0412094
  * Poutanen & Vilhu (1993) Poutanen, J. & Vilhu, O. 1993, A&A, 275, 337
  * Poutanen & Svensson (1996a) Poutanen, J., & Svensson, R. 1996a, ApJ, 470, 249
  * Poutanen et al. (1996b) Poutanen, J., Nagendra, K. N., & Svensson, R. 1996b, MNRAS, 283, 892
  * Pozdnyakov et al. (1983) Pozdnyakov, L. A., Sobol, I. M., Syunyaev, R. A. 1983, ASPRv, 2, 189
  * Psaltis & Johannsen (2012) Psaltis, D. & Johannsen, T. 2012, ApJ, 745, 1
  * Rauch & Blandford (1994) Rauch, K. P. & Blandford, R. D. 1994, ApJ, 421, 46
  * Ryan et al. (2015) Ryan, B. R., Dolence, J. C., & Gammie, C. F. 2015, ApJ, 807, 31
  * Ryan & Dolence (2020) Ryan, B. R. & Dolence, J. C. 2020, ApJ, 891, 118
  * Shcherbakov & Huang (2011) Shcherbakov, R. V., Huang, L., 2011, MNRAS, 410, 1052
  * Schnittman et al. (2006) Schnittman, J. D., Krolik, J. H., & Hawley, J. F. 2006, ApJ, 651, 1031
  * Schnittman & Krolik (2010) Schnittman, J. D., & Krolik, J. H. 2010, ApJ, 712, 908, doi: 10.1088/0004-637X/712/2/908
  * Schnittman & Krolik (2013) Schnittman, J. D., & Krolik, J. H. 2013, ApJ, 777, 11
  * Semel & López Ariste (1999) Semel, M. & López Ariste, A. 1999, A&A, 342, 201
  * Shcherbakov & Huang (2011) Shcherbakov, R. V. & Huang, L. 2011, MNRAS, 410, 1052
  * Spanier (1959) Spanier, J. 1959, Monte Carlo Methods and Their Application to Neutron Transport Problems (WAPD-195)
  * Stern et al. (1995) Stern, B. E., Begelman, M. C., Sikora, M., et al. 1995, MNRAS, 272, 291
  * Synge (1957) Synge, J. L. 1957, The Reletivistic Gas (Amsterdam: North-Holland).
  * Sunyaev & Titarchuk (1985) Sunyaev, R. A. & Titarchuk, L. G. 1985, A&A, 143, 374
  * Takahashi & Umemura (2017) Takahashi, R. & Umemura, M. 2017, MNRAS, 464, 4567. doi:10.1093/mnras/stw2479
  * Tsunetoe et al. (2020) Tsunetoe, Y., Mineshige, S., Ohsuga, K., et al. 2020, PASJ, 72, 32
  * Vincent et al. (2011) Vincent, F. H., Paumard, T., Gourgoulhon, E., et al. 2011, Classical and Quantum Gravity, 28, 225011
  * Vincent et al. (2020) Vincent, F. H., Wielgus, M., Abramowicz, M. A., et al. 2020, arXiv:2002.09226
  * Walker & Penrose (1970) Walker, M. & Penrose, R. 1970, Communications in Mathematical Physics, 18, 265
  * Weinberg (2008) Weinberg, S. 2008, Cosmology, by Steven Weinberg. ISBN 978-0-19-852682-7. Published by Oxford University Press, Oxford, UK, 2008.
  * Whitney (2011) Whitney, B. A. 2011, Bulletin of the Astronomical Society of India, 39, 101
  * Wienke (1985) Wienke, B. 1985, Astr. Ap., 152, 336
  * Yang & Wang (2013) Yang, X. & Wang, J. 2013, ApJS, 207, 6
  * Younsi et al. (2012) Younsi, Z., Wu, K., & Fuerst, S. V. 2012, A&A, 545, A13
  * Younsi & Wu (2013) Younsi, Z. & Wu, K. 2013, MNRAS, 433, 1054
  * Yuan et al. (2009) Yuan, Y.-F., Cao, X., Huang, L., et al. 2009, ApJ, 699, 722
  * Zane et al. (1996) Zane, S., Turolla, R., Nobili, L., et al. 1996, ApJ, 466, 871
  * Zhang et al. (2019) Zhang Wen-da, J. D., Krolik, J. H. 2019, ApJ, 777, 11