Using Eq. (169) (and recalling that $\left\langle{\dots}\right\rangle_{i}$ has
turned into $\left\langle{\cdots}\right\rangle_{\star}$, and that we have made
the replacement $\underline{\underline{K}}\to i\underline{\underline{G}}$), we
then have

$\displaystyle C(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\left\langle{x(t)x(t^{\prime})}\right\rangle_{\star},$
$\displaystyle G(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-i\left\langle{x(t)\widehat{x}(t^{\prime})}\right\rangle_{\star}.$
(180)

With the definition of $\left\langle{\cdots}\right\rangle_{\star}$, and
recalling that $Z_{1}[\underline{\psi}=\underline{0},\underline{h}]=1$ in Eq.
(173) for any $\underline{h}$, we find

$G(t,t^{\prime})=\frac{\delta}{\delta
h(t^{\prime})}\left\langle{x(t)}\right\rangle_{\star}.$ (181)

This means that the self-consistent problem for the macroscopic order
parameters $\underline{M}$, $\underline{\underline{C}}$ and
$\underline{\underline{G}}$ is given by

$\displaystyle\dot{x}(t)$ $\displaystyle=$ $\displaystyle
x(t)\left[1-x(t)+\Gamma\sigma^{2}\int
dt^{\prime}G(t,t^{\prime})x(t^{\prime})+\mu M(t)+\eta(t)+h(t)\right],$
$\displaystyle M(t)$ $\displaystyle=$
$\displaystyle\left\langle{x(t)}\right\rangle_{\star},$ $\displaystyle
C(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\left\langle{x(t)x(t^{\prime})}\right\rangle_{\star},$
$\displaystyle G(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{\delta}{\delta
h(t^{\prime})}\left\langle{x(t)}\right\rangle_{\star},$ (182)

where $\left\langle{\cdots}\right\rangle_{\star}$ is an average over
realisation of the process in the first line. We have used the first relation
in Eq. (169) to write down the self-consistency relation for $\underline{M}$.

We will discuss the interpretation of this effective single-species process
and the associated self-consistency relations further in Sec. 9.

##  8 Generating functional analysis for random Lotka–Volterra equations III:
Proper derivation of single effective-species process

[For the purposes of this chapter, we return to the situation in which the
$\underline{\psi}_{i}$ and $\underline{h}_{i}$ can (at least potentially) be
different for different values of $i$.]

The dodgy step of our vanilla derivation consisted in simply setting
‘averages’ of hatted quantities to zero, at a point in the calculation where
we could not be sure that these were averages over a bona fide stochastic
process. More precisely, we set $\underline{P}=\underline{0}$ and
$\underline{\underline{L}}=\underline{\underline{0}}$ at the beginning of Sec.
7.1. Secondly, from the vanilla derivation alone, we cannot be sure what
exactly the order parameters $\underline{M},\underline{\underline{C}}$ and
$\underline{\underline{G}}$ mean in terms of the original Lotka–Volterra
process. Of course it is natural to assume that
$C(t,t^{\prime})=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\overline{\left\langle{x_{i}(t)}\right\rangle}$
for example, where the overbar describes average over realisation of the
disorder in the original problem, and $\left\langle{\dots}\right\rangle$ an
average over initial conditions (if they are random). We’d also suspect that
$M(t)=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\overline{\left\langle{x_{i}(t)}\right\rangle}$,
and $G(t,t^{\prime})=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta}{\delta
h_{i}(t^{\prime})}\overline{\left\langle{x_{i}(t)}\right\rangle}$. But we have
not formally proved any of this. The main purpose of Sec. 8.1 is to show
precisely this. Along the way we also provide a proper justification for
setting $\underline{P}=\underline{0}$ and
$\underline{\underline{L}}=\underline{\underline{0}}$. This then results in an
effective single-species measure and associated single-species process, as
described in Sec. 8.2

###  8.1 Identification of order parameters

Our derivation can be divided into three steps:

  * •

In step 1, we define order parameters for the original Lotka–Volterra problem
(as opposed to the effective process resulting from the saddle-point
calculation). These original order parameters will be called
$\underline{m},\underline{p},\underline{\underline{c}},\underline{\underline{k}}$
and $\underline{\underline{\ell}}$.

  * •

In step 2 we then write these original order parameters as derivatives of the
disorder-averaged generating functional.

  * •

In step 3 we then show that
$\underline{m},\underline{p},\underline{\underline{c}},\underline{\underline{k}}$
and $\underline{\underline{\ell}}$ can be expressed in exactly the same way as
the order parameters
$\underline{M},\underline{P},\underline{\underline{C}},\underline{\underline{K}}$
and $\underline{\underline{L}}$ which we introduced when we carried out the
saddle-point calculation. Thus, both sets of order parameters are identical at
the saddle point.

####  8.1.1 Step 1: Definition of order parameters from the original
Lotka–Volterra problem

The objects
$\underline{M},\underline{P},\underline{\underline{C}},\underline{\underline{K}}$
and $\underline{\underline{L}}$ were introduced in Eqs. (155) and we have
shown that they can be written in the form of Eqs. (169) in the thermodynamic
limit.

Let us now define the following

$\displaystyle m(t)$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\overline{\left\langle{x_{i}(t)}\right\rangle},$
$\displaystyle c(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\overline{\left\langle{x_{i}(t)x_{i}(t^{\prime})}\right\rangle},$
$\displaystyle k(t,t^{\prime})$ $\displaystyle=$ $\displaystyle
i\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta}{\delta
h_{i}(t^{\prime})}\overline{\left\langle{x_{i}(t)}\right\rangle},$
$\displaystyle p(t)$ $\displaystyle=$ $\displaystyle
i\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta}{\delta
h_{i}(t)}\overline{\left\langle{1}\right\rangle},$
$\displaystyle\ell(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta}{\delta
h_{i}(t)\delta h_{i}(t^{\prime})}\overline{\left\langle{1}\right\rangle},$
(183)

where the angle brackets stand for an average over realisations of the
original Lotka–Volterra dynamics in Eq. (1) for given disorder. In our case
the only randomness of the dynamics (for fixed disorder) is in the potentially
random initial conditions. However, one could, in principle, also include
dynamic noise. The overbar describes an average over realisations of the
disorder, i.e., the interaction matrix elements $\alpha_{ij}$ in Eq. (1).

We note the obvious fact $\overline{\left\langle{1}\right\rangle}=1$ for all
choices of $h_{i}(\cdot)$. This means that $p(t)$ and $\ell(t,t^{\prime})$ are
both zero by construction. We will use this to show that $P(t)=0$ and
$L(t,t^{\prime})=0$.

We re-iterate that the objects in Eqs. (183) are constructed from the original
dynamics. This is clear immediately for $m(t)$ and $c(t,t^{\prime})$. These
two objects can, in principle, be measured from numerical simulations, modulo
the fact that simulations will necessarily be for finite $N$ and that any
measurements will have to rely on a finite number of realisation runs. For a
given (large) number of species $N$ one runs a large number of realisations of
the random Lotka–Volterra system (each realisation with a new draw of the
interaction coefficients, and with a new initial condition). The object
$\frac{1}{N}\sum_{i}\overline{\left\langle{x_{i}(t)x_{i}(t^{\prime})}\right\rangle}$
for example can then be obtained from these simulations, as an average of the
object $\frac{1}{N}\sum_{i}x_{i}(t)x_{i}(t^{\prime})$ over realisations of the
disorder and the initial condition. Measuring $k(t,t^{\prime})$ is a little
more complicated (one would have to run the dynamics for different
$\underline{h}_{i}$), but nonetheless, $\underline{\underline{k}}$ can (at
least in principle) be obtained from simulating the original Lotka–Volterra
system.

####  8.1.2 Step 2: Order parameters for the original problem as derivatives
of the disorder-averaged generating functional

We know that moments of the $\\{x_{i}(t)\\}$ can be expressed as derivatives
of the generating functional with respect to the $\\{\psi_{i}(t)\\}$,
evaluated at $\underline{{\mbox{\boldmath$\psi$}}}=0$. We start with a fixed
realisation of the disorder in mind, so the word ‘moment’ refers to an average
over initial conditions only at this point. For example

$\left\langle{x_{i}(t)x_{i}(t^{\prime})}\right\rangle=-\left.\frac{\delta^{2}Z[\underline{{\mbox{\boldmath$\psi$}}}]}{\delta\psi_{i}(t)\delta\psi_{i}(t^{\prime})}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0},$
(184)

where $Z[\underline{{\mbox{\boldmath$\psi$}}}]$ is the generating functional
of the system for the relevant fixed realisation of the interaction matrix
[Eq. (141)].

Taking derivatives of the disorder-averaged generating functional instead
implements an additional average over the interaction coefficients. We have

$\displaystyle m(t)$ $\displaystyle=$
$\displaystyle-i\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta}{\delta\psi_{i}(t)}\left.\overline{Z[\underline{{\mbox{\boldmath$\psi$}}}]}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0},$
$\displaystyle c(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta^{2}}{\delta\psi_{i}(t)\delta\psi_{i}(t^{\prime})}\left.\overline{Z[\underline{{\mbox{\boldmath$\psi$}}}]}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0},$
$\displaystyle k(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta^{2}}{\delta\psi_{i}(t)\delta
h_{i}(t^{\prime})}\left.\overline{Z[\underline{{\mbox{\boldmath$\psi$}}}]}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0},$
$\displaystyle p(t)$ $\displaystyle=$ $\displaystyle
i\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta}{\delta
h_{i}(t^{\prime})}\left.\overline{Z[\underline{{\mbox{\boldmath$\psi$}}}]}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0},$
$\displaystyle\ell(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta^{2}}{\delta
h_{i}(t)\delta
h_{i}(t^{\prime})}\left.\overline{Z[\underline{{\mbox{\boldmath$\psi$}}}]}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0}.$
(185)

####  8.1.3 Step 3: Equivalence to order parameters at the saddle point

We now use the saddle-point expression for $\overline{Z}$ to show that the
expressions in Eqs. (8.1.2) are the same as those in Eqs. (169) at the saddle
point. We recall the definition of the average
$\left\langle{\cdots}\right\rangle_{i}$ in Eq. (170), as well as the
definition of the object $\Omega$ in Eq. (163). We repeat the latter here

$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle
N^{-1}\sum_{i}\ln\bigg{[}\int{\cal D}\underline{x}{\cal
D}\widehat{}\underline{x}\,p_{0}[x(0)]\exp\left(i\int
dt~{}\psi_{i}(t)x(t)\right)$ (186) $\displaystyle\times\exp\left(i\int
dt~{}\widehat{x}(t)\left(\frac{\dot{x}(t)}{x(t)}-[1-x(t)]-h_{i}(t)\right)\right)$
$\displaystyle\times\exp\left(-i\int
dt~{}dt^{\prime}\left[\widehat{C}(t,t^{\prime})x(t)x(t^{\prime})+\widehat{L}(t,t^{\prime})\widehat{x}(t)\widehat{x}(t^{\prime})+\widehat{K}(t,t^{\prime})x(t)\widehat{x}(t^{\prime})\right]\right)$
$\displaystyle\times\exp\left(-i\int
dt~{}[\widehat{M}(t)x(t)+\widehat{P}(t)\widehat{x}(t)]\right)\bigg{]}.$

From Eq. (159) we have

$\overline{Z}[\underline{{\mbox{\boldmath$\psi$}}}]=\int
D[\underline{M},\widehat{}\underline{M}]{\cal
D}[\underline{P},\widehat{}\underline{P}]{\cal
D}[\underline{\underline{C}},\widehat{}\underline{\underline{C}}]{\cal
D}[\underline{\underline{K}},\widehat{}\underline{\underline{K}}]{\cal
D}[\underline{\underline{L}},\widehat{}\underline{\underline{L}}]\,\exp\left(N\left[\Psi+\Phi+\Omega+{\cal
O}(N^{-1})\right]\right),$ (187)

noting that only $\Omega$ contains the $\\{\psi_{i}(t)\\}$ and
$\\{h_{i}(t)\\}$.

Let us focus on the expression for $c(t,t^{\prime})$ in Eqs. (8.1.2) as an
example. We have

$\displaystyle c(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta^{2}}{\delta\psi_{i}(t)\delta\psi_{i}(t^{\prime})}\left.\overline{Z[\underline{\psi}]}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0},$
(188) $\displaystyle=$
$\displaystyle-\lim_{N\to\infty}\frac{1}{N}\sum_{i}\frac{\delta}{\delta\psi_{i}(t)}\int
D[\underline{M},\widehat{}\underline{M}]{\cal
D}[\underline{P},\widehat{}\underline{P}]{\cal
D}[\underline{\underline{C}},\widehat{}\underline{\underline{C}}]{\cal
D}[\underline{\underline{K}},\widehat{}\underline{\underline{K}}]{\cal
D}[\underline{\underline{L}},\widehat{}\underline{\underline{L}}]\,$
$\displaystyle\times\bigg{\\{}\left.N\frac{\delta\Omega}{\delta\psi_{i}(t^{\prime})}\exp\left(N\left[\Psi+\Phi+\Omega+{\cal
O}(N^{-1})\right]\right)\bigg{\\}}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0}$
$\displaystyle=$ $\displaystyle-\lim_{N\to\infty}\frac{1}{N}\sum_{i}\int
D[\underline{M},\widehat{}\underline{M}]{\cal
D}[\underline{P},\widehat{}\underline{P}]{\cal
D}[\underline{\underline{C}},\widehat{}\underline{\underline{C}}]{\cal
D}[\underline{\underline{K}},\widehat{}\underline{\underline{K}}]{\cal
D}[\underline{\underline{L}},\widehat{}\underline{\underline{L}}]\,\bigg{\\{}$
$\displaystyle\left.\left[N\frac{\delta^{2}\Omega}{\delta\psi_{i}(t)\delta\psi_{i}(t^{\prime})}+N^{2}\left(\frac{\delta\Omega}{\delta\psi_{i}(t)}\right)\left(\frac{\delta\Omega}{\delta\psi_{i}(t^{\prime})}\right)\right]\right.$
$\displaystyle\left.\times\exp\left(N\left[\Psi+\Phi+\Omega+{\cal
O}(N^{-1})\right]\right)\bigg{\\}}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0}.$

Following the rules of saddle-point integration (and using the normalisation
of the generating functional,
$\left.\overline{Z[\underline{\psi}]}\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0}=1$),
the integral is given by the the object
$\left.N\frac{\delta^{2}\Omega}{\delta\psi_{i}(t)\delta\psi_{i}(t^{\prime})}+N^{2}\left(\frac{\delta\Omega}{\delta\psi_{i}(t)}\right)\left(\frac{\delta\Omega}{\delta\psi_{i}(t^{\prime})}\right)\right|_{\underline{{\mbox{\boldmath$\psi$}}}=0}$,
evaluated at the saddle point. We therefore find

$\displaystyle c(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left[N\frac{\delta^{2}\Omega}{\delta\psi(t)\delta\psi_{i}(t^{\prime})}+N^{2}\left(\frac{\delta\Omega}{\delta\psi_{i}(t)}\right)\left(\frac{\delta\Omega}{\delta\psi_{i}(t^{\prime})}\right)\right]\right|_{\mbox{SP},\underline{{\mbox{\boldmath$\psi$}}}=0},~{}~{}~{}~{}$
(189)

where the subscript ‘SP’ indicates that this is to be evaluated at the saddle
point.

We next use

$\Omega=\frac{1}{N}\sum_{j}\ln Z_{1}[\underline{\psi}_{j},\underline{h}_{j}],$
(190)

with $Z_{1}[\underline{\psi}_{j},\underline{h}_{j}]$ given in Eq. (167).
Therefore (and leaving out the argument $\underline{h}_{i}$ inside $Z_{1}$),
we have identities such as
$N\frac{\delta\Omega}{\delta\psi_{i}(t)}=\frac{\delta\ln
Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t)}$ for any fixed $i$. We can
then write

$\displaystyle c(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left[\frac{\delta^{2}\ln\,Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t)\delta\psi_{i}(t^{\prime})}+\left(\frac{\delta\ln\,Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t)}\right)\left(\frac{\delta\ln\,Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t^{\prime})}\right)\right]\right|_{\mbox{SP},\underline{{\mbox{\boldmath$\psi$}}}=0}.$

The first term in the square brackets is

$\displaystyle\frac{\delta^{2}\ln\,Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t)\delta\psi_{i}(t^{\prime})}$
$\displaystyle=$
$\displaystyle\frac{\delta}{\delta\psi_{i}(t)}\frac{1}{Z_{1}[\underline{\psi}_{i}]}\frac{\delta
Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t^{\prime})}$ (192)
$\displaystyle=$
$\displaystyle-\underbrace{\frac{1}{Z_{1}[\underline{\psi}_{i}]^{2}}\frac{\delta
Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t)}\frac{\delta
Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t^{\prime})}}_{=\left(\frac{\delta\ln\,Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t)}\right)\left(\frac{\delta\ln\,Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t^{\prime})}\right)}+\frac{1}{Z_{1}[\underline{\psi}_{i}]}\frac{\delta^{2}Z_{1}[\underline{\psi}_{i}]}{\delta\psi_{i}(t)\delta\psi_{i}(t^{\prime})}.$

Therefore,

