2.3.1 The bispectrum

In contrast to spectral analysis (e.g. using the fast Fourier transform), which gives the distribution of variance of a sinusoidal signal with frequency, bispectral analysis describes the distribution of non-sinusoidality with frequency (King, 1996). The skewness of a cycle is described by the real part of the bispectrum, whereas the asymmetry of a cycle is deconvolved in the imaginary part of the bispectrum. The bispectrum shows nonconservative relative energy exchanges among frequencies of a single time series. Energy transfers in the bispectrum computed on the LR04 stack are expressed in ‰³ kyr², which constitutes a statistical measure of energy. (We refer to Sect. 2.4.1 for an explanation of conservative net energy exchanges and Sect. 4.2 for a discussion of the scaling of conservative to absolute energy exchanges). Energy transfers can only occur, and be analysed by the bispectrum, if both the frequencies (f) and phases (Phi, φ, in radians) of these non-sinusoidal cycles are coupled: i.e. $f_{\mathrm{1}}+f_{\mathrm{2}}=f_{\mathrm{3}}$ and ${\mathit{\phi }}_{\mathrm{1}}+{\mathit{\phi }}_{\mathrm{2}}={\mathit{\phi }}_{\mathrm{3}}$, respectively (Hagelberg et al., 1991). Thus, for any possible combination of three frequencies, the bispectrum assesses whether there is a coupling, and if so, whether energy exchanges occur. The transfer of energy in these so-called triad interactions is nonlinear, because the changes in cycle amplitudes (A) do not sum and have a currently unknown, and probably variable, coupling through time (i.e. $A_{\mathrm{1}}+A_{\mathrm{2}}\ne A_{\mathrm{3}}$; see Sect. 4.2). We note that the meaning of “cycle”, “frequency” and/or “periodicity” is different when referring to either the time or frequency domain. For example, a “single” cycle in a time series is often composed of multiple frequency components, which can be deconvolved by (bi)spectral analysis. However, for textual purposes, we use these words interchangeably, regardless of their reference to a specific domain.

The bispectrum is defined by Eq. (4):

where E[…] is the ensemble average of the triple product of complex Fourier coefficients C at the difference frequencies f₁, f₂, and their sumf₃ (i.e. $=f_{\mathrm{1}}+f_{\mathrm{2}}$), and the asterisk indicates complex conjugation (Hasselmann et al., 1963). In this study, we focus on the imaginary part of the bispectrum (hereafter often referred to as “bispectrum”, unless indicated otherwise), following studies on ocean waves (e.g. Herbers et al., 2000; de Bakker et al., 2015), because according to bispectral theory most energy transfers are associated with asymmetric cycle shapes (i.e. wave forms), and because the climate cycles of the Middle and Late Pleistocene have high amplitudes (Fig. 1) and are highly asymmetric (Fig. 2). Furthermore, conservative net energy exchanges between the coupled frequencies that are located in the imaginary part of the bispectrum can be computed if we assume a simple coupling coefficient between frequencies. Prior to bispectral analysis, we detrend the data and subsequently apply a tapering function using a Hann (a.k.a. Hanning) window that also mitigates spectral leakage (Hagelberg et al., 1991) and then multiply the data by an energy correction factor to adjust for the change in amplitude that results from the windowing.

2.3.2 Interpreting the bispectrum

The x and y axes in the bispectrum, which correspond to f₁ and f₂, respectively, are mirror images of each other and share a symmetry axis at x=y (Fig. 3) (Hasselmann et al., 1963). Therefore, only the positive, one-eighth part of the bispectrum is depicted, which is subdivided into 15 zones (see Sect. 2.4.2, Table A1). For triad interactions, only two outcomes exist: either two frequencies (by definition, f₁ and f₂) gain energy (though not necessarily in equal measures), and one frequency (i.e. f₃) loses energy, or f₁ and f₂ lose energy (ditto), and f₃ gains energy. Frequencies f₁ and f₂ are the so-called difference frequencies, while frequency f₃ is referred to as a sum frequency. Frequencies often participate in multiple interactions, and depending on the interaction, they act as either a difference frequency or a sum frequency. A particular frequency can thus receive energy through one interaction and simultaneously lose energy in another interaction.

