Quantifying nonergodicity in nonautonomous dissipative dynamical systems: An application to climate change

In nonautonomous dynamical systems, like in climate dynamics, an ensemble of trajectories initiated in the remote past defines a unique probability distribution, the natural measure of a snapshot attractor, for any instant of time, but this distribution typically changes in time. In cases with an aperiodic driving, temporal averages taken along a single trajectory would differ from the corresponding ensemble averages even in the infinite-time limit: ergodicity does not hold. It is worth considering this difference, which we call the nonergodic mismatch, by taking time windows of finite length for temporal averaging. We point out that the probability distribution of the nonergodic mismatch is qualitatively different in ergodic and nonergodic cases: its average is zero and typically nonzero, respectively. A main conclusion is that the difference of the average from zero, which we call the bias, is a useful measure of nonergodicity, for any window length. In contrast, the standard deviation of the nonergodic mismatch, which characterizes the spread between different realizations, exhibits a power-law decrease with increasing window length in both ergodic and nonergodic cases, and this implies that temporal and ensemble averages differ in dynamical systems with finite window lengths. It is the average modulus of the nonergodic mismatch, which we call the ergodicity deficit, that represents the expected deviation from fulfilling the equality of temporal and ensemble averages. As an important finding, we demonstrate that the ergodicity deficit cannot be reduced arbitrarily in nonergodic systems. We illustrate via a conceptual climate model that the nonergodic framework may be useful in Earth system dynamics, within which we propose the measure of nonergodicity, i.e., the bias, as an order-parameter-like quantifier of climate change.

Figure 1

Nonergodic mismatch pdfs $P({\delta }_{\tau }\left(t\right))$ based on the variable $\phi =y$, calculated over the numerical ensemble of $N={10}^6$ trajectories described in Appendix pp1. In each panel we compare three different values of the window length $\tau$. The time instants $t$ of observation are indicated above the panels and are taken from a stationary and a changing climate. In the insets we show the time dependence (3) of the parameter $F_0$ of (2), mark the time instants $t$ of observation by vertical gray dot-dashed lines, and indicate the different time windows (of length $\tau$) in the same color (shade of gray in print) as in the main plot. The ensemble size appears to realize the asymptotic limit in the sense that we do not find any considerable change in the graphs when plotting the results for $N=10000$, $20000$, or ${10}^6$. For the histogram, the bin size is 0.025.

Figure 2

The bias ${\mathcal{A}}_{\mu }\left({\delta }_{\tau }\right)$ as a “continuous” function of the window length $\tau$, calculated for the numerical ensemble of ${10}^6$ trajectories, where the nonergodic mismatch ${\delta }_{\tau }$ is based on the variable $y$, as in Fig. 1. The time instants $t$ of observation are year $-75$ (flat black line, stationary climate) and year 76 [magenta line (gray in print), changing climate].

Figure 3

The spread ${\sigma }_{\mu }\left({\delta }_{\tau }\right)$ as a function of the number of years, $\tau +1$, included in the temporal average, plotted on a doubly logarithmic scale, calculated for the numerical ensemble of ${10}^6$ trajectories, where the nonergodic mismatch ${\delta }_{\tau }$ is based on the variable $y$, as in Fig. 1. The time instants $t$ of observation are year $-75$ (black line, stationary climate) and year 76 [light blue line (gray in print), changing climate]. The dashed line is of slope $-1/2$ to guide the eye.

Figure 4

The ergodicity deficit ${\mathcal{A}}_{\mu }\left(\right|{\delta }_{\tau }\left|\right)$, compared to the bias ${\mathcal{A}}_{\mu }\left({\delta }_{\tau }\right)$ and the spread ${\sigma }_{\mu }\left({\delta }_{\tau }\right)$, as a function of the window length $\tau$, calculated for the numerical ensemble of ${10}^6$ trajectories, where ${\delta }_{\tau }$ is based on the variable $y$, as in Fig. 1. The time instants $t$ of observation are indicated in the panels and are taken from a stationary and a changing climate.

Figure 5

(a) Histogram over $y$ (coded by the brightness), and the ensemble average ${\mathcal{A}}_{\mu }\left(y\right)$ [red (gray in print) “continuous” line] as a function of time, taken at $t\mathrm{mod}T=T/4$ time instants, calculated over a numerical ensemble of ${10}^6$ trajectories. The bin size of the histogram is $0.025$. (b) Intersections of the snapshot attractor with the $(x,y)$ (i.e., the $z=0$) plane conditioned by $\dot{z}>0$ in the particular years indicated in the particular plots. These years are marked in panel (a) by black circles.

Figure 6

