An asymptotic-preserving dynamical low-rank method for the multi-scale multi-dimensional linear transport equation

We introduce a dynamical low-rank method to reduce the computational complexity for solving the multi-scale multi-dimensional linear transport equation. The method is based on a macro-micro decomposition of the equation. The proposed numerical method uses the low rank approximation only for the micro part of the solution. The time and spatial discretizations are done properly so that the overall scheme is second order accurate and asymptotic-preserving (AP); that is, in the diffusive regime, the scheme becomes a macroscopic solver for the limiting diffusion equation and is automatically low rank. We demonstrate the accuracy and efficiency of the proposed low rank method by a number of two-dimensional examples.


Introduction
The linear transport equation models particles such as neutrons or photons interacting with a background medium. This integro-differential equation is widely used in many science and engineering disciplines [3,4]. The linear transport equation belongs to the class of kinetic equations and is consequently posed in a five-dimensional phase space (3D in physical variable and 2D in angle or normalized velocity). This, in particular, implies that its full numerical simulation can be extremely expensive. The situation is further complicated if the scattering strength varies in space by several orders of magnitude, i.e., the equation is stiff in certain region of the domain and non-stiff elsewhere, then an explicit numerical scheme must resolve the smallest collision length.
Recently, a class of dynamical low rank methods has been introduced to solve high dimensional kinetic equations such as the Vlasov equation [8,9,12], the Boltzmann-BGK equation [7], and radiation transport equations [27,6]. Motivated by these advances, we develop an efficient numerical method to solve the multi-scale multi-dimensional linear transport equation. Employing a low-rank approximation 2 The linear transport equation and its macro-micro decomposition We are interested in the following time dependent linear transport equation in diffusive scaling: where f = f (t, x, v) is the distribution function of time t, position x = (x, y, z) ∈ Ω x ⊂ R 3 , and velocity v = (ξ, η, γ) ∈ S 2 which is confined to the unit sphere 1 . v denotes the integration over S 2 with respect to v. σ S (x) ≥ σ S min > 0 and σ A (x) ≥ 0 are the scattering and absorption coefficients, and G(x) is a given source term. Finally ε is the rescaled collision length, which can range between the kinetic regime ε ∼ O(1) and the diffusive regime ε ≪ 1.
The density ρ = 1 4π f v is defined as the angular average of f . In the limit ε → 0, ρ satisfies a diffusion equation which can be seen via the Chapman-Enskog expansion. Indeed, (2.1) can be written as On the other hand, taking 1 4π v of (2.1) yields which, upon substitution of (2.2), becomes with the diffusion matrix D given by (2.5) Therefore, as ε → 0 the limit of (2.1) is the diffusion equation (2. 6) In the macro-micro decomposition [21], we write f as where ρ is the macro part of the solution and g is the micro part. Note that g v = 0. Substituting (2.7) into (2.1) and taking 1 4π v , we obtain Subtracting (2.8) from (2.1) yields (2.9) The coupled system (2.8) and (2.9) is the macro-micro decomposition of the linear transport equation (2.1). In the limit ε → 0, we have from (2.9): which, when substituting into (2.8), yields the same diffusion equation (2.6). 3 The dynamical low-rank method for the linear transport equation We first constrain g(t, x, v) to a low rank manifold M such that where r is called the rank and the basis functions {X i } 1≤i≤r and {V j } 1≤j≤r are orthonormal: with ·, · x and ·, · v being the inner products on L 2 (Ω x ) and L 2 (S 2 ), respectively. With this low rank approximation, (2.8) becomes For (2.9), we write Equation (3.4), however, does not uniquely specify the dynamics of the low-rank factors X i , S ij , and V j . We therefore impose the following gauge conditions [18]: Let us emphasize that the resulting dynamics of g is independent of the specific gauge conditions chosen. However, using (3.5) is convenient as it allows us to easily obtain evolution equations in terms of the low-rank factors. To that end, we now project the right hand side of (3.4) onto the tangent space of M: where the orthogonal projector P g can be written as Using (3.7) and the gauge conditions we can in principle derive evolution equations for X i , S ij , and V j . However, this process requires inverting the matrix S = (S ij ). Since an accurate approximation mandates that S has small singular values, the resulting problem is severely ill-conditioned. Thus, we will use the projector splitting scheme introduced in [24]. For a corresponding mathematical analysis see [17]. This scheme has been extensively used in the literature, see e.g. [8,27,23], and extensions to various tensor formats have also been proposed [26,25]. The main idea is to split equation (3.6) into the following three subflows This is particularly convenient as for the first subflow V j is constant (in time), for the third subflow X i is constant, and for the second subflow both X i and V j are constant. Thus, we can write where After solving each subflow we use a QR decomposition to obtain X i and S ij from K j and S ij and V j from L i , respectively.