$\displaystyle c(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\frac{1}{Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}\frac{\delta^{2}Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}{\delta\psi_{i}(t)\delta\psi_{i}(t^{\prime})}\right|_{\mbox{SP},\underline{{\mbox{\boldmath$\psi$}}}=0}.$
(193)

Using the definition of $Z_{1}$ in Eq. (167) as well as that of the average
$\left\langle{\cdots}\right\rangle_{i}$ in Eq. (170), we can then take the
remaining derivatives, and find

$\displaystyle
c(t,t^{\prime})=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left\langle{x(t)x(t^{\prime})}\right\rangle_{i}\right|_{\mbox{SP},\underline{\psi}_{i}=0}$
(194)

Comparing with Eq. (169), we thus conclude that the quantity
$c(t,t^{\prime})$, defined from the original Lotka–Volterra problem, is the
same as the order parameter $C(t,t^{\prime})$ we introduced in the course of
the generating functional calculation, provided the expressions in Eq. (169)
are evaluated at the saddle-point and at $\mathbf{\underline{\psi}}=0$.
Analogous equivalences for the remaining order parameters can be demonstrated
along similar lines.

We can therefore summarise (all valid at the saddle point and for
$\mathbf{\underline{\psi}}=0$),

$\displaystyle\lim_{N\to\infty}M(t)$ $\displaystyle=$ $\displaystyle
m(t)=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left\langle{x(t)}\right\rangle_{i}\right|_{\mbox{SP},\underline{\psi}_{i}=0},$
$\displaystyle\lim_{N\to\infty}P(t)$ $\displaystyle=$ $\displaystyle
p(t)=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left\langle{\widehat{x}(t)}\right\rangle_{i}\right|_{\mbox{SP},\underline{\psi}_{i}=0}=0,$
$\displaystyle\lim_{N\to\infty}C(t,t^{\prime})$ $\displaystyle=$
$\displaystyle
c(t,t^{\prime})=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left\langle{x(t)x(t^{\prime})}\right\rangle_{i}\right|_{\mbox{SP},\underline{\psi}_{i}=0},$
$\displaystyle\lim_{N\to\infty}K(t,t^{\prime})$ $\displaystyle=$
$\displaystyle
k(t,t^{\prime})=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left\langle{x(t)\widehat{x}(t^{\prime})}\right\rangle_{i}\right|_{\mbox{SP},\underline{\psi}_{i}=0},$
$\displaystyle\lim_{N\to\infty}L(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\ell(t,t^{\prime})=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left\langle{\widehat{x}(t)\widehat{x}(t^{\prime})}\right\rangle_{i}\right|_{\mbox{SP},\underline{\psi}_{i}=0}=0.$
(195)

We have therefore demonstrated two relevant facts: (i) In the thermodynamic
limit and at the saddle point, the order parameters
$\underline{M},\underline{\underline{C}}$ and $\underline{\underline{K}}$,
initially introduced as bookkeeping variables in the generating-functional
calculation have an interpretation in terms of averages in the original
Lotka–Volterra problem. (ii) The order parameters $\underline{P}$ and
$\underline{\underline{L}}$ vanish at the saddle point.

###  8.2 Effective process for a representative species

The fields $\underline{\psi}_{i}$ have now served their purpose and we can set
them to zero. We also restrict the perturbation fields $h_{i}(t)$ to be
uniform across species, that is, replace $h_{i}(t)\to h(t)$. Recalling that
the only $i$-dependence of the average $\left\langle{\dots}\right\rangle_{i}$
is through the fields $h_{i}(t)$ and $\underline{\psi}_{i}$ [Eq. (170)] we
then find that the operation $\left\langle{\dots}\right\rangle_{i}$ becomes
independent of $i$. We simply write $\left\langle{\dots}\right\rangle_{*}$ for
this average, and include in the definition of
$\left\langle{\cdots}\right\rangle_{\star}$ that this is to be evaluated at
the saddle point.

We then have

$\left\langle{g[\underline{x},\widehat{}\underline{x}]}\right\rangle_{*}=\frac{\int{\cal
D}\underline{x}{\cal
D}\widehat{}\underline{x}\,g[\underline{x},\widehat{}\underline{x}]\,{\cal
M}[\underline{x},\widehat{}\underline{x}]}{\int{\cal D}\underline{x}{\cal
D}\widehat{}\underline{x}\,\,{\cal M}[\underline{x},\widehat{}\underline{x}]}$
(196)

with

$\displaystyle{\cal M}[\underline{x},\widehat{}\underline{x}]$
$\displaystyle=$ $\displaystyle p_{0}[x(0)]\,\times\exp\left(i\int
dt~{}\widehat{x}(t)\left(\frac{\dot{x}(t)}{x(t)}-[1-x(t)]-h(t)\right)\right)$
(197) $\displaystyle\times\exp\left(-i\int
dt~{}dt^{\prime}\left[\widehat{C}(t,t^{\prime})x(t)x(t^{\prime})+\widehat{L}(t,t^{\prime})\widehat{x}(t)\widehat{x}(t^{\prime})+\widehat{K}(t,t^{\prime})x(t)\widehat{x}(t^{\prime})\right]\right)$
$\displaystyle\times\left.\exp\left(-i\int
dt~{}[\widehat{M}(t)x(t)+\widehat{P}(t)\widehat{x}(t)]\right)\right|_{\mbox{SP}}.$

At the saddle point we have from Eq. (164) and using
$\underline{P}=\underline{0},\underline{\underline{L}}=\underline{\underline{0}}$,

$\displaystyle i\widehat{M}(t)$ $\displaystyle=$ $\displaystyle i\mu P(t)=0,$
$\displaystyle i\widehat{P}(t)$ $\displaystyle=$ $\displaystyle i\mu M(t),$
$\displaystyle i\widehat{C}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sigma^{2}L(t,t^{\prime})=0,$ $\displaystyle
i\widehat{K}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\Gamma\sigma^{2}K(t^{\prime},t),$ $\displaystyle
i\widehat{L}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sigma^{2}C(t,t^{\prime}).$ (198)

Using this, and again implementing the change of notation
$\underline{\underline{K}}\to i\underline{\underline{G}}$, we find

$\displaystyle{\cal M}[\underline{x},\widehat{}\underline{x}]$
$\displaystyle=$ $\displaystyle p_{0}[x(0)]\,\exp\left(i\int
dt~{}\widehat{x}(t)\left(\frac{\dot{x}(t)}{x(t)}-[1-x(t)]-\mu
M(t)-h(t)\right)\right)$ (199) $\displaystyle\times\exp\left(-\sigma^{2}\int
dt~{}dt^{\prime}\left[\frac{1}{2}C(t,t^{\prime})\widehat{x}(t)\widehat{x}(t^{\prime})+i\Gamma
G(t^{\prime},t)x(t)\widehat{x}(t^{\prime})\right]\right).$

Following the lines of Sec. 7.2, the single effective species measure ${\cal
M}[\underline{x},\widehat{}\underline{x}]$ again describes the effective
process

$\displaystyle\dot{x}(t)=x(t)\left[1-x(t)+\Gamma\sigma^{2}\int
dt^{\prime}G(t,t^{\prime})x(t^{\prime})+\mu M(t)+\eta(t)+h(t)\right].$ (200)

From Eq. (169) (and recalling that
$\lim_{N\to\infty}N^{-1}\sum_{i}\left\langle{\cdots}\right\rangle_{i}$ has
turned into $\left\langle{\cdots}\right\rangle_{\star}$), we also have the
self-consistency relations

$\displaystyle M(t)$ $\displaystyle=$
$\displaystyle\left\langle{x(t)}\right\rangle_{*},$
$\displaystyle\left\langle{\eta(t)\eta(t^{\prime})}\right\rangle_{\star}$
$\displaystyle=$
$\displaystyle\sigma^{2}C(t,t^{\prime})~{}\mbox{with}~{}C(t,t^{\prime})=\left\langle{x(t)x(t^{\prime})}\right\rangle_{*},$
$\displaystyle G(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-i\left\langle{x(t)\widehat{x}(t^{\prime})}\right\rangle_{*},$
(201)

and where the last relation can be written in the form
$G(t,t^{\prime})=\frac{\delta}{\delta
h(t^{\prime})}\left\langle{x(t)}\right\rangle_{\star}$.

###  8.3 Comments on the legitimacy of the vanilla approach

The main deficiency of the vanilla approach was to set the order parameters
$\underline{P}$ and $\underline{\underline{L}}$ to zero without proper
justification. Nonetheless, the vanilla method leads to the right effective
process. Once this process is established, the self-consistency relations for
$\underline{P}$ and $\underline{\underline{L}}$ result in the conclusion that
$\underline{P}$ and $\underline{\underline{L}}$ must vanish (because they are
averages of conjugates field over a bona fide stochastic dynamics, namely the
effective process). The vanilla approach effectively amounts to simply
asserting $\underline{P}=\underline{0}$ and
$\underline{\underline{L}}=\underline{\underline{0}}$, and then justifying
them retrospectively from the effective process. This is not incorrect, albeit
a little weaker than the bottom-up derivation from first principles in Sec. 8.

The second issue is that we have kept the source fields $\psi_{i}(t)$ general
throughout Sec. 7, even if setting $\underline{P}$ and
$\underline{\underline{L}}$ to zero is at best justified in the limit
$\underline{{\mbox{\boldmath$\psi$}}}\to\underline{\mathbf{0}}$. This is not a
serious problem though. In fact we are not using the source fields in Sec. 7,
so no harm would have been done if we had simply set them to zero from the
beginning. From the vanilla approach we would then not obtain a generating
functional, but instead the single effective particle measure ${\cal
M}[\underline{x},\widehat{}\underline{x}]$ in Eq (199). But of course the
generating functional and the single effective particle measure are really two
sides of the same coin, and describe the same effective process (the
generating functional is the Fourier transform of the probability measure in
the space of paths).

As an overall conclusion, it is therefore fair to say that the vanilla
approach is OK, but only if it is used responsibly and if the user is aware of
the subtleties.

##  9 Generating functional analysis for random Lotka–Volterra equations IV:
Dynamic mean-field process and its interpretation

We have now arrived at the dynamic mean-field description of the
Lotka–Volterra system. This description consists of the effective single-
species process

$\displaystyle\dot{x}(t)=x(t)\left[1-x(t)+\Gamma\sigma^{2}\int
dt^{\prime}G(t,t^{\prime})x(t^{\prime})+\mu M(t)+\eta(t)+h(t)\right],$ (202)

along with the self-consistency relations

$\displaystyle\sigma^{2}C(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\sigma^{2}\left\langle{x(t)x(t^{\prime})}\right\rangle_{*}=\left\langle{\eta(t)\eta(t^{\prime})}\right\rangle_{*},$
$\displaystyle G(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{\delta}{\delta
h(t^{\prime})}\left\langle{x(t)}\right\rangle_{*},$ $\displaystyle M(t)$
$\displaystyle=$ $\displaystyle\left\langle{x(t)}\right\rangle_{*}.$ (203)

The notation $\left\langle{\cdots}\right\rangle_{*}$ denotes an average over
realizations of the effective dynamics (202).

Technically, the first line in Eqs. (203) entails two different statements:
One is the actual self-consistency relation
$\sigma^{2}\left\langle{x(t)x(t^{\prime})}\right\rangle_{*}=\left\langle{\eta(t)\eta(t^{\prime})}\right\rangle_{*}$.
The second is the fact that the order parameter $C(t,t^{\prime})$
($=c(t,t^{\prime})$ in the original problem) is obtained from an average of
the object $x(t)x(t^{\prime})$ over realisations of the effective process,
i.e., $C(t,t^{\prime})=\left\langle{x(t)x(t^{\prime})}\right\rangle_{*}$.

We’ll now briefly talk about the mathematical structure of this self-
consistent problem, and then about its interpretation and relation to the
original Lotka–Volterra system.

###  9.1 Mathematical structure

First of all one needs to realise that Eqs. (202,203) are primarily a set of
equations determining the dynamical order parameters $M(t),C(t,t^{\prime})$
and $G(t,t^{\prime})$. They are the analog of Eq. (15) in the kinetic Ising
model that we discussed as a warm-up, and of Eqs. (135) and (136) in the toy
problem in Sec. 5.

Eqs. (202,203) do not provide explicit expressions for the order parameters.
Instead this is an implicit system, that needs to be solved self-consistently.
This means that the stochastic (integro-) differential equation for the
effective species abundance $x(t)$ involves the order parameters – the
quantities $M(t)$ and $G(t,^{\prime}t)$ appear on the right-hand side of Eq.
(202). The properties of the noise $\eta(t)$, also on the right-hand side of
the effective process, are determined by $\underline{\underline{C}}$, where
$C(t,t^{\prime})=\left\langle{x(t)x(t^{\prime})}\right\rangle_{*}$.

Conversely, Eqs. (203) indicate that $C(t,t^{\prime}),G(t,t^{\prime})$ and
$M(t)$ are to be obtained as averages over realisations of the effective
process. So we have here a Catch-22 situation, we need the order parameters to
construct realisations of the effective process, and at the same time, we need
these realisations to know what the order parameters are. This is what we mean
by a problem that has to be solved self-consistently.

How could one go about solving this set of equations? We’ll look for fixed
points of the effective dynamics in Sec. 10 – this can be done analytically.
Other than that, one would have to do this numerically. One possible method is
to proceed along the following steps:

  1. 1.

Make an initial guess for the functions $C(t,t^{\prime}),G(t,t^{\prime})$ and
$M(t)$.

  2. 2.

Use $C(t,t^{\prime})$ to generate realisations of the noise $\eta(t)$.

  3. 3.

Use these realisations, and the guesses for $G(t,t^{\prime})$ and $M(t)$ to
generate realisations of the process $x(t)$.

  4. 4.

From these realisations of $x(t)$, obtain new functions
$C(t,t^{\prime}),G(t,t^{\prime})$ and $M(t)$ via the self-consistently
relations.

  5. 5.

Repeat steps 2-4 until convergence (i.e., until the order parameters do not
change any more).

I would not say that these steps are easy, and I am not really going to go
into how one would actually do this (determining the response function
$G(t,t^{\prime})$ in step 4 for example is non-trivial). The point is here
just to explain the nature of the self-consistent problem. A procedure along
those lines is described in [16].

Alternatively, there is a celebrated method by Eissfeller and Opper [53]. This
is for discrete-time processes, and builds up the matrices $C(t,t^{\prime})$
and $G(t,t^{\prime})$ gradually, time step by time step. Realisations of
$x(t)$ and $\eta(t)$ are also constructed time step by time step in the
Eissfeller–Opper scheme.

###  9.2 Interpretation

This is a good point to briefly think about the interpretation of the
effective process (the mean field dynamics). Suppose we have found the order
parameters self-consistently. With those
$\underline{\underline{C}},\underline{\underline{G}}$ and $\underline{M}$ the
effective process in Eq. (202) then describes what a typical or representative
species ‘experiences’. The right-hand side is the change of abundance in time
of such a ‘typical’ species. The word ‘typical’ does not mean an ‘average’
over species. The right-hand side of the effective dynamics contains
randomness (the noise $\eta(t)$). This noise comes from the original problem –
the quenched disorder. For fixed disorder, different species in the original
problem will experience different right-hand sides of Eq. (1). The noise
$\eta(t)$ in the effective process reflects this.

To be a little more precise, the statement is the following. Eqs. (202,203),
if solved self-consistently, generate an ensemble of realisations $x(t)$ of
the effective process. The statistics of this ensemble is the same as that of
trajectories $x_{i}(t)$ one obtains in the original Lotka–Volterra problem
(for large $N$). We can phrase this slightly more formally as follows:

Equivalence of averages in original problem and at the level of single
effective species.
For any function $f[\underline{x}_{i},\widehat{}\underline{x}_{i}]$ we have
(at the saddle point and for
$\underline{{\mbox{\boldmath$\psi$}}}=\underline{\mathbf{0}}$),

$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\overline{\left\langle{f[\underline{x}_{i},\widehat{}\underline{x}_{i}]}\right\rangle}=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left\langle{f[\underline{x},\widehat{}\underline{x}]}\right\rangle_{i}\right|_{\mbox{SP},\underline{\psi}_{i}=0}.$
(204)

This is demonstrated in Appendix B.

We can say the same thing again in different words. Suppose I fix $\mu,\Gamma$
and $\sigma^{2}$, and I choose a sufficiently large $N$. I then draw a
realisation of the disorder, and integrate the Lotka–Volterra equations (1).
This coupled system produces trajectories
$\underline{x}_{1},\underline{x}_{2},\dots,\underline{x}_{N}$, one for each
species. One can therefore think of a distribution of single-species
trajectories. Modulo finite-size effects ($N$ will obviously have to be finite
for the numerical integration of the Lotka–Volterra equations) this
distribution is the same as the distribution of trajectories for the effective
process (at the self-consistent solution).

Figure 3: Illustration of the relation of the different sections on the
Lotka–Volterra problem.

###  9.3 Summary so far

In Sec. 6 we set up the generating functional for the Lotka–Volterra problem,
and then wrote the disorder-averaged generating functional in the form of a
saddle-point integral. The steps involved carrying out the disorder average,
and the introduction order parameters as book-keeping variables. We then
derived the saddle-point relations for the order parameters.

There are then two ways to proceed. Along the easier route (Sec. 7) we simply
assert that the order parameters $\underline{P}$ and
$\underline{\underline{L}}$ vanish, and we then arrive directly at the
effective dynamics. A more complete and longer route (Sec. 8) involves the
definition of order parameters at the level of the original problem,
identifying these with the book-keeping variables at the saddle point, and
explicitly showing that averages of hatted variables do indeed vanish (i.e.,
$\underline{P}=\underline{0}$ and
$\underline{\underline{L}}=\underline{\underline{0}}$).

Either way one arrives at the same dynamical mean field theory (Sec. 9).

Our next steps are now as follows: In Sec. 10 we make a fixed-point ansatz for
the effective process and derive analytical relations for the resulting order
parameters at the fixed point. In Sec. 11 we then study the stability of the
fixed point, and obtain a phase diagram for the Lotka–Volterra problem,
indicating for what combinations pf $\mu,\sigma^{2}$ and $\Gamma$ the system
is stable, and what types of instability occur where in parameter space.

The relation of these different sections to one another is also illustrated in
Fig. 3.

##  10 Fixed point analysis

The purpose of the perturbation fields $h_{i}(t)$ was to generate response
functions. The fields are not part of the actual model. Therefore, with the
exception of Sec. 10.2.3, we always set $h_{i}(t)\equiv 0$ throughout this
section.

###  10.1 Fixed-point ansatz and its consequences

####  10.1.1 Fixed-point ansatz for effective abundances

We now assume that the system reaches a stationary state and that this
stationary state does not depend on the initial condition (i.e., we assume
absence of long-term memory, see also [14]). The response function $G$ is then
a function of time differences only, i.e. $G(t,t^{\prime})=G(\tau)$, where
$\tau=t-t^{\prime}$. Causality dictates that $G(\tau<0)=0$.

We assume further that the dynamics reaches a fixed point. In the original
system this means that all $x_{i}$ tend to fixed-point values $x_{i}^{*}$ at
long times. At the level of the effective dynamics, all trajectories of the
single-species process tend to constant values $x(t)\to x^{*}$, but it is
important to realise that different realisations of the effective process will
have different fixed-point values. The quantity $x^{*}$ is therefore a static
(time-independent) but random quantity. For given model parameters
$\mu,\Gamma$ and $\sigma^{2}$ the distribution of $x^{*}$ corresponds to the
distribution of the $x_{i}^{*}$ in the original problem.

As a consequence the order parameter
$C(t,t^{\prime})=\left\langle{x(t)x(t^{\prime})}\right\rangle_{*}=\left\langle{(x^{*})^{2}}\right\rangle$
is constant (independent of $t$ and $t^{\prime}$). We write

$C(t,t^{\prime})\equiv q.$ (205)

####  10.1.2 The noise variable in the effective process becomes static

We next note that the relation

$\left\langle{\eta(t)\eta(t^{\prime})}\right\rangle_{*}=\sigma^{2}C(t,t^{\prime})$
(206)

in Eq. (203) means that
$\left\langle{\eta(t)\eta(t^{\prime})}\right\rangle_{*}=\sigma^{2}q$ for all
$t,t^{\prime}$. The only way this can be the case is if the random variable
$\eta(t)$ itself is independent of time. To see this, look at

$\displaystyle\left\langle{[\eta(t)-\eta(t^{\prime})]^{2}}\right\rangle_{*}$
$\displaystyle=$
$\displaystyle\left\langle{(\eta(t))^{2}}\right\rangle_{*}+\left\langle{(\eta(t^{\prime}))^{2}}\right\rangle_{*}-2\left\langle{\eta(t)\eta(t^{\prime})}\right\rangle_{*}$
(207) $\displaystyle=$ $\displaystyle
C(t,t)+C(t^{\prime},t^{\prime})-2C(t,t^{\prime}).$

If $C(t,t)$, $C(t^{\prime},t^{\prime})$ and $C(t,t^{\prime})$ are all equal to
$\sigma^{2}q$ for times $t,t^{\prime}$ in the stationary regime, then
$\left\langle{[\eta(t)-\eta(t^{\prime})]^{2}}\right\rangle_{*}=0$. If the
average of a square is zero, the object that is squared must be zero itself
with probability one. In other words, we find that $\eta(t)=\eta(t^{\prime})$
with probability one. Therefore, $\eta(t)$ is in effect a static random
variable, we write $\eta(t)=\eta^{*}$. The variance of $\eta^{*}$ is
$\sigma^{2}q$. The mean is zero (we already knew that
$\left\langle{\eta(t)}\right\rangle_{*}=0$ for all $t$).

####  10.1.3 Distribution of fixed-point abundances

We now introduce the fixed-point ansatz into the effective process in Eq.
(202). In doing this we note that $\int
dt^{\prime}G(t,t^{\prime})x(t^{\prime})$ becomes
$\int_{-\infty}^{t}dt^{\prime}G(t-t^{\prime})x^{*}$, where we have sent the
starting point of the dynamics to time $-\infty$ (this has the advantage that
any finite time $t$ is now in the fixed-point regime). We can then replace

$\int dt^{\prime}G(t,t^{\prime})x(t^{\prime})\to
x^{*}\int_{0}^{\infty}d\tau\,G(\tau)=x^{*}\chi,$ (208)

where we have defined

$\chi\equiv\int_{0}^{\infty}d\tau\,G(\tau).$ (209)

We then find from Eq. (202), setting $h(t)\equiv 0$ (as explained at the
beginning of Sec. 10),

$x^{*}\left[1-x^{*}+\Gamma\sigma^{2}\chi x^{*}+\mu M^{*}+\eta^{*}\right]=0.$
(210)

Keeping in mind that
$\left\langle{(\eta^{*})^{2}}\right\rangle_{*}=\sigma^{2}q$ we write
$\eta^{*}=\sqrt{q}\sigma z$ with $z$ a static Gaussian random variable of mean
zero and unit variance.

Eq. (210) always (for all $z$) has the solution $x^{*}=0$. The second
solution,

$x^{*}=\frac{1+\mu M^{*}+\sqrt{q}\sigma z}{1-\Gamma\sigma^{2}\chi}$ (211)

can only be physical for those values of $z$ for which this expression is non-
negative ($x^{*}$ is a fixed-point abundance, and must be non-negative). In
the following we use

$x(z)=\mbox{max}\left(0,\frac{1+\mu M^{*}+\sqrt{q}\sigma
z}{1-\Gamma\sigma^{2}\chi}\right),$ (212)

that is to say, we use the expression in Eq. (211) whenever this expression is
positive, and $x^{*}=0$ for values of $z$ for which the right-hand side of Eq.
(211) is negative (or equal to zero). The zero solution can be seen to be
unstable when the right-hand side of Eq. (211) is positive, see Sec. 11.2
below.

###  10.2 Self-consistency relations for order parameters

####  10.2.1 Summary of status so far

We have now achieved the following:

  1. 1.

We have made a fixed-point ansatz for the effective process.

  2. 2.

We have established the consequences for the macroscopic order parameters.
Only the quantities $q,M^{*}$ and $\chi$ (which are derived from the order
parameters) are relevant at the fixed point.

  3. 3.

We have expressed the fixed points of the effective process in terms of the
order parameters and a static noise variable $z$ (Gaussian of mean zero and
variance one) [Eq. (212)].

This almost completes our fixed-point simplifications of the self-consistent
problem, but not quite. We still have to express the quantities $M^{*},q$ and
$\chi$ as averages over the fixed-point abundances $x(z)$, or objects derived
from these abundances. This is what we will do next.

####  10.2.2 Self-consistent expressions for order parameters $M^{*}$ and $q$

The order parameter $M(t)$ is equal to $\left\langle{x(t)}\right\rangle_{*}$
at finite times, so it is clear that $M^{*}$ must be the average over the
fixed-point abundances. We have

$M^{*}=\left\langle{x(z)}\right\rangle_{*},$ (213)

where we note that, at the fixed-point, the average
$\left\langle{\cdots}\right\rangle_{*}$ is really an average over the Gaussian
random variable $z$. The functional form of $x(z)$ is given in Eq. (212).
Similarly, the quantity $q$ is the value the correlation function
$C(t,t^{\prime})=\left\langle{x(t)x(t^{\prime})}\right\rangle_{*}$ takes in
the fixed-point regime (where $C(t,t^{\prime})$ becomes independent of $t$ and
$t^{\prime}$ as discussed above). We therefore have

$q=\left\langle{x(z)^{2}}\right\rangle_{*}.$ (214)

####  10.2.3 Self-consistent relation for the integrated response $\chi$

[For the purposes of this subsection (and only this subsection) we allow non-
zero perturbation fields again. We impose $h_{i}(t)\equiv h(t)$ for all $i$.]

We have to work a little harder to express $\chi$ in terms of an average over
something to do with $x(z)$. First we remember that the response function
$G(t,t^{\prime})$ depends only on time differences in the stationary state,
and in absence of memory effects. We can write $G(\tau)=G(t,t-\tau)$, that is
$G(\tau)$ describes the average linear response of abundances at a time $t$ in
the stationary state to a (small) perturbation which is applied $\tau$ units
of time earlier,

$G(\tau)=\left.\frac{\delta}{\delta
h(t-\tau)}\left\langle{x(t)}\right\rangle_{*}\right|_{\underline{h}=\underline{0}}~{}\mbox{(with
$t$ in the stationary regime)}.$ (215)

The quantity $\chi$ is the integrated response function in the stationary
state, $\chi=\int_{0}^{\infty}d\tau\,G(\tau)$. It is helpful to briefly
imagine what this would look like in discrete time (with time step equal to
one). We’d have $\chi=\sum_{\tau=1}^{\infty}G(\tau)$, and $G(\tau)$ describes
the response of $x(t)$ to a perturbation at $t-\tau$ (with $t$ sufficiently
long so that all this happens in the stationary state). The point is that for
the discrete-time system

$\chi=\sum_{\tau>0}^{\infty}G(\tau)=\sum_{\tau>0}\left.\frac{\delta}{\delta
h(t-\tau)}\left\langle{x(t)}\right\rangle_{*}\right|_{\underline{h}=\underline{0}}~{}\mbox{(with
$t$ in the stationary regime)},$ (216)

so $\chi$ sums up the effects on $x(t)$ of a perturbation at time $t-1$, a
perturbation at time $t-2$, at $t-3$ and so on. So one would have to apply a
small perturbation at all times. The integrated response, in the fixed-point
regime, tells us how the fixed point value for $x$ is affected by a small
constant perturbation, acting at all times. This is illustrated in Fig. 4.

Figure 4: Illustration of the interpretation of the integrated response
$\chi=\sum_{\tau>0}G(\tau)$ in a discrete-time system. $G(\tau=1)$ describes
the average effect on $x(t)$ of a small perturbation one time step earlier
(upper panel), $G(2)$ is the average effect two time steps after a
perturbation has been applied, and $G(3)$ three steps after the perturbation.
Thus, $G(1)+G(2)+G(3)$ describes the linear response at a time $t$ when a
perturbation is applied simultaneously at times $t-1,t-2$ and $t-3$ (fourth
panel). The integrated response $\chi$ is then the response of the system if a
constant perturbation is applied at all earlier times. Now you may ask: OK, if
I perturb the system at some time $t_{1}$ and then I perturb the same system
again at some later time $t_{2}$, then surely the later perturbation at
$t_{2}$ acts on an already perturbed signal. Don’t I have to take this into
account somehow? Why are all derivatives in Eq. (216) to be taken at
$\underline{h}=\underline{0}$? The answer is that we are looking at linear
response, so we are not considering perturbations of perturbations. Instead we
are summing up the small effects of perturbations at different times
independently..

OK, so what we need to do it to apply a constant perturbation
$h(t^{\prime})\equiv h$. Instead of Eq. (210) we the find

$x^{*}\left[1-x^{*}+\Gamma\sigma^{2}\chi x^{*}+\mu
M^{*}+\eta^{*}{\color[rgb]{0,0,1}+h}\right]=0,$ (217)

and accordingly,

$x(z)=\mbox{max}\left(0,\frac{1+\mu M^{*}+\sqrt{q}\sigma
z{\color[rgb]{0,0,1}+h}}{1-\Gamma\sigma^{2}\chi}\right),$ (218)

instead of Eq. (212).

Following the above argument, we then have

$\chi=\left.\frac{d}{dh}\left\langle{x(z)}\right\rangle_{*}\right|_{h=0},$
(219)

with $x(z)$ given by Eq. (218). The quantity $1-\Gamma\sigma^{2}\chi$ will
always come out positive, so we proceed here assuming that this will be the
case. The second argument in the maximum in Eq. (218) is then positive if and
only if $z>-\frac{1+\mu M^{*}{\color[rgb]{0,0,1}+h}}{\sigma\sqrt{q}}$.
Therefore

$\left\langle{x(z)}\right\rangle_{*}=\int_{-\frac{1+\mu
M^{*}{\color[rgb]{0,0,1}+h}}{\sigma\sqrt{q}}}^{\infty}dz\,e^{-z^{2}/2}\frac{1+\mu
M^{*}+\sqrt{q}\sigma z{\color[rgb]{0,0,1}+h}}{1-\Gamma\sigma^{2}\chi}.$ (220)

From this we find

$\displaystyle\left.\frac{d}{dh}\left\langle{x(z)}\right\rangle_{*}\right|_{h=0}$
$\displaystyle=$
$\displaystyle\frac{1}{\sigma\sqrt{q}}\times\left(\mbox{integrand in
Eq.~{}(\ref{eq:chi_aux}) evaluated at $h=0$ and~{}}z=-\frac{1+\mu
M^{*}}{\sigma\sqrt{q}}\right)$ (221)
$\displaystyle+\frac{1}{1-\Gamma\sigma^{2}\chi}\int_{-\frac{1+\mu
M^{*}}{\sigma\sqrt{q}}}^{\infty}dz\,e^{-z^{2}/2}.$

The integrand in Eq. (220) vanishes when $h=0$ and $z=-\frac{1+\mu
M^{*}}{\sigma\sqrt{q}}$, therefore we are left with

$\chi=\frac{1}{1-\Gamma\sigma^{2}\chi}\int_{-\frac{1+\mu
M^{*}}{\sigma\sqrt{q}}}^{\infty}dz\,e^{-z^{2}/2}.$ (222)

####  10.2.4 Summary of self-consistency relations at the fixed point

[We now set $h(t)\equiv 0$ again.]

We had already established that $M^{*}=\left\langle{x(z)}\right\rangle_{*}$
and $q=\left\langle{(x(z))^{2}}\right\rangle_{*}$ at the fixed point. Only
non-zero $x(z)$ contribute to these averages. It is now helpful to introduce

$\Delta\equiv\frac{1+\mu M^{*}}{\sigma\sqrt{q}}.$ (223)

We then have

$\displaystyle M^{*}$ $\displaystyle=$
$\displaystyle\int_{-\Delta}^{\infty}Dz~{}x(z),$ $\displaystyle q$
$\displaystyle=$ $\displaystyle\int_{-\Delta}^{\infty}Dz~{}x(z)^{2},$ (224)

and

$\chi=\frac{1}{1-\Gamma\sigma^{2}\chi}\int_{-\Delta}^{\infty}Dz$ (225)

with the definition $Dz=e^{-z^{2}/2}/\sqrt{2\pi}$. In the integration range
$z>-\Delta$ we also have [from Eq. (212)]

$x(z)=\sqrt{q}\sigma\frac{\Delta+z}{1-\Gamma\sigma^{2}\chi}.$ (226)

Explicitly writing things out we therefore have

$\displaystyle\chi$ $\displaystyle=$
$\displaystyle\frac{1}{1-\Gamma\sigma^{2}\chi}\int_{-\Delta}^{\infty}Dz,$
$\displaystyle M^{*}$ $\displaystyle=$
$\displaystyle\sqrt{q}\sigma\frac{1}{1-\Gamma\sigma^{2}\chi}\int_{-\Delta}^{\infty}Dz~{}(\Delta+z),$
$\displaystyle 1$ $\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{(1-\Gamma\sigma^{2}\chi)^{2}}\int_{-\Delta}^{\infty}Dz~{}(\Delta+z)^{2}.$
(227)

Remark:
It isn’t particularly important for the further analysis, but nonetheless
useful to note that, in short form, the self-consistency relations can be
written as

$\displaystyle\chi$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{q}\sigma}\left\langle{\frac{\partial
x(z)}{\partial z}}\right\rangle_{*},$ $\displaystyle M^{*}$ $\displaystyle=$
$\displaystyle\left\langle{x(z)}\right\rangle_{*},$ $\displaystyle q$
$\displaystyle=$ $\displaystyle\left\langle{(x(z))^{2}}\right\rangle_{*}.$
(228)

with $x(z)$ given in Eq. (212).

Further remark (assumptions we have made for the fixed-point ansatz)
Eqs. (227) describe the system in the phase of unique stable fixed points.
What this means is that (i) the original Lotka–Volterra system approaches a
stable fixed point, and (ii) for any fixed realisation of the disorder this
fixed point does not depend on initial conditions (otherwise the system would
not be memory-less, and there could be different solutions for the order
parameters for different initial conditions). As we will describe below these
assumptions do not hold for all choices of the model parameters
$\mu,\sigma^{2}$ and $\Gamma$.

But even within the phase of unique stable fixed points some qualifications
are required. In this phase the above assumptions are valid in the
thermodynamic limit, i.e., for infinitely many species $N\to\infty$. In
practice, any numerical integration of the Lotka–Volterra equations (1) will
necessarily have to be carried out for a finite $N$. This means that in
simulations for parameters $\mu,\sigma^{2},\Gamma$ in the phase of unique
fixed points one may not always find a unique fixed point for all samples of
the disorder. Some samples might not converge to a fixed point, or perhaps one
finds multiple fixed points. But the fraction of these samples of the disorder
will go to zero for larger and larger system sizes $N$.

###  10.3 Further analysis

We now change integration variable $z$ into $-\mathpzc{z}$ (the symbol
$\mathpzc{z}$ is also pronounced ‘zet’). Doing this in each integral in Eqs.
(227) we have

$\displaystyle\chi$ $\displaystyle=$
$\displaystyle\frac{1}{1-\Gamma\sigma^{2}\chi}\int_{-\infty}^{\Delta}D\mathpzc{z},$
$\displaystyle M^{*}$ $\displaystyle=$
$\displaystyle\sqrt{q}\sigma\frac{1}{1-\Gamma\sigma^{2}\chi}\int_{-\infty}^{\Delta}D\mathpzc{z}~{}(\Delta-\mathpzc{z}),$
$\displaystyle 1$ $\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{(1-\Gamma\sigma^{2}\chi)^{2}}\int_{-\infty}^{\Delta}D\mathpzc{z}~{}(\Delta-\mathpzc{z})^{2}.$
(229)

These relations are for example given in Eqs. (9-11) of [3]. Together with

$\Delta=\frac{1+\mu M^{*}}{\sqrt{q}\sigma}$ (230)

this is a set of four equations for the four unknown quantities
$q,\chi,M^{*}$, and $\Delta$.

This can be reduced to a system of three equations for the three quantities
$q,M^{*},\chi$ by using Eq. (230) in Eqs. (10.3), but it will turn out
convenient for us not to eliminate $\Delta$ at this stage. This system is to
be solved as a function of the model parameters $\mu,\sigma^{2}$, and
$\Gamma$.

Along the way, we have made the assumption $1-\Gamma\sigma^{2}\chi>0$. This
can be checked retrospectively from the numerical solution. It is also
required self-consistently in the second relation in Eq. (10.3), as $M^{*}$ is
an average of non-negative abundances and therefore must be non-negative
itself, $M^{*}\geq 0$. (The integral in the second relation is clearly
positive, as the integrand is positive in the integration range.) We also note
that the first relation in Eq. (10.3) then implies $\chi>0$. This makes sense
intuitively, as the perturbation field in Eq. (202) has a positive effect on
the growth rate of $x$.

For future reference we also introduce

$\phi=\int_{-\Delta}^{\infty}Dz=\int_{-\infty}^{\Delta}D\mathpzc{z}.$ (231)

This means that $\phi$ is the probability that the standard Gaussian variable
$z$ takes values above $-\Delta$. This is the range in which $x(z)$ in Eq.
(212) is positive. Therefore $\phi$ is the fraction of surviving species in
the stable phase of the Lotka–Volterra system.

For future use, we further define

$w_{\ell}(\Delta)\equiv\int_{-\infty}^{\Delta}D\mathpzc{z}~{}(\Delta-\mathpzc{z})^{\ell},$
(232)

for $\ell=0,1,2$. Using the definition of the error function,

${\rm erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}dt\,e^{-t^{2}}$ (233)

and recalling
$D\mathpzc{z}=\frac{1}{\sqrt{2\pi}}e^{-\mathpzc{z}^{2}/2}d\mathpzc{z}$ we have
the following explicit expressions for the $w_{\ell}(\Delta)$,

$\displaystyle w_{0}(\Delta)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1+{\rm
erf}\left(\frac{\Delta}{\sqrt{2}}\right)\right],$ $\displaystyle
w_{1}(\Delta)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[e^{-\Delta^{2}/2}\sqrt{\frac{2}{\pi}}+\Delta\left(1+{\rm
erf}\left(\frac{\Delta}{\sqrt{2}}\right)\right)\right],$ $\displaystyle
w_{2}(\Delta)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+\Delta^{2}\right)\left[1+{\rm
erf}\left(\frac{\Delta}{\sqrt{2}}\right)\right]+\frac{1}{\sqrt{2\pi}}e^{-\Delta^{2}/2}\Delta.$
(234)

We also have

$\displaystyle w_{2}(\Delta)$ $\displaystyle=$ $\displaystyle
w_{0}(\Delta)+\Delta w_{1}(\Delta).$ (235)

###  10.4 Solution procedure

We now proceed to solve the system of Eqs. (10.3) in combination with Eq.
(230).

In principle one can do this as follows: Substitute Eq. (230) into Eqs.
(10.3), and then numerically solve the resulting system of three equations for
$q,\chi$ and $M^{*}$ for given $\sigma^{2},\mu$ and $\Gamma$. This is
possible, for example using a Newton–Raphson scheme. To do this one can use
the explicit forms of the integrals in Eqs. (10.3) given in Eqs. (10.3).
Nonetheless, the resulting Newton–Raphson problem shows some degree of
dependence on initial guesses for $q,M^{*}$ and $\chi$, and often turns out to
be rather unstable. So I do not recommend this ‘brute-force’ route.

Luckily, there is a much more efficient way of doing this – one can find a
parametric solution in closed form. The idea is that we fix $\mu$ and
$\Gamma$, but not $\sigma^{2}$. Then we use $\Delta$ as a parameter. This
means that we vary $\Delta$, and then obtain $\sigma^{2},q,\chi$ and $M$ as
functions of $\Delta$. In practice this produces a table of the following
type:

$\Delta^{(1)}$ | $(\sigma^{2})^{(1)}$ | $(M^{*})^{(1)}$ | $q^{(1)}$ | $\chi^{(1)}$
---|---|---|---|---
$\Delta^{(2)}$ | $(\sigma^{2})^{(2)}$ | $(M^{*})^{(2)}$ | $q^{(2)}$ | $\chi^{(2)}$
$\Delta^{(3)}$ | $(\sigma^{2})^{(3)}$ | $(M^{*})^{(3)}$ | $q^{(3)}$ | $\chi^{(3)}$
$\vdots$ | $\vdots$ | $\vdots$ | ⋮ | $\vdots$

Here, $\Delta^{(k)}$ is the $k$-th value of our parameter $\Delta$ that we use
in our scheme, and $(\sigma^{2})^{(k)}$, $(M^{*})^{(k)}$, $q^{(k)}$ and
$\chi^{(k)}$ are the corresponding values for $\sigma^{2},M^{*},q$ and $\chi$.
Once this table has been obtained, one can the ignore the first column (values
of $\Delta$) and plot, for example, $q$ as a function of $\sigma^{2}$.

To achieve the parametric solution we first square both sides of the first
relation in Eq. (10.3), and then divide by the expressions in the third
equation. We find

$\sigma^{2}\chi^{2}=\frac{w_{0}(\Delta)^{2}}{w_{2}(\Delta)}.$ (236)

From the first relation in Eq. (10.3) on the other hand we have

$\chi-\Gamma\sigma^{2}\chi^{2}=w_{0}(\Delta).$ (237)

Putting this together we find

$\chi=w_{0}(\Delta)+\Gamma\frac{w_{0}(\Delta)^{2}}{w_{2}(\Delta)}.$ (238)

For given $\mu$ and $\Gamma$ we have therefore found an explicit expression
for $\chi$ as a function of $\Delta$. We now use this in Eq. (236) and arrive
at

$\displaystyle\sigma^{2}$ $\displaystyle=$
$\displaystyle\frac{w_{0}(\Delta)^{2}}{w_{2}(\Delta)\chi^{2}}$ (239)
$\displaystyle=$ $\displaystyle\frac{w_{2}(\Delta)}{[w_{2}(\Delta)+\Gamma
w_{0}(\Delta)]^{2}}.$

For given $\mu$ and $\Gamma$, this is an explicit expression for $\sigma^{2}$
as a function of $\Delta$.

It remains to find $q$ and $M^{*}$ as functions of $\Delta$. To do this, it is
useful to introduce the following shorthand

$b\equiv 1-\Gamma\sigma^{2}\chi.$ (240)

Using Eqs. (238) and (239) we then have

$\displaystyle b$ $\displaystyle=$ $\displaystyle
1-\Gamma\frac{w_{2}}{(w_{2}+\Gamma
w_{0})^{2}}\left(w_{0}+\Gamma\frac{w_{0}^{2}}{w_{2}}\right)$ (241)
$\displaystyle=$ $\displaystyle 1-\Gamma\frac{w_{0}}{w_{2}+\Gamma w_{0}}$
$\displaystyle=$ $\displaystyle\frac{w_{2}}{w_{2}+\Gamma w_{0}},$

where we don’t write out the argument $\Delta$ of the $w_{\ell}$. For fixed
$\Gamma$, this is an explicit expression for $b$ as a function of $\Delta$.

From the second relation in Eq. (10.3) we then have

$\sigma\sqrt{q}=\frac{bM^{*}}{w_{1}},$ (242)

which, combined with $\Delta=\frac{1+\mu M}{\sigma\sqrt{q}}$, gives

$\displaystyle\frac{1}{M^{*}}$ $\displaystyle=$
$\displaystyle\frac{\Delta}{w_{1}}b-\mu$ (243) $\displaystyle=$
$\displaystyle\frac{\Delta}{w_{1}}\frac{w_{2}}{w_{2}+\Gamma w_{0}}-\mu.$

For given $\Gamma$ and $\mu$, this is an explicit expression for $1/M^{*}$ as
a function of $\Delta$ (again we have not written out the argument $\Delta$ of
the $w_{\ell}$ on the right-hand side).

Finally, we have

$q=\left(\frac{bM^{*}}{\sigma w_{1}}\right)^{2}.$ (244)

Given $w_{1}=w_{1}(\Delta)$ and that we know $b,\sigma^{2}$ and $M^{*}$ as
functions of $\Delta$, we therefore have $q$ as a function of $\Delta$ (for
fixed $\Gamma$ and $\mu$).

Summary:
For given $\Gamma$ and $\mu$, Eqs. (10.3,230) can be solved parametrically as
follows: Vary $\Delta$ (in practice for example in range from $-10$ to $10$).
For each value of $\Delta$ calculate $w_{0}(\Delta),w_{1}(\Delta)$ and
$w_{2}(\Delta)$ from Eq. (10.3). Then:

$\displaystyle\chi$ $\displaystyle=$ $\displaystyle
w_{0}+\Gamma\frac{w_{0}^{2}}{w_{2}},$ $\displaystyle\sigma^{2}$
$\displaystyle=$ $\displaystyle\frac{w_{2}}{(w_{2}+\Gamma w_{0})^{2}},$
$\displaystyle\frac{1}{M^{*}}$ $\displaystyle=$
$\displaystyle\frac{\Delta}{w_{1}}\frac{w_{2}}{w_{2}+\Gamma w_{0}}-\mu,$
$\displaystyle q$ $\displaystyle=$
$\displaystyle\left(\frac{w_{2}}{w_{2}+\Gamma w_{0}}\frac{M^{*}}{\sigma
w_{1}}\right)^{2},$ $\displaystyle\phi$ $\displaystyle=$ $\displaystyle
w_{0}.$ (245)

For fixed $\Gamma$ and $\mu$, this gives $M^{*},q,\phi$ and $\chi$ as a
function of $\sigma^{2}$ in parametric form.

###  10.5 Interpretation of the order parameters

We now briefly turn to an interpretation of these order parameters in the
fixed-point phase in terms of the original Lotka–Volterra system. The notation
$x_{i}^{*}$ indicates abundances at stable fixed points.

  * •

The quantity $\phi=w_{0}(\Delta)$ is the fraction of surviving species, i.e.,

$\phi=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\overline{\Theta(x_{i}^{*})},$ (246)

where $\Theta(\cdot)$ is the Heaviside function, $\Theta(x)=1$ for $x>0$, and
$\Theta(x)=0$ for $x\leq 0$. In the original Lotka–Volterra equations (1)
extinctions can never occur at finite times. Instead, an abundance $x_{i}$ can
only go to zero asymptotically at infinite times if it is initially positive.
Obviously we cannot integrate the Lotka–Volterra system to infinite times
numerically. To measure $\phi$ from numerical integration (in a system with
finite $N$ obviously), one therefore has to apply a threshold, $\vartheta$, at
the finite end-time, $t_{f}$, of the integration, and identify species $i$ as
surviving if $x_{i}(t_{f})\geq\vartheta$.

  * •

The quantity $M^{*}$ is the mean abundance per species,

$M^{*}=\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}\overline{x_{i}^{*}},$ (247)

i.e., the first moment of abundances. (The sum includes the species that have
died out.).

  * •

The quantity $q$ is the second moment of the distribution of abundances,

$q=\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}\overline{(x_{i}^{*})^{2}},$
(248)

