Introduction
Why have we detected no sign of intelligent life beyond the Earth? This fundamental question continues to challenge our deepest-held beliefs about humanity and our place in the Universe. Fermi's Paradox forces us to confront our Copernican assumptions about our lack of uniqueness with the lack of extraterrestrial intelligences (ETIs, see e.g. Brin Reference Brin1983; Ćirković Reference Ćirković2009). Its strongest formulation can be given as follows (Tipler Reference Tipler1980).
Imagine a civilization constructs an interstellar probe that is self-replicating. Such a probe would be able to produce a copy every time it visits a new star system. As each copy makes copies, the number of self-replicating probes (SRPs) grows exponentially, and every star in the Milky Way is explored on a timescale much, much shorter than its age. Estimates for this exploration timescale vary but are as short as ten million years (Nicholson and Forgan Reference Nicholson and Forgan2013), and perhaps shorter still.
Given that this timescale is much shorter than the age of the Earth, and only one ETI constructing SRPs is sufficient to produce this scenario, on balance we should expect to see an interstellar probe orbiting the Sun. And yet, we do not. How can this be resolved?
Among many possibilities, we can include solutions that require civilizations to be rare. However, as a single civilization is sufficient to swamp the galaxy in SRPs, we are effectively asking for humanity to be alone in the Universe.
It may well be the case that other intelligent beings exist and that their probes are en route, and may not arrive for several thousand or several million years. This demands that the biological timescale of most of the Milky Way is somehow correlated. Perhaps this is due to global regulation mechanisms, large-scale destructive events that reset the biological clocks of many civilizations simultaneously (Annis Reference Annis1999; Vukotic and Ćirković Reference Vukotic and Ćirković2007; Vukotic and Cirkovic Reference Vukotic and Cirkovic2008). However, there are no known astrophysical regulation mechanisms that are truly global. For example, when the Milky Way's central supermassive black hole enters an active accreting phase, the subsequent radiation output results in a regulation mechanism with a correlation length of order ten thousand light years (Balbi and Tombesi Reference Balbi and Tombesi2017). The Milky Way is much larger than this correlation length, and even a small handful of uncorrelated biospheres is fit to make this type of solution untenable.
Of course, it may well be the case that an interstellar probe is in the outer solar system, and we have not yet found it (Papagiannis Reference Papagiannis1978; Freitas Reference Freitas1983). While this is indeed plausible, the possible places for an interstellar probe to hide continues to decrease (see also Haqq-Misra and Kopparapu Reference Haqq-Misra and Kopparapu2012).
Other solutions suggest that SRPs themselves are rare. There have been several arguments put forth regarding the safety of SRPs as a technology. Sagan and Newman (Reference Sagan and Newman1983) suggest that civilizations will voluntarily refrain from building SRPs, for a host of reasons. For example, self-replication could result in encoding errors. These ‘mutations’ will propagate from generation to generation with unforeseen, unintended consequences. The consequences of this could be severe, e.g. a genocidal conversion of a species and its technology into probes. It has been argued that these risks would persuade intelligent beings to place a moratorium on SRPs. Such arguments are notorious for their anthropological assumptions, and the larger assumption that ‘civilisations’ are unified in purpose and execution – an assumption humanity repeatedly invalidates (for a limited set of examples, see e.g. Collins Reference Collins2008; Denning Reference Denning2011; Lempert Reference Lempert2014, also Forgan Reference Forgan2017). This heterogeneity crucially undermines the likelihood of a Galactic moratorium. Once technology achieves a certain threshold, a very small number of individuals in a single civilization become capable of building illegal technology and producing a species-ending moment, or indeed an SRP fleet (see e.g. Sotos Reference Sotos2017).
We will explore a variant of the unintended consequences of SRPs, the ‘Predator-Prey’ hypothesis (Chyba and Hand Reference Chyba and Hand2005). In this scenario, a subset of SRPs mutate into predators of other SRPs. For an SRP to make a copy of itself, it is likely that cannibalizing another SRP will be the most energy-efficient solution. In the standard description of this scenario, predators reduce the prey population until the available prey are exhausted, reducing the visibility of SRPs in the Milky Way.
However, we should also expect in this scenario that the spread of SRPs across interstellar space will be modulated by the non-trivial population dynamics of predator-prey systems, which are among the most well-studied fundaments of mathematical biology. The full implications of these population dynamics in SRPs are largely unstudied (although see Wiley Reference Wiley2011).
In this work, we apply the classic Lotka-Volterra formalism for predator-prey systems to an interstellar network. In the network, each star plays host to a predator-prey population, which transmits and receives both predator and prey from neighbouring star systems. We will consider under what circumstances predators can reduce the population significantly, as well as what circumstances permit stable, significant populations of prey to remain in existence.
Method
We model the Galactic Habitable Zone (GHZ) as a graph of $N_\ast$
stars. The stars are distributed in space, with their locations fixed. The semimajor axes of the stars a i around the Galactic Centre are exponentially distributed to simulate the Milky Way's surface density profile:
with the scale radius r S = 3.5 kpc (Ostlie and Carroll Reference Ostlie and Carroll1996). The minimum and maximum permitted radius of the stars is [7, 10] kpc respectively, following the Gowanlock et al. (Reference Gowanlock, Patton and McConnell2011) model of the GHZ. We assume a uniform eccentricity distribution, provided that the star's orbit prevents its closest approach to the Galactic Centre being smaller than the inner radius of the GHZ. We also restrict the inclination of the orbits so that they do not exceed 0.5 radians. The longitude of the ascending node, the argument of periapsis and the true anomaly are uniformly sampled in the range [0, 2π] radians. Note that once the stellar orbital parameters are determined, we fix the stellar positions throughout the calculation.
Each star provides a vertex to our graph, and we construct the graph G such that for a star i, any star j within a minimum distance R min of i is connected by an edge. We then define the minimum spanning tree T of this graph, and use this for computing predator/prey evolution (see Fig. 1).
The stellar network on which we conduct our calculations. Each star represents an individual system on which we integrate the Lotka-Volterra equations. Lines indicate edge connections between stars, which determines the rate of outflow/inflow of predators and prey onto each star (see text). This stellar network is the minimum spanning tree T of the graph G (see text for definition).
For each star i in the network, we solve the Lotka-Volterra equations for the local numbers of prey R i, and number of predators P i, with inflow and outflow rates of prey determined by the edges that connect each star to its neighbours. For logistic growth, the equations for star i are:
To model exponential growth, we can simply set the carrying capacity K to very large values (1030). We give a list of all variables with their definition in Table 1. The outflow of prey from star i to star j, O R,ij is calculated assuming a fixed outflow rate o R,i, a fixed probe speed v probe, and the distance between i and j, D ij:
A list of variables used in this paper, with their definition
As we begin our simulations with many stars containing zero prey initially, we demand that O R,ij be zero until a sufficient time has elapsed for the first prey to arrive at i to complete the journey to j. The total outflow is then
Inflows are computed similarly:
where all outflow/inflow terms are pairwise, i.e.
and an identical set of equations are used for predator outflow and inflow (O P,i, I P,i).
Results
Tests
Single star, exponential growth
We test the code by considering a single star with no inflow or outflow, to ensure that we retrieve the solution to the classic or ‘vanilla’ Lotka-Volterra equations for exponential growth (i.e. infinite carrying capacity), where we now drop the i subscript for clarity:
Table 2 shows the parameter values for this run, and Figure 2 shows the resulting behaviour. We see that the predator/prey populations are oscillatory, both with period 8.6 Myr, and out of phase. This represents a fixed locus in predator-prey space (right panel of Fig. 2). This locus can be determined by condensing the coupled equations into a single equation:
If we define the above as a Hamiltonian:
The Lotka-Volterra solution assuming exponential growth in a single star system with zero inflow/outflow.
Parameter values used for the results shown in Figure 2
We can use standard Hamiltonian analysis (see Murray Reference Murray2004) to obtain a single solution for the predator/prey population:
This defines an initial point on the R−P locus, with all other points on the locus defining a surface of constant H.
Single star, logistic growth
We repeat the calculation of the section ‘Single star, exponential growth’, where we now impose a prey carrying capacity K R = 20 (see Table 3). We recover the standard result, that the stabilizing effect of carrying capacity damps the oscillations in predator/prey populations, until an equilibrium is eventually found at late times (Fig. 3). This is represented in R−P space by a spiral with end point given by the equilibrium solution. Note that the equilibrium prey value is much less than the carrying capacity, which is a common outcome in models of this type. Depending on the input parameters, steady-state solutions are possible where either the prey or predators dominate the combined population.
The Lotka-Volterra solution assuming logistic growth in a single star system with zero inflow/outflow. The prey carrying capacity is 20.
Parameter values used for the results shown in Figure 3
Logistic growth, globally constant parameters
Moderate, constant outflow rate
We now consider a full stellar network, where each star possesses the same fixed values for all parameters (see Table 4). We place an initial population (R 1, P 1) on star 1, with all other stars hosting zero prey/predators initially. The global picture is similar to the single star case (section ‘Single star, logistic growth’) - the initial oscillatory phase is quickly damped towards equilibrium values. Again, we find that the total prey population is lower than the maximum permitted by carrying capacity ($N_\ast K_{R,i}=10^4$
), but the predator population is also constrained, and hence we find the prey population can be sustained at relatively high levels (provided they do not exhaust local resources for self-replication).
Parameter values used for the results shown in the section ‘Moderate, constant outflow rate’
Interrogating individual stars reveals that all systems assume damped oscillatory states similar to the previous section, all with equal oscillation periods for both the prey and predator populations (now 9.17 Myr). This is made possible by the relatively large outflow rate (o R,i = 10−3), which seeds a system with sensible initial quantities of prey/predators while remaining sufficiently small that the internal predator-prey dynamics dominate the population's evolution.
As the initial prey population requires time to traverse the stellar network, each predator-prey system begins operating at a slightly different initial time. The oscillations seen in the global population have a slightly increased wavelength compared to individual systems, due to the constructive interference of many oscillatory curves, each slightly out of phase with each other, ‘smearing’ the curve over a longer period range.
Low, constant outflow rate
We consider the effect of outflow rate by repeating the previous experiment, with the outflow parameter reduced to very low values (o R,i = 10−9). We find that this significantly alters both the global behaviour and the behaviour of individual systems (Figs 5 and 6).
The predator/prey population in the Galactic Habitable Zone, assuming logistic growth over 500 stars with globally fixed predator/prey growth and death rates (see Table 4).
As Figure 4, but with significantly reduced outflow.
The predator/prey populations of selected stars, where the outflow rate is constant and very low (o R,i = 10−9).
The reduced outflow rate ensures that a greater time interval is required for all stars to host predator-prey systems with sufficient quantities of each population. As a result, the initial evolution of all systems depends strongly on their local environment – the number of directly connected neighbours, and the quantities of prey/predators arriving from each neighbour. Once each star in the entire Galaxy is sufficiently populated, the individual systems are able to attain an equilibrium state, resulting in a global equilibrium with total predator/prey counts very similar to the previous example (where o R,i = 10−3). Note that while the behaviour of each individual system is markedly different (Fig. 6), the periodic behaviour of prey and predators in any given system remains tightly coupled.
This run demonstrates the highly time-dependent nature of this solution to Fermi's Paradox. Two quite distinct phases of evolution can be characterized – a relaxation phase which persists for the first 200–400 Myr, followed by an equilibrium phase that endures beyond 400 Myr. Given that human SETI searches span a time interval of <10−4 Myr, any constraint we can place on the total number of probes in a given star system will only be a brief snapshot. As a result, our ability to use observations of any kind to constrain any of the parameters of our model will be extremely limited indeed, even with a bona fide detection of alien SRPs.
Spatially varying prey growth rates
It is quite likely that the prey growth rate around individual stars will vary, perhaps due to the quality of resources available for self-replication. This is likely to be a function of system chemical composition (i.e. stellar metallicity), but also the degree of element differentiation and chemical processing experienced by asteroids and minor bodies orbiting said star. The location of the debris may also limit its usefulness to SRPs – if most of the ‘desirable’ raw material resides inside a deep gravitational potential well, this places energy constraints on the probe's manoeuvrability, further restricting its maximum replication rate.
We consider this possibility by rerunning the previous calculation (section ‘Moderate, constant outflow rate’), but now randomly sampling b R,i in the range [0.5, 2.0]. The resulting total population can be seen in Figure 7. The steady equilibrium of the previous section has disappeared. Periodogram analysis shows oscillatory behaviour over a range of periods, with principal period 8.77 Myr. However, there is significant amplitude spread around this principal period, as the oscillation period of an individual system (equivalently, the velocity of the system's trajectory around its constant H surface in R−P surface) will depend on its individual Lotka-Volterra parameters (b R,i, d R,i, b P,i, d P,i). The combination of a range of oscillation periods results in a ‘smearing’ of the total oscillatory period.
As Figure 4, where we now allow the prey growth rate to vary uniformly amongst the 500 stars in the range b R,i = [0.5, 2.0].
Spatially varying prey carrying capacity
In a similar vein to the previous section, the resources provided by local asteroid belts may support varying levels of prey populations. We might therefore expect that the prey carrying capacity will vary between individual stars. We investigate this possibility by resetting b R,i = 0.6 and randomly sampling K R,i from a uniform distribution in the range [5, 30]. The resulting total populations again assume an oscillatory state, but with no clear principal period, and reduced variations in the value of R and P (Fig. 8). The total prey population tends to remain at larger values than the previous case, mostly because the prey growth rates can now be set at a relatively large value.
As Figure 4, where we now allow the prey carrying capacity to vary uniformly amongst the 500 stars in the range K R,i = [5, 30].
When we consider individual star systems (Fig. 9), we can see that the individual predator populations are difficult to distinguish from the total predator population. However, we can see that the prey populations show a variety of periodicities, defined not only by the local K R,i but also the carrying capacity of its neighbouring stars. For example, the prey population in the top right panel of Figure 8 exhibits a visible periodicity of around 100 Myr, whereas the bottom right panel shows no obvious evidence of periodicity.
The predator/prey populations of selected stars, where the prey carrying capacity varies uniformly in the range K R,i = [5, 30].
Notably, in contrast to every other simulation conducted so far, computing periodograms for individual systems reveals that predator and prey populations no longer share the same overall periodicity. This decoupling is a direct consequence of the predator populations practising exponential growth, while the prey populations exercise logistic growth with varying K.
Spatially varying outflow rates
Finally, we consider the effect of varying outflow rates between star systems by running the model with o R,i uniformly sampled in the range (10−4, 10−3) (Figs 10 and 11). As with the previous section, this variation forces the predator and prey populations of individual stars to oscillate on differing periods, with the oscillation period being sensitive to the local o and the value of o for connected neighbours. However, Figure 10 shows that this decoupling at local scales is not evident at global scales. Periodogram analysis confirms that the predator and prey populations on Galaxy-wide scales oscillate with identical periodicity!
As Figure 4, where we now allow the outflow parameter o R,i to vary in the range [10−4, 10−3].
The predator/prey populations of selected stars, where the outflow parameter o R,i varies uniformly in the range [10−4, 10−3]..
Discussion
Comparison with other studies of predator-prey dynamics
Predator-prey dynamics has a rich history of study in mathematical biology. The Lotka-Volterra equations have been applied to a variety of cases, in particular considering how the system is generalized to more than two species, as an attempt to model the food chain of an ecosystem (for examples see Macarthur and Levins Reference Macarthur and Levins1967; Smale Reference Smale1976; Palamara et al. Reference Palamara, Zlatic, Scala and Caldarelli2011; Gavina et al. Reference Gavina, Tahara, Tainaka, Ito, Morita, Ichinose, Okabe, Togashi, Nagatani and Yoshimura2018).
At a basic level, one can incorporate spatial effects into the vanilla Lotka-Volterra system by recasting it as a set of reaction-diffusion equations (Cross and Hohenberg Reference Cross and Hohenberg1993), where one can think of predators and prey as two reactants combining to form a product, and both entities diffuse spatially with a diffusion constant D, e.g.:
Such systems can be stable in the absence of diffusion (D=0), only to become unstable when diffusion is added (Turing Reference Turing1952). Of course, these models are inherently symmetric, and cannot account for spatial heterogeneities.
Most modern attempts to incorporate geography into calculations of this sort assume a lattice configuration, upon which either analytic or probabilistic solutions can be obtained. For example, Frachebourg et al. (Reference Frachebourg, Krapivsky and Ben-Naim1996) considered the behaviour of predator-prey systems on a 1D lattice. For a two-species system, one can model the entire evolution in terms of interfaces separating species. Over time, species tend to organize into a mosaic of alternating domains, with the size of each domain increasing linearly with time. If the number of species exceeds 5, these domains can become ‘frozen-in’ (although this depends on whether the species chain is symmetric, i.e. can species 1 eat species 2 and species 5?).
Tomé and de Carvalho (Reference Tomé and de Carvalho2007) construct a probabilistic cellular automaton, inspired by the Lotka-Volterra equations, on a 2D lattice. They are able to show that self-sustained stable oscillations can be set up in the system, just as in the ‘vanilla’ Lotka-Volterra case. However, these oscillations are stable against changes in the initial conditions.Footnote 1 A more generalized version of this result was obtained by Rozhnova and Nunes (Reference Rozhnova and Nunes2010), who considered random networks with k neighbours per node (see also Ohtsuki and Nowak Reference Ohtsuki and Nowak2006).
As far as the author is aware, there are no examples of a coupled Lotka-Volterra system computed on spatial graphs/networks (although see Palamara et al. Reference Palamara, Zlatic, Scala and Caldarelli2011 for an example on a network of species, to resemble food webs). Lotka-Volterra systems usually model spatially continuous environments, although deliberately inserting heterogeneity into a continuous environment has been shown to remove the sustained oscillations, as we have found in our analysis (McLaughlin and Roughgarden Reference McLaughlin and Roughgarden1991; Täuber Reference Täuber2011).
Limitations of the analysis
The flow of predators/prey between star systems depends heavily on the distance between them. In our model, we have kept stellar positions fixed. If we allowed the stars to move, we can expect that this will result in quasi-periodic forcing of the flow rates. The periodicity of both predator and prey populations for a given star would be further modified, to accommodate both the star's motion and that of its neighbours.
We also assume that each star system has an unending supply of resources. While the growth of prey is generally limited by the local carrying capacity K, we have not considered the possibility that replication eventually exhausts the local supply of raw materials. One could model this rather simply as a non-constant K that decays with time. If resources could be exhausted sufficiently quickly, that might limit the number of probes overall, and provide a resolution to Fermi's Paradox. However, if one considers the number of probes that can be produced from the available debris mass in the solar system, the exhaustion timescale of a typical star system is likely to be much too long to be of concern to this analysis.
That being said, we might note that the quality of raw material can vary from system to system, as a function of local metallicity. We have attempted to model this by allowing K to vary between stars. However, our models allowed K to be effectively random. In practice, K should vary according to the metallicity gradient of the Galaxy (e.g. Bergemann et al. Reference Bergemann, Ruchti, Serenelli, Feltzing, Alves-Brito, Asplund, Bensby, Gruiters, Heiter, Hourihane, Korn, Lind, Marino, Jofre, Nordlander, Ryde, Worley, Gilmore, Randich, Ferguson, Jeffries, Micela, Negueruela, Prusti, Rix, Vallenari, Alfaro, Allende Prieto, Bragaglia, Koposov, Lanzafame, Pancino, Recio-Blanco, Smiljanic, Walton, Costado, Franciosini, Hill, Lardo, de Laverny, Magrini, Maiorca, Masseron, Morbidelli, Sacco, Kordopatis and Tautvaišienė2014). This uniform variation may result in similar spatial variations in system periodicities.
It is also worth noting that predator probes can also scavenge other predators for resources to self-replicate. Adding such ‘omnivorous’ probes to a star system could have important consequences – reducing the predator population in this way could allow prey populations to grow to larger values.
In any case, the above additions to the analysis will not affect the final result – the predator/prey solution to the SRP formulation of Fermi's Paradox does not significantly reduce the SRP population, and is therefore not a viable solution.
Conclusions
In this paper, we have considered a proposed solution to Fermi's Paradox regarding the growth and spread of self-replicating interstellar probes. It has been proposed that if some SRPs were to ‘mutate’ and begin predating other probes, this would reduce the total population of probes, ensuring that humanity would not see them.
We conduct simulations of predator-prey probe evolution using the Lotka-Volterra equations, amongst a connected network of stars in the GHZ. We find that traditional competition can result in oscillating behaviour for the predator/prey populations at a given star, as well as equilibrium solutions where both local and global populations tend to fixed values. The nature of the system behaviour depends sensitively on the birth and death rates of each species, as well as the local carrying capacity and the flow of species between star systems. In any case, we find that significant quantities of prey probes can persist throughout the Galaxy – admittedly less than the maximum permitted by carrying capacity, but still sufficiently large that this solution to Fermi's Paradox is weak at best, and in effect, not a solution at all.
In summary, the self-replicating probe formulation of Fermi's Paradox remains, in our view, one of the strongest and most testing formulations, and an important check on our assumptions regarding the number of intelligent species in the Milky Way.
Author ORCIDs
Duncan H. Forgan, 0000-0003-1175-4388
Acknowledgments
The author gratefully acknowledges support from the ECOGAL project, grant agreement 291227, funded by the European Research Council (ERC) under ERC-2011-ADG. This research has made use of NASA's Astrophysics Data System Bibliographic Services. The code used in this paper is available at https://github.com/dh4gan/lotka-volterra-probes
