Boundary layer analysis of nonlinear reaction-diffusion equations in a smooth domain

Abstract In this article, we consider a singularly perturbed nonlinear reaction-diffusion equation whose solutions display thin boundary layers near the boundary of the domain. We fully analyse the singular behaviours of the solutions at any given order with respect to the small parameter ε, with suitable asymptotic expansions consisting of the outer solutions and of the boundary layer correctors. The systematic treatment of the nonlinear reaction terms at any given order is novel along the singular perturbation analysis. We believe that the analysis can be suitably extended to other nonlinear problems.


Introduction
Nonlinear reaction-diffusion equations arise in many areas in systems consisting of interacting components. The equations describe, e.g., chemical reactions, pattern-formation, population dynamics, predator-prey equations, and competition dynamics in biological systems (see, e.g., [7,11,12,[31][32][33]39]). One can consider a typical form of systems of reaction-diffusion equations in the form u t = D∆u + g(u), (1.1) where g = g(u) describes a change or a local reaction of the state u and D represents a diffusion coefficient matrix. It is also possible that the reaction g may depend on the spatial domain variable x and of a derivative of u, i.e., g = g(x, u, ∇u).
In real applications like a fast reaction system, the magnitude of some coefficients in the diffusion matrix D is relatively small and hence the system can be singularly perturbed.
In this article, for the singular perturbation and boundary layer analysis aimed here, we consider the steady state system of (1.1) and study the following scalar nonlinear singularly perturbed problem which can serve as a guide for more general systems: Here, < ε ≪ , Ω is a general smooth domain, f = f(x, y) and g = g(u) are given smooth functions with g( ) = , g ὔ (u) ≥ λ > for all u ∈ ℝ. For small ε > , the solutions to (1.2) display thin sharp transition layers called boundary layers which are formed due to the discrepancies between the limit solutions when ε = (see (2.3) below) and the boundary conditions in (1.2). The discrepancies are inevitable because the limit problem (see (2.2) below) loses high order derivatives and hence in general its solutions cannot meet the boundary conditions. Then, the small diffusion term −ε∆u ε smoothes out the discrepancies, which leads to sharp transition boundary layers.
Another motivation of studying boundary layers is the vanishing viscosity problem in fluid dynamics, see, e.g., [2, 5, 6, 8, 13, 23-25, 27, 28, 30, 35-37]. The typical question is on the behaviour of the Navier-Stokes flows at small viscosity, i.e., the limit behaviour or convergence to Euler flows as the viscosity tends to zero. The boundary layers play a crucial role for connecting the Navier-Stokes and Euler flows and they also do so for the singular perturbation analysis in the nonlinear reaction-diffusion equations considered here.
An additional motivation comes from the computational aspects in numerical simulations. Due to the thin boundary layers, the computational meshes are classically refined near the boundary ∂Ω and this causes high cost in the simulations. Rather than refining meshes we propose to enrich with suitable boundary layer correctors the Galerkin or finite element space (or finite volume space). Then, we are able to use a coarse mesh and this reduces substantially the computational cost. See, e.g., [17, 18, 20-22, 38, 41] for the method of spaces enriched with boundary layer correctors. For singular perturbations analysis, see [15,26,42,44] and also the recent review article [14]. See other perspectives in singular perturbations and boundary layers in [3,4,9,10,16,19,29,43].
In what follows, we discuss the problems posed on a channel domain in Section 2 which is relatively easier to handle thanks to the simple geometry of the boundary. In Section 3, we cast the nonlinear reactiondiffusion equations in a general domain. We need to take into account the geometrical properties, like curvature, using the boundary fitted coordinates. Throughout this paper, we systematically handle the nonlinear term g along the singular perturbation analysis at any orders. This nonlinear treatment can apply to other nonlinear problems.
An exponentially small term, denoted e.s.t., is a function whose norm in all Sobolev spaces H s (Ω) is exponentially small with, for each s, a bound of the form c e −c /ε γ , c , c , γ > , with c i , γ depending possibly on s.
To give an idea on how to construct the boundary layers, for now we assume which, as we will see, reduces the boundary layer at x = L , so that only the boundary layer at x = persists. Thanks to (1.3) and (2.3), ≤ λ(u ) ≤ (g(u ) − g( ))u = fu , and hence u = at x = L and u (x, y) = u (x, y + L ). (2.5)

Boundary layer analysis at order ε
We now construct a zeroth order corrector to account for the discrepancy between u ε and u at x = . Formally, substituting u ε ∼ u + θ in (2.1) and subtracting (2.2) from (2.1), we find that Using the stretched variablex = x/ ε and dropping non-stiff small terms, we find the zeroth order corrector equation for θ : However, in general u ε − u does not satisfy the boundary condition in (2.1), and hence at the boundary x = , so we propose a boundary layer corrector θ satisfying Although θ is not known explicitly, unlike in many linear problems, we can derive pointwise estimates for θ .
We can also deduce some norm estimates.

