Quasi-Stationary Distributions and Resilience: What to get from a sample?

TL;DR

This paper investigates how to estimate the resilience of multi-species birth-and-death processes by linking macroscopic dynamical stability to microscopic fluctuations, especially in large but finite populations.

Contribution

It establishes relations between resilience and process fluctuations, providing estimators to assess resilience from observed data on large stochastic systems.

Findings

Relations between resilience and microscopic fluctuations are derived.

Estimators for resilience are developed for different time scales.

The approach connects stochastic process behavior with dynamical system stability.

Abstract

We study a class of multi-species birth-and-death processes going almost surely to extinction and admitting a unique quasi-stationary distribution (qsd for short). When rescaled by $K$ and in the limit $K\to+\infty$, the realizations of such processes get close, in any fixed finite-time window, to the trajectories of a dynamical system whose vector field is defined by the birth and death rates. Assuming that this dynamical has a unique attracting fixed point, we analyzed in a previous work what happens for large but finite $K$, especially the different time scales showing up. In the present work, we are mainly interested in the following question: Observing a realization of the process, can we determine the so-called engineering resilience? To answer this question, we establish two relations which intermingle the resilience, which is a macroscopic quantity defined for the dynamical system, and the fluctuations of the process, which are microscopic quantities. Analogous relations are well known in nonequilibrium statistical mechanics. To exploit these relations, we need to introduce several estimators which we control for times between $\log K$ (time scale to converge to the qsd) and $\exp(K)$ (time scale of mean time to extinction).

Equations

$$\underset{K\to +\infty }{\lim }{\text{mathds}P}_{\lfloor K{x}_0\rfloor }\left(\underset{0\le t\le t_{\text{scaleto}H3pt}}{\sup }\text{dist}\left(\frac{\text{ushort}{N}^K(t)}{K},\text{ushort}x(t)\right)>\epsilon \right)=0$$

$$\frac{\mathrm{d}\text{ushort}x}{\mathrm{d}t}=\text{ushort}B(\text{ushort}x)-\text{ushort}D(\text{ushort}x)$$

$$\text{ushort}X=\text{ushort}B-\text{ushort}D$$

$$M^{\ast }=D\text{ushort}X(\text{ushort}{x}^{\text{scaleto}\ast 3pt}).$$

$${\rho }^{\text{scaleto}\ast 2.8pt}=-\sup \{\mathrm{Re}(z):z\in \mathrm{Sp}(M^{\ast })\}$$

$${\mu }_p^K=\int n_p\mathrm{d}\nu {}_{\text{scaleto}K3.5pt}(\text{ushort}n),p=1,\dots ,d.$$

$$\Sigma^{{\scriptscriptstyle K}}_{p,q}(\tau)={\mathds{E}}_{\nu\!_{\scaleto{K}{2.8pt}}}\!\left[\big{(}N_{p}^{{\scriptscriptstyle K}}(\tau)-\mu_{p}^{{\scriptscriptstyle K}}\big{)}\big{(}N_{q}^{{\scriptscriptstyle K}}(0)-\mu_{q}^{{\scriptscriptstyle K}}\big{)}\right],\;p,q\in\{1,\ldots,d\}.$$

$$\Sigma^{{\scriptscriptstyle K}}(\tau)=\operatorname{e}^{\tau M^{*}}\Sigma^{{\scriptscriptstyle K}}(0)+\mathcal{O}\big{(}\sqrt{K}\,\big{)}.$$

$$e^{\tau M^{\ast }}={\Sigma }^K(\tau ){\Sigma }^K(0{)}^{-1}+\mathcal{O}\left(\frac{1}{\sqrt{K}}\right).$$

$$S_p^{\text{ushort}\mu }(T,K)=\frac{1}{T}{\int }_0^T{N}_p^K(s)\mathrm{d}s$$

$$S_{p,q}^C(T,\tau ,K)=$$

$$\displaystyle\qquad\frac{1}{T\!-\!\tau}\int_{0}^{T-\tau}\!\!\big{(}N_{p}^{{\scriptscriptstyle K}}(s+\tau)\!-\!S^{{\scriptscriptstyle\ushort{\mu}}}_{p}(T,K)\big{)}\big{(}N_{q}^{{\scriptscriptstyle K}}(s)\!-\!S^{{\scriptscriptstyle\ushort{\mu}}}_{q}(T,K)\big{)}\hskip 1.0pt\mathrm{d}s.$$

$$\displaystyle\left|{\mathds{E}}_{\ushort{n}}\big{[}S^{{\scriptscriptstyle\ushort{\mu}}}_{p}(T,K)\big{]}-\mu_{p}^{{\scriptscriptstyle K}}\right|\leq$$

$$\displaystyle\qquad C\big{(}K+\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}\big{)}\left(\frac{1+\log K}{T}+\operatorname{e}^{-c\,(\scaleto{\|}{7pt}\ushort{n}\scaleto{\|}{7pt}_{\scaleto{1}{3pt}}\wedge K)}+\,T\operatorname{e}^{-\,c^{\prime}K}\right)$$

$$e^{\tau M_{\mathrm{emp}}^{\ast }(T,\tau ,K)}=S^C(T,\tau ,K)S^C(T,0,K{)}^{-1}.$$

$$\rho^{\scaleto{*}{3pt}}_{\scriptscriptstyle{\mathrm{emp}}}(T,\tau,K)=-\sup\big{\{}\operatorname{Re}(z):z\in\mathrm{Sp}\big{(}M^{*}_{\scriptscriptstyle{\mathrm{emp}}}(T,\tau,K)\big{)}\big{\}}.$$

$$\big{|}\rho^{\scaleto{*}{2.8pt}}_{\scriptscriptstyle{\mathrm{emp}}}(T,\tau,K)-\rho^{\scaleto{*}{2.8pt}}\big{|}\leq\mathcal{O}(1)\left(\frac{K^{2}}{\sqrt{T}}+\frac{1}{\sqrt{K}}\right)$$

$$\big{|}\rho^{\scaleto{*}{2.8pt}}_{\scriptscriptstyle{\mathrm{emp}}}(T,\tau,K)-\rho^{\scaleto{*}{2.8pt}}\big{|}\leq\frac{\mathcal{O}(1)}{\sqrt{K}}.$$

$${\mathcal{D}}_{p,p}^K=K{B}_p(\text{ushort}{x}^{\text{scaleto}\ast 3pt})=K{D}_p(\text{ushort}{x}^{\text{scaleto}\ast 3pt}).$$

$$M^{*}\Sigma^{{\scriptscriptstyle K}}+\Sigma^{{\scriptscriptstyle K}}{M^{*}}^{\intercal}+2\mathcal{D}^{{\scriptscriptstyle K}}=\mathcal{O}\big{(}\sqrt{K}\,\big{)}.$$

$${\rho }^{\text{scaleto}\ast 2.8pt}=\frac{K(B(\text{ushort}{x}^{\text{scaleto}\ast 3pt})+D(\text{ushort}{x}^{\text{scaleto}\ast 3pt}))}{2{\Sigma }^K}+\mathcal{O}\left(\frac{1}{\sqrt{K}}\right).$$

$$S^{{\scriptscriptstyle\mathcal{D}}}(T,K)=\frac{1}{T}\big{(}\mathrm{number\;of\;births\;up\;to\;time}\;T\big{)}.$$

$$\left|{\mathds{E}}_{\ushort{n}}\big{[}S^{{\scriptscriptstyle\mathcal{D}}}(T,K)\big{]}-KB(\ushort{x}^{\scaleto{*}{3pt}})\right|$$

$$S^{{\scriptscriptstyle\Sigma}}(T,K)=\frac{1}{T}\int_{0}^{T}\big{(}N^{{\scriptscriptstyle K}}(s)-S^{{\scriptscriptstyle\ushort{\mu}}}(T,K)\big{)}\big{(}N^{{\scriptscriptstyle K}}(s)-S^{{\scriptscriptstyle\ushort{\mu}}}(T,K)\big{)}\hskip 1.0pt\mathrm{d}s.$$

$$\frac{\mathrm{d}\text{ushort}x}{\mathrm{d}t}=\text{ushort}B(\text{ushort}x)-\text{ushort}D(\text{ushort}x)=\text{ushort}X(\text{ushort}x).$$

$$\text{scaleto}\Vert 9pt\text{ushort}x\text{scaleto}\Vert {9pt}_{\text{scaleto}14.pt}=\sum _{j=\text{scaleto}14.5pt}^{\text{scaleto}d5pt}{x}_j,\text{scaleto}\Vert 9pt\text{ushort}x\text{scaleto}\Vert {9pt}_{\text{scaleto}23.5pt}=\sqrt{\sum _{j=\text{scaleto}14.5pt}^{\text{scaleto}d5pt}{x}_j^{\text{scaleto}24pt}}.$$

$$\text{ushort}B(\text{ushort}{x}^{\text{scaleto}\ast 3pt})-\text{ushort}D(\text{ushort}{x}^{\text{scaleto}\ast 3pt})=\text{ushort}X(\text{ushort}{x}^{\text{scaleto}\ast 3pt})=\text{ushort}0.$$

$$⟨\text{ushort}X(\text{ushort}x),(\text{ushort}x-\text{ushort}{x}^{\text{scaleto}\ast 3pt})⟩\le -\beta \text{scaleto}\Vert 9pt\text{ushort}x\text{scaleto}\Vert {9pt}_{\text{scaleto}23.5pt}\text{scaleto}\Vert 10pt\text{ushort}x-\text{ushort}{x}^{\text{scaleto}\ast 3pt}\text{scaleto}\Vert {10pt}_{\text{scaleto}24pt}^{\text{scaleto}24pt}.$$

$$\hat{B}(y)=\underset{\text{scaleto}\Vert 8pt\text{ushort}x\text{scaleto}\Vert {8pt}_{\text{scaleto}13pt}=y}{\sup }\sum _{j=\text{scaleto}14.5pt}^d{B}_j(\text{ushort}x),\hat{D}(y)=\underset{\text{scaleto}\Vert 8pt\text{ushort}x\text{scaleto}\Vert {8pt}_{\text{scaleto}13pt}=y}{\inf }\sum _{j=\text{scaleto}14.5pt}^d{D}_j(\text{ushort}x)$$

$$F(y)=\frac{\hat{B}(y)}{\hat{D}(y)}.$$

