An asymptotic preserving scheme for kinetic models with singular limit

We propose a new class of asymptotic preserving schemes to solve kinetic equations with mono-kinetic singular limit. The main idea to deal with the singularity is to transform the equations by appropriate scalings in velocity. In particular, we study two biologically related kinetic systems. We derive the scaling factors and prove that the rescaled solution does not have a singular limit, under appropriate spatial non-oscillatory assumptions, which can be verified numerically by a newly developed asymptotic preserving scheme. We set up a few numerical experiments to demonstrate the accuracy, stability, efficiency and asymptotic preserving property of the schemes.


Introduction
We consider the following type of kinetic equations, (1.1) Here, f ε = f ε (t, x, v) is the probability density function at time t ≥ 0, of space variable x = (x 1 , . . . , x d ) T ∈ Ω and velocity v = (v 1 , . . . , v d ) T ∈ R d , with spatial domain Ω = R d or T d . Q is the interaction operator, which can be nonlinear in f ε and nonlocal in x and v.
The main property of the interaction operator Q of our concern is that it has a mono-kinetic equilibrium, namely where δ is the Dirac delta distribution, and ρ(t, x) and u(t, x) = (u 1 (t, x), . . . , u d (t, x)) T are macroscopic density and velocity, respectively, satisfying Under this setup, one can formally let ε → 0 in (1.1) and obtain an asymptotic solution, lim ε→0 f ε (t, x, v) = ρ(t, x)δ v=u(t,x) , (1.2) which is the equilibrium of Q, and is singular in v.
In this paper, we shall focus on the following two interaction operators, both of which have interesting biological applications. The first model is called aggregation system, where the interaction operator is defined as The operator consists two parts. The first term describes pairwise attraction-repulsion interactions, where K is the interaction potential. A natural biological assumption is that the strength of the interaction depends on the distance between two agents: attraction in large distance and repulsion in short distance. Hence, K = K(r) is radial, and it is decreasing when r is small and increasing when r is large. The second term represents relaxation in velocity. This term is less biologically motivated, but plays a crucial role in deriving an interesting asymptotic limit [2]. In fact, the mono-kinetic asymptotic solution (1.2) is rigorously derived in [16], see also [10]. Furthermore, (ρ, u) satisfy ∂ t ρ + ∇ x · (ρu) = 0, u(t, x) = − Ω ∇ x K(x − y)ρ(t, y)dy. (1.4) This limiting system is realized as the aggregation equation which appears in various contexts related to biological aggregation models. The equation has been intensively studied in the recent decade, and we refer to [19,25] and references therein. The second model is called 3-zone system, where (1.5) The artificial relaxation term in (1.3) is replaced by an alignment term, which models pairwise interactions in the middle range. The alignment force, proposed by Cucker and Smale in [7], describes the so called flocking phenomenon that agents align their velocities to the neighbors. Here, φ(x) is the influence function which represents the strength of alignment between two agents. It naturally depends on the distance between the agents, and decreases when the distance becomes larger. We also assume that φ is bounded and Lipschitz. Without loss of generality, we take The kinetic representation of Cucker-Smale model is derived in [15], analyzed in [4,24], and numerically studied in [21,24]. We refer readers to [8,20] for discussions on Cucker-Smale dynamics with singular influence function. The interaction operator (1.5) combines long-range attraction, short-range repulsion and mid-range alignment. Such 3-zone interaction framework is proposed in [22]. It has been very successful in biological and ecological modeling, and it is widely used in computer animations. As ε → 0, the asymptotic limit of (1.1) with interaction (1.5) is rigorously derived in [9], where mono-kinetic asymptotes (1.2) is justified, with (ρ, u) satisfying (1.6) The wellposedness theory of the limiting system (1.6) is also established in [9], with the additional equality on momentum conservation. The system serves as a more biologically relevant substitute to the aggregation equation (1.4). The goal of this paper is to design a universal numerical scheme for (1.1) that solves the equation in both the kinetic regime when ε = O(1), and the fluid regime when ε → 0. This type of numerical schemes is called asymptotic preserving (AP) and was originally introduced in [17]. The commutative diagram on the left hand side of Figure 1 illustrates the AP property. A scheme f h ε that approximates the solution f ε with discretization parameters h is AP if its stability requirement on h is independent of ε, and if its limit f h when ε tends to zero consistently serves as an approximation of the limiting solution f . Therefore, the scheme can be automatically applied to the limiting equation simply by setting ε → 0. AP schemes have been very successful in solving kinetic equations with different types of hydrodynamic limits, see, e.g., [18] for a recent review of AP schemes. In the conventional kinetic equations and corresponding AP schemes, the limiting profile is usually given by a smooth Maxwellian distribution. Hence one can use fixed grid points in velocity discretization with a cutoff. The study of kinetic equations with non-Maxwellian equilibrium has received attentions recently. AP schemes have been designed for the kinetic equations with heavy-tail equilibrium [5,6,26,27].
The equilibrium of the alignment operator Q for our system (1.1), on the contrary, is given by a δ-distribution in velocity space. As ε becomes small, the solution f ε becomes more and more singular. This addresses a major challenge in designing AP schemes for (1.1) as its direct discretization can not achieve high accuracy and stability for small ε due to the fact that the limit solution is singular.
To overcome the difficulty, we apply a family of transformations T ε to the original system (1.1). As illustrated in Figure 1, f ε is mapped to a new function The aim is to find appropriate transformations so that the limiting solution g = lim ε→0 g ε is not singular, and thus an AP scheme for g ε can be designed without worrying about the singularity. Figure 1. AP scheme under transformation Since the singularity of f is a δ-distribution in velocity, a natural choice of the transformation T ε is a scaling in velocity. Velocity scaling methods have been used in various kinetic systems with singular time asymptotic limits, see, e.g., [1,11,12,21]. The heart of the matter in these methods is to find an appropriate scaling factor ω ε to ensure that the rescaled function g ε is not singular. In [1,11,12], the choice of ω ε is based on the self-similar behavior of the spatial homogeneous equations, in which the transport part in (1.1) is omitted. The scaling factor used in these works is proven to be optimal only for homogeneous systems with self-similar initial configurations. In [21], a new scaling factor is introduced for kinetic flocking systems and it is shown to be exact for spatially homogeneous systems with all smooth initial conditions.
In this paper, we present an AP scheme for (1.1) based on the velocity scaling method, where the transformation is given by (2.2), with a scaling factor similar to the one proposed in [21]. We study the asymptotic behavior of the rescaled function g ε , and provide sufficient conditions to ensure g ε is non-singular uniformly in ε. The result indicates that our choice of scaling factor captures the right scaling. Moreover, it implies that our numerical scheme is indeed asymptotic preserving.
The rest of the paper is organized as follows. In Section 2, we describe the velocity scaling method, and show that with our choice of the scaling factor, the scaling is exact for the spatial "homogeneous" aggregation and 3-zone systems. In Section 3, we discuss the asymptotic behavior of the full system (1.1), and prove that the rescaled profile is not singular under appropriate non-oscillatory conditions. In Section 4, we design AP schemes for the systems after velocity scaling, and discuss the AP property of (1.1). Finally, in Section 5, we provide numerical experiments to illustrate the performance of the new scheme.

