A Relativistic Formula for the Multiple Scattering of Photons

The relativistic effects of radiation field are essential in a variety of astrophysical flow and jets, e.g., active galactic nuclei, gamma-ray bursts, and supernovae. In a black hole accretion flow with super-Eddington accretion rate, the photons and the fluid are tightly coupled and the radiative transfer (RT) of the scattered photons in a photon-trapping region is of great significance (Ohsuga & Mineshige 2007; Ohsuga et al. 2009; McKinney et al. 2014; Takeo et al. 2020; Liska et al. 2022). In order to investigate the dynamical radiation effects in a relativistic flow, some past studies tried to directly solve the RT equation, i.e., Boltzmann equation, for photons (Beloborodov 2011; Jiang et al. 2016; Ohsuga & Takahashi 2016; Takahashi & Umemura 2017; Asahina et al. 2020; Asahina & Ohsuga 2022; Takahashi et al. 2022) and neutrinos (e.g., Nagakura et al. 2018; Akaho et al. 2021, 2023).

Despite the significant efforts of the past studies, even for the simplest cases assuming isotropic and elastic scattering of photons, the time-dependent radiation field of photons (or neutrinos) in a relativistic flow has not been hitherto solved exactly. In a relativistic flow, the relativistic boosting effect is significant in the diffusion process (Krumholz et al. 2007; Shibata et al. 2014). In particular, a precise description of the intermediate state between the free-propagating state and the diffusion state has not been possible to date, and lack of analytical understanding has prevented a fundamental understanding of these phenomena. In fact, we have confirmed that some of the approximate formulas proposed in the past contain a violation of the law of causality. Additionally, the relativistic diffusion problem is a well-known unsolved problem that has not been solved for many years (e.g., Dunkel et al. 2007). The probability density function (PDF) provided in this Letter, describing relativistic diffusion for photons, likely represents the first instance of an analytical solution for the relativistic diffusion problem involving these particles.

In this Letter, we shall pursue the analytic approach to describe the collective behavior of the repeatedly scattered photons in a relativistically boosted medium, and we investigate the time evolution of the scattering photons in the relativistic flow both in the rest frame (Section 2) and in the laboratory frame (Section 3). In the final section (Section 4), the conclusions are given.

According to relativistic kinetic theory (e.g., Lindquist 1966; Ehlers 1971; Sachs & Ehlers 1971; Israel 1972), the collective behavior of particles is described by an invariant distribution function ${ \mathcal F }({x}^{\mu },{p}^{\mu })$ where x^μ and p^μ are, respectively, the coordinate and the momentum of a particle. The particle number density flux N^μ is given by ${N}^{\mu }=\int {{dN}}^{\mu }=\int { \mathcal F }({x}^{\mu },{p}^{\mu }){p}^{\mu }{dP}$ and the particle number N(x^μ )dV in a volume element dV at a time slice surface (where t = t₀) is given by the projection of N^μ onto the vector volume element, i.e., $N({x}^{\mu }){dV}\,={dV}\int { \mathcal F }({x}^{\mu },{p}^{\mu })(-{\hat{u}}_{\mu }{p}^{\mu }){dP}$ where ${\hat{u}}_{\mu }$ is a time-like unit vector. In this study, the PDF in a time slice surface are used to describe the collective behavior of scattered photons. The PDF $P({x}^{\mu }){| }_{t={t}_{0}}$ giving the probability that particles exist in some region at the time slice surface by volume integration is given by $P({x}^{\mu }){| }_{t={t}_{0}}=N({x}^{\mu })/{N}_{\mathrm{all}}$ where N_all is the number of photons in three space V at t = t₀ calculated as ${N}_{\mathrm{all}}={\int }_{{V}_{t={t}_{0}}}N({x}^{\mu }){dV}$. Then, the PDF P(t, r) satisfies the normalization ${\int }_{{V}_{t={t}_{0}}}P({x}^{\mu }){dV}=1$. The time component of N^μ (x^μ ) as N⁰(x^μ ) = N_all P(x^μ ).

