The Volterra Integrable case.
Novel analytical and numerical results.

M.Scaliaa, O. Ragniscob, B.Tirozzic, F.Zullod

aDipartimento di Matematica, Università degli Studi “La Sapienza”
P.le A.Moro 2, I-00185 Roma

bDipartimento di Matematica e Fisica, Università degli Studi “Roma TRE”
(retired)
Via della Vasca Navale 84,00146 Roma

c Dipartimento di Fisica, Università degli Studi “La Sapienza”

P-le Aldo Moro 2, I-00185 Roma

d DICATAM, Università degli Studi di Brescia

via Branze, 38 - 25123 Brescia, Italy &

INFN, Milano Bicocca

Piazza della Scienza 3, Milano, 20126, Italy

e-mail<EMAIL_ADDRESS>

###### Abstract

In the present paper we reconsider the integrable case of the Hamiltonian
$N$-species Volterra system, as it has been introduced by Vito Volterra in
1937 and significantly enrich the results already published in the ArXiv in
2019. In fact, we present a new approach to the construction of conserved
quantities and comment about the solutions of the equations of motion; we
display mostly new analytical and numerical results, starting from the
classical predator-prey model till the general $N$-species model.

KEYWORDS: Generalized Volterra Model, Hamiltonian formulation, Integrable
Systems

##  1 Preface

This paper is devoted to the memory of one of the authors, namely our beloved
colleague and friend Massimo Scalia, who sadly died on December 10, 2023,
after a dramatic car accident occurred at a crossroad not far from the centre
of Rome. Last summer, Massimo had been deeply affected by a major personal
tragedy, the death of his beloved companion (Adele Vannini) who passed away
after a long illness. He was able to react to this tremendous shock, fully
devoting himself to the two fields of interest that marked his whole life,
politics and science. We could even say that they were not two different
interests, inasmuch as, all along his career, he tirelessly struggled to
establish a bridge between them. The idea of resuming the brilliant approach
to Mathematical Biology introduced by Vito Volterra [2], recently revived by
Giorgio Israel in his beautiful monographs [3], has in fact characterized
several of his recent papers, among which we just quote the most recent one
[4]. From [2], Massimo took the key idea that conflicting variables, one
economical and the other ecological, could be thought as forming a pair of
conjugated variables (in the Hamiltonian language) in the spirit of the
predator-prey model [5]. The analogy could indeed be pushed forward to involve
a larger number of competing species, so deriving an N-species generalization
of the simplest two-species model. Once recast in a Hamiltonian form, it was
then natural to ask whether also such an enlarged model  could enjoy the
property of complete integrability. In fact we discovered that Volterra
himself had already addressed and solved this question, that amounted to
require a specific (but simple) structure for the interaction matrix. The
present article focusses on the study of the N-species Volterra model in the
integrable case, and as such is a reformulation and a significant extension of
the preprints appeared in the ArXiv [6], whose contents are here largely
incorporated.

##  2 Introduction

The current twofold crisis, at a global level, of Economics and Environment
has been deeply investigated and reported by many authors, some of which have
deeply criticized the deafness of Economists with respect to the environmental
crisis, that at most has been assessed for its dramatic consequences on GDP
(see the well-known “Stern’s Report” [7]). The double crisis and its
entanglement would request models, even before than a theory, able to put
together economic and environmental variables in order to build a global
stationary state to rule present predicament in the perspective of a
sustainable scenario. The latter theme is not a news, several attempts having
been realized starting from the seventies towards a definition of “steady
state”([8], [9], [10]) but it could be useful to face the problem with other
scientific tools, as it has been recently proposed in [11], [12], [13], where
the leading idea is to put together pairs, each constituted by an economic
variable and an environmental one, that present a behavior “predator-prey”
type, as it is suggested by some of the most important pairs one can select
for the model. Already fifty years ago, a similar idea was applied to build
the Goodwin model, but with a pair of variables only economic to describe an
economic cycle [14], the so called “class struggle” model, where a variable
linked to the wage rate assumes the role of predator and the one giving
employment rate is the prey.

##  3 The Model

Here, aiming at giving a self-contained description of our approach, we resume
what two of the authors wrote in [6]. As is well known, the original idea by
Vito Volterra [5] was that of determining the evolution of a two species
biological system, the so-called “predator-prey” model, answering a question
raised by his son in law, the biologist Umberto d’Ancona [15], who was
wondering why the total catch of selachians (mostly sharks) was considerably
raising during World War 1, with respect to other more desirable kind of
fishes, in correspondence with the decrease of fishing activity [16]. To
answer that question, Vito Volterra constructed a dynamical system that
enabled him to identify the essential features of what was going on,
elucidating the properties entailing the existence of a stable equilibrium
configuration (and of periodic orbits in its neighborhood), and unveiling the
asymptotic behavior of the system under general initial conditions. He quickly
realized that the predator-prey model was just the simplest example in a large
class of biological, or rather ecological systems with pairwise interaction.
He was soon interested in understanding the mathematical properties of the $N$
species pairwise interacting model, and devoted a considerable effort to find
suitable Lagrangian and Hamiltonian formulations, with the final aim of
achieving a description where the deep analogy with the well established
theory of mechanical systems stemming from the Maupertuis minimal action
principle be made transparent. We would say that not the whole “Biological-
Mechanical” dictionary that he proposed in his famous paper (dating back to
1937), Principes de Biologie Mathematique [2], resisted the future
developments of both disciplines, and some of the notions he tried to
introduce look nowadays a bit artificial. But we believe that the core of his
derivation is still alive, as it has been witnessed by very wide spread
applications over about a century in many scientific research subjects, such
as Populations demography, Bio-physics, Biomedicine, Ecology and also
Economics. We notice that in [2] his main aim was the formulation of this
generalized model in a Hamiltonian language, with the purpose of elucidating
the algebraic conditions leading to a completely integrable model. So we think
that it could be worth recalling the key ingredients of Volterra’s approach,
and even emphasizing the role that the very special case of completely
integrable dynamics could play in the search for Stationary State models in
the economic-ecological framework. His general Predator-Prey model reads:

$\frac{dN_{r}}{dt}=\epsilon_{r}N_{r}+\sum_{s\neq r=1}^{N}A_{rs}N_{r}N_{s}$ (1)

In (1), we have set all the parameters introduced in [5] $\beta_{r}=1\forall
r$; $\epsilon_{r}$ are the natural growth coefficients of each species and
$A_{rs}$ are interaction coefficients between species $r$ and species $s$ that
account for the probability of encountering between two individuals. Moreover,
if we assume that the matrix $A$, whose elements are $A_{rs}$, is nonsingular,
then the system of equations defining equilibrium configurations (other than
the trivial one $(\\{N_{r}\\}=\\{0\\})$), namely

$0=\epsilon_{r}+\sum_{s=1}^{N}A_{rs}N_{s}^{(0)}$ (2)

has a unique solution, say $N_{r}^{(0)}$, $r=1,\cdots,N$. If, in addition,
according to [5], we require $A$ to be skew-symmetric and $N$ to be even, then
the eigenvalues of $A$ will be purely imaginary and complex conjugate in
pairs. What is more important, however, is that the roots of (2) be all
positive. As a necessary condition, some of “the spontaneous growth”
coefficients $\epsilon_{j}$ have to be negative. Of course, being skew
symmetric, the matrix $A$ will be singular in the case of an $odd$ number of
species.

###  3.1 Lagrangian and Hamiltonian formulation

As many other researchers of his time, Volterra was feeling more assured if a
phenomenon quantified by Mathematics could find an analogue with Mechanics,
that moreover allowed resorting to the powerful formalism and theorems of the
latter. The quantity of life of the species $r$, defined as
$X_{r}=\int_{0}^{t}N_{r}(t^{\prime})dt^{\prime}$ suggested to Volterra the
introduction of a biological, or rather $ecological$, Lagrangian $\Phi$,
defined as:

$\Phi=\sum_{r}\epsilon_{r}X_{r}+\sum_{r}X_{r}^{\prime}\log
X_{r}^{\prime}+\frac{1}{2}\sum_{rs}A_{rs}X_{r}^{\prime}X_{s}$ (3)

In terms of (3), (1) can be written as Euler-Lagrange equations

$\frac{d}{dt}\frac{\partial\Phi}{\partial
X_{r}^{\prime}}-\frac{\partial\Phi}{\partial X_{r}}=0$ (4)

yielding the ODEs

$X_{r}^{{}^{\prime\prime}}=(\epsilon_{r}+\sum_{s}A_{sr}X_{s}^{\prime})X_{r}^{\prime}$
(5)

which are just (1), up to the substitution $N_{r}=X_{r}^{\prime}$, the
$\prime$ denoting time-derivative. In the following we will change the
notations, by setting:

$X_{r}:=q_{r};~{}~{}X_{r}^{\prime}:=\dot{q}_{r}$ (6)

In other words, we turn to the standard notation of Lagrangian mechanics

###  3.2 From Lagrange to Hamilton

The linear momenta, canonically conjugated to the quantities of life, are
defined as

$p_{r}=\frac{\partial\Phi}{\partial\dot{q}_{r}}=\log\dot{q}_{r}-\frac{1}{2}\sum_{s}A_{rs}q_{s}$
(7)

whence

$\dot{q}_{r}=\exp(p_{r}+\frac{1}{2}\sum_{s}A_{rs}q_{s})$ (8)

Note that, in terms of the original ecological variables the expression
$\sum_{r}\epsilon_{r}q_{r}-\dot{q}_{r}$ takes the form:

$\sum_{r}\epsilon_{r}q_{r}-\dot{q}_{r}=\sum_{r}\epsilon_{r}\int_{0}^{t}dt^{\prime}N_{r}(t^{\prime})-N_{r}$
(9)

Through a transformation of Legendre type we define the Hamiltonian

${\mathcal{H}}=\Phi-\sum_{r}\dot{q}_{r}p_{r}=\sum_{r}\epsilon_{r}q_{r}-\dot{q}_{r}=\sum_{r}\epsilon_{r}q_{r}-\exp(p_{r}+\frac{1}{2}\sum_{s}A_{rs}q_{s})$
(10)

Volterra [2] showed that (1) can be written in the standard hamiltonian form

$\dot{q}_{r}=\frac{\partial{\mathcal{H}}}{\partial p_{r}}$ (11)
$\dot{p}_{r}=-\frac{\partial{\mathcal{H}}}{\partial q_{r}}$ (12)

It is easily seen (see again [4]) that the Hamiltonian system (10),(11)(12)
has the following $N$ independent non autonomous integrals of motion:

${\mathcal{H}_{r}}=\frac{p_{r}-\frac{1}{2}\sum_{s}A_{rs}q_{s}}{\epsilon_{r}}-t~{}~{}r=1,\cdots,N.$
(13)

whence one can select $N-1$ autonomous integrals by taking for instance
${\mathcal{H}_{1,r}}\equiv{\mathcal{H}_{r}}-{\mathcal{H}_{1}}$, and have a
complete set by adding the Volterra $N-$ species Hamiltonian (10).

A more modern approach to the Hamiltonian structure underlying the generalized
Volterra system can be found in an elegant paper by R.Loja Fernandes and
W.Oliva [17], where a Poisson morphism is established between the original
system, living in ${\mathbb{R}}^{N}$ and equipped with a quadratic Poisson
structure, and the one recast, after Volterra, in a canonical Hamiltonian form
and thus living in ${\mathbb{R}}^{2N}$. However, we remark that in [17] the
issue of complete integrability is not addressed.

On the contrary , in our opinion, the question whether there exists a special
form of the matrix elements $A_{rs}$ entailing involutivity of the complete
set of integrals of motion $({\mathcal{H}},{\mathcal{H}_{1,r}})$ is a relevant
one to ask. It turns out that this form has been found by Volterra himself [2]
and is the following:

$A_{rs}=\epsilon_{r}\epsilon_{s}(B_{r}-B_{s})~{}~{}~{}r,s=1,\cdots N$ (14)

where $N$ is the number of competing populations and the $B_{r}$ are real (and
positive) numbers. Clearly, (14) can be cast in the compact form:

$A=\\{A_{rs)}\\}=[B,\epsilon\otimes\epsilon]$ (15)

where $B$ $=diag(B_{1},\cdots,B_{N})$, and $\epsilon$ is the vector
$(\epsilon_{1},\cdots,\epsilon_{N})^{t}$, meaning that $A$ is the commutator
of a diagonal matrix with distinct entries and a rank one matrix

###  3.3 A note about the equilibrium conditions

So, if we require integrability, (15) shows that the invertibility of $A$ has
to be given up for $N>2$. Indeed, in the integrable case $Ker(A)$ has
codimension 2, and correspondingly its Range $R(A)$ is two-dimensional.
Accordingly, the equilibrium conditions (2) read:

$\sum_{s=1}^{N}\epsilon_{s}(B_{s}-B_{r})N^{0}_{s}=1$

where we have denoted by $N^{0}_{s}(s=1,\cdots,N)$ the equilibrium population
numbers, implying that the equilibrium configuration be defined as the
intersection of the two hyperplanes

$\displaystyle\sum_{s}\epsilon_{s}N^{0}_{s}=0$ (16)
$\displaystyle\sum_{s}\epsilon_{s}B_{s}N^{0}_{s}=1$ (17)

Consequently an $equilibrium$ $configuration$ can exist only if some of the
$\epsilon_{s}$ are negative. Moreover, a $unique$ equilibrium position exists
only for $N=2$; setting $\mu=B_{2}-B_{1}$, we get:

$N^{0}_{1}=-\frac{1}{\mu\epsilon_{1}};~{}~{}~{}~{}N^{0}_{2}=\frac{1}{\mu\epsilon_{2}}$
(18)

For instance, in the case $N=3$, we have a one parameter family of equilibrium
solutions, reading ($0<\rho<1$)

$\epsilon_{1}N^{0}_{1}=\frac{\rho}{B_{1}-B_{3}};\hskip
2.84544pt\epsilon_{2}N^{0}_{2}=\frac{1-\rho}{B_{2}-B_{3}}\hskip
2.84544pt;\epsilon_{3}N^{0}_{3}=-(\epsilon_{1}N^{0}_{1}+\epsilon_{2}N^{0}_{2})$
(19)

In (19), in order the equilibrium species populations be positive, we have to
require $(B_{1}-B_{2})\epsilon_{1}$, and $(B_{2}-B_{3})\epsilon_{2}$ to be
positive quantities, while the ratios $\frac{\epsilon_{1}}{\epsilon_{3}}$ and
$\frac{\epsilon_{2}}{\epsilon_{3}}$ have to be negative. We end the present
subsection by remarking that in all cases, whether they are integrable or not,
the N-species Volterra system enjoy a sort of box structure, being equipped
with a number of invariant submanifolds, obtaining when only a (possibly not
ordered) subset of species is alive, this fact being dictated just by the
initial conditions. So, for the two species case, the axes $N_{1}=0$ and
$N_{2}=0$ are invariant submanifolds, for the 3-species case we have the 6
invariant submanifolds given by the axes and by the planes $N_{i}=0$, and so
forth. In particular the 3-species case gives rise to three Volterra-Lotka
systems.

##  4 More on the integrable case. The Volterra approach

However, once realized that integrability implies the non-uniqueness of the
equilibrium configuration ($\forall$ $N>2$), what is needed to stress is that,
for the special form of the matrix $A$ given by (14), the integrals of motion
are still functionally independent. This is readily seen as (13) shows that
the linear dependence of those integrals upon the momenta $p_{r}$ is in no-way
affected by by the specific form of the matrix $A$ (while the requirement that
the $B_{r}$ be all distinct is mandatory!), so that the rank of the Jacobian
matrix constructed with the gradients of the integrals of motion with respect
to the canonical coordinates is maximal (namely $N$) whatever be that form.
So, the integrable version of the $N$-species Volterra system is again a
genuine hamiltonian system with $N$ degrees of freedom.

Here we write down explicitly the expression of the Hamiltonian and of the
integrals of motion in the integrable case, resuming what we sketched in
formulas (13),( 14),(15). Denoting now by ${\mathcal{H}}_{int}$ the
Hamiltonian (10) we have:

${\mathcal{H}}_{int}=\sum_{r=1}^{N}\epsilon_{r}q_{r}-\exp[p_{r}+(\epsilon_{r}/2)\sum_{s=1}^{N}\epsilon_{s}(B_{r}-B_{s})q_{s}]$
(20)

and

${\mathcal{H}}_{r}=p_{r}/\epsilon_{r}-(1/2)\sum_{s=1}^{N}\epsilon_{s}(B_{r}-B_{s})q_{s}-t\quad
r=1,\cdots,N$ (21)

so that

${\mathcal{H}}_{rl}\equiv{\mathcal{H}}_{r}-{\mathcal{H}}_{l}=p_{r}/\epsilon_{r}-p_{l}/\epsilon_{l}-\frac{1}{2}(B_{r}-B_{l})\sum_{s=1}^{N}\epsilon_{s}q_{s},\quad
s=1,\cdots,N$ (22)

The constants of motion (22) are mutually in involution. So we can take for
instance $l=1$ and get $N-1$ independent integrals of motion in involution.
The set can be completed by adding any function of the Hamiltonian. for
instance the Hamiltonian itself.

The above formulas clearly show that, in the integrable case, both the
Volterra Hamiltonian and the involutive family of integrals of motion depend
on the quantities of life only through the inner products
$<{\vec{\epsilon}}|Q>$ and $<{\vec{\epsilon}}|B|Q>$, where $|Q>$ is the vector
of components $q_{j}$, while by $|B|Q>$ we have denoted the vector of
components $B_{j}q_{j}$. However, even in the completely integrable case, we
did not succeed in reducing our problem to quadratures for $N$ larger than
$2$, although this possibility is a well known result in Classical Mechanics
([18],[19]).

Of course, if one comes back to the Lagrangian formulation, one will get
expressions involving the variables $q_{s}$ and $\dot{q}_{s}$. Let us also
notice that

$\exp(\epsilon_{r}{\mathcal{H}}_{r})=\exp[p_{r}-\epsilon_{r}/2\sum_{s=1}^{N}\epsilon_{s}(B_{r}-B_{s})q_{s}-\epsilon_{r}t]$
(23)

Obviously, one could choose $\exp({\mathcal{H}_{rl}})$ as an alternative
legitimate form for an involutive family of integrals of motion. In the
simplest nontrivial case, $N=2$, (20) reads (the subscript $V$ refers to
Volterra):

${\mathcal{H}}_{V}=\epsilon_{1}q_{1}+\epsilon_{2}q_{2}-\exp[p_{1}+(1/2)\epsilon_{1}\epsilon_{2}(B_{1}-B_{2})q_{2}]-\exp[p_{2}-(1/2)\epsilon_{1}\epsilon_{2}(B_{1}-B_{2})q_{1}]$
(24)

The above formula can be slightly simplified through the canonical
transformation (in fact, just a rescaling):

$p_{j}\to\tilde{p}_{j}=p_{j}/\epsilon_{j};\quad
q_{j}\to\tilde{q}_{j}=\epsilon_{j}q_{j}$ (25)

that maps (22) into:

${\mathcal{H}}_{rl}=\tilde{p}_{r}-\tilde{p}_{l}+(1/2)(B_{r}-B_{l})\sum_{j}\tilde{q}_{j}$
(26)

and (24) into:

${\mathcal{H}}_{V}=\tilde{q}_{1}+\tilde{q}_{2}-\exp[\epsilon_{1}(\tilde{p}_{1}+1/2)(B_{1}-B_{2})\tilde{q}_{2})]-\exp[\epsilon_{2}(\tilde{p}_{2}-(1/2)(B_{1}-B_{2})\tilde{q}_{1}]$
(27)

It might be worth to notice that in the integrable case the $q$-dependence of
the Hamiltonian (20) materializes just through the $global$ variables
$Q:=\sum_{s}\tilde{q}_{s}$ and $BQ=\sum_{s}B_{s}\tilde{q}_{s}$. Of course, the
same $Mean$ $Field$ property is shared, even more perspicuously, by the
dynamics expressed in term of the “ecological” variables $N_{r}$.

###  4.1 Integrability of the N=2 case via the canonical formalism

Let us slightly simplify the notations in (27), by setting $\mu:=B_{2}-B_{1}$,
and introducing the new canonical variables:

$P_{1}=\frac{1}{\sqrt{\mu}}(\tilde{p}_{1}+\frac{1}{2}\mu\tilde{q}_{2});~{}~{}~{}Q_{1}=\frac{1}{\sqrt{\mu}}(-\tilde{p}_{2}+\mu\frac{1}{2}\tilde{q}_{1})$
(28)

$P_{2}=\frac{1}{\sqrt{\mu}}(\tilde{p}_{1}-\mu\frac{1}{2}\tilde{q}_{2});~{}~{}~{}Q_{2}=\frac{1}{\sqrt{\mu}}(\tilde{p}_{2}+\
\frac{1}{2}\mu\tilde{q}_{1})$ (29)

In terms of these new variables, the first integral:

${\mathcal{H}}_{12}=\tilde{p}_{1}-\tilde{p}_{2}-(\mu/2)(\tilde{q}_{1}+\tilde{q}_{2})$
(30)

takes the form

${\mathcal{H}}_{12}=\frac{1}{\sqrt{\mu}}(P_{2}-Q_{2})$ (31)

while the two-particle Hamiltonian reads:

${\mathcal{H}}_{V}=\frac{1}{\sqrt{\mu}}[Q_{1}+Q_{2}+P_{1}-P_{2}]-\exp{\sqrt{\mu}}(\epsilon_{1}P_{1})-\exp[-{\sqrt{\mu}}(\epsilon_{2}Q_{1})]$
(32)

Inserting the first integral (31), on the level surface
${\mathcal{H}}_{12}=C$, up to an irrelevant additive constant we can finally
write:

${\mathcal{H}}_{V}=\frac{1}{\sqrt{\mu}}(Q_{1}+P_{1})-\exp{\sqrt{\mu}}(\epsilon_{1}P_{1})-\exp[-{\sqrt{\mu}}(\epsilon_{2}Q_{1})]$
(33)

It follows that, in terms of these new coordinates, the above Hamiltonian is a
“one-body” hamiltonian (integrable by definition), which is nothing but the
traditional Volterra-Lotka Hamiltonian. One may pervene to more elegant
formulas by defining:

$X\equiv\exp{\sqrt{\mu}}(\epsilon_{1}P_{1});~{}~{}~{}Y\equiv\exp[-{\sqrt{\mu}}(\epsilon_{2}Q_{1})]$

leading to:

$\displaystyle\dot{X}=\epsilon_{2}X(1-\mu\epsilon_{1}Y)$ (34)
$\displaystyle\dot{Y}=\epsilon_{1}Y(1+\mu\epsilon_{2}X)$ (35)

The (nontrivial and stable) equilibrium position is the pair
$(X^{(0)}=-\frac{1}{\mu\epsilon_{2}},Y^{(0)}=\frac{1}{\mu\epsilon_{1}})$
whence it follows that it can belong to the first quadrant only if the
coefficients of spontaneous growth have opposite sign, as it is natural if we
require that the predator species can survive only by ”eating” the prey one.
As is well known ([2],[16]), the equations for the orbits of (34,35) can be
written in closed form:

$(X\exp(-\mu\epsilon_{2}X))^{\epsilon_{1}}~{}=~{}K(Y\exp(-\mu\epsilon_{1}Y))^{\epsilon_{2}}$
(36)

where $K$ is a positive constant.

###  4.2 Volterra’s approach

Here we briefly sketch the procedure followed by Volterra to integrate the
N-species system in terms of the natural coordinates.

Volterra defines the quantities ${\mathcal{N}}:=\sum_{s}\epsilon_{s}N_{s}$ and
${\mathcal{M}}:=1-\sum_{s}\epsilon_{s}B_{s}N_{s}$, then rewriting the original
dynamical system (66) as:

$\dot{N}_{r}=\epsilon_{r}N_{r}(1+\sum_{s}(B_{r}-B_{s})\epsilon_{s}N_{s})$

or, in other terms:

$\dot{N}_{r}=\epsilon_{r}N_{r}(B_{r}{\mathcal{N}}+{\mathcal{M}})$

namely:

$(1/\epsilon_{r})d/dt\log N_{r}=(B_{r}{\mathcal{N}}+{\mathcal{M}})$

Note that the previous equation implies
$\sum_{r}(\dot{N}_{r}-\epsilon_{r}N_{r})=0$, which is just the conservation of
the hamiltonian (10). It is evident that the involutivity constraints on the
coefficients $A_{rs}$ entail a typical ${\it Mean~{}Field}$ dynamics. Each
species interacts with the others through the ${\it collective~{}variables}$
${\mathcal{N}}$ and ${\mathcal{M}}$. By taking two different values of the
index $r$ and subtracting, the variable ${\mathcal{M}}$ can be eliminated,
resorting to:

$\frac{(1/\epsilon_{r})d/dt\log N_{r}-(1/\epsilon_{s})d/dt\log
N_{s}}{B_{r}-B_{s}}={\mathcal{N}}$

By introducing the variables

$Y_{k}\equiv{(1/\epsilon_{k}})\log N_{k}$ (37)

we immediately get $N-2$ integrals of motion written in terms of the natural
variables.

We arrive at the linear formula

$\frac{1}{B_{r}-B_{j}}(Y_{r}-Y_{j})-\frac{1}{B_{s}-B_{j}}(Y_{s}-Y_{j})=C_{rs}$

The simplest nontrivial case is the $3-species$ case, where an elementary
calculation shows that the three equations above are indeed the same,
yielding:

$(B_{2}-B_{3})Y_{1}+(B_{3}-B_{1})Y_{2}+(B_{1}-B_{2})Y_{3}=const.$ (38)

We remind the form of the equations of motion for the integrable three species
case:

$\displaystyle\dot{N}_{1}=\epsilon_{1}N_{1}+\epsilon_{1}\epsilon_{2}\
(B_{1}-B_{2})N_{1}N_{2}+\epsilon_{1}\epsilon_{3}(B_{1}-B_{3})N_{1}N_{3}$ (39)
$\displaystyle\dot{N}_{2}=\epsilon_{2}N_{2}+\epsilon_{2}\epsilon_{1}(B_{2}-B_{1})N_{2}N_{1}+\epsilon_{2}\epsilon_{3}(B_{2}-B_{3})N_{2}N_{3}$
(40)
$\displaystyle\dot{N}_{3}=\epsilon_{3}N_{3}+\epsilon_{3}\epsilon_{1}(B_{3}-B_{1})N_{3}N_{1}+\epsilon_{3}\epsilon_{2}(B_{3}-B_{2})N_{3}N_{2}$
(41)

and add the expression of the integrals of motion for the system (39, 40, 41)
in terms of the numerosities and of the quantities of life. Taking into
account (37), formula ($\Ref{Thirdintegral}$) takes the form

$N_{1}^{(B_{2}-B_{3})/\epsilon_{1}}N_{2}^{(B_{3}-B_{1})/\epsilon_{2}}N_{3}^{(B_{1}-B_{2})/\epsilon_{3}}=I_{123}$
(42)

Of course, (42) has a meaning only for nonzero initial data, where it can be
written as well as:

$(N_{1}/N_{1}{(0)})^{(B_{2}-B_{3})/\epsilon_{1}}(N_{2}/N_{2}{(0)})^{(B_{3}-B_{1})/\epsilon_{2}}(N_{3}/N_{3}{(0)})^{(B_{1}-B_{2})/\epsilon_{3}}=1$
(43)

The Hamiltonian ${\mathcal{H}}$ is written as:

${\mathcal{H}}=\sum_{r=1}^{N}\epsilon_{r}\int_{0}^{t}N_{r}(t^{\prime})dt^{\prime}-N_{r}\
$ (44)

Implying that its value equals minus the total population at $t=0$. The
involutive integrals $\exp({\mathcal{H}_{r}})$ read:

$\exp({\mathcal{H}_{r}})=N_{r}\exp[-\sum_{s}(B_{r}-B_{s})\epsilon_{r}\epsilon_{s}\int_{0}^{t}N_{s}(t^{\prime})dt^{\prime}]$
(45)

so each of them equals the corresponding initial population. In the next
section, devoted to both analytical and numerical results, we will focus
attention on the $three$ $species$ integrable case, which escapes the general
analysis presented by Volterra in( [2]), focussed on the study of an even
number of species, and always assuming invertibility of the interaction matrix
$A$.

##  5 Numerical and analytical results: the surprising integrable three
species case and beyond

In this section, and in the following subsection, we propose a way to
construct the integrals of motion for the model under scrutiny which is
alternative to the one followed by Volterra, recalled in the previous section
and in ([6]). We will as well integrate numerically the equations of motion
with the help of $Mathematica$ package, namely the (39,40),41) equations, for
different choices of the relevant parameters. To simplify a little bit the
notations, we have set:

$B_{1}-B_{2}=\alpha;~{}~{}B_{2}-B_{3}=\beta$ (46)

whence $B_{1}-B_{3}=\alpha+\beta$. By starting from the equations of motion

$\displaystyle\dot{N}_{1}=\epsilon_{1}N_{1}+\epsilon_{1}\epsilon_{2}\alpha
N_{1}N_{2}+\epsilon_{1}\epsilon_{3}(\alpha+\beta)N_{1}N_{3}$ (47)
$\displaystyle\dot{N}_{2}=\epsilon_{2}N_{2}-\epsilon_{2}\epsilon_{1})\alpha
N_{2}N_{1}+\epsilon_{2}\epsilon_{3}\beta N_{2}N_{3}$ (48)
$\displaystyle\dot{N}_{3}=\epsilon_{3}N_{3}-\epsilon_{3}\epsilon_{1}(\alpha+\beta)N_{3}N_{1}-\epsilon_{3}\epsilon_{2}\beta
N_{3}N_{2}\ $ (49)

Let us look for an integral of motion written as

$e^{A(N_{1}(t)+N_{2}(t)+N_{3}(t))}N_{1}(t)^{m}N_{2}(t)^{n}N_{3}(t)^{k},$ (50)

for suitable constants $A,m,n,k$. By deriving (50) with respect to time,
dividing by (50) itself and making use of (39-41) we find

$\begin{split}&\epsilon_{1}N_{1}\left(A-n\alpha\epsilon_{2}-k(\alpha+\beta)\epsilon_{3}\right)+\epsilon_{2}N_{2}\left(A+m\alpha\epsilon_{1}-k\beta\epsilon_{3}\right)+\\\
&+\epsilon_{3}N_{3}\left(A+n\beta\epsilon_{2}+m(\alpha+\beta)\epsilon_{1}\right)+\epsilon_{1}m+\epsilon_{2}n+\epsilon_{3}k=0.\end{split}$
(51)

The coefficients of the $N_{i}^{\prime}s$, $i=1,2,3$, are all compatible each
other if $\epsilon_{1}m+\epsilon_{2}n+\epsilon_{3}k=0$, giving e.g.
$A=k\beta\epsilon_{3}-m\alpha\epsilon_{1}$. So, the following quantity

$e^{\left(k\beta\epsilon_{3}-m\alpha\epsilon_{1}\right)\left(N_{1}(t)+N_{2}(t)+N_{3}(t)\right)}N_{1}(t)^{m}N_{2}(t)^{n}N_{3}(t)^{k}=I_{m,n,k}$
(52)

is a conserved quantity if the parameters $(m,n,k)$ satisfies
$\epsilon_{1}m+\epsilon_{2}n+\epsilon_{3}k=0$. Notice that (42) can be
rewritten as:

$\frac{N_{1}^{\beta/\epsilon_{1}}N_{3}^{\alpha/\epsilon_{3}}}{N_{2}^{(\alpha+\beta)/\epsilon_{2}}}=I_{123},$
(53)

and this conserved quantity corresponds to the choice
$(m,n,k)=(\beta/\epsilon_{1},-(\alpha+\beta)/\epsilon_{2},\alpha/\epsilon_{3})$.

Of the three parameters $(m,n,k)$, only two are independent because of the
relation $\epsilon_{1}m+\epsilon_{2}n+\epsilon_{3}k=0$. It is possible to
write two different surfaces by choosing properly the values of the constants
$m,n,k$: if these surfaces intersect by defining a closed curve, the
corresponding motion defined by (39-41) will be periodic. Let us make an
example. If we take

$\alpha=1,\;\beta=2\;\epsilon_{1}=\epsilon_{2}=1,\;\epsilon_{3}=-1.$

Then, two conserved quantities are given by

$\frac{N_{1}^{2}}{N_{2}^{3}N_{3}}=I_{1},\quad\quad
e^{-2(N_{1}+N_{2}+N_{3})}N_{2}N_{3}=I_{2}.$ (54)

We give the plot of the closed orbit $(N_{1}(t),N_{2}(t),N_{3}(t))$
numerically obtained by taking the initial conditions
$N_{1}(0)=N_{2}(0)=N_{3}(0)=1$, the plot of $I_{1}$, of $I_{2}$ and finally
the three plots together: as it can be seen from the last figure the two
surfaces $I_{1}$ and $I_{2}$ intersect in a closed curve and the motion is
constrained on this curve.

Figure 1: Plot of the closed orbit $(N_{1}(t),N_{2}(t),N_{3}(t))$
corresponding to the initial conditions $N_{1}(0)=N_{2}(0)=N_{3}(0)=1$ Figure
2: Plot of the surface $I_{1}$ in (54) corresponding to
$N_{3}=\frac{N_{1}^{2}}{N_{2}^{3}}$. Figure 3: Plot of the surfaces $I_{2}$ in
(54) corresponding to $e^{-2(N_{1}+N_{2}+N_{3})}N_{2}N_{3}=e^{-6}$. Figure 4:
Plot of $N_{3}=\frac{N_{1}^{2}}{N_{2}^{3}}$ and
$e^{-2(N_{1}+N_{2}+N_{3})}N_{2}N_{3}=e^{-6}$ with the numerically obtained
closed orbit (in black)

Together with the geometrical picture given by the orbits, we have the
dynamical picture associated with the time-behavior of the system.

We will consider a few different cases. Case 1a:
$\epsilon_{1}=\epsilon_{2}=1;~{}~{}\epsilon_{3}=-1;\alpha=0.5;\beta=-2.5$

$N_{1}{(0)}=1,N_{2}{(0)}=1,N_{3}{(0)}=1,t\in[0,20]$
$\displaystyle\dot{N}_{1}=N_{1}+0.5N_{1}N_{2}-3N_{1}N_{3},$ (55)
$\displaystyle\dot{N}_{2}=N_{2}-0.5N_{1}N_{2}-2.5N_{2}N_{3}$ (56)
$\displaystyle\dot{N}_{3}=-N_{3}+3N_{1}N_{3}+2.5N_{3}N_{2}$ (57)

Case
2a.$\epsilon_{1}=\epsilon_{2}=1;~{}~{}\epsilon_{3}=-1,\alpha=-0.5,\beta=-0.5$

$N_{1}{(0)}=1,N_{2}{(0)}=1/2,N_{3}{(0)}=1/2,t\in[0,20]$

$\dot{N}_{1}=N_{1}-0.5N_{1}N_{2}-N_{1}N_{3}$
$\dot{N}_{2}=N_{2}+0.5N_{1}N_{2}-0.5N_{2}N_{3}$
$\dot{N}_{3}=-N_{3}+N_{1}N_{3}+0.5N_{3}N_{2}$

The shape of oscillation is somehow sensible to the initial conditions, we
show a few interesting cases

Figure 5: Oscillations of $N_{3}$ with initial conditions $111$,
$\alpha=-0.5$, $\beta=-0.5$ Figure 6: Oscillations of the $N_{2}$ with initial
conditions $111$,$\alpha=-0.5$, $\beta=-0.5$

[H]

Figure 7: Oscillations of the $N_{1}$ with initial conditions
$111$,$\alpha=-0.5$, $\beta=-0.5$

Further interesting case:

$\alpha=0.5,\beta=2.5$