Figure 3Bispectral zonation scheme. Bispectra of climate cycles can be subdivided into 15 zones that reflect unique combinations of three frequency bandwidths. Orange, purple, and blue lines represent the main astronomical frequencies of the eccentricity, obliquity, and precession cycles, respectively (this colour coding is consistent throughout the paper). The grey shaded areas represent suborbital periodicities. Difference frequencies f₁ and f₂ are plotted along the along the horizontal and vertical axes, respectively. Sum frequencies f₃ (i.e. f₁+f₂) are depicted as diagonal lines in the bispectrum. See also Table A1 for a list of bispectral zones and Table A2 for a list of triad types among the main astronomical frequencies in the bispectrum.

Download

To help the interpretation of the bispectrum, we depict the main astronomical frequencies in the bispectrum with coloured lines (Fig. 3). Vertical and horizontal lines correspond to the difference frequencies (f₁ and f₂, respectively), and diagonal lines correspond to the sum frequency (f₃). Junctions between lines highlight the locations in the bispectrum where interactions occur between three periodicities (i.e. triad interactions) of which at least one is an astronomically paced climate cycle, which can also interact with itself (Fig. 3; see also Table A2). No frequency axis is associated with sum frequency f₃, but its value can be read off at the crossing points with the x and y axes of the bispectrum or by summing the x and y coordinates of the difference frequency axes at any point along the diagonals (Fig. 3).

We follow the convention and define a triad interaction as negative or positive if f₃ either loses (blue colours) or gains (red colours) energy, respectively (Fig. 4). In written form, energy gains at a particular frequency are marked by an upward-pointing arrow (↑) and losses by a downward-pointing arrow (↓). Bispectral notation of triad interactions are given as B_f(f₁, f₂, f₃) in the frequency domain (Myr⁻¹) and as B_p(p₁, p₂, p₃) in the periodicity domain (kyr). For palaeoclimatological purposes, we will mainly focus on the periodicities for which the corresponding frequencies can easily be looked up in Table A2. We refer to the relevant part of the bispectrum as “Re” (i.e. real) and “Im” (i.e. imaginary), which are given in superscript (B^Re and B^Im, respectively). Thus, if we consider a negative triad interaction, located in the imaginary part of the bispectrum, between p₁, a 40 kyr obliquity-paced cycle (i.e. f₁=25.0 Myr⁻¹), and p₂, a ∼95 kyr eccentricity-paced cycle (i.e. f₂=10.5 Myr⁻¹), that results in a loss of energy at p₃, the ∼28 kyr periodicity (i.e. f₃=35.5 Myr⁻¹), this is given as $B_{\mathrm{p}}^{\mathrm{Im}}$(40↑, 95↑, 28↓) (i.e. $B_{\mathrm{f}}^{\mathrm{Im}}$(25.0↑, 10.5↑, 35.5↓)) (Fig. 4a).

Figure 4Bispectra of Pliocene and Pleistocene climate cycles. (a) Imaginary part of a bispectrum of the Middle and Late Pleistocene “∼110 kyr world”. Age range =668 to 0 ka. Resampling resolution =1 kyr⁻¹. Window length =668 kyr. Block length =668 kyr. Number of blocks =1. Number of merged frequencies =1. Frequency resolution =1.50 Myr⁻¹. Degrees of freedom = 2. Min value = −598 ‰³ kyr². Max value =51 ‰³ kyr². (b) Imaginary part of a bispectrum of the mid-Pleistocene transition “∼80 kyr world”. Age range =1318 to 650 ka. Resampling resolution =1 kyr⁻¹. Window length =668 kyr. Block length =668 kyr. Number of blocks =1. Number of merged frequencies =1. Frequency resolution =1.50 Myr⁻¹. Degrees of freedom =2. Min value $=-\mathrm{57}$ ‰³ kyr². Max value =73 ‰³ kyr². (c) Imaginary part of a bispectrum of the Pliocene and Early Pleistocene “40 kyr world”. Age range =3000 to 1000 ka. Resampling resolution =1 kyr⁻¹. Window length =2000 kyr. Block length =2000 kyr. Number of blocks =1. Number of merged frequencies =3. Frequency resolution =1.50 Myr⁻¹. Degrees of freedom =6. Min value $=-\mathrm{21}$ ‰³ kyr². Max value =23 ‰³ kyr². Note the different scaling of the z axes between panels (a), (b), and (c). These bispectra are computed on the LR04 stack (see methods section) (Lisiecki and Raymo, 2005a).

Download

