## Abstract

We study the long-time asymptotics of reaction–diffusion-type systems that feature a monotone decaying entropy (Lyapunov, free energy) functional. We consider both bounded domains and confining potentials on the whole space for arbitrary space dimensions. Our aim is to derive quantitative expressions for (or estimates of) the rates of convergence towards an (entropy minimizing) equilibrium state in terms of the constants of diffusion and reaction and with respect to conserved quantities. Our method, the so-called entropy approach, seeks to quantify convergence to equilibrium by using functional inequalities, which relate quantitatively the entropy and its dissipation in time. The entropy approach is well suited to nonlinear problems and known to be quite robust with respect to model variations. It has already been widely applied to scalar diffusion–convection equations, and the main goal of this paper is to study its generalization to systems of partial differential equations that contain diffusion and reaction terms and admit fewer conservation laws than the size of the system. In particular, we successfully apply the entropy approach to general linear systems and to a nonlinear example of a reaction–diffusion–convection system arising in solid-state physics as a paradigm for general nonlinear systems.

## 1. Introduction

Many physical, chemical and biological processes are driven, on the one hand, by diffusion, a random particle motion microscopically described by a Brownian stochastic process, and, on the other hand, by reactions representing instantaneous interactions between particles. Typical examples where both the mechanisms occur are chemical kinetics, population dynamics, flame propagation and combustion, movement of biological cells in plants and animals, and charge-carrier transport in semiconductors.

To reiterate classical mathematical modelling, let denote the concentration vector describing *N* interacting species, where denotes the position variable and *t*>0 time. Then, according to Fick's law, the diffusion part of the motion is described by the following linear parabolic evolution equation:where *D*=*D*(*x*, *t*, *U*) is a positive definite, symmetric diffusion matrix, in general, depending on position *x*, time *t* and on the concentration vector *U* itself.

On the other hand, given a reaction process in terms of a ‘local’ dynamical system of the form *U*_{t}=*F*(*x*, *t*, *U*), the interaction between both reaction and diffusion leads to the following reaction–diffusion (RD) systems:

If the two processes occur in a spatially bounded domain *Ω*, suitable boundary conditions have to be imposed, for instance, *Dirichlet conditions U*=*U*_{D} on ∂*Ω* or, with *n* denoting the outer normal, *Neumann conditions* on ∂*Ω* and also linear combinations of either conditions (*Robin conditions*). Note that the total number of particles is conserved if zero outflow, i.e. *U*_{N}=0, is prescribed through the boundary ∂*Ω*.

Also, RD equations posed on the whole *n*-dimensional space occur in many important applications. In these cases, a confinement of particles occurs often by means of an exterior potential acting as a trap. A classical example of a confined diffusion process is the Fokker–Planck equation describing the interaction of Brownian diffusion with a harmonic (i.e. quadratic) potential as confinement.

RD systems are nowadays considered a classical topic, going (at least) back to the pivotal works of, for example, Fisher (1937) and Kolmogorov & Petrovsky (1937). The mathematical literature on RD equations is vast and we shall refer to classical textbooks (e.g. Smoller 1983; Rothe 1984; Amann 1985; Murray 2003 and references therein) for typical questions such as the existence of solutions and their global boundedness, stability and asymptotics, travelling waves, geometry and topology of attracting sets and singular limits.

From very early on, the question of stability versus instability of linear and nonlinear RD systems has been a key issue to mathematicians, physicists and biologists. In particular in 1952, Turing (1952) first pointed out a diffusion-induced instability of stable homogeneous reaction systems in chemistry (well known as *Turing instabilities*). More generally, classical mathematical analysis of the long-time asymptotic behaviour involves linearized stability techniques, spectral theory, perturbation and invariant regions arguments, and the Lyapunov stability (Conway *et al*. 1978; Fitzgibbon *et al*. 1997).

The approach presented in this paper is different and motivated by the recent great progress in the understanding of long-time asymptotics of scalar linear and nonlinear diffusion and diffusion–convection equations due to the so-called *entropy approach*. The entropy method refers to the general idea of a *functional inequality relationship* between an *entropy functional* of a system and its monotone change in time, usually called the *entropy dissipation*. Such an entropy–entropy dissipation inequality entails convergence to an entropy minimizing equilibrium state, at first in entropy and further in *L*^{1} using Cziszár–Kullback–Pinsker-type inequalities (cf. Arnold *et al*. 2001; Carrillo *et al*. 2001; Otto 2001; Del Pino & Dolbeault 2002). The entropy approach is *per se* a nonlinear method avoiding any kind of linearization and capable of providing explicitly computable convergence rates. Moreover, being based on functional inequalities rather than particular differential equations, it has the advantage of being quite robust with respect to model variations.

Related to the context of this paper, the entropy approach has already been applied to semiconductor drift–diffusion–Poisson systems (Arnold *et al*. 2000) or drift–diffusion–reaction–Poisson systems on bounded domains in (Glitzky *et al*. 1996; Glitzky & Hünlich 1997), but with the drawback of a proof based on an indirect contradiction argument without control of the rates and constants. We also mention the paper (Degond *et al*. 1997) dealing with general cross-diffusion systems. It was first in (Desvillettes & Fellner 2006, 2008, submitted) that explicit exponential convergence to equilibrium was shown via entropy methods for nonlinear RD systems modelling reversible mass action kinetics of two, three and four species. A general framework for RD systems (long-time asymptotic convergence, convergence rates, etc.), however, is still lacking.

In order to demonstrate the entropy approach, we consider the homogeneous Neumann problem for the heat equation on a bounded domain Owing to the homogeneous Neumann conditions, all constants are stationary states and the equilibrium state is determined by the conservation of the initial masswhere |*Ω*| denotes the volume of the domain. Multiplying the heat equation by *u* and integrating by parts yields the entropy (free energy) dissipation equationwhere we have used the *H*^{1}(*Ω*)–Poincaré inequality with constant , which is the spectral gap of the homogeneous Neumann–Laplace operator. Hence, after integration in time, we obtain the sharp decay estimate

