Diffusion approximation in a radiative transfer model for astrophysical flows

In this work, we present the diffusion approximation model for radiative transfer when we deal with optically thick astrophysical flows. Since the initial model is high CPU time demanding when dealing with its numerical approximation, solving this simpler system can provide a low cost strategy for the simulation of radiative media. We then use a finite-volume algorithm coupled with an implicit scheme for radiative contributions to solve this simplified system. Numerical experiments in the one-dimensional and two dimensional cases are presented to validate our numerical strategy and to prove the relevance of this asymptotic model.


Introduction
The comprehension of realistic astrophysical phenomena has been the subject of many challenging studies for decades, with various applications in stellar physics such as the study of jets in the formation of young stars involving matter accretion, stellar winds and dynamics of supernova remnants. Within this framework, a wide complexity occurs for different reasons. First, these phenomena involve both hydrodynamics and radiative processes (see [16], [18], [25]). The latter ones are driven by the propagation of photons at the speed of light, which means that typical velocities completely differ between hydrodynamic and radiation scales. This causes a strong coupling between fluid and radiation since the behavior of the medium at a given point may depend on quantities which need to be evaluated on a large spatial zone (the physical effects are said to be nonlocal). Moreover, in the general case the analytical expression for the physical values under study is out of reach due to the complexity of the governing model. In this case, one can use numerical methods in order to simulate the physical phenomena as accurately as possible. Consequently, mathematical algorithms that are involved have to be consistent with the desired solutions (see [21]). In particular, they have to tolerate discontinuities that may spatially propagate (known as "shock waves") and they have to produce physically relevant results (namely entropic solutions). Finally, one has to implement robust methods in order to handle extreme ratios that may be encountered for typical quantities such as density, pressure, radiative energy, time and space scales.
In this paper, we deal with an asymptotic model for the radiative transfer known as the diffusion approximation for which the mean free path of photons is assumed to be small in terms of the characteristic length. In this case, radiative effects are considered as local, leading to a simplified model which has been the subject of various studies in the last decade (see [2], [3], [9], [16], [23], [26]). We then propose a numerical strategy based on finite volume methods that are widely used for the resolution of hyperbolic problems, with an implicit coupling between hydrodynamics and radiative processes.
This paper is organized as follows: in Section 2, the physical model is discussed, with a focus on the diffusion approximation. In Section 3, the numerical method used for the computation of approximate solution is detailed. Some numerical results obtained in dimensions 1 and 2 are presented in Section 4. Finally, discussion and perspectives are addressed in Section 5.

Physical model
We first intend to present the general physical model for radiative transfer as well as the optically thick approximation leading to the diffusion approximation.

Pure hydrodynamic model
We here recall the classical hydrodynamic description of a given gas where we look for the time evolution of the physical quantities such as density ρ, pressure p, velocity field u = (u, v) T and total energy E per volume unit, that will be assumed to depend on time t and space x. From the conservation of mass, momentum and total energy, the governing model is the classical Euler system In order to get a closed system, one needs the equation of state that enables to evaluate the pressure in terms of the hydrodynamic quantities. For ideal gases, one has where γ stands for the adiabatic index (for diatomic gas, γ = 1.4). Note that it is possible to express the total energy with respect to the temperature T as where c v denotes the heat capacity at constant volume. System (1) has been intensively studied both from a mathematical and numerical point of view. It is well-known to belong to the class of hyperbolic systems requiring to use numerical methods capable to compute discontinuous states as well as entropic solutions.