Velocity scaling method
In this section, we present the velocity scaling method for (1.1). We shall follow the storyline of [21] to derive a rescaled system.

2.1.
Exact rescaling on spatial "homogeneous" system. As the main driving force of the system towards singularity is the interaction operator Q, we first consider the spatial "homogeneous" system omitting the free transport part. We rescale the velocity variable by and the transformation T ε is defined as Here, ω ε = ω ε (t, x) is a scaling factor, and u ε = u ε (t, x) is the macroscopic velocity defined as The rescaled function g ε has the following properties. First, the macroscopic density of g ε is the same as the macroscopic density of f ε , Second, with the shift by u ε , the first moment of g ε in ξ is zero for all x, namely the profile g ε is nicely centered in ξ, An appropriate choice of scaling factor ω ε should produce a non-singular rescaled function g ε , which neither concentrates nor spreads out as ε approaches zero, namely max ξ |g ε (t, x, ξ)| ≤ G, and supp ξ g ε (t, x, ξ) ⊂ B R (0), ∀x ∈ Ω, ∀t ∈ [0, T ], (2.6) where G and R are finite and independent of ε. Here, B R (0) denotes for a ball centered at origin and has radius R in (R d , | · | ∞ ).
To choose an appropriate ω ε , we represent the dynamics of f ε by the triple (g ε , u ε , ω ε ). The term ∂ t f ε can then be expressed by and the interaction kernel Q is expressed as Here, A ε = A ε (t, x) and B ε = B ε (t, x) differ for different models. For aggregation system (1.4), For 3-zone system (1.6), (2.10) Combining (2.7) and (2.8), we obtain an evolution equation for g ε To describe the dynamics of ρ ε and u ε , we take zeroth and first moments of f ε in (2.1): which in turn implies that Therefore, taking into account (2.12) and defining the scaling factor ω ε in (2.11) as leads to ∂ t g ε = 0, and thus g ε remains unchanged in all time. In this case, we say that the rescaling is exact with factor ω ε .
Since the initial profile f 0 (x) does not depend on ε, it is easy to check that the initial triple (g 0 , u 0 , ω 0 ) is also independent to ε. Hence, the solution g(t, x, ξ) = g 0 (x, ξ) remains the same while ε varies and condition (2.6) is clearly satisfied as long as g 0 satisfies (2.6). Moreover, as the scaling is exact, we can easily reconstruct f ε as follows: In the previous works [11,12], the scaling factor was chosen to be ω ε (t, x) = T ε (t, x) −1/2 , where T ε was the temperature of the system, that is, It was also shown that such scaling factor was exact for self-similar initial data. A new scaling factor proposed in [21] takes advantage of the structure of the interaction operator and it is exact for all initial data.