To mimic the effect of a stable reaction term, we now add a linear absorption term with a constant rateClearly, this shifts the spectrum, and convergence to the unique equilibrium state *u*_{∞}=0 ( is no more conserved) follows from the sharp estimate

Thus, even in this most simplistic scalar example equation, diffusion and stable reaction do not ‘cooperate’ in the rate of decay to equilibrium (since constant states are not affected by diffusion). Nevertheless, full asymptotic cooperation between diffusion and stable reaction is observed for the intermediate asymptotic state ,

Later on, we shall see that linear systems combining ‘purely diffusive’ modes obeying conservation laws (such as the Neumann–Laplace problem) with ‘diffusive–reactive’ modes, which are governed by both diffusion and reaction (similar to the Neumann–Laplace absorption problem) are particularly interesting as far as (non-)cooperation of reaction and diffusion effects on the convergence rates are concerned. In §2, an explicitly computable 2×2 system with constant coefficients will demonstrate the complicated system-related interaction effects of diffusion and reaction. Such systems are typical in many physical and biological applications.

The aim of this paper is to start developing a framework for the quantitative analysis of the long-time asymptotics of stable reaction–diffusion–convection systems based on the entropy approach. In particular, for drift–diffusion–reaction systems on the whole space with confinement, the present paper provides—to our knowledge—the first attempt at this generality.

*Outline.* In §2, we outline the presented approach in the clearest possible way, namely for the systems posed on bounded domains with constant equilibrium states. In §3, we shall turn to the analysis of whole-space systems with confinement in each component. Both §§2 and 3 deal with linear systems. An example of a nonlinear 2×2 RD system with a quadratic reaction term and confinement will be considered in §4 to indicate a possible general strategy and the technical difficulties when dealing with nonlinear systems. However, a general theory for systems with nonlinear reaction is still being developed.

We summarize the main results contained in the paper as follows.

In theorem 2.14, we prove the exponential stability of constant equilibrium states for linear systems on bounded domains. The structural assumptions in this result prelude to those in (the more important) theorem 3.6, which deals with inhomogeneous systems on the whole spaces and is shortly resumed below. A significant aspect is that, unlike in theorem 3.6, in theorem 2.14 we can allow for degenerate diffusion.

A general result on the exponential stability of inhomogeneous equilibria for linear systems on the whole space with confinement is proven in theorem 3.6, which—loosely speaking (without going into the technical assumptions here)—reads as follows. Let the RD system be linear and appropriately symmetrizable, assume that the reaction is ‘stable’, i.e. the eigenvalues of the reaction matrix are either zero (of position-independent multiplicity) or uniformly negative on the bounded domain in , on which the system is posed, subject to homogeneous Neumann boundary conditions. Note that each zero eigenvalue of the reaction matrix corresponds to a conservation law. Then, for each

*L*^{2}initial state, there is a unique constant equilibrium state, to which the solution converges exponentially in*L*^{2}as time tends to infinity. The rate of decay can be characterized in terms of the coefficient matrices of the system. This characterization is typically more complex than a simple ‘combination’ of the diffusive and reactive rates.A refined convergence result for a linearized 2×2 drift–diffusion–recombination system with confining potentials is shown in theorem 3.9. Again, without going into technical details, the result states the following: let the linear RD system be posed on the whole space , with appropriate confinement in every component of the system. As above, let the reaction be stable, generating a certain number of conservation laws. Then, a unique integrable equilibrium state exists for all initial data. The solution of the initial-value problem converges exponentially to this equilibrium, again with a characterizable rate.

For a nonlinear drift–diffusion–recombination example system, we show exponential convergence towards the inhomogeneous equilibria in theorem 4.4. We consider a nonlinear drift–diffusion system with a charge-carrier recombination–generation term (subject to linear confinement) arising in solid-state physics and prove the exponential convergence of the carrier densities towards the inhomogeneous equilibria (theorem 4.4).

*Notation.* For the sake of clarity, we recall some basic notation and conventions in matrix calculus. The linear space of square *N*×*N* matrices is denoted by . Given a function with real coefficients, we denote its *j*th row by *A*^{j}(*x*) and its *k*th column by *A*_{k}(*x*). The components of an *N*-dimensional column vector will be denoted by *U*_{j}, *j*=1, …, *N*. Moreover, we adopt the following conventions.

The (vector valued) divergence operator acts

*row-wise*on matrices, i.e.The Jacobian matrix of a vector map

*has gradients on its rows*, i.e. .

The standard scalar product in will be denoted by *a*.*b*. We also recall the definition of the scalar product of two matrices as . For future use, we recall the formula with .

*Important note.* Owing to length limit, some of the proofs are omitted in the paper. They can be found in the electronic supplementary material.

## 2. Linear systems on bounded domains with constant equilibria

This section studies linear RD systems on bounded domains with zero outflow conditions, which feature unique, spatial homogeneous equilibrium solutions. However, we emphasize that this section is not intended to be merely an exercise on linear techniques. Our goal here is to capture (in the simplest possible setting) the system-related structural difficulties in determining the long-time asymptotics of linear(ized) RD systems by applying the entropy approach (as outlined in §1), which promises to be rather robustly extendable to nonlinear problems.

It is well known that RD systems may feature Turing instabilities, a phenomenon where stable steady states of a spatial homogeneous reaction system are rendered unstable when diffusion, sufficiently different in the components of the system, is added. However, as we shall see, Turing instabilities cannot occur in systems featuring an entropy functional that dissipates due to both reaction and diffusion until a unique entropy minimizing equilibrium state is reached.