We note that within this study, we almost exclusively focus on the imaginary part of the bispectrum. This approach contrasts with previous palaeoclimate investigations that used the bicoherence spectrum (see, e.g. Hagelberg et al., 1991; von Dobeneck and Schmieder, 1999; Wara et al., 2000; Huybers and Curry, 2006; Da Silva et al., 2018). In short, bicoherence is a measure of the coherency between the real and imaginary parts of the bispectrum (i.e. where both parts have strong interactions) and is marked by positive values only (see, e.g. Elgar and Sebert, 1989; de Bakker et al., 2014). Conversely, the real and imaginary parts of the bispectrum transform skewed and asymmetric cycle geometries into their composite frequency components, respectively, and both parts are characterised by positive as well as negative interactions.

2.3.3 Quantifying geometries using the bispectrum

Adopting the bispectral method to extract cycle geometries (Figs. 2, S2, and S3) (Elgar, 1987), skewness and asymmetry are computed from the biphase. The biphase reflects the ratio between the real and imaginary parts of the bispectrum. The biphase is defined by Eq. (5):

where n and l range from 1 to the Nyquist frequency N, with n > l and $n+l\le N$ (Elgar, 1987), and f_p refers to the primary frequency (a.k.a. natural or resonance frequency, or first harmonic). We note that positive asymmetry values correspond to a dominance of negative interactions (i.e. a loss of energy at f₃) in the imaginary part of the bispectrum, and vice versa, negative values for this geometry are indicative of an overall dominance of positive interactions (i.e. a gain of energy at f₃).

2.4.1 Integration of the spectrum and bispectrum

To quantify power spectral density, we apply a standard integration of the power spectrum. To obtain conservative net energy transfers per frequency (defined by the nonlinear source term S_nl) and hence to extract the total gain or loss in energy per climate cycle over a specific time interval (i.e. window), we integrate over the imaginary part of the bispectrum and multiply it by a coupling coefficient. For palaeoclimate time series, no coupling coefficient has been determined previously. Therefore, we make minimum assumptions and use a coupling coefficient that only corrects for frequency $W_{(f_{\mathrm{1}},f_{\mathrm{2}})}=(f_{\mathrm{1}}+f_{\mathrm{2}})$. A comparable correction for frequency is also part of the Boussinesq scaling that is used for ocean waves (e.g. Herbers and Burton, 1997; Herbers et al., 2000). The correction for frequency ($W_{(f_{\mathrm{1}},f_{\mathrm{2}})}$) insures energy conservation during triad interactions, because more energy is exchanged among the higher frequencies. Furthermore, this correction allows for qualitative interpretations of conservative net energy exchanges across the spectrum. If a more appropriate coupling coefficient for palaeoclimatic purposes were to be deduced, it could also scale the conservative net energy transfers of the bispectrum (already corrected using S_nl during integration) to absolute energy transfers that can be compared to (changes in) power spectral density (see Sect. 4.2). We express the integral of S_nl in terms of an integration, or summation, over the positive quadrant of the bispectrum alone or equivalently in sum and difference interactions following Eq. (6):

where the sum contributions are expressed by Eq. (7):

and the difference contributions are expressed by Eq. (8):

The sum contributions are obtained by integrating diagonally over the bispectrum, and the difference contributions are obtained by integrating along vertical and horizontal lines, perpendicular to the x and y axes, respectively (Fig. 3). We track the evolution of net energy transfers through time by integrating over bispectra that are determined for consecutive and partially overlapping windows of the time series (i.e. a moving window) (Figs. 5, 6, and S3). Integration can be performed including all frequencies (i.e. over the entire imaginary part of the bispectrum) or over specific zones only, in the event that subsets of interactions need to be further examined (Fig. 3, Table A1) (de Bakker et al., 2015, 2016). After integration over the bispectrum, either totally or zonally, energy transfers computed on the LR04 stack are expressed in ‰³. Blue colours indicate a loss of energy at a specific frequency; red colours indicate a gain.

Figure 5(a) Input → (b) “black box” climate → (c) output. (a) Power spectral density of mean summer insolation (21 June to 21 September) at 65^∘ N (Paillard et al., 1996; Laskar et al., 2004). (b) Conservative net energy transfers (i.e. after correction using the coupling coefficient) during the Pliocene and Pleistocene, computed by integrating over the entire imaginary part of the bispectrum (see methods section). Input data are the resampled LR04 stack (Lisiecki and Raymo, 2005a). (c) Power spectral density of the LR04 stack (Lisiecki and Raymo, 2005a). For all three panels, the window length is 668 data points (668 kyr), and the step size between partially overlapping windows is 50 data points (50 kyr). We used no blocks and did not merge frequencies. These settings yield a frequency resolution of 1.50 Myr⁻¹ and 2 degrees of freedom.

Download