$$\underset{\text{ushort}x\in \text{mathds}{R}_{\text{scaleto}+3.5pt}^{\text{scaleto}d3.5pt}}{\inf }\underset{1\le j\le d}{\inf }\frac{D_j(\text{ushort}x)}{{\sup }_{1\le \ell \le d}{x}_{\ell }}>\xi .$$

Full text

Keywords: birth-and-death process, dynamical system, engineering resilience, quasi-stationary distribution, fluctuation-dissipation relation, empirical estimators.

1. 1 Introduction and main results

2. 1.1 Context and setting

3. 1.2 Main results

4. 1.3 A “fluctuation-dissipation” approach

5. 1.4 Standing assumptions

6. 1.5 A numerical example

7. 1.6 Organization of the paper

8. 2 Time evolution of moments of the process and moments of the QSD

9. 2.1 Time evolution of moments starting from anywhere

10. 2.2 Moments estimates for the qsd

11. 3 Controlling time averages of the estimators

12. 4 Fluctuation and correlation relations

13. 4.1 Proof of Theorem 1.1

14. 4.2 Proof of Theorem 1.2

15. 5 Variance estimates for the estimators

16. A Proof of the two variance estimates

17. A.1 Starting from anywhere: proof of Proposition 5.1

18. A.2 Starting from the qsd: proof of Corollary 5.2

19. B Counting the number of births

20. C Gaussian limit for the rescaled qsd

Titre en français :

Distributions quasi-stationnaires et résilience : que peut-on obtenir des données ?

Résumé en français :

Nous étudions une classe de processus de naissance et mort avec plusieurs espèces dans la situation où l’extinction est certaine et la distribution quasi-stationnaire est unique. Si on fixe un intervalle de temps fini et qu’on normalise les réalisations d’un tel processus par un paramètre d’échelle $K$, elles deviennent arbitrairement proches, dans la limite $K\to+\infty$, des trajectoires d’un certain système dynamique dont le champ de vecteurs est défini à partir des taux de naissance et mort. Quand le système dynamique admet un seul point fixe attractif, nous avons précédemment analysé le comportement du processus pour des valeurs de $K$ finies et pour des temps finis, c’est-à-dire, le comportement intermédiaire entre les deux comportements limites évoqués ci-dessus. La question principale qui nous intéresse dans le présent article est la suivante : si on observe une réalisation du processus, pouvons-nous estimer la résilience au sens de l’ingénieur (engineering resilience) ? Pour répondre à cette question, nous démontrons deux relations entremêlant la résilience, qui est une quantité macroscopique définie pour le système dynamique sous-jacent, et les fluctuations du processus, qui sont elles des quantités microscopiques. De tels genres de relations sont bien connues en mécanique statistique hors d’équilibre. Afin d’exploiter ces relations nous introduisons plusieurs estimateurs empiriques que nous parvenons à contrôler pour des temps entre $\log K$, qui est l’échelle de temps pour observer la convergence vers la distribution quasi-stationnaire, et $\exp(K)$, qui est l’échelle du temps moyen d’extinction.

Mots clés :

processus de naissance et mort, systèmes dynamiques, résilience, distribution quasi-stationnaire, relation de fluctuation-dissipation, estimateurs empiriques.

1 Introduction and main results

1.1 Context and setting

The ability of an ecosystem to return to its reference state after a perturbation stress is given by its resilience, a concept pioneered by Holling. Resilience has several faces and multiple definitions [5]. In the traditional theoretical setting of dynamical systems, that is, differential equations, one of them is the so-called engineering resilience. It is concerned with what happens in the vicinity of a fixed point (equilibrium state) of the system, and is given by minus the real part of the dominant eigenvalue of the Jacobian matrix evaluated at the fixed point. It can also be defined as the reciprocal of the characteristic return time to the fixed point after a (small) perturbation. In this paper, we are interested in how to determine the engineering resilience from the data. But which data? The drawback of the notion of engineering resilience is that we do not observe population densities governed by differential equations. Instead, we count individuals which are subject to stochastic fluctuations. Can we nevertheless infer the resilience? The subject of this paper is to show that this is possible in the framework of birth-and-death processes which are, in a sense made precise below, close to the solutions of a corresponding differential equation, at certain time and population size scales.