A first order in time scheme
Our goal is to solve the coupled system (3.3) and (3.6) using the projector splitting integrator outlined in the previous section. We now proceed by deriving the evolution equations corresponding to the subflows given by equations (3.8)-(3.10). (3.12) • L-step: (3.14) Therefore, for the overall system, we can construct a simple first order in time scheme. Suppose at time step t n , we have (X n i , V n j , S n ij , ρ n ). To obtain the solution (X n+1 i , V n+1 j , S n+1 ij , ρ n+1 ) at t n+1 we proceed as follows: 1. K-step: Solve (3.12) for a full time step ∆t, update from (X n i , V n j , S n ij ) to (X n+1 i , V n j , S ij ) using ρ n . Specifically, given K n j = r i=1 X n i S n ij , we discretize (3.12) using a first order IMEX scheme (i.e., forward-backward Euler scheme) as where the term σ S K j is treated implicitly to overcome the stiffness induced by a small ε. We then perform the QR decomposition of K n+1 j to obtain the updated basis functions X n+1 i and the matrix 2. L-step: Solve (3.13) for a full time step ∆t, update from (X n+1 ij ) using ρ n . Specifically, given L n i = r j=1 S (1) ij V n j , we discretize (3.13) (similar to (3.12)) as follows (3.17) We then perform the QR decomposition of L n+1 i to obtain the updated basis V n+1 j and matrix S (2) ij : (3.18) 3. S-step: Solve (3.14) for a full time step ∆t, update from (X n+1 ij , we discretize (3.14) (similar to (3.12)) as follows kj . 4. ρ-step: Solve (3.3) for a full time step ∆t, update from ρ n to ρ n+1 using (X n+1 i , V n+1 j , S n+1 ij ). Specifically, given ρ n , we discretize (3.3) as For clarity, we will refer to the above scheme as the K-L-S-ρ scheme in the following.

AP property of the first order scheme
In this subsection, we analyze the AP property of the first order scheme introduced in the previous section. Our conclusion is summarized in the following proposition.
Remark 3.2. If, for a given initial value (X 0 i , S 0 ij , V 0 j ), one of the conditions ξ, η, γ ∈ span({V 0 j } r j=1 ) is not satisfied, we can simply add them to the approximation space. For example, if ξ ∈ span({V 0 j } r j=1 ), we considerX where h is an arbitrary function. We then orthogonalizeX 0 andṼ 0 (e.g. using the Gram-Schmidt process) and use the result as the initial value in our algorithm. This increases the rank to at most r + 3.
Proof. In the K-step, let ε → 0, we have from (3.15): Without loss of generality, we assume that the three components of ∇xρ n σ S : ∂xρ n σ S , ∂yρ n σ S and ∂z ρ n σ S are linearly independent 2 . Then after the QR decomposition of K n+1 j , the span of the new basis {X n+1 i } 1≤i≤3 would be the same as span{ ∂xρ n σ S , ∂yρ n σ S , ∂z ρ n σ S }. In other words, we can write In the L-step, let ε → 0, we have from (3.17): Since the matrix A := ( X n+1 i σ S X n+1 k x ) 1≤i≤r,1≤k≤r is symmetric positive definite (since σ S > 0), hence invertible (whose inverse, say, is matrix B = (b ki ) 1≤k≤r,1≤i≤r ), we have After the QR decomposition of L n+1 k , we can write (by a similar argument as above) We may write (3.26) as AS n+1 = C. Since the matrix A is invertible, we know that the matrix S n+1 is unique. We next claim that the S n+1 defined as 2 If they are linearly dependent, say, span{ ∂xρ n σ S , ∂yρ n σ S , ∂z ρ n σ S } = span{ ∂xρ n σ S }, then one just needs to replace the second and third components of X 0 by X n+1 2 and X n+1 3 and the same analysis carries over. satisfies (3.26), where the middle matrix is of size r × r, with −I 3×3 in the first 3 × 3 block and zero elsewhere. Indeed, using (3.22) and (3.25) we have On the other hand, substituting (3.28) into (3.20) gives which is the forward Euler scheme for the limiting diffusion equation (2.6).

Some other first order schemes and their AP property
From the operator splitting point of view, the previously introduced K-L-S-ρ scheme is certainly not the only first order scheme. In fact, one can switch the order of K, L, and S steps arbitrarily and still obtains a first order scheme. For example, the L-K-S-ρ scheme is also first order and preserves the same asymptotic limit as the K-L-S-ρ scheme (since the proof of Proposition 3.1 still holds if one switches the K and L steps). Nonetheless, for some other first order schemes, such as L-S-K-ρ, S-L-K-ρ, K-S-L-ρ, and S-K-L-ρ schemes, their AP property needs to be examined individually. Fortunately, as we will show in the following, by slightly different arguments these schemes all have the same asymptotic limit as the K-L-S-ρ scheme.
After the first two substeps (L-S or S-L), the span of the updated basis {V n+1 j } 1≤j≤r will contain v. After the substep K, one has K n+1 Substituting g n+1 into the last ρ step recovers (3.31).
After the first two substeps (K-S or S-K), one has where D 1 is an invertible r × r matrix. After the substep L, one has satisfies (3.34). Indeed, for such L k , one has which, upon projection onto the space spanned by {X n+1 i } 1≤i≤r , yields (3.34). On the other hand, substituting g n+1 into the last ρ step recovers (3.31).
Remark 3.3. The discussion in this subsection implies that one has the flexibility to choose the updating order of K, L and S, while still maintaining the AP property. This flexibility is crucial in designing second order schemes, where one needs to properly compose these steps to achieve high order as well as preserve the asymptotic limit.

A second order in time scheme and its AP property
We now extend the first order scheme to second order. Due to the operator splitting necessary in the low rank method, a straightforward application of the IMEX-RK scheme as used in [2,14] does not work (there a coupled system for ρ and g is solved simultaneously; in the present work ρ has to be "frozen" while updating g). In the following, we propose a scheme that maintains second order in both kinetic and diffusive regimes. It is a proper combination of the almost symmetric Strang splitting [10,11] and the IMEX-RK scheme.
Suppose at time step t n , we have (X n i , V n j , S n ij , ρ n ). To obtain the solution (X n+1 i , V n+1 j , S n+1 ij , ρ n+1 ) at t n+1 , we proceed as follows: 1. ρ-step: Solve (3.3) for a half time step ∆t/2, update from ρ n to ρ n+ 1 2 using (X n i , V n j , S n ij ).
2. K-step: Solve (3.12) for a half time step ∆t/2, update from (X n i , V n j , S n ij ) to (X ij ) using ρ n+ 1 2 .
More specifically, in step 1, we use the forward Euler scheme to discretize (3.3): In steps 2-7, we use a second order IMEX-RK scheme to discretize the system for K, L or S. Let us take step 2 for example, (3.38) whereÃ = (ã pq ),ã pq = 0 for q ≥ p and A = (a pq ), a pq = 0 for q > p are s × s matrices. Along with w = (w 1 , . . . ,w s ) T , w = (w 1 , . . . , w s ) T , they can be represented by a double Butcher tableau: Here we employ the ARS(2,2,2) scheme whose double tableau is given by Finally, in step 8, we use the midpoint scheme to discretize (3.3): Let us analyze the AP property of the above second order scheme. First, steps 2-4 (K-L-S) are (almost) the same as steps 1-3 in the first order K-L-S-ρ scheme (as discussed in Section 3.2), hence as ε → 0, one has Furthermore, steps 5-6 (S-L-K) are (almost) the same as steps 1-3 in the first order S-L-K-ρ scheme (as discussed in Section 3.3), hence as ε → 0, one has which is a second-order explicit RK scheme for the limiting diffusion equation (2.6). Therefore, the scheme is AP.
Remark 3.4. There are many other choices to construct the second order scheme by altering the order of K, L and S, as long as the steps 2-4 are symmetric with respect to steps 5-7. Note that the AP property is always guaranteed due to the flexibility in the first order scheme.

Fully discrete scheme
It remains for us to specify the discretization in the physical space and velocity space. This is the purpose of this section.

Velocity discretization
For the velocity space S 2 , we adopt the discrete velocity method 3 . The velocity points {v i } i=1,...,Nv and weights {w i } i=1,...,Nv are chosen according to the Lebedev quadrature on S 2 . Then all the integrals of the form F (v) v are approximated as (3.46)

Spatial discretization
For the physical space Ω x , we assume the third dimension is homogeneous and the domain is rectangular so that we consider x = (x, y) ∈ [a, b] × [c, d]. For simplicity, we assume periodic boundary condition.
In the following, we describe a second order finite difference method in space. We use simplified notations such as ρ k,l = ρ(x k , y l ), (K i ) k+ 1 2 ,l = K i (x k+ 1 2 , y l ) to denote numerical solutions evaluated at the corresponding grid points. We use the first order K-L-S-ρ scheme in time. The discussion for other time discretization methods is similar. Figure 1: The staggered grids. ρ is located at the red dots; g (hence {K i , X i } i=1,...,r ) is located at the blue diamonds.
• K-step Note that the system (3.15), in matrix form, can be written as where K n = [K n 1 , K n 2 , . . . , K n r ] T and 48) It is clear that the matrices V 1 and V 2 are not necessarily symmetric hence the system might not be hyperbolic. Therefore, to get a reasonable spatial discretization for (3.15), we propose to discretize the original equation (3.4) and then project the resulting scheme.
Specifically, we first discretize (3.4) as where ξ + = max(0, ξ), ξ − = min(0, ξ). A second order upwind operator is applied to the spatial derivatives of g and a central difference operator is applied to the spatial derivatives of ρ. More precisely, we use and Derivatives in y are defined similarly.
We then project the equation (3.49) onto the space spanned by {V j } 1≤j≤r , which yields (3.52) Here the scheme is given at the grid points (x k+ 1 2 , y l ). The scheme at the grid points (x k , y l+ 1 2 ) is similar.