2.4.2 Spectral and bispectral zones

To obtain more insight on the energy that is exchanged among eccentricity, obliquity, and precession cycles, we integrate the spectra and bispectra over separate zones (Fig. 3, Table A1). For bispectra, such a zonation approach was first applied in research concerned with nearshore ocean waves to distinguish between infragravity wave and sea-swell wave frequencies (de Bakker et al., 2015), and is adapted here to investigate the most important climate cycle bandwidths of the Pliocene and Pleistocene. The boundaries between the climate cycle zones are (arbitrarily) defined at frequencies of 17.8, 35.6, and 80.0 Myr⁻¹ (i.e. periodicities of ∼56.2, ∼28.1, and 12.5 kyr). In total, 15 bispectral zones are defined (Fig. 3, Table A1), each of which represents a unique combination of three main groups of astronomically paced climate cycles and suborbital periodicities (i.e. supraorbital frequencies) that can interact and exchange energy with one other. We note that many periodicities without primary astronomical origin are included in these “eccentricity”, “obliquity”, and “precession” zones. We focus here on the results of the eight most important zones, namely those involving astronomically paced climate cycles (Figs. 7, 8, 9, S4, and S5). To better highlight the roles of all nonlinear interactions involving eccentricity, obliquity, or precession, we recombine (i.e. sum) the individual zones (Fig. 10).

3.3.1 Power spectral density of insolation

The time-evolutive power spectra of insolation are marked by high spectral power at the precession and obliquity periodicities, and the near absence of spectral power at the eccentricity periodicities (Figs. 5a, 6, and 7a) (Hays et al., 1976; Berger, 1977). Furthermore, it is characterised by amplitude variability in spectral power at these periodicities, which is governed by the long-term (eccentricity) modulations of obliquity and precession. These amplitude modulations have durations of ∼180 kyr (not detected here due to window length) and 1.2, 2.4, and 4.8 Myr for obliquity, and 405 kyr and 2.4 Myr for eccentricity-modulated precession (Figs. 5a and 7a). Concurrent with the MPT, we document the sharpest decrease in total spectral power of the entire Pliocene and Pleistocene, which is largely due to a decrease in precession power from maximum values (Fig. 7a) (Laskar et al., 2004). This is indicative of a reduction in the perturbation of the Earth's system that results from a smaller difference in insolation forcing between the hemispheres (i.e. reduced seasonality).

Figure 7Integration of the spectra and bispectra per astronomical bandwidth (top right). (a) Totally and zonally integrated spectra of mean summer insolation (21 June to 21 September) at 65^∘ N. (b) Totally and zonally integrated bispectra of the LR04 stack. (c) Totally and zonally integrated spectra of the LR04 stack. Computational settings are as in Fig. 5.

Download

3.3.2 Energy transfers among climate cycles

Figures 5b and S3 depict the conservative net energy transfers among (astronomically forced) climate cycles as present in the LR04 stack. The net energy transfers show how energy is redistributed over the power spectrum among astronomically paced and non-astronomically paced climate cycles. They reveal strong increases in both energy gains (in red) and losses (in blue) toward the present, and from ∼1000 ka onward several consecutive shifts in energy transfers from lower to higher periodicities occur (Fig. 5b). Energy gains are highly concentrated at the 40 kyr periodicity from ∼2500 ka onward and are marked by a modest 1.2 Myr beat in amplitude (Figs. 5b and 7b). A second location of strong energy gains is the ∼80 kyr periodicity from ∼1000 ka onward (Fig. 6b), which at ∼700 ka shifts (Fig. 5b) to a dominant gain for the 95 kyr component (Fig. 6a). Further energy gains are located at a ∼200 kyr periodicity from ∼750 ka onward, a 50 to 60 kyr periodicity between ∼1500 and 500 ka, and (very weak) at the 24 kyr periodicity from ∼650 ka onward (Fig. 5b). Energy losses are also observed in the total integration over the bispectra (Figs. 5b, 6, and 7b). However, these losses are generally of smaller amplitude per cycle and are distributed across a greater number of periodicities, in comparison to the highly concentrated energy gains (Figs. 5b and 6). Despite this more dispersed pattern of energy losses, we highlight the more conspicuous decreases: the 40 kyr periodicity from ∼550 ka onward, the 30 to 36 kyr periodicities between 1500 and 1000 ka, the 30 to 26 kyr periodicities between ∼2500 and ∼2000 ka, the ∼28 kyr periodicity from ∼800 ka onward, the 24 kyr periodicity between ∼1000 and ∼800 ka, the 22 kyr periodicity from ∼700 onward, the 19 kyr periodicity from 2000 ka onward, and the 15 kyr periodicity intermittently from ∼1000 ka onward (Fig. 5b).