Radiative model
We now take into account the interaction between photons and the matter which involve absorption, emission and scattering processes (see [16]). When a beam of radiation passes through a material, a loss of energy occurs due to the opacity of the material and can be expressed in terms of extinction coefficient χ(t, x; n, ν) defined such that an element of material of length dl and cross section dS, oriented normal to a beam of radiation having specific intensity I = I(t, x; n, ν), propagating along n into solid angle dΩ in frequency band dν in an interval time dt removes an amount of energy δE = χ(t, x; n, ν)I(t, x; n, ν) dl dS dΩ dν dt.
The emission coefficient η(t, x; n, ν) of the material is defined such that the amount of energy released by a material element of length dl and cross section dS propagating along n into solid angle dΩ in frequency band dν in an interval time dt is δE = η(t, x; n, ν) dl dS dΩ dν dt.
In the same way, η can be separated into a thermal part η t and a scattering part η s such that η(t, x; n, ν) = η t (t, x; n, ν) + η s (t, x; n, ν).
The transfer equation is obtained computing the difference between the amount of energy between position x and x + ∆x and time t and t + ∆t which equals the difference between the amount of energy created by emission and the amount absorbed by the material. Consequently, one has I(t + ∆t, x + ∆x; n, ν) − I(t, x; n, ν) dS dt dΩ dν = η(t, x; n, ν) − χ(t, x; n, ν)I(t, x; n, ν) dl dS dt dΩ dν.
From the intensity I evaluated at a given frequency ν, it is possible to define radiative energy E, radiative flux F and radiative pressure P as the first moments of I, that is x; n, ν) n dΩ, From these expressions, we compute the total radiative quantities integrating all the elementary contributions over all frequencies Integrating the transport equation (3) for all possible frequencies and directions, one gets x; n, ν) − χ(t, x; n, ν)I(t, x; n, ν) dΩ dν := −cG 0 .
At local thermal equilibrium, the source terms arising in these equations can be expressed in the frame where particles are always at rest, the comoving frame (with subscript "0"), as involving the first radiative constant a R , the gas temperature T and Planck, energy and flux mean opacities defined by where B(ν, T ) is the Planck function. To simplify (G 0 , G), we assume that the radiation has a blackbody spectrum which means that E(ν) ∝ B(ν, T ), leading to κ E = κ P . In order to close the system, we use the relation P R = DE R where the Eddington tensor (see [11]) is given by The Eddington factorχ is obtained by minimizing the radiative entropy (see [4], [17]) and writesχ One obtains with this closure relationship the so-called M1 model for the description of the radiative transfer. Finally, one gets the system of radiation hydrodynamics by coupling this model with the Euler system Let us remark that in Eqs. (7), radiative quantities are calculated in the laboratory frame while source terms are expressed in the comoving frame. It is possible to switch from one frame to the other one using Lorentz transformations (see [15]). In the laboratory frame case, left-hand terms are unchanged but right-hand terms are modified as a consequence of fluid motion; indeed, one has to take into account the Doppler shift and the aberration of photons (see [14]). In the comoving frame case, right-hand terms are unchanged but advection terms appear in the left-hand side of equations (see [12]). Among all numerical codes designed for the simulation of radiation hydrodynamics, HADES code based on the M1 model has been developed in LUTH (see [13]). In the following section, we introduce an asymptotic regime leading to numerical routines used in HADES within the framework of local radiative effects.

Diffusion approximation
It is well-known that two extremal regimes can be observed in radiative transfer, leading to simplified models depending on the ratio between the mean free path λ of photons and the characteristic length L of the phenomenon under study. The first one is known as the optically thin case for which λ L. It is then possible to derive an Euler system involving a cooling function as a source term in the equation of the energy conservation.
In the second one referred as the optically thick case, we assume that λ L. Radiative effects are assumed to be local since emitted photons are absorbed almost instantaneously and do not come into play outside the considered spatial point. In this specific case, the radiative flux F R in the comoving frame expresses as where κ R is the Rosseland mean opacity defined by Since the radiation field is assumed to be isotropic, the radiative pressure rewrites . Keeping all terms up to O (u/c) and dropping at least some terms that are of order O (u/c) 2 in the streaming limit and the static and diffusion limits, one obtains the new system of radiation hydrodynamics in the diffusion approximation where radiative quantities are expressed in the comoving frame (see [9]) as where we have removed the subscript "0".

Numerical method
We now aim to solve numerically Syst. (8). For hydrodynamics, great efforts have been paid in order to derive robust numerical methods especially designed to tolerate discontinuities and high Mach numbers. We write (8) under the general form with We perform our computations on a cartesian bidimensional uniform spatial grid. Considering a time discretization with a prescribed time step δt, we then globally split the resolution of Eq. (9) using Strang algorithm (see [20]) to treat alternatively the homogeneous part and the nonhomogenous part at each time increment. Let us mention that S(U ) involves divergence terms since when considering the full first-order contribution, the eigenvalues of the Jacobian matrix cannot be easily calculated, which makes the numerical resolution quite delicate.

