AN ASYMPTOTIC PRESERVING SCHEME FOR KINETIC CHEMOTAXIS MODELS IN TWO SPACE DIMENSIONS

In this paper, we study two-dimensional multiscale chemotaxis models based on a combination of the macroscopic evolution equation for chemoattractant and microscopic models for cell evolution. The latter is governed by a Boltzmann-type kinetic equation with a local turning kernel operator which describes the velocity change of the cells. The parabolic scaling yields a non-dimensional kinetic model with a small parameter, which represents the mean free path of the cells. We propose a new asymptotic preserving numerical scheme that reflects the convergence of the studied micro-macro model to its macroscopic counterpart—the Patlak-Keller-Segel system—in the singular limit. The method is based on the operator splitting strategy and a suitable combination of the higher-order implicit and explicit time discretizations. In particular, we use the so-called even-odd decoupling and approximate the stiff terms arising in the singular limit implicitly. We prove that the resulting scheme satisfies the asymptotic preserving property. More precisely, it yields 2010 Mathematics Subject Classification. Primary: 92C17, 35Q20, 76M45; Secondary: 76P99, 65M06, 65M70.

1. Introduction.Chemotaxis is the movement of cells towards a chemical substance in a medium and it is often modeled by systems of PDEs.The classical PDE chemotaxis model is the Patlak-Keller-Segel (PKS) system [49,38,39,40], which is derived at the macroscopic level in terms of the cell density and chemoattractant concentration.In the two-dimensional (2-D) case, this model reads as where x = (x, y) ∈ Ω ⊂ R 2 are spatial variables and t is time, ρ(x, t) is the cell density and S(x, t) is the chemoattractant concentration, D and α are positive diffusion constants, χ is the chemotactic sensitivity constant, and the positive constants γ and β stand for the production and degradation rate of the chemoattractant, respectively.The constant τ determines the type of the system: It is parabolic-parabolic if τ = 1 and parabolic-elliptic for τ = 0.In recent years, several modifications of the PKS system have also been studied; see, e.g., [5,9,26,27,28,42,50] and references therein.
In order to describe the chemotaxis at the cellular (microscopic) level, a class of Boltzmann-type kinetic equations has been developed.A stochastic approach based on the velocity-jump process was introduced in [54] and was later used in the framework of kinetic chemotaxis models in [1,48,52].The velocity-jump process characterizes the movement in two phases, namely, run and tumble.During the run phase, the cells move (almost) linearly with constant speed and in the tumble phase, they reorient their motion with a new velocity and direction.The Boltzmann-type kinetic model reads as where f := f (x, t, v) is the probability density function (pdf) of cells at the position x with the velocity v = (u, v) ∈ V ⊂ R 2 at a given time t, and f := f (x, t, v ).In (3), T [S] is the turning kernel operator, which describes the velocity change from v to v at (x, t), that is, T [S] := T [S](x, t, v, v ) and T * [S] := T [S](x, t, v , v).Specifying the turning kernel T is a crucial point in the kinetic chemotaxis modeling, which will be discussed in the §2.Notice that the microscopic pdf f is related to the macroscopic cell density ρ in the following way: Applying the parabolic scaling to (3) yields the following non-dimensionalized kinetic equation (see, e.g., [7,23,47]): with the non-dimensional scaling parameter (mean-free path) ε and the new notations: The question of convergence (in the singular limit as ε → 0) of the kinetic model (5) to the PKS system (1), (2) has been extensively studied.More precisely, the global in time convergence was proven in the parabolic-elliptic case in [7].In the parabolic-parabolic case, only local convergence results were established; see [31].We also refer the reader to [47] for more results on the limiting process.
It is well-known that if the total number of cells is sufficiently large, a concentration phenomenon may occur and can be modeled by both the PKS (1), (2) and kinetic-chemotaxis ( 5), (2) systems.It should be observed that under the assumption that no-flux conditions are imposed at the boundary of the domain Ω, the total mass is conserved for both models, but the solution behavior depends on the value of M .For instance, the solution of the 2-D PKS system (1), ( 2) may develop δ-type singularities in finite time if M is larger than some critical value M c ; see, e.g., [10,20,22,24,4,46,33,16,21].Otherwise, the solution of ( 1), (2) exists globally in time.In the parabolic-elliptic case (τ = 0), the critical mass values are explicitly available, while this is not the case for the parabolic-parabolic system (τ = 1); see, e.g., [50].The kinetic-chemotaxis system (5), (2) exhibits a similar behavior, which depends, however, not only on the value of the initial mass M , but also on the choice of the specific kernel T ε ; see, e.g., [3,7].At the same time, the kinetic-chemotaxis system provides a more detailed description of the underlying cell dynamics and thus may be advantageous in a variety of applications.The goal of this paper is twofold.First, we develop an efficient and accurate asymptotic preserving (AP) numerical method for the kinetic-chemotaxis system (5), (2).Secondly, we use this AP method to experimentally study the blowup behavior of the kinetic chemotaxis model and compare it with the macroscopic PKS model for both parabolic-elliptic (τ = 0) and parabolic-parabolic (τ = 1) cases.
One of the difficulties in achieving an efficient approximation of ( 5), ( 2) is related to the fact that the studied system is stiff when 0 < ε << 1.If a naive numerical discretization is used, then one may need to take both spatial and temporal discretization parameters to be proportional to O(ε) or even O(ε 2 ) due to stability restrictions, which may become unaffordable for small ε.To overcome this difficulty, we develop an AP scheme, which yields a consistent approximation of the limiting macroscopic PKS system as ε → 0 and is stable on a coarse spatio-temporal grid with the mesh parameters being independent of ε.The concept of AP methods was introduced in [37,41] to solve kinetic transport equations in diffusive regimes.The main idea is based on a suitable operator splitting and odd-even decoupling as will be presented §3 below.This approach was afterwards generalized for a variety of kinetic models; see, e.g., [19,34,35,14,30,13,29,11,12,43] and reference therein.
The AP property of our numerical approach is achieved by implementing an operator splitting technique combined with an idea of the even-odd formulation; see, e.g., [36,6].We split the system (5), (2) into the following three parts: a non-stiff linear transport system, stiff nonlinear relaxation equations for the pdf, and the macroscopic equation (2) for the chemoattractant concentration.The linear transport system is then solved by a second-order upwind method, while the nonlinear relaxation part is solved exactly.Finally, the chemoattractant equation is treated either by fast Fourier transform (FFT) in the parabolic-elliptic case or by a semidiscrete centered-difference discretization along with an implicit-explicit (IMEX) time integration in the parabolic-parabolic case.We would like to stress that the operator splitting has to be designed carefully so that the kinetic and macroscopic parabolic/elliptic equations are coupled in the way that allows one to design an efficient numerical method; see §3.2.
The paper is organized as follows.In §2, we describe a particular choice of the turning kernel for the underlying kinetic-chemotaxis model.The numerical method is developed in §3 and its AP property is then proven in §4.Finally, in §5, the performance of the proposed AP scheme is illustrated on a number of numerical examples suggesting that our AP method may be a suitable tool to get a better insight to the blowup behavior of kinetic chemotaxis model in two-space dimension for both parabolic-elliptic and parabolic-parabolic cases.We refer the reader to a related study on blowup behavior [6], where only parabolic-elliptic one-and radially symmetric two-dimensional models have been considered.