We describe the coordinates x^μ and momentum k^μ of a photon as x^μ = (ct, r ) and k^μ = (ℏ ω/c, ℏ k ) where c is the speed of light and ℏ is the Planck constant divided by 2π. It is useful to introduce the normalization with the mean free path of photon ℓ₀ in the fluid rest frame (Rybicki & Lightman1979; Mihalas & Mihalas 2013). Here, we introduce ${x}_{* }^{\mu }=({t}_{* },{{\boldsymbol{r}}}_{* })\equiv {x}^{\mu }/{{\ell }}_{0}=({ct}/{{\ell }}_{0},{\boldsymbol{r}}/{{\ell }}_{0})$ and ${k}_{* }^{\mu }=({\omega }_{* },{{\boldsymbol{k}}}_{* })\,\equiv {k}^{\mu }{{\ell }}_{0}=({\hslash }\omega {{\ell }}_{0}/c,{\hslash }{\boldsymbol{k}}{{\ell }}_{0})$. We also define r_* ≡ ∣ r _*∣ and k_* ≡ ∣ k _*∣.

The PDF, P(t_*, r _*), of a scattered photon is a function where P(t_*, r _*)d r _* gives the probability that the scattered photon is in spatial region d r _* in the spatial hypersurface at time t_*. These scattered photons include photons that have experienced n times scatterings where n = 0, 1, 2, ⋯ . We have found the PDF, P(t_*, r _*), is expressed as the sum P(t_*, r _*) = ${\sum }_{n=0}^{\infty }{P}_{n}({t}_{* },{{\boldsymbol{r}}}_{* })$ where P_n (t_*, r _*) is a PDF in four-dimensional space-time that represents the distribution of the next scattering point of photons scattered n times, which satisfy the normalization ∫P_n (t_*, r _*)d⁴ x_* = 1 where d⁴ x_* is an invariant volume element of four-dimensional space-time. ⁷ In the following, we report the method and results of analytic calculation of P_n (t_*, r _*) under the assumptions of isotropic elastic scattering.

Consider a process in which, at time t_* = 0, a large number of photons are instantaneously emitted isotropically from a point in a static medium and spread out while repeatedly elastic scattering in the medium. Define the point of the photon emission as the coordinate origin, and denote the position vector in spatial coordinate as r _* and the distance from the origin as r_*(=∣ r _*∣). In this case, the PDF $P({x}_{* }^{\mu })$ at a time t_* for the scattered photons exhibits a spatially spherical symmetric distribution. Consequently, the PDF is a function of t_* and r_*, i.e., P(t_*, r_*) and P_n (t_*, r_*). The equation to be solved to derive the analytical formula for P_n (t_*, r_*) is obtained by a method similar to that used in previous studies of random flight (e.g., Rayleigh 1919; Hughes 1995).

If we set the initial position of a photon at the origin of the coordinates, the probability ${p}_{0}({x}_{* }^{\mu })$ is given by the delta function as ${p}_{0}({x}_{* }^{\mu })$ $=\,\delta ({x}_{* }^{\mu })$, and its Fourier transform is ${\hat{p}}_{0}({k}_{* })$ = 1. We denote the coordinate ${x}_{* }^{\mu }$ describing P_n as x_*n . The relation between P_n (x_*n ) and P_n−1(x_*n−1) is described by the convolution given as P_n (x_*n ) = ∫p(x_*n − x_*n−1)P_n−1(x_*n−1)d⁴ x_*n−1, where p(x_*n − x_*n−1) is a probability that a photon is scattered from x_*n−1 to x_*n . The PDF of an isotropic scattering p(t_*, r _*) is given by $p({t}_{* },{{\boldsymbol{r}}}_{* })=p({t}_{* },{r}_{* })=\delta ({t}_{* }-{r}_{* })\theta ({t}_{* }){e}^{-{r}_{* }}/(4\pi {r}_{* }^{2})$, and its Fourier transform is $\hat{p}({\omega }_{* },{k}_{* })\,=\arctan [{k}_{* }/(1-i{\omega }_{* })]/{k}_{* }$. Since the convolution theorem tells ${\hat{P}}_{n}({k}_{* })=\hat{p}({k}_{* }){\hat{P}}_{n-1}({k}_{* })$ (e.g., Folland 1992), we obtain ${\hat{P}}_{n}{({k}_{* })=\hat{p}{({k}_{* })}^{n+1}{\hat{p}}_{0}({k}_{* })=\hat{p}({k}_{* })}^{n+1}$. Then, the PDF P_n (t_*, r_*) is calculated by the inverse Fourier transform given by ${P}_{n}{({t}_{* },{r}_{* })=(2\pi )}^{-4}\int {e}^{{{ik}}_{* \mu }{x}_{* }^{\mu }}{\{\hat{p}({\omega }_{* },{k}_{* })\}}^{n+1}{d}^{4}{k}_{* }$. Finally, we can obtain the variable separation form of P_n (t_*, r_*) given by ${P}_{n}({t}_{* },{r}_{* })={t}_{* }^{n-3}{e}^{-{t}_{* }}{V}_{n}({r}_{* }/{t}_{* })\theta ({t}_{* })/(4\pi )$ where θ(t_*) is a step function and V_n (r_*/t_*) is an even function of r_*/t_*. Since the function P_n (t_*, r_*) contains r_* in the form r_*/t_*, the function P_n (t_*, r_*) is in variable separation form of variables t_* and v ≡ r_*/t_*. It should be noted that the function V_n (v) defined in the range 0 ≤ v ≤ 1, i.e., 0 ≤ r_* ≤ t_*, ensures that the formula for the PDF preserves causality. The derivation of the analytical expression for the function V_n (v) is explained below.