2.2.
Rescaling on the full system with free transport. We now apply the scaling argument to the full system (1.1). The presence of free transport destroys the self-similar structure of the spatial homogeneous system (2.1) and therefore it is in general impossible to find an exact scaling. We thus extend the idea of the new scaling factor to find a non-singuar rescaled function, in the sense of (2.6). The free transport term can be expressed in terms of (g ε , u ε , ω ε ) as follows, Adding this new contribution to (2.11), yields (2.14) Obtaining an exact scaling in this case would require finding a scaling factor ω ε that satisfies Since ω ε = ω ε (t, x) is independent on the velocity variable ξ, such ω ε does not exist. Instead, we take ω ε which satisfies We again set ω 0 (x) ≡ 1, namely we do not perform scaling at t = 0.
By taking the first two moments moments of (1.1), we deduce the dynamics of macroscopic density ρ ε and velocity u ε : where P ε is the pressure tensor defined as Note that equation (2.17) can be rewritten in the following non-conservative form: 19) and the two forms are equivalent in the non-vacuum region where ρ ε (x) > 0. Taking into account (2.15), (2.16) and (2.17), equation (2.14) can be rewritten as 20) or in the following equivalent conservative form: Unlike the spatial "homogeneous" system (2.1), the rescaled function g ε does change in time now, and it varies with different ε. To validate our choice of scaling factor for the full system, it is important to check that g ε satisfies (2.6) uniformly in ε, particularly when ε approaches zero.

