Geological and biological data

There are various types of geological and biological data which are used to study climate history or evolution. Widely used measures of biological change are as follows.

• • Species or genera diversity. This measures the number of different species or genera present at any one time. A species is the lowest level taxonomic rank; a genus is composed of many species (above this comes family, order, class, etc).

• • Extinction. This can be expressed either as the number (or fraction) of species or genera which become extinct in a time interval, or as a time series of delta functions marking epochs of mass extinction but without any measure of the extent of the extinction.

• • Origination. As with extinction, but for newly created species or genera.

The above are measured via the fossil record and have been recorded back to at least the beginning of the Phanerozoic eon, some 545 Myr BP (before present). This time marks a significant increase in the diversity of life on Earth, in particular the occurrence of hard-shelled animals. Diversity, extinction and origination are not necessarily correlated. For example, a large extinction concurrent with a large origination event would result in little change in the total diversity.

There are several tracers of geological activity.

• • Impact craters.

• • Iridium layers. It has been proposed that iridium – which is rare in the Earth's crust – could be delivered by a meteor or comet and then spread over a large area following the impact. This was proposed by Alvarez et al. (Reference Alvarez, Alvarez, Asaro and Michel1980) as evidence of an impact at the K-T (Cretaceous–Tertiary) boundary 65 Myr BP.

• • Flood basalts. Giant volcanic eruptions result in basalt lava covering large areas of land or sea crust. As the lava has low viscosity it forms layers rather than a classic volcano, although on giant scales it can form mountain ranges and plateaus. Examples are the Deccan traps in India (occurring around 65 Myr BP) and the Siberian traps (around 250 Myr BP).

• • Orogenic events (mountain building) and plate tectonics. These affect climate because of their influence on atmospheric and oceanic circulation as well as the formation of icecaps at the poles.

• • Geomagnetic reversals, the orientation of the field being preserved in some rocks.

• • ¹⁸O temperature proxy. ¹⁸O is heavier than the much more abundant oxygen isotope ¹⁶O, so evaporation of water leads to a partial separation of these (stable) isotopes. The degree of separation is temperature dependent: the warmer the water, the higher the ¹⁸O content of the evaporated water. Hence the relative abundance of ¹⁸O to ¹⁶O (written δ¹⁸O) may be used as proxy for the ocean temperature (Dansgaard Reference Dansgaard1964). When this water is precipitated out of the atmosphere as rain or snow it can remain permanently frozen in places such as Greenland and Antarctica. There is then a positive correlation between ¹⁸O and temperature, e.g. δ¹⁸O =0.7T−13.6, where δ¹⁸O is measured in parts per thousand and T in Celsius (the constants depend on geographical location). This has been recorded in several ice cores going back as far as 5.3 Myr BP (Lisiecki & Raymo 2005). Some of these cores, e.g. the Vostok and EPICA ice cores, have also measured δ¹⁸O of the air trapped in the ice, thus providing a measure of atmospheric temperatures. The precipitated water can also be taken up by land animals and δ¹⁸O later measured in their fossils, whereby there is also a positive correlation. In contrast, the oceans from which the water evaporated become depleted in ¹⁸O. This can be measured in marine fossils, in particular in the calcite shells of foraminifera (millimetre-sized creatures) on the ocean floor. In this case, δ¹⁸O is negatively correlated with the ocean temperature at the time the organism died (Epstein et al. Reference Epstein, Buchsbaum, Lowenstam and Urey1961). Such fossils have been used to trace climate back as far as 600 Myr BP (e.g. Veizer et al. (Reference Veizer1999), IPCC (Reference Solomon, Qin, Manning, Chen, Marquis, Averyt, Tignor and Miller2007)).

• • Sea level variations, driven either by local uplift (e.g. plate tectonics) or through changes in the amount of water locked in glaciers and ice caps.

• • Anoxic events. These are periods when oxygen is completely depleted from the oceans and so are implicated in mass extinctions. They are believed to be a consequence of changes in climate and ocean circulation.

In this section I will review several studies which claim to have detected periodicities in one or more of the above records. Periodicity detection is often controversial because there are many different ways to analyse time series and there are differences in opinion of what constitutes a significant detection. Yet it is appealing, because detection of a period makes it easier to associate extinction events to a (recurring) physical mechanism. However, periodicities in the data do not automatically imply periodicities in the driver; they may just reflect the characteristic response of the system to an impulsive trigger. More generally, triggers of extinction could be pulsed but the system may respond smoothly. On the other hand, a smoothly varying external driver could give rise to a smooth response in the system (e.g. species origination and extinction), but due to limited sensitivity (e.g. preservation in the fossil record) only those events which rise above some background level are observed. An underlying continuous but variable phenomenon is thus converted into an apparently pulsed one.

With these considerations in mind, let us look at claims of periodicities in the geological data.

Periods in climate data over the past few million years

Atmospheric temperature has been recorded in Antarctic and Greenland ice cores going back as far as 800 kyr BP at high temporal resolution (50–600 yr per interval). Figure 1 shows the temperature variation over the past 420 kyr extracted from the Vostok ice core (Petit et al. Reference Petit1999). We clearly see cyclic behaviour on a timescale of about 100 kyr. Figure 2 shows the temperature extending back even further, using the δ¹⁸O proxy measured in foraminifera. The 100 kyr oscillation period is still seen, but extends only back to about 700 kyr BP. Before that a shorter period clearly dominates; a periodogram analysis identifies it as 41 kyr, with lower power than the 100 kyr one (Muller & MacDonald Reference Muller and MacDonald2000). The relative power of peaks in the periodogram obviously depends on the time window selected.

Fig. 1.

Air temperature as measured in the Vostok ice core over the past 423 kyr from data in Petit et al. (Reference Petit1999).

Fig. 2.

Earth temperature variation over the past 3 Myr. Note that time decreases from left to right (the present is at t=0). Reproduced from Figure 1.6 from Muller & MacDonald (Reference Muller and MacDonald2000) with kind permission of R.A. Muller.

The consensus is that these periodic variations reflect global temperatures and are caused by variations in the solar irradiance (insolation) at the Earth. (Although these variations are obviously both positive and negative, they are often referred to as causing ice ages.) The question is then what causes the variation in the solar irradiance. A priori plausible astronomical mechanisms on timescales of a million years or less are intrinsic solar variability and variations in the eccentricity or inclination of the Earth's orbit about the Sun or in the orientation of the Earth's spin axis. These will be discussed in the subsections ‘Solar variability’ and ‘Variations in the Earth's orbit about the Sun’. Spiral arm and Galactic plane crossings as well as encounters with molecular clouds occur on much longer timescales, so can be ruled out a priori. There is limited data on pre-quaternary (before 2.6 Myr BP) climate, based for example on carbon and boron isotope ratios as proxies for CO₂, and δ¹⁸O and the Mg/Ca ratio in foraminifera as proxies for temperature (IPCC Reference Solomon, Qin, Manning, Chen, Marquis, Averyt, Tignor and Miller2007), Section 6.3. There are no obvious periods in these.

Periodicities in biological variation over the Phanerozoic

Many studies performed since the early 1980s have searched for periodicity in biological extinction or diversity data. These have frequently used the compendium of fossil marine animal genera compiled by the late J.J. Sepkoski.Footnote ¹

The 26 Myr period

Raup & Sepkoski (Reference Raup and Sepkoski1984) generated controversy in the mid 1980s when they claimed that there was a 26 Myr period in the extinction record. They used a database on 3500 marine animal taxonomic families spanning the 253–11 Myr BP (Mesozoic and Cenozoic eras). The data are expressed as the fraction of extinctions in 39 stages, i.e. a time series with 39 points. When presented in this way, these data show 12 maxima (although one is extremely small) which appear to occur quasi-periodically. These stages have a mean length of 6.2 Myr, but because the exact durations are uncertain, they did not normalize the extinction fraction by the duration. This introduces uncertainties into the amplitudes of the time series.Footnote ² They use the times of the 12 peaks as a set of delta functions as their raw data. Based on a Fourier and autocorrelation analysis they identify a significant (p<0.01) period between 27 and 35 Myr. (The data are similar to those from their later study, reproduced here in Fig. 3.)

Fig. 3.

Percentage of extinction among marine animal genera, reproduced from Raup & Sepkoski (Reference Raup and Sepkoski1986) (Figure 1b). Reprinted with permission from the AAAS.

Raup & Sepkoski use a time-domain technique for pulsed events to re-estimate the period. This technique has been used in other studies (e.g. Stothers (Reference Stothers1979), Raup (Reference Raup1985a)) so I describe it briefly. The basic idea is to create a synthetic sequence of events at a certain period and phase and to calculate a goodness-of-fit which measures the standard deviation between each event in the real data and the closest one in the synthetic sequence. This is repeated for a series of phases (e.g. 1, 2, …, P) up to the full period (P), yielding the best-fitting phase and corresponding goodness-of-fit for that period. This is repeated for a range of periods and the best-fitting period/phase identified. This is a type of phase dispersion minimization technique, similar in principle to several others in the literature (e.g. Stellingwerf (Reference Stellingwerf1978), Cincotta et al. (Reference Cincotta, Méndez and Núñez1995)). A measure of significance can be defined using Monte Carlo simulations with random data: from these we could calculate a probability distribution over the goodness-of-fit and so assign a p-value to the best-fit period/phase. Some applications instead ‘normalize’ the goodness-of-fit based on the uniform distribution, using the procedure described by (Stothers Reference Stothers1979). However (Lutz Reference Lutz1985), has shown that this overestimates the significance (see subsection ‘Geomagnetic reversals’).

One advantage of this technique is that because it identifies the nearest point in the synthetic series to each measured event, it is relatively insensitive to missing data. That is, it can identify a period in events even if some are missing. However, one must then provide an explanation for why some events are missing.

Using this technique (Raup & Sepkoski Reference Raup and Sepkoski1984), refine their 27–35 Myr period in the extinction data to 26 Myr. Although the authors acknowledge the significant uncertainty in both the dating and the amplitudes of their data, they ultimately claim this period to be significant. It gives rise to 10 cycles over the 250 Myr time span, such that the phased fit matches reasonably well to only eight of the 12 maxima (see Figure 1 of Raup & Sepkoski (Reference Raup and Sepkoski1984)). Raup & Sepkoski (Reference Raup and Sepkoski1986) revisited these data and revised their conclusions to include just eight significant extinction events, but still with a 26 Myr period (their Figure 1b is reproduced here as Fig. 3). They report a high formal significance for this, p<2×10⁻⁴, although a visual inspection of the data shows that the extinction peaks are not evenly spaced. They reduce their significance to p=0.05 once dating errors are considered. It is important to realise that p is just the probability of obtaining some data or statistic assuming that a particular null hypothesis (e.g. events drawn at random from a uniform distribution) is true. It is not the probability that the null hypothesis is true, nor is 1−p the probability that some particular (and untested) alternative hypothesis (“periodicity”) is true. This will be discussed more in the subsection ‘Hypothesis testing’.

The conclusion of Raup & Sepkoski was criticized and strongly undermined by Stigler & Wagner (Reference Stigler and Wagner1987) in a re-analysis of the data. Using the same time series analysis method, they agreed with the conclusion of Raup & Sepkoski that the null hypothesis (that this time series is random) can be confidently rejected (from Monte Carlo tests they derive a p-value of 0.006±0.001). However, rejection of this null hypothesis does not prove periodicity and Stigler & Wagner demonstrate that non-periodic time series (such as a Moving Average), when analysed with the method of Raup & Sepkoski, reveal significant periods. In other words, the time series method is not specific enough. (A similar conclusion was reached by Kitchell & Pena (Reference Kitchell and Pena1984), discussed in the subsection ‘Complex systems: non-periodicity and multiple mechanisms’.)