Nevertheless, even for those ‘entropically stabilized’ systems, it remains a non-trivial question as to how convergence to equilibrium can be quantified in terms of reaction rates and diffusion constants, in particular, as in general, no maximum principle holds. Let us discuss an explicitly computable *example in one dimension*.

We consider reversible mass action-type reactions with for two species and described by the concentrations *a* and *b* on the interval *Ω*=[0,1], with homogeneous Neumann boundary conditionswhere *d*_{a} and *d*_{b} are the positive diffusion constants and *l* and *k* are the reaction rates. The above system has a unique positive, constant equilibrium state (*a*_{∞},*b*_{∞}) that balances the reaction (i.e. ) and satisfies the conservation of mass for given non-negative initial data.Linearizing around the equilibrium states and denoting *u*=*a*−*a*_{∞} and *v*=*b*−*b*_{∞} leads to the linear system(2.1)where . By rescaling we are able to set, without loss of generality, , i.e.which can be solved explicitly in terms of a Fourier series. By denoting the even mirrored extension of a function , the eigenproblem *ϕ*=λ_{k}*ϕ* on [0,1] with homogeneous Neumann boundary conditions transforms into the eigenproblem on [−1,1] with periodic boundary conditions. The resulting orthonormed eigensystem (*λ*_{0}=0 with eigenfunction and *λ*_{k}=−(*kπ*)^{2}, for all *k*≥1) Fourier decomposes the mirrored concentrations and , i.e.The individual Fourier modes yield the eigenvalues with the negative trace and the non-negative determinant . In particular, we have *Δ*>0 iff *k*≥1 and then . Moreover, we have , and the convergence to equilibrium is determined by the eigenvalueswhere *μ*_{1}(0) corresponds to the mass conservation law and the rate of convergence is given by the minimum of |*μ*_{2}(0)| and |*μ*_{1}(1)|. One can check whetherand the rate of convergence depends, in a non-trivial fashion, on the reaction parameters, the diffusion constants and a geometric parameter of the domain. Note that the value of |*μ*_{1}(1)| involves both reaction and diffusion parameters, except in the case of equal diffusion constants *d*_{a}=*d*_{b}=*d*, when |*μ*_{1}(1)|=*dπ*^{2} reflects only the diffusion. In the degenerate cases of one vanishing diffusivity *d*_{a}=0 or *d*_{b}=0, we have exponential decay with a rate involving the other positive diffusivity and the reaction.

### (a) Statement of the problem and remarks

We consider linear systems in the symmetrized form(2.2)with , *N*≥1, denoting the unknowns depending on time *t*≥0 and position within a bounded domain *Ω* with a sufficiently smooth (e.g. Lipschitz) boundary. We further prescribe the initial datum(2.3)and assume zero-flux boundary conditions(2.4)where denotes the outer normal unit vector at the boundary *∂Ω*. The structural assumptions on system (2.2) are the following.

*Non-degenerate symmetrizer matrix S.*The matrix is constant, symmetric and strictly positive definite. More precisely, there exists a constant>0 such that for all . For future use, we also define*s*as the maximum eigenvalue of*s**S*.*Positive semi-definite diffusion matrix*. The matrix is symmetric with eigenvalues for an integer*d*<*N*and for all . Moreover, has a constant eigenvector orthonormed matrix*F*, i.e.(2.5)*Negative semi-definite reaction matrix*. The matrix is symmetric with eigenvalues for an integer and for all . Moreover, admits a constant orthogonal eigenvector matrix*E*, i.e.(2.6)*Compatibility of reaction and diffusion matrices*. The following condition holds:(2.7)

A major subclass of the system (2.2), which we have in mind, is the linear RD system of the form(2.8)where the matrices *D*(*x*) and *R*(*x*) are, although not necessarily symmetric, *simultaneously symmetrizable*, i.e. there exists a non-degenerate, symmetric, positive definite, constant matrix such that *SD*(*x*) and *SR*(*x*) are symmetric. Multiplication of the system (2.8) by the symmetrizer matrix *S* yields the symmetrized system (2.2) with and (see, for instance, (2.1) in the above example with the symmetrizer matrix ). We remark that the structural assumptions and the theorems presented in this section can be stated more straightforwardly in the symmetrized version (2.2).

The set of assumptions (A1)–(A4) could be relaxed by dropping the symmetry assumptions on and and by instead requiring that the symmetric parts and satisfy all other hypotheses set up above for and , respectively. In terms of the original formulation (2.8) of a RD system, this means that one can avoid requiring that *D*(*x*) and *R*(*x*) are simultaneously symmetrizable as stated in remark 2.1 by instead supposing that(2.9)The exponential convergence statement in theorem 2.14 still holds valid under assumptions (2.9), with the only problem of a more involved notation. Since the existence of a simultaneous symmetrizer *S* is guaranteed in all the applied models motivating our theory, we shall simplify the treatment by requiring (A1)–(A4) instead of (2.9). We shall prove in proposition 2.11 that (2.9) are the minimal assumptions needed to have a reasonable energy structure for (2.8).

The set of conditions (A1)–(A4) and its generalization (2.9) share similarities with the structural condition required in the study of nonlinear systems of conservation laws (cf. Friedrichs 1954; Shizuta & Kawashima 1985; Kreiss & Lorenz 1989). In particular, condition (A4) is reminiscent of the so-called Kawashima condition, which is well known to ensure, for instance, convergence to equilibrium for hyperbolic–parabolic systems by ensuring that the normal modes corresponding to the nonlinear convection part do not ‘meet’ the degeneracy in the parabolic part.

In the case of non-zero integers *d* and *c* in assumptions (A2) and (A3), the diffusion and reaction parts separately are not able to drive the system towards a unique homogeneous equilibrium state. In particular, the degeneracy of the reaction matrix implies that a certain number of integral quantities are conserved with respect to time (see lemma 2.9). Therefore, the system needs a partial amount of diffusion to be driven towards a unique homogeneous equilibrium. More precisely, the diffusion operator should be non-degenerate along the zero eigenvectors of the reaction matrix. This is assured by the orthogonality condition (A4).

