the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
The onset of chaos in nonautonomous dissipative dynamical systems: a loworder oceanmodel case study
Stefano Pierini
Mickaël D. Chekroun
Michael Ghil
A fourdimensional nonlinear spectral ocean model is used to study the transition to chaos induced by periodic forcing in systems that are nonchaotic in the autonomous limit. The analysis relies on the construction of the system's pullback attractors (PBAs) through ensemble simulations, based on a large number of initial states in the remote past. A preliminary analysis of the autonomous system is carried out by investigating its bifurcation diagram, as well as by calculating a metric that measures the mean distance between two initially nearby trajectories, along with the system's entropy. We find that nonchaotic attractors can still exhibit sensitive dependence on initial data over some time interval; this apparent paradox is resolved by noting that the dependence only concerns the phase of the periodic trajectories, and that it disappears once the latter have converged onto the attractor. The periodically forced system, analyzed by the same methods, yields periodic or chaotic PBAs depending on the periodic forcing's amplitude ε. A new diagnostic method – based on the crosscorrelation between two initially nearby trajectories – is proposed to characterize the transition between the two types of behavior. Transition to chaos is found to occur abruptly at a critical value ε_{c} and begins with the intermittent emergence of periodic oscillations with distinct phases. The same diagnostic method is finally shown to be a useful tool for autonomous and aperiodically forced systems as well.
Understanding the mechanisms that lead to the onset of chaos in dissipative dynamical systems is of fundamental importance both from a cognitive viewpoint and for the correct use of the mathematical models on which the systems are based. Chaos arises in such systems as a control parameter in the governing equations crosses a given threshold. A huge amount of work has been devoted to analyzing the transition to chaos in the framework of autonomous dynamical systems, i.e., in systems in which the external forcing and the coefficients do not depend on time. The various routes to chaos in autonomous dissipative systems – in the presence of timeindependent forcing — include perioddoubling cascades, intermittency and crisis, quasiperiodic routes, and global bifurcations (e.g., Nicolis, 1995; Hilborn, 2000; Ott, 2002; Tél and Gruiz, 2006; Strogatz, 2015).
Nonautonomous dissipative dynamical systems represent a crucial extension of autonomous systems for practical applications, since the external forcing in most real systems – whether deterministic, random or both – depends, typically, on time. Despite their importance, nonautonomous systems have received, until recently, less attention than autonomous systems. Transition to chaos induced by timedependent forcing has, nonetheless, been studied in several significant cases. A classical example is the Van der Pol oscillator (Van der Pol, 1920, 1926), in which chaotic relaxation oscillations emerge under the effect of an external periodic forcing. A few more recent examples in the climate sciences include (i) transition to chaos due to quasiperiodic forcing (Le Treut and Ghil, 1983; Ghil, 1994); (ii) modification of the autonomous transition by periodic forcing (Sushama et al., 2007); and (iii) the important contributions of Anna Trevisan and colleagues to the data assimilation problem for chaotic systems, in which the data stream can be seen essentially as a timedependent forcing (Trevisan and Uboldi, 2004; Carrassi et al., 2008a, b; Trevisan and Palatella, 2011).
The onset of chaos is analyzed here in the framework of nonautonomous systems, which has been received rapidly increasing attention recently in the context of climate dynamics (Ghil et al., 2008; Chekroun et al., 2011; Bódai and Tél, 2012; Bódai et al., 2013; Pierini, 2014; Drótos et al., 2015, 2017; Ghil, 2015, 2017; Pierini et al., 2016; Lucarini et al., 2017). Our study focuses on a fourdimensional nonlinear spectral ocean model (Pierini, 2011), which is subjected to periodic forcing, chosen as the simplest form of time dependence. Cases will be considered that are nonchaotic in the autonomous limit, so that the chaos that emerges in the system is strictly associated with the nonstationarity of the forcing.
The study makes use of ensemble simulations performed with many initial states distributed in a given subset of phase space, following the methodology of Pierini (2014) and Pierini et al. (2016). The overall idea is that the relevant information in the climate system must be derived from statistical analyses of an ensemble of different system trajectories, each corresponding to a different initial state, provided that the corresponding trajectories have converged to the system's timedependent attractor.
Such an attractor is called a pullback attractor (PBA; e.g., Ghil et al., 2008; Chekroun et al., 2011; Kloeden and Rasmussen, 2011; Carvalho et al., 2012) in the mathematical literature and a snapshot attractor (e.g., Romeiras et al., 1990; Bódai and Tél, 2012; Bódai et al., 2013) in the physical literature; it provides the natural extension to nonautonomous dissipative dynamical systems of the classical concept of an attractor that is fixed in time for autonomous systems. A global PBA is defined as a timedependent set 𝒜(t) in the system's phase space that is invariant under its governing equations, along with the equally timedependent, invariant measure μ(t) supported on this set, and to which all trajectories starting in the remote past converge (Arnold, 1998; Rasmussen, 2007; Kloeden and Rasmussen, 2011; Carvalho et al., 2012). In the deterministic case, it is understood that 𝒜(t) depends also on the particular forcing, say F(t), that is being applied, but this dependence is usually not kept track of in the notation. In the random case, the PBA is called a random attractor, and the dependence on the specific realization ω of the noise process is often included in the notation, as 𝒜(t,ω).
Pierini et al. (2016) rigourously proved that a weakly dissipative nonlinear model like the one used there and herein does possess a global PBA, subject to mild integrability conditions on the forcing. In the present study, the numerical approach used for the systematic investigation of the system's PBAs follows Pierini (2014) and Pierini et al. (2016). Further diagnostic tools will be introduced for the present periodicforcing setup, and a new diagnostic tool will also be proposed to monitor the onset of chaos in our nonautonomous system.
The paper is organized as follows. In Sect. 2, the mathematical model is described. In Sect. 3, the main properties of the autonomous system are summarized and an apparent paradox related to the sensitivity to initial states in the periodic regime is discussed. In Sect. 4, the results obtained for the periodically forced system are presented and discussed; the new crosscorrelationbased method specifically formulated to characterize the onset of chaos is introduced and applied to the specific case at hand. This method helps characterize the transition to chaos as the amplitude of the periodic forcing increases, as well as document the coexistence of local PBAs with chaotic and nonchaotic behavior within the model's global PBA. In Sect. 5, the same method is shown to be a useful tool also for autonomous and aperiodically forced systems. Finally, in Sect. 6 the results are summarized and conclusions are drawn. An appendix illustrates in greater detail the coexistence of local PBAs that are chaotic and nonchaotic in the setting of a periodically forced Van der Pol–Duffing oscillator.
The highly idealized model of the oceans' winddriven, doublegyre circulation (Ghil, 2017, and references therein) used in the present study is governed by the system of four nonlinear, coupled ordinary differential equations derived by Pierini (2011); Vannitsem (2014), Vannitsem and De Cruz (2014) and De Cruz et al. (2018) used this model to represent the ocean component in their loworder climate models. The author introduced such a loworder model to complement the process studies on the Kuroshio Extension's lowfrequency variability previously carried out with a much more detailed, primitive equation ocean model (e.g., Pierini, 2006; Pierini et al., 2009; Pierini and Dijkstra, 2009). The same loworder model was later used by Pierini (2014) and Pierini et al. (2016) to explore the PBAs of the system in various cases. Here we merely review the main aspects of the model; for all the technical details and parameter values, the interested reader should kindly refer to Pierini (2011).
The dynamics are governed by the evolution equation of potential vorticity in the quasigeostrophic approximation on the beta plane for a shallow layer of fluid, superimposed on an infinitely deep quiescent lower layer. Pierini (2006) found such a reducedgravity model to be a good approximation for process studies of the Kuroshio Extension's lowfrequency variability. The flow is described by the streamfunction ψ(x,t): like in the previous studies, ψ, the horizontal coordinates $\mathbf{x}=(x,y)$ and the time t are dimensionless, but the dimensional time will be plotted in all the time series presented in this study to emphasize the typical timescales of the oceanic phenomena under investigation.
A fourdimensional spectral model is obtained by expanding the streamfunction in a rectangular domain as follows:
The orthonormal basis i〉 is defined as follows:
where α is a real positive constant. This basis satisfies the freeslip boundary conditions along the borders of the rectangular domain; it also captures the oceanic flow's westward intensification thanks to the exponential factor first introduced in the twodimensional model of Jiang et al. (1995).
The four nonlinear coupled ordinary differential equations that govern the evolution of the vector $\mathbf{\Psi}\left(t\right)=\left({\mathrm{\Psi}}_{\mathrm{1}},{\mathrm{\Psi}}_{\mathrm{2}},{\mathrm{\Psi}}_{\mathrm{3}},{\mathrm{\Psi}}_{\mathrm{4}}\right)$ can be written (see Pierini, 2011) as
The coefficients of the nonlinear and linear terms in the equation are encapsulated by the rank3 and rank2 tensors J and L, respectively, and the forcing is represented by the vector w (see Pierini, 2011). The forcing w is obtained from a suitable doublegyre surface wind stress curl, while G is defined in the present paper to be periodic,
with period ${T}_{\mathrm{p}}=\mathrm{2}\mathit{\pi}/\mathit{\omega}$, while γ and ε are positive dimensionless parameters.
To construct the system's PBAs, ensembles of forward time integrations are carried out; each of these starts at t=0 from a different initial point contained in a given subset Ω of the model's fourdimensional phase space and ends at $t={T}_{\ast}=\mathrm{400}$ years; as shown in Figs. 3 and 9 below, T_{∗} is much greater than the spinup time in all cases. Following Pierini (2014) and Pierini et al. (2016), the fourdimensional hypercube Ω is defined as follows:
and the initial data are all chosen to satisfy Ψ_{1}=Ψ_{2} and Ψ_{3}=Ψ_{4}, i.e., to lie within a plane set embedded in Ω.
The ensembles will consist of 15 000 initial data at t=0 that are regularly spaced either in Ω or in a small subset thereof. For the sake of graphical representation, maps of various quantities will be plotted in the rectangle $\mathrm{\Gamma}\equiv \mathit{\{}\mathrm{70}\le {\mathrm{\Psi}}_{\mathrm{1}}\le \mathrm{150},\mathrm{150}\le {\mathrm{\Psi}}_{\mathrm{3}}\le \mathrm{120}\mathit{\}}\subset \mathrm{\Omega}$ that lies in the (Ψ_{1},Ψ_{3}) plane. In the discussion of the results, we will refer, for the sake of simplicity and concision, to the model's trajectories as being defined in the (Ψ_{1},Ψ_{3}) plane but, naturally, the actual trajectories evolve in the full fourdimensional phase space.
3.1 The autonomous model's attractors
We begin by analyzing some basic properties of the autonomous system that will be useful in the subsequent investigation. The bifurcation diagram of Fig. 1a shows the range of variability of Ψ_{1} vs. the forcing parameter γ. The value γ=1 corresponds to a global bifurcation that manifests itself by a sudden transition from a smallamplitude limit cycle to a relaxation oscillation with a much higher amplitude. The previous results in this respect (Pierini, 2011; Pierini et al., 2016) will be further bolstered by those in Sect. 5 herein (Figs. 14 and 15), which are based on the diagnostic tool proposed in Sect. 4.2.
Figure 1b shows the limit cycle in Γ arising from arbitrary initial data for γ=1.1, which corresponds to the red vertical line in the bifurcation diagram of panel (a). For γ=1.35 the attractor is chaotic (green line in panel a; see Sect. 5.1 further below). In this case, the map of the suitably normalized decimal logarithm of the probability density function (PDF) of the trajectories in Γ is plotted in Fig. 1c; it is defined by ${q}_{k}\left(t\right)={\mathrm{log}}_{\mathrm{10}}\left[\mathrm{1000}\times {n}_{k}\left(t\right)\right]$. Here n_{k} is the number of trajectories contained at time t in the kth cell belonging to the same regular grid of N square cells of width ΔΨ that is used in our ensemble simulations, with N=15 000 and ΔΨ=2, and it is plotted at $t={T}_{\ast}=\mathrm{400}$ years, i.e., after spinup.
In an autonomous dynamical system, the attractors do, by definition, not depend on time, i.e., an attractor is a geometric object in phase space that is fixed in time. However, any attractor that is not a fixed point – whether a limit cycle, torus or strange attractor – can contain timedependent trajectories. Such ensembles of trajectories arising from specific sets of initial states will be plotted to illustrate the attractors of the autonomous system studied herein.
Following Pierini et al. (2016), in Fig. 2a, b the attractors that correspond to the two cases in Fig. 1b, c are represented by the time evolution of ${P}_{{\mathrm{\Psi}}_{\mathrm{3}}}\left(t\right)={\mathrm{log}}_{\mathrm{10}}(\mathrm{1000}\times {p}_{{\mathrm{\Psi}}_{\mathrm{3}}})$, where ${p}_{{\mathrm{\Psi}}_{\mathrm{3}}}\left(t\right)$ is the PDF of localization of the Ψ_{3} variable; see Pierini et al. (2016) for technical details. The dense distribution of ${P}_{{\mathrm{\Psi}}_{\mathrm{3}}}\left(t\right)$ for γ=1.35, as seen in Fig. 2b, is clearly associated with the chaotic character of the flow, while the periodic distribution of ${P}_{{\mathrm{\Psi}}_{\mathrm{3}}}\left(t\right)$ that corresponds to γ=1.1 in Fig. 2a is due to the different phases that each trajectory attains on the limit cycle, depending on the initial point.
In Fig. 2c, d, the same attractors are characterized through the metric σ that was introduced by Pierini et al. (2016); this metric measures the mean divergence of trajectories over the total integration time T_{∗} and is defined as follows. The instantaneous Euclidean distance between two initially close trajectories is δ(t), and its normalized value is given by ${\mathit{\delta}}_{\mathrm{n}}\left(t\right)=\mathit{\delta}\left(t\right)/\mathit{\delta}\left(\mathrm{0}\right)$. Then σ is simply the average of δ_{n} over T_{∗},
Pierini et al. (2016) found the quantity σ to be a good indicator of the degree of sensitivity of the system's evolution with respect to the initial state during the phase of convergence to the attractor.
The determination of the PBAs of the periodically forced system and the application of the new qualitative and quantitative diagnostic methods proposed in Sect. 4 need an analysis of the behavior of trajectories that lie at t=0 on a given subset Ω of phase space, as is the case when calculating σ(X,Y) above. Thus, investigating the behavior of model trajectories as they emerge from Ω is the most unifying and distinctive feature of the present model study.
The map of σ in Fig. 2c reveals, in the autonomous case at hand, the same striking features found by Pierini et al. (2016) for the nonautonomous, aperiodicforcing case, namely the coexistence of extended regions of Γ with σ⩽1, shown by cold colors, and with σ>1, appearing as warm colors. In the first case, two trajectories that are initially close remain close at all times, as seen in Fig. 3a. In the second case, though, two trajectories that are initially close may attain a large phase difference once they have converged to the attractor (cf. Fig. 3b), while still remaining perfectly coherent.
In the chaotic case with γ=1.35, the warmcolor regions, in which σ>1, overwhelm the coldcolor regions, in which $\mathit{\sigma}\le \mathrm{1},$ cf. Fig. 2d. To illustrate the two types of behavior, Fig. 3c, d show the evolution of Ψ_{3}(t) of two initially nearby trajectories. If σ<1, as is the case near P_{1}, the two trajectories are virtually coincident (Fig. 3c). If, on the contrary, σ>1, as is the case near P_{2}, the two aperiodic signals lose their coherence (Fig. 3d). Finally, it is worth noting that, for simulations with sufficiently small γ (not shown), σ<1 everywhere.
A different and useful way of looking at these two types of behavior is to analyze the corresponding mixing properties of the flow in the model's phase space. To do so, one can make use of the system's entropy (Shannon, 1948):
Here Γ is decomposed into a regular grid of N square cells of width ΔΨ (with N=15 000 and ΔΨ=2, so that the grid corresponds to that of the initial data used in our ensemble simulations) and p_{k}(t) is the probability of localization in the kth cell at time t of the trajectories emanating at time t=0 from a given subset ϑ⊂Γ.
Figure 4a shows the intersection with the (Ψ_{1},Ψ_{3}) plane at t=300 years of 15 000 trajectories originating from the box ϑ_{1} that coincides with the ΔΨ×ΔΨ grid cell centered at P_{1}; the red dots correspond to the case γ=1.1 and the green dots to γ=1.35. Figure 4b shows ${S}_{{\mathit{\vartheta}}_{\mathrm{1}}}$ for the two cases; note that ${S}_{{\mathit{\vartheta}}_{\mathrm{1}}}\left(\mathrm{0}\right)=\mathrm{0}$ since all the initial states lie in the single cell ϑ_{1}, and thus p_{1}=1. The entropy of the periodic case γ=1.1, characterized by σ<1, oscillates between 0 and 1, with the final evolution limited to virtually a single cell over the limit cycle; the latter cell is enclosed in the red circle of Fig. 4a.
In the chaotic case γ=1.35, σ≤1 for 43 % of the points contained in ϑ_{1}, while σ>1 for the remaining points. The evolution of the former leads to the localized blue dots in Fig. 4a while the evolution of the latter leads to the green dots scattered over the strange attractor. The green line of Fig. 4b, giving ${S}_{{\mathit{\vartheta}}_{\mathrm{1}}}$ computed with all the trajectories, shows the gradual spreading of the initial points with σ>1.
Figure 4c, d show the same quantities for the initial ΔΨ×ΔΨ box ϑ_{2} centered at P_{2}. The chaotic case is similar to that for ϑ_{1}, but with a greater entropy; however, the periodic case differs in that now σ>1 – cf. Fig. 2c. Figure 4c shows that the asymptotic evolution of the very small ϑ_{2} covers a limited but significant part of the limit cycle, as seen by comparing this figure with Fig. 1b; the corresponding entropy in Fig. 4d eventually oscillates periodically between the values ${S}_{{\mathit{\vartheta}}_{\mathrm{2}}}\sim \mathrm{2.2}$–3.7. Figure 4 thus demonstrates clearly the usefulness of the metric σ in characterizing subsets of Γ and the effect of the control parameter γ.
Finally, it is worth stressing that, since the forcing is constant, the range of variability of the entropy in the chaotic case with σ>1 must tend to zero as the number of points tends to infinity. This tendency is clearly illustrated by the green line of Fig. 4d. However, the range of variability of ${S}_{{\mathit{\vartheta}}_{\mathrm{1}}}$ in the chaotic case (Fig. 4b) is still quite large after 400 years because, as pointed out above, the number of points with σ>1 contained in ϑ_{1} is relatively small.
3.2 An apparent paradox
We conclude the analysis of the autonomous system by discussing an apparent paradox. We have just seen that, in regions of Γ where σ>1, the trajectories for γ=1.1 exhibit sensitive phase dependence on initial data, as shown, for instance, by Fig. 3b, by the red dots in Fig. 4c and by the red curve in Fig. 4d. Sensitive dependence on initial data is usually associated with chaotic dynamics, but in this case the dynamics are periodic.
This paradox is resolved by noting that such sensitivity concerns only the phase of the periodic trajectories, as already noticed in the previous subsection and, in addition, it occurs only if the initial data lie outside the attractor, e.g, elsewhere on Γ; on the attractor, this phase sensitivity disappears, as we will show below. On the contrary, in the chaotic case γ=1.35, the sensitivity to initial data for trajectories with σ>1 always holds, off the attractor as well as on it. This is in excellent agreement with the chaotic character of the dynamics in the latter case.
We must show, therefore, that the trajectories are stable for γ=1.1 and unstable for γ=1.35, once they have settled onto the attractor. This distinction between the two cases can already be inferred from Figs. 3 and 4 but it is worth investigating the issue in greater detail. The usual quantitative approach relies on the computation of the leading finitetime Lyapunov exponent λ of each trajectory (e.g., Ott, 2002).
The results (not shown) are consistent with the assumption above, but the exponents are highly dependent on the time T_{λ} over which the finitetime exponents are computed, and on the amplitude of the perturbation superimposed on the reference trajectory at each time step T_{λ}. Moreover, the assumption of exponential divergence of chaotic trajectories is not fully met in our highly nonlinear framework, so that the transition between periodic and chaotic dynamics may actually occur at a value of λ that is not exactly equal to 0. A qualitative diagnostic method is instead illustrated in Fig. 5, and furthermore we propose an alternative quantitative method in Sect. 4.
Let us consider the points lying in the black rectangle shown in Fig. 4c at t=300 years: the corresponding evolution at four subsequent time instants, T_{Δ}=25 years, is shown in Fig. 5a for γ=1.1 (red dots). Note that the period of the orbits on the attractor is ${T}_{\ast}p=\mathrm{14.08}$ years, i.e., ${T}_{\mathrm{\Delta}}>{T}_{\ast}p$.
The stability of the trajectories under consideration is clearly demonstrated in Fig. 5a by the compact form and limited extent of the cluster: indeed, these points that start from t=300 years evolve anticlockwise around the attractor, covering it roughly 6 times during the interval 4T_{Δ}=100 years that separates the first snapshot from the last one. On the contrary, Fig. 5b shows that for γ=1.35 (green dots) the compact form of the initial cluster is lost already after a single 25year lapse of time. The trajectories will thus soon be scattered over the strange attractor, due to their divergence.
In conclusion, our autonomous system becomes chaotic for sufficiently large values of γ, e.g., for γ=1.35. The system's periodic regime spans a range of γ that includes the bifurcation at γ=1, which is apparent in Fig. 1a.
Moreover, we have shown that, when γ=1.1, regions of Γ exist within which the mean normalized distance σ between two initially nearby trajectories is larger than unity; see again Fig. 2c. In this case, despite the attractor's being a limit cycle, the trajectories leaving from such regions of phase space experience sensitive phase dependence on the initial data. Although sensitive dependence is typically associated with chaotic systems, this sensitivity is not in contradiction with the periodic character of the solutions: as a matter of fact, the trajectories under discussion are stable and the sensitive dependence disappears once the trajectories have converged onto the attractor.
The existence of regions with σ>1 for an autonomous periodic system is an important feature for the transition to chaos when the system is subjected to timedependent forcing: this issue will be discussed in the next section. Besides, the new diagnostic method introduced in Sect. 4.2 to help analyze transition to chaos in the nonautonomous case will be applied in Sect. 5.1 to the autonomous case.
For the idealized doublegyre model governed by Eq. (2), we have seen that the autonomous system given by ε=0 in Eq. (3) exhibits a limit cycle when γ=1.1; this limit cycle corresponds to the largeamplitude relaxation oscillation of Fig. 1b. However, Pierini (2014) showed that, when the same model, with the same value of γ, is subjected to periodic forcing with ε=0.2 and T_{p}=30 years, it exhibits chaotic, cyclostationary and cycloergodic behavior; see Figs. 2 and 3 therein and the related discussion. To understand the transition to deterministically chaotic behavior induced by the forcing, we will now apply in Sect 4.1 the methodology used in Sect. 3 to the attractors corresponding to γ=1.1 and T_{p}=30 years across the intervening parameter range $\mathrm{0}\le \mathit{\epsilon}\le \mathrm{0.2}$. In Sect. 4.2, the transition to chaos induced by the periodic forcing will be analyzed in greater detail through an additional method here developed explicitly for this purpose.
We know from the rigorous proof of Pierini et al. (2016, Appendix A) that our idealized ocean model possesses a global PBA, in the general case of timedependent forcing, whether periodic or aperiodic. PBAs are, in fact, timedependent mathematical objects that characterize the asymptotic behavior of a nonautonomous dissipative dynamical system (e.g., Chekroun et al., 2018, Fig. 2).
It is common, though, in the literature of periodically forced dynamical systems to study the asymptotic behavior of such a system by an iterated stroboscopic map $\mathit{\left\{}u\right(t)\to u(t+T\left)\mathit{\right\}}$, where u is the variable and T is the period. In the particular case in which the system's driver is periodic, so is the PBA; see Sect. 2.3.2 of Chekroun et al. (2018) for a rigorous proof. In contradistinction, a “normal” – i.e., forward rather than pullback – attractor visualized in the embedded space $\left(u\right(t),u(t+T\left)\right)$ and built by using points along a long trajectory $\mathit{\left\{}u\right(t):{t}_{\mathrm{0}}\le t\le {t}_{\mathrm{0}}+\mathcal{T}\mathit{\}}$ is static and therefore contains less information than the corresponding PBA, no matter how long the interval 𝒯 may be. Furthermore, PBAs built from ensembles of initial data allow us to visualize in one single picture the coexistence of different types of dynamical behavior in terms of disjoint PBAs; see, for instance, Fig. A2 in Appendix A herein. Besides, the PBA framework is useful for the visualization of fractal structures that arise when noise is superimposed to the periodic forcing. A stroboscopic map analysis may not easily reveal such fractal features; see Chekroun et al. (2011), as well as Sect. 3.4 of Chekroun et al. (2018).
4.1 The pullback attractors of the forced system
For all the above reasons, we now present the PBAs of our periodically forced ocean model. As already just mentioned, the PBAs of a periodically forced dissipative system are always periodic, but the system can be either chaotic or nonchaotic, depending on its parameter values. For the sake of simplicity, we will refer below to the PBAs of a chaotic and nonchaotic system, abbreviated as CPBAs and NPBAs, respectively.
For the autonomous case, the time evolution of P_{Ψ3}(t) and the map of σ(X,Y) – already shown in Fig. 2a and c, respectively – are again plotted in Fig. 6a and e for the sake of comparison. Figure 6d, h correspond to the reference CPBA studied by Pierini (2014). Note also that one of the two reference cases studied by Pierini et al. (2016) has the same values of γ and ε, but the latter parameter was multiplied by the aperiodic forcing shown in Fig. 17 below. Two intermediate cases that correspond to ϵ=0.05 and ϵ=0.1 are shown in Fig. 6b, f and c, g, respectively. In Sect. 4.2, we will show in greater detail that the transition to chaotic behavior occurs abruptly, when crossing a critical value ϵ_{c} that lies between the two intermediate values of 0.05 and 0.1; here we merely provide some qualitative arguments showing that, in fact, the case ϵ=0.05 is still periodic, while the case ϵ=0.1 is chaotic.
Before proceeding with the analysis of the results in Fig. 6, we recall that chaotic systems subjected to periodic forcing can be studied either by ensembles of trajectories – as done herein – or by stroboscopic averages using a single long trajectory, provided the assumption of cycloergodicity holds (e.g., Boyles and Gardner, 1983): the latter result extends the classical ergodicity property valid for strange attractors of autonomous systems (e.g., Eckmann and Ruelle, 1985) to chaotic, periodically forced systems. An example of this equivalence is shown in Fig. 3c, d of Pierini (2014) for the present model and for the parameter values corresponding to the P_{Ψ3}(t) and σ(X,Y) plotted in Fig. 6d, h herein (obviously, the trajectory used in that example was derived from a region with σ>1). However, the existence of regions with σ⩽1, as well as with σ>1, in the two chaotic cases (Fig. 6g, h) shows that the cycloergodicity assumption fails to hold for our idealized ocean model. Our system must, therefore, be investigated using the ensemble approach, which we pursue throughout this paper. In fact, had we only used the stroboscopic map method, we would never have discovered the existence of two types of local attractors, namely CPBAs and NPBAs, in our model.
Figure 7 provides further information on the four cases illustrated in Fig. 6. In Fig. 7a–d, the magenta dots represent the intersection with the (Ψ_{1},Ψ_{3}) plane at t=400 years of 15 000 trajectories, whose initial points (in blue) are evenly distributed in Γ at t=0; in addition, in Fig. 7e–h, the corresponding entropy S_{ϑ=Γ} is plotted as a function of time, along with the number n_{o} of cells that are occupied by at least one point. Clearly, the structure of the PBA snapshot in Fig. 7b, for ε=0.05, is very similar to that of the autonomous case in Figs. 1b and 7a, while the PBA snapshots plotted at t=400 years – for ε=0.1 and 0.2 in Fig. 7c and d, respectively – are quite different.
To understand this difference better, we focused in Fig. 8 on the subdomain ϑ_{2} of Γ that was defined in Fig. 4c and for which σ>1. The autonomous case has already been analyzed in Sect. 3.1: in fact, Fig. 8a, e are equivalent to Fig. 4c, d. In the case of ε=0.05, the same behavior is found, i.e., the sensitivity to initial data leads to only a compact subset of the attractor being covered; this implies the periodicity of the trajectories.
On the contrary, for ε=0.1 the intersection of the trajectories at t=400 years with the (Ψ_{1},Ψ_{3}) plane, shown by the magenta dots in Fig. 8c, is virtually indistinguishable from the one that appears in Fig. 7c, when the initial data are selected in the whole of Γ: this excellent match is an unequivocal sign of the mixing property of chaotic dynamics, as already discussed for the autonomous case γ=1.35 in connection with Figs. 4 and 5. That the same property holds for ε=0.2 in Fig. 8d is not surprising, since Pierini (2014) already recognized the model's chaotic behavior for this parameter value.
Finally, it is instructive to visualize Ψ_{3}(t) for a couple of trajectories that are very close at t=0, as done in Fig. 3 for the autonomous case. Figure 9 shows Ψ_{3}(t) for the four cases of Figs. 6–8 and for trajectories that emerge from ϑ_{2}. In Fig. 9a, b the two trajectories are periodic, but with a phase difference. In the two chaotic cases of ε=0.1, 0.2 in Fig. 9c, d, both trajectories are clearly aperiodic. In Fig. 9c, though, i.e., in the case that is closer to the transition, this aperiodicity is merely associated with a temporary shift in phase of an otherwise periodic signal within the intervals t≃40–120 years and t≃360–400 years.
The intermittent behavior seen in Fig. 9c appears – from many simulations that are not shown here – to be typical of chaotic solutions near the transition point ε_{c} and suggests a possible mechanism through which chaos is induced by an external periodic forcing. For values of ε just past ε_{c}, the model still tends to behave periodically, but the external forcing is sufficiently strong to entrain a trajectory occasionally into a nearby region, where the periodicity is preserved but the phase differs by a finite amount. Since these shifts are very sensitive to the initial data, the result is a chaotic trajectory characterized by separate intervals of periodic oscillations with a different phase. This mechanism also explains why the transition to chaos leads to a notable increase in the measure of the regions in Γ where sensitive dependence to initial data occurs; such an increase is visually obvious when comparing Fig. 6e, f with Fig. 6g, h. As ϵ increases further, the duration of the intervals of constant phase decreases, and the oscillations tend to become more genuinely aperiodic, as seen in Fig. 9d.
This behavior is similar to the intermittency found in autonomous dissipative systems (e.g., Manneville and Pomeau, 1979; Pomeau and Manneville, 1980), in which a trajectory switches back and forth from periodic to aperiodic oscillations provided a certain control parameter of the system – e.g., the amplitude of the steady, timeindependent forcing – crosses a given threshold. In our nonautonomous system, the amplitude of the periodic forcing ε plays a similar role. This transition to chaos induced by timedependent forcing appears, therefore, to be directly linked to the existence of regions in phase space in which sensitive dependence to initial data occurs in the limit of periodic solutions. Thus, the chaotic behavior merely due to the timedependent nature of the forcing can be traced back to the apparently paradoxical property of the autonomous system that was emphasized in Sect. 3.2. This striking observation deserves to be analyzed in greater depth in future studies.
4.2 Transition to chaos studied by a crosscorrelation method
Recognizing qualitatively whether a PBA is chaotic or not is relatively simple; e.g., this can be done through the heuristic arguments illustrated in Figs. 4 and 5 and through those outlined in the previous subsection and illustrated in Figs. 6–9. But how does one characterize the transition from periodic to chaotic dynamics as a control parameter, such as the amplitude ε of the periodic forcing, changes?
We have already discussed in Sect. 3.1 the limitations of using the mean finitetime Lyapunov exponents. Here we propose a new, simple and robust method that is particularly useful in our periodicforcing case, but can be applied also to any autonomous system and even to aperiodically forced systems; the latter situations will be addressed in Sect. 5. In Sects. 3 and 4.1, we have relied on the mixing properties of chaotic dynamics, as measured by the system's entropy S_{ϑ}, to recognize the occurrence of chaotic behavior. Now we rely on the emergence of aperiodic signals from a subset of Γ; this subset will necessarily be contained in the region where sensitive dependence on initial data occurs, i.e., where σ>1.
The most obvious approach would be to compute the power spectrum of each trajectory. Periodic signals can, however, be quite complex, as seen, for instance, in Fig. 9a, b; this complexity makes it quite difficult to identify a parameter whose value will distinguish, accurately and reliably, between periodic and chaotic dynamics, based solely on the Fourier spectra of a finite number of finitelength trajectories.
We propose a simpler alternative method that takes advantage of the ensemble simulations carried out to obtain the PBAs numerically. Let ${\stackrel{\mathrm{\u203e}}{\mathrm{\Psi}}}_{\mathrm{3}}$ and ζ_{3} be the mean and rootmean square values of ${\mathrm{\Psi}}_{\mathrm{3}}(X,Y,t)$, and consider the centered and normalized anomaly time series ${\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}}_{\mathrm{3}}(X,Y,t)$ and ${\stackrel{\mathrm{\u0303}}{\mathrm{\Psi}}}_{\mathrm{3}}\left({X}^{\prime},{Y}^{\prime},t\right)$, of Ψ_{3}, with
where (X,Y) and $\left({X}^{\prime},{Y}^{\prime}\right)$ are two points in Γ that are near to each other, and from which these two time series emerge at t=0. We can then compute the crosscorrelation between the two signals, after removing the initial transient, as
here ${T}_{\ast}=\mathrm{400}$ years is again the maximum integration time, and $T\le \mathit{\tau}\le T$, while T=50 years; once more, the following results are independent of T, provided it is sufficiently larger than the typical timescale of the phenomenon. Note also that, in the above definition, we have dropped the dependence of c on $({X}^{\prime},{Y}^{\prime})$ for the sake of conciseness.
Now, if σ<1, the two signals are periodic and virtually coincident, as seen, for instance, in Fig. 3a, c. Hence, defining the maximal crosscorrelation by
one will have $\mathrm{\Theta}(X,Y)\cong \mathrm{1}$, with Θ being attained at τ=0.
However, if σ>1, there are two possibilities:

either the PBA is not chaotic, in which case all couples (X,Y) and $({X}^{\prime},{Y}^{\prime})$ yield two periodic and virtually equal signals, apart from a finite phase difference, as seen, for instance, in Figs. 3b and 9a, b; in this case, again, Θ≅1, which will now occur at some lag τ≠0 that depends on the phase difference;

or the PBA is chaotic, in which case all couples yield two aperiodic and significantly different signals, as seen, for instance, in Figs. 3d and 9c, d; in this case, Θ will be substantially less than unity.
This alternative is illustrated in Fig. 10 for the four couples of trajectories plotted in Fig. 9a–d, all of which were initialized in ϑ_{2}, where σ>1.
It is then useful to analyze the maps of the parameter Θ. Figure 11a–d show Θ for the four cases of Figs. 6–9. In the two cases that we have already identified as nonchaotic, Θ varies within a range of values [0.95,1.1] that lies very close to unity, as expected; see Fig. 11a, b. There is only a small neighborhood of ${P}_{\mathrm{0}}=(\mathrm{16},\mathrm{83})$ in which Θ is very small: this is because P_{0} in the autonomous case is a fixed point; see Fig. 12.
The two cases that we have identified as chaotic appear here as Fig. 11c, d, and in them Θ exhibits in fact smaller values. These values lie in the range [0.2,0.8] for a large subset of the domain, where σ>1, as shown in Fig. 6g, h. Regions in which σ>1 but Θ≃1 are present as well, but the corresponding trajectories are nonetheless unstable once they have converged onto the PBA, because they will always pass sufficiently near trajectories that are chaotic, thanks to the mixing properties of the latter.
In summary, if Θ∼1 everywhere (yellow colors) we have an NPBA, whereas if Θ yields values that are sufficiently smaller than unity (grey colors) then we have a CPBA.
To summarize in a clear and simple way the information provided by the values of both σ and Θ, we introduce the integervalued parameter Φ, defined as follows:
Here Θ_{0} is a threshold value and Φ(X,Y) is plotted in Fig. 11e–h for Θ_{0}=0.8.
As an example of the usefulness of Φ, let us note that the σ maps in Fig. 6f and g are fairly similar, except for the more extended warmcolor region, where σ>1, in the second map. Recall, however, that the meaning of σ>1 is profoundly different if the system is chaotic, in which case mixing is present, as opposed to when it is not, in which case sensitive dependence to initial data concerns only the phase of the signal and is not accompanied by mixing. This ambiguity is resolved by the use of the step function Φ: if ϵ=0.05 (Fig. 11f), sensitive dependence to initial data, i.e., σ>1, yields the value of Φ=2 (green regions), since Θ≃1, which tells us that the system is not chaotic. On the contrary, if ϵ=0.1 (Fig. 11g), regions with Φ=2 (in red) appear within the green regions: this implies low Θ values and therefore chaos.
We conclude by analyzing the transition from NPBAs to CPBAs via a suitable function of the control parameter ε: this metric is provided by the average value 〈Θ〉_{Γ} of Θ over Γ. The graph of 〈Θ〉_{Γ}(ε) in Fig. 13 is obtained by performing many ensemble simulations of system trajectories with many distinct values of ε; the latter values are chosen to lie closer to each other, where the variation in 〈Θ〉_{Γ}(ε) is stronger. An abrupt transition from NPBAs, with 〈Θ〉_{Γ}≃1, to CPBAs occurs at ε_{c}≅0.09. Many additional analyses (not shown) for values just below and above ε_{c} confirm that this is in fact the critical value beyond which chaos sets in.
5.1 Application to the autonomous system
The diagnostic method proposed in Sect. 4.2 to monitor the transition from NPBAs to CPBAs in periodically forced systems relies on two properties: (i) in an NPBA all trajectories are periodic, and (ii) in a CPBA diverging aperiodic trajectories emerge from a subset of Γ, in which necessarily σ>1. Thus, the same crosscorrelationbased method can obviously be applied to an autonomous system as well. The method's application to the autonomous model studied in Sect. 3 will shed new light on the periodic vs. chaotic character of its solutions.
The graph in Fig. 14 shows 〈Θ〉_{Γ}(γ) and is obtained, like that of Fig. 13, by performing many ensemble simulations, each with a different value of γ, rather than ε, which equals zero in the present case. The first thing to notice is the sudden drop of 〈Θ〉_{Γ} at $\mathit{\gamma}={\mathit{\gamma}}_{\mathrm{c}}=\mathrm{1}$, where a global bifurcation separates smallamplitude limit cycles from largeamplitude relaxation oscillations, as shown in Pierini (2011) and in Sect. 3.1 here. In addition, Pierini (2012) investigated the stochastic version of this deterministic tipping point in the case of random forcing.
The chaotic nature of the attractor for γ_{c} is illustrated in Fig. 15. For γ=0.99, the Ψ_{3} time series exhibits the typical smallamplitude, purely periodic behavior studied in Sect. 3.1, while for γ=1 both small and largeamplitude oscillations occur irregularly in the same time series. The behavior at γ=1.01 illustrates the return to more regular behavior.
Thus, chaotic dynamics occurring in an extremely restricted γ range separates two different types of limit cycles. Figure 16 shows this dramatic transition in terms of σ: the chaotic nature of the flow for γ_{c}=1 is such that the warmcolored regions in which σ>1 overwhelm the coldcolored regions, as in Fig. 2d, where $\mathit{\gamma}=\mathrm{1.35}.$
For γ>1, the system is not chaotic – except for limited γ intervals centered at γ≃1.25 and γ≃1.335 – until a new abrupt drop of 〈Θ〉_{Γ} at $\mathit{\gamma}={\mathit{\gamma}}_{\mathrm{0}}=\mathrm{1.3475}$, shown by a dashed black line in Fig. 14. This drop signals the presence of chaotic attractors beyond γ_{0}; in particular, the chaotic case γ=1.35, shown by the solid green line and discussed in Sect. 3, lies just after this transition. It is worth noting that large fluctuations dominate the chaotic regime.
5.2 Application to an aperiodically forced system
Our crosscorrelation diagnostics have been shown to apply to both periodically forced and autonomous systems. Its validity, however, is even more general, since it extends to a large class of aperiodically forced systems as well. We choose the model setup of Pierini et al. (2016) to illustrate the latter possibility. A thorough analysis of this application is beyond the scope of the present study: we will therefore limit ourselves to analyzing the basic aspects of the problem and leave the details for a future investigation.
Pierini et al. (2016) considered the same system – governed by Eq. (2), and within the same parameter regime adopted here and in Pierini (2011). The forcing, though, was aperiodic and given by the following:
where ${\mathit{\epsilon}}^{\prime}>\mathrm{0}$ is a dimensionless coefficient and f(t) is a normalized, fixed realization of an Ornstein–Uhlenbeck process that has been smoothed to resemble multiannual windstress forcing of the midlatitude oceans' doublegyre circulation. Figure 17 shows G(t) for γ=1 and ${\mathit{\epsilon}}^{\prime}=\mathrm{0.2}$.
Figure 18 shows the evolution of two initially nearby trajectories emerging from P_{2}, along with the corresponding crosscorrelation, for γ=1.1, with ${\mathit{\epsilon}}^{\prime}=\mathrm{0.05}$ in panels (a–b) and ${\mathit{\epsilon}}^{\prime}=\mathrm{0.20}$ in panels (c–d); both cases have σ>1 and the corresponding time series of ${P}_{{\mathrm{\Psi}}_{\mathrm{3}}}$ are plotted in Fig. 4h, j of Pierini et al. (2016). The case ${\mathit{\epsilon}}^{\prime}=\mathrm{0.2}$ corresponds to the CPBA analyzed in detail by Pierini et al. (2016). The chaotic character of the solution is clearly visible from Fig. 18c; the crosscorrelation between the two signals is plotted in Fig. 18d and it is accordingly small.
On the contrary, the case ${\mathit{\epsilon}}^{\prime}=\mathrm{0.05}$ corresponds to an NPBA: the two signals in Fig. 18a develop a large phase difference after the initial transient, but are virtually identical and remain coherent at all times. Now, unlike in Figs. 3b and 9a, b, the two signals are not periodic because they are modulated by the aperiodic forcing, but this is irrelevant; in fact, the nonchaotic character of the solution can still be highlighted by the corresponding crosscorrelation in Fig. 18b, whose maximum value is Θ≃1, like in the autonomous and periodically forced case. Obviously, this is possible because the period of the modulated relaxation oscillation is much smaller than the timescale of the forcing; this is in fact the only condition required for the applicability of this diagnostic method to aperiodically forced systems.
We can therefore conclude that the parameter 〈Θ〉_{Γ} can be a valuable tool for monitoring the onset of chaos in aperiodically forced systems as well. For example, this diagnostic method can be applied to study the onset of chaos in systems that possess a drift mimicking global warming and other climate change scenarios (as done, for instance, in Drótos et al., 2015).
In this paper, we studied the transition from nonchaotic to chaotic PBAs in a nonautonomous system whose autonomous limit is nonchaotic, and in which, therefore, chaos is induced by the periodic forcing. The illustrative example chosen for this general problem was a loworder quasigeostrophic model of the midlatitude winddriven ocean circulation, subject to periodic forcing. The model was described and connected with previous work in Sect. 2.
We first investigated, in Sect. 3, the autonomous system, following up on the work of Pierini (2011), who obtained its governing equations and analyzed their solutions. Here, ensemble simulations based on a large number of initial data and the calculation of the system's entropy allowed us to determine novel and interesting features of the system subject to steady forcing.
To do so, we used the metric σ that was introduced by Pierini et al. (2016) and measures the time average of the distance between two trajectories that are very close at t=0 on a subset Γ of phase space. The analysis based on this metric yielded regions in Γ with σ values that can be either larger or less than unity. The spatial structure of these regions (cf. Fig. 2c, d here) is very similar to that of the nonautonomous case investigated by Pierini et al. (2016), as seen in Fig. 6 therein. This similarity suggests that the nonautonomous behavior of a dynamical system is profoundly influenced by the convergence properties of trajectories initialized off the attractor in the autonomous case; this finding, in turn, implies that ensemble simulations are very helpful in studying said properties.
Next we investigated, still in the autonomous case, the apparent paradox of regions with σ>1 coexisting in the periodic regime of γ=1.1 with the expected regions of σ<1. A large number of trajectories emanating from the small square box of Fig. 4c, for which σ>1, evolves into the extended red line shown in the same figure at a given time, t=300 years: this line belongs to the periodic attractor shown in Fig. 1b and the entropy evolution along it oscillates periodically (cf. Fig. 4d); but it is clearly distinct from the chaotic attractor apparent as the green cloud of points in the same figure panel. Further evidence for the stability of the trajectories in this periodic case with σ>1 is provided by Fig. 5a.
We conclude that sensitive phase dependence on initial data in a periodic regime may be present if the trajectories are initialized off the attractor, but that it disappears once the trajectories have converged onto the attractor. Clearly, generic sensitive dependence on initial data is, therefore, a necessary but not sufficient condition for chaotic behavior.
In Sect. 4, we studied the onset of chaos in the system subject to the periodic forcing given by Eq. (3) in the case of γ=1.1, in which the system is nonchaotic in the autonomous limit, so that chaos is induced by the periodic forcing. The PBAs were analyzed at first for four different sinusoidalforcing amplitudes, with $\mathit{\u03f5}=\mathrm{0},\mathrm{0.05},\mathrm{0.10}$ and 0.2, while using the same γ=1.1 and T_{p}=30 years.
We found that the first two cases, namely ϵ=0 and 0.05, are nonchaotic, while the other two, namely ϵ=0.10 and 0.20, are chaotic. A large number of trajectories emanating from the small square box where σ>1 in Fig. 8a–d evolves – depending on the value of ε – into two very different fixedinstant subsets of the model's PBA, with the snapshot taken at t=400 years, i.e., after convergence of the trajectories to the PBA. In the first two cases, this snapshot is a curvedline segment that belongs to the PBA, while in the latter two cases it covers the whole PBA, due to the typical mixing property of chaos.
An analysis of the trajectories, as shown in Fig. 9, for instance, indicates that the transition to chaos occurs via an intermittent emergence of periodic oscillations with different phases; see again Fig. 9c. We have shown that, for values of ϵ in the chaotic regime just above the transition, the periodic character of the system is still predominant, but the external forcing is now sufficiently strong to cause a trajectory's occasionally shifting to a phase space region in which a different phase prevails: since the shifts are very sensitive to the initial data, the result is a chaotic trajectory characterized by irregular jumps of the oscillatory solutions between distinct phases.
In Sect. 4.2, we introduced a novel diagnostic method for the study of the transition between nonchaotic and chaotic behavior as the amplitude ϵ of the periodic forcing increases. The method's basic idea is that in a nonchaotic regime any couple of initially nearby trajectories emerging from Γ remain coherent at all times, while in a chaotic regime aperiodic diverging trajectories emerge from a subset of Γ. Hence, a systematic recognition of the character of all the trajectories allows one to diagnose which of these two types of behavior occurs or whether the two actually coexist.
A simple and robust way to do this is to compute the crosscorrelation $c(X,Y,{X}^{\prime},{Y}^{\prime},\mathit{\tau})$ at lag τ between two initially nearby trajectories started at (X,Y) and $({X}^{\prime},{Y}^{\prime})$ in Γ, then compute its maximum value Θ(X,Y) over τ. If $\mathrm{\Theta}(X,Y)\simeq \mathrm{1}$ everywhere in Γ, then the system is nonchaotic and is therefore periodic under periodic forcing; on the contrary, if Θ is appreciably smaller than unity in some subset of Γ, the system is chaotic. The diagram of the average 〈Θ〉_{Γ}(ϵ) of Θ(X,Y) over Γ, as plotted in Fig. 13, reveals an abrupt transition to chaos at ϵ_{c}≃0.09.
This crosscorrelationbased method has also been applied to the autonomous system, in which the conditions required for the periodically forced case do apply as well. The diagram of 〈Θ〉_{Γ}(γ) in Fig. 14 reveals that chaotic dynamics occur at first within an extremely restricted range centered at the global bifurcation point $\mathit{\gamma}={\mathit{\gamma}}_{\mathrm{c}}=\mathrm{1}$, which separates smallamplitude, fairly smooth oscillations below γ_{c} from largeamplitude relaxation oscillations above it. A considerably broader range of chaotic behavior occurs in the autonomous case for values of γ greater than a threshold γ_{0}=1.3475.
We have then applied the crosscorrelationbased method to the aperiodically forced system studied by Pierini et al. (2016). Our results show that, in fact, this method can be applied to systems subject to aperiodic forcing when the system's intrinsic periodicity and the characteristic timescale of the external forcing are sufficiently well separated from each other, which is the case in Pierini et al. (2016); cf. Fig. 17 here. Once more, the crosscorrelations in Fig. 18c, d agree remarkably well with the character of the model trajectories in Fig. 18a, b.
Finally, the coexistence of local PBAs with chaotic vs. nonchaotic behavior within a global PBA – as first described by Pierini et al. (2016) in the aperiodicforcing case – was confirmed here for the periodically forced case; cf. Figs. 6 and 11. This situation was explored in greater depth in Appendix A for an even simpler, weakly dissipative nonlinear model, namely a Van der Pol–Duffing oscillator (e.g., Jackson, 1991), and additional references were given for this type of PBA bistability.
Overall, this paper provides additional insights into the complex and varied behavior that arises even in highly idealized atmospheric, oceanic and climate models from the interaction of nonlinear intrinsic dynamics with various types of external forcing. In addition, it stresses the importance of using the framework of nonautonomous dynamical systems and of their PBAs for a deeper understanding of this complexity and variety.
No data sets were used in this article.
The purpose of this Appendix is to provide further insight into the coexistence of local PBAs with quite different stability properties, as illustrated in the main text by Figs. 6 and 11. In complex autonomous systems, several local attractors may coexist for a given set of the system's parameters; each of these attractors possesses attracting sets of initial data that are typically separated by fractal boundaries (Grebogi et al., 1987). Basins of attraction with fractal boundaries have consequences for predictability: uncertainties in the initial state x_{0} may result in different types of dynamical behavior, depending on which basin x_{0} lies in; see, for instance, Roques and Chekroun (2011). When the number of the coexisting attractors is two, one speaks of bistability, and multistability refers colloquially to more than two coexisting attractors.
Our goal here is to illustrate how multistability of nonautonomous systems manifests itself unambiguously through the existence of disjoint local PBAs. In the case of periodically forced systems, such as those considered in this article, similar results can be inferred, of course, from the analysis of Poincaré maps. Nevertheless, the presence of other frequencies in the internal dynamics, the external forcing or the noise may render the analysis of Poincaré maps difficult, whereas the framework of PBAs naturally includes such additional levels of complexity (Chekroun et al., 2018). Pierini et al. (2016), for instance, already demonstrated the usefulness of this framework in the context of bistability.
Furthermore, the purpose of this Appendix is also to illustrate that multistability of local PBAs arises not only for our quasigeostrophic model, as discussed in the paper's main text. More generally, PBA coexistence occurs fairly often for externally forced systems, although a careful analysis of the flow's dependence on initial data may be required in practice in order to conclude on multistability.
A paradigm of multistability is provided by dissipative nonlinear systems that become Hamiltonian in the limit of vanishing dissipation, as is the case, for instance, in celestial mechanics (e.g., Ghil and Wolansky, 1992). In this situation, it is expected that the number of coexisting attractors exceeds any fixed bound in approaching this limit, as documented for various nearly integrable maps and flows; cf. Feudel et al. (1996), Zaslavsky and Edelman (2008), Celletti (2009), Pisarchik and Feudel (2014) and Dudkowski et al. (2016), and references therein.
To illustrate this multistability phenomenon, we consider the following periodically forced Van der Pol–Duffing oscillator, given by the secondorder nonlinear (ordinary differential equation) ODE:
where μ, b, F and ω are positive constants that determine the dynamical behavior of the system. This system is nonautonomous and its PBAs are analyzed hereafter in the $(x,\dot{x})$ plane. This nonlinear ODE arises in various applications such as in engineering, electronics, biology and neurology (Jackson, 1991; Kozlov et al., 1999; Kuznetsov et al., 2009). It combines the nonlinearity of the dissipation $\mathit{\mu}{\dot{x}}^{\mathrm{2}}x$, which characterizes the Van der Pol (1920) oscillator with that of the internal force −bx^{3}, which characterizes the Duffing (1918) oscillator.
Multistability was already numerically documented for Eq. (A1) by relying on Poincaré maps (Venkatesan and Lakshmanan, 1997; Dudkowski et al., 2016). In our calculations, we have followed Dudkowski et al. (2016) and assumed μ=0.2, $F=b=\mathrm{1.0}$ and ω=0.955. While this parameter regime does not correspond to the limit of vanishing dissipation that was mentioned above, it still allows for a coexistence of PBAs, given a careful choice of initial states.
The numerical protocol followed to analyze multistability for Eq. (A1) in terms of PBAs is described next. First, the initial data have been drawn uniformly in the two disjoint domains D_{1} and D_{2} of the $(x,\dot{x})$ plane, with
A total of 6000 initial data from each domain were propagated according to Eq. (A1). The ODE was integrated using a Runge–Kutta fourthorder method with a constant time step $\mathrm{\Delta}t={\mathrm{10}}^{\mathrm{2}}$, generating a total of 12 000 trajectories and keeping 10^{6} data points for each, after removal of the transient.
The majority of the initial data taken in the smaller domain D_{2} leads to a quasiperiodic orbit, while each of the 6000 initial data taken in D_{1} leads to a chaotic trajectory – as do a few “rare” initial data from D_{2}. An example of such a quasiperiodic trajectory is shown in blue in Fig. A1a, within the $(x,\dot{x})$ phase plane. This blue trajectory is superimposed upon a red, chaotic trajectory emanating from an initial point taken in D_{1}. The corresponding power spectra are shown in Fig. A1b, with the same blue and red color coding. The chaotic trajectory is clearly more diffuse within the phase plane than its quasiperiodic counterpart, and its power spectrum is quite a bit noisier. Both of these features are well known to be symptomatic of deterministic chaos (Eckmann and Ruelle, 1985).
By allowing the quasiperiodic trajectories that emanate from D_{2} to evolve up to t=2000, one obtains the set of blue points shown in Fig. A2. Somewhat surprisingly, this set does not form a closed curve: each blue dot in Fig. A2 actually corresponds to the state at t=2000 in the phase plane of a quasiperiodic orbit. One such orbit is represented in blue in Fig. A1a, after removal of the transient dynamics. Each blue dot in Fig. A2 corresponds to a different quasiperiodic orbit, whose frequency characteristics may change slightly from one blue dot to another. All these quasiperiodic orbits share, however, a spectral signature that resembles the one shown by the blue curve in Fig. A1b.
To illustrate further the distinction between quasiperiodic and chaotic orbits, the return maps for the minima of the x(t) variable have been computed. As is well known (e.g., Strogatz, 2015), if the return map contains just one point, the solution is periodic in time, with all minima having the exact same value, and the period of the oscillation can be estimated by calculating the time interval between two consecutive minima. If the return map contains continuouslooking curves that fill up with more and more points as the length of the orbit increases, the solution is quasiperiodic, while the presence of folds and selfsimilarity in the return map provides strong evidence for chaotic solutions. For the blue and red trajectories of Fig. A1a, we plot the corresponding return maps in Fig. A3a and b, respectively. The two plots clearly discriminate between the quasiperiodic nature of the former and the chaotic nature of the latter solution.
The initial data taken in D_{1}, when allowed to flow according to Eq. (A1), lead to a totally different local PBA that is formed by the red points shown in Fig. A2. Although the approximation of this local PBA shown here is relatively sparse, one can clearly discern the fact that its constitutive points are arranged according to a stretching and folding pattern that is typical of nonlinear, chaotic dynamics in the autonomous as well as in the nonautonomous setting (Chekroun et al., 2018).
The features that we find here to be exhibited by the local chaotic PBA are highly reminiscent of those that were obtained, in this deterministic case, by applying a standard Poincaré section analysis; cf. Dudkowski et al. (2016, Fig. 17c). In the presence of noise, though, the fine structure of the PBAs that results from stretching and folding in phase space is still captured by the PBA framework (Chekroun et al., 2018), whereas a Poincarémap approach would lead only to a cloud of points with no particular geometric structure. This statement was numerically illustrated by Chekroun et al. (2011) in their Fig. 7, by contrasting the upperright panel vs. the six lower panels of that figure.
Finally, we emphasize that it is not the disjointness of the two domains, D_{1} and D_{2}, that leads to the two distinct types of PBA, chaotic and quasiperiodic.^{1} Indeed, as mentioned earlier, even though the area of D_{2} is small, it still contains initial data whose evolution lands within the local PBA associated with chaos, i.e., with the other red points shown in Fig. A2.
SP devised the study, carried out the computations and wrote the first draft of the paper. MG helped integrate the work into the broader picture of applying the theory of nonautonomous and random dynamical systems to the study of climate variability and climate change. MDC devised the appendix to further illustrate the coexistence of local PBAs that are chaotic and nonchaotic, and carried out the corresponding computations. All three authors contributed to the writing of the paper.
Stefano Pierini is a member of the editorial board of the journal. The authors declare no conflict of interest.
This article is part of the special issue “Numerical modeling, predictability and data assimilation in weather, ocean and climate: A special issue honoring the legacy of Anna Trevisan (1946–2016)”. It is a result of a Symposium Honoring the Legacy of Anna Trevisan – Bologna, Italy, 17–20 October 2017.
Stefano Pierini would like to thank the “Università di Napoli Parthenope” for having supported his visit to the University of California at Los Angeles in August 2017 through grants D.R. 539 (2962016) and D.R. 953 (28112016, DSTE315) and for the partial support provided by the grant DSTE315B. Stefano Pierini and Michael Ghil gratefully acknowledge partial support from the MOMA project (PNRA1600196B) of the Italian P.N.R.A. This work has been partially supported by the Office of Naval Research (ONR) Multidisciplinary University Research Initiative (MURI) grants N000141210911 and N000141612073 (Mickaël D. Chekroun and Michael Ghil). Mickaël D. Chekroun also gratefully acknowledges support from the National Science Foundation (NSF) grants OCE1243175, OCE1658357 and DMS1616981.
This paper is dedicated to the memory of Anna Trevisan and to her
contributions to the applications of dynamical systems theory to the climate
sciences.
Edited by: Juan Manuel Lopez
Reviewed by: two anonymous referees
Arnold, L.: Random Dynamical Systems, Springer, Berlin, Germany, 1998. a
Bódai, T. and Tél, T.: Annual variability in a conceptual climate model: Snapshot attractors, hysteresis in extreme events, and climate sensitivity, Chaos, 22, 023110, https://doi.org/10.1063/1.3697984, 2012. a, b
Bódai, T., Károlyi, G., and Tél, T.: Driving a conceptual model climate by different processes: Snapshot attractors and extreme events, Phys. Rev. E, 87, 022822, https://doi.org/10.1103/PhysRevE.87.022822, 2013. a, b
Boyles, R. and Gardner, W. A.: Cycloergodic properties of discreteparameter nonstationary stochastic processes, IEEE T. Inform. Theory, 29, 105–114, 1983. a
Carrassi, A., Ghil, M., Trevisan, A., and Uboldi, F.: Data assimilation as a nonlinear dynamical systems problem: Stability and convergence of the predictionassimilation system, Chaos, 18, 023112, https://doi.org/10.1063/1.2909862, 2008a. a
Carrassi, A., Trevisan, A., Descamps, L., Talagrand, O., and Uboldi, F.: Controlling instabilities along a 3DVar analysis cycle by assimilating in the unstable subspace: a comparison with the EnKF, Nonlin. Processes Geophys., 15, 503–521, https://doi.org/10.5194/npg155032008, 2008b. a
Carvalho, A., Langa, J. A., and Robinson, J.: Attractors for InfiniteDimensional NonAutonomous Dynamical Systems, Springer, New York, USA, 2012. a, b
Celletti, A.: Periodic and quasiperiodic attractors of weaklydissipative nearlyintegrable systems, Regul. Chaotic Dyn., 14, 49–63, 2009. a
Chekroun, M. D., Simonnet, E., and Ghil, M.: Stochastic climate dynamics: Random attractors and timedependent invariant measures, Physica D, 240, 1685–1700, 2011. a, b, c, d
Chekroun, M. D., Ghil, M., and Neelin, J. D.: Pullback attractor crisis in a delay differential ENSO model, in: Advances in Nonlinear Geosciences, edited by: Tsonis, A., 1–33, Springer, available at: https://link.springer.com/chapter/10.1007/9783319588957_1, last access: 7 September 2018. a, b, c, d, e, f
De Cruz, L., Schubert, S., Demaeyer, J., Lucarini, V., and Vannitsem, S.: Exploring the Lyapunov instability properties of highdimensional atmospheric and climate models, Nonlin. Processes Geophys., 25, 387–412, https://doi.org/10.5194/npg253872018, 2018. a
Drótos, G., Bódai, T., and Tél, T.: Probabilistic concepts in a changing climate: A snapshot attractor picture, J. Climate, 28, 3275–3288, 2015. a, b
Drótos, G., Bódai, T., and Tél, T.: On the importance of the convergence to climate attractors, Eur. Phys. J.Spec. Top., 226, 2031–2038, 2017. a
Dudkowski, D., Jafari, S., Kapitaniak, T., Kuznetsov, N. V., Leonov, G. A., and Prasad, A.: Hidden attractors in dynamical systems, Phys. Rep., 637, 1–50, 2016. a, b, c, d
Duffing, G.: Erzwungene Schwingungen bei veränderlicher Eigenfrequenz und ihre technische Bedeutung, vol. 41/42 of Sammlung Vieweg, R. Vieweg & Sohn, Braunschweig, Germany, 1918. a
Eckmann, J.P. and Ruelle, D.: Ergodic theory of chaos and strange attractors, Rev. Modern Phys., 57, 617–656, 1985. a, b
Feudel, U., Grebogi, C., Hunt, B. R., and Yorke, J. A.: Map with more than 100 coexisting lowperiod periodic attractors, Phys. Rev. E, 54, 71, https://doi.org/10.1103/PhysRevE.54.71, 1996. a
Ghil, M.: Cryothermodynamics: the chaotic dynamics of paleoclimate, Physica D, 77, 130–159, 1994. a
Ghil, M.: A mathematical theory of climate sensitivity or, How to deal with both anthropogenic forcing and natural variability?, in: Climate Change: Multidecadal and Beyond, edited by: Chang, C. P., Ghil, M., M., L., and Wallace, J. M., 31–51, World Scientific Publ. Co./Imperial College Press, Singapore, 2015. a
Ghil, M.: The winddriven ocean circulation: Applying dynamical systems theory to a climate problem, Discrete Cont. Dyn.A, 37, 189–228, 2017. a, b
Ghil, M. and Wolansky, G.: NonHamiltonian perturbations of integrable systems and resonance trapping, SIAM J. Appl. Math., 52, 1148–1171, 1992. a
Ghil, M., Chekroun, M. D., and Simonnet, E.: Climate dynamics and fluid mechanics: Natural variability and related uncertainties, Physica D, 237, 2111–2126, 2008. a, b
Grebogi, C., Ott, E., and Yorke, J. A.: Chaos, strange attractors, and fractal basin boundaries in nonlinear dynamics, Science, 238, 632–638, 1987. a
Hilborn, R. C.: Chaos and Nonlinear Dynamics, Oxford University Press, Oxford, UK, 2000. a
Jackson, E. A.: Perspectives of Nonlinear Dynamics, Cambridge University Press, New York, USA, 1991. a, b
Jiang, S., Jin, F.F., and Ghil, M.: Multiple equilibria, periodic, and aperiodic solutions in a winddriven, doublegyre, shallowwater model, J. Phys. Oceanogr., 25, 764–786, 1995. a
Kloeden, P. E. and Rasmussen, M.: Nonautonomous Dynamical Systems, 176, Amer. Math. Soc., 2011. a, b
Kozlov, A., Sushchik, M., Molkov, Y., and Kuznetsov, A.: Bistable phase synchronization and chaos in system of coupled Van der PolDuffing oscillators, Int. J. Bifurcat. Chaos, 9, 2271–2277, 1999. a
Kuznetsov, A., Stankevich, N., and Turukina, L.: Coupled van der Pol–Duffing oscillators: Phase dynamics and structure of synchronization tongues, Physica D, 238, 1203–1215, 2009. a
Le Treut, H. and Ghil, M.: Orbital forcing, climatic interactions, and glaciation cycles, J. Geophys. Res., 88, 5167–5190, https://doi.org/10.1029/JC088iC09p05167, 1983. a
Lucarini, V., Ragone, F., and Lunkeit, F.: Predicting climate change using response theory: Global averages and spatial patterns, J. Stat. Phys., 166, 1036–1064, 2017. a
Manneville, P. and Pomeau, Y.: Intermittency and the Lorenz model, Phys. Lett. A, 75, 1–2, 1979. a
Nicolis, G.: Introduction to Nonlinear Science, Cambridge University Press, Cambridge, UK, 1995. a
Ott, E.: Chaos in Dynamical Systems, Cambridge University Press, Cambridge, UK, 2002. a, b
Pierini, S.: A Kuroshio Extension system model study: Decadal chaotic selfsustained oscillations, J. Phys. Oceanogr., 36, 1605–1625, 2006. a, b
Pierini, S.: Lowfrequency variability, coherence resonance, and phase selection in a loworder model of the winddriven ocean circulation, J. Phys. Oceanogr., 41, 1585–1604, 2011. a, b, c, d, e, f, g, h, i
Pierini, S.: Stochastic tipping points in climate dynamics, Phys. Rev. E, 85, 027101, https://doi.org/10.1103/PhysRevE.85.027101, 2012. a
Pierini, S.: Ensemble simulations and pullback attractors of a periodically forced doublegyre system, J. Phys. Oceanogr., 44, 3245–3254, 2014. a, b, c, d, e, f, g, h, i, j
Pierini, S. and Dijkstra, H. A.: Lowfrequency variability of the Kuroshio Extension, Nonlin. Processes Geophys., 16, 665–675, https://doi.org/10.5194/npg166652009, 2009. a
Pierini, S., Dijkstra, H. A., and Riccio, A.: A nonlinear theory of the Kuroshio Extension bimodality, J. Phys. Oceanogr., 39, 2212–2229, 2009. a
Pierini, S., Ghil, M., and Chekroun, M. D.: Exploring the pullback attractors of a loworder quasigeostrophic ocean model: the deterministic case, J. Climate, 29, 4185–4202, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y
Pisarchik, A. N. and Feudel, U.: Control of multistability, Phys. Rep., 540, 167–218, 2014. a
Pomeau, Y. and Manneville, P.: Intermittent transition to turbulence in dissipative dynamical systems, Commun. Math. Phys., 74, 189–197, 1980. a
Rasmussen, M.: Attractivity and Bifurcation for Nonautonomous Dynamical Systems, Springer, Berlin, Germany, 2007. a
Romeiras, F. J., Grebogi, C., and Ott, E.: Multifractal properties of snapshot attractors of random maps, Phys. Rev. A, 41, 784, https://doi.org/10.1103/PhysRevA.41.784, 1990. a
Roques, L. and Chekroun, M. D.: Probing chaos and biodiversity in a simple competition model, Ecol. Complex., 8, 98–104, 2011. a
Shannon, C. E.: The mathematical theory of communication, Bell Syst. Tech. J., 27, 379–423, 1948. a
Strogatz, S. H.: Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering, CRC Press, Boca Raton, FL, USA, 2015. a, b
Sushama, L., Ghil, M., and Ide, K.: Spatiotemporal variability in a midlatitude ocean basin subject to periodic wind forcing, Atmos. Ocean, 45, 227–250, https://doi.org/10.3137/ao.450404, 2007. a
Tél, T. and Gruiz, M.: Chaotic Dynamics, Cambridge University Press, Cambridge, UK, 2006. a
Trevisan, A. and Palatella, L.: On the Kalman Filter error covariance collapse into the unstable subspace, Nonlin. Processes Geophys., 18, 243–250, https://doi.org/10.5194/npg182432011, 2011. a
Trevisan, A. and Uboldi, F.: Assimilation of standard and targeted observations within the unstable subspace of the observationanalysisforecast cycle system, J. Atmos. Sci., 61, 103–113, https://doi.org/10.1175/15200469(2004)061<0103:AOSATO>2.0.CO;2, 2004. a
Van der Pol, B.: A theory of the amplitude of free and forced triode vibrations, Radio Rev., 1, 701–710, 1920. a, b
Van der Pol, B.: On relaxationoscillations, The London, Edinburgh and Dublin Phil. Mag. J. Sci., 2, 978–992, 1926. a
Vannitsem, S.: Stochastic modelling and predictability: analysis of a loworder coupled ocean–atmosphere model, Philos. T. Roy. Soc. A, 372, 20130282, https://doi.org/10.1098/rsta.2013.0282, 2014. a
Vannitsem, S. and De Cruz, L.: A 24variable loworder coupled ocean–atmosphere model: OAQGWS v2, Geosci. Model Dev., 7, 649–662, https://doi.org/10.5194/gmd76492014, 2014. a
Venkatesan, A. and Lakshmanan, M.: Bifurcation and chaos in the doublewell Duffing–Van der Pol oscillator: Numerical and analytical studies, Phys. Rev. E, 56, 6321, https://doi.org/10.1103/PhysRevE.56.6321, 1997. a
Zaslavsky, G. M. and Edelman, M.: Superdiffusion in the dissipative standard map, Chaos, 18, 033116, https://doi.org/10.1063/1.2967851, 2008. a
Other such quasiperiodic PBAs exist (not shown), whereas only one chaotic PBA seems to exist for the parameter regime analyzed herein.
 Abstract
 Introduction and motivation
 Model description
 The autonomous system
 The periodically forced system
 Further applications of crosscorrelation diagnostics
 Summary and conclusions
 Data availability
 Appendix A: Coexistence of pullback attractors in a Van der Pol–Duffing oscillator
 Author contributions
 Competing interests
 Special issue statement
 Acknowledgements
 References
 Abstract
 Introduction and motivation
 Model description
 The autonomous system
 The periodically forced system
 Further applications of crosscorrelation diagnostics
 Summary and conclusions
 Data availability
 Appendix A: Coexistence of pullback attractors in a Van der Pol–Duffing oscillator
 Author contributions
 Competing interests
 Special issue statement
 Acknowledgements
 References