For n = 0, P₀(t_*, r_*) = p(t_*, r_*) as denoted above. We can find the function V_n (v) obeys the equation vV_n (v) = S_n (v) for n = 1 and the differential equations dⁿ⁻¹{vV_n (v)}/dvⁿ⁻¹ = S_n (v) for n ≥ 2 where ${S}_{n}{(v)=({i}^{n-2}/\pi )\{{f}_{n}(-{iv})-(-1)}^{n+1}{f}_{n}({iv})\}$. Here, ${f}_{n}{(x)\equiv (\pi /2-{\tan }^{-1}x)}^{n+1}$. The solution of V_n (v)(n ≥ 2) is described as V₂(v) = W₂(v)/v and ${V}_{n}(v)\,={\sum }_{k=0}^{{k}_{\max }}{C}_{n}^{k}{v}^{2k}+{W}_{n}(v)/v$ (for n ≥ 3) where C_n ^k are constants, W_n (v) is a function satisfying W_n (0) = 0, and ${k}_{\max }=n/2-2$ ((n − 1)/2 − 1) for even (odd) integer n. Using the properties of the function P_n (t, r), the coefficient C_n ^k is computed as

$$\begin{eqnarray}&&{C}_{n}^{k}=\displaystyle \frac{2{\left(-1\right)}^{k}}{\pi (2k+1)!(n-3-2k)!}{\int }_{0}^{\infty }\displaystyle \frac{{\left({\tan }^{-1}X\right)}^{n+1}}{{X}^{n-1-2k}}{dX},\end{eqnarray} \tag{ 1 }$$

which can be analytically solved by integration by parts. With a new variable $y=2{\tanh }^{-1}v$, the function W_n (v) can be solved by n − 1 times integration by parts and we obtain the expression ${W}_{n}(v)={\sum }_{k\,=\,0}^{n+1}{c}_{n+1}^{k}{I}_{n-1}^{k}(y(v))$ where the coefficients c_n ^k is given as ${c}_{n}^{k}={(i\pi )}^{n-k-1}{2}^{-n}{}_{n}{C}_{k}\{{(-1)}^{n}-{(-1)}^{k}\}$ and ${I}_{n}^{k}(y)$ is solved by the recurrence relation given as

$$\begin{eqnarray}\begin{array}{rcl}{I}_{n}^{k}(y) & = & \left[{I}_{n-1}^{k}(\bar{y}){s}_{0}^{1}(\bar{y})h(\bar{y})+\displaystyle \sum _{i=1}^{n-2}{\left(-1\right)}^{i}{s}_{1}^{i}(\bar{y}){I}_{n-1-i}^{k}(\bar{y})\right.\\ & & {\left.+\,\displaystyle \sum _{i=0}^{k}{\left(-1\right)}^{i+n+1}{s}_{i+1}^{n-1}(\bar{y})\displaystyle \frac{{d}^{i}}{d{\bar{y}}^{i}}{\bar{y}}^{k}\right]}_{\bar{y}=0}^{\bar{y}=y},\end{array}\end{eqnarray} \tag{ 2 }$$

