ASSESSING THE INFLUENCE OF THE SOLAR ORBIT ON TERRESTRIAL BIODIVERSITY

Over the course of Earth's history, evolution has produced a wide variety of life. This is particularly apparent from around 550 Myr ago—the start of the Phanerozoic eon—when hard-shelled animals first appeared and were preserved in the fossil record. Since then we observe a general increase in the diversity of life, but with significant variation superimposed (Sepkoski et al. 2002; Rohde & Muller 2005; Alroy et al. 2008). The largest and most rapid decreases in biodiversity—defined here as the number of genera extant at any one time—are referred to as mass extinctions.

The cause of these variations in biodiversity in general, and mass extinctions in particular, has been the subject of intense study and speculation for over a century. Many mechanisms have been proposed for the observed variation, which we can place into four groups.

First, the variations are the result of inter-species interactions. Species compete for limited resources, and as one species evolves to compete in this struggle for survival, so other species will evolve too. This idea has been referred to as the "Red Queen hypothesis" (Van Valen 1973; Benton 2009; in reference to the Red Queen's race in Lewis Carroll's Through the Looking Glass, where Alice must run just to keep still). Recent ecological studies indicate that the interaction between species can minimize competition and enhance biodiversity (George & Hao 2009), although it is not obvious that these biotic factors are the main cause of large-scale patterns of biodiversity (Benton 2009; Alroy 2008).

Second, the environment changes with time, and species will evolve in response to this. This ideas is sometimes called the "Court Jester hypothesis" (Barnosky 2001; Benton 2009). Some of these (abiotic) geological changes are relatively slow, such as plate tectonics, atmospheric composition, and global climate (Sigurdsson 1988; Crowley & North 1988; Wignall 2001; Martí & Ernst 2005; Feulner 2009; Wignall et al. 2009). Others may be more rapid. Large-scale volcanism, for example, would inject dust, sulfate aerosols, and carbon dioxide into the atmosphere, resulting in a short-term global cooling, reduced photosynthesis, long periods of acid rain, and resulting ultimately in a long-term global warming (on a timescale of 10⁵ yr; Martí & Ernst 2005).

Third, extraterrestrial mechanisms could be involved, either through a direct impact on life or by changing the terrestrial climate. Variations in Earth's orbit (mostly its eccentricity) over ten to one hundred thousand year timescales are responsible for the ice ages (Hays et al. 1976; Muller 2000). Extraterrestrial mechanisms on longer timescales could also play a role. These include solar variability (Shaviv 2003; Lockwood 2005; Lockwood & Frohlich 2007), asteroid or comet impacts (Shoemaker 1983; Alvarez et al. 1980; Glen 1994), cosmic rays (Shaviv 2005; Sloan & Wolfendale 2008), supernovae (SNe), and gamma-ray burst (GRBs; Ellis & Schramm 1995; Melott & Thomas 2008; Domainko et al. 2013; for a review see Bailer-Jones 2009). For example, cosmic rays might influence Earth's climate if they play a significant role in cloud formation (through the formation of cloud condensation nuclei; Carslaw et al. 2002; Kirkby 2007). Secondary muons resulting from cosmic rays—as well as high-energy gamma rays from SNe—could kill organisms directly or damage their DNA (Thorsett 1995; Scalo & Wheeler 2002; Atri & Melott 2012).

Finally, the apparent variation may be, in part, the result of uneven preservation and sampling bias (Raup 1972; Alroy 1996; Peters 2005). Fossilization is relatively rare and some animals are more likely to be preserved than others. Furthermore, the degree of preservation of marine species (more common in the fossil record) depends on the amount of continental outcrop available at any time, and this depends on the sea level (Hallam 1989; Holland 2012). The number of species or genera living at any one time is not observed but must be reconstructed from the times at which species appear and disappear, which implies some kind of sampling or modeling. This can introduce a bias, although it may be diminished to some degree by various techniques (Alroy et al. 2001; Alroy 2010).

These four types of cause of biodiversity variation are not mutually exclusive. They probably all acted at some point and will also have interacted. For example, an asteroid impact could release so much carbon dioxide that long-term global warming has the biggest impact on biodiversity. Alternatively, a cool period would lower sea levels, leaving less continental shelf for the preservation of marine fossils, even though biodiversity itself may be unchanged.

Given the limited geological record, untangling the relevance of these different causes in most individual cases is difficult, if not impossible. More promising might be an attempt to identify the overall, long-term significance of these potential causes. The goal of this article is to do that for extraterrestrial phenomena.

Many of the astronomical mechanisms mentioned above are ultimately caused by the presence of nearby stars. Stars turn SNe, the source of gamma rays, and their remnants are a major source of cosmic rays (Koyama et al. 1995). Stars perturb the Oort Cloud, the main source of comets in the inner solar system (Rampino & Stothers 1984; García-Sánchez et al. 2001). Broadly speaking, when the Sun is in regions of higher stellar density, it is more exposed to extraterrestrial mechanisms of biodiversity change. In its orbit around the Galaxy (once every 200–250 Myr or so), the Sun's environment changes. For example, it oscillates about the Galactic plane with a (quasi) period of 50–75 Myr (Bahcall & Bahcall 1985), and in doing so moves through regions of more intense star formation activity in the Galactic plane. This is particularly true if the Sun crosses spiral arms (Gies & Helsel 2005; Leitch & Vasisht 1998), which it may do every 100–200 Myr or so.

Such changes in the solar environment have been used as the basis for many claims of a causal connection between the solar motion and mass extinctions and/or climate change. Typically, authors have identified a periodicity in the fossil record and then connected this to a plausible periodicity in the solar motion (Alvarez & Muller 1984; Raup & Sepkoski 1984; Davis et al. 1984; Muller 1988; Shaviv 2003; Rohde & Muller 2005; Melott & Bambach 2010; Melott et al. 2012). These comparisons are fraught with problems, however, some of which remained unmentioned by the authors. The first is the fact that the solar motion and past environment are poorly constrained by the astronomical data, so a wide range of plausible periods are permissible (Overholt et al. 2009; Mishurov & Acharova 2011), yet the coincident one is naturally chosen. Second comes the fact that the solar motion is not strictly periodic even under the best assumptions. Third, many of these studies have not performed a careful model comparison. Typically, they identify the best-fitting period assuming the periodic model to be true, but fail to accept that a non-periodic model might explain the data even better (Kitchell & Pena 1984; Stigler & Wagner 1987). In some cases, a significance test is introduced to exclude a specific noise model, but this is often misinterpreted, and the resulting significance overestimated. The reader is referred to Bailer-Jones (2009) for an in-depth review and references.

Here we attempt a more systematic assessment of the possible role of the solar orbit in modulating extraterrestrial extinction mechanisms. Our approach is new in a number of respects because we (1) do a numerical reconstruction of the solar orbit (rather than just assuming it to be periodic); (2) take into account the observational uncertainties in that reconstruction; (3) use models which predict the variation of the probability of extinction with time (rather than assuming that extinction events occur deterministically, for instance); (4) do proper model comparison (rather than using p-values in an oversimplified significance test); (5) compare not only the orbital model with the fossil record but also numerous reference (non-orbital) models, such as periodic, quasi-periodic, trend, periodic with constant background, etc.

Our method is as follows. Adopting a model for the distribution of mass in the Galaxy, we reconstruct the solar orbit over the past 550 Myr by integrating the Sun's trajectory back in time from the current phase space coordinates (position and velocity). This gives us a time series of how the stellar density in the vicinity of the Sun has varied over the Phanerozoic. We then assume that this density is approximately proportional to the terrestrial extinction probability (per unit time). That is, we adopt a non-specific kill mechanism linking the solar motion to terrestrial biodiveristy. This is naturally a strong and rather simple assumption, but it should be emphasized that we are interested in the overall plausibility of extraterrestrial phenomena rather than trying to identify a specific cause of individual extinction events. The resulting time series is then compared to several different reconstructions of the biodiversity record.

