1600-0870 Tellus A: Dynamic Meteorology and Oceanography 1600-0870 Stockholm University Press 10.16993/tellusa.25 VoR Original research paper On Equilibrium Fluctuations https://orcid.org/0000-0002-2308-6834 von Storch Jin-Song jin-song.von.storch@mpimet.mpg.de 1 2 Max Planck Institute for Meteorology, Hamburg, Germany Center for Earth System Research and Sustainability (CEN), University of Hamburg, Germany 20 09 2022 2022 74 2022 364 381 13 12 2021 30 08 2022 Copyright: © 2022 The Author(s) 2022 This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (CC-BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. See http://creativecommons.org/licenses/by/4.0/.

This paper considers a dynamical system described by a multidimensional state vector x. A component x of x evolves according to dx/dt = f(x). Equilibrium fluctuations are fluctuations of an equilibrium solution x(t) obtained when the system is in its equilibrium state reached under a constant external forcing. The frequencies of these fluctuations range from the major frequencies of the underlying dynamics to the lowest possible frequency, the frequency zero. For such a system, the known feature of the differential operator d()dt \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{{d( \cdot )} \over {dt}}}$ \end{document} as a high-pass filter makes the spectrum of f to vanish not only at frequency zero, but de facto over an entire frequency range centered at frequency zero (when considering both positive and negative frequencies). Consequently, there is a non-zero portion of the total equilibrium variance of x that cannot be determined by the differential forcing f. Instead, this portion of variance arises from many impulse-like interactions of x with other components of x, which are received by x along an equilibrium solution over time. The effect of many impulse-like interactions can only be realized by integrating the evolution equations in form of dx/dt = f(x) forward in time. This integral effect is not contained in, and can hence not be explained by, a differential forcing f defined at individual time instances.

low-frequency internal variability equilibrium variance and its decomposition in frequency domain spectrum at frequency zero differential forcing
1 Introduction

This paper studies the fluctuations of an equilibrium solution reached by a dynamical system under constant external forcing conditions. Applying this definition to climate, the fluctuations of our concern correspond to internal climate variability (ICV) resulting from internal processes without involving any time-varying external forcing. The existence of ICV has been known for a long time (see e.g. a review by Mitchell (1976)). A vast literature has been devoted to the description and the qualification of ICV, and to the understanding of mechanisms responsible for ICV. ICV has different appearances on different temporal and spatial scales. With respect to temporal scales, there was a great excitement generated by the discovery of James and James (1989) that internal dynamics of the atmosphere produce not only variations on the synoptic time scales of days but also variations on much longer time scales of about a few decades. Since then decadal and multi-decadal climate variations have drawn attention of many researchers. With respect to spatial and temporal scales, a variety of modes, ranging from atmospheric teleconnections (Panagiotopoulos et al., (2002) to the coupled modes, such as the El Niño Southern Oscilation (ENSO) (Rasmusson & Carpenter, 1982) in the tropical Pacific, and to the large-scale modes of extra-tropical SST, such as the Pacific Decadal Oscillation (PDO) (Newman et al., 2016) and the Atlantic Multidecadal Oscillation (AMO) (Kerr, 2000), has been investigated using both observations and model simulations.

Despite all the efforts, our understanding about ICV is still incomplete. An indication of the incompleteness can be found in a study on two atmospheric modes simulated by a coupled atmosphere-ocean general circulation model (GCM) (J.-S. von Storch, 1999a). These modes vary in phase with the two parts of the global atmospheric axial angular momentum M on time scales of months to a century. Because of the relation to M, and because M is completely determined by its torque F via dMdt=F \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{{dM} \over {dt}}} = F$ \end{document} , it was thought that the cause of the low-frequency variations of the two modes can be understood by examining the torque F. This turns out to be a wishful thinking. Against the expectation, the coherence and the phase spectra between M and F break down on time scales longer than a few months (bottom panel in Figure 1), inhibiting establishing any link between M and F on these scales. The result is irritating, since the angular momentum equation dMdt=F \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{{dM} \over {dt}}} = F$ \end{document} is one of those first principles in which our deterministic thinking is rooted.

Spectra of the global atmospheric axial absolute angular momentum M (top left) and the global torque F (top right), and the coherence (bottom left) and phase (bottom right) spectra between M and F. The spectra (adopted from (J.-S. von Storch 1999a) consist of a high-frequency part derived from 10-year daily data and a low-frequency part derived from 500-year monthly data, produced by a coupled atmosphere-ocean GCM. A further analysis using 200-year half-daily data produces essentially the same result (not shown). Despite the strong sampling errors, the same characteristics, namely the essentially white low-frequency plateau of M, a tendency for the spectrum of F to decrease with decreasing frequency at not too low frequencies before becoming white at the lowest frequencies, the coherent 90-phase relation at not too low frequencies and the break-down of the coherence at the lowest frequencies, are also found for the Lorenz solution shown in Figure 7.

Spectra of simulated global axial angular momentum and its torque

In this paper, we aim to clarify the mystery shown in Figure 1 and identify the mechanism responsible for ICV in a possibly general manner. Broadly speaking, there are, to my knowledge, two avenues that can lead to a more general explanation of climate variability. The first one emphasizes the role of nonlinear dynamics and is built on the dynamical systems approach. As reviewed by Dijkstra and Ghil (2005) for low-frequency variability of the ocean circulation, the observed ocean variability can be interpreted as emerging from the road of successive bifurcations leading through multiple equilibria to oscillatory and eventually chaotic solutions. What is less clear with this dynamical systems approach is the extent to which it can be applied to purely internally generated variability. A bifurcation behavior has to be initiated by varying a parameter of the considered system. For the ocean, the parameters studied in this context are wind stress and buoyancy fluxes at the sea surface. As both parameters are also the external forcing of the ocean, the distinction between internally and externally generated variability becomes blurred.

The second avenue emphasizes the stochastic nature and is culminated in the concept of stochastic climate models proposed by Hasselmann (1976). Following this concept, variations of a slow climate component are the response to continuous random excitation by fast climate components (such as the weather). The evolution equation of a slow climate component includes a stochastic forcing representing the fast components and takes the form of a stochastic differential equation. Such evolution equations suggest that the internal climate variability is stochastic and does not result from deterministic responses to a temporally varying external forcing. Despite the wide acceptance of the concept, it is still not definitely clarified what stochasticity means. The meaning of stochasticity is of fundamental importance: A physical phenomenon, such as ICV, cannot be governed by two different types of evolution laws, one formulated in terms of deterministic differential equations and the other in terms of stochastic differential equations. If we accept that all climate components, no matter fast or slow, result from deterministic laws (such as the conservation of angular momentum) that do not contain any stochastic element, a stochastic model can only be an approximation at most. To “upgrade” a stochastic model to something having the same status as the governing deterministic laws, one needs first to show that there is a physical “thing” that escapes the control of the governing deterministic evolution laws.

This papers shows that for dynamical systems of our concern such a “thing” exists. A dynamical system of our concern has the following three properties. First, the system described by the state vector x is governed by the Newtonian type of dynamics such that the temporal evolution of a component x of x is determined by

dxdt=f(x). \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\frac{{dx}}{{dt}} = f({\bf{x}}).$ \end{document}

d()dt \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{{d( \cdot )} \over {dt}}}$ \end{document} is the time derivative operator. The set of Eq.(1), formulated for all components of x, represents the deterministic laws governing the evolution of x. f = f(x) is a deterministic non-linear function of x, which contains a constant external forcing as well as a dissipation. A dissipation is needed to counteract the constant external forcing, so that the system reaches an equilibrium state under the constant external forcing, rather than runs away. Hereafter f is referred to as the differential forcing of x, also later when discrete version of Eq.(1) is considered. Eq.(1) represents the classical Newtonian dynamics in the sense that the time rate of change of a “momentum”, i.e. x in Eq.(1), equals the applied “force”, i.e. f(x) in Eq.(1).

Secondly, the system’s state vector x must be multidimensional so that the differential forcing of x, f(x), depends not only on x but also on other components of x. For such a multidimensional system, a differential forcing f of component x at time t describes the interaction of x with other components of x at time t.

Finally, when forced with a constant external forcing for a while, the system reaches an equilibrium state described by a time-varying equilibrium solution. An equilibrium solution fluctuates stationarily, non-periodically, and for ever into infinite time (when left alone). The fluctuations of an equilibrium solution are referred to as the equilibrium fluctuations. Equilibrium fluctuations have finite variances, referred to as the equilibrium variances. Equilibrium variances are determined by the constant external forcing, independent of the initial state, from which the system starts its “journey” toward the equilibrium state. For any component x(t) of an equilibrium solution x(t), x(t) at two time points become increasingly less well correlated with increasing separation between the two time points. As a result, all covariance functions of variables, that are functions of an equilibrium solution x(t) (e.g. the auto-covariance function of a component x of x(t), the auto-covariance function of the differential forcing f of x, and the cross-covariance function between x and f), decay fast in magnitude with increasing time lag.