Figure 8Disentangling “black box” climate. Conservative net energy transfers during the Pliocene and Pleistocene over specific zones in the imaginary part of the bispectrum (see methods section). Computational settings are as in Fig. 5b. (a) Zone 1. (b) Zone 2. (c) Zone 3. (d) Zone 4. (e) Zone 5. (f) Zone 6. (g) Zone 7. (h) Zone 8. See also Fig. 3 and Table A1.

Download

3.3.3 Power spectral density of climate cycles

Time-evolutive power spectral analysis of the LR04 record is characterised by high-amplitude variability at the 40 kyr obliquity-paced cycle, from approximately 4000 ka onward (Figs. 5c, 6, and 7c). In addition, distinct ∼80 and ∼110 kyr cycles of rapidly increasing amplitude are documented from ∼1200 ka onward, marking the MPT (Fig. 6b) and post-MPT (Fig. 6a) intervals, respectively (Lisiecki and Raymo, 2005a, 2007). We note that the 125 kyr component of the ∼110 kyr signal in the power spectra (Fig. 5c) is not present in the integration over the bispectrum (Fig. 5b) A relatively weak ∼200 kyr cycle is present from ∼750 ka onward (Fig. 5c). No strong precession response is present (Raymo and Nisancioglu, 2003; Raymo et al., 2006). The 40 kyr cycle in δ¹⁸O partially mimics the 1.2 Myr amplitude modulation of the obliquity cycle, especially from ∼3000 ka onward (Fig. 5c) (Lourens and Hilgen, 1997; Lisiecki and Raymo, 2007). We draw attention to a subtle “arc” in spectral power (absent in bispectral power), which starts with a strong 405 kyr cycle, at ∼3500 ka, that becomes gradually shorter in duration until it resonates with the 95 kyr eccentricity modulation of precession, at ∼2650 ka, and then culminates in several weaker 70 to 80 kyr cycles at ∼2000 ka (Fig. 5c; see also Fig. 1b) (Lisiecki, 2010; Meyers and Hinnov, 2010; Viaggi, 2018).

Figure 9Conservation of energy in triad interactions located in the imaginary part of Pliocene and Pleistocene bispectra. Conservativity of (a) the entire bispectrum (Fig. 5b) and (b) Zone 1, (c) Zone 2, (d) Zone 3, (e) Zone 4, (f) Zone 5, (g) Zone 6, (h) Zone 7, and (i) Zone 8, of the bispectrum (Fig. 8). Energy gains or losses of f₁ and f₂ (dashed lines) are computed using Eq. (8), those of f₃ (single coloured lines) following Eq. (7), and their differences (black lines) by Eq. (6).

Download

3.4.1 Zonal energy transfers

In Zone 1, from 1000 ka onward, we document a modest exchange of energy between the (eccentricity-related) periodicities in the range from ∼200 to 50 kyr (Fig. 8a). Notably, from ∼500 ka onward, periodicities of ∼80 to 50 kyr predominantly lose power, whereas those between ∼200 and ∼80 kyr mainly gain energy. Stronger exchanges are documented in Zone 2 (Fig. 8b). Especially the 40 kyr obliquity-paced cycle gains energy between 3500 and 900 ka in this zone, after which it loses energy from 900 ka onward. Also, the 50 to 60 kyr periodicity loses energy, especially from ∼600 ka onward. Between ∼1100 and ∼800 ka, a single periodicity of ∼87 kyr gains power. At ∼750 ka, this periodicity splits into an ∼80 kyr periodicity and a 95 kyr eccentricity-paced cycle, which both gain energy. In general, interactions in Zone 2 change from (weakly) positive (B^Im(E↓, E↓, O↑)) to negative (B^Im(E↑, E↑, O↓)) during the MPT (∼1000 ka), indicative of a change from escalating energy to higher frequencies, to a cascading from higher frequencies. We note that the energy exchanged through B^Im(E, E, O) interactions (Zone 2) is almost of equal amplitude to that of all other interactions combined. Zones 3 to 5 are mainly marked by increases in energy at the 40 kyr cycle, roughly from 2500 ka onward (Fig. 8c to e). A smaller gain is documented in these zones (especially in Zone 3) at periodicities between ∼70 and 50 kyr from ∼1300 ka onward, and in Zone 3 at the 95 kyr cycle from ∼650 ka onward. Furthermore, we document decreases in energy at the ∼30 to 27 kyr periodicities, the ∼28 to 22 kyr periodicities, and the main precession periodicities of 24, 22, and 19 kyr, for Zones 3, 4, and 5, respectively. The switchover during the MPT, at about 750 ka, of energy loss at the 24 kyr cycle to a loss at the ∼28 kyr periodicity is associated with B^Im(O, E, P) interactions (Zone 4) only. Zones 6 and 7 show weaker energy exchanges among the main astronomical periodicities. In Zone 6, the ∼27 and 24 kyr cycles show a small gain between ∼1000 and 700 ka, and from ∼700 ka onward, respectively. We document only very minimal direct fuelling of eccentricity-paced climate cycles by precession-paced climate cycles in this zone. Energy is lost in Zone 6 from 1500 ka onward at the 19 kyr cycle. In Zone 7, relatively weak gains are present at the main obliquity and precession cycles from 1500 ka onward. Energy losses of Zone 7 are located mainly at sub-precession periodicities of ∼15 to < 13 kyr (Fig. 8f and g). No energy exchanges of similar magnitude are documented in Zone 8 (Fig. 8h).

