#  Interactions Between Different Birds of Prey
as a Random Point Process

Gernot Akemann1, Nayden Chakarov2, Oliver Krüger2,
Adam Mielke3, Meinolf Ottensmann2, and Patricia Päßler1 1Faculty of Physics
and 2Faculty of Biology, Bielefeld University, P.O. Box 100131, D-33501
Bielefeld, Germany;
3DTU Compute, Richard Petersens Plads, bygning 324, 2800 Kgs. Lyngby, Denmark
<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

The two-dimensional Coulomb gas is a one-parameter family of random point
processes, depending on the inverse temperature $\beta$. It is proposed as a
simple statistical measure to quantify the intra- and interspecies repulsion
among three different highly territorial birds of prey, by comparing to the
spacing between their nests in the plane. Using data from the area of the
Teutoburger Wald over 20 years, we fit the nearest and next-to-nearest
neighbour spacing distributions between the respective nests of Goshawk, Eagle
Owl and the previously examined Common Buzzard to $\beta$ of the Coulomb gas.
Within each species, the repulsion measured in this way deviates significantly
from the Poisson process of independent points in the plane. In contrast, the
repulsion amongst each of two species is found to be considerably lower and
closer to Poisson. Methodologically we investigate the influence of the
terrain, of a shorter interaction range given by the two-dimensional Yukawa
interaction, and the statistical independence of the time moving average we
use for the yearly ensembles of occupied nests. A simple, approximate
analytical expression for the nearest neighbour spacing distribution derived
from non-Hermitian random matrix theory proves to be very useful, being valid
for the two-dimensional Coulomb gas close to Poisson at $\beta=0$.

##  1\. Introduction

The study of the statistics of random point processes in one (1D) and two
dimensions (2D) is a very active area of research, with many applications in
physics and other sciences. Examples in 2D, on which we will focus here,
include quantum optics [1], quantum chaos [2, 3, 4], condensed matter physics
[5], statistics [6] and ecology [7, 8, 9]. The goal of such an analysis is to
clarify, first, whether the points coming from experimental data or numerical
simulations are independent or not, and if not to quantify their correlations.
The former case is called Poisson point process, where the points are
distributed independently on a line in 1D, or in the plane in 2D. Closed form
expressions exist for arbitrary dimension $D$, cf. [4, Appendix A]. Many point
processes found in nature show correlations, and in particular repulsion
between points in a characteristic, universal way, such that simple models
from statistical mechanics apply. Before describing the point processes we use
in more detail, let us explain what kind of understanding could be expected
from such an approach to biological systems, which typically have a high
degree of complexity.

A central question in biology is to understand the distribution of animals or
plants in space and time. Furthermore, one would like to disentangle the
effect of direct interaction within one species, between different species,
and the effect of the environment. Such an understanding is of great value for
conservation etc. Competition is one of the most ubiquitous features of
ecology [10]. Animals compete for limited resources with individuals of their
own (intraspecific competition) or another species (interspecific
competition), see [11]. One of the easiest and hence most often used estimates
of competition is the nearest neighbour distance [7]. It has been shown to be
a useful measure, with the underlying assumption that a too near neighbour of
the same or another species ultimately leads to decreased reproduction or
survival [12, 13].

Two examples where 2D point processes have been applied in biology are the
spacial distribution of trees, see [8] for a review, and the distribution of
nests of birds of prey in space and time in our previous work [9]. Using data
from [14, 15, 13] on the annual locations of occupied nests of the Common
Buzzard in an area of 300 km2 in the Teutoburger Wald, Germany, cf. Section
2.1, we have been able to draw the following conclusions. The distribution of
their nests differs significantly from 2D Poisson statistics. We have been
able to quantify the repulsion between the locations of nests of these highly
aggressive birds of prey, by fitting the spacing distribution between nearest
(NN) and next-to-nearest neighbours (NNN) in radial distance to those of a
static 2D Coulomb gas of point charges in a confining (Gaussian) potential.
Such a point process is also called Gibbs ensemble, cf. [6]. The single fit
parameter used is the coupling strength $\beta$, representing the inverse
temperature in the 2D Coulomb gas. Details about this approach will be given
in Section 2.2. Moreover, with this simple parametrisation we have been able
to identify a change in time of the interaction strength measured by $\beta$,
over the observed period of 20 years, as a function of the increasing
population density. Let us emphasise that the fitted value of $\beta$ does not
have a direct biological meaning, but merely serves to quantify the repulsion,
in particular the deviation from 2D Poisson statistics at $\beta=0$.

A prime example from physics for such a transition from Poisson to correlated
random variables is the Bohigas–Giannoni–Schmit or quantum chaos conjecture
[18, 19, 20]. It relates to a further example for point processes that is
frequently used in physics [22], the eigenvalues of random matrices [23]. For
self-adjoint matrices these are real (1D), whereas non-Hermitian matrices have
complex eigenvalues in 2D. In both cases these eigenvalues are strongly
correlated random variables, originating from independently distributed matrix
elements in the classical Gaussian ensembles. The quantum chaos conjecture
states that eigenvalue statistics of the Hamiltonian of generic integrable
quantum systems will follow Poisson statistics, according to Berry–Tabor [21],
while that of fully chaotic Hamiltonians obeys random matrix statistics
(distinguished by their global symmetry under time reversal). Here, the NN
spacing distribution between eigenvalues has played an important role to
quantify the path from integrability to chaos in 1D, and ample numerical and
analytical evidence has been collected [22, 24]. This conjecture has been
extended to 2D for dissipative quantum systems [25], and evidence has been
presented, e.g. from the spectrum of the corresponding Liouville operator [2,
3, 4].

The repulsion of random matrix eigenvalues is very well understood. It follows
the logarithmic 2D Coulomb interaction, both for real eigenvalues in 1D and
complex ones in 2D, with only particular values for $\beta=1,2,4$ occurring in
1D, and $\beta=2$ in 2D 111This holds for the classical Ginibre ensembles, see
a further discussion in [2, 26].. The reason is that random matrix eigenvalues
represent determinantal or Pfaffian point processes and are thus integrable,
in the sense that all eigenvalue correlation functions are explicitly known
[23, 27]. Their universality in the limit of a large number of points $N$
(matrix dimension) has been proved in 1D, see [28] and references therein, and
in 2D, cf. [29] and [3] for the NN spacing distribution. This means that they
do not depend on the choice of a Gaussian distribution of matrix elements. The
quantum chaos conjecture has raised the obvious question of how to interpolate
between the Poisson point process, that can be viewed as a special case of
diagonal random matrices at $\beta=0$, and the specific $\beta$-values
occurring for random matrices. In 1D, heuristic, approximate descriptions
exist, e.g. the Brody distribution, see [22] for other proposals. In 2D, it
has been proposed [3] to directly use the 2D Coulomb gas at intermediate
values of $\beta$, and to determine the NN and NNN spacing distributions
numerically, that can then be compared with data. These findings are reviewed
in Section 2.2, in particular a corresponding approximate surmise for the 2D
NN spacing distribution valid for small values of $\beta$ [26], cf. Appendix
B. These NN and NNN spacings were then used for fits as a function of $\beta$,
for the spacing of eigenvalues of the Liouville operator of dissipative,
boundary driven systems in 2D [3], and in the annual spacing distributions of
nests of the Common Buzzard [9]. Here, every year represents a realisation of
the ensemble, and we have taken averages over all years, or windows of time
moving averages over several consecutive years, in order to be sensitive to
the time dependence.

The goal of this article is to go beyond the analysis of a single species [9].
For the three species of competing birds of prey living in the observed area,
the Common Buzzard, Goshawk and Eagle Owl, annual data for the locations of
their nests over the same area and period of time have been collected [15, 13,
16, 17]. Our goal is to try to quantify the interaction between each two
species, and to compare it to the interaction found within each species
individually. Therefore, in a first step we have repeated the analysis from
[9] for Goshawks and Eagle Owls individually, see Section 3.1. This includes
the dependence on the respective change in population, ranging between 1-30
pairs. Because the spacing distribution for Eagle Owls differs considerably,
whether they nest in the Teutoburger Wald (F) or in the adjacent plains (P),
we have further split them into 2 groups. Because of comparably low statistics
we have been unable to address the time dependence for the interaction between
different species so far, to be discussed in Section 3.2.

A further purpose of the present work is to analyse methodological issues in
Section 4, that came up in discussions when presenting the results of [9].
Because all 3 species discussed only nest in forested patches, a first
question addressed in Section 4.1 is whether the repulsion we observe within a
species via a fitted $\beta\neq 0$ is due to the (fractal) area of the forest,
or represents a true repulsion. We have generated a Poisson point process on
the given forest area in the monitored region, and made a fit to a real $D>0$,
for a $D$-dimensional Poisson NN spacing distribution. As we will see, this
makes the repulsion between nests even more pronounced, and is clearly not an
effect of the local area. In [9] the NN and NNN spacings distribution of nests
were found to be described by different (typically decreasing) values of
$\beta$, clearly indicating a shorter interaction range than 2D Coulomb.
Therefore we have investigated numerically the spacing distribution of the
Yukawa interaction in 2D in Section 4.2. It depends on 2 parameters, with a
shorter interaction than logarithmic at larger distance. It reduces to 2D
Coulomb in a limiting case and at small distance. Furthermore, the question of
"independence" of nests occupied in consecutive years is analysed in Section
4.3. Here, we quantify the reusage of nests within one or all species. A large
abundance of old nests exists in the area, and the precise locations of nests
in the monitoring process allows to quantify the respective percentage of
reusage.

For a comparison to the NN or NNN spacing distribution from the 2D Coulomb
gas, the mean density of the data has to be normalised to unity, a procedure
called unfolding. For 1D this is unique and easier than in 2D, cf. [22] vs.
[30, 3]. As an alternative quantity that does not require unfolding, spacing
ratios have been introduced in 1D [31] and 2D [4]. However, they are only
known analytically for Poisson and for random matrix statistics, the complex
Ginibre ensembles in 2D [4, 32]. In contrast to 1D [33], where the spacing
ratio was derived as a function of $\beta$, the interpolation in 2D is not so
clear, as we will see from our data analysis presented in Appendix A. All our
findings, open questions and possible future directions are discussed in
Section 5.