We remark that the assumption of a constant eigenvector matrix *E* in (A3) may be satisfied even in the case of a matrix with variable coefficients. Consider, as an example (with the notation introduced in (A1)–(A4)), the following 3×3 case:(2.10)with both *a*(*x*)−*b*(*x*) and *a*(*x*)+*b*(*x*) being uniformly negative.

The hypotheses of a constant eigenvector matrix *F* in (A2) is not needed in case the diffusion matrix *D*(*x*) is uniformly non-degenerate.

Before stating the theorems, we shall present two important examples of RD systems satisfying the set of assumptions (A1)–(A4). Both examples will be stated in the original form (2.8).

Let us consider the system(2.11)modelling transport and diffusion of charged particles in a semiconductor device combined with a recombination–generation mechanism called *band–trap capture and emission* (due to the presence of impurities; see Markowich *et al*. (1990, ch. 2.6) for further explanations). Here, *n* denotes the density of electrons, *p* denotes the density of holes and *n*_{tr} are the densities of occupied traps. The parameter *N*_{tr} represents the density of traps, therefore, *n*_{tr}<*N*_{tr}. For simplicity, we have neglected the effect of the self-consistent potential, therefore, coupling occurs only due to the recombination–generation terms. The equilibrium vectors satisfywhereas the uniqueness of the equilibrium state is achieved by imposingThe linearization of the system (2.11) around the unique equilibrium state gives the linear system(2.12)whereA symmetrizer for *R* is the diagonal matrixand the symmetrized reaction matrix readsThe matrix has the only zero eigenvector , which implies that assumption (A4) is satisfied. An elementary, but tedious, calculation shows that the two remaining eigenvalues *λ*_{2} and *λ*_{3} of are both strictly non-negative. Therefore, the symmetrized form of the linear system (2.12) satisfies all assumptions (A1)–(A4).

We consider a diffusive and reversible chemical reaction of the type on a bounded domain (cf. Desvillettes & Fellner 2008, submitted). Denote by *a*_{i}(*x*,*t*), *i*=1, …, 4, the concentration of the four reacting species and by *d*_{i}>0, *i*=1, …, 4, their respective diffusivity constants. Assuming mass action kinetics for the reactions, we obtain the systemwhere we have rescaled—without loss of generality—the reaction rates to 1. The stationary states consist of the unique set of positive constants that balance the reaction, i.e. , and satisfy the conservation laws of the systems. Linearization around those states produces the linearized reaction matrix *R* and the associated diagonal symmetrizer matrix *S*It is easy to check that all assumptions (A1)–(A4) are satisfied in this case.

### (b) Stationary states and their exponential stability

In this section, we prove exponential convergence to constant equilibrium states for linear systems on bounded domains. We shall use the symmetrized form (2.2) and work always in the framework set by assumptions (A1)–(A4). Note that the global existence of classical solutions for linear systems follows readily by standard fix-point and the Gronwall arguments (Rothe 1984; Amann 1985). A first natural problem is to detect and characterize stationary solutions to (2.2). In addressing these issues, we identify, in the following lemma, the eigenvectors with eigenvalue zero of the symmetrized reaction matrix as *conservation laws* of the system.

*Let U*(*x*,*t*) *be a classical solution to the system* (*2.2*). *Then, for every* *,* *the quantity*(2.13)*is conserved for all times t*≥0.

We compute the time derivative of *l*_{j}, for . Owing to (2.2) and (2.6), we havewhere we have used Green's theorem and the boundary condition (2.4). ▪

*For every fixed set of quantities* *,* *there exists a unique stationary solution U*^{∞} *to* (*2.2*) *such that*(2.14)*for all* . *Moreover, U*^{∞} *is constant and it satisfies* *with*

Let us multiply the stationary form of the system (2.2) by *U*^{∞}(*x*) and integrate over *Ω*. After integration by parts and in view of the boundary conditions, we obtainBut the term on the left-hand side is also non-negative due to assumption (A2),Consequently, both the above terms are equal to zero and we have . Let us then define(2.15)We haveSince *E* is non-degenerate, this implies . Therefore, the set of possible stationary states is the *d*-dimensional linear subspace of all *U*^{∞} with , which is equivalent to requiring(2.16)

On the other hand, condition (2.14) can be rewritten with (2.15) as(2.17)and, as the matrix *G* is constant and non-degenerate, there exists a unique *c*-dimensional constant vector satisfying (2.17). This completes the proof. ▪

For future use, we introduce, according to (2.15), the *normal-mode* vector variable(2.18)and the following *energy functional*(2.19)

In the following proposition, we will see that (2.19) is indeed a Lyapunov functional for (2.2), which will be used later on to prove, convergence to equilibrium for solutions to (2.2). Before that, let us detail further remark 2.2. Going back to the original form (2.8) of a linear RD system, we prove, in the following proposition, that a slightly weaker assumption than a simultaneous symmetrization is necessary to have a quadratic Lyapunov functional for (2.8). More precisely, we shall prove (under the assumption that and are constant, in order to avoid technicality) that conditions (2.9) are necessary for defined in (2.19) to be a Lyapunov functional for (2.8).

*Let R and D be constant matrices. Suppose the reaction matrix R has at least one left eigenvector V. Let S be a symmetric N×N matrix with constant coefficients. Then, the functional defined in* (*2.19*) *is a Lyapunov functional for* (*2.8*) *if and only if SD+DS and* *are non-negative definite*.

We first prove the ‘only if’ statement by contradiction. For a classical solution *U*(*t*) to (2.8), let us compute the evolution of the energy