An example of systems having the just stated properties is an atmospheric (or an oceanic or a coupled atmosphere-ocean) GCM. In this case, the dimension of x is determined by the grid used to discretize the governing equations of motions, which are of the type of Newtonian dynamics and can be transformed into a set of equations in form of Eq.(1). A component x of x is a prognostic variable at one grid point. The differential forcing f of x is a function not only of x but also of other variables at the same and the surrounding grid points. When forcing the GCM with a constant external forcing (e.g. from the sun) for some time, we obtain an equilibrium climate. When rerunning the same GCM under the same external forcing, we obtain the same climate characterized by the same equilibrium variances.

Another example is a gas as a classical many-particle system. Obviously, the governing dynamics are Newtonian and can be described by a set of equations in form of Eq.(1). A component x of x is the position or velocity of one particle. The differential forcing f of x is a function not only of x but also of positions and velocities of other particles. When putting a gas in contact with a heat reservoir as an external forcing, it reaches an equilibrium state characterized by a constant mean kinetic energy of the gas molecules. Independent of the initial states of the individual gas molecules, a gas brought in contact with the same heat reservoir reaches the same equilibrium state characterized by the same mean kinetic energy.

A final example is the 3-dimensional Lorenz model (Lorenz, 1963). The model with the parameters given in Lorenz (1963) is derived from a description of convective motions within a fluid layer. The external forcing, which represents the difference in temperature at the upper and lower boundaries of the layer, is merged into one of the model parameters. The Lorenz model is governed by three equations in form of Eq.(1). The differential forcing f of one component x is a function not only of x but also of the other components of the state vector x. With the same parameter (corresponding to the same external forcing), two equilibrium solutions starting from two different initial states fluctuate in similar manner (top row of Figure 2) and have the same variances. The auto-covariance functions of x (bottom panel in Figure 5) and f (top panel in Figure 6), as well as the cross-covariance function between x and f (bottom panel in Figure 6), decay fast in magnitude with increasing time lag.

Pieces of two equilibrium solutions (left and right column) obtained by integrating the Lorenz system from two different initial conditions. Top panel shows the three components (green, red, blue) of the two equilibrium solutions xs and the bottom one the respective differential forcing. For each record, the respective record mean is removed. Here and in all other figures, an equilibrium solution is obtained in two steps. The first step consists of a spin-up integration, that starts from a near-equilibrium state and lasts 5 × 103 time steps. A near-equilibrium state is a state obtained after a long integration (over 5 × 105 time steps), that is then perturbed by adding a small normally distributed random variable. In the second step, an equilibrium solution is produced by integrating the model from the end state of the spin-up integration. Throughout this paper, the Lorenz model is integrated using a Runge-Kutta method with a time step of 0.01.

Pieces of two equilibrium Lorenz solutions

To this end we note that the system of our concern is a simpler and idealized version of the reality. Take the climate as an example. The real climate, which is always subjected to some time-varying external forcing (e.g. to a changing atmospheric greenhouse gas concentration that leads to a changing radiative forcing of the system), is generally not in an equilibrium state. A climate subjected to a constant forcing is an artificially constructed and simpler version of the real climate. We believe that the real climate cannot be fully understood without first fulling understand this idealized climate.

We note further that even though some properties of equilibrium fluctuations resemble those produced by a chaotic dynamical system, we intend not to equal the systems of our concern with chaotic dynamical systems (Eckmann & Ruelle, 1985). The main reason is that the latter has been an acronym for properties associated with an attractor, which are not the concern of the present study.

In this paper we should use the Lorenz model to illustrate and support our ideas. We do so not only because the Lorenz model belongs to the dynamical systems of our concern, but also because of the limited computer resources. Given that the number of simulated years per day of a typical modern climate model is around 50 (Balaji et al., 2017), and given that we do need at least multicentury-long, if not millennia-long, simulations for studying the climate variability at the lowest frequencies, it is not feasible to produce a sufficiently large number of sufficiently long solutions of a climate model with the present computer facility. Such excessive computations, which are needed to avoid sampling errors, such as those found in Figure 1, can be carried out for the 3-dimensional Lorenz model.

The paper is organized as follows. After some preliminary notions introduced in Section 2, Section 3 discusses two requisites needed for the considerations in Section 4. The first one is a clean definition of spectrum (Section 3.1). The second is a specification of the low-frequency shape of spectrum for variables that are functions of an equilibrium solution x(t) (Section 3.2). We show that for these variables, whose covariance functions decay fast in magnitude with increasing time lag, the spectra of x and f and the cross-spectrum between x and f are white at the lowest frequencies. In Section 4, equipped with the spectrum definition identified in Section 3.1, we derive the known spectral relations between x and f arising from the differential operator in Eq.(1). Equipped with the low-frequency spectral shape identified in Section 3.2, we are able to point out that for component x having a white low-frequency spectral plateau, the variance related to the spectral plateau, that makes up a non-zero portion of the total equilibrium variance of x, cannot be determined by f, since f de facto ceases to fluctuate at near-zero frequencies. Section 5 discusses ideas suggesting why x fluctuates at the lowest frequencies, and what, if not the differential forcing f, determines the low-frequency spectral plateau of x. Conclusions and a discussion on the implication of this paper is provided in Section 6.

2 Preliminaries

Since a system described by Eq.(1) can generally only be solved numerically, we consider a discretized version of Eq.(1),

xsδtx(s1)δtδt=fx(sδt), \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\frac{{x\left( {s\delta t} \right) - x\left( {(s - 1)\delta t} \right)}}{{\delta t}} = f\left( {{\bf{x}}(s\delta t)} \right),$ \end{document}

where s ∈ ℤ is an integer and t = sδt indicates the time at the s-th time step. For the sake of simplicity, we set the integration time step δt to one and consider,

xsxs1=f(xs1)fs1. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${x_s} - {x_{s - 1}} = f({{\bf{x}}_{s - 1}}) \equiv {f_{s - 1}}.$ \end{document}

f is also known as the first difference of x. The limit δt → 0 turns out to be irrelevant and will be briefly discussed at a later point. Hereafter, for a clear distinction between discrete and continuous variables, we denote a discrete variable by a subscript, e.g. xs, and a continuous variable by an argument, e.g. x(t). For an easy indexing, we write a finite discrete record of a variable x as

{xs,N}=xssN,,1,0,1,,N, \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\{ {x_{s,N}}\} = \left\{ {{x_s}\left| {s \in \left\{ { - N, \cdots , - 1,0,1, \cdots ,N} \right\}} \right.} \right\},$ \end{document}

which contains total 2N+1 time steps. Here and hereafter, the subscript N denotes a quantity related to a finite record. In theory, N in Eq.(4) can go to infinity, since an equilibrium solution can, when left alone, be extended to infinity. In practice, N is always finite.

Two operators are needed. The first one is

E()limmEm(), \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $E( \cdot ) \equiv \mathop {\lim }\limits_{m \to \infty } {E_m}( \cdot ),$ \end{document}

where Em is an average over quantities related to m equilibrium solutions. The m equilibrium solutions correspond to m points in the respective phase space. They are obtained by integrating the dynamical system from m initial states that are different but belong to the same equilibrium state. Thus, E(∙) takes all reachable equilibrium solutions of the considered equilibrium state into account. Note that E(∙) is purely deterministic. The link of E(∙) to the expectation operator ɛ(·) used in probability theory is discussed in Appendix A.

The second operator is the Fourier transform F{} \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\cal F}\left\{ \cdot \right\}$ \end{document} , or the discrete Fourier transform when dealing with discrete variables. The Fourier transform F{} \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\cal F}\left\{ \cdot \right\}$ \end{document} of a discrete real-valued function α indexed by integer τ,F{α} \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\tau \in ,{\cal F}\{ \alpha \}$ \end{document} , is defined as

F{α}(ω)τ=ατe2πiτω, \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\cal F}\{ \alpha \} (\omega ) \equiv \sum\limits_{\tau = - \infty }^\infty {{\alpha _\tau }} {e^{ - 2\pi i\tau \omega }},$ \end{document}

which operates on infinite values of ατ and results in a function of ω defined on the real interval [–1/2, 1/2] (H. von Storch & Zwiers, 1999). If α is a solution of a dynamical system, τ indexes the time step. If α is a covariance function, τ indexes the time lag. F{α} \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\cal F}\{ \alpha \}$ \end{document} is generally a complex-valued function, but becomes real-valued, if ατ is symmetric about τ = 0. As pointed out by Priestley (1981), F{α} \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\cal F}\{ \alpha \}$ \end{document} does not always exist. A necessary condition for the existence of F{α} \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\cal F}\{ \alpha \}$ \end{document} is that α is absolutely summable, meaning