Figure 10Summed zonal integrations of the imaginary part of Pliocene and Pleistocene bispectra. Net energy transfers are computed by (a) summing Zones 1, 2, 3, 4, and 6 for triad interactions corresponding to the “eccentricity” bandwidth, (b) summing Zones 2, 3, 4, 5, and 7 for triad interactions corresponding to the “obliquity” bandwidth, and (c) summing Zones 4, 5, 6, 7, and 8 for triad interactions corresponding to the “precession” bandwidth. Panel (d) shows the corresponding summed zonal conservation of energy, with eccentricity in orange, obliquity in purple, and precession in blue.

Download

Overall, the zonal integrations are marked by negative triad interactions and thus reveal a cascade of energy from the (sub)precession bandwidth to the obliquity and (ultimately) eccentricity bandwidths through several successive interactions (Fig. 8). These analyses further show that many periodicities fulfil a dual role that often remains stable through time, namely they serve simultaneously and persistently as energy provider and receiver. For example, the 24 kyr precession-paced cycle gains energy through B^Im(P, E, P) and B^Im(P, O, P) interactions, as documented in Zones 6 and 7, meanwhile losing energy through B^Im(O, O, P) and B^Im(O, E, P) interactions, as can be observed in Zones 4 and 5. As an exception to this persistent behaviour of many lower periodicities in transferring energy to higher periodicities, the role of the 40 kyr obliquity-paced cycle evolves but only in B^Im(E, E, O) interactions (Zone 2). In Zones 3, 4, 5, and 7, the 40 kyr obliquity-paced climate cycle (predominantly) partakes in triad interactions as a difference frequency, and it continuously receives energy, mainly from precession. However, during the MPT (from ∼900 ka onward) in Zone 2, when the 40 kyr periodicity partakes in interactions solely as a sum frequency (f₃) of two eccentricity components (f₁+f₂), it starts providing energy to the eccentricity bandwidth. The net (i.e. total) result of these opposing roles for the 40 kyr obliquity-paced cycle after the MPT is that energy losses overtake the gains (Figs. 5b and S3). However, the initial decrease of the integrated “obliquity” zone, between ∼1000 and 600 ka, is linked to the ∼28 kyr periodicity, not the 40 kyr periodicity (Fig. 5b).

3.4.2 Energy conservation

Computations of energy conservation, hereafter also referred to as “conservativity”, indicate that both the total as well as the zonal energy gains in triad interactions are commensurate with losses, given our assumed coupling coefficient (Fig. 9). The dominance of negative interactions (i.e. positive values for f₁+f₂ and negative values for f₃) is indicative of a transfer within the Earth's climate–cryosphere system of insolation energy to lower frequencies. Energy is not well conserved in triad interactions documented by Zone 8. This is probably the result of the low number of interactions located in Zone 8 (Fig. 4) and the small area of this zone (Fig. 3).

3.4.3 Energy exchanges involving precession-, obliquity-, or eccentricity-paced climate cycles