##  2\. Description of Setup: Data and Point process

###  2.1. Biological description and data collection - birds of prey in the
Teutoburger Wald

In this subsection we describe the collection of data we use and the
geographical setup. The locations of occupied nests of three kinds of birds
were collected in late winter and early spring during the years 2000-2020 in
an area in and around the Teutoburger Wald close to Bielefeld. The three
species are the Common buzzard (Buteo buteo L.), Goshawk (Accipiter gentilis)
and Eagle Owl (Bubo bubo). The area of 300 km2 (8∘25’E and 52∘6’N) is located
in Eastern Westphalia. It consists of two 125 km2 grid squares and 50 km2 edge
areas, see Fig. 1 for illustration. The main habitat of these birds of prey is
the Teutoburger Wald and a cultivated landscape to the north and south of it.
This is a low mountain region of height up to 315 m above sea level, and we
will treat the area as approximately flat, that is two-dimensional. All three
species nest in forest patches. Their size varies from rows of trees to large
patches, of more than 10 km2. In total approximately 17% of the study site is
forested. The area has been intensively monitored for birds of prey, and the
resulting spatial data have been published and used before extensively [14,
15, 13].

Figure 1.  Left: Distribution of all observed occupied bird nests in the year
2020: Eagle Owls (large pink points), Goshawks (yellow points) and Common
Buzzard (dark blue points). There exist occupied nests outside the area with
marked points, which have not been monitored. The green shaded area marks the
approximate extend of the Teutoburger Wald based on the roads on either side,
cf. right plot. Eagle Owl nests within this green area are distinguished by a
black ring around the large pink points.
Right: A classification of the observed area in three types of terrain: City
(black), forest (green) and cultivated land in the plains (white). The lack of
population of birds within the smaller cities of Halle or Werther is well
visible on the left plot. The forest found in the North of the Teutoburger
Wald in that area consists of many small patches of irregular shape.

The forest patches were visited in March and April each year to look for
incubating birds, and if a nest was occupied in both months, the pair was
classified as breeding. The locations of these nests were marked in large
scale maps or using GPS-devices to monitor the spacial distribution of
occupied nests. In Fig. 1 left, a snapshot for the year 2020 is shown, with
occupied nests of the three species marked with coloured points, see caption.
The classification of the monitored area into forest, city or plains is shown
in Fig. 1 right.

The example year 2020 chosen has amongst the highest population for all three
kinds of birds. Already by eyesight several features emerge. Common Buzzards
are much more abundant than Goshawks and Eagle Owls. Whereas Common Buzzards
and Goshawks spread relatively evenly over the monitored area, the density of
Eagle Owls inside the Teutoburger Wald is much higher than outside. Because of
this difference in distribution of Eagle Owl nests, we decided to split the
group of Eagle Owls in two: those nesting in the forest area (F) in green in
the left plot (nests marked by a black rings around the pink dot), and those
in the plains (P) (nests marked by a pink dot only). Notice that the plains
were populated by Eagle Owls only after 2010, see Fig. 5 bottom right for
their yearly populations. The marking of the forest area in green in Fig. 1
left is of course approximate. In particular, one might argue whether the two
outliers of Eagle Owls South of the Teutoburger Wald should be counted as (F)
or (P), in view of the two large patches of forest in Fig. 1 right. Although
living in large forest patches, they are separated from the main forest by
cities and are thus counted as (P) here. Such a choice does not change the
quantitative picture, with the low statistics for Eagle Owls being difficult
anyhow.

A second feature is visible by eye: There are holes in the population
densities around cities (marked in black in Fig. 1 right). Although effects of
the edge of a density are known in point processes, see the discussion at the
end of Subsection 2.2 below, we will disregards them here and treat all points
as bulk points.

The third feature that is immediate from Fig. 1 left is that the mean spacing
between nests within one species is very different for Goshawks, Common
Buzzards, Eagle Owls (F) and Eagle Owls (P). It is approximately given by the
inverse density (for a constant density). While this mean in meters may encode
important biological information about the range of interaction within one (or
different) species, we are interested in comparing with (universal) features
from spatial statistics in simple 2D point processes. Such a comparison can
only be made quantitative, if the mean density is normalised, and the
fluctuations around this mean, as given for instance by the local spacing
distribution, is measured and compared with the predictions from (equally
normalised) point processes.

The procedure of normalising the mean density is called unfolding and very
well studied in applications of random matrix theory, see [22]. For data in
one spatial dimension it is unique and mostly straight forward, see [22, Sect.
3.2.1]. In two spatial dimensions it is more delicate, see [30, 3] for
different procedures. Here, we will apply the method proposed in [3]. Taking
for example the Common Buzzards only, each blue point in Fig. 1 left is
replaced by a Gaussian distribution of a certain width, the sum of which
defines the approximate mean density $\rho_{\rm ave}$ of points in 2D. The
measured spacing around point $z=(x_{0},y_{0})$ is then normalised by
rescaling with $\sqrt{\rho_{\rm ave}(x_{0},y_{0})}$, to achieve a spacing with
respect to a mean density of unity. This process is done for each year for
each of the 4 sets of birds, recalling that we have split the Eagle Owls into
2 sets, yielding an ensembles for each set. This can then be compared to
random point process to quantify the degree of repulsion, to be described in
the next subsection.

###  2.2. Random point processes - from Poisson to Coulomb gas in 2D

In this subsection we introduce the random point processes and their
corresponding nearest neighbour (NN) and next-to-nearest neighbour (NNN)
spacing distribution that we will use for a comparison to data for the
distribution of occupied nests of birds of prey.

We begin with the Poisson (Poi) random point process in two dimensions, cf.
Section 4.1 for general dimension $D$. It consists of $N$ independent,
uncorrelated points in the plane. The following normalised spacing
distributions in radial distance between NN and NNN are known in the limit of
a large number of points $N$:

(2.1) $\displaystyle p_{\mathrm{Poi}}^{\mathrm{(NN)}}(s)$ $\displaystyle=$
$\displaystyle\;\frac{\pi}{2}s\ {\rm e}^{-\pi s^{2}/4}\quad\sim s\ ,$ (2.2)
$\displaystyle p_{\mathrm{Poi}}^{\mathrm{(NNN)}}(s)$ $\displaystyle=$
$\displaystyle\;\frac{\pi^{2}}{8}s^{3}\ {\rm e}^{-\pi s^{2}/4}\,\sim s^{3}\ ,$

see e.g. [24, 4] for details of the derivation. The NN distribution has its
first moment normalised to unity to set the scale (for the NNN spacing the
first moment follows from this scale and is given by $3/2$). These spacing
distributions can be obtained by distributing $N$ uncorrelated points on a
disc of radius $R$. In the limit $N\to\infty$, after rescaling the mean
spacing between points as $\bar{s}=R\sqrt{\pi/(4N)}$, eqs. (2.1) and (2.2) are
obtained in units of $\bar{s}$. Despite resulting from uncorrelated points,
the 2D area measure in polar coordinates, $dxdy=sds\,d\theta$ for
$z=x+iy=se^{i\theta}$, leads to a linear (cubic) repulsion in the NN (NNN)
spacing distribution. We will not consider higher order spacings distributions
in the following. As we will see below, the situation of a uniform density on
a disc can also be obtained for a 2D Coulomb gas. Let us emphasise that the
above spacing distributions quantify local correlations amongst points at
distance $\sim 1/\sqrt{N}$, compared to global distances of the order of unity
on the (unit) disc.

Let us move to the 2D static Coulomb gas (Cou). For a finite number $N$ of
(charged) points it is given by the equilibrium distribution at inverse
temperature $\beta=(k_{B}T)^{-1}$, subject to the logarithmic long-range
interaction and a confining potential $V(z)$. The latter is chosen to be
Gaussian here for simplicity, $V(z)=|z|^{2}$,

(2.3) $\displaystyle\mathcal{P}_{\mathrm{Cou},\beta}(z_{1},\dots,z_{N})$
$\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}_{N,\beta}}\exp\left[\beta\sum_{{j,k=1;j<k}}^{N}\log|z_{k}-z_{j}|-\sum_{j=1}^{N}|z_{j}|^{2}\right],$

where $\mathcal{Z}_{N,\beta}$ is the normalising partition function. The $N$
points are represented by complex coordinates $z_{j}\in\mathbb{C}$,
$j=1,\ldots,N$, with the usual identification $\mathbb{C}\sim\mathbb{R}^{2}$.
Compared to standard conventions, where $\beta$ multiplies the entire
Hamiltonian including the confining potential, we have absorbed $\beta$ in
front of the potential by rescaling the coordinates $\beta V(z)\to V(z)$. The
rescaling is made in order to allow us to take the limit $\beta\to 0$ to reach
Poisson statistics, while maintaining a confining potential. Furthermore, the
charge is set to unity (or absorbed into $\beta$).

Such static Coulomb gases have been much studied in the recent mathematical
literature, and we refer to [34] for a recent review. After a further
rescaling of the potential $V(z)\to NV(z)$, in the large-$N$ limit the mean
density $\rho(z)$ of points condenses on the so-called droplet $S$. It is
given by the Laplacian of the potential,
$\rho(z)=\frac{2}{\beta}\partial_{z}\partial_{\bar{z}}V(z)$, Frostman’s
equilibrium measure. For general rotationally invariant potentials $V=V(|z|)$
the support $S$ is given by a disc of fixed radius, that can be rescaled to
the unit disc.

For the local correlations among points at distance of order $1/\sqrt{N}$ very
little is known for fixed $\beta>0$, apart from Poisson statistics at
$\beta=0$ (that extends to $\beta\sim 1/N$, cf. [35]) and the integrable case
$\beta=2$, when the point process (2.3) becomes determinantal, cf. [27].
Inspired by the Wigner surmise for the 1D Dyson gas, and its generalisation to
general $\beta$ based on a $2\times 2$ $\beta$-ensemble [36], a surmise (sur)
for the NN spacing distribution was derived from complex normal $2\times 2$
random matrices in [26]222Notice a typo in the normalisation constant in [26]:
$\alpha^{\beta}$ there should be replaced by $\alpha^{1+\beta/2}$ as here.