[H]

Figure 8: Oscillations of $N_{2}$ with initial conditions
$N_{1}(0)=1$,$N_{2}(0)=1/2$, $N_{3}(0)=1/2$,$\alpha=0.5$, $\beta=2.5$ Figure
9: Oscillations of $N_{1}$,$N_{3}$ species with initial conditions
$N_{1}(0)=1$,$N_{2}(0)=0$, $N_{3}(0)=1$, $N_{2}$ is zero, $\alpha=0.5$,
$\beta=2.5$

To conclude this subsection it sounds appropriate comparing the fully
nonlinear picture above described with the one corresponding to the linearized
evolution in the neighbourhood of (one suitable point) of the equilibrium
configuration. To this aim, Let us set

$N_{i}(t)=N_{i}^{0}+x_{i}(t),\quad i=1,2,3$ (59)

where $x_{i}$ are small quantities. At first order in the $x_{i}$, the set of
equations (1) becomes

$\begin{split}&\dot{x}_{1}=\frac{B_{1}-B_{2}}{B_{1}-B_{3}}\rho\epsilon_{2}x_{2}+\rho\epsilon_{3}x_{3}\\\
&\dot{x}_{2}=\frac{B_{2}-B_{1}}{B_{2}-B_{3}}(1-\rho)\epsilon_{1}x_{1}+(1-\rho)\epsilon_{3}x_{3}\\\
&\dot{x}_{3}=\left(\frac{B_{1}-B_{3}}{B_{2}-B_{3}}+\rho\frac{B_{2}-B_{1}}{B_{2}-B_{3}}\right)\epsilon_{1}x_{1}+\left(1-\rho\frac{B_{1}-B_{2}}{B_{1}-B_{3}}\right)\epsilon_{2}x_{2}\end{split}$
(60)

We need the eigenvalues of the matrix of the coefficients, given by

$M=\begin{pmatrix}0&\frac{B_{1}-B_{2}}{B_{1}-B_{3}}\rho\epsilon_{2}&\rho\epsilon_{3}\\\
\frac{B_{2}-B_{1}}{B_{2}-B_{3}}(1-\rho)\epsilon_{1}&0&(1-\rho)\epsilon_{3}\\\
\left(\frac{B_{1}-B_{3}}{B_{2}-B_{3}}+\rho\frac{B_{2}-B_{1}}{B_{2}-B_{3}}\right)\epsilon_{1}&\left(1-\rho\frac{B_{1}-B_{2}}{B_{1}-B_{3}}\right)\epsilon_{2}&0\end{pmatrix}$
(61)

The characteristic polynomial is given by:

$\lambda^{3}-\lambda^{2}\textrm{Tr}(M)+\frac{\lambda}{2}\left(\textrm{Tr}(M)^{2}-{Tr}(M^{2})\right)-\textrm{De}t(M)=0,$
(62)

and, since $\textrm{Tr}(M)=\textrm{Det}(M)=0$ it reduces to

$\lambda(\lambda^{2}-\frac{1}{2}\textrm{Tr}(M^{2}))=0.$ (63)

We see that an eigenvalue is always zero, the sum of the other two is zero. So
the three eigenvalues are $(0,\lambda_{1},-\lambda_{1})$. If $\lambda_{1}$ is
real, then the fixed point is unstable. If $\lambda_{1}$ is imaginary (it
happens if $\textrm{Tr}(M^{2})<0$) then the matrix $M$ must have three
independent eigenvectors in order the point to be stable [16]. The eigenvector
corresponding to the eigenvalue $0$ is
$(\frac{B_{2}-B_{3}}{\epsilon_{1}},\frac{B_{3}-B_{1}}{\epsilon_{2}},\frac{B_{1}-B_{2}}{\epsilon_{3}})$
whereas the other two must be independent. The point is not asymptotically
stable obviously. So, if $\textrm{Tr}(M^{2})<0$ one has periodic orbits, the
period being $\textrm{T}=\frac{2\sqrt{2}\pi}{\sqrt{|\textrm{Tr}(M^{2})|}}$. By
setting $B_{1}=B_{2}+\alpha$ and $B_{2}=B_{3}+\beta$, we get:

$\textrm{Tr}(M^{2})=2\epsilon_{3}\left(\epsilon_{2}+\rho(\epsilon_{1}-\epsilon_{2})\right)-2\alpha\rho(1-\rho)\frac{\epsilon_{1}(\epsilon_{2}-\epsilon_{3})\alpha+\epsilon_{3}(\epsilon_{2}-\epsilon_{1})\beta}{\beta(\alpha+\beta)}$

Recalling that, in all cases we have considered, we always choose
$\epsilon_{1}=1$ $\epsilon_{2}=1$, $\epsilon_{3}=-1$, the above expression
becomes:

$-2-4\frac{\alpha^{2}\rho(1-\rho)}{\beta(\beta+\alpha)}$

which is surely negative, $\forall 0<\rho<1$, provided
$\frac{\alpha}{\beta}\in(-1,+\infty)$. Once fulfilled, the condition will hold
for any point on the equilibrium line. Finally, let us compare the behaviour
of the nonlinear model and its linearization, given in subsection (3.3) (see
equations (60) and the discussion after these equations). We make two
numerical examples by taking the following values of the parameters:

$\epsilon_{1}=\epsilon_{2}=1,\;\;\epsilon_{3}=-1,\;\;\alpha=1,\;\;\beta=2,\;\;$
(64)

The equilibrium positions are given by

$N_{1}^{0}=\frac{\rho}{3},\;\;N_{2}^{0}=\frac{1-\rho}{2},\;\;N_{3}^{0}=\frac{1}{2}-\frac{\rho}{6}.$
(65)

We start with an initial condition close to the equilibrium point
$(N_{1}^{0},N_{2}^{0},N_{3}^{0})=(\frac{1}{6},\frac{1}{4},\frac{5}{12})$
corresponding to $\rho=\frac{1}{2}$: we choose
$(N_{1}(0),N_{2}(0),N_{3}(0))=(\frac{1}{6}+0.1,\frac{1}{4},\frac{5}{12})$. The
trajectory determined by the numerical solution of the non-linear system (1)
and the one determined by the analytical solution of the set of equation (60)
are compared in figure (10). Also, the evolution of the population $N_{1}(t)$
and the evolution of $N_{1}^{0}+x_{1}(t)$ is given in figure (11)

Figure 10: Plot of the closed orbit $(N_{1}(t),N_{2}(t),N_{3}(t))$ for the
system (1) corresponding to the initial conditions $N_{1}(0)=\frac{1}{6}+0.1$,
$N_{2}(0)=\frac{1}{4}$,$N_{3}(0)=\frac{5}{12}$ with
$\epsilon_{1}=\epsilon_{2}=1$, $\epsilon_{3}=-1$, $\alpha=1$, $\beta=2$ (in
black) and of
$(\frac{1}{6}+x_{1}(t),\frac{1}{4}+x_{2}(t),\frac{5}{12}+x_{3}(t))$, where
$(x_{1}(t),x_{2}(t),x_{3}(t))$ are solution of the system (60) with initial
conditions $x_{1}(0)=0.1$, $x_{2}(0)=0$,$x_{3}(0)=0$ (in red). Figure 11: Plot
of the evolution of $N_{1}(t)$ for the system (1) corresponding to the initial
conditions $N_{1}(0)=\frac{1}{6}+0.1$,
$N_{2}(0)=\frac{1}{4}$,$N_{3}(0)=\frac{5}{12}$ with
$\epsilon_{1}=\epsilon_{2}=1$, $\epsilon_{3}=-1$, $\alpha=1$, $\beta=2$ (in
black) and of $\frac{1}{6}+x_{1}(t)$, where $x_{1}(t)$ is the solution of the
system (60) with initial conditions $x_{1}(0)=0.1$, $x_{2}(0)=0$,$x_{3}(0)=0$
(in red).

Now we take an initial condition closer to the equilibrium point
$(N_{1}^{0},N_{2}^{0},N_{3}^{0})=(\frac{1}{6},\frac{1}{4},\frac{5}{12})$ with
respect to the previous one. We choose
$(N_{1}(0),N_{2}(0),N_{3}(0))=(\frac{1}{6}+0.01,\frac{1}{4},\frac{5}{12})$.
Again, the trajectory determined by the numerical solution of the non-linear
system (1) and the one determined by the analytical solution of the set of
equation (60) are compared in figure (12). Also, the evolution of the
population $N_{1}(t)$ and the evolution of $N_{1}^{0}+x_{1}(t)$ is given in
figure (13). As it can be seen by the figures the trajectories are really
close. Also, the period of the linear system
$\textrm{T}=\frac{2\sqrt{2}\pi}{\sqrt{|\textrm{Tr}(M^{2})|}}$ is a very good
approximation for the period of the system in these cases.

Figure 12: Plot of the closed orbit $(N_{1}(t),N_{2}(t),N_{3}(t))$ for the
system (1) corresponding to the initial conditions
$N_{1}(0)=\frac{1}{6}+0.01$, $N_{2}(0)=\frac{1}{4}$,$N_{3}(0)=\frac{5}{12}$
with $\epsilon_{1}=\epsilon_{2}=1$, $\epsilon_{3}=-1$, $\alpha=1$, $\beta=2$
(in black) and of
$(\frac{1}{6}+x_{1}(t),\frac{1}{4}+x_{2}(t),\frac{5}{12}+x_{3}(t))$, where
$(x_{1}(t),x_{2}(t),x_{3}(t))$ are solution of the system (60) with initial
conditions $x_{1}(0)=0.01$, $x_{2}(0)=0$,$x_{3}(0)=0$ (in red). Figure 13:
Plot of the evolution of $N_{1}(t)$ for the system (1) corresponding to the
initial conditions $N_{1}(0)=\frac{1}{6}+0.01$,
$N_{2}(0)=\frac{1}{4}$,$N_{3}(0)=\frac{5}{12}$ with
$\epsilon_{1}=\epsilon_{2}=1$, $\epsilon_{3}=-1$, $\alpha=1$, $\beta=2$ (in
black) and of $\frac{1}{6}+x_{1}(t)$, where $x_{1}(t)$ is the solution of the
system (60) with initial conditions $x_{1}(0)=0.01$, $x_{2}(0)=0$,$x_{3}(0)=0$
(in red).

The agreement between the “exact” nonlinear picture and the linearized one
looks extremely, and somehow astonishingly, good.

###  5.1 Extension to $n$ species

We recall that the integrable case of the Volterra equations with $n$ species
reads

$\frac{dN_{r}}{dt}=\epsilon_{r}N_{r}+\sum_{s\neq r=1}^{n}A_{rs}N_{r}N_{s}$
(66)

where the matrix of the interactions must be taken as

$A_{rs}=\epsilon_{r}\epsilon_{s}(B_{r}-B_{s})~{}~{}~{}r,s=1,\cdots n$ (67)

Let us look for an integral of motion of the type:

$e^{A\sum_{i}N_{i}(t)}\prod_{i=1}^{n}N_{i}(t)^{w_{i}},$ (68)

where $A$ and $w_{i}$ are $N+1$ suitable constants. By requiring that the
derivative of (68) with respect to time is zero and by using equations (66),
we get, after some manipulations:

$\sum_{k=1}^{n}\epsilon_{k}w_{k}=0,\quad\quad
A=-\sum_{k=1}^{n}B_{k}\epsilon_{k}w_{k},$ (69)

so that the quantities

$e^{-\sum_{k=1}^{n}B_{k}\epsilon_{k}w_{k}\left(\sum_{i}N_{i}(t)\right)}\prod_{i=1}^{n}N_{i}(t)^{w_{i}}=I_{1,..,n},$
(70)

are a parametric family of conserved quantities, depending on $n-1$ parameters
since the constraint $\sum_{k=1}^{N}\epsilon_{k}w_{k}=0$ must be satisfied.
Notice that equation (70) can be written as the following product of functions
of a single variable:

$\prod_{i=1}^{n}\frac{N_{i}^{w_{i}}}{e^{AN_{i}}}=I_{1,..,n}$ (71)

