Is the solar system's galactic motion imprinted in the Phanerozoic climate?

A new δ(18)O Phanerozoic database, based on 24,000 low-Mg calcitic fossil shells, yields a prominent 32 Ma oscillation with a secondary 175 Ma frequency modulation. The periodicities and phases of these oscillations are consistent with parameters postulated for the vertical motion of the solar system across the galactic plane, modulated by the radial epicyclic motion. We propose therefore that the galactic motion left an imprint on the terrestrial climate record. Based on its vertical motion, the effective average galactic density encountered by the solar system is 0.172±0.006stat±0.006sysM∘pc⁻³. This suggests the presence of a disk dark matter component.

The Milky Way is a barred spiral galaxy. As a consequence, the solar system revolves around the galaxy and carries out a small oscillation in the direction perpendicular to the dense galactic disk, with a period much shorter than its orbital period1. For very small amplitudes of oscillation, the disk-crossing period P1/2 depends solely on the density ρ0 in the galactic plane, via (with G being the gravitational constant, assuming a constant density disk, otherwise the period may be somewhat larger2). Previous determinations of the half-period range between 30 and 42 Ma234. The methodologies for measuring the density of the galactic disk can be classified into three categories. In the first approach, the known components are simply summed up, giving the total baryonic mass density at the galactic plane of about 0.09 to 45. The second approach is based on the vertical kinematics of stars6. Assuming that a given tracer population of stars is kinematically relaxed, there is then a relation between the vertical density of the tracer stars and their dispersion velocity. This enables derivation of the underlying vertical dependence of the gravitational potential, and from it the density, yielding values from 0.07 to 23. These two approaches are often used to estimate the amount of dark matter (DM) in the disk. It was argued that the stellar kinematics are inconsistent with a higher end estimate of the central mass density7, implying that DM in the disk may account for the deficient mass. Since the Hipparcos data set is presently the most extensive astrometric data available, analyses based on it should be considered the most reliable. They appear to converge towards lower values of 0.10 to 23, suggesting that the unobserved DM in the disk is at most the amount expected from a spherical or oblate DM halo, around 0.01 to . The DM density is important because its exact value is essential for searches of weakly interacting massive particles8.

The third approach to derive the density of the galactic disk is based on the periodicity of the solar systems vertical oscillation (VO). The periodic perturbation of the Oort cloud should increase the population of comets crossing the Earth pathway, potentially leaving an imprint in the terrestrial cratering record. Estimates for cratering periodicity range from 26 to 36 Ma (e.g., ref. 9), but the sparse statistics are disputed10.

Paleo-climatic records could serve as alternative chronometers if a physical mechanism exists to link the solar system's vertical motion to the terrestrial climate. Several suggested mechanisms could potentially do so: 1) perturbation of comets in the Oort cloud that, after disintegration into dust, could potentially cool the climate11, 2) collision with interstellar gas clouds12, 3) climate modulation via cosmic rays131415. In the latter case, the VO would translate into an oscillatory variation in the cosmic ray flux, because the cosmic rays density depends on the distance from the galactic plane16. In addition, the lower interstellar pressure would allow the heliosphere to puff up and increase the cosmic ray energy loses as they propagate into the solar system17.

Since the cosmic ray mechanism can potentially explain climate variations on longer time scales18, we will assume that the link does indeed operate and with it construct in what follows a model for the average climate. However, two points should be noted. First, any VO/climate correlation by itself does not prove that the climate link is through any particular mechanism. Second, the cosmic ray climate link is controversial and it should therefore not be taken for granted. For example, various criticisms were raised on the validity and implications of the long time scale correlations19, the cloud cover cosmic-ray correlation on the decade time scale20 and during Forbush decreases21, or whether the observed link between atmospheric ionization and the nucleation of small particles can affect the number density of cloud condensation nuclei22.

Here we use a new database of δ18O values over 488 Ma, which covers almost the entire Phanerozoic, and that transcend the aforementioned limitations. We will show that it exhibits a very clear 32 Ma periodicity.

The data