where the sum again includes species that have died out.

  * •

The quantity $\chi$ is a susceptibility, and measures how strongly species
abundances change if an external perturbation is applied. It is not easy to
measure this directly in simulations (but not impossible either).

Remark:
We re-iterate that Eqs. (245) were derived assuming that (i) the
Lotka–Volterra dynamics converges to a fixed point, and (ii) that this fixed
point is, for a given realisation of the disorder, independent of initial
conditions (we assumed the absence of memory effects). As we will see below,
this is not the case for all choices of the model parameters $\mu,\sigma^{2}$
and $\Gamma$. The procedure in Eqs. (245) only describes the outcome of the
system in the phase of unique stable fixed points.

###  10.6 Species abundance distribution

Figure 5: Species abundance distribution $P(x)$ [Eq. (250)]. Solid black curve
shows $p_{\rm surv}(x)$ in for $x>0$ [Eq. (251)]. The dashed line shows the
extension of this curve to negative $x$. The orange histogram is from
numerical integration of the Lotka–Volterra system. Here all abundances are
non-negative by construction. The peak at $x=0$ are extinct species.
Parameters are $\mu=-1$, $\Gamma=-0.5$ and $\sigma^{2}=1$. Numerical data is
from $400$ realisations of the disorder, and $N=200$. The numerical
integration is performed using an Euler scheme with adaptive time step.

