A Unified Gas-Kinetic Particle Method for Radiation Transport in an Anisotropic Scattering Medium

In this paper, a unified gas kinetic particle (UGKP) method is developed for radiative transfer in both absorbing and anisotropic scattering media. This numerical method is constructed based on our theoretical work on the model reduction for an anisotropic scattering system. The macroscopic solver of this method directly solves the macroscopic anisotropic diffusion equations, eliminating the need to solve higher-order moment equations. The reconstruction of macroscopic scattering source in the microscopic solver, based on the multiscale equivalent phase function we proposed in this work, has also been simplified as one single scattering process, significantly reducing the computational costs. The proposed method has also the property of asymptotic preserving. In the optically thick regime, the proposed method solves the diffusion limit equations for an anisotropic system. In the optically thin regime, the kinetic processes of photon transport are simulated. The consistency and efficiency of the proposed method have been validated by numerical tests in a wide range of flow regimes. The novel equivalent scattering source reconstruction can be used for various transport processes, and the proposed numerical scheme is widely applicable in high-energy density engineering applications.


Introduction
Thermal radiative transfer (TRT) equations, which describe the time evolution of radiative intensity and its interaction with the background material, have wide applications in astrophysics, atmospheric physics, inertial confinement fusion (ICF), high-temperature flow systems, plasma physics, etc. [1][2][3][4][5][6][7][8].TRT equations are intrinsically nonlinear due to the absorption-emission process, which renders the system difficult to solve [9][10][11].Furthermore, the intrinsic high dimensionality of TRT equations significantly increases the computational costs.Developing numerical methods with high accuracy and high efficiency has become an important topic over the past few decades.
Currently, three kinds of models are used to describe the radiation transport process: the kinetic model, the macroscopic moment model, and the machine learning model.The kinetic model, also known as the Boltzmann transport equation, comprehensively describes the photon transport, absorbing, emission, and scattering processes at a mesoscale level.It is the most fundamental equation of photon transport, with seven dimensions, i.e., the physical space, the velocity space, the frequency space, and time.However, it is computationally expensive to solve directly.The macroscopic moment model, in contrast, offers a means to reduce the high-dimensional Boltzmann transport equation to a low-dimensional equation based on the spherical harmonic function expansion [12].Nonetheless, further research is necessary to bolster the reliability of its closed function.The machine learning model, as a relatively new model, utilizes deep neural networks to develop efficient macroscopic models [13].However, the generalization ability and confidence of the macroscopic models generated by deep neural networks require additional validation.
Generally, the numerical methods for kinetic equations are categorized into the deterministic method and the stochastic method.The deterministic methods include the macroscopic moments method [14,15], and microscopic discrete ordinate (SN) method [16][17][18][19][20][21].The moment method expands the radiant intensity in a spherical harmonic function space, and the SN method discretizes the velocity space into grid points.For the stochastic method, the radiative transfer equation is solved by simulating the interactions of individual radiation particles with the background material.The implicit Monte Carlo (IMC) method, proposed by Fleck and Cummings [22], is a popular Monte Carlo method for solving the TRT equations.The deterministic method has an explicit formation, which could provide analytic solutions.However, it will suffer from the nonphysical ray effects.Conversely, the stochastic method does not suffer the ray effects and can be applied for high-dimensional phase space directly with good robustness.However, the stochastic process of individual particles will lead to statistical results.In addition, both the deterministic method and stochastic method are single-scale numerical methods, which are unsuitable for problems with complex regimes.
Recently, the multiscaled unified gas kinetic scheme (UGKS) method has been developed to address problems spanning optically thin and thick regimes [23][24][25][26].The UGKS method utilizes a finite volume formulation to solve the macroscopic energy equations, while the discrete ordinates method (DOM) is employed for solving the microscopic transport equation [27,28].The integral solution of the transport equation bridges the macroscopic and microscopic equations, which covers the physics from the free transport to the diffusion limit and makes the UGKS method asymptotic preserving (AP).Based on the framework of UGKS, different numerical schemes have been constructed, including the discrete unified gas kinetic scheme (DUGKS) method [29][30][31], the implicit unified gas kinetic scheme (IUGKS) method [32][33][34][35], the unified gas-kinetic particle (UGKP) method [36,37], and the unified gas-kinetic wave-particle (UGKWP) method [38][39][40].These methods share a similar framework with the UGKS method, making them asymptotic preserving as well.Furthermore, the Monte Carlo solver of the UGKP and UGKWP methods avoids the interdependence of complex unstructured mesh in microscopic transport equation solving, which makes the UGKP and UGKWP methods more suitable for problems with such mesh than the UGKS method.
In our prior work, we extended the UGKP method for the frequency-dependent radiative system, considering both absorption-emission and isotropic scattering processes under a general unstructured mesh [37].However, practical applications usually involve media with anisotropic scattering properties.Up until this point, both the UGKP and UGKWP methods have not considered the effects of the anisotropic scattering process yet.In the present work, we take the next step by extending the UGKP method for the anisotropic scattering system, building upon our previous theoretical foundations [41].The rest of this paper is organized as follows.In Section 2.1, we briefly introduce the UGKP method for the isotropic scattering system (hereinafter referred to as the isotropic-UGKP method).Then, we present the UGKP method for the anisotropic scattering system (hereinafter referred to as the anisotropic-UGKP method) in Section 2.2, based on the foundation laid out in Section 2.1.The asymptotic preserving property of the proposed anisotropic-UGKP method will be discussed in the following Section 2.3.The results of several numerical results will be presented in Section 3 to validate the proposed anisotropic-UGKP method.Finally, we conclude and discuss future work in Section 4.

