Non-technical Summary
Is it environment or life that drives macroevolution? A recent analysis of the massive paleobiology database argues that the answer depends on the timescale. At short timescales, less than 40 million years, it is the environment, at longer timescales, life can effectively adapt. Both the environment and life are scaling—they fluctuate over the full range of scales from millions to hundreds of millions of years (the megaclimate regime). In this paper, we present a simple model of this scaling “crossover” phenomenon. The model has some unusual features: it is fully random and is based on fractional (rather than classical integer-ordered) differential equations.
The model is driven by temperature (a proxy for the environment) and the turnover rate (a proxy for life); it has two exponents, a cross-over time and two correlations, yet it is able to reproduce not only the statistics of the temperature, diversity, extinction, origination, and turnover rates, but it also effectively reproduces the pairwise correlations between them, and this over the whole range of timescales. If forced deterministically, it gives the response to bolide impact or other sharp forcing events.
Introduction
Several centuries of paleontological research revealed that the evolution of life on Earth is characterized by high temporal complexity characterized by periods of sluggish and predictable evolution (Jablonski Reference Jablonski1986; Casey et al. Reference Casey, Saupe and Lieberman2021) with mass extinctions characterized by selectivity that is low or different in kind than in “background intervals” (Raup Reference Raup1992a, Reference Raup1994; Payne and Finnegan Reference Payne and Finnegan2007). There are also mass evolutionary radiations that are sometimes contemporaneous with mass extinctions (Cuthill et al. Reference Cuthill, Guttenberg and Budd2020). Moreover, the factors and modes of macroevolution apparently vary with time—for example, the Cambrian explosion or Ediacaran–Cambrian radiation and post-Cambrian evolution (Gould Reference Gould1990; Erwin Reference Erwin2011; Mitchell et al. Reference Mitchell, Harris, Kenchington, Vixseboxse, Roberts, Clark, Dennis, Liu and Wilby2019); environment (Jablonski et al. Reference Jablonski, Roy and Valentine2006; Miller and Foote Reference Miller and Foote2009; Kiessling et al. Reference Kiessling, Simpson and Foote2010; Boyle et al. Reference Boyle, Sheets, Wu, Goldman, Melchin, Cooper, Sadler and Mitchell2013; Spiridonov et al. Reference Spiridonov, Brazauskas and Radzevičius2015; Tomašových et al. Reference Tomašových, Jablonski, Berke, Krug and Valentine2015); and timescales (Van Dam et al. Reference Van Dam, Aziz, Sierra, Hilgen, W, Ostende, Lourens, Mein, Van Der Meulen and Pelaez-Campomanes2006; Spiridonov et al. Reference Spiridonov, Stankevič, Gečas, Šilinskas, Brazauskas, Meidla, Ainsaar, Musteikis and Radzevičius2017b; Crampton et al. Reference Crampton, Meyers, Cooper, Sadler, Foote and Harte2018; Beaufort et al. Reference Beaufort, Bolton, Sarr, Suchéras-Marx, Rosenthal, Donnadieu, Barbarin, Bova, Cornuault and Gally2022). Furthermore, macroevolution is strongly influenced by Earth's systems: geological, climatic, and paleoceanographic factors (Marshall et al. Reference Marshall, Webb, Sepkoski and Raup1982; Lieberman and Eldredge Reference Lieberman and Eldredge1996; Lieberman Reference Lieberman, Macintyre and Clegg2003; Saupe et al. Reference Saupe, Qiao, Donnadieu, Farnsworth, Kennedy-Asser, Ladant, Lunt, Pohl, Valdes and Finnegan2019; Carrillo et al. Reference Carrillo, Faurby, Silvestro, Zizka, Jaramillo, Bacon and Antonelli2020; Halliday et al. Reference Halliday, Holroyd, Gheerbrant, Prasad, Scanferla, Beck, Krause, Goswami, Prasad and Patnaik2020), but also by biotic interactions, which can translate into patterns that are apparent on extremely long timescales of tens to hundreds of millions of years (Vermeij Reference Vermeij1977, Reference Vermeij2019; Jablonski Reference Jablonski2008; Erwin Reference Erwin2012). There are also questions on the role of general stochasticity and path dependence/memory in evolutionary dynamics (Schopf Reference Schopf1979; Hoffman Reference Hoffman, Nitecki and Hoffman1987; Gould Reference Gould, Briggs and Crowther2001, Reference Gould2002; Cornette and Lieberman Reference Cornette and Lieberman2004; Erwin Reference Erwin2011, Reference Erwin2016). The question is: can we reconcile in a single simple model this multitude of hierarchically organized and causally heterogenous processes producing macroevolutionary dynamics, while maintaining simplicity and conceptual clarity? Here, we argue that we can.
The development of large, high temporal resolution databases—both of past climate indicators (Veizer et al. Reference Veizer, Ala, Azmy, Bruckschen, Buhl, Bruhn, Carden, Diener, Ebneth and Godderis1999; O'Brian et al 2017, Song et al. Reference Song, Wignall, Song, Dai and Chu2019; Grossman and Joachimski Reference Grossman and Joachimski2022) and of paleobiological information, such as the Paleobiology Database (PBDB; Alroy et al. Reference Alroy, Marshall, Bambach, Bezusko, Foote, Fürsich and Hansen2001, Reference Alroy, Aberhan, Bottjer, Foote, Fürsich, Harries, Hendy, Holland, Ivany and Kiessling2008) or NOW (Jernvall and Fortelius Reference Jernvall and Fortelius2002; Žliobaitė et al. Reference Žliobaitė, Fortelius and Stenseth2017; Žliobaitė Reference Žliobaitė2022)—is transforming our understanding of macroevolution. Time series are frequently long enough that they can be studied systematically, not just as chronologies to be compared with other chronologies, but as functions of temporal scale, that is, the behavior of their fluctuations as functions of duration (or equivalently, their behavior as functions of frequency).
Before attempting to understand processes at specific timescales, it is important to understand their context, that is, the dynamical regime in which they operate. Dynamical regimes are objectively defined by scaling; they are regimes over which fluctuations are scaling (see the review by Lovejoy [Reference Lovejoy2023]). By definition, a scaling regime is one in which fluctuations ΔT (in some quantity such as temperature) are of the form ΔT(Δt) ∝ Δt H, where Δt is duration, or “lag,” scale, and H is an exponent. If such power law relationships hold (typically for various statistics), then long- and short-duration events/fluctuations only differ quantitatively; they are qualitatively the same. This is because over such a regime, long-duration fluctuations, ΔT(λΔt), at scale, λ Δt ( λ > 1), are related to the shorter-duration fluctuations, ΔT(Δt), by: ΔT(λΔt) = λHΔT(Δt), that is, the fluctuations at different timescales differ only in their amplitudes, λH (with the equality understood in an appropriate stochastic sense). In addition, we can already distinguish the qualitatively different types of regime by the sign of the exponent H. H > 0 implies that fluctuations increase with scale and can appear nonstationary, whereas H < 0 implies that they decrease with scale, they appear to converge.
An important consequence for our understanding of deep-time biogeodynamics—here understood as the joint Earth–life systems—is the robustness of the “megaclimate” regime. The megaclimate regime was first discovered in benthic paleotemperatures (Lovejoy Reference Lovejoy2015), and it was characterized by “positive scaling” (a shorthand for H > 0) on the basis of long paleotemperature data from ocean core stacks (Veizer et al. Reference Veizer, Godderis and Francois2000; Zachos et al. Reference Zachos, Pagani, Sloan, Thomas and Billups2001), see also Lovejoy (Reference Lovejoy2013) for the shorter timescale weather, macroweather, and climate regimes (up to Milankovitch scales). This implies that the difference between temperatures typically becomes larger and larger at epochs separated by longer and longer intervals of time. Theoretically, megaclimate is the hypothesis that there is a unique (presumably highly nonlinear) biogeological dynamical regime that operates over timescales spanning the range ≈1 Myr to (at least) several hundred millions of years. This would be the consequence of a unique (albeit complex, nonlinear) underlying dynamic that is valid over this wide range of scales; presumably it involves a scaling (hence hierarchical) mechanism that operates from long to short durations. A consequence is the existence of a statistical scaling regimes (notably of paleotemperatures), empirically verified throughout the Phanerozoic. While its inner scale appears to be fairly robust at around 1 Myr, its outer scale (the longest duration over which it is valid) is not known, although it appears to be at least 300 Myr. The megaclimate regime implies that the underlying biology–climate dynamics are essentially the same over these timescales: that is, at long enough time scales the statistics are stationary (although up to the outer, longest, limiting scale of the regime, they may appear to be nonstationary).
The hypothesis that biology and the climate are linked and that climate is a crucial and defining variable in ecological and evolutionary turnovers (Vrba Reference Vrba1985, Reference Vrba1993; Eldredge Reference Eldredge, Crutchfield and Schuster2003; Lieberman et al. Reference Lieberman, Miller and Eldredge2007; Hannisdal and Peters Reference Hannisdal and Peters2011; Mayhew et al. Reference Mayhew, Bell, Benton and Mcgowan2012; Crampton et al. Reference Crampton, Cooper, Sadler and Foote2016; Spiridonov et al. Reference Spiridonov, Brazauskas and Radzevičius2016, Reference Spiridonov, Kaminskas, Brazauskas and Radzevičius2017a, Reference Spiridonov, Samsonė, Brazauskas, Stankevič, Meidla, Ainsaar and Radzevičius2020a,Reference Spiridonov, Stankevič, Gečas, Brazauskas, Kaminskas, Musteikis and Kaveckasb; Mathes et al. Reference Mathes, Van Dijk, Kiessling and Steinbauer2021) is hardly controversial. However, the scope and utility of the megaclimate notion would increase if it could be backed up by direct analysis of paleobiological series, particularly of extinction and origination rates. This has now been done. A recent paper (Spiridonov and Lovejoy Reference Spiridonov and Lovejoy2022), hereafter SL, found that genus-level extinction and origination rates exhibited scaling statistics over roughly the same range as the paleotemperatures confirming that the megaclimate includes these key macroevolutionary parameters (see Fig. 9, the left-hand side of the figure in the section “The Statistics of the Simulated Series Resampled at the Data Sampling Times,” for a plot of the data used in SL and modeled in this paper).
The shortest scale of SL's paleobiological time series was closer to ≈ 3 Myr (the average age resolution was 5.9 Myr), which corresponds to the durations of the shortest PBDB stages—a standard shortest time resolution for Phanerozoic-scale global biodiversity analyses (e.g., Alroy et al. Reference Alroy, Aberhan, Bottjer, Foote, Fürsich, Harries, Hendy, Holland, Ivany and Kiessling2008; Alroy Reference Alroy2010b), although note that some sub-age data are available at resolutions closer to 1 Myr). Systematic reviews and multiple case studies revealed that even variously defined (molecular, morphological, phylogenetic, and taxic) evolutionary rates universally exhibit negative temporal scaling (H < 0) behavior (Gingerich Reference Gingerich1993, Reference Gingerich2001, Reference Gingerich2009; Roopnarine Reference Roopnarine2003; Harmon et al. Reference Harmon, Pennell, Henao-Diaz, Rolland, Sipley and Uyeda2021; Spiridonov and Lovejoy Reference Spiridonov and Lovejoy2022), which suggests the universality of the temporal scaling—hence hierarchical—evolutionary dynamics. An inner megaclimate scale of ≈1 Myr was also proposed in Lovejoy (Reference Lovejoy2015) and is discussed in the nonspecialist book Weather, Macroweather and Climate (Lovejoy Reference Lovejoy2019). The scaling, and thus by implication dominance of timescale symmetric hierarchical interactions, was also detected on multimillion-year timescales in sedimentation rates/stratigraphic architecture (Sadler Reference Sadler1981), sea level (Spiridonov and Lovejoy Reference Spiridonov and Lovejoy2022), and dynamics of continental fragmentation (Spiridonov et al. Reference Spiridonov, Balakauskas and Lovejoy2022), which shows universality of the pattern in major Earth systems as well. Therefore, the time-scaling patterns of evolution and megaclimate overlap at the very wide range of temporal scales (from ≈106 to > 4 × 108 yr), which motivates the development of quantitative models that explicitly tackle and integrate these timescale symmetries.
If macroevolution and climate respect wide range scaling, then it may be possible to resolve a long-standing debate in macroevolution. In terms first posed by Van Valen (Reference Van Valen1973), we may ask: are evolutionary processes dominated by external factors—especially climate, the “Court Jester” (Barnosky Reference Barnosky2001; Benton Reference Benton2009)—or by life itself—the “Red Queen” (Van Valen Reference Van Valen1973; Finnegan et al. Reference Finnegan, Payne and Wang2008)? SL proposed a scaling resolution of the debate in which at scales below a critical transition time τ of ≈40 Myr, the climate process is dominant, but there is a crossover beyond which life (self-regulating by means of geodispersal and competition) is dominant. SL thus quantitatively concluded that at long enough timescales, the Red Queen ultimately overcomes the Court Jester. The scaling processes of the Earth system here are playing a double role (thus the Geo-Red Queen theory)—climate fluctuations growing with timescale cause perturbations in diversity to grow in size, but at the same time, at longer and longer timescales, fluctuating climates and plate tectonics cause the mixing and competitive matching of biota, thus effecting global synchronization that results in a crossover when an unstable and wandering diversity regime changes to a longer timescale fluctuation canceling or stabilizing regime (Spiridonov and Lovejoy Reference Spiridonov and Lovejoy2022).
Physicists use the term “crossover” as shorthand to describe analogous phenomena involving processes that are subdominant over one scale range but eventually become dominant at longer scales. However, such transitions are typically modeled by Markov processes such that the autocorrelations are exponential, implying that at the critical timescale, the transition between two regimes is fairly sharp. In SL, on the contrary, in keeping with the basic megaclimate scaling dynamics, the crossover was postulated to be a consequence of the interaction of two scaling processes, that is, the transition is a slow, power law decay of one and the slow emergence of another. An analogous scaling crossover phenomenon was found in phytoplankton, in which the competing scaling processes were phytoplankton growth (with turbulence) and a predator–prey process of zooplankton grazing (Lovejoy et al. Reference Lovejoy, Currie, Tessier, Claeredeboudt, Roff, Bourget and Schertzer2001).
SL argued that both macroevolution and climate respect wide-range statistical scaling, but that nevertheless, their quantitative and qualitative differences are significant, and these differences were the key to understanding the diversity (D) statistics that appeared to be involve a crossover between two different power laws. While temperature (T) fluctuations vary with timescale, Δt as ΔT(Δt) ≈ ΔtHT with HT ≈ 0.25, the corresponding laws for extinction (E) and origination (O) have HE ≈ HO ≈ −0.25. When H > 0, fluctuations grow with scale, such that the corresponding series tend to “wander” without any tendency to return to a well-defined value, and they appear “unstable.” On the contrary, when H < 0, successive fluctuations tend to have opposite signs, such that they increasingly cancel over longer and longer timescales, and they fluctuate around a long-term value, thus appearing “stable.”
To deepen our understanding, it is necessary to build a quantitative model of the interaction of climate and life. In recognition of the strongly nonlinear nature of evolutionary dynamics, numerous deterministic chaos models such as predator−prey models (e.g., Huisman and Weissing Reference Huisman and Weissing1999; Caraballoa et al. Reference Caraballoa, Colucci and Han2016) have been developed. Although extensions with some stochastic forcing exist, for example, in Vakulenko et al. (Reference Vakulenko, Sudakov and Mander2018), the stochasticity is a perturbing noise on an otherwise deterministic system. In paleontology, the model of exponential (unconstrained) proportional growth of diversity was historically popular (Stanley Reference Stanley1979; Benton Reference Benton1995) or expanded to include possible accelerations due to niche construction effects (a second-order positive feedback)—a hyperbolic model (Markov and Korotayev Reference Markov and Korotayev2007). These simple models of expansion were contrasted with single or coupled logistic models of resource-constrained competitive macroevolutionary dynamics, sometimes also supplemented with random perturbations that account for effects of mass extinctions (Sepkoski Reference Sepkoski1984, Reference Sepkoski, Jablonski, Erwin and Lipps1996); or implicitly hierarchical, and also competition-constrained, Gompertz models (Brayard et al. Reference Brayard, Escarguel, Bucher, Monnet, Brühwiler, Goudemand, Galfetti and Guex2009). However, such models assume that only a few degrees of freedom are important (typically fewer than 10), whereas the true number is likely to be astronomical. It is therefore logical to model the process in a stochastic framework (involving infinite dimensional probability spaces), where the primary dynamics are stochastic, using the scaling symmetry as a dynamical constraint. Therefore, there is growing recognition of stochastic models as essential tools for understanding macroevolutionary dynamics. Actually, some of the first models that tried to explain complexities of macroevolutionary dynamics were stochastic Markovian birth and death models (Gould et al. Reference Gould, Raup, Sepkoski, Schopf and Simberloff1977; Raup and Valentine Reference Raup and Valentine1983; Raup Reference Raup1985, Reference Raup1992b; Nee Reference Nee2006). Several recent applications of linear stochastic differential equations were used in causal inference of macroevolutionary drivers and competitive interactions between clades (Liow et al. Reference Liow, Reitan and Harnik2015; Reitan and Liow Reference Reitan and Liow2017; Lidgard et al. Reference Lidgard, Di Martino, Zágoršek and Liow2021).
Beyond the realism of implicitly involving larger numbers of degrees of freedom, stochastic models have the advantage that they may be linear, even though the corresponding deterministic models may be highly nonlinear. Also, by the simple expedient of using fractional-ordered differential equations rather than the classical integer-ordered ones, stochastic models can readily handle scaling, which is rarely explicitly considered in macroevolutionary analyses. Fractional differential equations provide a natural way of defining processes that vary over a wide range of timescales. Although, we forced the model with a nonintermittent Gaussian white noise in this paper, in future, this can easily be replaced by strongly intermittent (multifractal) processes that are presumably necessary to realistically model the extremely intermittent behavior, implied, for example, by mass extinctions or thermal climatic events, which are observed in macroevolution and megaclimate respectively. Scaling processes are characterized by a slow (power law) decay of the memory, much slower than an exponential rate, such that values of a time series are nontrivially correlated, even if they are separated by long time periods. Therefore, the dynamics of the system are potentially conditioned not only on the current state of the system but also on its distant past. This is exactly the property that is exploited in constructing state-of-the-art descriptive and predictive models of long-memory phenomena at shorter timescales (weeks to years), namely macroweather forecasts, based on the Scaling LInear Macroweather Model (SLIMM; Lovejoy et al. Reference Lovejoy, Del Rio Amador and Hébert2015) and StocSIPS (Del Rio Amador and Lovejoy Reference Del Rio Amador and Lovejoy2019, Reference Del Rio Amador and Lovejoy2021; Lovejoy and Del Rio Amador Reference Lovejoy and Del Rio Amador2023). The same basic model using deterministic climate forcings produces climate projections to 2100, notably with much lower uncertainty than classical general circulation model (GCM) approaches (Hébert et al. Reference Hébert, Herzschuh and Laepple2022; Lovejoy Reference Lovejoy2022b; Procyk et al. Reference Procyk, Lovejoy and Hébert2022).
The useful scaling property of fractional equations arises because they have impulse response functions (Green's functions)—and hence solutions—that are based on scaling (power laws) rather than the exponential Green's functions associated with integer-ordered differential equations. In general, fractional derivatives are simply convolutions with power laws of different orders, and convolutions with different domains of integration define different types of fractional derivatives. In the fractional equations discussed in this paper, the particularly simple Weyl fractional derivative is used; in the frequency domain, it simply corresponds to a power law filter. Finally, it could be noted that although fractional derivatives were discovered several centuries ago, there has been an explosion of interest in them in the last decade or so, and many books on the topic are now available (e.g., Oldham and Spanier Reference Oldham and Spanier1974; Miller and Ross Reference Miller and Ross1993; Hilfer Reference Hilfer2000; West et al. Reference West, Bologna and Grigolini2003; Petras Reference Petras2011; Baleanu et al. Reference Baleanu, Diethelm, Scalas and Trujillo2012; Klafter et al. Reference Klafter, Lim and Metzler2012; Owolabi and Atangana Reference Owolabi and Atangana2019).
In this paper, we therefore build a simple model for biodiversity (D) that can reproduce and explain SL's findings. The model is parsimonious: it has only two scaling drivers—the climate and life—and by construction it reproduces the observed scaling crossover at 40 Myr. Although the model has two basic exponents (climate and life) and two coupling coefficients between temperature and turnover and between turnover and diversity, a total of four basic parameters, it satisfactorily reproduces the fluctuation statistics of D, T, E, and O as well as the turnover (E+ = O + E) and difference E − = O − E over the range of ≈3 Myr to several hundred millions of years (the longest scales available). The six variables imply 15 pairwise correlations, and the correlations that are implied by the model are not single values between pairs of variables at a unique timescale, but each correlation is a function of the timescale indicating an agreement between the model and data over a wide range of scales.
In the model, the driver of the life processes is turnover. As with many other properties of groups of biological individuals, turnover can be defined at many levels. Here we use the taxic approach to modeling and analysis, such that we consider macroevolutionary turnover at the level of genera. Similarly, as in the case of organisms from different species in populations, turnover of genera in the biota shows total perturbations to a given diversity state. Namely perturbation by subtraction (extinction) and addition (origination), which represent cases of creative destruction and the destructive creation respectively, which could work as stabilizing or destabilizing factors of the global diversity depending on the circumstances (Cuthill et al. Reference Cuthill, Guttenberg and Budd2020). Beyond this, it explains the 15 pairs of scale by scale fluctuation correlations over the same observed range. The data are from the SL paper; they represent stage-level time series of Phanerozoic marine animal genera O and E (second-for-third origination and extinction proportion [Alroy Reference Alroy2015; Kocsis et al. Reference Kocsis, Reddin, Alroy and Kiessling2019] not standardized for the duration of stages), sample standardized using the shareholder quorum method (Alroy Reference Alroy2010a) genus diversity of Phanerozoic marine animals based on PBDB data (https://paleobiodb.org). The paleotemperatures (T) are also the same as in the SL paper, obtained from Song et al. (Reference Song, Wignall, Song, Dai and Chu2019). The rates used in the analyses are proportions, per-lineage rates, which is a natural way of describing processes working on a per capita basis (such as evolution).
The inclusion of the maximal amount of data and the taxonomic precision are the essential trade-offs of any taxic macroevolutionary study. With sufficiently well-preserved fossils, it is often possible to make accurate genus-level identification (relatively high accuracy and robustness of identification). On the other hand, if the taxon can be identified only to the family level, as is often a case for multi-element taxa, probably very few remains were preserved, or in the case of single skeleton taxa, the preservation is not adequate. Therefore, the decision to use higher-rank taxonomic data far removed from the species level could increase the risk of including more noise to the data. We chose the genus level as a compromise, which is the closest taxonomic level to the species level, while on the other hand, covering more occurrences than the fossil record resolved to the species level.
Because rates of macroevolution originations and extinctions can be estimated in many ways, and their properties (accuracy, and precision) can vary depending on the properties of the sampling process and the evolutionary process itself, we performed a sensitivity analysis using Sepkoski's genus-level marine animal data and Foote's per-lineage rates (see Appendix 2). Analysis shows that despite the differences in datasets, their completeness, differences in standardization (Sepkoski's data are not sample standardized), and differences in extinction and origination rates used, the basic results pertaining to the scaling are nearly identical in both cases. Therefore, in this paper, we only discuss the more complete and sample-standardized paleobiological time series derived from PBDB data. Currently, the PBDB data are the best source for multi-lineage global-scale analysis of evolutionary patterns, despite the presence of possible distortions related to the uneven spatial sampling, due to objective geological heterogeneities of the fossil record and various historical and socioeconomic factors that significantly impacted the study of ancient life (Raja et al. Reference Raja, Dunne, Matiwane, Ming Khan, Nätscher, Ghilardi and Chattopadhyay2022; Ye and Peters Reference Ye and Peters2023). Future work toward a more even representation of the global-scale data will certainly improve the accuracy of diversity estimates; this in itself would be an interesting test of the model.
As a final comment, we should note that the basic—simplest—stochastic crossover process is the fractionally integrated fractional relaxation noise (ffRn process), whose properties were only fully elucidated very recently (Lovejoy Reference Lovejoy2022a) in the context of long-term weather forecasts (Del Rio Amador and Lovejoy Reference Del Rio Amador and Lovejoy2021) and climate projections (Procyk et al. Reference Procyk, Lovejoy and Hébert2022). The new model has conceptual commonalities with the environmental “stress model” of M. Newman that attempted to replicate the scaling statistics of extinction intensities of marine biota (Newman Reference Newman1997; Newman and Palmer Reference Newman and Palmer2003). The model presented here is more sophisticated, as it ties the principal macroevolutionary variables—O and E—to a principal geophysical scaling process—the megaclimate—in producing realistic multi-timescale global dynamics of marine animal biodiversity, while keeping its conceptual simplicity in transparently using a few crucial parameters of time-scaling and correlations. The model is also implicitly hierarchical as implied by its scaling relations, and this is a desirable feature of a unified evolutionary theory (Eldredge Reference Eldredge1985, Reference Eldredge1989; Gould Reference Gould2002; Lieberman et al. Reference Lieberman, Miller and Eldredge2007).
The Model
The Equations
Wide Range Scaling
The SL picture is one where the extra-biological factors (“the climate”) are scaling and drive biodiversity from ≈1 Myr to ≈40 Myr, where the crossover occurs, followed by the domination of biotic regulation at the longer timescales, which is also enabled by global homogenization of biota at long timescales caused by plate tectonics and changes in climate zones (Geo-Red Queen dynamics). The endemism and peculiarities or contingencies in evolutionary innovation occur at all times and places. Because evolutionary innovations are inherently unpredictable and can change carrying capacities of ecosystems and communities (Erwin Reference Erwin2012), their effects at the global scale of measurement of diversity are destabilizing: they act as random shocks. Because all innovations appear in certain times and places, it takes a certain time for the dispersion of innovations and the re-equilibration of the global-scale biosphere (adaptation of other taxa) following perturbations. The principal agent of mixing of biotas is plate tectonic motion (which also separates biota at shorter timescales). Therefore, the only way biota can readjust and equilibrate at the global scale in the light of multiscale local and regional evolutionary innovations is by means of geodispersal mediated by competition and co-adaptation. Endemicity exists at all times, but the duration of endemicity of concrete faunas is limited by unstable Earth dynamics that work toward homogenization (the same faunas or floras cannot be endemic forever). Larger taxa either disperse to their maximal capacity and persist (and become effectively cosmopolitan) or they go extinct.
The Basic Diversity Equation
Based on this picture, we propose the following Fractional MacroEvolution Model (FMEM). First we describe the model, then we comment on it.
The basic diversity equation is:
$$\tau ^h\displaystyle{{d^h( {D-s_TT} ) } \over {dt^h}} + D = s_EE_ + ; \;\,\quad E_ + = O + E$$where h is a (fractional) order of differentiation whose physical interpretation is given shortly, τ is the crossover timescale (≈40 Myr), and E + = E + O is the turnover anomaly, that is, the deviation of the turnover rate from its long-term average (in the model, E + can therefore be either positive or negative). For reference, it is defined on the right, where O and E are the anomalies of the origination and extinction rates with respect to their long-term averages (the diversity D is a similarly defined anomaly with respect to its long-term average). Whereas D and E + are already nondimensional, the temperature anomaly T must be nondimensionalized, for example, by the standard deviation of its increments at some convenient reference scale, say 1 Myr. sT and sE are constants that are determined by the coupling between T and D (sT) and E + and D (sE; see eq. 17).
The Drivers
The basic drivers are the climate (T) and life (E +), themselves driven by Gaussian white noises γT, γE:
$$\matrix{ {\tau ^{\alpha + h}\displaystyle{{d^{\alpha + h}T} \over {dt^{\alpha + h}}} = \gamma _T} \cr {\tau ^\alpha \displaystyle{{d^\alpha E_ + } \over {dt^\alpha }} = \gamma _E} \cr } $$where α is the basic biology (extinction and origination rates) exponent (α ≈ 0.25 as deduced from SL's analysis), and h is the exponent difference (contrast) between the temperature and biology, from SL's analysis: h = 0.75 − α ≈ 0.5. As discussed later (eq. 7 or eq. 35), combined with the diversity equation (eq. 1), these determine D. Equations (1) and (2) specify the model dynamics; see Figure 1 for an overall schematic.
Figure 1. A schematic showing the way the various parts of the Fractional MacroEvolution Model (FMEM) fit together. The basic drivers are shown at the top; physical drivers are the temperature (T) and turnover rate (E +). These are shown at the right, because they have nontrivial properties, such that they are best modeled as the responses to more elementary causes—the temperature and turnover rate forcings (fT, fE+). In the paper, we primarily discuss the simple case that reproduces the Paleobiology Database (PBDB) statistics, where these are Gaussian white noises (fT = γT, fE+= γE+), However, deterministic forcings such as bolide impacts are also discussed, shown here with both T and E+ forced with a Dirac function of amplitude f0,T, f0,E+, respectively. More general forcings can be used and their responses can be obtained using the impulse response functions. The middle line shows how the T, E + drivers determine the diversity (D). Finally, to complete (close) the model, we need a diagnostic equation that enables us to determine E −, E, O; this is shown in the bottom line.
The derivatives are fractional; in this paper, we use the semi-infinite “Weyl” fractional derivatives. For the arbitrary function W(t), the ζ-ordered Weyl fractional derivative is defined as:
$$\matrix{ {\displaystyle{{d^\zeta W} \over {dt^\zeta }} = \displaystyle{1 \over {\Gamma ( {1-\zeta } ) }}\displaystyle{d \over {dt}}\int\limits_{-\infty }^t {{( {t-s} ) }^{-\zeta }W( s ) ds} ; \;} & {0 < \zeta < 1} \cr } $$where Γ is the usual gamma function, and s is an unimportant variable of integration. In this paper, the range of differentiation, 0 < ζ < 1, is adequate, but more generally, if ζ is outside the range, equation (3) can simply be combined with ordinary integer-ordered differentiation to obtain the required result. Because fractional derivatives (and their inverse, fractional integrals) are, as in equation (3)—generally convolutions, different fractional operators are defined on different ranges of integration for the convolutions. The semi-infinite Weyl derivatives are particularly easy to deal with, because they are simply power law filters in Fourier space, see equation (12) (for more information on fractional equations, see, e.g., Miller and Ross Reference Miller and Ross1993; Podlubny Reference Podlubny1999).
The γ's in equation (2) are Gaussian white noises; they are proportional to “unit” white noises, γ. Unit white noises have the properties:
$$\left\langle {\gamma ( {t_1} ) \gamma ( {t_2} ) } \right\rangle = \delta ( {t_1-t_2} ) ; \;\,\quad \matrix{ {\left\langle {\gamma^2} \right\rangle = 1; \;} & {\left\langle \gamma \right\rangle = 0} \cr } $$where the angle brackets indicate ensemble (statistical) averaging, and δ is the Dirac function. Equation (2) therefore implies that T, E + are fractional integrals of white noises. Depending on the values of the exponents, these are fractional Gaussian noises (fGns) and fractional Brownian motions (fBms) (Mandelbrot and Van Ness Reference Mandelbrot and Van Ness1968; see the later discussion on the small- and large-scale limits).
Completing the Model, the Diagnostic Equation
Equations (1) and (2) determine D, E +, T. However, for the model to determine E and O, we need a final equation for E −:
$$E_-= \tau _D^{} \displaystyle{{dD} \over {dt}}; \;\quad E_- = O-E$$This is just the differential form of the usual discrete-time definition of diversity: D j+1 = D j(1 + O j − E j), where j is a time index. τD is the discretization time, the basic resolution of the series. Equation (5) plays no role in the dynamics, conventionally, it defines D (see eq. 34 for its Fourier expression). Mathematically, equation (5) is thus a “diagnostic equation,” because it simply allows us to close (complete) the model by determining O, E:
$$\matrix{ {O = ( {E_ + + E_-} ) /2} \cr {E = ( {E_ + -E_-} ) /2} \cr } $$A schematic of the full model is given in Figure 1 showing its various parts, including the possibility of deterministic forcing discussed later (eq. 12).
Discussion
Diversity as an ffRn
The diversity model was written in a nonstandard way (eq. 1), because in this form, its basic behavior is transparent. When h > 0, the fractional term is the highest-order derivative; at high frequencies it therefore dominates the zeroth-order (D) term, such that at short lags, Δt < τ, diversity fluctuations ΔD ∝ ΔT such that D follows the temperature. However, at low frequencies (Δt > τ), the zeroth-order term dominates, and we have instead ΔD ∝ ΔE +. By inspection, the model therefore reproduces the crossover at lag τ, and the crossover will be scaling due to the scaling of T, E + (eq. 2). The mathematical structure of the model is clearer if we substitute the drivers in terms of their own Gaussian forcings γT, γE (eq. 2), rewriting equation 1 as:
$$\tau ^h\displaystyle{{d^hD} \over {dt^h}} + D = \tau ^{-\alpha }\displaystyle{{d^{-\alpha }} \over {dt^{-\alpha }}}( {s_T\gamma_T + s_E\gamma_E} ) $$ The term on the right-hand side represents the total forcing of the diversity. (Note: d −α/dt −α is a fractional integral order α: for Weyl derivative and integrals it is the inverse of the α-order derivative 
$d^\alpha /dt^\alpha $).
The linear combination of white noises s Tγ T + s Eγ E is also a white noise. The D equation is thus an order h-ordered fractional relaxation equation forced by an order α fractionally integrated white noise, that is, it is a “fractionally integrated fractional relaxation” process (ffRn; Lovejoy Reference Lovejoy2022a). The basic “unit” ffRn process Uh ,α(t) satisfies:
$$\left({\displaystyle{{d^{h + \alpha }} \over {dt^{h + \alpha }}} + \displaystyle{{d^\alpha } \over {dt^\alpha }}} \right)U_{\alpha , h} = \gamma $$where γ is the unit white noise defined earlier (eq. 4), and we have used the fact that for Weyl fractional derivatives, fractional differentiation and integration commute. If time is rescaled (t→t/τ), we see (from eq. 7) that D is proportional to U α,h. We note that if h = 1, the D equation (eq. 1) would be a classical relaxation equation, and if forced by a white noise (i.e., if α = 0), D would be a classical Ornstein-Uhlenbeck (OU) process. OU processes are thus the classical special cases of crossover phenomena involving high-frequency processes with exponential decorrelations (e.g., Markov processes) that lead to white noise behavior at low frequencies. They are currently the conventional approaches to the modeling and analysis of microevolutionary as well as macroevolutionary dynamics (Khabbazian et al. Reference Khabbazian, Kriebel, Rohe and Ane2016; Bartoszek et al. Reference Bartoszek, Glémin, Kaj and Lascoux2017; Liow et al. Reference Liow, Uyeda and Hunt2022).
Deterministic Behavior: Impulse Response Functions
The D process—the solution to equation (7)—is the response of the operator 
$\left({{{d^{h + \alpha }} \over {dt^{h + \alpha }}} + {{d^\alpha } \over {dt^\alpha }}} \right)$ to a white noise forcing. The general behavior of responses to linear operators is determined by their impulse response (Green's) functions G α,h that satisfy:
$$\left({\displaystyle{{d^{h + \alpha }} \over {dt^{h + \alpha }}} + \displaystyle{{d^\alpha } \over {dt^\alpha }}} \right)G_{\alpha , h} = \delta ( t ) $$(Lovejoy Reference Lovejoy2022a), where δ(t) is the Dirac (“delta”) function. G α,h can be expressed in terms of “generalized exponentials” or Mittag-Leffler functions e h,h+α as:
$$\matrix{ G_{\alpha , h} ( t ) = t^{h-1 + \alpha} e_{h, h + \alpha} ( \!\!-\!t^h ) = t^{\alpha -1} \sum\limits_{n = 1}^\infty ( {-1} ) ^{n + 1} \displaystyle{t^{nh} \over \Gamma ( \alpha + nh ) }; & t \ge 0 \hfill \cr{} & 0 < \alpha < 1/2; \quad 0 < h < 2 \hfill \cr G_{\alpha , h} ( t ) = 0; \hfill & {t < 0}\hfill\cr } $$
$$e_{a, b}( z ) = \sum\limits_{n = 0}^\infty {\displaystyle{{z^n} \over {\Gamma ( {an + b} ) }}} $$ At small t, the leading order term is therefore 
$G_{\alpha , h}( t ) \approx {{t^{\alpha -1 + h}} \over {\Gamma ( {\alpha + h} ) }}$. The large t (asymptotic) expansion is:
$$\matrix{ {G_{\alpha , h}^{} ( t ) = t^{\alpha -1}\mathop \sum \limits_{n = 0}^\infty \displaystyle{{{( {-1} ) }^n} \over {\Gamma ( {\alpha -nh} ) }}t^{{-}nh}; \;} & {t > > 1} \cr } $$ (Podlubny Reference Podlubny1999). Whereas the small t expansion is t α−1 times terms of positive powers of h, the large t expansion is in terms of t α−1 times terms in negative powers of h, with leading term 
$G_{\alpha , h}^{} ( t ) = {{t^{\alpha -1}} \over {\Gamma ( \alpha ) }}$. Unless h = 0, G α,h(t) therefore transitions between two different power laws. The special case h = 0 that applies to the temperature and turnover forcings (eq. 2), corresponds to the pure power law 
$G_{\alpha , 0}^{} ( t ) = {{t^{\alpha -1}} \over {\Gamma ( \alpha ) }}$. G α,h has the property that if it is (fractionally) integrated ζ times, the result is just G α+ζ,h. As explained in Appendix 1, G α,h is useful for numerical simulations.
At a typical highest resolution of global datasets with timescales ≈ of 1 Myr, G α,h gives the deterministic response to forcings that are effectively impulses at this scale, for example, a bolide strike (Alvarez et al. Reference Alvarez, Alvarez, Asaro and Michel1980; During et al. Reference During, Smit, Voeten, Berruyer, Tafforeau, Sanchez, Stein, Verdegaal-Warmerdam and Van Der Lubbe2022), supernova or gamma ray burst (Fields et al. Reference Fields, Melott, Ellis, Ertel, Fry, Lieberman, Liu, Miller and Thomas2020), or even much slower hyperthermal event such as the Paleocene–Eocene thermal maximum (Gingerich Reference Gingerich2006; McInerney and Wing Reference McInerney and Wing2011) or the Cenomanian–Turonian event (Eaton et al. Reference Eaton, Kirkland, Hutchison, Denton, O'Neill and Parrish1997; Meyers et al. Reference Meyers, Siewert, Singer, Sageman, Condon, Obradovich, Jicha and Sawyer2012; Venckutė-Aleksienė et al. Reference Venckutė-Aleksienė, Spiridonov, Garbaras and Radzevičius2018), extensive volcanic eruption episodes, or other short-timescale (below ≈ 1 Myr) stressors. Figure 2 shows a plot of impulse responses for temperature and turnover demonstrating their singular nature for the empirical parameters estimated in SL (α ≈ 0.25, h ≈ 0.5). These singular responses combine apparently contradictory features: on the one hand, the falloff at short times following the impulse is very sharp, whereas on the other hand, the decay is very slow at long times, so its effects take a long time to disappear. The sharpness feature is desirable; because mass extinctions, and potentially other episodes of dramatic biotic change, effectively represent periods of almost infinite turnover rate, it is near instantaneous or singular-like (Foote Reference Foote1994, Reference Foote2005). Indeed, the global stratigraphic stages and substages are based on the episodes of turnover. Figure 3 shows the impulse response for the diversity that has two different power law regimes: a high-frequency t α−1+h regime and a low-frequency t α−1 regime (the leading terms in eqs. 10, 11) with a transition at various scales τ indicated. As expected, the former (temperature dominance) corresponds to the singularity t −0.25 and the latter to the singularity t −0.75 (turnover dominance; both are shown in Fig. 2). Note that because the equations are linear, the impulse responses from the deterministic forcings shown in Figures 2 and 3 will be superposed onto the stochastic white noise–driven responses that we calculate in the section entitled “Numerical Simulations.”
Figure 2. The impulse (delta function) response 
$G_{\alpha , 0}^{} ( t ) = t^{\alpha -1}/\Gamma ( \alpha ) $ for fractional integrals of order α normalized for the same response after 1 Myr. The bottom corresponds to the turnover (E +) response α = ¼, and the top corresponds to the temperature (T) response with α = ¾. Notice the long-term effects. Due to causality, the impulse response is 0 for t < 0.
Figure 3. The impulse response G α,h(t/τ), with α = ¼, h = ½, corresponding to the diversity (D) response, for critical transition times τ = 1, 4, 16, 64, and 256 Myr (bottom to top). The empirical value is τ ≈ 40 Myr (SL). Due to causality, the impulse response is 0 for t < 0.
The other aspect of the singular impulse responses is their long time decays that are much slower than conventional exponential decays and can be thought of as a kind of memory that characterizes the responses to both deterministic and stochastic forcing. Our model thus predicts that there will be long-term consequences of bolide impacts or other catastrophic events. This is in accord, for example, with the findings of Krug et al. (Reference Krug, Jablonski and Valentine2009) and Krug and Jablonski (Reference Krug and Jablonski2012) that the Cretaceous–Paleogene (K-Pg) mass extinction caused by the effects of the Chicxulub asteroid impact changed long-term origination rates and their spatial distribution, a situation that persists today, 66 Myr after the event, in accord with this long-memory feature of the FMEM model. The slow decay in the response is also a desired property in modeling macroevolutionary dynamics, as it can replicate effects of phylogenetic inertia or “phylogenetic constraint” (Gould Reference Gould2002) and also inertia of persistence of geological structures that can affect dynamics of biodiversity for periods, eras, or even eons. The genetic composition is basically the material trans-generational memory of biological individuals of all levels of the Linnaean hierarchy (Eldredge Reference Eldredge, Jablonski, Erwin and Lipps1996). Therefore, extinction (or origination) of genera or taxa of higher ranks can change the functioning and the properties of the biosphere for hundreds of millions of years. The same goes for the formation of such structures as oceanic basins or continents with their myriads of possible configurations—their origins and disappearances impose new boundary conditions for geophysical and macroevolutionary dynamics for hundreds of millions of years (Nance Reference Nance2022; Spiridonov et al. Reference Spiridonov, Balakauskas and Lovejoy2022).
We could further note that the long memory can be exploited to make future predictions. This is because for Gaussian forcing (eq. 2), E + and the increments of T are long-memory fGn processes and—as discussed earlier—D is an ffRn. The predictability of the former has been mathematically studied in Gripenberg and Norros (Reference Gripenberg and Norros1996) and has been exploited for monthly and seasonal forecasts in Del Rio Amador (2019, 2021) and Lovejoy and Del Rio Amador (Reference Lovejoy and Del Rio Amador2023), and the predictability properties of ffRn processes have been studied in Lovejoy (2022).
Solving the Model
Readers only interested in results can skip ahead until the basic model equations (38, 39) that are given at the end of the section “Full model statistics”.
Fractional derivatives are generally convolutions (with power laws, eq. 3), hence different ranges of integration in the convolution yield different fractional derivatives. Most often (e.g., the Riemann-Liouville and Caputo fractional derivatives), the latter are taken from time = 0 to t, in which case the initial conditions are important and dealing with them is technically somewhat complex. In these cases, the main tool is the Laplace transform. Here, however, we consider statistically stationary white noise forcing that starts at time = −∞. In this case, we can use the “Weyl” fractional derivative (a convolution from −∞ to t, eq. 3) whose Fourier transform (“F.T.”) is particularly simple:
$$\displaystyle{{d^h} \over {dt^h}}\mathop {\mathop \leftrightarrow \limits^{F.T.} {( {i\omega } ) }^h}\limits^{} $$where ω is the Fourier conjugate variable, that is, the frequency (when h is an integer, the formula may be familiar to the reader). If we Fourier transform equations (1) and (2), we can write the model in matrix form as:
$$\underline {\widetilde{S}} ( \omega ) = ( {i\omega \tau } ) ^{-\alpha }\widetilde{{\underline{\underline F} }}( \omega ) \left({\matrix{ {\matrix{ 1 \cr 0 \cr {s_T} \cr } } & {\matrix{ 0 \cr 1 \cr {s_E} \cr } } \cr } } \right)\left({\matrix{ {{\widetilde{\gamma }}_T} \cr {{\widetilde{\gamma }}_E} \cr } } \right)$$
$$\widetilde{{\underline S }}( \omega ) = \left({\matrix{ {\widetilde{T}} \cr {\widetilde{{E_ + }}} \cr {\widetilde{D}} \cr } } \right)$$
$$\underline{\underline F} ( \omega ) = \left({\matrix{ {{( {i\omega \tau } ) }^{{-}h}} & 0 & 0 \cr 0 & 1 & 0 \cr 0 & 0 & {\displaystyle{1 \over {1 + {( {i\omega \tau } ) }^h}}} \cr } } \right)$$(the single underlining indicates a vector, the double underlining, a matrix, the Fourier transform is denoted with a tilde).
As noted earlier, the D forcing is a linear combination of white noises (eq. 7), such that the sum on the right-hand side of equations (7) and (13) is a correlated white noise. However, from the data (see fig. 5), we see that E + and T are themselves correlated. We therefore rewrite the model in terms of two statistically independent (
$\left\langle {\gamma_1^{} \gamma_2^{} } \right\rangle = 0$) unit (
$\left\langle {\gamma_1^2 } \right\rangle = \left\langle {\gamma_2^2 } \right\rangle = 1$) white noise drivers γ1, γ2:
$$\left({\matrix{ {\gamma_T} \cr {\gamma_E} \cr } } \right) = \left({\matrix{ {\sigma_T} & 0 \cr 0 & {\sigma_E} \cr } } \right)\left({\matrix{ 1 & 0 \cr {\rho_E} & {\sqrt {1-\rho_E^2 } } \cr } } \right)\left({\matrix{ {\gamma_1} \cr {\gamma_2} \cr } } \right)$$So that:
$$\matrix{ {\sigma _T^2 = \left\langle {\gamma_T^2 } \right\rangle ; \;} & {\sigma _E^2 = \left\langle {\gamma_E^2 } \right\rangle ; \;} & {\matrix{ {\rho _E^{} = \displaystyle{{\left\langle {\gamma_T^{} \gamma_E^{} } \right\rangle } \over {\sigma _T^{} \sigma _E^{} }}; \;} & {\left\langle {\gamma_T^{} } \right\rangle = \left\langle {\gamma_E^{} } \right\rangle = 0} \cr } } \cr } $$where σT is the standard deviation of γT, σE of γE, and ρE is the T, E + correlation. Equation (14) is the standard Cholesky decomposition for two correlated random variables, noises.
Fourier transforming equation (14) and using equation (13), we can write the model as:
$$\underline {\widetilde{S}} ( \omega ) = ( {i\omega \tau } ) ^{-\alpha }\underline{\underline F} ( \omega ) \underline{\underline \sigma } \underline{\underline \rho } \underline{\underline {\widetilde{\gamma }}} $$
$$\underline{\underline \sigma } = \left({\matrix{ {\sigma_T^{} } & 0 & 0 \cr 0 & {\sigma_E^{} } & 0 \cr 0 & 0 & {\sigma_D^{} } \cr } } \right); \;\quad \underline{\underline \rho } = \left({\matrix{ {\matrix{ 1 \cr {\rho_E} \cr {\rho_D} \cr } } & {\matrix{ 0 \cr {\sqrt {1-\rho_E^2 } } \cr {{\mathop{\rm sgn}} ( r ) \sqrt {1-\rho_D^2 } } \cr } } \cr } } \right); \;\quad \underline {\widetilde{\gamma }} = \left({\matrix{ {{\widetilde{\gamma }}_1} \cr {{\widetilde{\gamma }}_2} \cr } } \right)$$Where the parameters:
$$\matrix{ {\sigma _D^{} = s_T\sigma _T\sqrt {1 + 2\rho _Er + r^2} ; \;} & {r = \displaystyle{{s_E\sigma _E} \over {s_T\sigma _T}}} \cr {\rho _D = \displaystyle{{1 + r\rho _E} \over {\sqrt {1 + 2r\rho _E + r^2} }}} & {} \cr } $$depend on both the driver statistics (σT, σE, and ρE) and the model parameters sT, sE. While σD does parametrize the amplitude of the diversity fluctuations, unlike σT, σE (which must be ≥ 0), it is not a true standard deviation: if sT < 0, it will be negative. Similarly, we will see that ρD determines the D, E+, and T correlations but is not itself a correlation coefficient and depends on the sign of the ratio r.
Stochastic Response to White Noise Forcing
Scaling Processes: fGn and fBm
We are interested in the statistical properties of the solutions 
$\widetilde{{\underline S }}( \omega ) $. These can be expressed in terms of fGn, fBm, and ffRn processes. Before discussing the full statistics that include the cross-correlations, let us therefore discuss the statistics of each of the main variables individually.
Let us start with the scaling processes T, E + that are of the form:
$$\matrix{ {\mathop {\displaystyle{{d^{\alpha _X}X} \over {dt^{\alpha _X}}}}\limits^{} = \gamma \mathop \leftrightarrow \limits^{F.T.} {( {i\omega } ) }^{\alpha _X}\widetilde{X} = \widetilde{\gamma }} & {} \cr } $$For the statistics, we can determine the power spectrum:
$$\matrix{ {E_X( \omega ) = \left\langle {{\vert {\widetilde{X}} \vert }^2} \right\rangle = \displaystyle{1 \over {2\pi }}{\vert \omega \vert }^{-\beta _X}; \;} & {\beta _X = 2\alpha _X} \cr } $$where βX is the spectral exponent, and we have used the fact that the spectrum of a Gaussian white noise is flat:
$$\left\langle {{\vert {\widetilde{\gamma }( \omega ) } \vert }^2} \right\rangle = \displaystyle{1 \over {2\pi }}\left\langle {\gamma^2} \right\rangle = \displaystyle{1 \over {2\pi }}$$E X(ω) is thus the basic form of the T, E + spectra (not to be confused with the extinction rate E). From the Wiener-Khintchin theorem, the (real-space) autocorrelation function RX(Δt) is the inverse transform:
$$\eqalign{\matrix{&{R_X( {\Delta t} ) = \left\langle {X( t ) X( {t-\Delta t} ) } \right\rangle \propto \Delta t^{H_X}\mathop \leftrightarrow \limits^{F.T.} {\widetilde{R}}_X( \omega )} \cr&\hskip1pc= E_X( \omega ) = \left\langle {{\vert {\widetilde{{X( \omega ) }}} \vert }^2} \right\rangle \propto {\vert \omega \vert }^{-\beta _X}; \; \cr&\hskip-1.5pc{H_X = \displaystyle{{\beta _X-1} \over 2} = \displaystyle{{\alpha _X} \over 2}} \cr } }$$ The technical difficulty is that due to a low-frequency divergence, the inverse transform of pure power spectra (eq. 19) only converges for βX < 1 (i.e., αX < ½, HX < 0); this is the fGn regime appropriate for E +. Even here, R X(Δt) is infinite for Δt = 0. Because R X(0) is the variance, fGn processes are (like the white noise special case αX = 0) generalized functions that must be averaged (integrated) over finite intervals in order to represent physical processes. Averaging to yield a finite resolution process is adequate for βX > −1 (αX > −½, HX > −1) such that the fGn process is defined for −1 < βX < 1 (i.e., −½ < αX < ½, −1 < HX < 0). After X is averaged over a finite resolution τr, the result is Xτr with root-mean-square (RMS) statistics 
$\left\langle {X_{\tau_r}^2 } \right\rangle ^{1/2}\propto \tau _r^{H_x} $. Because HX < 0, the empirical statistics will be highly sensitive to the resolution τr.
When αX ≥ ½, the low-frequency divergences imply that the X(t) process is nonstationary (the process generally “wanders off” to plus or minus infinity). However, for 1 < βX < 3 (i.e., ½ < αX < 3/2, 0 < HX < 1; this is the range appropriate for T: HT ≈ 0.25, βT ≈ 3/2), its increments are (stationary) fGn processes; this regime defines the fBm process. Finally, because all physical scaling processes exist over finite ranges of scale, there will be finite outer (longest) timescale (smallest frequency) such that truncating the spectrum at low frequencies (as for the ffRn processes, see the section “Two Scaling Regimes: fRn and ffRn”) leads to an overall stationary process.
When analyzing paleo-series, it is convenient to analyze the statistics in real space, the main reason being that these are easier to interpret (the difficulty in interpretation is the cause of the recently discovered quadrillion error in climate temperature spectra [Lovejoy Reference Lovejoy2015]). An additional reason is that uniformly sampled paleo-series are typically not available: indeed, the geochronologies themselves are typically scaling (see Appendix 3 and Fig. A3.1 for more discussion). For the moment, the problem of spectral estimation from data with scaling measurement densities (i.e., scaling number of measurements per time interval) is an unsolved problem.
We have already noted that the autocorrelation functions are only adequate for HX < 0 (αX < ½, βX < 1), this is why, when 0 < HX < 1, it is conventional to define fluctuations using differences ΔX(Δt) = X(t − Δt) − X(t), which are stationary over this range. Differences avoid low-frequency divergences but will still have high-frequency divergences when HX < 0. To avoid problems at both the small scale (resolution dependencies) and large scales (nonstationarity), it is convenient to use Haar fluctuations. Over the interval Δt, the Haar fluctuation ΔX(Δt) is defined as the difference between the average of the first and second halves of the interval.
$$\matrix{ {{\left\langle {\Delta X{( {\Delta t} ) }^2} \right\rangle }^{1/2}\propto \Delta t^{H_X}\mathop {\leftrightarrow E_X( \omega ) }\limits^{} \propto \omega ^{-\beta _X}; \;} & {\matrix{ {-1 < H_X < 1} \cr {-1 < \beta _X < 3} \cr {\beta _X = 2H_X + 1} \cr } } \cr } $$(valid for Haar fluctuations). Over the indicated range of parameters, the Haar fluctuations are stationary and are independent of the resolution.
Comparing equations (7) and (2), we find:
$$\matrix{ {\matrix{ {{\left\langle {\Delta E_ + {( {\Delta t} ) }^2} \right\rangle }^{1/2}\propto \Delta t^{H_E}; \;} & {H_E = \alpha -\displaystyle{1 \over 2}} \cr } } \cr {\matrix{ {{\left\langle {\Delta T{( {\Delta t} ) }^2} \right\rangle }^{1/2}\propto \Delta t^{H_T}; \;} & {H_T = h + \alpha -\displaystyle{1 \over 2}} \cr } } \cr } $$Two Scaling Regimes: fRn and ffRn
From equations (8) and (9), the basic Fourier transforms of ffRn processes and their impulse responses are:
$$\!\!\!\!\!\!\!\!\matrix{ \widetilde{U}_{\alpha , h}( \omega ) = \displaystyle{\widetilde{\gamma} \over ( i\omega )^\alpha ( 1 + ( i\omega ) ^h )} ; \cr \widetilde{G}_{\alpha , h} ( \omega ) = \displaystyle{1 \over ( i\omega ) ^\alpha ( 1 + ( i\omega ) ^h )} ;}\quad 0 < \alpha < 1/2; \quad 0 < h < 2 $$The fractional relaxation noise (fRn) process is the special case where α = 0. The ffRn power spectrum is therefore:
$$E_{\alpha , h}( \omega ) = \left\langle {{\vert {{\widetilde{U}}_{\alpha , h}} \vert }^2} \right\rangle = \displaystyle{1 \over {2\pi {\vert \omega \vert }^{2\alpha }{\vert {1 + {( {i\omega } ) }^h} \vert }^2}}$$E α,h(ω) is thus the basic form of the D spectrum.
The full statistical properties of ffRn processes (including series expansions) are discussed in Lovejoy (2022); however, for our purposes, the low- and high-frequency scaling exponents are sufficient. For these, equation (25) yields:
$$\matrix{ {E_{\alpha , h}( \omega ) \propto {\vert \omega \vert }^{-\beta }; \;} & {\matrix{ {\beta _l = 2\alpha ; \;} & {\omega < < 1} \cr {\beta _h = 2( {\alpha + h} ) ; \;} & {\omega > > 1} \cr } } \cr } $$(“h” for high frequency, “l” for low frequency). To obtain the basic fluctuation statistics, it is sufficient to apply equation (22) to each regime separately. Indeed, more generally, “Tauberian theorems” (e.g., Feller Reference Feller1971) imply that if the spectrum is a power law over a wide enough range, then the corresponding (second-order) real-space statistics will also be scaling. Therefore:
$$\matrix{ {{\left\langle {\Delta U_{\alpha , h}{( {\Delta t} ) }^2} \right\rangle }^{1/2}} & {\matrix{ {\propto \Delta t^{H_l}; \;} & {\matrix{ {H_l = \alpha -\displaystyle{1 \over 2}; \;} & {\Delta t > > 1} \cr } } \cr {\propto \Delta t^{H_h}; \;} & {\matrix{ {H_h = \alpha + h-\displaystyle{1 \over 2}; \;} & {\Delta t < < 1} \cr } } \cr } } \cr } $$Using the empirical values α ≈ 0.25, h ≈ 0.5, we see that E + is a fractional Gaussian noise and that T is an fBm process. Also, we find (cf. eqs. 7, 27) that HD ≈ HT (Δt ≪ τ) and HD ≈ HE (Δt ≫ τ).
The Full-Model Statistics: Spectra, Correlations
The Basic Model
The model is linear and has stationary Gaussian (white noise) forcing (T, E +), therefore D, E −, E, O are also Gaussian, such that their statistics are determined by spectra and cross-spectra, or equivalently in real space (via the Wiener-Khintchin theorem), by the autocorrelations and cross-correlations:
$$R_{ij}( {\Delta t} ) = \left\langle {S_i( t ) S_j( {t-\Delta t} ) } \right\rangle \mathop \leftrightarrow \limits^{F.T.} \widetilde{R}_{ij}( \omega ) = \left\langle {{\widetilde{S}}_i( \omega ) \widetilde{S}_j^\ast ( \omega ) } \right\rangle $$(the diagonal terms are the spectra of the components: 
$\widetilde{R}_{ii}( \omega ) = E_i( \omega ) $, the asterisk indicates the complex conjugate). In matrix notation:
$$\eqalign{\underline{\underline {\widetilde{R}}} ( \omega ) &= \left\langle {\widetilde{{\underline S }}{\widetilde{{\underline S }}}^{T\ast }} \right\rangle = \vert {\omega \tau } \vert ^{{-}2\alpha }\underline{\underline F} ( \omega ) \underline{\underline \sigma } \underline{\underline \rho } \left\langle {\underline {\widetilde{\gamma }} {\underline {\widetilde{\gamma }} }^{T\ast }} \right\rangle \underline{\underline \rho } ^{T\ast }\underline{\underline \sigma } ^{T\ast }\underline{\underline F} ( \omega ) ^{T\ast } \cr &= \displaystyle{{{\vert {\omega \tau } \vert }^{{-}2\alpha }} \over {2\pi }}\underline{\underline F} ( \omega ) \underline{\underline \sigma } \underline{\underline \rho } \underline{\underline \rho } ^{T\ast }\underline{\underline \sigma } ^{T\ast }\underline{\underline F} ( \omega ) ^{T\ast }} $$where the superscript T indicates the transpose, the asterisk indicates the complex conjugate, and we have used:
$$\left\langle {\underline {\widetilde{\gamma }}.{\underline {\widetilde{\gamma }} }^{{\ast} T}} \right\rangle = \left\langle {\left({\matrix{ {{\widetilde{\gamma }}_1} \cr {{\widetilde{\gamma }}_2} \cr } } \right)\left({\matrix{ {{\widetilde{\gamma }}_1} & {{\widetilde{\gamma }}_2} \cr } } \right)} \right\rangle = \displaystyle{1 \over {2\pi }}\left({\matrix{ 1 & 0 \cr 0 & 1 \cr } } \right) = \displaystyle{1 \over {2\pi }}\underline{\underline 1} $$The key correlation matrix (from eq. 16) is:
$$\underline{\underline \rho } \underline{\underline \rho } ^{T\ast } = \left({\matrix{ 1 & {\rho_{TE}} & {\rho_{TD}} \cr {\rho_{TE}} & 1 & {\rho_{ED}} \cr {\rho_{TD}} & {\rho_{ED}} & 1 \cr } } \right)$$where
$$\matrix{ {\rho _{TE} = \rho _E; \;} & {\rho _{TD} = \rho _D; \;} & {\rho _{ED} = \rho _E\rho _D + {\mathop{\rm sgn}} ( r ) \sqrt {1-\rho _E^2 } \sqrt {1-\rho _D^2 } } \cr } $$and
$$\underline{\underline \sigma } \underline{\underline \rho } \underline{\underline \rho } ^{T\ast }\underline{\underline \sigma } ^{T\ast } = \left({\matrix{ {\sigma_T^2 } & {\rho_{TE}\sigma_T^{} \sigma_E^{} } & {\rho_{TD}\sigma_D\sigma_T^{} } \cr {\rho_{TE}\sigma_T^{} \sigma_E^{} } & {\sigma_E^2 } & {\rho_{DE}\sigma_E^{} \sigma_D} \cr {\rho_{TD}\sigma_D\sigma_T^{} } & {\rho_{DE}\sigma_E^{} \sigma_D} & {\sigma_D^2 } \cr } } \right)$$Completing the Model: The Diagnostic Equation for E−
Before writing down the final spectra, let us complete the system with the help of the diagnostic equation that allows us to determine E − from D (and hence E, O, eq. 6).
The Fourier transform of the diagnostic equation (eq. 5) is:
$$\widetilde{E}_- = \left({\displaystyle{{\tau_D^{} } \over \tau }} \right)( {i\omega \tau } ) \widetilde{D}$$Therefore the full system is:
$$\left({\matrix{ {\widetilde{T}} \cr {\widetilde{{E_ + }}} \cr {\widetilde{D}} \cr {\widetilde{{E_-}}} \cr } } \right) = ( {i\omega \tau } ) ^{-\alpha }\left({\matrix{ {{( {i\omega \tau } ) }^{{-}h}} & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & {\displaystyle{1 \over {1 + {( {i\omega \tau } ) }^h}}} & 0 \cr 0 & 0 & 0 & {\displaystyle{{i\omega \tau } \over {1 + {( {i\omega \tau } ) }^h}}} \cr } } \right)\left({\matrix{ {\sigma_T^{} } & 0 & 0 & 0 \cr 0 & {\sigma_{E_ + }^{} } & 0 & 0 \cr 0 & 0 & {\sigma_D^{} } & 0 \cr 0 & 0 & 0 & {\displaystyle{{\tau_D} \over \tau }} \cr } } \right)\left({\matrix{ {\matrix{ 1 \cr {\rho_E} \cr {\rho_D} \cr {\rho_D} \cr } } & {\matrix{ 0 \cr {\sqrt {1-\rho_E^2 } } \cr {{\mathop{\rm sgn}} ( r ) \sqrt {1-\rho_D^2 } } \cr {{\mathop{\rm sgn}} ( r ) \sqrt {1-\rho_D^2 } } \cr } } \cr } } \right)\left({\matrix{ {\gamma_1} \cr {\gamma_2} \cr } } \right)$$From this, we can find E, O:
$$\matrix{ {\widetilde{E} = \displaystyle{1 \over 2}( {{\widetilde{E}}_ + {-}{\widetilde{E}}_-} ) } \cr {\widetilde{O} = \displaystyle{1 \over 2}( {{\widetilde{E}}_ + { + } {\widetilde{E}}_-} ) } \cr } $$The explicit formulae for E ± are:
$$\matrix{ {{\widetilde{E}}_ + { \,=\, } {( {i\omega \tau } ) }^{-\alpha }\left[{\sigma_{E_ + }\left({\rho_E{\widetilde{\gamma }}_1 + \sqrt {1-\rho_E^2 } {\widetilde{\gamma }}_2} \right)} \right]} \cr {{\widetilde{E}}_-{ \,=\, } {( {i\omega \tau } ) }^{-\alpha }\displaystyle{{i\omega \tau _D} \over {1 + {( {i\omega \tau } ) }^h}}\left({\rho_D{\widetilde{\gamma }}_1 + {\mathop{\rm sgn}} ( r ) \sqrt {1-\rho_D^2 } {\widetilde{\gamma }}_2} \right)} \cr } $$The overall final statistics are:
$$\eqalign{& \underline{\underline {\widetilde{R}}} ( \omega ) = \left({\matrix{ {\left\langle {{\vert {\widetilde{T}} \vert }^2} \right\rangle } & {\left\langle {\widetilde{T}\widetilde{{E_ + ^\ast }}} \right\rangle } & {\left\langle {\widetilde{T}\widetilde{{D^\ast }}} \right\rangle } & {\left\langle {\widetilde{T}\widetilde{{E_-^\ast }}} \right\rangle } \cr {\left\langle {\widetilde{{E_ + }}\widetilde{{T^\ast }}} \right\rangle } & {\left\langle {{\vert {\widetilde{{E_ + }}} \vert }^2} \right\rangle } & {\left\langle {\widetilde{{E_ + }}\widetilde{{D^\ast }}} \right\rangle } & {\left\langle {\widetilde{{E_ + }}\widetilde{{E_-^\ast }}} \right\rangle } \cr {\left\langle {\widetilde{D}\widetilde{{T^\ast }}} \right\rangle } & {\left\langle {\widetilde{D}\widetilde{{E_ + ^\ast }}} \right\rangle } & {\left\langle {{\vert {\widetilde{D}} \vert }^2} \right\rangle } & {\left\langle {\widetilde{D}\widetilde{{E_-^\ast }}} \right\rangle } \cr {\left\langle {\widetilde{{E_-}}\widetilde{{T_{}^\ast }}} \right\rangle } & {\left\langle {\widetilde{{E_-}}\widetilde{{E_-^\ast }}} \right\rangle } & {\left\langle {\widetilde{{E_-}}\widetilde{{D^\ast }}} \right\rangle } & {\left\langle {{\vert {\widetilde{{E_-}}} \vert }^2} \right\rangle } \cr } } \right)\cr & = \vert {\omega \tau } \vert ^{{-}2\alpha }\left({\matrix{ {{\vert {\omega \tau } \vert }^{{-}2h}\sigma_T^2 } & {{( {i\omega \tau } ) }^{{-}h}\rho_{TE}\sigma_T^{} \sigma_E^{} } & {\displaystyle{{\rho_{TD}\sigma_D\sigma_T} \over {{( {i\omega \tau } ) }^h( {1 + {( {-i\omega \tau } ) }^h} ) }}} & {\displaystyle{{\rho_{TD}\tau_D\sigma_T{( {-i\omega \tau } ) }^{1-h}} \over {\tau ( {1 + {( {i\omega \tau } ) }^h} ) }}} \cr {{( {-i\omega \tau } ) }^{{-}h}\rho_{TE}\sigma_T^{} \sigma_E^{} } & {\sigma_E^2 } & {\displaystyle{{\rho_{ED}\sigma_E\sigma_D} \over {( {1 + {( {-i\omega \tau } ) }^h} ) }}} & {\displaystyle{{\rho_{ED}\tau_D\sigma_E( {-i\omega \tau } ) } \over {\tau ( {1 + {( {i\omega \tau } ) }^h} ) }}} \cr {\displaystyle{{\rho_{TD}\sigma_D\sigma_T} \over {{( {-i\omega \tau } ) }^h( {1 + {( {i\omega \tau } ) }^h} ) }}} & {\displaystyle{{\rho_{ED}\sigma_E\sigma_D} \over {( {1 + {( {i\omega \tau } ) }^h} ) }}} & {\displaystyle{{\sigma_D^2 } \over {{\vert {1 + {( {i\omega \tau } ) }^h} \vert }^2}}} & {\displaystyle{{\sigma_D^{} \tau_D( {-i\omega \tau } ) } \over {\tau {\vert {1 + {( {i\omega \tau } ) }^h} \vert }^2}}} \cr {\displaystyle{{\rho_{TD}\tau_D\sigma_T{( {i\omega \tau } ) }^{1-h}} \over {\tau ( {1 + {( {i\omega \tau } ) }^h} ) }}} & {\displaystyle{{\rho_{ED}\tau_D\sigma_E( {i\omega \tau } ) } \over {\tau ( {1 + {( {i\omega \tau } ) }^h} ) }}} & {\displaystyle{{\sigma_D^{} \tau_D( {i\omega \tau } ) } \over {\tau {\vert {1 + {( {i\omega \tau } ) }^h} \vert }^2}}} & {\displaystyle{{\tau_D^2 {\vert {\omega \tau } \vert }^2} \over {\tau_{}^2 {\vert {1 + {( {i\omega \tau } ) }^h} \vert }^2}}} \cr } } \right)} $$Using equations (36) and (37), the spectra of E, O can be determined:
$$\matrix{ {\left\langle {{\vert {\widetilde{E}} \vert }^2} \right\rangle = \displaystyle{1 \over 4}\left({\left\langle {{\vert {\widetilde{{E_ + }}} \vert }^2} \right\rangle + \left\langle {{\vert {\widetilde{{E_-}}} \vert }^2} \right\rangle -2\left\langle {\widetilde{{E_ + }}\widetilde{{E_-^\ast }}} \right\rangle } \right)\approx \displaystyle{1 \over 4}\left\langle {{\vert {\widetilde{{E_ + }}} \vert }^2} \right\rangle } \cr {\left\langle {{\vert {\widetilde{O}} \vert }^2} \right\rangle = \displaystyle{1 \over 4}\left({\left\langle {{\vert {\widetilde{{E_ + }}} \vert }^2} \right\rangle + \left\langle {{\vert {\widetilde{{E_-}}} \vert }^2} \right\rangle + 2\left\langle {\widetilde{{E_ + }}\widetilde{{E_-^\ast }}} \right\rangle } \right)\approx \displaystyle{1 \over 4}\left\langle {{\vert {\widetilde{{E_ + }}} \vert }^2} \right\rangle } \cr } $$The far-right approximation can be seen from equation (37) using the fact that τD is the resolution of the series, such that for the full range of empirically accessible frequencies, we have ωτD < 1. In addition, because τ > τD, the factor |ωτ D/(1 + (iωτ)h)| < <1.
The Properties of the Model
Scaling Properties
High- and Low-Frequency Exponents
Given that both the model equations and the corresponding model statistics are both readily expressible in the spectral domain, it is tempting to empirically evaluate the model directly using spectra. Unfortunately, although spectra are commonly used, they suffer from various technical issues such as spectral leakage (the “smearing out” of spectral variance across a range of frequencies), effects that are important whenever strong spectral peaks are present. Problematic peaks occur not only in quasi-periodic signals (where most of the effort has been made [e.g., Springford et al. Reference Springford, Eadie and Thomson2020]) but also scaling signals, especially when the spectral exponent β (eq. 19) is significant (β = 0 for white noise [Hébert et al. Reference Hébert, Rehfeld and Laepple2021]). As we discuss in Appendix 3, the challenges for estimating spectra are much greater when the data are sampled at irregular intervals (as they typically are in paleobiology), and in particular, when the chronology itself (i.e., the temporal measurement density) is scaling (Appendix 3, Fig. A3.1). Indeed, the problem of how best to estimate the spectra of data with scaling chronologies (i.e., with holes over a wide range of scales) is still unsolved, and many nonstandard approaches give large biases.
Even if one can carefully handle the technical aspects, spectra still have the difficulty that their interpretation remains problematic. This was dramatically illustrated when Lovejoy (Reference Lovejoy2015) discovered that the iconic (and still frequently cited) spectrum of Mitchell (Reference Mitchell1976) was in error by a factor of 1015, an error that had not been noticed for four decades and was probably a consequence of spectra that are often not plotted with amplitude units at all or with undetected errors in the units. If Mitchell's spectrum had been accurate, two consecutive million year average Earth temperatures would only have differed by about 10 microKelvins (μK)—yet this patently false implication was not noticed because the units of the spectrum (K2yr) were not intuitive, whereas an RMS Haar fluctuation of 10 μK would have been obviously problematic. Even spectral updates as recent as 2020 are in error by a factor 1011 (see the review by Lovejoy [2023]). The comparison of the fluctuation analyses in this paper with those of the spectra (Appendix 3) highlights the power of fluctuation analysis when applied to irregularly sampled data.
Fortunately, over scaling ranges, the fluctuation and spectral statistics are both scaling with exponents as indicated in the previous section. Therefore, to interpret the statistics (eqs. 38, 39) in real space, it suffices to use the fact that Fourier scaling implies real-space scaling and to use the above relations between real-space and Fourier scaling exponents (eq. 22). In matrix form, the spectral exponents are therefore:
$$\beta _h = \left({\matrix{ {2( {\alpha + h} ) } & {2\alpha + h} & {2( {\alpha + h} ) } & {2( {\alpha + h} ) -1} \cr {2\alpha + h} & {2\alpha } & {2\alpha + h} & {2\alpha + h-1} \cr {2( {\alpha + h} ) } & {2\alpha + h} & {2( {\alpha + h} ) } & {2( {\alpha + h} ) -1} \cr {2( {\alpha + h} ) -1} & {2\alpha + h-1} & {2( {\alpha + h} ) -1} & {2( {\alpha + h-1} ) } \cr } } \right)$$
$$\beta _l = \left({\matrix{ {2( {\alpha + h} ) } & {2\alpha + h} & {2\alpha + h} & {2\alpha + h-1} \cr {2\alpha + h} & {2\alpha } & {2\alpha } & {2\alpha -1} \cr {2\alpha + h} & {2\alpha } & {2\alpha } & {2\alpha -1} \cr {2\alpha + h-1} & {2\alpha -1} & {2\alpha -1} & {2( {\alpha -1} ) } \cr } } \right)$$(The elements correspond to T, E +, D, and E −, left to right, top to bottom). Using the relationship between H and β (eq. 22), the high and low frequency (here small and large times, t) have exponents:
$$H_h = \left({\matrix{ {\alpha + h-\displaystyle{1 \over 2}} & {\alpha + \displaystyle{{h-1} \over 2}} & {\alpha + h-\displaystyle{1 \over 2}} & {\alpha + h-1} \cr {\alpha + \displaystyle{{h-1} \over 2}} & {\alpha -\displaystyle{1 \over 2}} & {\alpha + \displaystyle{{h-1} \over 2}} & {\alpha + \displaystyle{h \over 2}-1} \cr {\alpha + h-\displaystyle{1 \over 2}} & {\alpha + \displaystyle{{h-1} \over 2}} & {\alpha + h-\displaystyle{1 \over 2}} & {\alpha + h-1} \cr {\alpha + h-1} & {\alpha + \displaystyle{h \over 2}-1} & {\alpha + h-1} & {\alpha + h-\displaystyle{3 \over 2}} \cr } } \right)$$While at low frequencies, large Δt (i.e., large lags) we have:
$$H_l = \left({\matrix{ {\alpha + h-\displaystyle{1 \over 2}} & {\alpha + \displaystyle{{h-1} \over 2}} & {\alpha + \displaystyle{{h-1} \over 2}} & {\alpha + \displaystyle{h \over 2}-1} \cr {\alpha + \displaystyle{{h-1} \over 2}} & {\alpha -\displaystyle{1 \over 2}} & {\alpha -\displaystyle{1 \over 2}} & {\alpha -1} \cr {\alpha + \displaystyle{{h-1} \over 2}} & {\alpha -\displaystyle{1 \over 2}} & {\alpha -\displaystyle{1 \over 2}} & {\alpha -1} \cr {\alpha + \displaystyle{h \over 2}-1} & {\alpha -1} & {\alpha -1} & {\alpha -\displaystyle{3 \over 2}} \cr } } \right)$$We should add here that because E, O are linear combinations of E +, E −, their exponents will the maximum of those of E +, E −, so that:
$$\matrix{ {H_{h, E_ \pm } = \max \left({\alpha -\displaystyle{1 \over 2}, \;\alpha + h-\displaystyle{3 \over 2}} \right) = \alpha -\displaystyle{1 \over 2}; \;} & {h < 1} \cr {H_{l, E_ \pm } = \max \left({\alpha -\displaystyle{1 \over 2}, \;\alpha -\displaystyle{3 \over 2}} \right) = \alpha -\displaystyle{1 \over 2}} & {} \cr } $$We see that for the physically relevant parameters, H = α − ½ = −0.25 for both E, O, over the whole range (close to the data, see SL and Fig. 4).
Figure 4. This shows the Phanerozoic marine animal macroevolutionary analysis of the six series discussed in this paper; D, T, O, E are replotted from SL. The dashed lines show the theory slopes (eq. 44) with transition at Δt ≈ 40 Myr, i.e., log10Δt ≈ 1.6.
To get a concrete idea of the implications of model, let us use the rough empirical estimates from SL of α = 0.25, h = 0.5. Plugging these values into equations (41) and (42), we obtain:
$$ \eqalign{H_h &= \left({\matrix{ {0.25} & 0 & {0.25} & {-0.25} \cr 0 & {-0.25} & 0 & {-0.5} \cr {0.25} & 0 & {0.25} & {-0.25} \cr {-0.25} & {-0.5} & {-0.25} & {-0.75} \cr } } \right)\cr H_l &= \left({\matrix{ {0.25} & 0 & 0 & {-0.5} \cr 0 & {-0.25} & {-0.25} & {-0.75} \cr 0 & {-0.25} & {-0.25} & {-0.75} \cr {-0.5} & {-0.75} & {-0.75} & {-1.25} \cr } } \right)}$$(again, for T, E+, D, and E −, left to right, top to bottom). We can see that the Haar fluctuations will be useful for all the series over the whole range of frequencies/scales, the only exception being ΔD(Δt) at long lags (Hl < −1, lower right corner of the Hl matrix with Hl < −1). In this case, the Haar fluctuations “saturate,” and the spurious (limiting) value Hl = −1 is obtained.
With these results, we can comment on the issue of the stationarity/nonstationarity of the statistics. In the real world, scaling regimes only exist over finite ranges of scale, they have finite inner and outer limits. However, the outer scaling limits in mathematical series may be infinite. In this case, scaling series with low frequency Hl < 0 are stationary, whereas if Hl > 0, they will be nonstationary (although even here, if Hl < 1, the increments of the series will be stationary). From equation (44), we see that the only correlation with Hl > 0 is for the temperature; hence all the evolutionary responses (including cross-correlations) are stationary. However, even the temperature will be stationary at timescales beyond the outer limit of the scaling regime (which we argue is at least of the order of the length of the Phanerozoic). That is why we prefer the term “wandering” when Hl > 0: for such scaling processes, the nonstationarity may be spurious, the series only appears to be nonstationary over the observed (finite) range of scales. On a more fundamental note, no empirical series is stationary or nonstationary, the latter are properties of models, not data.
Normalized Correlations
The cross-spectra and cross-covariances (eq. 38) can be used to determine the normalized correlations that were estimated in SL:
$$\rho _{\,jk}( {\Delta t} ) = \displaystyle{{\left\langle {\Delta S_j( {\Delta t} ) \Delta S_k( {\Delta t} ) } \right\rangle } \over {{\left\langle {\Delta S_j{( {\Delta t} ) }^2} \right\rangle }^{1/2}{\left\langle {\Delta S_k{( {\Delta t} ) }^2} \right\rangle }^{1/2}}}$$(using Haar fluctuations). However, from equations (41) and (42), we find that their exponents (whether at high or low frequencies) are 2H jk − (H jj + H kk) = 0, that is, they are not power laws and only vary at sub–power law rates, and they are therefore nontrivial (i.e., they are significant) over the whole range of Δt. Because there are six series (T, E, D, O, E +, E −) there are 15 pairs whose fluctuation correlations may be determined over the observed range of 3 ≈ < Δt ≈ < 400 Myr; see Figure 5. The key correlations are those that correspond to the model parameters: ρ E = ρ TE, ρ D = ρ TD, see equation (32). We can already see that the correlations are quite noisy, a consequence of the low resolution and variable sampling of the series. To make a proper model–data comparison, we therefore turn to numerical simulations.
Figure 5. The (normalized) pairwise correlations of the 15 pairs of the six series as functions of lag. Several of these are reproduced from SL.
Numerical Simulations
The Statistics of the Simulated Series
The model has two fundamental exponents (α, h), two basic correlations (ρE = ρTE+, ρD = ρTD), and a crossover timescale τ. The third correlation (ρDE) is a derived parameter (eq. 32). In addition, there are two amplitude factors, σT, σE, but these will depend on the nondimensionalization/normalization of the series; on log-log plots, they correspond to an up–down shift, and on (normalized) correlation plots, the normalization eliminates them; they will not be considered further.
We used the results of SL to fix the values α = 0.25, h = 0.5, τ = 32 Myr. (This is the nearest power of 2 to the slightly larger—but only roughly estimated—value τ = 40 Myr in SL. Also, due to the scaling of the model, the more significant parameter is log τ and log2 40 = 5.3, not far from log2 32 = 5.) This leaves the only unknown parameters as the TE and TD correlations (ρE = ρTE+, ρD = ρTD; Fig. 5).
Before comparing the model directly to the (noisy) data, we checked that we were able to numerically reproduce the theoretically expected behavior. The basic modeling technique is to use convolutions with various (impulse response) Green's functions; this is detailed in Appendix 1, but follows the methods described in Lovejoy (2022). The main numerical problems are the small scales that have singular power law filters that are not trivial to discretize, and there are some (easier to handle) long-time (low-frequency) issues.
Rather than attempting to rigorously determine optimum parameters, as indicated earlier, we fixed the exponents α = 0.25, h = 0.5, and the crossover scale τ = 32 Myr. With guidance of the Figure 5 correlations for ρTE +, ρTD and some numerical experimentation, we took ρE = 0.5, ρD = −0.1 (hence ρTE + = 0.5, ρTD = −0.1, ρDE + = −0.9; i.e., the sign of r was taken as negative, eq. 32). We then performed simulations at a resolution of 250 kyr resolution, with simulation length of ≈4 Gyr (214 = 16,384 points), shown in Figure 6. This single very large realization has statistics that are close to those of an infinite ensemble. We postpone a discussion of the significance of the correlations to the section “The Statistics of the Simulated Series Resampled at the Data Sampling Times.”
Figure 6. The previous 214 simulation degraded from ¼ Myr resolution to 1 Myr. Curves normalized by their standard deviations and then offset by 5 units in the vertical for clarity.
According to the model (see the diagonal elements in eq. 44), the only series with positive low-frequency scaling exponent (Hl > 0) is the temperature (Hl = 0.25), which indeed shows “wandering” behavior (second from the bottom in Fig. 6); from the figure, one can see its long-range correlations as low-frequency undulations. This is also true for D, but only up to the crossover scale (≈32 Myr), after which consecutive 32 Myr intervals tend to cancel (Hl < 0, eq. 44). The other series are, on the contrary, “canceling” (Hl < 0, Hh < 0) especially E − (eq. 44). We can also visually make out some of the correlations, but this is clearer at lower resolution, as discussed later.
On these simulations, we can check that the theoretical scaling is obeyed; this was done using Haar fluctuations; see Figure 7 where the theory slopes (from eqs. 43, 44) are shown as reference lines. Note that because the Haar analysis “saturates” at H = −1, the low-frequency Hl = −1.25 value for E − (eq. 44, lower right-hand diagonal element) yields a slope of −1 (not −1.25); however, the other slopes are accurately estimated. Note that the theory–simulation agreement is not perfect, mostly because the theory is for the average statistics over an infinite ensemble, whereas Figure 7 is from a single—albeit large—simulation.
Figure 7. Simulation 214 =16,384 points with theoretical slopes indicated. The transition scale τ is 27 = 128 units, indicated by dashed vertical lines. The could represent a series modeled at 250 kyr resolution with a total simulation length of 4 Gyr and with crossover at 32 Myr. Due to its length, this simulation has statistics that are close to the ensemble averaged statistics. The parameters are: α = 0.25, h = 0.5, ρE = ρTE = 0.5, ρD = ρTD = −0.1 (with derived DE correlation ρDE = −0.9).
We can also work out the 15 correlations as functions of lag (Fig. 8). The figure shows the model parameters ρTE + = 0.5 (= ρE = 0.5), ρTD = –0.1 (= ρD = –0.1) as solid black reference lines, and the derived correlation ρDE + = –0.9 (eq. 32) as a dashed reference lines. Also shown are dashed theory lines for the TE, TO correlations (predicted to be equal to equal to TE + at long lags, eq. 39) and the DE, DO correlations (predicted to be equal to DE + at long lags, see eq. 39). We can see that the correlations approach the theoretical correlations at large lags, although the results are somewhat noisy.
Figure 8. The 15 pairwise correlations from the 214 realization in Fig. 7. Only two of the correlations were prescribed, and this only at a single resolution; the rest are consequences of the model, the two exponents α, h, and the crossover time τ = 27 (shown as short dashed vertical lines). The two prescribed correlations (DT, TE +) are shown as solid horizontal lines, and the derived correlations are shown as dashed lines (DE + from DT, TE +, eq. 32) and then TE, TO (predicted to be equal to equal to TE + at long lags, eq. 39) and DE, DO (predicted to be equal to DE +, at long lags, eq. 39). Note that these are from a single realization of the process, not the ensemble average. In addition, the statistics of some are fairly sensitive to irregularly sampled (and small size) of the empirical data; compare with Fig. 11.
The Statistics of the Simulated Series Resampled at the Data Sampling Times
Before making more effort at parameter fitting and comparing the model to data, it is important to take into account the small number of empirical data points and their irregular sampling (the corresponding geochronology itself turns out to be scaling; see Appendix 3, Fig. A3.1). Figure 9 shows the result for a simulation with the same parameters, but with a 1 Myr temporal resolution (right-hand side), resampled at the same times as the data (left-hand side). Because the model and data are only expected to have similar statistics, the detailed “bumps” and “wiggles” are unimportant, but one can nevertheless make out realistic-looking variability, including correlations between the series. Note that the model respects causality, so when there is a large extinction event, the curve is asymmetric with a rapid upturn being followed by a slower downturn (however, we have followed convention such that the present is at the left and the past at the right).
Figure 9. Model–simulation comparison with all series normalized by their standard deviations. The simulation was at 1 Myr resolution, and it was sampled at the same (irregular) times as the data (84 points over the last 500 Myr). Each curve was displaced by 5 units in the vertical for clarity. Due to causality, the series are asymmetric, with time running from right to left. The simulation is on the right.
We can now consider the fluctuation scaling and correlation statistics on the resampled series and compare them with both the data and the results from the same simulations but at a regular 1 Myr resolution (Fig. 10; see Fig. A2.1 for the E, O plots for the alternative Sepkoski stages database; Appendix 2, the E, O scaling and correlations are nearly identical). The figure shows a log-log plot of the RMS fluctuations as a function of the lag. To make the comparison, they were normalized by their standard deviations, but this is somewhat arbitrary, such that the up–down displacement (corresponding to a different nondimensionalization/normalization) is unimportant. To judge the realism of the model, the appropriate comparison is between the shapes of the resampled model output (red) and the data (black). We can see that the two are fairly close, although both model and data are noisy due to the small number of points and the irregular sampling. The agreement must be assessed not only by allowing for (relative) vertical shifts, but also by noting that the scales on the top D, T comparisons are such that the fluctuations vary only over a small range (for the data, factors of ≈1.7 for D and ≈2 for T) for lags varying over range of about a factor 100. In comparison, the E+, E −, E, O ranges are closer to factors of 10. Aside from this, these basic fluctuation statistics are fairly close to the data.
Figure 10. A comparison of the RMS Haar fluctuations for the 1 Myr resolution simulations discussed earlier (brown), from the simulation resampled at the data times (red), and from the data (black), these two irregularly sampled series are shown in Fig. 9. The relative vertical offsets of the curves are not significant; they correspond to specific normalizations/nondimensionalizations. We see that in general, the resampling at the data times (red) yields a closer fit to the data (black) than the analysis of the simulation at uniform (1 Myr) intervals; this is especially true for E −, O, E, E +. FMEM, Fractional MacroEvolution Model.
Figure 11. The pairwise correlations from the same three series as in fig. 10 with the same color codings: i.e., data (black), brown the simulation at a uniform 1 Myr resolution, and (red), the simulation resampled at the data times. The resampling notably improves the correlations for DE+, DO, DE, E+E −, E −O, OE, and to a lesser extent the OE, TE + comparisons; for the others, it is about the same. FMEM, Fractional MacroEvolution Model. The solid and dashed horizontal lines are the same as those in figure 8.
Figure 10 also gives important information about the effect of the sampling: compare the resampled (red) and uniformly sampled analyses (brown). The resampling is particularly important for E+, E −, E, O, although the effects are mostly at small lags for E+, E, O but at large lags for E −. This information should prove useful in interpreting a variety of real-world extinction and origination data.
Finally, we can compare the 15 pairwise correlations (Fig. 11). Again, to judge the realism of the model, compare the red and black correlations. Although—as expected—these are fairly noisy, we see that the agreement is quite good, significantly, it is generally much better than the agreement between the uniformly sampled correlations (brown curves) and data (black). By comparing the red (resampled) and brown (uniformly sampled) correlations, we see that the resampling is especially important for the DE+, DO, DE, E+E −, E −O, OE, correlations and to a lesser extent, the OE, TE + comparisons; for the others, it is about the same. We could note the successful prediction that the E+E, E+O, OE correlations should be ≈1 and the E −E correlations should be ≈ −1. Interestingly, the prediction that the E −O correlations should be ≈ −1 (eq. 39) is verified with the uniform sampling (i.e., it is indeed a property of the model), yet the resampling (red in the lower left graph in Fig. 11) makes it > 0 and aligns it closely with the observations. In other words, when the pure model predictions are poor (brown vs. black), there are many instances where the effects of nonuniform sampling are particularly strong, such that overall the model explains the data fairly well: overall 6 fluctuation plots (Fig. 10) and 15 correlations (Fig. 11) with 5 adjustable parameters (α, h, τ, ρE, ρT).
Discussion of the Model and Physical Significance of the Correlations
The model was motivated by an attempt to model the diversity process as a scaling crossover phenomenon with wandering climate (paleotemperature) and stabilizing life (turnover) scaling drivers. In the course of the model development, it became clear both theoretically (due to the definition of the diversity, eq. 5) and empirically that rather than E, O being fundamental, it was rather the turnover E + that is fundamental (indeed, the E + and E − statistics are quite different (Figs. 4, 10), and the E+E − correlations are nearly zero (Figs. 5, 11). In any event, the model predicted that E, O would follow the E + statistics (eq. 39; Figs. 4, 10, and the E+E and E+O correlations in Figs. 4, 11).
A more counterintuitive finding concerns the correlations. To start with, the model specifies that the diversity is primarily driven by the temperature up until the crossover scale, yet the temperature and diversity are negatively correlated over the entire range! Although at any given time lag, the DT correlation is small (−0.1), it means that there is a (weak) tendency for the diversity fluctuations to decrease when temperature fluctuations increase and vice versa, but this is not enough to offset the overall temperature control of the diversity that implies that consecutive temperature fluctuations tend to add up (HT = 0.25 > 0), and this is a stronger overall effect.
There is an additional more subtle effect. Consider that at each scale, the imposed TE + correlation is moderate and positive (ρTE+ = 0.5), and together, ρTD (the temperature diversity correlation) and ρTE+ (the temperature turnover correlation with r < 0, eq. 32) imply that at each lag, DE + is negatively correlated (reaching the theory value ρDE+ ≈ −0.9, at long lags; see the DE + correlation, the brown curve in Fig. 11). As the turnover E + also drives the diversity (eq. 1), at each scale, we thus have a tendency for T and E + fluctuations to increase (or decrease) together but for D and E + (and hence T and D) to have opposite tendencies. The overall result is that the weak anticorrelation of D with T and D with E + at any fixed scale is still dominated by the stronger effect of T fluctuations growing with scale and dominating the E + driver at lags < τ.
We could remark that ρTE+ = 0.5 > 0 indicates a tendency for temperature changes to “stimulate” the turnover: periods of increasing temperatures tending to be associated with increasing turnovers, and periods of decreasing temperatures with decreasing turnovers. Also there is a strong anticorrelation between D and E + (ρDE+ ≈ −0.9) that indicates that turnover decreases with diversity (note that this anticorrelation seems to nearly disappear after the nonuniform sampling; see Fig. 11, second in the top row). However over the range of scales that E + dominates dynamics of D (i.e., Δt > τ, as HE + ≈ −0.25 < 0), successive E + fluctuations tend to cancel, and on long timescales, the latter effect is dominant, such that HD = HE + ≈ −0.25—this is a scaling region of biotic self-regulation.
Discussion and Conclusions
The driver of macroevolutionary biodiversity has famously been reduced to a dichotomy between life and the environment: the metaphor of Red Queen versus Court Jester (Van Valen Reference Van Valen1973; Barnosky Reference Barnosky2001). Using genus-level time series from the PBDB, Spiridonov and Lovejoy (Reference Spiridonov and Lovejoy2022; SL) systematically analyzed fluctuations in extinction (E) and origination (O) rates, biodiversity (D), and paleotemperatures (T) over the Phanerozoic. They did this as a function of timescale from the shortest (≈3 Myr) to longest lags available (≈400 Myr), and their analysis included the correlations of the fluctuations at each scale. They concluded that T, E, O—the basic climate and life parameters—showed evidence of wide range scaling, supporting the hypothesis that over this range, there is a single biogeological “megaclimate” (Lovejoy Reference Lovejoy2015) regime with no fundamental timescale. However, they found that D followed the T fluctuations up until a critical time τ ≈ 40 Myr, whereas at longer timescales, it followed life (E, O): D was a scaling crossover phenomenon. At the shorter timescales Δt < τ, following the temperature, the D scaling exponent HD ≈ +0.25 (i.e., > 0), indicating that fluctuations tended to grow with scale, leading to “wandering” behavior. In contrast, for time lags Δt > τ, following E, O, its scaling exponent was HD ≈ −0.25 (i.e., < 0), hence successive fluctuations tended to cancel, resulting in long-time stabilization of diversity by life.
The model assumes, first, that in the Phanerozoic, dynamical processes occur over a wide range of timescales, and second, that the basic driving processes (here T, E +) are also scaling, such that over a wide range, they define no characteristic timescale. This is the simplest assumption, yet it is compatible with rich behavior. The critical scale for the diversity (estimated at ≈40 Myr) is simply the scale at which one process (life/turnover) starts to dominate the other (climate/temperature).
The Phanerozoic evolution of life on our planet is full of contrasting patterns, quantitative and qualitative dynamical transitions, yet it may be compatible with an underlying scale symmetry. For example, transitions related to era boundaries include changes in complexity of marine animal communities in post-Paleozoic times (Wagner et al. Reference Wagner, Kosnik and Lidgard2006), the increase in longevities of animal genera in the same post-Paleozoic period (Miller and Foote Reference Miller and Foote2003), or changes in the structure of bivalve communities after the K-Pg mass extinction (Aberhan and Kiessling Reference Aberhan and Kiessling2015), yet each of these processes occur over wide ranges of scale.
In addition, it could be argued that the dynamics of biodiversity should be subdivided into these smaller (substage) time intervals, potentially extending the scaling range to shorter timescales (≈1 Myr or less). Important boundaries that modified macroevolutionary dynamics at global scales can also be found at period and subperiod boundaries: apparently the appearance of calcifying plankton at the beginning of the Jurassic (or so called “post-Conodontozoic”; Ferretti et al. Reference Ferretti, Bancroft and Repetski2020) should have had a significant effect in stabilizing hyperthermal events and thus decreasing the volatility of biospheric evolution (Wignall Reference Wignall2015; Eichenseer et al. Reference Eichenseer, Balthasar, Smart, Stander, Haaga and Kiessling2019). In addition, presumably the appearance, global spread, and deep impact on the climate of terrestrial plants in the Silurian is highly significant (Lenton et al. Reference Lenton, Dahl, Daines, Mills, Ozaki, Saltzman and Porada2016). There are also shorter timescale transitions—such as transitions toward more volatile dynamics of zooplankton between the Early and Middle Ordovician toward the Late Ordovician (Crampton et al. Reference Crampton, Cooper, Sadler and Foote2016).
Such changes are thus numerous and add to the plausibility that such diverse and sometimes sharply varying processes may be scaling. Therefore, the subdivision of the Phanerozoic into shorter time intervals in search of local scaling laws—in our view at least at this stage of the understanding and data availability—is rather premature and arbitrary and could preclude the search of overarching commonalities that can hide in this apparent heterogeneity. Moreover, scaling multifractal processes are more general than the Gaussian ones considered here; they are intermittent, typically involving occasional large changes in intensity (e.g., apparent transitions), occurrences of extremes. Physically, this could correspond to local dominance of controlling mechanisms (which can be hierarchically modulated by feedbacks from biota or external forcing) at all timescales where scaling regimes apply.
In this study, as well as in the study which first described scaling properties of Phanerozoic marine genus-level evolution, SL analyzed the processes and time series as a homogenous stationary system. Although this is an assumption that the statistical properties are constant over the eon, it is nevertheless consistent with high variability and sudden transitions (i.e., volatile and intermittent behavior). Although the exact FMEM model presented here was not intermittent, intermittency and its associated sharp transitions could easily be introduced simply by replacing the (Gaussian) white noise forcings by multifractal ones.
To clarify our ideas, to better understand the geobiodynamics, and to better understand and quantify the limitations, biases, and other data issues, we proposed the simple FMEM to reproduce the observations. It is a model of macroevolutionary biodiversity driven by paleotemperature (the climate proxy) and the turnover rate (E + = O + E), the life proxy. To fit with basic empirical scaling statistics and theoretical ideas about the macroclimate regime (form timescales of roughly 1 Myr to at least 400 Myr), these drivers were taken to be scaling with climate dominating at short timescales and life at long timescales. Therefore, FMEM suggests a possible way to combine into a single stochastic framework both: (1) the destabilizing geophysical (and possibly astrophysical) processes (Raup Reference Raup1991, Reference Raup1992a,b; Melott and Bambach Reference Melott and Bambach2014; Fields et al. Reference Fields, Melott, Ellis, Ertel, Fry, Lieberman, Liu, Miller and Thomas2020) with (2) the stabilizing, density-dependent, and self-regulating biotic processes. The model is specified by a simple parametrization based on two scaling exponents and two pairwise correlations (between T and E + and between T and D).
Because the FMEM model is linear, it can be used to model a superposition of stochastic and deterministic processes such as bolide impacts: the responses to the deterministic (external) and stochastic (internal) forcings simply add. As the impulse response function is a (singular) power law, a short impulse on the system will generate a sharp initial response followed by a long power law decay (i.e., much slower than the classical exponentials), decaying with somewhat different exponents depending on whether the impulse is on the temperature or turnover or both (see Fig. 2).
The model has two unusual characteristics: first, it is stochastic, such that the crossover from climate to life dominance is thus a scaling (power law) not standard exponential (i.e., Markov process–type) transition. Stochastic models involve infinite dimensional probability spaces, they are therefore natural model types in systems with huge numbers of degrees of freedom. We believe that they are intrinsically more realistic than strongly nonlinear but deterministic chaos-type models (including those that are deterministic but are perturbed by noises). When the intermittency is strong, scaling stochastic models must be nonlinear (e.g., multifractal cascade processes), and this can easily be included in further model improvements—the Gaussian forcing (γ1, γ2, eq. 14) need only be replaced by a multifractal one. Here, intermittency is neglected, and linear stochastic equations with Gaussian white noise forcings are used (linear stochastic models can often be used even when the underlying dynamics are strongly nonlinear).
The other unusual FMEM characteristic is that it is a system of fractional differential equations. Unlike the familiar integer-ordered differential equations that typically have exponential impulse response functions (Green's functions), fractional equations typically have power law response functions and are natural ways to model scaling processes. These impulse response functions are physical models of bolide impacts and similar nearly instantaneous processes, and we discussed some implications.
The model is also highly parsimonious with two scaling exponents and a crossover time τ determined by the PBDB data as analyzed SL. These determined the basic scaling characteristics of the six series: T, E+, D, E − (= O − E), O, E. In addition, the model had two correlations that were specified: those between T and E + and between T and D. From these, the other 13 pairwise correlations (out the 15 possible pairs of the six series) were implicitly determined and were compared with the data.
The fractional derivatives were of the Weyl type such that their Fourier transforms are simply power laws. Because the system is ultimately forced by two Gaussian white noises, only the second-order statistics (i.e., the spectra and correlation functions) are needed, and these are easily obtained: the basic solutions are ffRns that were recently introduced (Lovejoy Reference Lovejoy2022a). In future, more realistic intermittent (multifractal) forcings could be used instead of the Gaussian white noise. Beyond exhibiting the full solution to the equations with a full statistical characterization, we then implemented the model numerically, first verifying the model against the theoretically predicted behavior. By producing simulations at 1 Myr resolution, we are able to resample the output at the same irregular sampling times as the PBDB. The statistical characteristics of the results (the six scaling curves showing the fluctuations as functions of timescale) plus the 15 pairwise correlations as functions of timescale are all quite close to the data, and in several cases, the agreement could be clearly attributed to the limitations, biases, and so on of the data. In particular, this was the case of the DE+, DO, DE, E+E −, E −O, OE correlations that are much closer to the data following the irregular sampling than with the original model outputs uniformly sampled at 1 Myr resolution.
Given the model's simplicity, it thus is remarkably realistic. This is fortunate, as until higher-resolution (global-scale) time series become available (e.g., Fan et al. Reference Fan, Shen, Erwin, Sadler, MacLeod, Cheng and Hou2020), more complex models may not be warranted. In any case, the model is able to help explain some subtle points about the interaction of different correlated series that are also strongly self-correlated over wide ranges of timescales, and this with quantitatively and qualitatively different scaling behaviors (“wandering” vs. “canceling”/self-stabilizing).
Acknowledgments
We would like to thank the editors L. H. Liow and M. J. Hopkins and the reviewers, C. Marshall and three anonymous reviewers, for many helpful suggestions that significantly strengthened and clarified the article. The research project was supported by grant S-MIP-21-9 “The Role of Spatial Structuring in Major Transitions in Macroevolution” and National Science and Engineering Research Council grant 222858 “Scaling and Multifractals in Geosystems.” This paper is Paleobiology Database official publication no. 470.
Competing Interests
The authors declare no competing interests.
Appendix 1. Numerical Simulations
Because the model is linear, the obvious simulation method is to use Fourier techniques. The main problem is that the small scales have singular power law filters that are not trivial to discretize, there may also be some long time (low-frequency) issues. A convenient way is to use techniques developed for simulating ffRn processes discussed in Lovejoy (2022). ffRn processes can be simulated by convolving Gaussian white noises with the ffRn Green's function G α,h (eqs. 9, 10). A somewhat better numerical technique is to use the step-response Green's function (= G α+1,h; it is the smoother—and hence easier to handle—integral of G α,h), followed by a numerical differentiation.
Appendix 2. Sepkoski's Genus Compendium Data and Macroevolutionary Metrics
2.1. The Rates
To explore the robustness of the results based on the PBDB data, we used the Sepkoski database of stratigraphic ranges for marine animal genera (Sepkoski Reference Sepkoski2002, as implemented in an online tool; Peters Reference Peters2013). We used the option of including all available fossilized Metazoan phyla. This yielded 19,713 resolved to substage-resolution stratigraphic-range records of the genera. Based on this record, the standing diversities and Foote's per-lineage per-substage extinction (q) and origination (p) rates (Foote Reference Foote2000; Peters and Foote Reference Peters and Foote2002) were calculated in the web implementation:
$$q = {-}\ln \left({\displaystyle{{N_{bt}} \over {N_{bl} + N_{bt}}}} \right)\eqno( A2.1) $$
$$p = {-}\ln \left({\displaystyle{{N_{bt}} \over {N_{\,ft} + N_{bt}}}} \right)\eqno( A2.2) $$Here, N bt is the number of taxa crossing bottom and the top of the interval, N bl is the number of taxa crossing bottom of the interval and disappearing in the interval, and N ft is the number of taxa first appearing in a given interval and crossing the upper stratigraphic boundary. Here and in the case of PBDB data, rates were not normalized for the durations of intervals, because most of the evolutionary turnover (as was shown in previous studies) is concentrated at or near the stratigraphic boundaries, which most often are defined by such events (Foote Reference Foote2000, Reference Foote2005; Payne and Heim Reference Payne and Heim2020).
Averages of substage boundary ages were used for assigning the diversity, extinction, and origination rate ages. The analyzed substage record consists of 154 points, which started at 535 and ended at 0.9 Myr. The average resolution of the record is 3.5 Myr, being higher than the standard PBDB binning. According to the stratigraphic-range nature of Sepkoski's compendium data (vs. occurrence-based records in the PBDB), no sampling standardization was performed.
Figure A2.1. Scaling of extinction (black) and origination (brown) rates estimated from Sepkoski's genus-level compendium data. Red line shows time-scaling of correlations between extinctions and originations (for the correlations, the values are not logarithmic: 0 represents no correlation, 1 represents maximum correlation). This analysis is similar to that used in SL but on the Sepkoski stages database (for a detailed comaprison, compare with the bottom left and right of Fig. 9). The bottom curves are the root-mean-square (RMS) origination and extinction rate fluctuations with reference line H = −0.25 (as for the Paleobiology Database [PBDB]). The correlations (red, use the same numerical scale as the fluctuations but are linear, varying up to nearly a maximum of 1 indicated by the horizontal dashed line) behave similarly to the PBDB correlations: they are low until about 40 Myr, after which they are high. Approximately at the crossover scale of 40 Myr, macroevolutionary rates functionally track each other, which results in negative scaling of diversity at the longest timescales—therefore showing the same pattern as sample standardized PBDB data for Phanerozoic marine animals and using second-for-third macroevolutionary rates (SL).
2.2. Scaling of Rates in Sepkoski's Genus Compendium Data
The analysis of Foote's rates calculated from the Sepkoski's marine animal compendium data shows an essentially identical pattern of scaling to the second-for-third rates calculated from PBDB data on marine animal genera (Fig. A2.1). Despite the different kinds of data—taxic ranges versus subsamples of occurrences—and differences in the methodologies of estimations in extinction and origination rates, the scaling exponent in both cases was identical H = –0.25, which shows that stabilizing behavior and the timescales of relaxation are essentially identical in both cases. Moreover, the transition to synchronization of extinction and origination rates (high correlations) is also achieved at long timescales after the crossover timescale of 40 Myr. Therefore, despite different extinction metrics sometimes showing quite distinct patterns (Foote Reference Foote1994, Reference Foote2000; Alroy Reference Alroy2015), the patterns of scaling and the statistical association between extinction and origination rates is robust either to the kind of data chosen and the metrics used. Additionally, the traditional taxic range based Sepkoski's data shows the same qualitative timescale-dependent pattern of transition between incoherent dynamics at shorter timescales (low correlations of E and O) toward stable dynamics (high correlations between E and O) at timescales > 40 Myr.
Appendix 3. Comparing Spectral and Fluctuation Analyses
The FMEM model as well as its corresponding statistical properties are most easily expressed in the Fourier domain; it is therefore natural and tempting to empirically evaluate it via spectral rather than fluctuation analysis. Indeed, had the PBDB been uniformly sampled, this would have been straightforward enough; standard discrete Fourier transforms could be used with standard windowing functions to reduce spectral leakage. This would have been quite adequate for determining scaling regimes and exponents. In certain cases, for example, where accurate determination of frequencies of spectral peaks are required, such as astrophysical applications—multitaper methods (MTM; Thomson Reference Thomson1982) could be used for further precision (as advocated, e.g., in Springford et al. [Reference Springford, Eadie and Thomson2020]).
But analyzing paleo-series poses unique problems for data analysis. Not only are the chronologies irregular, but there is increasing evidence that they are in fact themselves typically scaling! This means that the number of data points taken per time interval—the measurement density ρ(t)—is scaling, that is, the density fluctuations obey 〈Δρ(Δt)q〉 ∝ Δt ξ(q) (see Fig. A3.1 for the PBDB fluctuations). This is the chronological (temporal) counterpart of the spatial scaling of measuring stations reported in Lovejoy et al. (Reference Lovejoy, Schertzer and Ladoy1986); see also Lovejoy and Schertzer (Reference Lovejoy and Schertzer2013: chapter 3). Figure A3.1 also shows that the scaling of the measurement densities is non-Gaussian; this is because ξ(1) ≈ −0.15 but ξ(2) ≈ −0.35, such that ξ(2) ≠ ξ(1). A full characterization of scaling measurement densities will be discussed in another publication. For now, the key point is that its scaling implies potential biases in spectral exponents, and this is an area that has largely been neglected in recent work on spectral analysis.
To illustrate the difficulties of using spectral instead of fluctuation analysis, we used the Lomb-Scargle method (Lomb Reference Lomb1976) combined with a standard Hann window to reduce spectral leakage. Figure A3.2 shows the resulting spectra for the extinction and origination rates, and Figure A3.3 shows the diversity and temperature spectra. These can be compared with the RMS Haar fluctuations (black curves in Fig. 9 or the bottom curves in Fig. A2.1). The spectrum is sufficiently noisy that if the straight reference lines—whose slopes were inferred from the real-space fluctuation analyses (using β = 1 + 2H)—had not been added, the scaling would not have been obvious. Indeed, it frequently happens that in order to properly interpret paleo-spectra, performing an initial Haar fluctuation analysis is useful—even necessary—for their interpretation.
Numerical experiments show that in spite of the windowing, the low frequencies are still very sensitive to spectral leakage and so are strongly biased. For example, Hébert et al. (Reference Hébert, Rehfeld and Laepple2021) examined this issue numerically by determining the Lomb-Scargle spectrum of scaling processes and found that there were large spectral biases and that these increased as the spectral exponent β increased. This is unsurprising, because as β increases, the low frequencies that are most affected by leakage contain a larger and larger fraction of the total spectral variance. Yet Hébert et al. (Reference Hébert, Rehfeld and Laepple2021) underestimated the problem: the numerically simulated data gaps that they used were confined to a fairly narrow range of timescales: they did not consider the effect of scaling chronologies—that is, sampling irregulariites over a wide range of scales.
Although spectral leakage is a known—and serious—limitation of Lomb-Scargle spectra, Springford et al. (Reference Springford, Eadie and Thomson2020) argue that the problem can be alleviated using MTM (Thomson Reference Thomson1982). The problem with this proposal is that there is no theoretical basis for using the MTM on irregular data: the “Slepian” tapers (window functions) are no longer orthogonal, such that averaging the spectrum of many windowed tapers can make spectral biases worse, not better, a finding that we will elaborate in another publication. An alternative suggested by Hébert et al. (Reference Hébert, Herzschuh and Laepple2022) is to fill the gaps with linearly interpolated values and use a regular FFT combined with the MTM method on the resulting uniformly spaced series. Unfortunately, this interpolation method also leads to quite serious biases (especially at high frequencies). This is because the scaling exponent H is the maximum possible order of differentiation, and linear interpolations (i.e., of order 1) have H = 1, whereas the series in our case have H ≈ 0.25 (temperature) and H ≈ −0.25 (origination and extinction), such that the linear segments filling the gaps are far too smooth. The result is an uncontrolled mixing of data with different H values. A possible alternative might be to first estimate H using another method—for example, Haar fluctuations—and then to use fractal interpolation (e.g., Navascués et al. Reference Navascués, Pacurar and Drakopoulos2022) to do the interpolation using the correct exponent to fill in the gaps. Unfortunately, even this method may not be used if the measurement density is multifractal (i.e., the H exponent is only a partial characterization of a scaling process).
It should be stressed that real-space (fluctuation) analysis is advantageous, as the interpetation of the fluctuations is straightforward (recall the discussion of the missing quadrillion in the “The Drivers” section), and in addition, it is easy to handle missing data (irregular, including scaling chronologies) and also to determine correlations between series each having uneven chronologies—an advantage that we fully exploited here. In principle, these correlations could be studied by using cross-spectral analysis, but for single realizations (i.e., single series), this requires first breaking the series into segments, and this is more difficult when the chronologies are irregular.
Figure A3.1. The Haar fluctuations of the measurement density for the Paleobiology Database (PBDB) used here. The first two moments (mean, mean square) are shown with reference lines indicating the corresonding scaling exponents (logarithmic slopes). Over the range of ≈10–300 Myr, the chronologies themselves (i.e., ρ(t)) do not have well-defined resolutions; they are scaling.
Figure A3.2. Extinction (brown) and origination rates (red) obtained by the Lomb-Scargle method using a Hanning window. The frequencies higher than (4 Myr)−1 were not shown, as they are at higher frequencies than the mean resolution.
Figure A3.3. Diversity (red) and temperature (gray) spectra. The reference slopes are the theoretical slopes corresonding to H = −0.25 (spectral exponent β = 0.5), and H = 0.25 (spectral exponent β = 1.5).
Open access