• L-step and S-step
One can add spatial discretization to (3.17) and (3.19) directly. First of all, we approximate the inner product x by a midpoint rule: Then we approximate the spatial derivatives of ρ and X i at (x k+ 1 2 , y l ) and (x k , y l+ 1 2 ) by Derivatives in y are treated similarly.

AP property of the fully discrete scheme
Similar to the semi-discrete case, in the limit ε → 0, the K-L-S steps yield g n+1 which, when substituting into (3.56), give This is an explicit standard 5-point finite difference scheme applied to the limiting diffusion equation (2.6) at grid points (x k , y l ). The limiting scheme at grid points (x k+ 1 2 , y l+ 1 2 ) can be considered similarly. Therefore, the fully discrete scheme is also AP.

A Fourier analysis of the low-rank structure of the solution
In this section, we analyze the behavior of the solution to the linear transport equation by performing a simple Fourier analysis. Our focus is in the kinetic regime because the rank is already proved to be small in the diffusive regime.
For simplicity, we consider the 1D slab geometry x ∈ [0, 2π] with periodic boundary condition, and v ∈ [−1, 1] (so v = 1 −1 · dv). Also we assume σ A = G = 0. Then the macro-micro system of the linear transport equation reads: (4.1) Projecting the above system onto the Fourier space of x yields whereρ k (t),ĝ k (t, v) andσ k are the Fourier coefficients of ρ, g and σ S , respectively. For a constant σ S we haveσ 0 = σ S ,σ k = 0, k = 0, (4.3) and the system (4.2) reduces to i.e., ρ(0, x) and g(0, x, v) are band-limited, then the latter solutions will remain in the same frequency range. In this case the solution is clearly low-rank. This analysis is similar to the analysis conducted in [12], where it was shown that for the linearized Vlasov-Maxwell equation the solution remains low rank if it is initially in a form similar to (4.5). However, the present situation is different in the sense that if we have a non-constant σ S , as is commonly the case in practice, then even if the initial value is in that form additional Fourier modes are excited gradually with time. This is as far as we can go with such an argument.
However, it should not be taken to imply that performing a low-rank approximation is necessarily futile in such a situation. In fact, the dynamical low-rank integrator makes no assumptions that the space dependence has to take the form of a finite number of Fourier modes; this is purely an artifact of the analysis done here. Hence, just because we have an infinite number of Fourier modes does not necessarily imply the solution can not be captured by a low-rank scheme. In fact, from the numerical tests in the next section, we can see that the rank of the solution in the kinetic regime when σ S is not constant can be rather intricate, but that often relatively small ranks are sufficient to obtain an accurate approximation to the dynamics of interest.