Asymptotic behavior
This section is devoted to studying the asymptotic behavior of equations (2.15)-(2.20), as ε → 0. The goal is to understand whether g ε is non-singular under the proposed rescaling when ε is small. The result also supports the AP property of the numerical scheme that will be discussed in Section 4 below.
We denote and R ε (t) be the smallest number such that We also recall that g ε is non-singular if condition (2.6) is satisfied and hence we shall show that G ε (t) and R ε (t) are bounded independent of ε, for all t ∈ [0, T ], under appropriate assumptions.
3.1. Non-oscillatory assumptions. We start our discussion with two assumptions on the solution triple (g ε , u ε , ω ε ). The first one is a spatially non-oscillatory assumption on the rescaled function g ε , 3) implies non-oscillatory bounds on macroscopic quantities. Indeed, for density ρ ε , we have For pressure P ε , we have the following estimate and condition (3.3) implies If condition (3.3) is violated, then g ε becomes more oscillatory when ε gets smaller, in which case one can not expect to design AP numerical scheme for g ε . The second assumption is the Lipchitz apriori bound on the macroscopic velocity It should be observed that taking ε → 0 in (2.16) and (2.19), one can formally obtain the limiting system, for which condition (3.7) is also satisfied. The argument has been rigorously proved in [16] for the aggregation system (1.4) and in [9] for the 3-zone system (1.6). For both systems, the limiting velocity u is Lipschitz globally in time, under suitable regularity assumptions on kernels K and φ, and thus satisfies (3.7). The regularity for the limiting system does not imply, however, that (3.7) holds uniformly in ε. In fact, the convergence of f ε to ρδ v=u is only weak- * in measure. This does not rule out the possibility of oscillation in x as ε → 0.
In [3,23], it has also been proven that condition (3.7) is satisfied when the system (2.16), (2.19) is considered in the pressureless regime and does not depend on g ε , i.e., and subject to subcritical initial data. Moreover, the subcritical region becomes larger when ε gets smaller. Therefore, (3.7) is satisfied uniformly for ε ∈ [0, ε 0 ] if the initial profile u 0 lies in the subcritical region of the system (3.8) with ε = ε 0 . The result for pressureless system (3.8) can be easily extended to the general dynamics (2.19) when the pressure term is Lipschitz bounded uniformly in ε. Indeed, from (3.5) and (3.6), we know P ε /ρ ε and ∇ x P ε /ρ ε are uniformly bounded by R 2 ε . This together with the estimate on ω ε (see Section 3.2) implies boundedness of the pressure term in (2.19). The Lipschitz bound can also be obtained by the additional non-oscillatory assumption |∇ ⊗2 x g ε (t, x, ξ)| g ε (t, x, ξ). We omit the proof and redirect the reader to [3,23] for relavant discussions.
We have thus argued that condition (3.7) holds under appropriate setup if g ε is non-oscillatory in x. Therefore, both (3.3) and (3.7) are considered as spatially nonoscillatory assumptions on g ε . Given these two assumptions, we are going to prove that g ε is non-singular, uniformly in ε.

3.2.
The scaling factor. In this section, we use the evolution equation (2.15) to estimate both the scaling factor ω ε and its gradient ∇ x ω ε . To this end, we begin with the following proposition.
Proposition 3.1. Assume A ε is bounded below by a positive constant c that is independent of ε, then ω ε (t, ·) L ∞ (Ω) tends to 0 as ε → 0.
Proof. Consider a flow map X ε (t, x) such that (3.9) Along each characteristic path, we have which in turns yields Collecting all paths, we obtain which vanishes as ε → 0.
Since lim ε→0 f ε is singular, a correct rescaling has to have a factor ω ε vanishes as ε → 0. This is true for our choice of ω ε .
Remark 3.2. Under appropriate settings, the lower bound assumption on A ε is valid for both the aggregation system (2.9) and 3-zone system (2.10). Indeed, for the aggregation system, A ε ≡ 1, while for the 3-zone system, A ε can be estimated by provided φ is lower bounded by φ min > 0. Note that ρ ε (t, ·) L 1 (Ω) = 1 due to mass conservation, and the fact that f ε is a probability distribution. The assumption on φ can be further relaxed (see e.g. [23]). We omit the details.
Proof. We start with the estimate on ∇ x A ε and obtain from (2.10): where the first inequality is due to non-oscillatory condition (3.4). Here, we recall our assumption that φ is bounded and φ L ∞ = 1. We now estimate ∇ x ω ε by applying operator ∇ x to equation (2.15). Once again we consider the flow map (3.9) and obtain the following equation along each characteristic path: This implies As ω 0 (x) ≡ 1, we obtain M(0, x) ≡ 0. Given any t ∈ [0, T ], we combine all paths and use (3.7) (3.10), to obtain which also vanishes as ε → 0.
Remark 3.3. It follows from Proposition 3.3 that if the rescaled function g ε does not spread out, it will not concentrate either. In particular, for the aggregation system (2.9) there is no concentration regardless of the size of the support of g ε since ∇ x ω ε ≡ 0 in this case.
We are left to prove that g ε does not spread out. The growth of the support of g ε is equivalent to the spread of the characteristic paths in (3.11), whose dynamics implies the following estimate on R ε (t): From the non-oscillatory bounds (3.5)-(3.7), we obtain The estimate has the form d dt where a ε can be determined by Propositions 3.1 and 3.2. The last inequality (3.12) allows us to prove that under appropriate assumptions on ε and R(0), the function R ε (t) in (3.2) is bounded globally in time.
Proof. We first consider aggregation system (2.9). From the estimate (3.10) and the fact that ∇ x ω ε = 0, we have We denote S ε (t) := R ε (t) exp(− c ε t) and obtain from (3.12) the dynamics of S ε : We now take ε 0 = c 2C 2 and R 0 = C 2 C 1 and observe that since c − C 2 ε > 0, S ε has an invariant region 0, c − C 2 ε C 1 ε and the following inequality holds: Therefore, S ε (t) ≤ c − C 2 ε C 1 ε for all t ≥ 0 and we conclude with the bound Next, we turn our attention to the 3-zone system (2.10). By propositions 3.1 and 3.2, we have The extra exponential growth e C 2 t in (3.13) is due to the estimate of ∇ x ω ε . It can be controlled by the exponential decay provided C 2 < c ε . In fact, we have Using the same argument as the aggregation system, we obtain the following bound: It should be pointed out that the two bounds obtained above are not uniform in ε. Uniform bounds can only be achieved up to a finite time, provided that a ε (t) is uniformly bounded.
Proof. For the aggregation system (2.9), a ε (t) is uniformly bounded by C 1 and from which is a Ricatti-type first order ODE. Therefore, there exists a finite time T = T (C 1 , C 2 , R(0)) > 0, such that R ε (t) remains finite in [0, T ]. Since T does not depend on ε, the bound is uniformly in ε.
For the 3-zone system (2.10), we use the estimate (3.13) to obtain The 1 ε term can be controlled by the exponentially decay, namely, for all ε ∈ [0, ∞), which in turns implies The right hand side in the last inequality is an increasing function in t. Therefore, we conclude that a ε (t) is bounded by , which does not depend on ε, and thus according to (3.12), R ε (t) is uniformly bounded in ε for any finite time t ∈ [0, T ].
Putting everything together, we prove that g ε is non-singular. It provides a strong support that our choice of ω ε captures the right scaling.