(2.4) ${p}^{\rm(NN)}_{{\rm
sur},\beta}(s)=\frac{2\alpha^{1+\beta/2}}{\Gamma[1+\beta/2]}\,s^{1+\beta}\exp[-\alpha
s^{2}]\sim s^{1+\beta},$

where $\alpha={\Gamma[(3+\beta)/2]^{2}}/{\Gamma[1+\beta/2]^{2}}$. Its
behaviour at small values $s\to 0$ is as expected heuristically from (2.3),
with one power from the radial measure and a power $\beta$ from the
Vandermonde determinant. Unfortunately, it is well known [25] from the
integrable case at $\beta=2$ in 2D, cf. (2.9) below, that here $N=2$ is not a
good approximation to the large-$N$ limit. However, in the limit $\beta\to 0$
eq. (2.4) exactly reproduces the NN spacing of the Poisson distribution (2.1).
This characteristic is opposite to the Wigner surmise in 1D, which becomes
more accurate for increasing values of $\beta$, rather than for $\beta\to 0$.

Figure 2. Left: The approximate, surmised NN spacing distribution
${p}^{\rm(NN)}_{{\rm sur},\beta_{\rm eff}}(s)$ from (2.4) and (2.5), varying
from $\beta=0$ (red) to $\beta=1$ (violett) in steps of 0.1. The maximum
increases and moves from left to right. Right: At $\beta=2$ the surmise (full
line) is no longer close to the exact answer from the Ginibre ensemble (2.9),
truncated at $N=20$ with rescaled, normalised first moment (dashed line).

In order to improve the approximation of the surmise (2.4) in 2D to larger
values of $\beta$ (and with exact results for $N>2$ being unavailable), an
effective $\beta_{\rm eff}$ was introduced in [26], by fitting a third order
polynomial in $\beta$ to the numerically determined spacing of the 2D Coulomb
gas in the range of $\beta\in[0,3]$:

(2.5) $\beta_{\rm eff}(\beta)=2.108\beta-0.190\beta^{2}+0.030\beta^{3}.$

Consequently in ${p}^{\rm(NN)}_{{\rm sur},\beta_{\rm eff}}(s)$ the
normalisation constant has to change too, $\alpha(\beta)\to\alpha(\beta_{\rm
eff})=\alpha_{\rm eff}$, to ensure a normalised spacing and first moment equal
to unity. This leads a reasonable approximation up to $\beta\approx 0.5$, with
a standard deviation of up to $\sigma=2.6\cdot 10^{-2}$. See Appendix B for
plots in this range, and Fig. 2 right for a comparison at $\beta=2$, where the
surmise clearly fails. In [26] a more detailed comparison is presented,
including higher $\beta$-values and their Kolmogorov-Smirnov distances.
Obviously, the fit using an effective $\beta_{\rm eff}$ spoils the
heuristically expected proportionality at very small argument $\sim
s^{1+\beta_{\rm eff}}$. However, the limit $\beta\to 0$ is still reproduced
exactly. Being based on $2\times 2$ matrices, there is no approximate
prediction for the NNN spacing possible.

In our comparison to data in the next section we will use the cumulative
distribution function (CDF) of (2.4), in order to avoid any dependence on the
choice of histograms. It is given by

(2.6) $\displaystyle{E}^{\rm(NN)}_{{\rm sur},\beta_{\rm eff}}(s)$
$\displaystyle=$ $\displaystyle 1-\frac{\Gamma\left(1+\frac{\beta_{\rm
eff}}{2},\alpha_{\rm eff}s^{2}\right)}{\Gamma\left(1+\frac{\beta_{\rm
eff}}{2}\right)},$

where $\Gamma(n+1,x)=\int_{x}^{\infty}t^{n}e^{-t}dt$ is the upper incomplete
Gamma function.

For completeness we give the analytical result for the NN [25] and NNN [38]
spacing distribution for $\beta=2$, that follows from the Ginibre ensemble
[37]. Here, the 2D Coulomb gas (2.3) has a representation in terms of complex
eigenvalues of complex non-Hermitian random matrices with Gaussian
distribution. In this case, the logarithmic interaction term can be written in
terms of the modulus square of the Vandermonde determinant,

(2.7)
$\Delta_{N}(z_{1},\ldots,z_{N})=\det[z_{i}^{j-1}]_{i,j=1}^{N}=\prod_{k>l}^{N}(z_{k}-z_{l})\
,$

of the $N$ eigenvalues $z_{j}$. Hence the point process is determinantal, and
all complex eigenvalue correlation functions are explicitly known at finite
$N$. The spacing distributions can be derived from the limiting CDF or gap
probability $E_{{\rm Gin}}(s)$, to find an eigenvalue at the origin and the
closest non-zero complex eigenvalue at radial distance $s$,

(2.8) $E_{{\rm
Gin}}^{\rm(NN)}(s)=\prod_{j=1}^{\infty}\frac{\Gamma(1+j,s^{2})}{j!}\ .$

For finite $N$ the product extends only to $N-1$. The limiting spacing
distributions at infinite matrix dimension follow from differentiation,
compare [25] for NN and [38] for NNN:

(2.9) $\displaystyle p^{\mathrm{(NN)}}_{\mathrm{Gin}}(s)$ $\displaystyle=$
$\displaystyle E_{{\rm
Gin}}^{\rm(NN)}(s)\sum_{j=1}^{\infty}\frac{2s^{2j+1}{\rm
e}^{-s^{2}}}{\Gamma(1+j,s^{2})}\sim s^{3},$ (2.10) $\displaystyle
p^{\mathrm{(NNN)}}_{\mathrm{Gin}}(s)$ $\displaystyle=$ $\displaystyle E_{{\rm
Gin}}^{\rm(NN)}(s)\sum_{j,k=1;k\neq
j}^{\infty}\frac{\gamma(1+j,s^{2})}{\Gamma(1+j,s^{2})}\frac{2s^{2k+1}{\rm
e}^{-s^{2}}}{\Gamma(1+k,s^{2})}\sim s^{5},$

where $\gamma(1+k,s^{2})=\int^{s^{2}}_{0}t^{k}{\rm e}^{-t}{\rm d}t$ is the
lower incomplete Gamma function, and we give again the behaviour at $s\to 0$.
The products converge very rapidly and both spacing distributions are
normalised to unity. For simplicity, we have given the expressions above where
the first moment $\bar{s}_{1}$ of the NN spacing is not yet normalised to
unity. It can only be determined numerically, for a product truncated at
sufficiently high value of, say $N\approx 20$, or larger. The NN with first
moment of unity is then obtained by rescaling
$\hat{p}^{\mathrm{(NN)}}_{\mathrm{Gin}}(y)=\bar{s}_{1}p^{\mathrm{(NN)}}_{\mathrm{Gin}}(\bar{s}_{1}y)$,
and the correspondingly rescaled NNN spacing reads
$\hat{p}^{\mathrm{(NNN)}}_{\mathrm{Gin}}(y)=\bar{s}_{1}p^{\mathrm{(NNN)}}_{\mathrm{Gin}}(\bar{s}_{1}y)$.

The spacing distributions (2.9) and (2.10) hold throughout the bulk of the
spectrum inside the supporting unit disc and are known to be highly universal.
That is they hold for a much larger class of confining potentials $V(|z|)$
than Gaussian, in general random normal matrix ensembles [29].

In the Ginibre ensemble at $\beta=2$, the complex eigenvalue correlations at
the edge of the droplet have also been investigated. They are also universal
and agree with the correlations at an inner edge, that is a droplet with a
hole, that can be studied in the induced Ginibre ensemble [39]. In our data
there is no outer edge in Fig. 1, it is simply given by the area of
observation and there are birds nesting outside the area containing dots.
However, we do observe inner edges, as the cities of Bielefeld, Halle or
Werther show up as holes in Fig. 1. We have not been able to address such edge
correlations, mainly due to lack of statistics, but also because it is not so
clear where to draw the boundaries of the inner edges. In our analysis we have
thus treated all points as bulk points.

In the next section we will use both to the approximate NN spacing
${p}^{\rm(NN)}_{{\rm sur},\beta_{\rm eff}}(s)$ from the surmise [26], and the
numerically determined NN and NNN spacing from [3] as functions of $\beta$, in
order to quantify the repulsion within one species, and between different
species.

##  3\. Interaction of species

###  3.1. Interaction within one species: Common Buzzard, Goshawk and Eagle
Owl

Let us describe the procedure initiated in [9] for Common Buzzards, to
quantify the supposedly repulsive interaction between these birds of prey, by
fitting the NN and NNN spacing distributions between occupied nests. Each year
is observed for each species and treated as an ensemble. In principle, for
each year a value of $\beta$ can be determined for each species, as it is done
in Fig. 3 left column for the NN spacing distribution in the year 2020 as an
example. As it can be seen, the quality of the fit decreases rapidly with the
amount of data available, in particular for the Eagle Owls. We show both the
fit using the Coulomb gas (full line) with the spacings determined numerically
in [3] are in discrete steps of size $0.1$ for $\beta$, and the explicit
formula (2.4) for the surmise (dashed line), with $\beta$ varying
continuously. The fits are done using the CDF, compare (2.6), to avoid any
dependence on the choice of binning. To guide the eye we nevertheless show the
spacing distribution compared to histograms of the data. The fitted
$\beta$-values vary over a wide range. For example, in the year 2020 in Fig. 3
left column, we find the strongest repulsion among Goshawks ($\beta=3$),
followed by the Common Buzzard ($\beta=1.1$) down to Eagle Owls (F) with
$\beta=0.7$. Apart from the latter plot with very low statistics, the surmise
and Coulomb gas fit agree well. Notice that the continuous respectively
discrete fit in $\beta$ adds to the discrepancy between the two, see Appendix
B for a comparison at equal values. We find that all fitted values are far
from 2D Poisson statistics at $\beta=0$. Let us emphasise, that the fitted
parameter $\beta$ does not have a biological meaning, but that it is the
relative strength that allows us to compare between different species of
birds, or their interaction in the next Subsection 3.2.

2020 2011-2020 2000-2020

Eagle Owls (F) Goshawks Buzzards