τ=|ατ|<. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\sum\limits_{\tau = - \infty }^\infty | {\alpha _\tau }| < \infty .$ \end{document}

A variable that is a function of an equilibrium solution is not absolutely summable. The spectrum that decomposes the variance of this variable in frequency domain must be defined as the Fourier transform of the auto-covariance function of this variable.

3 Two requisites

This section seeks for equilibrium solutions of a system of our concern, a) a general definition of spectrum, and b) a specification of the general low-frequency shape of spectra of variables that are functions of an equilibrium solution.

3.1 Spectrum definition

Throughout this subsection and if not specified otherwise, x denotes a variable related to an equilibrium solution, i.e. x can be a component of an equilibrium solution or a differential forcing. We first concentrate on the definition of auto-spectrum and extend the definition to cross-spectrum at the end of this subsection.

The spectrum of x decomposes the variance of x in frequency domain. In order to concentrate on variance, we consider throughout this paper only anomaly variables obtained by subtracting the respective mean from the variable. We assume, and as can be expected from many Newtonian type systems, that the variance is ergodic, i.e.

x2¯=E(x2) \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\overline {{x^2}} = E({x^2})$ \end{document}

with

x2¯limN x2¯N, x2¯N 1(2N+1)s=NNxs2, \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\overline {{x^2}} \equiv \mathop {\lim }\limits_{N \to \infty } \;{\left( {\overline {{x^2}} } \right)_N},{\rm{ }}{\left( {\overline {{x^2}} } \right)_N} \equiv {\rm{ }}\frac{1}{{(2N + 1)}}\sum\limits_{s = - N}^N {{{\left( {{x_s}} \right)}^2}} ,$ \end{document} E(x2)limm Em(x2),    Em(x2)1mi=1mxsi2. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $E({x^2}) \equiv \mathop {\lim }\limits_{m \to \infty } \;{E_m}({x^2}),\;\;\;\;{E_m}({x^2}) \equiv \frac{1}{m}\sum\limits_{i = 1}^m {{{\left( {x_s^i} \right)}^2}} .$ \end{document}

xs in Eq.(9a) is taken from a finite record {xs,N} indexed as in Eq.(4), with s running from –N to N. xsi \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $x_s^i$ \end{document} in Eq.(9b) indicates the value of x associated with the i-th equilibrium solution at an arbitrary time point s. Since x is a function of an equilibrium solution, the exact value of s in Eq.(9b) is irrelevant. In words, ergodicity in variance means that the variance obtained by averaging xs2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $x_s^2$ \end{document} over a sufficiently large number of time points equals the variance obtained by averaging xs2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $x_s^2$ \end{document} over a sufficiently large number of ensemble points. Figure 3 shows that for the three Lorenz components, the variance is indeed ergodic.

Variances (x2¯)N \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${(\overline {{x^2}} )_N}$ \end{document} (as defined in Eq.(9a)) obtained from a single Lorenz equilibrium solution via time average (solid) and variances Em(x2) (as defined in Eq.(9b)) obtained by averaging an ensemble of Lorenz equilibrium solutions (dashed), for the components of the 3-dimensional Lorenz state vector x (green, red, blue). Abscissa shows 2N with 2N + 1 being the length of the equilibrium solution used for time average, and m being the number of equilibrium solutions used for ensemble average.

Ergodicity in variance

The ergodicity in variance induces one to think that the spectrum is also ergodic, and can be derived by applying the discrete Fourier transform either to an auto-covariance function defined using a time average, or to an auto-covariance function defined using an ensemble average. The former is given by

cτlimNcτ,N \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${c_\tau } \equiv \mathop {\lim }\limits_{N \to \infty } {c_{\tau ,N}}$ \end{document}

with

cτ,N1(2N+1|τ|)s=N+|τ|Nxsτxs,  for τ0,1,,2N1(2N+1|τ|)s=NN|τ|xsτxs,  for τ2N,,1  \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${c_{\tau ,N}} \equiv \left\{ {\begin{array}{*{20}{c}} {\frac{1}{{(2N + 1 - |\tau |)}}\sum\nolimits_{s = - N + |\tau |}^N {{x_{s - \tau }}{x_s},{\rm{ for }}\ \tau \in \left\{ {0,1, \cdots ,2N} \right\}} }\\ {\frac{1}{{(2N + 1 - |\tau |)}}\sum\nolimits_{s = - N}^{N - |\tau |} {{x_{s - \tau }}{x_s},{\rm{ for }}\ \tau \in \left\{ { - 2N, \cdots , - 1} \right\}} \;} \end{array}} \right.$ \end{document}

and the latter by

γτExsτxs=limmγτ,m, \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\gamma _\tau } \equiv E\left( {{x_{s - \tau }}{x_s}} \right) = \mathop {\lim }\limits_{m \to \infty } {\gamma _{\tau ,m}},$ \end{document}

with

γτ,mEm(xsτxs). \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\gamma _{\tau ,m}} \equiv {E_m}({x_{s - \tau }}{x_s}).$ \end{document}

The two lines in Eq.(11) differ only with respect to the values taken by the running index s. Appendix B (below Eq.(B14)) explains why 1/(2N + 1 – |τ|) is chosen as the normalization factor in Eq.(11), even though other factors are used in statistics literature (H. von Storch & Zwiers, 1999). As in Eq.(9a), xs in Eq.(11) is taken from a finite record {xs,N}. As in Eq.(9b), the choice of s in Eq.(13) is irrelevant for x that is a function of an equilibrium solution.

Both cτ and γ| result from an average of two-time products xs–txs. Similar to the variance that can be obtained by averaging the product xs2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $x_s^2$ \end{document} either over time points or over ensemble points, an auto-covariance at lag τ can be obtained by averaging the two-time product xs–τxs either over time points following Eq.(11) and Eq.(10), or over ensemble points following Eq.(13) and Eq.(12). However, there is a big difference between the two averages. The former is an average over different numbers of two-time products depending on τ, since the number of available two-time products in a single record {xs,N} equals 2N – |τ| + 1 and hence decreases with increasing |τ|. On the contrary, the latter is an average over the same number of two-time products, regardless of the value of τ.

Assuming that auto-covariance function exists, either the limit given in Eq.(10) or the limit given in Eq.(12) must converge. If for any value of τ, the two limits converge to the same value, the auto-covariance function is ergodic. However, cτ,N with large values of τ cannot converge in the limit N → ∞ due to the decrease in the number of available two-time products with increasing time lag. For any value of N in the limit N → ∞, cτ,N at time lags close to N is always averaged over a limited number of two-time products. For τ = 2N, cτ,N reduces to the “raw” two-time product x–NxN without being subjected to any average! Thus, it is the lack of available two-time products at large time lags that makes it impossible for cτ,N to converge to zero with increasing N, even when cτ,N would converge to zero, if there were sufficiently many two-time products available for averaging. To ensure the convergence, we must increase the available two-time products, which is only possible by employing more equilibrium solutions.

For the Lorenz system, Figure 4 shows how cτ,N behaves with increasing N; and Figure 5 shows how γτ,m behaves with increasing m. Since both cτ,N and γτ,m are symmetric about τ = 0, only positive lags are shown. At small time lags with τ ≤ 1000 (left columns in Figures 4 and 5), cτ,N converges with increasing N towards γτ,m obtained with the largest value of m: cτ,N obtained with 2N = 106 (bottom left panel in Figure 4) is visually identical to γτ,m obtained with m = 106 (bottom left panel in Figure 5). However, we need – for defining a spectrum – the values of cτ,N and γτ,m at all τ with τ ∈ ℤ, since the Fourier transform is an operation on auto-covariance function at all τ with τ ∈ ℤ (see Eq.(6)). When considering all time lags (right columns in Figure 4 and Figure 5), the behavior of cτ,N with increasing N is clearly different from the behavior of γτ,m with increasing m. While γτ,m becomes increasingly smooth and the magnitudes of γτ,m at large time lags converge to zero with increasing m, cτ,N at large time lags remains erratic without any sign of convergence, reflecting the above described lack of two-time products available for time averaging.

Auto-covariance functions cτ,N (as defined in Eq.(11)), obtained by averaging for each value of τ the two-time products xs–τxs available from a single equilibrium solution of length 2N + 1 with 2N = 2 × 103 in the first, 2N = 104 in the second, and 2N = 105 in the third, and 2N = 106 in the last row, for the three components of the Lorenz state vector x (green, red, blue). The left column shows cτ,N at small time lags. The right column shows cτ,N at all time lags.

Auto-covariance functions of components x derived from time average

Auto-covariance functions γτ,m (as defined in Eq.(13)), obtained by averaging for each value of τ the two-time products xs–τxs from an ensemble of equilibrium solutions with the ensemble size being m = 102 in the first, m=104 in the second, m = 105 in the third, and m = 106 in the last row, for the three components of the Lorenz state vector x (green, red, blue). The left column shows γτ,m at small time lags. The right column shows γτ,m at all time lags.