Let us now describe our framework. We consider a population made of $d$ species interacting with one another. Suppose that the state of the process at time $t$, which we denote by $\ushort{N}^{{\scriptscriptstyle K}}(t)=(N_{\scaleto{1}{4.5pt}}^{{\scriptscriptstyle K}}(t),\ldots,N_{\scaleto{d}{4.5pt}}^{{\scriptscriptstyle K}}(t))$, is $\ushort{n}=(n_{\scaleto{1}{4pt}},\ldots,n_{i},\ldots,n_{\scaleto{d}{4pt}})\in{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, where $n_{i}$ is the number of individuals of the $i$th species. Then the rate at which the population increases (respectively decreases) by one individual of the $j$th species is $\scaleto{K}{6pt}B_{j}(\scaleto{\frac{\ushort{n}}{K}}{11pt})$ (respectively $\scaleto{K}{6pt}D_{j}(\scaleto{\frac{\ushort{n}}{K}}{11pt})$), where $K$ is a scaling parameter. Under the assumptions we will make, the process goes extinct, i.e., $\ushort{0}$ is an absorbing state, with probability one. There are two limiting regimes for the behavior of this process. The first one is to fix $K$ and let $t$ tend to infinity, which leads inevitably to extinction. The second one consists in fixing a time horizon and letting $K$ tend to $+\infty$, after having rescaled the process by $K$. In this limit, the behavior of the rescaled process is governed by a certain differential equation. More precisely, given any $0<t_{\scaleto{H}{4pt}}<+\infty$ and any $\varepsilon>0$ and $x_{0}\in\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, we have

[TABLE]

where $\textup{dist}(\cdot,\cdot)$ is the Euclidean distance in $\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, and $\ushort{x}(t)$ is the solution of the differential equation in $\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$

[TABLE]

with initial condition $\ushort{x}_{0}$. We refer to [4, Chapter 11] for a proof. We use the notations $\ushort{x}=(x_{\scaleto{1}{4pt}},\ldots,x_{\scaleto{d}{4pt}})$, $\ushort{B}(\ushort{x})=(B_{\scaleto{1}{4pt}}(\ushort{x}),\ldots,B_{\scaleto{d}{4pt}}(\ushort{x}))$, and so on and so forth. We will make further assumptions (see Subsection 1.4) on the birth and death rates to be in the following situation. The vector field

[TABLE]

has a unique attracting fixed point $\ushort{x}^{\scaleto{*}{3pt}}$ (lying in the interior of $\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$). We denote by $M^{*}$ its differential evaluated at $\ushort{x}^{\scaleto{*}{3pt}}$, namely

[TABLE]

We then define the (engineering) resilience as

[TABLE]

where \mathrm{Sp}\big{(}M^{*}\big{)} denotes the spectrum (set of eigenvalues) of the matrix $M^{*}$. Under our assumptions, we have $\rho^{\scaleto{*}{2.8pt}}>0$.

We can now formulate more precisely the goal of this paper. Given a finite-length realization of the process $(\ushort{N}^{{\scriptscriptstyle K}}(t),t\leq T)$, with large, but finite $K$, we want to build an estimator for $\rho^{\scaleto{*}{2.8pt}}$. To this end, we need a good understanding of the behavior of the Markov process $(\ushort{N}^{{\scriptscriptstyle K}}(t))$ in an intermediate regime between the two limiting regimes described above. This was done in a previous work of ours [3], and this can be roughly summarized as follows. All states $\ushort{n}\neq\ushort{0}$ are transient and $\ushort{0}$ is absorbing, hence the only stationary distribution is the Dirac measure sitting at $\ushort{0}$. The mean time to extinction behaves like $\exp(\Theta(K))$. (We recall Bachmann-Landau notations below.) If we start in the vicinity of the state $\ushort{n}^{\scaleto{*}{3pt}}=\lfloor K\ushort{x}^{\scaleto{*}{3pt}}\rfloor$, that is, if the initial state has its coordinates of size of order $K$, then either the process wanders around $\ushort{n}^{\scaleto{*}{3pt}}$ or it gets absorbed at $\ushort{0}$. More precisely, there is a unique quasi-stationary distribution (qsd, for short) $\nu\!_{\scaleto{K}{3.5pt}}$ which describes the statistics of the process conditioned to be non-extinct before time $t$. Without this conditionning, the law of the process at time $t$ is well approximated by a mixture of the Dirac measure at $\ushort{0}$ and the qsd $\nu\!_{\scaleto{K}{3.5pt}}$, for times $t\in\left[cK\log K,\exp(\Theta(K))\right]$, where $c>0$, in the sense that the total variation distance between them is exponentially small in $K$, provided that $K$ is large enough. We will rely on these results that will be recalled precisely later on. We will also need to prove further properties.

1.2 Main results

To estimate the engineering resilience $\rho^{\scaleto{*}{2.8pt}}$, we will establish a matrix relation involving $M^{*}$. Let $\ushort{\mu}^{{\scriptscriptstyle K}}=(\mu_{1}^{{\scriptscriptstyle K}},\ldots,\mu_{d}^{{\scriptscriptstyle K}})$ be the vector of species sizes averaged with respect to $\nu\!_{\scaleto{K}{3.5pt}}$, that is,

[TABLE]

For each $\tau\geq 0$, define the matrix

[TABLE]

In Section 4.1, we will prove the following result.

Theorem 1.1**.**

For all $\tau\geq 0$ we have

[TABLE]

Some comments are in order. If $\tau$ is equal to, say, $1/K$ then the estimate becomes useless. More generally, if $\tau$ is too small then $\operatorname{e}^{\tau M^{*}}$ is too close to the identity matrix. Moreover, we will show later on that $\ushort{\mu}^{{\scriptscriptstyle K}}$ and $\Sigma^{{\scriptscriptstyle K}}(\tau)$ are of order $K$. Hence the estimate becomes irrelevant if $\tau$ becomes proportional to $\log K$. Indeed, without knowing the constant of proportionality, $\operatorname{e}^{\tau M^{*}}\Sigma^{{\scriptscriptstyle K}}(0)$ can be of the same order than the error term.

Before proceeding further, we recall the following classical Bachmann-Landau notations.

Notations**.**

Given $a\in\real$, the symbol $\mathcal{O}(K^{a})$ stands for any real-valued function $f(K)$ such that there exists $C>0$ and $K_{0}>0$ such that, for any $K>K_{0}$, $|f(K)|\leq CK^{a}$. Note in particular that $\mathcal{O}(1)$ will always mean a strictly positiveconstant that we don’t want to specify. Sometimes, we will also use the symbol $\Theta(K^{a})$ stands for any real-valued function $f(K)$ such that there exist $C_{1},C_{2}>0$ and $K_{0}>0$ such that, for any $K>K_{0}$, $C_{1}K^{a}\leq f(K)\leq C_{2}K^{a}$. One can naturally generalize $\Theta(K^{a})$ to vector-valued functions. For instance, for $\ushort{n}\in\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$ we write $\ushort{n}=\Theta(K^{a})$ if $n_{i}=\Theta(K^{a})$ for $i=1,\ldots,d$.

Relation (1.2) allows to determine $M^{*}$. Indeed, we have

[TABLE]

This formula suggests that in order to estimate $M^{*}$, we need estimators for $\Sigma^{{\scriptscriptstyle K}}(0)$ and $\Sigma^{{\scriptscriptstyle K}}(\tau)$. Given a finite-length realization of \big{(}\ushort{N}^{{\scriptscriptstyle K}}(t),0\leq t\leq T\big{)} up to some time $T>0$, we define estimators for $\mu_{p}^{{\scriptscriptstyle K}}$ and $\Sigma^{{\scriptscriptstyle K}}_{p,q}(\tau)$, for $0\leq\tau<T$, $p,q\in\{1,\ldots,d\}$, $K\in\mathbb{N}^{*}$ by

[TABLE]

and

[TABLE]

Under suitable conditions on $\ushort{n}$, $K$ and $T$, $S^{{\scriptscriptstyle\ushort{\mu}}}(T,K)$ well approximates $\ushort{\mu}^{{\scriptscriptstyle K}}$. More precisely, we will prove an estimate of the following form (see Theorem 3.4 for a precise statement)

[TABLE]

for every $\ushort{n}\in{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, $p=1,\ldots,d$, where $c,c^{\prime},C$ are positive constants. We use the notation $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}=\sum_{i=1}^{d}n_{i}$. Let us comment on this bound. Roughly speaking, it can only be useful if $T$ is much smaller than $\exp(\mathcal{O}(1)K)$ if $\ushort{n}$ is, say, of order $K$. For instance, suppose that, for $K$ large enough, we want the bound to be $\Theta(K^{-a})$, for some $a>0$. One can check that this is possible if $\ushort{n}=\Theta(K)$ and $T=\Theta(K^{a+1}\log K)$. (Note in particular that, in this situation, we have a consistent estimator when $K\to+\infty$.) However, when $T$ becomes $\exp(\mathcal{O}(1)K)$ or larger, we know that {\mathds{E}}_{\ushort{n}}\big{[}S^{{\scriptscriptstyle\ushort{\mu}}}_{p}(T,K)\big{]}\approx 0, because with high probability, at this time scale the process is absorbed at $\ushort{0}$. This is the manifestation of the fact that the only stationary distribution is the Dirac measure at $\ushort{0}$. Consistently, our bound becomes very bad in that regime.

An estimate of the same kind holds for $S^{{\scriptscriptstyle C}}(T,\tau,K)$ which well approximates $\Sigma^{{\scriptscriptstyle K}}(\tau)$ in the appropriate ranges.

Remark 1.1**.**

It is possible to use discrete time instead of continuous time in the above averages. Indeed the key results (in particular Proposition 3.3) are obtained for discrete times.

We can now define the empirical matrix $M^{*}_{\scriptscriptstyle{\mathrm{emp}}}(T,\tau,K)$ by

[TABLE]

We will see later on that, in appropriate regimes, $S^{{\scriptscriptstyle C}}(T,0,K)$ is near $\Sigma^{K}(0)$ and $S^{{\scriptscriptstyle C}}(T,\tau,K)$ is near $\Sigma^{K}(\tau)$ (see Propositions 5.4 and 5.6). The matrix $\Sigma^{K}(0)$ is invertible as a covariance matrix of a non-constant vector and is $\Theta(K)$ (see Proposition 2.9). Then (1.2) implies that $\Sigma^{K}(\tau)$ is invertible and the same holds for $S^{{\scriptscriptstyle C}}(T,\tau,K)$. These remarks imply that the matrix $M^{*}_{\scriptscriptstyle{\mathrm{emp}}}$ is well defined.

We define the empirical resilience by

[TABLE]

Our main result (Theorem 5.7) is then the following.

Theorem**.**

For $\tau=\Theta(1)$, $\ushort{n}=\Theta(K)$ (initial state) and $0<T\ll\exp(\Theta(1)K)$, and $K$ large enough, we have

[TABLE]

with a probability larger than $1-1/K$. In particular, if $T\geq K^{5}$, we have

[TABLE]

Several comments are in order. The dependence on the initial state $\ushort{n}$ is somewhat hidden and involved in the fact that the estimates hold “with a probability larger than $1-1/K$”. Indeed, the estimate of this probability results from Chebychev inequality and variance estimates in which the process is started in $\ushort{n}$. What the symbol $\ll$ precisely means is not mathematically defined. It means that we need to consider $T$ “much smaller than something exponentially big in $K$”. Indeed, since we do not control explicitly the various constants appearing in exponential terms in $K$, we have to consider $T$ which varies on a scale smaller than $\exp(\Theta(1)K)$, for instance \exp\big{(}\Theta(1)\sqrt{K}\,\big{)}. The reader is invited to step through the proof of Theorem 5.7 for the more precise, but cumbersome bound we obtain.

1.3 A “fluctuation-dissipation” approach

The above estimator for the engineering resilience, based on (1.3), is valid for any $d$. In the case $d=1$ (only one species), we have another, simpler, estimator based on a “fluctuation-dissipation relation”. This relation is in fact true for any $d$ and of independent interest. Let $\mathcal{D}^{{\scriptscriptstyle K}}$ be the $d\times d$ diagonal matrix given by

[TABLE]

We have the following result. We write $\Sigma^{{\scriptscriptstyle K}}$ instead of $\Sigma^{{\scriptscriptstyle K}}(0)$, and the transpose of a matrix $M$ is denoted by $M^{\intercal}$.

Theorem 1.2**.**

We have

[TABLE]

This relation is proved in Section 4.2. For background on fluctuation-dissipation relations in Statistical Physics, we refer to [7, sections 2-3]. Note that the matrix $\Sigma^{{\scriptscriptstyle K}}$ is symmetrical, but in general the matrix $M^{*}$ is not (see [3]). Note also that each term in the left-hand side of (1.7) is of order $K$, as we will see below.

If $\Sigma^{{\scriptscriptstyle K}}$ and $\mathcal{D}^{{\scriptscriptstyle K}}$ are known, the matrix $M^{*}$ is not uniquely defined, except for $d=1$ (see for example [8]). For $d=1$, (1.7) easily gives the resilience since it becomes a scalar equation:

[TABLE]

Remark 1.2**.**

The quantity $K(B(\ushort{x}^{\scaleto{*}{3pt}})+D(\ushort{x}^{\scaleto{*}{3pt}}))=2KB(\ushort{x}^{\scaleto{*}{3pt}})$ is the average total jump rate $K\nu\!_{\scaleto{K}{3.5pt}}(B(\ushort{n}/K)+D(\ushort{n}/K))$ up to $\mathcal{O}(1)$ terms. This follows from a Taylor expansion of $B(\ushort{n}/K)+D(\ushort{n}/K)$ around $\ushort{x}^{\scaleto{*}{3pt}}$, Theorem 2.6 and Proposition 2.7.

In the case $d=1$, an estimator for $\mathcal{D}^{{\scriptscriptstyle K}}$ is

[TABLE]

In Section 5, we establish a bound for

[TABLE]

which depends on $T,K$ and $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}$, and is small in the relevant regimes. The estimator we use for $\Sigma^{{\scriptscriptstyle K}}$ is

[TABLE]

Again, we can control how well this estimator approximates $\Sigma^{{\scriptscriptstyle K}}$. This provides another estimator for $\rho^{\scaleto{*}{3pt}}$, with a controlled error.

1.4 Standing assumptions

The two (regular) vector fields $\ushort{B}(\ushort{x})$ and $\ushort{D}(\ushort{x})$ are given in $\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$. We assume that their components have second partial derivatives which are polynomially bounded. Obviously, we suppose that $B_{j}(\ushort{x})\geq 0$ and $D_{j}(\ushort{x})\geq 0$ for all $j=1,\ldots,d$ and $\ushort{x}\in\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$. A dynamical system in $\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$ is defined by the vector field $\ushort{X}(\ushort{x})=\ushort{B}(\ushort{x})-\ushort{D}(\ushort{x})$, namely

[TABLE]

For $\ushort{x}\in\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, we use the following standard norms:

[TABLE]

We now state our hypotheses.

1. H.1

The vectors fields $\ushort{B}$ and $\ushort{D}$ vanish only at $\ushort{0}$. 2. H.2

There exists $\ushort{x}^{\scaleto{*}{3pt}}$ belonging to the interior of $\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$ (fixed point of $\ushort{X}$ ) such that

[TABLE] 3. H.3

Attracting fixed point: there exist $\beta>0$ and $R>0$ such that $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}^{\scaleto{*}{3pt}}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}<R$, and for all $\ushort{x}\in\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$ with $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}<R$,

[TABLE] 4. H.4

The fixed point $\ushort{0}$ of the vector field $\ushort{X}$ is repelling (locally unstable). Moreover, on the boundary of $\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, the vector field $\ushort{X}$ points toward the interior (except at $\ushort{0}$). 5. H.5

Define

[TABLE]

and for $y>0$, let

[TABLE]

We assume that there exists $\,0<L<R\,$ such that $\sup_{y>L}F(y)<\scaleto{1/2}{10pt}$ and $\lim_{y\to+\infty}F(y)=0.$ 6. H.6

