Equilibrium and surviving species in a large Lotka–Volterra system of differential equations

Abstract

Lotka–Volterra (LV) equations play a key role in the mathematical modeling of various ecological, biological and chemical systems. When the number of species (or, depending on the viewpoint, chemical components) becomes large, basic but fundamental questions such as computing the number of surviving species still lack theoretical answers. In this paper, we consider a large system of LV equations where the interactions between the various species are a realization of a random matrix. We provide conditions to have a unique equilibrium and present a heuristics to compute the number of surviving species. This heuristics combines arguments from Random Matrix Theory, mathematical optimization (LCP), and standard extreme value theory. Numerical simulations, together with an empirical study where the strength of interactions evolves with time, illustrate the accuracy and scope of the results.

Appendices

A Simulation details

Simulations were performed in Python. All the figures and the code are available in Clenet (2022).

Simulations on the properties of surviving species are performed in two different ways. The theoretical solutions are obtained resolving numerically the system of equations of heuristics 1. We use a solver (cf. scipy.optimize) to find a local minimum of the function defined by the system of equations (a modification of the Powell hybrid method). The empirical solutions are computed using a Monte Carlo experiment. We simulate a large number of matrix matrix B, we resolve the associated LCP problem using the Lemke’s algorithm. Then, we use the LCP solution to calculate the properties of the surviving species: proportion of survivors, etc. Finally, we make an average on the ensemble of experiments. The Lemke algorithm is implemented in the lemkelcp package and can be found on Lamperski (2019). The dynamics of the Lotka–Volterra are achieved by a Runge–Kutta of order 4 (RK4) implemented in the code.

B Proof of Theorem 2

We have

$$\begin{aligned} I- B +I -B^T = 2I - (B+B^T) = 2I - \left( \frac{A+A^T}{\alpha \sqrt{n}}+\frac{2\mu }{n} \varvec{1}\varvec{1}^*\right) . \end{aligned}$$

Notice that $2I- (B+B^T)$ is positive definite iff the top eigenvalue of $B+B^T$ is lower than 2:

$$\begin{aligned} \lambda _{\max }(B+B^T) < 2 \,. \end{aligned}$$

(23)

We first focus on the random part $(A+A^T)/\alpha$ which is a symmetric matrix with independent ${{\mathcal {N}}}(0,2/\alpha ^2)$ entries above the diagonal (note that the distribution of the diagonal entries is different from the off-diagonal entries, with no asymptotic effect). In this case, it is well known that the largest eigenvalue of the normalized matrix (or equivalently its spectral norm since the matrix is symmetric) a.s. converges to the right edge of the support of the semi-circle law (see (Bai and Silverstein 2010, Th. 5.2)):

$$\begin{aligned} \lambda _{\max }\left( \frac{A+A^T}{\alpha \sqrt{n}} \right) \xrightarrow [n\rightarrow \infty ]{a.s.} \frac{2 \sqrt{2}}{\alpha }\,. \end{aligned}$$

(24)

In the centered case ($\mu = 0$), condition (23) occurs if $\alpha > \sqrt{2}$.

We now consider the general case where $\mu \ne 0$. Notice that the rank-one perturbation matrix $P = \frac{2\mu }{n} \varvec{1}\varvec{1}^*$ admits a unique non zero eigenvalue $2\mu$. Denote by $\check{A}=\frac{A+A^T}{\alpha \sqrt{n}}$. We are interested in the top eigenvalue of the symmetric matrix $\check{A}+P$. Based on a result by Capitaine et al. (Capitaine et al. 2009, Th. 2.1), we have:

$$\begin{aligned} \lambda _{\max }(\check{A}+P) \xrightarrow [n \rightarrow \infty ]{a.s}\left\{ \begin{array}{ll} 2\mu +\frac{1}{\alpha ^2 \mu } &{}\text {if } \mu > \frac{1}{\sqrt{2} \alpha }\,, \\ \frac{2 \sqrt{2}}{\alpha } &{}\text {else.} \end{array}\right. \end{aligned}$$

This result is illustrated in Fig. 11.

Fig. 11

Spectrum (histogram) of the Hermitian random matrix $B+B^T$ ($n=1000$, $\alpha =\sqrt{2}$). The solid line represents the semi-circular law. In Fig. 11a, $\mu =0$ and there is no oulier. In Fig. 11b, $\mu =1.5$ and one can notice the presence of an eigenvalue outside the bulk of the semi-circular law. The dashed line indicates its theoretical value

Assume first that $\mu \le \frac{1}{\alpha \sqrt{2}}$ (corresponding to zone ${{\mathcal {C}}}$ in Fig. 2), then $\lambda _{\max }( \check{A}+P) \xrightarrow [n\rightarrow \infty ]{a.s.} \frac{2\sqrt{2}}{\alpha }$, which is strictly lower than 2 (cf. condition (23)) if $\alpha >\sqrt{2}$. Hence $\lambda _{\max }(\check{A} +P)$ is eventually strictly lower than 2 under this condition.

Assume now that $\mu >\frac{1}{\alpha \sqrt{2}}$ (corresponding to zone ${{\mathcal {B}}}$ in Fig. 2), then

$$\begin{aligned} \lambda _{\max }( \check{A}+P) \xrightarrow [n\rightarrow \infty ]{a.s.} 2\mu +\frac{1}{\alpha ^2 \mu }\,. \end{aligned}$$

We are interested in the conditions for which $2\mu +\frac{1}{\alpha ^2 \mu }<2$ or equivalently

$$\begin{aligned} 2\alpha ^2 \mu ^2 - 2\alpha ^2 \mu +1<0\,. \end{aligned}$$

(25)

An elementary study of the polynomial $\xi (X)= 2\alpha ^2 X^2 - 2\alpha ^2 X +1$ yields that $\xi$’s discriminant is positive if $\alpha >\sqrt{2}$,

$$\begin{aligned} \xi (\mu ^{\pm })=0\quad \Leftrightarrow \quad \mu ^{\pm }= \frac{1}{2} \pm \frac{1}{2} \sqrt{1- \frac{2}{\alpha ^2}}\, \end{aligned}$$

and $\xi \left( \frac{1}{\alpha \sqrt{2}}\right) <0$, so that $\frac{1}{\alpha \sqrt{2}}\in (\mu ^-, \mu ^+)$. In particular condition (25) is fulfilled if

$$\begin{aligned} \mu \in \left( \frac{1}{\alpha \sqrt{2}}\,\ \frac{1}{2}+\frac{1}{2} \sqrt{1-\frac{2}{\alpha ^2}}\right) . \end{aligned}$$

Under this condition, (25) is fulfilled and a.s. $\limsup _{n\rightarrow \infty }\lambda _{\max }(\check{A}+P) < 2$, which completes the proof: we can then rely on Theorem 1 to conclude.

C Construction of the heuristics

We first discuss Heuristics 1 and establish Equations (12), (13) and (14).

1.1 C.1 Equation (12)

We first recall a result on order statistics of a Gaussian sample. Consider a family $(Z_k)_{k\in [n]}$ of i.i.d. random variables ${{\mathcal {N}}}(0,1)$ and the associated order statistics

$$\begin{aligned} Z_1^*\ \le \ Z_2^*\ \le \cdots \le \ Z_n^*\,. \end{aligned}$$

Consider an index $\lfloor n\alpha \rfloor \in [n]$ where $\alpha \in (0,1)$ is fixed, then the typical location of $Z^*_{ \lfloor n\alpha \rfloor }$ is $\Phi ^{-1}(\alpha )$:

$$\begin{aligned} Z^*_{ \lfloor n\alpha \rfloor }\simeq \Phi ^{-1}(\alpha )\quad \text {as}\quad n\rightarrow \infty \,, \end{aligned}$$

(26)

see for instance (Smirnov 1949; Balkema and De Haan 1978).

Let $\varvec{x}^*$ be the equilibrium of (1) and consider the random variable

$$\begin{aligned} \check{Z}_k= \sum _{i\in {{\mathcal {S}}}} B_{ki} x_i^* = (B\varvec{x}^*)_k. \end{aligned}$$

We assume that asymptotically the $x_i^*$’s are independent from the $B_{ki}$’s, an assumption supported by the chaos hypothesis, see for instance Geman and Hwang (1982). Denote by $\mathbb {E}_{\varvec{x}^*}=\mathbb {E}(\,\cdot \mid \varvec{x}^*)$ the conditional expectation with respect to $\varvec{x}^*$. Notice that conditionally to $\varvec{x}^*$, the $\check{Z}_k$’s are independent Gaussian random variables, whose two first moments can easily be computed, see Sect. 1 below for the details:

$$\begin{aligned} \mathbb {E}_{\varvec{x}^*} \check{Z}_k= \mu \,\hat{p}\, \hat{m}\quad \text {and}\quad \text {var}_{\varvec{x}^*}(\check{Z}_k)= \frac{\hat{p}\hat{\sigma }^2}{\alpha ^2} \,. \end{aligned}$$

Notice that the fact that $\mathbb {E}_{\varvec{x}^*}$ and $\text {var}_{\varvec{x}^*}(\check{Z}_k)$ only depend on $\hat{p}, \hat{\sigma }$ and $\hat{m}$ which are (supposedly) converging quantities supports the idea that $\check{Z}_k$ is unconditionally a Gaussian random variable with moments:

$$\begin{aligned} \mathbb {E} \check{Z}_k= \mu \,p^*\, m^*\quad \text {and}\quad \text {var}(\check{Z}_k)= \frac{p^* (\sigma ^*)^2}{\alpha ^2} \,, \end{aligned}$$

where $p^*,m^*,\sigma ^*$ are resp. the limits of $\hat{p}, \hat{m}, \hat{\sigma }$. We now introduce the standard Gaussian random variables $(Z_k)_{k\in [n]}$ where

$$\begin{aligned} Z_k= \frac{\check{Z}_k - \mathbb {E} \check{Z}_k}{\sqrt{\text {var}(\check{Z}_k)}} = \alpha \frac{\check{Z}_k - \mu \,p^*\, m^*}{\sigma ^* \sqrt{p^*}}\,. \end{aligned}$$

Consider the equilibrium $\varvec{x}^*=(x_k^*)_{k\in [n]}$. If $k\in {{\mathcal {S}}}$, that is $x_k^*>0$, we have

$$\begin{aligned} 1- x_k^* + (B\varvec{x}^*)_k=0 \quad \Rightarrow \quad 1+ (B\varvec{x}^*)_k>0 \,. \end{aligned}$$

This identity has two implications:

$$\begin{aligned} x_k^*= 1 + (B\varvec{x}^*)_k\quad \text {and}\quad 1+ (B\varvec{x}^*)_k>0 \qquad \text {if}\ k\in {{\mathcal {S}}}\,. \end{aligned}$$

Relying on the representation $(B\varvec{x}^*)_k=\check{Z}_k$, we obtain the representation

$$\begin{aligned} x_k^* = 1+(B\varvec{x}^*)_k \ =\ 1+\mu \,p^* \,m^*+\frac{\sigma ^*\sqrt{p^*}}{\alpha } Z_k \qquad \text {if}\quad k\in {{\mathcal {S}}}. \end{aligned}$$

(27)

and the condition:

$$\begin{aligned} 1+(B\varvec{x}^*)_k \ =\ 1+\mu \,p^* \,m^*+\frac{\sigma ^*\sqrt{p^*}}{\alpha } Z_k\, >\, 0. \end{aligned}$$

If $k\notin {{\mathcal {S}}}$ then

$$\begin{aligned} 1+ (B\varvec{x}^*)_k\ =\ 1+\mu \,p^*\,m^* +\frac{\sigma ^*\sqrt{p^*}}{\alpha } Z_k\ \le \ 0 \end{aligned}$$

by the non invadability condition. Otherwise stated,

$$\begin{aligned} \left\{ \begin{array}{lcll} Z_k&{}\le &{} -\frac{\alpha (1+\mu \, p^* m^*)}{\sigma ^*\sqrt{p^*}} &{}\text {if}\ k\notin {{\mathcal {S}}},\\ Z_k&{}>&{} -\frac{\alpha (1+\mu \, p^* m^*)}{\sigma ^*\sqrt{p^*}} &{}\text {if}\ k\in {{\mathcal {S}}}.\\ \end{array} \right. \end{aligned}$$

Considering the order statistics of the $Z_k$’s we obtain:

$$\begin{aligned} Z_1^*\le \cdots \le Z_{\varvec{i}}^*\le -\frac{\alpha (1+\mu \, p^* m^*)}{\sigma ^*\sqrt{p^*}} \le Z_{\varvec{i}+1}^*\le \cdots \le Z_n^*\,. \end{aligned}$$

Now, there are exactly $n-|{{\mathcal {S}}}|=n(1-\hat{p})$ indices before the threshold corresponding to the components of $\varvec{x}^*$ equal to zero. In particular, index $\varvec{i}=n(1-\hat{p})$ corresponds to the value

$$\begin{aligned} Z^*_{\varvec{i}}\simeq -\frac{\alpha (1+\mu \,p^*\,m^*)}{\sigma ^*\sqrt{p^*}}. \end{aligned}$$

Relying on (26), we finally obtain

$$\begin{aligned} \Phi ^{-1}(1-\hat{p}) = -\frac{\alpha (1+\mu \,p^* \, m^*)}{\sigma ^*\sqrt{p^*}}. \end{aligned}$$

It remains to replace $\hat{p}$ by its limit $p^*$ to obtain (12).

1.2 C.2 Details on Eq. (12): moments of $\check{Z}_k$

We compute hereafter the conditional mean and variance of $\check{Z}_k=(B\varvec{x}^*)_k$ with respect to $\varvec{x}^*$. We rely on the following identities:

$$\begin{aligned} \mathbb {E} B_{ki} =\frac{\mu }{n},\quad \mathbb {E} (B_{ki})^2 = \frac{1}{\alpha ^2 n } + \frac{\mu ^2}{n^2} \simeq \frac{1}{\alpha ^2 n},\quad \mathbb {E} B_{ki}B_{kj} = \frac{\mu ^2}{n^2}\quad (i\ne j)\,. \end{aligned}$$

We first compute the conditional mean:

$$\begin{aligned} \mathbb {E}_{\varvec{x}^*}(\check{Z}_k)= & {} \sum _{i\in [n]} \mathbb {E}(B_{ki}) x_i^* = \sum _{i\in {{\mathcal {S}}}} \mathbb {E}(B_{ki}) x_i^* = \frac{\mu }{n}\sum _{i\in {{\mathcal {S}}}} x_i^*, = \mu \frac{|{{\mathcal {S}}}|}{n} \frac{1}{|{{\mathcal {S}}}|}\sum _{i\in {{\mathcal {S}}}} x_i^*, \\= & {} \mu \, \hat{p}\, \hat{m}\,. \end{aligned}$$

We now compute the second moment:

where the approximation in (a) follows from the fact that

$$\begin{aligned} \frac{1}{|{{\mathcal {S}}}|^2} \sum _{i,j\in {{\mathcal {S}}}}x_i^*x_j^* = \frac{1}{|{{\mathcal {S}}}|^2} \sum _{i \ne j}x_i^*x_j^* + {{\mathcal {O}}} \left( \frac{1}{|{{\mathcal {S}}}|} \right) \,. \end{aligned}$$

We can now compute the variance:

$$\begin{aligned} \text {var}_{\varvec{x}^*} \left( \check{Z}_k\right) = \mathbb {E}_{\varvec{x}^*}\left( \check{Z}_k^2\right) - \left( \mathbb {E}_{\varvec{x}^*}\check{Z}_k\right) ^2 \ =\ \frac{\hat{p}\, \hat{\sigma }^2}{\alpha ^2}\,. \end{aligned}$$

1.3 C.3 Equation (13)

Our starting point is the following generic representation of an abundance at equilibrium (either of a surviving or vanishing species):

$$\begin{aligned}&x_k^* = \left( 1+\mu \,p^* m^* +\frac{\sigma ^*\sqrt{p^*}}{\alpha } Z_k \right) \varvec{1}_{\{Z_k> - \delta ^*\}},\\&\quad = \left( 1+\mu \,p^* m^* \right) \varvec{1}_{\{Z_k> - \delta ^*\}} +\left( \frac{\sigma ^*\sqrt{p^*}}{\alpha } Z_k \right) \varvec{1}_{\{Z_k> - \delta ^*\}}. \end{aligned}$$

Summing over ${{\mathcal {S}}}$ and normalizing,

$$\begin{aligned} \begin{aligned} \frac{1}{|{{\mathcal {S}}}|}\sum _{k \in \mathcal {S}}x_k^{*}&=(1+\mu \, p^* m^*)\frac{1}{|{{\mathcal {S}}}|}\sum _{k \in \mathcal {S}}\varvec{1}_{\{Z_k> -\delta ^*\}}+\frac{\sigma ^*\sqrt{p^*}}{\alpha }\frac{1}{|{{\mathcal {S}}}|}\sum _{k \in \mathcal {S}}Z_k\varvec{1}_{\{Z_k> -\delta ^*\}}, \\ \hat{m}&{\mathop {=}\limits ^{(a)}} (1+\mu \, p^* m^*)+\frac{\sigma ^*\sqrt{p^*}}{\alpha }\frac{n}{|{{\mathcal {S}}}|}\frac{1}{n}\sum _{k \in [n]}Z_k\varvec{1}_{\{Z_k> -\delta ^*\}}, \\ \hat{m}&{\mathop {\simeq }\limits ^{(b)}}(1+\mu \,p^* m^*)+\frac{\sigma ^*\sqrt{p^*}}{\alpha } \frac{1}{\mathbb {P}(Z> -\delta ^*)}\mathbb {E}(Z\varvec{1}_{\{Z> -\delta ^*\}}), \\ \hat{m}&\simeq (1+\mu \, p^* m^*)+\frac{\sigma ^*\sqrt{p^*}}{\alpha }\mathbb {E}(Z \mid Z > -\delta ^*). \end{aligned} \end{aligned}$$

where (a) follows from the fact that $|{{\mathcal {S}}}| = \sum _{k \in \mathcal {S}}\varvec{1}_{\{Z_k > -\delta ^*\}}$ (by definition of ${{\mathcal {S}}}$), (b) from the law of large numbers $\frac{1}{n} \sum _{k\in [n]} Z_k \varvec{1}_{\{Z_k>-\delta \}} \xrightarrow [n\rightarrow \infty ]{} \mathbb {E}Z \varvec{1}_{\{Z>-\delta \}}$ and $\frac{|{{\mathcal {S}}}|}{n} \xrightarrow [n\rightarrow \infty ]{} \mathbb {P}(Z>-\delta ^*)$ with $Z\sim {{\mathcal {N}}}(0,1)$. It remains to replace $\hat{m}$ by its limit $m^*$ to obtain (13).

Equation (14) can be obtained similarly.

1.4 C.4 Equation (14)

As for the proof of (13), we start from the generic representation of $x_k^*$:

$$\begin{aligned} x_k^*\ {}&= \left( 1+\mu \, p^* m^* +\frac{\sigma ^*\sqrt{p^*}}{\alpha } Z_k\right) \varvec{1}_{\{Z_k> -\delta ^*\}},\\&= \left( 1+\mu \,p^* m^*\right) \varvec{1}_{\{Z_k> -\delta ^*\}}+\frac{\sigma ^*\sqrt{p^*}}{\alpha }Z_k\varvec{1}_{\{Z_k > -\delta ^*\}}\, . \end{aligned}$$

Taking the square, we get

$$\begin{aligned} x_k^{*2}&= \left( 1+\mu \,p^* m^*\right) ^2\varvec{1}_{\{Z_k> -\delta ^*\}}\\&\quad +2(1+\mu \,p^* m^*)\frac{\sigma ^*\sqrt{p^*}}{\alpha }Z_k\varvec{1}_{\{Z_k> -\delta ^*\}} +\frac{(\sigma ^*)^2p^*}{\alpha ^2}Z^2_k\varvec{1}_{\{Z_k > -\delta ^*\}}\,. \end{aligned}$$

Summing over ${{\mathcal {S}}}$ and normalizing, we get

$$\begin{aligned}&\frac{1}{|{{\mathcal {S}}}|}\sum _{k \in \mathcal {S}}(x_k^*)^{2} = (1+\mu \, p^* m^*)^2\frac{1}{|{{\mathcal {S}}}|}\sum _{k \in \mathcal {S}}\varvec{1}_{\{Z_k> -\delta ^*\}} \\&\quad +2(1+\mu \,p^* m^*)\frac{\sigma ^*\sqrt{p^*}}{\alpha } \frac{1}{|{{\mathcal {S}}}|}\sum _{k \in \mathcal {S}}Z_k\varvec{1}_{\{Z_k> -\delta ^*\}}\\&\quad +\frac{(\sigma ^*)^2 p^*}{\alpha ^2}\frac{1}{|{{\mathcal {S}}}|}\sum _{k \in \mathcal {S}}Z^2_k\varvec{1}_{\{Z_k > -\delta ^*\}}\,. \end{aligned}$$

Finally, we conclude by replacing the empirical means by their limits

$$\begin{aligned} \frac{1}{|{{\mathcal {S}}}|}\sum _{k \in \mathcal {S}}Z^i_k\varvec{1}_{\{Z_k> -\delta ^*\}}= & {} \mathbb {E}(Z^i\mid Z>-\delta ^*)\,\quad i=1,2\,. \end{aligned}$$

and get

$$\begin{aligned} \hat{\sigma }^2&= (1+\mu \,p^* m^*)^2+2(1+\mu \,p^* m^*)\frac{\sigma ^*\sqrt{p^*}}{\alpha }\mathbb {E}(Z \mid Z> -\delta ^*)\\&\quad +\frac{(\sigma ^*)^2p^*}{\alpha ^2}\mathbb {E}(Z^2 \mid Z > -\delta ^*)\,. \end{aligned}$$

It remains to replace $\hat{\sigma }$ by its limit $\sigma ^*$ to obtain (14).

D Density of the distribution of the persistent species

Assume that $x^*>0$, and let $f=\mathbb {R}\rightarrow \mathbb {R}$ be a bounded continuous test function. We have

$$\begin{aligned} \mathbb {E}f(x_k^*)= & {} \textbf{E}\left[ f\left( 1+ \frac{\sigma ^*\sqrt{p^*}}{\alpha } Z_k+\mu \,p^* m^*\right) \ \bigg |\ Z_k> -\delta ^* \right] \,, \\= & {} \int _{-\infty }^{\infty }f\left( 1+\mu \,p^* m^*+\frac{\sigma ^*\sqrt{p^*}}{\alpha }u\right) \frac{\varvec{1}_{\{u>-\delta ^*\}}}{1-\Phi (-\delta ^*)}\frac{e^{-\frac{u^2}{2}}}{\sqrt{2\pi }}du\,, \\= & {} \int _{0}^{\infty }f(y) e^{-\frac{1}{2}\left( \frac{\alpha }{\sigma ^* \sqrt{p^*}}y-\delta ^*\right) ^2} \frac{\alpha }{\sqrt{2\pi } \Phi (\delta ^*)\,p^*\, \sigma ^*}\,dy\,, \end{aligned}$$

hence the density of $x^*_k$.

1.1 D.1 Theoretical estimation of the diversity index

Recall that $|\mathcal {S}| = n \hat{p}$ is the number of surviving species and that

$$\begin{aligned} p_i = \frac{x_i}{\sum _{j\in {{\mathcal {S}}}} x_j} \end{aligned}$$

is the frequency of (surviving) species i.

To find a theoretical estimate of Hill number of order 1, we proceed by expansion and set

$$\begin{aligned} p_i = \frac{1}{|\mathcal {S}|}+\delta _i,\quad |\delta _i|\ll \frac{1}{|S|} \end{aligned}$$

where $\delta _i$ represents the deviation of species i from the standard frequency if all surviving species have the same frequency. Notice that $\sum _{i\in \mathcal {S}}\delta _i=0$.

$$\begin{aligned} H' \ =\ - \sum _{i\in \mathcal {S}} p_i \log (p_i) = - \sum _{i\in \mathcal {S}} \left( \frac{1}{|\mathcal {S}|}+ \delta _i \right) \log \left( \frac{1}{|\mathcal {S}|}+\delta _i \right) . \end{aligned}$$

We use the Taylor-Young formula of order 2 to decompose the log:

$$\begin{aligned} \log \left( \frac{1}{|\mathcal {S}|}+\delta _i \right)= & {} \log \left( \frac{1}{|\mathcal {S}|}\right) +|\mathcal {S}|\delta _i-\frac{\delta _i^2|\mathcal {S}|^2}{2}+\delta _i^3\varepsilon (\delta _i) \,, \\\approx & {} \log \left( \frac{1}{|\mathcal {S}|}\right) +|\mathcal {S}|\delta _i-\frac{\delta _i^2|\mathcal {S}|^2}{2} .\\ H'\approx & {} - \sum _{i\in \mathcal {S}} \left( \frac{1}{|\mathcal {S}|}+\delta _i \right) \left( \log \left( \frac{1}{|\mathcal {S}|}\right) +|\mathcal {S}|\delta _i-\frac{\delta _i^2|\mathcal {S}|^2}{2}\right) \,, \\= & {} - \sum _{i\in \mathcal {S}}\left[ \frac{1}{|\mathcal {S}|}\log \left( \frac{1}{|\mathcal {S}|}\right) +\delta _i-\frac{\delta _i^2|\mathcal {S}|}{2}+\delta _i \log \left( \frac{1}{|\mathcal {S}|} \right) + |\mathcal {S}|\delta _i^2-\frac{\delta _i^3 |\mathcal {S}|^2}{2} \right] \,, \\= & {} \log (|\mathcal {S}|) - \sum _{i\in \mathcal {S}} \frac{\delta _i^2|\mathcal {S}|}{2} +\sum _{i\in \mathcal {S}} \frac{\delta _i^3|\mathcal {S}|^2}{2}. \end{aligned}$$

Notice that $\sum _{i=1}^{|\mathcal {S}|} \frac{\delta _i^3|\mathcal {S}|^2}{2}$ is negligible since $|\delta _i|\ll |S|^{-1}$. The term 1 corresponds to the maximum value that the Shannon diversity index can take if $|\mathcal {S}|$ are present in the system. It remains to develop the second term of the r.h.s.

$$\begin{aligned} -\frac{1}{2}\sum _{i\in \mathcal {S}}\delta _i^2|\mathcal {S}|= & {} -\frac{|\mathcal {S}|}{2}\sum _{i\in \mathcal {S}} \left( \frac{x_i}{\sum _{j\in \mathcal {S}} x_j}-\frac{1}{|\mathcal {S}|} \right) ^2\, \\= & {} -\frac{|\mathcal {S}|}{2} \sum _{i\in \mathcal {S}} \left( \frac{x_i^2}{(\sum _{j\in \mathcal {S}}x_j)^2}-\frac{2}{|\mathcal {S}|}\frac{x_i}{\sum _{j\in \mathcal {S}} x_j}+\frac{1}{|\mathcal {S}|^2} \right) \, \\= & {} -\frac{|\mathcal {S}|}{2} \sum _{i\in \mathcal {S}} \left( \frac{x_i^2}{(\sum _{j\in \mathcal {S}}x_j)^2}\right) +\frac{1}{2}\, \\= & {} -\frac{|\mathcal {S}|}{2}\frac{\sum _{i\in \mathcal {S}}x_i^2}{|\mathcal {S}|^2 (\frac{1}{|\mathcal {S}|}\sum _{j\in \mathcal {S}}x_j)^2}+\frac{1}{2}\,\\= & {} -\frac{1}{2}\frac{\frac{1}{|\mathcal {S}|}\sum _{i\in \mathcal {S}}x_i^2}{ (\frac{1}{|\mathcal {S}|}\sum _{j\in \mathcal {S}}x_j)^2}+\frac{1}{2}\, \\= & {} -\frac{1}{2} \frac{\hat{\sigma }^2}{(\hat{m})^2}+\frac{1}{2}\, \\= & {} -\frac{1}{2}\left( \frac{\hat{\sigma }^2}{\hat{m}^2}-1 \right) . \end{aligned}$$

Finally the Hill number of order 1 can be computed as:

$$\begin{aligned} e^{H'}\approx & {} e^{\log (|\mathcal {S}|)-\frac{|\mathcal {S}|}{2}\sum _{i=1}^{|\mathcal {S}|}\delta _i^2}\,, \\\approx & {} |\mathcal {S}|\left( 1-\frac{|\mathcal {S}|}{2}\sum _{i=1}^{|\mathcal {S}|}\delta _i^2\right) \ =\ |\mathcal {S}|\left( 1-\frac{1}{2} \frac{\hat{\sigma }^2}{(\hat{m})^2}+\frac{1}{2}\right) \ =\ \frac{|\mathcal {S}|}{2}\left( 3- \frac{\hat{\sigma }^2}{(\hat{m})^2}\right) . \end{aligned}$$

Replacing $|{{\mathcal {S}}}|$ by $np^*$ and $\hat{\sigma }$ and $\hat{m}$ by their limits, we get the desired result:

$$\begin{aligned} e^{H'}\approx & {} \frac{np^*}{2}\left( 3- \frac{(\sigma ^*)^2}{(m^*)^2}\right) . \end{aligned}$$