*Case 1*. Suppose first that *SD*+*DS* is not non-negative definite. Then, there exists a vector such that . Let us define the vector function for a scalar smooth function *f*. Then, we haveWe then set . We consider the solution to (2.8) with initial datum *U*_{0}. At time *t*=0, we haveSince the functional is unbounded on the set , it is always possible to find a function such that the right-hand side above is strictly positive. Therefore, is increasing at *t*=0 and this proves the assertion.

*Case 2*. Suppose that is not non-positive. Let such that . Let us consider the constant initial datum . Trivially, one has at time *t*=0, which implies that the energy is increasing. Vice versa, it is straightforward to see that (*SD*+*DS*) and (*SR*+*R*^{T}*S*) being non-negative definite is a sufficient condition for to be non-increasing in time. ▪

Note that symmetrization is not the only way to show stability for a RD system. Consider, for instance, the two-dimensional systemwhere , *d*>0 andSuppose that the eigenvalues *r*_{11} and *r*_{22} of *R* are real, distinct and negative, and let *E* be the eigenvector matrix of *R*, such that . Let us also suppose that , ensuring that the symmetric part of of *R* has two eigenvalues of distinct sign (easy to check). In this example where the diffusion matrix is a multiple of the identity, a diagonalization method is more efficient than symmetrization. Indeed, the new variable *V*:=*EU* satisfies the uncoupled systemwhich implies that *V* decays to zero exponentially. The symmetrization method in the alternative sense (2.9) would also lead to the same result (by a more involved computation) by using the symmetrizing matrix , where *s*_{1} and *s*_{2} are the positive constants satisfying , which is possible if and only if as .

For the sake of completeness, we recall a well-known example of a linear 2×2 RD system featuring Turing instabilities (cf. Grindrod 1991), i.e. a system of the formwith *U*=(*u*,*v*), , *d*_{1} and *d*_{2} positive and having two distinct non-negative real eigenvalues. We shall point out that Turing instabilities can occur only if *R* is not symmetric and *d*_{1}≠*d*_{2}, which partly motivates the symmetrization setting (A1)–(A3) as a reasonable structural framework (besides the fact that such a framework fits certain target applications). Applying a Fourier transform gives the following system for the vector :The eigenvalues of *R* being negative implies, in particular, that *r*_{11}+*r*_{22}<0. By computing explicitly the eigenvalues of *B*(ξ), one easily realizes that exponential growth of one of the two normal modes on a closed non-empty *ξ* interval may occur only if , and this condition can be satisfied only if *d*_{1}≠*d*_{2} owing to the original condition *r*_{11}+*r*_{22}<0 on the trace of *R*. Imposing *R* to be symmetric implies , and therefore *r*_{11} and *r*_{22} both have a negative sign, which is in contradiction with the above condition (*r*_{11}*d*_{2}+*r*_{22}*d*_{1})>0.

We now state the main result of this section.

*Let U*(*x*,*t*) *be a classical solution to* (*2.2*) *with initial datum* *satisfying*(2.20)*for some given* . *Then, for all times t*≥0*,* *the following estimate holds*(2.21)*where* ; *is the optimal constant in the Poincaré inequality on the domain Ω;* *L and C*_{F} *are constants depending on the matrices E and F*, *respectively; the constants* *s**,* *μ**and* *r**are defined in assumptions* (*A1*)*–*(*A3*)*;* *and δ is an arbitrary positive small constant*, . *In particular**,* *exponentially fast in L*^{2} *as* .

In fact, exponential convergence in all *H*^{s} norms holds by interpolation between the above *L*^{2} decay and *H*^{s} *a priori* bounds, which are propagated by the linear system.

The proof relies on the estimate of the evolution of the relative energy functional . A standard symmetrization argument, the use of a vector-valued version of the Poincaré inequality and a careful decomposition of the quadratic form via the normal modes are its main ingredients. The complete version of the proof can be found in the electronic supplementary material. ▪

## 3. Confined linear systems on with integrable equilibria

In this section, we deal with linear reaction–convection–diffusion systems posed on the whole of . The main difficulty with respect to §2 consists of the presence of a *confinement matrix* and the spatial inhomogeneity of the integrable stationary states. As we shall see in this case, a theory applicable to a significant class of examples can only be developed when the symmetrizer for the reaction matrix has variable coefficients.

Instead, the systems under consideration in this section read as(3.1)with , *d*≥1, *t*≥0 and . We prescribe the initial datum(3.2)As explained in example 3.3 below, the system (3.1) generalizes a model of *N* reacting–diffusing species convected by *N* external potentials. In the sequel, we shall first develop a general theory for *N*×*N* systems including the exponential asymptotic stability result in theorem 3.6. Then, we shall focus on the 2×2 case, where improved results can be obtained under less restrictive assumptions.

### (a) General theory for N species

In this section, we shall study the system (3.1) coupled with the initial datum (3.2), posed on the whole of . We shall assume the structural hypotheses as follows.

*Essentially convex and diagonal confinement.*The confinement matrix is of the form and such that there exist two functions and a constant*K*>0 such that and, for all*j*'s,(3.3)Moreover, we shall assume that the functions are*L*^{∞}perturbations of uniformly convex functions for all*j*=1, …,*N*.*Positive definite diffusion matrix.*The diffusion matrix*D*(*x*) is symmetric and positive definite. Moreover, the following inequality is satisfied:(3.4)for a certain positive constant*d*_{0}.*Confinement compatible reaction matrix.*The confinement matrix*S*(*x*) is a symmetrizer for the reaction matrix*R*(*x*), i.e.Moreover, the symmetrized reaction matrix has eigenvaluessatisfying , for all and there exist two positive constants*C*_{1}and*C*_{2}such that(3.5)for all and for all .*Bounded eigenvector matrix E of the reaction matrix*. The orthogonal eigenvector matrix*E*(*x*) of satisfiesMoreover,*E*(*x*) has uniformly bounded coefficients , , such that for all*i*and*j*.*Conservation laws.*The reaction matrix*R*(*x*) has*d*constant left zero eigenvectors, i.e. there exist such that

