An accurate reaction-diffusion limit to the spherical-symmetric Boltzmann equation

We resolve a long standing question regarding the suitable effective diffusion coefficient of the spherically-symmetric transport equation, which is valid at long times. To that end, we generalize a transport solution in three dimensions for homogeneous media, to include general collisional properties, including birth-death events and linearly anisotropic scattering. This is done by introducing an exact scaling law relating the Green function of the pure-scattering case with the general collision case, which is verified using deterministic and Monte-Carlo simulations. Importantly, the effective diffusion coefficient is identified by inspecting the transport solution at long times.

Introduction.The study of three-dimensional (3D) time-dependent transport phenomena is at the heart of many fields in physics, chemistry and biology.For example in astrophysics, when a massive star explodes, there is a blast of energy that is released in a relatively focused spot [1,2].Similarly, in inertial confinement fusion a large focused release of energy occurs as a result of fusion reactions [3][4][5].In these cases, X-ray photons are propagating fast toward the medium due to the radiative transfer mechanism.Correctly modeling the propagation of these heat waves is a prerequisite for analyzing the outgoing signals, e.g., from a supernovae explosion.Timedependent transport phenomena also appear in electron transport in hot plasmas [6], modeling coda waves of local earthquakes [7], the distribution of neutrons inside a nuclear reactor [8][9][10], and optical properties of biological tissues in biology [11][12][13][14].Notably, in the latter, the balance between scattering and absorption inside the medium are of great interest [15].Even in modeling the transport of bacteria that diffuse in a medium and may reproduce or die, a transport equation can be defined to model the population dynamics [16].
The exact particle propagation is modeled via the Boltzmann (transport) equation for the probability density function (PDF) P ( r, t, Ω) per unit volume d 3 r and direction d Ω, at time t.Notably, in most of the fields specified above, the governing equation of transport is of non-interacting particles, which are aptly described by the linear Boltzmann equation.When the medium is highly scattering, the exact solution of this equation includes a diffusive region in the bulk, and ballistic 'tails' due to particles that almost do not undergo collisions [17][18][19].However, obtaining the solution of the transport equation in the general case is not amenable analytically, and numerically it is highly time consuming.Notably, in a one-dimensional (1D) geometry, where the particles are free to travel only forward and backward, the linear transport equation reproduces the Telegrapher's equation, which can be solved analytically [17,18,[20][21][22].
It is well known that at sufficiently long times, the exact transport solution tends to a reaction-diffusion solution via the central limit theorem, due to multiple scattering events.The central part of the distribution func-tion has a Gaussian shape with increasing width, which corresponds to an effective diffusion coefficient.Naturally, applying a diffusion approximation to the transport equation gives rise to non-physical tails of particles at r > vt [23].Yet, even in the bulk region where the diffusion approximation holds, it remains an open question of which diffusion coefficient to use to accurately reproduce the bulk solution of the full transport equation.
Here, we present a simple yet rigorous derivation of the accurate diffusion coefficient, valid in the spherical symmetric case, in the presence of an arbitrary set of processes, including birth (or branching) and death (or absorption) events, that may occur upon an interaction between an agent and the medium.This diffusion limit is crucial, e.g., in modeling radiative transfer in the non-LTE (Local Thermodynamic Equilibrium) regime, when photon emission dominates absorption [2], which occurs in various astrophysical scenarios, such as shock breakout phase during supernovae, corona of accretion disks and nebular phase of stellar explosions, to name a few.
Our analysis is based on the generalization of the solution of Paasschens [18], which was obtained for homogeneous infinite medium and pure isotropic scattering.While in [35] the solution was extended to the case of linear anisotropic scattering, we here generalize the solution to hold for an arbitrary collision scenario.We first derive an exact scaling relation (in space and time) of the 3D transport solution between the pure-scattering and general-collision case, which holds for the general anisotropic case as well, and find the Green function so-lution for a point source term.The scaling relation is verified via two different numerical schemes: (i) a solution of the deterministic time-dependent equation with the discrete ordinates method (S N method) [36], and (ii) stochastic Monte-Carlo (MC) simulations mimicking the different probabilistic events, see supplemental material (SM), Sec.I.This scaling relation allows us to extend the analytical solution of Paasschens for the general-collision case, from which we infer the diffusion coefficient by computing its long time asymptotics.
Spherical Symmetric Transport Solution.In the case of no external forces and assuming mono-energetic particles, the equation that governs the physical process is the linear Boltzmann transport equation.For a sphericalsymmetric setting, the transport equation takes the following form [8,18,25,36] (see SM, Sec.II): Here the PDF, P (r, t, µ), is a function of the radius r, time t, and the direction of the particle with respect to r, µ ≡ Ω • r, where Ω is the direction of the particle before the collision.In addition, we have assumed that the interaction between particles is negligible, and the medium is locally in equilibrium, such that the Boltzmann collision term can be modeled using the definitions of the mfps for different interactions (or alternatively cross sections).These are the mfp of absorption, ℓ a , and of scattering, ℓ s , where we denote by ℓ t the total mfp due to all physical events.In general, the scattering event is a function of the angular deflection during scattering, i.e., the angle between the direction of the particle before the scattering Ω and after Ω′ [37].As a result, the macroscopic cross section for scattering has the form Σ , where f (µ 0 ) is a normalized distribution of the deflection angle cosine µ 0 ≡ Ω • Ω′ .The function f is problem dependent, and is commonly taken in the Henyey-Greenstein form [38] (see SM, Sec.IVa) [39].Therefore, the collision term in Eq. ( 1) can be written explicitly as a "gain-loss" term using the above mfps definitions.
Additional notations in Eq. ( 1) include the particle velocity v and an external source Q ext (r, t, µ), which may be a function of space, direction and time.Finally, c(r, t) represents the mean number of particles that are emitted in an interaction, including sources [8-10, 23-25, 40-43], for example due to cell division, and satisfies [44]: Here P(r, t) = 1  2)], see legend, in the isotropic case.Our analytical solution, Eq. (5) (dashed lines) is compared with the numerical solution using the SN (dash-dotted lines) and MC (solid lines) methods.The noisy solution at r ≃ 0 is due to the small volumes of the numeric cells near the origin.r and time.In the case of no internal source and homogeneous media, 0 c 1 becomes the ratio of the total and scattering mfps, . If the medium involves internal sources such as birth events proportional to P( r, t), one has c = (ℓ −1 Here n is the mean number of particles created in a birth event and ℓ b is the birth mfp.Clearly, in systems with c < 1 (c > 1) the expected number of particles decays (exponentially grows) in time.For example, in astrophysics, the emitted black-body energy density B and radiation energy density E allow us to define c as While in LTE, E ≈ B and thus ω eff ≈ 1, in non-LTE, when the matter is hotter than the radiation, B > E and ω eff > 1, and vice versa [2].Naturally, when c > 1, the exponential growth is eventually arrested by nonlinear effects.
Henceforth, we will refer to the cases without and with internal sources as the sub-scattering and superscattering problem, respectively.We will also focus on the case of infinite homogeneous medium; i.e., we assume c(r, t) = c is a constant.Importantly, we perform our calculations below using an external point-like source at the origin: Q ext (r, t) = δ(r)δ(t)/4πr 2 , where δ(x) is the Dirac delta function.In this case, the calculated population density will serve as the Green function, which will allow the subsequent calculation of the population density for any space-and time-dependent external source term by the Green's convolution integral.
In Fig. 1 we compare our analytical solution to Eq. ( 1), see Eq. (5) below, to numerical solutions using the S N and MC methods.Here we show P(r, t) as a function of the rescaled position r = r/ℓ t , at different rescaled times t = tv/ℓ t , for various values of c, 0.5 c 2. Both the analytical and numerical solutions show ballistic behavior near r = vt.Yet, the ballistic region shrinks and the diffusive behavior becomes dominant as time advances.
We now solve Eq. ( 1) is the general collision scenario, i.e., c = 1.To do so, we identify a scaling relation between the PDF, P (r, t, µ; c), in the general c = 1 case, and the PDF, P (r, t, µ; 1), in the pure-scattering, c = 1, case.Previously, it was shown that for c 1, the scaling rela- tion reads [18,30,35]: P (r, t, µ; c) = e −vt/ℓa P (r, t, µ; 1), with c = ℓ −1 s /(ℓ −1 a +ℓ −1 s ).We propose the following exact scaling relation for the Boltzmann equation in spherical symmetry (1), for a general value of c (see SM, Sec.III): where ℓ t denotes the total mfp.In Fig. 2 we numerically verify the scaling relation by plotting P(r, t) using the S N and MC methods, for various values of c, both in the ballistic and diffusion regions.We now use this scaling relation to find the solution for the population density P(r, t; 1), for c = 1.Here, our starting point is the result in the pure-scattering case, c = 1, obtained by Paasschens [18].This result assumes a homogeneous infinite medium, and was obtained exactly for 2D and 4D geometries using a Fourier transforms in space and time.Interestingly, this method is less useful in the 3D case as the inverse transform cannot be found analytically.Instead, Paasschens used an interpolation between the 2D and 4D solutions [18], which yields P(r, t; 1) in the 3D case: where Θ(r) is the Heaviside step function and Γ(x) is the gamma function.Notably, see below, at long times Eq. ( 4) converges to the diffusion solution with D = vℓ t /3 and c = 1.Solution (4) can be extended to the linear anisotropic scattering case, upon using ℓ tr instead of ℓ t , where ℓ −1 tr ≡ ℓ −1 t (1 − g), and g = μ0 denotes the average scattering angle cosine [35].Here, ℓ tr is called the transport mfp [25].However, we show that this solution is valid only at long times (see SM Sec.IVa).
With the scaling relation (3) we have found, and using Eq. ( 4) for c = 1, the general solution for the population density, P(r, t; c) for c = 1 reads (see SM, Sec.IVb): This solution coincides with the Paasschens solution for c = 1 and g = 0.It can be shown to be valid for g ≪ 1 at all times; yet, it becomes valid for all values of g, −1 < g < 1, at long times (see SM, Sec.IV).In Fig. 1 our analytical solution [Eq. ( 5)] is shown to excellently agree with the numerical results of both S N and MC, for all values of c, in the isotropic case.In Fig. 3 we test the accuracy of Eq. ( 5) for the general anisotropic case for c = 1.Here, we plot the PDF for the dimensionless spatial position r = 3, as a function of t for isotropic (g = 0), forward (g = 0.3) and backward (g = −0.3)scattering.Importantly, we have checked that our results exactly reproduce earlier benchmark results for these g values [9].While for g = 0 our solution coincides with Eq. ( 4), and is thus accurate at all times, for g = 0 the picture is different.For forward scattering (g = 0.3) Eq. ( 5) yields a good agreement only at late times, whereas for backward scattering (g = −0.3),Eq. ( 5) yields a good agreement also for earlier times, see Fig. 3.This is because forward scattering increases anisotropy, thereby delaying the validity of the diffusion assumption.Conversely, backward scattering causes particles to return to the origin, thereby increasing isotropy.It is also demonstrated that for the various values of g, the diffusion solution yields solutions that propagate faster than the particle speed at early times, and converge to the exact solution at long times.
Determining the Diffusion Coefficient at Long Times.We now identify the effective diffusion coefficient by computing the late-time asymptotics of Eq. (1).At late times, t ≫ v/ℓ t , the PDF is close to being isotropic and we can use the so-called diffusion assumption, P (r, t, µ) ≈ P(r, t) + 3µP 1 (r, t), with P 1 (r, t) ≡ 1 −1 µdµP (r, t, µ) be-ing the flux.Integrating Eq. ( 1) once over dµ and once over µdµ yields two equations for the zeroth and first moments of the PDF.After some algebra we obtain: Note that while Eq.(6a) is exact, Eq. ( 6b) is an approximated Fick's law, using the diffusion assumption, and neglecting the ∂P 1 (r, t)/∂t term compared to (v/ℓ t )P 1 (t, r) [45].Substituting Eq. (6b) in Eq. (6a) yields a of reaction-diffusion equation [46,47], with a diffusion constant of D tr ≡ vℓ tr /3, where ℓ tr was defined above [8,12,13,24,25].The solution of Eqs. ( 6) yields: Yet, this solution yields an incorrect scaling relation P diff (r, t; c) = P diff (r, t; 1)e −(1−c)vt/ℓt and thus, violates the scaling in Eq. ( 3) for c = 1 [12,13].This is also evident by comparing Eq. ( 7) to numerical results, revealing that the approximation given by the first two moment leading to Eqs. ( 6) breaks down at c = 1 (see Fig. 4).
To find the correct diffusion coefficient for any c = 1, we use Eq. ( 5) in the limit of t → ∞, which reads: This solution coincides with the diffusion solution (7), but with a modified diffusion coefficient: D ′ tr ≡ D tr /c = vℓ tr /(3c), which is exact, even when the classical diffusion assumption is invalid.This result can also be obtained directly by applying the scaling relation to Eq. ( 7) for c = 1: P diff (cr, ct; c) = c 3 e −(1−c)v/ℓtr P diff (cr, ct; 1).Notably, while our result holds for any value of c, in the special case of c 1, the diffusion constant becomes D = vℓ s /3, and provides a rigorous basis for previous semi-empirical studies that obtained a similar result [2,9,12,[26][27][28][29][30].Finally, one can also derive a diffusion equation directly from the transport equation (1) at t → ∞, without explicitly assuming the diffusion assumption, which yields the correct diffusion coefficient in the close vicinity of c ≃ 1 (see SM, Sec.V).
In Fig. 4 (a-c) we compare the probability density P(r, t; c) given by Eq. ( 8) as function of r, with numerical solutions of Eq. (1) at late times, using the S N and MC methods, for different values of g and c.The figure shows that the naïve choice of D = vℓ t /3 yields large errors.On the other hand, using the accurate diffusion constant D = vℓ t /(3c(1 − g)) yields an excellent agreement indicating that the proposed general-collision diffusion coefficient accurately reproduces the solution to the transport equation [Eq.( 5)] at late times, for any c or g.A complementary view to the spatial profiles is to use the value at the center of the Gaussian, P(0, t → ∞) as a function of time [as can be seen from Eqs. ( 7) and ( 8)], which can be numerically found by fitting the S N solution.In Fig. 4(d) we plot P(0, t)e (1−c)vt/ℓt and compare the different diffusion coefficients from Fig. 4(a).Again, using D = vℓ tr /(3c) yields an excellent agreement with the S N results.However, the naïve choice of D = vℓ t /3 yields a highly inaccurate slope.
Discussion.A long-standing question in the field of statistical physics, is the identification of the diffusion coefficient in the diffusion limit of the 3D transport equation, for an arbitrary set of reactions, i.e., c = 1, and weakly-anisotropic scattering.Here we identify the correct diffusion coefficient by deriving a generalized solution for the point source Green function for the 3D spherical-symmetric transport equation.This was done by introducing an exact scaling relation between the general-collision (c = 1) and pure-scattering (c = 1) solutions.We verified our analytical results using the S N method and Monte-Carlo simulations.Unlike previous results [2,9,12,[26][27][28][29][30], the diffusion coefficient we have derived is accurate for any c, i.e., for arbitrary birthdeath reactions.These results may provide important insight into fields such as photon diffusion in tissue optics and radiative transfer in non-equilibrium astrophysics.

2 1− 1 1 FIG. 1 .
FIG.1.The population density P(r, t) at different times for various values of c [Eq. (2)], see legend, in the isotropic case.Our analytical solution, Eq. (5) (dashed lines) is compared with the numerical solution using the SN (dash-dotted lines) and MC (solid lines) methods.The noisy solution at r ≃ 0 is due to the small volumes of the numeric cells near the origin.

2 FIG. 2 .
FIG.2.A numerical verification of the scaling relation [Eq.(3)] for various values of c.The solid lines represent the SN (left) and MC (right) results for the population densities, calculated with a specific value of c, while the dashed lines are the rescaled results with c = 1, using Eq.(3).

3 FIG. 3 .
FIG.3.The population density P(r, t) versus the normalized time vt/ℓ, at r = 3 (right), for various values of anisotropy g.Shown are the results obtained from the extended Paasschens solution (solid lines), SN simulations (dashed lines) and the accurate diffusion[9] (dotted lines).

4 FIG. 4 .
FIG. 4. (a)The population density P(r, t) at long times (vt/ℓ = 100), for c = 0.6 and g = 0.2, versus the normalized space coordinate vt/ℓ.Shown are the full transport solutions obtained from SN solutions (blue curves), MC simulations (red curves) and the full analytical solution (5) (magenta), which are compared with the results of the naïve diffusion theory (black) and the accurate diffusion coefficient (8) (green curves).(b) The population density P(r, t) at long times (vt/ℓ = 100), for different values of c (g = 0.2), exact (SN , MC and analytic) versus the accurate diffusion solution.(c) Same for different values of g (c = 0.6).(d) A log-log plot of the population density at the origin P(0, t), for c = 0.6 and g = 0.2, as a function of the normalized time vt/ℓ.Here, results of the naïve (green) and accurate (red) diffusion solutions are compared with the SN solution (blue).