A significant source of uncertainty in the reconstructed solar orbit is the current phase space coordinates (or "initial" conditions). We therefore sample over these to build up a set of (thousands of) possible solar orbits and compare each of these with the data. The comparison is done by calculating the likelihood of the data for each orbit. Rather than finding the single most likely orbit, we calculate the average likelihood over all orbits. This is important because it properly takes into account the uncertainties (whereas selecting the single most likely orbit would ignore them entirely). Indeed, these initial conditions can be considered as the six parameters of this orbital model (for a fixed Galactic mass distribution). We are, therefore, averaging the likelihood for this model over the prior plausibility of each of its parameters. This average—or marginal—likelihood is often called the "evidence." This is just the standard, Bayesian approach to model assessment, which avoids the various flaws of hypothesis testing (Jeffreys 1961; Winkler 1972; MacKay 2003), but unfortunately it has seen little use in this field of research.

The next step is to compare this evidence with that calculated for various reference models (sometimes called "noise" or "background" models, depending on the context). One such model is a purely sinusoidal model, parameterized by an amplitude, period, and phase. We generate a large number of realizations of the model for different combinations of the parameters, calculate the likelihood of the data for each, and average the results. This averaging plays the crucial role of accommodating the complexity of the model. A complex model with lots of parameters can often be made to fit an arbitrary data set well. That is, it will give a high maximum likelihood. But this does not make it a good model, precisely because we know that it could have been made to fit any data set well! Such models are highly tuned, so while the maximum likelihood fit may be very good, a small perturbation of the parameters results in poor predictions. Unless supported by the data very well, such models are less plausible. A simpler model, in contrast, may not give such an optimal fit, but it is typically more robust to small perturbations of the model parameters or the data, so it gives a good fit over a wider portion of the parameter space. The model evidence embodies and quantifies this tradeoff, which is why it—rather than the maximum likelihood—should be used to compare models.

We have selected four data sets for our study. The first two are compilations of the variation of extinction rate over time, from Rohde & Muller (2005) and Alroy et al. (2008). In the latter, we use the extinction rate standardized to remove the sampling bias. Both report a magnitude as a function of time. The second two data sets just record the time of mass extinction events. Here we take the times of the "big 5" mass extinctions and 18 mass extinctions identified by Bambach (2006) based on Sepkoski's earlier work (Sepkoski & Raup 1986). Each mass extinction is represented as a (normalized) Gaussian on the time axis, the mean representing the best estimate of the date of the event and the standard deviation the uncertainty. We refer to these as "discrete" data sets, as they just list the discrete dates at which events occur (we do not use any magnitude information). The two rate data sets we therefore refer to as "continuous" (even though in practice the rates are also recorded at discrete time points).

The paper is arranged as follows. We first introduce the data sets. In Section 3, we explain the modeling approach, how we calculate the likelihoods and evidence, and how we compare the models. In Section 4, we describe how we reconstructed the solar orbit, and also quantify the degree of periodicity typically present (as a strict periodicity has often been assumed in the past). In Section 5, we calculate and compare the evidences for the various models and data sets and test the sensitivity of the results to the model parameters and uncertainties in the data. We conclude in Section 6.

Sepkoski and others have identified five extinction events to be "mass extinctions" (Sepkoski & Raup 1986), often referred to as the "big five." Other studies have identified different candidates for these, or have identified a "big N" for some other value of N. For example, Bambach et al. identify three mass extinction events as being globally distinct (Bambach et al. 2004). Here we adopt a set of 18 mass extinction events (or B18) selected by Bambach (2006) using an updated Sepkoski genus-level database. They are consistently identifiable in different biodiversity data sets and when using different tabulation methods. The second of our discrete data sets is the "big five" as identified from among the B18. Other choices of events are of course possible and our results will, in general, depend on this choice (although, as we will see, the results are rather consistent).

The times and durations of the events are listed in Table 1. The time, τ, is the midpoint between the start age and the end age of the substage in which the extinction occurred, and the substage duration, d, is the difference between these. The geological record does not resolve the extinction event, so the extinction presumably took place more rapidly than this substage duration. In that case, τ is our best estimate of the true (but unknown) time, t, at which the extinction occurred, and d is a measure of our uncertainty in this estimate. Uncertainty is represented by probability, so we interpret an "event" as the probability distribution P(τ|t, d). This is the probability that we would measure the event time as τ, given t and d. We could represent this as a rectangular ("top-hat") distribution of mean t and width d, but this assigns exactly zero probability outside the substage duration, which implies certainty of the start and end ages. Even though their (relative) ages have uncertainties which are less than the event duration, we nonetheless accommodate some uncertainty in these start and end ages by considering each event to be a Gaussian distribution with mean t, and standard deviation σ equal to the standard deviation of the rectangular distribution, which is $\sigma = d/\sqrt{12}$. This Gaussian distribution is broader than the corresponding rectangular distribution. Note that the intensity of the mass extinction is not taken into account. A Gaussian is normalized, so its peak value is determined by its standard deviation (Figure 1).¹

Zoom In Zoom Out Reset image size

The discrete data sets are naturally biased in that they only select periods of high extinction rate. It may well be that extraterrestrial phenomena are only relevant in causing (or contributing to) mass extinctions, but a priori it is natural to ask how the overall extinction rate varies. The extinction rate, E(t), is the fraction of genera which go extinct in a stratigraphic substage divided by its duration. This is directly proportional to the variation of extinction probability per unit time. For one of our extinction rate data sets we use the linearized and interpolated data set constructed by Rohde & Muller (2005) as reported in Bambach (2006). We denote this as RM. The other data set is the "three-timer" extinction rate from the Paleobiology database² (Alroy et al. 2008). The data are binned into 48 intervals averaging 11 Myr in duration. The counts are derived from 281,491 occurrences of 18,541 genera within 42,627 fossil collections. We use the data set processed using their subsampling method in order to reduce the sampling bias and denote this A08. Both data sets are reported as lists of extinction rates at specific times, {r_j, τ_j}. These two continuous data sets are plotted in the lower row of Figure 1.

The goal of this work is to compare how well various models predict the paleontological data sets. We do this by calculating, for a given data set, the Bayesian evidence for each model. If the models are equally probable a priori, then the one with the highest evidence is the best predictor of the data. This does not exclude the possibility that there exists a better model which we have not yet tested, but it at least allows us to conclude that the lower evidence models are neither appropriate nor sufficient explanations of the phenomenon. (Indeed, we never assume a model is "true," just that it is better than the alternatives.)

The modeling approach is described in full by Bailer-Jones (2011a, 2011b), so it will only be outlined here. Let D denote the paleontological time series and M the model. Examples of M are the orbital model or the periodic model. Each model has a set of parameters, θ. The likelihood of the model, P(D|θ, M), is the probability of obtaining D from model M with its parameters set to some specific values of θ. Normally we do not know the exact values of these parameters, and the data—being noisy and imperfectly fit by the model—do not determine them exactly either. We therefore average the likelihood over all possible values of θ, weighting each by how plausible that value of θ is. This weighted average is the evidence. This weighting is given by the prior probability distribution, P(θ|M). In the case of the orbital model, where θ is the current phase space coordinates of the Sun, P(θ|M) is determined by the observational uncertainties in these measurements. Mathematically, the evidence is

$$\begin{equation} P(D|M)=\int _\theta P(D|\theta,M) P(\theta |M) d\theta. \end{equation} \tag{ 1 }$$

It gives the probability of getting the data from the model, regardless of the specific values of the parameters, i.e., it measures how well the model explains the data. The absolute value of the evidence is not of interest, so we generally deal with the ratio of two evidences for two models, known as the Bayes factor. The evidence is a far better measure of the suitability of a model than is the p-value (Jeffreys 1961; Kass & Raftery 1995; Winkler 1972; MacKay 2003; Bailer-Jones 2009).

It is worth stressing again that by averaging over the parameters, the evidence is not sensitive to model complexity per se. This is in contrast to the likelihood at the best-fitting parameters (the maximum likelihood): a more complex (flexible) model will always fit the data better, and so will always deliver a higher maximum likelihood. The evidence reports the average likelihood, so it will only increase if the extra complexity gives a net benefit over the plausible parameter space. The model complexity does not then need to be considered separately in some ad hoc way.