Suppose now are at a point $\mu,\sigma^{2},\Gamma$ in the phase with unique
fixed points. We can then determine the order parameters $q,M^{*}$ and $\chi$
from Eqs. (245). The species abundance distribution (distribution of the
$x_{i}^{*}$’s) can then be deduced from Eq. (212), which we repeat here,

$x(z)=\mbox{max}\left(0,\frac{1+\mu M^{*}+\sqrt{q}\sigma
z}{1-\Gamma\sigma^{2}\chi}\right),$ (249)

and where we recall that $z$ is a static Gaussian random variable with mean
zero and variance one.

Thus the abundance $x(z)$ is zero whenever the second argument inside the
maximum is zero or negative. Consequently, the distribution of abundances has
a delta peak at zero, plus a part which describes the surviving species,

$P(x)=(1-\phi)\delta(x)+p_{\rm surv}(x)$ (250)

The weight of the delta-peak at $x=0$ is $1-\phi$, and we have
$\int_{0}^{\infty}dx~{}p_{\rm surv}(x)=\phi$.

The functional form of $p_{\rm surv}$ is a clipped Gaussian. More
specifically,

$p_{\rm
surv}(x)=\frac{1}{\sqrt{2\pi\Sigma^{2}}}\exp\left(-\frac{(x-x_{S})^{2}}{2\Sigma^{2}}\right)\Theta(x),$
(251)

with the Heaviside function $\Theta(x)$, and

$\displaystyle x_{S}$ $\displaystyle=$ $\displaystyle\frac{1+\mu
M^{*}}{1-\Gamma\chi\sigma^{2}},$ $\displaystyle\Sigma^{2}$ $\displaystyle=$
$\displaystyle\frac{q\sigma^{2}}{(1-\Gamma\sigma^{2}\chi)^{2}}.$ (252)

These can be read off directly from Eq. (249), as the mean and variance of the
object $\frac{1+\mu M^{*}+\sqrt{q}\sigma z}{1-\Gamma\sigma^{2}\chi}$, keeping
in mind that $z$ is a standard Gaussian random variable.

The quantities $x_{S}$ and $\Sigma^{2}$ are the mean and variance respectively
of an unclipped Gaussian, and they are not to be confused with the moments of
the clipped Gaussian $p_{\rm surv}(\cdot)$. That is to say in general we have
$M=\int_{0}^{\infty}dx\,p_{\rm surv}(x)x\neq x_{S}$ and
$q=\int_{0}^{\infty}dx\,p_{\rm surv}(x)x^{2}\neq\Sigma^{2}$. Using the
definitions of $x_{S}$ and $\Sigma^{2}$, and the relation $\phi=w_{0}(\Delta)$
with $\Delta=(1+\mu M^{*})/(\sigma\sqrt{q})$ one can directly check that
$\int_{0}^{\infty}dx~{}p_{\rm surv}(x)=\phi$.

An example is shown in Fig. 5. This is for specific parameters $\mu=-1$,
$\Gamma=-0.5$ and $\sigma^{2}=1$. The solid curve shows $p_{\rm surv}(x)$ for
$x>0$. The figure also shows $p_{\rm surv}(x)$ for $x<0$ (dashed black line).
This is the part of the Gaussian that gets clipped, i.e., the fraction of
extinct species ($x^{*}=0$) is equal to they shaded grey area under $p_{\rm
surv}$ at negative $x$. The orange histogram is from numerical integration of
the Lotka–Volterra system. More precisely the data is from $400$ independent
realisations of the interaction matrix, and for a system size of $N=200$
species. The numerical integration is performed using an Euler scheme with
adaptive time step.

###  10.7 Test against numerical integration of the Lotka–Volterra dynamics

We test the predictions of the theory against results from a numerical
integration of the Lotka–Volterra equations in Fig. 6. We show both the mean
abundance $M^{*}$ and the fraction of surviving species $\phi$. Results are
only shown in the phase of unique stable fixed points, as the theory is
otherwise not applicable.

Figure 6: Static order parameters in the unique fixed point phase. We show the
mean abundance $M^{*}$ and the fraction of surviving species $\phi$ as a
function of the variance of interactions $\sigma^{2}$. The parameter $\mu$ is
set to $\mu=-1$, the value of $\Gamma$ is as indicated in the legend. Solid
lines are from the theory in the unique fixed-point phase [Eqs. (245)]. The
markers are from a numerical integration of the Lotka–Volterra system
($N=200$, averaged over $100$ realisations of the disorder). In the numerical
integration surviving species are identified as those with $x_{i}>0.01$ at the
end of the integration.

##  11 Stability analysis and phase diagram

We will now ask when the validity of the self-consistent equations (227), or
equivalently, Eqs. (10.3) breaks down. In principle this could happen in
several different ways. For example the equations could simply not have any
physically meaningful solutions at all for some combinations of
$\mu,\sigma^{2}$ and $\Gamma$. Or the Lotka–Volterra system could show memory,
i.e., there could be multiple physically meaningful solutions for
$C(t,t^{\prime}),G(t,t^{\prime}),M(t)$, depending on initial conditions for
example.

As we will see, the relevant ways in which the assumptions of unique fixed
points break down in the Lotka–Volterra system are: (i) an onset of divergent
abundances, and (ii) a linear instability. These have been identified for
example in [1, 2, 3, 16]. We now discuss the two types of transition in turn.

###  11.1 Onset of divergent solutions

The onset of divergent abundances is signalled by a divergent mean abundance
$M^{*}$. That is, if we start in the stable phase and then modify parameters
$\mu,\sigma^{2}$ and $\Gamma$, we may approach combinations of these
parameters at which $M^{*}$ diverges.

