DIFFUSIVE LIMITS OF 2D WELL-BALANCED SCHEMES FOR KINETIC MODELS OF NEUTRON TRANSPORT

. Two-dimensional dissipative and isotropic kinetic models, like the ones used in neutron transport theory, are considered. Especially, steady-states are expressed for constant opacity and damping, allowing to derive a scattering 𝑆 -matrix and corresponding “truly 2D well-balanced” numerical schemes. A first scheme is obtained by directly implementing truncated Fourier–Bessel series, whereas another proceeds by applying an exponential modulation to a former, conservative, one. Consistency with the asymptotic damped parabolic approximation is checked for both algorithms. A striking property of some of these schemes is that they can be proved to be both 2D well-balanced and asymptotic-preserving in the parabolic limit, even when setting up IMEX time-integrators: see Corol-laries 3.4 and A.1. These findings are further confirmed by means of practical benchmarks carried out on coarse Cartesian computational grids.


Dissipative (2 + 1)-dimensional kinetic models
Consider a damped kinetic model, where (x) ≥ 0 is the opacity, for x = (, ) ∈ R2 and v = (cos , sin ).In the special case where (x) ≡ 1, the model studied in [25] is clearly recovered, hence our intention hereafter is to extend these findings toward the more general model (1.1).We shall be especially interested in its "four-stream approximation" [14,28] with diagonal microscopic velocities, }︂ , which yield the following system, (in standard notation) :=  + +  − +  + +  − .This discrete kinetic model, within a convenient rescaling of variables, relaxes towards the damped heat equation, and rotating the axes for simplicity, it comes Adding the four microscopic balance laws, and subtracting the first (second) with the third (fourth) ones, This last equation implies that, formally, | ± −  ± | = () and (1.3) holds: see Section 5 in [28] for related rigorous results.

Plan of the paper
Following the roadmap proposed in [25], we intend to study numerical approximations of (1.2) endowed with both 2D well-balanced (WB) and asymptotic-preserving (AP) properties.To this end, two distinct numerical processes will be introduced: the first one, given in Section 3, strongly relies on the data of two-dimensional steady-states for (1.1).Such steady-states are derived in Section 2, following original ideas given in [5,6]: in both papers, it is shown that solutions of stationary elliptic equations yield, thanks to Laplace transforms, microscopic steady-states of related kinetic models.A second numerical scheme is proposed in Section 4, being essentially an exponential modulation of the one given in [25] for the special case  ≡ 0. Its WB and AP properties are studied, in the light of what was previously done in Section 3. A main difference between both schemes is that only the second one appears to be unconditionally positivity-preserving; numerical tests displayed in Section 5 reveal that it is slightly more diffusive, though, especially in kinetic regime with stiff parameters.Concluding remarks are given in Section 6.

Dissipative isotropic scattering with adsorption
In the original model (1.1), the expression of the stationary regimes is obtained (by passing to the limit  → +∞ in the method of characteristics), see e.g.[12], page 34 as a Laplace transform (again, letting  = 1), see [6]: Again, another integration in v produces a Fredholm integral equation on , 1 The Fourier coordinates "diagonalize" the discrete scattering operator.
in which the circular integral on (x − v) will reduce to the pointwise value (x) if  solves a convenient elliptic differential operator [16,34].Yet, the diffusive approximation suggests the "modified Helmholtz problem", which solution satisfies (see e.g.[3], Sect.2.4 for more details on Green and modified Bessel functions), By inserting (2.8) into (2.6), one gets that the constant  (different of the 1d case, see [8]) must satisfy: Being a solution to (2.7), the macroscopic density  (expressed in polar coordinates and separating variables) rewrites as a Fourier-Bessel series in any disk of radius  > 0, so that the corresponding stationary microscopic density  (x, v) follows by computing Laplace transforms of Bessel functions because trigonometric functions depend only on  (and v = (cos , sin ) ∈ S 1 ).
As  = 1, an interesting relation comes out, so that the general term proceeds by the summation formula (2.10), For instance, since  2 = 1 −  2 = (1 + ), we get the expression of the first terms in (2.9):
Remark 2.1.Both (2.4) and (2.12) express steady-states only in a special case where x = (cos , sin ) and v = ±(cos , sin ) are colinear.This is enough to build up the -matrix (3.3) and to derive the corresponding 2D-WB scheme in the form (3.4).However, this is less general than what was found in equation (3.1) of [25], where explicit steady-states were given for arbitrary arguments (x, v) ∈ R 2 × S 1 .

