A domain decomposition method for a two-scale transport equation with energy flux conserved at the interface, Kinetic and Related Models

When a linear transport equation contains two scales, one diffusive 
and the other non-diffusive, it is natural to use a domain 
decomposition method which couples the transport equation with a 
diffusion equation with an interface condition. One such method was 
introduced by Golse, Jin and Levermore in [11], where an 
interface condition, which is derived from the conservation of 
energy density , was used to construct an efficient non-iterative 
domain decomposition method. 
   
In this paper, we extend this domain decomposition method to 
diffusive interfaces where the energy flux is conserved . Such 
problems arise in high frequency waves in random media. New 
operators corresponding to transmission and reflections at the 
interfaces are derived and then used in the interface conditions. 
With these new operators we are able to construct both first and 
second order (in terms of the mean free path) non-iterative domain 
decomposition methods of the type by Golse-Jin-Levermore, which will 
be proved having the desired accuracy and tested numerically.

1. Introduction. The one dimensional transport equation is Ψ(t, x, µ) is the energy density of waves (or particles) at position x and time t, and µ = cos θ where θ is the angle between the x-axis and the moving direction of waves (or particles). σ(x, µ, µ ′ ) > 0 is the scattering cross-section at x which is usually symmetric in µ and µ ′ , c(x, µ, µ ′ ) is emission rate at x, i.e. the average number of where µ 1 and µ 2 are respectively the cosines of incidence and transmission angles. µ 3 , the cosine of reflection angle, satisfies µ 3 = −µ 1 by the reflection law, see Figure  1.
In [13], we developed an efficient numerical method for such a problem. The basic idea, following those of [14,15], was to build the interface transmission and reflection into the numerical flux. In this paper, we extend the method into the case where the mean free path (or the dimensionless Knudsen number) has two different scales, one is the transport scale-where the Knudsen number is of O(1)-and the other diffusive scale -where the Knudsen number is very small. In the diffusive regime the transport equation (1.1) can be replaced by the computationally more efficient, macroscopic diffusion equation. Through a domain decomposition method that couples the transport equation with the diffusion equation, we obtain a multiscale method which is much more efficient than the approach using the full transport equation (1.1) in the entire domain.
To better illustrate the basic idea, in this paper, we will concentrate on the steady solution Ψ(x, µ) of (1.1) in the case of isotropic scattering and slab geometry for The method in the present paper can be generalized to more general boundary conditions, but will not be pursued here. For more physical background of (1.1), and (1.3)-(1.4), see for example [7]. We are concerned with the numerical computation of (1.3)- (1.4) in the case where the scattering cross-section σ have different orders of amplitude in different subdomains of (x L , x R ). Such situations arise when the background medium is made of (very) different materials. Let the small parameter ǫ be the nondimensional mean free path, i.e. the ratio of the average distance between two successive collisions of a particle with the background medium to the size of the domain x R − x L . Specifically, we consider the background medium made of two different materials with an interface located at x M ∈ (x L , x R ), and denote the background density and sound speed to be ρ 1 , v 1 in (x M , x R ) and ρ 2 , v 2 in (x L , x M ). At x M , the scattering cross-section σ and emission rate c have discontinuities, given by the parameter ǫ: where we assume for some constants γ * and γ * ; and Notice that these assumptions exclude the case where c takes the value 1, namely, the purely scattering case.
Since one can define a new space variable by the so-called "Optical thickness" in the radiative transfer [6], , there is no loss of generality in assuming σ to be piecewise constant as in (1.5).
The mean-free path is small in the region (x M , x R ), hence the solution Ψ ǫ is µ-isotropic (independent of µ) to leading order, and governed by the diffusion approximation of (1.3): see Lemma 2.7 below. We name the domain (x M , x R ) as "the diffusion region". Contrarily, Ψ ǫ is µ-dependent in the domain (x L , x M ), which will be named as "the transport region".
In [11] (followed by the numerical work [22]), Golse, Jin and Levermore developed an efficient domain decomposition method for the two-scale transport equation (1.3)-(1.5) with energy density conserved at the interface. In their method, the transport and diffusion equations are decoupled to the leading order, thus a first order (in ǫ) domain decomposition completely avoids any iteration. To obtain a second order accuracy in ǫ, a correction in the interface condition is needed, which requires consequently just one iteration.
Nevertheless, in many other applications of the transport equation, like the propagation of energy density for waves in heterogeneous media with weak random fluctuation in the high frequency regime [21], the energy flux is conserved at the interface. It is known that the radiative transfer equation (1.1) is the weakly coupling limit (or semiclassical or high frequency limit) of linear high frequency waves through weakly and highly fluctuating random (for example, stationary Gaussian) media, which can be derived by the Wigner transform for the underlying wave or Schrödinger type equations and then taking the high frequency limit, see for example [21,18]. In such problems, the energy flux is conserved at the interface [1] where the local index of refraction contains discontinuities, and the waves satisfy the Snell's Law of Refraction (1.2). Note that the conservation of energy density considered in [11] also conserves energy flux when the local index of refraction is continuous.
For other models of radiative transfer or radiation-matter coupling (especially in time-dependent case and in moving medium), energy flux conservation is of essential importance for simulating the energy transport, exchange and redistribution in globe system of physics [5,16,20].  In this paper, we develop a domain decomposition method which satisfies the energy flux conservation at the interface. Our new contribution is to derive some new operators to account for (specular or diffuse) scattering at the interface, which are needed in the interface conditions used to couple the transport and the diffusion equations via a domain decomposition method of [11]. We then give rigorous L 1 error estimates, similar to the L 2 -estimate in [11]. Numerical examples are used to show the performance of the method.
The result of this paper differs from [11] in that the scattering is different at the interface, thus different interface behavior exhibits which requires different interface conditions and a slightly different convergence theory. It differs from [13] in that here we couple a transport with a diffusion equation while in [13] two transport equations-with discontinuous cross-sections-are coupled.
The paper is organized as follows. In section 2, we present some basic transmission and reflection operators used in [1], and then derive several new operators that are needed for the interface conditions. We then use these operators in both the first and second order (in ǫ) domain decomposition methods in section 3 and prove their L 1 convergence. The numerical discretizaton is given in section 4 and some numerical experiments are conducted in section 5 to show the performance of the methods. We conclude the paper in section 6.
2. The interface conditions. In this section, we derive the interface conditions needed for the domain decomposition method. Several new operators will be introduced for the interface conditions. The properties of these new operators will also be studied.
, i, j = 1, 2 map the functions defined on Γ j + onto the functions on Γ j − , where the spaces L p µ are defined by L p µ (X) = L p (X; |µ|dµ) and R ii corresponds to reflection, while R ij ,i = j describes the transmission across the interface. See Figure 2 for an illustration of these operators.
Since the energy flux is conserved across the interface at x = x M , the reflection and transmission operators must satisfy the following condition: . The usual Fresnel reflection-transmission coefficients are given by [17,8]: where ρ 1 , v 1 and ρ 2 , v 2 are the background density and sound speed of the diffusion region (x M , x R ) and the transport region (x L , x M ) respectively, and v 1 , v 2 , µ 1 and µ 2 satisfy the Snell's Law of Refraction (1.2). For example, if one considers the case that the waves at the interface undergoes specular reflection and diffuse transmission, then the operators R ij in the interface condition (2.1) can be expressed in terms of F ij as (2.4) For physical and analytical reasons, we assume the operators R ij satisfy the following assumptions (a part of Hypothesis 2.1 in [1]): 1. The operators R ij are bounded on L ∞ (Γ j + ) and on L 1 µ (Γ j + ), and the operators R ii are nonamplifying, i.e. R ii

