Interactions between different birds of prey as a random point process

The two-dimensional (2D) Coulomb gas is a one-parameter family of random point processes, depending on the inverse temperature β. Based on previous work, it is proposed as a simple statistical measure to quantify the intra- and interspecies repulsion among three different highly territorial birds of prey. Using data from the area of the Teutoburger Wald over 20 years, we fit the nearest-neighbour and next-to-nearest neighbour spacing distributions between the respective nests of the goshawk, eagle owl and the previously examined common buzzard to β 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 2D Yukawa interaction, and the statistical independence of the time moving average we use for the yearly ensembles of occupied nests. We also check that an artificial random displacement of the original nest positions of the order of the mean level spacing quickly destroys the repulsion measured by β > 0. A simple, approximate analytical expression for the nearest-neighbour spacing distribution derived from non-Hermitian random matrix theory proves to be very useful.


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,5], condensed matter physics [6], statistics [7] and ecology [8,9,10].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.[5,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 [11].Animals compete for limited resources with individuals of their own (intraspecific competition) or another species (interspecific competition), see [12].One of the easiest and hence most often used estimates of competition is the nearest neighbour distance [8].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 [13,14].
Two examples where 2D point processes have been applied in biology are the spacial distribution of trees, see [9] for a review, and the distribution of nests of birds of prey in space and time in our previous work [10].Using data from [15,16,14] on the yearly locations of occupied nests of the Common Buzzard in an area of 300 km 2 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.[7].The single fit parameter used is the coupling strength β, 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 β, over the observed period of 20 years, as a function of the increasing population density.Let us emphasise that the fitted value of β does not have a direct biological meaning, but merely serves to quantify the repulsion, in particular the deviation from 2D Poisson statistics at β = 0.
Currently, we do not have an underlying microscopic model that would lead to a NN spacing distribution between the nests of birds of prey of (approximately) 2D Coulomb type.In 1D, a model for the spacing between parked cars and birds on a power line has been developed in [17], leading to a Gamma-distribution.It differs from Hermitian random matrix theory with a logarithmic Coulomb interaction in 1D, to be described below.Ref. [17] is based on independently distributed symmetric spacings and provides a link to the perception of space by humans and birds.Because the absolute ordering into consecutive spacings between cars (birds) in 1D is absent in 2D, a simple generalisation of such an approach is not possible.Instead, our argument for using the spacing distribution of the 2D Coulomb gas is based on universality, as it has been observed to apply in several very different examples in physics: in dissipative quantum spin chains [2,4,5], the dissipative kicked rotor [2,3] and Quantum Chromodynamics with quark chemical potential [18].This is of course a strong assumption for the biological system, but it has worked surprisingly well for the Common Buzzard [10].Other authors make different strong mathematical assumptions when modelling point processes from biology, for example that they are determinantal [7].The data are then fitted by a presumed correlation kernel, which relates to the two-point function.A plethora of models have been formulated to explain intra-and interspecific competition and density dependence in ecological time series, and they commonly make a number of ecological assumptions.Here, we test how powerful a simple approach from statistical mechanics might be, that does not carry any ecological assumptions with it.
A prime example from physics for a transition from Poisson to correlated random variables is the Bohigas-Giannoni-Schmit or quantum chaos conjecture [19,20,21].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 eigenvalues are real (1D), whereas non-Hermitian matrices have complex eigenvalues in 2D.In both cases the 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 [24], 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 evidence has been collected, cf.[22] and references therein.It is based for example on 740 spacings in the Sinai billiard, and 1726 spacings from nuclear scattering data.A complete analytical picture based on a semi-classical expansion has taken over 25 years to be developed, cf.[25] and references therein.The quantum chaos conjecture has been extended to 2D for dissipative quantum systems [26], and so far numerical evidence has been presented, e.g. from the spectrum of the corresponding Liouville operator [2,4,5].
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 β = 1, 2, 4 occurring in 1D, and β = 2 in 2D1 .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,28].Their mathematical universality in the limit of a large number of points N (matrix dimension) has been proved rigorously in 1D, see [29] and references therein, and in 2D, cf.[30] and [4] 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 β = 0, and the specific β-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 [4] to directly use the 2D Coulomb gas at intermediate values of β, and to determine the NN and NNN spacing distributions numerically, that can then be compared with data.This is reviewed in Section 2.2, in particular a corresponding approximate surmise for the 2D NN spacing distribution valid for small values of β [27], cf.Appendix B. These NN and NNN spacings were then used for fits as a function of β, for the spacing of eigenvalues of the Liouville operator of dissipative, boundary driven systems in 2D [4], and in the yearly spacing distributions of nests of the Common Buzzard [10].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 extend the analysis of a single species [10] to two further species, and to quantify their respective interactions.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 [16,14,31,32].Our goal is to try to quantify the interspecies interaction between each two species, and to compare it to the intraspecies interaction found within each species individually.Therefore, in a first step we have repeated the analysis from [10] for Goshawks and Eagle Owls individually, see Section 3.1.This includes the dependence on the respective change in population over 20 years, ranging between 12-29 and 1-23 pairs, respectively.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.Due to a comparably low statistics we have been unable to address the time dependence for the interspecific interaction between different species so far, to be discussed in Section 3.2.
While this study might be perceived as rather limited with regard to the sample sizes and time spans involved, each species-specific time series with complete coverage of over 20 years in a study area of 300 km 2 is already amongst the largest data sets for any predatory bird [15], and the inclusion of three species from the same area over such a time-span and detail is probably unique in the world.Over the course of the study, 3377 breeding attempts of Common Buzzard have been documented, over 423 Goshawk breeding attempts and 174 Eagle Owl breeding attempts, cf.Appendix D.
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 [10].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 β ̸ = 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 [10] the NN and NNN spacing distributions of nests were found to be described by different (typically decreasing) values of β, 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. [18,4].As an alternative quantity that does not require unfolding, spacing ratios have been introduced in 1D [33] and 2D [5].However, they are only known analytically for Poisson and for random matrix statistics, the complex Ginibre ensembles in 2D [5,34].In contrast to 1D [35], where the spacing ratio was derived as a function of β, the interpolation in 2D is not so clear, as we will see from our data analysis presented in Appendix A. In Appendix C we check that a random displacement of the observed positions of the nest by the order of the mean level spacing destroys the repulsion in the NN spacing and quickly leads to a Poisson distribution.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 km 2 (8 • 25'E and 52 • 6'N) is located in Eastern Westphalia.It consists of two 125 km 2 grid squares and 50 km 2 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 km 2 .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 [15,16,14].
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 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.The extent of the monitored area is roughly 12 × 25 km.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.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 and Appendix D Table 2 for their annual 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, cf.Table 2.
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).
The mean spacing in kilometers, e.g.m B = 0.63 km for the Common Buzzard and m G = 2.63 km for the Goshawk in 2020, is very different and encodes important biological information about the range of interaction within one (or different) species.Goshawks are among the most territorial of all European birds of prey and have large territory sizes of commonly between 5 and 50 km 2 and regularly make foraging trips of 5-10 km distance [36].Interactions between neighbouring territories are rather aggressive and intruders into the territory are regularly killed [36].As the mean nearest neighbour distance in our study site m G = 2.63 km is very low compared to other study sites due to our high population density [15], regular antagonistic interaction between the territorial pairs is exceedingly likely.However, 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 [18,4] for different procedures.Here, we will apply the method proposed in [4], which goes as follows.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 ρ ave of points in 2D.The measured spacing around point z = (x 0 , y 0 ) is then normalised by rescaling with ρ 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.

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 nextto-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 will mainly discuss the spacing distributions in the following.However, for fits we will exclusively use the cumulative distribution function (CDF), which is independent of the choice of binning into histograms for the data.
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 : see e.g.[25,5] 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 → ∞, after rescaling the mean spacing between points as s = R π/(4N ), eqs.(2.1) and (2.2) are obtained in units of s.Despite resulting from uncorrelated points, the 2D area measure in polar coordinates, dxdy = sds dθ for z = x + iy = se iθ , 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 constant uniform density on a disc can also be obtained for a 2D Coulomb gas.Let us emphasise that the above spacing distributions quantify (the absence of) local correlations amongst points at distance ∼ 1/ √ 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 β = (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, where Z N,β is the normalising partition function.The N points are represented by complex coordinates z j ∈ C, j = 1, . . ., N , with the usual identification C ∼ R 2 .Compared to standard conventions, where β multiplies the entire Hamiltonian including the confining potential, we have absorbed β in front of the potential by rescaling the coordinates βV (z) → V (z).The rescaling is made in order to allow us to take the limit β → 0 to reach Poisson statistics, while maintaining a confining potential.Furthermore, the charge is set to unity (or absorbed into β).
Such static Coulomb gases have been much studied in the recent mathematical literature, and we refer to [37] for a recent review.Other classes of repulsion are possible such as a power law decay also called Riesz gas, cf.[38] for mathematical results and references therein.Due to the universality of the 2D Coulomb gas mentioned in the introduction we will stick to this class of interactions.After a further rescaling of the potential V (z) → N V (z), in the large-N limit the mean density ρ(z) of points condenses on the so-called droplet S. It is given by the Laplacian of the potential, ρ(z) = 2 β ∂ z ∂ 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/ √ N very little is known analytically for fixed β > 0, apart from Poisson statistics at β = 0 (that extends to β ∼ 1/N , cf. [39]) and the integrable case β = 2, when the point process (2.3) becomes determinantal, cf.[28].Inspired by the Wigner surmise for the 1D Dyson gas, and its generalisation to general β based on a 2 × 2 β-ensemble [40], a surmise (sur) for the NN spacing distribution was derived from complex normal . Its behaviour at small values s → 0 is as expected heuristically from (2.3), with one power from the radial measure and a power β from the Vandermonde determinant.Unfortunately, it is well known [26] from the integrable case at β = 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 β → 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 β, rather than for β → 0. In order to improve the approximation of the surmise (2.4) in 2D to larger values of β (and with exact results for N > 2 being unavailable), an effective β eff was introduced in [27], by fitting a third order polynomial in β to the numerically determined spacing of the 2D Coulomb gas in the range of β ∈ [0, 3]:

Consequently in p
(NN) sur,β eff (s) the normalisation constant has to change too, α(β) → α(β eff ) = α eff , to ensure a normalised spacing and first moment equal to unity.This leads a reasonable approximation up to β ≈ 0.5, with a standard deviation of up to σ = 2.6 • 10 −2 .See Appendix B for plots in this range, and Fig. 2 right for a comparison at β = 2, where the surmise clearly fails.In [27] a more detailed comparison is presented, including higher β-values and their Kolmogorov-Smirnov distances.Obviously, the fit using an effective β eff spoils the heuristically expected proportionality at very small argument ∼ s 1+β eff .However, the limit β → 0 is still reproduced exactly.Being based on 2 × 2 matrices, there is no approximate prediction for the NNN spacing possible.
2 Notice a typo in the normalisation constant in [27]: α β there should be replaced by α 1+β/2 as here.
In our comparison to data in the next section we will use the CDF of (2.4) for the fitting procedure, in order to avoid any dependence on the choice of histograms.It is given by dt is the upper incomplete Gamma function.For completeness we give the analytical result for the NN [26] and NNN [42] spacing distribution for β = 2, that follows from the Ginibre ensemble [41].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, 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 Gin (s), to find an eigenvalue at the origin and the closest non-zero complex eigenvalue at radial distance s, For finite N the product extends only to N − 1.The limiting spacing distributions at infinite matrix dimension follow from differentiation, compare [26] for NN and [42] for NNN: dt is the lower incomplete Gamma function, and we give again the behaviour at s → 0, too.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 s1 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 ≈ 20, or larger.The NN with first moment of unity is then obtained by rescaling p(NN) Gin (s 1 y), and the correspondingly rescaled NNN spacing reads p(NNN) (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 [30].
In the Ginibre ensemble at β = 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 [43].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 (nests) as bulk points.
In the next section we will use both to the approximate NN spacing p sur,β eff (s) from the surmise [27], and the numerically determined NN and NNN spacing from [4] as functions of β, in steps of 0.1, and we refer to [4] for details about the numerics.In particular around 200000 spacings have been used there to generate these numerical distributions.

Interaction of species
3.1.Interaction within one species: Common Buzzard, Goshawk and Eagle Owl.Let us describe the procedure initiated in [10] 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 β 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.This has motivated us to consider time moving averages in the other columns, to be discussed below.In all figures we show both the fit using the Coulomb gas (full line), with the spacings determined numerically in [4] in discrete steps of size 0.1 for β, and the explicit formula (2.4) for the surmise (dashed line), with β 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 β-values vary over a wide range.For example, in the year 2020 in Fig. 3 left column, we find the strongest repulsion among Goshawks (β = 3), followed by the Common Buzzard (β = 1.1) down to Eagle Owls (F) with β = 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 β 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 β = 0. Let us emphasise, that the fitted parameter β 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.Without the successful application to the Common Buzzard in [10] we would not have embarked onto an analysis of data from other birds of prey, with a considerably lower statistics.The inference of biological information from blood samples of the Common Buzzard with low statistics has been analysed in [32].
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), cf.Table 2 in Appendix D. Between Common Buzzards and Goshawks there is a factor of 10 in abundance, and typically 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 [10] we choose a time moving average of 5 years, that gave a better resolution of the time dependence of the fitted β as a measure of repulsion.In [10] we also presented the fits for  2 in Appendix D for details.Shown are histograms and spacing distribution, but we emphasise that the fits were obtained using the cumulative distributions.Clearly when moving from individual years to time moving averages, the number of spacings for Goshawks and Eagle Owls (F) leads to an improved statistics and quality of the fits.
all individual years (compare Fig. 4 left in [10]).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 β-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 β-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 β value goes down for the Goshawks over the years, see Fig. 4. Notice that their population is peaked in 2012.For the Eagle Owls (F) we obtain a fitted value β = 0.1, being very close to 2D Poisson in this window of time average.Also here the averaged fitted β 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 β for all three species, as shown in Fig. 3 right column.The strongest repulsion is again amongst Goshawks (β = 0.8), followed by Common Buzzards (β = 0.6), and almost equally among Eagle Owls (F) with β = 0.5.The particularly high β for intraspecific repulsion in Goshawk makes a lot of ecological sense, as the Goshawk is one of the most territorial species in Europe [44].The Goshawk preys on larger prey and needs exclusive access to hunting areas [45].Therefore, intruding individuals are always met with high aggression and are regularly killed by the territory owners [46,44,47], cf.Fig. 1 for their relatively large spacings.All values are significantly above 2D Poisson statistics (β = 0).Because of the strong growth in population for Common Buzzards, and for Eagle Owls (F), see Figs. 4 and 5, respectively, and Table 2 Appendix D, the later years are weighted much higher here in the average over all years, because they contribute to more spacings per year.
In the following we will discuss our findings for the β-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 β-values obtained from NN and NNN fits do not agree in general, sometimes lying systematically below (Common Buzzards, Eagle Owls (F)) 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 [10] 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 β and γ) the phenomena of obtaining two different values for β from NN and NNN remains, see the discussion there.Therefore, we kept the simpler 2D Coulomb interaction as a measure for repulsion.
In Figs. 4 and 5 we do not give error bars for the fitted values of β.In our previous paper [10] we made an attempt to estimate the error by fitting β to approximately the same number of points from a 2D 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 β.Because in the ensembles of Goshawks and Eagle Owls (F) 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 [10].The NN values vary in range between β = 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 β = 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 more 2D 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 β-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.2. Notice the lack of nests in the plains (right) until 2010.
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.
Finally let us discuss the findings for the Eagle Owls.Because of the lowest statistics, see Fig. 5 and Table 2, about a factor of 2-10 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 β 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 β = 0.In contrast, the NNN values remain consistent with β = 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 Common Buzzards, Eagle Owls (F) show a decreasing with short range interaction, despite an increase or stabilisation in population, respectively.It be noted that for relatively few spacings of the order 20, the is strongly peaked, often leading to larger values for the fitted β.For a better resolution the first NN value at β = 3 is suppressed in Fig. 5 top left.
time dimension of the strength the repulsion also reflects the species' biology as well as population trends in the study area.In Common Buzzard, we see an increase of β 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 [14].Why β 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 effectively 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 β ≈ 0 for the fit, this may not yet indicate a complete absence of repulsion.

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 2D 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 2D 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 β-value to the spacings between species A and B.
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 to 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 subsection.The respective sums of all spacings per species over the entire period can be found in Table 2, last row.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.  2, where for Eagle Owls we distinguish (F) and (P) in the top row and bottom row, respectively.For the comparison with Eagle Owls (F) we have 174 spacings (top left and middle), and for Eagle Owls (P) we have 40 spacings (bottom left and middle).In the comparison Goshawk to Buzzard (right column) we have 423 spacings available.When the best fit to the 2D Coulomb gas gives β = 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.
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: β = 0.6 for Common Buzzards, β = 0.8 for Goshawks, and β = 0.5 for Eagle Owls (F).Clearly, the repulsion measured in β 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 β = 0.15, followed by that between Goshawks and Eagle Owls (F), and Common Buzzards and Eagle Owls (P) at β = 0.1 each.Do these low values of β 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 β-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 [12,48].The effect of Eagle Owls is rather different and has been shown to be clearly negative for both Common Buzzard and Goshawk [16,14].These interactions are, however, very dynamic at a very 14], and hence simplistic 2D Coulomb gas picture employed here is perhaps at its limits in detecting 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 the dimension of roughly 12 × 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.[5], and we simply analytically continue to real D > 0 here: It is normalised with first moment equal to unity.The repulsion ∼ s D−1 originates entirely from the (fractional) area measure.
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 ≫ 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 For the a Kolmogorov-Smirnov fit to the CDF of eq.(4.1), 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 β = 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 β 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 β ≥ 0, we encounter the situation that the fit to β = 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 β 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 β.For instance, when finding β NN > β NNN we expect an interaction weaker than Coulomb, or of shorter range.Likewise, for β NN < β 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 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 [51] for a standard work on quantum field theory, which we follow here for the derivation.
The Yukawa potential can be defined in D (integer) dimensions3 by taking the inverse Fourier transformation of the propagator with mass m and coupling constant g: The integral can be performed using polar coordinates, and in D = 2 we obtain It only depends on the radial distance r = |⃗ 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 [52, 6.532].It has a logarithmic singularity at the origin, If we define a length scale by γ = 1/m, we obtain a one-parameter deformation of the 2D Coulomb interaction.Namely, for fixed distance r and large scale γ ≫ (m ≪ 1) we are back the logarithmic repulsion (shifted by a constant).At fixed γ, 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 β, we use the following shifted potential (4.3): A comparison between the two potentials at fixed β = 1 is shown in Figure 8 left.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 [10].That is, we use a Metropolis-Hastings algorithm [54,53] 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 after −V before ) if it leads to higher energy, otherwise it is rejected.(Note that by including β 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 γ < ∞, as the global spectrum is non-uniform.

Potential
Nearest Neighbour Next-to-Nearest Neighbour Distance While the potential V Yukawa (r) differs considerably for the parameter values chosen, the spacing distribution converges rather rapidly to the one of the Coulomb potential at γ → ∞, 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 for all fixed values of β for deformed that behave only locally as a logarithm [55].
The fitted value of the scale γ could in principle be translated into an actual correlation length, but the effect of lowering the parameter γ is qualitatively similar to lowering β, and only becomes noticeable at relatively small values of γ 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 β-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 β 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 β 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 β and γ.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.In this subsection we do not split the group of Eagle Owls into two.
The amount of reusage of occupied nests from one year to the next year varies between 30% (Common Buzzards), 45% (Goshawks) and 65% (Eagle Owls (F)+ (P)) within one species, on average.The reuse of any nest from any previous year is about 10-20% higher.The frequent reuse of nest sites over years and even decades has been documented previously for the study area and is typical of avian predators [48].Habitat heterogeneity in many dimensions is commonly observed in these species [44], and individual preferences for certain habitat features are also common [48], leading to frequent reuse of nests.In other species, the same nest has been in use for centuries [49] and even millennia [50].
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 picture, 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 [10] 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.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 [48].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 [48,16].
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 is 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.It is approximately constant in Fig. 9 and largely fluctuates around the mean.As we have seen in Fig. 4 (and Table 2) for the Common Buzzard, their population has grown from 63 to 267 pairs, and from a single pair to 23 for all Eagle Owls in total, see Fig. 5 and Table 2.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.

Dicussion and Open Questions
This study has shown that a simple, one-parameter dependent approach from statistical without any biological assumptions, has the potential to explain important features in this threespecies guild of avian predators.We have shown that the intraspecific repulsion measured through β clearly deviates Poisson statistics.The relative ranking in repulsion strength amongst the three species, averaged over all has been found compatible with experience from ecology.
The correlation between β 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 63 to 267 pairs has been paralleled by a clear rise in β, 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 β.In contrast for the Eagle Owls, with the population increasing from a single pair to over 23, we have seen a steep decrease in β for NN and no repulsion between NNN.This relation is not clear to us and may be due to the limitation of our approach, or an artefact due to very low statistics.The change in β 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 (F) 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 β = 0.15 − 0.1 above Poisson at β = 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 (P) 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 our approach (and perhaps present statistics).
We have shown that a simple formula based on 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 measures.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.Finally, as a cross check we have made a numerical experiment for two examples of year with medium to low statistics.Introducing a random displacement of the original locations of the nests by an amount comparable to the mean spacings between the nests, the displacement destroys the repulsion we have seen in fitting β to the NN spacing, and rapidly leads to a Poisson distribution.
It would be very valuable to find an underlying microscopic model that leads at least approximately to the kind of spacing distribution of a 2D Coulomb gas, as the one we have applied in our comparison with the data.We feel that this example of transfer between disciplines has been an insightful exercise.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.Each moment is normalised.The angular and radial integrals decouple in all cases.The moments for complex Ginibre are found in [34] and reproduced here in Table 1 1.Numerically determined moments in the complex Ginibre ensemble from [34], rounded to 3 digits (middle column), and the corresponding values for Poisson from (A.5) [5].
For the 2D Coulomb gas with intermediate values 0 < β < 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.
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 β-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 Common Buzzards Goshawks Figure 10.moments from complex spacing ratios in (A.5) for Poisson and from Table 1 for complex Ginibre [34] 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.
Poisson and Ginibre and then moving clearly below the Ginibre value.In comparison, in Section 3.1 for the fitted β-values for NN and NNN spacing distribution for Goshawks, see Figure 4 right, a clear decent from β-values from about 1.5 down for 0.5 for both NN and NNN was detectable.
In the comparison of β-values within one species and amongst different species, it was important to have even approximate values for β, 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.Table 2.The number of nests per each species per year, including the split of the Eagle owls into (F) and (P) in the last two columns.The sum of all nests per species over all years is given in the last row.

Figure 1 .
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.The extent of the monitored area is roughly 12 × 25 km.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.

Figure 3 .
Figure 3.Comparison of the fitted β-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 β (dashed line).We present examples for fits for a single year 2020 (left column), an average over 10 consecutive years from 2011-2020 (middle column), and over all years (right column).For the Common Buzzards (top row) we have 228 spacings in 2020, 2187 spacings in 2011-2020, and 3377 spacings in all years.The corresponding numbers of spacings for Goshawks (middle row) are 19, 226 and 423 respectively, and for Eagle Owls (F, bottom row) they read 16, 116 and 174, cf.Table2in Appendix D for details.Shown are histograms and spacing distribution, but we emphasise that the fits were obtained using the cumulative distributions.Clearly when moving from individual years to time moving averages, the number of spacings for Goshawks and Eagle Owls (F) leads to an improved statistics and quality of the fits.

Figure 4 .Figure 5 .
Figure 4. Top row: Fit of β 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).The number of spacings per fit ranges between 1050 and 2420 for the Buzzards, and between 170 and 290 for the Goshawks, as can be seen from the plots in the 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.

6 .Figure 6 .
Figure 6.2D Coulomb gas and surmise fits for the β-values of the NN spacing distributions (full respectively dashed line) between all pairs of species for all years, cf.Table2, where for Eagle Owls we distinguish (F) and (P) in the top row and bottom row, respectively.For the comparison with Eagle Owls (F) we have 174 spacings (top left and middle), and for Eagle Owls (P) we have 40 spacings (bottom left and middle).In the comparison Goshawk to Buzzard (right column) we have 423 spacings available.When the best fit to the 2D Coulomb gas gives β = 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.

Figure 7 .
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).

Figure 8 .
Figure 8. Left: Comparison of the 2D Coulomb potential (γ = ∞, dashed) and the the 2D Yukawa potential (4.6) varying from γ = 0.01 (bottom red full line) to γ = 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 γ.The maximum increases with γ and moves from left to right towards the rightmost curve, given by the 2D Coulomb potential.

1 0
This leads to the following expression for the first moments in the Poisson case: drr r cos(θ)ρ Poi (r, θ) = 0 .(A.5) If a new configuration is generated from an old 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, as observed for the Eagle Owls.

Table
below for completeness.