It should be noted that the uncertainty in the dating is at least 6 Myr, which is 23% of the 26 Myr period claimed by Raup & Sepkoski. They reasonably argued that dating errors would blur out a real periodic signal, making it even harder to find, so that this cannot be used as an argument against the significance of a detection. But the situation is not quite so simple, because uncertainties are not necessarily random. An important phenomenon is the ‘Signor–Lipps’ effect: because fossilization is rare, the last appearance of a fossil in a stratigraphic layer predates the final extinction, so a species/genus event may appear to go extinct in an earlier stage than it really did. (The effect is further enhanced by dating errors.) Stigler & Wagner showed that on account of the variable stage length in this data set, it is prone to showing an apparent period of around 26 Myr when analysed with this time series analysis method. (The value of 26 Myr is a consequence of the typical duration of the stages.) Raup & Sepkoski (Reference Raup and Sepkoski1988) responded to this, but in their reply Stigler & Wagner (Reference Stigler and Wagner1988) stuck to and extended their original criticism.

It is worth asking what is the probability of detecting a period at all in the presence of dating errors? This was examined by Heisler & Tremaine (Reference Heisler and Tremaine1989). They used a Monte Carlo technique to simulate time series with dating errors from a nominal time series with eight equally spaced points, emulating the eight extinction events of Raup & Sepkoski (Reference Raup and Sepkoski1986). In each time series, points are randomly jittered by an amount chosen at random from a Gaussian distribution with zero mean and standard deviation Δ. Heisler & Tremaine showed that if Δ is more than 13% of the period (4.6 Myr for a 26 Myr period) then the period cannot be detected with a confidence above 0.90. If Δ is increased to 23% (6 Myr, the average dating error) then the probability of detecting a period drops to 0.55. Heisler & Tremaine therefore conclude that the period detected by Raup & Sepkoski is either due to a statistical fluke or to a biasing of the data.Footnote ³ This conclusion seems unfounded, because a 50:50 chance of detecting the period is reasonable odds. Interestingly, if the time errors are only 7% of the period, then the period can still be detected with 100% confidence (Figures 1 and 2 of Heisler & Tremaine (Reference Heisler and Tremaine1989)). This might lead one to speculate that the dating uncertainties in the Raup & Sepkoski data have been overestimated, thus reconciling the two analyses. However, this is not the case, because the dating uncertainties are intrinsic to the placing of the extinctions in broad and variable-length (average of 30 Myr long) stages. If the actual extinctions were not distributed at random within these, then this would be a case of a fluke which Heisler & Tremaine refer to.

As we shall discuss in subsection ‘Reconciling a period detection with its non-detectability’, neither the Raup & Sepkoski rejection of a null hypothesis of randomness, nor the Heisler & Tremaine argument that a period could not be confidently detected, are sufficient to draw a conclusion about the existence of a period.

A period of around 26 Myr has been identified by other authors. Rampino & Caldeira (Reference Rampino and Caldeira1992) compiled data on seven different types of geological events including mass extinctions, orogenic events, sea-floor spreading and flood basalt volcanism. The data are expressed as 80 delta functions spanning the past 260 Myr. To create a contiguous data set more convenient for power spectral analysis they first smooth these data using a moving average (both Gaussian and top-hat functions are tried). Power spectra are calculated from the Fourier transform of the autocorrelation function using a Tukey window. They identify several significant periods when analysing just the individual data sets, but a peak at 26 Myr stands out when they combine all seven sets (as well as in some individual sets). The significance is tested by seeing how many out of 1000 pseudo-random time series give a power equal to or greater than that detected at 26 Myr. They claim that these tests reject the null hypothesis at p=1×10⁻⁴, although how such a small p can be achieved from just 1000 tests is unclear (as it would imply 0.1 tests showed an equal or higher power). Moreover, the probability that the random data show a peak with this power or more at any period over 10–65 Myr is 0.045, and this is surely the more relevant figure.

What is interesting about this study is that they combine data on different types of event. But are the data independent? After all, geological changes are used to define the boundaries in the geological timescale, which in turn is used to date the events. Rampino & Caldeira nonetheless suggest that the periodicity is real and could reflect regularities in core/mantle dynamics (‘pulsation tectonics’). Noting also that a similar period has been claimed for the impact cratering record, they suggest that the two may not be independent, i.e. impacts could cause geological change. Napier (Reference Napier1988) performed a similar study, combining data on impact craters, mass extinctions (from Raup & Sepkoski (Reference Raup and Sepkoski1984)) and the geological events from Rampino & Caldeira (Reference Rampino and Caldeira1992). He identifies a period of 27±1 Myr and links this to asteroid/comet impacts (or ‘bombardment episodes’) driven by Galactic plane passages.

The 62 Myr period

Rohde & Muller (Reference Rohde and Muller2005) examined the variation of the biological diversity in the Sepkoski (Reference Sepkoski, Jablonski and Foote2002) database, recalibrated according to the timescale published by the International Commission on Stratigraphy (ICS 2004, see Gradstein et al. (Reference Gradstein, Ogg and Smith2005)). The data – both the complete data set and one pruned of single-occurrence genera and ‘poorly dated’ genera – suggest a periodicity which is visually stronger after the removal of a slowly varying component with a third-order polynomial (see Fig. 4). A power spectral analysis reveals a strong peak at 62±3 Myr and a weaker peak at 140±15 Myr.

Fig. 4.

Figure 1 from Rohde & Muller (Reference Rohde and Muller2005), showing the variation in genera diversity before (a) and after (b) removing poor data. (c) shows (b) after removal of a third-order polynomial trend (with the 62 Myr period sinusoid overplotted). (d) shows the same having also subtracted the 62 Myr fit (with a 140 Myr period sinusoid overplotted). (e) is the power spectrum showing two alternative backgrounds (W and R) for the significance calculation. Note that time decreases from left to right (the present is at t=0). Reprinted by permission from Macmillan Publishers Ltd (Nature 2005).

When one talks about the ‘significance’ of the detection of an event, such as a peak in a power spectrum, one is normally concerned with the probability that the event could have been caused by a ‘random’ process. This is typically the noise associated with some background (what we would measure in the absence of a signal). However, often the noise process is not well understood, so there is no unique definition of ‘random’: we could draw events at random from a uniform or a Gaussian distribution, for example. In classical hypothesis testing, one tries to reject the null hypothesis that the signal was produced by some specified background.

In the case of power spectral analysis it is usual to calculate the significance of peaks in a power spectrum by comparing it to the power spectrum of the background process. Rohde & Muller use two different backgrounds: a complete randomization of the detrended data (R); and a randomization in consecutive stages (of length 27 Myr), which preserves the short-term correlations (W); see Figure 4. W reflects the fact that if the data are periodic then nearby points are not independent. Both are examples of ‘red noise’, in which the background is higher at lower frequencies. W of course has more power than R at high frequencies (and less at low), with the choice of stage length to use for W affecting its power spectrum (in the limit of small stages it will look like the R power spectrum; this illustrates that there is no simple definition of ‘random background’ as a null hypothesis). The probability of these achieving a power as high as that seen in the 62 Myr period anywhere in the spectrum is quoted as <0.0013 for R and 0.01 for W. These decrease to around 1×10⁻⁴ when we only consider the probability of getting a peak at 62 Myr; that is, if we suspect a peak at this period (before we look at the data) and examine the evidence for a peak only at this period. This might be valid if one has independent evidence for a peak at this period and if one ignores any other peaks. However, there is no prior reason to suspect a peak at 62 Myr, so the higher – yet still significant – p-values should be used. The 140 Myr peak is only significant at p<0.01 if one limits the significance test to this period (otherwise p=0.71 and 0.13 for R and W, respectively), and then only for background model W. We should not regard this as a significant detection.

The plot of the fitted 62 Myr sinusoidal curve to the data (Fig. 4) shows reasonable agreement in period and phase with the detrended data for much of the Phanerozoic, although the fit is quite poor over the past 150 Myr. This is not inconsistent with finding a significant period, because even if this is a statistically significant detection, it does not automatically follow that it explains a significant part of the total variance in the data. All we have found is that it describes significantly more variance than the background process, and more than any other single period. The analysis has not ruled out that some other process may explain as much or even more of the variance in the data. Indeed, we see in Fig. 4b that the variation due to the long-term trend swamps the periodic variations, and may reflect a far more important aspect of biodiversity change.

Rohde & Muller stress that the 62 Myr period is only found when using the ICS 2004 timescale. They also ask whether this is a period in the true diversity or just in the observed diversity. One subsequent study (Smith & McGowan Reference Smith and McGowan2005) has suggested that much of the signal is a selection effect (see the subsection ‘How reliable are the data, how appropriate are the methods?’). Incidentally, Rohde & Muller were not the first to suggest a 62 Myr period. Thirty years earlier, Thomson (Reference Thomson1976) suggested there was such a period in the number of genera in various groups of lower vertebrates and certain invertebrates, although this very short article presents just a superficial analysis.

Methodologically the Rohde & Muller result has been criticized by Omerbashich (Reference Omerbashich2006) on three accounts: (1) interpolation of the time series to give uniform spacing (increasing the original 167 data points to 2170 points; this is necessary for using the fast Fourier transform, FFT); (2) the polynomial detrending prior to taking the power spectrum; (3) the use of zero padding in the Fourier transform. Interpolation can cause problems, e.g. by introducing short-scale correlations in the data, but is acceptable as long as one realizes that it corresponds to a smoothing at high frequencies. Detrending is a standard procedure and seems valid, provided one recognizes that it assumes that whatever causes the long-term trend is independent of whatever causes the periodic variations. Zero-padding the time series to lengthen the time span is a common technique with the FFT to reach the maximum possible resolution of the power spectrum. Note that this does not (cannot) add information, so one should estimate the maximum resolution supported by the data and take care not to overinterpret the results.

Omerbashich reanalysed the original 167 points from Rohde & Muller using the Lomb–Scargle (LS) periodogram (Lomb Reference Lomb1976; Scargle Reference Scargle1982) and without applying any detrending. (LS works with unevenly sampled data and does not need zero padding. It goes under various names, including Gauss–Vaníček spectral analysis.) He found the 62 Myr period to be insignificant at the p=0.01 level, and instead reports significant periods at 194, 140 and (according to his Figure 2a) 77 Myr. However, this difference appears to be the result of not detrending, rather than the data preprocessing and spectral analysis technique used: Cornette (Reference Cornette2007) found that a LS analysis of the original 167 data points, detrended, still gives rise to a significant period at 62 Myr. In contrast, Cornette found no period to be significant via either LS or FFT if the data are not detrended. This suggests that the detrending increases the sensitivity to the periods reported. Cornette also concludes that the additional 140 Myr period reported by Rohde & Muller is not statistically significant. (This is not surprising, because this period is a quarter of the total time window, and thus sensitive to the exact detrending applied.)

Lieberman & Melott (Reference Lieberman and Melott2007) also reanalysed the Rohde & Muller biodiversity data and also detect the 62 Myr period as significant (see also Lieberman & Melott (Reference Lieberman and Melott2009)). They further analysed selected time windows of the biodiversity data: while the 62 Myr periodicity is strong over the period 520–150 Myr BP, there is no evidence for a significant period over the last 150 Myr, which, as mentioned above, is not surprising when one inspects the time series fit to the data (Fig. 4). Thus any astronomical mechanism invoked to explain the 62 Myr period must also explain why it did not significantly influence biodiversity in the past 150 Myr. This could imply there are multiple mechanisms at work. Lieberman & Melott (Reference Lieberman and Melott2007) also examine extinction and origination data and find a significant period at 27 Myr, although they caution that this might be an artefact of the data. Melott & Bambach (in preparation) also find a 62 Myr period in other geological and biodiversity measures (private communication).

Many of the published articles on periods (26 Myr or 62 Myr) in biological diversity are based on the same database or updated versions thereof by Sepkoski. Melott (Reference Melott2008) therefore looked for periods in the Paleobiology Database, a more recent compilation from many sources where particular attention has been made to correct for the preservation bias.Footnote ⁴ He finds a period of 63 Myr (using correlation analysis, FFT and LS), consistent with the 62±3 Myr period found by Rohde & Muller (Reference Rohde and Muller2005), further checking that the two had the same phase (to within 1.6 Myr, i.e. within the measurement uncertainty). Two other statistically significant peaks are disregarded as artefacts of detrending and variations in the interval length. There is no figure showing how the fit corresponds to the original data, so it is unclear whether there is a good fit over the whole Phanerozoic. Alroy (Reference Alroy2008) looked for periodicity on data selected from the same database. He finds no significant autocorrelation and says that the power spectrum (Figure S2 in the online supplement) is consistent with white noise: there is no significant power around 26 Myr, although he suggests weak evidence for a significant period ‘somewhat longer than 62 Myr’. (The definition of ‘significance’ depends on a 95% confidence band around the median power, the origin of which is not clear.)