Figure 3. Comparison of the fitted $\beta$-values from the NN spacing
distribution from 2D Coulomb in steps of $0.1$ (full line), respectively from
the surmise (2.4) and (2.5) with continuous $\beta$ (dashed line), for Common
Buzzards (top row), Goshawks (middle row), and Eagle Owls (F, bottom row). We
present examples for fits for a single year 2011 (left column), an average
over 10 consecutive years from 2011-2020 (middle column), and over all years
(right column). Shown are histograms and spacing distribution, whereas the
fits were obtained form the cumulative distributions.

The population of pairs of Common Buzzards and Goshawks per year is shown as
red crosses in Fig. 4 bottom left, respectively right, and for Eagle Owls in
Fig. 5 bottom in the forest (F, left) and plains (P, right). Between Common
Buzzards and Goshawks there is a factor of 10 in abundance, and another factor
of 2 between Goshawks and Eagle Owls, with the latter being completely absent
in the plains up to the year 2010. To remedy the resulting low statistics, we
have introduced a time moving average, averaging over 10 consecutive ensembles
(years) for all three species. The correspondingly averaged population is
shown as a black line in the population plots in Figs. 4 and 5, respectively.
For comparison, because of the better statistics for the Common Buzzards, in
[9] we choose a time moving average of 5 years, that gave a better resolution
of the time dependence of the fitted $\beta$ as a measure of repulsion. In [9]
we also presented the fits for all individual years (compare Fig. 4 left in
[9]). This result is very noisy and makes it difficult to see an overall trend
in time, which is why we do not present such plots here.

The fitted $\beta$-value for such time averaged NN spacing distributions for
each species is shown in Fig. 3 middle column, choosing the average over the
period 2011-2020 as an example, that includes the single year 2020 from the
left column. It is striking to see that in all three cases the ensemble
average leads to a much lower $\beta$-value. This is especially true for the
Goshawks, which now show a repulsion comparable to the Common Buzzards. There
is a clear trend that the fitted $\beta$ value goes down for the Goshawks over
the years, see Fig. 4. Notice that their population is peaked in 2011 and
2012. For the Eagle Owls (F) we obtain a fitted value $\beta=0.1$, being very
close to 2D Poisson in this window of time average. Also here the averaged
fitted $\beta$ goes down over time, see Fig. 5 top left.

Finally, if we abandon any time resolution and make an ensemble average over
all available years, we obtain closer values for $\beta$ for all three
species, as shown in Fig. 3 right column. The strongest repulsion is again
amongst Goshawks ($\beta=0.8$), followed by Common Buzzards ($\beta=0.6$), and
almost equally among Eagle Owls (F) with $\beta=0.5$. The particularly high
$\beta$ for intraspecific repulsion in Goshawk makes a lot of ecological
sense, as the Goshawk is one of the most territorial species in Europe [40].
The Goshawk preys on larger prey and needs exclusive access to hunting areas
[41]. Therefore, intruding individuals are always met with high aggression and
are regularly killed by the territory owners [42, 40, 43]. All values are
significantly above 2D Poisson statistics ($\beta=0$). Because of the strong
growth in population for Common Buzzards, and Eagle Owls (F), see Figs. 4 and
5 respectively, the later years are weighted much higher here in the average
over all years, because they contribute to more spacings per year.

$\beta$-fit Buzzards $\beta$-fit Goshawks

Common Buzzard population Goshawk population

Figure 4. Top row: Fit of $\beta$ for the NN spacing for Common Buzzards
(left) and Goshawks (right) from the Coulomb gas (blue crosses) and surmise
(circles), as well as for the NNN spacings distribution (red crosses). We use
a time moving average of 10 years, where the middle year is indicated on the
x-axis (e.g. 2004+1/2 for the period 2000-2009). Bottom row: Population of
birds per year (red crosses) for Common Buzzards (left) and Goshawks (right),
over the entire period of time 2000-2020. The resulting average populations
using a time moving average of 10 years is shown in comparison (black line),
where again the middle year is chosen as a label, as in the top plots. For a
better comparison these points are connected by the line of the moving
average.

In the following we will discuss our findings for the $\beta$-values fitted to
the time moving average. Two general observations can be made. First, the
values obtained using the surmise agree quite well with the NN spacings form
the 2D Coulomb gas, apart from the years with low statistics for the Eagle
Owls (F). In particular, they follow the same trend as seen from the NN
Coulomb gas fits and thus can be used as a very simple approximate, but
analytical tool. Second, the $\beta$-values obtained from NN and NNN fits do
not agree in general, sometimes lying systematically below (Common Buzzards,
Eagle Owls) or above (Goshawks) the NN value. This discrepancy may be used to
compare the observed interaction range to that of the Coulomb interaction. The
same observation was already made in [9] and motivated us to introduce and
study an interaction with shorter range, the 2D Yukawa interaction, see
Section 4.2. It turned out, however, that also with such a two-parameter
family (with $\beta$ and $\gamma$) the phenomena of obtaining two different
values for $\beta$ from NN and NNN remains, see the discussion there.
Therefore, we kept the simpler 2D Coulomb gas model as a measure for
repulsion.

In Figs. 4 and 5 we do not give error bars for the fitted values of $\beta$.
In our previous paper [9] we made an attempt to estimate the error by fitting
$\beta$ to approximately the same number of points from a Coulomb gas
simulation, that is $N=200$ for the Common Buzzards. This estimate (which is
not related to the data of occupied nests) seemed to overestimate the
fluctuations, e.g. when comparing to a linear fit for the trend in $\beta$.
Because in the ensembles of Goshawks and Eagle Owls the number of nests per
year is smaller by a factor of 10, we refrain from giving such rough estimates
for the errors.

For the Common Buzzards the following picture emerges from Fig. 4 left, that
was already found in [9]. The NN values vary in range between $\beta=0.3$ and
0.8 over the observed period. Although they still fluctuate considerably for
this time moving average, an increasing trend can be identified. In contrast,
the NNN values remain close to $\beta=0$ for half of the observed period, and
then jump approximately to the NN values from the mid-year 2010 onwards. In
comparison, the population growth is approximately linear from around 100
pairs on average to above 200 (10 year average given by the black curve in
Fig. 4 bottom left). Apparently below a certain population threshold the
interaction range is rather limited to NN, and then increases in range to be
Coulomb like up to NNN. Notice that because of unfolding the data, the mean
spacing is the same in all years. Thus the trivial effect of having a smaller
spacing for a larger density is removed here.

For the Goshawks we observe in Fig. 4 right that there is a clear trend for
the fitted $\beta$-value to go down from approximately 1.2 to 0.4 for NN and
from 1.6. to 0.5 for NNN. In contrast to the Common Buzzards, the NNN values
are systematically above the NN values, apart from the last middle year. This
means that the interaction range is longer than 2D Coulomb, and the repulsion
decreases over time. The population development is also different from the
Common Buzzards when comparing the bottom rows in Fig. 4. Notice the
difference in abundance of a factor of 5-10 between the two species. Although
the Goshawk population fluctuates between 12 and 29 pairs, the average only
slightly goes up from 17.5 to 22.5, stabilises and then slightly goes down
again. Thus a different factor than the Goshawk population seems to be at work
here. If it lies in the interaction between the two other species, it will be
answered (negatively) in the next Subsection 3.2.

$\beta$-fit Eagle Owls (F) $\beta$-fit Eagle Owls (P)

Eagle Owl (F) population Eagle Owl (P) population

Figure 5. The same plots as in Fig. 4 for Eagle Owls (F) left column, and for
the Eagle Owls (P) right columns. Top row: corresponding $\beta$-fits for time
moving averages (F) left, and a fit over all years (P) right. Values above
$\beta=2$ are not shown. Bottom row: corresponding populations of Eagle Owls.
Notice the lack of nests in the plains (right) until 2010.

Finally let us discuss the findings for the Eagle Owls. Because of the lowest
statistics, see Fig. 5, about a factor of 2 in abundance below the population
of Goshawks, this is very difficult, especially for the much smaller group
nesting in the plains (P). For those nesting in the forest (F), the $\beta$
value from NN clearly decreases from large values above 1.5 down to 0.2-0.1,
which is almost comparable with 2D Poisson at $\beta=0$. In contrast, the NNN
values remains consistent with $\beta=0$ over the entire period. A comparison
with the population development from 5 pairs on average to 12 shows an
approximate linear increase on average, see Fig. 5 bottom left. The effect is
thus opposite to the Common Buzzards, the Eagle Owls (F) and Goshawks show a
decreasing repulsion, with short range interaction, despite an increase or
stabilisation in population, respectively. It should be noted that for very
few spacings of the order of unity, the distribution is close to a delta peak,
thus leading to very high values for the fitted $\beta$. For a better
resolution the first NN value at $\beta=3$ is suppressed in Fig. 5 top left.

The time dimension of the strength of the repulsion also reflects the species’
biology as well as population trends in the study area. In Common Buzzard, we
see an increase of beta over time, most likely due to the increase in the
population density until carrying capacity might have been approached in the
last five years. The strength of intraspecific competition is expected to
increase with increasing population density. With regard to the Goshawk, a
decreasing beta over time is mirrored by a decrease in population density over
the last five years, most likely due to displacement by Eagle Owls [13]. Why
$\beta$ shows a decrease for Eagle Owls, the population of which has increased
very rapidly over the last 20 years, is difficult to explain. It could be that
the population is still not approaching carrying capacity and hence
progressively smaller repulsions measured through NN spacings are observed.

One remark regarding 2D Poisson is in order here. In Subsection 4.1 we ask the
question, if the scattered patches of forest (see Fig. 1) where all birds
prefer to nest, does not introduce an effective lower dimension than 2. The
effective reduction we find there, by generating a Poisson point process on
the forest patches only, is from $D=2$ to approximately $D=1.66$, see
Subsection 4.1 for more details. We could therefore conclude, that even when
finding $\beta\approx 0$ for the fit, this may not yet indicate a complete
absence of repulsion.

###  3.2. Interaction among two species