3.2.1. Discrete Data Sets

As explained in Section 2.1, each event in a discrete data set can be considered as a Gaussian distribution. For event j this gives the probability that the mass extinction occurs at time τ_j, given that the true time is t_j and the uncertainty in our measurement is σ_j. That is,

$$\begin{equation} P(\tau _j | \sigma _j, t_j) = \frac{1}{\sqrt{2\pi }\sigma _j} \, e^{-(\tau _j - t_j)^2/2\sigma _j^2}. \end{equation} \tag{ 2 }$$

In order to compare the measurement of this event with the predictions of a model, we calculate the likelihood. This is given by an integral over the unknown true time

$$\begin{eqnarray} P(\tau _j | \sigma _j, \theta, M)&=&\int _{t_j} P(\tau _j | \sigma _j, t_j, \theta, M) P(t_j | \sigma _j, \theta, M) {d} t_j \nonumber \\ &=&\int _{t_j} P(\tau _j | \sigma _j, t_j) P(t_j | \theta, M) {d}t_j. \end{eqnarray} \tag{ 3 }$$

The first term in the integral—Equation (2)—is sometimes called the measurement model. The second term is the prediction of the time series model, i.e., the probability (per unit time) that a mass extinction occurs at time t_j. Our time series models are, therefore, stochastic in the sense that they do not attempt to predict when mass extinctions occurred, but rather how the probability of occurrence of a mass extinction varies over time. Note that the likelihood is just measuring the degree of overlap between the data and the model predictions, averaged over all time.

Equation (3) gives the likelihood for a single event. Assuming all events are measured independently,³ then the likelihood for all the data is just the product of the event likelihoods:

$$\begin{equation} P(D|\theta, M) = \prod _j P(\tau _j| \sigma _j,\theta, M), \end{equation} \tag{ 4 }$$

where D = {τ_j}. For the sake of the likelihood and evidence calculation we do not consider the {σ_j} as data, although they are of course measured. That is because D is defined as just those quantities that are predicted by the measurement model.

3.2.2. Continuous Data Sets

Mathematically, the likelihood for the continuous data is very similar to the discrete case, but the interpretation is different.

Consider the measurement model, P(τ|σ, t), in Equation (2) (we consider just one event, so we drop the subscript j). We have interpreted this as the probability of a discrete extinction event being measured at τ, but we could equivalently interpret it as the probability density (i.e., probability per unit time) of extinction at time τ. Now, rather than characterizing the probability density as a Gaussian with mean t and standard deviation σ, we could consider an arbitrary function, characterized by a series of top-hat functions, {p_i} (a histogram), each top-hat characterized by a center t_i, height r_i, and width δ_i. We can then replace P(τ|σ, t) with ∑_ip_i(τ|t_i, δ_i), where

$$\begin{eqnarray} p_i(\tau | t_i, \delta _i) = \left\lbrace \begin{array}{@{}ll@{}}r_i & {\rm when} & t_i - \delta _i/2 < \tau < t_i + \delta _i/2\\ \ms 0 & {\rm otherwise.} \end{array} \right. \end{eqnarray} \tag{ 5 }$$

We now see that

$$\begin{equation} \lim _{\delta _i \rightarrow 0} \sum _i p_i(\tau | t_i, \delta _i) = E(t), \end{equation} \tag{ 6 }$$

i.e., we get a continuous function of the variation of the extinction probability with time t (as τ and t become equivalent in this limit). The likelihood is then

$$\begin{equation} P(D|\theta, M) = \int _{t} E(t) P(t | \theta, M) {d}t. \end{equation} \tag{ 7 }$$

In practice, we characterize E(t) using the extinction rate, r_j, tabulated at each time τ_j, which is equivalent to assuming that extinction rate is constant over the substage (or that we have zero uncertainties on the measured times).

We can actually apply this interpretation to the discrete data sets too. In both cases, the data provide the variation of extinction probability (per unit time) as a function of time, or something proportional to that. The proportionality constant is irrelevant because we keep the data fixed when comparing different models using the evidence.

The time series models appear in the equation for the likelihood in the form P(t|θ, M), i.e., the extinction probability (per unit time) as a function of time predicted by model M at parameters θ. It is important to realize that this probability density function (PDF) over t is normalized, i.e., integrates over all time to unity. This is key to the model comparison because a model which assigns a lot of probability to extinctions at some particular time must necessarily assign lower probability elsewhere. This follows because we are not trying to model the absolute value of the extinction rate, but just its relative variations.

In addition to specifying the functional form of the models, we must also specify the prior probability distribution of the model parameters, P(θ|M). This describes our prior knowledge of the relative probability of different parameter settings. For example, given the timescale in the data, we are not interested in models with timescales less than a few million years or more than a few hundred million years. It is often difficult to be precise about priors, and the evidence and therefore Bayes factors often depend on the choice. The choice of prior must therefore be considered part of the model (e.g., "periodic model with permissible periods between 10 and 100 Myr" is distinct from "periodic model with permissible periods between 50 and 60 Myr"). We investigate the sensitivity of the results to changes in the prior in Section 5.3.

Table 2 summarizes the functional form of the models, which are now briefly described. Figure 2 plots examples of some of these models. The range of the data is taken to be 0–550 Myr BP.

• Uniform. Constant extinction PDF over the range of the data. This has no parameters.
• RB/RNB. Random model in which a set of N times is drawn at random from a uniform distribution extending over the range of the data. A Gaussian with standard deviation σ = 10 Myr is assigned to each of these, and then a constant B added before normalizing. This is the RB model. The RNB ("random no background") model is just the special case of B = 0, which produces a model which is similar to our discrete data. In practice, we fix N and B and calculate the evidence by averaging over a large number of realizations of the model. Specifically, when modeling the B5 and B18 data sets, we fix B = 0 and N = 5, 18, respectively.
• PB/PNB. Periodic model of period T and phase β (model PNB). There is no amplitude parameter because the model is normalized over the time span of the data. Adding a background B to this simulates a periodic variation on top of a constant extinction probability (model PB).
• QPM. A quasi-periodic model in which the phase is a sinusoid with amplitude A_Q, period T_Q, and phase θ (it becomes the same as the PNB model if A_Q = 0).
• SP. A monotonically increasing or decreasing nonlinear trend in the extinction PDF using a sigmoidal function characterized by the steepness of the slope, λ, and the center of the slope, t₀. In the limit that λ becomes zero, the model becomes a step function at t₀, and in the limit of very large λ becomes the Uniform model.
• SSP. Combination of SP and PNB.
• OM(P)/SOM(P). The orbital/semi-orbital model with/without spiral arms, defined in Section 4.5.

Zoom In Zoom Out Reset image size

Table 2. The Mathematical Form of the Time Series Models and Their Corresponding Parameters

Model Name	Pu(t|θ, M)	Parameters
Uniform	1	None
RNB/RB	$\sum _{n=1}^{N}\mathcal {N}(t; \mu _n,\sigma)$+B	σ, N, B
PNB/PB	1/2{cos [2π(t/T + β)] + 1}+B	T, β, B
QPM	1/2{cos [2πt/T + AQcos (2πt/TQ) + β] + 1}	T, β, AQ, TQ
SP	$[1+e^{(t-t_0)/\lambda }]^{-1}$	λ, t0
SSP	PNB+SP	T, β, λ, t0
OM(P)/SOM(P)	$n(\overrightarrow{r_\odot }(t),\overrightarrow{v_\odot }(t))$	$\overrightarrow{r_\odot }(t=0), \overrightarrow{v_\odot }(t=0)$

Note. Time t increases into the past and P_u(t|θ, M) is the unnormalized extinction probability density predicted by the model.

Download table as:  ASCIITypeset image

Except for the orbital model, we adopt a uniform prior over all model parameters over the range specified in Table 3.

Table 3. Range of Parameters Adopted in the Model Prior Parameter Distributions