In summary, there have been claims of a periods in the fossil record of 26 Myr spanning 250–10 Myr and of 62 Myr spanning 520–150 Myr (formal uncertainties in the periods are 1–3 Myr). The former (and older) claim has been heavily criticized on grounds of data and methodology and I believe now has little credibility (more on this in the subsection ‘Issues in time series analysis’). The 62 Myr period is on a stronger footing, although less time has passed for this to be re-analyzed/criticized and many of the issues I will discuss later apply.

Periodicities in the geological record over the Phanerozoic

Complex systems: non-periodicity and multiple mechanisms

Although there is a suggestion of periodicity in the biological and geological records, a single period is often not a good fit, at least not over long timescales. We saw this already for climate data over the past 3 Myr (Fig. 2), where there is a clear change in the period. Attempting to derive a single period persistent over the whole time series may give a meaningless result, depending on the method used. Rampino (Reference Rampino1998) has suggested that a better periodic fit to extinction and cratering data is obtained by allowing a variable period. Alternatively one could allow phase shifts with a constant period. Of course, both of these would require a causal mechanism. Later we will see that the orbit of the Sun about the Galaxy naturally accommodates a variable period of its vertical oscillations through the Galactic disk, either due to spiral arms or the R-dependence of the potential.

There are many studies in the literature which conclude that biodiversity does not vary periodically. In his review, Bambach (Reference Bambach2006) argues that in the Sepkoski genus-level database there are only three events which have a level of extinction significantly higher than the background, thus challenging the claim that mass extinctions are common enough to even attempt to describe them as periodic. (Some have made claim of a ‘big five’ of mass extinctions, others just a single ‘big one’. This depends on how we define ‘big’ or ‘significant’.) Based on the diverse impact of lower intensity extinctions on biological diversity, he concludes that these are unlikely to have a single cause. This is an important point, because if several mechanisms are at work, then even if one or more of these are periodic, attempts to derive a single period will give spurious results. Bambach also concludes that the variability and the peaks in the extinction rate are not a result of incompleteness in the fossil record.

Many authors have argued for the existence of multiple extinction mechanisms. White & Saunder (Reference White and Saunder2005) note that large impacts and mass volcanism have occurred much more frequently in the past 300 Myr than have mass extinctions, so neither of these in isolation could have caused the mass extinctions. They present a statistical analysis which shows that a few random coincidences of these two mechanisms is not unlikely, and could explain the frequency of mass extinctions. In support of this, Hallam (Reference Hallam2004) notes that there are many large impact craters for which no contemporary mass extinction has been identified. Arens & West (Reference Arens and West2008) likewise argue that a simultaneity of volcanism and impacts (or, more generally, what they call ‘press’ and ‘pulse’ events) is necessary to explain mass extinctions.

It is well known that apparent periodicities in data can be produced by non-periodic processes. Chaotic phenomena, for example, can exhibit quasi-periodic signals. It is not even obvious that extraordinary events such as mass extinctions require an external driver (whether extraterrestrial or not, or periodic or not). From an analysis of the correlation between biodiversity and extinction, in particular the observation that a high rate of extinction tends to follow times of large biodiversity, Alroy (Reference Alroy2008) concludes that at least some variation in the fossil record can be explained by purely ecological interactions (e.g. predation and competition). However, this diversity–extinction correlation is weak (ρ=0.44). Numerical models have also shown how, in such interdependent systems, small perturbations can produce large extinctions and large perturbations may result only in small extinctions (e.g. Plotnick & McKinney (Reference Plotnick and McKinney1993)). The implication is that we should not necessarily look for a linear correlation between cause and effect. It has further been suggested that evolution may be a self-organized critical phenomenon arising from the interaction between species (e.g. Kauffman & Johnsen (Reference Kauffman and Johnsen1991)). Such dynamical systems can show highly nonlinear responses, with extreme events an almost inevitable consequence of the system dynamics. Evidence supporting this idea is the power-law distribution of the magnitude of extinction events and taxa lifetimes (Newman Reference Newman1997), although the data are noisy and other fits are possible. Stanley (Reference Stanley1990) has suggested that apparent periodicities in biodiversity are a result of delayed ecosystem recovery from mass extinction and thus a feature of the system rather than any driver.

Kitchell & Pena (Reference Kitchell and Pena1984) fit various time-domain models to the extinction data of Raup & Sepkoski (Reference Raup and Sepkoski1984) and found that a stochastic autoregressive (AR) model provides a good fit (just as Stigler & Wagner (Reference Stigler and Wagner1987) found that a moving average process is a good model for these data; see the subsection ‘The 26 Myr period’). In this AR model, each value of the time series is a linear combination of the previous five events. The system retains a memory and can display pseudo-periodic behaviour, even though there is no external driving force. Kitchell & Pena find that the AR model fits the data better than either a continuous periodic variation or a model with periodic impulses. Indeed, they claim that the lack of strict periodicity and the large variations in extinction magnitude rule out a deterministic explanation of the data, which in turn suggests that biodiversity is a dynamic phenomenon. Of course, just because the AR model fits the data well, this does not “rule out” the periodic model – we must compare the posterior probabilities of the models – and does not rule out the presence of an external driver. Yet the onus is on us to show that the periodic model is considerably more likely, and the good AR fit argues a single external driver being the sole or even dominant cause.

How reliable are the data, how appropriate are the methods?

Geological data are far from perfect. One of the most significant issues is that of dating. There exist different dating scales and one of the more robust results reported – the 62 Myr periodicity of Rohde & Muller (Reference Rohde and Muller2005) (subsection ‘The 62 Myr period’) – was only found when the Sepkoski (Reference Sepkoski, Jablonski and Foote2002) data were recalibrated on the ICS 2004 timescale. There may be good reasons to assume that this is the most accurate timescale to date, but that was presumably thought of earlier timescales too. In addition to such calibration uncertainties, there are significant random errors in dating extinctions or impact craters, which can extend to tens of Myr.

Beyond dating, there is the problem of the incompleteness of the geological record. Concerning the impact record, presumably not all craters have been found, some have eroded or been subducted into the Earth's interior, and the amount of land surface available for impacts has varied over time. Similarly, the efficiency of fossilization depends on environment and species, thus giving rise to biases. Fossilization is anyway rare, and as we already discussed in the subsection ‘The 26 Myr period’ it gives rise to the Signor–Lipps effect. In addition to possibly producing spurious periods under some circumstances (as already discussed), this effect tends to smooth peaks in the extinction record, making extinctions appear more gradual than they really were (e.g. Hallam (Reference Hallam2004), Chapter 3). Fossils have also not been searched for uniformly over the globe, introducing a geographical bias (compounded by the fact that the continental plates have migrated with time).

One must also ask how reliably the recorded variable measures the phenomenon of interest. By measuring the size of an extinction through the number of genera which survive, we put genera in which almost all species died on an equal footing with genera which were unaffected. A different measure will be obtained if families or the total biodiversity are used instead. All of this is measured via fossils, but to what extent do the data represent variations in biological life rather than variations in their preservation in the fossil record? If we are open to the idea of a mechanism which has a widespread impact on life, then we should be equally open to the idea of a mechanism which can affect preservation. Peters & Foote (Reference Peters and Foote2002) showed that much of the observed variability in the marine fossil record can be explained by variations in the amount of rock exposed for fossil preservation, rather than variations in the extinction rate itself (see also Smith (Reference Smith2007)). Smith & McGowan (Reference Smith and McGowan2005) give evidence suggesting that this is the source of the biodiversity variation reported by Rohde & Muller (Reference Rohde and Muller2005) (see the subsection ‘The 62 Myr period’). They note that long-lived taxa, which are less affected by this sampling bias, do not show periodic variations. Rock exposure variations could be due to changes in sea level, themselves plausibly driven by plate tectonics or climate change (locking ice into glaciers and ice caps). Of course, one must still explain why short-lived taxa and/or sea levels appear to vary periodically and whether preservation issues interfere with periods on all timescales.

Our variable of interest is often measured indirectly via a proxy, which must be calibrated. With ice cores, for example, one measures δ¹⁸O as a function of ice depth, yet δ¹⁸O must be converted to temperature (using calibration terms which depend on geographical location) and depth must be converted to time. Our proxies/samples may also be contaminated. For example, Royer et al. (Reference Royer, Berner, Montañez, Tabor and Beering2004) showed that the δ¹⁸O measured in marine carbonates has to be corrected for water pH. Patterson & Smith (Reference Patterson and Smith1987) claimed that 75% of the extinctions in families reported by Raup & Sepkoski (Reference Raup and Sepkoski1984, Reference Raup and Sepkoski1986) are artificial extinctions introduced by taxonomic definitions. This includes ‘pseudo-extinctions’ of species which do not die out but rather evolve into something else.

In terms of the techniques, there are many ways of analysing a time series. Many appear appropriate, but they can produce different results. There are decisions of how to calculate the significance, what to accept as significant (and so report anything at all), what filtering to use, what data to omit, etc. Several authors (e.g. Grieve et al. (Reference Grieve, Sharpton, Rupert and Goodacre1988), Napier (Reference Napier1988), Lieberman & Melott (Reference Lieberman and Melott2007), Melott (Reference Melott2008)) have found statistically significant periods in data which they then argue away (perhaps quite reasonably) on other grounds.

Is our choice of data appropriate? If extinction data are essentially a contiguous time series, does it make sense to focus on selected large extinction events and try to explain just these? This would rule out testing mechanisms which are capable of explaining both low and high amplitude events. For example, Alroy (Reference Alroy2008) find that the distribution of extinction and origination rates of marine invertebrates over the whole Phanerozoic can be fit well with a log normal (after detrending). This does not necessarily imply a single common mechanism. But do we need to invoke occasional catastrophic triggers if in fact large amplitude events can be explained as the tail of a more mundane process?

This summary of the issues is not intended as blanket scepticism, but rather as a reminder that there are choices to be made in data analysis. It is not as objective as it sometimes appears to be. Time series analysis is a complicated business: there are decisions of which data to include, what preprocessing to do, which methods to use and which significance tests to apply. We should question the assumptions and identify the uncertainties in order to examine what impact they have on potentially far-reaching conclusions.

How independent are the studies?

Numerous articles have identified supposedly significant periodicities with a period in the range 25–33 Myr. This has led many authors to speculate an astrophysical cause, partly because of the lack of a plausible terrestrial mechanism for such a periodicity. Despite various criticisms, the sheer number of studies converging on a similar values is noteworthy. But are they independent? First, many studies use the same data sets and dating system, so are subject to the same systematic errors and sampling biases. Second, some studies use the same methods, some of which have been demonstrated to have deficiencies. Third, what is taken as ‘evidence’ or ‘significance’ is often inadequate. Fourth, many other studies find different periods or no periods at all.

A phenomenon worth noting is the ‘band wagon’ effect, in which the presence of a published value biases authors' analyses (consciously or otherwise) and the conclusions they choose to publish towards confirming that published value. Studies with results lying far from the current trend may not be published at all. In the current context, the presence of a mechanism which has a period near to one which the data could support may play a similar role. The band wagon effect is described in the context of distance measures to the Large Magellanic Cloud by Schaefer (Reference Schaefer2008), who notes that the estimates are far more consistent with each other than expected based on the reported uncertainties in the individual estimates.

Hypothesis testing