In this subsection we quantify the repulsion between all combinations of two
different species of birds, while keeping the two groups of Eagle Owls (F) and
(P) separate, see Figure 1 left. Following the Coulomb gas analogy, it would
be perhaps natural to associate to each species a different charge, according
to their observed (average) repulsion strength. However, such a multi-
component Coulomb gas would depend on the ratio of the different charges,
which may change over the years. If one charge is very abundant, we may
consider the others to be screened, but this is also not always the case. The
multi-parameter fit is therefore not easy to do, and we stick to the simple
one-parameter fits for each pair of species instead, see Figure 6. We thus fit
one $\beta$-value to the spacings between species $A$ and $B$.

vs. Common Buzzards  vs. Goshawks

Eagle Owls (P) Eagle Owls (F)

Goshawks

vs. Common Buzzards

Figure 6. 2D Coulomb gas and surmise fits for the $\beta$-values of the NN
spacing distributions (full respectively dashed line) between all pairs of
species for all years, where for Eagle Owls we distinguish (F) and (P) in the
top left row and bottom left row, respectively. When the best fit to the 2D
Coulomb gas gives $\beta=0$ (bottom middle and right), we fit instead to the
Poisson distribution with $D<2$ from Eq. (4.1) (dashed line), with the
resulting value for $D$ given in the inset.

In order to avoid an over counting of spacings, we always go through the
number of points (nest) of the less abundant species between the two, and find
its NN of the more abundant species in a given year, defining an ensemble.
Because in all years the number of pairs of Common Buzzards is larger than the
pairs of Goshawks which in turn is larger than the number of pairs of Eagle
Owls ((F) or (P)), the ordering is clear. For example, for the Common Buzzard
- Goshawk interaction we go through all Goshawk nests in a given year and find
its NN Common Buzzard nest. Our statistic is thus limited by the population of
Goshawks or Eagle Owls, respectively. For that reason, we choose to make an
average over all ensembles of 21 years in this section. Furthermore, we
restrict ourselves to NN spacings only. The unfolding of the data has been
made using all occupied nests from all three species per year, to get an
approximate mean global density. This certainly represents a simplification,
but does not introduce any extra bias.

The following picture emerges from Fig. 6. All fitted values are closer to
$\beta=0$, with the largest repulsion observed between Common Buzzards and
Eagle Owl (F) at $\beta=0.15$. Thus the 2D Coulomb gas and surmise values are
very close. The NN spacing between Goshawks and Eagle Owls (P) is not well
fitted by Poisson at $D\leq 2$, nor by a 2D Coulomb gas at $\beta>0$.

Because we average over all years, we have to compare with the values in Fig.
3 right column, for the repulsion within each species for all years:
$\beta=0.6$ for Common Buzzards, $\beta=0.8$ for Goshawks and $\beta=0.5$ for
Eagle Owls (F). Clearly, the repulsion measured in $\beta$ is much weaker
between different species than within one species. The strongest interspecies
repulsion measured in this way is between Common Buzzards and Eagle Owls (F)
with $\beta=0.15$, followed by that between Goshawks and Eagle Owls (F), and
Common Buzzards and Eagle Owls (P) at $\beta=0.1$ each. Do these low values of
$\beta$ close to zero indicate, that there is almost no repulsion between
different species, as in a random Poisson point process in 2D? As already
mentioned in the previous subsection, we investigate in Subsection 4.1 if the
random point process on the forest patches reduces the dimension, finding
$D=1.66$ as effective dimension. This seems to indicate that even low
$\beta$-values close to zero represent a certain repulsion. For the
interaction between Common Buzzards and Goshawks we (perhaps coincidentally)
find that the best fit is obtained by the Poisson distribution (4.1) close to
this effective dimension at $D=1.64$, and in that sense there is no repulsion
also according to this measure between the two. The plot with the lowest
statistics between Goshawks and Eagle Owls (P) could be interpreted in the
same way, beside the low quality of the fit.

The finding of a small repulsion between the three species is in line with
empirical findings that intraspecific competition is commonly stronger than
interspecific competition in avian predators [11, 44]. The effect of Eagle
Owls is rather different and has been shown to be clearly negative for both
Common Buzzard and Goshawk [15, 13]. These interactions are, however, very
dynamic at a very small scale [15, 13], and hence the simple models employed
here are perhaps at their limits with regard to detect these interactions.

##  4\. Methodology - Variations of the random point process

###  4.1. Poisson point process in varying dimension $D$

In this section we investigate if the 2D Poisson point process used so far
indeed describes the situation of nests placed as independent random variables
in the plane. It has been observed that all three species of birds of prey
invariably breed in forest patches. A look on Fig. 1 right thus poses the
question, if the forest represents a lower dimensional, fractal domain in the
full two-dimensional area. Notice that we completely ignore the elevation of
the terrain, by treating it as two-dimensional. It varies from about 70 m to
300 m in height above sea level, compared to the dimension of roughly
$12\times 25$ km extent of the monitored area.

We will try to answer this question by generating a Poisson point process
solely on the forested area, determine its NN spacing distribution
numerically, and compare it to the analytic result for the NN spacing
distribution of a Poisson point process for general dimension $D$, by fitting
$D$ as a free parameter. This distribution is well known for integer dimension
$D$, see e.g. [4], and we simply analytically continue to real $D>0$ here:

(4.1) $p_{\rm
Pois,D}^{\rm(NN)}(s)=D\left(\frac{\Gamma(1/D)}{D}\right)^{D}s^{D-1}\exp\left[-\left(\frac{\Gamma(1/D)}{D}\right)^{D}s^{D}\right]\sim
s^{D-1}\ .$

It is normalised with first moment equal to unity. The repulsion $\sim
s^{D-1}$ originates entirely from the (fractional) area measure.

Figure 7.  Left: NN spacing distribution of the Poisson point process (4.1)
for various dimensions from $D=2$ (red) down to $D=1$ (green), in steps of
0.1. The maximum is moving from right to left when going from $D=2$ to $D=1$.
Right: Fit of the dimension $D$ in (4.1) (solid curve) to the random point
process generated on the forested area as described in the main text
(histograms). From the fit we obtain an effective dimension of $D=1.66$. For
comparison, the NN spacing distribution in $D=2$ is also shown (dashed curve).

The Poisson point process on the forest is generated as follows. A certain
number $N$ of points is distributed independently in the monitored area. If
they fall onto a green forest patch in Fig. 1 right, they are accepted, else
they are rejected. Note that we accept all areas with trees, i.e., not just
the main forest depicted as a green band in Figure 1. In the limit $N\gg 1$
with a fixed area, we would recover a collection of two-dimensional patches,
given the number of points per patch is large enough. In order to capture
roughly the same length scale as the distances between occupied nests, we stop
after having $N=200$ accepted points, which is about the number of Common
Buzzard pairs in the entire area. For the fit we use a Kolmogorov-Smirnov fit
to the cumulative distribution of eq. (4.1), see Figure 7 for the result.

We find that the Poisson point process generated as described indeed leads to
slight reduction of dimension from $D=2$ to an effective dimension of
$D=1.66$. Consequently, we may conclude that a fit of data to a Coulomb gas
with $\beta=0$, or equivalently the Poisson point process with $D=2$, still
reflects a small repulsion, compared to the process on the forested area.
Consequently this makes the repulsion found for small $\beta$ more pronounced.

In Subsection 3.2, where the repulsion between different species is quantified
by fitting to the NN spacing of the 2D Coulomb gas with $\beta\geq 0$, we
encounter the situation that the fit to $\beta=0$ is apparently still not
satisfactory, because the maximum of the NN data is further to the left. In
that cases we rather fit the dimension $D$ of the Poisson NN distribution
(4.1), that has a maximum further to the left, see Fig. 6 and the discussion
there.

###  4.2. Varying correlation length: 2D Yukawa interaction

In the previous Section 3 we have seen that fitting $\beta$ independently to
the NN and NNN spacing distribution may lead to different values. This
indicates that the repulsion between the nests cannot always be described by a
2D Coulomb interaction at a single inverse temperature $\beta$. For instance,
when finding $\beta_{\rm NN}>\beta_{\rm NNN}$ we expect an interaction weaker
than Coulomb, or of shorter range. Likewise, for $\beta_{\rm NN}<\beta_{\rm
NNN}$ we expect an interaction stronger than Coulomb, or of longer range on
that scale.

This motivates us to study a Coulomb-like interaction in this subsection,
where the interaction range can be varied: The Yukawa potential $V_{\rm
Yukawa}$ in $D=2$ dimensions. We note in passing that the Yukawa interaction
arises in particle physics from the scattering of distinguishable Fermions
(points) in the non-relativistic limit, see [47] for a standard work on
quantum field theory, which we follow here for the derivation.

The Yukawa potential can be defined in $D$ (integer) dimensions333In most
textbooks only $D=3$ is considered, including [47]. by taking the inverse
Fourier transformation of the propagator with mass $m$ and coupling constant
$g$:

(4.2) $V_{\rm
Yukawa}^{(D)}(\vec{x})=\int\frac{dq^{D}}{(2\pi)^{D}}\frac{-g^{2}}{|\vec{q}|^{2}+m^{2}}e^{i\vec{q}\cdot\vec{x}},\quad\vec{x}\in\mathbb{R}^{D}.$

The integral can be performed using polar coordinates, and in $D=2$ we obtain

(4.3) $\displaystyle V_{\rm Yukawa}^{(2)}(\vec{x})$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{4\pi^{2}}\int_{0}^{\infty}dqq\int_{0}^{2\pi}d\theta\frac{e^{iqr\cos\theta}}{q^{2}+m^{2}}$
$\displaystyle=$
$\displaystyle-\frac{g^{2}}{2\pi}\int_{0}^{\infty}dqq\frac{J_{0}\left(qr\right)}{q^{2}+m^{2}}$
$\displaystyle=$ $\displaystyle-\frac{g^{2}}{2\pi}K_{0}\left(mr\right)\ .$

It only depends on the radial distance $r=|\vec{x}|$. Here, we used the
standard integral representation of the Bessel function of the first kind
$J_{0}(y)$, and $K_{0}(y)$ denotes the modified Bessel function of the second
kind. The last equality follows from [48, 6.532]. It has a logarithmic
singularity at the origin,

(4.4) $K_{0}\left(mr\right)\sim-\log(mr/2)-C,\quad\mbox{for}\ mr\to 0,$