Model Name	Range of Prior
Uniform	None
RNB/RB	σ = 10 Myr, N ∈ {5, 18}, $B\in \lbrace 0,\frac{1}{\sqrt{2\pi }\sigma }\rbrace$
PNB/PB	10 < T < 100, 0 < β < 2π, B ∈ [0, 1]
QPM	10 < T < 100, 0 < β < 2π, 0 < AQ < 0.5, 200 < TQ < 500
SP	−100 < λ < 100, 100 < t0 < 500
SSP	10 < T < 100, 0 < β < 2π, −100 < λ < 100, 100 < t0 < 500
OM(P)/SOM(P)	Initial conditions (see Table 5)

Notes. Except for OM(P)/SOM(P), a uniform prior for all parameters for all models is adopted which is constant inside the range shown, and zero outside. The prior PDF of parameters in the OM(P)/SOM(P) model is Gaussian and specified by the uncertainties in the initial conditions.

Download table as:  ASCIITypeset image

The integral in Equation (1) is a multidimensional integral over the parameter space, and cannot be calculated analytically. As in Bailer-Jones (2011a), we estimate it using a Monte Carlo method, by drawing model parameters at random from the prior distribution and calculating the likelihood at each. If the set of N parameter draws is denoted {θ}, then Equation (1) can be approximated as the average likelihood,

$$\begin{equation} P(D | M) \, \simeq \, \frac{1}{N} \sum _{\theta } P(D | \theta, M). \end{equation} \tag{ 8 }$$

In the following simulations, we adopt N = 10,000 unless noted otherwise.

We adopt an analytic Galaxy potential, Φ_G(R, z) (described in cylindrical coordinates), comprising three components:

$$\begin{equation} \Phi _G=\Phi _b+\Phi _h+\Phi _d. \end{equation} \tag{ 9 }$$

The first two components are spherically symmetric distributions which represent the bulge and halo using Plummer's model (Plummer 1911),

$$\begin{equation} \Phi _{b,h}=-\frac{{GM}_{b,h}}{\sqrt{R^2+z^2+b_{b,h}^2}}, \end{equation} \tag{ 10 }$$

and the third is an axisymmetric disk according to Miyamoto & Nagai (1975):

$$\begin{equation} \Phi _d=-\frac{{GM}_d}{\sqrt{R^2+\left(a_d+\sqrt{z^2+b_d^2}\right)^2}}. \end{equation} \tag{ 11 }$$

In Equations (10) and (11), R is the galactocentric radius and z is the distance from the Galactic midplane. M_b, M_h, and M_d are the mass of the bulge, halo, and disk, respectively. a_d and b_d are the scale length and height (respectively) of the disk, and b_d and b_h are the scale lengths of the bulge and halo, respectively. We adopt the numerical values for these as given in Table 4 (see García-Sánchez et al. 2001). In addition to the these three components, we introduce into some of the models a time-dependent potential due to the Galactic spiral arms, denoted Φ_s(R, ϕ, z, t). It is defined below in Section 4.3.

There is, of course, significant uncertainty not only in the value of the parameters of the potential but also in the functional form of our model. It is no doubt a simplification of the true potential of the Galaxy, so specific numerical values of quantities such as orbital periods and amplitudes inferred should not be interpreted too literally. Below, we investigate the sensitivity of our model comparison to changes in the potential.

To calculate the motion of a body through the potential from given initial conditions, we solve Newton's equations of motion, which in cylindrical coordinates are

$$\begin{eqnarray} \ddot{R} -R\dot{\phi }^2 &=& -\frac{\partial \Phi }{\partial R} \nonumber \\ R^2 \ddot{\phi }+2R\dot{R}\dot{\phi } &=& -\frac{\partial \Phi }{\partial \phi } \nonumber \\ \ddot{z} &=& -\frac{\partial \Phi }{\partial z}. \end{eqnarray} \tag{ 12 }$$

We solve these equations by numerical integration using the lsoda method implemented in the R package deSolve, with a time step of 0.1 Myr.

The initial conditions are the current phase space coordinates (three spatial and three velocity coordinates) of the Sun. These are derived from observations with a finite accuracy, so our initial conditions are Gaussian distributions, with mean equal to the estimated coordinate and standard deviation equal to its uncertainty (Table 5). In order to calculate an orbit, we draw the initial conditions at random from these prior distributions, and a large number of draws give us a sampling of orbits which will be used later (e.g., in the evidence calculations).

Table 5. The Current Phase Space Coordinates of the Sun, Represented as Gaussian Distributions, and Used as the Initial Conditions in our Orbital Model

	R	VR	ϕ	$\dot{\phi }$	z	Vz
	(kpc)	(kpc Myr−1)	(rad)	(rad Myr−1)	(kpc)	(kpc Myr−1)
Mean	8.0	−0.01	0	0.0275	0.026	0.00717
Standard deviation	0.5	0.00036	0	0.003	0.003	0.00038

Download table as:  ASCIITypeset image

We derive our initial conditions from a number of sources: the distance to the Galactic center comes from astrometric and spectroscopic observations of the stars near the black hole of the Galaxy (Eisenhauer et al. 2003). The Sun's displacement from the galactic plane is calculated from the photometric observations of classical Cepheids by Majaess et al. (2009). The Sun's velocity is calculated from Hipparcos data by Dehnen & Binney (1998).

The model for the spiral arms is described by their geometry and their gravitational potential. However, for the arm crossing periodicity test in the next section, we ignore the mass of the spiral arms when calculating the solar orbit and only consider their location. Likewise, in one class of variants of the orbital model, OM and SOM (defined later), we ignore the arms entirely (for both the orbit and stellar density calculations). This is done so that we can see the additional effect of the arms, the form and mass of which are poorly determined by current observations.

The geometric model comprises two logarithmic spiral arms, the positions of which in circular coordinates, (R, ϕ), are given by

$$\begin{equation} \phi _s(R)=\alpha \log (R/R_{\min })+\phi _{\min }, \end{equation} \tag{ 13 }$$

where α is a winding constant, R_min is the inner radius, and ϕ_min is the azimuth at that inner radius. The radius of the spiral arm ranges from R_min to R_max. Of the various arm models offered by Wainscoat et al. (1992), we selected the main two spiral arms, 1 and 1', with ϕ_min given by Vanhollebeke et al. (2009) and other parameters given by Wainscoat et al. (1992; see Table 6). Their location in the plane of the Galaxy is shown in the left panel of Figure 3. The arms rotate rigidly with a constant angular velocity (pattern speed) of Ω_p = 20 km s⁻¹ kpc⁻¹ (Martos et al. 2004; Drimmel 2000).

Zoom In Zoom Out Reset image size

The gravitational potential of the arms is described by the (first term of the) analytic potential of Equation (8) of Cox & Gómez (2002). This is

$$\begin{eqnarray} \Phi _s &=& -\frac{4\pi GH}{K_1 D_1} \rho _0 e^{-\frac{R-r_0}{R_s}} \nonumber\\ &&\times \cos (N[\phi -\phi _s(R,t)])\left[{\rm sech}\left(\frac{K_1 z}{\beta _1}\right)\right]^{\beta _1}, \end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray*} K_1 &=& \frac{N}{R \sin (\alpha)},\\ \beta _1 &=& K_1 H(1+0.4K_1H),\\ D_1&=& \frac{1+K_1H+0.3(K_1H)^2}{1+0.3K_1H},\\ \phi _s(R,t)&=&\phi _s(R)+\Omega _p t. \\ \end{eqnarray*}$$

The amplitude of the spiral density distribution is

$$\begin{equation} \rho _A(R, z) = \rho _0 e^{-\frac{R-r_0}{R_s}} {\rm sech}^2(z/H), \end{equation} \tag{ 15 }$$

where R_s is the radial length of the drop-off in the density amplitude of the arms, ρ₀ is the midplane arm density at the fiducial radius r₀, and H is the scale height of the spiral density. This is modulated by a sinusoidal pattern in ϕ to give the overall density which corresponds approximately⁴ to the above potential:

$$\begin{equation} \rho _s(R,\phi, z, t) = \rho _A(R,z)\cos [N(\phi -\phi _s(R,t))], \end{equation} \tag{ 16 }$$

