  There is evidence that ice age cycles are paced by astronomical forcing,
suggesting some kind of synchronisation phenomenon. Here, we identify the type
of such synchronisation and explore systematically its uniqueness and
robustness using a simple paleoclimate model akin to the van der Pol relaxation
oscillator and dynamical system theory. As the insolation is quite a complex
quasiperiodic signal involving different frequencies, the traditional concepts
used to define synchronisation to periodic forcing are no longer applicable.
Instead, we explore a different concept of generalised synchronisation in terms
of (coexisting) synchronised solutions for the forced system, their basins of
attraction and instabilities. We propose a clustering technique to compute the
number of synchronised solutions, each of which corresponds to a different
paleoclimate history. In this way, we uncover multistable synchronisation
(reminiscent of phase- or frequency-locking to individual periodic components
of astronomical forcing) at low forcing strength, and monostable or unique
synchronisation at stronger forcing. In the multistable regime, different
initial conditions may lead to different paleoclimate histories. To study their
robustness, we analyse Lyapunov exponents that quantify the rate of convergence
towards each synchronised solution (local stability), and basins of attraction
that indicate critical levels of external perturbations (global stability). We
find that even though synchronised solutions are stable on a long term, there
exist short episodes of desynchronisation where nearby climate trajectories
diverge temporarily (for about 50 kyr). (...)