where $C\approx 0.577$ is the Euler–Mascheroni constant, and vanishes
exponentially for large distance as

(4.5) $K_{0}\left(mr\right)\sim\sqrt{\frac{\pi}{2mr}}\
e^{-mr},\quad\mbox{for}\ mr\to\infty.$

If we define a length scale by $\gamma=1/m$, we obtain a one-parameter
deformation of the 2D Coulomb interaction. Namely, for fixed distance $r$ and
large scale $\gamma\gg 1$ ($m\ll 1$) we are back to the logarithmic repulsion
(shifted by a constant). At fixed $\gamma$, however, the interaction range of
the 2D Yukawa potential is much shorter and decays exponentially, while still
being logarithmic at short distances.

In order to match with the point process (2.3) of the 2D Coulomb gas at a
given inverse temperature $\beta$, we use the folloqing shifted potential
(4.3):

(4.6) $\displaystyle V_{\rm Yukawa}^{(2)}(r)$ $\displaystyle=$
$\displaystyle-\beta\left(K_{0}\left(r\gamma^{-1}\right)-\log(2\gamma)+C\right)\
,$

identifying $\frac{g^{2}}{2\pi}=\beta$. A comparison between the two
potentials at fixed $\beta=1$ is shown in Figure 8 left.

Potential Nearest Neighbour Next-to-Nearest Neighbour

Distance

Figure 8. Left: Comparison of the 2D Coulomb potential ($\gamma=\infty$,
dashed) and the the 2D Yukawa potential (4.6) varying from $\gamma=0.01$
(bottom red full line) to $\gamma=10$ (top blue full line), see also inset in
the right plot for the colour coding. Middle and right: numerically determined
NN and NN spacing distributions, respectively, for the same values of
$\gamma$. The maximum increases with $\gamma$ and moves from left to right
towards the rightmost curve, given by the 2D Coulomb potential.

The NN (NNN) spacing distributions shown in Fig. 8 middle (right) have to be
determined numerically again, as it was done for the 2D Coulomb interaction in
[9]. That is, we use a Metropolis-Hastings algorithm [50, 49] with (4.6) as
the potential. The points are initialised independently, and then perturbed
iteratively. Each perturbation is always accepted if it leads to lower energy,
and accepted with probability $e^{-(V_{\rm after}-V_{\rm before})}$ if it
leads to higher energy, otherwise it is rejected. (Note that by including
$\beta$ in Equation (4.6), we have made it dimensionless.) After a number of
iterations, the points are considered to be a sample of the potential. On
average, each point is perturbed 100 times. This is considered enough, as
increasing the amount of iterations does not change the NN spacing
distribution. Unfolding is necessary when computing the NN spacing
distribution of the Yukawa potential for $\gamma<\infty$, as the global
spectrum is non-uniform.

While the potential $V_{\rm Yukawa}^{(2)}(r)$ differs considerably for the
parameter values chosen, the spacing distribution converges rather rapidly to
the one of the Coulomb potential at $\gamma\to\infty$, in particular for NN.
The reason is that it probes local correlations among the points (which are
closer for NN than for NNN), which seems to be rather robust under
deformations of the potential. As an aside, such a universality has been
proven rigorously in 1D for all fixed values of $\beta$ for deformed
potentials that behave only locally as a logarithm [51].

The fitted value of the scale $\gamma$ could in principle be translated into
an actual correlation length, but the effect of lowering the parameter
$\gamma$ is qualitatively similar to lowering $\beta$, and only becomes
noticeable at relatively small values of $\gamma$ for the NN. Additionally, in
order to compare our data with the spacing distribution we first need to
unfold. This makes a translation back into a real correlation distance that
may depend on the terrain (local density) rather cumbersome. The two-parameter
fit with the Yukawa potential does not reconcile the discrepancy between
different $\beta$-values obtained for the NN and NNN distributions, and we
therefore did not pursue the approach further. In other words, adding a length
scale does not allow us to fit both NN and NNN with the same parameter values.
Recall that the actual value of $\beta$ does not have a biological meaning. It
rather serves as a relative measure in comparing different species, and the
sole purpose of adding a length scale was to compensate for differences in the
fitted $\beta$ for NN and NNN distributions. Fitting both NN and NNN at the
same time would not improve this, as we would merely observe an average of the
two distributions. Because of the relatively small amount of available data,
especially for the Eagle Owls, it is probably not correct to weight the
contributions from NN and NNN equally in determining $\beta$ and $\gamma$.
However, any other choice of weighting them easily becomes arbitrary.

###  4.3. Independence of points: Reuse of nests

One of the key assumptions in deriving the Poisson NN and NNN spacing
distributions in (2.1) and (2.2) is the independence of points in this
process. In this subsection we investigate if in the absence of correlation
(repulsion), which we found in some data sets by agreeing with Poisson, it is
justified. The reason for this question is that some of the nests are reused,
not only within one species but also between different species, and we will
present data about this fact below. Such an analysis is possible because all
old and new nests were marked precisely with GPS coordinates, up to an error
of about 10 m.

Buzzards Goshawks Eagle Owls
Probability of Reuse

Figure 9. Percentage of nests that have been reused per year by each species.
The black crosses indicate the fraction of nests that were reused from the
previous year by the same species, i.e. Common Buzzards reusing a Common
Buzzard nest in the left plot. The black full line indicates the mean over the
entire period of observation. The red circles indicate the fraction of nests
that one species reused from all nests occupied by any of the 3 species in any
previous year. i.e. a Common Buzzard reusing a nest from a Common Buzzard,
Goshawk or Eagle Owl from any previous year in the left plot. The red dashed
line again indicates the mean value. Notice that both percentages may
sometimes coincide.

The amount of reusage of occupied nests from one year to the next year varies
between 30% (Common Buzzards), 45% (Goshawks) and 65% (Eagle Owl) within one
species, on average. The reuse of any nest from any previous year is about
10-20% higher. The frequent re-use of nest sites over years and even decades
has been documented previously for the study area and is typical of avian
predators [44]. Habitat heterogeneity in many dimensions is commonly observed
in these species [40], and individual preferences for certain habitat features
are also common [44], leading to frequent re-use of nests. In other species,
the same nest has been in use for centuries [45] and even millennia [46].

What is the consequence for our data analysis? In Section 3 we treated every
set of occupied nests per year as an independent ensemble, in order to compare
it to a statistical mechanics model, that represents an equilibrium
configuration. The above data show that this is a simplification, that is only
partly justified. The fact that a substantial fraction of nests is reused,
both for the repulsion within one species or among two different species,
means that not all NN and NNN positions are new.

In order to improve our statistics we have then made an ensemble average over
10 consecutive years (in [9] over 5 years). We have not analysed if after two
or more years the positions become more independent or randomised, but it
seems reasonable to assume that. Our situation could be compared to ensemble
generation in statistical mechanics. If a new configuration is generated from
an old one, e.g. with the Metropolis–Hastings algorithm used to generate the
2D Coulomb and 2D Yukawa ensembles in this paper, one has to monitor the
autocorrelation time before a new configuration can be accepted as
independent. In the generation of points, we could simply let enough
iterations go by, but for the observed data we cannot wait one or more years,
until the next set of occupied nests is more independent, without losing much
information (and statistics). This is because the repulsion within one or more
species can change over time, as we have observed.

There are ecological reasons why certain locations of nests are more popular
than others, which is why they get reused. Plus it takes time and effort to
build a new nest. Certain nest sites are associated with higher reproductive
success, independent of the individuals using these nests [44]. The proximity
of high quality hunting habitat, less exposure to predators and parasites, and
less disturbance from humans are just some of the key factors [44, 15].

The fact that there is a large number of used nests in the observed area seems
to guarantee at least to some extent, that a kind of "random sampling" can
take place between consecutive years. There an estimate on the total number of
500-600 reusable nests in the area for birds of prey.

As a final remark, there seems to be no correlation between the population
growth and the percentage of reusage of nests, which is approximately constant
in Fig. 9. As we have seen in Fig. 4 for the Common Buzzard, their population
has grown from 50 pairs to over 250, and from a single pair to over 20 for all
Eagle Owls in total, see Fig. 5. The Goshawk population has also grown but
more slowly, and then declined again. Together with the high number of empty
nests available each year this could be argued to be in favour of a certain
independence of the locations of nests each year, despite the high fraction of
reusage.

##  5\. Dicussion and Open Questions

This study has shown that a simple, one-parameter model from statistical
mechanics, without any biological background, nevertheless has surprisingly
well explained important features in this three-species guild of avian
predators. We have shown that the intraspecific repulsion measured through
$\beta$ clearly deviates from Poisson statistics. The relative ranking in
repulsion strength amongst the three species, averaged over all years, has
been found compatible with experience from ecology.

The correlation between $\beta$ and the time dependence of the population have
been found plausible for at least two of the three species analysed. For the
Common Buzzard the steep rise in population from about 50 to 250 pairs has
been paralleled by a clear rise in $\beta$, in particular in the NNN spacing,
showing also an increased interaction range. For the Goshawk the population
has been more stable around 20 pairs, with a decline after a peak around 10
years ago. Both NN and NNN show a clear decrease in the fitted $\beta$. In
contrast for the Eagle Owls, with the population increasing from a single pair
to over 15, we have seen a steep decrease in $\beta$ for NN and no repulsion
between NNN. This relation is not clear to us and may be due to the limitation
of the model, or an artefact due to very low statistics. The change in $\beta$
over time might hence be a fruitful further topic for investigation.

At the same time, the results also clearly indicate some limitations when it
comes to interspecific interactions. We have found the repulsion between Eagle
Owls in the forest on the one hand, and the Common Buzzard respectively
Goshawk on the other hand to be significantly smaller than the respective
intraspecific repulsion, which is plausible. Both interspecies repulsions have
values of $\beta=0.15-0.1$ above Poisson at $\beta=0$. Taken together with the
effect of the terrain, effectively lowering the dimension of the Poisson
distribution, this makes the repulsion found more pronounced and significant.
For the Eagle Owls in the plains the same order of repulsion was only observed
with the Common Buzzard, and no repulsion was seen for the Goshawk, neither
among the latter two. In view of the continuing growth of population in Common
Buzzards and Eagle Owls, notably populating the plains for the latter, and the
population of Goshawks being under pressure, this seems to be a limitation of
the model (and present statistics).