Unlike in §2, in the present case, we are not including possible degeneracy of the diffusion matrix *D*(*x*) in the assumptions. This choice is motivated mostly by the sake of avoiding further restrictions in the set of assumptions and further complications in the notation. However, we are confident that it would be possible to allow degeneracies of *D*(*x*) and still achieve a unique, asymptotically stable equilibrium by some sort of orthogonality condition in the spirit of (2.7). This appears clearly in the estimate (3.8) of the electronic supplementary material, where the last term on the right-hand side (which has a non-negative sign and it is not needed in the estimate) would not appear in the case of a suitably chosen degeneracy of *D*(*x*).

The convexity assumption in (B1) assures that each one of the *L*^{1} densities , *j*=1, …, *N*, satisfies the *weighted Poincaré inequality*for all such that . Weighted Poincaré inequalities on the whole space are the subject of a very wide literature (e.g. Arnold *et al*. 2001; Arnold & Dolbeault 2005; Arnold *et al*. 2007). More refined conditions could be prescribed involving the diffusion matrix *D*(*x*) (see the above references).

The above assumptions generalize the main class of examples we have in mind, arising from the applications, i.e. the following linear system of *N* convection–reaction–diffusion equations with inhomogeneous equilibria.

N×N reaction–drift–diffusion system on . Let us consider the system(3.6)where , , *t*≥0; *R*^{j} denotes the *j*th row of a reaction matrix *R*; and *V*_{j} is a *confining potential* acting on the *j*th species *u*_{j}. A simple computation shows that the system (3.6) can be written in the above form (3.1) withand *R*(*x*) is any reaction matrix such that *S*(*x*)*R*(*x*) is symmetric. Assumption (B1) can be fulfilled by assuming that all potentials *V*_{j}(*x*) are all continuous and convex with equal asymptotic behaviour of the tails as . However, this assumption on the tails is often too restrictive. Section 3*b* will show how to relax it in a 2×2 example system.

For future use, we write the system (3.1) in the symmetrized form(3.7)

The proof of the following lemma is straightforward (see the proof of lemma 2.9).

*Let U*(*x*, *t*) *be a classical solution to the system* (*3.1*) *such that U decays rapidly at* . *Then, for every* *,* *the quantity**is conserved for all times t*≥0.

*For every fixed set of real numbers* *,* *there exists a unique stationary solution U*^{∞}(*x*) *to* (*3.1*) *such that*(3.8)*for all* . *Moreover, U*^{∞}(*x*) *has the form*(3.9)*for a certain constant vector* .

The proof of this proposition follows a strategy similar to that proposed in the proof of proposition 2.10. It can be found in the electronic supplementary material. ▪

We introduce the following variables:Similar to §2, we shall consider the energy functional

*Let U*(*x*, *t*) *be a classical solution to the system* (*3.1*) *with initial datum U*_{0}*,* *such that U decays rapidly at* . *Suppose*(3.10)*for fixed quantities* . *Let the stationary state U*^{∞} *be uniquely determined by the l*_{j} *quantities as in* *proposition 3.5*. *Then, there exists a fixed constant* >0 (*depending on the structural assumptions* (*B1*)–(*B5*)) *such that*(3.11)*for all t*≥0.

The proof of the present theorem follows a strategy similar to that of theorem 2.14. It can be found in the electronic supplementary material. ▪

### (b) A typical 2×2 example case

In this section, we shall detail a 2×2 whole-space drift–diffusion–recombination system, which is, in fact, the linearization of the nonlinear semiconductor model (4.1) to be studied in §4.(3.12)where(3.13)which implies, in particular, that . As stated in remark 3.2 in §3*a*, the set of assumptions (3.13) ensures the validity of a suitable weighted Poincaré inequality. The reaction rate shall denote a continuous function (further assumptions on *F* will be specified below) and is typically thought to be of *Shockley–Read–Hall* type(3.14)for some positive, bounded below functions .

With the vector variable , the above system (3.12) is of the form (3.1), with(3.15)and the matrix *R*(*x*)(3.16)For the present 2×2 system, we are able to follow a different line of the estimates compared with §3*a* and obtain an improved result, under less restrictive assumptions than the set (B1)–(B5).

Note that the considered example is a good representant for a general two-species theory. Indeed, one can easily verify that the set of assumptions on the reaction matrix in (B1)–(B5) forces the above choice (3.16) for *R* in the 2×2 case, except that one could have a left eigenvector corresponding to the zero eigenvalue different from (1, −1). A more general left eigenvector can be obtained by multiplying one of the two equations in (3.12) by an arbitrary constant, which does not imply any significant change in the computations below. Moreover, it is easily seen that, in order to have the matrix *SR* symmetric, *R* must be of the form (3.16) up to a scalar multiplier. Our result in the 2×2 case will require less restrictive assumptions on the scalar multiplier *F* than those inferred by the set of assumptions (B1)–(B5). Finally, we could easily generalize the following computations to the case of a more general *D* satisfying assumption (B1). However, in order to simplify the presentation, we shall consider *D*=*S*^{−1} as in (3.15).

We shall use the entropy (or free energy) functionalwhich coincides with in the vector notation, consistently with the strategy developed in §§1 and 2. After introducing the variables and and multiplying equation (3.12) with the vector (*z*_{1}, *z*_{2}), we obtainwithand the measures and . Denoting further the normalized measureswe estimate the Fisher information terms using a weighted Poincaré inequality (cf. Arnold *et al*. 2001)The constants *P*_{1} and *P*_{2} are the whole-space Poincaré constants with respect to *ξ*_{1} and *ξ*_{2}, respectively; therefore, they depend on *V*_{1} and *V*_{2}. Then, for a suitable constant *C*>0, we are looking for the following entropy–entropy dissipation estimate:under the constraint for the conservation of mass(3.17)

