A kinetic reaction model: decay to equilibrium and macroscopic limit

We propose a kinetic relaxation-model to describe a generation-recombination reaction of two species. The decay to equilibrium is studied by two recent methods for proving hypocoercivity of the linearized equations. Exponential decay of small perturbations can be shown for the full nonlinear problem. The macroscopic/fast-reaction limit is derived rigorously employing entropy decay, resulting in a nonlinear diffusion equation for the difference of the position densities.


Introduction.
1.1. The model. We consider the system where f and g depend on position x ∈ T 3 , the three dimensional torus of volume one, on velocity v ∈ R 3 , and on time t ≥ 0. They represent the phase space densities of chemical reactants A and B, which are produced (with nonnegative velocity profiles χ 1 and χ 2 , respectively) by the decomposition of a substance C. The density of the substance C is not subject of our study and is assumed to be fixed. On the other hand the substances A and B can recombine to form C and thus be eliminated from our system. Similar models have been used for generation and recombination of electron-hole pairs in semiconductors [5,6]. The probability of the reaction is depending on the position density of the reaction partner. We consider the system (1) subject to initial conditions
The last line will be needed for an L 2 averaging lemma with the weight 1/χ j . The largest value to be expected for the exponent is θ = 1, which is achieved, e.g., for Gaussian distributions, the prototypical examples for the χ j , but more generally also for χ j (v) ≤ c(1 + |v| 2 ) −k with k > 1.
Note that, at least formally, the mass difference is conserved: as can be seen by subtraction of the two equations and subsequent integration. This is to be expected since by the reaction molecules of A and B are created and destroyed pairwise. We introduce the unique constant ρ ∞ > 0, such that and expect convergence as t → ∞ of solutions of (1), (2) to the steady state satisfying ρ f∞ = ρ ∞ and ρ g∞ = 1/ρ ∞ . This is supported by the decay properties of the entropy functional which decreases as long as (f, g) is different from (ρ(x)χ 1 (v), where the superscript denotes evaluation at v . Among these functions (f ∞ , g ∞ ) is the only solution of (1) with the same mass difference as the initial data. Spectral stability of the equilibrium will be investigated by linearization: where for simplicity the perturbations have again been denoted by f and g, now satisfying Rigorous results for kinetic equations for chemically reacting species with nonlinear reaction models are scarce in the literature. An example is an existence result [14] for a model with quadratic reaction terms under rather weak natural assumptions on the initial data, where stability is based on entropy decay. Since existence of solutions is not the main focus of this work, we shall make rather strong assumptions on the data, with the consequence that the existence and uniqueness proof in the following section is based on weighted L ∞ estimates and rather straightforward.
The analysis of the decay to equilibrium is complicated by the fact that the entropy dissipation (6) vanishes not only for the equilibrium, but on a larger set of local equilibria. If exponential decay to equilibrium can still be proven, the system is called hypocoercive [16]. The authors have been involved in the development of two different abstract procedures for proving hypocoercivity for linear equations [9,13], both based on the construction of suitable Lyapunov functionals (or modified entropies), whose dissipation functionals have appropriate coercivity properties. The method of [9] is based on a slightly tilted, weighted L 2 -norm, while [13] works in a H 1 setting and can be extended to higher regularity. In Section 2 we show that both methods are applicable to a linearized version of (1). Since the estimates of the existence result already provide neutral stability of the equilibrium, the decay results can be extended to a local asymptotic stability result with exponential convergence for the full nonlinear model. The decay rates proven by both methods can, in principle, be computed explicitly. Complete formulas would however be rather complicated, whence we did not attempt a comparison. An essential difference between the methods is the weaker assumptions on initial data in [9]. On the other hand, the method of [13] has the potential to provide strong convergence properties including derivatives.
In Section 1.3 the macroscopic/fast-reaction limit is carried out formally, leading to a nonlinear diffusion equation for the difference of the position densities of the reactants. Similar results have been derived for reaction-diffusion systems [2,7,12] and for coagulation-fragmentation models [3,4]. A rigorous justification of the limit is the subject of Section 3. It is based on an analysis of the entropy dissipation functional (6) and adapts the procedure of [15], where compactness is obtained from an averaging lemma in weighted L 2 -spaces. We prove a slightly generalized version compared to [15].
We note that with the torus we chose the simplest geometric setting. Natural modifications include bounded domains with specular reflection boundary conditions or whole space problems with confining potentials. We conjecture that our results can be extended to these situations, however with considerably more technical effort for the latter (see, e.g., the hypocoercivity results with confining potentials in [9]).

1.2.
Existence of solutions. The entropy decay relation (6) would suggest an existence result for initial data with bounded entropy. Such a result for a similar problem has been proven in [14]. The main ingredients are entropy inequalities, weak L 1 compactness and velocity averaging. These ideas might be transferable to our situation. However, for our purposes we need more information on the solutions. Under stronger assumptions on the initial data, a global existence result can be proved easily.

LUKAS NEUMANN AND CHRISTIAN SCHMEISER
Then the initial value problem (1), (2) has a unique global mild solution (f, Proof. The mild formulation of the equation for f is given by and an analogous equation holds for g. It is easily seen that the set of (f, g) defined by the estimates of the theorem is mapped into itself by the right hand sides. This provides the a priori estimate needed for the continuation of a local solution constructed by Picard iteration.
1.3. Formal macroscopic limit. In this section we formally derive a macroscopic limit of (1). The limit will be made rigorous in Section 3. Since by (3) the mean velocities of the equilibrium distributions vanish, we adopt a diffusive (or parabolic) scaling t → t/ε 2 and x → x/ε and derive We substitute the ansatz Balancing the leading order terms gives ρ g0 f 0 = χ 1 and ρ f0 g 0 = χ 2 . This is equivalent to the existence of ρ 0 (x, t) such that

Due to (3) the solvability condition
where the second terms on the right hand side constitute the general solution of the homogeneous problem. Note that ρ f1 = ρ 1 , ρ g1 = −ρ 1 /ρ 2 0 . Now we substitute this into the limit of the conservation equation where we have introduced the positive definite symmetric matrices This can be written as the nonlinear diffusion equation for the new unknown where we have introduced the diffusion matrix The unknown m is the zeroth order approximation of the difference of the position densities of f and g.

2.
Long time properties. In this section we study decay to equilibrium for solutions of (1) and of the linearized problem (7), (8). In order to estimate the decay towards the equilibrium quantitatively we introduce the L 2 scalar product, weighted with the steady state measure, Throughout this article we denote by · w the norm induced by this scalar product. At this point we introduce some notation for function spaces. The subscripts t, x, v will indicate spaces of functions of t ∈ (0, ∞), x ∈ T 3 , v ∈ R 3 , and the superscripts χ 1 , χ 2 will indicate the weights 1/χ 1 , 1/χ 2 . The superscript w (meant to abbreviate w = (χ 1 , χ 2 )) will be used for spaces of pairs of functions, where the weight 1/χ 1 is used for the first component and 1/χ 2 for the second. Thus, ·, · is a scalar product on . The orthogonal projection onto the null space of the linearized collision operator that is the space of local equilibria, is given by Straightforward computations yield called microscopic coercivity in [9]. This gives a quantitative estimate of the decay towards the local equilibrium introduced by the linearized collision operator. In the spatially homogeneous situation such an estimate is enough to prove exponential decay to equilibrium for the linearized equation.
For spatially non homogeneous situations we expect the densities to become constant as these are the only local equilibria that also annihilate the transport part of the equation. The complete relaxation mechanism can be seen as a combination of local relaxation in the velocity direction by the collision operator and an interplay between mixing by the transport operator and confinement in our bounded spatial domain. In the following we will study two recent methods ( [9,8] and [13]) to estimate the decay in the spatially non-homogeneous situation. Both rely on properties of the linearized collision operator that we will verify in the sequel. The microscopic coercivity property (13) is needed in both approaches.
2.1. Coercivity in a weighted L 2 -space. In this section we apply the abstract convergence theory of [9]. This approach uses a modified L 2 entropy functional to quantify the mixing effect of the transport. We start by writing the linearized equation (7) in the abstract form with F = (f, g) T , with the transport operator and with the linearized collision operator L given in (11). The following result for abstract linear ODEs in Hilbert spaces has been proven in [9].
Then there exist positive constants C and λ, depending on λ m , λ M , C M , such that the semigroup generated by L − T satisfies Remark 1. The approach of [9] relies on the Lyapunov (or modified entropy) functional with a small positive constant δ. Its time derivative along solutions of (14) is The first term on the right hand side suggests that the microscopic coercivity assumption (16) is one of the necessary ingredients. Since the operator ATΠ can be interpreted as the application of the map z → z 1+z to (TΠ) * TΠ, the second condition (17), called macroscopic coercivity, is also plausible. As a consequence of (16) and (17), the sum of the first two terms in the entropy dissipation (20) is coercive: Theorem 2.1 will be applied with L given by (11), T given by (15), and with H = L 2,w x,v equipped with the scalar product (10) and the norm · w . The projection Π is then given by (12). It is easily seen that L is symmetric and T is antisymmetric. Because of (13), i.e. λ m = min{ρ ∞ , 1/ρ ∞ }, it remains to verify (17)-(19). A straightforward calculation shows that (17) is equivalent to (8). Thus, by the positive definiteness and since the application of Π involves an integration with respect to v, the assumptions (3) imply (18). The linear collision operator L is easily seen to be bounded. For the verification of (19) it is thus sufficient to prove the boundedness of AT or, equivalently, of its adjoint (AT) * = −TA * = −T 2 Π(1 + (TΠ) * TΠ) −1 . The equation

LUKAS NEUMANN AND CHRISTIAN SCHMEISER
The norm T 2 ΠG w is equivalent to the L 2 (T 3 )-norm of ∇ 2 x u G , whose boundedness in terms of the L 2 (T 3 )-norm of u F (and therefore in terms of F w ) is a consequence of elliptic regularity. This proves (19) and completes the verification of the assumptions of Theorem 2.1.
Since the maximum principle estimates of Theorem 1.1 already imply a stability (but not asymptotic stability) result for the nonlinear problem, the decay result can be extended to a local result for the nonlinear case by the same method.
with positive constants C and λ.
Proof. We start by writing the problem in terms of the unknown F = (f −ρ ∞ χ 1 , g− χ 2 /ρ ∞ ) T . Then we proceed as above producing the entropy decay relation (20), however with LF replaced by This difference is given by and, by Theorem 1.1, ρ f is close to ρ ∞ and ρ g close to 1/ρ ∞ , yielding with a small constant γ. Therefore the entropy dissipation for the nonlinear equation is a small perturbation of the entropy dissipation of the linearized problem, which does not destroy its coercivity.

2.2.
Coercivity in H 1 . When studying coercivity of the collision operator we saw that in L 2 w the operator provides coercivity only with respect to the velocity distribution. The strategy in [13] is to transfer some of this dissipation effect in the velocity to the spacial variable by using the mixing effect of the transport operator. This method was mainly inspired by discussions with C. Villani (see [16]) and results by Y. Guo (see for example [11] or the references in [13]).
The regularizing effect of the transport operator can be quantified by looking at the time evolution of mixed derivative terms. The weights are in principle the same as in the previous section, however on the level of derivatives it is much simpler to use a standard L 2 -norm and rescale the equation to get rid of the weights. This also leads to a much cleaner notation and facilitates frequent applications of partial integration. To make the distinction from the weighted norms more obvious we denote the standard L 2 -norm in x and v by · L 2 .
We rescale the linearized equation by replacing and derive the rescaled equation We denote the standard L 2 scalar product in x and v and acting component-wise on the entries of F = (f, g) T by (·, ·) and write the orthogonal (with respect to the L 2 scalar product) projection on the kernel ofL as The rescaled collision operator naturally inherits the microscopic coercivity estimate with the same constant. Now we are in the position to calculate the time evolution of derivatives. For any ν > 0 we have, using the divergence theorem and the periodic boundary conditions in x, where the gradients have to be applied component-wise to F = (f, g) T . We used integration by parts in the equality and Cauchy-Schwarz as well as the fact thatL acts locally in x, in the estimate. A simple computation shows thatL is bounded in L 2 and thus the second term in the right hand side of (23) can be compensated by which is again a consequence of locality ofL and the microscopic coercivity. Thus by combining (23) and (24) the transport operator provides some control of the x-derivatives, if the last term in (23) can be controlled by coercivity properties of the collision operator in the velocity derivatives (cf. (28)). We introduce the Sobolev space H 1 with the (unweighted) norm where again the norms as well as the derivatives are understood component-wise in the entries of F = (f, g) T .
Following [13], the Lyapunov functional is constructed, which is, for δ small enough, equivalent to the square of the H 1 -norm but uses the mixing of derivatives to get a coercivity estimate also for the spatially non-homogeneous situation. Exponential decay of H 1 [F ] is then derived from and convergence in H 1 follows from the equivalence mentioned above. Since the essential mixing effect of the transport operator is quantified by means of (23) this approach is suitable only for proving coercivity in spaces of differentiable functions. On the one hand this is a restriction of the method, which, on the other hand, can also be employed to prove convergence of higher derivatives. This allows to use embedding theorems and lends itself to studying nonlinear equations in the perturbative regime. In the model we study here, however, we have good a priory bounds already and thus there is no need to go beyond first order derivatives as we will see in the proof of Theorem 2.5.
To control remainder terms in the time evolution of H 1 [F ] and quantify the coercivity of the linearized reaction operator on the level of velocity-derivatives structural assumptions are needed. We will verify these in the sequel. To get decay of the velocity derivatives it is essential that the linearized reaction operator can be split,L = K − Λ, into a "loss" part which is (in our case trivially) coercive, and a "gain" part which is regularizing in v as long as √ χ 1 and √ χ 2 are regular. Indeed, for ∇ v acting component-wise in the two functions, Cauchy-Schwarz together with Young inequality lead to where C depends on Since the method relies on H 1 type estimates, coercivity in H 1 of the loss part is necessary. In our case the same estimate as in (26) results in where in general negative terms of lower order derivatives are allowed but not needed in this case. Properties (26)-(28) together with (22) and the boundedness ofL ensure that we can use the main theorem from [13] to derive convergence for the linearized problem. If resorting to the equation before rescaling the result can be formulated as follows: Then the solution (f, g) of the linearized problem (7) subject to initial conditions f (t = 0) = f I , g(t = 0) = g I exists globally and converges exponentially to the equilibrium distribution. For T 3 R 3 (f I − g I )dv dx = 0 the equilibrium is zero and we have f / √ χ 1 (·, ·, t) 2 H 1 (dv dx) + g/ √ χ 2 (·, ·, t) 2 where C and τ depend on the L 2 bound onL and the constants in the estimates (22), (26)-(28). The constant C also depends on the norm of the initial data.
Convergence in higher order Sobolev spaces can be derived straightforwardly provided the estimates (27), (28) can be generalized to higher order derivatives, as is easily verified for our model. This feature is useful mainly in applying the results to the nonlinear system in a perturbative setting. Control of the bilinear contribution in the interaction is given by applying the chain rule, the Hölder inequality, and using Sobolev embedding to lower the exponents in the norm to two again (see [13] for details).
Here however we want to give a stronger result -in the sense that less regularity is necessary -by using the a priori bounds of Theorem 1.1.
Theorem 2.5. Let the assumptions of Theorem (2.4) hold. Moreover let f I and g I satisfy the assumptions of Theorem 1.1 with γ 1 and γ 2 small enough. Then the solution of (1) with initial data (f I , g I ) exists globally in time and converges exponentially to the unique stationary state, more precisely with positive constants C and τ showing the same dependencies as in Theorem 2.4 and τ also depending on γ 1 and γ 2 .
Proof. After rescaling the new unknown f obeys the bound (cf. Theorem 1.1) and in the same way we obtain for the rescaled g F = (f, g) T is the solution of the rescaled equation with the bilinear contribution given by Now using (25) the time evolution of the modified H 1 norm can be estimated by for some positive constant C representing the equivalence of H 1 and the squared H 1 -norm. Using the multiplicative structure of B(F, F ) we infer from the bounds (29), (30), the normalization of χ 1 and χ 2 as well as Cauchy-Schwarz, that with a constant γ that becomes arbitrarily small as γ 1 and γ 2 decrease. Together with (31) and the equivalence of H 1 with the square of the H 1 -norm this proves exponential decay of H 1 [F ]. Rewriting the result in the original scaling leads to the statement of the theorem.
Proof. We represent f by the Fourier transform with respect to t and by the Fourier series with respect to x: with the lattice T 3 * dual to the torus T 3 , implying For each λ > 0, we introduce a smooth, nonnegative real function ψ λ (z) ≤ 1, satisfying ψ λ (z) = 0 for |z| ≤ λ and ψ λ (z) = 1 for |z| ≥ 2λ. Now we estimate, using , With the optimal choice λ = |ξ| θ/(2+θ) , we obtain completing the proof.
Proof. Because of the boundedness of ρ f and ρ g and of Lemma 3.1, f and the function h := satisfy the assumptions of Lemma 3.2 with χ = χ 1 (after even extension to t < 0).
With an analogous argument for g we obtain ρ f , ρ g ∈ L 2 t H θ/(2+θ) x uniformly in ε.
The conservation law the observation and Lemma 3.1 imply ρ f − ρ g ∈ H 1 t H −1 x which, after interpolation (see Lemma 3.4 below) with the averaging result gives t,x uniformly in ε .
As a consequence, for each 0 ≤ a < b and compact K ⊂ R 3 , a subsequence of ρ f − ρ g converges strongly in L 2 ((a, b) × K) as ε → 0. Since the same is true for √ ρ f ρ g → 1, it also holds for ρ f and ρ g individually as a consequence of the L ∞ bounds. Another application of Lemma 3.1 completes the proof of the convergence statement.
For the derivation of the limiting problem, we pass to the limit in (33) in the distributional sense, denoting the weak limits of R 3 vf ⊥ dv and R 3 vg ⊥ dv, which exist because of (34), by J 1 and J 2 , respectively. Now we multiply the equation for f by v/ε and integrate with respect to v obtaining By the uniform-in-ε boundedness of R 3 v ⊗ vf dv (consequence of Theorem 1.1 and (3)) and by the strong convergence of ρ g , we can pass to the limit, leading to the desired equation for J 1 . For J 2 we proceed analogously.
For completeness we outline the proof of the interpolation lemma used in the proof of Theorem 3.3 with γ = θ 2+θ . Lemma 3.4.