where N is the number of arms and ϕ_s(R) is given in Equation (13). The values of the relevant parameters are given in Table 6.

In previous studies of the impact of astronomical phenomena on the terrestrial biosphere, it has frequently been assumed that the solar motion shows strict periodicities in its motion perpendicular to the Galactic plane, and sometimes also with respect to spiral arm crossings. We investigate this here using our numerical model.

For each orbit k, we calculate the intervals between successive crossings, {Δtⁱ}_k (separately for midplane and spiral arm crossings), where i indexes the crossing. We then calculate the sample mean and sample standard deviation of these intervals for each orbit:

$$\begin{eqnarray} \overline{\Delta t_k} &=& \frac{1}{N_k-1}\sum _{i=1}^{N_k-1} \Delta t_k^i \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} \sigma _k &=& \sqrt{\frac{1}{N_k-2} \sum _{i=1}^{N_k-1} \left(\Delta t_k^i-\overline{\Delta t_k}\right)^2} \, \end{eqnarray} \tag{ 18 }$$

where N_k is the number of crossings in the kth orbit. To assess the periodicity of the crossing intervals, we define the degree of aperiodicity as

$$\begin{equation} a_k=\sigma _k/\overline{\Delta t_k}. \end{equation} \tag{ 19 }$$

An orbit with a = 0 is strictly periodic.

We investigated the variation of the aperiodicity of the solar orbit with the six parameters (initial conditions). This parameter space is too large to report on extensively here, but we find that the aperiodicity is most sensitive to R(t = 0) and $\dot{\phi }(t=0)$. In the following, we vary these initial conditions individually by drawing 10⁴ samples from the corresponding initial condition distribution. (Larger sample sizes did not alter the results significantly.) We simulate the solar orbits using the arm-free potential, Φ_G.

4.4.1. Midplane Crossings

Some earlier studies claimed that Galactic midplane crossings trigger increases in terrestrial extinction due to an enhanced gamma-ray or cosmic-ray flux or due to a larger perturbation of the Oort Cloud. These are directly related to the increased stellar density and increased occurrence of star-forming regions. The larger tidal forces are postulated to enhance the disruption of the Oort Cloud (Rampino & Stothers 2000; Matese et al. 1995), and the higher density of massive stars—and thus high-energy radiation as well as increased SN rate—raises the average flux the Earth is exposed to. The periodicity of the Sun's vertical motion—not least its period, phase, and the assumed stability of this period—is central to these claims. We examine these using our model.

The results of varying just the initial galactocentric radius of the Sun, R(t = 0), are shown in Figure 4. We see in the top right panel that about 90% orbits have an aperiodicity less than 0.1. In the lower two panels, we see how a varies with the value of the initial condition and with the average crossing interval.

Zoom In Zoom Out Reset image size

The aperiodicity is 0.002 (nearly strict periodicity) at $\overline{\Delta t}=61.8$ Myr. This corresponds to a 1:1 resonance between the vertical motion and the radial motion. Its value is close to a period in the biodiversity data of 62 ± 3 Myr claimed by Rohde & Muller (2005). Little should be read into this coincidence, however, as there is no good (i.e., independent) reason to select the specific initial condition that leads to this period over any other. Moreover, changing the parameters of the Galactic potential—which is not very well known—changes this period. (For example, if we increase the mass of the Galactic halo, then the values of $\overline{\Delta t}$ are decreased.) The other minimum in the aperiodicity in the bottom panels is 0.02 at $\overline{\Delta t}=100.2$ Myr. This corresponds to an approximately circular orbit in the midplane. If we set V_R(t = 0) = 0, V_z(t = 0) = 0 and z = 0, then this solar orbit would be strictly circular.

The cumulative curve (top right panel of Figure 4) makes a sharp turn at a' = 0.1. This is because of a sudden decrease in the number of orbits with large aperiodicities. Similarly, the discontinuities in the lower panels are caused by changes in the (small) number of discrete plane crossings which occur for different aperiodicity ranges.

If we now vary the initial condition $\dot{\phi }(t=0)$ instead, then the periodicity test gives very similar results: we find a nearly strict periodicity at $\overline{\Delta t}=60$ Myr and another minimum in the aperiodicity at about 100 Myr. That means the nearly strict periodicity is mainly determined by a combination of R(t = 0) and $\dot{\phi }(t=0)$.

In summary, we see that the majority of the simulated orbits (90%) are quite close to periodic (a ⩽ 0.1) in their motion vertical to the midplane, although strict periodicity essentially never occurs.

4.4.2. Spiral Arm Crossings

Regarding spiral arms as regions of increased star formation activity and stellar density, the mechanisms of mass extinction considered for midplane crossing could likewise be applied to spiral arm crossings, and have been by some authors (Leitch & Vasisht 1998; Gies & Helsel 2005). However, such studies have oversimplified the solar motion by failing to take into account the considerable uncertainties in the current phase space coordinates of the Sun, and thus in its plausible orbits. Some studies have even claimed a connection between spiral arm crossings and the terrestrial biosphere after having fit the solar motion to the geological data, but such reasoning is clearly circular.

We examine here the periodicity of spiral arm crossings (although we note that some studies in the literature claiming a spiral arm-extinction link just consider the crossing times and do not claim a periodic crossing). The crossing intervals are longer than with the midplane, so we only include in our analyses models in which there are at least three arm crossings. We assume that the arms have indefinite vertical extent, so that a crossing on the x–y plane is always a true encounter. In reality, the Sun might pass over or under the arms, thus reducing the overall relevance of spiral arm crossings to terrestrial extinction.

Figure 5 shows the result of this analysis for the 7407 orbits (out of the original sample of 10,000) that exhibit at least three arm crossings. The cumulative probability (top right panel) shows that about 40% of the orbits have an aperiodicity larger than 0.2. In other words, it is not very likely that the solar orbit and spiral arms are so tuned to give periodic crossings. The lower two graphs show how a varies with R(t = 0) and $\overline{\Delta t_k}$. The numerous gaps in these plots are a consequence of the fact that not all orbits for certain ranges of R(t = 0) had at least three arm crossings, and so were removed from the analysis. We see, therefore, that the crossing interval is very sensitive to R(t = 0).

Zoom In Zoom Out Reset image size

Note that we have neglected the mass of the arms in the orbital calculations. When we include it, the values of aperiodicity increase and there is an even less clear dependence of a on R(t = 0) or $\overline{\Delta t_k}$.

In summary, we find it unlikely that spiral arm crossings are even close to periodic. If the pattern speed of the spiral arms has not been constant in the past 550 Myr, or if the pattern itself has not been stable, then this conclusion is strengthened further.

4.5.1. Derivation of the Extinction Rate from the Stellar Density Variation

As outlined in Section 1, various astronomical mechanisms for biological extinction have been identified, including comet impacts (from Oort Cloud perturbation), gamma rays (from SNe or GRBs), and cosmic rays (from SN remnants; Ellis & Schramm 1995; García-Sánchez et al. 2001; Gies & Helsel 2005). The intensity of all of these depends on the local stellar density. If we consider a general mechanism involving flux from nearby stars, then the flux from a single star is proportional to f/d², where f is the relevant surface flux and d the distance. The sum of this over the whole relevant volume of space around the Sun is proportional to the total intensity and thus the extinction probability (per unit time).

Let us assume that the extinction rate, E, is linearly proportional to the flux, and that the number density of relevant stars is proportional to the total stellar number density (stars per unit volume), n. Because the density of spiral arms is much less than the density of the other components, we consider at first only the time-independent density arising from the halo, disk, and bulge. The density is calculated from the corresponding potential (defined in Section 4.1) using Poisson's equation.

In an axisymmetric cylindrical coordinate system, the extinction rate at the Sun is then

$$\begin{eqnarray} E(R_\odot, z_\odot) &=& C \int \!\! \int \frac{n(R,z)}{d^2} \,R\,dR\,dz\nonumber \\ &=&C \int \!\! \int \frac{n(R,z)}{(R-R_\odot)^2+(z-z_\odot)^2}\,R\,dR\,dz,\nonumber \\ \end{eqnarray} \tag{ 20 }$$