*Consider measurable functions z*_{i}, *i*=1,2*, such that* (*3.17*) *holds*. *Let* *be integrable with respect to the measure* *and satisfy*(3.18)*Then,**holds, provided that*(3.19)*with K*_{1} *defined in* (*3.22*) *and where* 0<*ϵ*<1 *can be chosen in order to maximize the constant K*.

In contrast to theorem 3.6, the strategy of the proof presented here uses the entropy dissipation coming from the reaction to construct a lower bound of the entropy dissipation in terms of the relative entropy. This is done crucially for the special case of constant *z*_{i} using the conservation law and the integrability of *F* w.r.t. the measure . To obtain the entropy–entropy dissipation estimate for general *z*_{i}, we employ an ansatz around constant states and take a sufficiently small fraction of the reaction–dissipation such that the remainders are controlled by the dissipation coming from the diffusion. To control the remainders, the bounds (3.18) are required naturally. Nevertheless, for general systems such as those stated in theorem 3.6, this approach seems at least clumsy to formalize.

*Step 1*. We consider the particular situation when *z*_{1} and *z*_{2} are constant. We then calculate the relative entropy as(3.20)while for the entropy dissipation we have(3.21)Thus, using the conservation law to express one in terms of the other, we find, by comparing (3.20) with (3.21) (symmetrized in terms of and ), that(3.22)

*Step 2*. We employ the ansatz , *i*=1,2 with . Expanding the relative entropy yields(3.23)On the other hand, expanding the entropy dissipation, we estimate the contribution coming from the reaction using Young's inequality for *ϵ*<1Then, using (3.18), we continue to estimate again with Young's inequality*Step 3*. By the previous steps, we may now estimate for a constant *K*<1and the statement of lemma 3.7 follows. ▪

We have therefore proven the following.

*Let n*_{∞}*,* *p*_{∞}*,* *V*_{1} *and V*_{2} *satisfy* (*3.13*) *and let assumption* (*3.18*) *on F hold. Let* *be a classical solution to the system* (*3.12*) *with initial datum* (*u*_{0}, *v*_{0}) *having finite energy* . *Then,**with the constants K**,* *K*_{1} *and ϵ as in* *lemma 3.7*.

Comparing the result of theorem 3.9 with that in theorem 3.6, we find the improvements as follows.

The exponential rate of convergence is much easier to compute.

We do not need the two potentials

*V*_{1}and*V*_{2}to have equivalent behaviour of the tales at , as was implicitly required in (B1) (cf. example 3.3).The assumptions (3.18) on

*F*are less restrictive than those implicitly inferred by (B3). In particular, (3.18) match with the Shockley–Read–Hall reaction rate (3.14).

## 4. A reaction–diffusion model with nonlinear reaction

In this section, we study the nonlinear model system arising in semiconductor and plasma physics(4.1)where *n* and *p* model two species of charged particles subject to confinement and to a recombination–generation mechanism *R*(*n*, *p*). We suppose non-negative initial dataand the following assumptions.

The confining potentials

*V*_{n}and*V*_{p}satisfyMoreover, is finite for*i*∈{*n*,*p*}, and we define and introduce the related measures .The recombination–generation term is of the form for a constant

*δ*>0, which—without loss of generality—shall be rescaled as*δ*=1. The scalar function is assumed to be such that(4.2)for constants*A*_{1}>0,*A*_{2},*A*_{3}≥0, which includes the typical Shockley–Read–Hall form (3.14).

The nonlinear system (4.1) models recombination–generation mechanisms in semiconductor theory when neglecting the influence of the self-consistent potential. The physical much more relevant and mathematically more interesting case of system (4.1) coupled to a Poisson equation for a self-consistent field is nevertheless beyond the scope of this paper and will be the subject of an upcoming article (Di Francesco *et al*. in preparation). For the existence of global solutions of the system (4.1) (with self-consistent potential), we refer to (Wu *et al*. 2008). We recall that the initial mass is conserved for all *t*>0(4.3)Given *M* as fixed, the equilibrium *n*=*n*_{∞}, *p*=*p*_{∞} is uniquely determined by(4.4)

The *relative entropy E*=*E*(*n*, *p*) of the system (4.1) with respect to the equilibrium (4.4),(4.5)dissipates (i.e. ) with the entropy dissipation(4.6)Note that the first two integrals in the entropy dissipation quantify how drift and diffusion tend to make the ratios *n*/*n*_{∞} and *p*/*p*_{∞} constant, while the third reaction-type integral drives their product to be equal to 1. Hence, the conservation of mass (4.3) determines that the entropy dissipation vanishes only for the equilibrium (4.4).

In order to prove the main result of this section, we need a technical result concerning with a uniform *L*^{∞} bound for the solution (*n*(*t*), *p*(*t*)) to (4.1). This issue is addressed in the electronic supplementary material.

### (a) Exponential convergence to equilibrium

In the following, we show exponential convergence (with constants that can all be made explicit) towards the unique equilibrium states *n*_{∞} and *p*_{∞} as defined in (4.4). In addition to assumptions (NL1) and (NL2), we will suppose that:

The confining potentials are equal and is—without loss of generality—normalized with .

There exists a constant lower bound , and moreover due to the bounds provided in the electronic supplementary material, .

*Notation.* We rewrite the drift–diffusion fluxes *J*_{n} and *J*_{p} in terms of their stationary state *μ*, i.e.and *J*_{n} and *J*_{p} vanish iff *n* and *p* are proportional to *μ*. Moreover, we introduce the shortcutsand observe that and due to .

*Let n* and *p be the non-negative functions in* *satisfying the conservation law* *,* *as given in* (*4.3*). *Suppose assumptions* (*NL1*)*,* (*NL3*) *and* (*NL4*) *hold*. *Then, the following inequality holds for a constant K depending only on the stated quantities*(4.7)