Homogeneous part
Setting S(U ) = 0 in (9), we deal with the hyperbolic system We use a Strang splitting to solve (10) separately in x and y directions and we perform a finite volume approach.
Defining respectively δt and δx as the time step and the space step, and noting t n = nδt and x j±1/2 = (j±1/2)δx, we calculate the average value U n j of the solution of (10) at time t n and over the cell For this aim, we use the MUSCL-Hancock scheme (see [24]) which is second-order accurate in time and space. The evaluation of numerical fluxes requires to compute the solution of the Riemann problem, that is from an initial data consisting in two constant states separated by an initial discontinuity, for linearized equations. We then use HLLC solver (see [22]) to compute the flux F hllc between each interface x j+ 1 2 for which we add the new component for E R . The classical form of the finite volume scheme expresses as where the boundary flux at the interface separating constant states In (11), S L , S * and S R are the wave speeds corresponding respectively to the eigenvalues u − c, u and u + c of the system. Moreover, the two intermediate states U * L and U * R are given by and where u * = S * . For both intermediate states U * L and U * R , the fifth components are obtained using Rankine-Hugoniot jump conditions (see [8], [19]) as for classical Euler equations with HLL type Riemann solvers.

Nonhomogeneous part
We now deal with the reduced system that is treated using a fully time implicit discretization, from the state U n = (ρ n , (ρu) n , (ρv) n , E n , E n R ) T resulting from the previous step. Since the first component of the total system is left unchanged, we have ρ n+1 = ρ n . Furthermore, we choose to use (2) to compute the total energy. In such a way, we obtain an implicit system in which u, T and E R are sought at time t n+1 . Since we perform a dimensional splitting, we here describe our algorithm when only x derivative operators are considered. The time discretized system writes In (14), the unknown quantities continuously depend on spatial location (x, y). Since this is a fully nonlinear system due to kinetic and thermal contributions, we seek a linearly implicit approximation of the nonlinear terms. First, we have dropping off the quadratic contribution in u n+1 −u n . Similarly, one finds that (T n+1 ) 4 −3(T n ) 4 +4(T n ) 3 T n+1 and plugging these expressions into (14) gives a simplified system that turns out to be linear in u n+1 , T n+1 and E n+1 R . After some computations, we find that E R solves that takes the general form M E n+1 R = f . Once this linear equation is solved, it is possible to compute u n+1 as u n+1 = u n − δt 3ρ n ∂ x E n+1 R and T n+1 as the solution of the linear equation involving a diagonal matrix. We finally update the value of total energy E n+1 as E n+1 = ρ n c v T n+1 + 1 2 ρ n (u n+1 ) 2 . We now spatially discretize the two operators ∂ x cδt 3κ n R ∂ x and β n ∂ x in (15) on the one-dimensional uniform grid x 1 , . . . , x N . The first one is approximated at x i using the formula . The second one is discretized using δx that enables to mimic the propagation across characteristic curves of the transport operator β∂ x for any sign of the speed. Consequently, we express β n ∂ x as Taking the complete two-dimensional case, we would obtain pentadiagonal contributions that require more computational efforts for the numerical implementation. This explains the reason why we performed an alternate direction splitting that consists of one-dimensional iterative resolutions. Since the system that has to be inverted is tridiagonal but not symmetric, we use a LU decomposition algorithm. Once again, the second-order Strang strategy is used, giving a good compromise between computational efficiency and numerical accuracy.

Numerical results
In this Section, we present various numerical simulations that have been performed in order to validate our numerical code. We first present a shock-tube test in the context of radiative flows in strong equilibrium regime, that appears as a generalization of the usual hydrodynamic benchmark test. We then focus on a radiative shock raised by a left-sided piston. We also simulate subcritical and supercritical shocks and finally present twodimensional simulations, the first one in a radial geometry and the second one with an initial data constant on four quadrants.

Shock-tube test
We begin with a test that has been presented in [26]. This is a purely one-dimensional test for an ideal gas with adiabatic index γ = 5/3 and mean molecular weight µ = 1 in the relationship between the pressure and the temperature p = ρk B T /(µm H ) where k B is the Boltzmann contant and m H the hydrogen mass. A discontinuity initially separates two constant states with the following values: ρ L = 10 −2 kg m −3 , T L = 1.5 × 10 6 K, u L = 0 m s −1 and ρ R = 10 −2 kg m −3 , T R = 3 × 10 5 K, u R = 0 m s −1 . Consequently, the initial density is uniform and the contrast between left and right temperatures is 5. In this test, it is assumed an initial thermal equilibrium in the gas, which means that E r = a R T 4 . We impose constant opacities set to κ P = 10 8 m −1 and κ R = 10 10 m −1 . In this regime, radiative pressure cannot be neglected. We plot in Figures 1 and 2 profiles of density, velocity, total pressure p + P rad and radiative energy that have been computed with the following numerical parameters: we considered 400 gridcells for the normalized domain [0, 1] with a initial discontinuity localized at the center of this domain. The Courant number in our simulations has been taken equal to 0.9. Transparency  Numerical results that have been obtained are in good agreement with the ones presented in [26]. Furthermore, the structure of the solution is similar to the one found in a purely hydrodynamic regime. In particular, a depletion occurs for the density as seen in the left-hand side of Figure 1. A shock wave is generated and propagates with velocity v s ≈ 300 km s −1 . Let us mention that the elementary waves are magnified and propagate faster in this radiative context than the hydrodynamic case.