Two-dimensional well-balanced scheme
Hereafter, a uniform Cartesian grid is used, with ∆ = ∆ and ∆ > 0 a time-step.Standard notation is

Derivation of the dissipative 𝑆-matrix
Following the same procedure as Section 3.1 of [25], a 2D -matrix,   = M  −1 , is deduced from the expression of kinetic (both incoming and outgoing) steady-states (2.12), where , a 4 × 4 matrix which columns are orthogonal to each other, and The Fourier polynomial (2.12) of both incoming/outgoing states yields an analogue of (2.5), which gives back (2.5) when  → 0, so that  → 1, because In order to express the -matrix, the inverse of  is needed, , so that, by setting  0 =  − 0 / + 0 , and so on.A discrete version of (2.3) is deduced, (see Fig. 2) where each function   () rewrites as,
Proof.Under the CFL restriction, the marching scheme (3.4) is a convex combination, and preserves positivity as soon as all the entries of   are nonnegative, hence  small enough.The dissipation of  1 and  ∞ norms is a consequence of the symmetry of   where both lines and columns add to less than one.For consistency, when  ≪ 1, by (3.2), it comes ) (like conservative case  = 0), so that both  1 ,  2 behave like ,  in equation (3.4) of [25], and which yields that (3.4) is consistent with (1.2) as  → 0. To establish the "2D well-balanced" property, notice first that, in a given cell centered in  − 1 2 ,  + 1 2 a sufficient condition for a state to be invariant by (3.4) is, in the prescribed computational cell.Assume now that there locally exists a truncated steady-state of the form (2.12), f (σ(cos , sin ), ±(cos , sin )),  ∈ (0, 2), which samples for   =  4 +   2 ,  ∈ {0, 1, 2, 3} match all the 8 values located around the computational cell, In that case, by construction, the -matrix (3.3) maps the 4 "incoming values" (v • x < 0) into 4 "outgoing values" (v • x > 0) which all correspond to f , so that , and this exactly means that the discrete "4 velocities" restriction of f given by (2.12) is left invariant by (3.4).
Remark 3.2.Oppositely, if the opacity is stiff,  ≫ 1 (local hydrodynamic limit), the sign of the entries in the -matrix isn't obvious.These expressions allow to get necessary conditions for positivity in (3.4): for instance, in order to get  0 −  2 ≥ 0 in the limit of infinite stiffness  ≫ 1 with  ≃ 1,  ≃ 0, one should have Under the same assumptions, if (x), (x) are constant in the whole computational domain, then any truncated steady-state of the form (2.12) restricted to four velocities, Proof.Given  > 0, if all parameters are constant, any steady-state f (x, v) can be split into "patches" of the form (2.12) inside the disks   of radius  and centered in )︁ , ,  ∈ Z 2 (see Fig. 2) inside which the proof of Theorem 3.1 applies.