2.
Local turning kernel.In this section, we select a specific turning kernel T ε , which models the reorientation process of the cells and will be used in the kinetic equation (5).To this end, we first write the formal asymptotic expansion (see, e.g., [7,6,23,47,15]): Here, the leading term T 0 [S] = F (v) > 0 is the bounded velocity distribution at the equilibrium, which satisfies the following assumptions: The coefficient of the second term in (7), T 1 [S], describes the new favorable direction of the cells and we consider the positive taxis towards the chemoattractant.Throughout this paper, we employ the turning kernel operator introduced in [3], which represents a small perturbation of the equilibrium state: where a + := max(a, 0).Consequently, the kinetic equation ( 5) becomes: We note that the use of this local kernel may result in solutions that blow up in finite time.On the other hand, alternative turning kernels have also been studied; see, e.g., [7,31,32].It should be observed that the solution properties depend on the choice of the kernel.For instance, it was shown in [7] that the kinetic model ( 5), (2) with certain global kernel operators has a global solution for any initial mass.On the other hand, if the local kernel operator (9) is used, the solution may blow up.Indeed, as it has been proven in [3] for the parabolic-elliptic case, if the total mass M > M c = 16π |V | , then the solution blows up, while if M < m c = 0.806π |V | , then a classical solution exists globally in time.In [25], a different global kernel that reflects the cell detection and response to an external signal in finite size radius was proposed.In this case, the limiting PKS-type system contains a nonlocal sensing term and the existence of the bounded global solution was proven.
We also note that in the parabolic-parabolic case, the criteria for blow up or global existence are not explicitly known.In this paper, we develop an AP numerical method for the kinetic chemotaxis model (10), (2) and then use it to numerically investigate possible blowup scenarios in the parabolic-parabolic case.
3. Numerical method.In this section, we present an AP scheme for the kinetic equation (10) coupled with the chemoattractant concentration equation (2).To this end, in §3.1, we first rewrite the system (5), (2) in the form convenient for numerical simulations using an even-odd formulation; see, e.g.[36,6].Then, in §3.2, we implement the Strang splitting approach, [53], by separating stiff and non-stiff parts of the system.In this setup, the stiff subsystem is solved exactly as described in §3.2.1, while the non-stiff subsystem becomes a system of linear transport equations, which is solved by a second-order upwind method presented in §3.2.2.