where $A$ is defined in (69). The function $f(x)=\frac{x^{w}}{e^{ax}}$ plays a
crucial role in a quite simple proof that the orbits of the original two-
populations Volterra model is periodic (see e.g. [16]). Indeed, for $x$
positive, it has a maximum in $x=\frac{w}{a}$ and decreases to zero by going
to $x=0$ and $x=\infty$, similar to a gaussian bell shaped curve. Let us look
at equation (71) supposing to vary just the value of the parameter
$I_{1,..,n}$ (the other parameters being fixed):

$f_{1}(N_{1})\cdot f_{2}(N_{2})\cdot\ldots\cdot f_{n}(N_{n})=I_{1,...,n},\quad
f_{i}(N_{i})\doteq\frac{N_{i}^{w_{i}}}{e^{AN_{i}}}$ (72)

Each of the functions $f_{i}(N_{i})$ has a maximum in $N_{i}=\frac{w_{i}}{A}$
equal to $M_{i}=\left(\frac{w_{i}}{Ae}\right)^{w_{i}}$: it follows that if

$I_{1,...,n}>\prod_{i=1}^{n}M_{i}$ (73)

then the equation (72) has no real solution. Also, if

$I_{1,...,n}=\prod_{i=1}^{n}M_{i}$ (74)

then equation (72) has just one real solution, given by
$N_{i}=\frac{w_{i}}{A}$, $i=1,...,n$. Let us suppose now that

$I_{1,...,n}<\prod_{i=1}^{n}M_{i}.$ (75)

We set

$I_{1,...,n}=\lambda\prod_{i=2}^{n}M_{i},\quad\lambda<M_{1}$ (76)

and look at the equation

$f_{1}(N_{1})\cdot f_{2}(N_{2})\cdot\ldots\cdot
f_{n}(N_{n})=\lambda\prod_{i=2}^{n}M_{i},\quad\lambda<M_{1}$ (77)

The equation $f_{1}(N_{1})=\lambda$ has just two real solutions, let us call
them $N_{1}^{-}$ and $N_{1}^{+}$, since $f_{1}(N_{1})$ increases from zero to
$M_{1}$ and then decreases to zero at infinity. Clearly one has
$N_{1}^{-}<\frac{w_{1}}{A}$ and $N_{1}^{+}>\frac{w_{1}}{A}$. It follows that
if $N_{1}=N_{1}^{-}$ or $N_{1}=N_{1}^{+}$ equation (79) has just one real
solution (the other $N_{i}$ being given by $N_{i}=\frac{w_{i}}{A}$,
$i=2,...,n$). When $N_{1}<N_{1}^{-}$ or $N_{1}>N_{1}^{+}$ equation (79) has no
solution, since one has

$f_{2}(N_{2})\cdot\ldots\cdot
f_{n}(N_{n})=\frac{\lambda}{f_{1}(N_{1})}\prod_{i=2}^{n}M_{i}>\prod_{i=2}^{n}M_{i}.$
(78)

Finally it remains the case $N_{1}\in(N_{1}^{-},N_{1}^{+})$. In this case,
having fixed the value of $N_{1}$, we can repeat the above reasoning on
$N_{2}$, since now we have

$f_{2}(N_{2})\cdot\ldots\cdot
f_{n}(N_{n})=\frac{\lambda}{f_{1}(N_{1})}M_{2}\prod_{i=3}^{n}M_{i}=\hat{\lambda}\prod_{i=3}^{n}M_{i},\quad\hat{\lambda}<M_{2}$
(79)

and there will be two real solutions to the equation
$f_{2}(N_{2})=\hat{\lambda}$, $N_{2}^{-}<\frac{w_{2}}{A}$ and
$N_{2}^{+}>\frac{w_{2}}{A}$. By repeating the same until the last $N$, we see
that equation (79), for $N_{i}>0,i=1,...,n$, represents a compact closed
surface, isomorphic to the n-sphere. In the following we plot the surface in
the case $n=3$. We take

$\frac{N_{1}^{2}}{e^{2N_{1}}}\frac{N_{2}^{3}}{e^{2N_{2}}}\frac{N_{3}^{4}}{e^{2N_{3}}}=0.001.$
(80)

The maximum of the function
$\frac{N_{1}^{2}}{e^{2N_{1}}}\frac{N_{2}^{3}}{e^{2N_{2}}}\frac{N_{3}^{4}}{e^{2N_{3}}}$
equals to $54e^{-9}\sim 0.0067$ so equation (80) possesses solutions. The
corresponding surfaces is plotted in figure (14)

Figure 14: Plot of
$\frac{N_{1}^{2}}{e^{2N_{1}}}\frac{N_{2}^{3}}{e^{2N_{2}}}\frac{N_{3}^{4}}{e^{2N_{3}}}=0.001.$

##  6 Concluding remarks

A general comment is needed here, on the comparison between our results and
those displayed on pages (245-250) of [2], where the three basic properties of
fluctuations in conservative associations are stated and explained. There,
first of all Volterra distinguishes between even and odd number of species. In
the odd case in fact the coefficient matrix $A$, being skew symmetric, has to
be singular: Volterra concludes that equilibrium states will be impossible,
and the number of individuals in some species will grow indefinitely or go to
zero, and on the long run. only an even number of species will survive.
Curiously enough, Volterra did not take into account the fact that complete
integrability could change, even drastically, the above scenario. In fact, in
section 5 we have shown that even for an odd number of species, in spite of
the singular nature of the matrix $A$, there might bounded trajectories and
even $periodic$ orbits.

It is somehow surprising that the results we have derived do not really rely
upon the hamiltonian formulation of integrable systems; as a matter of fact,
even the results we have got for the two species case could have been derived
without resorting to the hamiltonian formulation. Moreover, ubiquitous
features in the mathematical picture of integrable system, such as Lax
representation, do not appear in our description.

A further question, more important in view of the possible applications, is
related to the construction of a richer model, still enjoying the complete
integrability property, but involving a larger number of conjugated pairs of
variables, all of them being significant in the economic-ecological approach.
Work is progress in all the outlined directions.

Acknowledgments

F.Z. wishes to acknowledge the support of Università degli Studi di Brescia;
INFN, Gr. IV - Mathematical Methods in NonLinear Physics and ISNMP -
International Society of Nonlinear Mathematical Physics. O.R. and F.Z. wish to
acknowledge the support of GNFM-INdAM.

## References

  * [1]
  * [2] V.Volterra, Principes de Biologie Mathematiques, Acta biotheoretica (Leiden,3, 6-39)(1937)
  * [3] G.Israel, Vito Volterra, Book on Mathematical Biology (1931), Landmark Writings in Western Mathematics, 1640-1940, edited by I. Grattan-Guinness, Amsterdam, Elsevier, 2005, pp. 936-944.
  * [4] M.Scalia,A.Angelini,F.Farioli,GF Mattioli,O.Ragnisco,M.Saviano, An Ecology and Economy coupling model. A global stationary state model for a sustainable economy in the Hamiltonian formalism in Ecological Economics ,172,106497(2020)
  * [5] V.Volterra, Lecons sur la Theorie Mathematique de la lutte pour la vie, Gauthier-Villars Paris (1927).
  * [6] M.Scalia, O.Ragnisco The Volterra integrable case, ArXiv 1903 03595 [nlin.si], 2019.
  * [7] N. Stern, Stern Review on the Economics of the Climate Change, report presented by Sir Nicholas Stern, Head of the Government Economic Service and Advisor to the Government on the economics of climate change and development, to the Prime Minister and the Chancellor of the Exchequer,http://webarchive.nationalarchives.gov.uk/20100407172811 , http:// www.hm-treasury.gov.uk/stern review report.htm (2006).
  * [8] N. Georgescu-Roegen, The entropy law and the economic process, Harvard University Press, Cambridge MA (1971)
  * [9] N. Georgescu-Roegen, The Steady State and Ecological Salvation. A thermodynamic Analysis, Bioscience, 27, 266-70 (1977)
  * [10] H.Daly, Steady-State Economics, Island Press, Washington, D.C.(1977); Steady-State Economics. Second edition with new essays, Island Press, Washington, D.C. (1991)
  * [11] O.Ragnisco and M.Scalia, A Stationary State Model for a Sustainable Economy, Preprint on Research Gate (2020).
  * [12] A.Angelini, F.Farioli, G.Mattioli and M.Scalia, Le due crisi: crisi del capitalismo e crisi ambientale. Una soluzione sostenibile? (Parte II), Culture della sostenibilità , 17, 83-97 (2016)
  * [13] M.Scalia, A.Angelini, F.Farioli, G.Mattioli and M.L.Saviano, The chariots of Pharaon at the Red sea: the crises of capitalism and environment. A modest proposal towards sustainability, I Quaderni di Culture della Sostenibilità, Num. Spec. I Sem. anno IX, 4-64 (2016).
  * [14] R.M.Goodwin, A growth Cycle, in C.H. Feinstein, Editor, Socialism, Capitalism and Economic Growth, Cambridge University Press (1967).
  * [15] G.Israel, A.M. Gasca, Letters between Umberto D’Ancona and Vito Volterra in The Biology of Numbers , Birkhauser Basel, 130-200 (2002).
  * [16] M.Braun, Differential Equations and Their Applications, Predator-Prey problems,(pp. 441-449), Springer Berlin (1992).
  * [17] R. L. Fernandes and W. M. Oliva, Hamiltonian dynamics of the Lotka-Volterra equations, In International Conference on Differential Equations (Lisboa, 1995) (pp. 327-334). World Sci. Publ., River Edge, NJ. (1998).
  * [18] V.I. Arnold, Mathematical Methods of Classical Mechanics (2nd Edition), Springer Berlin (1980).
  * [19] R.Abraham and J.Marsden Foundation of Mechanics, Taylor and Francis Ltd (2020).
  * [20] H. Goldstein, C.P. Poole and J.L. Safko, Classical Mechanics (3rd Edition), Addison-Wesley, New York ( 2001).


KCL-PH-TH-2021-68

# Footprints of population III stars in the gravitational-wave background

Katarina Martinovic  Theoretical Particle Physics and Cosmology Group, Physics
Department,
King’s College London, University of London, Strand, London WC2R 2LS, UK
Carole Périgois  LAPP, CNRS, 9 Chemin de Bellevue, 74941 Annecy-le-Vieux,
France  Tania Regimbau  LAPP, CNRS, 9 Chemin de Bellevue, 74941 Annecy-le-
Vieux, France  Mairi Sakellariadou  Theoretical Particle Physics and Cosmology
Group, Physics Department,
King’s College London, University of London, Strand, London WC2R 2LS, UK

(September 2021)

###### Abstract

We investigate detection prospects of the gravitational-wave background (GWB)
that originates from the merging of compact objects formed by the collapse of
population III stars. Younger population I/II stars lead to a GWB in the
LIGO/Virgo frequency band at the inspiral phase, while population III stars
would likely show up at the later merger and ringdown phases. We show that,
using a network of third-generation detectors, we may be able to separate a
population I/II signal from a population III one, provided we can subtract
individual coalescence events. A detection of a population III GWB could
reveal important information, such as the average redshifted total mass.

Introduction—We have witnessed a rapid expansion of gravitational-wave (GW)
astrophysics in the last decade due to the success of the Advanced LIGO and
Virgo GW detectors Aasi et al. (2015); Acernese et al. (2015) at uncovering
signals from numerous compact binary coalescences (CBCs) KAGRA:2013rdx . The
most recent Advanced LIGO/Virgo observing run, O3, presented us with dozens of
new merger events and has significantly expanded the stellar graveyard
LIGOScientific:2020ibl . Despite the increase in the number of detections, we
are yet to observe with confidence an event that would suggest that the
progenitor compact objects are remnants of the oldest stars in the Universe
LIGOScientific:2020kqk  \- the theoretically-postulated population III (pop
III) stars (GW190521 could be a potential candidate Pacucci:2017jbp ). Pop III
stars are thought to have formed at high redshifts and as such have low
metallicity compared to the more recently formed, population I/II (pop I/II)
stars Bromm:1999du ; Bromm:2003vv ; Yoshida:2003rw . These old stars have
hitherto evaded sky surveys deSouza:2011ea ; Bowler:2016qen ;
2019BAAS…51c.449W , and their detection remains an objective for upcoming
experiments, such as the James Webb Space Telescope Gardner:2006ky .