They are linear, and positive in the sense that
The reflection operators are not totally reflecting, that is, there exist positive constants η j and subsets I j ∈ Γ j + of positive measure such that, for all nonnegative functions h of support supp(h) ⊂ I j , Note that these assumptions will guarantee Lemma 2.1 below, and allow for a sufficiently wide class of operators (Remark 2.2 in [1]).

2.2.
The response operator and the Chandrasekhar H-function. The domain decomposition method needs the following two basic quantities for the transport equation. First we define the Chandrasekhar H-function via an integral equation ( The density of particles emerging from a purely scattering half-space is Note that the response operator R also appears in the diffusion region as shown in Figure 2.
The study of Chandrasekhar's H-function and of the response operator R is referred to [9] and [10]. Another presentation given by stochastic processes can be found in [4].

Derivation of the new interface operators.
Here we derive the operators S and T in Figure 2, which is needed in for the interface conditions of the domain decomposition methods.
Let I be the identity operator, then by one has The well-definedness of the operators S and T , which depends on the invertibility of the operator I − R 11 R, is given in the next lemma by Bal and Ryzhik (Corollary 4.4 in [1]). First, we give a basic property of a positive linear operator.
Lemma 2.2. If P is a positive linear operator, then for any function Proof: Notice that For future use, we need the following three properties of the operators. Lemma 2.3. The operators R 11 , R 21 and R satisfy the following property: where L p µ means L p (·, |µ| dµ). Note that (2.12) and (2.13) come from equation (4.7) of [1], and (2.14) comes from equation (5.13) of [1].
for some constants σ * , σ * and c * and for each x ∈ [x L , x M ]. Consider the transport equation The following positivity lemma was due to Bardos [3]: The idea of the proof for this lemma mostly follows those in [11] except here we work on L 1 rather than L 2 . Since σ(x) ≥ σ * > 0, c(x) ≤ c * < 1, due to Lemma 2.4, we have f n → 0 in

Proof
. By (2.15), this will imply µ∂ x f n → 0 in which is contradicting to (2.19). Therefore the second inequality is strict.
Note that the second equality comes from (2.14) in Lemma 2.3, and the last equality comes from (2.2) by takeing i = 2 and j = 1. Because S is a linear positive operator, by Lemma 2.2 we have Therefore, S L 1 µ = 1, and S is a contraction map.
For convenience, in the next lemma, we restate some results of Proposition 4.1 in [11] about the convergence of the diffusion approximation with boundary conditions.
In the diffusive region x ∈ [x M , x R ], the scaled transport equation is Let Ψ ǫ (x + M , ·) be the right limit of Ψ ǫ (x, ·) at the interface x = x M , and Lemma 2.7. Error estimates for the diffusion approximation: • At order O(ǫ), define N ǫ 0 as the solution of

22)
where R is defined in (2.6) and D 0 is a positive constant.
• At order O(ǫ 2 ), define N ǫ 1 as the solution of where λ is called the extrapolation length, given as follows: where D 1 is a positive constant. Here in the first inequality of (2.21) and (2.24) we only consider x M < x ′ M < x < x ′ R < x R to avoid the complication due to the boundary layers near x = x M and x = x R which does not change the results (see [11]).
3. The coupling algorithms. In this section we revise the coupling algorithm in [11] with the new interface operators we derived in the previous section.