where $h(y)=1+\cosh y$ and the function ${s}_{n}^{m}(y)$ is defined by ${s}_{0}^{1}(y)=\tanh (y/2)/h(y)$, ${s}_{n}^{m}(y)={\int }_{0}^{y}{s}_{n-1}^{m}({y}_{1}){{dy}}_{1}$, and ${s}_{0}^{m}(y)\,={s}_{1}^{m-1}(y)/h(y)$. The coefficients C_n ^k and the functions ${s}_{n}^{m}(y)$ can be solved analytically by Mathematica to obtain an analytical expression for the function V_n (v). We obtain the analytic results of V_n (v) as, e.g., V₁(v) = (2/v)t(v), ${V}_{2}(v)={-({\pi }^{2}/4)f(v)+3f(v)t(v)}^{2}-(3/v){f}_{2}(v)+(6/v){f}_{+}(v)t(v),$ and

$$\begin{eqnarray*}\begin{array}{rcl}{V}_{3}(v) & = & -\displaystyle \frac{{\pi }^{2}}{2}f(v)-\displaystyle \frac{{\pi }^{2}({v}^{2}+1)}{2v}t(v)+6f(v)t{\left(v\right)}^{2}\\ & & +2{vf}{\left(v\right)}^{2}t{\left(v\right)}^{3}+\displaystyle \frac{{\pi }^{2}}{2}\{{f}_{+}(v)+{f}_{-}(v)\}\\ & & +\,\displaystyle \frac{12}{v}{f}_{+}(v)t(v)-12{f}_{+}(v)t{\left(v\right)}^{2}+6{f}_{3}(v)\\ & & +\,6(2t(v)-{v}^{-1}){f}_{2}(v),\end{array}\end{eqnarray*}$$

where $t(v)={\tanh }^{-1}v$, f(v) = 1/v − 1, ${f}_{\pm }(v)=\mathrm{log}2/(1\pm v)$, f_n (v) = Li[(v − 1)/(v + 1)] where Li_n is the polylogarithm of order n. Similarly, we obtained analytical expressions for V_n (v) (n = 4, 5, ⋯ , 10), but as n grows, the expressions become longer and it is not practical to calculate P(t_*, r_*) using these closed forms of the analytic expression of V_n (v). Instead, when we compute P(t_*, r_*), we use the series expansion of V_n (v) for large n given as ${V}_{n}(v)={\sum }_{k=0}^{\infty }{C}_{n}^{k}{v}^{2k}$ where the coefficients of the series expansion are calculated by Equation (1) and the expansion coefficient of S_n (v), which can be analytically obtained.

When we calculate P(t_*, r_*) from the sum of P_n (t_*, r_*), we should set the range of n. If t is large (t_* ≳ 10²), since the major portion of P_n resides in the range of $n-\sqrt{n}\lt {t}_{* }\lt n+\sqrt{n}$, we set the range of n as t_* − c_* $\sqrt{{t}_{* }}\lt n\lt {t}_{* }\,+\,$ c_* $\sqrt{{t}_{* }}$ where the value of c_* is set to achieve sufficient numerical precision (c_*≳3). When t_* is much larger (t_* > 10³), P can be treated as a perturbation of the diffusion approximation P_diff (defined below). In this case, the expansion coefficients of P are obtained by extrapolation with the precomputed P at some time t_* and P_diff. In this study, we use 20th-order interpolation for this calculation. On the other hand, when t_* is not large (t_* ≲ 10²), the sum of P_n can be calculated directly.