Pop III stars may solve some of the puzzles in black hole formation, as well
as help understanding the early epochs of the Universe such as reionisation
and galaxy evolution 1984ApJ…277..445C ; Yoshida:2003ab ; Trenti:2009cj .
Numerical simulations show that these primordial stars could have led to the
formation of super-massive black holes at high redshifts 2018MNRAS.475.4104C ;
Hirano:2015wxa ; Heger:2002by . Mergers of such heavy remnants would appear in
the millihertz frequency range explored by future space-based detectors such
as LISA. The scope of this study, however, is detection prospects of
terrestrial detector networks, and we therefore focus on models that predict a
pop III signal in the LIGO/Virgo frequency range. The contribution to the
gravitational-wave background (GWB)111Often referred to as the stochastic
gravitational-wave background Rosado:2011kv . from a superposition of
unresolved pop III-seeded CBCs has been explored in several studies
Belczynski:2001uc ; Suwa:2007du ; Inayoshi:2016hco ; Belczynski:2016ieo ;
Ng:2020qpk . They show significant deviation of a pop III star signals from a
pop I/II stars signal due to different mass and redshift distributions
Ng:2020qpk . The GWB is comprised of many sources, of astrophysical or
cosmological origin, but we expect the CBC signal to be the foreground to all
sources KAGRA:2021kbb . In this study, for the first time, we consider the
possibility of separating pop I/II and pop III GWB contributions. Numerous
models suggest the total CBC background is dominated by pop I/II. However, pop
III can be uncovered using subtraction techniques and studying the residual
backgrounds Regimbau:2014uia ; Sachdev:2020bkk . As the sensitivity of
detectors increases and GW interferometers see more individual CBC events, a
pop III residual background emerges as the dominant signal over pop I/II
residual background. We first study how to detect the GWB from pop III stars,
and in the case of a successful detection, we explore subsequent implications,
namely information about masses and redshifts of the population.

The paper is organised as follows: in Sec. II we present the pop III models
and their resulting GWB in different detector networks, highlighting
characterisation of the residual backgrounds. We then introduce Bayesian
analysis used in Sec. III, and discuss search filters we consider for pop III
stars. Sec. IV is an implications study in case of a detection of a pop III
signal. We use StarTrack (ST) simulation data Belczynski:2005mr  and apply our
detection methods, ultimately showing consistency of our implications analysis
with the underlying population. We select the ST data since it is the most
recent extensive catalogue of merging binaries from pop I, II and III stars
that lead to a GWB in the LIGO/Virgo frequency range Perigois:2020ymr .

Population III GWB— The GWB is defined as the superposition of GWs from all
unresolved sources. It is characterised by the dimensionless parameter
$\Omega_{\rm GW}$ Allen:1997ad , expressed as the ratio of GW energy density
per logarithmic frequency bin $\text{d}\rho_{\rm GW}/\text{d}\,\text{ln}(f)$,
normalised by the critical energy density of the Universe, $\rho_{\rm
c}=(3H_{0}^{2}c^{2})/(8\pi G)$:

$\Omega_{\rm GW}(f)=\frac{1}{\rho_{\rm c}}\frac{\text{d}\rho_{\rm
GW}(f)}{\text{d}\text{ln}(f)},$ (1)

with $H_{0}=67.9$ km s-1 Mpc-1 Planck:2018lbu .

Here we concentrate on the CBC contribution to the GWB, namely from pop I/II
stars and the theoretical pop III stars. One can express the quantity
$\Omega_{\rm GW}$ in terms of CBC source parameters $\theta$ (masses and
spins), as Perigois:2020ymr  :

$\Omega_{\rm GW}(f)=\frac{f}{\rho_{c}H_{0}}\int\text{d}\theta
p(\theta)\int_{0}^{z_{\rm
up}(\theta)}\text{d}z\frac{R(z;\theta)\frac{\text{d}E_{\rm
GW}(f_{s};\theta)}{\text{d}f_{s}}}{(1+z)E_{z}(z)},$ (2)

where $p(\theta)$ is the probability distribution of the source parameters,
${\text{d}E_{\rm GW}}/{\text{d}f_{\rm s}}$ in the energy density emitted by a
single source at a redshift $z$ with parameters $\theta$, $f_{\rm s}$ is the
emitted frequency in the source frame $f_{\rm s}=f(1+z)$ and $z_{\rm
up}(\theta)$ is the maximal redshift at which a compact binary with parameters
$\theta$ can form. The factor $(1+z)$ in the denominator converts the merger
rate $R(z,\theta)$ from the source to the detector frame, and $E_{z}(z)$
accounts for the considered cosmology, i.e. the expansion history of the
Universe, $E_{z}(z)=\sqrt{\Omega_{\rm m}(1+z)^{3}+\Omega_{\Lambda}},$ with
$\Omega_{\rm m}=0.31,\Omega_{\Lambda}=0.69$ Planck:2018lbu .

For the total population, when all sources are included in the background, the
merger rate $R(z;\theta)$ at redshift $z$ for sources with parameters $\theta$
is given in the source frame, per unit of comoving volume and time. It is
derived from the star formation rate, corrected by the time delay between the
birth of the progenitors and the merger of the compact objects Barack:2018yly
. To calculate the residual background when individually detected sources are
removed, one has to multiply the total rate by a factor
$1-\epsilon(z,\theta)$, where the efficiency $\epsilon(z,\theta)$ is the
probability for a source at redshift $z$ with parameters $\theta$ to be
detected, integrated over inclination, polarisation and position in the sky
(see Regimbau:2014nxa ).

In the case of binary neutron stars and neutron star-black hole mergers, we
only consider the inspiral phase and assume that the emission of GWs stops at
the last stable orbit. For binary black holes, we consider the three different
regimes of the coalescence (inspiral, merger and ringdown phase) given by the
corresponding phenomological waveforms Ajith:2009bn  calculated for circular
orbits. The energy density is Regimbau:2008nj :

$\displaystyle\frac{\text{d}E_{\rm GW}}{\text{d}f_{\rm s}}(f_{\rm s})$
$\displaystyle=\frac{5(G\pi)^{2/3}\mathcal{M}_{\rm
c}^{5/3}F_{\iota}}{12}f_{\rm s}^{-1/3}$ (3)
$\displaystyle\times\begin{cases}(1+\sum^{3}_{i=2}\alpha_{i}\nu^{i})^{2}&\text{if
}f_{s}<f_{\rm merg}\\\ f_{\rm s}w_{\rm
m}(1+\sum^{2}_{i=1}\epsilon_{i}\nu^{i})^{2}&\text{if }f_{\rm merg}\leq f_{\rm
s}<f_{\rm ring}\\\ f_{\rm s}^{1/3}w_{\rm r}\mathcal{L}^{2}(f_{\rm s},f_{\rm
ring},\sigma)&\text{if }f_{\rm ring}\leq f_{\rm s}<f_{\rm cut}\end{cases}$

where $f_{\rm merg}$, $f_{\rm ring}$ and $f_{\rm cut}$ are the frequencies at
the start of merger, start of ringdown and end of emission in the source
frame, respectively. The chirp mass of the binary is a combination of the
individual masses of the compact objects, $\mathcal{M}_{\rm
c}=\frac{(m_{1}m_{2})^{3/5}}{(m_{1}+m_{2})^{1/5}}$. $\mathcal{L}(f_{\rm
s},f_{\rm ring},\sigma)$ is the Lorentzian function centered at $f_{\rm
ring}$, with width $\sigma$, and $w_{m}$, $w_{r}$ are the normalisation
constants ensuring the continuity between the three phases. The factors
$\epsilon_{i}$ and $\alpha_{i}$ and the frequencies $f_{\rm merg}$, $f_{\rm
ring}$ and $f_{\rm cut}$ follow from analytical waveforms detailed in
Ajith:2009bn  and depend on the symmetric mass ratio
$\eta=(m_{1}m_{2})/(m_{1}+m_{2})^{2}$ of the progenitors’ masses, $m_{1}$ and
$m_{2}$, and the effective spin of the system
$\chi=[(m_{1}\vec{s}_{1}+m_{2}\vec{s}_{2})/(m_{1}+m_{2})]\vec{L}/L$. Once the
spectrum of $\Omega_{\rm GW}$ is calculated, we can estimate the corresponding
signal-to-noise ratio (SNR) for a given network of $N$ detectors
Romano:2016dpx :

$\text{SNR}=\frac{3H_{0}^{2}}{10\pi^{2}}\sqrt{2T}\left[\int_{0}^{\infty}\text{d}f\>\sum_{I=1}^{N}\sum_{J>I}\frac{\gamma_{IJ}^{2}(f)\Omega_{\rm
GW}^{2}(f)}{f^{6}P_{J}(f)P_{J}(f)}\right]^{1/2}\,,$ (4)

with $T$ the observational time, $P_{I}$ and $P_{J}$ the one-sided power
spectral noise densities of detectors $I$ and $J$, and $\gamma_{IJ}$ the
normalised isotropic overlap reduction function characterising the distance
and the relative orientation between $I$ and $J$ for sources isotropically
distributed in the sky.

In this work, we consider the StarTrack model FS1 for pop III
Belczynski:2016ieo . Assumptions about characteristics of the initial binary
pop III stars are discussed in 10.1093/mnras/stv2629 . The FS1 model assumes
that pop III stars were formed in large gas clouds with a star formation rate
that peaks at redshift $z\sim 12$, while the star formation rate for pop I/II
stars peaks at $z\sim 2$ (see Fig. 4 in Belczynski:2016ieo ). Even though pop
III stars are less abundant than pop I/II, the ST model FS1 considers a rather
optimistic ratio of pop III to pop I/II stars. The corresponding background
and its detectability have been calculated Perigois:2020ymr  using a catalogue
of sources rather than the analytical expression in Eq. 2. The residual
background catalogue is obtained by subtracting all sources individually
detected by the interferometer network. For each source $k$ we calculate the
individual SNR $\rho^{k}$ assuming optimal-matched filtering and uncorrelated
gaussian noise in the detectors as follows:

$\left(\rho^{k}\right)^{2}=\sum_{I=1}^{N}4\int_{f_{i,\min}}^{f_{i,\max}}\frac{\left|\tilde{H}^{k}\right|^{2}}{P_{I}(f)}\
df,$ (5)

where

$\tilde{H}^{k}=F_{+,I}(f,\Theta^{k},\psi^{k})\tilde{h}^{k}_{+}(f)+F_{\times,I}(f,\Theta^{k},\psi^{k})\tilde{h}^{k}_{\times}(f),$
(6)

with $F_{+,I}$ and $F_{\times,I}$ the antenna factors of detector $I$ for
polarisations $+$ and $\times$ that depend on source inclination $\Theta^{k}$
and position in the sky $\psi^{k}$, while $\tilde{h}_{+}^{k}$ and
$\tilde{h}_{\times}^{k}$ are the Fourier transforms of the gravitational
waveforms of the source $k$. A residual catalogue is computed by removing all
sources with $\rho^{k}>12$.

If pop III exists, its signal will be superposed with a pop I/II signal. In
the case of a dominant pop I/II signal, the pop III signal will remain hidden
underneath it, and one can only place upper limits on the amplitude of the pop
III contribution to $\Omega_{\rm GW}(f)$. If, however, a pop III signal is the
dominant one, then we could detect deviations from the 2/3 CBC power law, and
even getting insight on the mass and redshift distribution of pop III stars.
We explore the last scenario by considering two terrestrial networks of third-
generation (3G) detectors: (i) Einstein Telescope (ET) at the Virgo site, and
(ii) ET at the Virgo site with two Cosmic Explorers (CE) at the LIGO Hanford
and Livingston sites.

