Quantitative Propagation of Chaos in the bimolecular chemical reaction-diffusion model

TL;DR

This paper proves that a stochastic particle system modeling bimolecular chemical reactions converges to a mean field PDE in the large particle limit, providing a quantitative propagation of chaos result using a novel probabilistic approach.

Contribution

The paper introduces a new probabilistic proof for quantitative propagation of chaos in a reaction-diffusion particle system, extending previous entropy-based methods.

Findings

Convergence of particle system to reaction-diffusion PDE as particle number grows

Quantitative estimates on the rate of propagation of chaos

Development of a new martingale-based proof technique

Abstract

We study a stochastic system of $N$ interacting particles which models bimolecular chemical reaction-diffusion. In this model, each particle $i$ carries two attributes: the spatial location $X_t^i\in \mathbb{T}^d$, and the type $\Xi_t^i\in \{1,\cdots,n\}$. While $X_t^i$ is a standard (independent) diffusion process, the evolution of the type $\Xi_t^i$ is described by pairwise interactions between different particles under a series of chemical reactions described by a chemical reaction network. We prove that in the large particle limit the stochastic dynamics converges to a mean field limit which is described by a nonlocal reaction-diffusion partial differential equation. In particular, we obtain a quantitative propagation of chaos result for the interacting particle system. Our proof is based on the relative entropy method used recently by Jabin and Wang \cite{JW18}. The key ingredient of the relative entropy method is a large deviation estimate for a special partition function, which was proved previously by technical combinatorial estimates. We give a simple probabilistic proof based on a novel martingale argument.

Equations

$$S_k+S_l\stackrel{R}{\to }{S}_{k^{\prime }}+S_{l^{\prime }},k,l,k^{\prime },l^{\prime }\in \{1,\cdots ,n_s\}.$$

$$S_1+S_2\to 2S_2.$$