Lemma 2.2.
There exists a constant c > , independent of ε, such that Proof. The second estimate of (2.9) directly follows from (2.7). To obtain the first estimate, we introducē θ = −u ( , y)e −x/ ε δ(x), where δ(x) is a smooth cut-off function with δ(x) = for x ∈ [ , L / ] and δ(x) = for x ∈ [ L / , ∞). We observe that for a.e. y ∈ ℝ, where the last inequality follows from the mean value theorem and the L -estimate ofθ xx ,θ , and θ . This implies the lemma.
We can also obtain the lower bound of |θ (x, y)|.

Boundary layer analysis at arbitrary order ε n , n ≥
Outer expansion. We now consider the higher order outer expansions u ε ∼ ∑ ∞ j= ε j u j . Substituting in (2.1) and using (2.2), we formally write (2.14) Dropping O(ε n+ ) terms, we have We identify at the order O(ε j ), j = , , . . . , n, and find We then obtain, e.g., More generally, we recursively obtain To construct the higher order correctors, we assume, for simplicity, that f is infinitely flat at x = L , i.e., using the multi-index notation α = (α , α ) with |α| = α + α . (2.18) This implies that the u j , j ≥ , are infinitely flat at x = L , that is Thus, we only have boundary layers at x = corresponding to u j .
Correctors. We now proceed with the determination of the correctors. Substituting We subtract (2.14) from (2.19) to obtain (2.20) We first need to handle the nonlinear term to identify the quantities of order ε j and this is discussed below.
Thus, we may write the multi-index notations as (2.23) Hence, using the multi-index notations in (2.23), we formally write For the analysis below, we estimate the truncation error corresponding to the expansion (2.24).
Lemma 2.6. There exists a constant C > , independent of ε, such that and the multi-index notations are given in (2.23).
Proof. We first note that the G n given above can be written as Thanks to the multinomial theorem, we observe that On the other hand, we find from Taylor's theorem that Here, ξ is between (u + θ ) and ∑ n j= ε j (u j + θ j ) and ξ is between u and ∑ n j= ε j u j . The lemma follows by observing that |R | ≤ Cε n+ .
We now define the boundary layer correctors θ j at order O(ε j ). From (2.20) and (2.24), using the stretched variablex = x/ ε at each order O(ε j ), j = , , . . . , we identify Dividing by ε j , rearranging terms in the latter equation at each order O(ε j ), and using the fact that we rewrite (2.28) as We supplement the boundary condition on θ j , for each j = , , . . . , by Remark 2.7. We note that the corrector equations for θ j , j ≥ , are all linear and this allows us to directly apply the maximum principle. Differentiating the equations in y, the maximum principle also holds for ∂ m θ j ∂y m (x, y) for m ≥ , j ≥ .
Proof. We use the maximum principle to prove the lemma. Let L be the linear operator given by We introduce a barrier function Ψ = C exp(− λ ε x), where C will be chosen later. We use the mathematical induction on m starting from the case j = . By (2.7), we already have |θ (x, y)| ≤ c exp(− λ ε x). We then see that by the mean value theorem, for some η between (u + θ ) and θ . We also find Since g ὔ (u + θ ) ≥ λ and |g ὔὔ (η)||θ ||u y | is bounded onΩ, we can find a positive constant C ὔ such that By the boundary conditions of θ , we obtain where C = max(|u ( , y)|, C ὔ ). The maximum principle implies that We suppose by induction that for k ≤ (m − ), m ≥ , there exists a positive constant C k satisfying We then find that with some multivariate polynomials h k and η k between (u + θ ) and u . Since g and u are smooth, there exists C ὔ m such that We infer from the boundary conditions for θ that and using the maximum principle we obtain that where P (θ ) = θ and P k (θ ) = ∑ α +⋅⋅⋅+(m+ )α m+ =k ∏ m+ i= (∂ i y θ ) α i for k ≥ with for some multivariate polynomial h k and η k between (u + θ ) and u . Since g and u j are smooth and h k are polynomials, we can find a positive constant C ὔ m , by the result for θ , such that By the boundary condition on θ , we also have We find from (2.35) and (2.36) that for m ≥ , by the maximum principle where C m = max(C ὔ m , | ∂ m ∂y m u ( , y)|). We now suppose by induction that for k ≤ (j − ), j ≥ , there exists a positive constant C jm such that To prove (2.31) at order k = j, differentiating (2.29) in y we note that the first and third terms of the right-hand side of (2.29) are similarly estimated as for the case θ by (2.37). We thus estimate the second term there.
Observing that for k ≥ , it suffices to show that for any m ≥ , To prove this, we note that and using the factorization a n − b n = (a − b) ∑ n i= a n−i b i− , where α = (α , . . . , α j− ). Differentiating (2.40) in y, thanks to the mean value theorem, the left-hand side of (2.39) can be written as the sum of the products of θ k and their derivatives in y for k ≤ (j − ). We then conclude, by assumption (2.37), that (2.39) holds true.
We hence proved the following convergence theorem.