The studies discussed above are concerned with assessing evidence for periodicity in time series data. This is an example of hypothesis testing. Probably the most common approach to this now in use is that developed by Fisher (Reference Fisher1925). The general idea is to define a null hypothesis, a model for producing the data in the absence of the effect one is investigating. We then calculate the probability that the data (or rather a statistic based on them) are predicted by this model. If this probability is low, we ‘reject’ the null hypothesis at this confidence level, which suggests that some alternative hypothesis may be more likely to explain the data. As a simple example, imagine we have two groups of people, one which has received special training, the other not, and we wish to assess the impact of the training on scores in an exam. The typical null hypothesis states that the training makes no difference. If we assume that the scores of the individuals in each group are distributed according to a Gaussian, we would perform a t-test to examine whether the means of the two groups differ by a significant amount, where ‘significant’ is referred to the pooled standard deviation of the two groups. The result is the ‘p-value’ or P(D|H ₀), the probability of observing the data we did, D, given that the null hypothesis, H ₀, is true. (This is sometimes called the ‘likelihood’ of the hypothesis.)

In this simple example the null hypothesis (H ₀; means equal) is the complement of an alternative (H _A; means unequal) which is of interest. The two hypotheses cover all possibilities, so the alternative hypothesis is implicit in the definition of the null. This is not generally the case however. In the case of assessing the significance of peaks in a power spectrum of a (zero mean) time series (e.g. the biodiversity data shown in Fig. 4), we could define a null hypothesis as a time series with the same sampling drawn from a Gaussian distribution with mean zero and standard deviation, σ, set equal to the standard deviation in the original data. We can then measure its power spectrum (or rather the probability distribution of the power at any period via Monte Carlo simulations) and calculate the p-value of this null hypothesis. If p is small (typically we require p<0.01) we ‘reject’ the null hypothesis. Many people automatically assume that this therefore ‘accepts’ the alternative hypothesis. But this is not the case, because the alternative here is not the complement of the null. Indeed, we have not even specified (let alone tested) the alternative! All we have done is assign a low probability to one specific null hypothesis. There might be other null hypotheses which predict the data with a higher probability. There are certainly other ways to specify the null (a different value of σ, non-Gaussian noise, retain short-term correlations, etc).

Ideally we try and reject several different null hypotheses, as Rohde & Muller did with their W and R background models. However, the model of interest (periodic variability) is not tested in orthodox hypothesis testing. The null hypothesis may be constructed using some properties of the measured data (e.g. the time sampling), but the data themselves are not tested.

Orthodox hypothesis testing has the further curiosity that it assesses the probability of getting data that are not observed. In the t-test example above, we do not actually calculate the probability of getting the data under the null hypothesis: the probability of observing a specific value from a continuous distribution is infinitesimally small. Instead, we calculate the probability that the means differ by the measured amount or more. Likewise with the power spectra: we would calculate (via Monte Carlo simulations) the probability that the null hypothesis can produce the observed power or more. But why should we be interested in the probability of observing data we never actually saw? This issue goes to the heart of criticisms of hypothesis testing (and the limitation of ‘proof by contradiction’ or ‘falsifiability’), which have been discussed extensively elsewhere (e.g. Berger (Reference Bergern2003), Jaynes (Reference Jaynes2003), Christensen (Reference Christensen2005), Section 1.4 of Sober (Reference Sober2008)).

A yet more important point about hypothesis testing is the interpretation of the p-value. If we calculate a low p-value for a null hypothesis H ₀, then we have found that P(D|H ₀) is low. This does not mean that P(H _A|D) is high! (Some articles nonetheless interpret 1−p as the probability that the alternative is true. This is wrong.) Strictly we cannot even ‘reject’ H ₀; for this we would need to know P(H ₀|D) which is not the same as P(D|H ₀). We can illustrate this with a simple example. Imagine you draw a card at random from a deck of cards and it is the ace of spades. This is D. The probability of drawing this, assuming it is a normal deck (H ₀), is P(D|H ₀)=1/52. But you are unlikely to tell me that the probability that the deck of cards is normal [P(H ₀|D)] is 1/52. Or take a more extreme example: a particular person wins the lottery with a chance of 1 in 10⁸. Would you claim that the probability that the lottery is fair is 10⁻⁸? We have a similar problem with scientific data: the probability of getting the data we actually observe is very small under almost any hypothesis (vanishingly so with continuous variables).Footnote ⁵ The point is that even a low P(D|H ₀) or p-value may provide more support for the null hypothesis than for any other alternative. We cannot know how low p should be in order to reject the null hypothesis (lottery is fair) without knowing the p-value of alternative hypotheses (lottery is rigged, lottery is fair but someone bought all of the tickets, etc). If a data point lies 10σ from the mean of a Gaussian, we can only say the data did not come from that Gaussian if we accept there is an alternative origin. How likely we think there to be an alternative is quantified using the prior probability of the model. This is the probability that the model is true independent of (before using) the specific data D (Sivia Reference Sivia1996; Jaynes Reference Jaynes2003; Gregory Reference Gregory2005).Footnote ⁶ Only if the prior probability of the alternative is very small might this measurement give evidence in favour of the Gaussian origin.

The solution to this dilemma is to do direct model comparison, that is, to compare the likelihoods, P(D|H _i), for two explicit models (their ratio is called the ‘odds ratio’). If we give the two models equal prior probabilities, P(H _i), then the one with the largest likelihood better predicts the data. Yet this still does not give the posterior probability, P(H _M|D), for the hypothesis of interest (H _M, e.g. periodicity at some period). We can only calculate this if we know the complete set of alternative models, H _i. We can then use Bayes' theorem (which follows from the basic axioms of probability)

(1)

The Bayesian approach to hypothesis testing is to explicitly test (calculate P(D|H _i)) all plausible models for the data, including the model of interest. Only in this way can we calculate the quantity we are actually interested in, P(H _M|D). This avoids having to calculate the probability of observing data we might have seen but did not. Although this approach overcomes the limitations of orthodox hypothesis testing, it presents a new problem, namely the need to specify all plausible alternative models (all those which do not have a vanishingly small prior probability). In most real-world problems it is almost impossible to define all plausible alternatives. That is, the model space is incomplete. But at least this encourages us to define and test as many plausible alternatives as we can think of.Footnote ⁷

The fundamental difference between orthodox hypothesis testing and the Bayesian approach is that the former tests the probability of a single null hypothesis whereas the latter always compares the probabilities of two or more alternative hypotheses. We can only ‘reject’ a hypothesis if a better candidate is available: the unlikeliness of the data itself is not enough. This limitation of orthodox hypothesis testing I shall refer to as incomplete inference. Rejecting a null hypothesis may be a useful first step, but only if it fills a large part of the hypothesis space (Fig. 5). If, in practice, we interpret a low p-value as evidence against the null, then it is because we implicitly assume the alternatives to have low prior probabilities this idea carries over to time series analysis. We can find the power spectrum of any data set (the Fourier transform is just a basis function projection). However, the peaks are only interesting if they cannot be produced by any other non-periodic model which is plausibly responsible for the data. There are methods of significance estimation for time series analysis which do not depend on p-values for the probability of unobserved data (e.g. Gregory & Loredo (Reference Gregory and Loredo1992), Sturrock (Reference Sturrock2008)). These typically result in lower significance estimates than do p-value estimates.