There exists $y_{0}>0$ such that $\int_{y_{0}}^{\infty}\widehat{D}(y)^{-1}\hskip 1.0pt\mathrm{d}y<+\infty$ and $y\mapsto\widehat{D}(y)$ is increasing on $\left[y_{0},+\infty\right[$. 7. H.7

There exists $\xi>0$ such that

[TABLE] 8. H.8

Finally, we assume that

[TABLE]

(By $\partial_{x_{j}}$ we mean the partial derivative with respect to $x_{j}$.)

Assumptions H.5 and H.6 ensure that the time for “coming down from infinity” for the dynamical system is finite. Together with H.3, this also implies that $\ushort{x}^{\scaleto{*}{3pt}}$ is a globally attracting stable fixed point. More comments on these assumptions can be found in [3].

1.5 A numerical example

We consider the two-dimensional vector fields

[TABLE]

and

[TABLE]

where all the coefficients are positive. This is a model of competition between two species of Lotka-Volterra type. We have taken

[TABLE]

Assumptions H.1 and H.4 are easily verified numerically. Assumptions H.5 and H.6 are satisfied because $\widehat{B}(y)\leq(a+b+e+f)y$ and $\widehat{D}(y)\geq(c\wedge h)y^{\scaleto{2}{4pt}}/4$. Concerning H.2, we checked numerically that there is a unique fixed point inside the positive quadrant, namely $\ushort{x}^{\scaleto{*}{3pt}}=(0.3567,1.4855)$. It remains to check H.3, namely that

[TABLE]

where

[TABLE]

We first checked that the numerator $N(\ushort{x})=\langle\ushort{X}(\ushort{x}),(\ushort{x}-\ushort{x}^{\scaleto{*}{3pt}})\rangle$ is negative and vanishes only at $\ushort{0}$ and $\ushort{x}^{\scaleto{*}{3pt}}$. It is easy to check that $N(\ushort{x})<0$ for $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}$ large enough. We have verified numerically that the only solutions of the equations $\partial_{x_{1}}N=\partial_{x_{2}}N=0$ in the closed positive quadrant are $\ushort{x}^{\scaleto{*}{3pt}}$ and $\ushort{z}=(0.1739,0.4361)$, with $N(\ushort{z})=-0.2852$, thus this is a negative local minimum. This implies that $N(\ushort{x})<0$ in the closed positive quadrant, except at $\ushort{0}$ and $\ushort{x}^{\scaleto{*}{3pt}}$ where it vanishes. This implies that $R\leq 0$ in the closed positive quadrant. It is easy to check that

[TABLE]

This implies that $R<0$ except perhaps at $\ushort{0}$ and $\ushort{x}^{\scaleto{*}{3pt}}$. Near $\ushort{0}$ we have by Taylor expansion

[TABLE]

and, since the vector $D^{\intercal}\ushort{B}(0)\,\ushort{x}^{\scaleto{*}{3pt}}$ has positive components, there exists $\varrho>0$ such that for all $\ushort{x}\in{}_{\scaleto{+}{5pt}}^{\scaleto{2}{4pt}}$

[TABLE]

If $\ushort{y}=\ushort{x}-\ushort{x}^{\scaleto{*}{3pt}}$ is small, we have by Taylor expansion (since $\ushort{X}(\ushort{x}^{\scaleto{*}{3pt}})=\ushort{0}$)

[TABLE]

One can check numerically that the two real eigenvalues of the symmetric matrix

[TABLE]

are strictly negative, the largest being numerically equal to $-0.786$. This completes the verification of hypothesis H.3.

Illustrating standard experiments on populations of cells or bacteria, we have chosen $K=10^{5}$ and simulated a unique realization of the process with $T=100$ which contains about $5.10^{7}$ jumps (cell divisions or deaths). The resilience computed from the vector field is numerically equal to $0.547$. We have computed $\rho^{\scaleto{*}{3pt}}_{\scriptscriptstyle{\mathrm{emp}}}(100,1,10^{5})$. The relative error, that is $|\rho^{\scaleto{*}{3pt}}_{\scriptscriptstyle{\mathrm{emp}}}(100,1,10^{5})-\rho^{\scaleto{*}{3pt}}|/\rho^{\scaleto{*}{3pt}}$, is equal to $0.022$.

Note that the situation we are interested in is completely different from standard statistical approach where one can repeat the experiments.

1.6 Organization of the paper

In Section 2, we will study the time evolution of the moments of the process and we will prove moment estimates for the qsd. In Section 3, we will obtain control on the large time behavior of averages for the process. In Section 4, we will prove the relations (1.2) and (1.7). In Section 5, we will apply these relations to obtain approximate expressions of the engineering resilience in terms of the covariance matrices for the qsd. From the results of Section 3, we will deduce variance bounds for the estimators (1.4), (1.5) and (1.8), starting either in the qsd or from an initial condition of order $K$.

2 Time evolution of moments of the process and moments of the QSD

2.1 Time evolution of moments starting from anywhere

The generator $\mathscr{L}_{{\scriptscriptstyle K}}$ of the birth and death process $\ushort{N}^{{\scriptscriptstyle K}}=(\ushort{N}^{{\scriptscriptstyle K}}(t),t\geq 0)$ is defined by

[TABLE]

where $\scaleto{\ushort{e}^{(\ell)}}{11.5pt}=(0,\ldots,0,1,0,\ldots,0)$, the $1$ being at the $\ell$-th position, and $f:{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}\to\real$ is a function with bounded support. We denote by $(S^{{\scriptscriptstyle K}}_{t},t\geq 0)$ the semigroup of the process $\ushort{N}^{{\scriptscriptstyle K}}$ acting on bounded functions, that is, for $f:{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}\to\real$, we have

[TABLE]

For $A>1$, let

[TABLE]

Notice that we will use either $\scaleto{\|}{9pt}\hskip 0.5pt\!\cdot\!\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}$ or $\scaleto{\|}{9pt}\hskip 0.5pt\!\cdot\!\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}$. They are of course equivalent but one can be more convenient than the other, depending on the context. We have the following result.

Theorem 2.1**.**

There exists a constant $C_{\mathrm{(\ref{borneexp})}}>0$ such that for $K$ large enough, the operator group $S^{{\scriptscriptstyle K}}_{1}$ extends to exponentially bounded functions and

[TABLE]

Proof.

Introduce the function $G_{{\scriptscriptstyle K}}$ defined on $[\hskip 1.0pty_{0},+\infty)$ by

[TABLE]

Assumption H.6 implies that $G_{\!{\scriptscriptstyle K}}$ is well defined and decreasing on $[y_{0},+\infty)$. We can define its inverse function on $(0,s_{0}]$ for $s_{0}>0$ small enough (independent of $K$). Take $0<\eta\leq s_{0}\wedge\frac{1-\,\operatorname{e}^{-1}}{4}$. Then there is a unique positive function $y_{{\scriptscriptstyle K}}$ defined by

[TABLE]

Note that $y_{{\scriptscriptstyle K}}(s)\geq y_{0}$ and $\lim_{s\downarrow 0}y_{{\scriptscriptstyle K}}(s)=+\infty$. Let

[TABLE]

Note that

[TABLE]

Using the Lipschitz continuity of $\widehat{D}$ (and then its differentiability almost everywhere) and (2.3), we obtain

[TABLE]

We now consider the function

[TABLE]

to which we apply Itô’s formula at time $t\wedge\mathcal{T}_{A}$. We get

[TABLE]

We have

[TABLE]

Note that

[TABLE]

It follows from H.5 that there exists a number $\zeta>y_{0}$ such that if $y>\zeta$, then $F(y)<(2\operatorname{e})^{-1}$.

If $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}<\zeta K$ we get

[TABLE]

For $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}\geq K(\zeta\vee y_{{\scriptscriptstyle K}}(t))$ we have

[TABLE]

since $\dot{\varphi}_{{\scriptscriptstyle K}}(t)=\eta K\widehat{D}(y_{{\scriptscriptstyle K}}(t))\varphi_{{\scriptscriptstyle K}}(t)$ and $\widehat{D}(\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}/K)\geq\widehat{D}(y_{{\scriptscriptstyle K}}(t))$.

Finally, for $\zeta K\leq\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}<Ky_{{\scriptscriptstyle K}}(t)$ we get

[TABLE]

We deduce that

[TABLE]

The result follows by letting $A$ tend to infinity and by monotonicity. ∎

We deduce moment estimates for the process which are uniform in the starting state, and in time, for times larger than $1$.

Corollary 2.2**.**

For all $t\geq 1$, the semi-group $(S_{t})$ maps functions of polynomially bounded modulus in bounded functions. In particular, for all $q\in\mathds{N}$, we have

[TABLE]

Proof.

We have

[TABLE]

since for all $x\geq 0$, $x^{q}\operatorname{e}^{-x}\leq q^{q}\operatorname{e}^{-q}$. Inequality (2.4) follows from Hölder’s inequality and Theorem 2.1. Let us now consider $t>1$. From the Markov property and by using the previous inequality, we deduce that

[TABLE]

The proof is finished. ∎

For times $t$ less than $1$, the moment estimates depends on the initial state.

Proposition 2.3**.**

For each integer $q$, there exists a constant $c_{q}>0$ such that for all $K>1$, $t\geq 0$ and $\ushort{n}\in{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$

[TABLE]

Proof.

We have only to study the case $t<1$, the other case being given in (2.4). We prove the result for $q$ even, namely $q=\scaleto{2}{6pt}q^{\prime}$. The result for $q$ odd follows from Cauchy-Schwarz inequality. Letting

[TABLE]

we have

[TABLE]

Using H.5 and the equivalence of the norms, we see that there exists a constant $c_{q^{\prime}}>0$ such that if $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}>c_{q^{\prime}}K$

[TABLE]

Moreover, we can take $c_{q^{\prime}}$ large enough such that for all $\ushort{n}$

[TABLE]

Applying Itô’s formula to $f_{q^{\prime}}$ we get as in the proof of Theorem 2.1

[TABLE]

(Recall that $\mathcal{T}_{A}$ is defined in (2.2).) The result follows by letting $A$ tend to infinity. ∎

2.2 Moments estimates for the qsd