$$\displaystyle\left\{\begin{array}[]{ll}X^{i}_{t}&=X_{0}^{i}+\displaystyle\int_{0}^{t}\sigma(\Xi_{s}^{i})dB^{i}_{s},\\ \Xi_{t}^{i}&=\Xi_{0}^{i}+\displaystyle E^{i}\left(\frac{1}{N}\sum_{j=1,j\neq i}^{N}\int_{0}^{t}\mathbbm{1}\{(\Xi_{s-}^{i},\Xi_{s-}^{j})=(1,2)\}\Phi(X_{s}^{i}-X_{s}^{j})ds\right),\\ &\qquad\mbox{ for }1\leq i\leq N.\end{array}\right.$$

$${\mu }_N(t)=\frac{1}{N}\sum _{i=1}^N{\delta }_{(X_t^i,{\Xi }_t^i)}$$

$$\displaystyle\left\{\begin{array}[]{ll}\partial_{t}u&=\displaystyle\frac{\sigma(1)^{2}}{2}\Delta u-(\Phi*w)u,\\ \partial_{t}w&=\displaystyle\frac{\sigma(2)^{2}}{2}\Delta w+(\Phi*w)u,\end{array}\right.$$

$$(\Phi \ast u)(x)={\int }_{{\mathbb{T}}^d}\Phi (x-y)u(y)dy.$$

$$\displaystyle\left\{\begin{array}[]{ll}\dot{u}&=-\lambda uw,\\ \dot{w}&=\lambda uw,\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ll}U_{t}&=U_{0}-\displaystyle E\left(\lambda\int_{0}^{t}U_{s-}W_{s-}ds\right),\\ W_{t}&=W_{0}+\displaystyle E\left(\lambda\int_{0}^{t}U_{s-}W_{s-}ds\right),\end{array}\right.$$

$$d{X}_t^i=\sigma (t,X_t^i)d{B}_t^i+\frac{1}{N}\sum _{j=1,j\ne i}^{N}F(X_t^i-X_t^j)dt,1\le i\le N.$$

$$\underset{N\ge 2}{\sup }\mathbb{E}\left[\exp \left(\frac{1}{N}\sum _{i,j=1,i\ne j}^{N}f(Y_i,Y_j)\right)\right]<\infty ,$$

$$\sum _{k=1}^n{\alpha }_k{S}_k\stackrel{R}{\to }\sum _{k=1}^n{\beta }_k{S}_k.$$

$$\mathfrak{C}\subset \{{\mathfrak{e}}_k+{\mathfrak{e}}_l:k,l\in \mathfrak{S}\}.$$

$$k+l\stackrel{R}{\to }{k}^{\prime }+l^{\prime },(k,l)\stackrel{R}{\to }(k^{\prime },l^{\prime }),R:(k,l)\to (k^{\prime },l^{\prime }).$$

$$k\le l,k^{\prime }\le l^{\prime }.$$

$$R^{-}=\{k,l\},R^{+}=\{k^{\prime },l^{\prime }\}.$$

$$\Pi =\mathbb{X}\times \mathfrak{S},$$

$$dm=dx\otimes d\#$$

$${\mathbf{y}}_N$$

$$y_k$$

$${\mathbf{x}}_N$$

$${\mathbf{\xi }}_N$$

$$m^N:=m^{\otimes N}=\stackrel{N}{\overbrace{m\otimes \cdots \otimes m}}.$$

$$⟨\phi ,\rho ⟩$$

$$\sigma :\mathfrak{S}\to {\mathbb{R}}^{+},$$

$${\Phi }_R\in L^1(\mathbb{X})\cap L^{\infty }(\mathbb{X}),\text{mbox}with{\Phi }_R(x)={\Phi }_R(-x)\text{mbox}forallx\in \mathbb{X},$$

$$\displaystyle\chi_{k}(\xi)=\left\{\begin{array}[]{ll}1&\mbox{if }\xi=k,\\ 0&\mbox{else}.\end{array}\right.$$

$${\chi }_R^{-}(\xi ,{\xi }^{\prime })={\chi }_k(\xi ){\chi }_l({\xi }^{\prime }),{\chi }_R^{+}(\xi ,{\xi }^{\prime })={\chi }_{k^{\prime }}(\xi ){\chi }_{l^{\prime }}({\xi }^{\prime }).$$

$${\chi }_R^{-}(k,l)=1,{\chi }_R^{-}(l,k)=0,{\chi }_R^{+}(k^{\prime },l^{\prime })=1,{\chi }_R^{+}(l^{\prime },k^{\prime })=0.$$

$$X_t^i=X_0^i+{\int }_0^t\sigma ({\Xi }_s^i)d{B}_s^i,1\le i\le N.$$

$${\tilde{E}}_R^{ij}(t)$$

$$\frac{1}{N}\cdot {\chi }_R^{-}({\Xi }_{t-}^i,{\Xi }_{t-}^j)\cdot {\Phi }_R(X_t^i-X_t^j).$$

$${\Xi }_t^i$$

Full text

1. Introduction

In this paper, we consider a class of stochastic interacting particle systems modeling a chemical reaction-diffusion process. In this model, there are $N$ particles, indexed by $i\in\{1,\cdots,N\}$. Each particle carries two attributes: a location $X_{t}^{i}\in\mathbb{T}^{d}$ (the $d$-dimensional torus), and a chemical type $\Xi_{t}^{i}\in\{S_{1},S_{2},\cdots,S_{n}\}=\{1,\cdots,n_{s}\}$, where $n_{s}$ is the number of distinct chemical species. As time $t$ progresses, $X_{t}^{i}$ diffuses as a standard Brownian motion in space independently, with the speed of diffusion depending on its type $\Xi_{t}^{i}$. The type $\Xi_{t}^{i}$ changes in time according the pairwise chemical interactions between different particles, with transition rates depending on their locations, types, and a set of bimolecular chemical reactions $R_{1},R_{2},\cdots,R_{n_{r}}$ of the form

[TABLE]

Specifically, a pair of particles $i,j$ with types matching the input of a reaction $R$ (i.e. $(\Xi_{t}^{i},\Xi_{t}^{j})=(k,l)$ or $(l,k)$ in (1.1)) may react and instantly change to types $(k^{\prime},l^{\prime})$ or $(l^{\prime},k^{\prime})$ at a random time that depends on their spatial location $X_{t}^{i},X_{t}^{j}$. Our main result (Theorem 2.2) shows that in the limit $N\to\infty$, the empirical measure of the particles converges in a suitable sense to the solution of a nonlinear system of reaction-diffusion equations.

The precise description of the stochastic system will be given in the next section, using the notations from chemical reaction network theory. For the time being, let us consider a simple special case, which involves only two types of particles $S_{1},S_{2}$ (and hence $n_{s}=2$), and a single irreversible chemical reaction (hence $n_{r}=1$)

[TABLE]

As mentioned earlier each particle diffuses independently in space with diffusivity depending on its type. A reaction, which turns a type-$S_{1}$ particle into a type-$S_{2}$ particle, happens at a random time, with a rate depending on its location relative to the type-$S_{2}$ particles. The reaction (1.2) is irreversible, in the sense that once a particle turns into type-$S_{2}$, it can never turn back to type-$S_{1}$ again. Hence, we have the following stochastic system:

[TABLE]

where $\sigma:\{1,2\}\to\mathbb{R}^{+},\Phi\in(L^{1}\cap L^{\infty})(\mathbb{T}^{d})$ are prescribed diffusion coefficients and a reaction rate kernels, $\{B^{i}_{t}\}_{1\leq i\leq N}$ are independent standard $d$-dimensional Brownian motions, and $\{E^{i}(t)\}_{1\leq i\leq N}$ are independent unit Poisson jump processes. Here, if the type $\Xi_{t}^{i}$ of a particle $i$ initially is 1 (hence a type-$S_{1}$ particle), it then turns into 2 at a rate of $\frac{1}{N}\sum_{j:\Xi_{t}^{j}=2}\Phi(X_{t}^{i}-X_{t}^{j})$. The scaling $\frac{1}{N}$ here is critical in deriving the mean field limit. When the type $\Xi_{t}^{i}$ turns into 2, it remains a type-2 forever, because the rate is then zero (as $\mathbbm{1}\{(\Xi_{s}^{i},\Xi_{s}^{j})=(1,2)\}=0$).

The objective of this paper is to establish the convergence, as $N\to\infty$, of the empirical measure of the process (which is a time-dependent random measure on $\mathbb{T}^{d}\times\{1,2,\cdots,n\}$)

[TABLE]

to the solution $\bar{\rho}$ of a mean field limit equation, which is a deterministic system of $n$ non-local reaction-diffusion equations, provided that the initial distribution $\mu_{N}(0)$ converges to $\bar{\rho}(0)$ in some appropriate sense. For instance, in the special case (1.2), the limiting system of equations is

[TABLE]

where the components $u,w:[0,\infty)\times\mathbb{T}^{d}\to\mathbb{R}^{+}$ represent the distribution of type-1,2 particle respectively, and “$*$” denotes the convolution operator

[TABLE]

The limit system in the general case, $n\geq 2$, is given below in (2.30).

In the case that the initial distribution is well-mixed (constant density), the mean-field limit coincides with a mass action system [17]. In this case, the normalized concentrations $\{\bar{\rho}_{k}(t)\}_{k=1}^{n}$ for the $n$ chemical species do not depend on $x$, and they satisfy a system of $n$ ordinary differential equations. In the special case (1.2), this is a simple Lotka-Volterra system, with $(u,w)=(\bar{\rho}_{1},\bar{\rho}_{2})$, and some $\lambda>0$:

[TABLE]

(cf. (1.10)). There is also a stochastic counterpart of these systems in which case the total number $\{M_{k}(t)\}_{t\geq 0}$ of each species $k\in\{1,\cdots,n\}$ is a counting process satisfying a certain coupled stochastic system described by appropriate time-changed Poisson processes, whose rates depend on the current configuration $\{M_{k}(t)\}_{1\leq k\leq n}$. Again, for the special case (1.2), the stochastic system is given by (with $(U,W)=(M_{1},M_{2})$)

[TABLE]

where $E$ is a unit rate Poisson process (cf. (1.6)). With an appropriate scaling (in the rate of reactions), it can be shown that the process $\{M^{N}_{k}(t)\}_{1\leq k\leq n}$ described by the stochastic system converges, in a certain sense, to the solution $\{\bar{\rho}_{k}(t)\}_{1\leq k\leq n}$ of a mass action ODE system as $N\to\infty$ [25, 3]. This convergence result (at least for the case of bimolecular reactions), is a special case of our main result, when the initial distribution is well-mixed spatially (i.e., $X^{i}_{t}$, conditioned on $\Xi^{i}_{t}$, is a uniform random variable).

Models for chemical reactions with spatial diffusion have been studied for more than a century, going back to the work of Smoluchowski [34] (see also [8, 1]). Nevertheless, there are relatively few works that make a mathematically rigorous and quantitative connection between stochastically interacting particles (as a microscopic model) and systems of reaction-diffusion equations, as a mean field limit – this is the motivation for our work. De Massi, Ferrari, and Lebowitz [6] derived a scalar reaction-diffusion equation as the limit of a Glauber-type spin system on a lattice $\mathbb{Z}^{d}$. This result was generalized by Durrett and Neuhauser [10] to allow for more than two states/types, leading to a system of reaction-diffusion equations for reactions of the form $S_{k}\to S_{j}$, with rate depending on the density of other types. The reaction-diffusion limit is useful for studying phase transitions in the underlying stochastic model. See [5] for related results in the spatially-discrete setting. In a continuum setting, Oelschläger [30] analyzed a system of diffusing particles in which each particle may give birth $(S_{k}\to 2S_{k})$, die ($S_{k}\to\emptyset$), or change its type $(S_{k}\to S_{j})$ at rates that depend on the density of other types, leading to a system of reaction-diffusion equations in the infinite population limit. See [7, 23, 24, 12] for related works involving scalar reaction-diffusion equations. Whether in the spatially discrete or continuous setting, most of these models involve only one particle changing its type at a time, although the rate of the reaction may depend on the types and locations of other particles. The recent paper [2] studies the mean field limit of a leader-follower dynamics, which models transitions between two labels (followers and leaders). The mean field limit obtained there involves transport and reaction, but without diffusion. Moreover, the reaction rate depends only on the global state of the system, but not on specific locations. Compared to these works, the stochastic systems that we study allow for two particles to change type simultaneously with location-dependent reaction rate, in a reaction of the form $S_{1}+S_{2}\to S_{3}+S_{4}$, for example. However, the total number of particles is conserved in our systems.

The derivation of macroscopic equations from interacting particle systems is a classical topic in mathematical physics dating back to Maxwell and Boltzmann. A prototype microscopic model that has been studied a lot in the literature is the following McKean-Vlasov system of stochastic differential equations:

[TABLE]

For a Lipschitz continuous force $F$ and a diffusion coefficient $\sigma$, this system has a mean field limit which is a nonlocal drift-diffusion equation (usually called the McKean-Vlasov equation); see e.g. [35] and [14] for a proof. An important concept in studying the particle system (1.14), and others like it, is propagation of chaos, due to Kac [20] and McKean [28, 29] (see also the classical monograph by Sznitman [35]), which roughly means that as $N\to\infty$, any finite group of the particles are asymptotically chaotic (independent) if they are initially chaotic (independent); see Definition 2.1 for a precise statement. In recent years, some progress has been made in proving mean field limits of McKean-Vlasov systems with singular interaction forces or kernel; we refer the interested reader to [13, 18, 15, 19, 33]. Among these results, we specifically mention the recent work by Jabin and Wang [19], who investigate the mean field limit for stochastic systems (1.14) with $F\in W^{-1,\infty}$ which in particular includes the Biot-Savart kernel. They obtain a quantitative estimate for the propagation of chaos in terms of the relative entropy between the $k$-marginals of the joint distribution of $N$ particles and the $k$-tensorized distribution of the mean field limit. The relative entropy method was initiated by Yau [41] in the proof of hydrodynamic limits of Ginzburg-Landau Models, where the relative entropy of the current density with respect to some local Gibbs equilibrium serves as a controllable norm whose convergence to zero implies the existence of macroscopic limit; see [21, 36] for more discussion about this method in interacting particle models from statistical physics. The recent work [19] demonstrated the potential power of the relative entropy method in the study of mean field limits of singular stochastic dynamics. Compared to [41], the relative entropy approach used in [19] (and in the present article) is based on the idea of controlling directly the relative entropy between the density of the particle system and that of the mean field limit, instead of the local equilibrium measure of the particle system.

Our strategy of proving propagation of chaos in the chemical reaction model is similar to that of Jabin and Wang [19]. Specifically, we show the propagation of chaos by proving an explicit estimate for the relative entropy between the joint law of particle system and the tensorized law of the mean field limit. The key ingredient of the proof is a large deviation inequality of the form

[TABLE]

where $\{Y_{i}\}_{i\in\mathbb{N}}$ are i.i.d. random variables, and $f$ is a $L^{\infty}$ function satisfying some appropriate cancellation condition (see Lemma 4.3). Note that this inequality requires weaker regularity condition than the usual condition in obtaining large deviation bounds through Varadhan’s lemma (c.f. [37, Theorem 2.5] and [9, Theorem 1.2.1]) since the function $f$ needs not be continuous. The inequality (1.15) (with a weaker assumption on $f$) was proved by Jabin and Wang (see [19, Theorem 4]) using combinatorial estimates. One of the contributions of this work is to provide a short probabilistic proof, assuming $f$ is bounded. Our proof relies on identifying a martingale structure in the exponential of (1.15) which allows us to conclude the estimate from a sharp Marcinkiewicz-Zygmund inequality for a martingale difference sequence.

One may wonder why we only restrict our attention to the case of bimolecular reactions (1.1) (two inputs, two outputs), instead of a more general reaction structure. First, although one may still discuss the mean field limit for general reaction structure, the notion of the propagation of chaos only applies to those systems that conserve the total number of particles in time. Since our work is based on this notion, it excludes the case of those reactions with uneven inputs and outputs. It is also worthwhile to consider chemical reactions taking $m$ inputs, $m$ outputs for $m\geq 3$. However, proving a mean field limit for these models requires a delicate estimate of the quantity (1.15), which now involves $m$-body interactions. This turns out to be highly-nontrivial and will be investigated in a future work.

The rest of the paper is organized as follows. In the coming section, we begin with a brief introduction to general notions of chaos, chemical reaction networks and some notations, then we set up the $N$-particles system and state the main result of its propagation of chaos property (Theorem 2.2). In Section 3, we will determine the infinitesimal generator $\boldsymbol{\mathcal{L}}_{N}$ of the process $(\boldsymbol{X}^{N}_{t},\boldsymbol{\Xi}^{N}_{t})=\{(X_{t}^{i},\Xi_{t}^{i})\}_{1\leq i\leq N}$ defined by the $N$-particles system, and its dual $\boldsymbol{\mathcal{L}}_{N}^{*}$. The proof of the main theorem is presented in Section 4, which is based on the idea of establishing appropriate differential inequality for the normalized relative entropy and then applying Grönwall lemma. The large deviation inequality and a lemma used in the proof of this main theorem will be proved in Section 5. An appendix, briefly discusses the well-posedness and regularity of the mean field limit equation and the Fokker-Planck equation for $N$-particles system, is also included in Section 6.

Acknowledgement

The work of James Nolen was partially funded through grant DMS-1351653 from the National Science Foundation.

2. Settings and the Main Result

2.1. Bimolecular chemical reaction network

To describe the model, we adopt the notations from the theory of chemical reaction networks (see e.g. [3]), so let us now state the definition of general chemical reaction networks. A chemical reaction network is a triplet $(\mathfrak{S},\mathfrak{C},\mathfrak{R})$ given as follows.

1. (1)

$\mathfrak{S}=\{1,2,\cdots,n\}=\{S_{1},\cdots,S_{n}\}$ is a finite set of chemical species. 2. (2)

$\mathfrak{C}=\{C_{1},C_{2},\cdots,C_{n_{c}}\}\subset\mathbb{N}_{0}^{n}$ is a finite collection of chemical complexes. For instance, $C=\sum_{k=1}^{n}\alpha_{k}\mathfrak{e}_{k}\in\mathbb{N}_{0}^{n}$ (where $\{\mathfrak{e}_{k}\}_{k=1,\cdots,n}$ is the unit vector in $\mathbb{R}^{n}$, $\alpha_{k}\in\mathbb{N}_{0}$) represents the complex formed by a total of $\alpha_{k}$ type-$S_{k}$ particles for each $1\leq k\leq n$. 3. (3)

$\mathfrak{R}=\{R_{1},R_{2},\cdots,R_{n_{r}}\}\subset\mathfrak{C}\times\mathfrak{C}$ is a finite collection of chemical reactions. For instance, $R=(C,C^{\prime})\in\mathfrak{R}$ denotes the chemical reaction between complexes $C\to C^{\prime}$. If $C=\sum_{k=1}^{n}\alpha_{k}\mathfrak{e}_{k},C^{\prime}=\sum_{k=1}^{n}\beta_{k}\mathfrak{e}_{k}$, then the reaction can be written in the form of a chemical equation

[TABLE]

The complex $C$ (the left hand side of the above) will be called the input of the reaction $R$, whereas $C^{\prime}$ (the right hand side) is its output.

In this work, we restrict our attention to the case of a bimolecular reaction network, which is a particular case of a chemical reaction network, whose complexes $C\in\mathfrak{C}$ are of bimolecular form, namely,

[TABLE]

This means that all reactions $R\in\mathfrak{R}$ have the form $S_{k}+S_{l}\xrightarrow{R}S_{k^{\prime}}+S_{l^{\prime}},$ for some $k,l,k^{\prime},l^{\prime}\in\mathfrak{S}$. In this case, we use the notation

[TABLE]

It is possible that $k=l$ and/or $k^{\prime}=l^{\prime}$, so that the two input molecules and/or the two output molecules are of the same species (e.g. (1.2)). Throughout this work, for every chemical reaction equation $R\in\mathfrak{R}$ written as above, to avoid ambiguity we always assume the ascending order for both input and output, that is,

[TABLE]

For $R\in\mathfrak{R}$, we denote $R^{-}$ (resp. $R^{+}$) the input (resp. output) of the reaction. For instance, a chemical reaction of the form (2.1), we have

[TABLE]

We remark that the number of particles, and hence the total mass, are conserved in every bimolecular reaction.

2.2. Notations

We now introduce some general notations being used in the work. Throughout this paper, we will consider stochastic processes in the state space

[TABLE]

where $\mathbb{X}=\mathbb{T}^{d}$ (the $d$-dimensional torus), and the discrete set $\mathfrak{S}=\{1,2,\cdots,n\}$ represents the type of particles. Though the present work is mainly based on the case of a torus, we sometimes will also consider the whole space setting, i.e., $\mathbb{X}=\mathbb{R}^{d}$. We set the variable $y=(x,\xi)\in\Pi$, with $x\in\mathbb{X}$, $\xi\in\mathfrak{S}$, and denote

[TABLE]

the canonical measure on $\Pi$, that is, the product of Lebesgue measure on $\mathbb{T}^{d}$ (or $\mathbb{R}^{d}$), and the counting measure $\#$ on $\{1,2,\cdots,n\}$. For $N\geq 1$, denote variables

[TABLE]

Specifically, the boldface symbols ${\boldsymbol{y}}_{N},{\boldsymbol{x}}_{N},{\boldsymbol{\xi}}_{N}$ are for elements in the $N$-fold product spaces $\Pi^{N},\mathbb{X}^{N},\mathfrak{S}^{N}$, while the normal typeset symbols $y,x,\xi$ denote elements in $\Pi,\mathbb{X},\mathfrak{S}$. Likewise, the boldface uppercase symbols (e.g., $\boldsymbol{X}_{t}^{N},\boldsymbol{\Xi}_{t}^{N}$) denote processes on the $N$-fold product spaces, while the normal uppercase symbols (e.g., $X_{t},\Xi_{t}$) denote processes on $\mathbb{X},\mathfrak{S}$. Also, let $m^{N}$ be the $N$-fold product of the measure $m$, which is a measure on $\Pi^{N}$:

[TABLE]

Finally, for $N\geq 1,p\in[1,\infty]$, and $k\geq 0$, denote the spaces of measures and functions on $\Pi^{N}$ as below:

• •

$\mathcal{M}(\Pi^{N})$ - the family of Borel measures on $\Pi^{N}$;

• •

$\mathcal{P}(\Pi^{N})$ - the family of probability measures on $\Pi^{N}$;

• •

$C_{0}(\Pi^{N})$ - the space of continuous functions on $\Pi^{N}$, vanishing at infinity;

• •

$L^{p}(\Pi^{N})$ - the space of all $L^{p}$- integrable functions on $\Pi^{N}$ (w.r.t. $dm^{N}$);

• •

$C_{0}^{k}(\Pi^{N})$ - the space of functions $f$ so that $f(\cdot,{\boldsymbol{\xi}}_{N})\in C_{0}^{k}(\mathbb{X}^{N})$ for all ${\boldsymbol{\xi}}_{N}\in\mathfrak{S}^{N}$;

• •

$W^{k,p}(\Pi^{N})$ - the space of functions $f$ so that $f(\cdot,{\boldsymbol{\xi}}_{N})\in W^{k,p}(\mathbb{X}^{N})$ for all ${\boldsymbol{\xi}}_{N}\in\mathfrak{S}^{N}$.

Also, define a bilinear form $\langle\cdot,\cdot\rangle:C_{0}(\Pi^{N})\times\mathcal{M}(\Pi^{N})\to\mathbb{R}$ by

[TABLE]

2.3. The stochastic chemical reaction-diffusion system

Fix dimension $d\geq 1$ and a bimolecular chemical reaction network $(\mathfrak{S},\mathfrak{C},\mathfrak{R})$, with $\mathfrak{S}=\{1,2,\cdots,n\}$. Fix a function

[TABLE]

to represent the diffusion coefficient for each type of particle. For each reaction $R\in\mathfrak{R}$ of the form (2.1), we associate a non-negative function

[TABLE]

which will be called the reaction kernel of the chemical reaction $R$. The reaction kernel $\Phi_{R}$ will be used to define the rate at which the reaction $R$ between two particles occurs. One typical choice of reaction kernel, as introduced in [8], could be a cut-off function supported on finite ball, i.e $\Phi_{R}(x)=\chi_{\{|x|\leq r\}}$ for some $r>0$. The symmetry assumption (2.6) is not essential to us and it is only used for the purpose of simplifying expressions in the mean field limit; see Section 3.3 for more details.

Next, we introduce some useful indicator functions. Given a fixed type $k\in\mathfrak{S}$, we denote by $\chi_{k}:\mathfrak{S}\to\{0,1\}$ the indicator

[TABLE]

For $R\in\mathfrak{R}$ with $(k,l)\xrightarrow{R}(k^{\prime},l^{\prime})$, let $\chi_{R}^{\pm}:\mathfrak{S}^{2}\to\{0,1\}$ be the indicator function

[TABLE]

That is, $\chi^{-}_{R}(\xi,\xi^{\prime})$ (resp. $\chi^{+}_{R}$) indicates the event when $(\xi,\xi^{\prime})$ matches with the input $R^{-}$ (resp. output $R^{+}$) of the reaction $R$. We stress here these indicators are order sensitive when $k\neq\ell$ and $k^{\prime}\neq\ell^{\prime}$, in the sense that

[TABLE]

We next introduce the “random ingredients” for the model. Let $i\in\mathbb{N}$ label particles involved in the modeled chemical reaction-diffusion process. For each particle $i$, we associate the following independent random variables and stochastic processes to it:

1. (1)

a $\Pi$-valued random variable $(X_{0}^{i},\Xi^{i}_{0})$; 2. (2)

a standard Brownian motion $\{B^{i}_{t}\}_{t\geq 0}$.

For each $R\in\mathfrak{R}$ and *ordered *pair $(i,j)\in\mathbb{N}^{2}$ with $i\neq j$, we associate also

1. (3)

an independent unit rate Poisson process $\{E^{ij}_{R}(t)\}_{t\geq 0}$;

These Poisson processes will be used to define the counts, up to a certain time, of the type-$R$ reactions that happens between the (ordered) pair of distinct particles $i$ and $j$. The collection of processes $\{B_{t}^{i}\}_{t\geq 0},\{E^{ij}_{R}(t)\}_{t\geq 0}$, for $i,j\in\mathbb{N}$ and $R\in\mathfrak{R}$ are assumed to be independent. Let $(\Omega,\mathcal{F},\mathbb{P})$ be the probability space on which these random variables and processes are defined.

Fix $N\gg 1$ large and consider the stochastic system of $N$-particles described as follows. For each $i\in\{1,\cdots,N\}$, we consider the process $\{(X_{t}^{i},\Xi^{i}_{t})\}_{t\geq 0}$ on the state space $\Pi=\mathbb{X}\times\mathfrak{S}$, with $X^{i}_{t}\in\mathbb{X}$ representing the location of the $i$-th particle at time $t\geq 0$, and $\Xi^{i}_{t}\in\mathfrak{S}=\{1,2,\cdots,n\}$ representing its type. Each particle $i$ diffuses in the spatial domain $\mathbb{X}$ independently, with diffusion coefficient $\sigma(\Xi^{i}_{t})$, depending on its current type $\Xi^{i}_{t}$. Therefore, $X_{t}^{i}$ satisfies the SDE

[TABLE]

Next consider the type process $\{\Xi_{t}^{i}\}$ for a given particle $i$, which is a pure jump process. The type has initial value $\Xi_{0}^{i}$. To describe the evolution, we introduce the reaction counter process $\tilde{E}^{ij}_{R}(t)$ for a given reaction $R\in\mathfrak{R}$ and the (ordered) pair of particles $(i,j)$ (with $i\neq j$), which is the following time-change of the Poisson process $E^{ij}_{R}(t)$:

[TABLE]

Specifically, this is a counting process that counts the number of type-$R$ reactions occurring between the pair $(i,j)$ up to time $t$. The rate of these reactions depends on the relative locations $X_{t}^{i},X_{t}^{j}$, through the reaction kernel $\Phi_{R}$, and is given by

[TABLE]

The scale $\frac{1}{N}$ is the virtue of the mean field interaction, so that the pairwise interaction diminishes as the number of particles $N\to\infty$. Observe that the rate is nonzero only when $\chi_{R}^{-}(\Xi_{t-}^{i},\Xi_{t-}^{j})=1$ and when $\Phi_{R}(X_{t}^{i}-X_{t}^{j})\neq 0$. In other words, a reaction of type-$R$ can happen between the pair $(i,j)$ only when the type $(\Xi_{t-}^{i},\Xi_{t-}^{j})$ matches with the input $R^{-}$ of the reaction (i.e., $\Xi_{t-}^{i}=k,\Xi_{t-}^{j}=l$, see (2.7)), and when their locations $X_{t}^{i},X_{t}^{j}$ are sufficiently close (if $\Phi_{R}$ has a compact support).

As time progresses, a jump $\Xi_{t-}^{i}\to\Xi_{t}^{i}$ occurs when a reaction $(k,l)\xrightarrow{R}(k^{\prime},l^{\prime})$ between particle $i$ and some other particle $j\neq i$ takes place. Namely, a jump in $\Xi_{t}^{i}$ happens at the time $t$ when either a reaction of type-$R$ between the pair $(i,j)$, or $(j,i)$ occurs, and hence, when either reaction counter $\tilde{E}_{R}^{ij}(t)$ or $\tilde{E}^{ji}_{R}(t)$ jumps. When that happens, the type $(\Xi_{t-}^{i},\Xi_{t-}^{j})$ instantly turns into the output $R^{+}$ of the reaction $R$: either $(k^{\prime},l^{\prime})$ or $(l^{\prime},k^{\prime})$. To avoid ambiguity, we enforce the rule of assignment that if the reaction of type-$R$ occurs between the pair $(i,j)$ of particles, the type $(\Xi_{t-}^{i},\Xi_{t-}^{j})$ turns into $(k^{\prime},l^{\prime})$ in order (recall the convention (2.2) that $k^{\prime}\leq l^{\prime}$). This also means, if a reaction occurs between the pair $(j,i)$ instead, and hence $(\Xi_{t-}^{j},\Xi_{t-}^{i})=(k,l)$, it also turns into the same output: $(\Xi_{t}^{j},\Xi_{t}^{i})=(k^{\prime},l^{\prime})$. Specifically, if $(\Xi_{t-}^{i},\Xi_{t-}^{j})$ has value $(k,l)$ (resp. $(l,k)$), then the reaction $\tilde{E}^{ji}_{R}$ (resp. $\tilde{E}^{ij}_{R}$) will not fire, as the rate $\chi_{R}^{-}(\Xi_{t-}^{j},\Xi_{t-}^{i})=0$ (resp., $\chi_{R}^{-}(\Xi_{t-}^{i},\Xi_{t-}^{j})=0$ see (2.8)). Therefore, with the rule of assignment, when a reaction of type-$R$ happens, it always turns types $k\to k^{\prime}$, $l\to l^{\prime}$ in ascending order. In total, the evolution of the type process is described by the following stochastic integrals against the reaction counter processes $\tilde{E}^{ij}_{R}$

[TABLE]

In summary, the dynamics of the process $(\boldsymbol{X}^{N}_{t},\boldsymbol{\Xi}^{N}_{t})=\{(X^{i}_{t},\Xi_{t}^{i})\}_{1\leq i\leq N}$ is described by the following SDE system:

[TABLE]

Observe that $\Xi_{t}^{i}$ takes values in $\{1,\cdots,n\}$. The process $t\mapsto\int_{0}^{t}(k^{\prime}-\Xi_{s-}^{i})d\tilde{E}^{ij}_{R}(s)$ is piecewise constant, making a jump of size $(k^{\prime}-k)$ only at times when $\Xi^{i}_{t}=k$. With probability one, the jumps in $\tilde{E}^{ij}_{R}(t)$ occur at distinct times.

We can also write this equation in vector form. Let $\boldsymbol{B}_{t}^{N}=(B_{t}^{1},\cdots,B_{t}^{N})$, which is the standard Brownian motion on $\mathbb{X}^{N}$, and $\Sigma_{N}(\boldsymbol{\Xi}_{t}^{N})$ be the $(N\times N)$-diagonal matrix with diagonal entries $\sigma(\Xi_{t}^{i})$, for $1\leq i\leq N$. Introduce also the function $\Theta^{ij}_{R}:\mathfrak{S}^{N}\to\mathfrak{S}^{N}$ by

[TABLE]

That is, $\Theta_{R}^{ij}({\boldsymbol{\xi}}_{N})$ is obtained by changing only the $i$- and $j$-coordinate of ${\boldsymbol{\xi}}_{N}$ from $\xi_{i},\xi_{j}$ to the output $k^{\prime},l^{\prime}$ of the reaction $R$ respectively, while leaving the other coordinates unchanged. For later use, we also introduce $\tilde{\Theta}_{R}^{ij}:\mathfrak{S}^{N}\to\mathfrak{S}^{N}$ the “reverse” of the map $\Theta_{R}^{ij}$, which turns the $i,j$-coordinates into the input $k,l$ of $R$:

[TABLE]

Then the vector form of (2.13) is written as follows:

[TABLE]

*Remarks. * Before ending this subsection, let us make one comment about the rule of assignment. By our construction, whenever a reaction $R:(k,l)\to(k^{\prime},l^{\prime})$ takes place, it always turns $k\to k^{\prime},l\to l^{\prime}$, following the ascending order of the input and output. Perhaps a more realistic situation will be turning the input into a prescribed choice of output. For instance, when considering a reaction $R:(1,2)\to(2,3)$, it is more natural to have $1\to 3$, while 2 remains $2\to 2$, instead of following the ascending order $1\to 2,2\to 3$. To be precise, for those reactions $R:(k,l)\to(k^{\prime},l^{\prime})$ with distinct inputs $k\neq l$, we prescribe a surjective assignment $\theta_{R}:\{k,l\}\to\{k^{\prime},l^{\prime}\}$. Whenever a reaction of type-$R$ occurs, it turns $k\to\theta_{R}(k)$, and $l\to\theta_{R}(l)$. If $k=l$, we then switch the types into $(k^{\prime},l^{\prime})$ or $(l^{\prime},k^{\prime})$, each with probability $\frac{1}{2}$. This is equivalent to first fix an order of the right hand side (the output) of (2.1) for each $R$, namely,

[TABLE]

(the order of the output does not matter when $k=l$). Then a reaction turns the input $(k,l)$ into the output, following the order of the right hand side above. One may construct the dynamics with this rule, and of course the result of propagation of chaos still holds, except with a minor difference in the mean field limit equation (2.34) (more precisely, the definition of $\bar{T}_{R}^{+}$ in (2.36)).

2.4. Notions of entropy and chaos.

Before stating our main results, we define some notions of entropy and chaos. For a generic Polish space $\Pi$, let $\mathcal{P}(\Pi)$ be the family of probability measures on $\Pi$. For $N\geq 2$, we denote by $\mathcal{P}_{sym}(\Pi^{N})$ the set of symmetric probability measures on the product space $\Pi^{N}$, that is the set of laws of exchangeable $\Pi^{N}$-valued random variables. Given a symmetric probability density $\rho_{N}\in\mathcal{P}_{sym}(\Pi^{N})$, let us denote by $\rho_{k,N}\in\mathcal{P}(\Pi^{k})$, $1\leq k\leq N$, the $k$-marginal of $\rho_{N}$. For two probability measures $\mu,\nu\in\mathcal{P}(\Pi)$, the relative entropy of $\mu$ with respect to $\nu$ is defined as

[TABLE]

If $\mu$ and $\nu$ have densities $f,g$ with respect to a measure $m$ on $\Pi$, then

[TABLE]

For two symmetric probability measures $\rho_{N},\bar{\rho}_{N}\in\mathcal{P}_{sym}(\Pi^{N})$ with $k$-marginals $\rho_{k,N}$ and $\bar{\rho}_{k,N}$, we also define the following normalized relative entropies

[TABLE]

Definition 2.1** (Kac’s chaos).**

Let $\rho_{N}\in\mathcal{P}_{sym}(\Pi^{N})$ be a sequence of symmetric probability measures on $\Pi^{N},N\geq 1$ and let $\bar{\rho}\in\mathcal{P}(\Pi)$. We say $\rho_{N}$ is $\bar{\rho}$-Kac’s chaotic if one of the following three equivalent conditions holds:

(i) the sequence of two marginals $\rho_{2,N}\rightharpoonup\bar{\rho}\otimes\bar{\rho}$ (converges weakly) as $N\rightarrow\infty$;

(ii) for all $k\geq 1$ fixed,

[TABLE]

(iii) the law $\hat{\rho}_{N}$ of the empirical measure $\mu^{N}_{\boldsymbol{X}^{N}}$ associated to $\boldsymbol{X}^{N}=\{X^{i}\}_{i=1}^{N}\sim\rho_{N}$ satisfies

[TABLE]

The first notion of chaos given in Definition 2.1 (ii) was defined by Kac [20, Section 3]. Sznitman [35] proved the equivalence of the three formulations above. We also refer to [16, Theorem 1.2] for more quantitative statements about Kac’s chaos.

The primary goal of this paper is to prove that the law of the chemical reaction system (2.24) is Kac’s chaotic for all $t\geq 0$. We will achieve this goal by proving a quantitative bound on the normalized relative entropy between the joint distribution of the $N$-particles and the tensorized law of the mean field limit; see Theorem 2.2.

2.5. Main result

The process $\{(\boldsymbol{X}_{t}^{N},\boldsymbol{\Xi}_{t}^{N})\}_{t\geq 0}$ defines a Markov process on the state space $\Pi^{N}$, whose generator $\boldsymbol{\mathcal{L}}_{N}$ and adjoint $\boldsymbol{\mathcal{L}}_{N}^{*}$ will be determined in the next section. Thus, the joint probability distribution $\rho_{N}(t)$ (for $t\geq 0$) of the process $(\boldsymbol{X}_{t}^{N},\boldsymbol{\Xi}^{N}_{t})$, which is a probability measure on $\Pi^{N}$, satisfies the Fokker-Planck equation

[TABLE]

This equation is in fact a linear parabolic system. From the theory of parabolic equations, $\rho_{N}$ admits a smooth (in ${\boldsymbol{x}}_{N}$) density w.r.t. the measure $m^{N}$, provided that the initial data $\rho(0)$ admits a nonnegative $L^{1}$-density; see Section 6, specifically, Proposition 6.1. Abusing notation, we will use the same symbol $\rho_{N}(t,\cdot)$ to denote its density.

Our main objective is to show, as $N\to\infty$, the joint distribution $\rho_{N}$ converges, in some sense, to the tensorization $\bar{\rho}_{N}=\bar{\rho}^{\otimes N}$ of the solution of the mean field limit equation $\bar{\rho}$. The mean field limit equation of this stochastic system, which will be formally derived in Section 3, is the following nonlocal reaction-diffusion system

[TABLE]

where $\bar{\rho}=(u_{1},\cdots,u_{n})$, namely, $\bar{\rho}(t,x,\xi)=u_{\xi}(t,x)$ for $1\leq\xi\leq n$, and $\delta_{\xi l}$ denotes the Dirac delta, and the operator $\Delta$ is the Laplacian in the spatial coordinate $x$. The initial data $\bar{\rho}_{0}=(u_{0,1},u_{0,2},\cdots,u_{0,n})$ is a probability density (with respect to $m$) on $\Pi$, namely,

[TABLE]

Alternatively, $\bar{\rho}$ satisfies the equation

[TABLE]

where $\bar{T}$ is the (nonlinear) operator

[TABLE]

and $\bar{T}_{R}^{\pm}$, for a reaction $R\in\mathfrak{R}$ with $(k,l)\xrightarrow{R}(k^{\prime},l^{\prime})$, is given by

[TABLE]

This nonlinear parabolic system (2.30) is a regularized version of the local reaction-diffusion system ((2.30) with the Dirac delta measure $\lambda_{R}\delta_{0}$, $\lambda_{R}>0$, in place of $\Phi_{R}$). The system is in fact globally well-posed for every nonnegative initial data $\bar{\rho}_{0}\in L^{1}(\Pi)$ from (2.31). Moreover, these solutions are nonnegative and regular, with the total mass $\int_{\Pi}\bar{\rho}(t,y)dm(y)$ conserved in time. Specifically, solutions are regular in the sense that

[TABLE]

Thus, the solution $\{\bar{\rho}(t,\cdot)\}_{t\geq 0}$ is a time-dependent probability density on $\Pi$. The proof of these results (well-posedness, regularity) will be presented in the appendix (Section 6).

Let $\bar{\rho}_{N}=\bar{\rho}^{\otimes N}$ be the tensorized law of $\bar{\rho}$, namely,

[TABLE]

Again, our goal is to show that the distribution $\rho_{N}$ of $(\boldsymbol{X}_{t}^{N},\boldsymbol{\Xi}^{N}_{t})$ converges to $\bar{\rho}_{N}$ in some appropriate sense. To this end, we consider the renormalized relative entropy between $\rho_{N},\bar{\rho}_{N}$:

[TABLE]

(cf. (2.25)). The main result of this work, stated as follows, shows that this quantity vanishes uniformly on $[0,T]$ for any $T>0$ as $N\to\infty$, provided $W_{N}(0)\to 0$.

Theorem 2.2**.**

Let $d\geq 1$, $(\mathfrak{S},\mathfrak{C},\mathfrak{R})$ be a bimolecular chemical reaction network, and $\mathbb{X}=\mathbb{T}^{d}$. Let $\sigma$, $\{\Phi_{R}\}_{R\in\mathfrak{R}}$ be the diffusion coefficients and reaction kernels from (2.5), (2.6). For $N\geq 2$, let $\rho_{N}$ be the law of the process $\{(\boldsymbol{X}_{t}^{N},\boldsymbol{\Xi}_{t}^{N})\}$ described by (2.13) with $\rho_{N}(0)\in\mathcal{P}_{sym}(\Pi^{N})$. Let $\bar{\rho}_{N}=\bar{\rho}^{\otimes N}$, where $\bar{\rho}$ is the unique global solution to the mean field limit equation (2.34) with initial data $\bar{\rho}_{0}$ from (2.31) (guaranteed by Proposition 6.1).

If $\bar{\rho}_{0}\in L^{\infty}(\Pi)$, $\inf_{(x,\xi)\in\Pi}\bar{\rho}_{0}(x,\xi)>0$, then for every $T>0$, there exists a constant $\Lambda_{T}>0$ depending on $\bar{\rho}_{0}$ such that the following estimate holds with $\|\Phi\|_{L^{\infty}}=\sum_{R\in\mathfrak{R}}\|\Phi_{R}\|_{L^{\infty}}$:

[TABLE]

As a direct consequence of the main theorem, we have the propagation of chaos property for the $N$-particle system (2.13).

Corollary 2.3** (Propagation of chaos).**

Assume the same settings as Theorem 2.2. If $\mathcal{H}_{N}(\rho_{N}\|\bar{\rho}_{N})(0)\rightarrow 0$ as $N\rightarrow\infty$, then for each $t\in[0,T]$, $\rho_{N}(t)$ is $\bar{\rho}(t)$-Kac’s chaotic (c.f. Definition 2.1).

Proof.

First, by Lemma 3.3, the probability distribution $\rho_{N}(t)$ is symmetric for any $t\in[0,T]$. Therefore, thanks to the monotonicity [39] of the normalized relative entropy:

[TABLE]

and the Csiszár-Kullback-Pinsker inequality [38, Remark 22.12]:

[TABLE]

one obtains by (2.38) and assumption that for any $k\geq 1$,

[TABLE]

*Remarks. *

1. With additional assumptions imposed for the mean field limit $\bar{\rho}$, one may in fact extend the main result to the whole space setting $\mathbb{X}=\mathbb{R}^{d}$; for example, the following: it holds for some $T,C_{T}>0$ that:

[TABLE]

Although unnatural, this comparability assumption is critical in establishing the bound (2.38) (as $\Lambda_{T}$ depends on $C_{T}$). Indeed, in the case of torus $\mathbb{X}=\mathbb{T}^{d}$, the condition (2.40) alone, which is weaker than the assumption $\inf_{(x,\xi)\in\Pi}\bar{\rho}_{0}(x,\xi)>0$ imposed by Theorem 2.2 (as a consequence of the maximum principle), is sufficient to imply the main result. Of course, to extend the result to $\mathbb{R}^{d}$, one may need more assumptions for $\bar{\rho}$, for instance, some appropriate regularity and decay at infinity condition for $\bar{\rho}$, so that Lemma 4.1 holds. We are by no means to list down these precise assumptions, and we leave it to the interested reader.

2. Since the relative entropy bound (2.38) from Theorem 2.2 has an explicit dependence in the $L^{\infty}$-norm of reaction kernels, one may use that to obtain the convergence to mean field limits for local reaction-diffusion equations. That is, we consider the stochastic system (2.13) with the reaction kernels $\Phi_{R}$ scaled according the the number of particles $N$. More precisely, given $r>0$, let $\Phi_{R}^{r}\in L^{1}(\mathbb{X})$ be the $L^{1}$-rescaled of the reaction kernel $\Phi_{R}$, namely,

[TABLE]

Let $\{(\boldsymbol{X}^{N,r}_{t},\boldsymbol{\Xi}^{N,r}_{t})\}_{t\geq 0}$ satisfy the system (2.24) with $\Phi_{R}^{r}$ in place of $\Phi_{R}$, $\rho_{N}^{r}(t)$ be its distribution, and $\bar{\rho}^{r}$ be the solution of (2.34), with $\Phi_{R}^{r}$ taking the role of $\Phi_{R}$ in (2.36). Then the bound (2.38) implies

[TABLE]

provided that the comparability condition (2.40) holds with $\bar{\rho}^{r}$ in place of $\bar{\rho}$, for some constant $C_{T}$ independent of $r>0$, namely,

[TABLE]

From this bound, if we choose $r=r_{N}$, so that

[TABLE]

(specifically, $r_{N}\gg(\log N)^{-1/d}$), then the propagation of chaos holds for the dynamics $\{(\boldsymbol{X}^{N,r_{N}}_{t},\boldsymbol{\Xi}^{N,r_{N}}_{t})\}_{t\geq 0}$, assuming $\mathcal{H}_{N}(\rho_{N}^{r}\|(\bar{\rho}^{r})^{\otimes N})(0)\leq O(N^{-1})$. In particular, if the initial $N$ particles are sampled independently from the law $\bar{\rho}^{r}(0)$, then $\rho_{N}^{r}(0)=(\bar{\rho}^{r})^{\otimes N}(0)$ and $\mathcal{H}_{N}(\rho_{N}^{r}\|(\bar{\rho}^{r})^{\otimes N})(0)=0$. In this case, as $N\to\infty$ the empirical measure (1.7) converges to a solution $\bar{\rho}$ of the local chemical reaction-diffusion system (2.30) with $\Phi_{R}=\delta_{0}$ (the Dirac delta measure at origin). Of course, all the claims made here are based on the assumptions of (2.41), and the $L^{1}$ convergence of $\bar{\rho}^{r}\to\bar{\rho}$ as $r\searrow 0$. We leave the justification of these assumptions and other details to the interested reader. See also [12] for discussion of this issue in the context of different particle systems.

3. Infinitesimal Generator of the Process $\{(\boldsymbol{X}_{t}^{N},\boldsymbol{\Xi}_{t}^{N})\}$

The main objective of this section is to determine the infinitesimal generator $\boldsymbol{\mathcal{L}}_{N}$ of the process $\{(\boldsymbol{X}_{t}^{N},\boldsymbol{\Xi}_{t}^{N})\}_{t\geq 0}$, and subsequently, its adjoint operator $\boldsymbol{\mathcal{L}}_{N}^{*}$ with respect to the inner product $\langle\cdot,\cdot\rangle$ in $L^{2}(\Pi^{N})$. Throughout this section, $N\geq 1$ will be fixed. For the sake of notational simplicity, we will suppress the subscript $N$ for ${\boldsymbol{y}}={\boldsymbol{y}}_{N}=({\boldsymbol{x}}_{N},{\boldsymbol{\xi}}_{N})=({\boldsymbol{x}},{\boldsymbol{\xi}})\in\Pi^{N}$, $\rho=\rho_{N}$, and the superscript for $\boldsymbol{X}_{t}=\boldsymbol{X}^{N}_{t},\boldsymbol{\Xi}_{t}=\boldsymbol{\Xi}^{N}_{t}$ in the computations involved.

3.1. Generator of the process

The generator of the process $(\boldsymbol{X}_{t},\boldsymbol{\Xi}_{t})$ has two components, corresponding to continuous diffusion (in ${\boldsymbol{x}}$) and to jumps in the ${\boldsymbol{\xi}}$ coordinate (discrete change of type). For $\varphi\in C_{0}^{2}(\Pi^{N})$ ($C^{2}$ in ${\boldsymbol{x}}$), let $\boldsymbol{\Delta}_{N}$ denote the spatial diffusion operator (in ${\boldsymbol{x}}$) on $\Pi^{N}$:

[TABLE]

Next, recall $\Theta^{ij}_{R}$ from (2.17), and define linear operators $\boldsymbol{S}_{N,R}^{ij},\boldsymbol{S}_{N,R},\boldsymbol{S}_{N}:C_{0}(\Pi^{N})\to C_{0}(\Pi^{N})$

[TABLE]

Proposition 3.1**.**

$(\boldsymbol{X}_{t}^{N},\boldsymbol{\Xi}_{t}^{N})$* defined by the $N$-particles system (2.13) is a Markov process on $\Pi^{N}$ with infinitesimal generator $\boldsymbol{\mathcal{L}}_{N}=\boldsymbol{\Delta}_{N}+\boldsymbol{S}_{N}$, where $\boldsymbol{\Delta}_{N},\boldsymbol{S}_{N}$ are from (3.1), (3.2).*

Proof.

Fix a smooth function $\varphi\in C_{0}^{2}(\Pi^{N})$ ($C^{2}$ in ${\boldsymbol{x}}$), and consider $\varphi(\boldsymbol{X}_{t},\boldsymbol{\Xi}_{t})$. Applying Itô’s formula (see e.g. [22, Chapter 9.3]) to the jump-diffusion equation (2.13) (see also (2.24)), we have

[TABLE]

where $M_{t}=\sum_{i=1}^{N}\int_{0}^{t}\sigma(\Xi^{i}_{s})\nabla_{x_{i}}\varphi(\boldsymbol{X}_{s},\boldsymbol{\Xi}_{s})\cdot dB^{i}_{s}$ is a martingale with $M_{0}=0$. Now we take expectation for each term of the identity above. Since $\rho$ is the law of the process $(\boldsymbol{X}_{t},\boldsymbol{\Xi}_{t})$, the expectation of the left hand side is given by

[TABLE]

The above also holds for the first term from the right hand side of (3.3) (with $t=0$). As for the second term $G_{t}$ of (3.3), it again follows by the definition of $\rho$ that

[TABLE]

Now consider the expectation of the last term $H_{t}$ in (3.3), which is given by

[TABLE]

Consider $L_{R}^{ij}(t)$ for a fixed $1\leq i,j\leq N$ with $i\neq j$ and reaction $(k,l)\xrightarrow{R}(k^{\prime},l^{\prime})$ from $\mathfrak{R}$. Observe that the integrand $\varphi(\boldsymbol{X}_{s},\Theta^{ij}_{R}(\boldsymbol{\Xi}_{s-}))-\varphi(\boldsymbol{X}_{s},\boldsymbol{\Xi}_{s-})$ is left-continuous (and hence predictable), and the integrator $\{\tilde{E}^{ij}_{R}(s)\}_{s\geq 0}$, by (2.13), has compensator

[TABLE]

Therefore, the expectation of the process $L_{R}^{ij}(t)$ equals to that of the integrand $\varphi(\boldsymbol{X}_{s},\Theta_{R}^{ij}(\boldsymbol{\Xi}_{s-}))-\varphi(\boldsymbol{X}_{s},\boldsymbol{\Xi}_{s-})$ integrating against the compensator $A_{R}^{ij}(s)$. That is,

[TABLE]

Summing up $\mathbb{E}L^{ij}_{R}(t)$ for $R\in\mathfrak{R}$, and $1\leq i,j\leq N$ with $i\neq j$, we find that $\mathbb{E}H_{t}$ from (3.3) is given by

[TABLE]

Combining all the computations for the expectation of (3.3), we conclude

[TABLE]

which shows that the generator of $\{(\boldsymbol{X}_{t},\boldsymbol{\Xi}_{t})\}_{t\geq 0}$ is $\boldsymbol{\mathcal{L}}_{N}$. ∎

3.2. The adjoint equation

Now we turn our attention to the adjoint operator of $\boldsymbol{\mathcal{L}}_{N}$. The main result of this subsection is stated as follow.

Proposition 3.2**.**

The adjoint of generator $\boldsymbol{\mathcal{L}}_{N}$ in the $L^{2}(\Pi^{N})$-sense is given by $\boldsymbol{\mathcal{L}}_{N}^{*}=\boldsymbol{\Delta}_{N}+\boldsymbol{S}_{N}^{*}$, where

[TABLE]

The law $\rho_{N}(t)$ of the process is the unique strong solution to the Fokker-Planck equation with initial data $\rho_{N}(0)\in L^{1}(\Pi^{N})$:

[TABLE]

Moreover, $\rho_{N}$ has the following regularity: for any $p\in[1,\infty)$

[TABLE]

If additionally $\rho_{N}(0)\in W^{2,\infty}(\Pi^{N})$, then $\rho_{N}\in C([0,\infty);W^{2,p}(\Pi^{N}))\cap C^{1}([0,\infty);L^{p}(\Pi^{N}))$.

Proof.

Observe that $\boldsymbol{\Delta}_{N}$ is self-adjoint (w.r.t. inner product $\langle\cdot,\cdot\rangle$ in $L^{2}(\Pi^{N})$), and thus $\boldsymbol{\Delta}_{N}^{*}=\boldsymbol{\Delta}_{N}$. We now compute the dual of the operator $\boldsymbol{S}_{N}$. Recall the definition from (3.2). To determine the dual of $\boldsymbol{S}_{N}$, we begin with that of $\boldsymbol{S}_{N,R}^{ij}$. Fix a pair of functions $\varphi,\psi\in C_{0}(\Pi^{N})$. Then

[TABLE]

From here, we claim that the following identity holds:

[TABLE]

Indeed, to verify this identity, by permuting the indices it suffices to check it for $(i,j)=(1,2)$. Write ${\boldsymbol{\xi}}=(\xi_{1},\xi_{2},\hat{\boldsymbol{\xi}})$, with $\hat{\boldsymbol{\xi}}=(\xi_{3},\xi_{4},\cdots,\xi_{N})\in\mathfrak{S}^{N-2}$. Then by the definitions of $\chi_{R}^{-},\Theta_{R}^{12},\tilde{\Theta}_{R}^{12}$ (see (2.7), (2.17), (2.21)), and recall also $R\in\mathfrak{R}$ is so that $(k,l)\xrightarrow{R}(k^{\prime},l^{\prime})$, we have

[TABLE]

Using the identity we just established, it follows

[TABLE]

Thus,

[TABLE]

Summing these operators up for $1\leq i,j\leq N$ with $i\neq j$, and $R\in\mathfrak{R}$ yields (3.5).

Existence, uniqueness, and regularity of solutions to (3.6) follows by the standard theory of parabolic systems. In order not to disrupt the flow of presentation, we postpone this part of the proof to the appendix (Proposition 6.1). By the theory of Markov processes, the law $\rho_{N}(t)$ of the process satisfies the forward equation (3.6) [11]. ∎

Despite being straightforward from our construction, for the sake of completeness we give the proof to the preservation of exchangeability for the dynamics (2.13) before ending this subsection.

Lemma 3.3**.**

If $\rho_{N}(0)\in\mathcal{P}_{sym}(\Pi^{N})$, then $\rho_{N}(t)\in\mathcal{P}_{sym}(\Pi^{N})$ for all $t>0$. Equivalently, the particles $\boldsymbol{Y}^{N}_{t}:=(\boldsymbol{X}^{N}_{t},\boldsymbol{\Xi}^{N}_{t})$ are exchangeable for any $t>0$ if the initial particles $\boldsymbol{Y}^{N}_{0}$ are exchangeable.

Proof.

Lemma 3.3 is an immediate consequence of Proposition 3.2, the uniqueness of the PDE (3.6) correspondent to the initial data, and the symmetry of the operator $\boldsymbol{S}_{N}^{*}$ from (3.5). Specifically, for any permutation $\tau$ on $\{1,\dots,N\}$ the operator $\boldsymbol{S}_{N}^{*}$ satisfies

[TABLE]

Here we are using $\tau\psi$ to denote the action of $\tau$ on a function $\psi:\Pi^{N}\to\mathbb{R}$ by $(\tau\psi)({\boldsymbol{y}}_{N})=\psi(\tau{\boldsymbol{y}}_{N})$, where $\tau{\boldsymbol{y}}_{N}$ denote the state variables with permuted indices:

[TABLE]

To see why the symmetry property (3.8) holds, observe that

[TABLE]

where here we suppress the subscript $N$ for ${\boldsymbol{y}},{\boldsymbol{x}},{\boldsymbol{\xi}}$. Recalling the definition (2.21), it is easy to see that for any pair of indices $(i,j)$ with $i\neq j$, and if $(m,\ell)=(\tau^{-1}(i),\tau^{-1}(j))$, then

[TABLE]

Consequently, (3.9) is equivalent to

[TABLE]

3.3. Formal derivation of the mean field limit

We now formally derive the mean field limit equation (2.34) using the result established earlier; a rigorous justification of the mean field limit is carried out in the next section.

If $\mu_{N}(t)\in\mathcal{P}(\Pi)$ is the empirical measure (1.7) and $\phi\in C_{C}^{\infty}(\Pi)$ is any test function, then we claim that

[TABLE]

where for the rest of this section $\langle\cdot,\cdot\rangle$ denotes the bilinear form on $C_{0}(\Pi)\times\mathcal{M}(\Pi)$ (i.e. (2.4) with $N=1$). Recall that the nonlinear operator $\bar{T}$ is defined at (2.35). Therefore, if $\mu_{N}(t)\rightarrow\bar{\rho}(t)dm$ in the appropriate sense as $N\to\infty$, where $\bar{\rho}(t)\in C(\Pi)$ is some smooth deterministic function, then formally passing to the limit in (3.11) we obtain the weak formulation of (2.34):

[TABLE]

The relation (3.11) may be derived as follows. From (3.4) we know that for any test function $\varphi\in C_{0}^{2}(\Pi^{N})$,

[TABLE]

In particular, if we choose $\varphi$ of the form $\varphi({\boldsymbol{x}},{\boldsymbol{\xi}})=N^{-1}\sum_{i=1}^{N}\phi(x_{i},\xi_{i})$ with $\phi\in C_{C}^{\infty}(\Pi)$, then the left side reads

[TABLE]

Moreover,

[TABLE]

Consider $\mathbb{E}[\boldsymbol{S}_{N}\varphi(\boldsymbol{X}_{s},\boldsymbol{\Xi}_{s})]$. Recall the definition of $\boldsymbol{S}_{N}=\sum_{R\in\mathfrak{R}}\boldsymbol{S}_{N,R}$ from (3.2), so consider $\mathbb{E}[(\boldsymbol{S}_{N,R}\varphi)(\boldsymbol{X}_{s},\boldsymbol{\Xi}_{s})]$ for a fixed $R:(k,l)\to(k^{\prime},l^{\prime})$ from $\mathfrak{R}$. By the definition of $\Theta_{R}^{ij}$,

[TABLE]

Consider $I_{1}=I_{1}^{R}$. Let us include the diagonal terms $(i=j)$, which is at most $CN^{-1}$, with constant $C$ depending on $\|\Phi_{R}\|_{L^{\infty}}$, $\phi$, into the double summation $I_{1}$. Noting also by definition that $\chi_{R}^{-}(\xi_{i},\xi_{j})=\chi_{k}(\xi_{i})\cdot\chi_{l}(\xi_{j})$, one then has

[TABLE]

($\#$ denotes the counting measure on $\mathfrak{S}$) where in obtaining the second and third equality above we used the fact that

[TABLE]

In the same manner, one can derive similar expressions for the remaining three terms in (3.14), specifically for $I_{3}$ we have

[TABLE]

Note that to derive the expression above for $I_{3}$ we also used the symmetry of the kernel $\Phi_{R}$ defined in (2.6), without which the convolution above would have been replaced by $\Phi_{R}^{-}\ast\mu_{N}$ where $\Phi_{R}^{-}(\cdot)=\Phi_{R}(-\cdot)$. We make this symmetry assumption to simplify expressions of the mean field limit (2.30). Summing up $I_{1},I_{3}$, and using the notation $\bar{T}_{R}^{+}$ from (2.36), we have

[TABLE]

Here, $\bar{T}_{R}^{\pm}$ from (2.36), initially defined on $C_{0}(\Pi)\to C_{0}(\Pi)$, can be extended naturally to a map from the space of signed measures on $\Pi$ to itself (as $\Phi_{R}$ is bounded integrable). The same computation for $I_{2},I_{4}$ gives

[TABLE]

and hence by the definition $\bar{T}_{R}=\bar{T}_{R}^{-}+\bar{T}_{R}^{+}$ (see (2.35))

[TABLE]

Summing up $R\in\mathfrak{R}$ in (3.14) and recalling $\bar{T}$ from (2.35), we arrive at

[TABLE]

Combining this with (3.12), (3.13), we have shown (3.11). This concludes the formal derivation of the mean field limit equation.

4. Proof of the Main Result (Theorem 2.2)

To prove Theorem 2.2, we will compare the joint distribution of $N$ particles with the tensorized law $\bar{\rho}_{N}:=\bar{\rho}^{\otimes N}$ of the mean field limit in terms of their relative entropy. Recall that the joint distribution $\rho_{N}(t)$ of the process $\{(X^{i}_{t},\Xi_{t}^{i})_{i=1,\cdots,N}\}_{t\geq 0}$, by Proposition 3.2, satisfies the Fokker-Planck equation (3.6). Before we estimate the relative entropy, let us first identify the PDE satisfied by the tensorized law $\bar{\rho}_{N}$. Recall that the mean field limit $\bar{\rho}$ satisfies the system (2.34) with initial data (2.31) (see Proposition 6.1 for the existence and uniqueness result). The tensorized law

[TABLE]

then solves the PDE system

[TABLE]

where $\boldsymbol{\Delta}_{N}$ is the diffusion operator from (3.1), and $\bar{\boldsymbol{T}}_{N}$ is as follows, with $\bar{T}$ from (2.35):

[TABLE]

Recall the normalized relative entropy of $\rho_{N}$ with respect to $\bar{\rho}_{N}$:

[TABLE]

Our main objective is to establish a differential inequality for this quantity, then invoke Grönwall lemma to obtain an estimate in terms of $N$ and $W_{N}(0)$. Taking the time derivative of $W_{N}(t)$ and using equations (3.6), (4.1), we have for all $t>0$ that

[TABLE]

where

[TABLE]

All functions $\rho_{N},\bar{\rho}_{N}$ involved above are evaluated at time $t>0$. The calculations above can be rigorously justified since the solutions have sufficient regularity, namely, for any $p\in[1,\infty)$,

[TABLE]

and $\rho_{N},\bar{\rho}_{N}>0$ are bounded away from zero when $t>0$. (See Propositions 3.2, 6.1 and 6.3.) We proceed with estimating the quantities $D(t)$ and $G(t)$.

4.1. Estimating $D(t)$, Dissipation of Relative Entropy

$D(t)$ is in fact nonpositive, due to the diffusive nature of the operator $\boldsymbol{\mathcal{L}}_{N}^{*}$. Specifically, it is due to the following lemma.

Lemma 4.1**.**

Let $\rho_{N},\bar{\rho}_{N}$ be given in Theorem 2.2, with $\mathcal{H}(\rho_{N}\|\bar{\rho}_{N})(0)<\infty$. Then the integral (4.4) is finite and nonpositive for all $t>0$.

Proof of Lemma 4.1.

The strict positivity condition for $\bar{\rho}_{0}$ in the main theorem is not essential in the proof, and we will prove the lemma with the weaker assumption of $\bar{\rho}_{0}$ being propagated in the sense defined before Proposition 6.3. We will first show the following integral is finite and nonpositive, for any appropriate pair of densities $\rho,\tilde{\rho}\in L^{1}(\Pi^{N})$:

[TABLE]

Let $\{u_{t},\tilde{u}_{t}\}_{t\geq 0}$ be the solutions of (3.6), with initial data $u_{0}=\rho,\tilde{u}_{0}=\tilde{\rho}$, namely,

[TABLE]

The proof of (4.5) then breaks into two parts: first, we show that the relative entropy between two solutions $W(t)=\mathcal{H}(u_{t}\|\tilde{u}_{t})$ is non-increasing w.r.t. the time variable $t\geq 0$, then show that the integral in (4.5) is the derivative of $W(t)$ evaluating at $t=0$, and thus must be non-positive. The precise condition for $\rho,\tilde{\rho}$ will be given later.

We start with verifying the claim that $W(t)$ is non-increasing, for any initial data $\rho,\tilde{\rho}\in L^{1}(\Pi^{N})$ with $W(0)=\mathcal{H}(\rho\|\tilde{\rho})<\infty$. In fact, this monotonicity property of relative entropy, often referred to as the data processing inequality [40], is well-known in physics and information theory community, and we will provide a short proof to it. Let $(\mathcal{X},\mu)$ be a measure space and consider two probability densities $p(x,y),\tilde{p}(x,y)\in L^{1}(\mathcal{X}^{2},\mu^{\otimes 2})$. Denote by $p(y|x),\tilde{p}(y|x)$ their conditional densities, and by $p(x),\tilde{p}(x)$ the $x$-marginals. With a slight abuse of notation, we denote the averaged conditional relative entropy between $p(y|x)$ and $\tilde{p}(y|x)$ by

[TABLE]

By the chain rule [9, Theorem C.3.1] of the relative entropy, we have

[TABLE]

Now return to the solutions $\{u_{t},\tilde{u}_{t}\}_{t\geq 0}$ defined earlier, at (4.6). Recall from Proposition 3.1 that they are probability densities of the Markov process $\{\boldsymbol{Y}_{t}=(\boldsymbol{X}_{t},\boldsymbol{\Xi}_{t})\}_{t\geq 0}$ defined by (2.13), with initial distribution $\rho,\tilde{\rho}$ respectively. Given $t,s\geq 0$, let $u_{t,s},\tilde{u}_{t,s}$ denote the joint probability density of $(\boldsymbol{Y}_{t},\boldsymbol{Y}_{s})$. Similarly, let $u_{t|s},\tilde{u}_{t|s}$ be the conditioned densities of $\boldsymbol{Y}_{t}$ given $\boldsymbol{Y}_{s}$. By the chain rule, for any $t,h\geq 0$,

[TABLE]

Since $u,\tilde{u}$ are defined by the same Markov process (with the same transition kernel), their condition densities coincide, i.e. $u_{t+h|t}=\tilde{u}_{t+h|t}$. As a result, the second term of the left hand side vanishes. By the non-negativity of relative entropy, it follows from the last equation that

[TABLE]

Therefore, $t\mapsto W(t)$ is nonincreasing.

We next proceed with verifying that the integral (4.5) is given by $W^{\prime}(0)$. Indeed, formally differentiating $W(t)$ respect to $t$ and using (4.6) and nonincreasing of $W(t)$ yield

[TABLE]

In conclusion, if $\rho,\tilde{\rho}\in L^{1}(\Pi^{N})$ are so that $t\mapsto W(t)=\mathcal{H}(u_{t}\|\tilde{u}_{t})$ is (right-hand) differentiable at $t=0$, then (4.5) holds.

Now let us apply (4.5) to conclude Lemma 4.1. Fix $t_{0}>0$ and consider the integral $D(t_{0})$ from (4.4). Using (4.5) with $\rho_{N}(t_{0}),\bar{\rho}_{N}(t_{0})$ taking the roles of $\rho,\tilde{\rho}$, we then have $D(t_{0})\leq 0$, provided that the correspondent $W(t)$ is differentiable at $t=0$. So it remains to check this differentiability condition.

First, from the definition, $\tilde{u}_{t}$ is the solution of (3.6) with initial data $\tilde{u}_{0}=\bar{\rho}_{N}(t_{0})=\bar{\rho}^{\otimes N}(t_{0})$, where $\bar{\rho}$ is the solution of (2.34). By Propositions 6.1 and 6.3 (and under the assumptions of Theorem 2.2), we know that $\bar{\rho}(t_{0})\in W^{2,p}(\Pi)$ for all $p\in[1,\infty)$, $\inf_{\Pi}\bar{\rho}(t_{0})>0$, and hence $\bar{\rho}_{N}(t_{0})\in W^{2,p}(\Pi^{N}),\inf_{\Pi^{N}}\bar{\rho}_{N}(t_{0})>0$. Thus, by Propositions 6.1 and 6.3 again, we have for any $p\in[1,\infty)$ and $\delta>0$ that

[TABLE]

Next, since $u_{t}$ is the solution of (3.6) with $u_{0}=\rho_{N}(t_{0})$, we have then $u_{t}=\rho_{N}(t_{0}+t)$. By the regularity of solutions (Proposition 6.1), we have $u_{t}\in C^{1}([0,\infty);L^{p}(\Pi^{N}))$ for any $p\in[1,\infty)$. Moreover, by Proposition 6.3, the function $u_{t}$ has the following property: for any ${\boldsymbol{\xi}}\in\mathfrak{S}^{N}$

[TABLE]

Let $\mathcal{A}\subset\mathfrak{S}^{N}$ be the set of all ${\boldsymbol{\xi}}$ so that the latter of the above holds. Then

[TABLE]

Since $u_{t},\tilde{u}_{t}\in C^{1}([0,\infty);L^{p}(\Pi^{N}))$ are uniformly bounded above, and below from [math] over a time interval $[0,\delta]\times\mathbb{X}^{N}$ on the set $\mathcal{A}$, this follows $\chi_{\mathcal{A}}\log(\frac{u_{t}}{\tilde{u}_{t}})\in C^{1}([0,\delta);L^{p}(\Pi^{N}))$ for any $p\in[1,\infty)$. This implies $W(t)$ is differentiable at $t=0$, which finishes the proof of Lemma 4.1. ∎

4.2. Estimating $G(t)$

Next consider $G(t)$ from (4.3). In the coming computation, we suppress the time variable $t$ for $\rho_{N},\bar{\rho}_{N}$, and the subscript $N$ for ${\boldsymbol{y}}_{N},{\boldsymbol{x}}_{N},{\boldsymbol{\xi}}_{N},\bar{\boldsymbol{T}}_{N},\boldsymbol{S}_{N}^{*}$ and so on, since $N$ will be fixed. To further ease our notation, we also introduce functions $(u_{1},\cdots,u_{n})=\bar{\rho}$, that is,

[TABLE]

Recall that these functions satisfy the mean field limit equation (2.34). Recall also the quantity $G(t)$ is the expectation, against $\rho_{N}$, of the function $\bar{\rho}_{N}^{-1}(\boldsymbol{T}-\boldsymbol{S}^{*})\bar{\rho}_{N}$. We begin with considering $\bar{\rho}_{N}^{-1}\bar{\boldsymbol{T}}\bar{\rho}_{N}$. By the definitions of $\bar{\boldsymbol{T}},\bar{T}$ from (4.2), (2.35), we have

[TABLE]

where, for a given reaction $(k,l)\xrightarrow{R}(k^{\prime},l^{\prime})$ from $\mathfrak{R}$,

[TABLE]

We can symmetrize (4.7) by introducing an extra summation over $j=1,\cdots,N$, so that

[TABLE]

On the other hand, by the definition of $\boldsymbol{S}^{*}_{N}$ from (3.5) and the tensorized law $\bar{\rho}_{N}$,

[TABLE]

where

[TABLE]

Now define

[TABLE]

By (4.9), (4.10), we have

[TABLE]

and so the quantity $G(t)$ is given by

[TABLE]

Introduce the quantity

[TABLE]

Then the diagonal term $G_{2}(t)$ is simply bounded by

[TABLE]

As for $G_{1}(t)$, to set up a differential inequality for Grönwall lemma, we need the following inequality which essentially follows from the variational characterization for relative entropy.

Lemma 4.2** ([19, Lemma 1]).**

Let $N\geq 1$ and $\rho,\bar{\rho}$ be two probability measures on the space $\Pi^{N}$. For every $\eta>0$ and $\Psi\in L^{\infty}(\Pi^{N})$, it holds

[TABLE]

Applying this lemma to $G_{1}(t)$, with $\Psi=N^{-2}\sum_{i,j=1,i\neq j}^{N}f_{t}(y_{i},y_{j})$, $\rho_{N}(t),\bar{\rho}_{N}(t)$ in place of $\rho,\bar{\rho}$, and with some $\eta>0$ (depending on $f_{t}$, c.f. (4.12)) to be determined later, we have

[TABLE]

To this end, we need an estimate of the exponential moment on the right hand side. This can be achieved by using the following large deviation inequality.

Lemma 4.3**.**

Let $(\Pi,\bar{\rho})$ be a probability space, $\{Y_{1},Y_{2},\cdots\}$ be a sequence of i.i.d. $\Pi$-valued random variables with common distribution $\bar{\rho}$, and $f\in L^{\infty}(\Pi^{2},\bar{\rho}^{\otimes 2})$ be a bounded measurable function satisfying the following marginal mean zero conditions:

[TABLE]

Then there exists a constant $\eta=\eta(f)>0$ such that

[TABLE]

Specifically, $\eta(f)$ can be chosen to be $2\sqrt{2}e\|f\|_{L^{\infty}(\Pi^{2},\bar{\rho}^{\otimes 2})}$.

As mentioned earlier, Lemma 4.3 was first proved by Jabin and Wang [19] (see Theorem 4 therein), and we will provide a shorter proof in the coming section, after completing the proof of the main result. To finish the estimation, we apply this lemma to the second integral of (4.14). Prior to that, we first verify the marginal mean zero condition (4.15), required by the lemma, for the function $f_{t}(y,y^{\prime})$.

Lemma 4.4**.**

The function $f_{t}(y,y^{\prime})$ defined in (4.12) satisfies the marginal mean zero conditions (4.15):

[TABLE]

Proof.

Recall from the definition (4.12) that $f_{t}(y,y^{\prime})=\sum_{R\in\mathfrak{R}}F_{t}^{R}(y,y^{\prime})$, where $F_{t}^{R}(y,y^{\prime})=A_{t}^{R}(y^{\prime})+\hat{A}_{t}^{R}(y)-B^{R}_{t}(y,y^{\prime})$, with $A_{t}^{R},\hat{A}_{t}^{R},B_{t}^{R}$ from (4.8), (4.11). The lemma is a direct consequence of the following four identities:

[TABLE]

[TABLE]

Let us begin with the first one. By (4.8),

[TABLE]

This establishes the first identity. Swapping the role of $k$ and $l$, and replacing $l^{\prime}$ by $k^{\prime}$ in the computation above yields the second one. Now consider the third identity. Again by the definition of $B_{t}^{R}$ (see (4.11)), we have

[TABLE]

The last step is due to $\sum_{\xi=1}^{n}\chi_{R}^{+}(\xi,\xi^{\prime})=\chi_{l^{\prime}}(\xi^{\prime})$, $\sum_{\xi=1}^{n}\chi_{R}^{-}(\xi,\xi^{\prime})=\chi_{l}(\xi^{\prime})$ (see (2.7)). Again, swapping the role of $k,l$ gives the fourth identity.

From these calculations, it follows that

[TABLE]

and similarly for $\int_{\Pi}f_{t}(y,y^{\prime})\bar{\rho}(t,y^{\prime})dm(y^{\prime})=0$. ∎

Applying Lemma 4.3 to the integral from (4.14) with $\bar{\rho}(t),f_{t}$ in place of $\bar{\rho},f$, and recall the quantity (4.13), it follows

[TABLE]

4.3. Conclusion

Substituting $D\leq 0$ and the estimates for $G_{1},G_{2}$ into (4.3), we have

[TABLE]

for some universal constant $C>0$ and $K_{t}$ from (4.13). Let us now get a bound for $K_{t}$, assuming $\bar{\rho}$ satisfies the comparability condition (2.40). Note that if $\mathbb{X}=\mathbb{T}^{d}$, then the solution $\bar{\rho}(t)$ of the system (2.34) is locally bounded and strictly positive, because of the assumption of Theorem 2.2 $\bar{\rho}_{0}\in L^{\infty}(\Pi),\inf_{y\in\Pi}\bar{\rho}(0,y)>0$ (see Propositions 6.1, 6.3). Therefore, for every $T>0$, the condition (2.40) holds for some constant $C_{T}>0$ depending on $\bar{\rho}_{0}$. By (4.8), (4.11), and Young’s inequality, it follows

[TABLE]

Hence, we have

[TABLE]

where recall that $\|\Phi\|_{L^{\infty}}=\sum_{R\in\mathfrak{R}}\|\Phi_{R}\|_{L^{\infty}}<\infty$, which is bounded (as $\mathfrak{R}$ is a finite set). Now applying Grönwall lemma to (4.17), we establish the bound

[TABLE]

The desired estimate from (2.38) follows, if we set $\Lambda_{T}=C\tilde{C}_{T}$. This finishes the proof of the main theorem.

5. The Large Deviation Inequality (Proof of Lemma 4.3)

The main objective of this section is to prove the large deviation inequality (4.16) from Lemma 4.3. Before proceeding to its proof, let us first make a few comments about the result. In Lemma 4.3 we do not require that $f$ is continuous. In fact, if $f$ is bounded and continuous on $\Pi^{2}$, then the classical large deviation principle of empirical measures (see e.g. [4]) would imply that

[TABLE]

where the constant $m$ is characterized by

[TABLE]

The last line above follows from the mean zero condition (4.15). From above one obtains that $m=0$ if $\eta$ is chosen large enough, which proves

[TABLE]

Despite being a weaker estimate compared to (4.16), this qualitative result alone is sufficient to conclude the propagation of chaos of the particle system (at least for the case of reaction kernels $\Phi_{R}$ being continuous), except without an explicit bound. The quantitative estimate (4.16) was proved by Jabin and Wang [19] under the assumption that $f$ satisfies

[TABLE]

for some constant $C>0$, which is a weaker assumption than $\|f\|_{L^{\infty}}<\infty$. Their proof relies on some sophisticated combinatorial analysis by taking account of cancellations due to the mean zero condition (4.15). We provide a simple probabilistic alternative. Our proof to Lemma 4.3 relies on a characterization of the exponential moment, given as follows:

Lemma 5.1**.**

Let $Z$ be a random variable satisfying the bound for some $\gamma>0$:

[TABLE]

Then it holds

[TABLE]

Proof.

Let $\eta=2e\gamma$. Expanding $e^{\eta^{-1}Z}$ by Taylor series, and using the assumption,

[TABLE]

The only inequality above follows by the simple inequality $k^{k}(k!)^{-1}\leq e^{k}$, which follows from the Stirling’s approximation. ∎

In the proof of (4.16) we use a Marcinkiewicz-Zygmund type inequality for martingales, which, roughly speaking, bounds the $L^{p}$-norm of a martingale by the root sums squared of the $L^{p}$-norm of its martingale increments. Specifically, we will use the following sharp version of inequality (5.1), due to Rio [32]. When $p=2$, (5.1) holds as an equality, which is due to the Itô isometry for discrete martingales. Also, the constant $p-1$ on the right hand side of (5.1) is known to be sharp. We point out also the sharpness of this estimate, specifically the growth rate as $p\to\infty$, plays a critical role in our argument.

Lemma 5.2** ([32, Theorem 2.1]).**

Let $p\geq 2$, and $\{X_{k}\}_{k\geq 1}$ be a sequence of $L^{p}$-martingale differences with respect to a filtration $\{\mathcal{F}_{k}\}_{k\geq 0}$. That is, for each $k\geq 1$, $X_{k}$ is $\mathcal{F}_{k}$-measurable, $X_{k}\in L^{p}$, and $\mathbb{E}[X_{k}|\mathcal{F}_{k-1}]=0$. Let $S_{n}=\sum_{k=1}^{n}X_{k}$ (which is a martingale). Then for all $n\geq 1$ we have

[TABLE]

Now we are ready to prove Lemma 4.3.

Proof of Lemma 4.3.

Let $\{Y_{j}\}_{j\in\mathbb{N}}$ be a sequence of $\Pi$-valued i.i.d. random variables, and $f\in L^{\infty}(\Pi^{2},\bar{\rho}^{\otimes 2})$ as in the lemma. For $N\geq 2$, denote the random variables

[TABLE]

By Lemma 5.1, the statement (4.16) follows if we show the following uniform-in-$N$ bound is valid:

[TABLE]

We first write $M_{N}$ as the sum of martingale differences. Indeed, we have

[TABLE]

($D_{1}$ is set to be zero). Observe here that $\{D_{k}\}_{k=1,\cdots,N}$ forms a $L^{p}$ martingale difference w.r.t. the filtration $\{\mathcal{F}_{k}=\sigma(Y_{1},\cdots,Y_{k})\}_{k\geq 0}$ ($\mathcal{F}_{0}$ is set to be the trivial $\sigma$-algebra) for any $p\in[1,\infty)$ (and hence $\{M_{k}\}_{1\leq k\leq N}$ is an $L^{p}$-martingale). Indeed, each $D_{k}\in L^{p}$ because $f$ is bounded. Moreover, the marginal mean zero condition (4.15) implies

[TABLE]

Using Lemma 5.2 (with $N,M_{N},D_{k}$ in place of $n,S_{n},X_{k}$), we establish the following bound for all $p\geq 2$:

[TABLE]

Now consider $\|D_{k}\|_{L^{p}}$ for each $1\leq k\leq N$. We again write $D_{k}$ as a sum of martingale difference, namely,

[TABLE]

Note that $\|B^{k}_{j}\|_{L^{p}}\leq 2\|f\|_{L^{\infty}(\Pi^{2},\bar{\rho}^{\otimes 2})}$ for each $k,j$, and $\{B_{j}^{k}\}_{1\leq j\leq k-1}$ again forms a sequence of martingale differences w.r.t. the filtration $\{\tilde{\mathcal{F}}_{j}\}_{0\leq j\leq k-1}$, where $\mathcal{F}_{j}=\sigma(Y_{k},Y_{1},\cdots,Y_{j})$. Specifically, $\mathbb{E}[B_{j}^{k}|\tilde{\mathcal{F}}_{j-1}]=0$ follows directly from the marginal mean zero condition (4.15). So we may again apply Lemma 5.2 (with $k-1,D_{k},B^{k}_{j}$ in place of $n,S_{n},X_{k}$) to obtain

[TABLE]

Inserting this estimate into (5.3), it follows

[TABLE]

Now we return to $A_{N}(f)=N^{-1}M_{N}$. Using the $L^{p}$-bound for $M_{N}$ that we just established, we have for all $k\geq 2$ that

[TABLE]

For $k=1$, we have $|\mathbb{E}A_{N}(f)|=0$ by (4.15). This concludes the proof of (5.2). ∎

6. Appendix: Well-posedness and Regularity of Semilinear Parabolic Systems

In this appendix, we provide a brief discussion on the systems (3.6) and (2.34), particularly, the well-posedness of the correspondent Cauchy problem, and regularity of solutions. These results follow by the classical theory of semigroups and the standard construction of solutions to ODEs using the contraction mapping theorem. Though elementary, for the sake of completeness we also present their proofs. For a detailed discussion, we point the reader to, for instance, [31].

To state a result that is applicable for both the (linear) Fokker-Planck equation (3.6) on $\Pi^{N}$ and the (nonlinear) mean field limit system (2.34) on $\Pi$, we will consider a general form of equation. Let $\mathbb{X}=\mathbb{T}^{D}$ or $\mathbb{R}^{D}$, with spatial dimension $D\in\mathbb{N}$, $G$ be a finite set of indices, and denote $\Gamma=\mathbb{X}\times G$. Denote also the variables $y=(x,\xi)\in\mathbb{X}\times G$, and the measure $dm=dx\otimes d\#$ on $\mathbb{X}\times G$, where $\#$ denotes the counting measure on the index set $G$. For $p\in[1,\infty]$, denote the Banach space and norm

[TABLE]

and $X_{p}^{+}\subset X_{p}$ be the (closed) subset of all nonnegative $(L^{1}\cap L^{p})(\Gamma)$-functions. Similarly, for $k\geq 0$ and $p\in[1,\infty]$, we denote $W^{k,p}(\Gamma)$ the Banach space of all functions $f$ on $\Gamma$ so that $f(\cdot,\xi)\in W^{k,p}(\mathbb{X})$ for every $\xi\in G$, with the norm $\|f\|_{W^{k,p}(\Gamma)}=\sum_{\xi\in G}\|f(\cdot,\xi)\|_{W^{2,p}(\mathbb{X})}$.

We consider the evolution equation on $\mathbb{X}\times G$, given by

[TABLE]

where $A$ is the $D$-dimensional elliptic operator on $\mathbb{X}\times G$ given by

[TABLE]

and $T^{\pm}:L^{1}_{loc}(\Gamma)\to L^{1}_{loc}(\Gamma)$ (possibly nonlinear). If regarding $\xi\in G$ as an index, one may view the forward equation (6.1) as a parabolic system with $n=|G|$ equations in variables $(t,x)\in(0,\infty)\times\mathbb{X}$. Namely, for each $\xi\in G$, the component $u_{\xi}(t,x)=\rho(t,x,\xi)$ satisfies the parabolic equation

[TABLE]

where $A_{\xi}$ is the elliptic operator on $\mathbb{X}$ given by (6.2) with $\xi\in G$ fixed.

Throughout this section we assume the following hypothesis for the maps $T^{\pm}$: it holds for some constant $C>0$ that

1. (T1)

(Boundedness and local Lipschitz continuity). For every $p\in[1,\infty)$ and $\rho,\hat{\rho}\in X_{p}$, it holds

[TABLE] 2. (T2)

(Nonnegativity of $T^{+}$). $T^{+}$ maps nonnegative functions to nonnegative functions; 3. (T3)

(Pointwise bound of $T^{-}$). For every $\rho\in L^{1}(\Gamma)$, we have

[TABLE] 4. (T4)

(Mass conservation). For every $\rho\in L^{1}(\Gamma)$, it holds

[TABLE]

Both forward equations (3.6) and (2.34) are special cases of (6.1). For the Fokker-Planck equation (3.6) on $\Pi^{N}$, we have the spatial domain $\mathbb{X}=\mathbb{T}^{dN}$ or $\mathbb{R}^{dN}$, the index set $G=\mathfrak{S}^{N}$ (and so $D=dN$, $\Gamma=\Pi^{N}$), $A=\boldsymbol{\Delta}_{N}$ from (3.1), and the maps $T^{\pm}$, by (3.5), are given by

[TABLE]

Since in this case $T^{\pm}$ are bounded linear in $X_{p}$ (as $\Phi_{R}\in L^{\infty}$), (T1)–(T4) hold with $C>0$ depending on $N$ and $\|\Phi\|_{L^{\infty}}=\sum_{R\in\mathfrak{R}}\|\Phi_{R}\|_{L^{\infty}}$. ((T4) is a straightforward computation from the definition (3.5) of $\boldsymbol{S}_{N}^{*}$.) For the limit system (2.34), we have $\mathbb{X}=\mathbb{T}^{d}$ or $\mathbb{R}^{d}$, $G=\mathfrak{S}$ (hence $D=d$, $\Gamma=\Pi$), and $A=\frac{1}{2}\sigma(\xi)^{2}\Delta$. The maps $T^{\pm}$ are given by $T^{\pm}=\sum_{R\in\mathfrak{R}}\bar{T}_{R}^{\pm}$ from (2.35). Conditions from (T2)–(T4) can be easily verified. For (T1), the boundedness condition (with $C$ depending on $\|\Phi\|_{L^{\infty}}$) is a direct consequence of Young’s inequality and the definition of $\bar{T}_{R}^{\pm}$ from (2.36). To show the local Lipschitz bound (6.4), it suffices to show (6.4) for a fixed $\bar{T}_{R}^{\pm}$, $R\in\mathfrak{R}$. For any $\bar{\rho},\hat{\rho}\in X_{p}$, by (2.36) it holds for all $(x,\xi)\in\Pi$ that

[TABLE]

Taking $\|\cdot\|_{X_{p}}$-norm for the above, and applying Young’s inequality, we have

[TABLE]

Hence, in this case $T^{\pm}$ also satisfy (T1)–(T4).

Notice that, for any $p\in[1,\infty)$, $A$ is the generator of a contractive (analytic) semigroup on $X_{p}=(L^{1}\cap L^{p})(\Gamma,dm)$, and denote $\{e^{tA}\}_{t\geq 0}$ the correspondent semigroup. This semigroup is in fact mass-conserving, in the sense that

[TABLE]

Moreover, by the estimate of heat potentials, the semigroup has the following $L^{p}\to L^{r}$ bound for any $1\leq p\leq r\leq\infty$:

[TABLE]

We say a $X_{p}$-valued function $\rho\in C([0,t_{0});X_{p})$ for some $t_{0}\in(0,\infty]$ is a mild solution to the equation (6.1) if

[TABLE]

It is a local solution if $t_{0}<\infty$, and a global solution if $t_{0}=\infty$. Notice also if a mild solution possesses higher regularity, namely, $\rho\in C((0,t_{0});W^{2,p}(\Gamma))\cap C^{1}((0,t_{0});X_{p})$, then it is a strong solution. If it is $C^{1}$ in time, $C^{2}$ in space (with respect to $x\in\mathbb{X}$), then it is a classical solution.

The main result of this section, which applies to both Fokker Planck equation (3.6) and mean field limit equation (2.34), is stated as follows.

Proposition 6.1**.**

Let $A$ be from (6.2) and $T^{\pm}$ satisfy (T1)–(T4). For every nonnegative $L^{1}(\Gamma)$ initial data $\rho_{0}$, the Cauchy problem (6.1) admits a unique global mild solution. The solution is nonnegative, and mass-conserving, in the sense that

[TABLE]

The solution has the following regularity for any $p\in[1,\infty)$:

[TABLE]

If additionally, for some $p\in[1,\infty)$ $\rho_{0}\in L^{p}(\Gamma)$, then $\rho\in C([0,\infty);L^{p}(\Gamma))$; if $\rho_{0}\in W^{2,p}(\Gamma)$, then $\rho\in C([0,\infty);W^{2,p}(\Gamma))\cap C^{1}([0,\infty);L^{p}(\Gamma))$.

Proof.

Let us first address the local well-posedness of the problem (6.1) with general $X_{p}$ initial data with any $p\in[1,\infty)$. Since the maps $T^{\pm}$ are assumed to be locally Lipschitz in $X_{p}$, by Picard-Lindelöf theorem for Banach-valued ODEs, a unique local mild solution $\rho\in C([0,t_{0});X_{p})$ exists for every $\rho_{0}\in X_{p}$.

We will first prove existence, uniqueness, non-negativity and mass conservation of a global mild solution $\rho\in C([0,\infty);X_{p}^{+})$ for every initial data $\rho_{0}\in X_{p}^{+}$, $p\in[1,\infty)$. Proposition 6.1 is then the special case $p=1$. Given $\rho_{0}\in X_{p}^{+}$, let $\rho\in C([0,t_{0});X_{p})$ be the mild solution of (6.1), guaranteed by the local well-posedness result earlier, and we now show $\rho$ is nonnegative. To this end, we consider (6.1) with $T^{+}(\rho)$ replaced by any $f\in C([0,\infty);X_{p}^{+})$, that is,

[TABLE]

We will first show that the above admits a unique nonnegative local mild solution. Existence and uniqueness of a local solution $\mu\in C([0,t_{0});X_{p})$ follows by the local Lipschitz bound (6.4) for $T^{-}$, and the same argument earlier. To prove its non-negativity, observe that for each $\xi\in G$, by (T3), (6.3) with $f$ in place of $T^{+}\rho$, and $f\geq 0$, $u_{\xi}(t,x)=\rho(t,x,\xi)$ is a super-solution of the following semilinear equation:

[TABLE]

$\beta(t)$ is finite on $[0,t_{0})$ because $\mu\in C([0,t_{0});L^{1}(\Gamma))$. By the comparison principle for semilinear equations ($v\equiv 0$ is a solution to above), we have $u_{\xi}(t)\geq 0$, i.e., $\mu(t)\in X_{p}^{+}$, for all $t\in[0,t_{0})$. In conclusion, we have shown that the local solution $\mu$ of (6.9) is in the space $C([0,t_{0});X^{+}_{p})$, if the initial data $\rho_{0}\in X_{p}^{+}$.

Now let $Y_{t_{0}}=C([0,t_{0}];X_{p}^{+})$ for $t_{0}>0$ small, and $\mathcal{P}(\cdot;\rho_{0}):Y_{t_{0}}\to Y_{t_{0}}$ be the map sending $f\in Y_{t_{0}}$ to the solution $\mu\in Y_{t_{0}}$ of (6.9). Subsequently, define $\mathcal{Q}(g;\rho_{0})=\mathcal{P}(T^{+}(g);\rho_{0})$, which again is a map from $Y_{t_{0}}\to Y_{t_{0}}$, because $T^{+}(g)\in X_{p}^{+}$ by (T2). Again, by choosing $t_{0}>0$ sufficiently small, using the Lipschitz bound from (T1) for $T^{+}$, one can show that $\mathcal{Q}(\cdot;\rho_{0})$ is contractive in $Y_{t_{0}}$, and hence has a unique fixed point $\mu\in Y_{t_{0}}$. Since the fixed point also satisfies the equation (6.1), by uniqueness $\mu$ agrees with the unique (local) mild solution $\rho$ of (6.1). Thus, $\rho\in C([0,t_{0}];X_{p}^{+})$. From here, the conservation of mass (6.8) holds for all $t\in[0,t_{0}]$, as a direct consequence of (6.5), (6.7) and (T4).

To show that the solution can be extended globally, it suffices to establish a global bound for the $\|\cdot\|_{X_{p}}$-norm of solutions. By the bound from (T1), (6.7), contractive nature of the semigroup $\{e^{tA}\}_{t\geq 0}$ and conservation of mass, a mild solution $\rho(t)$ to (6.1) satisfies

[TABLE]

where $C^{\prime}=2C\max\{1,\|\rho_{0}\|_{L^{1}(\Gamma)}\}$. Applying Grönwall lemma, we then have $\|\rho(t)\|_{X_{p}}\leq\|\rho_{0}\|_{X_{p}}e^{C^{\prime}t}$. With this bound, the solution can be extended indefinitely in time, and thus a unique global mild solution exists. To conclude, we have just shown if $\rho_{0}\in X_{p}^{+}$ for some $p\in[1,\infty)$, then (6.1) admits a unique global mild solution $\rho\in C([0,\infty);X_{p}^{+})$.

Now we address the regularity issue. From now on let us assume $\rho_{0}\in X_{1}^{+}$, and let $\rho\in C([0,\infty);X_{1}^{+})$ be the unique global mild solution guaranteed by the previous existence and uniqueness result. We now prove that $C((0,\infty);X_{p}^{+})$ for any $p\in[1,\infty)$. Using an induction argument, this follows by the following claim: if $\rho_{N}\in C((0,\infty);X_{p}^{+})$ for some $p\geq 1$, and $r>p$ is such that $p^{-1}-r^{-1}=D^{-1}$ (recall $D$ is the spatial dimension of $\mathbb{X}$), then $\rho_{N}\in C((0,\infty);X_{r}^{+})$. Indeed, regarding $\rho_{N}$ as a solution of (6.1) starting at $t=\delta_{0}$ for some fixed $\delta_{0}>0$, by (T1), (6.7) and the $L^{p}\to L^{r}$ bound (6.6), we have for all $t>0$ that

[TABLE]

with some constant $C=C(D,p,r,\|\rho_{0}\|_{L^{1}(\Gamma)})$. Now regarding $\rho_{N}$ as a solution of the forward equation (6.1) starting at $t=2\delta_{0}$, with initial data $\rho_{N}(2\delta_{0})\in X_{r}^{+}$, the previous existence and uniqueness result then implies $\rho_{N}\in C([2\delta_{0},\infty);X_{r}^{+})$. Since $\delta_{0}>0$ is arbitrary, it follows $\rho_{N}\in C((0,\infty);X_{r}^{+})$, which is our claim.

We have shown $\rho_{N}\in C((0,\infty);L^{p}(\Gamma))$ for any $p\in[1,\infty)$. By (T1) we have $T^{\pm}(\rho)\in C((0,\infty);L^{p}(\Gamma))$. By (6.3), since the forcing $T^{\pm}(\rho)$ is in $L^{p}(\mathbb{X})$, the standard parabolic regularity results implies $u_{\xi}(t,x)=\rho(t,x,\xi)\in C((0,\infty);W^{2,p}(\mathbb{X}))\cap C^{1}((0,\infty);L^{p}(\mathbb{X}))$ for every $\xi\in G$, e.g. see Theorem 7.22 of [27] or Section IV.3 of [26]. Hence, it follows $\rho_{N}\in C((0,\infty);W^{2,p}(\Gamma))\cap C^{1}((0,\infty);L^{p}(\Gamma))$, for any $p\geq 1$. Of course, if $\rho_{N}(0)\in W^{2,p}(\Gamma)$ initially, we may replace the open time interval $(0,\infty)$ above by the closed one $[0,\infty)$. This finishes the proof.

∎

Finally, we give a proof to the strict positivity of solutions correspondent to positive initial data (6.1), for the case of torus $\mathbb{X}=\mathbb{T}^{D}$.

Proposition 6.2**.**

Let $\mathbb{X}=\mathbb{T}^{D}$, $\rho_{0}\in L^{1}(\Gamma)$ with $\rho_{0}(x,\xi)\geq 0$, and let $\rho$ be the corresponding solution of (6.1). For any $\xi\in\{1,\dots,n\}$, if $\rho_{0}(\cdot,\xi)$ is positive on a set of positive measure, then $\inf_{(t,x)\in[t_{0},t_{1}]\times\Pi}\rho(t,x,\xi)>0$ for any $0<t_{0}<t_{1}$. Moreover, if $\inf_{(x,\xi)\in\Gamma}\rho_{0}(x,\xi)>0$, then $\inf_{(t,x,\xi)\in[0,t_{1}]\times\Pi}\rho(t,x,\xi)>0$ for any $t_{1}>0$.

Proof.

As shown in the previous proof, for each $\xi\in G$, $u_{\xi}(t,x)=\rho(t,x,\xi)$ is a super-solution of (6.10). By the conservation of mass, $\beta(t)=\beta_{0}$ for some $\beta_{0}>0$ depending on $\|\rho_{0}\|_{L^{1}(\Gamma)}$. Hence, we have

[TABLE]

The comparison principle implies

[TABLE]

where $H_{\xi}$ is the heat kernel on $\mathbb{X}=\mathbb{T}^{D}$ for the uniformly elliptic operator $A_{\xi}$. If $\rho_{0}(\cdot,\xi)$ is positive on a set of positive measure, then $H_{\xi}(t,x)*\rho_{0}(\cdot,\xi)$ is strictly positive on any compact subset of $(0,\infty)\times\mathbb{T}^{D}$. Moreover, if $\rho_{0}(\cdot,\xi)$ is bounded below by $\epsilon>0$, then $H_{\xi}(t,x)*\rho_{0}\geq\epsilon$ also holds. ∎

Now we give an improved positivity result using the particular structure of the forward equations (3.6) and (2.34). Given a subset of species $V_{0}\subset\mathfrak{S}$, we inductively define an increasing family of sets $V_{n}\subset\mathfrak{S}$ by

[TABLE]

(Recall the notation (2.3) for $R^{\pm}$.) That is, $V_{n+1}$ is the set of all chemical species (types) that are the products of reactions having inputs only from $V_{n}$. Then, we define the set

[TABLE]

We call this set $\overline{V_{0}}\subset\mathfrak{S}$ the closure of $V_{0}$ under the reaction network dynamics. Finally, we say that an initial density $\bar{\rho}_{0}(x,\xi)\geq 0$ on $\mathbb{X}\times\mathfrak{S}$ is propagating if $\overline{V_{0}}=\mathfrak{S}$ when

[TABLE]

Proposition 6.3**.**

Let $\bar{\rho}$ be a probability density on $\Pi$ satisfying (2.34) with initial condition $\bar{\rho}_{0}(x,\xi)$. Let $V_{0}=V_{0}(\bar{\rho}_{0})$ be defined by (6.13). Then for all $\xi\in\overline{V_{0}}$, $\inf_{(t,x)\in C}\bar{\rho}(t,x,\xi)>0$ for any compact set $C\subset(0,\infty)\times\mathbb{T}^{d}$. Moreover, for all $\xi\in\mathfrak{S}\setminus\overline{V_{0}}$, $\bar{\rho}(t,x,\xi)=0$ holds for all $t\geq 0$, $x\in\mathbb{T}^{d}$.

Proof.

By Proposition 6.2, we know that for any compact set $C\subset(0,\infty)\times\mathbb{T}^{d}$

[TABLE]

holds for all $\xi\in V_{0}$. Now, proceeding inductively, suppose that (6.14) holds for all $\xi\in V_{n}$, for some $n\geq 0$, with $V_{n}$ defined via (6.12). Then if $\xi^{\prime}\in V_{n+1}$, the definition of operator $\bar{T}^{+}_{R}$ in (2.36) implies

[TABLE]

holds for any compact set $C\subset(0,\infty)\times\mathbb{T}^{d}$. As in (6.11), this implies that for any compact $C^{\prime}\subset(0,\infty)\times\mathbb{T}^{d}$, there is a constant $c_{1}>0$ such that

[TABLE]

holds for all $(t,x)\in C^{\prime}$. This and the maximum principle implies that the condition (6.14) also holds for $\xi^{\prime}$. Since (6.14) holds for all such $\xi^{\prime}\in V_{n+1}$, we conclude by induction on $n$ that (6.14) holds for all $\xi\in\overline{V_{0}}$. ∎