where R and z are the galactocentric radius and height above the midplane, respectively, for some star, R_☉ and z_☉ are the corresponding (time-varying) coordinates of the Sun, and C is a constant. Note that the stellar number density, n, is proportional to the corresponding stellar density, ρ. Defining the distance from a star to the Sun as r ≡ R − R_☉ and Z ≡ z − z_☉, the extinction rate is

$$\begin{eqnarray} E(R_\odot,z_\odot) &=& C \int \!\! \int \frac{n(R_\odot +r,z_\odot +Z)}{r^2+Z^2}(R_\odot +r) drdZ.\nonumber \\ \end{eqnarray} \tag{ 21 }$$

The flux from a star falls off as 1/d², but we can truncate this integral at some upper distance because at some point the flux is too weak to influence the terrestrial biosphere. We take d_th = 50 pc as an upper limit.⁵ This is much smaller than the scale length of the disk and comparable to the scale height of the disk (see Table 4), so we can approximate n(R_☉ + r, z_☉ + Z)(R_☉ + r) by n(R_☉, z_☉ + Z)R_☉. The integral then becomes

$$\begin{equation} E(R_\odot,z_\odot)\simeq {CR}_\odot \int _{-d_{\rm th}}^{d_{\rm th}}\int _{-d_{\rm th}}^{d_{\rm th}}\frac{n(R_\odot,z_\odot +Z)}{r^2+Z^2} drdZ. \end{equation} \tag{ 22 }$$

Integrating over r gives

$$\begin{equation} E(R_\odot,z_\odot)\simeq 2CR_\odot \int _{-d_{\rm th}}^{d_{\rm th}}n(R_\odot,Z+z_\odot)\frac{\arctan (d_{\rm th}/Z)}{Z}\,dZ. \end{equation} \tag{ 23 }$$

The geometric factor (arctan(d_th/Z))/Z drops close to zero at about Z = 25 pc, and does so much more rapidly than the stellar density term, which follows the vertical profile of the disk (which has a much larger scale height of 250 pc). Thus to a reasonable degree of approximation, we can set n(R_☉, Z + z_☉) ≃ n(R_☉, z_☉) in this integral. The integral is then just over the geometric factor, which gives some constant (dependent on d_th, but of no further interest). Thus we are left with

$$\begin{equation} E(R_\odot,z_\odot)\simeq C^{\prime }R_\odot n(R_\odot,z_\odot), \end{equation} \tag{ 24 }$$

for some constant C'. For the solar motion, the relative variation of R_☉ is less than that of n(R_☉, z_☉), so we have

$$\begin{equation} E(R_\odot,z_\odot) \propto n(R_\odot, z_\odot). \end{equation} \tag{ 25 }$$

In other words, the extinction rate is just proportional to the stellar density at the location of the Sun. The approximations in Equations (22)–(25) still hold when we include the low-density spiral arms defined in Section 4.3, in which case we must also introduce the explicit dependence on azimuth and time:

$$\begin{equation} E(R_\odot, \phi _\odot, z_\odot, t) \propto n(R_\odot,\phi _\odot,z_\odot,t). \end{equation} \tag{ 26 }$$

In the above model, we assumed that the extinction rate is proportional to d⁻², i.e., the influence falls off like a flux on the surface of a sphere. We could generalize this dependence to be d^−k/2 for k ⩾ 0 in order to reflect other mechanisms, e.g., tidal effects.

In order to test the validity of the above approximations, we compare in Figure 6 the extinction rate as given by Equation (20) (by numerical integration) with the stellar number density n(R_☉, z_☉). We plot over ranges of R_☉ from 5 kpc to 10 kpc and z_☉ from −0.5 kpc to 0.5 kpc, in accordance with the ranges covered by the simulated solar orbits. We normalize the extinction rate (and the stellar density) by setting its integral over R_☉ and z_☉ to be unity. In the upper row of Figure 6, the difference between the stellar density and the extinction rate reaches a maximum in the midplane (z_☉ = 0); this is on account of the relatively large density gradient at z_☉ = 0. The maximum difference is only about 10% of the peak value of stellar density for all values of k. In the lower row, the largest difference is at the lower limit of R_☉. Note that the value of k has very little impact.

Zoom In Zoom Out Reset image size

In practice, most of the simulated orbits spend most of their time in the region 7 < R_☉/kpc < 9 and −0.3 < z_☉/kpc < 0.3, where the differences between local stellar density and extinction rate variation are even smaller. Thus to within a few percent, the stellar density at the Sun is a good predictor of the extinction rate. The time variation of this density is the time series model and forms the basis for what we refer to as the "orbital models," the forms of which we now define.

4.5.2. Definition of OM(P) and SOM(P)

The orbital model "OM" is the orbital model which does not include the spiral arm at all, neither in the gravitational potential (for calculation of the orbits) nor in the stellar density (for the extinction rate calculation). The orbital model OMP does include the spiral arm in both senses. Thus both OM and OMP are internally self-consistent.

Once normalized, E(t) is just the quantity P(t|θ, M) in Section 3.3 (and it is normalized to give the unit integral over the span of the data). The parameters of OM and OMP are the initial conditions of the orbit, and the corresponding priors are the Gaussian distributions summarized in Table 5. Thus, one orbit calculated from one draw of the initial conditions allows us to calculate one likelihood for these models (for given data set). Repeating this and averaging the resulting likelihoods gives the evidence for that orbital model (see Section 3.4).

For both of these models we consider four variations, labeled 1–4, according to which initial conditions we vary (and therefore sample over to build up the set of orbits).

In addition to these models, we define the "semi-orbital model," SOM. This is derived from the OM simply by subtracting from the predicted extinction rate a constant value, h, and setting all resulting negative values to zero. Here we simply set h to be the minimum value of the extinction rate (see Figure 7). This is intended to model the situation in which the flux causing the extinction must rise above some threshold before it has an effect. (We might consider this as an adaption of life to the extraterrestrial flux background.) In analogy to OM, SOM excludes the spiral arm. SOMP is SOM with the spiral arm potential and density included. Once again, we will consider four varieties according to which initial conditions are varied.

Zoom In Zoom Out Reset image size

We now calculate the Bayesian evidence (Equation (1)) for the various models for each data set. This is done by sampling from prior probability distributions of the model parameters (P(θ|M), Table 3), calculating the likelihood (Equation (4) for discrete time series, Equation (7) for continuous time series), and then averaging these for that model and data set (Equation (8)).

To calculate the evidences for the RNB and PNB models for the B5 and B18 data sets, we adopt a Monte Carlo sample size of 10⁶. In all other cases, we use a sample size of 10⁴. Larger sample sizes did not alter the estimated evidence significantly.⁶ This sample size is given according to the sensitivity test of the evidence to the sample size for all models and all data sets. This test shows that the evidence estimated from 10⁴ draws in the prior distribution is close to the real evidence when there is a background either in the model or in the data set.

As the absolute value of the evidence is not of interest, we report the ratio of evidence, the Bayes factor. Here we report Bayes factors with respect to the Uniform model. We regard a model as being significantly better than another when its evidence exceeds that of the other by a factor of 10 (Jeffreys 1961; Kass & Raftery 1995). Note that it is only meaningful to compare evidences—and therefore Bayes factors—for a fixed data set.