Auto-covariance functions of components x derived from ensemble average

Auto-covariance functions of differential forcing f of component x (top), and cross-covariance functions between x and f (bottom), estimated as γτ,mf \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _{\tau ,m}^f$ \end{document} and γτ,mxf \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _{\tau ,m}^{xf}$ \end{document} following Eq.(13), for the three Lorenz components (green, red, blue). The left column shows these functions at short time lags; whereas the right column shows them at all considered time lags. γτ,mxf \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _{\tau ,m}^{xf}$ \end{document} is calculated using x and f of the same equilibrium solution. The average Em(·) is performed for an ensemble of size m = 106. Note that the cross-covariance function γτxf \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _\tau ^{xf}$ \end{document} is not exactly anti-symmetric according to Eq.(24).

Covariance functions derived from ensemble average

Given that only γτ,m converges for all τ ∈ ℤ, whereas cτ,N cannot converge for large values of τ ∈ ℤ, we are left with no other choice than defining spectrum as the Fourier transform of γ obtained by averaging over all reachable equilibrium solutions. Another way to see this is to recall that a necessary condition for the Fourier transform of an infinite series to exist is that the series is absolutely summable. cτ is not absolutely summable, since the magnitude of cτ remains large at large time lags due to the lack of two-time products available for averaging. Generally, the spectrum obtained from one equilibrium solution is erratic and differs from the spectrum obtained from another equilibrium solution.

Hereafter, a spectrum, denoted by Γ(ω), will be defined as the Fourier transform of γτ

Γ(ω)F{γ}=τ=γτei2πωτ=limmτ=γτ,mei2πωτ. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\Gamma (\omega ) \equiv {\cal F}\{ \gamma \} = \sum\limits_{\tau = - \infty }^\infty {{\gamma _\tau }} {e^{ - i2\pi \omega \tau }} = \mathop {\lim }\limits_{m \to \infty } \sum\limits_{\tau = - \infty }^\infty {{\gamma _{\tau ,m}}} {e^{ - i2\pi \omega \tau }}.$ \end{document}

Since Γ(ω) is symmetric, Γ(ω) over non-negative frequencies is considered throughout this paper. Appendix B shows that Γ(ω) of a variable x decomposes the variance of x as defined in Eq.(8). Appendix B also shows that Γ(ω) can be equally expressed as

Γ(ω)=limNE(Sω,N)=limNlimmEm(Sω,N), \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\Gamma (\omega ) = \mathop {\lim }\limits_{N \to \infty } E({S_{\omega ,N}}) = \mathop {\lim }\limits_{N \to \infty } \mathop {\lim }\limits_{m \to \infty } {E_m}({S_{\omega ,N}}),$ \end{document}

where Sω,N is the sample spectrum (periodogram) defined in terms of Fourier coefficients used to Fourier decompose a finite record {xs,N} (see Eq.(B5)). In this paper, all spectra of the Lorenz system are estimated following Eq.(15), rather than Eq.(14), since Γ(ω) can be better approximated using Em(Sω,N) with a moderate value of m than using auto-covariance function γτ,m with the same value of m. The situation comes about since the Fourier transform is sensitive to the values of an auto-covariance function at large time lags. Even when an auto-correlation function would decay to zero when sufficiently averaged, an estimate of this auto-correlation function, obtained by averaging over limited number of two-time products, has generally a error structure such that it varies slowly around zero at large time lags (Bartlett, 1946; Bloomfield, 1976). To accurately estimate spectrum as the Fourier transform of an auto-covariance function, this error structure must be sufficiently suppressed, which requires considering a large value of m. Figure 5 shows that γτ,m is still clearly non-zero at large time lags when employing an ensemble of size m = 104 (second panel in the right column), and only converges visually to zero for m larger than 105 (third and fourth panel in the right column) . In contrast to that, Appendix B (Figure B1) shows for the Lorenz components that Em(Sω,N) converges fairly fast with m and we are able to obtain a stable spectrum estimate already with m = 103.

To further reinforce our conclusion that spectrum should be defined by taking all equilibrium solutions into account, Appendix C shows for the Lorenz system that the sample spectrum Sω,N, that is defined for a single equilibrium solution, does not converge with increasing N. The convergence of Sω,N in the limit N → ∞ is needed to obtain a spectrum that decomposes the variance as defined in Eq.(9a).

Before leaving this subsection, we extend the definition given in Eq.(14) for an auto-spectrum to the definition of a cross-spectrum between two variables, e.g. the cross-spectrum between a component x and its differential forcing f, Γxf (ω). Γxf (ω) is defined in analogy to Eq.(14) as the Fourier transform of the cross-covariance function

γτxfExsτfs. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _\tau ^{xf} \equiv E\left( {{x_{s - \tau }}{f_s}} \right).$ \end{document}

Similar to Eq.(15), Γxf (ω) is identical to the double limit of the sample cross-spectrum Sω,Nxf \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $S_{\omega ,N}^{xf}$ \end{document} . The latter is defined in analogy to Eq.(B5) in terms of Fourier coefficient of x multiplied with the complex conjugate of the Fourier coefficient of f.

3.2 Specification of low-frequency spectral shape

A spectrum Γ(ω) defined in Eq.(14) has the spectral value at ω = 0

Γ(0)=τ=γτ. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\Gamma (0) = \sum\limits_{\tau = - \infty }^\infty {{\gamma _\tau }} .$ \end{document}

Γ(0) is linked to the low-frequency shape of Γ(ω). In this subsection, Γ(ω) can be the auto-spectrum of a variable (e.g. a component x or its differential forcing f), or the cross-spectrum between two variables (e.g. between x and f). The variables are functions of an equilibrium solution whose properties are specified in Section 1. For such variables, covariance functions γτ decay fast in magnitude with increasing time lag. This is illustrated for the Lorenz system in the bottom panel of Figure 5 for the auto-covariance functions of x, γτx \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _\tau ^x$ \end{document} , and in Figure 6 for the auto-covariance functions of f, γτf \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _\tau ^f$ \end{document} , and for the cross-covariance functions between x and f, γτxf \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _\tau ^{xf}$ \end{document} .

For a decaying covariance function, there exists a finite time lag τ, hereafter referred to as the cut-off time lag, and a negligibly small positive ϵ, such that

|γτ|<ϵ  for all τ>τ. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $|{\gamma _\tau }| < \;\;{\rm for}\ {\rm all}\;\tau > {\tau ^\dagger }.$ \end{document}

We consider ϵ being negligibly small when Γ(ω) in Eq.(14) can be well approximated by

Γ(ω)τ=ττγτei2πωτ. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\Gamma (\omega ) \simeq \sum\limits_{\tau = - {\tau ^\dagger }}^{{\tau ^\dagger }} {{\gamma _\tau }} {e^{ - i2\pi \omega \tau }}.$ \end{document}

Knowing the value of τ that allows the above approximation, we can define (when considering only posivite frequencies) the frequency interval

Ω=[0,ω0)ω: ω 12πτ, \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\Omega = [0,{\omega _0}) \equiv \left\{ {\omega :\;\omega \ll {\rm{ }}\frac{1}{{2\pi {\tau ^\dagger }}}} \right\},$ \end{document}

with

ω0   for  ωΩ. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\omega \simeq 0\;\;\;{\rm{for}}\;\;\omega \in \Omega .$ \end{document}

Ω covers the low-frequency end of the entire non-negative frequency interval [0,1/2]. The exact value of ω0, the right end of Ω, and with that the exact width of Ω, depend on τ, which in turn depends on the accuracy of Eq.(19). For γτ whose magnitude converges sufficiently fast to zero, τ is finite for a negligibly small ϵ. The width of Ω is not zero. The near-zero frequencies ω ϵ Ω are physically meaningful for equilibrium solutions, which can extend to infinite times and fluctuate on infinite timescales.

For ω ϵ Ω and τ , we have 2πωτ < 2πωτ ≪ 1, so that e–i2πωτ in Eq.(19) can be well approximated by one, since all higher-order terms in the Taylor expansion of e–i2πωτ are negligibly small relative to one. We obtain

Γ(ω)τ=ττγττ=γτ=Γ(0),   for ωΩ. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\Gamma (\omega ) \simeq \sum\limits_{\tau = - {\tau ^\dagger }}^{{\tau ^\dagger }} {{\gamma _\tau }} \simeq \sum\limits_{\tau = - \infty }^\infty {{\gamma _\tau }} = \Gamma (0),\;\;\;{\rm{for}}\;\omega \in \Omega .$ \end{document}

Thus, the low-frequency part of Γ(ω) is essentially white and equal the respective spectral value at frequency zero, Γ(0).

For a spectrum whose low-frequency part is white, Γ(0) not only represents the spectral value at a single frequency, ω = 0, but also determines the white spectral level at the lowest frequencies and the variance related to this white low-frequency spectral level.