Setting $1/M^{*}=0$ in the third relation in Eqs. (245) we find that this
occurs when

$\mu=\frac{\Delta}{w_{1}(\Delta)}\frac{w_{2}(\Delta)}{w_{2}(\Delta)+\Gamma
w_{0}(\Delta)}.$ (253)

We also have from the second relation in Eqs. (245),

$\sigma^{2}=\frac{w_{2}(\Delta)}{[w_{2}(\Delta)+\Gamma w_{0}(\Delta)]^{2}}.$
(254)

For a given choice of $\Gamma$, Eqs. (253,254) provide a parametric
description of a line in the $(\mu,\sigma^{2})$-plane (the parameter being
$\Delta$), telling us where the onset of divergence of the mean divergence
$M^{*}$ occurs. At (and beyond) this line a non-zero fraction of species
abundances diverge, and the fixed-point assumptions no longer hold.

###  11.2 Linear stability analysis

The Lotka–Volterra system also shows a linear instability. This was first
identified in a similar replicator system in [13], and we proceed along the
lines of this reference. The analysis is divided into six steps.

####  11.2.1 Step 1: Add small perturbation to effective single-species
dynamics

Assume we choose parameters such that the system is in the phase of unique
stable fixed points. We then add white noise $\varepsilon\xi_{i}(t)$ to the
Lotka–Volterra equations, i.e., we have

$\dot{x}_{i}=x_{i}\left(1-x_{i}+\sum_{j}\alpha_{ij}x_{j}+\varepsilon\xi_{i}(t)\right).$
(255)

We assume $\left\langle{\xi_{i}(t)}\right\rangle=0$, and
$\left\langle{\xi_{i}(t)\xi_{j}(t^{\prime})}\right\rangle=\delta_{ij}\delta(t-t^{\prime})$.