The bias ${\mathcal{A}}_{\mu }\left({\delta }_{\tau }^{\left(2\right)}\right)$, the spread ${\sigma }_{\mu }\left({\delta }_{\tau }^{\left(2\right)}\right)$ and the ergodicity deficit ${\mathcal{A}}_{\mu }\left(\right|{\delta }_{\tau }^{\left(2\right)}\left|\right)$ from the nonergodic mismatch ${\delta }_{\tau }^{\left(2\right)}$ of the second cumulant based on the variable $x$ of the model (1, 2, 3), as functions of the window length $\tau$. The time instants $t$ of observation are indicated in the panels and are taken from a changing climate. Note that ${\mathcal{A}}_{\mu }\left({\delta }_{\tau }^{\left(2\right)}\right)$ is not identically zero either in year 76, which indicates nonergodicity for this year, too. The numerical ensemble, consisting of ${10}^6$ trajectories, is the one described in Appendix pp1.

Figure 7

The time dependence of the histograms (a) $P({\delta }_{\tau }\left(t\right))$ of the variable $y$ and (b) $P({\delta }_{\tau }^{\left(2\right)}\left(t\right))$ of the variable $x$, with a fixed window length of $\tau =72$ years, calculated over the numerical ensemble of $N={10}^6$ trajectories. In any single panel we compare different time instants $t$ of observation. The bin size is 0.025. Note that the lines for the two first time instants almost overlap.

Figure 8

The time dependence of a few characteristic statistical measures of the natural measure as indicated in the legend (derived from the numerical ensemble of ${10}^6$ trajectories). The blue line (the darkest gray in print) of the upper plot $[{\mathcal{A}}_{\mu }({\delta }_{\tau }\left(t\right))$ and ${\mathcal{A}}_{\mu }({\delta }_{\tau }^{\left(2\right)}\left(t\right))$ in panels (a) and (b), respectively] is the difference of the magenta (gray in print) and the black lines (the two lines of the lower plot) in view of (B1). Panel (a) concerns the average in the variable $y$, and panel (b) concerns the second cumulant in the variable $x$. The window length is $\tau =72$ yr. The time instants considered in Fig. 7 are marked by vertical gray dot-dashed lines. Observe that ${\sigma }_{\mu }({\delta }_{\tau }^{\left(n\right)}\left(t\right))$ and ${\mathcal{A}}_{\mu }(\left|{\delta }_{\tau }^{\left(n\right)}\left(t\right)\right|)$, $n=1,2$, lie very close to each other for $t<0$.

Figure 9

Results concerning the nonergodic mismatch ${\delta }_{\tau }$ based on the variable $y$ in the Lorenz 84 model driven by a double ramp, calculated for a numerical ensemble of $5\times {10}^5$ trajectories. (a)–(c) The ergodicity deficit ${\mathcal{A}}_{\mu }\left(\right|{\delta }_{\tau }\left|\right)$, compared to the magnitude $|{\mathcal{A}}_{\mu }\left({\delta }_{\tau }\right)|$ of the bias and to the spread ${\sigma }_{\mu }\left({\delta }_{\tau }\right)$, as a function of the window length $\tau$. The time instants $t$ of observation are indicated in the panels. (d) The time dependence of a few characteristic statistical measures as indicated in the legend. The window length is $\tau =72$ yr. The driving function $F_0\left(t\right)$ is indicated in the bottom plot, in orange (gray in print). The time instants considered in panels (a)–(c) are marked by vertical gray dot-dashed lines.

Figure 10

Same as Fig. 9 for the Lorenz 84 model driven by the quasiperiodic signal (D1), calculated for a numerical ensemble of $5\times {10}^5$ trajectories.

Figure 11

Same as Fig. 9 for the Hénon map (D2, D3) driven by a double ramp, calculated for a numerical ensemble of $5\times {10}^5$ trajectories.

Figure 12

Same as Fig. 9 for the Hénon map (D2, D3) driven by the quasiperiodic signal (D4), calculated for a numerical ensemble of $5\times {10}^5$ trajectories.

Figure 13

The pdf $P(d_{\tau }\left(t\right))$ of the nonergodic mismatch $d_{\tau }\left(t\right)$ (E4) based on the variable $\phi =y$ of (1)–(3), calculated over a numerical ensemble of $N=10000$ trajectories. In each panel we compare three different values of the window length $\tau$. The time instants $t$ of observation are as in Fig. 1. The bin size is 0.025.

Figure 14

The bias ${\mathcal{A}}_{\mu }\left(d_{\tau }\right)$, the spread ${\sigma }_{\mu }\left(d_{\tau }\right)$, and the nonergodic mismatch ${\mathcal{A}}_{\mu }\left(\right|d_{\tau }\left|\right)$ as “continuous” functions of the window length $\tau$, calculated for a numerical ensemble of $10000$ trajectories, where the nonergodic mismatch $d_{\tau }$ is based on the variable $y$, as in Fig. 13. The time instants $t$ of observation are indicated in the panels.