Boundary layer analysis from the Keller-Segel system to the aggregation system in one space dimension

This article is devoted to the discussion of the boundary layer which arises from the one-dimensional parabolic elliptic type Keller-Segel system to the corresponding aggregation system on the half space case. The characteristic boundary layer is shown to be stable for small diffusion coefficients by using asymptotic analysis and detailed error estimates between the Keller-Segel solution and the approximate solution. In the end, numerical simulations for this boundary layer problem are provided.


1.
Introduction. Chemotaxis is the oriented movement of a cell, a bacteria or an organism toward or away from a chemical stimulus. It helps them find food or avoid poisons. In mathematical biology, many of the models for chemotaxis are partial differential equations, [23]. Among them, the so-called Keller-Segel (KS) type systems [32,33,45] are the basic models for chemotaxis. They are used to describe the collective motion of cells or the evolution of the density of bacteria, where both diffusion and aggregation appear within the system. The parabolic elliptic KS system is written as where ρ(x, t) is the density of cells, c(x, t) is the chemoattractant, ε is the diffusion coefficient and χ is the chemoattractant sensitivity. Without loss of generality, we take χ = 1. Since the work by Jäger and Luckhaus, [30] in 1992, it has been widely studied in applied mathematics. Horstmann in 2003 gave detailed review on such systems in [24,25].
The system has two important properties, namely, the conservation of mass From the entropy-entropy dissipation relation, one can see that the total entropy decays in time. However the entropy is the difference between the quantities for the diffusion and for the aggregation, which hints clearly that the system allows both global existence and the finite time blow up of the solution. In the one dimensional case, the effect of the diffusion is stronger than that of the aggregation, therefore one can get the global existence of solutions [39,44]. In the two dimensional case, due to the Logarithmic potential, the effects of the diffusion and aggregation are "competing" at the same level. In [21], the authors obtained the critical mass 8π (with ε = 1). More precisely, if the initial mass is smaller than 8π, the solution exists globally; otherwise the solution blows up in finite time. Besides that, the two dimensional case was extensively studied by several authors, for example, [13,18,26,39,40,41,42]. For the higher dimensional case, the aggregation from Newtonian potential is "stronger" than the diffusion. The solution exists only for "small" initial data (i.e. ρ I L n 2 is small) and the blow-up could easily happen for a big class of initial data [8,9,30]. Many authors also tried to use a porus media type degenerate diffusion ∆ρ m , m > 1 (which means that the diffusion coefficient depends also on the density) to compete with the aggregation from the singular potential 1/|x| n−2 , [1,7,10,16,17,19,27,29,34,36,50,52,53,54]. Among them, the mass critical exponent 2 − 2/n (cf. [52]) and the energy critical exponent 2n/(n + 2) (cf. [16]) have been obtained. We refer more detailed discussions on the critical exponents to the review paper [57]. Now we look at the corresponding aggregation equation, i.e., the diffusion coefficient ε = 0, ρ t = −div(ρ∇c), (1.4) −∇c = ∇K * ρ, where K(x) = C(n)1/|x| n−2 is the fundamental solution of the Poisson equation. This problem still keeps the property of mass conservation (1.2). In [4], by considering the mass contained inside a ball of radius r, M (r, t) = B(r) ρ(x, t)dx, the authors transformed the system into the Burgers equation on a half line and proved that the radially symmetric solution blows up in finite time.
The aggregation problem (1.4), with a different interaction potential K(x), arises from several physical phenomena such as [20,37,38] for 2D vortex dynamics, [11,12,55] for flocking and swarming models in biology, [2,15,56] for the interacting models in granular media. In [35], it was shown that the local and global existence of the solution to (1.4) depends on the regularity of K. For "acceptable" potentials, the solution exists locally. One typical example of "acceptable" potential is K(x) = e −|x| . We refer the detailed definition of the "acceptable" potential to [35]. If furthermore the potential is in W k+2,∞ (R n ), the solution exists then globally. One example is K(x) = e −|x| 2 . For different choices of K, many of the works have shown that the smooth solutions generate finite time blow-up, see the references in [3,5,14,22,31]. In [14], under structured conditions on the interaction potential K, the global existence of weak measure solutions was proved. For more detailed results on aggregation equations, one can check the review paper [6]. If the potential has a special form K(x) = |x| α , 2 − n ≤ α < 2, smooth solutions are known to develop singularities in finite time. Furthermore, in [4,22], the critical space L p for initial data has been introduced, with which the solution immediately develops a Delta singularity.
In the current paper, we focus on the vanishing viscosity limit of the problem (1.1) in the one space dimensional case. By taking ε → 0 formally in (1.1), the limiting problem is (1.4), where K(x) is the Newtonian potential. For radially symmetric solutions, with small ε > 0, we expect that the solution behavior of (1.1) near the origin is mainly dominated by the solution of (1.4). While for the O( √ ε) order, an internal correction term has to be added. We will clarify this observation in the rest of the paper in the one dimensional case. The radially symmetric multidimensional case will be left for future discussion. The motivation of this work is to investigate the behavior of the solution to the problem (1.1) when the diffusion effect is relatively smaller than the effect of the aggregation. From the modeling point of view, it corresponds to the group of cells (or bacteria) which are more active in response to the chemical stimulus than to their own random walks. From another point of view, for fixed small ε 1, it costs a lot of energy and time to realize the numerical solution to (1.1). An asymptotic approximation of the solution with a detailed error estimate would be more helpful to understand the behavior of the solution itself.
For radially symmetric solutions in the one dimensional case, the problem is reduced to the analysis of the boundary layer (for small ε) to the following half line problem with homogeneous Neumann boundary condition, where the initial data ρ I ∈ L ∞ (R + ) ∩ L 1 (R + ). The formal limiting problem is In fact, both problems (1.5) and (1.6) have been studied in the literature: the wellposedness of (1.5) for fixed ε in [39,44], the blow-up behavior of the solution in [11] by transforming the problem into a Burger's equation. However, in order to study the boundary layer, we need detailed estimates forρ. Therefore we reprove the solvability of (1.6) by using the method of characteristics and additionally give a series of estimates forρ. Our goal in this paper is to study the vanishing viscosity limit of (1.5). Before we state our main result in this article, a short overview of the general vanishing viscosity problem with boundary layers will be given in this paragraph. Except in very simple situations, the viscous equations with small parameters are too difficult to solve. At the beginning of the 20th century, Prandtl proposed a division of the problem to make it more tractable. It was hypothesized that the viscous effects are only important in the thin region adjacent to the solid boundary, which is called the boundary layer phenomenon, and as long as the layer remains thin, a number of approximations can be constructed so that the equations are greatly simplified. To study the structure of the viscous solution, one can perform by combining the inviscid solution for the outer flow with a viscous boundary layer solution to the flow near the surface of the object. Although this method turns out to be powerful in many practical cases, it is still worth carrying out a general and systematical rigorous proof in mathematics. Since there is a loss of boundary conditions when the viscosity coefficient ε goes to zero and hence a boundary layer appears. Thus, it is commonly believed that the solutions for the viscous parabolic equations cannot be uniformly close to those for the inviscid hyperbolic equations. There are immense literatures on this aspect of the theory, see [47,49,58,62] and references therein. Most of the initial boundary value problems for fluid dynamics systems arise in various domains, and the boundary conditions for these problems are chosen according to the physical properties in each situation. There are two distinguished classes of boundary conditions, which yield substantially different behaviors, depending on whether the boundary is characteristic or not. For example, the in-flow and out-flow boundary conditions make the boundary be noncharacteristic, and in this case, there are boundary layers of size ε which are stable when the amplitude is small, see [48,51,59]. From the physical point of view, for Navier-Stokes equations, a natural assumption is that the boundary is impermeable; this makes the boundary be characteristic. More generally, for a one dimensional system where f (x, t) is the source term, we say that the boundary is characteristic for some state u(t, x) if det A(u(t, x)) = 0. If there is a loss of boundary condition in the vanishing viscosity limit process, it is pointed out in Prandtl's theory that characteristic boundary condition will develop the so-called characteristic boundary layers, with thickness of O( √ ε), and many nonlinear phenomena may occur in the layers. In this case, the leading boundary layer functions appear in the O(1)-term of the asymptotic approximate solution, and they satisfy a set of nonlinear Prandtl type equations, see details in [47,60]. From the mathematical point of view, there is no satisfying theory for the wellposedness and the justification of the Prandtl boundary layer equations, and it still remains open for general cases.
Here we study the asymptotic behavior between the solution ρ ε to the system (1.5) and the solutionρ to the inviscid system (1.6) on the time interval before blow up. Using (1.5) 2 and (1.6) 2 , the system (1.5) and (1.6) can be changed into a single equation (2.1)-(2.2) and (3.1)-(3.2) with a nonlocal term x 0ρ dξ or x 0 ρ ε dξ being a coefficient. Thus, the boundary {x = 0} is characteristic. The boundary layer equation in such a case is not an ordinary differential equation, but a partial differential equation. We use the two-scale asymptotic expansion to obtain the inner and boundary layer equations. More precisely, we build up an approximate solution whereρ(x, t) is the solution to (1.6), ρ 1 (x, t) and ρ 1 b (x/ √ ε, t) are higher correction terms of outer and boundary layer expansions respectively, which will be discussed at length in Section 3.1. It will be clear from our analysis that the O(1)-term boundary layer functions ρ 0 b is identically 0. Consequently, ρ 1 b as a leading boundary layer profile of order √ ε satisfies a linear equation (3.9). The main part of the proof is an energy estimate of the error equations given in Section 3.2. It is easy to see that a viscosity of size ε is not sufficient to stabilize boundary layer of size of O( √ ε), that is, the usual way of using the Poincare-type inequality to control the first-order term in the noncharacteristic case is not sufficient due to the viscosity term. Since the nonlinearity in the viscosity term may cause some difficulties, we perform a modified weighted energy estimate inspired by [28,61] on the error term and get the stability result. This eventually implies where C is a positive constant which is independent of ε. Moreover, for any given The arrangement of this paper is the following. In Section 2, we give the solvability of the aggregation problem (1.6) and develop a-priori estimates for the solution, which are going to be used in the stability estimates for the approximate solution. In Section 3, we carry out the boundary layer analysis and give the convergence rate estimates. In Section 4, we compare the behaviour of approximate solution, obtained from boundary layer analysis in Section 3, with its corresponding Keller-Segel and aggregation solution by the numerical simulations.
2. Solvability of the aggregation system. As mentioned in the introduction, in the multidimensional case, the existence of a unique radially symmetric solution for aggregation equation with Newtonian potential was obtained in Section 4 of [4]. The methods that have been used in [4] were to reduce the aggregation equation to the classical inviscid Burgers equation for m(r) = B(r) ρ(x)dx, the mass contained inside a ball of radius r. However, the existence result and solution properties obtained for m(r) in that article can not be directly applied in our current work. Therefore, in this section, we reprove the solvability and finite time blow up result of the aggregation system in one space dimension without using the transformation m(r). The methods are based on characteristics and a modified Picard iteration. More importantly, we arrive at a series of estimates which are useful in the next section for the asymptotic analysis.
Notice thatc x = − x 0ρ (ξ, t)dξ, the aggregation system (1.6) can be reduced into a differential-integral equation, The characteristics of this equation satisfy Along the characteristics, the solution of (2.1) can be represented as With these equivalent representation of the solution, we arrive at the following existence theorem of this section.
, which satisfies the following equation Step 1. Design an iteration scheme.
) in the following iteration processes.
We first estimate I n (t) = sup x0≥0 |x n+1 (x 0 , t) − x n (x 0 , t)| by using the definition of x n in the following Next we get an iteration for the second term on the right hand of (2.9).
where for fixed ξ,x 0 andx 0 are defined by By using the mean value theorem and Fubini theorem, we have where we have used the definition of x n and the notation ξ θ , µ θ,s > 0, i.e., For any fixed ξ ≥ 0, due to the fact that x n−1 (x 0 , s) = x n−2 (x 0 , s), we havê which implies that Hence, for s ≤ δ 2 ρ I ∞ , we have Therefore, for s ≤ δ 2 ρ I ∞ , (2.10) can be bounded by Noticing that a combination of both iterations in (2.9) and (2.12) gives the following estimates, where n 2 = n 2 when n is even and n 2 = n−1 2 when n is odd. In the end, after all the iterations, we have Therefore, by standard argument, we know that there exists a unique x(x 0 , t) which is defined for all (x 0 , t) ∈ R + × [0,T ] such that With the help of the iteration in (2.12), one can easily obtain that there existsρ(·, t) such that ρ n (·, t) →ρ(·, t) in L 1 , ∀t ∈ [0,T ].
Step 3. Unique solvability of the problem within the time interval [0,T ]. We need to show that the limit (x(x 0 , t),ρ(x, t)) which is obtained from the last step is the unique solution of the problem (2.5), or (x(x 0 , t),ρ(x, t)) is the solution for (2.3) and (2.4).
Since ρ I is continuous, we can take the limit in (2.7) to get the existence of the solution for (2.1). The uniqueness of the solution can be easily obtained by using the same idea as what has been used in (2.9) and (2.10).
Step 4. The existence time interval which has been proved in the last step is In the following by using the standard energy estimates, one can obtain the regularity ofρ with the additional assumption ρ I ∈ C 3 and (1 + x 2 ) d α ρ I dx α ∈ L ∞ [0, T ]; L 2 (R + ) for 0 < T < T * = 1/ ρ I ∞ with α = 1, 2, 3. We take derivative ∂ x directly on both side of (2.1), multiply it by (1 + x 4 )∂ xρ and integrate it over [0, ∞), By using integral by parts and Young's inequality
3. Boundary layer analysis. In this section, we investigate the boundary layer behavior between (1.5) and (1.6) as ε → 0. Actually, the Keller Segel system (1.5) can be rewritten into and its corresponding inviscid problem is exactly (2.1). As mentioned in the introduction, the existence of the solution for the one dimensional Keller-Segel problem on the whole space has already been proved. However, for completeness we still give a-priori estimates of this half space viscous problem here. Multiplying the viscous equation by ρ ε and integrating over R + , we have For any fixed ε, Thus by Gronwall's inequality, ∀τ > 0, where the constant C depends on ρ ε L 1 , ε and τ . Moreover, set w ε = ∂ x ρ ε and then w ε satisfies Similarly, we can obtain here the constant C also depends on ρ ε L 1 , ε and τ . We aim to study the asymptotic relation between ρ ε andρ as ε → 0 for 0 ≤ t ≤ T < T * = 1/ ρ I ∞ , i.e. before the blow up time of the aggregation problem (1.6). Since there is a loss of boundary condition, the boundary layer may appear at x = 0.