Estimates of CBC contributions to $\Omega_{\rm GW}$ from ST simulations
suggest that pop III signal is lost in the pop I/II foreground. For 2G
detector networks – even by including LIGO-Hanford, LIGO-Livingston, Virgo,
LIGO-India, and KAGRA – pop III is practically invisible and its contribution
to the global SNR is negligible, as it is shown in Perigois:2020ymr . However,
3G detectors such as ET and CE, may reveal a pop III background. The future
detectors will have unprecedented sensitivity and they will be able to
discover a great number of individual CBCs, thereby reducing the GWB
originating from unresolved CBCs. For ET+2CE, we uncover pop III after the
subtraction of individually resolved merger events. This follows because
subtraction methods are less efficient to detect the high redshift and low
frequency pop III CBCs. Being more difficult to resolve, binaries from pop III
persist, resulting in a large contribution to the residual CBC background in
3G detectors.

We compare in Fig. 1 the total and residual background for the two 3G
networks: ET (top) and ET+2CE (bottom). It confirms that the pop III
contribution in ET has a very small impact on the combined residual background
from pop I/II and pop III, while in ET+2CE the pop III residual background
clearly dominates for frequencies below $\sim$ 20Hz. In addition, Fig. 1 shows
a change in the shape of the background: The peak frequency of pop III changes
slightly while the slope characterising the end of emission decreases
dramatically when we remove individually detected sources.

Figure 1: Total and residual GWB of ET (top) and ET+2CE (bottom) detector
networks. The pop I/II and pop III contributions are shown in green and red,
respectively, with the combined residual signal shown in black.

To demonstrate the impact of subtraction of resolved CBCs on the population,
we show in Fig. 2 the probability density of the total redshifted mass,
$M_{\rm{tot}}^{z}=(1+z)(m_{1}+m_{2})$, and the merger rate $R(z)$ as a
function of redshift, between the whole catalogue and the residual one for
ET+2CE.

Figure 2: Comparison between the total (blue) and ET+2CE residual (orange)
catalogue for redshifted total mass distributions (top) and merger rates
(bottom).

Clearly, the sources remaining in the residual catalogue are the ones with the
highest redshift, affecting the total corrected mass distribution which is in
turn responsible for the changes in the GWB spectrum. We will estimate the
ET+2CE residual pop III parameters by filtering the corresponding background
and performing a Bayesian analysis.

Detection method— The stochastic pipeline takes strain data $\tilde{s}_{I,J}$
from detectors, $I,J$, and constructs cross-correlation statistics using
optimal filters Romano:2016dpx :

$\hat{C}_{IJ}(f)=\frac{2}{T}\frac{{\rm
Re}[\tilde{s}_{I}^{*}(f)\tilde{s}_{J}(f)]}{\gamma_{IJ}(f)S_{0}(f)},$ (7)

with $T$ the duration of the run, and $\gamma_{IJ}(f)$ the normalised overlap
reduction function as defined in Eq. 4. The estimator is normalised with
$S_{0}(f)=(3H_{0}^{2})/(2\pi^{2}f^{3})$ leading to
$\langle\hat{C}_{IJ}(f)\rangle=\Omega_{\rm GW}(f)$. We assume correlated noise
not to be a limiting factor to our detector sensitivity and consider all noise
to be gaussian. The variance is

$\sigma_{IJ}^{2}(f)\approx\frac{1}{2T\Delta
f}\frac{P_{I}(f)P_{J}(f)}{\gamma_{IJ}^{2}(f)S_{0}^{2}(f)}.$ (8)

Let us construct a gaussian log-likelihood,

$\displaystyle p(\hat{C}_{IJ}(f|\theta))$ $\displaystyle\propto$
$\displaystyle\exp\left[-\frac{1}{2}\sum_{IJ}\left(\frac{\hat{C}_{IJ}(f)-\Omega_{\rm
GW}(f|\theta)}{\sigma_{IJ}(f)}\right)^{2}\right],$

where $\Omega_{\rm GW}(f|\theta)$ represents the GWB model with parameters
$\theta$. This allows us to estimate the model parameters by finding the best-
fit to the cross-correlation data and minimising the likelihood function. Note
that we have made the simplifying assumption that the log-likelihood of a
detector network is the sum of log-likelihoods of the individual baselines. To
compare models and find which ones fit data better, we perform model selection
with Bayes factors. Bayes factor,
$\mathcal{B}^{\mathcal{M}_{1}}_{\mathcal{M}_{2}}$, is defined as the ratio of
evidences of model $\mathcal{M}_{1}$ to model $\mathcal{M}_{2}$, and if large
and positive, demonstrates preference for $\mathcal{M}_{1}$ over
$\mathcal{M}_{2}$.

Typically, we model a CBC signal as $\Omega_{\rm GW}(f)=\Omega_{\rm
ref}\leavevmode\nobreak\ (f/f_{\rm ref})^{2/3}$, with $f_{\rm ref}$ = 25 Hz.
This is because the CBCs detected so far have low masses that would lead to an
inspiral signal in the low-frequency range. This can be seen in Fig. 1 where
the total GWB from pop I/II and III in 3G detectors are presented. The pop III
spectrum shows clear deviation from a 2/3 power law, because these further
away stars will lead to more redshifted frequencies and therefore be detected
in their merger and ringdown phases. We test search filters different from a
2/3 power law to investigate if the deviation from pop I/II signal can be
identified in a parameter estimation study. Motivated by the shape of the
residual pop III signal in Fig. 1, we consider the following filters:

  * •

power law with varying spectral index (PL)

$\Omega_{\rm GW}^{\rm PL}(f)=\Omega_{\rm ref}\leavevmode\nobreak\ (f/f_{\rm
ref})^{\alpha}$ (10)

  * •

broken power law (BPL)

$\Omega_{\rm GW}^{\rm BPL}(f)=\begin{cases}\Omega_{\rm peak}(f/f_{\rm
peak})^{\alpha_{1}}&\text{for $f\leq f_{\rm peak}$},\\\ \Omega_{\rm
peak}(f/f_{\rm peak})^{\alpha_{2}}&\text{for $f>f_{\rm peak}$}.\end{cases}$
(11)

  * •

smooth BPL

$\Omega_{\rm GW}^{\rm SBPL}(f)=\Omega_{\rm peak}\leavevmode\nobreak\ (f/f_{\rm
peak})^{\alpha_{1}}\leavevmode\nobreak\ [1+(f/f_{\rm
peak})^{\Delta}]^{(\alpha_{2}-\alpha_{1})/\Delta}.$ (12)

  * •

triple BPL

$\Omega_{\rm GW}^{\rm TBPL}(f)=\begin{cases}\Omega_{\rm peak}(f/f_{\rm
peak}^{(1)})^{\alpha_{1}}&\text{for $f\leq f_{\rm peak}^{(1)}$},\\\
\Omega_{\rm peak}(f/f_{\rm peak}^{(1)})^{\alpha_{2}}&\text{for $f_{\rm
peak}^{(1)}<f\leq f_{\rm peak}^{(2)}$},\\\ k\,\Omega_{\rm peak}(f/f_{\rm
peak}^{(2)})^{\alpha_{3}}&\text{for $f>f_{\rm peak}^{(2)}$},\end{cases}$ (13)

where $k=(f_{\rm peak}^{(2)}/f_{\rm peak}^{(1)})^{\alpha_{2}}$ ensures
continuity of the piecewise function.

The priors for each model’s parameters can be found in the Appendix. If any of
the filters above are preferred over a 2/3 filter, this could be an indication
of the presence of a pop III signal.

Implications— In the case of a detection, we examine whether we can constrain
the mass/redshift distribution from the optimal search parameters. Following
the GWB expression (Eq. (2)), we see that the parameters impacting the
background shape are the redshift-dependent merger rate and the black holes’
mass distribution. To understand how these population characteristics relate
to model parameters, such as peak frequency and slope, we generate multiple
spectra. We make simplifying assumptions about our progenitors by assuming
spinless, equal-mass binaries Gerosa:2020bjb . We fix the merger rate and vary
the intrinsic mass input, observing how the shape of the GWB spectrum changes.
The results we find, however, change with a different choice of merger rate,
as described in the Appendix. This is because there is a degeneracy between
the effects that merger rate and mass distribution have on the GWB
Inayoshi:2016hco . We thus study the dependence of $\Omega_{\rm GW}$ on the
redshifted total mass of the population, $M^{z}_{\rm tot}=(1+z)(m_{1}+m_{2})$,
which is related to the merger rate, and find a relationship between the mass
and the peak frequency of the spectrum.

We generate GW spectra with a merger rate from ST, varying the redshifted
total mass, and we find an agreement (within 10%) between redshifted ringdown
frequency and the peak of the spectrum, see Table 1. We obtain the same
agreement if we use the merger rate from Inayoshi:2016hco , suggesting once
more that an estimate of the peak frequency can be used to constrain the
average redshifted total mass of the population. This relationship, therefore,
holds independently of the model used for the evolution of the pop III
binaries.

$M_{\rm{tot}}^{z}$ | $f_{\rm peak}$ | $f^{z}_{\rm ring}$ | % difference
---|---|---|---
100 | 166.2 | 165.8 | 0.20
200 | 83.7 | 82.9 | 1.0
300 | 56.5 | 55.3 | 2.1
400 | 43.1 | 41.4 | 3.9
500 | 35.2 | 33.2 | 5.7
Table 1: Variation of the peak of GWB spectra with a change in redshifted
total mass. We find agreement between the peak frequency and the redshifted
ringdown frrequency.

Results— We simulate one year of observation time with the ET+2CE network,
taking the CBC background from the ST catalogue. We find the best-fit models
to the residual GWB that remains after subtracting the individual sources. A
model selection study shows preference for other filters over a 2/3 PL, see
Table 2. The models with a broken power law shown in the last 3 rows are
clearly favoured over a single power law model. However, we do not observe a
great increase in Bayes factor for the smooth and triple BPL over just a BPL.
Therefore, we conclude that a BPL filter is sufficient for a pop III GWB
search.

model, $\mathcal{M}$ | ln $\mathcal{B}^{\mathcal{M}}_{2/3}$
---|---
PL | 29 000
BPL | 46 000
smooth BPL | 47 000
triple BPL | 46 000
Table 2: Log Bayes factor of pop III filters compared to the 2/3 power law
filter.

Already with a varying-index PL search we deduce that the 2/3 filter is not
appropriate, since the estimated power law index is $\alpha=-0.6$, as observed
in the corner plot in Fig. 3. The more intricate filters, however, fit the
$\Omega_{\rm GW}$ spectrum well and capture the presence of the peak
successfully. In order to understand the redshifted mass distribution of pop
III creating the GWB, we investigate the peak frequency of the signal. We
obtain a good estimate of the peak frequency using a BPL search filter as
shown in Fig. 4, $f_{\rm peak}=15.4$ Hz. The redshifted ringdown frequency
that matches the peak frequency, $f^{z}_{\rm ring}=15.4$ Hz, corresponds to
$M^{z}_{\rm tot}=1076M_{\odot}$. The ST redshifted mass distribution shown in
the top panel of Fig. 2, gives an average redshifted total mass of the
residual population, $<M^{z}_{\rm tot}>=1121M_{\odot}$. Therefore, our
estimate of average $M^{z}_{\rm tot}$ agrees well with the true value.
Finally, the estimated $M^{z}_{\rm tot}$ can be depicted as a curve in the
redshift-intrinsic total mass plane since $M^{z}_{\rm tot}=(1+z)M_{\rm tot}$,
see Fig. 5. Note that we have included a 10% uncertainty error in matching of
the ringdown and the peak frequency for consistency with our findings in Table
1.

Figure 3: Varying-$\alpha$ PL fit to residual GWB spectrum of pop I+II+III
from the ST simulation. We see that the $\alpha$ estimate is not 2/3 which
would be expected for pop I/II. Figure 4: BPL fit to residual GWB spectrum of
pop I+II+III from the ST simulation. The peak frequency is estimated to
$f_{\rm peak}=15.4$ Hz. Figure 5: $M^{z}_{\rm tot}=1076M_{\odot}$ constraint
shown in the $M_{\rm tot}$-z plane, including 10% error bars accounting for
the uncertainty of the estimate of $M^{z}_{\rm tot}$ from $f_{\rm peak}$.