Let us first recall (cf. [3]) that, under the assumptions of Section 1.4, there exists a unique qsd $\nu\!_{\scaleto{K}{3.5pt}}$ with support ${{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}\backslash\{\ushort{0}\}$. Further, starting from the qsd, the extinction time is distributed according to an exponential law with parameter $\lambda_{0}(K)$ satisfying (Theorem 3.2 in [3])

[TABLE]

where $d_{1}>d_{2}>0$ are constants independent of $K$. Recall also that for all $t>0$,

[TABLE]

where

[TABLE]

Finally, for all $f$ in the domain of the generator

[TABLE]

with the notation

[TABLE]

We use several notations from [3] that we now recall. Let

[TABLE]

For $\ushort{x}\in\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$ and $r>0$, $\mathcal{B}(\ushort{x},r)$ is the ball of center $\ushort{x}$ and radius $r$. We consider the sets

[TABLE]

where $\rho_{(4.2)}>0$ is a constant defined in [3, Corollary 4.2]. Note that since $\ushort{n}^{\scaleto{*}{3pt}}$ is of order $K$, we have $\Delta\subset\mathcal{D}$ for $K$ large enough. The first entrance time in $\Delta$ (resp. $\mathcal{D}$) will be denoted by $T_{\!\scaleto{\Delta}{5pt}}$ (resp. $T_{\mathcal{D}}$).

We first prove that the support of the qsd is, for large $K$, almost included in $\mathcal{D}$. (This will be important to control moments later on.)

Proposition 2.4**.**

There exists a constant $c_{\mathrm{(\ref{borneexterieur})}}>0$ such that for all $K$ large enough

[TABLE]

Proof.

We first recall two results from [3]. From Lemma 5.1 in [3], there exist $\gamma>0$ and $\delta\in(0,1)$ such that for all $K$ large enough

[TABLE]

By Sublemma 5.8 in [3], there exist two constants $C>0$ and $c>0$ such that for all $K$ large enough, and for all $t>0$

[TABLE]

Now, for $q\in\mathds{N}\backslash\{0\}$ define

[TABLE]

We will first estimate \sup_{\ushort{n}}{\mathds{P}}_{\!\ushort{n}}\big{(}\ushort{N}^{{\scriptscriptstyle K}}(t_{q})\in\mathcal{D}^{c},T_{\ushort{0}}>t_{q}\big{)}. Note that $\ushort{N}^{{\scriptscriptstyle K}}(t_{q})\in\mathcal{D}^{c}$ implies $T_{\mathcal{D}^{c}}\leq t_{q}$. We distinguish two cases according to whether $\ushort{n}\in\Delta$ or $\ushort{n}\in\Delta^{\!c}\backslash\{\ushort{0}\}$.

Let $\ushort{n}\in\Delta$. It follows from (2.10) that

[TABLE]

Now let $\ushort{n}\in\Delta^{\!c}\backslash\{\ushort{0}\}$. We have

[TABLE]

Using the strong Markov property at time $T_{\scaleto{\Delta}{5pt}}$ and (2.10) we obtain

[TABLE]

We bound the second term recursively in $q$.

[TABLE]

where we used the strong Markov property at time $t_{q-1}$ and (2.9). This implies

[TABLE]

Therefore

[TABLE]

Taking $q=\lfloor K\rfloor$ we conclude that there exists a constant $c^{\prime}>0$ such that for $K$ large enough

[TABLE]

This implies

[TABLE]

but by (2.6)

[TABLE]

and the result follows from (2.5). ∎

Corollary 2.5**.**

For each $q\in\mathds{N}$, there exists $C_{q}>0$ such that for all $K$ large enough

[TABLE]

Proof.

It follows at once from (2.6) (at time $1$) and Theorem 2.1 that

[TABLE]

for $K$ large enough. We have

[TABLE]

We use Hölder inequality to get

[TABLE]

The first result follows from (2.11) and Proposition 2.4. The second estimate follows from the first one, and the bound $\sup_{\ushort{n}\in\mathcal{D}}\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}\leq\mathcal{O}(1)K$. ∎

We now estimate centered moments.

Theorem 2.6**.**

For each $q\in{\mathds{Z}}_{+}$, there exists $C_{q}>0$ such that for all $K$ large enough

[TABLE]

Proof.

The proof consists in a recursion over $q$. The bound is trivial for $q=0$. For $q\in\mathds{N}$ define the function

[TABLE]

where

[TABLE]

Recall that $\scaleto{\ushort{e}^{(j)}}{11.5pt}$ is the vector with $1$ at the $j$th coordinate and [math] elsewhere. From the trivial identity

[TABLE]

it follows that

[TABLE]

Indeed, applying the trinomial expansion to (2.12), we obtain

[TABLE]

Observe that if $p_{1}\leq q-2,p_{1}+p_{2}+p_{3}=q$ and then $2p_{1}+p_{2}=p_{1}+q-p_{3}\leq 2q-2-p_{3}\leq 2q-2$, since $p_{3}\geq 0$. This implies that

[TABLE]

It follows that

[TABLE]

where

[TABLE]

To get this bound, we used the fact that

[TABLE]

Using (1.10) we get

[TABLE]

where

[TABLE]

Integrating the equation (2.13) with respect to $\nu\!_{\scaleto{K}{3.5pt}}$ and using (2.7), (2.14), (2.15) and Proposition 2.4, we obtain

[TABLE]

Observing that $\nu\!_{\scaleto{K}{3.5pt}}(f_{0})\leq 1$, it follows by recursion over $q$ that, for each integer $q$, there exists $C^{\prime}_{q}>0$ such that, for all $K$ large enough, $\nu\!_{\scaleto{K}{3.5pt}}(f_{q})\leq C^{\prime}_{q}K^{q}$. Finally we have

[TABLE]

since $\mathcal{D}\subset\mathcal{D}_{1}$. The result follows using the previous estimate and Corollary 2.5. ∎

The next result gives a more precise estimate for the average of $\ushort{n}$ (instead of an error of order $\sqrt{K}\,$).

Proposition 2.7**.**

We have

[TABLE]

where $\ushort{\mu}^{{\scriptscriptstyle K}}$ is defined in (1.1). Moreover, since $\|\ushort{n}^{\scaleto{*}{3pt}}-K\ushort{x}^{\scaleto{*}{3pt}}\|_{\scaleto{2}{4pt}}=\mathcal{O}(1)$, we have

[TABLE]

Proof.

Define the functions

[TABLE]

By Taylor expansion and the polynomial bounds on $\ushort{B}$ and $\ushort{D}$ we get

[TABLE]

for some positive integer $p$ independent of $K$. Using Cauchy-Schwarz inequality, identity (2.7), Corollary 2.5 and Proposition 2.4 we get

[TABLE]

From Proposition 2.4, Theorem 2.6 and (2.5) we get

[TABLE]

The result follows from the invertibility of the $d\times d$ matrix $(\partial_{m}B_{j}(\ushort{x}^{\scaleto{*}{3pt}})-\partial_{m}D_{j}(\ushort{x}^{\scaleto{*}{3pt}}))$ which follows from H.3. The other inequalities follow immediately. ∎

Corollary 2.8**.**

For all $K>0$, we have

[TABLE]

Proof.

Combine Proposition 2.7 and Theorem 2.6. ∎

We now show that $\Sigma^{{\scriptscriptstyle K}}$ is indeed of order $K$.

Proposition 2.9**.**

There exist two strictly positive constants $c_{\mathrm{(\ref{borninf})}}$ and $c_{\mathrm{(\ref{borninf})}}^{\prime}$ such that for all $K$ large enough, the matrix $\Sigma^{{\scriptscriptstyle K}}$ satisfies

[TABLE]

for the order among positive definite matrices, $\mathrm{Id}$ being the identity matrix, and, in particular,

[TABLE]

Proof.

We denote by $\widetilde{\Sigma}^{{\scriptscriptstyle K}}$ the positive definite matrix

[TABLE]

By (2.16) we have

[TABLE]

Let $\ushort{v}$ be a unit vector in d. We have

[TABLE]

From Lemma 5.3 in [3] there exists a constant $c>0$ such that for all $K$ large enough and all $\ushort{n}\in\Delta$,

[TABLE]

where $U_{\Delta}$ is the uniform distribution on $\Delta$. Therefore

[TABLE]

and we get

[TABLE]

The result follows. ∎

3 Controlling time averages of the estimators

For $T>0$, we define the time average of a function $f:{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}\to\real$ by

[TABLE]

The goal of this section is to obtain a control of $|S_{f}(T,K)-\nu\!_{\scaleto{K}{3.5pt}}(f)|$ for a suitable class of functions.

We recall the following result from [3, Theorem 3.1].

Theorem 3.1** ([3]).**

There exist $a>0$, $K_{0}>1$ such that, for all $t\geq 0$ and for all $K\geq K_{0}$, we have

[TABLE]

It is also proved in [3] that, for a time much larger than $\log K$ and much smaller than the extinction time (which is of order $\exp(\Theta(1)K)$), the law of the process at time $t$ is close to the qsd. The accuracy of the approximation depends on the initial condition. This suggests to study the distance between the law of the process at time $t$ and the qsd as a function of the initial condition, $K$ and $t$. This will result from (3.2) if {\mathds{P}}_{\!\ushort{n}}\big{(}T_{\ushort{0}}\leq t\big{)} can be estimated. In fact we prove a more general result.

Lemma 3.2**.**

For $\gamma\geq 0$, define \tau_{\gamma}=\inf\big{\{}t\geq 0:\|\ushort{N}^{{\scriptscriptstyle K}}(t)\|_{\scaleto{1}{4pt}}\leq\gamma K\big{\}}. There exist $\delta>0$, $\alpha>0$ and $C>0$ such that for all $\ushort{n}\in{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, $K\geq 1$, $0\leq\gamma\leq 1\wedge\frac{\alpha}{\scaleto{\|}{8pt}\hskip 0.5pt\ushort{x}^{\scaleto{*}{3pt}}\scaleto{\|}{8pt}_{\scaleto{1}{3pt}}}$ and $t\geq 0$, we have

[TABLE]

where

[TABLE]

Taking $\gamma=0$ in (3.3), we get

[TABLE]

Proof.

It follows from H.1 and H.3 (using Taylor’s expansion of $\ushort{X}(\ushort{x})$ near $\ushort{0}$) that there exists $\alpha_{0}\in(0,R)$ (where $R$ was introduced in Assumption H.3) such that for all $\ushort{x}\in\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$ satisfying $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}\leq\alpha_{\scaleto{0}{4pt}}$ we have

[TABLE]

