Bispectra of climate cycles show how ice ages are fuelled

Abstract. The increasingly nonlinear response of the climate–cryosphere
system to insolation forcing during the Pliocene and Pleistocene, as
recorded in benthic foraminiferal stable oxygen isotope ratios (δ18O), is marked by a distinct evolution in ice-age cycle frequency,
amplitude, phase, and geometry. To date, very few studies have thoroughly
investigated the non-sinusoidal shape of these climate cycles, leaving
precious information unused to further unravel the complex dynamics of the
Earth's system. Here, we present higher-order spectral analyses of the LR04
δ18O stack that describe coupling and energy exchanges among
astronomically paced climate cycles. These advanced bispectral computations
show how energy is passed from precession-paced to obliquity-paced climate
cycles during the Early Pleistocene (from ∼2500 to ∼750 ka) and ultimately to eccentricity-paced climate cycles during the
Middle and Late Pleistocene (from ∼750 ka onward). They also
show how energy is transferred among many periodicities that have no primary
astronomical origin. We hypothesise that the change of obliquity-paced
climate cycles during the mid-Pleistocene transition (from ∼1200 to ∼600 ka), from being a net sink into a net source of
energy, is indicative of the passing of a land-ice mass loading threshold in
the Northern Hemisphere (NH), after which cycles of crustal depression and
rebound started to resonate with the ∼110 kyr eccentricity
modulation of precession. However, precession-paced climate cycles remain
persistent energy providers throughout the Late Pliocene and Pleistocene,
which is supportive of a dominant and continuous fuelling of the NH ice ages
by insolation in the (sub)tropical zones, and the control it exerts on
meridional heat and moisture transport through atmospheric and oceanic
circulation.



5
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 (f1 and f2 respectively), and diagonal lines correspond to the sum frequency (f3). 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 10 frequency f3, but its value can be read off at the crossing points with the xand 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).

Quantifying geometries using the bispectrum
Adopting the bispectral method to extract cycle geometries (Fig. 2, Supp. Fig. 2, and Supp. Fig. 3) (Elgar, 1987), skewness 15 and asymmetry are computed from the biphase. The biphase reflects the ratio between the real and imaginary parts of the bispectrum, which (recall) correspond to skewed and asymmetric cycle shapes, respectively. The biphase is defined by Eq.
where n and l range from 1 to the Nyquist frequency N, with n > l and n + l ≤ (Elgar, 1987), and fp refers to the primary 20 frequency (a.k.a. natural or resonance frequency, or first harmonic).

Integration of the spectrum and bispectrum
To quantify power spectral density, we apply a standard integration of the power spectrum. To obtain conservative, net 25 energy transfers per frequency (defined by the nonlinear source term Snl), 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 with a coupling coefficient. For palaeoclimate time series, no coupling coefficient has been determined Clim. Past Discuss., https://doi.org/10.5194/cp-2019-43 Manuscript under review for journal Clim previously. Therefore, we make minimum assumptions and use a coupling coefficient that only corrects for a frequency of W(f1,f2) = (f1 + f2). This correction insures energy conservation during triad interactions and allows for qualitative interpretations. Furthermore, if a more appropriate coupling coefficient for palaeoclimatic purposes can be found, it could also scale the net energy transfers of the bispectrum (already corrected using Snl during integration) to absolute energy transfers that are equal to power spectral density (see Section 4.3). We express the integral of Snl in terms of an integration, 5 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): 10 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 15 partially overlapping windows of the time series (i.e., a moving window) (Fig. 5, Fig. 6, and Supp. Fig. S3). Integration can be performed including all frequencies (i.e., over the entire imaginary part of the bispectrum), or over specific zones only, in case subsets of interactions need to be further examined (Fig. 3, Table A1) (de Bakker et al., 2015;de Bakker et al., 2016).