From the point of view of the statistical mechanics model, we have shown that
a simple formula based on a surmise from random matrix theory allows to
quantify and discuss the intra- and interspecies NN repulsion as well, instead
of using the numerically generated spacings distributions from a full 2D
Coulomb gas, that also include NNN. Methodologically, we have quantified the
effect of the terrain. In particular, it is not responsible for the repulsion
observed. We have also discussed the 2D Yukawa interaction and spacing ratios
as alternative models. The extensive reusage of nests, especially among Eagle
Owls, has put the assumption of independence of ensembles as realised by
consecutive years under scrutiny. However, at least from a random matrix
perspective we know that there, spectral statistics shows a certain robustness
or universality under such perturbations.

We feel that this example of model transfer between disciplines has been an
insightful exercise for both disciplines involved. It remains to be seen
whether interspecific interactions might be better captured with another
decade of data with increased statistics, and with all three species at their
ecological carrying capacity.

Acknowledgements: We acknowledge funding by the German Research Foundation
(DFG) through grant CRC 1283/2 2021 "Taming uncertainty and profiting from
randomness and low regularity in analysis, stochastics and their applications"
(GA and PP), and CRC/Transregio 212 "A Novel Synthesis of Individualisation
across Behaviour, Ecology and Evolution: Niche Choice, Niche Conformance,
Niche Construction (NC)3" (OK). AM would like to thank Kurt Mielke for
improving the speed of the Metropolis-Hastings code and providing a server for
it, which made both the many parameter combinations of the Yukawa potential
and larger $N$ feasible. We thank Michael Baake for useful discussions and
comments.

##  Appendix A Alternative to Unfolding: Complex Spacing Ratios

In this appendix we briefly discuss complex spacing ratios as introduced in
[4] and compare a subset of our data to several moments of these. However, the
outcome is not conclusive, as we will see.

Complex spacing ratios have become a popular tool in analysing two-dimensional
data sets, compared to NN or NNN spacing distributions in radial direction.
The advantage compared to the latter is that, first, no unfolding is
necessary, which is often more difficult in two compared to one dimension.
Second, it also gives both angular and radial information about the
interaction between the points. It is defined as follows. Suppose we have an
eigenvalue at $z_{k}$, with NN at $z_{k}^{\rm NN}$, and NNN at $z_{k}^{\rm
NNN}$, in radial distance. The complex spacing ratio is then defined as

(A.1) $u_{k}=\frac{z_{k}^{\rm NN}-z_{k}}{z_{k}^{\rm NNN}-z_{k}}\ .$

The goal is then to determine the probability distribution
$\rho(u=re^{i\theta})=\rho(r,\theta)$ of this spacing ratio, in the limit of a
large number of particles $N\to\infty$ in the bulk of the spectrum, for a
given point process. For the Poisson process, it is known to be flat on the
unity disc [4],

(A.2) $\rho_{\rm Poi}(u)=\frac{1}{\pi}\Theta(1-|u|)\ .$

For the complex Ginibre ensemble explicit $(N-1)$-fold integral representation
are known [4], as well as approximate expressions that converge very rapidly
[32].

There are several disadvantages in our situation, however. First, the amount
of data we have available is rather small - even if very large for biological
standards. This makes even a qualitative comparison to 2D plots very hard, if
not impossible.

The difficulty to make 2D-fits was anticipated in [4], and thus alternatively
integrals over the angle or the radius where proposed,

(A.3) $\rho(r)=\int_{0}^{2\pi}d\theta r\rho(r,\theta)\
,\quad\rho(\theta)=\int_{0}^{1}drr\rho(r,\theta)\ ,$

and then to consider moments thereof. Here, we assume that the limiting
support of complex eigenvalues is normalised to the unit disc, as in Section
2.2 for the general 2D Coulomb gas. Analytic and approximate formulas have
been worked out both for the Poisson case, and the complex Ginibre ensemble,
which we reproduce here for completeness. Inserting the expression (A.2) for
the Poisson ensemble into the definition (A.3), we obtain [4]

(A.4) $\rho_{\rm Poi}(r)=2r\ ,\quad\rho_{\rm Poi}(\theta)=1/(2\pi)\ .$

This leads to the following expression for the first moments in the Poisson
case:

(A.5) $\displaystyle\begin{split}\langle r\rangle_{\rm
Poi}=&\int_{0}^{1}dr\,r\rho_{\rm Poi}(r)\ =\ \frac{2}{3}\ ,\\\ \langle
r^{2}\rangle_{\rm Poi}=&\int_{0}^{1}dr\,r^{2}\rho_{\rm Poi}(r)\ =\
\frac{1}{2}\ ,\\\ \langle\cos(\theta)\rangle_{\rm
Poi}=&\int_{0}^{2\pi}d\theta\cos(\theta)\rho_{\rm Poi}(\theta)\ =\ 0\ ,\\\
\langle\cos^{2}(\theta)\rangle_{\rm
Poi}=&\int_{0}^{2\pi}d\theta\cos(\theta)^{2}\rho_{\rm Poi}(\theta)\ =\
\frac{1}{2}\ ,\\\ \langle r\cos(\theta)\rangle_{\rm Poi}=&\
\int_{0}^{2\pi}d\theta\int_{0}^{1}drr\,r\cos(\theta)\rho_{\rm Poi}(r,\theta)\
=\ 0\ .\end{split}$

Each moment is normalised. The angular and radial integrals decouple in all
cases. The moments for complex Ginibre are found in [32] and reproduced here
in Table 1 for completeness.

limiting moment | Ginibre | Poisson
---|---|---
$\langle r\rangle$ | 0.739 | 2/3
$\langle r^{2}\rangle$ | 0.581 | 1/2
$\langle\cos(\theta)\rangle$ | -0.247 | 0
$\langle\cos(\theta)^{2}\rangle$ | 0.450 | 1/2
$\langle r\cos(\theta)\rangle$ | -0.189 | 0
Table 1. Numerically determined moments in the complex Ginibre ensemble from
[32], rounded to 3 digits (middle column), and the corresponding values for
Poisson from (A.5) [4].

For the 2D Coulomb gas with intermediate values $0<\beta<2$, to our knowledge
no such relations have been worked out, in particular if the corresponding
values are in between the Poisson and Ginibre values. Furthermore, the effect
of the interaction range not being of Coulomb type on the moments is not
clear. These are crucial for our quantitative comparison.

In Figure 10 we compare the moments of the complex spacing rations from
Poisson and complex Ginibre ensemble to those of the spacings among nests of
Common Buzzards (left) and Goshawks (right), for which we have the most data.
As in Section 3 we use a time moving average.

Common Buzzards Goshawks

Figure 10.  The moments from complex spacing ratios in (A.5) for Poisson and
from Table 1 for complex Ginibre [32] are compared to the corresponding
moments for Common Buzzards (left) and Goshawks (right). For a better
comparison we have put the moments where the Ginibre (Poisson) value is higher
in the top (bottom) row.

In the comparison Fig. 10 for the Common Buzzards (left) the top row with the
radial moments shows a clear tendency over time, where the value for the
moments increases from above Poisson towards Ginibre. This is consistent with
the findings from Section 3.1 for the fitted $\beta$-values for NN and NNN
spacing distribution for Common Buzzards, see Figure 4 left, where an increase
was seen for NN (NNN) from about 0.4 to 0.7 (0.0 to 0.8). In contrast, the
angular moments in the bottom row are less conclusive. While the left and
right plots are close to Poisson and show a decrease towards the negative
value for Ginibre, the middle plot is scattered somewhat between Poisson and
Ginibre and rather decreasing towards Poisson.

For the Goshawks the picture is even less clear from the complex spacing
ratios. While the radial moments in the top row show values above Ginibre
(which is above Poisson), with a slight tendency going down towards the
Ginibre value, the bottom plots with angular moment (left and right plot) are
closer to Poisson, with no clear trend. The bottom middle plot is quite
scattered, starting in between Poisson and Ginibre and then moving clearly
below the Ginibre value. In comparison, in Section 3.1 for the fitted
$\beta$-values for NN and NNN spacing distribution for Goshawks, see Figure 4
right, a clear decent from $\beta$-values from about 1.5 down for 0.5 for both
NN and NNN was detectable.

In the comparison of $\beta$-values within one species and amongst different
species, it was important to have even approximate values for $\beta$, that
quantify the repulsion. With the predictions from moments of complex spacing
rations such a quantitative comparison does not seem to be easily possible, at
least not with the predictions we have at hand. Second, we cannot exclude an
influence of the spherical geometry on the predicted complex spacing ratios, a
situation we do not have for our observed data. This seems to be indicated by
the unclear trend form the angular momenta. The spacing distributions we use
in the main text are local quantities which a priori seem to be less sensitive
to the geometry.

##  Appendix B Comparison surmise and Coulomb gas NN spacing for small $\beta$

Figure 11.  Comparison of the surmise (2.4) with (2.5) (dashed red line), and
the numerically determined NN spacing distribution from the 2D Coulomb gas [3]
(full blue line).

In this appendix we illustrate the deviation between the numerically
determined NN spacing distribution from [3] and the surmise [26] given in
(2.4) together with (2.5), for small values of $\beta$. The deviation
quantified in terms of the standard deviation and the Kolmogorov-Smirnov
distance can be found in [26, Table 1.].

The fits in $\beta$ from the 2D Coulomb gas to the data presented in Section 3
are done in steps of $0.1$ in $\beta$, except for an intermediate step at
$\beta=0.15$ for small $\beta$, see Fig. 11 top middle. In contrast, the fits
to the surmise in Section 3 are made with continuous values for $\beta$, as it
becomes apparent form the values given in the respective insets. This leads to
an additional deviation between the two, which is why we display them at equal
values for comparison here.