For $\alpha\in(0,\alpha_{0}]$ and $\delta>0$ to be chosen later on, we define

[TABLE]

It is easy to verify that if $\langle\ushort{n},\ushort{x}^{\scaleto{*}{3pt}}\rangle>\alpha\,K+\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}^{\scaleto{*}{3pt}}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}$ we have

[TABLE]

If $\alpha K-\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}^{\scaleto{*}{3pt}}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}\leq\langle\ushort{n},\ushort{x}^{\scaleto{*}{3pt}}\rangle\leq\alpha K+\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}^{\scaleto{*}{3pt}}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}$ we have

[TABLE]

For $\langle\ushort{n},\ushort{x}^{\scaleto{*}{3pt}}\rangle\leq\alpha K-\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}^{\scaleto{*}{3pt}}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}$, we have $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}\leq\langle\ushort{n},\ushort{x}^{\scaleto{*}{3pt}}\rangle/\zeta\leq\alpha K/\zeta$, where $\zeta$ is defined in (3.4), and

[TABLE]

where the function $\mathfrak{g}$ is defined by

[TABLE]

We have

[TABLE]

From the differentiability of the vector fields $\ushort{B}$ and $\ushort{D}$ and using (3.6), it follows that there exists a constant $\Gamma>0$ such that, for all $0\leq s\leq 1$ and $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}<\alpha_{0}$ we have

[TABLE]

Therefore we can choose $\delta>0$ and $0<\alpha<\alpha_{0}$ such that

[TABLE]

Therefore, for all $\ushort{n}$, we have

[TABLE]

For $\tilde{\gamma}>0$ (independent of $K$), we define

[TABLE]

We apply Ito’s formula to $\psi$ to get

[TABLE]

We have

[TABLE]

hence

[TABLE]

Then

[TABLE]

Therefore

[TABLE]

To conclude, observe that

[TABLE]

for $\tilde{\gamma}=\gamma\,\scaleto{\|}{9pt}\hskip 0.5pt\ushort{x}^{\scaleto{*}{3pt}}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}$ because for all $\ushort{n}\in{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$,

[TABLE]

and $\|\ushort{N}^{{\scriptscriptstyle K}}(\tau_{\gamma})\|_{\scaleto{1}{4pt}}\leq\gamma K$. ∎

We have the following result.

Proposition 3.3**.**

For all bounded functions $h:{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}\to\real$, $t\geq 0$, $\ushort{n}\in{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, and $K>K_{0}$, we have

[TABLE]

where $\alpha,\delta$ and $\zeta$ are defined in Lemma 3.2, and $a$ and $K_{0}$ are defined in Theorem 3.1.

Proof.

From the bound (3.2) we get

[TABLE]

This implies

[TABLE]

using (3.5). ∎

We now extend Proposition 3.3 to more general functions. For $q\in{\mathds{Z}}_{+}$, we define the Banach space $\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,q}$ by

[TABLE]

We have the following result for time-averages of functions in $\mathscr{F}\!\!_{{\scriptscriptstyle K}}$.

Theorem 3.4**.**

For all $K>K_{0}$, $f\in\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,q}$, $T>0$, and $\ushort{n}\in{{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, we have

[TABLE]

where $\alpha,\delta$ and $\zeta$ are defined in Lemma 3.2, and $\lambda_{0}(K)$ is defined in (2.5).

Remark 3.1**.**

One can check that if one modifies slightly the definition of the time average (3.1) by integrating from $1$ to $T+1$, then one can remove the term $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}^{\scaleto{q}{4pt}}$ from the previous estimate.

Proof.

For $f\in\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,q}$, Corollary 2.5 gives

[TABLE]

By Proposition 2.3 we have

[TABLE]

Hence for $T\leq 1$ we get

[TABLE]

For $T>1$, we have by the Markov property that

[TABLE]

where we set

[TABLE]

By Corollary 2.2, the function $g$ is bounded and

[TABLE]

Applying Proposition 3.3 to $g$ thus gives

[TABLE]

Integrating over $s\in[0,T-1]$ yields

[TABLE]

Using Lemma 3.5 (stated and proved right after this proof), we finally obtain

[TABLE]

This finishes the proof of the theorem. ∎

We used the following lemma in the previous proof.

Lemma 3.5**.**

For $f\in\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,q}$ and $g$ defined in (3.8) we have

[TABLE]

Proof.

We write

[TABLE]

Since $\nu\!_{\scaleto{K}{3.5pt}}$ is a qsd, it follows by Cauchy-Schwarz inequality that

[TABLE]

where we used Corollaries 2.2 and 2.5 and the fact that under $\nu\!_{\scaleto{K}{3.5pt}}$ the law of $T_{\ushort{0}}$ is exponential with parameter $\lambda_{0}(K)$. The lemma is proved. ∎

4 Fluctuation and correlation relations

4.1 Proof of Theorem 1.1

Let

[TABLE]

For $1\leq i\leq d$, let $f_{i}(\ushort{n})=\langle\ushort{n}-\ushort{n}^{\scaleto{*}{3pt}},\scaleto{\ushort{e}^{(i)}}{11.5pt}\rangle$. We have, since $B_{i}(\ushort{x}^{\scaleto{*}{3pt}})=D_{i}(\ushort{x}^{\scaleto{*}{3pt}}),\;1\leq i\leq d$, and $\ushort{n}^{\scaleto{*}{3pt}}/K-\ushort{x}^{\scaleto{*}{3pt}}=\mathcal{O}(1)/K$,

[TABLE]

As in the previous proof, we split the integrals according to whether $\ushort{N}^{{\scriptscriptstyle K}}(t)\in\mathcal{D}$ or $\ushort{N}^{{\scriptscriptstyle K}}(t)\in\mathcal{D}^{c}$. Using Cauchy-Schwarz inequality, Corollary 2.5, and the fact that $\nu\!_{\scaleto{K}{3.5pt}}$ is a qsd, the second contribution is exponentially small in $K$. In the first contribution, we use Taylor expansion around $\ushort{x}^{\scaleto{*}{3pt}}$. The error terms are bounded by

[TABLE]

Now we use Cauchy-Schwarz inequality, Theorem 2.6 and that $\nu\!_{\scaleto{K}{3.5pt}}$ is a qsd to obtain

[TABLE]

Since $M^{*}$ has a spectrum contained in the open left half-plane by H.3, we integrate the equation

[TABLE]

from [math] to $\tau$ using the method of constant variation and obtain

[TABLE]

We arrive at the desired relation by using (2.16).

4.2 Proof of Theorem 1.2

Recall that

[TABLE]

We will first do the proof with the following matrix instead of $\Sigma^{{\scriptscriptstyle K}}$:

[TABLE]

On the one hand we have by (2.7)

[TABLE]

By Theorem 2.6 and (2.5) the right-hand side of this equation is exponentially small in $K$. On the other hand, using formula (2.1) we have

[TABLE]

We split each integral by separating integration over $\mathcal{D}$ (defined in (2.8)) and $\mathcal{D}^{c}$. Inside $\mathcal{D}^{c}$, we apply Corollary 2.5 and use the assumption that $B$ and $D$ are polynomially bounded. Inside $\mathcal{D}$, we use Taylor’s formula around $\ushort{x}^{\scaleto{*}{3pt}}$ for the functions $B_{i}(\ushort{n}/K)-D_{i}(\ushort{n}/K)$, and $B_{i}(\ushort{n}/K)+D_{i}(\ushort{n}/K)$. We also use that $B_{i}(\ushort{x}^{\scaleto{*}{3pt}})=D_{i}(\ushort{x}^{\scaleto{*}{3pt}}),\;1\leq i\leq d$, and \ushort{n}^{\scaleto{*}{3pt}}/K-\ushort{x}^{\scaleto{*}{3pt}}=\mathcal{O}\big{(}\frac{1}{K}\big{)}. The error terms are then bounded by

[TABLE]

respectively. Using Theorem 2.6, both bounds are of order $\sqrt{K}$. We obtain

[TABLE]

which can be written in the more compact form

[TABLE]

where $\mathcal{D}^{{\scriptscriptstyle K}}$ is the diagonal matrix of averages birth (or death) rates. To finish the proof, it remains to replace $\widetilde{\Sigma}^{{\scriptscriptstyle K}}$ by $\Sigma^{{\scriptscriptstyle K}}$. This is done by using (2.17).

Remark 4.1**.**

Note that each term on the left hand side is of order $K$, see Corollary 2.8.

Remark 4.2**.**

We will see in Appendix C that the qsd $\,\nu_{K}\,$ around $\ushort{n}^{\scaleto{*}{3pt}}$ is well approximated at scale $\sqrt{K}$ by a Gaussian distribution. Dividing out (4.1) by $2K$ and taking the limit $K\to\infty$, we recover Relation (C.1), as expected from Theorem C.1.

5 Variance estimates for the estimators

It is straightforward to apply Theorem 3.4 to $S^{{\scriptscriptstyle\ushort{\mu}}}(T,K)$, $S^{{\scriptscriptstyle C}}(T,\tau,K)$, $S^{{\scriptscriptstyle\mathcal{D}}}(T,K)$, and $S^{{\scriptscriptstyle\Sigma}}(T,K)$, which are defined respectively in (1.4), (1.5), (1.8), and (1.9). This gives the bound (1.6) on $S^{{\scriptscriptstyle\ushort{\mu}}}(T,K)$ anounced in Section 1. The bounds for the other estimators all have the same structure. We will not state them.