*Step 1*. Upper bound for the relative entropy. We recall the logarithmic Sobolev inequality for the normalized reference measure d*μ*=*μ* d*x* corresponding to the strictly convex confining potential *μ*=*e*^{−V} (see Holley & Stroock (1987) to also permit *L*^{∞}-bounded perturbations)with *σ*_{n} as in assumption (NL1) and an analogue inequality holds for *p* with respect to *σ*_{p}. Thus, splitting the entropy density as , we obtainEmploying an elementary inequality for (Desvillettes & Fellner 2008) yields further that(4.8)Note that as a particular benefit of employing the logarithmic Sobolev inequality individually to the two species, we may estimate the *L* log *L* structure of the entropy on the level of the integrals and , for which bounds are sufficient. Alternatively, estimating directly the entropy density would require and bounds, and exponential convergence can be obtained in the case of slowly growing *L*^{∞} *a priori* estimates (Desvillettes & Fellner 2008).

*Step 2*. Mass exchange due to recombination and the conservation law. We quantify the recombination mass exchange between *n* and *p* under the constraint of the conservation law (4.3), i.e. in terms of a two-dimensional functional inequality involving the integrated quantitiesPrecisely, we show that for all , there exists a constant *C* such that(4.9)holds under the conservation law constraint(4.10)We remark that inequality (4.9) (together with (4.12) below) yields the entropy–entropy dissipation estimate (4.7) in the special situations when *n*/*μ* and *p*/*μ* are constant.

To prove (4.9) we factorize, for instance, with respect to *y*_{2}using the constraint (4.10). Hence, again using (4.10), the inequality (4.9) rewrites aswhich holds for a constant *C* as the fractionis bounded on [0,∞)^{2}, as can be checked by distinguishing the cases . Thus, altogether from (4.8) and (4.9), we have that(4.11)

*Step 3*. Lower bound for the entropy dissipation. Denoting the third reactive term of the entropy dissipation (4.6) *D*_{r}, we apply the elementary inequality , assumption (NL4) and Jensen's inequality to estimate(4.12)Hence, comparing (4.12) with (4.11) we are—qualitatively speaking—left to interchange the square root with the integrals. We therefore estimated some fraction (leaving the rest aside for the logarithmic Sobolev inequality) of the drift–diffusion terms of the entropy dissipation (4.6) by Poincaré's inequalityand analogue with *p* and *δ*_{p}. Expanding (4.12) in terms of *δ*_{n} and *δ*_{p}, we exploit that the first-order integrals vanish and estimate with Young's inequalitywhere by Cauchy–Schwartz and analogue for *p*. We remark that, since the entropy dissipation vanishes at the equilibrium in the second order in terms of *δ*_{n} and *δ*_{p}, the cancellation of the first-order integrals (in general ) is crucial for the proof and precisely the point where we are forced to make assumption (NL3).

To quantify what qualitatively was called interchanging square roots with integration, we observe (and again analogue for *p*) that(4.13)and is unbounded if and only if vanishes. We are therefore led to distinguish two cases. At first, for an *ϵ*>0 to be set below, we suppose the case and , and expand in terms of (4.13)and the first term on the right-hand side corresponds to (4.11). Thus, we obtain the entropy–entropy dissipation estimate (4.7) by choosing a constant large enough such that the remainder terms of order and are compensated by the Poincaré lower bounds of the drift–diffusion dissipation terms (after leaving some fraction aside to employ the logarithmic Sobolev inequality).

The other cases, e.g. , are in terms of *n*: in these far-from-equilibrium cases, we chose and find therefore thatHence, since the second term on the right-hand side of (4.11) is trivially bounded in terms of a constant , we obtain the stated entropy dissipation provided for a constant large enough. ▪

The above proof of lemma 4.1 obviously does not yield a sharp convergence rate. Nevertheless, the proof requires only natural *L*^{1} bounds (which follow very generally from the entropy decay), a non-degeneracy assumption on the reaction rate (NL4), which, as a consequence of the homogeneity of the nonlinear reaction term *R*(*n*, *p*), includes also the confining measure *μ* and suitable confinement assumptions (NL1) and (NL3). It remains an open problem to prove exponential (or algebraic) convergence without (NL3) or if (NL4) would be relaxed to *F*≥*C*_{F}>0.

An entropy–entropy dissipation estimate for the system (4.1) coupled to a self-consistent Poisson equation is the subject of a forthcoming paper (Di Francesco *et al*. in preparation).

Let *n and p be the solutions to the system* (*4.1*) *subject to non-negative initial data n*_{I} and *p*_{I}. *Suppose assumptions* (*NL1*)–(*NL4*) *hold. Then, the solution converges exponentially* (*with explicitly computable constants*) *to the unique equilibrium state* (*4.4*) *and the following estimate holds*:(4.14)*with a constant K depending explicitly on the initial masses* and *,* *the initial entropy E*(*n*_{I}, *p*_{I})*,* *the uniform convexity constant of the confinement potential V and the lower bound in* (*NL4*).

By the uniform *L*^{1} bound provided in the electronic supplementary material and the entropy–entropy dissipation estimate lemma 4.1, the results follow directly by integratingand the Csizsár–Kullback inequality(4.15) ▪

## Acknowledgments

This article was supported by the KAUST Investigator Award 2008 and the Wolfson Research Merit Award (Royal Society) of P.A.M. M.D.F. was partially supported by the Italian research project ‘Modelli iperbolici non lineari in fluido dinamica’ (INdAM GNAMPA 2008).

## Footnotes

↵† On leave of absence from the Faculty of Mathematics, University of Vienna, Nordbertstr. 15, 1090 Wien, Austria.

- Received May 27, 2008.
- Accepted July 17, 2008.

- © 2008 The Royal Society