The results are shown in Table 7. For the reference models we evaluate the evidence by sampling over all their model parameters, but in the case of OM and SOM we sample over just some of the parameters (initial conditions), keeping the others fixed, in order to investigate the impact of the different parameters. As shown in Section 4.4, the periodicity of the solar orbit is most sensitive to the initial conditions R(t = 0) and $\dot{\phi }(t=0)$. We therefore calculate the evidence for the OM (and SOM) models with four different sets of initial conditions being varied: R(t = 0) only; $\phi \dot{(}t=0)$ only; {R(t = 0) and $\dot{\phi }(t=0)\rbrace$; {R(t = 0), V_R(t = 0), $\dot{\phi }(t=0)$, and V_z(t = 0)}. In all cases, we fix ϕ(t = 0) and z(t = 0), the former because it has no impact on the solar motion in this axisymmetric potential, and the latter because the uncertainty in the current z position of the Sun has a limited impact on the subsequent orbit. To assess the effect of the spiral arm perturbation on the Bayes factors, we have selected four perturbed orbital models, OMP1, OMP4, SOMP1, and SOMP4, to compare with corresponding unperturbed orbital models.

For the B5 data set, the Bayes factors of all time series models relative to the Uniform model are less than 10. Thus, none of these models are a significantly better explanation of the data. One model, RNB, has a Bayes factor of less than 0.1, indicating that we can discount this one as being an unlikely explanation. Given that the Uniform model is the simplest model of the set, the principle of parsimony suggests that we should be satisfied with it as an explanation. This does not deny the possibility that some other model shows significantly higher evidence. After all, we can only ever make claims about models which we explicitly test.

The B18 data set includes more extinction events than the B5 data set, and not surprisingly it discriminates more between the models (the Bayes factors show a larger spread). (These results are also shown graphically in the upper panel of Figure 8.) The OM models are favored somewhat more than the other models—e.g., the Bayes factor of OM3 to SP is 0.63/0.019 = 33—although again no model is favored significantly more than the Uniform model. In contrast, several models are significantly disfavored (RNB, SP, OMP1, OMP4, SOM1, SOM2, SOMP1, and SOMP4). In particular, the perturbed orbital models, including OMP1, OMP4, SOMP1, and SOMP4, are less favored by the data than their corresponding unperturbed orbital models. All the other perturbed orbital models (not listed in Table 7) also have lower Bayes factors than the unperturbed orbital models.

Zoom In Zoom Out Reset image size

For the two continuous time series, RM and A08, the difference between the evidences of all of the models is not significant.⁷ Our broad conclusion is that no model significantly outperforms the Uniform model on any of the data sets. On the contrary, a few can be "rejected" on the grounds of significantly lower evidence. Recall that all of these models are predicting the extinction probability (per unit time). In terms of the discrete data sets, the Uniform model just means that the mass extinction events occur randomly in time. We find this to be no less probable than a periodic or quasi-periodic variation of the probability, or a monotonic trend in the probability, etc. In terms of the continuous data sets, we obviously do not believe that the Uniform model is a good explanation of the clearly apparent variations in the extinction rate (see Figure 1). But the analysis does tell us that this is no worse an explanation than the more complex models of the variation considered, such as periodic, orbital-model based, etc. Clearly, there must be yet other models which could explain the data even better. This may explain why previous authors have found an apparent periodicity in the data: the periodic model can explain the data to some degree, but actually no better than simpler models.

We have seen that the evidence hardly discriminates between any of the models on the continuous data sets, and only between some of them on the discrete data sets. (This is by no means inevitable. In other problems, the evidence can vary enormously between models.) This means that on average over their parameter space, the models differ little in their predictions. It is nonetheless interesting to see how the likelihood varies over the parameter space. (We would do this in particular to find the best-fitting parameters, although these are only meaningful if the overall model has been identified as the best explanation of the data.) We focus here mainly on the PNB and OM models for the B18 data. We again normalize the likelihood for a model by dividing it by the likelihood of the Uniform model to form the likelihood ratio. As the latter model has no parameters, it is likelihood constant and equal to its evidence. We denote the maximum value of the likelihood ratio as MLR.

Figure 9 shows how the likelihood varies over the two-dimensional space formed by the two parameters, period and phase, of the PNB model. There is significant variation. We see numerous local maxima, the largest likelihoods being around {T/Myr, β/rad} = {60, 4.5}. However, these maxima are rather narrow, so once the (much lower) likelihood in the other (equally plausible) regions are taken into account, the overall evidence for the model is not particularly high. If we are interested in the variation of likelihood with just the period, then we can marginalize this diagram over phase, and plot with respect to period, thereby forming a (Bayesian) periodogram (lower panel). We see a clear peak around 60 Myr. This is coincident with the period of 62 ± 3 Myr identified by Rohde & Muller (2005). It is tempting (but incorrect) to associate this peak value of the likelihood with the periodic model as a whole, and use it to claim a larger evidence for the periodic model. Certainly there is a degree of arbitrariness in the prior parameter distribution—in this case a uniform distribution—and narrowing this range around this peak would clearly increase the evidence. For example, if we truncate the period range from its current value of [10, 100] Myr to [50, 80] Myr, then the Bayes factor relative to the Uniform model increases from 0.62 to 1.5. This is a rather modest increase, but we could increase it to a significantly high value with an even narrower prior. However, we may not use the data to find the best-fitting parameters and then claim that we should only consider the model near to these. We would need some other reasoning or independent data to make such a selection. (The Rohde & Muller (2005) time series is not independent of B18 because both are based on the same paleontological data.) We do not see how, a priori, we could limit the plausible periods of periodicity to something as narrow as 50–80 Myr, let alone the much narrower range required to favor PNB significantly over other models. In the extreme limit of an infinitesimal region around the maximum likelihood, we end up doing model comparison using the maximum likelihood. Just out of interest, these values are shown in Table 7 and plotted in Figure 8. If we were to use this (incorrect) metric, then PNB and some of its variants have significantly higher likelihood than the Uniform model and several of the other models (although barely more than a factor of 10 above the random model, RB). Another way of seeing why this is the wrong approach was already discussed in Section 3.1: by focusing on the best fits, we simply favor the more complex model. We could always define a more flexible model and so fit even better.

Zoom In Zoom Out Reset image size

We labor this point because many of the claims for a periodicity in biodiversity data have made use of a maximum likelihood approach (of which χ² is a special case) or something equivalent. We must instead use the evidence for model comparison. (Maximum likelihood may be used for estimating the best parameters once we have established that we have the best model.) If the periodic model were in fact the true one, then of course only one period and phase would be true. In that case, the likelihood around these values would be so high as to result in a large evidence even when averaging over the broader parameter space (see simulations in Section 4.1 of Bailer-Jones 2011a for a demonstration).

Incidentally, the fact that we find a dominant period at all in the B18 data set is actually not that unlikely. The (Bayesian) periodogram of samples drawn at random from a uniform distribution often exhibits a period which has a likelihood larger than that of the true Uniform model (see Section 4.2 of Bailer-Jones 2011a). In other words, it is often possible to explain a random data set with some period, which is just a testament to how flexible the periodic model is.

Moving on from the periodic model, we show in Figure 10 the likelihood distribution for OM1 and OM2, i.e., where we vary the initial conditions R(t = 0) and $\dot{\phi }(t=0)$, respectively. The likelihood ratio varies by a factor of up to 10⁴, but its absolute value is never more than about two. That is, no value of the initial conditions gives a model much more favorable than the Uniform model, whereas some are far less favorable. If we had lower uncertainties on these phase space coordinates of the Sun, then we might be able to conclude something more definitive. For example, if R(t = 0) = 7.5 kpc, then the OM model would be even less favored. We see a similar degree of variability for the other OM models and data sets listed in Table 7.

Zoom In Zoom Out Reset image size

We performed a similar analysis for the other model, but for the sake of space report only the maximum likelihood ratios in Table 7.

For the B5 data set, the RNB model has the highest maximum likelihood ratio, yet its evidence was the lowest. This indicates that while one particular instance of RNB fits the data well, overall it is a poor model.

For the RM and A08 data sets, the evidence is very similar for all models. This means that the data are not able to discriminate between these models very well: they are equally good (or bad). However, the time series analysis model used here is not best suited to these data sets. These data can be better interpreted as valued time series, ones in which we have an extinction magnitude attached to each time (both, in general, with uncertainties), rather than a time variable probability of extinction. Bailer-Jones (2012) has extended the present model in order to work with such data sets; the results of its application to RM and A08 will be reported in a future publication.

In summary, we find that for none of the data sets is any model particularly favored over the simple uniform one. In particular, there is no evidence to suggest that the orbital model (with or without a background extinction level) is a particularly good explanation of the data.

The evidence is of course sensitive to the choice of prior parameter distribution, and often we have little reason to make a very specific choice. Here we test the sensitivity of the evidence to this, as well as to the parameters of the Galactic potential used in the orbital models, and to the age uncertainties for the discrete data sets.

The age uncertainties in the discrete data are taken into account by the likelihood function. However, it is often difficult to estimate uncertainties, and we additionally made a plausible, but not unique, translation of the estimated duration of a stratigraphic substage in order to estimate uncertainty (which is the standard deviation of a Gaussian for each event; see Section 2.1). To see how this affects our results, we scale the age uncertainties in the B18 data set by a constant factor of 1/4, 1/2, 2, and 4. For each of these modified data sets, we calculate the evidences for the models (S)OM1–4 and the Uniform model and recalculate the Bayes factors relative to the Uniform model. These are plotted in the top panel of Figure 11. The Bayes factors change by just a few percent, so a precise age uncertainty is not necessary.

Zoom In Zoom Out Reset image size

As a second test, we scale in the same way the uncertainties of the initial conditions of the orbital models (i.e., we change the width of the prior parameter distribution). The results are shown in the bottom panel of Figure 11. The change in evidence for any particular model is larger than in the previous case, but in most cases is less than a factor of five (except for the models including the spiral arms). Moreover, the absolute value of the Bayes factor remains below one.

Our results are therefore also insensitive to considerable imprecision in the uncertainties in the phase space coordinates of the Sun. This (and the previous conclusion) is also true for the B5 data set.

As a third sensitivity test, we allow the number of simulated random events, N, and the standard deviation of each event, σ, in the RNB and RB models to vary. (Earlier, we fixed N = 5 when drawing models for the B5 data set and N = 18 for the B18 data set.) For the B18 data set, we find that a larger number of peaks or larger standard deviation in the RNB model produces a significantly larger Bayes factor (see Table 8), although it is still below unity. The RB model shows much less sensitivity to N and σ. We similarly recalculate the evidence for the other models in response to various perturbations of their priors, as also listed in Table 8. The resulting Bayes factors do not change by more than a factor of 10 in any case, and often by much less.

Table 8. The Bayes Factors (Relative to the Uniform Model) on the B18 Data Set for Models with Priors Varied

Models	Varied Prior	BF
	None	0.62
	50 < T < 80	1.5
PNB
	10 < T < 200	0.31
	10 < T < 400	0.15
	None	0.80
PB	B = 1/2	0.88
	B = 2	0.97
	None	0.85
	0 < AQ < 1/4	0.62
QPM	0 < AQ < 1	0.54
	100 < TQ < 300	0.61
	100 < TQ < 500	0.58
	None	0.00050
	σ = 5 Myr	3.7 × 10−11
RNB	σ = 20 Myr	0.026
	N = 9	0.0085
	N = 36	0.29
	None	0.40
	σ = 5 Myr	0.73
	σ = 20 Myr	0.25
RB	$B=\frac{1}{2\sqrt{2\pi }\sigma }$	0.13
	$B=\frac{2}{\sqrt{2\pi }\sigma }$	0.45
	N=9	0.84
	N=36	0.20
	None	0.019
	−200 < λ < 0	0.070
SP
	0 < λ < 200	0.10
	10 < t0 < 500	0.047
	None	0.18
	10 < T < 200	0.17
SSP
	−200 < λ < 200	0.37
	10 < t0 < 500	0.27

Note. Each prior is varied individually (listed in the middle column) with the other fixed at their canonical values.

Download table as:  ASCIITypeset image

Finally, we test the sensitivity of the Bayes factors to changes in the parameters of the Galaxy model (canonical values listed in Table 4). The results are shown in Table 9 for the B18 data set. If we double the mass of the halo, for example, then the evidence for the OM and SOM models changes by no more than a factor of three. Some other changes produce smaller effects, some larger, but not more than by a factor of five (and note that a change in a factor of two of the scale lengths is beyond what is consistent with observed data). Changes in the parameters of models which include spiral arms (the OMP and SOMP models) can produce larger changes in the Bayes factors. However, most significantly, none of these changes produce a Bayes factor greater than one. In other words, none of these changes result in the orbital model becoming a better explanation for the paleontological than the Uniform model.

In summary, we find that the evidence for most models is not particularly sensitive to the age uncertainties, Galaxy model parameters, or reasonable changes to the prior parameter distributions.

Here we investigate how well our analysis can discriminate between models, by using simulated data which have been drawn from one of these models. Given sufficient data, such a discrimination will always be possible to some threshold Bayes factor, but here we are interested in the case where the data have similar properties (in particular, sparsity) to the real data we have been using.

We investigate this by simulating a number of time series from the RNB, OM1, and PNB models (which we will refer to as "generative models" when used in this way to distinguish from their use to calculate the evidence on given data). For each generative model, we fix the parameters to certain values, then sample 18 events from the resulting P(t_j|θ, M) to give a simulated time series, to which we then attach the measured age uncertainties. We generate 10 time series in this way (and below we average the Bayes factors over these and report that). For the OM1 and PNB generative models, we repeat this at 10 different values of the solar initial radius parameter (OM1) or period parameter (PNB) in the generative model (RNB has no parameters). We repeat the whole process a second time but using simulated age uncertainties drawn from a log normal distribution with standard deviation and mean calculated from the measured age uncertainties. We finally repeat the process a third time for the OM1 and PNB generative models, but now drawing data to have the same time sampling as our continuous data sets (for which age uncertainties are not used; see Section 3.2.2). (We do not do this for the RNB model as it predicts discrete events.)

For each simulated data set we calculate Bayes factors for the RNB, OM1, and PNB models relative to the Uniform model. For the data generated from the RNB model (with age uncertainties taken from the data), the Bayes factors for the three models are as follows: 0.18 for RNB; 0.57 for OM1; 0.52 for PNB. (We get almost identical values when the age uncertainties were drawn at random.) Thus, no model—not even the true one—is favored over the Uniform model (although none is significantly rejected either). This is not that surprising, however, because with only 18 events, and with the evidence effectively averaging the predicted times of events from the RNB model over all time, the Uniform and RNB models end up with similar predictive power. This is unavoidable because with the RNB model, we cannot decide in advance where the events are: we must average over all possibilities.

The results for applying the models to the data generated from the OM1 and PNB models are shown in Figure 12, where the horizontal axis shows how the Bayes factor varies with the one parameter which is varied in these generative models. The top row shows the results for data drawn from the OM1 model, for the discrete (B18-like) data (left) and the continuous data (right). We see that the (true) OM1 model is not significantly favored in either case (Bayes factor always less than 10), although not disfavored either (Bayes factor more than 0.1). In particular, the continuous data show no discriminative power.

Zoom In Zoom Out Reset image size

The lower two rows show conceptually the same thing, but now for data drawn from the PNB model for two different values of the phase parameter (the two rows). Here we see that at least for longer periods, the PNB model is generally correctly identified (on the basis of a large Bayes factor) when using the discrete data sets. Yet the continuous data still show no discriminative power.

This difference between discrete and continuous data sets is not unexpected. In the former, the likelihood is a product of the likelihood for many events, each of which is the convolution of the event with the model. In the latter, the likelihood is just the result of a single convolution of a continuous model over continuous data. This is not the best approach for modeling continuous data. A better choice is the recently developed continuous time series modeling method described in Bailer-Jones (2012), which will be used in future work.

Clearly, one could perform many more tests with more simulated time series, varying different parameters in the generative models, and with different permutations of the values of the fixed parameters. No doubt there are parts of parameter space where some models are favored over others, particularly, if we adopted more informative priors. Thus, while these results on simulated data give some check on the discriminative power of the method and data, they should not be overinterpreted to say anything too general. Nonetheless, the tests we have done confirm what we concluded based on the analysis of the real data. Specifically, while our analysis of the real data does not allow us to claim evidence in favor of the orbital-based models, it also cannot rule out these models. This is due partly to the lack of predictive power of the data and partly to the large flexibility (or broad prior parameter space) of the models. Better constraints on the solar orbit would help reduce the latter.