Figure 1 shows the PDF of scattering photons in a static medium from the time t_* = 10⁻³ to 10⁷. For t_* ≲ 10, the spikelike peak can be seen at r_* = t_*, which corresponds to the peaks held by P₁(t_*, r_*) and P₂(t_*, r_*). In other words, this peak represents the traces of photons propagating at the speed of light. We call this period the isotropic free-propagation state (F). On the other hand, in case t_* ≳ 10, the maximum of the distribution lies at the center of r_* = 0. As time passes (t_* ≫ 1), the PDF P(t_*, r_*) of the scattering photons approaches the distribution P_diff(t_*, r_*) of the diffusion approximation, i.e., ${P}_{\mathrm{diff}}({t}_{* },{r}_{* })={\{3/(4\pi {t}_{* })\}}^{3/2}\exp \{-3{r}_{* }^{2}/(4{t}_{* })\}$. Thus, we call this period the diffusion state (D). The period around t_* ∼ 1–10 is in the intermediate state (I) between the two states (F and D). That is, it has a local maximum of the PDF at the center of r_* = 0 and also a spikelike peak at r_* = t_*. Results of Monte Carlo (MC) simulations are shown as dots in Figure 1, confirming that our analytic formulas reproduce the simulation results. In the calculations in Figure 1, using the analytical solution can accelerate the computation by a factor of 10² ∼ 10⁵ compared to MC simulations. In addition, in the case of MC, it is difficult to accurately obtain values near the maximum of the PDF. Therefore, the analytical solution is superior to the MC calculation in terms of computation time, computational accuracy, and computational feasibility.

Zoom In Zoom Out Reset image size

We can also calculate the photon number density flux N^μ in the rest frame analytically. The time component of this flux is ${N}^{0}({x}_{* }^{\mu })={N}_{\mathrm{all}}P({x}_{* }^{\mu })$ and the flux satisfies the particle number conservation law, ∇_μ N^μ = 0. Then, the spatial component of the flux Nⁱ (i = 1, 2, 3) is calculated as Nⁱ = N_all Q(t_*, r_*) where Q(t_*, r_*) is defined as

$$\begin{eqnarray}&&Q({t}_{* },{r}_{* })=\displaystyle \frac{{r}_{* }}{{t}_{* }}\left\{P({t}_{* },{r}_{* })+\displaystyle \sum _{n=0}^{\infty }{c}_{n}({t}_{* }){X}_{n}\left(\displaystyle \frac{\sqrt{n}{r}_{* }}{{t}_{* }}\right)\right\}.\end{eqnarray} \tag{ 3 }$$

Here, ${c}_{n}({t}_{* })={n}^{3/2}{t}_{* }^{n-3}{e}^{-{t}_{* }}({t}_{* }-n)/n!$ and ${X}_{n}(u)\equiv {u}^{-3}{\int }_{0}^{u}{\tilde{u}}^{2}{U}_{n}(\tilde{u})d\tilde{u}$ where ${U}_{n}(u)\equiv {V}_{n}(u/\sqrt{n})n!/(4\pi \sqrt{{n}^{3}})$. The function X_n (u) can be calculated analytically as the closed analytic form or as the series expansion form as ${X}_{n}(u)={\sum }_{k=0}^{\infty }{C}_{k}{u}^{2k}/(2k+3)$ when U_n (u) is expanded as ${U}_{n}(u)={\sum }_{k=0}^{\infty }{C}_{k}{u}^{2k}$.

The photon number density flux in the laboratory frame N^μ is calculated as ${N}^{\mu }={{\rm{\Lambda }}}_{\nu }^{\mu }N{{\prime} }^{\nu }$ where ${{\rm{\Lambda }}}_{\nu }^{\mu }$ is a matrix representing the Lorentz transformation and $N{{\prime} }^{\nu }$ is the photon number density flux in the rest frame. The PDF in the laboratory frame is calculated as $P({x}_{* }^{\mu })={N}^{0}/N$ where N⁰ is the time component of the photon number density flux and N is the number density at a constant time surface in the laboratory frame, which is calculated by the integration of N⁰ in three-dimensional space in the laboratory frame.

Figure 2 shows the cross section in the x_*–z_* plane of the PDF calculated by the analytic formula in the laboratory frame. Figure 3 shows the time evolution of the mean of r_*, ${\overline{r}}_{* }$, in the laboratory frame defined by the following three-dimensional integration in space, ${\overline{r}}_{* }({t}_{* })={\int }_{3{\rm{d}}\,\mathrm{space}}{r}_{* }\,P({t}_{* },{{\boldsymbol{r}}}_{* }){d}^{3}{{\boldsymbol{r}}}_{* }$, where P(t_*, r _*) is the PDF at the constant t_* surface in the laboratory frame. Figure 3 shows the results of both the analytical formula (solid line) and the MC simulation (black circle).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 2, the case of Γβ = 0 (the first column) reflects the distribution described in Figure 1. For Γβ = 0, approximately, ${\overline{r}}_{* }={t}_{* }$ (i.e., light-speed propagation) for t_* ≲ t_a , and ${\overline{r}}_{* }=4\sqrt{{t}_{* }/(3\pi )}$ (diffusive expansion) for t_* ≳ t_a , where the time t_a is the boundary between the free-propagating and diffuse states, calculated as t_a = 16/(3π). In Figure 3, we can confirm that the analytical formulas and the MC simulation results are consistent for all periods.