4 Frailty of the differential forcing at near-zero frequencies

With the spectrum definition given in Eq.(14), we are now in the position to translate Eq.(3) formulated in time domain into spectral relations between x and f formulated in frequency domain. Such a relation can take the form of the cross-spectrum between x and f, Γxf, or a relation between the spectrum of x, Γx, and the spectrum of f, Γf. We derive the spectral relations following H. von Storch and Zwiers (1999, chapter 11.4.4). In their case, in which the spectral relation between a stochastic process Xt and its first difference Xt–Xt–1 is considered, a covariance function is indisputably defined in terms of an expectation operator ɛ(·) that involves joint probability functions (s. Appendix A). In our case, in which a dynamical system in its equilibrium state is considered, we have the choice of defining a covariance function either as cτ using just one equilibrium solution, or as γτ using all reachable equilibrium solutions. Once the choice of γτ is made, for the reasons given in the last section, the derivation can be carried out in the same way as in H. von Storch and Zwiers (1999) by replacing ɛ(·) by E(·).

We obtain, by shifting s in Eq.(3) to s–τ, multiplying the shifted Eq.(3) with Eq.(3), and then applying E(·),

γτf=2γτxγτ1xγτ+1x; \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _\tau ^f = 2\gamma _\tau ^x - \gamma _{\tau - 1}^x - \gamma _{\tau + 1}^x;$ \end{document}

and by multiplying Eq.(3) with xs–1–τ, and applying E(·),

γτxf=γτ+1xγτx. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _\tau ^{xf} = \gamma _{\tau + 1}^x - \gamma _\tau ^x.$ \end{document}

Applying Fourier transform to γx, γf, and γfx in Eq.(23) and Eq.(24) yields

Γf(ω)=Gf(ω) Γx(ω) \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\Gamma ^f}(\omega ) = {G^f}(\omega )\;{\Gamma ^x}(\omega )$ \end{document} Γxf(ω)=Gxf(ω) Γx(ω) \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\Gamma ^{xf}}(\omega ) = {G^{xf}}(\omega )\;{\Gamma ^x}(\omega )$ \end{document}

with the proportionality factors

Gf(ω)=(2sin(πω))2,   Gxf(ω)=(ei2πω1), \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${G^f}(\omega ) = {(2\sin (\pi \omega ))^2},\;\;\;{G^{xf}}(\omega ) = ({e^{i2\pi \omega }} - 1),$ \end{document}

where the following equality is used

τ=γτ±1ei2πωτ=e±i2πωτ=γτ±1ei2πω(τ±1)=e±i2πωτ=γτei2πωτ. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\sum\limits_{\tau = - \infty }^\infty {{\gamma _{\tau \pm 1}}} {e^{ - i2\pi \omega \tau }} = {e^{ \pm i2\pi \omega }}\sum\limits_{\tau = - \infty }^\infty {{\gamma _{\tau \pm 1}}} {e^{ - i2\pi \omega (\tau \pm 1)}} = {e^{ \pm i2\pi \omega }}\sum\limits_{\tau = - \infty }^\infty {{\gamma _\tau }} {e^{ - i2\pi \omega \tau }}.$ \end{document}

Thus, both Γf (ω) and Γxf (ω) are proportional to the spectrum of x, Γx (ω).

Γxf(ω) is a complex-valued function and can be described by its amplitude, known as the amplitude spectrum Axf(ω), or the coherence and phase spectrum, Cxf(ω) and Pxf(ω). Cxf(ω) is a normalized version of the amplitude of Γxf(ω), defined as (Axf(ω))2/(Γx(ω>f(ω)). Pxf(ω) equals the arctangent of the ratio of the imaginary to the real part of Γxf(ω). Given Γxf(ω) in Eq.(26), these spectra take the following forms

Axf(ω)=2sin(πω)Γx(ω) \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${A^{xf}}(\omega ) = 2\sin (\pi \omega ){\Gamma ^x}(\omega )$ \end{document} Pxfω =± π/2+πω, \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${P^{xf}}\left( \omega \right)\; = \pm {\rm{ }}\pi /2 + \pi \omega ,$ \end{document} Cxf(ω) =1,    for ω0. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${C^{xf}}(\omega ){\rm{ }} = 1,\;\;\;\;{\rm{for}}\;\omega \ne 0.$ \end{document}

Since the period of the tangent function is π, Pxf (ω) is not unique and given by two expressions in Eq.(30). Since Axf(ω) and Cxf(ω) are symmetric and Pxf(ω) is anti-symmetric, only positive frequencies need to be considered.

For sufficiently small ω satisfying 2πω ≪ 1, Gf(ω) and Gxf() can be well approximated by

Gf(ω)=(2πω)2,   Gxf(ω)=i2πω, \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${G^f}(\omega ) = {(2\pi \omega )^2},\;\;\;{G^{xf}}(\omega ) = i2\pi \omega ,$ \end{document}

Eq.(25), Eq.(26), Eq.(29) and Eq.(30) reduce, for 2πω ≪ 1, to

Γf(ω) (2πω)2 Γx(ω), \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\Gamma ^f}(\omega ) \simeq {\rm{ }}{(2\pi \omega )^2}\;{\Gamma ^x}(\omega ),$ \end{document} Γxf(ω)i2πω Γx(ω), \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\Gamma ^{xf}}(\omega ) \simeq i2\pi \omega \;{\Gamma ^x}(\omega ),$ \end{document} Axf(ω)2πωΓx(ω), \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${A^{xf}}(\omega ) \simeq 2\pi \omega {\Gamma ^x}(\omega ),$ \end{document} Pxf(ω) ±π/2. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${P^{xf}}(\omega ){\rm{ }} \simeq \pm \pi /2.$ \end{document}

The above spectral relations are known (H. von Storch & Zwiers, 1999; Jenkins & Watts, 1968; J.-S. von Storch, 1999b). They result from the first difference operator (·)s – (·)s–1 as a high-pass filter, characterized by the proportionality factor Gf(ω) (also known as the gain function). J.-S. von Storch (1999b) discussed the case when Γx(ω) is white at low frequencies. However, these relations are generally valid for solutions of any dynamical system that is governed by evolution equations in form of Eq.(1), independent whether the solution fluctuates stationary or not, and regardless how covariance functions of the solution behave with increasing time lag. In such a general case, Γx(ω) can have any spectral shape at the lowest frequencies, including having a power law Γx(ω) ∝ ω–n with n ≥ 2. Γf(0) would then be finite and non-zero if n = 2, and infinite if n > 2.

The situation is different when considering equilibrium solutions of a dynamical system of our concern specified in Section 1. In this case, all three covariance functions, γτx,γτf \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _\tau ^x,\gamma _\tau ^f$ \end{document} and γτxf \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\gamma _\tau ^{xf}$ \end{document} , have a finite cut-off time lag. Denote the largest of the three by τ. We can define a low-frequency interval Ω having a non-zero width (as it is done in Eq.(20)), such that for ω ϵ Ω, i.e. for ω ⋍ 0 but not strictly ω = 0, Γx(ω), Γf(ω), and Γxf(ω) are essentially white. Γx(ω) in Eq.(33)–Eq.(35) can then be replaced by a finite Γx(0) for ω ϵ Ω, noting that Γx(0) must be finite since the respective equilibrium variance is finite. Eq.(33)–Eq.(35) can then be replaced by

Γf(ω) (2πω)2 Γx(0)0, \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\Gamma ^f}(\omega ) \simeq {\rm{ }}{(2\pi \omega )^2}\;{\Gamma ^x}(0) \simeq 0,$ \end{document} Γxf(ω) 0, i2πω Γx(0)0, \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\Gamma ^{xf}}(\omega ) \simeq {\rm{ }}\left( {0,\;i2\pi \omega } \right)\;{\Gamma ^x}(0) \simeq 0,$ \end{document} Axf(ω) 2πωΓx(0)0. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${A^{xf}}(\omega ) \simeq {\rm{ }}2\pi \omega {\Gamma ^x}(0) \simeq 0.$ \end{document}

At frequency ω = 0, we have exactly

Γf(0)=0,  Γxf(0)=0,  Axf(0)=0. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\Gamma ^f}(0) = 0,\;\;{\Gamma ^{xf}}(0) = 0,\;\;{A^{xf}}(0) = 0.$ \end{document}

Thus, different from Eq.(33) which represents only a relation between Γx(ω) and Γf(ω) at the lowest non-zero frequencies, Eq.(37) specifies the spectral values of Γf(ω) at frequency ω ϵ Ω: f de facto ceases to fluctuate at frequency ω ϵ Ω due to the finite constant spectral level of Γx (and the smallness of ω), even though for any non-zero ω ϵ Ω, Γf(ω) is numerically not exactly zero and decreases with decreasing ω at the rate of ω2. Because f de facto ceases to fluctuate at frequency ω ϵ Ω, Γxf(ω) and Axf(ω) are de facto zero for ω ϵ Ω, even though for any non-zero ω ϵ Ω, Γxf(ω) is numerically not exactly zero and whose amplitude decreases with decreasing ω at the rate of ω, and the normalized amplitude spectrum, the coherence spectrum Cxf(ω), is one at all non-zero frequencies.