At the level of the effective process, this corresponds to adding white noise
$\xi(t)$ of zero mean and unit variance
[$\left\langle{\xi(t}\right\rangle=0,\left\langle{\xi(t)\xi(t^{\prime})}\right\rangle=\delta(t-t^{\prime})$],

$\displaystyle\dot{x}(t)=x(t)\left[1-x(t)+\Gamma\sigma^{2}\int
dt^{\prime}G(t,t^{\prime})x(t^{\prime})+\mu
M(t)+\eta(t)+\varepsilon\xi(t)\right].$ (256)

We now study fluctuations $y(t)$ about a fixed point of Eq. (175).
Specifically, we assume that we run the unperturbed system ($\varepsilon=0$)
to a fixed point and then imagine we ‘switch’ the perturbation $\xi(t)$ on
after the fixed point has been reached. As a consequence the abundance $x$
will be perturbed about the fixed point. We write

$x(t)=x^{*}+\varepsilon y(t).$ (257)

In the fixed-point phase and when there are no perturbations, the noise
$\eta(t)$ in the effective process becomes a static random variable at large
times. This is no longer the case in the perturbed system, and we denote the
resulting additional term in the self-consistent noise by $\varepsilon v(t)$,
i.e.,

$\eta(t)=\eta^{*}+\varepsilon v(t).$ (258)

We then expand both sides of Eq. (256). The leading-order term
($\varepsilon=0$) is simply the fixed-point relation [Eq. (210)],

$\displaystyle 0=x^{*}\left[1-x^{*}+\Gamma\sigma^{2}\chi x^{*}+\mu
M^{*}+\eta^{*}\right].$ (259)

Comparing terms of linear order in $\varepsilon$ we find,

$\displaystyle\varepsilon\dot{y}(t)$ $\displaystyle=$
$\displaystyle\varepsilon y(t)\left[1-x^{*}+\Gamma\sigma^{2}\chi x^{*}+\mu
M^{*}+\eta^{*}\right]$ (260) $\displaystyle+x^{*}\left[-\varepsilon
y(t)+\Gamma\sigma^{2}\int dt^{\prime}\\{G(t,t^{\prime})\varepsilon
y(t^{\prime})\\}+\varepsilon v(t)+\varepsilon\xi(t)\right].$

####  11.2.2 Step 2: Independence of random variables

In our analysis, we assume that

  1. 1.

The random variables $\xi(t)$ are independent of the static random variable
$x^{*}$ (equivalently, $\xi(t)$ is independent of $\eta^{*}$).

  2. 2.

The random variables $v(t)$ are independent of the static random variable
$x^{*}$.

  3. 3.

The random variables $\underline{v}$ and $\underline{\xi}$ are independent.

The first statement is clear, as $x^{*}$ is the fixed point of the effective
process in absence of the perturbation $\xi(t)$, which is only added after the
fixed point has been reached.

The assumption that $v(t)$ and $x^{*}$ are independent from one another is a
little more intricate to justify. It helps to remember that the noise
$\underline{\eta}=\eta^{*}+\varepsilon\underline{v}$ in the effective process
reflects the disorder in the original Lotka–Volterra problem, and ultimately
results from the term $\sum_{j\neq i}\alpha_{ij}x_{j}$ in the equation for
$\dot{x}_{i}$. This is also transparent from the derivation of the effective
process via the cavity method in Sec. 12. As such, the noise variable $v(t)$
results from the fluctuations of $\sum_{j\neq i}\alpha_{ij}x_{j}(t)$ in the
dynamics for $x_{i}$ away from the fixed point value $\sum_{j\neq
i}\alpha_{ij}x_{j}^{*}$. Within a linear approximation it is then sensible to
assert that $v(t)$ is independent of $x^{*}$.

The argument for the third item is similar. The term $\underline{v}$ results
from the disordered interactions of one given species with all other species.
The noise $\underline{\xi}$ acts directly on that focal species itself and,
within the linear approximation, has nothing to do with the other species.
Therefore $\underline{v}$ is independent of $\underline{\xi}$.

We now proceed to linearise Eq. (256) in $\varepsilon$, i.e. in $y(t),v(t)$
and $\xi(t)$. We first do this for effective species with zero abundance (step
3), and then subsequently for effective species with a positive abundance
(step 4).

####  11.2.3 Step 3: Further analysis for species with vanishing fixed-point
abundance

Linearising Eq. (256) about a fixed point $x^{*}=0$ one has from Eq. (260) and
writing again $\eta^{*}=\sqrt{q}\sigma z$,

$\frac{dy(t)}{dt}=y(t)\left[1+\mu M^{*}+\sqrt{q}\sigma z\right].$ (261)

Within our ansatz, the object in the square brackets is negative for fixed
points at zero [see Eq. (212), and recall that $1-\Gamma\sigma^{2}\chi>0$]. We
conclude that $y(t)\to 0$ for $t\to\infty$. We note explicitly that the
presence of the noise $\varepsilon\xi(t)$ makes no difference for $\dot{y}$ to
linear order in $\varepsilon$, i.e., Eq. (261) is valid (to linear order in
$\varepsilon$) for the model with noise $\varepsilon\xi(t)$ in Eq. (256). This
implies the following: Imagine we run the unperturbed Lotka–Volterra system
for a particular realisation of the disorder, and a particular species $i$
goes extinct at fixed point of this system. If we then switch on the noise and
perturb the fixed point, then species $i$ will still remain extinct (provided
the noise strength is small enough for the linear approximation to work).

In essence the result in this step means, that once a species has gone
extinct, the system is stable against re-invasion of that species.

Remark:
We also note at this point that Eq. (261) indicates that a zero fixed point is
not stable if the object in the square brackets on the right-hand side is
positive. This retrospectively justifies the use the non-zero solution $x^{*}$
in Eq. (212) when $1+\mu M^{*}+\sqrt{q}\sigma z>0$.

####  11.2.4 Step 4: Further analysis for non-zero fixed-point abundances

For non-zero $x^{*}$ we have $1-x^{*}+\Gamma\sigma^{2}\chi x^{*}+\mu
M^{*}+\eta^{*}=0$ [see Eq. (210)], and hence Eq. (260) turns into

$\displaystyle\frac{dy(t)}{dt}$ $\displaystyle=$ $\displaystyle
x^{*}\left[-y(t)+\Gamma\sigma^{2}\int_{0}^{t}G(t,t^{\prime})y(t^{\prime})dt^{\prime}+v(t)+\xi(t)\right].$
(262)

Recalling $G(t,t^{\prime})=G(t-t^{\prime})$ in the fixed-point phase, we
transform into Fourier space (Fourier transforms are indicated by tildes), and
obtain

$\displaystyle\frac{i\omega\widetilde{y}(\omega)}{x^{*}}$ $\displaystyle=$
$\displaystyle\left(\Gamma\sigma^{2}\widetilde{G}(\omega)-1\right)\widetilde{y}(\omega)+\widetilde{v}(\omega)+\widetilde{\xi}(\omega).$
(263)

This leads to

$\widetilde{y}(\omega)=\frac{\widetilde{v}(\omega)+\widetilde{\xi}(\omega)}{\frac{i\omega}{x^{*}}+\left(1-\Gamma\sigma^{2}\widetilde{G}(\omega)\right)},$
(264)

and we remember that this is valid only for extant species (species that have
not gone extinct at the fixed point of the unperturbed system).

The noise $\underline{\xi}$ has mean zero, and within the linear approximation
in Eq. (262) it is then reasonable to conclude that
$\left\langle{y(t)}\right\rangle_{y}=0$ and
$\left\langle{v(t)}\right\rangle_{v}=0$, for any fixed $x^{*}$.

####  11.2.5 Step 5: Self-consistency relations

The self-consistency relation
$\left\langle{\left\langle{\eta(t)\eta(t^{\prime})}\right\rangle}\right\rangle_{*,v}=\sigma^{2}\left\langle{\left\langle{x(t)x(t^{\prime})}\right\rangle}\right\rangle_{*,y}$
translates into

$\displaystyle\left\langle{(\eta^{*})^{2}}\right\rangle_{*}+\varepsilon\left\langle{\left\langle{\eta^{*}v(t)}\right\rangle}\right\rangle_{*,v}+\varepsilon\left\langle{\left\langle{\eta^{*}v(t^{\prime})}\right\rangle}\right\rangle_{*,v}+\varepsilon^{2}\left\langle{v(t)v(t^{\prime})}\right\rangle_{*,v}$
(265) $\displaystyle=$
$\displaystyle\sigma^{2}\left(\left\langle{(x^{*})^{2}}\right\rangle_{*}+\varepsilon\left\langle{\left\langle{x^{*}y(t)}\right\rangle}\right\rangle_{*,y}+\varepsilon\left\langle{\left\langle{x^{*}y(t^{\prime})}\right\rangle}\right\rangle_{*,y}+\varepsilon^{2}\left\langle{y(t)y(t^{\prime})}\right\rangle_{*,y}\right)$

We have used the notation $\left\langle{\cdots}\right\rangle_{*}$ for averages
over $\eta^{*}$ or $x^{*}$, and
$\left\langle{\left\langle{\cdots}\right\rangle}\right\rangle_{*,y}$ stands
for a combined average over $x^{*}$ (or $\eta^{*}$) and $y$, and similarly for
$\left\langle{\left\langle{\cdots}\right\rangle}\right\rangle_{*,v}$.

We know that
$\left\langle{(\eta^{*})^{2}}\right\rangle_{*}=\sigma^{2}\left\langle{(x^{*})^{2}}\right\rangle_{*}$
(self-consistency relation in the earlier fixed-point analysis). The linear
terms in $\varepsilon$ in Eq. (265) vanish because $v(t)$ and $y(t)$ each have
zero average for any fixed $x^{*}$. Thus we are left with the requirement

$\left\langle{v(t)v(t^{\prime})}\right\rangle_{*,v}=\sigma^{2}\left\langle{y(t)y(t^{\prime})}\right\rangle_{*,y}.$
(266)

Given the independence of $\underline{v}$ from $x^{*}$ we can drop the average
over $x^{*}$ on the left. We therefore have

$\left\langle{v(t)v(t^{\prime})}\right\rangle_{v}=\sigma^{2}\left\langle{y(t)y(t^{\prime})}\right\rangle_{*,y}.$
(267)

####  11.2.6 Step 6: Criterion for linear instability

We now proceed to calculate
$\left\langle{|\widetilde{y}(\omega)|^{2}}\right\rangle_{*,y}$. This is the
Fourier transform of the correlation function
$\left\langle{y(t+\tau)y(\tau)}\right\rangle_{*,y}$ in the stationary state
($t\to\infty$), see also Sec. 3.3. We recall that $\underline{y}$ is a
perturbation about $x^{*}$, and that the statistics of $y(t)$ (equivalently,
those of $\widetilde{y}(\omega)$) depend on the fixed-point value $x^{*}$.
Specifically $\lim_{t\to\infty}y(t)=0$, when $x^{*}=0$, whereas
$\widetilde{y}(\omega)$ is given by Eq. (264) when $x^{*}>0$. In other words

$|\widetilde{y}(\omega)|^{2}=\left\\{\begin{array}[]{cl}\left|\frac{\widetilde{v}(\omega)+\widetilde{\xi}(\omega)}{\frac{i\omega}{x^{*}}+\left(1-\Gamma\sigma^{2}\widetilde{G}(\omega)\right)}\right|^{2}&\mbox{for
surviving species},\\\ {}\hfil&{}\\\ 0&\mbox{for extinct
species}.\end{array}\right.$ (268)

We now carry out an average on both sides. This involves a three-fold
averaging procedure, over fixed-point abundances $x^{*}$, realisations of the
noise $\underline{v}$ and over realisations of $\underline{\xi}$. Using the
independence properties from Step 2, we have

$\displaystyle\left\langle{|\widetilde{y}(\omega)|^{2}}\right\rangle_{*,y}$
$\displaystyle=$ $\displaystyle\left(\int dx^{*}\,p_{\rm
surv}(x^{*})\left|\frac{i\omega}{x^{*}}+\left(1-\Gamma\sigma^{2}\widetilde{G}(\omega)\right)\right|^{-2}\right)$
(269)
$\displaystyle\times\left(\left\langle{|\widetilde{v}(\omega)|^{2}}\right\rangle_{v}+\left\langle{|\widetilde{\xi}(\omega)|^{2}}\right\rangle_{\xi}\right).$

Next we use the self-consistency relation in Eq. (267), which translates into

$\left\langle{|\widetilde{v}(\omega)^{2}|}\right\rangle_{v}=\sigma^{2}\left\langle{|\widetilde{y}(\omega)|^{2}}\right\rangle_{*,y}.$
(270)

Thus, we find

$\displaystyle\left\langle{|\widetilde{y}(\omega)|^{2}}\right\rangle_{*,y}$
$\displaystyle=$ $\displaystyle\left(\int dx^{*}\,p_{\rm
surv}(x^{*})\left|\frac{i\omega}{x^{*}}+\left(1-\Gamma\sigma^{2}\widetilde{G}(\omega)\right)\right|^{-2}\right)$
(271)
$\displaystyle\times\left(\sigma^{2}\left\langle{|\widetilde{y}(\omega)|^{2}}\right\rangle_{*,y}+\left\langle{|\widetilde{\xi}(\omega)|^{2}}\right\rangle_{\xi}\right).$

Following [13] we focus on $\omega=0$. We have $\widetilde{G}(0)=\chi$, and
the term inside the integral over $x^{*}$ becomes independent of $x^{*}$.
Noting that $\int_{0}^{\infty}dx^{*}p_{\rm surv}(x^{*})=\phi$, we then arrive
at

$\left\langle{|\widetilde{y}(0)|^{2}}\right\rangle_{*,y}=\phi\left(\Gamma\sigma^{2}\chi-1\right)^{-2}\left[\sigma^{2}\left\langle{|\widetilde{y}(0)|^{2}}\right\rangle_{*,y}+1\right].$
(272)

Eq. (272) can be re-written as

$\left\langle{|\widetilde{y}(0)|^{2}}\right\rangle_{*,y}=\frac{\phi}{\left(\Gamma\sigma^{2}\chi-1\right)^{2}-\phi\sigma^{2}}.$
(273)

This indicates that $\left\langle{|\widetilde{y}(0)|^{2}}\right\rangle_{*,y}$
diverges when $\phi\sigma^{2}=(1-\Gamma\sigma^{2}\chi)^{2}$. When
$\phi\sigma^{2}>(1-\Gamma\sigma^{2}\chi)^{2}$, the right-hand side is
negative, leading to a contradiction (and therefore a breakdown of the fixed-
point ansatz). One has $\phi\sigma^{2}<(1-\Gamma\sigma^{2}\chi)^{2}$ in the
stable phase, consistent with a well-defined (positive) quantity
$\left\langle{|\widetilde{y}(0)|^{2}}\right\rangle$.

The onset of linear instability is thus signalled by the condition
$\phi\sigma^{2}=(1-\Gamma\sigma^{2}\chi)^{2}$. This in turn leads to
$\Delta=0$ in Eqs. (10.3). To see this we substitute
$\phi\sigma^{2}=(1-\Gamma\sigma^{2}\chi)^{2}$ into
$1=\frac{\sigma^{2}}{(1-\Gamma\sigma^{2}\chi)^{2}}\int_{-\infty}^{\Delta}D\mathpzc{z}(\Delta-\mathpzc{z})^{2}$
[the third relation in Eq. (10.3)], and find

$\phi=\int_{-\infty}^{\Delta}D\mathpzc{z}~{}(\Delta-\mathpzc{z})^{2}.$ (274)

On the other hand we also have $\phi=\int_{-\infty}^{\Delta}D\mathpzc{z}$.
Comparing the two expressions gives $\Delta=0$, and hence $\phi=1/2$.

Using this in the first and third relations in Eqs. (10.3), we have,
respectively,

$\displaystyle 2\chi$ $\displaystyle=$
$\displaystyle\frac{1}{1-\Gamma\sigma^{2}\chi},$ $\displaystyle 2$
$\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{(1-\Gamma\sigma^{2}\chi)^{2}}.$ (275)

From this we find $\chi^{2}=1/(2\sigma^{2})$. Using $\chi>0$ (which we know
from Sec. 10.3) we have $\chi=1/(\sqrt{2}\sigma)$. Substituting this in the
first relation in Eq. (11.2.6) in turn leads to

$\sigma_{c}^{2}(\Gamma)=\frac{2}{(1+\Gamma)^{2}}.$ (276)

This is the linear instability condition reported for example in [1, 2, 3].

###  11.3 Stability diagram

The stability diagram is shown in Fig. 7 in the $(\mu,\sigma^{2})$ plane. On
the left I show the lines of divergent abundance (solid) and of the onset of
linear instability (dashed) for different values of $\Gamma$. The linear
instability sets in at $\sigma^{2}=2/(1+\Gamma^{2})$ [Eq. (276)], the onset of
diverging abundance is given by Eqs. (253,254). On the right-hand panel I
indicate the typical behaviour of the system in the different phases.

Figure 7: Stability diagram for the random Lotka-Volterra system in the
$(\mu,\sigma)$ plane. Solid lines mark the onset of diverging abundances
[$M\to\infty$, Eqs. (253,254)], dashed lines are the onset of linear
instability [Eq. (276)]. The dotted lines also show the condition of diverging
abundances. However, the fixed-point equations are no longer valid in this
part of the phase diagram, as the linear instability has already set in. On
the left we show the phase diagram for different values of the symmetry
parameter $\Gamma$. On the right we fix $\Gamma=-0.5$ and indicate the typical
behaviour observed in the different parts of the phase plane (see text for
further details).

There are three different phases:

Unique stable fixed point phase. The system converges to one single fixed
point for any fixed realisation of the disorder. That is to say, once the
interaction matrix has been drawn, the Lotka–Volterra dynamics always goes to
the same stable fixed point, regardless of the initial conditions (provided
all $x_{i}$ are positive at the beginning). This is the phase described by our
fixed-point theory [for example Eqs. (10.3) and (245)]. This phase is realised
in the part of the phase diagram below the linear instability line (dashed)
and to the left of the (solid) line marking the onset of diverging abundance.

Diverging abundance. Starting in the phase with unique stable fixed points and
increasing the mean of interactions $\mu$ (keeping $\sigma^{2}$ fixed) the
system eventually crosses into the phase of diverging abundance. Within the
stable fixed point phase $M$ remains finite and there is a unique stable fixed
point. However, as one approaches the solid line in the stability diagram from
the left, the fixed point abundances gradually increase, and eventually the
fixed point ‘escapes’ to infinity, so that $M^{*}\to\infty$.

Multiple fixed points or persistent dynamics. In the part of the phase diagram
labelled ‘multiple fixed points or persistent dynamics’ there are several
possible types of behaviour. For example, it is possible that the system never
really settles down, the abundances keep fluctuating indefinitely. This is
what I mean by ’persistent dynamics’ (the dynamics persists forever). It is
also possible that the system has many stable fixed points for a given
realisation of the random interactions. Which one is reached, depends on
initial conditions. This is what we mean by ‘multiple fixed points’. I add
that this phase is sometimes also called ‘multiple attractor’ (MA) phase, see
e.g. [1, 2]. The ‘multiple fixed points or persistent dynamics’/MA phase is
usually observed in a relatively small part of the phase diagram just above
the linear instability line (dashed). If the variance of interactions
($\sigma^{2}$) is increased further, abundances typically diverge at some
point. In the phase diagram we indicate the boundary between the ‘multiple
fixed points or persistent dynamics’ phase and the phase with diverging
abundances by a dotted line. This line is obtained from Eqs. (253,254).
However, it is important to realise that the linear instability line has been
crossed at this point (the dotted line is above the dashed line in the phase
diagram). So technically, our fixed-point equations (10.3) do not hold any
longer. In this sense, the dotted line is only an approximation. Starting from
the multiple-attractor phase, the line indicates ‘roughly’ where abundances
start to diverge.

Figure 8: Sample trajectories. The figure shows single trajectories of the
system in different phases of the stability diagram. On the right: stable
phase ($\mu=-1,\Gamma=-0.5,\sigma^{2}=1$). Middle: Oscillatory behaviour
($\mu=-1,\Gamma=-0.5,\sigma^{2}=10$). Right: Irregular, potentially chaotic
dynamics ($\mu=-0.1,\Gamma=-0.5,\sigma^{2}=10$). Each plot is for a single
realisation with $N=200$ species. Grey dots in the background show all species
abundances. A selection of trajectories $x_{i}(t)$ is highlighted.

Examples of possible trajectories are shown in Fig. 8. For further examples,
see also [54].

Further remark:
The phase diagram is obtained from the analysis in the thermodynamic limit. I
stress again that not all realisations for finite $N$ show the same behaviour
in a given part of the phase diagram. For example, it is possible that
individual samples diverge in the phase labelled ‘unique fixed point’, or that
some samples converge to single stable fixed points above the line marking the
linear instability. The stability diagram describes the behaviour for
$N\to\infty$, which we would also expect to be dominant (but not fully
exclusive) for sufficiently large $N$.

It is possible to systematically classify the different behaviours of
realisations from a numerical integration of the Lotka–Volterra system. One
can then systematically identify the dominant behaviour in the different parts
of the stability diagram. Details can be found for example in [15, 18].

##  12 Cavity approach to dynamic mean-field theory

In this section we describe an alternative derivation of the effective single-
species process. This is based on the ‘cavity method’. For background material
see e.g. [5]. The method has been used for the study of Lotka–Volterra systems
for example in [1, 2, 16]. A an overview of applications of the cavity method
to models in community ecology can be found in [55].

Assume we have a solution $x_{1}(t),\dots,x_{N}(t)$ for the original
Lotka–Volterra problem

$\dot{x}_{i}=x_{i}\left(1-x_{i}+\sum_{j}\alpha_{ij}x_{j}+h_{i}(t)\right),$
(277)

for given parameters $\mu,\Gamma,\sigma^{2}$ and a fixed large, but finite
$N$. Now we introduce an additional species $i=0$. To do this we need
additional interaction coefficients $\alpha_{i,0}$ and $\alpha_{0,i}$
($i=1,\dots,N$. These are chosen with the same statistics as the $\alpha_{ij}$
($i,j=1,\dots,N$). We have included perturbation fields $h_{i}(t)$. The
reasons for this will become clear below.

The introduction of the new species will change the dynamics for
$x_{1},\dots,x_{N}$, and we have

$\frac{d}{dt}{\widetilde{x}_{i}}(t)=\widetilde{x}_{i}(t)\left(1-\widetilde{x}_{i}(t)+\sum_{j=1}^{N}\alpha_{ij}\widetilde{x}_{j}(t)+\alpha_{i0}x_{0}(t)+h_{i}(t)\right),$
(278)

We have used the notation $\widetilde{x}_{i}$ for species abundances in the
system with the additional species. (I am aware that we used tildes for
Fourier transforms in earlier sections. Unfortunately the selection of
possible diacritics that can be added to letters is finite. The risk of
confusion should be minimal though, throughout this section tildes are not
used for Fourier transforms.) We do not give $x_{0}$ a tilde, as it is clear
that $x_{0}$ only exists in the system with species $0$ added.

Given that the effect of adding one single species to a system with $N\gg 1$
existing species is small, we can work in the linear-response regime. We then
have

$\widetilde{x}_{i}(t)=x_{i}(t)+\int_{0}^{t}dt^{\prime}\frac{\delta
x_{i}(t)}{\delta h_{i}(t^{\prime})}\alpha_{i0}x_{0}(t^{\prime}).$ (279)

The time-evolution of $x_{0}$ on the other hand is given by

$\frac{d}{dt}x_{0}(t)=x_{0}(t)\left[1-x_{0}(t)+\sum_{i=1}^{N}\alpha_{0,i}\widetilde{x}_{i}(t)+h_{0}(t)\right].$
(280)

Using Eq. (279) the dynamics for $x_{0}$ can be written as

$\displaystyle\frac{d}{dt}x_{0}(t)$ $\displaystyle=$ $\displaystyle
x_{0}(t)\left[1-x_{0}(t)+\sum_{i=1}^{N}\alpha_{0,i}x_{i}(t)+h_{0}(t)\right.$
(281)
$\displaystyle+\left.\sum_{i=1}^{N}\alpha_{0,i}\alpha_{i,0}\int_{0}^{t}dt^{\prime}\frac{\delta
x_{i}(t)}{\delta h_{i}(t^{\prime})}x_{0}(t^{\prime})\right].$

We first deal with the last term on the right-hand side. This term can be
approximated as follows for large $N$,

$\sum_{i=1}^{N}\alpha_{0,i}\alpha_{i,0}\frac{\delta x_{i}(t)}{\delta
h_{i}(t^{\prime})}\approx\Gamma\sigma^{2}G(t,t^{\prime})$ (282)

This is justified because
$\overline{\alpha_{0,i}\alpha_{i,0}}=\Gamma\frac{\sigma^{2}}{N}+\left(\frac{\mu}{N}\right)^{2}$
and $G(t,t^{\prime})=\frac{1}{N}\sum_{i=1}^{N}\frac{\delta x_{i}(t)}{\delta
h_{i}(t^{\prime})}$.

Making this replacement in Eq. (281) we arrive at

$\displaystyle\frac{d}{dt}x_{0}(t)$ $\displaystyle=$ $\displaystyle
x_{0}(t)\left[1-x_{0}(t)+\Gamma\sigma^{2}\int_{0}^{t}dt^{\prime}\,G(t,t^{\prime})x_{0}(t^{\prime})+\sum_{i=1}^{N}\alpha_{0,i}x_{i}(t)\right],$
(283)

where we have set the perturbation fieldf $h_{0}(t)$ to zero (the field is no
longer required, now that we have generated the response function
$\underline{\underline{G}}$).

Next we look at the remaining term $\sum_{i}\alpha_{0,i}x_{i}(t)$ in Eq.
(283), noting that $\alpha_{0,i}$ has mean $\mu/N$ and variance
$\sigma^{2}/N$, and can therefore be written in the form
$\alpha_{0,i}=\frac{\mu}{N}+\frac{\sigma}{\sqrt{N}}u_{0,i}$, with a Gaussian
random variable $u_{0,i}$ of mean zero and variance one, and with no
correlations between $u_{0,i}$ and $u_{0,j}$ for $i\neq j$. We then have

$\sum_{i=1}^{N}\alpha_{0,i}x_{i}(t)=\frac{\mu}{N}\sum_{i}x_{i}(t)+\frac{\sigma}{\sqrt{N}}\sum_{i}u_{0,i}x_{i}(t).$
(284)

Writing $M(t)=\frac{1}{N}\sum_{i}x_{i}(t)$, we therefore find that the
Gaussian random variable $\sum_{i}\alpha_{0,i}x_{i}(t)$ has mean $\mu M(t)$.
The variance and correlations in time of
$\frac{\sigma}{\sqrt{N}}\sum_{i}u_{0,i}x_{i}(t)$ are,

$\displaystyle\frac{\sigma^{2}}{N}\overline{\left(\sum_{i}u_{0,i}x_{i}(t)\right)\left(\sum_{i}u_{0,i}x_{i}(t^{\prime})\right)}$
$\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{N}\sum_{i,j=1}^{N}\overline{u_{0,i}u_{0,j}}x_{i}(t)x_{j}(t^{\prime})$
(285) $\displaystyle=$
$\displaystyle\frac{\sigma^{2}}{N}\sum_{i=1}^{N}x_{i}(t)x_{i}(t^{\prime}),$
$\displaystyle=$ $\displaystyle\sigma^{2}C(t,t^{\prime})$

We have used the fact that $u_{0,i}$ and $u_{0,j}$ are uncorrelated for $i\neq
j$, and that the variance of each $u_{0,i}$ is $\sigma^{2}/N$, so that
$\overline{u_{0,i}u_{0,j}}=\delta_{ij}\sigma^{2}/N$. In the last step, we have
then written
$C(t,t^{\prime})=\frac{1}{N}\sum_{i=1}^{N}x_{i}(t)x_{i}(t^{\prime})$.

Putting things together, we have found that we can write

$\sum_{i}\alpha_{0,i}x_{i}(t)=\mu M(t)+\eta_{0}(t),$ (286)

where $\eta_{0}(t)$ is Gaussian noise with mean zero, and with correlations in
time given by

$\left\langle{\eta_{0}(t)\eta(t^{\prime})_{0}}\right\rangle=\sigma^{2}C(t,t^{\prime}).$
(287)

Using this, we recover the effective dynamics in Eq. (202), namely,

$\displaystyle\frac{d}{dt}x_{0}(t)$ $\displaystyle=$ $\displaystyle
x_{0}(t)\left[1-x_{0}(t)+\mu
M(t)+\Gamma\sigma^{2}\int_{0}^{t}dt^{\prime}\,G(t,t^{\prime})x_{0}(t^{\prime})+\eta_{0}(t)\right].$
(288)

This derivation [Eqs. (284) and (285) in particular] highlights the origin of
the coloured noise in the effective process. The term $\eta_{0}(t)$ in the
dynamcics for $x_{0}$ results from the disordered interaction
$\sum_{i}\alpha_{0,i}x_{i}(t)$ of species $0$ with all other species
$i=1,\dots,N$. We used this in Sec. 11.2.2 for the linear stability analysis
of fixed points of the Lotka–Volterra system.

The reader may legitimately ask why we went through dozens pages of material
to introduce the generating-functional method, when the same effective process
can be derived in about two pages using the cavity method. This is partly a
matter of taste. Personally I’d argue that the generating-functional approach,
while longer, is clearer and leaves less space for ambiguity. Once the initial
problem is defined, the sequence of steps one needs to go through is well
defined, and in a sense ‘fail safe’. It is like an locomotive that you put on
the tracks. Once you have done this, it is just a matter of having enough
‘steam’ to plough ahead. The cavity path seems a little more cumbersome, for
example the replacements after Eq. (282) and those near Eq. (285), while
sensible, are not justified to the same level as the steps in the generating-
functional calculation leading to the effective single-species process. The
generating functional method can also be used to derive further (non-Gaussian)
properties of the surviving community, and its stability, see [17].

##  13 Summary and discussion

In summary, we have presented a detailed step-by-step guide through the
generating-functional analysis for a Lotka–Volterra system with random
interactions. We first reviewed the idea of mean-field theory (Sec. 2) and
some of the central mathematical tools (Sec. 3). Then we discussed the concept
of generating functions and functionals (Sec. 4). In Sec. 5 we then looked at
a relatively simple toy model, with no disorder. Sections 6 to 9 then
contained the actual generating-functional analysis of the Lotka–Volterra
problem. We first derived the saddle-point equations, and then proceeded along
the ‘vanilla’ route, before we presented a more complete derivation of the
dynamic mean-field theory for the Lotka–Volterra system (the effective single-
species process and associated self-consistency relations). Sections 10 and 11
then focused on a fixed-point ansatz and a stability analysis of this fixed
point respectively. From this we were able to establish the stability diagram.
In Sec. 12 finally, we presented an alternative derivation of the effective
dynamics, using the so-called ‘cavity’ method.

I hope the notes are sufficiently self-contained so that beginners do not have
to compile the relevant pieces from a number of different sources. I also hope
that the length of this paper is not too daunting. I tried to be as clear as I
could about the structure of the document as a whole (how the different
sections relate to one another), and about the contents of each section. I
also tried to balance pedagogical efficiency with analytical rigour (the
rigour of a theoretical physicist, not that of a pure mathematician). One of
my main aims was to tell the reader at any point of the calculation what the
assumptions and simplifications are that lead to a certain result, what
earlier results are used to make the next step, and to make sure that it is
possible keep track the flow of the overall argument.

There is much that is not said in these notes. As mentioned in the
introduction I don’t really motivate why one would be interested in the
disordered Lotka–Volterra model (I point to the relevant literature though).
Further, I restrict the notes to the ‘standard’ generalised Lotka–Volterra
equations (1). There are many extensions one could study and that have been
studied in the literature, either with generating functionals or with other
methods. This includes more structured interaction matrices and networks [56,
19, 18, 20], intrinsic noise [24, 57], MacArthur type resource models [58, 59]
and the study of the detailed properties of the interactions between surviving
species [1, 17]. One should also note the complementary physics literature
using equilibrium methods (e.g. the replica approach) to Lotka–Volterra
systems with symmetric interactions (see e.g. [23, 60, 61]). Exciting recent
work uses dynamic mean field theory to study the unstable phase [27].

As a final note, I’d like to re-iterate that the generating-functional
approach to disordered systems has been used and continues to be used in a
range of different areas, including spin glass physics [5], random-matrix
theory [62], neural networks [10, 51], game theory and related topics [9, 52,
63]. While these notes focus on the Lotka–Volterra model, the main steps in
particular in Sections 6-9 can be transferred without much change to other
disordered dynamical system. Because of this I hope that these notes are
useful to a broader set of readers.

## Acknowledgements

I am grateful to David Sherrington and Ton Coolen who introduced me to
generating functionals in the first place (more than twenty years ago during
my PhD). I thank Andrea De Martino, Reimer Kühn and in particular Joseph W.
Baron for collaboration and discussions which sharpened my understanding of
the method. I also thank (in alphabetical order) Sebastian Castedo, Alaina
Cockerell, Mark Crumpton, Joshua Holmes, Thomas Jun Jewell, Carlotta Nunzi,
Lyle Poley, Enrique Rozas, Christopher Ryder, and Laura Sidhom for useful
comments on earlier versions of these notes.

##### Funding information

I acknowledge partial financial support from the Agencia Estatal de
Investigación and Fondo Europeo de Desarrollo Regional (FEDER, UE) under
project APASOS (PID2021-122256NB-C21, PID2021-122256NB-C22), and the María de
Maeztu programme for Units of Excellence, CEX2021-001164-M funded by
MICIU/AEI/10.13039/501100011033.

## References

  * [1] G. Bunin,  _Interaction patterns and diversity in assembled ecological communities_ ,  arXiv preprint arXiv:1607.04734 (2016),  https://arxiv.org/abs/1607.04734.
  * [2] G. Bunin,  _Ecological communities with lotka-volterra dynamics_ ,  Phys. Rev. E 95, 042414 (2017).
  * [3] T. Galla,  _Dynamically evolved community size and stability of random lotka-volterra ecosystems_ ,  EPL (Europhysics Letters) 123(4), 48004 (2018).
  * [4] J. Hofbauer and K. Sigmund,  _Evolutionary games and population dynamics_ ,  Cambridge university press (1998).
  * [5] M. Mézard, G. Parisi and M. Virasoro,  _Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications_ , vol. 9,  World Scientific Publishing Company, London (1987).
  * [6] C. De Dominicis,  _Dynamics as a substitute for replicas in systems with quenched random impurities_ ,  Physical Review B 18, 4913 (1978).
  * [7] H. Sompolinsky and A. Zippelius,  _Relaxational dynamics of the edwards-anderson model and the mean-field theory of spin-glasses_ ,  Phys. Rev. B 25(11), 6860 (1982).
  * [8] T. R. Kirkpatrick and D. Thirumalai,  _p-spin-interaction spin-glass models: Connections with the structural glass problem_ ,  Phys. Rev. B 36(10), 5388 (1987).
  * [9] A. C. C. Coolen,  _The Mathematical Theory of Minority Games: Statistical Mechanics of Interacting Agents)_ ,  Oxford University Press, Oxford, UK (2005).
  * [10] A. C. Coolen, P. Sollich and R. Kühn,  _Theory of Neural Information Processing Systems_ ,  Oxford University Press, Oxford, UK (2005).
  * [11] A. Altland and B. D. Simons,  _Condensed Matter Field Theory_ ,  Cambridge University Press, Cambridge (2010).
  * [12] J. A. Hertz, Y. Roudi and P. Sollich,  _Path integral methods for the dynamics of stochastic and disordered systems_ ,  Journal of Physics A: Mathematical and Theoretical 50(3), 033001 (2016).
  * [13] M. Opper and S. Diederich,  _Phase transition and 1/f noise in a game dynamical model_ ,  Physical Review Letters 69(10), 1616 (1992).
  * [14] T. Galla,  _Random replicators with asymmetric couplings_ ,  Journal of Physics A: Mathematical and General 39(15), 3853 (2006).
  * [15] L. Sidhom and T. Galla,  _Ecological communities from random lotka-volterra dynamics with nonlinear functional response_ ,  Physical Review E 101, 032101 (2019).
  * [16] F. Roy, G. Biroli, G. Bunin and C. Cammarota,  _Numerical implementation of dynamical mean field theory for disordered systems: application to the lotka-volterra model of ecosystems_ ,  Journal of Physics A: Mathematical and Theoretical 52, 484001 (2019).
  * [17] J. W. Baron, T. J. Jewell, C. Ryder and T. Galla,  _Breakdown of random-matrix universality in persistent lotka-volterra communities_ ,  Phys. Rev. Lett. 130, 137401 (2023).
  * [18] J. I. Park, D.-S. L. Lee, S. H. Lee and H. J. Park,  _Incorporating heterogeneous interactions for ecological biodiversity_ ,  preprint arXiv:2403.15730 https://arxiv.org/abs/2403.15730.
  * [19] L. Poley, T. Galla and J. W. Baron,  _Interaction networks in persistent lotka-volterra communities_ ,  preprint arXiv:2404.08600 https://arxiv.org/abs/2404.08600.
  * [20] F. Aguirre-López,  _Heterogeneous mean-field analysis of the generalized lotka-volterra model on a network_ ,  preprint arXiv:2404.11164 https://arxiv.org/abs/2404.11164.
  * [21] S. Azaele and A. Maritan,  _Large system population dynamics with non-gaussian interactions_ ,  preprint arXiv:2306.13449 https://arxiv.org/abs/2306.13449.
  * [22] A. C. C. Coolen, L. Holmberg and J. E. Barrett,  _Principles of Survival Analysis_ ,  https://nms.kcl.ac.uk/ton.coolen/allnotes/POSA.pdf.
  * [23] G. Biroli, G. Bunin and C. Cammarota,  _Marginally stable equilibria in critical ecosystems_ ,  New Journal of Physics 20(8), 083051 (2018).
  * [24] A. Altieri, F. Roy, C. Cammarota and G. Biroli,  _Properties of equilibria and glassy phases of the random lotka-volterra model with demographic noise_ ,  Physical Review Letters 126(25), 258301 (2021).
  * [25] M. Barbier, C. de Mazancourt, M. Loreau and G. Bunin,  _Fingerprints of high-dimensional coexistence in complex ecosystems_ ,  Phys. Rev. X 11, 011009 (2021).
  * [26] I. Akjouj, M. Barbier, M. Clenet, W. Hachem, M. Maïda, F. Massol, J. Najim and V. C. Tran,  _Complex systems in ecology: a guided tour with large lotka–volterra models and random matrices_ ,  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 480, 20230284 (2024).
  * [27] T. Arnoulx de Pirey and G. Bunin,  _Many-species ecological fluctuations as a jump process from the brink of extinction_ ,  Phys. Rev. X 14, 011037 (2024).
  * [28] R. M. May,  _Stability and complexity in model ecosystems_ ,  In _Stability and Complexity in Model Ecosystems_. Princeton university press (2019).
  * [29] S. Allesina and S. Tang,  _The stability–complexity relationship at age 40: a random matrix perspective_ ,  Population Ecology 57(1), 63 (2015).
  * [30] R. M. May,  _Will a large complex system be stable?_ ,  Nature 238(5364), 413 (1972).
  * [31] R. M. May,  _Stability in multispecies community models_ ,  Mathematical Biosciences 12(1), 59 (1971).
  * [32] S. Allesina and S. Tang,  _Stability criteria for complex ecosystems_ ,  Nature 483(7388), 205 (2012).
  * [33] D. Gravel, F. Massol and M. A. Leibold,  _Stability and complexity in model meta-ecosystems_ ,  Nature Communications 7, 12457 (2016).
  * [34] T. Namba,  _Multi-faceted approaches toward unravelling complex ecological networks_ ,  Population Ecology 57(1), 3 (2015).
  * [35] J. Grilli, T. Rogers and S. Allesina,  _Modularity and stability in ecological communities_ ,  Nature Communications 7(1), 1 (2016).
  * [36] T. Gibbs, J. Grilli, T. Rogers and S. Allesina,  _Effect of population abundances on the stability of large random ecosystems_ ,  Physical Review E 98, 022410 (2018).
  * [37] L. Stone,  _The feasibility and stability of large complex biological networks: a random matrix approach_ ,  Scientific Reports 8(1), 1 (2018).
  * [38] Y. V. Fyodorov and B. A. Khoruzhenko,  _Nonlinear analogue of the may-wigner instability transition_ ,  Proceedings of the National Academy of Sciences 113(25), 6827 (2016).
  * [39] G. Barabás, M. J. Michalska-Smith and S. Allesina,  _Self-regulation and the stability of large ecological networks_ ,  Nature Ecology & Evolution 1(12), 1870 (2017).
  * [40] T. Gross, L. Rudolf, S. A. Levin and U. Dieckmann,  _Generalized models reveal stabilizing factors in food webs_ ,  Science 325(5941), 747 (2009).
  * [41] P. Landi, H. O. Minoarivelo, Å. Brännström, C. Hui and U. Dieckmann,  _Complexity and stability of ecological networks: a review of the theory_ ,  Population Ecology 60(4), 319 (2018).
  * [42] C. Song and S. Saavedra,  _Will a small randomly assembled community be feasible and stable?_ ,  Ecology 99(3), 743 (2018).
  * [43] J. M. Yeomans,  _Statistical Mechanics of Phase Transitions_ ,  Clarendon Press (1992).
  * [44] M. Evans,  _MP4/P4 Statistical Physics_ ,  https://www2.ph.ed.ac.uk/~mevans/sp/.
  * [45] R. Wong,  _Asymptotic approximations to integrals_ ,  Academic Press, New York and London (1989).
  * [46] N. G. de Bruijn,  _Asymptotic Methods in Analysis_ ,  Dover books, Amsterdam (2010).
  * [47] M. Evans,  _MP4 Methods of Mathematical Physics_ ,  https://www2.ph.ed.ac.uk/~dmarendu/MOMP/.
  * [48] V. Kaplunovsky (2011), https://web2.ph.utexas.edu/~vadim/Classes/2011f/saddle.pdf.
  * [49] K. Jacobs,  _Stochastic Processes for Physicists: Understanding Noisy Systems_ (2010).
  * [50] C. De Dominicis,  _Techniques de renormalisation de la théorie des champs et dynamique des phénomènes critiques_ ,  Journal de Physique Colloques 37, C1 (1976).
  * [51] A. C. C. Coolen,  _Chapter 15 statistical mechanics of recurrent neural networks ii — dynamics_ ,  In F. Moss and S. Gielen, eds., _Neuro-Informatics and Neural Modelling_ , vol. 4 of _Handbook of Biological Physics_ , pp. 619–684. North-Holland (2001).
  * [52] T. Galla and J. D. Farmer,  _Complex dynamics in learning complicated games_ ,  Proceedings of the National Academy of Sciences 110(4), 1232 (2013).
  * [53] H. Eissfeller and M. Opper,  _New method for studying the dynamics of disordered spin systems without finite-size effects_ ,  Physical Review Letters 68, 2094 (1992).
  * [54] F. L. Paixà,  _Stochastic and deterministic approaches to generalised random Lotka-Volterra communitiesk_ ,  Master’s thesis, Institute for Cross-Disciplinary Physics and Complex Systems IFISC (UIB-CSIC) (2020), https://ifisc.uib-csic.es/en/publications/stochastic-and-deterministic-approaches-to-general/.
  * [55] M. Barbier and J.-F. Arnoldi,  _The cavity method for community ecology_ ,  bioRxiv (2017),  https://www.biorxiv.org/content/early/2017/06/23/147728.
  * [56] L. Poley, J. W. Baron and T. Galla,  _Generalized lotka-volterra model with hierarchical interactions_ ,  Phys. Rev. E 107, 024313 (2023).
  * [57] F. Larroya and T. Galla,  _Demographic noise in complex ecological communities_ ,  Journal of Physics: Complexity 4, 025012 (2023).
  * [58] T. Yoshino, Yoshimiand Galla and K. Tokita,  _Statistical mechanics and stability of a model eco-system_ ,  Journal of Statistical Mechanics: Theory and Experiment 20007, P09003 (2007).
  * [59] M. A. Advani, G. Bunin, and P. Mehta,  _Statistical physics of community ecology: a cavity solution to macarthur’s consumer resource model_ ,  Journal of Statistical Mechanics: Theory and Experiment 2018, 033406 (2018).
  * [60] V. Ros, F. Roy, G. Biroli and G. Bunin,  _Quenched complexity of equilibria for asymmetric generalized lotka–volterra equations_ ,  Journal of Physics A: Mathematical and Theoretical 56(30), 305003 (2023).
  * [61] V. Ros, F. Roy, G. Biroli, G. Bunin and A. M. Turner,  _Generalized lotka-volterra equations with random, nonreciprocal interactions: The typical number of equilibria_ ,  Phys. Rev. Lett. 130, 257401 (2023).
  * [62] J. W. Baron,  _A path integral approach to sparse non-hermitian random matrices_ ,  preprint arxiv:2308.13605 (2023),  https://arxiv.org/abs/2308.13605.
  * [63] J. Garnier-Brun, M. Benzaquen and J.-P. Bouchaud,  _Unlearnable games and "satisficing” decisions: A simple model for a complex world_ ,  preprint arxiv:2312.12252 https://arxiv.org/abs/2312.12252.

##  Appendix A Derivation of the effective process for the Lotka–Volterra
problem without division by species abundances

In this appendix we describe how the generating-functional analysis would have
to be modified in order to avoid having to divide both sides of the
Lotka–Volterra equations (1) by $x_{i}$.

Instead of Eq. (141) we start from

$\displaystyle Z[{\mbox{\boldmath$\psi$}}]$ $\displaystyle=$
$\displaystyle\int{\cal D}\underline{\mathbf{x}}{\cal
D}\widehat{}\underline{\mathbf{x}}\exp\Bigg{(}i\sum_{i}\int
dt\Bigg{[}\widehat{x}_{i}(t)\Bigg{(}\dot{x}_{i}(t)-x_{i}(t)[1-x_{i}(t)+\sum_{j\neq
i}\alpha_{ij}x_{j}(t)+h_{i}(t)]\Bigg{)}\Bigg{]}\Bigg{)}$ (289)
$\displaystyle\times\exp\left(i\sum_{i}\int dt~{}x_{i}(t)\psi_{i}(t)\right),$

We now introduce the following short-hand

$f_{i}[\mathbf{x}(t)]=1-x_{i}(t)+\sum_{j\neq i}\alpha_{ij}x_{j}(t)+h_{i}(t).$
(290)

This can be introduced into the generating functional via suitable delta
functions, and we find

$\displaystyle Z[{\mbox{\boldmath$\psi$}}]$ $\displaystyle=$
$\displaystyle\int{\cal D}\underline{\mathbf{x}}{\cal
D}\widehat{}\underline{\mathbf{x}}~{}{\cal D}\underline{\mathbf{f}}{\cal
D}\widehat{}\underline{\mathbf{f}}~{}\exp\Bigg{(}i\sum_{i}\int
dt\Bigg{[}\widehat{x}_{i}(t)\Bigg{(}\dot{x}_{i}(t)-x_{i}(t)f_{i}(t)\Bigg{)}\Bigg{]}\Bigg{)}$
(291) $\displaystyle\times\exp\left(i\sum_{i}\int
dt~{}\widehat{f}_{i}(t)\left\\{f_{i}(t)-[1-x_{i}(t)+\sum_{j\neq
i}\alpha_{ij}x_{j}(t)+h_{i}(t)]\right\\}\right)$
$\displaystyle\times\exp\left(i\sum_{i}\int dt~{}x_{i}(t)\psi_{i}(t)\right),$

The disorder is now contained in the term $\exp\left(i\sum_{i\neq j}\int
dt~{}\widehat{f}_{i}(t)\alpha_{ij}x_{j}(t)\right)$. This is the equivalent of
the expression in Eq. (143). The calculation then proceeds following the same
steps as in Sec. 6.4, with the replacement
$\widehat{x}_{i}(t)\to\widehat{f}_{i}(t)$.

Following these steps, one ends up with the following disorder-averaged
generating functional of the same form as in Eq. (159),

$\overline{Z}[\underline{{\mbox{\boldmath$\psi$}}}]=\int
D[\underline{M},\widehat{}\underline{M}]D[\underline{P},\widehat{}\underline{P}]D[\underline{\underline{C}},\widehat{}\underline{\underline{C}}]D[\underline{\underline{K}},\widehat{}\underline{\underline{K}}]D[\underline{\underline{L}},\widehat{}\underline{\underline{L}}]\exp\left(N\left[\Psi+\Phi+\Omega+{\cal
O}(N^{-1})\right]\right).$ (292)

The expressions for $\Psi$ and $\Phi$ are as before [see Eqs. (160) and
(161)]. We have

$\displaystyle\Psi$ $\displaystyle=$ $\displaystyle i\int
dt~{}[\widehat{M}(t)M(t)+\widehat{P}(t)P(t)]$ (293) $\displaystyle+i\int
dt~{}dt^{\prime}\left[\widehat{C}(t,t^{\prime})C(t,t^{\prime})+\widehat{K}(t,t^{\prime})K(t,t^{\prime})+\widehat{L}(t,t^{\prime})L(t,t^{\prime})\right],$

and

$\displaystyle\Phi$ $\displaystyle=$ $\displaystyle-\frac{1}{2}\sigma^{2}\int
dt~{}dt^{\prime}\left[L(t,t^{\prime})C(t,t^{\prime})+\Gamma
K(t,t^{\prime})K(t^{\prime},t)\right]$ (294) $\displaystyle-\mu\int
dt~{}M(t)P(t).$

The contribution $\Omega$, however, now takes a slightly different form,

$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle
N^{-1}\sum_{i}\log\bigg{[}\int{\cal D}\underline{x}_{i}{\cal
D}\widehat{}\underline{x}_{i}{\cal D}\underline{f}_{i}{\cal
D}\widehat{}\underline{f}_{i}p_{0}[x_{i}(0)]\exp\left(i\int
dt~{}\psi_{i}(t)x_{i}(t)\right)$ (295) $\displaystyle\times\exp\left(i\int
dt~{}\widehat{x}_{i}(t)\left(\dot{x}_{i}(t)-x_{i}(t)f_{i}(t)\right)\right)\exp{\left(i\int
dt~{}\widehat{f}_{i}(t)\\{f_{i}(t)-[1-x_{i}(t)+h_{i}(t)]\\}\right)}$
$\displaystyle\times\exp\left(-i\int
dt~{}dt^{\prime}\left[\widehat{C}(t,t^{\prime})x_{i}(t)x_{i}(t^{\prime})+\widehat{L}(t,t^{\prime})\widehat{f}_{i}(t)\widehat{f}_{i}(t^{\prime})+\widehat{K}(t,t^{\prime})x_{i}(t)\widehat{f}_{i}(t^{\prime})\right]\right)$
$\displaystyle\times\exp\left(-i\int
dt~{}[\widehat{M}(t)x_{i}(t)+\widehat{P}(t)\widehat{f}_{i}(t)]\right)\bigg{]}.$

Variation with respect to the unhatted order parameters gives the same results
as before [Eq. (164)],

$\displaystyle i\widehat{M}(t)$ $\displaystyle=$ $\displaystyle i\mu P(t),$
$\displaystyle i\widehat{P}(t)$ $\displaystyle=$ $\displaystyle i\mu M(t),$
$\displaystyle i\widehat{C}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sigma^{2}L(t,t^{\prime}),$ $\displaystyle
i\widehat{K}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\Gamma\sigma^{2}K(t^{\prime},t),$ $\displaystyle
i\widehat{L}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sigma^{2}C(t,t^{\prime}).$ (296)

However, the variation with respect to hatted order parameters now leads to

$\displaystyle M(t)$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}N^{-1}\sum_{i}\left\langle{x(t)}\right\rangle_{i},$
$\displaystyle P(t)$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}N^{-1}\sum_{i}\left\langle{\widehat{f}(t)}\right\rangle_{i},$
$\displaystyle C(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}N^{-1}\sum_{i}\left\langle{x(t)x(t^{\prime})}\right\rangle_{i},$
$\displaystyle K(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}N^{-1}\sum_{i}\left\langle{x(t)\widehat{f}(t^{\prime})}\right\rangle_{i},$
$\displaystyle L(t,t^{\prime})$ $\displaystyle=$
$\displaystyle\lim_{N\to\infty}N^{-1}\sum_{i}\left\langle{\widehat{f}(t)\widehat{f}(t^{\prime})}\right\rangle_{i},$
(297)