As shown in Figure 2, for the other cases, Γβ = 0.1 (the second column) and 1.0 (the third column), the distribution of the PDF is biased in the +z_*-direction due to the boost effect at all the times. For t_* ≲ 1, we see that most of the photons are concentrated near z_* = t_*. That is, for Γβ > 0 and t_* ≲ 1, they propagate at the speed of light anisotropically in the boost direction, in contrast to the isotropic propagation at the speed of light for Γβ = 0. When t_* ≲ 1, we can see the anisotropic ballistic motion (denoted as (B) in Figure 2) whose mean velocity is nearly equal to the speed of light. Figure 3 also confirms that the propagation is at the speed of light when t_* ≲ 1. For t_b ≲ t_*, the scattered photons are in a state of dynamic diffusion (denoted as (D²) in Figure 2), where the time t_b is the boundary between the diffusion and the dynamic diffusion and is calculated as ${t}_{b}={(4/(3\pi {\beta }^{2})){[1+(1+2{\beta }^{2}/{\rm{\Gamma }})}^{1/2}]}^{2}$. That is, as can be seen in Figure 2, they expand diffusively while being boosted in the z_*-direction. In this case, the center of the distribution moves approximately according to ${\overline{r}}_{* }=\beta {t}_{* }-8\beta /(3\pi {\rm{\Gamma }})$. This can be seen in both Figures 2 and 3. In Figure 2, when Γβ ≳ 0.1, we can see the intermediate state (I) between the two states (B and D²) around t_* ∼ 1.

In this study, under the assumption of isotropic and elastic scattering, we analytically describe the PDF $P({x}_{* }^{\mu })$ in a time slice surface and the PDF ${P}_{n}({x}_{* }^{\mu })$ in space-time that describes the collective behavior of scattering photons in a relativistic fluid. The derived equations reproduce the results of the MC simulations and succeeds in smoothly describing the intermediate state between the free-streaming (F) and the diffusion (D) state when Γ ∼ 0 or between the anisotropic ballistic motion (B) and the dynamic diffusion (D²). Although we considered a point source represented by the delta function in this study, we believe that the distribution of photons emitted from spatially spread out or temporally continuous radiation sources can be represented by the superposition of the present analytical solution. In the future, the authors intend to implement this study in the general relativistic radiative transfer simulations and extend it to more general case of photon diffusion. Based on the present study, some of the authors of this Letter are working on calculations of the distribution function ${ \mathcal F }({x}_{* }^{\mu },{k}_{* }^{\mu })$ in phase space for the scattering photons.

We are deeply grateful to the anonymous referee for the insightful and constructive comments, which have significantly improved the quality of the manuscript. The authors are grateful to H. Itoh, H. Takahashi, and T. Kawashima for valuable discussions. This work was supported by JSPS KAKENHI grant Nos. 15H03638 (M.U.), 16K05302 (R.T. and M.U.), 22J10256 (M.M.T.), 18K03710, 21H04488 (K.O.), 19H00697 (M.U. and R.T.), 21H01132 (R.T., M.U., K.O., and Y.A.), 18K13581, 23K03445 (Y.A.), 22KJ0382 (M.T.), and 22K03686 (N.K.). This research was also supported in part by Multidisciplinary Cooperative Research Program in CCS, University of Tsukuba. This work was also supported by MEXT as Program for Promoting Researches on the Supercomputer Fugaku (Structure and Evolution of the Universe Unraveled by Fusion of Simulation and AI, JPMXP1020240219; K.O., black hole accretion disks and quasi-periodic oscillations revealed by general relativistic hydrodynamics simulations and general relativistic radiation transfer calculations, JPMXP1020240054; K.O.), and by Joint Institute for Computational Fundamental Science (JICFuS, K.O.).