3.1.
Even-odd formulation.In this section, we follow [36,6] and introduce new variables r 1 , j 1 , r 2 and j 2 by considering the so-called even-odd formulation.We assume that v ∈ From now on, we consider v ∈ V + only and rewrite equation (10) as the system of four equations obtained by substituting (u, v), (−u, −v), (u, −v), and (−u, v) into (10): where f (±u, ±v) is used instead of f (x, t, ±u, ±v) for the sake of simplicity.We then define the new variables with a one-to-one correspondence between them and f : It is instructive to point out that the macroscopic cell density ρ can be obtained from ( 4) and ( 12) in terms of the new variables r 1 and r 2 : Substituting ( 13) into (11), yields the following system: Since the left-hand sides of the second and fourth equations in ( 15) include stiff terms with the 1 ε 2 coefficients, we add and subtract u(r 1 ) x − v(r 1 ) y and u(r 2 ) x + v(r 2 ) y from the second and fourth equations, respectively, so that we finally obtain the following system for r 1 , j 1 , r 2 and j 2 : in which, all of the stiff terms are moved to the right-hand side.

3.2.
Operator splitting.In order to develop an efficient numerical method, we implement the operator splitting, in which the left-hand (non-stiff) and right-hand (stiff) sides of the system ( 16) are treated separately.To this end, we first introduce the vector W := (r 1 , j 1 , r 2 , j 2 ) T and write the system ( 16), (2) in the following form: where and We then implement the splitting approach by considering the following two subsystems: We note that in the subsystem (19), only the W variable is evolved in time while S remains unchanged there.Assuming that the solution at time t is available, we evolve it to the next time level using an operator splitting algorithm, [44,45,53], of either the first order: or the second-order: Here, L 1 and L 2 stand for numerical solution operators for the stiff, (19), and non-stiff , (20), subsystems, respectively.
Remark 1.It should be observed that the order of the operators in ( 21) and ( 22) is interchangeable.
Remark 2. For simplicity of presentation, we will describe the numerical schemes used to solve each one of the subsystems (19) and (20) in the context of the firstorder operator splitting (21).We also use the first-order splitting in the proof of the AP property of the proposed scheme.However, in all of the numerical experiments reported in §5, we implement the second-order operator splitting (22) to decrease the effect of splitting errors on the computed solutions.
Before proceeding with the description of numerical methods for the subsystems ( 19) and ( 20), we consider a computational domain, Ω×V + , where ] of size ∆x∆y with the cell centers (x i , y j ) = (x i− 1 2 +∆x/2, y j− 1 2 +∆y/2); i = 1, . . ., N x , j = 1, . . ., N y .We also introduce a uniform grid in the velocity domain V + consisting of N θ grid points: We also denote by 3.2.1.L 1 : Numerical solution of the stiff subsystem (19).In this section, we present a numerical scheme for the stiff subsystem (19).We start by solving the equations for r 1 and r 2 , (r keeping in mind that the chemoattractant concentration S satisfies τ S t = 0.It is instructive to point out that not only S, but also the macroscopic cell density ρ does not change in time during this substep.Indeed, from ( 14) and ( 24) we obtain which follows directly from (8) and the definition of J in (10).We assume that the numerical solution is available at time level t = t n and denote by Taking into account that ρ * i,j = ρ n i,j and S * i,j = S n i,j , we obtain the following semidiscrete approximations for (r 1 ) * i,j,k and (r 2 ) * i,j,k from (24): ) Here, we approximate the derivatives (S x ) n i,j and (S y ) n i,j with the central differences, and compute |∇S n i,j | = ((S x ) n i,j ) 2 + ((S y ) n i,j ) 2 .The linear ODEs in (26) are then solved exactly in time to obtain where . We now solve the equations for j 1 and j 2 : Equipped with the updated values of (r 1 ) * i,j,k and (r 2 ) * i,j,k , we use the central differences to compute substitute them into ( 28) and arrive at the following semi-discrete approximations for (j 1 ) * i,j,k and (j 2 ) * i,j,k : Finally, we solve the linear ODEs (30) exactly to obtain (j 1 ) * i,j,k and (j 2 ) * i,j,k ) T (for the sake of brevity, we omit the precise formulae).
3.2.2.L 2 : Numerical solution of the non-stiff subsystem (20).In this section, we describe the numerical solution operator L 2 and denote by .
As one can see from (20), the linear hyperbolic system with the constant coefficient matrices A 1 and A 2 given by ( 18) is, in fact, decoupled from the chemoattractant concentration equation.The latter is either the Poisson equation α∆S = βS − γρ, (32) if τ = 0, or the parabolic equation if τ = 1.Therefore, we will first evolve the solution of (31) to obtain W n+1 i,j,k , from which the values of the macroscopic density ρ n+1 i,j will be calculated and used in either (32) or (33) to compute S n+1 i,j .Upwind method for (31).We begin with the derivation of a second-order semidiscrete upwind approximation for the system (31).To this end, we first introduce the matrix which is used to simultaneously diagonalize the matrices A 1 and A 2 .We then define a change of variables U := Q −1 W and rewrite the system (31) in the diagonal form: where B 1 and B 2 are the following diagonal matrices: Next, we split B 1 and B 2 into the sum of non-negative and non-positive definite matrices: ), and rearrange (34) in an equivalent form: Multiplying the last equation by Q, we recover the original system (31): where We note that A + 1 and A + 2 are non-negative definite, while A − 1 and A − 2 are nonpositive definite so that and one can easily design upwind finite-difference schemes for the system (35).According to the upwind approach, we introduce the secondorder forward and backward finite-difference approximations for the spatial derivatives in (35): which are then used to construct the following second-order semi-discrete upwind scheme for (31): The system of time dependent ODEs (36) should be numerically integrated in time using a stable and sufficiently accurate ODE solver.For example, using the firstorder forward Euler method, a fully discretization of ( 36) can be written in the flux form as follows: It is important to stress that according to the definitions in ( 12), both r 1 and r 2 (and hence ρ, see ( 14)) should be positive, which is not guaranteed unless the scheme (37), ( 38) is used with a very small (possibly impractical) timestep ∆t.We therefore implement a draining timestep technique, which was introduced in [2].To this end, we denote by where m = 1, 2, (•) + := max(•, 0), and δ is a small positive number taken to be δ = 10 −14 in the numerical examples reported in §5.We then replace the first (m = 1) and third (m = 2) equations in (37) with where ∆t We note that using ( 39)-( 41) we, in fact, limit the fluxes through the boundaries of spatial cell (i, j) to prevent too much mass living the cell yielding unphysical negative values of r 1 or r 2 .It can be easily verified that now (r 1 ) n+1 i,j,k ≥ 0 and (r 3 ) n+1 i,j,k ≥ 0, and thus ρ n+1 i,j , which is computed from ( 14) using the midpoint rule, is non-negative, namely, Remark 3.Even though, we show here only one forward Euler step, in all of our computations below, we solve the semi-discrete system (36) using the three-stage third-order strong stability preserving (SSP) Runge-Kutta method: see, e.g., [18,17,51].Since SSP methods consist of a convex combination of forward Euler steps, the computed values of (r 1 ) n+1 i,j,k and (r 2 ) n+1 i,j,k , as well as ρ n+1 i,j , are still guaranteed to be non-negative.
Spectral methods for (32) and (33).Equipped with the point values of the macroscopic density, ρ n+1 i,j at time level t = t n+1 , we implement the spectral method to update the values of the chemoattractant concentration S. To this end, we remind the reader that the Neumann boundary conditions are imposed for both S and ρ.Therefore, the discrete Fourier coefficients S ,m (t) and ρ ,m (t) can be computed from the available point values S i,j (t) and ρ i,j (t), respectively, using the fast cosine Fourier transform and the solution at time t can be approximated by Substituting ( 43) into (2) yields In the elliptic (τ = 0) case, the values S n+1 ,m are immediately computed from In the parabolic (τ = 0) case, equation ( 44) can be solved exactly on the interval [t n , t n+1 ]: t n ρ ,m (s)e −ω ,m (s−t n −∆t) ds.
We then approximate the integral in the last equation using the trapezoidal rule to obtain Finally, we use the inverse fast cosine Fourier transform to compute the point values {S n+1 i,j } out of the set of the discrete Fourier coefficients { S n+1 ,m }.