Without the assumption (2.17)
If we consider a general smooth function f , i.e., if we remove the assumptions (2.17), we also expect similar boundary layers at x = L . Let us denote similar boundary layers θ j at x = by θ j l , j ≥ , given in (2.29). Similarly, we define the boundary layers at x = L , denoted by θ j r , which satisfy equations (2.29) but with different boundary conditions, i.e., Here, where R is given in (2.42)-(2.43) and .
We now note that R is exponentially small. Indeed, we find that, for all α ≥ , which is an exponentially small term. Then, the convergence analysis similarly follows as in the above, from which we infer the following theorem.

General domains
We now return to the case of a general smooth domain, where equation (1.2) is posed: Here, Ω is a general smooth domain, f = f(x, y) and g(u) are given smooth functions with The outer solutions u j , j ≥ , are the same as in (2.16). That is, in Ω, we have However, the boundary layers appear in the direction normal to the curved boundary. Thus, the boundary fitted coordinates, i.e., the normal and tangential components along the boundary, are necessary to devise the boundary layer correctors. Here, we consider smooth boundaries ∂Ω which are parametrised by an arclength η and we assume that (X(η), Y(η)) ∈ ∂Ω is a regular curve, i.e., the tangent vector T = (X ὔ , Y ὔ ) ̸ = for the arclength ≤ η < L is measured counterclockwise, where L is the length of the boundary ∂Ω. Hence, we may assume that it has unit speed, i.e., (X ὔ ) + (Y ὔ ) = (see [34]).

Boundary layer analysis at orders ε and ε
Unlike the channel domain (2.1), we will need to introduce a corrector θ to settle the error from the curvature κ(η).
In general, to find appropriate asymptotic expansions for the boundary layers, we preform the following expansions near the boundary ξ = . Subtracting (2.41) from (3.9) and observing that ε is the thickness of the boundary layers, as indicated in previous sections, which suggests to use the stretched variableξ = ξ/ ε, we find that To address the terms at all the orders of ε in the boundary layers, we have to resolve the effect of curvature κ(η).
Thanks to the mean value theorem, we observe that D α (g(u +θ ) − g(u )) = finite sum of products of D βθ with β ≤ (l − , m), (3.21) and we can thus inductively prove (3.19) for any l ≥ as long as the case for l = is proved.
We integrate over (ξ , ∞), to find that where Let μ = p λ, where p is a constant independent of ε to be chosen later. From estimates (3.19), we note that |θ ξ (t, η)| ≤ c exp(− λt/ ). Using the same argument as in the proof of Lemma 3.3, we only need to show that for < p < , Following then the same procedure as in the proof of Lemma 3.3 with the boundary conditions (3.16), we can obtain (3.24) forθ . For l ≥ , differentiating (3.15) inξ and using estimates (3.19), the lemma is inductively proved.
We now find the norm estimate in L and H forθ andθ in the next lemma. We now introduce the analogue of Theorem 2.3.

Theorem 3.7. Assume that f is a general smooth function and Ω is a general smooth domain. Then, there exists a positive constant c > such that
Proof. We define w = u ε − u − θ − εθ , then we find by (3.14) and by (3.15), where + g(u +θ ) + g ὔ (u +θ ) εθ − g(u +θ + εθ ). (3.33) We note that thanks to the Taylor expansion, we have g(u +θ + εθ ) − g(u +θ ) = g ὔ (u +θ ) εθ + T , On the other hand, We then find from (3.36) and the mean value theorem that for some ζ between u ε and (u + θ + εθ ). We then conclude

Boundary layer analysis at arbitrary orders ε n and ε n+ , n ≥
Similarly as in (2.20), we formally write −ε∆ ∞ j= ε j (θ j + εθ j+ ) + g ∞ j= ε j (u j + θ j + εθ j+ ) − g ∞ j= ε j u j = . (3.37) Because of the one-dimensional nature of these boundary layers near the boundary in the direction normal to the boundary, we now introduce the boundary fitted coordinates. We transform the Laplacian ∆ as in (3.7) and (3.8).
Multiplying, respectively, (3.14) by ε , (3.15) by ε , (3.48) by ε j and (3.49) by ε j+ , from j = to j = n, and adding these resulting equations we find that where I n k are described as in (3.40).