Numerical results
In this section, we present several numerical examples to demonstrate the accuracy and efficiency of the proposed low rank method. In all examples, we consider a two-dimensional square domain in physical space, i.e. x = (x, y) ∈ [a, b] 2 and periodic boundary conditions. Note that in some of the examples (e.g., the line source problem), one has to choose a large number of points in the angular direction to obtain a reasonable solution (for both the full tensor and low rank methods). This is the well-known drawback of the discrete velocity or collocation method. If a Galerkin rather than collocation approach is adopted, one could potentially use fewer discretization points (or bases). As the focus of the paper is on the low rank method, we leave the comparison of different angular discretizations to a future study.

(5.2)
Let the scattering and absorption coefficients be σ S = 1, σ A = 0, then the source term G is given by We use this source term and the initial condition ρ(t = 0, x, y) and g(t = 0, x, y, ξ, η, γ) as input for our low rank method and compute the solution up to a certain time. Note that the source term here depends also on time and velocity, hence the scheme needs to be modified accordingly to take into account this dependency. We omit the details. We consider both the first order scheme in Section 3.1 and the second order scheme in Section 3.4, coupled with the second order spatial discretization described in Section 3.5.2. We always take N v = 590 Lebedev quadrature points on the sphere S 2 [1]. Since we know a priori the rank of the exact solution g is 1, we fix r = 5 in the low rank method which is certainly sufficient to obtain an accurate solution.
We vary the spatial size ∆x and the value of ε, and evaluate the error at t = 0.1 as Since the proposed schemes are AP, we expect them to be stable under a hyperbolic CFL condition when ε ∼ O(1) and a parabolic CFL condition when ε ≪ 1. Specifically, we consider three kinds of CFL conditions: mixed CFL condition ∆t ∼ c 1 ∆x 2 + c 2 ε∆x, hyperbolic CFL condition ∆t ∼ ∆x, and parabolic CFL condition ∆t ∼ ∆x 2 .
The results of the first order (in time) scheme are shown in Figure 2. Under the mixed CFL condition, we expect to see first order convergence in the kinetic regime (ε ∼ O(1)) and second order in the diffusive regime (ε ≪ 1), which is clearly observed in Figure 2 (left). Under the parabolic CFL condition, we always expect second order convergence, which is also clear in Figure 2 (right). For the second order (in time) scheme, we don't expect order higher than two in the diffusive regime since ∆t ∼ ∆x 2 and the error behaves as O(∆t 2 + ∆x 2 ) = O(∆x 4 + ∆x 2 ). Hence we only test its performance in the kinetic regime (ε ∼ O(1)) under the hyperbolic CFL condition, where ∆t ∼ ∆x and the error is O(∆t 2 + ∆x 2 ) = O(∆x 2 ). The result is shown in Figure 3, where we can clearly see the uniform second order accuracy of the scheme (in contrast to the first order scheme).