Radiative shock in the laboratory frame
We now focus on the formation of a radiative shock based on a laboratory experiment in a medium of krypton for which γ = 5/3. Initially, the gaz is at rest with density ρ 0 = 0.168 kg m −3 and temperature T 0 = 1 eV. Furthermore, the gas is assumed to be initially in equilibrium with radiation, which implies that E R,0 = a R T 4 0 . To raise the shock, we define an artificial piston with uniform density ρ p = 10 3 ρ 0 , velocity v p = 20 km s −1 and temperature T p = T 0 . We perform a parallel plan simulation using 2.000 × 5 cells with space step δx = δy = 10 −6 m, where the piston initially fills 20 cells in x-space. Since we do not take into account ionization processes in our model, we roughly estimate it by dividing per 10 the mean molecular weight of krypton, which gives µ = 8.3798. Moreover, we assume that opacities are constant and equal to the space step. A comparison between results obtained with the M1 model and the diffusion approximation model at time t = 40 ns is shown on Figure 3. One can observe a very good quantitative and qualitative accordance between the two models, either for density and radiative energy. In this simulation, the maximal Mach number is found to be around 8. Let us notice that computational time is approximately 12 hours for a parallelized computation of 16 processors for the M1 model with HADES code and 28 minutes for a sequential computation using the diffusion approximation model. It has to be pointed out that using the M1 model, the CFL condition restricted significantly the time step in the hyperbolic part due to the speed of waves which depends on the speed of light. Consequently, it is worth dealing with the diffusion approximation model for optically thick flows rather than the full radiative transfer model for which computation times are much greater.

Subcritical and supercritical shocks
This benchmark has been introduced by [5] and was used to validate several numerical codes for radiation hydrodynamic models (see [1], [2], [6], [7]). Depending on temperatures T − and T + respectively observed at left-hand side and right-hand side of the shock, two cases occur: if the fluid velocity is sufficiently low, then T + T − and the shock is said to be subcritical, whereas for large enough fluid velocity, then T − ∼ T + and the shock is said to be supercritical. In both cases, a temperature overshoot is found at the shock location. We have here simulated a cold medium with γ = 1.4 and µ = 1 in the initial configuration ρ 1 = 7.78 × 10 −7 kg m −3 , T 1 = 10 K for two different velocities: u 1 = −6 × 10 3 m s −1 for the subcritical shock and u 1 = −2 × 10 4 m s −1 for the supercritical shock. Once again, the media is initially in equilibrium with radiation. The opacity has been taken constant and equal to κ = 3.1 × 10 −8 m −1 . Our computations have been performed on the spatial domain 7 × 10 8 m discretized with 300 cells with a CFL number equal to 0.1. We impose a reflexive boundary condition for the velocity at the left-hand side of the computational domain. Snapshots of the solutions are  presented in Figure 4 at different times. Curves are plotted in the moving frame (that is with respect to x − u 1 t) in order to make a comparison with the ones obtained in the literature. The total CPU time was 9.4 s for the first simulation and 8.7 s for the second one. In both cases, our results are in good agreement with the ones obtained in the cited papers. For the subcritical shock, the increasing of T − produces a flux of order σ R T 4 2 that preheats the material ahead the front shock and locally generates a temperature overshoot. For the supercritical case, the overshoot temperature is also observed at the interface between pre-shock and post-shock states with similar temperatures.

