## Abstract

Evidence from power spectra of deep-sea oxygen isotope time series suggests that the climate system of Earth responds nonlinearly to astronomical forcing by frequency modulating eccentricity-related variations in insolation. With the help of a simple model, it is shown that frequency modulation of the approximate 100,000-year eccentricity cycles by the 413,000-year component accounts for the variable duration of the ice ages, the multiple-peak character of the time series spectra, and the notorious absence of significant spectral amplitude at the 413,000-year period. The observed spectra are consistent with the classic Milankovitch theories of insolation, so that climate forcing by 100,000-year variations in orbital inclination that cause periodic dust accretion appear unnecessary.

Understanding the climate of the recent geological past is of importance, because finding out how Earth's environment has transformed the heat input from the sun into climate variations can help predict future global climate change. Most of our knowledge about climatic variations in the Plio-Pleistocene (the last 5.2 million years) is extracted from time series of oxygen isotope ratios, a proxy for global ice volume generally known as the δ^{18}O records (Fig. 1A). The records show that during the last million years, Earth has experienced at least 10 major glaciations, which according to the astronomical theory of the ice ages (1) are the consequence of secular variations in insolation caused by changes in Earth's orbital eccentricity, axial tilt, and longitude of perihelion (Earth's closest approach to the sun). The theory finds support in the fact that the spectra of the δ^{18}O records contain some of the same frequencies as the astronomical variations (2–4), but a satisfactory explanation of how the changes in orbital eccentricity are transformed into the ∼100-ky (1 ky = 1000 years) quasi-periodic fluctuations in global ice volume indicated by the data has not yet been found (5).

For instance, a fundamental difficulty is to understand the notable absence in the δ^{18}O data of a significant response to the 413-ky component of the orbital eccentricity, whose signal power is of the same order of magnitude as the 95-ky component. This has been called the “400-ky problem” by Imbrie and Imbrie (6).

One of the well-known features of the δ^{18}O time series from deep-sea sediments is the slight but evident change in the duration of consecutive glacial periods (7), which oscillates from about 120 to 80 ky. For instance, the interval between the last two interglacials is 120 ky, while 400 thousand years ago (ka) the interglacial interval was nearly 80 ky, with three successive interglacials occurring in less than 200 ky. This is reflected clearly in the spectrogram of Fig. 1B (moving window Fourier transform) of Ocean Drilling Program (ODP) Site 806 (8), a representative paleoclimate record. The power of the eccentricity band (∼0.01 cycles/ky) in the spectrogram occurs along an approximately sinusoidal curve (as indicated by the dashed curve overlay), apparently because of periodic switching of main power from ∼120 ky to ∼80 ky every 400 ky or so. This fact suggests frequency modulation (FM) of the ∼100-ky eccentricity signals by the longer-period, 413-ky component. Such frequency modulation is entirely similar to electronic modulation of a high-frequency carrier by a low-frequency modulating signal, as used in FM radio and television broadcasting (9).

In addition to the characteristic sinusoidal shape of the spectrogram, frequency modulation of a single frequency carrier widens the spectral band by the addition of sidebands. Sidebands are spectral peaks distributed symmetrically on both sides of a carrier's peak at intervals equal to integer multiples of the modulating frequency (10). The appearance of orderly distributed sidebands in the power spectra of all the δ^{18}O data analyzed here further supports the presence of frequency modulation. For instance, the power spectrum of Site 806 (Fig. 1C) shows evidence of frequency modulation as the observed 75- and 123-ky peaks coincide with the predicted sidebands of a 95-ky carrier frequency-modulated by a 413-ky signal, and the 85-, 107-, and 143-ky peaks coincide with the predicted sidebands of a 95-ky carrier frequency-modulated by the 826-ky subharmonic of the 413-ky signal (11).

Another important prediction of the FM hypothesis is that the power spectrum of an FM signal does not contain the modulating signal as a spectral peak, but as the separation interval between sidebands (9), a fact that is entirely consistent with the well-known absence of a distinctive 413-ky spectral peak in paleoclimate data.

It is important to point out that the spectra in Fig. 1are obtained from an orbitally tuned (12) time series. Orbital tuning is routinely applied to raw δ^{18}O time series in order to establish a common time scale for all records regardless of location and sedimentation rate differences. However, the “clock” used to develop the chronology is the rhythm of any of the orbital periodicities, including eccentricity. Hence, it is possible that orbital tuning introduces bias or even spectral peaks not present in the data, and thus it is useful and prudent to analyze and compare both tuned and untuned records. In their original form the orbitally untuned records are of δ^{18}O variation versus depth; the depth scale is converted to time using known published ages (determined by biostratigraphic or magnetic methods, or both) at one or two anchor points (13), and assuming a constant sedimentation rate for the rest of the sediment core. In spite of its obvious drawbacks, this assumption produces spectra similar enough to the tuned data to permit easy cross-identification of spectral peaks. As shown in Figs. 2 through 4, sidebands separated from the main spectral power as predicted by the FM hypothesis are a common feature in both tuned and untuned spectra from different sites, and as discussed below (Fig. 4), the effect of tuning appears to reduce, rather than to increase, the number of recognizable spectral peaks in the eccentricity band.

Since what follows hinges on the interpretation of power spectra, it is important to note that because of the strong age dependence of the eccentricity band (Fig. 1B), resolving individual spectral peaks depends critically on the window length chosen to calculate the spectrum. This becomes an important issue, because the apparent presence of a single power peak at 100 ky (Fig. 3) in several δ^{18}O recordings has been interpreted as due to variations in Earth's orbital inclination, which force the planet to periodically cross a ring of extraterrestrial dust, thus causing the ice ages of the last million years (14). This conclusion is, however, based on power spectra computed from time windows spanning just the last 600 ky, an interval half the duration of typical eccentricity signals, which have significant power from 0 to about 1200 ka (Fig. 1B). As shown in Fig. 3, a 600-ky window limits the resolving power in the eccentricity band so much that important details are lost, and it is thus plausible that the presence of the 100-ky single spectral peak is due simply to merging of the 95-, 107-, and 125-ky eccentricity-related spectral peaks because the time window is too short, rather than from periodic changes in orbital inclination (15, 16).

If frequency modulation indeed occurs, it should be possible to replicate the observed variation in ice age duration seen in the δ^{18}O records using an appropriate FM formula. A simple mathematical form of a frequency-modulated signal is(1)where Ω = 2π*f*
_{c}, ω = 2π*f*
_{m} and *f*
_{c}and *f*
_{m} are the frequencies of the carrier and the modulating signals respectively, ɛ = Δ*f*
_{c}/*f*
_{m} is the modulation index (17), and Δ*f*
_{c} is the frequency deviation, a measure of the bandwidth produced by the modulation (9, 10). Since the predominant frequency shifts from ∼120 to ∼80 ky, which is also the observed shift in ice age duration, ɛ ≈ 0.7 (moderate modulation) for *f*
_{m} = 1/413 cycle/ky.

Trigonometric functions of the form in Eq. 1 can be represented by an infinite series of Bessel functions of increasing order, such as
(2)which clearly indicates that the Fourier spectrum of a FM signal such as given by Eq. 1 is composed of a central spectral peak at the carrier frequency*f*
_{c} and additional spectral peaks (sidebands) at the frequencies *f*
_{c} ±*nf*
_{m} for *n* = 1, 2, 3 … , so that the modulating frequency *f*
_{m} = ω/2π appears in the power spectrum only as the interval between spectral peaks, but not as a spectral peak itself (9). A plausible model to simulate the δ^{18}O records should thus contain the 125-, 100-, and 95-ky harmonics as carrier signals and the 413-ky and, as explained below, a 826-ky subharmonic as the modulating signals. The following model is possibly the simplest
(3)where *t* is the time in kiloyears. The constants *a* through *c *are adjustable parameters, with *a* typically three times greater than *b* or *c*, reflecting the relative strength of the eccentricty components. The term ɛ′ is the modulating index for the subharmonic, and since ɛ ∼ 1 then ɛ′ ∼ 2, since*f*
_{m} is halved. In spite of the simplicity of Eq. 3, its spectrum (Fig. 4), reproduces every spectral peak in the power spectra of the tuned and untuned records, including the prominent combination tone at 107 ky (Figs. 2 and 4) which occurs only if the additional modulation by a 826-ky subharmonic is included in the model, as in Eq. 3. Frequency modulation as described by Eq. 3 is nonlinear in the sense that the effect produced by two modulating frequencies is not equal to the sum of the signals produced separately by each (18).

The 107-ky spectral peak is here interpreted as a combination tone (1/107 = 1/95 – 1/826) produced by nonlinear response of the climate system. Although the precise mechanism by which this may occur in the real climate system is not resolved here, an explanation is attempted below. On the other hand, Eq. 3 uses superposition of each modulated carrier to construct the total signal, and so in this sense the response is assumed linear.

Figure 4 shows that even when a spectrum is distorted by tuning, the sidebands can still be identified and the general structure of the spectrum understood. The distortion occurs because tuning in this case targeted the eccentricity frequencies (19), thus enhancing the amplitude of the 95-ky eccentricity carrier to the expense of its sidebands, as becomes apparent when comparing Fig. 4B with Fig. 4A. Accordingly, to compensate for tuning and simulate the untuned spectra required just a 20% increase in ɛ′ in Eq. 3.

It should be noted that the model spectra shown in Fig. 4 account not only for all the periods present, but also for those not present. As mentioned before, FM theory predicts that the modulating frequencies appear only as intervals between sidebands, and not as spectral peaks. Hence, if the role of the 413-ky eccentricity component is indeed to frequency-modulate the higher frequency components, we can expect to find signal power at 413 ky in all the δ^{18}O records, though not in the form of a spectral peak, but rather as the interval between sidebands (clear examples are Figs. 1C, 3A, and 4). Thus, the apparent absence of a 413-ky spectral peak in Earth's climate response, also known as “the 400-ky problem” (6), may not be a problem after all, but the logical consequence of the fact that the 1/413 ky^{−1} signal acts as a frequency modulator.

The simulated FM signal Φ_{FM}(*t*), shown in Fig. 5 compared to selected δ^{18}O time series, successfully reproduces the subtle change in time interval between consecutive ice ages. For times before 1 Ma, the fit deteriorates, because the period of the ice ages becomes dominated by the 41-ky obliquity signal and the eccentricity carriers gradually disappear (Fig. 1, A and B).

How does the climate system frequency-modulate the astronomical forcing? A natural place to search for an analogy is the mechanism of the resonant circuit of an electrical oscillator that generates FM signals as its capacitance or inductance are made to vary (9). The simplest analogue paleoclimatic mechanism is a conceptual representation of the climate system as a (potentially) resonant oscillator (20), where the linear extent of the ice sheet, *L*, satisfies an equation of the form(4)where Ω = (1/*C*
_{L}
*C*
_{T})^{1/2}is the natural angular frequency of the climate model,*C*
_{T} the ocean's thermal overturning time, and*C*
_{L} the ice cap's relaxation time. There are at least two ways to include external forcing in Eq. 4: with an additive term of the form *F* cos ω_{0}
*t*, which introduces the possibility of resonance when ω_{0} → Ω, and/or with a time-varying parameter (parametric excitation), assuming that Ω changes with the dimension of the ice sheet (21) through C_{L}, which is proportional to a critical dimension of the ice sheet (22).

If *C*
_{L} is of the form*C*
_{L}(1 + Δ*C*
_{L}/*C*
_{L} cos ω*t*) where Δ*C*
_{L} ≪*C*
_{L} and ω is the modulating frequency, Eq. 4becomes(5)which is the Mathieu equation (23), whose periodic solutions are frequency-modulated sinusoids (24) of the form(6)where A and φ are integration constants. The solution to Eq. 6 is of the same form as each of the terms of Eq. 3used to simulate the δ^{18}O records. Nonlinear forms of Eq. 5 are known to develop instabilities, bifurcations, and period-doubling cascades into chaos (25).