The original database of some 4500 δ18O measurements from low-Mg calcitic fossils23, a reflection of Phanerozoic climate24, was shown to encompass periodicities of about 145 and 32 Ma25. This database was later updated26 and once more prior to the present study, while keeping its original structure. It now includes over 24000 δ18O values. It lists separately data for the deep sea (>300 m depth) and surficial waters, the latter separately for the low-, mid- and high-latitudes and, where possible, also for the subtropical realms. Due to biological evolution, and the resulting diversity of fossil phyla, none of the subsets covers the entire 542 Ma of the Phanerozoic time span. This necessitates splicing of the partial records. The entire δ18O and δ13C database and its reduction procedures are described in the Methods below, and elaborated in the Supplementary Materials.

The secular δ18O trend homogenized into 1 Ma bins is reproduced in Fig. 1. Visual inspection of this figure, spectral (Supplementary Fig. S3) and wavelet analyses (Supplementary Figs. S4, S5) show that all subsets (tropical to high-latitude realms) exhibit a similar 30 to 40 Ma oscillation, confirming its robustness and global validity. This permits the splicing procedure over the almost entire Phanerozoic, providing corrections for offsets relative to the low-latitude subset, the latter covering most of the studied time span, are taken into account. We tested 5 possible splicing combinations and they all yielded similar outcomes (see Supplementary Fig. S5), with only small quantitative differences. All subsequent discussion is therefore based on a “master” set that splices the two most complete subsets, the mid-latitude one for the <200 Ma interval and the low-latitude one for the time interval >200 Ma. The detrended and high pass filtered δ18O data of this “master” (ML200) set are reproduced in fig. 2.

Figure 1

The Gaussian filtered δ18O data, separated into four groups.

The green line is the low-latitude, blue the mid-latitude, red the high-latitude, and the black line the deep sea subset. The latter three subsets were shifted to minimize the χ2 between them and the low-latitude subset (see Supplementary Materials for details). Note that the low-latitude data show a warming for the past 15 Ma while the three other subsets exhibit cooling. Note also the data gaps around 110 and 210 Ma. The dotted vertical lines denote time intervals used for splicing the different combinations of subsets.

Figure 2

The linearly detrended and high pass filtered ML200 δ18O dataset (in red) for Fourier modes shorter than 49 Ma.

This is because we do not model processes on longer time scales. The simulated VO motion of the solar system in the galaxy (blue) has a secondary frequency modulation caused by the epicyclic motion of the solar system that generates slightly shorter VO periods around 130 Ma and 300 Ma and longer ones in between. Because the vertical potential changes adiabatically with the epicyclic motion, the vertical amplitude is larger when the period is longer. The shaded region denotes the 95% confidence range for the measured δ18O obtained from the finite number of data points in each bin and the variance in the data.

The model

To test whether the apparent δ18O oscillations are related to the galactic VO of the solar system, we need to construct a model and compare it with the above data. Such a model requires two components. One is integrating the galactic orbit, and in particular, the VO, while the second component is linking the distance from the galactic place to the terrestrial climate.

For the orbit integration, we assume that there is a negligible drop of the density with distance from the galactic plane, implying that the vertical oscillation is harmonic. A realistic density fall will typically introduce a 10% error2. Namely, the density fitted describes an average density experienced by the solar system, not the actual density at the galactic plane. We also assume that the density falls exponentially with radial distance, such that it can be written as: For an R = 3.5 ± 0.5 kpc decay length1, a 650 pc radial epicyclic motion of the solar system (RO) would give rise to a 17 ± 3% variation in the density (See §S7 for a summary of the kinematic parameters). This variation should be higher if any flaring of the disk is present. For comparison, the second largest modulation, due to the spiral arms, is expected to result in only a 5 to 10% density variation2728.

The orbit can then be integrated to give the VO and a RO, as described in Supplementary Materials S3. One of the more interesting aspects of this orbit is the coupling between the average disk density and the vertical oscillation. When the solar system is closer to the galactic centre, the density is higher and the vertical oscillation is both faster and with a smaller amplitude than when it is further out.

The main hypothesis in the model is that the average global temperature depends on the distance from the galactic plane. The lowest order expansion of any symmetric relation should be second order, that is, The constants will depend on the actual physical link between z and the terrestrial temperature, if such a link exists. Again we note that we do not assume any particular model to be responsible for this link, however, we assume for practical reasons that it is the cosmic ray climate link13 since it provides us with quantitative predictions for the actual temperature variations expected.

We expect the cosmic rays to diffuse out of the galactic plane, giving a drop with height. Since the cosmic ray halo is much larger than the typical 100 pc amplitude of the VO29, the lowest order expansion should be adequate. Because the oscillations are nearly harmonic, the T(z) relation associated with eq. 2 would then be close to sinusoidal, with an amplitude and frequency that vary with the RO oscillations.