Radial test
We also performed a radial simulation on the cartesian domain [−8, 8] × [−8, 8] in a medium characterized by γ = 5/3, µ = 1 and that is assumed to be initially at rest with uniform density ρ 0 = 10 −8 kg m −3 . We impose the Gaussian radiative energy E R (x, y) = exp(−(x 2 + y 2 )) and assuming that the gas is in equilibrium with radiation imply that initially E = ρk B (E R /a R ) 1/4 /(γ − 1). The numerical value of the density ρ 0 has been chosen in such a way that E ∼ E R . We run the simulation with two sets of opacities κ P = κ R = 10 3 m −1 and κ P = κ R = 10 5 m −1 , for a 400 × 400 spatial discretization. The CFL number has been taken equal to 0.9 and the simulation has been performed until time t = 10 −3 s.
No analytical solution is available for this test. However, some of tendencies observed are in line with the ones expected. In Figure 5 and Figure 6, we clearly notice that waves propagate with larger velocities in the radiation case than in the hydrodynamics one. More precisely, velocities are enhanced for smaller opacities. This may suggest that this velocity does not only result from the hyperbolic part of the system where opacities do not come into play. In Figure 6, it is observed that the total pressure becomes larger in the central part of the domain which can be explained by the initial contribution of the Gaussian radiative pressure in this region. For smaller opacities, a homogenization of the total pressure occurs as the consequence of the diffusion term in the equation governing the radiative energy, which increases when the Rosseland opacity decreases. The same effect is encountered for the temperature.   Figure 6. Trace on y = 0 for different κ P = κ R = κ: κ = 10 5 m −1 (triangle), κ = 10 3 m −1 (asterisk) and pure hydrodynamic case (disk). Left: total pressure. Right: temperature.

Quadrant test
We finally perform a two-dimensional computation inspired by [10] where the initial data consists in constant states (see Figure 7 We deal with an overdensity and overpressure between zones 1, 3 and zones 2, 4 with symmetric velocity field (an horizontal and a vertical one for zones 2 and 4 and a diagonal one for the zone 3). The radiative energy is initialized to E R = a R T 4 for a gas characterized by γ = 1.4 and µ = 1. One can refer to Figure 7 (left) for a better understanding of the initial profile. Using this set of initial values, the resolution of the pure hydrodynamic problem exhibits shock waves between the different states with a nonconstant one that develops inside an increasing central domain (see Figure 7, right). Note that the interfaces between zones (1, 2), (2, 3), (3,4) and (4,1) are propagating with speeds given by Rankine-Hugoniot one-dimensional values. The solution is computed using 400 × 400 gridcells and the CFL number has been taken equal to 0.9. In Figure 8 are viewed the plots obtained for the density, first in a purely hydrodynamic configuration (left) and compared with the results in the diffusion regime (right), considered with the same initial states. In both cases, it can be observed the formation of the central structure delimited with circular shock waves. It means that even in the diffusion case, the general structure of the solution mimics the one in the absence of radiation. However, the density found in the diffusion case involves waves that once again propagate faster than in the hydrodynamic case. Moreover, the solution exhibits two dense regions in the poles of the structure in presence of radiation. This accumulation is also present when viewing the profile of radiative energy (see Figure 9, left). A slight density overshoot and a depletion effect at the center of the intermediate zone can be noticed. The plot of the trace of the density in the diagonal line that is presented in Figure 9 (right) for different values of opacities clearly confirms this effect. Qualitatively, this appears to be the same effect as the one observed in one-dimensional simulations: once again, the radiation magnifies both the values and the speed of propagation of density fronts. Let us finally mention that the cartesian crossed-line structure that is observed around the intermediate state in Figure 8 (right) is not a consequence of the directional splitting used in the method: the same simulation performed with a π/4-angle rotation has given the same result. Indeed, these secondary waves are generated by the radiation processes.

Discussion and perspectives
In this work, we have investigated the optically thick approximation that turns out to be a simplification of the classical M1 model for radiative transfer which is well-known to be costly in terms of numerical resolution. This reduced system of equations involves both an hyperbolic part as well as a diffusion term in the equation governing the radiative energy. A numerical code has been designed with a splitting algorithm that treats differently each contribution, involving a totally implicit scheme for the handling of the radiative part. Numerical simulations have shown a good compromise between accuracy of computed results and CPU time. In particular, benchmark tests have given the expected solutions and the same profiles as the ones obtained with the M1 model are recovered. Furthermore, all our simulations have shown the influence of radiation on the hydrodynamics, whether on the acceleration of wave velocities or in the hydrodynamic quantities. In a forthcoming study, this model will be simulated in more realistic astrophysical configurations where numerical results could be compared with experimental data.