Asymptotic preserving schemes
Now we design a asymptotic preserving scheme to solve (1.1). To avoid the singularity limit, we use velocity scaling method and express the solution f ε by the rescaled function g ε , together with the scaling factor ω ε and macroscopic velocity u ε . We have shown in Theorem 3.6 that under our proposed rescaling, g ε is non-singular uniformly in ε. Therefore, we proceed to design AP schemes for the rescaled system, where singularity is no longer an obstacle. 4.1. AP schemes for the rescaled systems. Let us recall the dynamics of the solution triple (g ε , u ε , ω ε ) and rewrite equations (2.15), (2.17) and (2.20) in the numerical friendly conservative representations, with ρ ε (t, x) defined in (2.4) and satisfying the continuity equation (2.16).
To obtain an AP scheme for system (4.1), we introduce an increasing sequence 0 < t 0 < t 1 · · · < t n · · · of times with uniform time step ∆t = t n+1 − t n and denote by q n the value of any unknown quantity q at time t n , i.e., q n (·) ≈ q(t n , ·). The canonical first order in time explicit-implicit time discretization for (4.1) reads: where the non-stiff fluxes are treated explicitly and the stiff terms are treated implicitly. To evolve the solution in time, we first compute g n+1 ε from the first equation in (4.2), which is fully explicit as g ε is non-singular, and its dynamics does not explicitly depend on ε. Then, ρ n+1 ε is obtained from the integration of g n+1 ε in ξ coordinate. Next, we use an implicit solver to compute ρ n+1 ε u n+1 ε from the second equation. Noting that the operator ρ ε B ε is a symmetric operator on ρ ε u ε , one can simply apply a conjugate-gradient method. Finally, ω n+1 ε can be obtained easily from the third equation since A n+1 ε only depends on ρ n+1 ε and hence can be computed explicitly. One can derive a second order time discretization scheme by applying, say, a backward differentiation formula (BDF) on the time derivative, an extrapolation on the explicit terms and a fully implicit solver on the stiff terms. We omit the details here and refer the reader to [13,14,26].
A fully discrete scheme should be obtained by consistent spatial and velocity discretizations, for instance, by using a finite volume method thanks to the conservative structure of the equations; see, e.g., [21] for the references. Importantly, since g ε is non-singular, the discretizations are independent of ε.
We summarize the entire procedure of the proposed numerical approach for solving (1.1). Given initial data f 0 , we set ω 0 ≡ 1, compute u 0 by (2.3) and g 0 by performing velocity scaling transformation T ε in (2.2). Then, we evolve the dynamics (4.1) on (g ε , u ε , ω ε ) using appropriate AP scheme, for instance (4.2), until a target time t. Finally, we apply the inverse transformation T −1 ε to obtain the solution f ε at time t. Note that T −1 ε has an explicit form which is easy to implement numerically.