Numerical Methods
In this section, we present our proposed UGKP method for a gray radiative transfer system with an anisotropic scattering process.The following formulas are presented in a two-dimensional Cartesian space.The angle direction is denoted by is the cosine value of the angle between the propagation direction Ω and the z-axis, and γ ∈ [0, 2π) is the angle between the projection vector of Ω onto the xy-plane and the x-axis.Because of the symmetry of angular distribution in the two-dimensional Cartesian case, we only need to consider ζ > 0.

Isotropic-UGKP Method
In this subsection, we provide a brief introduction to the isotropic-UGKP method for the gray radiative transfer system.

Gray Radiative Transfer System with Isotropic Scattering
The gray radiative transfer equations describe the transport of radiation and its energy exchange with material.Under the assumption of the local thermal equilibrium (LTE), the time-dependent gray radiative transfer equations in the absence of both external and internal sources can be written in the following scaled form: ( where I t, r, Ω is the radiation intensity, which depends on time t, spatial variable r, and angular variable Ω; T(t, r) is the material temperature; φ = acT 4 with the radiation constant a and the speed of light c and ρ = Id Ω; σ a ( r, T) is the absorption coefficient; σ s ( r, T) is the scattering coefficient; ε > 0 is the Knudsen number; L ε a and L ε s are two parameters depending on ε; U m ( r, t) is the material energy density; and C V > 0 is the heat capacity.Taking the angular integration of the radiation transport equation in (1), we obtain the following macroscopic system: where ΩI = ΩId Ω, and β = ∂φ/∂U m = 4acT 3 /C v .

Macroscopic Solver for the Isotropic-UGKP Method
Equation ( 2) can be discretized on both spatial and time variables based on the finite volume method.On the unstructured mesh illustrated in Figure 1, a conservative finite volume numerical scheme for the macroscopic system (2) has the following form: where ∆t = t n+1 − t n is the time interval; ρ n+1 j and φ n+1 j are the cell averaged value at time t n+1 in cell j; and Φ n+1 j,k is the macroscopic flux across the cell edge k.The ρ n+1 j , φ n+1 j , and Φ n+1 j,k are defined as: where V j is the volume of mesh j; p k,m is the center of edge k; l k is the length of edge k; and n k is the unit normal vector of edge k.The construction of Φ n+1 j,k is the key to the isotropic-UGKP method.To calculate the macroscopic flux Φ n+1 j,k in Equation ( 4), one must determine the radiation intensity in the vicinity of the center of edge k by solving the following equation: Equation ( 5) is formatted under the local coordinate system on edge center p k,m , with x and y as the local orthogonal coordinates, and µ and ξ as the local normal and tangential directions.The integral solution of Equation ( 5) is: where λ = c(L ε a σ a + L ε s σ s )/ε.Then, the macroscopic flux is evaluated by The macroscopic quantities φ and ρ in the integral solution (6) can be reconstructed by a piecewise linear polynomial around the edge k as follows: with φ n+1 j,k and ρ n+1 j,k as the value at edge k.The spatial derivatives in (8) are calculated by: where j denotes the neighboring cell which has the common edge k with cell j.The φn+1 k,1 , φn+1 k,2 , ρn+1 k,1 , and ρn+1 k,2 , as the macroscopic quantities on the two vertexes (p k,1 and p k,2 ) of edge k, are calculated as the average macroscopic quantities of those cells that shared the common vertex.The projected lengths in Equation ( 9) are given by: The time derivatives in (8) are given by: With the reconstruction of Equations ( 8)- (11), the macroscopic flux Φ n+1 j,k is evaluated by: with the effective diffusion coefficients in Φ n+1,macro j,k as: The Φ n+1,micro j,k of the macroscopic flux (12) contributes from microscopic photons with free transport, while the Φ n+1,macro j,k accounts for the contribution from macroscopic diffusive flux.Therefore, the macroscopic flux (12) is a multiscale numerical flux, which bridges the microscopic radiant flux and macroscopic diffusive flux.
Upon calculating the Φ n+1,micro j,k in Equation ( 12) by the microscopic solver (which will be discussed in the next subsection), a coupled nonlinear system for the macroscopic quantities φ and ρ is established.This system is solved by the following source iteration Algorithm 1.

Microscopic Solver for the Isotropic-UGKP Method
Following the first-order expansion of the macroscopic quantities, Equation ( 6) can be written as: Equation ( 15) means that the radiation intensity at t n+1 is a combination of the free transport and diffusion radiation intensity whose relative contribution is determined by the e −λ∆t and 1 − e −λ∆t term.With the macroscopic quantities ρ n+1 j and φ n+1 j acquired through the macroscopic solver in algorithm 1, the microscopic radiation intensity at t n+1 can be recovered based on Equation (15).The e −λ∆t I transport (t n+1 , x , y ) in (15) contributes from freely streamed photons in the time step t ∈ [t n , t n+1 ], which can be further tracked in the next time step.Thus, the re-emitted and scattered photons are re-sampled from 1 − e −λ∆t I di f f usion n+1 j in (15) only.The numbers of the re-sampled photons are determined by: where E j,re−emitted and E j,scattered are the total energy of re-emitted and scattered photons of cell j, respectively; and w re f is the reference weight of photons.Then, the re-emitted and scattered photons are sampled in each mesh isotropically (pointed out by the 1/2π in Equation ( 15)) and uniformly.The photon-tracking process in the isotropic-UGKP method is similar to the traditional Monte Carlo method, except that those "interacted" photons will be immediately removed from the simulation.The subsequent behavior of these "interacted" photons is evaluated by the macroscopic solver in algorithm 1 instead of the meticulous tracking involved in the traditional Monte Carlo method.During particle tracking, the Φ n+1,micro j,k in Equation ( 12) is calculated by: where ∑ i stands for the summation of all photons i that are transported across the edge k, and w i is the weight of MC particle i.
The microscopic solver of the isotropic-UGKP method is summarized by the following Algorithm 2.
Algorithm 2 Photon re-sampling and tracking algorithm for microscopic radiation intensity evolution.
1: Photon re-sampling: 2: for all mesh cell do 3: Re-sample the re-emitted and scattered photons isotropically and uniformly in mesh cell.if Photon leaks out of the system or Photon has "interacted" with material then 12: Remove the photon from the simulation.end if 16: end for

Anisotropic-UGKP Method
In this subsection, we present the anisotropic-UGKP method for the gray radiative transfer system.We begin by highlighting the distinction between the isotropic and anisotropic systems and subsequently discuss the modifications made to transition from the isotropic-UGKP method to the anisotropic-UGKP method.

The Difference between Isotropic and Anisotropic System
For an anisotropic system, the scaled form of the time-dependent gray radiative transfer equations under the same condition in Section 2.1.1 are: where S = p Ω → Ω I Ω d Ω represents the scattering source of anisotropic scattering.The p Ω → Ω here is the scattering phase function, which depends on the incoming direction Ω and the outgoing direction Ω.Compared with (1) and (18), the only difference between isotropic and anisotropic scattering systems is the scattering source term, changing from ρ/2π into S.
As previously discussed in our theoretical work [41], the scattering phase function can be approximated by a finite series of polynomials: with C * l as the corresponding coefficients of polynomial Ω • Ω l , and θ representing the scattering angle between Ω and Ω.Based on (19), the scattering source S maintains a relationship with the moments of radiation intensity I through: where is the lth-order moments of I, and . The necessary modifications of the macroscopic and microscopic solver for the anisotropic-UGKP method, because of the scattering source changing from ρ/2π in an isotropic system to S in an anisotropic system, are discussed in the following subsections.

Macroscopic Solver for the Anisotropic-UGKP Method
In our prior work, we proved that the pure anisotropic scattering radiative transfer equation satisfies the diffusion equation in the optically thick limiting regime [41].The diffusion coefficient of pure anisotropic scattering systems is related to the average cosine of the scattering angle cosθ by: Therefore, the diffusion coefficient of systems containing both absorption-emission and anisotropic scattering processes can be derived through a similar method: Hence, the macroscopic equations and macroscopic flux of the anisotropic systems remain the same structure as those in the isotropic system in (3) and (12).But the effective diffusion coefficients in macroscopic flux should be adopted: which will approach diffusion coefficient (22) within the optically thick regimes.Then, the same source iteration method in Section 2.1.2is utilized to solve the macroscopic equations in the anisotropic-UGKP method.

Microscopic Solver for the Anisotropic-UGKP Method
The integral solution of radiation intensity around the center of edge k in the anisotropic system (18) follows the same formation as in (6) expect that the scattering source is changed into S: Therefore, the following formation can be derived under the first-order expansion of the macroscopic quantities: Similar to (15), the re-emitted and scattered photons are re-sampled from 1 − e −λ∆t I di f f usion n+1 j only.The re-emitted photons are re-sampled in a similar manner in the isotropic-UGKP method.The scattered photons, however, need to be sampled with the same angular distribution as S n+1 j .
Based on Equation ( 20), S n+1 j can be written as: which is related to moments J (l) n+1 j with 0 ≤ l ≤ L. Thus, one common way is constructing the UGKP method for all these moments above.In doing so, the macroscopic system must contain the macroscopic equation of all moments' components, and the corresponding macroscopic flux of each macroscopic equation need also be calculated.
As a result, the calculating costs will increase exponentially with the highest order L increasing (since the total components' number of J (l) is 3 l ).What is more, S n+1 j has a complex relationship with moments, making it challenging to directly sample the angular distribution from S n+1 j .
To overcome these issues, we propose a revised phase function as an equivalent to the calculation and angular sampling of S n+1 j , drawing on our prior theoretical work [41].This equivalent phase function, which will be derivated and discussed in Section 2.2.4 in detail, incorporates a factor f applied to all odd-order coefficients of the original phase function (19): The factor f varies between 0 and 1 in the diffusion and free transport limit regimes, respectively.
With (27), we introduce modifications to the photon tracking and re-sampling processes within the microscopic solver of the anisotropic-UGKP method.In the anisotropic-UGKP method, photons are freely transported in the time step t ∈ [t n , t n+1 ] first.At the end of time step t n+1 , photons are then sampled to be freely transported (with probability e −λ∆t in Equation ( 25)), absorbed (with probability 1 − e −λ∆t 1 in Equation ( 25)), or scattered (with probability 1 − e −λ∆t 1 λ cL ε s σ s ε in Equation ( 25)), which is based on the current mesh's cross-sections.The freely streamed and re-emitted photons are treated in a similar way as in the isotropic-UGKP method.The scattered photons, in contrast, will be redirected into a new direction, with the scattering angle sampled from (27).Since photons are freely transported during the photon tracking, the Φ n+1,micro j,k term of macroscopic flux is also adapted: With this adapted microscopic solver, there is no need to calculate the macroscopic fluxes and solve the macroscopic equations of high-order moments, significantly reducing the computational costs.The scattered photons are also easily re-sampled with ( 27) compared with directly sampling the angular distribution from S n+1 j .This adapted microscopic solver is summarized by the following Algorithm 3: Algorithm 3 Photon re-sampling and tracking algorithm for the adapted microscopic radiation intensity evolution.Sample the photon's status: freely transported, absorbed, or scattered.Remove the photon from the simulation.Re-sample the re-emitted photons isotropically and uniformly in mesh cell.if Photon leaks out of the system then 22: Remove the photon from the simulation.

23:
end if 24: end for

Derivation and Discussion of the Equivalent Multiscale Phase Function
In this subsection, we commence with the derivation of the equivalent multiscale phase function (27).Then, we address the normalization condition, non-negativity condition, and multiscale property of ( 27).
Proposition 1.The discretized scattering source S n+1 j satisfies: with p eq Ω → Ω as the equivalent phase function (27).The I transport here stands for the freely transported radiation intensity during a time step [t n , t n+1 ].
Proof of Proposition 1.Based on the definition of S and Equation ( 20), the discretized scattering source S n+1 j equals: With Equation ( 25), the components of moments J (l) have the form: where S (l) is the lth-order moments of S, and S (l) l) .The even-order moments have the following relationship for all flow regimes [41]: In the diffusion limit, we have [37,41]: Thus, the components of moments J (l) can be written as: Therefore, Equation ( 29) can be proved: where a reverse derivation of Equation ( 20) is used.Based on Equation (35), the angular distribution of scattering source S n+1 j can be equivalent to one single scattering process with (27), saving the computational costs.Proposition 2. The equivalent multiscale phase function (27) satisfies the normalization condition: Proof of Proposition 2. For the original phase function (19), it satisfies the normalization condition: What is more, the normalization condition of the phase function is only dependent on the even-order coefficients C * l=even .Hence, with identical even-order coefficients, the equivalent multiscale phase function satisfies the normalization condition as the original one.Proposition 3. The equivalent multiscale phase function (27) satisfies the non-negativity condition: Proof of Proposition 3. To substantiate the non-negativity condition (38), the range of factor f needs to be determined first.Taking L ε a = L ε s , the factor f has the following form: where ω = σ s /(σ a + σ s ) is the scattering albedo with 0 ≤ ω ≤ 1, and −1 < cosθ < 1.For all cases below with ω = 0, the factor f satisfies 0 ≤ f ≤ 1: For the special case with ω = 0, the scattering process does not exist, while another special case with cosθ = 0 stands for a phase function without any odd-order terms.Hence, the odd-order terms of ( 27) are not considered in either of these two cases.
In the optically thin regime, the equivalent multiscale phase function (27) degenerates to the original phase function (19); 2.
In the optically thick regime, the equivalent multiscale phase function (27) preserves the isotropy of the scattering source.
Proof of Proposition 4. As indicated in Equation ( 40), the factor f ranges from 0 to 1 in the diffusion and free transport limit regimes, rendering the equivalent phase function ( 27) a multiscale phase function.The factor f equals 1 in the optically thin regime as λ → 0. Thus, the equivalent multiscale phase function (27) aligns with the original phase function (19).With λ increasing, more frequent collision processes lead to a reduction in the anisotropy of the macroscopic scattering source, which is signified by the factor f decreasing.In the optically thick regime as λ → ∞, the factor f tends toward 0, causing all odd-order coefficients in (27) to approach 0. Therefore, the equivalent multiscale phase function (27) transforms into a phase function without any odd-order terms, preserving the isotropy of the macroscopic scattering source [41].

Summary for the Anisotropic-UGKP Method
For the anisotropic-UGKP method, the microscopic solver furnishes the free-streaming numerical flux for solving the macroscopic equations, and the solution provided by the macroscopic solver serves as the closure source for the microscopic intensity.The macroscopic and microscopic solver are closely coupled for solving the macroscopic energy and microscopic intensity.The anisotropic-UGKP method is concisely summarized as following Algorithm 4.

Algorithm 4
The algorithm for the anisotropic UGKP method.Stream all particles by the microscopic solver, and calculate the free-streaming flux (28).

4:
Apply the macroscopic solver to evolve the macroscopic quantities.

5:
Re-sample the re-emitted and scattered source photons and recover the microscopic radiation intensity with Equation (25).6: end for

Asymptotic Analysis
The asymptotic preserving (AP) property is important for the construction of a multiscale scheme.For the free transport and single-temperature diffusive regimes, the asymptotic preserving analysis is similar to that in our previous work [37].For brevity, we omit the repetition of those details here.In this subsection, we focus on demonstrating the asymptotic preserving property of the proposed anisotropic-UGKP method in the other diffusion regimes.
Proposition 5.The proposed anisotropic-UGKP method preserves the two-temperature diffusion system: Proof of Proposition 5.For the two-temperature diffusion regime or non-equilibrium diffusion limit, the Knudsen number ε approaches 0, with L ε a = ε and L ε s = 1/ε.The effective diffusion coefficients in (23) have the following orders: Then, the macroscopic numerical flux ( 12) for an anisotropic system has the following form: Substituting the macroscopic numerical flux (46) into the macroscopic Equation (3), we have Equation ( 47) is a standard nine-point scheme for the diffusion limit equation in ( 44) for an anisotropic system independent of the parameter ε.The second equation of ( 44) is also independent of the parameter ε under L ε a = ε and L ε s = 1/ε.Thus, the convergence of Equation ( 44) is inherently satisfied.Furthermore, the equivalent multiscale phase function (27) also maintains the isotropy of the angular distribution of macroscopic scattering source in the diffusion regime.The above discussions affirm that the proposed anisotropic-UGKP method preserves the two-temperature non-equilibrium diffusion limit.Proposition 6.For systems with absorption-emission and scattering processes equally dominated, the proposed anisotropic-UGKP method preserves the following diffusion equation: Proof of Proposition 6.For systems with absorption-emission and scattering processes equally dominated, the Knudsen number ε approaches 0 and L ε a = L ε s ∼ 1/ε.The effective diffusion coefficients in (23) have the following orders: (49) Therefore, the macroscopic numerical flux (49) has the following form: Apply the asymptotic analysis to the macroscopic Equation ( 3), and we have the O ε −2 equation: Therefore, the macroscopic numerical flux (50) can be rewritten as: Coupling the macroscopic numerical flux (50) and the macroscopic Equation (3), we derive the O ε 0 equation: It shows that Equation (53) becomes a standard nine-point scheme for the diffusion limit Equation (48) coupled with (51).The equivalent multiscale phase function (27) keeps the isotropy of the angular distribution of the macroscopic scattering source as well.Hence, the above discussions confirm that the proposed anisotropic-UGKP method preserves the diffusion limit equation for systems with absorption-emission and scattering processes equally dominated.

Numerical Tests
Although our anisotropic-UGKP method is proposed for thermal radiative transfer problems, it is also suitable for neutron transport problems with pure scattering (since they have the same formation).Therefore, we present three numerical examples to validate the proposed anisotropic-UGKP method in this section, including a 1D neutron problem and two 2D radiation heat transfer problems.In the following examples, the unit of the length is taken to be centimeter (cm), the mass unit is gram (g), the time unit is nanosecond (ns), the temperature unit is kilo electronvolt (keV), and the energy unit is 10 9 Joules (GJ).Under the above units, the speed of light is 29.98 cm/ns, and the radiation constant a is 0.01372 GJ/(cm 3 keV 4 ).

One-Dimensional (1D) Neutron Problem
In this 1D neutron problem, we consider the following initial condition and two different cross-sections: σ, x ∈ [0, 0.35) ∪ (0.65, 1.35) ∪ (1.65, 2], with the computational domain as x ∈ [0, 2].This neutron problem is a pure scattering problem with ω = σ s /(σ a + σ s ) = 1.To verify the AP property of the anisotropic-UGKP method, we investigate two scenarios: one with σ = 1 and the other with σ = 10 3 in (54).These cases represent the free transport and diffusion limit, respectively.The neutron speed is taken as 1, and the simulation time is 1 and 100 for the transport and diffusion limit.
For case B with σ = 10 3 , the cross-section spans a vast range from 0 to 10 5 , covering a large range from the free transport regime to the diffusive regime.Thus, the computing is very challenging.The isotropic scattering case is simulated first to verify the codes, whose results are given in Figure 2 and compared with the reference results in [42].Subsequently, the anisotropic scattering case with the following F1 phase function is conducted to verify the anisotropic-UGKP method: with C as the normalization factor.Figure 2 also displays a comparison between the simulation results and reference results in [42] for this F1 scattering case.As can be seen, good agreement is observed for isotropic and F1 scattering cases in both kinetic and diffusion regimes.The cross-section at x = 1 is perpetually 0 for case B, switching the anisotropic-UGKP method to full MC simulation around x = 1, and therefore leading to the statistical fluctuations around x = 1.

Two-Dimensional (2D) Radiation Heat Transfer Problem
In this 2D radiation heat transfer problem, a square with the side length of 1 enclosed by four boundaries is considered, as illustrated in Figure 3.The bottom wall is kept hot with a non-dimensional temperature of T 1 = 1, while all other walls are kept cold with T 0 = 0.The medium is also kept cold with T 0 = 0 initially.Given the focus of this work on the anisotropic-UGKP method, most of the cases considered here are pure scattering (with scattering albedo ω = 1), where the scattering process is the most significant.The isotropic scattering case is also simulated here to verify the codes, while four anisotropic cases with different phase functions are explored to verify the proposed anisotropic-UGKP method.The four anisotropic phase function varying with scattering angle θ changing are illustrated in Figure 4, and their Legendre expansion coefficients are listed in Table 1.As can be seen, the anisotropic phase functions vary significantly with the scattering angle θ changing, which leads to different transfer behavior in this problem.The normalized net heat flux Q y /(acT 4 1 /4) in the y-direction and radiation intensity G/acT 4  1 along the centerline (x = 0.5) are compared with the reference results in [43].Figures 5 and 6 show the results for isotropic and different anisotropic phase function cases, respectively.It can be clearly observed that the computational results from our anisotropic-UGKP method align closely with the reference results for all cases.Furthermore, the influence of the anisotropic scattering on the energy transfer is also clearly shown: the forward scattering phase function transports more radiation heat than the isotropic case, which leads to a flatter distribution of radiation intensity; meanwhile, less radiation heat is transferred with the backward scattering phase function compared with the isotropic case, causing a steeper distribution of radiation intensity.These results verify the accuracy of the anisotropic-UGKP method in modeling 2D anisotropic-scattering models.In addition, the F2 scattering case with varying cross-sections ranging from 0.01 to 100 is simulated.The results of the normalized net heat flux in the y-direction along the centerline (x = 0.5) with different cross-sections are presented in Figure 7 compared with the reference results from [43] (σ = 0.01∼10) and [44] (σ = 40∼100).As can be seen, the results of the anisotropic-UGKP method agree well with the reference results for all cases.Furthermore, Figure 8 displays the normalized net heat flux in the y-direction along the centerline (x = 0.5) with different scattering albedos ω, which also shows good agreement with the reference results in [43].With ω = 0, the energy absorption reduces the radiative heat transfer, which leads to the steepest distribution of Q y /(acT 4  1 /4).The absorption process becomes less important with ω increasing, causing a flatter distribution.The above results confirm the capability of the present anisotropic-UGKP method in simulating the absorbing and anisotropic scattering systems of all regimes.

Two-Dimensional (2D) Radiation Heat Transfer Problem with Collimated Incidence
We also take a 2D radiation heat transfer problem with collimated incidence to further test the present anisotropic-UGKP method.In this problem, the simulation domain is similar to the problem above.The difference lies in that the bottom wall is also kept cold, and a collimated beam is normally incident through the top wall with intensity I c , as illustrated in Figure 9.The scattering albedo is also set to be ω = 1 for most of the cases considered here, and the same anisotropic phase functions are simulated.The normalized net heat flux Q y /|I c | in the y-direction and radiation intensity G/|I c | along the centerline (x = 0.5) for different anisotropic phase functions are shown in Figure 10 compared with the reference results in [45].The computational results from our anisotropic-UGKP method agree well with the reference results for all anisotropic cases.In addition, the normalized net heat flux values under different cross-sections and different albedos are given in Figure 11 and Figure 12, respectively.The results of the anisotropic-UGKP method also show good agreement with the reference results from [45].Hence, the above results clearly show the capability of the present anisotropic-UGKP method for solving radiative transfer problems with absorbing and anisotropic scattering media in all regimes.

Conclusions
In this work, we developed the anisotropic-UGKP method based on our prior research.The macroscopic solver of the proposed anisotropic-UGKP method directly solves the macroscopic diffusion equations of the anisotropic scattering system.In addition, we proposed a revised multiscale phase function as an equivalent to the calculation and reconstruction of macroscopic scattering source S, covering from the free transport limit to the diffusion limit regimes.With this multiscale equivalent phase function, the scattered photon re-sampling process in the microscopic solver of the proposed anisotropic-UGKP method is simplified as one single scattering process.Therefore, the proposed anisotropic-UGKP method does not need to calculate the macroscopic fluxes and solve the macroscopic equations of high-order moments, significantly reducing computational costs.Furthermore, the anisotropic-UGKP method also has the asymptotic preserving (AP) property in both optically thin and thick regimes.A set of radiative transport problems, including 1D neutron problems and 2D radiation heat transfer problems, are conducted to validate the proposed anisotropic-UGKP method.All results agree well with the reference results from prior studies, showing the accuracy and efficiency of our anisotropic-UGKP method across a broad range of computational conditions.The scheme and code will be applied in the high-energy density physics engineering applications.

Figure 1 .
Figure 1.A cell j of the generalized quadrilateral mesh with c j as the cell center and p k,m as the center of edge k.The length of the edge k is l k .The two vertexes of edge k are p k,1 and p k,2 , and the unit normal and tangential vector are n k and τ k , respectively [37].

10 :until
Photon leaks out of the system or Photon has "interacted" with material or Photon reaches the end of time step 11:

13 :
else if Photon reaches the end of time step then 14:Continually simulate the photon in the next time step.15:

4 :if 5 : 6 :
Photon is freely transported then Continually simulate the photon in the next time step.else if Photon is absorbed then 7:

8 :
else if Photon is scattered then 9:Sample the scattering angle from(27) and redirect the photon.10: end if 11: end for 12: for all mesh cell do 13:

20 :until
Photon leaks out of the system or Photon reaches the end of time step 21:

1 :
Initialize the flow field and sample photon particles.2: for Simulation time less than the final time do 3:

Figure 3 .
Figure 3. Schematic of 2D radiation heat transfer problem with a hot wall.

Figure 4 .
Figure 4.The four anisotropic phase functions varies with scattering angle θ changing in the 2D radiation heat transfer problem.

Figure 5 .
Figure 5.The normalized net heat flux Q y /(acT 4 1 /4) in the y-direction and radiation intensity G/acT 4 1

Figure 9 .
Figure 9. Schematic of 2D radiation heat transfer problem with collimated incidence.

Table 1 .
The Legendre expansion coefficients of the four anisotropic phase functions used in the 2D radiation heat transfer problem.