Consistency with damped diffusive limits in 2D
In order to study its diffusive limit, we shall again decompose    in a form directly inspired by [24,25], By analogy with coefficients defined in Section 4.1 of [25], and neglecting 2 2 ℐ 2 (. ..), ) and the decomposition of    =  0 +   1 follows: Accordingly, an IMEX scheme emerges from treating implicitly the stiff terms in (1)/.
Yet, (3.14) is reached by shifting indexes like in (3.11), along with -adding both first and last two equations, so as to get the preservation of averages, , and  +1 , =   , ; -subtracting both first and last two equations while defining macroscopic fluxes, -and finally taking advantage of the elementary relation, )︂ .
Remark 3.5.The assumptions of Corollary 3.4 are similar to the ones of Theorem 3.1.Indeed, to ensure that, in a given cell, holds, it is sufficient that there exists a local steady-state of the form (2.12) for the rescaled equation By sending formally  → 0, having divided by 4  , coefficients split into along with (3.11), which yields, Yet, invoking (3.15) in both (3.12) and (3.13) yields diffusive fluxes and dissipative terms, 2 2 2 )︀ )︁ .
We choose to set up the horizontal/vertical discrete Laplace operator with 2 )︀ )︁ , so that another "diagonal" Laplace operator, involving   −   = () (for fine grid) appears.Clearly, both f and g solve similar damped diffusion equations, if both (x) and (x) are smooth.Accordingly, if initial data are "well-prepared", they are likely to stay so for later times because the Maxwellian gap f  , − g  , satisfies a parabolic linear equation, too.As  = f + g, assuming f , = g , for all ,  ∈ Z 2 brings the limiting scheme, )︀ Assuming that the parameters aren't stiff, meaning ∆ √  ≪ 1, so the aforementioned asymptotic scheme is consistent, as ∆ → 0 with the diffusion approximation (1.3).Moreover, assuming that (x), (x) are positive constants, the scheme (3.16) rewrites as follows, )︀ , and under the fine-grid approximation (3.2),   +   ≥ 0. Consequently, an elementary manner to make it unconditionally positivity-preserving is to treat implicitly all the dissipative terms, the ones multiplied by   .

Positive, composite two-dimensional scheme
Consider again the following kinetic system with adsorption (1.1), in the particular case where σ, κ ∈ (0, 1) are positive constants.Thus, the change of variable  =   σκ allows to convert (1.1) into the conservative equation (formerly studied in [25]), In the general case, (x) and (x) vary, so this computation cannot be made on a global scale.However, by assuming that  and  are frozen in each disc of radius  = ∆/ √ 2 and centered in )︁ , see Figure 2, this trick can still be applied locally and this is enough for our purposes.Accordingly, by means of an exponential modulation of the scheme given in equation (3.5) of [25], a time-marching process for (1.2) can be deduced, where )︁ , and the -matrix S(),  ≥ 0 corresponds to the limit  → 0 in (3.3) and reads (see [25], Eq. (3.4)):
Proof.Under the CFL restriction, the scheme (4.1) realizes a nonnegative combination; moreover, all the entries in the symmetric matrix S((1 − )) are nonnegative and positivity is preserved.Yet, assume  is small enough so that the CFL restriction gives a time-step small enough to linearize the exponential modulation, Hence, since  ≪ 1 and ∆ ≪ 1, the time-marching scheme (4.1) is approximated by , and this expression is consistent with system (1.2),

Numerical steady-states
Since it involves an exponential modulation in the time variable, the ability of (4.1) to preserve continuous steady-states is not obvious.Accordingly, denote outgoing and incoming states, respectively, as follows: Accordingly, the scheme (4.1) rewrites simply , so that Yet, assume ∆ is small enough to linearize, , thus another scattering matrix appears, An interesting question is to relate it with the "exact" -matrix in (3.4), at least for small values of .Notice that at this stage, we can assume that  = 1 without loss of generality thanks to a simple rescaling of .First, the entries of the -matrix in [25] are the limits when  → 0 of (3.4), recalled in (2.5): Yet, in S  , the 's must be multiplied by 1 − , so its entries agree with (3.5): This computation shows that the steady-states of (4.1), at least when  ≪ 1 (so exponential and Bessel functions behave like polynomials), are consistent with the "exact ones" preserved by (3.4) as long as necessary linearizations are licit.

Consistency with the damped diffusion limit
The diffusive scaling of (1.2) corresponds to (3.6) at the discrete level; thus, applying the exponential modulation to the IMEX scheme given in equation (4.2) of [25] with shorthand notation where the scattering -matrix in (4.1) splits into Maxwellian and diffusive parts, .