As an example, if the cosmic ray density at the maximum amplitude z is 30% smaller than at the plane, it would translate into a 2 ± 1°C temperature variation30 (with the plane being colder). This would correspond to without ice-volume effects, or more with.

The free parameters include the present vertical location z and velocity v,0, the average density ρ0 and its radial gradient ρ0/R, the amplitude of the radial oscillation a and time of perigalacticon t, as well as the parameters A18 and B18 appearing in eq. 2. Since many of the parameters have independent observations, the best fit parameters in any viable model should be consistent with the observations. Note also that there is a degeneracy between and B18, as well as between the radial density gradient and the amplitude of the radial motion.

Results

We optimize these free model parameters by minimizing the residuals between the predicted model and the detrended “master” set (fig. 2). Because the data is noisy, with many local minima, we employ a genetic algorithm for minimization. The best fit gives an average VO of nearly 32 Ma, with a 175 Ma frequency modulation by RO. We use the bootstrap method on the different spliced datasets to estimate the statistical and systematic errors, as described in §S4 of the Supplementary Materials. For the exact fit parameters see Table 1.

Table 1

Summary of the paleoclimate derived galactic parameters

Parameter		Paleoclimate based value*	Astronomically based value
VO period		31.9 ± 0.6stat ± 0.6sys Ma†	30–42 Ma ‡
Density at the plane	Epicyclic center:	0.172 ± 0.006stat ± 0.006sys	—
	Locally:	0.208 ± 0.018stat ± 0.015sys	0.007−0.25 ‡
VO phase (last plane crossing)		3.9 ± 1.2stat ± 0.8sys Ma	2.8 ± 1.2 Ma§
RO Period		175 ± 9stat ± 12sys Ma	172 ± 7 Ma¶
RO Phase (next perigalacticon)		29 ± 17stat ± 8sys Ma	21 ± 7 Ma
Average δ18O range††		1.03 ± 0.06stat ± 0.11sys ‰

*The nominal values are obtained after having shifted the old geological chronology to the recent one45. The systematic error estimated is the sum of the systematic error obtained by fitting all data combinations of δ18O subsets and from the updated geological chronology. The statistical error is estimated using a bootstrap method.

†This value corresponds to the average oscillation period experienced by the solar system, where the epicyclic center is about 0.5 kpc further out than the present location.

‡Different kinematically based density estimates yield a wide range of periods, with the recent Hipparcos based determinations giving densities at the lower end of this range. See Supplementary Material S6 for a detailed discussion.

§Astronomical estimate for the last crossing is based on the solar system being 15–30 pc North of the galactic plane, moving northwards at 7 ± 1 km/sec46.

¶Using Hipparcos based determination of the Oort constants47, plus an offset due to the location of the epicyclic center.

††The average variation in δ18O between maximum warmth and maximum coldness, as obtained in the model fit, corresponds to ≤4°C.

The plot of δ18O as a function of the modeled vertical height z in the galactic disk (fig. 3) demonstrates that the average δ18O(z) can be adequately described with a quadratic form in z that is assumed in the model fit. The figure also shows that these δ18O oscillations are present in all latitudinal subsets, again confirming the robustness and global validity of the observations. The total δ18O variation is about 1‰, roughly 4°C, or less if the dynamics of ice caps accounts for some of the oscillations in the isotope signal. The impact of the secondary modulation of VO by RO, which explains the 30–40 Ma range of values around the primary 32 Ma frequency in the wavelet analysis (Supplementary Figs. S3–S5), is illustrated in fig. 4.

Figure 3

The “master” set δ18O data (Fig. 2) plotted against the modeled vertical location in the galaxy.

Normalization is relative to the maximal z0 of the recent maximal excursion around 20 Ma ago. The small orange dots are the actual 1 Ma binned data. The thick blue error bars are the averages of the combined data binned to 10 equal vertical bins. The red curve is a parabolic best fit: δ18O/[‰] = 0.42 − 0.09(z/z0) − 0.76(z/z0)2. The additional points are a coarser binning of the latitudinally separated subsets (colours are the same as in Fig. 1). Because of the poor coverage, the latitudinal data subsets cannot be high pass filtered. Instead, they have the low pass component of the “master” set removed. For legibility, an offset of −1‰ is applied. The fact that a similar vertical dependence appears in all 4 unrelated subsets indicates that the phenomenon is a global one.