To highlight the separate roles of the precession, obliquity, and eccentricity bandwidths during nonlinear interactions within the climate–cryosphere system, we recombine (i.e. sum) their respective bispectral zones of Fig. 8. To do so, we include those zones that contain either one, two, or three frequency components of a particular bandwidth (Figs. 3 and 10). The general pattern of an energy cascade is robust, but it becomes clearer that precession cycles do not contribute much to the fuelling of eccentricity-paced climate variability directly (Fig. 10a–c). The 40 kyr periodicity during the Middle and Late Pleistocene gains energy from precession-dominated interactions but loses energy in obliquity- and eccentricity-dominated interactions (Fig. 10a–c). A comparison of conservativity of the recombined zones that include at least one precession, obliquity, or eccentricity component indicates that approximately similar amounts of energy are exchanged in (triad) interactions involving obliquity as in those involving eccentricity (Fig. 10d).

4.1.1 Age uncertainty and signal-to-noise ratio

By combining more than 50 individual benthic foraminiferal δ¹⁸O records, the LR04 stack simultaneously decreases the age uncertainty and increases the ratio of signal-to-noise considerably in comparison to any record individually (Lisiecki and Raymo, 2005a). Reduced age uncertainties and improved signal-to-noise ratios of stacked records may be obtained through the computation of averaged sedimentation rates in combination with statistical “binning” (Lisiecki and Raymo, 2005a), a probabilistic approach using a hidden Markov model (Ahn et al., 2017) or, as is the case in this reanalysis of the (1 kyr resampled/binned) LR04 stack, by additional statistical “smoothing” (Chaudhuri and Marron, 1999). For successful application of advanced bispectral analysis techniques, good age control and (very) high signal-to-noise ratios are needed, because the distributions of only small amounts of signal over time and frequency are considered (Hasselmann et al., 1963; King, 1996). For consecutive computational steps that use bispectra (i.e. the bispectrum and the integration over the bispectrum), increasingly better constrained and more noise-free data are required, because each step zooms in further on these small amounts of signal that are associated with non-sinusoidal geometries. The bispectra of climate cycles presented here show clear interactions (Fig. 4), and the integrations of the bispectrum reveal qualitatively consistent patterns of energy transfers for analysis with varying window lengths (Figs. 5b and S3). Therefore, we argue that the LR04 stack has sufficiently accurate age calibrations and a high enough signal-to-noise ratio for the successful application of advanced bispectral analysis.

4.1.2 Resolvability in time and frequency domains

Similar to the power spectrum, the bispectrum is marked by a trade-off between resolution in the time versus frequency domains (Fig. S3). Bispectral results become more significant, and can yield greater degrees of freedom, if a larger number of similarly shaped cycles, or waveforms, can be included for analysis (i.e. a lengthening of the window). For example, studies on natural nearshore waves or those on flume-generated waves (periodicities of seconds to minutes) use hour-long stable time series with high sampling resolutions and wave numbers (de Bakker et al., 2014). Bispectral analyses of such data sets are well resolved in the frequency domain. Furthermore, they permit the computation of robust confidence levels on the results and the selection of statistical settings that yield high degrees of freedom (e.g. Elgar and Guza, 1985; de Bakker et al., 2015). However, for the Pliocene and Pleistocene record, similar statistical objectives are unattainable, because no two climate cycles are alike and because of the low number of ∼110 kyr cycles in the Pleistocene (King, 1996; Lisiecki, 2010). As a result of this absence of “ergodicity” (i.e. a stable waveform through space and time), which characterises the unique Pliocene–Pleistocene climatic evolution as captured by the LR04 stack, relatively short window lengths are preferred to obtain well-resolved results in the time domain. The Hann windowing function visually clarifies the interactions in both the bispectra as well as the integrations over the bispectrum. Similarly to previous bispectral studies on the palaeoclimate archive (Hagelberg et al., 1991; King, 1996), we cannot define confidence levels for the bispectral analysis of the LR04 stack, because the number of cycles per window is too small. However, the bispectral couplings among climate cycles documented here are robust for varying bispectral settings (Figs. 5, 8, and S3 to S5).

4.1.3 Quantifying cycle geometries of non-ergodic records