That Γf(ω) is de facto zero has a strong physical implication. As it is shown in Figure 7a for the Lorenz system, the spectrum of component x, Γx(ω), is non-zero at the lowest frequencies. Figure B2 further shows how Γx(ω) converges to a non-zero Γx(0), as ω approaches zero. With Γx(0) ≠ 0 and Ω having a non-zero width, the spectral plateau of Γx(ω) over the frequency range Ω makes up a non-zero portion of the total equilibrium variance of x. This portion of the total equilibrium variance of x, hereafter referred to as Γ(0)-variance of x, cannot be generated by fluctuations of the differential forcing f of x, since f de facto ceases to fluctuate at ω ϵ Ω.

Spectra of component x (top left), spectra of the differential forcing f of x (top right), coherence (bottom left) and phase (bottom right) spectra between x and f for the Lorenz components (green, red, and blue). All spectra are derived as the average of m = 1000 sample spectra or sample cross-spectra, each derived from an equilibrium solution of length 2N = 107. A sample cross-spectrum is derived from x and f obtained from the same equilibrium solution. Following Eq.(B5), the sample spectra of x and f are defined as (2N+1)x^ωx^ω* \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $(2N + 1){\hat x_\omega }\hat x_\omega ^*$ \end{document} and (2N+1)f^ωf^ω* \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $(2N + 1){\hat f_\omega }\hat f_\omega ^*$ \end{document} , and the sample cross-spectrum between x and f is defined as (2N+1)x^ωf^ω* \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $(2N + 1){\hat x_\omega }\hat f_\omega ^*$ \end{document} , where x^ω \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\hat x_\omega }$ \end{document} and f^ω \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\hat f_\omega }$ \end{document} are, respectively, the Fourier coefficients used to Fourier decompose finite records of x and f. The ensemble averaged cross-spectra is used to obtain the coherence and phase spectrum between x and f, following the usual definition of these spectra. The colored vertical lines, which are considered in Section 6, indicate the frequency ω* at which the spectrum of x starts to deviate more than 5% from its low-frequency plateau. ω* is determined from running averaging the spectra shown in a) to further reduce the remaining sample variations.

Spectra of Lorenz components x and their differential forcings

The above consideration applies only for solutions with infinite length. In practice where only finite records are available., complications arise. Figure 7 shows for the three Lorenz components, estimates of the spectra of x (top left), the spectra of differential forcing f (top right), and the coherence and phase spectra between x and f (bottom). All three Lorenz components reveal similar spectral behavior. The spectra of f have a maximum at frequencies around several 10–2 cycles per time step (cpts), and decrease thereafter with decreasing frequency, with the rate of decrease being close to ω2, consistent with Eq.(37). The decrease of spectrum of f is stopped at frequencies around 10–5–10–4 cpts. Thereafter and down to the smallest resolved frequency, the spectra of f are white. The coherence spectrum increases from zero at the highest frequencies, at which the spectral values of f are extremely small, to one at intermediate frequencies, consistent with Eq.(31). The phase at these intermediate frequencies is –90°, consistent with Eq.(30). At frequencies lower than 10–5–10–4 cpts, at which the spectra of f stop to decrease with decreasing frequency, the coherence and –90°-phase relation completely break down.

The low-frequency behaviors shown in Figure 7 are closely linked to the frailty of f. Appendix D shows that a spectrum estimated using finite records is subjected to truncation errors. Appendix D shows further that the spectrum of truncation error is white at the lowest frequencies. As the spectrum of f decreases with ω and becomes smaller than the truncation error at sufficiently low frequencies, the truncation error dominates. The spectrum of f approaches the white error spectrum. At these low frequencies, the coherence breaks down and the phase deviates strongly from -90°. For the Lorenz system, the frailty of f at the lowest frequencies can be more clearly seen by plotting the spectra in linear (rather than log) scale. Figure 8 shows Γf(ω) and Axf(ω) in linear scale. Since Γf(ω) at the lowest frequencies, at which f is weaker than the white truncation errors, is at least five to six orders of magnitude smaller than Γf(ω) at frequencies around 10–2 cycles per time step (cpts), Γf(ω) and Axf(ω) at the lowest frequencies appear to be zero relative to Γf(ω) and Axf(ω) at frequencies around 10–2 cpts.

Γf (same as in Figure 7b) and amplitude spectrum Axf(ω) between component x and its differential forcing f for the Lorenz components (green, red, blue), estimated as in Figure 7, but plotted in linear scale.

Spectral maxima of differential forcings

Note that the truncation error is not critical for the spectrum of x at the lowest frequencies. Relative to the high low-frequency spectral plateau of x, the truncation error can always be made negligibly small by considering a not too small N.

Before leaving this section, we point out that the spectral behaviors at the lowest frequencies found for the Lorenz system (Figure 7) resemble those found for the global axial angular momentum M and the global torque F simulated by a coupled GCM (Figure 1). In both cases, albeit that the spectra obtained from the GCM are much noisier, the behaviors of both models are characterized by a) an essentially white low-frequency plateau of spectrum of M or x, b) a tendency of decreasing spectral level of F or f with decreasing frequency before the decrease is stopped at the lowest frequencies by the truncation error, and c) the lack of coherence between M and F or between x and f at the lowest frequencies. The similarity suggests that these spectral features are general properties of a Newtonian type dynamical system subjected to similar truncation errors.

Generally, spectra are derived using limited data in practice. The resulting spectrum estimate is subjected both to sampling error and to truncation error. The former arises when applying the average Em to an ensemble of too small size. When no ensemble is available, one often uses an average over several chunks of a record, or a smoothing in frequency domain. For instance, the spectra in Figure 1 were obtained by cutting the available simulation into 5 pieces. The truncation error on the other hand emerges when estimating near-zero spectral values at the lowest frequencies using records of finite length. The spectra in Figure 1 are contaminated by both types of errors. The spectra shown in Figure 7, which are derived by using ensembles of size 1000, are essentially free of sampling error and mostly only contaminated by the truncation error.

We point out another side of the problem arising from using finite records, often encountered in climate applications. Generally, a prognostic variable describing a climate in an equilibrium state should reveal a low-frequency spectral plateau. Using records with limited length impedes fully resolving such spectral plateaus. For this reason, white low-frequency plateaus of prognostic climate variables (in an equilibrium climate state) are not always perceived.

Finally, we point out that the time derivative operator d()dt \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{{d( \cdot )} \over {dt}}}$ \end{document} filters out low-frequency spectral signals in the same way as the first difference (·)s – (·)s–1. When translating Eq.(1) into a spectral relation in frequency domain, we find that Γf is proportional to Γx in a way similar to that given in Eq.(25), but with the proportionality factor Gf(ω) = (2sin(πω))2 being replaced by (2πω)2, which approaches Gf(ω) = (2sin(πω))2 with decreasing ω. Thus, the result concerning the low-frequency parts of spectra obtained from Eq.(3) applies also to that obtained from Eq.(1). For this reason, d()dt \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{{d( \cdot )} \over {dt}}}$ \end{document} is not further considered.

5 Discussions

Fluctuations of a component x of an equilibrium solution at near-zero frequencies, that contribute to a non-zero Γ(0)-variance, cannot be determined by the fluctuations of its differential forcing f, since f de facto does not fluctuate at the lowest frequencies. But what else, if not f, generates the fluctuations of x at near-zero frequencies? Here we do not pretend to answer this question rigorously. Instead, we discuss a few ideas which we think are crucial for addressing this question. The overall outcome of the discussion – that the fluctuations of x at near-zero frequencies result from an integral effect – seems to be reasonable and self-evident, as these fluctuations definitely exist when the dynamical system of our concern is integrated forward in time.

We start by considering for a given dynamical system, the major frequencies (or the major timescales) of the underlying dynamics, defined as the frequencies (timescales), at which the spectra of f have their maxima. For the Lorenz model (Figures 7b and 8a), the most pronounced spectral peak is found for the differential forcing of the third component (blue line) located at the frequency near 1.42 × 10–2 cpts, which corresponds to a timescale of about 70 time steps. For the first and second component (green and red), the spectra of f are less peaky. There is a maximum located at frequencies slightly lower than 1.42 × 10–2 cpts and a maximum at frequencies around 2 × 10–2 cpts. We hence consider the frequencies of a few 10–2 cpts as the major frequencies of the Lorenz dynamics. The maxima of the spectra of x (Figure 7a) are located at similar frequencies as the spectra of f, although the spectra either further increase with decreasing frequency or drops only slightly at frequencies lower than the major frequencies. Consistent with these spectrum maxima, the time series of x and f (Figure 2) fluctuate visually only on the major timescales of the underlying dynamics. For x, whose spectrum has a white low-frequency plateau, the total low-frequency variance is evenly distributed over all low frequencies, making the fluctuations at an individual low frequency extremely weak and impossible to spot (even when sufficiently long records are used to resolve these low frequencies). Generally, we assume that the major timescale of a dynamical system of our concern is finite.