Test with Gaussian initial value
In this test case, we consider a smooth Gaussian initial condition: with zero absorption coefficient and source term σ A = G = 0.

Constant scattering coefficient σ S
We first consider σ S ≡ 1 and focus on the AP property of the proposed method. Therefore, we set ε = 10 −6 and compare our first order low rank method with the reference solution obtained by integrating (3.58), which solves the limiting diffusion equation directly. In the low rank method, we use N x = N y = 128, N v = 590 Lebedev quadrature points on S 2 , and time step ∆t = 0.1∆x 2 + 0.1ε∆x, and fix the rank as r = 5. In solving the diffusion equation, we use N x = N y = 512 and time step ∆t = 0.75∆x 2 .
The solutions at t = 0.1 are shown in Figure 4, where they match very well. As the theory predicts, in the limiting diffusive regime, the solution g should be become rank-2. To confirm this, we track the singular values of the matrix S, see Figure 5. Clearly, the effective rank is 2 (two singular values are above the threshold of 10 −5 , which is on the order of the spatial error).

Variable scattering coefficient σ S
We then set ε = 0.01 (an intermediate regime) and consider a spatially dependent scattering coefficient whose profile is shown in Figure 6. This is a challenging test as σ S (x,y) ε varies in a large range [0.1, 100]. Our aim here is to investigate the rank dependence of the low rank method and its performance compared with the full tensor method. Specifically, we compare the first order low rank method with the first order IMEX method that solves the macro-micro decomposition of the linear transport equation directly [21] (referred to as the full tensor method in the following). We use the same spatial mesh, same CFL condition ∆t = 0.1 min(σ S )∆x 2 + 0.1ε∆x, and same N v = 2702 Lebedev quadrature points on S 2 for both methods. In the low rank method, we choose different ranks from 20 to 120.
The comparison of the low rank solution and the full tensor solution on a 256 × 256 mesh at different times is shown in Figure 7 (top). We can see that the low rank solution matches well with the full tensor solution except for rank r = 20. To quantitatively understand the rank dependence, we compute the difference of two solutions on the same mesh as follows and track how this evolves in time under certain fixed ranks r ranging from 20 to 120. The results are shown in Figure 7 (bottom). The common trend is that once the rank is increased to a certain level, the difference saturates. This is because then the spatial error dominants. Also it is clear that the rank of the solution in this problem increases gradually with time.
In addition, we record the computational time needed to compute the solution to t = 0.012 for both methods on an i7-8700k @3.70 GHz CPU in Figure 8. The speedup of the low rank method is significant, especially for a large number of spatial points N x .

