[1]Francesco Caravelli

These authors contributed equally to this work.

[1]T-4, T-Division, Los Alamos National Laboratory,
Bikini Atoll rd, Los Alamos, 87544, NM, USA

2]CCS-3, CCS-Division, Los Alamos National Laboratory,
Bikini Atoll rd, Los Alamos, 87544, NM, USA

# On the combinatorics of Lotka-Volterra equations

<EMAIL_ADDRESS>Yen Ting Lin<EMAIL_ADDRESS>*  [

###### Abstract

We study an approach to obtaining the exact formal solution of the 2-species
Lotka - Volterra equation based on combinatorics and generating functions. By
employing a combination of Carleman linearization and Mori-Zwanzig reduction
techniques, we transform the nonlinear equations into a linear system,
allowing for the derivation of a formal solution. The Mori-Zwanzig reduction
reduces to an expansion which we show can be interpreted as a directed and
weighted lattice path walk, which we use to obtain a representation of the
system dynamics as walks of fixed length. The exact solution is then shown to
be dependent on the generator of weighted walks. We show that the generator
can be obtained by the solution of PDE which in turn is equivalent to a
particular Koopman evolution of nonlinear observables.

###### keywords:

Lotka-Volterra equations, Generating functions, Koopman evolution, formal
solution, Carleman linearization, Mori-Zwanzig

###### Contents

  1. 1 Introduction
  2. 2 Formal solution method
    1. 2.1 Linearization
    2. 2.2 The Mori-Zwanzig formalism for linear system
  3. 3 Simple example: $\dot{x}=x^{2}$
    1. 3.1 Carleman linearization and Mori-Zwanzig reduction
    2. 3.2 Series summation
    3. 3.3 Generating function of directed lattice walks
    4. 3.4 Closed form PDE for the generating function
    5. 3.5 Connection to the Koopman’s representation
    6. 3.6 Properties of the formal solution
  4. 4 The Lotka-Volterra equations
    1. 4.1 Carleman linearization of LV
    2. 4.2 Mori-Zwanzig reduction
    3. 4.3 Directed walks on lattice
    4. 4.4 Example of lattice coefficient calculation
    5. 4.5 Lotka-Volterra via generating functions
    6. 4.6 Comments on Lotka-Volterra with N species and hyperlattices
  5. 5 Conclusions

##  1 Introduction

The two species Lotka-Volterra equations are the archetypal nonlinear
equations, exhibiting both laminar and oscillatory behavior. These are also
known as predator-prey equations and are a pair of coupled differential
equations used to model the interactions between two species in an ecosystem.
Developed independently by Alfred J. Lotka and Vito Volterra in the early 20th
century [1, 2, 3, 4], these equations provide a mathematical framework for
understanding the dynamics of predator-prey relationships and to model the
nonlinearity of co-dependent species. This model has experienced some renewed
interest in later years because random multi-species models exhibit marginally
stable equilibria [5]. The techniques used to study these problems come from
disordered systems, which have shed some light on a century-old problem. While
the LV equations have been extensively studied, and various analytical and
numerical methods have been employed to find solutions and analyze their
behavior, we believe that the results and combinatorial structure presented in
this paper have never been explored in the past [6, 7].

The Lotka-Volterra equations assume two main populations: the predator and the
prey. They describe how the populations of these species change over time
based on their interactions and the availability of resources. In the most
general form, the equations are as follows:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle\alpha x+\beta xy$
(1) $\displaystyle\frac{dy}{dt}$ $\displaystyle=$ $\displaystyle\gamma
y+\delta xy$ (2)

where $x$ and $y$ represent the real-valued population densities of the prey
(e.g. rabbits) and predator (e.g. wolves) species, respectively. The
quantities $dx/dt$ and $dy/dt$ represent the rates of change of the population
densities over time. The parameters $\alpha>0$ and $\gamma<0$ represent the
per-capita growth rate of the prey and predator population in the absence of
predation, respectively (e.g. rabbits reproduce indefinitely in the absence of
wolves, while wolves die off in the absence of rabbits). The nonlinear terms
are controlled by the parameters $\beta$ and $\delta$: $\beta<0$ and
$\delta>0$ represent the rate of change of prey and predator population (per
unit of prey population, per unit predator population) due to predation.
Above, one often assumes the ecologically meaningful parameter regime:
$\alpha>0$, $\beta<0$, $\gamma<0$, and $\delta>0$, where oscillatory dynamics
can occur as a result of feedback, instead of unstable growth or complete
decay. This paper is predominantly focused on the two species case, but we
will later discuss some features of higher dimensional LVs as well in view of
the combinatorial analysis of this paper.

The LV equations illustrate the complex interplay between predators and prey
in an ecosystem. They demonstrate how fluctuations in one population can
influence the dynamics of the other population. These equations have been used
to study various ecological systems, such as the interaction between hares and
lynx [8] or bacteria and Bacteriophages [9], and provide insights into
population dynamics, cycles, and stability [10]. While the LV equations offer
a simplified representation of predator-prey dynamics, they provide a valuable
tool for understanding the fundamental principles of ecological interactions
and have contributed to developing more sophisticated ecological models (e.g.,
[11, 12]). This is also why understanding the structure underlying the two-
species model can provide insights into more complex ecological networks [13].

The present manuscript is the result of an attempt to highlight certain
features of coupled and nonlinear ODE that have been neglected. Using the
Lotka-Volterra (LV) equations as our model, we study the relationship between
quantities that are combinatorial in nature, such as generating functions of
weighted walks, the Koopman operator, and how these two are related formal
exact solution of the two species LV.

The relationship between Koopman evolution [14, 15, 16], Mori-Zwanzig methods
[17, 18, 19], and ordinary differential equations (ODEs) lies at the
intersection of dynamical systems theory and statistical mechanics [20]. The
Koopman operator is a mathematical tool used in the study of dynamical
systems. It describes the evolution of observables (functions) over time
without explicitly solving the underlying ODEs. The Koopman operator provides
a linear representation of the dynamics, enabling the analysis of complex
systems often through spectral methods and linear algebra techniques. It is
particularly useful for systems with high-dimensional state spaces. On the
other hand, Mori-Zwanzig methods are a set of mathematical techniques used in
statistical mechanics to derive effective equations of motion for a system
with many degrees of freedom [21, 22, 23, 24, 25, 26, 27]. The main idea is to
systematically eliminate the “fast” degrees of freedom (unresolved variables)
to obtain a reduced model that only considers the “slow” (resolved) or
relevant variables. This reduced model is often in the form of ODEs and
captures the long-term behavior of the system while discarding fast
oscillations or fluctuations. The Koopman operator and Mori-Zwanzig methods
are related in the context of understanding the dynamics of complex systems
[28, 29]. The Koopman operator provides an abstract and linear perspective on
the system’s evolution [30, 31], while Mori-Zwanzig methods focus on obtaining
reduced models for high-dimensional systems in terms of ODEs. The Koopman
operator can be used to study the dynamics of observables, and by exploiting
its spectral properties, one can identify relevant slow modes that capture the
system’s long-term behavior. The Koopman operator has been used extensively in
fluid dynamics [32], where it has been linked to the Dynamic Mode
Decomposition [33, 34]. By constructing a reduced model using Mori-Zwanzig
methods based on these slow modes, one can effectively obtain a system of ODEs
that approximates the original dynamics.

As we show in this paper for two sets of differential equations, the
additional insight that we gain in this manuscript is that the Koopman
evolution is indeed also the same operator that generates certain weighted and
directed lattice walks, which directly enter into the formal solution of these
equations. These lattice walks can be thought to be associated with monomials
in the initial conditions.

The paper is organized as follows. We first introduce in section 2 the two key
techniques that we use to formally solve the equations, given by a Carleman
linearization [35, 36] first in Sec. 2.1, followed by a Mori-Zwanzig reduction
in Sec. 2.2. We then first solve the well-known equation $\dot{x}=x^{2}$ in
Sec 3.1-3.2-3.3, showing that the solution has an interpretation in terms of a
one-dimensional lattice walk, and show that the solution depends on the
generator of lattice walks. We derive a closed-form PDE for these generators
in Sec. 3.4, with an analysis of the utility of this representation in Sec.
3.6. We then focus on the Lotka-Volterra equations in Sec. 4. We introduce the
variables used for the Carleman linearization in Sec. 4.1, and perform the
Mori-Zwanzig reduction in 4.2, introducing the expression of the formal
solution of the LV. In Sec. 4.3 we provide the interpretation of the formal
solution in terms of weighted lattice walks, we introduce the lattice
coefficients, showing that we can in principle use Monte Carlo techniques to
obtain the solution in Sec. 4.4. In Sec. 4.5 we show that the solution can be
analyzed in terms of the generating function of lattice walks, and provide a
closed-form quasi-linear PDE for the solution of the memory kernel of the
solution, which is the generating function lattice walks, and analyzing the
formal solution of the PDE in terms of the Lagrange-Charpit method. In Sec.
4.6 we provide a brief analysis of the N-species Lotka-Volterra equations.
Conclusions follow.

##  2 Formal solution method

We begin by providing the main ingredients which allow to formally write the
solution of the LV equations. We first discuss Carleman linearization, and
then discuss the Mori-Zwanzig formalism.

###  2.1 Linearization

Carleman linearization [35] is a powerful technique used to transform
nonlinear differential equations into linear ones [36], facilitating their
analysis and solution. This method involves introducing additional variables
and expressing the original nonlinear equations as a linear system with
respect to these new variables. By doing so, Carleman linearization enables
the use of well-established linear techniques to investigate the behavior and
properties of the system. In this paper, we utilize Carleman linearization to
obtain lattice path expansions, which serve as the foundation for deriving the
exact formal solution of the 2-species Lotka-Volterra equation. Consider a
system of differential equations of the form

$\frac{d\vec{x}}{dt}=\vec{f}(\vec{x}).$ (3)

We can linearize this system by considering a set of observables given by the
terms of observables of the form $r_{\vec{a}}=\prod_{i}x_{i}^{a_{i}}$ with
$a_{i}\in\mathbb{N}$. Then, from eqn. (3), expanding $f_{i}(\vec{x})$
(assuming they are analytical), one also writes the infinite chain of
equations in terms of $r_{i}$. If $f$ is time-invariant, then the result
equation $\dot{\vec{r}}=A\vec{r}$ is also time-invariant, and formally the
solution for this system is known. However, the chain of equations is
infinite-dimensional. We use instead a variant of this technique, which is the
Mori-Zwanzig formalism. However, the key property of quadratic equations such
as Lotka-Volterra, as we show below, is that they are amenable to a
combinatorial treatment.

###  2.2 The Mori-Zwanzig formalism for linear system

Let us consider the simplest example of Mori-Zwanzig coarse-graining [21]. We
have a linear ODE system of the form

$\displaystyle\frac{d}{dt}\vec{x}=A\vec{x}$ (4)

with $\vec{x}(0)=\vec{x}_{0}$. We wish to express the system’s dynamics in a
generalized Langevin formalism. We then assume that our observables are a
subset of the vector $\vec{x}$, e.g.

$\displaystyle\vec{x}(t)=\begin{pmatrix}\vec{y}(t)\\\ \vec{z}(t)\end{pmatrix}$
(5)

with $\vec{x}\in\mathcal{R}^{N}$, while $\vec{y}\in\mathcal{R}^{N-m}$ and
$\vec{z}\in\mathcal{R}^{m}$. We partition the matrix $A$ in blocks, so that

$\displaystyle A=\begin{pmatrix}A_{rr}&A_{ru}\\\ A_{ur}&A_{uu}\end{pmatrix}$
(6)

where $A_{rr}$ is $(N-m)\times(N-m)$, $A_{ru}$ is $N\times(m\times m)$,
$A_{ur}$ is $m\times(N-m)$ and $A_{uu}$ is of size $(m\times m)$. Then, we can
write

$\displaystyle\frac{dy}{dt}$ $\displaystyle=$ $\displaystyle
A_{rr}\vec{y}+A_{ru}\vec{z}$ (7) $\displaystyle\frac{dz}{dt}$ $\displaystyle=$
$\displaystyle A_{ur}\vec{y}+A_{uu}\vec{z}.$ (8)

In the blocks above, $r$ stands for resolved variables, while $u$ stands for
unresolved variables. Our observables, thus the resolved variables, are the
$y$ components. We thus write the formal expression for $z$, which we assume
to be unresolved, e.g.

$z(t)=e^{A_{uu}t}\vec{z}_{0}+\int_{0}^{t}ds\ e^{A_{uu}(t-s)}A_{ur}\vec{y}(s).$
(9)

We now plug this expression into the first set of equations, obtaining

$\displaystyle\frac{dy}{dt}$ $\displaystyle=$ $\displaystyle
A_{rr}\vec{y}+\int_{0}^{t}ds\
A_{ru}e^{A_{uu}(t-s)}A_{ur}\vec{y}(s)+A_{ru}e^{A_{uu}t}\vec{z}_{0}$ (10)

We then identify the noise term $F(t)=A_{ru}e^{A_{uu}t}\vec{z}_{0}$ and the
kernel operator, given by

$\displaystyle K(t-s)=A_{ru}e^{A_{uu}(t-s)}A_{ur}.$ (11)

from which then we obtain

$\displaystyle\frac{d\vec{y}}{dt}$ $\displaystyle=$ $\displaystyle
A_{rr}\vec{y}+\int_{0}^{t}ds\ K(t-s)\vec{y}(s)+\vec{F}(t),$ (12)

which is in the generalized Langevin equation. Thus, our problem reduces to a
simpler problem to analyze than the generalized Langevin equation, and the
memory kernel can be written down explicitly.

##  3 Simple example: $\dot{x}=x^{2}$

The techniques introduced in the previous sections can be showcased on a
particular ODE, whose solution has many of the key characteristics of the LV
equations. In particular, after having applied the Carleman linearization and
the Mori-Zwanzig formalism, we will see that the solution of this equation
obtained via power series resummation can be shown to be connected to the
generator of (weighted) lattice walks in one dimension. We consider the
nonlinear ODE

$\displaystyle\frac{dx}{dt}=x^{2},\ x(0)=x_{0}.$ (13)

The solution is given by $x(t)=x_{0}(1-x_{0}t)^{-1}$. To see this, we can
write

$\frac{1}{x^{2}}\frac{d}{dt}x=\frac{d}{dt}(-\frac{1}{x})=1$ (14)

and thus

$\displaystyle\frac{1}{x_{0}}-\frac{1}{x(t)}=t$ (15)

which, inverting, gives

$\displaystyle x(t)=\frac{x_{0}}{1-x_{0}t}.$ (16)

###  3.1 Carleman linearization and Mori-Zwanzig reduction

We now use the Carleman linearization, we have that the observables are the
Taylor powers, given by

$\displaystyle r_{k}=x^{k}.$ (17)

It is not hard to see that the equation above can be written in the form

$\displaystyle\frac{d}{dt}\vec{r}=A\vec{r}$ (18)

with $A_{ij}=i\delta_{i+1,j}$. We will treat $r_{1}$ as the resolved
observable, and the rest $r_{k}$, $k\geq 2$ as the unresolved observables.

We now use eqn. (10). First, note that $A_{rr}=A_{ur}=0$. Then, our equation
for the observable $r_{1}$ reads

$\displaystyle\frac{dr_{1}}{dt}=\frac{dx}{dt}=A_{ru}e^{A_{uu}t}\vec{z}_{0}.$
(19)

and then

$\displaystyle x(t)=x_{0}+A_{ru}\int^{t}ds\ e^{A_{uu}s}\vec{z}_{0}.$ (20)

In the equation above, $\vec{z}_{0}=(x_{0}^{2},x_{0}^{3},\cdots)$. The vector
$A_{ru}$ is infinite and looks like $A_{ru}=(1,0,\cdots,0)$ while $A_{uu}$ is
a square super-diagonal matrix, where the super diagonal is infinite, and
given by $(2,3,4,\cdots).$ We must then calculate exp$(A_{uu}t)$.

###  3.2 Series summation

We use the expression

$\displaystyle e^{A_{uu}t}=\sum_{k=0}^{\infty}\frac{t^{k}}{k!}A_{uu}^{k},$
(21)

for which we need to calculate an expression for the first row of the matrix,
as

$\displaystyle
A_{ru}\sum_{k=0}^{\infty}\frac{t^{k}}{k!}A_{uu}^{k}\vec{z}_{0}=\sum_{j}(e^{A_{uu}t})_{1j}x_{0}^{j+1}.$
(22)

Note that

$\displaystyle(A_{uu})_{ij}=(i+1)\delta_{i+1,j}$ (23)

Then,

$\displaystyle(A_{uu}^{r})_{ab}$ $\displaystyle=$
$\displaystyle\sum_{k_{1}\cdots
k_{r-1}}(A_{uu})_{ak_{1}}(A_{uu})_{k_{1}k_{2}}\cdots(A_{uu})_{k_{r-1}b}$ (24)
$\displaystyle=$ $\displaystyle(a+1)(a+2)\cdots(a+r)\delta_{a+r+1,b}$
$\displaystyle=$ $\displaystyle\frac{(a+r)!}{a!}\delta_{a+r,b}$

Then,

$\displaystyle(e^{A_{uu}t})_{ab}=\sum_{r=0}^{\infty}\frac{t^{r}}{r!}(A_{uu}^{r})_{ab}=\sum_{r=0}^{\infty}\frac{(a+r)!t^{r}}{a!r!}\delta_{a+r,b}$
(25)

from which we get then

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle
x_{0}+\int^{t}ds\sum_{r=0}^{\infty}\frac{(r+1)!s^{r}}{1!r!}\sum_{b=1}^{\infty}\delta_{r+1,b}x_{0}^{b+1}$
(26) $\displaystyle=$ $\displaystyle
x_{0}+\sum_{r=0}^{\infty}\frac{(r+1)!t^{r+1}}{r!(r+1)}x_{0}^{r+2}$ (27)
$\displaystyle=$
$\displaystyle\sum_{r=0}^{\infty}x_{0}^{r+1}t^{r}=x_{0}\sum_{r=0}^{\infty}x_{0}^{r}t^{r}=\frac{x_{0}}{1-x_{0}t}.$
(28)

which is exactly the solution of $x^{\prime}=x^{2}$ with $x(0)=x_{0}$.The
solution above was derived already in [37] along similar lines. We show now
that there is a completely different approach to obtaining the solution.

Figure 1: The path graph of the $\dot{x}=x^{2}$ model

###  3.3 Generating function of directed lattice walks

Here we wish to show that we can express the solution in terms of the
generating function of a walk on a path graph.

We simply define a walk as a weighted sequence of nodes. To see this, we first
introduce the path graph, as illustrated in Fig. 1. The path graph is a
directed graph, whose nodes are observables $r_{k}:=x^{k}$ and weighted edges
are the interactions between the observables, i.e., the weight between nodes
$r_{k}$ and $r_{k+1}$ is $k$ (because $\dot{r_{k}}=kr_{k+1}$). We further
define $\mathcal{N}^{k}_{(i,j)}$, $i,j\in\mathbb{N}$, as the product of
weights along the (only) path of length $k$ connecting nodes $r_{i}$ and
$r_{j}$. However, since our graph is directed, we must have
$\mathcal{N}^{k}_{(i,j)}=\delta_{k,|i-j|}\mathcal{M}(i,j)$, where

$\mathcal{M}(i,j):=\frac{(j-1)!}{(i-1)!}.$ (29)

We now make a connection between the path graph and the key term in Eq. (22),
$\left(A_{uu}\right)^{k}_{1,j}$:

$\left(A_{uu}\right)^{k}_{1,j}=\sum_{i_{1}}\ldots\sum_{i_{k}}\left(A_{uu}\right)_{1,i_{1}}\ldots\left(A_{uu}\right)_{i_{k},j}.$
(30)

Because of the super-diagonal structure of the matrix $A_{uu}$ (see Eq. (23)),
each of the $A_{uu}$ in the above product can be interpreted as a weight of a
single step moving from a node $r_{j}$ to the next $r_{j+1}$. As such, the
whole sum can be interpreted as the product of the weights along a path with
length $k$ connecting nodes $r_{2}$ and $r_{j+1}$ (note that the indices of
unresolved observables begin with $2$:

$\left(A_{uu}\right)^{k}_{1,j}=\delta_{k,j-1}\mathcal{M}(2,j+1)=\delta_{k,j-1}j!.$
(31)

Then, Eq. (22) can be expressed straightforwardly

$\displaystyle
A_{ru}\sum_{k=0}^{\infty}\frac{t^{k}}{k!}A_{uu}^{k}\vec{z}_{0}{}$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{j=1}^{\infty}\frac{t^{k}}{k!}\left(A_{uu}^{k}\right)_{j}x_{0}^{j+1}$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{j=1}^{\infty}\frac{t^{k}}{k!}\delta_{k,j-1}\mathcal{M}(2,j+1)x_{0}^{j+1}=\sum_{k=0}^{\infty}\frac{t^{k}}{k!}\mathcal{M}(2,k+2)x_{0}^{k+2}.$
(32)

We now define a generating function [38] for the weighted walks on the path
graph:

$G(s,x):=\sum_{k=0}^{\infty}\frac{s^{k}x^{k+2}}{k!}\mathcal{M}(2,k+2),$ (33)

which can be calculated analytically, and we explicitly computed it in Sec.
3.2. One can see that the Taylor expansion of $G$ in the parameter $s$
controls the length of the walk, while the Taylor expansion in $x$ controls
between which nodes of the graph the walk occurs, as the terms monomials
$s^{k}x^{k+2}$ multiply the number of walks $\mathcal{N}(2,k+2)$ between node
$2$ and $2+k$. Plugging these in the definition Eq. (29), we obtain:

$G(s,x)=x^{2}\sum_{k=0}^{\infty}\frac{(sx)^{k}}{k!}(k+1)!=x^{2}\sum_{k=0}^{\infty}(k+1)(sx)^{k}=\frac{x^{2}}{(1-sx)^{2}}.$
(34)

We can then interpret the solution of the differential equation as the time
integrand of the generating function of these weighted walks evaluated at
$x=x_{0}$:

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle
x_{0}+\int_{0}^{t}G(x_{0},s)ds$ (35) $\displaystyle=$ $\displaystyle
x_{0}+\frac{tx_{0}^{2}}{1-tx_{0}}=\frac{x_{0}}{1-tx_{0}}$ (36)

which is the expression we had obtained before, but interpreted now in terms
of the generating function.

###  3.4 Closed form PDE for the generating function

We now establish that the generating function satisfies a Partial Differential
Equation (PDE), which is what we will use for the Lotka–Volterra equation
below. Applying $\partial_{x}$ and $\partial_{s}$ to the generating function
Eq. (33), we obtain:

$\displaystyle\partial_{x}G(s,x)={}$ $\displaystyle
x\sum_{k=0}^{\infty}(k+2)\frac{t^{k}x^{k}}{(k-1)!}\mathcal{M}(2,k+2)=x^{3}(k+2)\sum_{k=0}^{\infty}\frac{t^{k}x^{k}}{k!}\mathcal{M}(2,k+2),$
(37) $\displaystyle\partial_{s}G(s,x)={}$ $\displaystyle
x^{2}\sum_{k=1}^{\infty}\frac{t^{k-1}x^{k}}{(k-1)!}\mathcal{M}(2,k+2)=x^{3}\sum_{k=0}^{\infty}(k+2)\frac{t^{k}x^{k}}{k!}\mathcal{M}(2,k+2).$
(38)

In the derivation above, we have used the identity

$\displaystyle\mathcal{M}(j,k+3)=\left(k+2\right)\mathcal{M}(j,k+2).$ (39)

Clearly, the above equations suggest the closed-form PDE

$\displaystyle\partial_{t}\big{(}G(x,t)\big{)}=x^{2}\partial_{x}G(x,t).$ (40)

The solution of the differential equation can be derived using the method of
characteristics. We impose the Lagrange-Charpit equality

$\displaystyle\frac{dt}{-1}=\frac{dx}{x^{2}}=\frac{dG}{0}$ (41)

From the first two, we get

$\displaystyle-t+\Phi=-1/x.$ (42)

From the latter, we get $dG=0\rightarrow G=\Gamma$. Imposing the Cauchy
surface condition that $\Gamma$ depends on $\Phi$,
$G=\Gamma(\Phi)=f(\frac{1-tx}{x})$. Thus, $G(x,t)=f(\frac{1-tx}{x})$. We know
however that $G(x,0)=x^{2}$, and thus $f(\cdot)=(\cdot)^{-2}$. Then,
$G(x,t)=x^{2}/(1-xt)^{2}$ which is the function we found earlier while summing
the series explicitly. Then, the interpretation of the solution of Eq. (35)
can also be applied to the case of other nonlinearly coupled equations.

###  3.5 Connection to the Koopman’s representation

The formal solution Eq. (35) can be derived from yet another method, which is
arguably the most formal and general derivation. Such a derivation is tightly
connected to Koopman’s representation [14, 15], with which one aims to
prescribe the evolutionary equation of the observables. Let us consider a
general observable $g:\mathbb{R}^{1}\rightarrow\mathbb{R}^{1}$. In the Koopman
picture, the states do not move and stay at $x_{0}$, the initial condition,
but the observable $g$ becomes time-dependent, denoted by
$g_{t}:=\mathcal{K}_{t}g$ where $\mathcal{K}_{t}$ is the finite-time Koopman
operator. We write a field $\psi(x,t)$ satisfying the following equation

$\displaystyle\partial_{t}\psi(x,t)={}$ $\displaystyle\mathcal{L}\psi(x,t),$
(43a) $\displaystyle\psi(x,0)={}$ $\displaystyle g(x),$ (43b)

where $\mathcal{L}$ is the (backward) generator of the process,
$\mathcal{L}=x^{2}\partial_{x}$. The above PDE and the associated initial
condition are often referred to as the Liouville equation [39, 40], albeit
this nomenclature is not standard in non-equilibrium statistical physics111In
statistical physics, the Liouville equation is the forward equation of the
dynamics, where the state variables are time-dependent and the observables are
static, which is the adjoint of Eq. (43). We note that the structure of the
PDE is exactly the one derived earlier.

The solution of

$\psi(x_{0},t)=\left(\mathcal{K}_{t}g\right)(x_{0})\equiv
g\left(x(t;x_{0})\right)$ (44)

contains all the information of the ODE system, with an arbitrary $g$. For
example, the formal state solution $x(t)$ can be obtained by setting
$g(x):=x^{2}$, as will be seen below. One can also consider the indicator
function $g(x):=\delta(x-x^{\prime})$ to probe if the solution is at a query
point $x^{\prime}$ at an arbitrary time $t>0$.

Interestingly, from this Koopman viewpoint, the PDE above is often solved
implicitly in terms of the ODE solution. In general, solving the backward PDE
(43) is a challenging task. For the $\dot{x}=x^{2}$ model, we have the ODE
solution

$x(t)=\frac{x_{0}}{1-x_{0}t}.$ (45)

In line with the method of characteristics, we can “pull-back” by express the
$x(t)$ as a function of the initial condition $x_{0}$:

$\psi(x,t)=g\left(\frac{x}{1-tx}\right).$ (46)

One can show that the PDE is satisfied:

$\displaystyle\partial_{t}\psi(x,t)={}$ $\displaystyle
g^{\prime}\left(\frac{x}{1-tx}\right)\frac{x^{2}}{\left(1-tx\right)^{2}}$
(47a) $\displaystyle x^{2}\partial_{x}\psi(x,t)={}$ $\displaystyle
x^{2}g^{\prime}\left(\frac{x}{1-tx}\right)\left[\frac{1}{1-tx}+\frac{xt}{\left(1+tx\right)^{2}}\right]$
(47b) $\displaystyle=$ $\displaystyle
g^{\prime}\left(\frac{x}{1-tx}\right)\frac{x^{2}}{\left(1-tx\right)^{2}}.$
(47c)

To connect to the solution in the previous calculation, we note that the
formal solution of the ODE can be written as

$x(t;x_{0})=x_{0}+\int_{0}^{t}x^{2}\left(s\right)\,ds,$ (48)

With this expression, it is clear that the observable of interest should be
$g(x):=x^{2}$, leading to the solution

$\psi(x,t)=g\left(\frac{x}{1-tx}\right)=\left(\frac{x}{1-tx}\right)^{2},$ (49)

which is exactly the generating function Eq. (33).

###  3.6 Properties of the formal solution

We note that this approach allows us to write explicitly the properties of the
solution as a function of the initial conditions. For instance, for the
example above, we have

$\displaystyle\partial_{x_{0}}x(t)$ $\displaystyle=$ $\displaystyle
1+\int_{0}^{t}\partial_{x_{0}}G(x_{0},t)dt$ $\displaystyle=$ $\displaystyle
1+\frac{1}{x_{0}^{2}}\int_{0}^{t}\partial_{t}G(x_{0},t)dt$ $\displaystyle=$
$\displaystyle 1+\frac{G(x_{0},t)-G(x_{0},0)}{x_{0}^{2}}$ $\displaystyle=$
$\displaystyle
1+\frac{G(x_{0},t)-x_{0}^{2}}{x_{0}^{2}}=\frac{G(x_{0},t)}{x_{0}^{2}}$ (51)

which can be promptly checked using the solution. It is also interesting to
note that there is another way of deriving this PDE, knowing the expression
for the solution. Using Eqn. (3.6), and from the fundamental theorem of
calculus, we can equate

$\displaystyle x^{2}=x^{\prime}(t)=G(x,t)$ (52)

which we note to be also satisfied. We know from the exact solution that
$G(x,t)=x_{0}^{2}/(1-x_{0}t)^{2}=x(t)^{2}$. Now note that from the expression
above, we can obtain

$\displaystyle\partial_{t}G=2xx^{\prime}=2x^{3}=x^{2}\partial_{x}G,$ (53)

which is exactly the PDE we derived using the generating function method.

We wish to show that a similar approach is also possible for the Lotka-
Volterra equations, where however the combinatorial structure is slightly more
involved, and where the solution is only formal.

##  4 The Lotka-Volterra equations

After having discussed at length the techniques used for the simpler toy model
of eqn. (13), let us consider the case of the Lotka-Volterra equation (LV).
The solution of the quadratic case shares in fact many similarities with the
formal solution we obtain below.

First, we begin by writing the LV equations in the form

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle\alpha x+\beta
xy\equiv f(x,y)$ $\displaystyle\frac{dy}{dt}$ $\displaystyle=$
$\displaystyle\gamma y+\delta xy\equiv g(x,y).$ (54)

In the typical formulation of the LVs, $\alpha>0$, $\beta<0$, $\gamma<0$, and
$\delta>0$.

###  4.1 Carleman linearization of LV

As for the case of the toy model, we perform a Carleman linearization, where
now the variables are the following monomials

$\displaystyle r_{(a,b)}=x^{a}y^{b}.$ (55)

We have then

$\displaystyle\frac{dr_{(a,b)}}{dt}$ $\displaystyle=$ $\displaystyle
ax^{a-1}(\beta xy+\alpha x)y^{b}+by^{b-1}(\delta xy+\gamma y)x^{a}$ (56)
$\displaystyle=$ $\displaystyle b\delta x^{a+1}y^{b}+a\beta
x^{a}y^{b+1}+x^{a}y^{b}(a\alpha+b\gamma)$ $\displaystyle=$ $\displaystyle
a\beta r_{(a,b+1)}+b\delta r_{(a+1,b)}+(a\alpha+b\gamma)r_{(a,b)}$

then, if we consider the combined index $(a,b)$, it follows that the Lotka-
Volterra equations, the linearized LV equation becomes then

$\displaystyle\frac{dr_{(0,1)}}{dt}$ $\displaystyle=$ $\displaystyle\delta
r_{(1,1)}+\gamma r_{(0,1)}$ $\displaystyle\frac{dr_{(1,0)}}{dt}$
$\displaystyle=$ $\displaystyle\beta r_{(1,1)}+\alpha r_{(1,0)}$
$\displaystyle\frac{dr_{(1,1)}}{dt}$ $\displaystyle=$ $\displaystyle\beta
r_{(1,2)}+\delta r_{(2,1)}+(\alpha+\gamma)r_{(1,1)}$
$\displaystyle\frac{dr_{(1,2)}}{dt}$ $\displaystyle=$ $\displaystyle\beta
r_{(1,3)}+2\delta r_{(2,2)}+(\alpha+2\gamma)r_{(1,2)}$
$\displaystyle\frac{dr_{(2,1)}}{dt}$ $\displaystyle=$ $\displaystyle 2\beta
r_{(2,2)}+\delta r_{(3,1)}+(2\alpha+\gamma)r_{(2,1)}$ (57)
$\displaystyle\vdots$

It is important to stress that the one above is an exact representation of the
Lotka-Volterra system. The LV equations can represent a path graph, which is
the directed two-dimensional lattice shown in Fig. 2.

###  4.2 Mori-Zwanzig reduction

Because of the linear relationship between the variables, we can write the
equations in the form

$\displaystyle\frac{d}{dt}\vec{r}=A\vec{r},$ (58)

where the vector $(\vec{r})_{(a,b)}(t)=x^{a}(t)y^{b}(t)$, and where the matrix
$A$ is an upper triangular matrix that can be expressed as a weighted directed
graph like the one shown in Fig. 2. Let us call
$\vec{v}=(r_{(1,0)},r_{(0,1)})\equiv(x,y)$ and
$\vec{z}=(r_{(1,1)},r_{(2,1)},r_{(1,2)},\cdots)$, so that
$\vec{r}^{t}=(\vec{v}^{t},\vec{z}^{t})$. Using equation (10), we have

$\displaystyle\frac{d\vec{z}}{dt}$ $\displaystyle=$ $\displaystyle
A_{ru}e^{A_{uu}t}\vec{z}$ (59)

which is decoupled from the variables $r_{(1,0)}$ and $r_{(0,1)}$.

Figure 2: The matrix $A$ from the Carleman linearization in transition between
the variables $r_{(a,b)}$ for a portion of the wedge containing the variables
$x=r_{(1,0)}$ and $y=r_{(0,1)}$. The variable $r_{(0,0)}$ are constants and
thus $dr_{(0,0)}/dt=0$. The pink variables are not involved in the lattice
calculations.

The matrix $A$, whose first two rows and columns correspond to $r_{(0,1)}$ and
$r_{(1,0)}$, given the ordering $r_{(0,1)}$, $r_{(1,0)}$, $r_{(1,1)}$,
$r_{(2,1)}$, $r_{(1,2)}$, can be written in the form:

$\displaystyle A=\left[\begin{array}[]{cc|cccc}\gamma&0&\delta&0&0&\cdots\\\
0&\alpha&\beta&0&0&\cdots\\\ \hline\cr
0&0&(\alpha+\gamma)&\delta&\beta&\cdots\\\ 0&0&0&(\alpha+2\gamma)&0&\cdots\\\
0&0&0&0&(2\alpha+\gamma)&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right]$ (66)

From the structure of the matrix above, we see that the solution can be
written in the form

$\displaystyle\frac{d\vec{v}}{dt}$ $\displaystyle=$ $\displaystyle
A_{rr}\vec{v}+A_{ru}e^{A_{uu}t}\vec{z}_{0}$

The equation above is a linear system of 2 equations with a non-homogeneous
forcing. The set of linearly coupled differential equations above can be
solved analytically, with a solution of the form

$\vec{v}(t)=e^{A_{rr}t}\vec{v}_{0}+e^{A_{rr}t}\int_{0}^{t}ds\
e^{-sA_{rr}}A_{ru}e^{A_{uu}s}\vec{z}_{0}.$ (67)

The simplicity of the solution is only apparent. While $A_{rr}$ is a two-
dimensional diagonal matrix, $A_{uu}$ is a semi-infinite matrix. While semi-
infinite, the elements of the matrix $(e^{A_{uu}s})_{(a,b)(c,d)}$ can be
calculated combinatorially, analyzing
$e^{A_{uu}s}=\sum_{k=0}^{\infty}(A_{uu})^{k}\frac{s^{k}}{k!}$. We first
discuss the interpretation of this matrix.

###  4.3 Directed walks on lattice

Let us then consider $(A_{uu})^{k}_{(a,b),(c,d)}$. Since the graph
representing $A$ is directed, we must have $c>a\geq 0$, $d>b\geq 0$ in our
notation. In particular, we must have $|a-c|+|b-d|\leq k$. If
$\alpha=\gamma=0$, then the inequality becomes an inequality, and this becomes
a problem of counting paths on the lattice from $(a,b)$ to $(c,d)$. Let us now
formalize the combinatorics of such a process. We call $\mathcal{R}=c-a$ and
$\mathcal{D}=d-b$. Then, in order for a directed walk to start from $(a,b)$
and reach $(c,d)$ in $k$ steps, no matter the order, if we go down
$\mathcal{D}$ and right $\mathcal{R}$ times, we will definitely get to
$(c,d)$. We can stop anytime then, for a total amount of
$k-\mathcal{U}-\mathcal{D}=\mathcal{S}$ times, for a walk of $k$ steps. Walks
can then be mapped as permutations of a vector of moves Down, Right, and Stop
of the form

$\displaystyle\left(\underbrace{D\cdots D}_{\mathcal{D}\
\textit{times}}\underbrace{R\cdots R}_{\mathcal{R}\
\text{times}}\underbrace{S\cdots S}_{\mathcal{S}\ \text{times}}\right)$ (68)

Given a configuration, and depending on the edges it goes through, then we
will have a numerical coefficient $C_{\sigma}$ which depends on the
permutation222In fact, in this representation, one has to divide by a symmetry
factor corresponding to the number permutations corresponding to the same
path. If for instance, we have $(DDRRSS)$ the sequence with which we perform
the first or second down move is actually independent. Thus, we must divide a
path like the one above by $1/2^{3}.$, a coefficient
$\beta^{\mathcal{R}}\delta^{\mathcal{D}}$ and a coefficient $\rho$ which
depends on the $\alpha$ and $\gamma$ values and product of the weights along
the path. We note that in the limit $\alpha,\gamma\rightarrow 0$, $\rho$
becomes $1$ if $\mathcal{R}+\mathcal{D}=k$ otherwise, e.g. the walk has to
start from $(a,b)$ and reach the final $(c,d)$ point in exactly $k$ moves. We
can then write

$(A^{k}_{22})_{(a,b)(c,d)}=\sum_{\sigma\in\mathcal{S}_{k}}C^{k}_{(c,d)}(\sigma)\beta^{c-a}\delta^{d-b}\rho(\sigma)\theta_{k}(\sigma)(c-a)\theta(d-b)$
(69)

where is the Heaviside function, defined as $\theta(x)=1$ for $x\geq 0$ and
zero otherwise and $\mathcal{S}_{k}$ is the group of permutations of $k$
elements. We then have

$(e^{A_{uu}s})_{(a,b)(c,d)}=\sum_{k=0}^{\infty}\sum_{\sigma\in\mathcal{S}_{k}}C^{k}_{(c,d)}(\sigma)\beta^{c-a}\delta^{d-b}\rho_{k}(\sigma)\theta(c-a)\theta(d-b)\frac{s^{k}}{k!}.$
(70)

Let us now write eqn. (67) explicitly. This is a 2-component equation, while
$\vec{z}_{0}=(x_{0}y_{0},x_{0}^{2}y_{0},x_{0}y_{0}^{2},\cdots)$ is the vector
of the appropriate powers of initial conditions. First, we note that
$e^{-sA_{rr}}$ is a diagonal matrix with elements $e^{-s\gamma}$,
$e^{-s\alpha}$, corresponding to $(0,1)$, $(1,0)$. Similarly,
$A_{ru}e^{A_{uu}s}$ selects only the first two rows of $e^{A_{uu}s}$ with the
appropriate weights. Then, we have

$e^{-sA_{rr}}A_{ru}e^{A_{uu}s}\vec{z}_{0}=\sum_{c,d}\begin{pmatrix}\delta
e^{-s\gamma}(e^{A_{uu}s})_{(1,1),(c,d)}(\vec{z}_{0})_{(c,d)}\\\ \beta
e^{-s\alpha}(e^{A_{uu}s})_{(1,1),(c,d)}(\vec{z}_{0})_{(c,d)}\\\ \end{pmatrix}$
(71)

Then, if we define the diagonal matrix $\text{diag}(\delta e^{-s\gamma},\beta
e^{-s\alpha})$ and

$\eta(s,x_{0},y_{0})=\sum_{c,d}(e^{A_{uu}s})_{(1,1),(c,d)}(\vec{z}_{0})_{(c,d)}$
(72)

we get

$e^{-sA_{rr}}A_{ur}e^{A_{uu}s}\vec{z}_{0}=\eta(s)\begin{pmatrix}\delta
e^{-s\gamma}&0\\\ 0&\beta e^{-s\alpha}\end{pmatrix}\vec{1}.$ (73)

and thus the solution can be written as

$\displaystyle\begin{pmatrix}y(t)\\\
x(t)\end{pmatrix}=\begin{pmatrix}e^{\gamma t}y_{0}\\\ e^{\alpha
t}x_{0}\end{pmatrix}+\int_{0}^{t}ds\ \eta(s,x_{0},y_{0})\begin{pmatrix}\delta
e^{\gamma(t-s)}\\\ \beta e^{\alpha(t-s)}\end{pmatrix}$ (74)

We can see from the equation above that for $t=0$ we recover the initial
conditions.

From a perturbative standpoint, we see that in the equation above we need to
calculate integrals of the form

$\displaystyle I_{m}(t)=\int_{0}^{t}\eta(s,x_{0},y_{0})e^{-ms}ds$ (75)

for arbitrary $m$. If we write

$\displaystyle\eta(s)=\sum_{c,d\geq
1}\sum_{k=0}^{\infty}\frac{s^{k}}{k!}\sum_{\sigma\in\mathcal{S}_{k}}C^{k}_{(c,d)}(\sigma)\beta^{c-1}\delta^{d-1}\rho_{k}(\sigma)(\vec{z}_{0})_{(c,d)}.$
(76)

We see that the expansion can be reduced to expressions of the form

$\displaystyle I^{k}_{m}(t)=\int_{0}^{t}s^{k}e^{-ms}ds.$ (77)

Furthermore, $I^{k}_{m}(t)$ can be expressed as

$I^{k}_{m}(t)=\frac{\Gamma(k+1)-\Gamma(k+1,mt)}{m^{k+1}},\text{ if }m>0,$ (78)

where $\Gamma(k+1)=(k+1)!$ and $\Gamma(k+1,mt)=\int_{t}^{\infty}e^{-s}s^{k}\
ds$ is the incomplete Gamma function. Note that for small values of $t$,
$I^{k}_{m}(t\rightarrow 0^{+})={(mt)^{k}}/{(1+k)}+O(t^{k+1})$, and thus at
small values of $t$ the expansion also corresponds to a time Taylor expansion.
Instead, in the limit $m\rightarrow 0$, we have from the properties of the
$\Gamma$ functions that

$\displaystyle\lim_{m\rightarrow 0}I_{m}^{k}(t)=\frac{t^{1+k}}{1+k}.$ (79)

which is consistent with setting $m=0$ in eqn. (78) and integrating $s^{k}$.

The summation over $c$ and $d$ can be swapped with the summation over $k$ and
$\sigma$, from which we get

$\displaystyle\Gamma_{k}=\sum_{c=1}^{\infty}\sum_{d=1}^{\infty}\sum_{\sigma\in\mathcal{S}_{k}}\delta^{d-1}y_{0}^{d}\beta^{c-1}x_{0}^{c}C^{k}_{(c,d)}(\sigma)\rho_{k}(\sigma)$
(80)

from which we can write

$\displaystyle y(t)$ $\displaystyle=$ $\displaystyle e^{\gamma t}y_{0}+\delta
e^{\gamma t}\Omega_{\gamma}(t)$ (81) $\displaystyle x(t)$ $\displaystyle=$
$\displaystyle e^{\alpha t}x_{0}+\beta e^{\alpha t}\Omega_{\alpha}(t)$ (82)

where

$\Omega_{m}(t)=\sum_{k=0}^{\infty}\Gamma_{k}I^{k}_{m}(t)$ (83)

which is the final expression for the lattice paths expansion solution of the
Lotka-Volterra equation. We define the parameters $\Gamma_{k}$ as the lattice
coefficients. However, we wish to stress that Eq. (78) is true only for $m>0$.
Because conventional parametrization demands $\gamma<0$, one must rely on Eq.
(77) to evaluation of $I^{k}_{\gamma}(t)$. It is clear that
$I^{k}_{\gamma}(t)$ is not converging for larger values of $t$ for any
$k\in\mathbb{Z}_{\geq 0}$. Consequently, this formal method by lattice
coefficients is unfortunately of little use, if not only formally.

This is now a good moment to stop and focus on the meaning of these
expressions. Equation (74), albeit formal, provides an interpretation of the
exact solution of the 2-species Lotka-Volterra equations in terms of directed,
and weighted, walks on lattice paths. As a result, the function $\eta(s,x,y)$
can in principle be estimated via a Monte Carlo approach, with random lattice
walks stopping at $(a,b)$ in order to estimate the contributions of order
$x_{0}^{a}y_{0}^{a}$ in the initial condition. This is the first result of
this paper. Yet, admittedly this is a quite cumbersome approach. Although in
the next section, we provide examples for the calculation of the lattice
coefficients, it is worth mentioning that the purpose of the rest of the paper
will be to show that $\eta$ is the solution of a particular differential
equation.

###  4.4 Example of lattice coefficient calculation

The solutions of eqn. (81)-(82) are, as a matter of fact, written explicitly
in terms of coefficients that need to be evaluated at every order $k$ of the
lattice expansion. However, unless an exact expression for $\Omega_{m}(t)$ can
be found, this is only formal for the time being. In fact, the symmetric group
$\mathcal{S}_{k}$ on a set of $k$ elements has order $k!$. It is then
immediate to see that the complexity of each order increases dramatically.
However, we still think that this approach can lead to a way of formalizing
this solution.

To show how this expansion works, let us consider the first few orders. We
define the boundary set $\mathcal{B}_{k}$ of the nodes in the lattice that can
be reached in $k$ steps. For $k=0$, we identify $\mathcal{S}_{0}=\emptyset$.
In this case $\Gamma_{0}=r_{(1,1)}=x_{0}y_{0}$. At $k=1$,
$\mathcal{S}_{1}=\\{Id\\}$, the boundary of nodes that can be reached after
one step is composed of $\mathcal{B}_{1}=\\{r_{(1,2)},r_{(2,1)}\\}$. It is
then easy to see that $\sigma=Id$, $\Gamma_{1}=\beta x_{0}y_{0}^{2}+\delta
x_{0}^{2}y_{0}$. The first non-trivial case is $k=2$. The boundary of the
nodes that can be reached is
$\mathcal{B}_{2}=\\{r_{(1,2)},r_{(2,1)},r_{(2,2)}\\}$. The set
$\mathcal{S}_{2}=\\{\sigma_{1}=\\{1,2\\},\sigma_{2}=\\{2,1\\}\\}$.

For the nodes $r_{(1,2)},r_{(2,1)}$ the walk has to stop once at least and go
right and down respectively, while for $r_{(2,2)}$ once right and down at
least.

Thus, the transitions are given by

  1. 1.

for $r_{(1,2)}=x_{0}y_{0}^{2}$, our moves must be $\sigma_{1}(R,S)=(R,S)$ or
$\sigma_{2}(R,S)=(S,R)$, thus $C_{2}(\sigma_{1})=C_{2}(\sigma_{2})=1$, while
$\rho_{2}(\sigma_{1})=(\alpha+\gamma)$,
$\rho_{2}(\sigma_{2})=(\alpha+2\gamma)$, and $\beta^{c-1}\delta^{d-1}=\beta$;

  2. 2.

for $r_{(2,1)}=x_{0}^{2}y_{0}$ we have $\sigma_{1}(D,S)=\\{D,S\\}$ and
$\sigma_{2}(D,S)=(S,D)$, thus $C_{2}(\sigma_{1})=1=C_{2}(\sigma_{2})=1$, while
$\rho_{2}(\sigma_{1})=\alpha+\gamma$ and
$\rho_{2}(\sigma_{2})=2\alpha+\gamma$, with $\beta^{c-1}\delta^{d-1}=\delta$;

  3. 3.

for $r_{(2,2)}$ we must have $\sigma_{1}(R,D)=(R,D)$ and
$\sigma_{2}(R,D)=(D,R)$, and $C_{2}(\sigma_{1})=C_{2}(\sigma_{2})=2$, and
$\rho_{2}=1$ and $\beta^{c-1}\delta^{d-1}=\beta\delta$.

These correspond then to

$\displaystyle\Gamma_{2}$ $\displaystyle=$
$\displaystyle\beta\big{(}(\alpha+\gamma)+(\alpha+2\gamma)\big{)}x_{0}y_{0}^{2}$
$\displaystyle+\delta\big{(}(\alpha+\gamma)+(2\alpha+\gamma)\big{)}x_{0}^{2}y_{0}$
$\displaystyle+\big{(}(\beta)(2\delta)+(2\beta)(\delta)\big{)}x_{0}^{2}y_{0}^{2}$
$\displaystyle=$
$\displaystyle\beta(2\alpha+3\gamma)x_{0}y_{0}^{2}+\delta(3\alpha+2\gamma)x_{0}^{2}y_{0}+4\beta\delta
x_{0}^{2}y_{0}^{2}.$

We then get that the expansion up to order $k=2$ is given by

$\displaystyle y(t)$ $\displaystyle=$ $\displaystyle e^{\gamma t}y_{0}+\delta
e^{\gamma t}\Big{(}x_{0}y_{0}\frac{1-e^{-\gamma
t}}{\gamma}+\Gamma_{1}I_{\gamma}^{1}(t)+\Gamma_{2}I_{\gamma}^{2}(t)\Big{)}$
$\displaystyle+O(t^{3})$ $\displaystyle x(t)$ $\displaystyle=$ $\displaystyle
e^{\alpha t}x_{0}+\beta e^{\alpha t}\Big{(}x_{0}y_{0}\frac{1-e^{-\alpha
t}}{\alpha}+\Gamma_{1}I_{\alpha}^{1}(t)+\Gamma_{2}I_{\alpha}^{2}(t)\Big{)}$
$\displaystyle+O(t^{3})$

which provides an example of the form of the solution. Then up to order two,
the solution above versus the numerically obtained one is shown in Fig. 3. We
see that this expansion is problematic at a perturbative level, because of its
slow convergence, and for $\beta<0,\gamma<0$, it has an alternating sign which
does not improve the convergence either. Albeit in the following we focus on
deriving an exact equation for $\eta$, the approach described above allows in
principle to study the solution of LV equations using Monte Carlo techniques.

Figure 3: Analytical at order $k$ versus forward Euler method with $dt=0.01$
(solid lines), for the parameters $\alpha=1$, $\beta=-0.5$, $\gamma=-0.5$,
$\delta=2$. We see that the expansion captures the short time behavior of the
system.

###  4.5 Lotka-Volterra via generating functions

The analysis of the previous section has highlighted the importance of
deriving exact equations for $\eta$. As we saw in the plots above, the
perturbative expansion is slowly converging. We thus need a better expression
for the solution.

We now know however that $\eta(t,x,y)$ in eqn. (72) is exactly the generator
function of walks of length $k$ from $(1,1)$ to $(c,d)$, and the coefficient
corresponding to $t^{k}/k^{2}(x^{c}y^{d})$ is the weighted walk of length $k$.
Similar to what we introduced earlier for the simpler model in eqn. (29) for
the model of eqn. (13), we now have a weighted walk coefficient
$\mathcal{N}^{k}_{(x_{1},y_{1}),(x_{2},y_{2})}$ which is the sum of all
possible walks from vertex $(x_{1},y_{1})$ to vertex $(x_{2},y_{2})$ on a
lattice. In our case, since we are interested in walks of length $k$ starting
always from node $(1,1)$, we can simply define
$\mathcal{N}^{k}_{(1,1),(a,b)}=N^{k}_{(a,b)}$.

We can then write, similarly to what we did for the simpler ODE of eqn. (13),
a recursion relation of the form

$N^{k}_{(a,b)}=(\alpha a+b\gamma)N^{k-1}_{(a,b)}+a\beta
N^{k-1}_{(a,b-1)}+b\delta N^{k-1}_{(a-1,b)},$ (85)

where $N^{k}_{a,b}$ is the product of the weights along a path of length $k$
connecting $(1,1)$ and $(a,b)$ (the dependence of $(1,1)$ is suppressed for
brevity). It is clear that we have the following conditions: (1)
$N^{0}_{(1,1)}=1$, (2) $N^{k}_{(0,1)}=N^{k}_{(1,0)}=0$, $k\in\mathbb{Z}_{\geq
0}$, and (3) $N^{r}_{(a,b)}=0$ if $r\neq a+b-2$, $a,b\in\mathbb{N}$. Then, the
generating function can be written in the form

$\displaystyle\eta(t,x,y)=\sum_{a,b=1}^{\infty}\sum_{k=0}^{\infty}N^{k}_{(a,b)}\frac{t^{k}}{k!}x^{a}y^{b}.$
(86)

We multiply the left and right of the recursion relation Eq. (85) by
$t^{k-1}x^{a}y^{b}/{(k-1)!}$ and sum over $a$, $b$, and $k$ from $1$ to
$\infty$. Let us analyze this term by term. On the left-hand side, we obtain
$\sum_{a,b=1}^{\infty}\sum_{k=1}^{\infty}t^{k-1}x^{a}y^{b}N_{(a,b)}^{k}/{(k-1)!}=\partial_{t}\eta(t,x,y).$
On the right hand side, the summation over $k$ is exactly right to give the
$t$ dependence, and we focus on the $a,b$ summations.

$\displaystyle\alpha:$
$\displaystyle\sum_{a,b=1}^{\infty}\sum_{k=1}^{\infty}\frac{t^{k-1}}{(k-1)!}ax^{a}y^{b}N_{(a,b)}^{k-1}=x\partial_{x}\eta(t,x,y)$
$\displaystyle\gamma:$
$\displaystyle\sum_{a,b=1}^{\infty}\sum_{k=1}^{\infty}\frac{t^{k-1}}{(k-1)!}bx^{a}y^{b}N_{(a,b)}^{k-1}=y\partial_{y}\eta(t,x,y)$
$\displaystyle\delta:$
$\displaystyle\sum_{a,b=1}^{\infty}\sum_{k=1}^{\infty}\frac{t^{k-1}}{(k-1)!}{b}x^{a}y^{b}N_{(a-1,b)}^{k-1}=xy\partial_{y}\eta(t,x,y)$
$\displaystyle\beta:$
$\displaystyle\sum_{a,b=1}^{\infty}\sum_{k=1}^{\infty}\frac{t^{k-1}}{(k-1)!}{a}x^{a}y^{b}N_{(a,b-1)}^{k-1}=x{y}\partial_{x}\eta(t,x,y),$

from which we obtain the PDE

$\displaystyle\partial_{t}\eta(t,x,y)$ $\displaystyle=$
$\displaystyle\Big{(}(\alpha x+\delta xy)\partial_{x}+(\gamma y+\beta
yx)\partial_{y}\Big{)}\eta(s,x,y)$ (87) $\displaystyle=$
$\displaystyle\Big{(}f(x,y)\partial_{x}+g(x,y)\partial_{y}\Big{)}\eta(t,x,y)$
$\displaystyle\eta(0,x,y)$ $\displaystyle=$ $\displaystyle xy$ (88)
$\displaystyle\eta(t,0,y)$ $\displaystyle=$ $\displaystyle\eta(t,x,0)=0$ (89)

where the boundary conditions are derived as follows. Since
$N^{k}_{(1,0)}=N^{r}_{(0,1)}=0$, it follows that $\eta(s,0,y)=\eta(s,x,0)=0$.
Since $t=0$ must imply that the only surviving term is $N^{0}_{(1,1)}xy$ with
$N^{0}_{(1,1)}=1$, $\eta(0,x,y)=xy$ must be supplied as initial condition.

The PDE can also be derived from the formal Koopman’s picture, analogous to
the technique presented in Sec. 3.5. To get to the same form as the PDE above,
we first rewrite the ODE to absorb the linear terms

$\displaystyle\frac{d}{dt}\left(e^{\alpha t}x(t)\right)={}$
$\displaystyle\beta e^{\alpha t}x(t)\,y(t),$ (90a)
$\displaystyle\frac{d}{dt}\left(e^{\gamma t}x(t)\right)={}$
$\displaystyle\delta e^{\gamma t}x(t)\,y(t).$ (90b)

To this end, we can see that both terms involve $x(t)y(t)$, hinting that the
key observable is $g(x,y):=xy$. As such, the backward PDE (Eq. (43)) is

$\displaystyle\partial_{t}\eta(x,y,t)={}$ $\displaystyle\left[\left(\alpha
x+\beta xy\right)\partial_{x}+\left(\gamma y+\delta
xy\right)\partial_{y}\right]\eta(x,t),$ (91a) $\displaystyle\eta(x,y,0)={}$
$\displaystyle g(x,y)=xy,$ (91b)

and once $\eta$ is solved (recall that
$\eta(x_{0},y_{0},t)=x(t;x_{0})\,y(t;y_{0})$, we can integrate Eq. (90):

$\displaystyle x(t)={}$ $\displaystyle
x(0)+\beta\int_{0}^{t}e^{\alpha(t-s)}\eta(x_{0},y_{0},t)\,ds,$ (92a)
$\displaystyle y(t)={}$ $\displaystyle
y(0)+\gamma\int_{0}^{t}e^{\gamma(t-s)}\eta(x_{0},y_{0},t)\,ds.$ (92b)

which is exactly the same as the solution obtained from combining Carleman
linearization, Mori-Zwanzig formalism, and generating function of the directed
weighted graph.

To see that the equation above is correct, let us consider inserting
explicitly the formal solution (Eq. (74)) into the ODE. We have

$\displaystyle\alpha x+\beta xy$ $\displaystyle=$ $\displaystyle
x^{\prime}(t)$ (93) $\displaystyle=$ $\displaystyle\alpha e^{\alpha
t}x_{0}+\beta\left[\alpha\int_{0}^{t}\eta(x_{0},y_{0},t)e^{\alpha(t-s)}ds+\eta(x_{0},y_{0},t)\right]$

We now note that on the right-hand side, we have

$\alpha x+\beta xy=x^{\prime}(t)=\alpha x(t)+\beta\eta(x_{0},y_{0},t)$ (94)

from which we obtain the identity

$\eta(x_{0},y_{0},t)\equiv\eta(x(t),y(t))=x(t)y(t)$ (95)

we now take the derivative with respect to time, obtaining

$\displaystyle\partial_{t}\eta(x,y)$ $\displaystyle=$ $\displaystyle
x^{\prime}(t)y(t)+x(t)y^{\prime}(t)$ $\displaystyle=$ $\displaystyle(\alpha
x+\beta xy)y(t)+(\gamma y+\delta xy)x(t)$

Now note that we can write

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle\partial_{y}\eta$ (96)
$\displaystyle y(t)$ $\displaystyle=$ $\displaystyle\partial_{x}\eta$ (97)

From which we obtain

$\displaystyle\partial_{t}\eta=(\alpha x+\beta xy)\partial_{x}\eta+(\gamma
y+\delta xy)\partial_{y}\eta$ (98)

Which is the expression obtained using the generating function method.

The solution of the PDE above can be solved (formally, and in principle) on a
Cauchy surface determined by constants of integration emerging from the
Lagrange-Charpit equations, as the above is a quasi-linear equation, exactly
as in the case of the simpler example described earlier. However, our attempts
to solve the equation analytically lead only to a complicated implicit form
that we decided to omit in the present manuscript.

###  4.6 Comments on Lotka-Volterra with N species and hyperlattices

Such methodology for the generation of exact solutions of Lotka-Volterra
equations can be generalized to higher dimensions. We consider then the set of
coupled differential equations of the form

$\frac{dx_{i}}{dt}=\alpha_{i}x_{i}+x_{i}\sum_{j}^{N}\beta_{ij}x_{j}$ (99)

where the matrix $\beta$ is zero on the diagonal. Similarly to what we had
done in the previous section, we consider the observables
$r_{\vec{a}}=\prod_{i=1}^{N}x_{i}^{a_{i}}$, and focus on their time
derivatives

$\displaystyle\frac{d}{dt}r_{\vec{a}}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N}a_{j}x_{j}^{a_{j-1}}\frac{dx_{j}}{dt}\prod_{r\neq
j}x_{r}^{a_{r}}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N}a_{j}x_{j}^{a_{j-1}}(\alpha_{j}x_{j}+x_{j}\sum_{k}^{N}\beta_{jk}x_{k})\prod_{r\neq
j}x_{r}^{a_{r}}$ $\displaystyle=$
$\displaystyle(\sum_{j}a_{j}\alpha_{j})r_{\vec{a}}+\sum_{j=1}^{N}a_{j}\sum_{k}^{N}\beta_{jk}r_{\vec{a}+\hat{x}_{k}}$
$\displaystyle=$
$\displaystyle\tilde{\alpha}_{a}r_{\vec{a}}+\sum_{k=1}^{N}\tilde{\beta}_{k}r_{\vec{a}+\hat{x}_{k}}.$
(101)

As a result, we obtain a lattice equation in higher dimensions. Equation (68)
has to be generalized, as now we have $N$ dimensions for the lattice path. It
is known that in general this construction can be expressed in a tensorial
representation [41].

The key difference between the $2-$species case and the higher dimensional one
is that the lattice depends on the initial point, which makes the analysis
slightly more complicated. For $3-$species, the initial points of the lattice
paths are given by $x=r_{(1,0,0)},y=r_{(0,1,0)}$ and $z=r_{(0,0,1)}$,
analogously to $x=r_{(1,0)}$ and $y=r_{(0,1)}$ in Fig. 2 in the $2-$species.
It can be promptly seen that, since $\dot{x}$ has the monomials $x$,$xy$ and
$xz$ on the right hand side, $\dot{y}$ has the monomials $y$,$yx$ and $yz$ and
$\dot{z}$ instead $z$, $zx$ and $zy$, the lattice expansion implies

$\displaystyle r_{1,0,0}$ $\displaystyle\rightarrow$ $\displaystyle
r_{1,0,0},r_{1,1,0},r_{1,0,1}$ (102) $\displaystyle r_{0,1,0}$
$\displaystyle\rightarrow$ $\displaystyle r_{0,1,0},r_{1,1,0},r_{0,1,1}$ (103)
$\displaystyle r_{0,0,1}$ $\displaystyle\rightarrow$ $\displaystyle
r_{0,0,1},r_{1,0,1},r_{0,1,1}$ (104)

From this we see that the lattice walks for the resolved variables begin at
three sublattices starting at $r_{0,1,1}$, $r_{1,0,1}$ and $r_{1,1,0}$ which
do not coincide as in the two-dimensional case. For comparison, in the $2-$
species case we had

$\displaystyle r_{1,0}$ $\displaystyle\rightarrow$ $\displaystyle
r_{1,0},r_{1,1}$ (105) $\displaystyle r_{0,1}$ $\displaystyle\rightarrow$
$\displaystyle r_{0,1},r_{1,1}$ (106)

which implies that we could write the solution in terms of a single generating
function $\eta$ corresponding to walks from $r_{1,1}$ to $r_{a,b}$. However,
since the lattice paths for a higher number of species are directed, the
sublattices resulting from these paths do not coincide in this case. Thus, in
the multi-species case, one has to resort to multiple generating functions. In
the $3-$species case, these are $\eta_{xy}$, $\eta_{xz}$ and $\eta_{yz}$,
related to walks from $r_{1,1,0}$, $r_{1,0,1}$ and $r_{0,1,1}$ to a generic
point $r_{a,b,c}$ respectively. Although by symmetry these three must be
related, we see immediately that the hyperlattices require a separate
discussion that goes beyond the purpose and space limitations of this note. We
can see however that our methodology can also be applied here, with due
attention.

##  5 Conclusions

This paper presents a novel approach to analyzing the Lotka-Volterra equations
by introducing a lattice path expansion, preceded by a Carleman linearization
and by a Mori-Zwanzig reduction which we used to derive the formal solution.
We first analyzed the simpler model $\dot{x}=x^{2}$, which although had been
solved before, as we have shown her the solution has an interpretation in
terms of lattice path expansion. As far as we know, such lattice structure
employed in this study was not previously known or utilized in the simpler
model examined in the paper, nor had it been applied to the Lotka-Volterra
equations. Our paper elucidates the interplay between the formal solutions of
the ODEs, the generating functions of the lattice walks, and how these can be
interpreted in terms of Koopman evolution.

This approach provides a powerful method for understanding the Lotka-Volterra
equations. Indeed, it came as a posteriori surprise that the generators of
these lattice path expansions have also, in turn, an interpretation in terms
of Koopman evolution operators, applied to a particular observable.

Some critical comments are in order. Although we could not provide an exact
solution to the two-dimensional Lotka-Volterra equations, what we think to be
an important contribution of this paper is the connection between
combinatorial methods and the solution of nonlinear ODE, and in turn the
relationship between the generator of walks on graphs to the Koopman evolution
of observables. We think of the Carleman and Mori-Zwanzig techniques as a
method to obtain the formal solution that in principle can be applied, clearly
with harder work, to a large class of ODE that can be expressed in terms of
monomials of the dynamical variables. The formal solution of the LV equations
should be seen only as that, as it is written in terms of the generator of
lattice walks on the graph. Solving the PDE of such a generator is as hard, if
not more difficult, than solving the ODE directly. In future works we will
focus on the numerical analysis of the PDE for the LV equation, trying to find
connections between the transition to oscillatory regime in the Lotka-Volterra
ODEs and to properties of the related PDE.

In future work, we will focus on the numerical integration of eqn. (98) trying
to elucidate these connections.

Acknowledgments

The work of FC and YTL was carried out under the auspices of the NNSA of the
U.S. DoE at LANL under Contract No. DE-AC52-06NA25396, and in particular
support from LDRD via 20220063DR, 20230338ER, and 20230627ER.

## References

  *   * Lotka [1925] Lotka, A.J.: Elements of physical biology. Nature 116 (1925)
  * Lotka [1926a] Lotka, A.J.: Fluctuations in the abundance of a species considered mathematically. Nature 119(2966), 558–560 (1926)
  * Lotka [1926b] Lotka, A.J.: Contribution to the theory of periodic reactions. Journal of Physical Chemistry 14(3), 271–274 (1926)
  * Volterra [1928] Volterra, V.: Variations and fluctuations of the number of individuals in animal species living together. Journal de Conseil International pour l’Exploration de la Mer 3(1), 3–51 (1928)
  * Biroli et al. [2018] Biroli, G., Bunin, G., Cammarota, C.: Marginally stable equilibria in critical ecosystems. New Journal of Physics 20(083051) (2018)
  * Iannelli and Pugliese [2014] Iannelli, M., Pugliese, A.: An Introduction to Mathematical Population Dynamics. Springer, Cham (2014)
  * Baigent [2016] Baigent, S.A.: Lotka-Volterra Dynamics: An Introduction. World Scientific, ??? (2016)
  * Gilpin [1973] Gilpin, M.E.: Do hares eat lynx? The American Naturalist 107(957), 727–730 (1973) https://doi.org/10.1086/282870
  * Maslov and Sneppen [2017] Maslov, S., Sneppen, K.: Population cycles and species diversity in dynamic Kill-the-Winner model of microbial ecosystems. Scientific Reports 7(1), 39642 (2017) https://doi.org/10.1038/srep39642
  * Strogatz [2000] Strogatz, S.H.: Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering, 1st edn. CRC press, ??? (2000)
  * Martcheva [2015] Martcheva, M.: Ecological Context of Epidemiology, pp. 249–279. Springer, Boston, MA (2015). https://doi.org/10.1007/978-1-4899-7612-3_10
  * Grunert et al. [2021] Grunert, K., Holden, H., Jakobsen, E.R., Stenseth, N.C.: Evolutionarily stable strategies in stable and periodically fluctuating populations: The rosenzweig–macarthur predator–prey model. Proceedings of the National Academy of Sciences 118(4), 2017463118 (2021) https://doi.org/10.1073/pnas.2017463118
  * Grilli et al. [2017] Grilli, J., Barabás, G., Michalska-Smith, M.J., Allesina, S.: Higher-order interactions stabilize dynamics in competitive network models. Nature 548(7666:210) (2017)
  * Koopman and Neumann [1932] Koopman, B.O., Neumann, J.v.: Dynamical systems of continuous spectra. Proceedings of the National Academy of Sciences 18(3), 255–263 (1932) https://doi.org/10.1073/pnas.18.3.255
  * Koopman [1931] Koopman, B.O.: Hamiltonian systems and transformation in Hilbert space. Proceedings of the National Academy of Sciences 17(5), 315–318 (1931) https://doi.org/10.1073/pnas.17.5.315
  * Kamb et al. [2020] Kamb, M., Kaiser, E., Brunton, S.L., Kutz, J.N.: Time-delay observables for koopman: Theory and applications. SIAM Journal on Applied Dynamical Systems 19(2), 886–917 (2020) https://doi.org/10.1137/18M1216572
  * Mori [1965] Mori, H.: Transport, collective motion, and Brownian motion. Progress of theoretical physics 33(3), 423–455 (1965)
  * Zwanzig [1973] Zwanzig, R.: Nonlinear generalized Langevin equations. Journal of Statistical Physics 9(3), 215–220 (1973)
  * Zwanzig [2001] Zwanzig, R.: Nonequilibrium Statistical Mechanics. Oxford University Press, Oxford (2001)
  * Evans and Morriss [2008] Evans, D.J., Morriss, G.P.: Statistical Mechanics of Nonequilibrium Liquids, 2nd ed edn. Cambridge University Press, Cambridge (2008)
  * Hijón et al. [2010] Hijón, C., Español, P., Vanden-Eijndenc, E., Delgado-Buscalioni, R.: Mori–zwanzig formalism as a practical computational tool. Faraday Discussions 144 (2010)
  * Darve et al. [2009] Darve, E., Solomon, J., Kia, A.: Computing generalized Langevin equations and generalized Fokker-Planck equations. Proceedings of the National Academy of Sciences 106(27), 10884–10889 (2009) https://doi.org/10.1073/pnas.0902633106
  * Lin and Lu [2021] Lin, K.K., Lu, F.: Data-driven model reduction, wiener projections, and the koopman-mori-zwanzig formalism. Journal of Computational Physics 424, 109864 (2021) https://doi.org/10.1016/j.jcp.2020.109864
  * Li et al. [2015] Li, Z., Bian, X., Li, X., Karniadakis, G.E.: Incorporation of memory effects in coarse-grained modeling via the mori-zwanzig formalism. The Journal of Chemical Physics 143(24), 243128 (2015) https://doi.org/10.1063/1.4935490 https://doi.org/10.1063/1.4935490
  * Wang et al. [2019] Wang, S., Li, Z., Pan, W.: Implicit-solvent coarse-grained modeling for polymer solutions via mori-zwanzig formalism. Soft Matter 15, 7567–7582 (2019) https://doi.org/10.1039/C9SM01211G
  * Parish and Duraisamy [2017] Parish, E.J., Duraisamy, K.: Non-markovian closure models for large eddy simulations using the mori-zwanzig formalism. Phys. Rev. Fluids 2, 014604 (2017) https://doi.org/10.1103/PhysRevFluids.2.014604
  * Falkena et al. [2019] Falkena, S.K.J., Quinn, C., Sieber, J., Frank, J., Dijkstra, H.A.: Derivation of delay equation climate models using the Mori-Zwanzig formalism. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 475(2227), 20190075 (2019) https://doi.org/10.1098/rspa.2019.0075
  * Lin and Lu [2021] Lin, K.K., Lu, F.: Data-driven model reduction, Wiener projections, and the Koopman-Mori-Zwanzig formalism. Journal of Computational Physics 424, 109864 (2021) https://doi.org/10.1016/j.jcp.2020.109864
  * Lin et al. [2021] Lin, Y.T., Tian, Y., Livescu, D., Anghel, M.: Data-Driven Learning for the Mori–Zwanzig Formalism: A Generalization of the Koopman Learning Framework. SIAM Journal on Applied Dynamical Systems 20(4), 2558–2601 (2021) https://doi.org/10.1137/21M1401759
  * Mauroy and Mezic [2016] Mauroy, A., Mezic, I.: Global Stability Analysis Using the Eigenfunctions of the Koopman Operator. IEEE Transactions on Automatic Control 61(11), 3356–3369 (2016) https://doi.org/10.1109/TAC.2016.2518918
  * Wanner and Mezić [2020] Wanner, M., Mezić, I.: Robust approximation of the stochastic koopman operator. arXiv preprint arXiv:2011.00078 (2020)
  * Arbabi and Mezić [2017] Arbabi, H., Mezić, I.: Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the koopman operator. SIAM Journal on Applied Dynamical Systems 16(4), 2096–2126 (2017) https://doi.org/10.1137/17M1125236 https://doi.org/10.1137/17M1125236
  * Schmid [2010] Schmid, P.J.: Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 656, 5–28 (2010) https://doi.org/10.1017/S0022112010001217
  * Schmid et al. [2011] Schmid, P.J., Li, L., Juniper, M.P., Pust, O.: Applications of the dynamic mode decomposition. Theoretical and Computational Fluid Dynamics 25(1-4), 249–259 (2011)
  * Carleman [1932] Carleman, T.: Application de la théorie des équations intégrales linéaires aux systèmes d’équations différentielles non linéaires. Acta Math. 59, 63–87 (1932) https://doi.org/10.1007/BF02546499
  * Kowalski and Steeb [1991] Kowalski, K., Steeb, W.-H.: Nonlinear Dynamical Systems And Carleman Linearization. World Scientific, Singapore (1991)
  * [37] Lin, Y.-T., Lowrie, R.B., Aslangil, D., Subası, Y., Sornborger, A.T.: Koopman von neumann mechanics and the koopman representation: A perspective on solving nonlinear dynamical systems with quantum computers. arXiv:2202.02188
  * Wilf [2005] Wilf, H.S.: Generatingfunctionology, 3rd edn. Academic Press, Boston, MA (2005)
  * Chorin et al. [2000] Chorin, A.J., Hald, O.H., Kupferman, R.: Optimal prediction and the mori–zwanzig representation of irreversible processes. Proceedings of the National Academy of Sciences 97(7), 2968–2973 (2000) https://doi.org/10.1073/pnas.97.7.2968
  * Alexandre J. Chorin et al. [2002] Alexandre J. Chorin, Hald, O.H., Kupferman, R., Chorin, A.J., Hald, O.H., Kupferman, R.: Optimal prediction with memory. Physica D: Nonlinear Phenomena 166(3-4), 239–257 (2002) https://doi.org/10.1016/S0167-2789(02)00446-3
  * Forets and Pouly [2017] Forets, M., Pouly, A.: Explicit error bounds for carleman linearization. arXiv:1711.02552 (2017)