Conclusions— GWs emitted from CBCs formed by pop III stars could represent a
promising detection channel of the first stars formed in the Universe in a
very low-metallicity environment. Although 2G detectors are unable to detect
the GWB from pop III stars, we have demonstrated that future GW
interferometers could lead to a detection. Employing the ET+2CE 3G detector
network, we could successfully subtract enough individual merger events to
unravel the elusive pop III GWB. Subtraction methods are more effective for
nearby sources, so the residual pop I/II signal can become sub-dominant to the
pop III residual. The 2/3 power law approximation of GWB energy density for
the CBC background breaks in this case due to the higher redshifted masses of
pop III. Instead of the 2/3 power law, the GWB spectrum peaks in the low-
frequency LIGO/Virgo range.

A model selection study showed that we could observe the peak caused by the
unresolved merger and ringdown GWs from pop III stars. With a good estimate of
the peak frequency, we could even deduce the redshifted total mass
distribution of the residual population. Taking data from the recent ST binary
sources catalogue, we have demonstrated the effectiveness of our Bayesian
analysis combined with implications of retrieved parameters. One should note
that we have derived the relationship between peak frequency and redshifted
total mass assuming equal-mass and non-spinning binaries. Future work should
involve relaxing these assumptions. Additionally, one could further
investigate the connection between the negative slope of the $\Omega_{\rm GW}$
spectrum and pop III properties such as merger rate and mass distribution.

Acknowledgments—The authors acknowledge access to computational resources
provided by the LIGO Laboratory supported by National Science Foundation
Grants PHY-0757058 and PHY-0823459. This paper has been given LIGO DCC number
2100337. K.M. is supported by King’s College London through a Postgraduate
International Scholarship. M.S. is supported in part by the Science and
Technology Facility Council (STFC), United Kingdom, under the research grant
ST/P000258/1.

Numerous software packages were used in this paper. These include matplotlib
Hunter (2007), numpy van der Walt et al. (2011), scipy Virtanen et al. (2020),
bilby Ashton et al. (2019), dynesty Speagle (2020), PyMultiNest Buchner, J. et
al. (2014).

## Appendix

Priors—All of the models have the same log-uniform prior distribution for the
GWB amplitude ranging between $10^{-13}$ and $10^{-5}$. As for the remaining
parameters, we use

  * •

PL: $\alpha=\mathcal{N}(0,3.5)$,

  * •

BPL: $\alpha_{1}=U(2/3,5/3),\alpha_{2}=U(0,-8),f_{\rm peak}=U(10,100)$,

  * •

smooth BPL: $\alpha_{1}=U(2/3,5/3)$, $\alpha_{2}=U(0,-8)$, $f_{\rm
peak}=U(10,100),\Delta=U(0,10)$,

  * •

triple BPL: $\alpha_{1}=U(2/3,5/3)$,
$\alpha_{2}=U(0,-8),\alpha_{3}=\delta(2/3),f_{\rm peak}^{1}=U(10,100),f_{\rm
peak}^{2}=U(10,100).$

For models with a break frequency we use a uniform prior for the first power
law index between 2/3 and 5/3, since this represents the inspiral/merger
regime of the CBC. Triple BPL has the third spectral index fixed to
$\alpha_{3}=2/3$ since we expect the inspiral phase of pop I/II signal to
dominate at higher frequencies.

Intrinsic mass distribution— We study the relation between peak frequency of
pop III GWB spectrum and the mass distribution of the sources. We fix the
merger rate as a function of redshift to be the one of ST. For total mass,
$M_{\rm{tot}}=m_{1}+m_{2}$ varying between 10 and 90 $M_{\odot}$, we generate
$\Omega_{\rm{GW}}$ spectra and record the frequency at which the spectra are
maximum. We then find a best-fit curve for the data,

$f_{\rm peak}=f_{0}\left(\frac{90M_{\odot}}{M_{\rm tot}}\right)\rm Hz,$

with $f_{0}=12.8$ Hz see Fig. 6. However, changing the merger rate to the one
from Inayoshi:2016hco , we find a different best fit curve, with $f_{0}=53.7$
Hz, implying that the intrinsic mass may be difficult to extract from the
estimate of the spectrum peak. We find more promising results if we study the
redshifted total mass and its relation to peak frequency, as described in the
main text.

Figure 6: StarTrack merger rate evolved, equal-mass binaries. We find a
relationship between peak frequency and total intrinsic mass of the merger.
This is a model-dependent statement.

## References

  * Aasi et al. (2015) J. Aasi, B. P. Abbott,et. al., Class. Quant. Grav. 32, 074001 (2015), URL https://doi.org/10.1088/0264-9381/32/7/074001.
  * Acernese et al. (2015) F. Acernese et al. (Virgo), Class. Quant. Grav. 32, 024001 (2015), eprint 1408.3978.
  * (3) B. P. Abbott et al. [KAGRA, LIGO Scientific, Virgo and VIRGO], Living Rev. Rel. 21, no.1, 3 (2018) eprint 1304.0670.
  * (4) R. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. X 11, 021053 (2021), eprint 2010.14527.
  * (5) R. Abbott et al. [LIGO Scientific and Virgo], Astrophys. J. Lett. 913, no.1, L7 (2021) eprint 2010.14533.
  * (6) F. Pacucci, A. Loeb and S. Salvadori, Mon. Not. Roy. Astron. Soc. 471 no.1, L72-L76 (2017), eprint 1706.09892
  * (7) V. Bromm, P. S. Coppi and R. B. Larson, Astrophys. J. Lett. 527, L5-8 (1999) eprint 9910224.
  * (8) V. Bromm and R. B. Larson, Ann. Rev. Astron. Astrophys. 42, 79-118 (2004) eprint 0311019.
  * (9) N. Yoshida, T. Abel, L. Hernquist and N. Sugiyama, Astrophys. J. 592, 645-663 (2003) eprint 0301645.
  * (10) R. S. de Souza, N. Yoshida and K. Ioka, Astron. Astrophys. 533, A32 (2011) eprint 1105.2395.
  * (11) R. A. A. Bowler, R. J. McLure, J. S. Dunlop, D. J. McLeod, E. R. Stanway, J. J. Eldridge and M. J. Jarvis, Mon. Not. Roy. Astron. Soc. 469, no.1, 448-458 (2017) eprint 1609.00727.
  * (12) Windhorst Rogier et al, Bulletin of the American Astronomical Society 51 3, 449 (2019), eprint 1903.06527
  * (13) J. P. Gardner, J. C. Mather, M. Clampin, R. Doyon, M. A. Greenhouse, H. B. Hammel, J. B. Hutchings, P. Jakobsen, S. J. Lilly and K. S. Long, et al. Space Sci. Rev. 123, 485 (2006), eprint 0606175.
  * (14) Carr, B. J. and Bond, J. R. and Arnett, W. D., Astrophys. J. 277 445-469 (1984).
  * (15) N. Yoshida, V. Bromm and L. Hernquist, Astrophys. J. 605, 579-590 (2004), eprint astro-ph/0310443.
  * (16) M. Trenti and M. Stiavelli, Astrophys. J. 694 879-892 (1984), eprint 0901.0711.
  * (17) S. Chon and T. Hosokawa, and N. Yoshida, Mon. Not. Roy. Astron. Soc. 475 no.3 4104-4121 (2018), eprint 1711.05262.
  * (18) S. Hirano, T. Hosokawa, N. Yoshida, K. Omukai and H. W. Yorke, Mon. Not. Roy. Astron. Soc. 448, no.1, 568-587 (2015), eprint 1501.01630.
  * (19) A. Heger, C. L. Fryer, S. E. Woosley, N. Langer and D. H. Hartmann, Astrophys. J. 591, 288-300 (2003), eprint 0212469.
  * (20) P. A. Rosado, Phys. Rev. D 84, 084004 (2011) eprint 1106.5795.
  * (21) K. Belczynski, V. Kalogera and T. Bulik, Astrophys. J. 572, 407-431 (2001) eprint 0111452.
  * (22) Y. Suwa, T. Takiwaki, K. Kotake and K. Sato, Astrophys. J. Lett. 665 L43-L46 (2007), eprint 0706.3495.
  * (23) K. Inayoshi, K. Kashiyama, E. Visbal and Z. Haiman, Mon. Not. Roy. Astron. Soc. 461 no.3, 2722-2727 (2016), eprint 1603.06921.
  * (24) K. Belczynski, T. Ryu, R. Perna, E. Berti, T. L. Tanaka and T. Bulik, Mon. Not. Roy. Astron. Soc. 471, no.4, 4702-4721 (2017), eprint 1612.01524.
  * (25) K. K. Y. Ng, S. Vitale, W. M. Farr and C. L. Rodriguez, Astrophys. J. Lett. 913, no.1, L5 (2021) eprint 2012.09876.
  * (26) R. Abbott et al. [KAGRA, Virgo and LIGO Scientific], Phys. Rev. D 104, no.2, 022004 (2021), eprint 2101.12130.
  * (27) T. Regimbau, D. Meacher and M. Coughlin, Phys. Rev. D 89, no.8, 084046 (2014) eprint 1404.1134.
  * (28) S. Sachdev, T. Regimbau and B. S. Sathyaprakash, Phys. Rev. D 102, no.2, 024051 (2020) eprint 2002.05365.
  * (29) K. Belczynski, V. Kalogera, F. A. Rasio, R. E. Taam, A. Zezas, T. Bulik, T. J. Maccarone and N. Ivanova, Astrophys. J. Suppl. 174, 223 (2008) eprint 0511811.
  * (30) C. Périgois, C. Belczynski, T. Bulik and T. Regimbau, Phys. Rev. D 103 no.4 043002 (2021) eprint 2008.04890.
  * (31) B. Allen and J. D. Romano, Phys. Rev. D 59, 102001 (1999), eprint 9710117.
  * (32) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A8 (2020), eprint 1807.06210.
  * (33) L. Barack, V. Cardoso, S. Nissanke, T. P. Sotiriou, A. Askar, C. Belczynski, G. Bertone, E. Bon, D. Blas and R. Brito, et al. Class. Quant. Grav. 36, no.14, 143001 (2019) doi:10.1088/1361-6382/ab0587 [arXiv:1806.05195 [gr-qc]].
  * (34) T. Regimbau, K. Siellez, D. Meacher, B. Gendre and M. Boër, Astrophys. J. 799, no.1, 69 (2015) eprint 1410.2739. Ajith:2009bn
  * (35) P. Ajith, M. Hannam, S. Husa, Y. Chen, B. Bruegmann, N. Dorband, D. Muller, F. Ohme, D. Pollney and C. Reisswig, et al. Phys. Rev. Lett. 106, 241101 (2011), eprint 0909.2867.
  * (36) T. Regimbau and V. Mandic, Class. Quant. Grav. 25, 184018 (2008) eprint 0806.2794.
  * (37) T. Ryu, T.L. Tanaka, R. Perna Monthly Notices of the Royal Astronomical Society 456 1 223-238 (2015) eprint 1509.05427
  * (38) J. D. Romano and N. J. Cornish, Living Rev. Rel. 20, no.1, 2 (2017) eprint 1608.06889.
  * (39) D. Gerosa, S. Vitale and E. Berti, Phys. Rev. Lett. 125, no.10, 101103 (2020) doi:10.1103/PhysRevLett.125.101103 [arXiv:2005.04243 [astro-ph.HE]].
  * Hunter (2007) J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * van der Walt et al. (2011) S. van der Walt, S. C. Colbert, and G. Varoquaux, Computing in Science Engineering 13, 22 (2011).
  * Virtanen et al. (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al., Nature Methods 17, 261 (2020).
  * Ashton et al. (2019) G. Ashton et al., Astrophys. J. Suppl. 241, 27 (2019), eprint 1811.02042.
  * Speagle (2020) J. S. Speagle, Mon. Not. Roy. Astron. Soc. 493, 3132 (2020), eprint 1904.02180.
  * Buchner, J. et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., Hsu, L., Rangel, C., Brightman, M., Merloni, A., Salvato, M., Donley, J., and Kocevski, D., A&A 564, A125 (2014), URL https://doi.org/10.1051/0004-6361/201322971.