In this section we prove two variance estimates for any time average $S_{f}(T,K)$ with $f\in\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,q}$. In the first one, one starts from anywhere in ${{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, while in the second one the starting distribution is the qsd. Recall that $S^{{\scriptscriptstyle\Sigma}}(T,K)=S^{{\scriptscriptstyle C}}(T,0,K)$. We will only give the proofs of these estimates for $S^{{\scriptscriptstyle\Sigma}}(T,K)$, since manipulating $S^{{\scriptscriptstyle C}}(T,\tau,K)$ is cumbersome but otherwise the proofs are the same.

Proposition 5.1**.**

There exist strictly positive constants $\delta^{\prime},\zeta^{\prime},\alpha^{\prime},\theta^{\prime},C^{\prime}$ and $K_{0}\geq 2$ such that, for all $K\geq K_{0}$, $f\in\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,q}$ (see Definition 3.7), $T\geq 0$, and $\ushort{n}\neq\ushort{0}$, we have

[TABLE]

where $c_{q}$ was defined in Proposition 2.3.

One can use Chebyshev inequality to bound {\mathds{P}}_{\!\ushort{n}}\big{(}\big{|}S_{f}(T,K)-\nu\!_{\scaleto{K}{3.5pt}}(f)\big{|}>\delta\big{)}, for any $\delta>0$.

The proof of Proposition 5.1 is postponed to Appendix A. The previous estimate, as well as all the estimates we will give below, have the same behaviour in their dependence in $K$, $\ushort{n}$ and $T$. They display the qualitative behaviour that we met several times:

The bounds are not useful for $K$ too small. 2. 2.

If $K$ is large, the bounds are not useful if $\ushort{n}$ is small (order one) because the process can be absorbed at $\ushort{0}$ in a time of order one with a sizeable probability. 3. 3.

Finally, for $K$ large and $\ushort{n}$ of order $K$, the time $T$ must be large enough (polynomial in $K$ in our bounds) but not too large (less than an exponential in $K$ because the process can reach the origin with high probability in such large times).

Integrating the previous estimate with respect to the qsd, we get the following control.

Corollary 5.2**.**

There exist two positive constants $C^{\prime\prime}>0$ and $\theta^{\prime\prime}$ such that for all $K\geq K_{0}$, for all $f\in\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,q}$ and for all $T\geq 0$, we have

[TABLE]

where $K_{0}$ is as in the previous proposition, $c_{q}$ is defined in Proposition 2.3, and $C_{q}$ is defined in Corollary 2.5.

Observe that the previous inequality is only useful in the range $0\leq T\leq\operatorname{e}^{\theta^{\prime\prime}\!K}$. The proofs of the two previous estimates are postponed to Appendix A.

We now apply the previous results to our estimators.

Proposition 5.3**.**

For all $1\leq p\leq d$, we have

[TABLE]

and

[TABLE]

Proof.

The proof follows by applying Proposition 5.1 and Corollary 5.2 to the functions $f(\ushort{n})=n_{j}$, $1\leq j\leq d$, which belong to $\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,\scaleto{1}{4pt}}$. ∎

Proposition 5.4**.**

For $1\leq p,p^{\prime}\leq d$ and for all $\ushort{n}\neq\ushort{0}$, we have

[TABLE]

and

[TABLE]

Proof.

The proof follows by applying Proposition 5.1 and Corollary 5.2 to the functions $f(\ushort{n})=n_{p}n_{p^{\prime}}$, $1\leq p,p^{\prime}\leq d$, which belong to $\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,\scaleto{2}{4pt}}$. ∎

Proposition 5.5**.**

There exist positive constants $\widetilde{C},\tilde{\theta}$, $\tilde{\delta}$, $\tilde{\zeta}$ and $\tilde{\beta}$ such that for all $K\geq 2$, $T>0$ and $1\leq\ell\leq d$,

[TABLE]

where

[TABLE]

and $A_{\ell}>0$, $q_{\ell}\in\mathds{N}$, are such that, for all $\ushort{x}\in\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$,

[TABLE]

*The existence of $A_{\ell}$ and $q_{\ell}$ follows from the assumptions on $B$. The constants $C_{q_{\ell}}$ and $c_{q_{\ell}}$ are defined in Corollary 2.5 and Lemma 2.3, respectively.

We also have*

[TABLE]

[TABLE]

where

[TABLE]

Proof.

First observe that

[TABLE]

where $\mathcal{N}^{{\scriptscriptstyle K}}_{\ell}(0,T)$ is defined in Appendix B. By assumption, the function $f_{\ell}(\ushort{n})=K^{q_{\ell}}\,B_{\ell}\left(\frac{\ushort{n}}{K}\right)\in\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,q_{\ell}}$. Let $\mathfrak{m}$ be any probability measure on ${{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$ having all its moments finite. We apply Theorem 3.4 to the function $f_{\ell}$, and then using integration against $\mathfrak{m}$ we get

[TABLE]

[TABLE]

We now apply the identity in Proposition B.1 and divide by $K^{q_{\ell}-1}$. We obtain

[TABLE]

We now estimate

[TABLE]

The second integral is bounded from above by $\mathcal{O}(1)/K$ using the polynomial bound on $B_{\ell}$ and the first estimate in Corollary 2.5. For the first integral we use Taylor expansion around $\ushort{x}^{\scaleto{*}{3pt}}$ to first order, then Cauchy-Schwarz inequality, and finally Theorem 2.6 for $q=1$. Therefore we obtain

[TABLE]

Now we apply the estimate in Proposition B.1 to obtain

[TABLE]

For the first term we use either Corollary 5.2 or Proposition 5.1. For the second term we use (5.1). For the third and last term we apply Theorem 3.4, integrate with respect to $\mathfrak{m}$ and use (5.2). To finish the proof, we replace $\mathfrak{m}$ by either $\delta_{\ushort{n}}$ or $\nu\!_{\scaleto{K}{3.5pt}}$. ∎

Recall that $B_{p}(\ushort{x}^{\scaleto{*}{3pt}})=D_{p}(\ushort{x}^{\scaleto{*}{3pt}})$, $1\leq p\leq d$.

Proposition 5.6**.**

Under the assumptions of Proposition 5.1 and Corollary 5.2, we have, for all $1\leq p,p^{\prime}\leq d$, and $\tau\geq 0$,

[TABLE]

and

[TABLE]

Proof.

The proof requires some simple modifications of the proofs of Propostions 5.1 and 5.2. This is left to the reader. ∎

Remark 5.1**.**

If one modifies slightly the definition of the estimator by integrating from time $1$, then, in the four previous propositions, one can replace the factor $(\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}+K)$ by $K$, and the factor $(\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}^{\scaleto{2}{3.6pt}}+K^{\scaleto{2}{4pt}})$ by $K^{\scaleto{2}{4pt}}$.

Recall that we defined in Section 1 an empirical matrix $M^{*}_{\scriptscriptstyle{\mathrm{emp}}}(T,\tau,K)$ by

[TABLE]

and an empirical resilience by

[TABLE]

From the above results one can derive various statistical estimates for the difference between $\rho^{\scaleto{*}{3pt}}_{\scriptscriptstyle{\mathrm{emp}}}(T,\tau,K)$ and $\rho^{\scaleto{*}{2.8pt}}$. We have the following result which was stated at the end of Section 1.2. As already mentioned, we use the symbol $\ll$ which is not rigorously defined to formulate a more transparent bound. The reader can easily step through the proof to get a more precise, but rather cumbersome bound. Let us also note that the dependence on the initial state $\ushort{n}$ is related to the part “with a probability larger than $1-1/K$” of the statement. Indeed, the estimate of this probability results from Chebychev inequality and variance estimates in which the process is started in $\ushort{n}$.

Theorem 5.7**.**

For $\tau=\Theta(1)$, $\ushort{n}=\Theta(K)$ (initial state) and $0<T\ll\exp(\Theta(1)K)$, and $K$ large enough, we have

[TABLE]

with a probability higher than $1-1/K$. In particular, if $T\geq K^{\scaleto{5}{4.5pt}}$, we have

[TABLE]

Proof.

It follows from Propositions 5.4 and 5.6 and the standing assumptions that, with a probability higher that $1-1/K$, we have

[TABLE]

and

[TABLE]

($\|\cdot\|$ stands for any matrix norm on $\mathds{R}^{d\times d}$ since they are all equivalent.) We now use Theorem 1.1 and Proposition 2.9 to obtain

[TABLE]

The result follows since $\tau$ is of order one. ∎

Appendix A Proof of the two variance estimates

A.1 Starting from anywhere: proof of Proposition 5.1

It is enough to prove the result for $\|f\|_{{\scriptscriptstyle K}\!,q}=1$. We have

[TABLE]

Step 1 is to estimate the contribution of the range $0\leq t_{1}\leq t_{2}\leq 1$. Using Cauchy-Schwarz inequality and Proposition 2.3 we get

[TABLE]

Step 2 is to estimate the contribution in the range $0\leq t_{2}-1\leq t_{1}\leq t_{2}$. This implies that $T>1$. We have using again Proposition 2.3

[TABLE]

**Step 3

**(1) Using the Markov property and the definition of $g$ (see (3.8)) we have

[TABLE]

Let us first write

[TABLE]

as the sum of $J_{1}(\ushort{n})$ and $J_{2}(\ushort{n})$ where

[TABLE]

and

[TABLE]

We further decompose $J_{1}(\ushort{n})$ as $J_{1,1}(\ushort{n})+J_{1,2}(\ushort{n})$ where

[TABLE]

and

[TABLE]

Since $\ushort{0}$ is an absorbing state, we have for all $\ushort{n}\neq\ushort{0}$ that

[TABLE]

(2) We start by estimating $J_{2}(\ushort{n})$. Since $\ushort{0}$ is an absorbing state, we have

[TABLE]

Note that $g(\ushort{0})={\mathds{E}}_{\ushort{0}}[f(\ushort{N}^{{\scriptscriptstyle K}}(1))]=f(\ushort{0})$. Since we are going to use Lemma 3.2, we write $J_{2}(\ushort{n})=J_{2,1}(\ushort{n})+J_{2,2}(\ushort{n})$ where

[TABLE]

and

[TABLE]

We first estimate $J_{2,1}(\ushort{n})$. Using (3.9), Lemma 3.2 with $\gamma=0$, and since $f$ belongs to $\mathscr{F}\!\!_{{\scriptscriptstyle K}\!,q}$ (see (3.7)), we have

[TABLE]

where we used Proposition 2.3 for the second inequality.

We now estimate $J_{2,2}(\ushort{n})$ by splitting it as $J_{2,2,1}(\ushort{n})+J_{2,2,2}(\ushort{n})$ where

[TABLE]

and

[TABLE]

where

[TABLE]

Proceeding as before we get

[TABLE]

We used Lemma 3.2 with $\gamma=0$.

We now handle $J_{2,2,2}(\ushort{n})$.

Note that $\gamma\leq 1\wedge\frac{\alpha}{\|\ushort{x}^{\scaleto{*}{3pt}}\|_{\scaleto{1}{3pt}}}$. We proceed as before with $f$ and $g$, and we use Lemma 3.2 with