Given the finiteness of the major timescale of the underlying dynamics, it becomes even more striking that x fluctuates at zero and almost-zero frequencies. We approach this issue by considering a Brownian particle embedded in fluid molecules. The differential forcing of the Brownian particle describes the collisions of the particle with its surrounding fluid molecules. The major timescale of the underlying dynamics is the collision timescale. With each individual collision, the Brownian particle receives an impulse signal, which is a function of time that is only non-zero within a very short time interval and zero everywhere outside the interval. An impulse signal has an essentially white spectrum. Thus, a collision received by a Brownian particle makes the particle to fluctuate a little at almost all frequencies, especially the lowest frequencies.

The Brownian motion example can be extended to a dynamical system such as a GCM or the Lorenz model. A particular interaction of a component x with the other components, as described by the differential forcing f of x, produces a signal in x that is non-zero within a time interval, but zero outside the time interval. Since the length of the interval, which corresponds to the major timescale of the underlying dynamics, is always negligibly short relative to the infinite timescale, the signal can be considered as an impulse signal. Similar to a collision received by the Brownian particle, an interaction with other components received by x makes x to fluctuate a little at almost all frequencies, especially the lowest frequencies.

We emphasize that different from a delta function, which has a peak of zero duration, an impulse-like signal received by x has a “peak” lasting over a non-zero duration. Even for a collision of a Brownian particle with a molecule, which is completed within a time interval that is too short to be recognized by the human eye, the duration of the collision is, relative to the major timescale of the differential forcing that describes the collisions, not zero. Figure 9 shows typical spectra of impulse-like signals with a peak of non-zero duration. Different from the spectrum of a delta function, which is white and non-zero over the entire frequency domain, the spectrum of a signal of non-zero duration has a white low-frequency plateau. We emphasize also that x can only receive a full impulse-like signal, after x has completed the interaction with the other components, or in other words, after x has evolved along a solution over the duration of the interaction. Consequently, the low-frequency fluctuations arising from an impulse-like interaction received by x is not compatible with the low-frequency fluctuations of the differential forcing f, which is a function of every time instance during the interaction.

Spectra of two impulse-like signals of length 1 × 107 + 1. The signals are not zero at the central n = 11 (red) and n = 101 (blue) time steps, but zero at other 1 × 107 + 1 – n time steps. The values of the signal at the central n time steps are obtained from an AR1-process with the AR1 coefficient being 0.9. The spectra shown are averages over 100 sample spectra obtained from 100 similarly generated impulse-like signals. The spectra are white and non-zero over the low-frequency interval [0,ωs], with ωs (vertical lines) being dependent on n. The smaller the value of n, the larger the frequency ωs, the wider the interval [0,ωs]. As n → 1, the spectrum of an impulse-like signal approaches the spectrum of a delta function, which is white and equals the lag-zero auto-covariance (the only non-zero value of the auto-covariance function) at all frequencies according to Eq.(14). This n-dependent spectral behavior is independent of the details of the AR1-process used to produce the non-zero values of the impulse-like signal at the central n time steps. ωs is identified as the frequency at which the spectrum deviates 1% from its value at the lowest resolved frequency.

Spectra of impulse-like signals

Along an equilibrium solution, component x receives not only just one interaction with other components of x, but many interactions. While each impulse-like interaction received by x makes x to fluctuate a little at the lowest frequencies, many impulse-like interactions can produce the non-zero Γ(0)-variance of x. We refer the net effect of many interactions of x with other components of x along a solution as the integral effect of all other components on component x. The net integral effect cannot be represented as a differential forcing defined at individual time instances, but can only be obtained by integrating f forward in time.

Still unclear is, as an equilibrium solution continues into an infinite time, what stops the integral effect to include more impulse-like interactions, thereby keeping the Γ(0)-variance and the respective total equilibrium variance finite. Purely integrating (i.e. simply adding) interactions of x with other components along a solution resembles a random walk that integrates a white noise over time. Such a pure integration makes the variance to increase with increasing time, and making an equilibrium state unsustainable. Obviously, there must be a dissipation, or a negative feedback as suggested by Hasselmann (1976), that counteracts the pure integration of interactions of x with other components.

Hasselmann identified the required dissipation as a part of the underlying dynamics, i. e. as a part of the differential forcing f. Such a dissipation dissipates variations of x at every time instance. We know from the example of Brownian motion that there is another type of dissipation: As a Brownian particle moves among fluid molecules, it experiences a friction that results from numerous collisions. This friction does not exist in the governing equations formulated in terms of differential forcing, since within a time interval δt → 0 that defines the differential operator ddt() \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{d \over {dt}}}( \cdot )$ \end{document} , a Brownian particle receives only isolating collisions. It appears that for a system in its equilibrium state reached under a constant external forcing, there exist two types of dissipation: One appears as a part of differential forcing and counteracts the constant external forcing at every time instance, which is needed to prevent the system to run away under the constant external forcing. The other appears as a part of the integral effect that does not exist in δt → 0 and counteracts the continuous integration of impulse-like interactions over time. As the impulse-like interactions make x to fluctuate at low frequencies, the latter type of dissipation is linked to (low-frequency) fluctuations, in line with the dissipation-fluctuation theorem (DFT) (Risken, 1989, chapter 7). DFT states that the response of a system to a change in its external forcing, which depends sensitively on the dissipation processes inside the system, is closely related to the natural fluctuations occurring without the change in the external forcing. Clearly, this dissipation issue and its link to DFT need further investigations.

Finally, we point out that the differential forcing f considered in this paper differs from a mechanical force fmech considered in classical mechanics in two aspects. First, when exerting a mechanical force fmech on an object, the force is equally received by every part of the object. For a system considered in this paper, although all components of x are subjected to the same external forcing, the differential forcings of any two different components are generally different. The external forcing is not equally received by all components of x. Secondly, a mechanical force fmech exerted on an object can be long-lasting, whereas a differential forcing considered in this paper is always short living. In the simplest case in which fmech is constant, the momentum of the object would be constantly accelerated, and the variance of the momentum would be infinite. In terms of spectral relation given in Eq.(33), which holds both for solutions of a system considered in this paper and for the evolution of the momentum caused by a mechanical force fmech, the spectrum Γfmech(ω) \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\Gamma ^{{f_{mech}}}}(\omega )$ \end{document} of fmech would be a line spectrum with a non-zero value at ω = 0, and the spectrum of the momentum of the object would be proportional to Γfmech(ω)(2πω)2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\Gamma ^{{f_{mech}}}}(\omega ){(2\pi \omega )^{ - 2}}$ \end{document} and hence infinite at ω = 0. Different from that, a differential forcing considered in this paper does not have any energy (variance) at frequency zero, and de facto does not have any energy at all near-zero frequencies (see Figure 8a). Such a forcing is ephemeral. The ephemerality of f is reflected in decaying correlation functions, a feature not existing in classical mechanics.

The ephemerality is linked to the multidimensionality of the dynamical system considered in this paper. For a multidimensional system, even though the form of f(x) as a function of x does not change over time, the contributions of different components of x to f(x) generally change over time. In case of a Brownian motion, the contributions to f(x) at different times stem from collisions of the Brownian particle with different molecules. Because of this constant change in contributors to f(x), f cannot be long-lasting.

An ephemeral differential forcing f of x generates fluctuations in two different ways: one via f at each time instance; the other via the ephemerality of f related to the decaying and reawakening of f over time – due to the constant external forcing and the dissipation contained in f. The former makes x to fluctuates at the frequency of f, in the same way as in classical mechanics that the momentum variations at a frequency result from the variations of the mechanical force at the same frequency. The latter, which can only be realized by integrating the dynamical system forward over time, makes x to fluctuate at the lowest frequencies via impulse-like interactions of x with other components of x.

6 Concluding remarks

This paper considers a dynamical system described by a multidimensional state vector x such that the temporal evolution of a component x of x is governed by dxdt=f(x). \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{{dx} \over {dt}}} = f({\rm{x}}).$ \end{document} . Under a constant external forcing, the system reaches an equilibrium state characterized by equilibrium fluctuations that have well-defined finite equilibrium variances and decaying covariance functions. We show that an equilibrium variance can only be decomposed in frequency domain when taking all equilibrium solutions of the considered equilibrium state into account. We show further that for a system of our concern, the differential operator d()dt \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{{d( \cdot )} \over {dt}}}$ \end{document} , together with the white low-frequency spectral plateau of x, constrain the differential forcing f of x such that f de facto does not fluctuate at near-zero frequencies. The variance related to the fluctuations of x at near-zero frequencies can hence not be generated by fluctuations of f. Instead, it results from an integral effect of many impulse-like interactions of x with other components of x occurring along an equilibrium solution over time. This integral effect cannot be represented as a differential forcing and only emerges when integrating the system forward in time. The forward integration sets the arrow of time toward an equilibrium state.