This is the same as in Eq. (169), but with the replacement
$\widehat{}\underline{x}_{i}\to\widehat{}\underline{f}_{i}$.

Following the steps of Sec. 8.1 one shows that $\underline{P}=\underline{0}$
and $\underline{\underline{L}}=\underline{\underline{0}}$ at the saddle point,
and in the limit
$\underline{{\mbox{\boldmath$\psi$}}}\to\underline{\mathbf{0}}$. Using this
and the remaining saddle-point relations in Eq. (296) one finds (after setting
$h_{i}(t)\equiv h(t)$, and replacing $\underline{\underline{K}}\to
i\underline{\underline{G}}$),

$\left\langle{g[\underline{x},\widehat{}\underline{x},\underline{f},\widehat{}\underline{f}]}\right\rangle_{*}=\frac{\int{\cal
D}\underline{x}{\cal D}\widehat{}\underline{x}{\cal D}\underline{f}{\cal
D}\widehat{}\underline{f}\,g[\underline{x},\widehat{}\underline{x},\underline{f},\widehat{}\underline{f}]\,{\cal
M}[\underline{x},\widehat{}\underline{x},\underline{f},\widehat{}\underline{f}]}{\int{\cal
D}\underline{x}{\cal D}\widehat{}\underline{x}{\cal D}\underline{f}{\cal
D}\widehat{}\underline{f}\,{\cal
M}[\underline{x},\widehat{}\underline{x},\underline{f},\widehat{}\underline{f}]}$
(298)