Construction of approximate solution to viscous problem.
3.1.1. Inner expansion. Away from the boundary {x = 0}, the viscous solution to (3.1)-(3.2) is expected to be approximated by the following regular series Direct calculation gives with the following initial conditions For simplicity we only write the remainder term as R in = R in (ε, ρ 0 , ρ 1 ), and the estimate of R in can be obtained by the properties of ρ 0 and ρ 1 , which will be given later. It is easy to see (3.3) is exactly the equation (2.1), then ρ 0 is exactlyρ and its regularity has already been obtained from theorem 2.1 that Proof. For convenience, we denote by C a positive constant depending on T, ρ 0 ∞ , We multiply (3.4) by (1 + x 2 )ρ 1 and integrate over (0, +∞) to obtain where I i , i = 1, 2, 3, 4 can be estimated as follows.
where we have used Hölder's inequality. It is easy to see and It then follows from Gronwall's inequality that Then taking ∂ j x , j = 1, 2 on (3.4) and using the similar arguments give Therefore the regularity of ρ 1 can be obtained by the Sobolev embedding theorem.
Together with Theorem 2.1 , we have Boundary layer expansion. We use the following two-scale expansion to approximate the viscous solution uniformly up to the boundary: where ρ in is the inner solution constructed in the previous subsection. Then In view of and R bd ∼ O(ε). We use the same expansion in the boundary condition to have ). Therefore, ρ 0 b (y, t) satisfies the following initial boundary value problem It is easy to see that ρ 0 b ≡ 0 is the unique solution to the above initial boundary value problem. Then the initial boundary value problem of ρ 1 b reduces to (1 + y 2k )|∂ 2α y ρ 1 b | 2 dy < ∞ with α = 1, 2. Therefore y∂ y ρ 1 b ∞ < ∞. Proof. We multiply the equation (3.9) by (1+y 2 )ρ 1 b , and integrate over R + to obtain 1 2
3.2.1. The error equation. Set φ = ρ ε − ρ a . Then by the property of ρ ε and ρ a , φ solves There exists an absolute positive constant C 1 , such that for some uniform constant C 2 .
Proof. Multiplying φ on both sides of (3.14) and integrating over R + give where we have used the boundary condition (3.15). By (3.12) and (3.13), we have Gronwall's inequality implies that To justify the pointwise estimate of φ, we need to derive the L 2 -estimate on the spatial derivative ∂ x φ.
Proof. Let ψ = ∂ x φ. Then ψ solves Similar estimates as in Lemma 3.3 lead to and Gronwall's inequality, we have In view of the structure of the approximate solution ρ a , we can easily get (1.7).
4. Numerical simulations. In this section, we illustrate the behavior of the approximate solution to the Keller-Segel system with respect to its corresponding aggregation problem. The presented numerical simulations focus on the two-scale expansion of the approximate solution near the boundary, where the boundary layer may appear due to the loss of boundary condition. In this context, the simulations of the Keller-Segel system (1.5), its corresponding aggregation problem (1.6) and the approximate solution by solving the inner expansion (3.3)-(3.4) and boundary layer expansion (3.9) are performed. All simulations are restricted to the unit interval Ω = [0, 1] with the appropriate boundary condition. We employ an equidistant spatial discretization using a finite difference method with an explicit forward evolution in time. We discretize the xt− plane by choosing a mesh size h := ∆x and a time step k := ∆t, and define the discrete mesh points (x i , t n ) by x i = ih, i = 0, 1, 2, . . . , N, t n = nk, n = 0, 1, . . . .
The discretization of the Keller-Segel equation (3.1) together with the boundary conditions (3.2) will be expressed as (4.2b) Similarly, the discretization of the aggregation equation (2.1) can be written as f (x, t n )dx is evaluated by the trapezoidal rule. The integrand function f will be ρ andρ in the equations (4.2) and (4.3), respectively.
The approximate solution consists of two parts, namely, the inner expansion and the boundary layer expansion. The inner expansion for the approximate solution can be computed from equations (3.3) and (3.4). These equations have exactly the same structure as the aggregation equation (2.1), therefore their discretized version will be close to the equation (4.3). Now, only the boundary layer expansion remains to compute in order to get the approximate solution, as expressed in the equation (3.7). To get the boundary layer expansion of the approximate solution, we scale the domain Ω by √ ε and then solve the equations (3.8) and (3.9) on the scaled domain. So, the spatial discretization will be transformed as y i = ih ε , i = 0, 1, . . . , N. In this way, we get ρ bd (y, ε , t) corresponding to each inner value expansion ρ in (x, t). Since the first part of the boundary layer expansion ρ 0 b (y, t) is identically zero, only the term ρ 1 b ( x √ ε , t) needs to be solved from the equation (3.9). Then, the discretization of equation (3.9) will be with the initial condition  Example. The implementation of the described discretization has been performed with the initial condition ρ I (x) = 1 − x. In Figure 1, the time evolution of the Keller-Segel (viscous) solution, its corresponding aggregation (inviscid) solution and the two-scale series approximation (3.7) of the Keller-Segel solution are shown.
One may observe that near the boundary x = 0, the approximate solution is close to the Keller-Segel solution, whereas the aggregation solution blows up in finite time. We perform the test for different value of ε and observe the error between the Keller-Segel solution and its approximate solution, as given in equation (3.16). We take the spatial mesh of size h = 0.001 for these simulations. The size of time step k has been chosen such that k = ∆t < 0.5h 2 ε . The time integration has been performed for k = 5 × 10 −6 , which is found to be optimal for our numerical discretization. The respective error bounds for the different value of ε are presented in Table 1.
We have also compared the computation time between the simulation of the Keller-Segel equation and its approximate solution obtained from the boundary layer analysis. We expect that for smaller values of ε the computation of the approximate solution could be performed on larger time scales. This means that the simulation of the approximate solution will be computationally efficient compared to the solution of Keller-Segel model. We perform the simulation of Keller-Segel equation on a fix time step of size k = 5 × 10 −6 . Then, the time interval is scaled by a factor k f ac to simulate the approximate solution. The elapsed time is given in the Table 2  The grid convergence of the implementation has also been examined. We illustrate the grid convergence index (GCI) [43,46] procedure for this implementation. GCI yields the percentage of absolute error between the actual computed value and the asymptotic numerical value. We will consider the error bound between viscous and its approximate solution, as the computed value to observe the GCI. The value of GCI shows that how much the solution would change with a further refinement of the grid. A small value of GCI indicates that the computation is within the asymptotic range, therefore grid refinement may not improve the solution.
For a chosen grid of size h, we compute the error bound and then compute the relative error of the error bound (4.8) on halving the grid, given as For different grid spacing, the simulation has been performed for ε = 0.01 and time step of size k = 5 × 10 −6 . The respective values of E h and its relative error are given in the Table 3.  Table 3. Grid convergence.
The order of convergence is observed from these results, p = log e log 2 = 0.9471. (4.10) Now, the GCI with respect to the fine grid solution can be computed as where r = 2 is the refinement factor of the grid and F s = 1.25 is a factor of safety [43]. Similarly, the GCI with respect to coarse grid is Subsequently, the decisive indicator, GCI coarse /(2 p × GCI f ine ) = 0.9 is close to one. It indicates that the solution behave well in the asymptotic range.