The O(ǫ) coupling.
• Step 1. Solve the steady transport problem where S is the decomposition operator defined in (2.10). • Step 2. Using the density Φ 0 provided in Step 1, solve the diffusion problem where H is the Chandrasekhar function defined in (2.5).
Remark 3.1. Notice that Θ 0 is independent of µ, thus there will be boundary layers of thickness O(ǫ) at the interface inside the diffusive region, which are not captured by this coupling method.
The procedure has the same advantage as the one in [11]. The transport region and the diffusion region are completely decoupled at the O(ǫ) approximation so that the solution does not depend on ǫ. Actually there is a flux of order O(ǫ) at the interface from the diffusion region to the transport region, which requires the coupling of these two regions to get the next order of accuracy O(ǫ 2 ).

3.2.
The O(ǫ 2 ) coupling. In order to obtain a second order scheme, one needs another coupling procedure.
• Step 1. Solve Step 2. With the results of Φ 1 , we continue to compute where λ is given in (2.23).

3.3.
Convergence of the two coupling algorithms. In this subsection, we establish the convergence of the two coupling algorithms introduced in the previous subsection. The proof follows those in [11] except that here we work on the L 1 space rather than the L 2 space in [11]. We will skip many details of the part of the proof which resembles those in [11], and only give the parts where things are slightly different. Then and, for each be the left limit of Ψ ǫ (x, ·) at the interface x = x M , and For µ ∈ (0, 1], define φ ǫ (µ) : . Notice that (2.22) in Lemma 2.7 will imply Ψ + 1 = R(Ψ − 1 )+O(ǫ) actually in (2.7). Therefore the boundedness of the operator (I − R 11 R) −1 , thanks to Lemma 2.1, will lead to φ ǫ L ∞ ≤ C 0 ǫ, where C 0 is a positive constant.
The remaining proof resembles that of Theorem 3.2 in [11] except that the space L 2 µ used in [11] is replaced by L 1 µ . The details are omitted here.
4. Numerical algorithms. We choose the same numerical schemes as those used in [22] except the treatment of the interface conditions.
First, we use the Gauss quadrature to approximate the collision operator. The Gaussian quadrature points in [0, 1] are given by µ m , with the corresponding weight A m , for m = 1, · · · , M , satisfying Then, we discretize the operators R ij and R i , i, j = 1, 2, numerically by the Gauss quadrature. Note that the Gauss quadrature will provide a higher accuracy in the computation of the Chandrasekhar H-function. Since the operators are linear at the interface, the discrete forms of the operators will be constant matrices. Denote the constant matrix P ij to be the discrete form of R ij and the constant matrix Q to be the discrete form of R, then the discrete form of S in (2.10) is simply given by S = P 21 Q(I − P 11 Q) −1 P 12 + P 22 . Similarly, we denote the discrete form of T in (2.10) as and the discrete form of (I − R 11 R) −1 R 11 as In addition, one should be careful when discretizing R 12 and R 21 , because they map the discretized function values on µ m , m = 1, · · · , M , into the ones that are not one of µ k 's. Therefore one needs to do interpolation during the discretization. Note that this interpolation is not needed when discretizing the operator S, although it involves with R 12 and R 21 , since it maps the discretized function values on µ m , m = 1, · · · , M , into the ones on the same µ m 's. Next, we use the upwind scheme for the transport equation, and center difference for the diffusion equation.
Step 1. In domain [x L , x M ], set • Using the Gauss quadrature again, we obtain the discrete boundary conditions Φ m 0 = F L (µ m ), m = 1, · · · , M, where S ml is the (m, l)-component of the matrix S in (4.1).
with the boundary condition where U lk is the (l, k)-component of the matrix U in (4.2).
Step 1. Observing that (3.3) differs from (3.1) only on the right hand side of the equations, so one just needs a discretization of ∂ x Θ 0 (x M ), and we use the second order formula Step 2. We use fictitious points to discretize the diffusion boundary conditions in order to gain a second order accuracy in the diffusion region (to be the same order of accuracy as the interior of the domain [x M , x R ]). Using then the left boundary condition is discretized as where U lk and V lk are the (l, k)-components of the matrix U in (4.2) and the matrix V in (4.3) respectively. Using the first equation of (3.4) at x M , By combining these two equations, one gets the second order approximation of the diffusion boundary condition at x M : Similarly, the diffusion boundary condition at x R is discretized by In the interior of the domain, the same second order centered difference scheme is used as in (4.4). We take ∆x = ∆µ = 0.005 and ∆t = ǫ∆x/5 to numerically simulate the evolutionary transport equation (1.1) with the boundary conditions (1.4) for any time t, the two-scale isotropic coefficients (1.5) and the interface condition (2.1). We run the codes with the initial condition Ψ(t = 0, x, µ) = 1.0 until the final time T = 30 to make sure we have got the steady solution. Note that the mesh size of the domain decomposition method is ten times as big as the one used in the direct simulation of the evolutionary two-scale transport equation.
• O(ǫ) results. Figure 3 (left) shows that when ǫ is getting smaller, the steady solution Ψ ǫ is getting closer to the O(ǫ) numerical solution Ψ 0 . Besides, Table  1 gives the L 1 errors of the average energy density between the steady solution Ψ ǫ and the O(ǫ) numerical solution Ψ 0 for different ǫ values, in which one can see that when the ǫ is reduced by a half, the error is almost reduced by a half, which coincides with the first order convergence given by Theorem 3.2. Moreover, we also give a 2d-plot of the O(ǫ) numerical solution Ψ 0 in Figure  3 (right). • O(ǫ 2 ) results. The 2d plots of the steady solution Ψ ǫ and the O(ǫ 2 ) numerical solution Ψ 1 ǫ at ǫ = 0.10 are given in Figure 4. Figure 5 gives the comparison of the average energy density between the O(ǫ) numerical solution Ψ 0 , the O(ǫ 2 ) numerical solution Ψ 1 ǫ and the steady solution Ψ ǫ at ǫ = 0.10 and ǫ = 0.20. The error plots are shown in Figure 6.
• Efficiency. To show the efficiency, we do the numerical simulation for the case ǫ = 0.0001, where the "true" solution is approximated by the O(ǫ 2 ) algorithm using ∆x = 0.01 which is good enough based on Theorem 3.3 and the above numerical simulation. As we can see in Figure 7, the O(ǫ) decoupling algorithm using the coarse mesh size ∆x ≫ ǫ still provides good approximation results compared to the "true" solution, especially in the diffusive region [0, 1].    6. Conclusion. By introducing the new interface operators (2.10), we extend the domain decomposition method by Golse, Jin and Levermore in [11] to the diffusive interfaces where the energy flux is conserved. It inherits the merit of the method proposed in [11], which decouples the transport and diffusion equations to the leading order and completely avoids any iteration. By an O(ǫ) correction in the interface condition, a second order accuracy is obtained, which requires only one iteration. The desired accuracy for this domain decomposition method is proved and tested numerically. Besides, we also show the numerical efficiency of this method in the case where the mean free path is pretty small in the diffusion region.
In the future, we will extend the method to high space dimensions. This can be done by implementing a similar interface condition in the direction normal to the interface, as was done in a related case of elastic wave in [12].