Spectral and bispectral zones 20
To obtain more insight in 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 and sea-swell 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 periodicities of ~56.2, ~28.1 and 12.5 kyr (i.e., 25 frequencies of 17.8, 35.6 and 80.0 Myr −1 ). 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 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 (Fig. 7, Fig. 8, Fig. 9, 30 Supp. Fig. S4 and Supp. Fig. 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). Clim. Past Discuss., https://doi.org/10.5194/cp-2019-43 Manuscript under review for journal Clim. Past Discussion started: 24 April 2019 c Author(s) 2019. CC BY 4.0 License.

80↓, 40↑) and A
?@ (40↓, ∞↓, 40↑), which also indicate a gain for 40 kyr cycles by transfers from numerous other frequencies. Overall, the Late Pliocene to Early Pleistocene bispectrum depicts strong energy gains at the obliquity paced 40kyr cycle through many triad interactions, some of which are with other astronomically-paced climate cycles.
The second example is a bispectrum across the mid-Pleistocene transition (i.e., the "~80-kyr world"), and shows several 5 negative interactions where either p1 or p2 is equal to ~80 kyr (Fig. 4b) However, we observe an overall gain of energy at the 70-to 80-kyr periodicities, and a dual role for the 40-kyr cycles, which operate as both energy receiver and supplier simultaneously. We note that the ~80-kyr periodicity is not an obvious triad of the primary astronomical cycles (Table A2), apart from the junction at A ?@ (95↑, 405↑, 77↓), where no interactions are documented in this particular bispectrum. The 80-kyr cycle thus starts as a lower harmonic cycle of the 40-kyr cycle, which in turn is a first harmonic of obliquity forcing, as well as a "combination tone" of precession cycles, Table S2. This lower 15 harmonic origin for the 80-kyr cycle is supported by the strong negative interaction at A ?@ (80↑, 80↑, 40↓).

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 ( Fig. 5a, Fig. 6, and Fig. 7a) (Hays et al., 1976;Berger, 1977). Furthermore, it is characterized by variability in spectral power governed by the long term 10 (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 ( Fig. 5a and Fig. 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 (Fig. 7a). This is indicative of a reduced perturbation of the Earth's system that results from a smaller difference in insolation forcing between 15 the hemispheres.

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-20 paced and nonastronomically-paced climate cycles. They reveal strong increases in both energy gains (in red) and losses (in blue) toward the present, and from ~1,000 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 cycle from ~2,500 ka onward and are marked by a modest 1.2 Myr beat in amplitude ( Fig. 5b and Fig. 7b). A second location of strong energy gains is the ~80-kyr cycle from ~1,000 ka onward (Fig. 6b), which at ~700 ka shifts (Fig. 5b) to a dominant gain for the 95-kyr component (Fig.  25 6a). Further energy gains are located at a ~200-kyr cycle from ~750 ka onward, a 50-kyr to 60-kyr cycle between ~1500 and 500 ka, and (very weak) at the 24-kyr cycle from ~650 ka onward (Fig. 5b). Energy losses are also observed in the total integration over the bispectra (Fig. 5b, Fig. 6 and Fig. 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 ( Fig. 5b and Fig. 6). Despite this more dispersed pattern of energy losses, we highlight the more conspicuous decreases: the 40-kyr cycle from ~550 ka onward, the 30-to 36-kyr cycles between 1,500 and 1,000 ka, the 30-to 26-kyr cycles between ~2,500 and ~2,000 ka, the ~28-kyr cycle from ~800 ka onward, the 24-kyr cycle between ~1,000 and ~800 ka, the 22-kyr cycle from ~700 onward, the 19-kyr cycle from 2,000 ka onward, and the 15-kyr cycle intermittently from ~1,000 ka onward ( Fig. 5b).  (Lisiecki andRaymo, 2005, 2007). A relatively weak ~200-kyr cycles is present from 10 ~750 ka onward (Fig. 5c). No strong precession response is present (Raymo and Nisancioglu, 2003;Raymo et al., 2006). The 40-kyr cycle in d 18 O partially mimics the 1.2 Myr amplitude modulation of the obliquity cycle, especially from ~3,000 ka onward ( Fig. 5c) (Lourens and Hilgen, 1997;Lisiecki and Raymo, 2007). We draw attention to a subtle "arc" in spectral power, which begins with a strong 405-kyr cycle, at ~3.500 ka, that becomes gradually shorter in duration until it resonates with the 95-kyr eccentricity modulation of precession, at ~2,650 ka, and then culminates in several weak 70 to 80 kyr cycle 15 at ~2,000 ka ( Fig. 5c, see also Fig. 1b) (Lisiecki, 2010;Meyers and Hinnov, 2010;Viaggi, 2018).

Zooming in on energy transfers
To better understand the different roles of specific (bandwidths of) climate cycles in nonlinear interactions, we determine conservative net energy transfers per bispectral zone (Fig. 8, Fig. S4, and Fig. S5). These zones reflect the "layers" that make 20 up total energy exchanges as depicted in Figure 5b. We assess if triad interactions are conservative of energy for each of the bispectral zones and for the entire bispectrum ( Fig. 9). Lastly, we recombine the bispectral zones again for frequency bandwidths involving precession-, obliquity-or eccentricity-paced climate cycles to highlight their influences in nonlinear interactions ( Fig. 10).

Zonal energy transfers
In Zone 1, from 1,000 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 resampled) LR04 stack, by statistical "smoothing" (Chaudhuri and Marron, 1999). For successful application of advanced bispectral analyses techniques, good age control and (very) high signal to noise ratios are needed, because the distributions of only small amounts of variance 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 noise-freer data is required, because each step zooms in further on these small amounts of variance. 5 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 ( Fig. 5b and Supp. Fig. 1).
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.

Resolvability in time and frequency domains
Like the power spectrum, the bispectrum is marked by a trade-off between resolutions in the time versus frequency domain ( Fig. S1). 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 15 time series with high sampling resolutions and wave numbers (de Bakker et al., 2014). Bispectra analyses of such datasets 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). 20 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. 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 are too small. However, the bispectral couplings among climate cycles 25 documented here are reproducible for varying bispectral settings (Fig. 5, Fig. 8, and Supp. Fig. 3 to Supp. Fig. 5).

Quantifying cycle geometries of nonergodic records
In contrast to the (qualitatively) reproducible bispectral results (see previous section), absolute geometry values are more strongly affected by the nonergodic nature of the Pliocene-Pleistocene record. Choices of (i) window length, (ii) method of 30 detrending the time series, and (iii) the application of a windowing function or not, can have profound effects on reproducibility (van Peer, 2018). These, somewhat arbitrary choices affect both ways of quantifying geometries similarly (i.e., using the central moments or higher order spectra (Fig. 2, Supp. Fig. 2, and Supp. Fig. S3) (van Peer, 2018). The sensitivity of geometry computations to data processing techniques was not previously appreciated in full, and may have resulted in over-optimistic confidence levels for cycle geometry values computed on an Oligocene and early Miocene 5 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, like those present the LR04 stack of the Middle and Late Pleistocene, are much less sensitive to data processing methods, and have been independently reproduced on at least three occasions ( Fig. 2) (Hagelberg et al., 1991;King, 1996;Lisiecki and Raymo, 2007). Remarkably, kurtosis values between ~3,000 and ~350 ka even change sign between computations that do apply a Hann window (resulting in 10 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 kurtosis. This greater sensitivity of kurtosis to windowing is not surprising considering that fourth moment statistics (and the trispectrum) zoom in on even smaller amounts of variance in comparison to the third central moments and the bispectrum (Collis et al., 1998). 15

Coupling coefficient
For ocean waves, the relationship between the integration of the bispectrum and absolute nonlinear energy transfers is determined by a coupling coefficient, which is based on the second-order Boussinesq theory (Herbers and Burton, 1997;Herbers et al., 2000). This Boussinesq approximation describes how energy is transferred from the primary ocean waves to 20 higher and lower frequency components by near-resonant nonlinear triad interactions. A coupling coefficient is needed (i) to ensure the conservation of energy within a single triad interaction (i.e., all energy lost at a certain frequency is compensated by an energy gain at the other two frequencies, and vice versa), and (ii) to obtain absolute energy transfers, which are directly comparable with the changes in the power spectrum when the waves propagate toward the coast (Herbers et al., 2000;de Bakker et al., 2014). Here, we apply a similar approach to palaeoclimate time series, and demonstrate that energy 25 conservation can be ensured using a simple coupling coefficient that multiplies the energy transfers at each specific sum frequency after integration, with that specific frequency (f3) (see Section 2.4.1). Comparison of energy gains and losses, for both the integration over the entire imaginary part of the bispectrum, as well as for each bispectral zone separately, shows that sum ( KL \ , see Eq. 7) and difference ( KL ] , see Eq. 8) interactions balance each other, which is indicative of energy conservation, and qualitatively trustworthy signals in energy transfers (Fig. 9). However, the second step toward absolute 30 energy transfers is not developed here (see qualitative comparison of spectral and bispectral results in Fig. 6), because a more fundamental understanding is needed of the many climatologic and oceanographic, biologic, sedimentologic and lithologic processes that affect the (globally integrated) benthic foraminiferal d 18 O record.

Revisiting the problem of the ice ages 5
The common denominator of (i) the 40-kyr problem of the Pliocene and Early Pleistocene, (ii) the ~110-kyr problem of the Middle and Late Pleistocene, and (iii) the (nonlinear) origins of many climate cycles without primary astronomical pacemaker (e.g., the ~28-kyr and 80-kyr cycles), is the mismatch in the distribution of power spectral density between benthic foraminiferal d 18 O and records and insolation (Raymo and Nisancioglu, 2003;Lisiecki, 2010;Lourens et al., 2010).
Furthermore, (iv) during the mid-Pleistocene transition, this discrepancy increases significantly (Lisiecki and Raymo, 2005;10 Clark et al., 2006;Huybers, 2007), by the shifting of spectral power of climate records from the 40-kyr, via the ~56-kyr and the ~80-kyr, to the ~110-kyr periodicity (Rial and Anaclerio, 2000;Lourens et al., 2010;Viaggi, 2018). In light of the bispectral results obtained in this study, we revisit the problem of the ice ages (Agassiz, 1840;Milankovitch, 1941;Hays et al., 1976), in search for the most parsimonious explanation(s) for the origins and evolution of climate cycles with and without astronomical periodicities. The answer that we distil from the bispectra of Pliocene and Pleistocene climate cycles, 15 in greater clarity than before (Hagelberg et al., 1991;Rial and Anaclerio, 2000;Lourens et al., 2010), is that these problems (i.e., i to iv) are intrinsically linked by the many, increasingly nonlinear responses of Earth's climate-cryosphere system to insolation forcing. Namely, the bispectra confirm (Hays et al., 1976;Hagelberg et al., 1991), that energy is transferred from the lower to the higher periodicities from at least 2,500 ka onward. More importantly, they show (i) how this happens, i.e., through which nonlinear triad interactions (Fig. 4), and (ii) how their involvement in transferring energy evolves trough time 20 (Fig. 5, Fig. 8, and Fig. 10), especially during MPT.

The origins of climate cycles: toward a better mechanistic understanding
To advance the understanding of the origins of Pliocene and Pleistocene climate cycles, we explore the potential climatic and cryospheric processes that cause the documented nonlinear triad interactions and energy transfers. Similar to the difficulties 25 in linking spectral peaks of benthic foraminiferal d 18 O records, which describe the integrated effect of deep water temperatures and cryosphere (Hays et al., 1976;Elderfield et al., 2012;Rohling et al., 2014;Crowhurst et al., 2018), to particular mechanisms, it is not straightforward to identify (singular) causes for bispectral peaks (Fig. 4). Should we attempt to link every individual triad interaction to a specific climatic or cryospheric process? Or, alternatively, should we search for

The cryospheric obliquity motor
Second, we hypothesize that the role reversal of obliquity paced climate cycles (from net energy receiver to net energy provider) across the MPT, in both ?@ (E, E, O)-interactions (Fig. 4) and total net energy transfers ( Fig. 5b and Fig 8b), is linked to a resonance of (auto-) cycles of crustal sinking and a delayed rebound with ~110-kyr (allo-) cycles of eccentricity modulated precession. In this hypothesis, crustal sinking is caused by the passing of a land-ice mass-loading threshold, and 5 the delayed response, is caused by the inertia of bulky, continental-sized ice sheets, which causes a hysteresis-pathway for glacial-interglacial cycles (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 with the ~110-kyr eccentricity paced cycles, mainly the ~95-kyr component (Fig. 6a) (Hagelberg et al., 1991;Tziperman et al., 2006;Abe-Ouchi et al., 2013). We speculate that the resonance of crustal sinking and a delayed isostatic 10 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).

Climatic and tectonic boundary conditions
If our interpretations of the bispectra are indeed supportive of an astronomically-paced monsoonal/AMOC "push" of (heat and) moisture, followed by a resonant lithospheric-asthenospheric "pull" through a cooling of the NH polar region after an initial phase of ice growth, for sufficiently large glacial land-ice masses to form, then they would highlight the importance of dynamic topography and the shape (i.e., extent, volume, and especially height) of the NH ice sheets in determining the 20 pattern and geometry of climate-cryosphere variability observed for the Middle and Late Pleistocene (Roe, 2006;Bintanja and Van de Wal, 2008;Creveling et al., 2015). However, these mechanisms do not explain why large-scale NH glaciations only developed during the Late Pliocene and Pleistocene, and not during earlier time-periods with comparable astronomical configurations. The Pliocene-Pleistocene climatic-cryospheric evolution can therefore not be seen outside a context of longterm, multi-Myr changes in climatic and tectonic boundary conditions. The "arc" in spectral power (~3,500-~2,000 ka) (Fig.  25 5c), which is largely absent in bispectral power (Fig. 5b), and which precedes the strong glaciations of the MPT and Late Pleistocene, may be indicative of these changing boundary conditions. We speculate that this arc reflects a change in the response time of the cryosphere, potentially through glacial landscaping by earlier (smaller) glacial cycles (Clark and Pollard, 1998;Tabor and Poulsen, 2016) or through changing tectonic/geographic boundary conditions and associated oceanic circulation patterns (Haug and Tiedemann, 1998;Kender et al., 2018). These developments could have 30 preconditioned the Earth's system for the full-blown bipolar glacial cycles of the Middle and Late Pleistocene. These cycles were amplified by variability in radiative forcing on astronomical time scales and set against a background of an all-Cenozoic-time low "# $ %&' (Shackleton, 2000;Beerling and Royer, 2011;Martinez-Boti et al., 2015;Chalk et al., 2017).

Outlook
We present considerable advances in the understanding of the origins of many Pliocene and Pleistocene climate cycles by 5 applying advanced bispectral analysis to the palaeoclimate archive. However, several lines of research remain to be explored further. Here we list those that we think need further attention in the future: • For nearshore ocean waves, energy transfers are fully described by the imaginary part of the bispectrum (Norheim et al., 1998). Further supporting evidence for the importance of the imaginary part in describing most/or all energy exchanges comes from spectral analysis of synthetically skewed and asymmetric time series (not shown). These 10 spectra indicate that the higher harmonic peaks of skewed cycles are about an order of magnitude smaller than those of asymmetric cycles, especially for higher harmonic peaks (3f, 4f, …) than the first higher harmonic (2f). However, for palaeoclimate studies, a more comprehensive understanding of the real part of the bispectrum may be still be important for obtaining a(n) (even) more complete description of nonlinear energy exchanges (e.g., see the ~28-kyr cycle), even though their contribution is probably very small (Supp. Fig. 1 and Supp. Fig. 2). 15 • Although we obtain conservative energy exchanges with our simple coupling coefficient, other studies are needed to scale these exchanges to the power spectrum (Fig. 6). Coupled climate-cryosphere models or energy balance models may be used to derive or approximate coupling coefficients for palaeoclimates that are based on "first principles" (i.e., a physicochemical understanding of the Earth system). Alternative coupling coefficients may also change how the bispectrum redistributes energy from climate cycles with low to those with high periodicities. 20 However, we note that energy conservation should be respected for every triad interaction in the bispectrum.
• The proposed relationships between nonlinear triad interactions as documented in bispectra of climate cycles and specific processes of the climate-cryosphere system (i.e., monsoons, isostatic rebound), remain speculative. More data (analysis) and/or (isotope-enabled) modelling studies are needed to support these hypotheses.
• There is a limited amount of information that can be extracted from studying a single climate record. Bispectral 25 analysis of other (Pliocene-Pleistocene) records may shed more light on the dynamics of Earth's palaeoclimate system. Not only in elucidating the couplings among astronomical cycles, but also across the spectral continuum, and in relation to 1/f "red" noise and/or the power-law scaling of energy exchanges between frequencies (see Section 2.4.1. and Section 4.2) (Hasselmann, 1976;Bak et al., 1987Bak et al., , 1988Bak, 1996;Pelletier, 1998;Huybers and Curry, 2006). Examples of interactions that could be studied in this way are those between astronomical and sub-30 astronomical (e.g., millennial) climate cycles (Hagelberg et al., 1994;Wara et al., 2000;Da Silva et al., 2018), and among cycles with centennial, decadal and annual durations (King Hagelberg and Cole, 1995). We emphasize that Clim. Past Discuss., https://doi.org/10.5194/cp-2019- other applications of bispectral analysis techniques to the palaeoclimate archive should be limited to records with good age control and high signal-to-noise ratios (see Section 4.1.1). Therefore, we argue that for these purposes data stacks are preferred.

Conclusions 5
In this interdisciplinary study, we present a new, higher-order spectral analysis and interpretation of Pliocene and Pleistocene climate cycles as present in the well-established LR04 globally-averaged benthic-foraminiferal d 18 O stack. These advanced bispectral analysis techniques were initially developed by the nearshore ocean wave community, and are adapted here to make them more applicable for paleoclimatic research purposes. They show, more thoroughly than before, the complexities of the Pliocene-Pleistocene climatic-cryospheric evolution, namely that (i) a nonlinear energy exchanges resulted in the 10 distinct asymmetric cycle geometries of the Middle and Late Pleistocene, (ii) energy is transferred across the power spectrum, among astronomical and nonastronomical climate cycles, (iii) precession-paced climate cycles fuelled both the 40kyr and 110-kyr ice age cycles of the Pliocene and Pleistocene, (iv) the role of obliquity-paced climate cycles changes from being a net sink into a net source of energy during the MPT (i.e., in ?@ (E, E, O)-interactions), and (v) after the MPT these obliquity-paced climate cycles helped fuel the ~110-kyr eccentricity-paced climate cycles of the Middle and Late 15 Pleistocene.
We postulate two distinct fuelling mechanisms for the ice ages to explain the energy cascade of negative interactions from climate cycles with high frequencies to those with low frequencies: (i) a continuous Pliocene-Pleistocene climatic "precession motor" and (ii) a Middle-to-Late Pleistocene cryospheric "obliquity motor". We interpret the dominant 20 precession fuelling of obliquity paced climate cycles (i.e., the precession motor) as evidence of a low latitude, (sub-) tropically-driven climate system, in which the monsoons (and oceanic meridional overturning circulation) transport heat and moisture to higher latitudes. We argue that the role-reversal of obliquity-paced climate cycles during the MPT is indicative of the passing of a land-ice mass-loading threshold, after which (auto-) cycles of crustal sinking and (delayed) rebound start to resonate, i.e., become frequency-and phase-coupled, with the ~110-kyr (allo-) cycles of eccentricity-modulated 25 precession. This resonance may have caused a shift of sensitivity to (obliquity-paced) insolation changes at increasingly higher latitudes, especially during terminations.
Despite the fact that the evidence for energy transfers agrees well with previously published mechanisms, their unprecedentedly detailed description here is an important step forward toward a more comprehensive solution of the problem 30 of the ice ages, because they largely reconcile the mismatch in power spectral density between records that capture the combined effects of deep-sea temperatures and land-ice volumes and those of insolation. Furthermore, we can now state with greater certainty than before, that the geometry (i.e., asymmetry) of the Pliocene and Pleistocene ice ages is in very close agreement with the classical Milankovitch Theory of astronomical climate forcing and its derivatives/superlatives that give greater weight to Earth's internal nonlinear responses. If the bispectral energy transfers documented here indeed track the entire redistribution of power over the spectrum, then the 40-kyr problem of the Pliocene and Early Pleistocene, and the ~110-kyr problem of the Middle and Late Pleistocene are reduced to finding appropriate coupling coefficients that describe 5 the relevant climatic and cryospheric mechanisms-a promising outlook.

Acknowledgments
We thank X. Bertin for inviting D. L. to present this work at an early stage at the Université de La Rochelle. We are grateful to T. E. van Peer for making us aware of the difficulties in reproducing geometries of nonergodic records. We thank L. E.   Absolute amplitude modulations and filters of the ~110-kyr, 40 kyr, and ~21-kyr cycles. (d) Phase evolution of the ~110-kyr cycles with respect to eccentricity (Laskar et al., 2011a;Laskar et al., 2011b), computed over a 750 kyr sliding window. 5 Blackman-Tukey phase estimates are only shown when coherent (Paillard et al., 1996). We focus the phase interpretation on the time interval with high amplitude ~110-kyr cycles (transparent area). Phase relationship of d 18 O to obliquity is not shown because it is assumed during the astronomical tuning process. Phase to precession is not shown because globally averaged   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 5 grey shaded areas represent suborbital periodicities. Difference frequencies f1 and f2 are plotted along the along the horizontal and vertical axes, respectively. Sum frequencies f3 (i.e., f1 + f2) 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.   September 21) at 65°N (Paillard et al., 1996;Laskar et al., 2004). (b) Conservative net energy transfers during the Pliocene and Pleistocene, computed by integrating over the entire imaginary part of the bispectrum (see Methods). Input data is the resampled LR04 stack (Lisiecki and Raymo, 2005). (c) Power spectral density of the LR04 stack (Lisiecki and Raymo, 5 2005). 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 −1 , and 2 degrees of freedom.      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 f1 and f2 (dashed lines) are computed using Eq. 8, those of f3 (single coloured lines) following Eq. 7, and their difference (black lines) by Eq. 6. 5  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 5 zonal conservativities with eccentricity in orange, obliquity in purple, and precession in blue.  Table A1. Bispectral zones defined in this study. These zones correspond to  Table A2. Bispectral triads among climate cycles with astronomical frequencies. One, two, or three astronomical frequencies can partake in triad interactions, and we refer to these coordinates in the bispectrum as single, double or triple junctions, respectively. These triads reflect positions in the frequency-frequency domain where nonlinear energy transfers among astronomically paced climate cycles are likely to occur. They correspond to the crossing points of the coloured lines 5 in Fig. 3. In case of a "single junction", an astronomically paced cycle exchanges energy with itself, or with a frequency that is only slightly offset. Single junctions where f1 = f2 constitute the first higher harmonic interactions. Double or triple junctions are "combination tones" of two or three different astronomically-paced cycles (Hagelberg et al., 1991;King, 1996;Rial and Anaclerio, 2000).