Figure 4

A raster plot of the detrended 1 Ma binned master set of δ18O data.

The vertical axis spans the Phanerozoic. The horizontal axis is the time folded over a 32 Ma period. For convenience, the horizontal axis is duplicated and shifted by 32 Ma. The blue and red circles (connected by dashed lines) are the modeled galactic plane crossings (blue) and the maximal excursions from the plane (red), respectively. An unmodulated 32 Ma signal would appear as a vertical line over the entire 490 Ma interval. The apparent modulation due to the radial epicyclic motion (RO) of the solar system is expressed as sinusoidal variations of this line. Xs denote bins with insufficient measurements. The disk radii and colour correspond to the detrended and high pass filtered δ18O signal, as given by the scale on the right (in ‰).

Statistical significance of the 32 Ma signal

In order to assess the statistical significance of the 32 Ma signal, we adopt the Multitaper Method (MTM) often used when studying climate signals. The method and the specific toolkit used are described at length in ref. 31. It allows us to estimate the amount of structure at different spectral bands and to assess the statistical significance by estimating the amount of noise level as a function of frequency expected for the null hypothesis.

Fig. 5 describes the MTM spectrum. Because there is no particular preference for any noise model, we conservatively assess the noise level as a function of frequency under the assumption of an AR(1) model31. This modeling allows for more noise to occur naturally at low frequencies, ensuring that the significance is not artificially high.

Figure 5

The Multitaper Spectrum of the DML75 dataset, using 3 tapers.

Without a good model for the noise estimate, the noise level is conservatively estimated assuming an AR(1) model, while allowing for more noise at low frequencies than a “white” noise model would. See ref. 31 for a description of the method. The dominant low frequency peaks are at f ~ 0.006 cycle/Ma and f ~ 0.03 cycle/Ma which are also clearly evident in the wavelet and Fourier analyses (see §S2.5 of the supplementary material). Significant higher frequencies are also present. However, they are absent from the wavelet and Fourier analyses. It is therefore not clear whether they correspond to real climate variations or to artifacts in the dataset.

The spectrum reveals two statistically significant peaks at low frequencies, around 150 Ma and around 30 Ma, with a significance that is more than 99.9% and 99.99% respectively. If we assume a “locally white” model for the null hypothesis instead of an AR(1) process, the significance is 99.99% per frequency bin for both peaks. (The significance for pure “white noise” will be much higher).

Since we expect the galactic half period to be 30 to 42 Ma34, and since the frequency resolution is about 0.001 Ma−1, there are roughly 10 frequency bins in which the galactic signal can appear. Thus, the probability that the given frequency band will randomly yield a signal which is as large as that observed in any of the bins, is about 10−3.

Table 1 reveals that the nominal phases of the geological 32 Ma signal and the astronomical data are consistent within 1.1 Ma of each other (with typical errors of 1.2 to 1.5 Ma on each measurement). Thus, there is a about a 15% probability that the two unrelated signals will have a phase difference which is as small as the one observed or within the combined error. This implies that the probability that there will be a random fluctuation giving a galactic like oscillation (i.e., within 30 to 42 Ma period) and a phase which is consistent with the astronomical data is about 1.5 × 10−4.

The probability that this random fluctuation would also have a secondary modulation that is consistent both in phase and period with the radial epicyclic motion is smaller by at least an order of magnitude. Thus the statistical significance of the 32 Ma signal and its secondary modulation is of order .

Discussion

Given the consistency between the vertical and radial oscillations and the paleoclimate data, and the low probability that it could be mimicked by random fluctuations, we conclude with high confidence that the terrestrial temperature has a component which is quadratic in the distance from the galactic plane. Although this can be naturally explained through the cosmic ray climate link13, the observations by themselves do not prove it.