In contrast to the relatively more robust bispectral results (see previous section), absolute geometry values are more strongly affected by the non-ergodic nature of the Pliocene–Pleistocene record. Choices of (i) window length, (ii) method of detrending the time series, and (iii) the application of a windowing function or not can have profound effects on the robustness of the results (Tim E. van Peer, personal communication, 2018). These somewhat arbitrary choices affect both ways of quantifying geometries similarly (i.e. using the central moments or higher-order spectra; Figs. 2, S2, and S3) (Tim E. van Peer, personal communication, 2018). The sensitivity of geometry computations to data processing techniques was not previously appreciated in full and may have resulted in overoptimistic confidence levels for cycle geometry values computed on an Oligocene and early Miocene climate record (Liebrand et al., 2017). We note that the uncertainty in skewness and asymmetry becomes larger for cycles with values near the mean (=0). The relatively high values of about −1 or +1, similar to those computed on the LR04 stack of the Middle and Late Pleistocene, are much less sensitive to data processing methods and have been independently computed on at least three occasions (Fig. 2) (Hagelberg et al., 1991; King, 1996; Lisiecki and Raymo, 2007). Remarkably, excess kurtosis values between ∼3000 and ∼350 ka even change sign between computations that do apply a Hann window (resulting in leptokurtic cycles; Fig. 2) and those that do not (yielding platykurtic cycles that reach values of −1.0 during the Pleistocene; not shown). By comparison, skewness and asymmetry values appear much more robust to various methods of data processing than values of excess kurtosis. This greater sensitivity of excess kurtosis to windowing is not surprising considering that fourth-moment statistics (and the trispectrum) zoom in on even smaller amounts of signal in comparison to the third central moments and the bispectrum (Collis et al., 1998).

5.2.1 The climatic precession motor

First, we speculate that the fuelling of the ice ages by (mainly) the precession periodicities is largely linked to the low- to midlatitude monsoons (Rutherford and D'Hondt, 2000; Berger et al., 2006; Bosmans et al., 2015a). The (sub)tropical zones, where most insolation is received, constitute Earth's heat engine, and their monsoons function as both a source of moisture and as a meridional teleconnection to build up large polar (land-)ice volumes. The interactions between NH land-ice volumes and the African and Asian monsoons are well documented for large parts of the Pliocene and Pleistocene (deMenocal, 1995; Sun et al., 2006; Cheng et al., 2016). The monsoonal (i.e. atmospheric) redistribution of heat and moisture was most probably aided by orbital-scale variability in the strength of Atlantic (i.e. oceanic) meridional overturning circulation (AMOC) (Lisiecki et al., 2008; Ziegler et al., 2010; Elderfield et al., 2012). Based on the bispectral results, we infer that during the Pliocene and Early Pleistocene this predominantly monsoonally driven precession motor fuels the 40 kyr obliquity-paced ice-age cycles, aided by more linear climatic–cryospheric responses resulting from variability in insolation at this periodicity, either zonally/cross-equatorially integrated (Huybers and Tziperman, 2008; Bosmans et al., 2015b) or locally at higher latitudes. The precession motor is especially strongly expressed in B^Im(O, O, P) interactions (Zone 5, Fig. 8e). However, uncertainty remains over the exact moisture sources for the Pleistocene land-ice volume fluctuations (e.g. Werner et al., 2001; de Boer et al., 2012).

5.2.2 The cryospheric obliquity motor

Second, we hypothesise that the role reversal of obliquity-paced climate cycles (from net energy receiver to net energy provider) across the MPT, in both B^Im(E, E, O) interactions (Fig. 4) and total net energy transfers (Figs. 5b and 8b), is linked to a resonance of (auto)cycles of crustal sinking and a delayed rebound with ∼110 kyr long (allo)cycles of eccentricity-modulated precession. In this hypothesis, crustal sinking is triggered by the passing of a land-ice mass loading threshold, and the delayed response to both eccentricity-modulated precession and obliquity is caused by the inertia of bulky, continental-sized ice sheets, which also results in a nonlinear positive ice-albedo feedback and hence a hysteresis pathway for glacial–interglacial cycles (Calov and Ganopolski, 2005; Abe-Ouchi et al., 2013). This obliquity motor evolves during and after the MPT from a lower harmonic of two to three obliquity-paced cycles (Fig. 6b) (Ridgwell et al., 1999; Huybers and Wunsch, 2005) into a resonance (i.e. phase coupling/locking or synchronisation) with the ∼110 kyr eccentricity-paced cycles, mainly the ∼95 kyr component (Fig. 6a) (Hagelberg et al., 1991; Tziperman et al., 2006; Crucifix, 2012; Abe-Ouchi et al., 2013). We speculate that the resonance of crustal sinking and a delayed isostatic rebound of the lithosphere–asthenosphere also shifted the sensitivity of the climate–cryosphere system to insolation forcing at higher latitudes, where the relative influence of obliquity compared to precession is greater than at lower latitudes. We conjecture that this shift in sensitivity to obliquity forcing was especially enhanced during glacial terminations (Huybers and Wunsch, 2005; Huybers, 2011; Konijnendijk et al., 2015).