4.2.
Asymptotic preserving property. Now, we verify the AP property of our numerical scheme. Recall the limiting system of (1.1) as ε → 0 satisfies mono-kinetic asymptotes (1.2) f (t, x, v) = ρ(t, x)δ v=u(t,x) , with macroscopic quantities (ρ, u) satisfying The goal is to check that f n ε converges to f n as ε → 0, at the discrete level. As discussed in Section 3, the spatial non-oscillatory conditions (3.3) and (3.7) play an important role and guarantee that g ε is non-singular, in the sense of (2.6). This argument, stated in Theorem 3.6, can be extended to semi-discrete or fully discrete dynamics with appropriate choices of discretizations.
We shall consider the first order scheme (4.2) as an example. The semi-discrete version of Theorem 3.6 implies that g n ε is non-singular, namely max ξ |g n ε (x, ξ)| ≤ G, and supp ξ g n ε (x, ξ) ⊂ B R (0), ∀x ∈ Ω, (4.5) where G, R are constants independent of ε, if for all x ∈ Ω, ξ ∈ R d , where C 1 , C 2 are constants which do not depend on ε. Assuming (4.6) holds, we, first, check f n ε converges to a mono-kinetic profile as ε → 0. It is enough to show that for any given x ∈ Ω, the size of supp v f n ε (x, v) tends to 0 as ε → 0. From (4.3) and (4.5), we obtain A semi-discrete version of proposition 3.1 implies that ω n ε (x) → 0 as ε → 0, which finishes the proof. Indeed, Here, A n ε ≥ c due to Remark 3.2.
Next, we show that the macroscopic quantities (ρ n ε , u n ε ) converges to (ρ n , u n ), which solves the semi-discrete version of the limiting system (4.4): To this end, we integrate the g n ε equation in (4.2) with respect to ξ to obtain ρ n+1 ε − ρ n ε ∆t + ∇ x · (ρ n ε u n ε ) = 0. Clearly, the limiting system as ε → 0 is the first equation in (4.7).
For the second equation in (4.7), we rewrite the u n ε equation in (4.2) as follows where the right hand side is of order O(ε), thanks to the non-oscillatory condition (4.6). Taking the limit ε → 0, we obtain ρ n+1 B n+1 = 0. It should be observed, that the AP property can be also verified for full discrete schemes. Detailed discretization can be found in, e.g. [21].
Note that the discrete non-oscillatory conditions (4.6) can be monitored during numerical simulations. Practically, instead of monitoring the oscillation on g n ε (x, ξ) for all x ∈ Ω, ξ ∈ R d , we only need to keep track of the oscillation for ρ and P , namely the discrete version of conditions (3.4) and (3.6). The AP property is guaranteed to hold as long as there is no violation of the following assumptions where C is a constant which does not depend on ε.

Numerical experiments
In this section, we demonstrate the performance of the proposed schemes on a number of numerical examples. We note that the velocity scaling method in Section 2 and the resulting AP scheme in Section 4 are dimension independent. For simplicity, numerical simulations are performed on a 1-D by 1-D phase space, with periodic spatial domain Ω = T = [−π, π]. In particular, we consider the computation domain (x, ξ) ∈ [−π, π] × [−6, 6], and pick initial data such that R in (2.6) is much smaller than 6. So, the solution will vanish at the boundary. Unless otherwise specified, we always take N x = 128 and N ξ = 64 grid points in the phase space. We take ∆t = ∆x/20 to satisfy the CFL condition, where ∆x is spatial mesh size.
In this section, we focus on the 3-zone system (1.5). The aggregation system (1.3) can be solved similarly. The alignment kernel φ of the 3-zone system is given by In the Examples 1-3 below, the interaction is modeled by the Morse potential K(x) = −e −|x|/2 + e −|x| .

Example 1 -Validation of the assumptions.
The first test is to check whether the spatial non-oscillatory assumption (4.8) is valid for a typical initial value problem of (1.1). The rescaled system (4.1) is numerically solved subject to the initial data where M(ξ) = 1 √ 2π e −ξ 2 /2 . We track the time evolution of max x |∇xρε| ρε , max x |∇xPε| ρε and max x |∇ x u ε |, for different values of ε. The results shown in Figure 2 suggest that the assumption (4.8) is valid and the bounds are uniform with respect to ε.