In addition, it should be noted that although a galactic driver can naturally explain a stable ~32 Ma cycle, there are terrestrial processes that could drive climate variations on the ~32 Ma time scale as well. The most prominent is probably mantle convection periodically producing plumes that result in large volcanic eruptions/igneous provinces. These eruptions will in turn add aerosols and carbon dioxide to the oceanic-atmospheric system32 and either cool or warm the climate. The influence of the Earth mantle connection could also drive changes or even reversals in Earth's magnetic field3334, which would modulate the atmospheric ionization. The advantage of the proposed galactic forcing over a terrestrial driver is that it will produce a gradual (sinusoidal) fluctuations as found for most of the δ18O record with a relatively steady periodicity, while volcanic forcing will more likely produce random abrupt perturbations followed by gradual relaxations to climate base levels (e.g., “volcanic winters”35). The non-gradual behaviour would also characterize other galactic forcings, such as periodic asteroid/comet bombardments triggered by the galactic tides36.

As a side note, there is also evidence for a ~ 60 Ma periodicity in the literature373839404142. For example, Meyers & Peter38 find a 56 ± 3 Ma period in the sedimentation record, which can be explained as an orogenic oscillator or as oscillatory dynamics in the mantle convection. Thus, these ~60 Ma periodicities are probably unrelated to the 32 Ma cycle discussed here, unless there is a very large north-south asymmetry relative to the galactic plane, and the 60 Ma cycle can be shown to be in phase with double the 32 Ma cycle discussed here.

The high accuracy of the estimated VO periodicity also enables estimating the density at the galactic plane, and may potentially point to an interesting conundrum. Our effective density (for an equivalent homogeneous disk) is at the radial center of the epicyclic motion, and at the present galactic radius (near perigalacticon, see Table 1). These values are inconsistent with the density of obtained by analyzing the dynamics of A & F stars in the Hipparcos catalogue4. Some of this inconsistency can be explained given that the different numbers do not reflect the same density. For example, the different model parameters may have changed over the past 500 Ma due to diffusion in parameter space, such that the 32 Ma reflects the average witnessed by the solar system over the past 500 Ma4344, and not the exact parameters today. In addition, we show in the Supplementary Materials that the Hipparcos based determination has a large unaccounted systematic uncertainty because the local stellar population is not dynamically relaxed. It is large enough to explain the whole inconsistency. It is worth to note in this respect that the GAIA star catalogue will be published in the near future. Since it will be much more extensive than the Hipparcos catalogue, both in terms of number of stars and their distance from the solar system, we expect under the above model to see convergence between the astronomically derived quantities and their respective geological counterparts. On the other hand, a more extensive geological database could provide additional evidence in the form of a second “frequency modulation” that corresponds to the density variations associated with the spiral arm passages. These variations are 2–3 times smaller than the RO and therefore harder to detect.

The density at the galactic plane should also be compared to the observed baryonic matter and the estimated halo dark matter densities, which only contribute and respectively (see ref. 4 and references therein). This discrepancy may imply the presence of a large unaccounted component that is not the usual halo dark matter.

Methods

Raw data

The first step in the present work was to expand the previously available database of δ18O and δ13C measurements for low-magnesium calcite shells of marine fossils for the Phanerozoic (last 542 Ma). This present database is expanded by almost 30% relative to the previous compilation26, to about 24,000 measurements (Table S1), with brachiopods, foraminifera and belemnites accounting for the bulk of the samples. It includes data from more than 20 new references that were published prior to July 2011 (see the 2nd sheet of Datafile S1). More information about the compilation is available in the Supplementary Materials. The database itself can be downloaded as Data File S1.

Homogenized time series

The first step in the data reduction process was to construct temporally homogeneous time series at equidistant 1 Ma intervals. This was carried out by Gaussian filtering the raw data of each habitat (tropical to high-latitude) separately. The second step was to derive the average δ18O offsets between the low altitude and each of the other habitats, by comparing the average δ18O over periods with sufficient data in each habitat pair. Because no habitat has a complete coverage of the last 490 Ma, the third step was to combine the data sets by splicing different habitats together, once the aforementioned offsets are applied to the time series. This step produced 5 different sets spanning the last 490 Ma. The different sets were all used in the analysis, with the variations giving a handle on possible systematic errors. The nominal set was the ML200 which includes the midlatitude data over the last 200 Ma, and the low-altitude data before that. It includes only one splice and avoids the recent low-latitude data that show atypical low δ18O values for the last 15 Ma, if compared to all the habitats. Although the raw data is based on the 2004 geological chronology, the time scale was adjusted to the new 2012 chronology45.

Author Contributions

The geochemical data analyses and compilation was carried out by A.P. and J.V. The astrophysical modeling was done by N.J.S. The manuscript was written jointly by N.J.S., A.P. and J.V.