As an aside, there are of course many examples in science where doubt has been cast on a theory on the basis of a poor fit to data, without their being a concrete alternative defined at the time. Examples include the perihelion precession of Mercury (which couldn't be satisfactorily explained by Newtonian mechanics), or the ultraviolet spectra of black bodies (which were incorrectly predicted by classical physics). In essence, a low p has sometimes been essential for motivating the search for alternative hypotheses, even though a formal approach to deciding how low p should be before searching for alternatives (rather than repeating the experiment or reassessing the measurement errors) presumably has been rarely adopted. Nonetheless, a proper statistical approach is still required in order then to compare the alternatives.

Fig. 5.

Schematic illustration of hypothesis testing. The hatched box shows the full hypothesis space. The left-hand (blue) box indicates the space covered by a certain periodic hypothesis and the two right-hand (red) boxes represent two alternative hypotheses for random models. In orthodox hypothesis testing for periodicity, one of the random model hypotheses will be considered unlikely if P(D|H _random1) (the p-value) is low, even though this does not imply P(H _random1|D) is low. Moreover, this alone cannot imply that P(D|H _periodic) is high, because generally there are other untested hypotheses such as H_random2.

Reconciling a period detection with its non-detectability

In the subsection ‘The 26 Myr period’ I discussed the detection of a significant periodicity in extinction data by Raup & Sepkoski (Reference Raup and Sepkoski1984, Reference Raup and Sepkoski1986) plus the criticism by Heisler & Tremaine (Reference Heisler and Tremaine1989) that, given the large dating errors, the probability of detecting a true period in these data is very small. How can we reconcile these two claims? It turns out that the two articles are testing different hypotheses. Let D represent the detection of a period in some time series data by some method. Raup & Sepkoski calculate the probability, P(D|H _r), of detecting the period under the null hypothesis that the data were generated by some random process, H _r. Because this probability is low they reject H _r. Heisler & Tremaine, by examining the recovery rate of a period in simulated time series with timing errors, evaluate P(D|H _p), where H _p is the hypothesis that the process which produced the data is periodic. Because it is low they say one is unlikely to detect the period, hence the apparent period must be caused by something else (e.g. noise). Yet neither study is sufficient to decide whether the data give evidence for H _p, because neither calculate P(H _p|D).

H_r is defined by the Monte Carlo method with which Raup & Sepkoski generate random time series, from which they calculate their p-value. They use the low value of p to (1) reject H_r and (2) infer that H_p is true. As described in the previous section, these two inferences make the additional assumption that H _r is the only alternative to H _p. This is not true, because we could have defined other random processes to calculate the p-value. If we are nonetheless generous to Raup & Sepkoski and assume that H _r is the only alternative to H _p then we can write (the horizontal bar means ‘not’). The two hypotheses under consideration are then related by

(2)

where P(D) is the probability of detecting the period at all (under either hypothesis). We are interested in . This is related to the above quantities via Bayes' theorem

(3)

Substituting Eqn 2 into this gives

(4)

If we have no reason to prefer H _p over , then we could set their prior probabilities to be equal. They then cancel out of the above equation leaving

(5)

Thus, in order to decide whether the data, D, favour H _p over , we must examine the odds ratio, . H _p is favoured if this ratio is less than 1, and Eqn 5 gives us a measure of confidence in this in terms of an actual probability. This is the standard Bayesian approach to testing two hypotheses with equal priors (e.g. Gregory (Reference Gregory2005)). In contrast, working with only or only P(D|H _r) are examples of what I called incomplete inference.

What does this calculation reveal for the analyses of Raup & Sepkoski and Heisler & Tremaine? It is not easy to convert their quoted confidence levels into consistent probabilities because they use different methods for detecting periods and for assigning significances. Furthermore, Raup & Sepkoski (Reference Raup and Sepkoski1986) look for evidence of periodicity at a range of periods, for which there is a higher probability of detecting something significant than at a specific period. So what follows is only approximate.

Raup & Sepkoski (Reference Raup and Sepkoski1986) identify the 26 Myr period because it crosses a 99.9% confidence level in their non-parametric test, implying p=0.001. The confidence for a period at any period must be lower (they are implicitly using prior information in identifying 26 Myr with the stated confidence). They also describe Monte Carlo simulations of 500 random time series (from H _r), of which 23 gave rise to a significant (p<0.001) period at some value. So really lies somewhere between 0.001 and 23/500=0.046 depending on what we are testing.

On the other hand, the simulations by Heisler & Tremaine (Reference Heisler and Tremaine1989) deal with the detectability of a periodic signal at a single frequency. They find that (only) 55% of their error-perturbed time series result in this period being detected, implying that P(D|H _p)=0.55.

Putting these two values, and P(D|H _p)=0.55, into Eqn 5 gives P(H _p|D)=0.998. On this basis, the strong evidence for a period at 26 Myr outweighs the somewhat modest probability of detecting the period at all. Had we used instead of 0.001 (which is more realistic, as we are interested in any period), we would get P(H _p|D)=0.92, which is still in favour of a period but much less confident. More significantly, there are neglected hypotheses (H _p is not the only alternative to H _r) so there should be additional terms on the right-hand side of Eqn 2 and thus in the denominator of Eqn 5. As these terms are always positive they would decrease P(H _p|D). Neglecting hypotheses leads us to overestimate the confidence in H _p. As discussed in the subsection ‘The 26 Myr period’, Stigler & Wagner (Reference Stigler and Wagner1987) showed that non-periodic models do produce a significant peak with the analysis method of Raup & Sepkoski, so important alternative models have been neglected and P(H _p|D) is certainly overestimated here.

There is another important issue, namely that of the priors. We have so far assumed (in Eqn 4) that the two hypotheses being tested have equal prior probabilities. That is, we have assumed that the unconditional probability of getting a period at 26 Myr (plus/minus some bin width) is 0.5. But why would we want to assume a large probability for a specific period before we have even seen the data? It seems more reasonable to give equal priors to the hypotheses ‘periodic’ (for any period) and ‘non-periodic’. This must lower the prior for a specific period by a large amount, thus reducing the posterior probability for that period. Ultimately, the Raup & Sepkoski analysis is insufficient to provide significant evidence for H _p.

There are other subtle issues, and at some point the quality and quantity of the data may not justify a much more thorough analysis. The main point of this discussion was to convince the reader that assessing evidence for periodicity is not a trivial matter, and is far from being concluded by a low p-value for some specific null hypothesis.

Minor body impacts, Oort cloud perturbation and comet capture

There are hundreds of large impact craters on the Earth (Shoemaker Reference Shoemaker1983).Footnote ⁸ It is widely accepted that the mass extinction 65 Myr ago at the K-T boundary was caused at least in part by an asteroid impact, as evidenced by a global iridium deposit (Alvarez et al. Reference Alvarez, Alvarez, Asaro and Michel1980) and the identification of the Chicxulub crater in Yucatan, Mexico (Hildebrand et al. Reference Hildebrand, Penfield, Kring, Pilkington, Camargo, Jacobsen and Boynton1991). Some contest this view, and it is possible that some of the extinction was caused by major volcanism which occurred at the time in central India (the Deccan traps). Evidence put forward in support of this is an apparent increase in extinction prior to the impact. Indeed, there has been a significant debate between proponents of volcanic and impact causes of mass extinctions Glen (Reference Glen1994).Footnote ⁹ Both Hallam (Reference Hallam2004) and Alvarez (Reference Alvarez2003) argue that while there is strong evidence for a giant impact having caused the K-T extinction (shocked quartz, tektites, iridium, an appropriate crater), such evidence is lacking for other mass extinctions.

These considerations aside, it is clear that impacts of large asteroids or comets have occurred many times and can cause widespread extinction. The mechanism is either the violence of the impact itself (blast, fires, earthquakes, tsunamis) or changes in the climate. Concerning the latter, stratospheric dust and sulfates released by the impact (as well as soot from fires) would remain in the atmosphere for a year or so and result in severe global cooling, a similar consequence to massive volcanic eruptions. Carbon may also be injected into the atmosphere, and combined with CO₂ from fires could lead to a longer-term (10⁵ yr) global warming. A very large impact could even eject the atmosphere. It has been estimated that the Chicxulub crater (diameter 180 km) was caused by an impactor with a kinetic energy of 10⁸ Mt TNT equivalent (1 Mt TNT=4.2×10¹⁵ J) (Toon et al. Reference Toon, Turco and Covey1997). This may have released enough dust to make the atmosphere so opaque that photosynthesis stopped and animals had insufficient light to forage for food. Assuming this was an asteroid with a typical density of 2500 kg m⁻³ and impacted with a relative velocity of 15 km s⁻¹, it would have had a diameter of 10–15 km. For comparison: the Tunguska object and the impactor at Meteor Crater in Arizona both had a kinetic energy of 10–15 Mt TNT and diameters of around 50 m; Krakatoa exploded with about 50 Mt TNT (Shoemaker Reference Shoemaker1983; Toon et al. Reference Toon, Turco and Covey1997); the Hiroshima nuclear bomb had a yield of 13 kt TNT (it had an efficiency of just 1%); the most powerful nuclear weapons tested by the Americans and Soviets had yields of up to 50 Mt TNT (Garwin & Charpak (Reference Garwin and Charpak2001); the largest weapons in current arsenals are ‘only’ around 1 Mt TNT). Impacts with an energy of 10⁴–10⁵ Mt TNT (corresponding to comets or asteroids about 1 km in size) become significant on a global scale and are estimated to occur once every 300 000 years or so (Toon et al. Reference Toon, Turco and Covey1997; Chapman Reference Chapman2004).

The main sources of potential impactors are Near-Earth Asteroids (such as the Atens, Apollos and Amors) and comets Shoemaker (Reference Shoemaker1983). Ongoing surveys for near-earth asteroids have either already detected or will soon detect almost all down to sizes of about 1 km (e.g. Morrison (Reference Morrison2003), Harris (Reference Harris2008)). In contrast, the survey completeness down to Tunguska-sized objects (50 m, which are estimated to impact once every thousand years or so) is just a few percent.

The major (and perhaps sole) source of long-period comets entering the inner Solar System is the Oort cloud. This is composed of minor bodies which orbit the Sun with aphelion distances of around 30 000–100 000 AU (cf. 30 AU for the approximate radius of Neptune's orbit) and are believed to be a remnant of the formation of the Solar System. As the size of the Oort cloud is similar in magnitude to the distance to the nearest star (Proxima Centauri, at 270 000 AU), the Oort cloud could be perturbed by the passage of nearby stars or Giant Molecular Clouds (GMCs). This could increase the frequency with which comets are kicked into the inner Solar System and potentially hit the Earth. Such a perturbation would probably release many comets, creating a comet shower in the inner Solar System and many impacts on the Earth spread over a few Myr. This is consistent with claims for some mass extinctions being drawn out over a similar timescale. It could also explain the subsequent discovery that the iridium feature at some extinction boundaries is not a sharp spike but has ‘shoulders’, the finite width of the feature reflecting the stochastic distribution of impacts (Glen (Reference Glen1994), p. 70). Even without an actual impact, it has also been suggested that dust from comets could enter the Earth's atmosphere and affect its climate (Hoyle & Wickramasinghe 1978; Torbett Reference Torbett1989; Shaviv Reference Shaviv2003).

Given that the stellar density increases towards the Galactic midplane and to a lesser extent inside spiral arms, both have been suggested as triggers for the perturbations (Napier & Clude Reference Napier and Clube1979; Napier Reference Napier1988). The perturbation could also kick comets away from the Sun, and as the other stellar system presumably has its own Oort cloud, our Sun may capture its comets (Clube & Napier Reference Clube and Napier1982a). (This has implications for the interpretation of dates of impact material in the Solar System: it may have come from another stellar system.) Clube & Napier (Reference Clube and Napier1982b) speculate that large impacts change the angular momentum of the Earth's core and mantle which could in turn trigger geomagnetic reversals and plate tectonic activity. In other words, different phenomena as recorded in the geological record may have a causal connection, which may in turn be related to mass extinctions or climate change.

Heisler & Tremaine (Reference Heisler and Tremaine1989) have suggested that the Galactic tide is probably a more significant source of gravitational perturbations than the passage of GMCs, in which case we would expect no relation between impacts and disk plane crossing or spiral arm passages. Wickramasinghe & Napier (Reference Wickramasinghe and Napier2008) estimated the flux of comets due to perturbations from the Galactic tide and molecular clouds. They find that the flux increases about an order of magnitude above the background rate on timescales of 25–35 Myr, which is consistent with what they describe as a weak periodicity in the cratering record of 36 Myr.

The perturbation/impact mechanism is certainly a plausible one for causing mass extinctions. The relevance depends on the number and mass of comets in the Oort cloud as well as the size and frequency of the perturbing effects. Large impacts could cause significant devastation and wipe out species, but whether they actually would depends on the complex reaction of ecosystems. The climate would recover relatively quickly (decades, unless significant amounts of CO₂ are released) so the effect on the biosphere is an impulse on geological timescales.

Cosmic rays

It has been suggested that Galactic cosmic rays could have an impact on the Earth's climate via cloud formation (e.g. Shaviv & Veizer (Reference Shaviv and Veizer2003), Carslaw et al. (Reference Carslaw, Harrison and Kirkby2002), Shaviv (Reference Shaviv2005), Kirkby (Reference Kirkby2007)). The basic argument is: (1) cosmic rays cause ionization in the troposphere; (2) these ions act as nucleation sites for water droplets which form clouds; (3) low altitude clouds contribute a net negative radiative forcing (cooling). Hence an increased cosmic ray flux would cause global cooling.

The cosmic ray–cloud/climate mechanism has many uncertainties. While cosmic rays are an important source of ionization in the atmosphere, it is not yet clear whether they are an important source of nucleation compared to neutral molecules (e.g. Jørgensen & Hansen (Reference Jørgensen and Hansen2000)). Even if they are, these nuclei must first grow (by a factor of a million in volume) by condensation and coagulation into cloud condensation nuclei (CCN) before becoming a source of cloud formation. The mechanism by which this growth proceeds remains uncertain. Partly for this reason, it has not yet been established whether the observed amplitudes of variation in the cosmic ray flux induces sufficient variation in the CCN density (Kirkby Reference Kirkby2007). Even if cloud drop nucleation does increase with cosmic ray flux, this does not necessarily translate into a larger areal coverage of clouds: it could rather increase their height or optical depth (Jørgensen & Hansen Reference Jørgensen and Hansen2000; Carslaw et al. Reference Carslaw, Harrison and Kirkby2002). Another important point is that clouds have both a cooling effect (by reflecting sunlight) and a warming one (by reradiating thermal radiation back to the Earth). Although it seems likely that low altitude clouds (<3 km) produce a net negative radiative forcing (e.g. IPCC (Reference Solomon, Qin, Manning, Chen, Marquis, Averyt, Tignor and Miller2007) Section 8.6.3.2, Kirkby (Reference Kirkby2007)), they will only contribute a net cooling when they occur over land or sea which has a lower albedo than the cloud. Over arctic regions, snow and ice provide a strong cooling by reflecting sunlight. Low altitude cloud cover here will reduce this effect and thus contribute a net warming.

Possible cloud formation is not the only impact of cosmic rays. The ions generated could set up a global atmospheric electric current which itself may have other atmospheric effects (Carslaw et al. Reference Carslaw, Harrison and Kirkby2002). Furthermore, the muons created by high energy cosmic rays from supernovae or gamma-ray bursts could kill organisms directly or damage their DNA.

The Earth is exposed to a continuous flux of Galactic cosmic rays (with energies of a few to a few tens of GeV), a large portion of which is believed to originate in shock fronts in supernova remnants (e.g. Lockwood (Reference Lockwood, Haigh, Lockwood and Giampapa2005)). The cosmic ray flux reaching the Earth is modulated by solar activity via the interaction of the solar wind with the Earth's magnetosphere (see the subsection ‘Solar variability’). Thus solar variability on a timescale of years is a potential mechanism for cosmic ray induced climate variability. Cosmic rays are also emitted from the Sun itself, generated at the shock fronts of explosive events on the Sun's surface (e.g. flares, coronal mass ejections) and typically have energies below 1 GeV. In addition, a nearby supernova would generate a large and potentially lethal burst of cosmic rays (subsection ‘Supernovae and gamma-ray bursts’). The progenitors of core-collapse supernovae are short-lived, massive stars in star forming regions. As these (and their remnants) are concentrated towards the Galactic plane and in spiral arms, this has motivated some researchers to look for evidence of a correlation between climate/extinction and the solar motion on timescales of tens to hundreds of million years (see the section ‘Solar motion through the Galaxy’).

What evidence is there for cosmic rays affecting climate on geological timescales? Shaviv & Veizer (Reference Shaviv and Veizer2003) compare a δ¹⁸O temperature proxy constructed by Veizer et al. (Reference Veizer1999) over the Phanerozic with the cosmic ray flux inferred by Shaviv (Reference Shaviv2003) from meteorites and three atmospheric CO₂ proxies. They find that the cosmic ray flux – but not the CO₂ level – correlates with the temperature record and so conclude that CO₂ has much less effect on global temperatures than other research has shown. Rahmstorf et al. (Reference Rahmstorf2004) refute this result on three grounds. First, they claim that the purported correlation is largely a result of several arbitrary ‘adjustments’ to the data. Second, they note the tenuousness of the adopted method of inferring variations in the cosmic ray flux from meteorites (see the subsection ‘Spiral arm crossings’). Third, they report work by Royer et al. (Reference Royer, Berner, Montañez, Tabor and Beering2004) which shows that the temperature calibration of the δ¹⁸O proxy from Veizer et al. (Reference Veizer1999) must be corrected for sea water pH. When this is done, the correlation between δ¹⁸O and cosmic ray flux vanishes. Royer et al. (Reference Royer, Berner, Montañez, Tabor and Beering2004) acknowledge that cosmic rays may have some influence on climate, but that this is probably minor on multimillion year timescales compared to the effect of CO₂.

The cosmic ray–climate link is particularly controversial because some researchers have claimed it explains a significant part of post-industrial global warming. The palaeontological record and a supposed influence on it by astronomical phenomena have been adopted to support this by claiming that if cosmic rays are relevant to climate on hundred million year timescales then they must be relevant on decadal timescales. Although this is a logical non sequitur – we must at least consider the amplitudes of the effects – it is still relevant to ask what evidence there is for a link on other timescales.

A specific claim is that cloud cover over the past few decades correlates well with both the cosmic ray flux and the solar activity (e.g. Friis-Christensen & Lassen (Reference Friis-Christensen and Lassen1991), Svensmark & Friis-Christensen (Reference Svensmark and Friis-Christensen1997), Marsh & Svensmark (Reference Marsh and Svensmark2000)). The idea is that solar activity affects the solar wind, which in turn provides the Earth with some shielding against Galactic cosmic rays. Those reaching the Earth interact further with the Earth's magnetosphere, resulting in a latitude dependence of the cosmic ray flux. Many of these claims have been rebutted (e.g. Jørgensen & Hansen (Reference Jørgensen and Hansen2000), Laut (Reference Laut2003), Damon & Laut (Reference Damon and Laut2004), Sloan & Wolfendale (Reference Sloan and Wolfendale2008)). For example, the data and data analysis supporting the claims of Friis-Christensen & Lassen (Reference Friis-Christensen and Lassen1991) and Svensmark & Friis-Christensen (Reference Svensmark and Friis-Christensen1997) have been strongly criticized by Laut (Reference Laut2003) and Damon & Laut (Reference Damon and Laut2004). After correcting for apparent flaws in the methodology, they show that there is no link between cloud cover and cosmic rays. Jørgensen & Hansen (Reference Jørgensen and Hansen2000) note that the reported correlation is relatively weak and that better established mechanisms of observed events provide a more plausible explanation (e.g. El Niño–Southern Oscillation, volcanism). Damon & Laut (Reference Damon and Laut2004) further report on new data covering 1992–2003 which confirm no correlation between cosmic ray flux and total (all altitudes) global cloud cover. Furthermore, before 1994 changes in the cloud cover lag behind changes in the cosmic ray flux by six months, whereas from what we know of the cosmic ray mechanism the lag should be not much more than one day (Laut Reference Laut2003).

Using two independent estimates of low altitude cloud cover from the International Satellite Cloud Climatology Project over the period 1983–1999, Kristjánsson et al. (Reference Kristjánsson, Staple, Kristiansen and Kaas2002) show that there is either no significant correlation or even a negative correlation between cosmic rays and low altitude cloud cover. On the other hand, Pallé & Butler (Reference Pallé and Butler2002) summarize arguments for and against a cosmic ray–cloud link over the past 50–120 years and conclude that a lack of quality data does not allow us to ‘totally dismiss the link between [Galactic cosmic rays] and cloudiness’.

Erlykin et al. (Reference Erlykin, Sloan and Wolfendale2009) show that there is a common variation between the cosmic ray flux in the Earth's atmosphere, sunspot number, solar irradiance and global average surface temperature over the past 50 years (1956–2002) which varies on a timescale of twice the 11-year solar cycle. (Only two ‘cycles’ are seen so we cannot call this ‘periodic’. The analysis uses an 11-year moving average centered on the point of interest. While this smooths shorter timescale variations, it does not eradicate them and, importantly, by using a symmetric smoothing window no phase shift is introduced.) The temperature, solar irradiance and sunspot number variations are in phase, whereas the cosmic ray flux lags behind by 2–4 years, so it cannot be a cause of the temperature variation. They further estimate that the direct impact of cosmic rays on the radiative forcing of the Earth is less than +0.07° C since 1956, so contributing less than 14% of global warming.

In summary, there is no strong evidence that cosmic rays have significantly influenced climate on geological timescales. There is now a broad consensus among scientists that, although cosmic rays may have an effect on cloud microphysics (via a mechanism only poorly understood) and thus on climate, the data show that they had at most a minor influence on post-industrial global warming (IPCC Reference Solomon, Qin, Manning, Chen, Marquis, Averyt, Tignor and Miller2007). This is further reinforced by the fact that we have much more support for an alternative hypothesis for climate change, namely CO₂ and other greenhouse gases. This consensus is not always reflected proportionally in the media.

Supernovae and gamma-ray bursts

Supernovae release large amounts of energy in the form of hard x-rays and cosmic rays. These could cause widespread extinction through at least three mechanisms. First, the radiation can kill organisms (on the hemisphere facing the blast) via direct cell destruction or damage to DNA. Second, ionizing radiation creates nitric oxide (NO) in the upper atmosphere, destroying ozone. A single blast could leave the ozone layer depleted for hundreds of years, exposing life to harmful solar UV radiation which can damage DNA even in water at depths of several metres (Ruderman Reference Ruderman1974). Third, NO₂ formed from the NO is a strong absorber of visible radiation from the Sun leading to global cooling (Thomas et al. Reference Thomas2005). As just discussed, cosmic rays may also affect climate via cloud formation. Tanaka (Reference Tanaka2006) estimates that a supernova within 12–15 pc would increase the flux of 10–100 GeV cosmic rays by a factor of 4–8, although how this translates to condensation nuclei is unclear. Apart from the initial blast, cosmic rays are also emitted from the supernova remnant for millions of years thereafter by shock wave processes.

Ellis & Schramm (Reference Ellis and Shramm1995) examined the possibility that nearby supernovae could have caused mass extinctions on Earth. Using a supernova rate of 0.1 per year in our Galaxy and an average stellar density of 1 pc⁻³, they estimate that a supernova would occur within 10 pc of the Sun every 240 Myr or so (this is only an order of magnitude estimate). Assuming that core-collapse of massive stars are the dominant cause of Galactic supernovae, then the probability of a nearby supernova increases as the Solar System crosses spiral arms, although not by a lot because the space density of supernovae may not be significantly larger in spiral arms (see the subsection ‘Spiral arm crossings’).There will also be a change in the local space density of massive stars (and thus supernovae) as the Sun oscillates vertically through the Galactic plane. The magnitude of variation depends on the scale height for massive stars and the amplitude of the motion. It we adopt 100 pc for the former with an exponential profile (Da-li & Zi Reference Da-li and Zi2008) and 70 pc for the latter (Gies & Helsel Reference Gies and Helsel2005), then the maximum increase in density is only exp(70/100)=2.

Gamma-ray bursts would have a similar effect on the Earth as supernovae, but could be effective out to distances of several kpc. Thomas et al. (Reference Thomas2005) have modelled their impact on the atmosphere and biosphere is some detail. Melott et al. (Reference Melott, Lieberman, Laird, Martin, Medvedev, Thomas, Cannizzo, Gehrels and Jackman2004) suggest that gamma-ray bursts occur at a rate which could cause two or more mass extinctions on the Earth every billion years, and single out in particular the late Ordovician event.

Solar variability

The total solar radiation reaching the top of the Earth's atmosphere (the solar irradiance) is 1367 W m⁻² averaged over the orbit. Convection and magnetic activity in the Sun's atmosphere (Sun spots) result in a variation of the solar irradiance by 0.1% over the 11 year solar cycle. This amplitude is too small to cause the ice ages on the Myr timescale discussed earlier. Lockwood (Reference Lockwood, Haigh, Lockwood and Giampapa2005) states that periodic variations in the Earth's climate on timescales of a decade or less are mostly smoothed out by the atmosphere–ocean coupling and large thermal capacity of the oceans. (This would therefore include the annual 7% variation in solar irradiance due to the eccentricity of the Earth's orbit.) However, UV and shorter wavelength variations are larger (Shaviv Reference Shaviv2005). Although UV is mostly absorbed in the stratosphere, the variations could be propagated down through the atmosphere. Furthermore, the intensity of spots varies on longer cycles, and longer durations with fewer spots (such as the Maunder Minimum around 1645–1715) have been associated with climate change. According to Lockwood & Fröhlich (Reference Lockwood and Fröhlich2007), the evidence suggests that solar variability has had an impact on climate over the past few centuries.

In addition to variations in the electromagnetic flux, there are variations in the cosmic rays (solar protons and electrons) emitted from the Sun. These charged particles form the solar wind and are responsible for the heliosphere, which modulates the flux of Galactic cosmic rays reaching the Earth. An increase in the solar activity ‘strengthens’ the heliosphere and so lowers the cosmic ray flux reaching the Earth. The cloud mechanisms described in the subsection ‘Cosmic rays’ have been invoked by some authors to provide the connection between solar activity and Earth climate. Lockwood & Fröhlich (Reference Lockwood and Fröhlich2007) show that data on the potentially relevant phenomena of the Sun – total solar irradiance, solar magnetic flux, cosmic ray flux from neutron counts, sun spot number – since 1985 show variations in the direction opposite to that required for them to be responsible for the recent rise in global temperatures. According to G. Feulner (private communication, April 2009), terrestrial temperature variations over the 11-year solar cycle are detectable, and, moreover, that they can be reproduced within climate models by the variation in the solar irradiance alone, without having to invoke effects of cosmic rays on the cloud cover.

Variations in the Earth's orbit about the Sun

The orbit of the Earth about the Sun is nominally an ellipse, but it is perturbed by the gravitational force of other bodies in the Solar System. This causes both the inclination of the Earth's orbital plane and the eccentricity of its orbit to vary. The Sun and Moon also impart torques on the Earth which cause the Earth's spin axis to precess and its obliquity (angle with respect to the orbital axis) to vary.

These perturbations have been modelled very accurately using classical mechanics. By expressing them as a series expansion we can identify terms with different periods and amplitudes. The eccentricity of the Earth's orbit varies between almost 0 and 0.05 (it is currently 0.017 and decreasing) with dominant periods of 95, 125 and 400 kyr with relative amplitudes of 1.2, 1.0 and 1.7 respectively. Together these account for 90% of the signal. (These and other figures in this paragraph are taken from Muller & MacDonald (Reference Muller and MacDonald2000).) The average distance of the Earth from the Sun depends on the eccentricity, so variations in this translate into variations in the average annual solar irradiance.Footnote ¹⁰ In contrast, variations in precession, obliquity and orbital inclination only affect the geographical distribution of the flux. These can nonetheless induce an ice age by preventing winter ice from melting in the summer, and subsequent changes in tree cover and sea ice enhance this via a positive feedback (IPCC (Reference Solomon, Qin, Manning, Chen, Marquis, Averyt, Tignor and Miller2007), Chapter 6). The precession period of the Earth's spin axis is 25.8 kyr, but as this is modulated by eccentricity variations, the precession parameter (which is relevant to variation of the solar irradiance) shows periods at 19, 22 and 24 kyr. The axis of the Earth's orbital plane also precesses. The orbital inclination varies between 0° and 4° at a dominant period of about 70 kyr. Because the directions of both the orbital axis and the spin axis vary with different periods and amplitudes, then the angle between them (the obliquity) varies too. The result is that the extent of the tropics varies between 22.1 and 24.5° (it is currently 23.5° and decreasing) with a dominant period of 41 kyr (and weaker periods at 29 and 53 kyr).

The similarity between periods in the climate record over the past 3 Myr (discussed in the subsection ‘Periods in climate data over the past few million years’) and periods in the perturbation terms of the Earth's orbit over a similar timescale, has lead many to suggest a causal connection (e.g. Haye et al. (Reference Haye, Imbrie and Shackleton1976)). This is sometimes called the Milankovitch theory of the ice ages, although others proposed it before him and others have modified it since. Its broad formulation is generally accepted to explain the occurrence of recent ice ages (e.g. IPCC (Reference Solomon, Qin, Manning, Chen, Marquis, Averyt, Tignor and Miller2007)), although some of the details are debated. Part of the evidence comes in the form of very narrow peaks in the climate power spectrum. Narrow peaks imply a non-dissipative process (Muller & MacDonald Reference Muller and MacDonald2000). This in turn is taken to imply a mechanism which is similarly non-dissipative, such as planetary orbits. One may argue that the simultaneity of ice ages in the northern and southern hemispheres argues against this, but because the northern hemisphere possesses two thirds of the terrestrial land surface, it is the northern hemisphere insolation which is responsible for triggering the ice ages.

In summary, while Earth axis and orbit variations may well explain some climate change on timescales of a few tens to a few hundreds of kyr, they seem not to offer an explanation for change or periodicity in climate on longer timescales.

Other mechanisms

Encounters of the Solar System with interstellar clouds could produce global cooling as dust from the cloud lowers the solar irradiation. Shaviv (Reference Shaviv2003) has coupled this with the idea that a bow shock from the cloud reduces the size of the heliosphere to below 1 AU and thus exposes the Earth to a greater flux of Galactic cosmic rays. Contrarily, Hoyle & Lyttleton (Reference Hoyle and Lyttleton1939) suggested that matter from an interstellar cloud falling into the Sun could raise the Sun's luminosity (via the release of gravitational energy) and that this in turn could trigger an ice age via increased precipitation. McCrea (Reference McCrea1975) took up this idea and suggested that sufficient matter cloud be provided by the Sun's passage through dust lanes at the edge of spiral arms.

Another suggestion is that the Sun is in a wide binary system, with a faint M dwarf or later-type companion in a long period orbit (Davis et al. Reference Davis, Hut and Muller1984; Whitmire & Jackson Reference Whitmire and Jackson1984). If the companion (colourfully named ‘Nemesis’ in the article of Davis et al. (Reference Davis, Hut and Muller1984)) is in a highly elliptical orbit, then perihelion passages at around 30 000 AU would be close enough to perturb the Oort cloud and eject large numbers of comets toward the inner Solar System (see subsection ‘Minor body impacts, Oort cloud perturbation and comet capture’). This idea was originally proposed as a mechanism to explain the claimed 26 Myr periodicity in extinction from Raup & Sepkoski (Reference Raup and Sepkoski1984). To achieve an orbital period of this order, then with the required eccentricity of around 0.7 (or more) the companion would need a semi-major axis of order 10⁵ AU. Davis et al. (Reference Davis, Hut and Muller1984) estimated that such a companion could perturb the orbits of some 10⁹ Oort comets, of which 25 would hit the Earth. However, the semi-major axis for the proposed companion is unusually large for binary systems (0.5 pc), leading some authors to suggest that the Galactic tide or close encounters with stars or interstellar clouds would unbind the system on a timescale of a Gyr (Torbett & Smoluchowski (Reference Torbett and Smoluchowski1984)). Even if it remains bound, these perturbations are likely to make the orbit unstable: one calculation predicts that a 26 Myr orbital period would vary by 10–20% over 250 Myr (Hut Reference Hut1984). Interestingly, this could actually speak in favour of this mechanism if extinctions are demonstrated to be quasi-periodic (rather than strictly periodic or non-periodic). No candidate for Nemesis has been found in deep, all-sky surveys.

Motion perpendicular to the Galactic plane

There have been many attempts to model the motion of the Sun perpendicular to the plane. Using a model for the gravitational potential of the disk, Bahcall & Bahcall (Reference Bahcall and Bahcall1985) derive a quasi-harmonic motion with a period ranging from 52 to 74 Myr, the range reflecting uncertainties in both the distribution of dark matter (which dominates) as well as the current vertical position and velocity of the Sun. The maximum displacement from the plane varies from 49 to 93 pc, with the motion being significantly non-harmonic beyond about 40 pc (because not all the matter is concentrated in the plane, but rather falls off exponentially).

Svensmark (Reference Svensmark2006) takes a fundamentally different approach. Instead of fitting the motion of the Sun using astronomical data, he fits it to δ¹⁸O temperature proxy data (subsection ‘Geological and biological data’) from the past 200 Myr, using the ad hoc assumption that the temperature of the Earth, T, varies as the square of the distance, z, of the Sun from the Galactic midplane. Thus δ¹⁸O becomes a proxy for z ². Adopting the harmonic oscillator model for the z motion, T then varies with twice the oscillation frequency. As his δ¹⁸O sample shows a period of around 30 Myr, this translates to a vertical oscillation period of 60 Myr. (The disk potential model allows for two perturbations due to two spiral arm crossings, so the vertical period is not strictly periodic. See the subsection ‘Spiral arm crossings’.)

There is a methodological problem with this approach, however, because it assumes that the Earth gets hotter the further it is from the Galactic plane. This, in turn, is based on the assumption that the mean global temperature is controlled by cosmic rays via cloud cover, as discussed in the subsection ‘Cosmic rays’. However, because the solar motion is derived from the δ¹⁸O (temperature) data, we cannot then use this to claim that the solar motion fit provides support for the cosmic ray model. We could only do this if we test the assumption T∝z ² by comparing the δ¹⁸O data with an independent determination of the solar motion. Otherwise we have circular reasoning. Note also that, because the assumption is not based on a physical model – it is chosen as the simplest form which gives the required symmetry of T(z) – z amplitudes cannot be derived from δ¹⁸O. It should also be mentioned that the phase of the fit (i.e. the current z coordinate of the Sun, taken as 9±4 pc north of the plane) is an input to the model, not a prediction. The most we can conclude is that the selected δ¹⁸O data appear to show a period which is of order half that of other determinations of the solar vertical period. However, this study does not really lend support to the idea that an increased flux of cosmic rays from the Galactic plane causes a net cooling of the Earth.

Svensmark is by no means the only author to use geological data to try and infer solar motion. There were several efforts in the 1980s, such as the work of Shuter & Klatt (Reference Shuter and Klatt1986), to try to constrain the dark matter content of the disk in this way (see also Bahcall & Bahcall (Reference Bahcall and Bahcall1985) and Wickramasinghe & Napier (Reference Wickramasinghe and Napier2008)). Given the significant uncertainties in the geological record, plus the considerable doubt about the existence of (stable) periods, this is a rather suspect approach.

Medvedev & Melott (Reference Medvedev and Melott2007) also constructed an astronomical model in which extinctions are caused by cosmic rays. They are motivated by the apparent 62±3 Myr period in biodiversity variation from Rohde & Muller (Reference Rohde and Muller2005), discussed in the subsection ‘Periodicities in biological variation over the Phanerozoic’. This is close to the period of the solar vertical oscillations in many dynamical models and thus twice the period of Galactic plane crossings. If there is a causal connection, then extinctions cannot be associated with star forming regions, spiral arms or anything else concentrated in the plane, but would have to be associated with some plane asymmetry. They note that the maxima in Rohde & Muller's analysis coincide with the Sun being near to its northern-most displacement, which is the direction to the Virgo supercluster. Based on these considerations, they propose a model in which a Galactic bow shock is produced by the Galactic wind and the motion of the Galaxy towards Virgo. This shock is a source of cosmic rays. The Galactic magnetic field shields the Sun from these to an extent which varies with distance from the Galactic midplane: there is considerably less shielding (perhaps five times the cosmic ray flux on the Earth) when the Sun is at its northern-most displacement.

They test this using a model of the solar motion from Gies & Helsel (Reference Gies and Helsel2005) to predict the cosmic ray flux at the Earth over the Phanerozoic and compare it with Rohde & Muller (Reference Rohde and Muller2005) results. Minima in diversity phase well with predicted cosmic ray maxima, perhaps implying that cosmic ray-induced cooling causes mass extinctions (although they do not pinpoint a specific extinction mechanism). Although the data interpretation involves many assumptions and the model is rather sketchy, the authors argue that the 62 Myr period in the biodiversity data demands an extragalactic explanation.

Many other papers have invoked plane crossings to explain mass extinctions or climate variations with periods around 25 to 33 Myr (e.g. Raup & Sepkoski (Reference Raup and Sepkoski1984), Raup & Sepkoski (Reference Raup and Sepkoski1986), Rampino & Caldeira (Reference Rampino and Caldeira1992), Napier (Reference Napier1988)) as discussed in the subsections ‘The 26 Myr period’ and ‘Periodicities in the geological record over the Phanerozoic’. Two frequently invoked triggers are an increased supernova rate (subsection ‘Supernovae and gamma-ray bursts’) and increased rate of comet impacts (section ‘Minor body impacts, Oort cloud perturbation and comet capture’). However, the increases may not be that large. For example, the local supernova rate at midplane is perhaps only double that at maximum distance from the plane (subsection ‘Supernovae and gamma-ray bursts’), so the exposure to cosmic rays (for example) from supernovae and their remnants would also only change by this amount. Whether this is a physically important difference depends on the details of the climate change or extinction model. Some authors have proposed a connection between the impact cratering record, mass extinctions and Galactic plane crossings (subsection ‘Impact cratering’), but one claimed period of around 13 Myr for cratering does not fit here.

While the estimates of the parameters of the solar orbit have a wide range, there is a reasonable consensus that the Sun is currently near the midplane (z≃15 pc compared to a vertical oscillation amplitude of around 70 pc). The Sun therefore recently underwent a midplane passage, yet the lack of evidence for a major extinction within the past 15 Myr has led some to suggest that plane crossings cannot be the cause of mass extinctions (e.g. Leitch & Vasisht (Reference Leitch and Vasisht1998)). We could argue that not every midplane crossing will cause a mass extinction (or ice age), but when we permit ourselves to pick and choose we cannot look for periodic phenomena, and indeed could fit almost any model to the geological data.

In the discussions of sections ‘Evidence from geological time series’ and ‘Issues in time series analysis’, I concluded that there is no good evidence for periodicities in climate or mass extinction in the range 25–33 Myr. This seems to rule out a relevant periodic influence of Galactic plane crossings, and other studies looking for a specific connection have not provided convincing evidence of a causal link. The recent 62 Myr period in biodiversity appears to be on a stronger footing, yet its identification is also subject to the limitation of ‘incomplete inference’ discussed in the subsection ‘Hypothesis testing’, and is not without criticism (see the subsection ‘How reliable are the data, how appropriate are the methods?’). It remains an interesting suggestion to be explored further in other data sets. However, given that the periodicity is only significant over the interval 520–150 Myr BP, a connection to the solar z-motion is not obvious. The link to the uncertain cosmic ray–cloud connection and even more uncertain bow shock model are additional, unresolved steps in the causal chain.

Spiral arm crossings

Spiral arms are sites of increased star formation in the Galactic disk, characterized in external galaxies by the excess of massive young blue stars and giant molecular clouds. The origin of spiral arms is debated. The classical model originating in the work of Lindblad (Reference Lindblad1938) and Lin & Shu (Reference Lin and Shu1964) explains them as gravitational in origin, whereby they can be understood as density waves propagating through the disk. In this model the spiral pattern is fixed and rotates as a rigid structure with angular velocity Ω_p (the pattern speed). In contrast, stars and gas in the disk rotate differentially, so they will pass in and out of the arms even if on exact circular orbits. Within the corotation radius the stars overtake the spiral arms and beyond it are overtaken by the arms. Only at the corotation radius do stars on circular orbits not move relative to the spiral pattern. The position of this radius is disputed, and various studies have placed it at Galactocentric radii ranging from 3 to 16 kpc, some also placing it near to the Sun (e.g. Dias & Lépine (Reference Dias and Lépine2005)). Different investigations have likewise inferred a wide range of values of the pattern speed, with estimates ranging from 10 to 30 km s⁻¹ kpc⁻¹ (e.g. Gies & Helsel (Reference Gies and Helsel2005)). There is not yet a consensus from observations that our Galaxy really does have a global spiral pattern (grand-design spiral arms) rather than just many short arms segments a few kpc long. Even with the former there is debate whether there are four or two arms. A recent compilation and synthesis of the literature is given by Vallée (Reference Vallée2008).

Several studies have implicated passages of the Sun through spiral arms in mass extinctions and/or climate change. Possible mechanisms are similar to those posited in the case of disk plane crossings, namely encounters with supernovae, gravitational perturbations of the Oort cloud and comet capture. These studies all assume that there is a grand-design spiral structure, although they differ in their assumptions about the shape, size and speed of the spiral pattern and number of arms. Some turn the problem around and assume that extinctions or ice ages were caused by arm crossings and use this to try and constrain the spiral structure.

Gies & Helsel (Reference Gies and Helsel2005) use the model for the disk gravitational potential from Dehnen & Binney (Reference Dehnen and Binney1998a) and for the solar motion derived from Hipparcos data by Dehnen & Binney (Reference Dehnen and Binney1998b) to derive the motion of the Sun over the past 500 Myr. Adopting a four-arm model of the spiral structure (with an additional local arm segment) they determine when passages occurred for various values for the pattern speed. They compare these with the midpoints of four ice age epochs. (They offer several different times for the occurrence of the ice ages from different sources, with discrepancies of up to tens of Myr.) Using a value for the spiral arm pattern speed of Ω_p=20 km s⁻¹ kpc⁻¹ (in the middle of the range of published estimates), they show that there is little correlation between arm crossings and the ice age midpoints. Testing a range of pattern speeds, as well as values for R ₀ and the disk scale length in the gravitational potential model of the disk, they obtain a better coincidence when adopting a pattern speed of Ω_p=14.4 km s⁻¹ kpc⁻¹ (Fig. 7). (They keep the solar angular motion fixed at in their simulations; it is the relative angular speed of the Sun to the spiral pattern, , which is relevant.) This gives arm crossing midpoints at 80, 156, 310 and 446 Myr BP. However, there are three parameters which can be varied in their model (plus four alternative ice age dating schemes are considered) whereas the quality of the fit is judged based on just four ice age epochs. In other words, the model has considerable degrees of freedom to be constrained by little data, so it should not be hard to obtain a good fit even if the model parameters are constrained in their range. The ice ages are long and the spiral arms wide (they adopt 0.75 kpc), so some coincidence is almost inevitable.

Fig. 7.

Path of the Sun relative to the spiral pattern in the Galaxy as reconstructed by Gies & Helsel (Reference Gies and Helsel2005) using=14.4 km s⁻¹ kpc⁻¹. The Galactic centre is at (0,0) and the Sun currently at (0,8.5) (the stars and arms rotate clockwise). The diamonds mark time intervals of 100 Myr, crosses selected mass extinctions and the thick solid line icehouse intervals. Reproduced from Gies & Helsel (Reference Gies and Helsel2005) with kind permission of D.R. Gies and the AAS.

Gillman & Erenler (Reference Gillman and Erenler2007) examine the temporal distribution of numerous geological markers (including extinctions and impact craters) over the past 700 Myr. By wrapping these at a period of 180 Myr, they notice that the events fall into three reasonably well-separated ‘zones’. If this period is associated with the period between arm passages, then the zones correspond to different phases of the orbit between passages. A fit with a time-domain model (a generalized linear model) gives a period of 175.96±0.43 Myr. (This surprisingly small uncertainty is the formal fitting error as reported: it does not explicitly account for the uncertainties in the geological dating.) Adopting a four-arm model of the Galaxy with rigid rotation, they identify this figure as the period between arm crossings, which implies Ω_p=18.2 km s⁻¹ kpc⁻¹. However, the geological events are not concentrated around arm crossings. As their Figure 2 shows, the events are spread across the whole phase of the purported crossing period, so the relevance of the arms is not obvious. The importance of the period of 176 Myr seems to be that it divides the geological events into the three zones, but neither the statistical nor astronomical significance of these is clear. There is no significance analysis of this or any other period, so its uniqueness is hard to assess. A connection with the spiral structure of this model would also require the spiral pattern to have rotated rigidly and not evolved for the past 700 Myr (four arm crossings).

Shaviv (Reference Shaviv2003) studied the correlation between cosmic ray fluxes, ice age epochs and spiral structure. He suggests that there is a causal connection and, depending on the data used, derives a periodicity for spiral arm crossings of between 134±22 and 163±50 Myr, with a best fit of the climate proxy data to a spiral arm model yielding 143±5 Myr (using a pattern speed of Ω_p=16 km s⁻¹ kpc⁻¹). The fit indicates a lag of the midpoint of glaciations behind arms crossings of 33±20 Myr. He further concludes that the data indicate the spiral pattern to have been stable over the past billion years, although this assumes the cosmic ray/spiral arm mechanism to be responsible for the climate record (and this seems incompatible with the Myr time lag). This work is based on interpreting the observed clustering in meteorite cosmic ray exposure ages (as measured by the ⁴¹K/⁴⁰K isotope ratio) as evidence for variability in the cosmic ray flux. This is a rather indirect method with a chain of assumptions, and it has been argued that these data are consistent with no clustering and thus no significant cosmic ray variability (Jahnke Reference Jahnke2005).

The model for the solar motion of Svensmark (Reference Svensmark2006) discussed in the subsection ‘Motion perpendicular to the Galactic plane’ is one dimensional, but it uses a potential for the disk which includes two perturbations, considered to be spiral arm crossings. The increased mass of the spiral arms leads to an acceleration of the vertical motion, with the epoch, duration and amplitude of the perturbation being free parameters in his fit to the Earth temperature proxy data. (Only two perturbations are included based on the prior evidence of the times of spiral arm crossings.) The inferred times of the arm crossings are 31 and 142 Myr BP. However, as we already discussed in the subsection ‘Motion perpendicular to the Galactic plane’, because the terrestrial temperature is assumed to depend on distance from the Galactic plane, this does not provide independent evidence that the solar motion (let alone cosmic rays) triggers climate change.

An arm crossing per se does not automatically imply the Solar System will experience anything fundamentally different. Torbett (Reference Torbett1989) notes that the stellar density inside spiral arms is a factor of only about 1.1 higher than outside the arms. Svensmark (Reference Svensmark2006) quotes values from the literature of 1.5–1.8 and 1.5–3 for external galaxies. Scoville & Sanders (Reference Scoville, Sanders and Smoluchowski1986) estimate the probability of having a ‘close’ encounter with a GMC inside a spiral arm as only about 0.1 per arm crossing. If crossings occurred on average every 100 Myr, there is a probability of 0.9⁵=0.59 of no close GMC encounter at all in the past 500 Myr. Of course, this depends on what one regards as an ‘encounter’, and as already discussed, the low likelihood, P(D|H), of a hypothesis (without assessing the alternatives) is insufficient to reject it (see the subsection ‘Issues in time series analysis’). Leitch & Vasisht (Reference Leitch and Vasisht1998) estimate the number of supernovae encountered during a spiral arm passage at 0.5, assuming that the supernova only has an effect on the Earth if it passes within 10 pc (Ellis & Schramm Reference Ellis and Shramm1995). This is just an order-of-magnitude figure (not a probability) from a simple volume calculation using the supernova rate (1/30 yr⁻¹ in the whole Galaxy), progenitor lifetime (10 Myr) and scale height (100 pc), plus the spiral arm length, width and pattern speed. The true average rate could easily differ from this by an order-of-magnitude or more. Moreover, GMCs, star formation and supernovae also occur in the disk outside of spiral arms, so even taking the ‘star formation region encounter hypothesis’ to be true, we may not expect a very high correlation of geological events with spiral arm transits.

All of these attempts to associate climate change or extinctions with spiral arm crossings are very sensitive to the exact morphology and pattern speed of the spirals arm. Yet the spiral structure of our Galaxy is poorly known: there is still a debate over whether it has a four-arm or two-arm structure, and estimates of the pattern speed vary by a factor of three (12–30 km s⁻¹ kpc⁻¹ in Table 3 of Shaviv (Reference Shaviv2003); see also Perryman (Reference Perryman2009), Section 9.7). Coupled to this is the uncertainty in the corotation radius of the pattern: the nearer it lies to the Sun, the lower the relative velocity of the Sun with respect to the arms and so the lower the frequency of arm crossings. Moreover, the studies described above assume that the spiral pattern has rotated with a fixed angular speed and with a fixed pattern over hundreds of millions of years. Yet some N-body simulations predict that spiral arms are unstable, showing significant changes in their structure in less than a rotation period (e.g. Sellwood & Carlberg (Reference Sellwood and Carlberg1984)).

Another implicit assumption is that the Sun's motion can be described by a smooth gravitational potential and has not experienced any close encounter for hundreds of Myr. There are thousands of known GMCs in the disk which could give the Sun an additional acceleration. The clouds would have now dispersed so it would be almost impossible to reconstruct these events even with accurate stellar kinematics. Whatever the source, Wielen (Reference Wielen1977) found that the dispersion velocity of stars increases with time due to orbital diffusion. This produces a velocity change of order 10 kms⁻¹ in a single orbit of the Sun around the Galaxy.

In summary, it is likely that the Sun has crossed the spiral arms up to a handful of times over the Phanerozoic. However, the specific conclusions of the cited studies of a connection to climate change or mass extinctions are very sensitive to the very uncertain structure, kinematics and evolution of the spiral arms, and to the uncertainties in the Galactic potential used to reconstruct the motion of the Sun. It is premature to draw a connection, let alone use any apparent one to provide support for a specific extinction/climate change mechanism such as cosmic rays or Oort cloud perturbation. However, we can draw at least one useful conclusion from these studies: as the arms are unlikely to show a perfectly symmetric structure and fixed rotation relative to the Sun, it is unlikely that arm crossings are periodic, and almost certainly not frequent enough to explain the 26 or 62 Myr periods (or values near these) discussed in the subsection ‘Periodicities in biological variation over the Phanerozoic’. Moreover, using these studies we can identify what better astronomical data we need to make progress. Some studies attempted to use the geological record to formally fit (or loosely constrain) models for the spiral structure. In my opinion, given the uncertainties in the data and the large degree of freedom in the models (or the use of numerous assumptions), this is, at best, inconclusive. Although it is frustrating that the astronomical data are not yet adequate, it is nonetheless essential that independent astronomical data be used to model astronomical phenomena, which only then are compared with the geological data.

Radial variations

Goncharov & Orlov (Reference Goncharov and Orlov2003) claim that 13 mass extinction events have a ‘repetition interval’ of 183±3 Myr. It is not clear what they mean by this, because the events are not periodic with this period (their extinction dates range from 2 to 469 Myr BP, much less than 12×183 Myr) and the events (their Figure 1) do not cluster around this period. They suggest that extinctions could be described by the radial motion of the Sun in its orbit around the Galactic centre: they quote a model in which R ₀ varies from 7.96 to 8.15 kpc with a period of 183 Myr. However, a superposition of the variation of R ₀ over time in this model with these extinction events (their Figure 3) shows only a very slight clustering of events at the apocentre and pericentre positions. They note that if only a subset of the events are retained then the clustering is better. However, if one arbitrarily select events then I suspect it would be possible to make a number of periods ‘fit’ the data. Without a systematic analysis of the significance of the clustering at a range of periods and the sensitivity to the data retained, it is not possible to draw any firm conclusions from this work. Of course, it is a priori possible that a subset of extinction events do show periods and are caused by a periodic mechanism. But unless one can find a signature in the data to decide which events to select independently of the period of the mechanism, the data themselves cannot be taken to support this hypothesis. The hypothesis ‘some of the events are periodic’ is simply too general and too flexible to have much predictive power and thus to get much support from the data.