The significance of the present study depends on the amount of the total equilibrium variance that results from the integral effect. We assume that this variance equals the Γ(0)-variance, i.e. the variance that is associated with the white low-frequency spectral plateau of x. For the Lorenz system, we estimate the Γ(0)-variance in terms of ν defined as

ν=0ω*Γx(ω)dω01/2Γx(ω)dω×100%. \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} $\nu = \frac{{\int_0^{{\omega _*}} {{\Gamma ^x}} (\omega )d\omega }}{{\int_0^{1/2} {{\Gamma ^x}} (\omega )d\omega }} \times 100\% .$ \end{document}

ω* is the frequency that for ω > ω* the spectrum of component x, Γx(ω), deviates from its white low-frequency plateau. For the three Lorenz components, ω* is indicated by the colored vertical lines in Figure 7a). ν amounts to 13.86%, 11.48%, 3.79% of the total variance, respectively, suggesting that only a small portion of the equilibrium variance is determined by the integral effect. We speculate that given a component x of a dynamical system, the amount of Γ(0)-variance depends on the “inertia” of the component. A more inert component needs to interact with other components a greater number of times before it starts to fluctuate. The integral effect that results from a larger number of interactions of x with other components leads to a higher white low-frequency plateau of x, and with that a larger Γ(0)-variance of x. It is conceivable that the Γ(0)-variance of a Brownian particle is essentially identical to its entire equilibrium variance.

This paper has strong implication for inert components of a dynamical system. The equilibrium variance of an inert component, which is largely determined by the integral effect, can no longer be deciphered by considering processes described by the differential forcing at individual time instances. If we are not content with estimating the total equilibrium variance by integrating the governing equations in form of dxdt=f(x). \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{{dx} \over {dt}}} = f({\bf{x}}).$ \end{document} forward in time, but wish to cast the integral effect in an evolution equation, we need to find a new framework that quantifies and represents the integral effect. Such a framework can take the form of a stochastic difference equation, which represents the role of the integral effect in generating variance in terms of a stochastic forcing, and at the same time, contains a dissipation counteracting the continuous generation of variance by the stochastic forcing. Here I deliberate use the word “difference”, as oppose to “differential”, to emphasize that the integral effect is an effect not expressible in term of the differential operator ddt() \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{d \over {dt}}}( \cdot )$ \end{document} . Within an infinitesimal time interval associated with the operator ddt() \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} ${\textstyle{d \over {dt}}}( \cdot )$ \end{document} , the integral effect does not exist; we have only the differential forcing at every time instance. The integral effect is the physical “thing” brought up in Section 1 that escapes the underlying deterministic dynamics represented by the differential forcing f. In this sense, a framework that quantifies the integral effect based on stochastic difference equations is not an approximation of the underlying deterministic differential equation, but represents processes not represented by the differential forcing. The subject of such a framework are fluctuations at near-zero frequencies only, rather than fluctuations at the major frequencies of the underlying dynamics.

The concept of stochastic climate models proposed by Hasselmann (1976) is an attempt toward formalizing such a framework. If such a framework is established for the dynamical system of our concern (e.g. for the climate), we would have for any one component x of the state vector x two types of evolution equations: One describes the detailed interaction of x with other components occurring at every time instance in terms of the differential forcing f(x) of x. The other describes the generation and the dissipation of the portion of the equilibrium variance of x that arises from the integral effect occurring along a solution over time. The total equilibrium variance of x can be obtained by integrating the first type of evolution equations, but cannot be fully explained by the dynamics represented by only one of the two types of evolution equations.

Data accessibility statements

All data and subsequently the figures are generated by Matlab-scripts archived in the publication repository of the Max Planck Society under http://hdl.handle.net/21.11116/0000-0009-9904-6.

The additional file for this article can be found as follows:

Appendix A to D. DOI: https://doi.org/10.16993/tellusa.25.s1

Acknowledgement

I thank Christian Reick, Cathy Hohenegger, Eduardo Zorita, Hans von Storch and Jochem Marotzke for helpful discussions. My thanks go also to the two anonymous reviewers for their constructive comments.

Competing Interests

The authors have no competing interests to declare.

Balaji, V, Maisonnave, E, Zadeh, N, Lawrence, BN, Biercamp, J, Fladrich, U, Wright, G, et al. 2017. Cpmip: Measurements of real computational performance of earth system models in cmip6. Geoscientific Model Development, 10(1): 1934. Retrieved from https://gmd.copernicus.org/articles/10/19/2017/. DOI: 10.5194/gmd-10-19-2017 Bartlett, M. 1946. On the theoretical specification of sampling properties of auto correlated time series. J.R. Stat. Soc., B8: 2741. DOI: 10.2307/2983611 Bloomfield, P. 1976. Fourier Analysis of Time Series: An Introduction. Wiley. Dijkstra, HA and Ghil, M. 2005. Low-frequency variability of the large-scale ocean circulation: A dynamical systems approach. Reviews of Geophysics, 43: RG3002. DOI: 10.1029/2002RG000122 Eckmann, J and Ruelle, D. 1985. Ergodic theory of chaos and strange attractors. Reviews of Modern Physics, 57(3): 617656. DOI: 10.1103/RevModPhys.57.617 Gardiner, C. 2004. Handbook of Stochastic Methods for physics, chemistry, and natural sciences. Springer. DOI: 10.1007/978-3-662-05389-8 Hasselmann, K. 1976. Stochastic climate models Part I. Theory. Tellus, 28(6): 473485. DOI: 10.1111/j.2153-3490.1976.tb00696.x James, S and James, PM. 1989. Ultra-low-frequency variability in a simple atmo spheric circulation model. Nature, 342: 5355. DOI: 10.1038/342053a0 Jenkins, GM and Watts, DG. 1968. Spectral analysis and its applications. Holden-Day. Kerr, RA. 2000. A North Atlantic climate pacemaker for the centuries. Science, 288: 19841985. DOI: 10.1126/science.288.5473.1984 Lasota, A and Mackey, M. 1994. Chaos, fractals, and noise: Stochastic aspects of dynamics. In Applied mathematical sciences, 97: 472. 2nd ed. New York: Springer. DOI: 10.1007/978-1-4612-4286-4 Lorenz, E. 1963. Deterministic nonperiodic flow. Journal of Atmospheric Sciences, 20(2): 130141. DOI: 10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2 Mitchell, J. 1976. An overview of climatic variability and its causal mechanisms. Quaternary Research, 6(4): 481493. Retrieved from https://www.sciencedirect.com/science/article/abs/pii/0033589476900211. DOI: 10.1016/0033-5894(76)90021-1 Newman, M, Alexander, MA, Ault, TR, Cobb, KM, Deser, C, Lorenzo, ED, Smith, CA, et al. 2016. The pacific decadal oscillation, revisited. Journal of Climate, 29(12): 43994427. Retrieved from https://journals.ametsoc.org/view/journals/clim/29/12/jcli-d-15-0508.1.xml. DOI: 10.1175/JCLI-D-15-0508.1 Panagiotopoulos, F, Shahgedanova, M and Stephenson, D. 2002. A review of Northern Hemisphere winter-time teleconnection patterns. Journal de Physics IV, 12: 2747. DOI: 10.1051/jp4:20020450 Priestley, M. 1981. Spectral analysis and time series. Academic Press. Rasmusson, E and Carpenter, T. 1982. Variations in tropical sea surface temperature and surface wind fields associated with the Southern Oscillation- El Niño. Monthly Weather Review, 110: 354384. DOI: 10.1175/1520-0493(1982)110<0354:VITSST>2.0.CO;2 Risken, H. 1989. The Fokker-Planck equation. In (p. 472). Springer. DOI: 10.1007/978-3-642-61544-3 von Storch, H and Zwiers, FW. 1999. Statistical analysis in climate research. Cambridge University Press. DOI: 10.1007/978-3-662-03744-7_2 von Storch, J-S. 1999a. The reddest atmospheric modes and the forcings of the spectra of these Modes. Journal of the Atmospheric Sciences, 16141626. DOI: 10.1175/1520-0469(1999)056<1614:TRAMAT>2.0.CO;2 von Storch, J-S. 1999b. What determines the spectrum of a climate variable at zero frequency? Journal of Climate, 21242127. DOI: 10.1175/1520-0442(1999)012<2124:WDTSOA>2.0.CO;2