Example 2 -Consistency test.
In this example, we verify that the solution to the rescaled system (4.1) is consistent with the original system (1.1). The original system is integratedsimulations in time by the forward Euler method, while the rescaled system is evolved by the AP scheme (4.2). The original system is very difficult to solve for small ε and long time, due to the fact that the solution f ε is approaching a singular delta function in velocity space. Hence, we taksimulationse ε = 1 and run the simulations until the final time t = 0.7 in this test. N v = 512 points are used in solving the original system (compare with N ξ = 64 for the rescaled system).
The following initial condition for the original system (1.1) is used, which is equivalent to the rescaled system (4.1) solved subject to the initial condition Time snapshots of the density ρ 1 (x) and the macroscopic velocity u 1 (x) at different time t are compared in the top of Figure 3. The solutions to different systems are almost identical, demonstrating that the rescaled system is consistent with the original system. In the bottom of Figure 3, we show the distributions f 1 (x, v) (left, solved from the original system) and g 1 (x, ξ) (right, solved from the rescaled system) at time t = 0.7. As one can see, f 1 (x, v) is getting concentrated in the velocity space, making it difficult to simulate with fixed grid points. In contrast, g 1 (x, ξ) has a finite support in the rescaled velocity space ξ and a fixed grid in ξ can be used for the simulation. 5.3. Example 3 -Asymptotic preserving test. Now we test the AP property of the scheme (4.2). More specifically, we compare the solutions of (4.2) with vanishing ε to the solution of the limiting system (1.6). We use the following initial data g 0 (x, ξ) = ρ 0 (x)M(ξ), ρ 0 (x) = 0.01 + e −20x 2 , u 0 (x) = 0, ω 0 (x) = 1, for the scheme (4.2). The limiting system (1.6) with initial condition (ρ 0 , u 0 ) is wellposed with momentum conservation condition We refer to [9] for analysis and numerical schemes for the limiting system. The comparison of the density ρ ε (x) and macroscopic velocity u ε (x) at time t = 1 is given in Figure 4. Different ε's are used for the scheme (4.2). The results clearly demonstrate that as ε vanishes, the solution obtained from (4.2) approach the solution to the limiting system, demonstrating the AP property of (4.2).  which describe two groups of agents in the same location moving to opposite directions.
The strength of interactions between agents are characterized by the value of ε. In Figure 5, we take ε = 1, hence a weak interaction is used. Time snapshots of the distribution g 1 (x, ξ), the density ρ 1 (x), the momentum ρ 1 (x)u 1 (x) and the scaling factor ω 1 (x) at different times are provided. It can be observed that the two groups continue moving toward opposite directions and eventually are separated from each other. The scaling factor ω 1 (x) decays to 0 uniformly in x. The alignment begins to dominate after a long time simulation, driving the momentum ρ 1 u 1 to zero.
In Figure 6, we plot the solution of same problem with a strong interaction by taking ε = 10 −4 . The effects of alignment and attraction/repulsion are much stronger than the free transport. The alignment plays a role in two aspects. First, it pushes ω ε to 0 immediately, describing all agents in the same location moving with the same velocity. This makes the two groups stick together. Second, after a long time, the alignment drives the momentum ρ ε u ε to zero for all x, hence forming a flocking pattern. The attraction/repulsion determines the shape of this pattern. In Figure 6, we also include the stationary solution (see e.g. [9]) of the limiting system (1.4) and note that it agrees very well with the long time profile of the aggregation system.  Figure 5. Example 4: Time snapshots of the solution to the aggregation system. From left to right: the distribution g ε (x, ξ), the density ρ ε (x), the momentum ρ ε (x)u ε (x) and the scaling factor ω ε (x). In this test ε = 1.  Figure 6. Example 4: Time snapshots of the solution to the aggregation system. From left to right: the distribution g ε (x, ξ), the density ρ ε (x), the momentum ρ ε (x)u ε (x) and the scaling factor ω ε (x). In this test ε = 10 −4 . The stationary solution ρ and ρu of the limiting system (1.6) is illustrated by red dashed lines in the last row.