[TABLE]

to get

[TABLE]

(3) Let us now estimate $\left|J_{1,2}(\ushort{n})-\nu\!_{\scaleto{K}{3.5pt}}(f)^{\scaleto{2}{4pt}}\right|$ for all $\ushort{n}\neq\ushort{0}$. We have

[TABLE]

(3)-(i) By Theorem 3.1 and since $\ushort{N}^{{\scriptscriptstyle K}}(t_{1})\neq\ushort{0}$, we have

[TABLE]

Hence, using Proposition 2.3, we get for all $\ushort{n}\neq\ushort{0}$

[TABLE]

(3)-(ii) We have

[TABLE]

[TABLE]

Define $0<\gamma^{\prime}<\beta$ by

[TABLE]

We split the right hand side in two terms:

[TABLE]

The first term is estimated using the growth property of $f$, Lemma 3.2, and Cauchy-Schwarz inequality, namely

[TABLE]

To deal with the second term, we observe using Lemma 3.2 and Proposition 2.3 that, if $\|\ushort{N}^{{\scriptscriptstyle K}}(t_{1})\|_{\scaleto{1}{4pt}}>\gamma^{\prime}K$, then

[TABLE]

Now

[TABLE]

(3)-(iii) Let us now prove that for all $\ushort{n}\neq\ushort{0}$,

[TABLE]

For $0\leq t_{1}\leq 1$, using Proposition 2.3 we obtain

[TABLE]

We now deal with $t_{1}>1$. The Markov property gives

[TABLE]

where

[TABLE]

is a function bounded by $\mathcal{O}(1)K^{q}$. For $\ushort{n}\neq\ushort{0}$, we use Theorem 3.1 and Corollary 2.2 to get

[TABLE]

Since $\nu\!_{\scaleto{K}{3.5pt}}$ is the qsd, we have

[TABLE]

Using Corollary 2.5, Lemma 3.2 and the properties of $f$ we obtain

[TABLE]

and (A.1) is proved.

(3)-(iv) Let us note that

[TABLE]

Proposition 2.3 and Lemma 3.5 give

[TABLE]

(3)-(v) Collecting the informations given in the four previous estimates, we obtain a precise estimation of $\left|J_{1,2}(\ushort{n})-\nu\!_{\scaleto{K}{3.5pt}}(f)^{\scaleto{2}{4pt}}\right|$ for all $\ushort{n}\neq\ushort{0}$.

(3)-(vi) We have

[TABLE]

Collecting the above relevant estimates we obtain that there exist $\delta^{\prime},\zeta^{\prime},\beta^{\prime},\theta^{\prime}$ (all being positive and independent of $K$) such that

[TABLE]

Now we have

[TABLE]

The final result for $T\geq 1$ follows by collecting all estimates. For $T<1$ the bound follows directly from Proposition 2.3.

A.2 Starting from the qsd: proof of Corollary 5.2

The result follows from Proposition 5.1 by integrating over $\ushort{n}$ with respect to the qsd. More precisely, we have

[TABLE]

The integrals of the first and third terms with respect to the q.s.d are estimated using Corollary 2.5. We deal with second term:

[TABLE]

The third integral is estimated using the fact that the integrand is exponentially small in $K$. The second integral is estimated using the first estimate in Corollary 2.5. We finally deal with the first integral. If $\ushort{n}\in\mathcal{D}$ then $\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{1}{4.pt}}\geq\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}\geq\scaleto{\|}{9pt}\hskip 0.5pt\ushort{n}^{\scaleto{*}{3pt}}\scaleto{\|}{9pt}_{\scaleto{2}{3.5pt}}\scaleto{/2}{9pt}$. If $\{\|\ushort{n}\|_{\scaleto{1}{4pt}}<\beta^{\prime}K/\zeta^{\prime}\}\cap\mathcal{D}\neq\emptyset$, on this set we have \operatorname{e}^{-\delta^{\prime}\big{(}\zeta^{\prime}\frac{\|\ushort{n}\|_{\scaleto{1}{3pt}}}{K}\wedge\,\beta^{\prime}\big{)}K}\leq\operatorname{e}^{-\delta^{\prime}\zeta^{\prime}\frac{\|\ushort{n}^{\scaleto{*}{2.2pt}}\|_{\scaleto{2}{3pt}}}{2}} (exponentially small in $K$). The estimate follows.

Appendix B Counting the number of births

Denote by $\mathcal{N}^{{\scriptscriptstyle K}}_{\ell}(t_{1},t_{2})$ the number of births of species of type $\ell$ between the times $t_{1}$ and $t_{2}$ ($1\leq\ell\leq d$, $0\leq t_{1}\leq t_{2}$).

Proposition B.1**.**

For any probability measure $\mathfrak{m}$ on ${{\mathds{Z}}}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$, we have

[TABLE]

and

[TABLE]

Proof.

Recall that the generator of the process is given in (2.1). Let us now give a pathwise representation of the process. We introduce $d$ independent point Poisson measures $M_{\ell}(\!\hskip 1.0pt\mathrm{d}s,\!\hskip 1.0pt\mathrm{d}\theta)$ on $\mathds{R}_{\scaleto{+}{4.5pt}}^{\scaleto{2}{4.5pt}}$ with intensity $\hskip 1.0pt\mathrm{d}s\hskip 1.0pt\mathrm{d}\theta$. We define the $d$-dimensionnal càd-làg process $(N_{t},t\in\mathds{R}_{\scaleto{+}{4.5pt}})$

[TABLE]

Then the number of births of species of type $\ell$ occuring between the times $t_{1}$ and $t_{2}$ is given by

[TABLE]

Using the Markov property we get at once the first identity.

We now establish the estimate. Indeed

[TABLE]

By the $\mathds{L}^{2}$-isometry for jump processes (see [6, Formula (3.9) p.62]), we have

[TABLE]

This finishes the proof. ∎

Appendix C Gaussian limit for the rescaled qsd

We have the following theorem of independent interest. A part of this theorem partially generalizes a result obtained in [2] for models involving a single species ($d=1$). Recall that $\ushort{n}^{\scaleto{*}{3pt}}=\lfloor K\ushort{x}^{\scaleto{*}{3pt}}\rfloor$.

Theorem C.1**.**

For all $K>1$, define the measure $\mathfrak{a}_{K}$ on the Borel $\sigma$-algebra of d by

[TABLE]

Then $(\mathfrak{a}_{{\scriptscriptstyle K}})_{{\scriptscriptstyle K}}$ converges weakly to the centered Gaussian measure with covariance matrix

[TABLE]

where $\mathcal{B}^{*}$ is the diagonal matrix with entries $B_{\ell}(\ushort{x}^{\scaleto{*}{3pt}})=D_{\ell}(\ushort{x}^{\scaleto{*}{3pt}})$. The matrix $\mathcal{S}$ is also the unique symmetric solution of the (Lyapunov) equation (fluctuation-dissipation relation)

[TABLE]

Remark C.1**.**

We have

[TABLE]

This follows by dividing out equation (1.7) by $K$, letting $K$ tend to infinity, and using the uniqueness of the (symmetric) solution of (C.1).

Proof.

By Theorem 2.6, the family of measures $(\mathfrak{a}_{{\scriptscriptstyle K}})_{{\scriptscriptstyle K}}$ is tight. For $\ushort{p}\in{}^{d}$ define

[TABLE]

It follows also from Theorem 2.6 that the family of functions $(H_{{\scriptscriptstyle K}})$ is uniformly bounded in $C^{2}$. We will prove that

[TABLE]

This will entail that there is only one weak accumulation point for $(\mathfrak{a}_{{\scriptscriptstyle K}})_{{\scriptscriptstyle K}}$. The proof will be the consequence of Prokhorov Theorem [1]. Using (2.7) and (2.5), we have

[TABLE]

We also have

[TABLE]

Using Taylor expansion, and the moments estimates and the polynomial bounds on $B_{\ell}$ and $D_{\ell}$ (and $B_{\ell}(\ushort{x}^{\scaleto{*}{3pt}})=D_{\ell}(\ushort{x}^{\scaleto{*}{3pt}})$) we obtain

[TABLE]

We conclude that every accumulation point $\breve{H}$ of $(H_{{\scriptscriptstyle K}})_{{\scriptscriptstyle K}}$ is bounded in $C^{1}$, satisfies $\breve{H}(\ushort{0})=1$, and is a solution of the equation

[TABLE]

Then (C.2) follows from Lemma C.2 (stated and proved right after this proof) with $A=M^{*}$. ∎

Lemma C.2**.**

Let $(B_{j})$ be $d$ strictly positive numbers and $A$ a real $d\times d$ matrix such that $\mathrm{Sp}(A)\subset\{z\in\mathds{C}:\operatorname{Re}(z)<0\}$. Then there exists a unique $C^{1}(\mathds{R}^{d},\mathds{R})$ function $H$ satisfying $H(\ushort{0})=1$ and

[TABLE]

This function is given by

[TABLE]

where

[TABLE]

where $\mathcal{B}$ is the diagonal matrix with entries $(B_{j})$. The matrix $S$ is also the unique symmetric solution of the equation

[TABLE]

Proof.

We use the method of characteristics. For all $\ushort{p}\in\mathds{R}^{d}$, we define the function $\ushort{p}(s),s\geq 0$ as the solution of

[TABLE]

Let

[TABLE]

Let $H$ be a solution of (C.3). It is easy to check that for all $\ushort{p}\in\mathds{R}_{\scaleto{+}{5pt}}^{\scaleto{d}{5pt}}$ and $s\in\real$

[TABLE]

Integrating from [math] to $u$ yields

[TABLE]

From the spectral properties of $A$ we get

[TABLE]

Therefore

[TABLE]

and

[TABLE]

Finally we get from the spectral properties of $A$

[TABLE]

This finishes the proof of the lemma. ∎

Acknowledgements: We thank the two anonymous referees for fruitful comments and suggestions.

Funding: S. M. has been supported by the Chair “Modélisation Mathématique et Biodiversité” of Veolia Environnement-Ecole Polytechnique-Museum national d’Histoire naturelle-Fondation X. P. C. and S. M. warmly thank the Basal Conicyt CMM AFD170001 project. J.-R. C. and P. C. also acknowledge the hospitality of the Instituto de Física de San Luis Potosí.