with the single effective-species measure

$\displaystyle{\cal
M}[\underline{x},\widehat{}\underline{x},\underline{f},\widehat{}\underline{f}]$
$\displaystyle=$ $\displaystyle p_{0}[x(0)]\,\exp\left(i\int
dt~{}\widehat{x}(t)\left[\dot{x}(t)-x(t)f(t)\right]\right)$
$\displaystyle\times\exp\left(i\int
dt~{}\widehat{f}(t)\left(f(t)-[1-x(t)+h(t)]\right)\right)$
$\displaystyle\times\exp\left(-\sigma^{2}\int
dt~{}dt^{\prime}\left[\frac{1}{2}C(t,t^{\prime})\widehat{f}(t)\widehat{f}(t^{\prime})+i\Gamma
G(t^{\prime},t)x(t)\widehat{f}(t^{\prime})\right]\right).$

This describes the process

$\dot{x}(t)=x(t)f(t),$ (300)

where

$f(t)=1-x(t)+\Gamma\sigma^{2}\int
dt^{\prime}G(t,t^{\prime})x(t^{\prime})+\eta(t)+h(t),$ (301)

and where $\eta(t)$ is Gaussian noise with mean zero and
$\left\langle{\eta(t)\eta(t^{\prime})}\right\rangle_{\star}=\sigma^{2}C(t,t^{\prime})$.
Eqs. (300) and (301) can be combined to give

$\dot{x}(t)=x(t)\left[1-x(t)+\Gamma\sigma^{2}\int
dt^{\prime}G(t,t^{\prime})x(t^{\prime})+\eta(t)+h(t)\right],$ (302)

which is the effective process in Eq. (202).

From Eqs. (297) we have the self-consistency relations
$M(t)=\left\langle{x(t)}\right\rangle_{\star}$, and
$C(t,t^{\prime})=\left\langle{x(t)x(t^{\prime})}\right\rangle_{\star}$. Noting
that a derivative of averages $\left\langle{\cdots}\right\rangle_{*}$ with
respect to $h(t^{\prime})$ now ‘brings down’ a factor
$-i\widehat{f}(t^{\prime})$ we also have
$G(t,t^{\prime})=-iK(t,t^{\prime})=-i\left\langle{x(t)\widehat{f}(t)}\right\rangle_{\star}=\frac{\delta}{\delta
h(t^{\prime})}\left\langle{x(t)}\right\rangle_{\star}$.

##  Appendix B Equivalence of averages over the effective dynamics and
averages in the original Lotka–Volterra problem

In this appendix we demonstrate the validity of the identity in Eq. (204),
which we repeat here for convenience,

$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\overline{\left\langle{f[\underline{x}_{i},\widehat{}\underline{x}_{i}]}\right\rangle}=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\left\langle{f[\underline{x},\widehat{}\underline{x}]}\right\rangle_{i}\right|_{\mbox{SP},\underline{\psi}_{i}=0}.$
(303)

The object on the left corresponds to an observable in the original
Lotka–Volterra problem. The object on the right is a quantity to be evaluated
from the effective single-particle measure at the saddle-point of the
generating functional calculation.

It is perhaps not immediately clear in what sense the object on the right is
an observable in the original problem. First we notice that
$\left\langle{\cdots}\right\rangle$ is an average over realisations of the
original dynamics for fixed disorder, and the overbar stands for the average
over the disorder. So this is manifestly for the original Lotka–Volterra
problem. The issue is the occurence of the variables $\widehat{x}_{i}$, which
don’t have a direct physical meaning. What we mean by the notation on the left
is

$\displaystyle\overline{\left\langle{f[\underline{x}_{i},\widehat{}\underline{x}_{i}]}\right\rangle}$
$\displaystyle=$
$\displaystyle\lim_{\underline{{\mbox{\boldmath$\psi$}}}\to\underline{\mathbf{0}}}\Bigg{\\{}\int{\cal
D}\underline{\mathbf{x}}{\cal
D}\widehat{}\underline{\mathbf{x}}\,p_{0}[\mathbf{x}(0)]\,f[\underline{x}_{i},\widehat{}\underline{x}_{i}]\exp\left(i\sum_{i}\int
dt\,\psi_{i}(t)x_{i}(t)\right)$
$\displaystyle\times\overline{\exp\Bigg{(}i\sum_{i}\int
dt\Bigg{[}\widehat{x}_{i}(t)\Bigg{(}\frac{\dot{x}_{i}(t)}{x_{i}(t)}-[1-x_{i}(t)+\sum_{j\neq
i}\alpha_{ij}x_{j}+h_{i}(t)]\Bigg{)}\Bigg{]}\Bigg{)}}\Bigg{\\}}$

(careful not to confuse the overbar with a fraction).

Now, any reasonable function
$f[\underline{x}_{i},\widehat{}\underline{x}_{i}]$ can be expanded into a
power series in terms of the $x_{i}(t)$ and $\widehat{x}_{i}(t)$, so we can
assume that $f$ is a polynomial, with individual terms for example of the form
$x_{i}(t_{1})x_{i}(t_{2})^{2}\widehat{x}_{i}(t_{3})^{2}\widehat{x}_{i}(t_{4})^{3}$.
We can then make the replacement

$f[\underline{x}_{i},\widehat{}\underline{x}_{i}]\to
f\left[-i\frac{\delta}{\delta\underline{\psi}_{i}},\frac{\delta}{-\delta\underline{h}_{i}}\right],$
(305)

i.e.,

$\displaystyle\overline{\left\langle{f[\underline{x}_{i},\widehat{}\underline{x}_{i}]}\right\rangle}$
$\displaystyle=$
$\displaystyle\lim_{\underline{{\mbox{\boldmath$\psi$}}}\to\underline{\mathbf{0}}}\Bigg{\\{}\int{\cal
D}\underline{\mathbf{x}}{\cal
D}\widehat{}\underline{\mathbf{x}}\,p_{0}[\mathbf{x}(0)]\,f\left[-i\frac{\delta}{\delta\underline{\psi}_{i}},\frac{\delta}{-\delta\underline{h}_{i}}\right]\exp\left(i\sum_{i}\int
dt\,\psi_{i}(t)x_{i}(t)\right)$
$\displaystyle\times\overline{\exp\Bigg{(}i\sum_{i}\int
dt\Bigg{[}\widehat{x}_{i}(t)\Bigg{(}\frac{\dot{x}_{i}(t)}{x_{i}(t)}-[1-x_{i}(t)+\sum_{j\neq
i}\alpha_{ij}x_{j}+h_{i}(t)]\Bigg{)}\Bigg{]}\Bigg{)}}\Bigg{\\}}.$

For example
$\overline{\left\langle{x_{i}(t_{1})x_{i}(t_{2})^{2}\widehat{x}_{i}(t_{3})^{3}\widehat{x}_{i}(t_{4})^{2}}\right\rangle}$
turns into $(-i)^{5}\frac{\delta^{5}}{\delta h_{i}(t_{3})^{3}\delta
h_{i}(t_{4})^{2}}\overline{\left\langle{x_{i}(t_{1})x_{i}(t_{2})^{2}}\right\rangle}$,
or alternatively into

$\displaystyle\overline{\left\langle{x_{i}(t_{1})x_{i}(t_{2})^{2}\widehat{x}_{i}(t_{3})^{3}\widehat{x}_{i}(t_{4})^{2}}\right\rangle}$
$\displaystyle=$
$\displaystyle\lim_{\underline{{\mbox{\boldmath$\psi$}}}\to\underline{\mathbf{0}}}\Bigg{\\{}\int{\cal
D}\underline{\mathbf{x}}{\cal
D}\widehat{}\underline{\mathbf{x}}\,p_{0}[\mathbf{x}(0)]\frac{(-i)^{8}\delta^{8}}{\delta\psi_{i}(t_{1})\delta\psi_{i}(t_{2})^{2}\delta
h_{i}(t_{3})^{3}\delta h_{i}(t_{4})^{2}}$ (307)
$\displaystyle\times\exp\left(i\sum_{i}\int
dt\,\psi_{i}(t)x_{i}(t)\right)\overline{\exp\left(\cdots\right)}\Bigg{\\}}$
$\displaystyle=$
$\displaystyle\frac{(-i)^{8}\delta^{8}}{\delta\psi_{i}(t_{1})\delta\psi_{i}(t_{2})^{2}\delta
h_{i}(t_{3})^{3}\delta
h_{i}(t_{4})^{2}}\overline{Z}[\underline{{\mbox{\boldmath$\psi$}}}]\Big{|}_{\underline{\psi}=\underline{\mathbf{0}}}.$

The statement in Eq. (303) is limited to observables which can be written as a
sum over species ($N^{-1}\sum_{i}\dots$) of quantities that relate to single
species only. For example, the object
$N^{-1}\sum_{i}\frac{\delta^{2}}{h_{i}(t_{1})h_{i}(t_{2})}\overline{\left\langle{x_{i}(t_{3})^{2}}\right\rangle}$
qualifies, but
$N^{-2}\sum_{ij}\overline{\left\langle{x_{i}(t_{1})x_{j}(t_{2})}\right\rangle}$
does not (or at least, it is not covered by our derivation).

We therefore start from

$A=\lim_{N\to\infty}\left.\frac{1}{N}\sum_{i}{\cal A}_{\rm
ind}\left(\frac{\delta}{\delta\underline{\psi}_{i}},\frac{\delta}{\delta\underline{h}_{i}}\right)\overline{Z}[\underline{{\mbox{\boldmath$\psi$}}},\underline{\mathbf{h}}]\right|_{\underline{\psi}=\underline{\mathbf{0}}},$
(308)

assuming we are interested in the thermodynamic limit. The object ${\cal
A}_{\rm ind}$ is to be understood as an ‘operator’, composed of derivatives
with respect to $\underline{\psi}_{i}$ and/or $\underline{h}_{i}$. The
subscript ‘ind’ indicates that the operation ${\cal A}_{\rm ind}$ acts on
individual species only. We write ${\cal A}_{i}\equiv{\cal A}_{\rm
ind}\left(\frac{\delta}{\delta\underline{\psi}_{i}},\frac{\delta}{\delta\underline{h}_{i}}\right)$
in the following to keep the notation compact. We also introduce the operator

${\cal A}=\lim_{N\to\infty}\frac{1}{N}\sum_{i}{\cal A}_{i},$ (309)

so that we have

$A={\cal
A}\,\overline{Z}[\underline{{\mbox{\boldmath$\psi$}}},\underline{\mathbf{h}}]\big{|}_{\underline{\psi}=\underline{\mathbf{0}}}.$
(310)

From Eq. (159) we have

$\overline{Z}[\underline{{\mbox{\boldmath$\psi$}}}]=\int
D[\underline{M},\widehat{}\underline{M}]{\cal
D}[\underline{P},\widehat{}\underline{P}]{\cal
D}[\underline{\underline{C}},\widehat{}\underline{\underline{C}}]{\cal
D}[\underline{\underline{K}},\widehat{}\underline{\underline{K}}]{\cal
D}[\underline{\underline{L}},\widehat{}\underline{\underline{L}}]\,\exp\left(N\left[\Psi+\Phi+\Omega+{\cal
O}(N^{-1})\right]\right),$ (311)

and only $\Omega$ contains the $\\{\psi_{i}(t)\\}$ and $\\{h_{i}(t)\\}$. From
Eqs. (163) and (167), we further have

$\Omega=\frac{1}{N}\sum_{j}\ln\,Z_{1}[\underline{\psi}_{j},\underline{h}_{j}].$
(312)

Thus,

$\exp(N\Omega)=\prod_{j}Z_{1}[\underline{\psi}_{j},\underline{h}_{j}].$ (313)

Applying the operator ${\cal A}_{i}$ (for a given fixed $i$) we then have

$\displaystyle{\cal A}_{i}\exp(N\Omega)$ $\displaystyle=$ $\displaystyle{\cal
A}_{i}\left(\prod_{j}Z_{1}[\underline{\psi}_{j},\underline{h}_{j}]\right)$
(314) $\displaystyle=$ $\displaystyle\left(\prod_{j\neq
i}Z_{1}[\underline{\psi}_{j},\underline{h}_{j}]\right){\cal
A}_{i}Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]$ $\displaystyle=$
$\displaystyle\left(\prod_{j}Z_{1}[\underline{\psi}_{j},\underline{h}_{j}]\right)\frac{{\cal
A}_{i}Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}{Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}$
$\displaystyle=$ $\displaystyle\exp(N\Omega)\frac{{\cal
A}_{i}Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}{Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}.$

Thus,

$\displaystyle{\cal A}\overline{Z}[\underline{{\mbox{\boldmath$\psi$}}}]$
$\displaystyle=$ $\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\int
D[\underline{M},\widehat{}\underline{M}]{\cal
D}[\underline{P},\widehat{}\underline{P}]{\cal
D}[\underline{\underline{C}},\widehat{}\underline{\underline{C}}]{\cal
D}[\underline{\underline{K}},\widehat{}\underline{\underline{K}}]{\cal
D}[\underline{\underline{L}},\widehat{}\underline{\underline{L}}]\,{\cal
A}_{i}\exp\left(N\left[\Psi+\Phi+\Omega+{\cal O}(N^{-1})\right]\right)$ (315)
$\displaystyle=$ $\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\int
D[\underline{M},\widehat{}\underline{M}]{\cal
D}[\underline{P},\widehat{}\underline{P}]{\cal
D}[\underline{\underline{C}},\widehat{}\underline{\underline{C}}]{\cal
D}[\underline{\underline{K}},\widehat{}\underline{\underline{K}}]{\cal
D}[\underline{\underline{L}},\widehat{}\underline{\underline{L}}]\,\Bigg{\\{}\frac{{\cal
A}_{i}Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}{Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}$
$\displaystyle\hskip 80.00012pt\times\exp\left(N\left[\Psi+\Phi+\Omega+{\cal
O}(N^{-1})\right]\right)\Bigg{\\}}.$

Using the rules of saddle-point integration (and the fact that the generating
functional is normalised), we arrive at

$\displaystyle{\cal A}\overline{Z}[\underline{{\mbox{\boldmath$\psi$}}}]$
$\displaystyle=$
$\displaystyle\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left.\frac{{\cal
A}_{i}Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}{Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]}\right|_{\underline{\psi}=\underline{\mathbf{0}},\mbox{\footnotesize
SP}}.$ (316)

We next use the definition of $Z_{1}$ in Eq. (167) as well as that of the
average $\left\langle{\cdots}\right\rangle_{i}$ in Eq. (170). The operator
${\cal A}_{i}$, applied to $Z_{1}[\underline{\psi}_{i},\underline{h}_{i}]$,
‘brings down’ the same powers of $\underline{x}$ and $\widehat{}\underline{x}$
as those of $\underline{x}_{i}$ and $\widehat{}\underline{x}_{i}$ in the
original definition of $f[\underline{x}_{i},\widehat{}\underline{x}_{i}]$ [see
e.g. the text just after Eq. (B)]. We therefore have

$A=\lim_{N\to\infty}\frac{1}{N}\sum_{i}\left\langle{f[\underline{x},\widehat{}\underline{x}]}\right\rangle_{i},$
(317)

which is what we wanted to show.