A Relativistic Formula for the Multiple Scattering of Photons

We have discovered analytical expressions for the probability density function (PDF) of photons that are multiply scattered in relativistic flows, under the assumption of isotropic and inelastic scattering. These expressions characterize the collective dynamics of these photons, ranging from free-streaming to diffusion regions. The PDF, defined within the light cone to ensure the preservation of causality, is expressed in a three-dimensional space at a constant time surface. This expression is achieved by summing the PDFs of photons that have been scattered n times within four-dimensional space-time. We have confirmed that this formulation accurately reproduces the results of relativistic Monte Carlo simulations. We found that the PDF in three-dimensional space at a constant time surface can be represented in a separable variable form. We demonstrate the behavior of the PDF in the laboratory frame across a wide range of Lorentz factors for the relativistic flow. When the Lorentz factor of the fluid is low, the behavior of scattered photons evolves sequentially from free propagation to diffusion, and then to dynamic diffusion, where the mean effective velocity of the photons equates to that of the fluid. On the other hand, when the Lorentz factor is large, the behavior evolves from anisotropic ballistic motion, characterized by a mean effective velocity approaching the speed of light, to dynamic diffusion.

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 pro-posed 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 PDF provided in this paper, describing relativistic diffusion for photons, likely represents the first instance of an analytical solution for the relativistic diffusion problem involving these particles.
In this paper, 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 ( §2) and in the laboratory frame ( §3).In the final section ( §4), the conclutions are given.

DISTRIBUTION OF SCATTERING PHOTONS IN A STATIC MEDIUM
According to relativistic kinetic theory (e.g., Lindquist 1966;Ehlers 1971;Israel 1972;Sachs & Ehlers 1971), the collective behavior of particles is described by an invariant distribution function F(x µ , p µ ) where x µ and p µ are, respectively, the coordinate and the momentum of a particle.The particle number density flux N µ is given by N µ = dN µ = F(x µ , p µ )p µ dP and the particle number N (x µ )dV in a volume element dV at a time slice surface (where t = t 0 ) is given by the projection of N µ onto the vector volume element, i.e., N (x µ )dV = dV F(x µ , p µ )(−û µ p µ )dP where ûµ is a timelike unit vector.In this study, the probability density function (PDF) in a time slice surface are used to describe the collective behavior of scattered photons.The PDF P (x µ )| t=t0 giving the probability that particles exist in some region at the time slice surface by volume integration is given by P (x µ )| t=t0 = N (x µ )/N all where N all is the number of photons in three space V at t = t 0 calculated as N all = Vt=t 0 N (x µ )dV .Then, the PDF The time component of N µ (x µ ) as N 0 (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 ℓ 0 in the fluid rest frame (Rybicki & Lightman 1991;Mihalas & Mihalas 2013).Here, we introduce x µ * = (t * , r * ) ≡ x µ /ℓ 0 = (ct/ℓ 0 , r/ℓ 0 ) and k µ * = (ω * , k * ) ≡ k µ ℓ 0 = (ℏωℓ 0 /c, ℏkℓ 0 ).We also define r * ≡ |r * | and k * ≡ |k * |.
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 * ) = ∞ n=0 P n (t * , r * ) where P n (t * , r * ) is a PDF in four-dimensional spacetime that represents the distribution of the next scattering point of photons scattered n times, which satisfy the normalization P n (t * , r * )d 4 x * = 1 where d 4 x * is an invariant volume element of four-dimensional spacetime. 1In 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, 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 µ * ) 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 n=0 Pn(t * , r * ) corresponds to the series expansion of the function P (t * , r * ) in terms of scattering number n.Specifically, the analytic solution derived below satisfies the normalization Pn(t * , r * ) 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 µ * ) is given by the delta function as p 0 (x µ * ) = δ(x µ * ), and its Fourier transform is p0 (k * ) = 1.We denote the coordinate x µ * 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 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 0 (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 for even (odd) integer n.Using the properties of the function P n (t, r), the coefficient C k n is computed as 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 solved by the recurrence relation given as where h(y) = 1+cosh y and the function s m n (y) is defined by s 1 0 (y) = tanh(y/2)/h(y), s m n (y) = y 0 s m n−1 (y 1 )dy 1 and s m 0 (y) = s m−1 1 (y)/h(y).The coefficients C k n and the functions s m n (y) can be solved analytically by Mathematica to obtain an analytical expression for the function V n (v).We obtain the analytic results of , and where ] 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 form 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 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 2 ), since the major portion of P n resides in the range of n − √ n < t * < n + √ n, we set the range of n as where the value of c * is set to achieve sufficient numerical precision (c * ≳ 3).When t * is much larger (t * > 10 3 ), 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 pre-computed 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 2 ), 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 −3 to 10 7 .For t * ≲ 10, the spike-like peak can be seen at r * = t * , which corresponds to the peaks held by P 1 (t * , r * ) and P 2 (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 diff (t * , r * ) = {3/(4πt * )} 3/2 exp{−3r 2 * /(4t * )}.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 spike-like 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 2 ∼ 10 5 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.

DISTRIBUTION OF SCATTERING PHOTONS IN A RELATIVISTIC FLOW
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 µ * ) = N all P (x µ * ) and the flux satisfies the particle number conservation law, ∇ µ N µ =0.Then, the spatial component of the flux ).The function X n (u) can be calculated analytically as the closed analytic form or as the series expansion form as The photon number density flux in the laboratory frame N µ is calculated as N µ = Λ µ ν N ′ ν where Λ µ ν is a matrix representing the Lorentz transformation and N ′ ν is the photon number density flux in the rest frame.The PDF in the laboratory frame is calculated as P (x µ * ) = N 0 /N where N 0 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 0 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 * , r * , in the laboratory frame defined by the following three-dimensional integration in space, r * (t * ) = 3d space r * P (t * , r * )d 3 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).
In Figure 2, the case of Γβ = 0 (the first column) reflects the distribution described in Figure 1.For Γβ = 0, approximately, r * = t * (i.e., light-speed propagation) for t * ≲ t a , and r * = 4 t * /(3π) (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 Monte Carlo 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 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 2 ) in Figure 2], where the time t b is the boundary between the diffusion and the dynamic diffusion and is calculated as 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 r * = βt * − 8β/(3πΓ).This can be seen in both Figure 2 and Figure 3.In Figure 2, when Γβ ≳ 0.1, we can see the intermediate state (I) between the two states (B and D 2 ) around t * ∼ 1.