•
Oppositely to (3.7) and (3.8), here, both coefficients   and   have the same limit, so that the asymptotic scheme will have a 5-points stencil, like the one of [25].An index-shift in (4.3) yields: )︀ Summing former equations and passing to the limit in the coefficients   ,   gives the asymptotic scheme, which is clearly consistent with (1.3).

Rigorous diffusive limit for constant parameters
For ( , ) a real sequence, we shall use the notations In order to make the above formal computation rigorous, we first derive estimates for the IMEX scheme (4.4).
Lemma 4.4.Let us assume that the CFL parabolic condition holds.Then, the scheme is nonnegative.Moreover, if the initial data are bounded in  1 ∩  ∞ , then the sequences   ±, , and   ±, , are bounded in  ∞ , and we have the dissipation inequality Proof.Inverting the block-diagonal matrix in the left hand side in (4.4), we get )︀ )︂ , and similar expressions for   +,+1 , and   −,+1
The second condition being more restrictive than the first one, we only have to verify this latter condition, i.e.
Constant coefficients.We assume now that (x) and (x) are constant and we study rigorously the diffusive limit when  → 0. We recall the definitions Then the scheme (4.4) rewrites )︀ Lemma 4.5.Let ,  ∈ R + be constants; under the parabolic CFL condition (4.7), the scheme (4.9) is TVD.
Proof.As in the proof of Lemma 5.4 from [25], inverting the block-diagonal matrix in the left hand side in (4.9) gives .Under the condition (4.7), each coefficient in the left hand side are nonnegative.Then by linearity, and after using a triangle inequality, we obtain with similar expressions for ⃒ ⃒ ⃒ Adding all these expressions, we obtain )︁ ]︂ .
Summing over  and , we get after shifting the indexes, Theorem 4.6 (Asymptotic preserving property).Let us assume that the parabolic stability condition holds and that initial data, independent of , are smooth enough so that, (recall notations in (4.6)) for some constant  ≥ 0. )︀ which satisfy: Moreover, we have, for all  ∈ N * , By summing the equations (4.11) and (4.12), we deduce the following result: Corollary 4.7.Let us denote   = 1 4 (f  + g  ).Under the same assumptions as in Theorem 4.6, we have and Proof.By the same token as in the proof of Lemma 4.5, we deduce that the sequences (︀   ±, , )︀ , and )︀ are Cauchy sequences with respect to  in ℓ 1 .Thus, we deduce that, when  → 0, they converge to some limit sequences denoted respectively (︀  ±, , )︀ , and )︀ .We may pass into the limit  → 0 in (4.9), we deduce , where we recall the definition of the matrix Thus, the limit verifies  +,+1 From the expressions of   and   in (4.8), we deduce that   → 1  and   → 1  when  → 0. Then passing into the limit in (4.13) and (4.14), recalling that  ±, , = 1 2 f  , and  ±, , = 1 2 g  , , we obtain (4.11) and (4.12).Then, denoting   , = f  , − g  , , we get straightforwardly from (4.11) and (4.12) and From the assumptions on the initial data in Theorem 4.6, we have ‖∆f 0 ‖ 1 + ‖∆g 0 ‖ 1 ≤ .Then, for all  ∈ N, we have We may inject this latter inequality into (4.15).Taking the absolute value, under the condition (4.10), and summing over  and , we deduce for some nonnegative constant .Applying a discrete Gronwall inequality, we get

Numerical results
All the practical tests displayed hereafter are conducted on the unit square coarsely gridded with 32 × 32 points in order to demonstrate both the robustness and accuracy of the time-marching schemes.Results are displayed on Figure 3: despite the stiffness of the benchmark, results are quite similar, except for (4.1) showing slightly more numerical viscosity (compare the macroscopic densities on the first row of Fig. 3).Both time-marching algorithms succeed in stabilizing correctly, being a consequence of 2D well-balanced properties.
Results at time  = 0.15 are given in Figure 4: both schemes appear to be quite similar.Especially, the Maxwellian gap f − g is practically identical for both algorithms, the one produced by (3.11) being slightly bigger.Macroscopic densities obtained from (3.11) are 5% higher compared to (4.3), which might confirm that the exponential modulation is endowed with a slightly higher numerical dissipation.