## References

  * [1] R. Bardenet, A. Feller, J. Bouttier, P. Degiovanni, A. Hardy, A. Ranc̨on, B. Roussel, G. Schehr, and C. I. Westbrook, From point processes to quantum optics and back. Preprint arXiv:2210.05522, 2022.
  * [2] R. Hamazaki, K. Kawabata, N. Kura, M. Ueda, Universality classes of non-Hermitian random matrices. Phys. Rev. Research 2, 023286 (2020) [arXiv:1904.13082].
  * [3] G. Akemann, M. Kieburg, A. Mielke, T. Prosen, Universal signature from integrability to chaos in dissipative open quantum systems. Phys. Rev. Lett. 123, 254101, 1–6 (2019) [arXiv:1910.03520].
  * [4] L. Sá, P. Ribeiro, and T. Prosen, Complex Spacing Ratios: A Signature of Dissipative Quantum Chaos. Phys. Rev. X 10, 021019 (2020). [arXiv:1910.12784].
  * [5] S. Torquato and F. H. Stillinger, Local density fluctuations, hyperuniformity, and order metrics. Phys. Rev. E, 68(4), 041113 (2003).
  * [6] F. Lavancier, J. Møller, and E Rubak, Determinantal point process models and statistical inference. Royal Stat. Soc. B 77, 853–877 (2015) [arXiv:1205.4818].
  * [7] P. J. Clark and F. C. Evans, Distance to nearest neighbour as a measure of spatial relationships in populations. Ecology 35, 445-453 (1954).
  * [8] R. Law, J. Illian, D. F. Burslem, G. Gratzer, C. Gunatilleke, and I. Gunatilleke, Ecological information from spatial patterns of plants: insights from point process theory. J. Ecology 97(4), 616–628 (2009).
  * [9] G. Akemann, M. Baake, N. Chakarov, O. Krüger, A. Mielke, M. Ottensmann, and R. Werdehausen, Territorial behaviour of buzzards versus random matrix spacing distributions. J. Theoretical Biology 509, 110475:1–7 (2021) [arXiv:2003.090204].
  * [10] M. Begon, C. R. Townsend, Ecology: From individuals to ecosystems. 5th Edition. Wiley-Blackwell, Hoboken, New Jersey, USA, 2021.
  * [11] O. Krüger, Interactions between Common Buzzard Buteo buteo and Goshawk Accipiter gentilis: trade-offs revealed by a field experiment. Oikos 96, 441–452 (2002).
  * [12] N. Chakarov, and O. Krüger, Mesopredator release by an emergent superpredator: a natural experiment of predation in a three level guild. PloS one 5(12), e15229 (2010).
  * [13] A. K. Müller, N. Chakarov, H. Heseker, and O. Krüger, Intraguild predation leads to cascading effects on habitat choice, behaviour and reproductive performance. J. Animal Ecolology 85, 774–784 (2016).
  * [14] O. Krüger, N. Chakarov, J.T. Nielsen, V. Looft, T. Grünkorn, B. Struwe-Juhl, and A. P. Møller, Population regulation by habitat heterogeneity or individual adjustment? J. Animal Ecology 81, 330–340 (2012).
  * [15] N. Chakarov, M. Boerner, and O. Krüger, Fitness in Common Buzzards at the cross–point of opposite melanin-parasite interactions. Functional Ecolology 22, 1062–1069 (2008).
  * [16] N. Chakarov, M. Pauli and O. Krüger, Immune responses link parasite genetic diversity, prevalence and plumage morphs in common buzzards. Evolutionary Ecology 31, 51–62 (2017).
  * [17] A.K Mueller, N. Chakarov, O. Krüger O and J.I. Hoffman, Long-term effective population size dynamics of an intensively monitored vertebrate population. Heredity 117, 290–299 (2016).
  * [18] O. Bohigas, M. J. Giannoni, and C. Schmit, Characterization of chaotic quantum spectra and universality of level fluctuation laws. Phys.Rev. Lett. 52, 1 (1984).
  * [19] O. Bohigas, M. J. Giannoni, and C. Schmit, Spectral properties of the Laplacian and random matrix theories J. Physique Lettres 45, 1015–1022 (1984).
  * [20] G. Casati, F. Valz-Gris, and I. Guarneri, Connection between quantization of nonintegrable systems and statistical theory of spectra. Lett. Nuovo Cimento Soc. Ital. Fis. 28, 279 (1980).
  * [21] M. V. Berry and M. Tabor, Level clustering in the regular spectrum. Proc. Roy. Soc. A 356, 375 (1977).
  * [22] T. Guhr, A. Müller-Groeling, and H. A. Weidenmüller, Random-matrix theories in quantum physics: Common concepts. Phys.Rep. 299, 189–425 (1998) [arXiv:cond-mat/9707301].
  * [23] M. L. Mehta. Random Matrices, volume 142 of Pure and Applied Mathematics. Elsevier/Academic Press, Amsterdam, 3rd edition, 2004.
  * [24] F. Haake, Quantum Signatures of Chaos, 3rd Edition, Springer, Heidelberg, 2010.
  * [25] R. Grobe, F. Haake, H.-J. Sommers, Quantum distinction of regular and chaotic dissipative motion. Phys. Rev. Lett. 61, 1899–1902 (1988)
  * [26] G. Akemann, A. Mielke, and P. Päßler, Spacing distribution in the 2D Coulomb gas: Surmise and symmetry classes of non-Hermitian random matrices at non-integer $\beta$. Phys. Rev. E 106, 014146 (2022) [arXiv:2112.12624].
  * [27] P.J. Forrester, Log-Gases and Random Matrices. Princeton University Press, Princeton, 2010.
  * [28] A. Kuijlaars, Universality, chapter 6 in G. Akemann, J. Baik, and P. Di Francesco (eds.), The Oxford Handbook of Random Matrix Theory, 1st ed., Oxford University Press, Oxford (2011) [arXiv:1103.5922v2].
  * [29] Y. Ameur, H. Hedenmalm, and N. Makarov, Fluctuations of eigenvalues of random normal matrices. Duke Math. J. 159, 31–81 (2011) [arXiv:0807.0375].
  * [30] H. Markum, R. Pullirsch, and T. Wettig, Non-Hermitian random matrix theory and lattice QCD with chemical potential. Phys. Rev. Lett. 83, 484–487 (1999) [arXiv:hep-lat/9906020].
  * [31] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007) [arXiv:cond-mat/0610854].
  * [32] I.G. Dusa, and T. Wettig, Approximation formula for complex spacing ratios in the Ginibre ensemble. Phys. Rev. E 105, 044144 (2022) [arXiv:2201.05104].
  * [33] Y. Atas, E. Bogomolny, O. Giraud, P. Vivo, and E. Vivo, Joint probability densities of level spacing ratios in random matrices, J. Phys. A: Math. Theor. 46, 355204 (2013) [arXiv:1305.7156].
  * [34] S. Serfaty, Microscopic description of Log and Coulomb gases, in: Borodin, A., Corwin, I, Guionnet, A. (eds.), Random Matrices, AMS, Providence, RI, 341–387 (2019). [arXiv:1709.04089].
  * [35] G. Lambert, Poisson statistics for Gibbs measures at high temperature. Annales de l’Institut Henri Poincaré, Probabilités et Statistiques. 57(1) Institute of Mathematical Statistics (2021) [arXiv:1912.10261].
  * [36] P. Vivo, and S.N. Majumdar, On invariant $2\times 2$ $\beta$-ensembles of random matrices Physica A 387, 4839 (2008) [arXiv:0706.2476].
  * [37] J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices. J. Math. Phys. 6, 440–449 (1965).
  * [38] G. Akemann, M. J. Phillips, and L. Shifrin, Gap probabilities in non-Hermitian random matrix theory. J. Math. Phys. 50, 063504, 1–32 (2009) [arXiv:0901.0897].
  * [39] J. Fischmann, W. Bruzda, B. A. Khoruzhenko, H.-J. Sommers, and K. Zyczkowski. Induced Ginibre ensemble of random matrices and quantum operations. J. Phys. A 45(7), 075203:31 (2012) [arXiv:1107.5019].
  * [40] O. Krüger, and J. Lindström, Habitat heterogeneity affects population growth in goshawk Accipiter gentilis. J. Animal Ecolology 70, 173–181 (2001).
  * [41] O. Krüger, and U. Stefener, Nahrungsökologie und Populationsdynamik des Habichts (Accipiter gentilis) im östlichen Westfalen. Vogelwelt 117, 1–8 (1996)
  * [42] S. Cramp and K.E.L Simmons, Vol. II: Hawks to bustards, Handbook of the birds of Europe, the Middle East and North Africa : the birds of the Western Palearctic, Oxford University Press, Oxford, 1980.
  * [43] O. Krüger, Long-term demographic analysis in goshawk Accipiter gentilis: the role of density-dependence and stochasticity. Oecologia 152, 459–471 (2007).
  * [44] O. Krüger, Analysis of nest occupancy and nest reproduction in two sympatric raptors: common buzzard Buteo buteo and goshawk Accipiter gentilis. Ecography 25 523–532 (2002).
  * [45] D.H. Ellis, T. Craig, E. Craig, S. Postupalsky, C.T. LaRue, R.W. Nelson, D.W. Anderson, C.J. Henny, J. Watson, B.A. Millsap, J.W. Dawson, K.L Cole, E.M. Martin, A. Margalida and P. Kung, Unusual raptor nests around the world. J. Raptor Research 43, 175–198 (2009).
  * [46] K.K. Burnham, W.A. Burnham, and I. Newton, Gyrfalcon Falco rusticolus post-glacial colonization and extreme long-term use of nest sites in Greenland. Ibis 151, 514–522 (2009).
  * [47] M.P. Peskin, and D.V. Schroeder, An Introduction To Quantum Field Theory. CRC Press, Boca Raton, 1995.
  * [48] A. Jeffrey and D. Zwillinger (eds.) Gradshteyn and Ryzhik’s Table of Integrals, Series, and Products. Academic Press, New York, sixth edition, 2000.
  * [49] W.K. Hastings, Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97–109 (1970).
  * [50] D. Chafaï, and G. Ferré, Simulating Coulomb gases and log-gases with hybrid Monte Carlo algorithms, J. Stat. Phys. 174, 692–714 (2019) [arXiv:1806.05985].
  * [51] M. Venker, Particle Systems with Repulsion Exponent $\beta$ and Random Matrices. Electron. Commun. Probab. 18(83), 1–12 (2013) [arXiv:1209.3178].