CONCLUSIONS
In this study, under the assumption of isotropic and elastic scattering, we analytically describe the PDF P (x µ * ) in a time slice surface and the PDF P n (x µ * ) in spacetime that describe 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 2 ).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 paper are working on calculations of the distribution function F(x µ * , k µ * ) in phase space for the scattering photons.
We are deeply grateful to the anonymous referee for their insightful and constructive comments, which have significantly improved the quality of the manuscript.

Figure 2 .
Figure 2. The probability density function P (t * , r * ) of scattering photons in the laboratory frame of a relativistic flow with Γβ = 0.0 (the first column), 0.1 (the second column) and 1.0 (the third column) at the time t * = 0.001, 1.0, 10.0, 100.0 and 1000.0 (from left to right).In each plot, the abscissa and the ordinate represent x * -axis and z * -axis, respectively, and the outermost circle (dotted circle) shows the location of r * = t * .The symbols (F, B, I, D, D 2 ) at the bottom-right corner of each plot represents the state of the scattering photons (see, texts).The maximum and minimum values, Pmax and Pmin, of the color bar are set to show the distribution of the probability density function.The circle with the radius of √ t * and the center at (x * , z * ) = (0, βt * ) are shown at the time t * = 10.0, 100.0 and 1000.0 (white dotted circle).

Figure 3 .
Figure 3.The time evolution of the mean radial distance r * of the probability distribution of the scattering photons in the laboratory frame of the fluid with Γβ = 0.0, 10 −4 , 10 −3 , 10 −2 , 0.1, 1.0 and 10 (from top to bottom).The results of the formula (solid line) and the MC simulations (blank circle) are shown.The formula representing the state of the scattering photons are shown: r * ∼ t * (isotropic free-streaming or anisotropic ballistic motion), r * ∝ √ t * (diffusion, i.e., diffusive expansion) and r * ∼ βt * +const.(dynamic diffusion, i.e., boosted diffusion).