Two-material test
The two-material test models a domain with different materials with discontinuities in material cross sections and source term. It is a slight modification of the lattice benchmark problem for linear transport equation. Here we choose the computational domain as [0, 5] 2 with the absorption coefficient σ A and scattering coefficient σ S given as in Figure 9. The source term is given by  Figure 9: Section 5.3: two-material test. Profiles of absorption coefficient σ A (left) and scattering coefficient σ S (right). Each square block in the computational domain is a 0.5 × 0.5 square. In the left figure, yellow square blocks represent that σ A = 10 and for the rest blue region σ A = 0; in the right figure, blue square blocks represent that σ S = 0 and for the rest yellow region σ S = 1.
We set ε = 1 and compare the first order low rank method with the first order full tensor method. For both methods, we choose N x = N y = 250, N v = 2702 Lebedev quadrature points on S 2 , and same mixed CFL condition ∆t = 0.1 min(σ S )∆x 2 + 0.1ε∆x. The initial condition is given by We test different ranks from 40 to 300 in the low rank method and compare it with the full tensor solution.
The error and computational time are reported in Figure 10. It is clear that at around rank r = 150, the spatial error dominates and increasing the rank further will have no gain in solution accuracy. Moreover, at r = 150, the efficiency of the low rank method is clearly better than the full tensor method. We then fix r = 150 and plot both the low rank solution and full tensor solution at t = 1.7 in Figure 11, where a good match is obtained. In addition, we consider another scenario with ε = 0.1. The same parameters are used as in the case of ε = 1, except we set the rank r = 100 in the low rank method (because we expect the rank of the solution to decrease as ε decreases). The solutions of the low rank method and full tensor method at time t = 0.6 are shown in Figure 12, where we again observe good agreement. An optimal (and possibly smaller) rank can be determined similarly as in Figure 10, we omit the result.

Line source test
We finally consider the line source test which is another important benchmark test for the linear transport equation. Here we approximate the initial delta function via (5.5) with a much smaller ς 2 = 4 × 10 −4 . σ S = 1 and σ A = G = 0. We set ε = 1 and compare the first order low rank method with the full tensor method. For both methods, we choose the computational domain as [−1.5, 1.5] 2 with N x = N y = 150, N v = 5810 Lebedev quadrature points on S 2 , and the same mixed CFL condition ∆t = 0.025∆x 2 + 0.025ε∆x. We fix the rank as r = 600 in the low rank method. The density profiles of both methods at time t = 0.7 are shown in Figure 13. We can see that the solutions match well.
We would like to mention that this is a difficult problem compared to the cases considered previously. Many more points are need on the sphere to get a reasonable solution. Nevertheless, there are still oscillations in the solution (for both the full tensor and the low rank method). This is a well-known  artifact in the S N method. In addition, we found that a higher rank and a more stringent CFL condition is needed in the low rank method. We believe part of the reason are the numerical oscillations, which can be tempered by applying a proper filter or using a positivity-preserving scheme. We refer to [20], and references therein, for more details.

Conclusion
We have introduced a dynamical low-rank method for the multi-scale multi-dimensional linear transport equation. The method is based on a macro-micro decomposition of the equation and uses the low rank approximation only for the micro part of the solution. The key feature of the proposed scheme is that it is explicitly implementable, asymptotic-preserving in the diffusion limit, and maintains second order in both kinetic and diffusive regimes. A series of numerical examples in 2D including some wellknown benchmark tests have been performed to validate the accuracy, efficiency, rank dependence, and AP property of the proposed method. Some interesting ongoing and future work includes adaptive rank selection and the theoretical investigation of rank dependence of the solution in the kinetic regime.