AP property.
As it was mentioned in the Introduction, the solutions of the studied kinetic-chemotaxis model are expected to converge to the corresponding solutions of PKS system as ε → 0. In this section, we show that the proposed numerical scheme for (10), (2) provides a consistent discretization of (1), (2) in the limiting ε → 0 case.In other words, the numerical method is AP.For the simplicity of the presentation, we prove the AP property for the first-order splitting (21) and note that a straightforward extension to the second-order splitting ( 22) can be derived.We first observe that when ε → 0 the equations in (27) reduce to We then substitute (47) into (30) and derive the following formulae in the ε → 0 limit: where (ρ x ) n i,j and (ρ y ) n i,j are obtained from ( 29) and ( 42) and equal to Next, we consider the first and third equations in the semi-discrete upwind scheme (36), which after the forward Euler time discretization read as Substituting ( 47) and ( 48) into (49), adding the above two equations and multiplying by 2 yield where we have used the following notations: . We now multiply (50) by v 0 ∆θ, sum it over all k, and use (42) to obtain which can be used to show that (51) provides a consistent approximation of (1).Example 1-Parabolic-elliptic system.In this example taken from [6], we consider the system (10), (2) in the parabolic-elliptic (τ = 0) and non-stiff (ε = 1) regime with α = γ = 1, β = 0, and F (v) ≡ 1/(2π).The system is solved on the domain Ω = [−2, 2] × [−2, 2] and subject to the following Gaussian-shaped initial data: where M is a total mass.According to the theoretical results in [3], there are two critical mass thresholds: if M > M c = 8, the solution blows up in finite time, while if M < m c = 0.403, a global classical solution exists.It is still, however, unclear whether the solution remains bounded if m c ≤ M ≤ M c .
We investigate the solution behavior by computing the ratio ||ρ|| ∞ /M for different values of M .We run the computations on a uniform grid with N x = N y = 128 until final time T = 6.The results are presented in Figure 1 for M = 1, 5, 7, 8 and 9.As one can see, when M = 1 the maximum density decays for all times while for other values of M the solution exhibits an initial growth.For M = 5 and M = 7, the maximum density decays at later times and the solution clearly remains bounded.At the same time, for M = 9, which is above the critical threshold M c = 8, the maximum density increases and eventually saturates.In fact, this solution blows up and its maximum saturation phenomenon is attributed to the fact that the magnitude of finite-difference approximations of δ-type singularities is always proportional to 1/(∆x∆y).The blowup is also confirmed by the data presented in the   Example 2-Parabolic-parabolic system: Blowup at the center.In this example, we consider the system (10), (2) in the parabolic-parabolic (τ = 1) and stiff (ε 1) regime with α = β = γ = 1, and F (v) ≡ 1/(2π).The system is solved on the domain Ω = [−1/2, 1/2] × [−1/2, 1/2] and subject to the following initial data: where M is a total mass.Example 2a.We first consider the case when the initial chemoattractant concentration is a Gaussian-shaped function given by S(x, y, 0) = 500e −50(x 2 +y 2 ) . ( We study the time evolution of maximum density for three different values of the initial mass, M = 1, 8 and 11.In Figure 3, we plot ||ρ|| ∞ as a function of time for four consecutive meshes.We conclude that the behavior of the solution depends on the total mass: if M is large enough, the maximum norm of the cell density grows rapidly and saturates after blowup due to the finite-difference approximation limitation as discussed in Example 1.However, for smaller M , the maximum density first increases and then decreases without blowing up. In this case, the spiky structure at the center of the computational domain develops much slower that in Example 2a, since S is not concentrated at the center initially.This seems to affect the value of the critical mass.For instance, when M = 8, the solution does not blow up in contrast to Example 2a; see Figure 5 (left), where we plot ||ρ|| ∞ as a function of time for three consecutive meshes.When larger values of M are considered, the solution blows up as expected; see Figure 5.It should be pointed out, however, that one can observe rapid change in the solution magnitudes for both M = 9.5 and M = 11, but these changes occur after the blowup times.The latter can be estimated as t = 0.24 (for M = 9.5) and t = 0.058 (for M = 11), since at these times the ratio of ||ρ|| ∞ computed on the  We conclude this example with an experimental convergence study.To this end, we take a small final time T = 0.01 and measure the experimental convergence rate for the density in the L ∞ -norm by computing is the estimated relative L ∞ -error, and ρ N is the solution computed on a mesh with N = N x = N y .
The results presented in Table 1 clearly demonstrate second order of accuracy of the proposed AP scheme.where M stands for a total mass.According to the analytical results in [20] and the numerical simulations in [8], the solution of the PKS model with the corresponding initial and boundary conditions moves towards the upper right corner of the computational domain and blows up there.In view of these results, it is interesting to numerically investigate whether the solution of the kinetic chemotaxis model behaves similarly.
In Figures 6 -8, we plot the density computed on the uniform grid with N x = N y = 128 for M = 3, 7 and 11.The time evolution is shown on each of these figures (from left to right).As one can see, in all of the cases, even when M = 3, the solution blows up though the blowup time is much smaller for larger values of M .One can also clearly see that the solutions propagate towards the upper corner of the computational domain.However, when M = 11, the solution blows up much earlier than it would reach the corner.We have also performed the same numerical experiments with ε = 10 −8 and the obtained results were very similar.This indicated that for larger initial masses the solutions of the kinetic chemotaxis model may not converge to the corresponding PKS solution.considered.The Patlak-Keller-Segel model for chemotaxis can be recovered in the singular limit of the kinetic model, when the mean free path converges to 0. It is well-known that the solutions of both the kinetic and Patlak-Keller-Segel models may blow up in finite time, when the initial mass is larger than a critical mass.The critical values, however, are different for macroscopic and kinetic models.
Computing these blowing up solutions requires a reliable and efficient approximation of singular functions or distributions.Our numerical method is based on the so-called even-odd decoupling followed by the Strang splitting and a suitable combination of the evolution of the chemoattractant, relaxation and transport steps for the evolution of cell density.Using the splitting allows one to exactly solve the stiff relaxation system whereas the transport system is linear and is numerically solved using an upwind approach.The macroscopic equations for the chemoattractant are treated by the spectral method.
We have proven that as far as the global solution exists, our numerical scheme is asymptotic preserving and yields a consistent approximation of the Patlak-Keller-Segel model when the mean free path converges to 0. Our numerical experiments indicate that in the blowing up regime, the solutions of the kinetic chemotaxis and Patlak-Keller-Segel systems may behave differently.To the best of our knowledge, this does not contradict the theoretical results since after the blowup time the solutions can only be understood as measure-valued solutions for which no rigorous analysis is available.
In the future works, we intend to study the kinetic model with global turning kernels and to analyze whether there is a global weak solution for both parabolic and elliptic chemoattractant equation for general two dimensional initial data.

5 .
Numerical results.In this section, we test the proposed AP scheme on several numerical examples and also study the behavior of the solutions of the kineticchemotaxis system (10), (2) in the ε → 0 regime.In all of the examples below, we take v 0 = 1 and N θ = 32.The parameter ε = 1 in Example 1 and ε = 10 −5 in all other examples.
Figure 2 (right), where we plot the time evolution of ||ρ|| ∞ computed on three consecutive meshes: as one can see, at time T = 6 the value of ||ρ|| ∞ increases by a factor of four as the grid is refined.This behavior of ||ρ|| ∞ is clearly different from the one observed in the case of M = 7 < M c shown in Figure 2 (left).

Figure 6 .
Figure 6.Example 3: The displacement of the density for M = 3.

Figure 7 .
Figure 7. Example 3: The displacement of the density for M = 7.

Figure 8 .
Figure 8. Example 3: The displacement of the density for M = 11.