Conclusion and outlook
Two numerical schemes endowed with both 2D well-balanced (WB) and asymptotic-preserving (AP) properties, extending the one previously in [25], were studied in this paper.The first one involves a -matrix directly built from the expression of exact steady-states for (1.1), namely truncated Fourier-Bessel series (2.12).The resulting scheme was proved to be 2D-WB; a novelty is that, in Corollary 3.4, it is also proved that its IMEX reformulation (3.11) is endowed with similar properties, too.A drawback of (3.4) is that it isn't unconditionally positivity-preserving.Accordingly, a simpler scheme (4.1), relying on an exponential modulation of the one given in [25], was proved to be positivity-preserving, but endowed with weaker well-balancing features.Its IMEX reformulation (4.3) is again shown to be consistent with asymptotic damped diffusion behavior, see (4.5).Such diffusive limit is even rigorously established by means of convenient estimates.The next stepping stone consists in applying the same program to two-dimensional kinetic models rendering biased velocity redistribution, like in semiconductor and chemotaxis dynamics modeling, see [7].Particular equations belonging to this class of models are still simple enough to admit steady-states which can be again expressed as Fourier-Bessel series, hence permitting to derive 2D-WB numerical schemes; in such cases, a drift-diffusion equation takes the place of (2.7).This somehow extends the one-dimensional ideas of [23,24] toward problems in higher dimensions.A more challenging question is about the extension of these 2D-WB schemes to kinetic models involving more than 4 microscopic velocities.In such case, two distinct possibilities coexist: -in a first one, the disks inside which the elliptic equation is solved in order to derive the -matrix by means of a Laplace transform should grow bigger (see Fig. 5), like for instance  = ∆ with a 8-velocity model.The corresponding 8 × 8 -matrix will keep trace of steady-states made of Fourier polynomials involving 8 terms (instead of 4 as in (2.12)), so that 8 "incoming states" will be needed; -a second one can be drawn by following the canvas of so-called "tailored schemes" (see [26,27]).Namely, at each point of coordinates (, ) ∈ Z 2 , the stationary kinetic equation is solved in a disk of radius  = ∆, see Figure 6 so as to derive "outgoing states" at its center, instead of along its boundary (as done here by means of (2.3) and in [7,25] as well).This asks for a slightly different -matrix, but may result into an easier way to process kinetic models with more than 4 discrete velocities.
Concerning 3D extensions, the theoretical ideas of [5,6] still apply but the "Poisson kernel" of the elliptic equation may have a different form (see e.g.[29]) because of the supplementary dimension; then cubature formulas may be needed in order to take full advantage of kinetic steady-states, again expressed as Fourier polynomials, and involving spherical harmonics, see e.g.[31].
Appendix A. Second order IMEX "Midpoint rule" scheme Following [9,32], other IMEX time-integrators, possibly high-order can be substituted to (3.10).A well-known example is the second order in time "midpoint IMEX rule", which in our notation, rewrites as follows: 2  -the first matches the IMEX scheme (3.10)where ∆ is changed into ∆/2, hence if .
Accordingly, the second order midpoint IMEX scheme is 2D well-balanced inside its domain of stability.
In order to illustrate the behavior of this second-order IMEX integrator, the same example as Section 5.2 was set up, with  = 10 −2 though.Results are displayed on Figure A.1: the mass is slightly bigger for the midpoint rule IMEX scheme, but the main differences show up in the Maxwellian gaps f − g which have a quite different shape, despite both have roughly the same amplitude.The second-order IMEX scheme isn't stable when  ≪ ∆.

Figure 6 .
Figure 6.Disk of radius  = ∆ with "outgoing states" at its center in a "tailored" framework.