π-theorem generalization of the ice-age theory

Analyzing a dynamical system describing the global climate variations requires, in principle, exploring a large space spanned by the numerous parameters involved in this model. Dimensional analysis is traditionally employed to deal with equations governing physical phenomena to reduce the number of parameters to be explored, but it does not work well with dynamical ice-age models, because, as a rule, the number of parameters in such systems is much larger than the number of independent dimensions. Physical reasoning may, however, allow us to reduce the number of effective parameters and apply dimensional analysis in a way that is insightful. We show this with a specific ice-age model (Verbitsky et al., 2018), which is a low-order dynamical system based on ice-flow physics coupled with a linear climate feedback. In this model, the ratio of positive-to-negative feedback is effectively captured by a dimensionless number called the “V number”, which aggregates several parameters and, hence, reduces the number of governing parameters. This allows us to apply the central theorem of the dimensional analysis, the π theorem, efficiently. Specifically, we show that the relationship between the amplitude and duration of glacial cycles is governed by a property of scale invariance that does not depend on the physical nature of the underlying positive and negative feedbacks incorporated by the system. This specific example suggests a broader idea; that is, the scale invariance can be deduced as a general property of ice age dynamics if the latter are effectively governed by a single ratio between positive and negative feedbacks.

1 Introduction

Mathematical modeling of Pleistocene ice ages using astronomically forced spatially resolving models of continental ice sheets, the ocean, and the atmosphere has always been and remains a computational challenge. Therefore, though higher-resolution models (e.g., Abe-Ouchi et al., 2013) and models of intermediate complexity (e.g., Verbitsky and Chalikov, 1986; Chalikov and Verbitsky, 1990; Gallée et al., 1991; Ganopolski et al., 2010) are gaining popularity, it has been argued for a long time that significantly less computationally demanding dynamical models may provide just as much insight as the models with more degrees of freedom (Saltzman, 1990). However, even though the computational load for solving dynamical equations is minimal, the work and number of experiments needed for spanning the full parameter space is easily overwhelming. Analyzing a dynamical system of ice ages is thus, in principle, a difficult task. In mathematical physics, the method of dimensional analysis (e.g., Barenblatt, 2003) has been traditionally employed to take advantage of symmetry or invariance principles and, as a result, to reduce the number of effective parameters. It has not been applied to low-order models of the Pleistocene climate, because in such models the number of governing parameters is much larger than the number of independent dimensions. Indeed, the number of independent dimensions in a dynamical system does not exceed the number of variables (it may be smaller if some variables have the same or dependent dimensions), to which one adds time, which is always present in a dynamical system. For example, the dynamical system of Saltzman and Verbitsky (1993) described the evolution of four variables: ice volume (in cubic meters), CO₂ concentration (in parts per million), ocean temperature (in degrees Celsius), and bedrock depression (in meters). The number of independent dimensions, including time, was thus four. This system had 18 parameters, including the amplitude and the period of the external forcing. In such a case, the π theorem (Buckingham, 1914) – the tenet of dimensional analysis – is of little help in simplifying the analysis and effectively providing physical insight, because, even in the dimensionless form, the system would still contain 14 (18 − 4) dimensionless groups.

In Verbitsky et al. (2018), we derived a dynamical model of the Pleistocene climate from the scaled conservation equations of viscous non-Newtonian ice and combined them with an equation describing the evolution of the climate temperature. The work was motivated by the prospect of delivering a low-order, parsimonious approach to the problem of understanding glacial–interglacial cycles. The state of the ice–climate system is summarized by a 3-dimensional vector: glaciation area S (in square meters), ice sheet basal temperature θ (in degrees Celsius), and climate temperature ω (in degrees Celsius). The number of independent dimensions, including time, is thus three. However, despite our effort to be parsimonious in the physical description, the model includes 12 parameters, which is still much larger than the number of independent dimensions. As now we may have nine (12 − 3) dimensionless groups; this is obvious progress relative to the Saltzman and Verbitsky (1993) model but not enough for an effective use of the π theorem. The situation changed dramatically when we discovered that the dynamical properties of the system are largely defined by the dimensionless V number incorporating eight model parameters and measuring the ratio of climate system positive feedback over the negative feedback from the ice sheet itself. At once, seven parameters are effectively eliminated, and using the π theorem became an attractive prospect. We first applied the π theorem reasoning to investigate the propagation of millennial forcing into ice-age dynamics (Verbitsky et al., 2019a) and found that the millennial forcing introduces a disruption, i.e., shifts the system equilibrium point, and this disruption is proportional to the second degree of the forcing period.

In this paper we will apply this approach systematically to all model variables. This will allow us to demonstrate that, in the model, glacial area and climate temperature are scale invariant in the orbital frequencies domain (in the case of the climate temperature – even beyond this domain) and observe that this property does not depend on the specific physical nature of the climate system feedbacks. This observation is important. The empirical analysis of paleoclimate series shows that there is a rich spectral content and points to the existence of “spectral slopes” (e.g., Huybers and Curry, 2006; Lovejoy and Schertzer, 2013). Lovejoy and Schertzer (2013) evoke some generic process, such as the principle of “cascades” which is tightly linked to the concept of scale invariance of the equations. For example, the scale invariance of fluid-dynamics equations is exploited to provide inferences about spectral slopes of turbulent flows. However, to our knowledge, there is no available theory supporting scale invariance in regimes associated with glacial–interglacial dynamics. Yet, paleoclimate simulations with more sophisticated models, including the seminal paper by Abe-Ouchi et al. (2013) and the simulations with CLIMBER provided by Ganopolski et al. (2010), tend to focus on the response of the ice-sheet climate system to orbital forcing and discuss the respective amplitudes of the 100, 41, and 21–23 kyr periods, but none discuss the slope of the power spectrum down to the millennium scale. Therefore, we believe that our research will provide at least some important elements that should help us to bridge both approaches

Accordingly, our paper is structured as follows. First, we will briefly recapture equations, parameters, and dimensions of the Verbitsky et al. (2018) model. Then we will remind readers of the essence of the π theorem, apply it to all model variables, and discuss its implications.

2 A dynamical model of Pleistocene glacial rhythmicity

The nonlinear dynamical model of the global climate system (Verbitsky et al., 2018) is derived from the scaled equations of ice sheet thermodynamics, combined with a linear feedback equation involving an effective “temperature”, which describes the climate state outside the ice region.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{4}}{\mathrm{5}}}{\mathit{\zeta }}^{-\mathrm{1}}{S}^{\mathrm{3}/\mathrm{4}}\left(a-\mathit{\epsilon }{F}_{\mathrm{S}}-\mathit{\kappa }\mathit{\omega }-c\mathit{\theta }\right)\text{(2)}& {\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}t}}={\mathit{\zeta }}^{-\mathrm{1}}{S}^{-\mathrm{1}/\mathrm{4}}\left(a-\mathit{\epsilon }{F}_{\mathrm{S}}-\mathit{\kappa }\mathit{\omega }\right)\left\{\mathit{\alpha }\mathit{\omega }+\mathit{\beta }\left[S-S_{\mathrm{0}}\right]-\mathit{\theta }\right\}\text{(3)}& {\displaystyle \frac{\mathrm{d}\mathit{\omega }}{\mathrm{d}t}}={\mathit{\gamma }}_{\mathrm{1}}-{\mathit{\gamma }}_{\mathrm{2}}\left[S-S_{\mathrm{0}}\right]-{\mathit{\gamma }}_{\mathrm{3}}\mathit{\omega }\end{array}$$

The model variables and their dimensions are defined as follows: S (in square meters) is the glaciation area, θ (in degrees Celsius) is the basal ice sheet temperature, and ω (in degrees Celsius) is the effective global climate temperature. The third equation implicitly accounts for the effect of the response of CO₂ concentration, along with other radiative feedbacks.

Model parameters along with their dimensions are as follows: ζ (measured in meters to the 1/2 power) is the “shape” factor of the ice sheet; a (in meters per second) is the characteristic rate of snow precipitation; F_S is normalized mid-July insolation at 65^∘ N (Berger and Loutre, 1991); ε (in meters per second) is the amplitude of the external forcing; κ (in meters per second per degree Celsius) and c (in meters per second per degree Celsius) are sensitivity parameters, describing, correspondingly, climate temperature and basal sliding impacts into ice-sheet mass balance; the dimensionless coefficient α describes basal temperature sensitivity to global climate temperature changes; coefficient β (in degrees Celsius per square meter) defines basal temperature dependence on ice sheet dimensions; S₀ (in square meters) is a reference glaciation area; and γ₁ (in degrees Celsius per second), γ₂ (in degrees Celsius per square meter per second), and γ₃ (in per second) define climate temperature evolution, $\mathrm{1}/{\mathit{\gamma }}_{\mathrm{3}},$ being a time constant. If the forcing is periodic, then we may consider that the system dynamics are described by an additional parameter: the forcing period T (in seconds). Thus we have a system of three variables, three (including time) independent dimensions, and 12 parameters. The system Eqs. (1)–(3) is not sensitive to initial conditions, and, therefore, we do not include the latter into the list of parameters.

Physical reasoning and numerical experiments (Verbitsky et al., 2018) led us to the suggestion that the system response is essentially determined by the V number, measuring a balance between positive and negative model feedbacks as follows:

$$\begin{array}{}\text{(4)}& V={\displaystyle \frac{\mathrm{1}}{\mathit{\beta }}}\left(\mathit{\alpha }+{\displaystyle \frac{\mathit{\kappa }}{c}}\right)\left({\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{3}}}}-{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{1}}}{S_{\mathrm{0}}{\mathit{\gamma }}_{\mathrm{3}}}}\right).\end{array}$$

Here parameter β is a measure of ice-sheet negative feedback. The term $(\mathit{\alpha }+\mathit{\kappa }/c)({\mathit{\gamma }}_{\mathrm{2}}/{\mathit{\gamma }}_{\mathrm{3}}-{\mathit{\gamma }}_{\mathrm{1}}/{\mathit{\gamma }}_{\mathrm{3}}/S_{\mathrm{0}})$ measures the climate system positive feedback (Verbitsky et al., 2018).

If we assume that the V number effectively captures the behavior of the model with respect to the eight parameters included in its definition, then the number of parameters is effectively reduced to five: V, ζ, a, ε, and T. We assume further that parameter ζ in Eqs. (1)–(2) is a constant, thus assuming an invariant relationship between ice thickness H and glaciation area S ($H=\mathit{\zeta }{S}^{\mathrm{1}/\mathrm{4}}$; Verbitsky et al., 2018). We also note that the V number has been assembled using components of the steady-state solution of the system Eqs. (1)–(3) (Verbitsky et al., 2018). Obviously, parameter ζ, as a multiplier, is not part of this steady-state solution. Therefore our hypothesis that the V number defines the behavior of the model in fact also includes the assumption that the impact of the parameter ζ on the system behavior, at the reference value, is weak. As a result, we end up with the assumption that the response of the system to external forcing is essentially determined by no more than four parameters: V, a, ε, and T. We will now learn how to profit from this advantage.

3 Dimensional analysis of model variables

3.1 Period of the system response to the external forcing, P

We previously noticed (Verbitsky et al., 2018) that with weak climate positive feedback (V∼0), the system, exhibits fluctuations in response to the astronomical forcing with a dominating period of about 40 kyr, which may arise as either a direct response to obliquity or a doubled-period response to the forcing associated with climatic precession (2×20 kyr). When the climate positive feedback intensifies such that V∼0.75 and external forcing is strong, the system evolves with a doubled obliquity period. We can therefore assume that the period of the system response to the external forcing, P, is a function of the V number, the amplitude of the external forcing, ε, and the period of the external forcing, T. We thus begin with the most general hypothesis:

$$\begin{array}{}\text{(5)}& P=\mathit{\psi }\left(V,a,\mathit{\epsilon },T\right).\end{array}$$

It is at this stage that the π theorem intervenes. Specifically, it stipulates that a physical relationship should not depend on a system of units, and therefore, in the dimensionless form, the number of dimensionless arguments is equal to the total number of the governing parameters minus the number of governing parameters with independent dimensions (Buckingham, 1914). If we select dimensions of ε and T as independent dimensions, then application of the π theorem to the Eq. (5) gives us the following:

$$\begin{array}{}\text{(6)}& {\displaystyle }P/T=\mathrm{\Psi }\left(V,\mathit{\epsilon }/a\right)\text{(7)}& {\displaystyle }P=T\mathrm{\Psi }\left({\mathrm{\Pi }}_{\mathrm{1}},{\mathrm{\Pi }}_{\mathrm{2}}\right);{\mathrm{\Pi }}_{\mathrm{1}}=V,{\mathrm{\Pi }}_{\mathrm{2}}=\mathit{\epsilon }/a.\end{array}$$

Figure 1 presents a sketch of what the function $\mathrm{\Psi }\left(V,\mathit{\epsilon }/a\right)$ may look like, qualitatively. The underlying idea is that the Pleistocene history of the climate system may be understood as a trajectory in the $\left[V,\mathit{\epsilon }/a\right]$ space (Crucifix and Verbitsky, 2019). The shape and location of the period-doubling domain Ψ=2 is expected to depend on the forcing period.

Figure 1 A typical illustrative $\mathrm{\Psi }\left(V,\mathit{\epsilon }/a\right)$ function. Red arrow represents hypothetical trajectory of the Pleistocene history of the system, which includes the following: from doubled precession periods of the early Pleistocene to doubled obliquity periods of the late Pleistocene.

It is interesting that Fig. 1 is consistent with a similar map produced by a conceptual model built on a completely different principle, i.e., the simple oscillator type model of Daruka and Ditlevsen (2016). In both cases, the obliquity period doubling requires relatively intense external forcing in combination with the relatively high V number (or reduced damping in the case of Daruka and Ditlevsen, 2016). This similarity implies that the importance of the V number for a climate system dynamics may extend well beyond the Verbitsky et al. (2018) model.

3.2 Amplitude of the glacial area variations, Ś

We begin again with the most general hypothesis. We suggest that the amplitude of glacial area variations Ś is a function of the V number, the characteristic rate of snow precipitation, a, the amplitude of the external forcing ε, and the period of the system response P as it is described by Eq. (7). The relationship between the period of the response and that of the forcing may therefore be nontrivial. It means that the system response may exhibit original forcing periods or multiples of them.

$$\begin{array}{}\text{(8)}& \mathit{Ś}=\mathit{\phi }\left(V,a,\mathit{\epsilon },P\right)\end{array}$$

If the hypothesis Eq. (8) is true, then, taking dimensions of ε and P as independent dimensions and using the π theorem, we obtain:

$$\mathit{Ś}/\left({\mathit{\epsilon }}^{\mathrm{2}}{P}^{\mathrm{2}}\right)=\mathrm{\Phi }\left(V,\mathit{\epsilon }/a\right),$$

and finally:

$$\begin{array}{}\text{(9)}& \mathit{Ś}={\mathit{\epsilon }}^{\mathrm{2}}{P}^{\mathrm{2}}\mathrm{\Phi }\left({\mathrm{\Pi }}_{\mathrm{1}},{\mathrm{\Pi }}_{\mathrm{2}}\right).\end{array}$$

Neither Π₁ nor Π₂ contains P. Equation (9) therefore implies that, at constant amplitude of the external forcing ε, the amplitude of glacial area variations is scale invariant with a frequency slope equal to 2. Figure 2 (Ś with reference parameter values) presents a numerical test of the hypothesis Eq. (8) and of its implication Eq. (9). Here, we measure the system response to single-sinusoid forcings of constant amplitude and periods T varying from 5 to 50 kyr. The system responds to this forcing with periods P ranging from 5 to 100 kyr, because forcing periods T of 40 and 50 kyr produce response periods P of 80 and 100 kyr, correspondingly. It can be seen that the Ś-amplitude frequency slope, β_a, is close to 2 (i.e., β_a=1.8) for periods between 30 and 100 kyr. It means that the amplitude of glacial area variations is scale invariant in the orbital domain.

Figure 2 The system response to a single-sinusoid external forcing of constant amplitude and different periods: (1) Ś with reference parameter values; (2) $\stackrel{\mathrm{´}}{\mathit{\omega }}$ with reference parameter values; (3) Ś with intensive climate temperature and weak albedo positive feedbacks; (4) $\stackrel{\mathrm{´}}{\mathit{\omega }}$ with intensive climate temperature and weak albedo positive feedbacks; (5) Ś with weak climate temperature and intensive albedo positive feedbacks; (6) $\stackrel{\mathrm{´}}{\mathit{\omega }}$ with weak climate temperature and intensive albedo positive feedbacks; (7) Ś with intensive climate temperature positive and ice-sheet basal temperature negative feedbacks; (8) $\stackrel{\mathrm{´}}{\mathit{\omega }}$ with intensive climate temperature positive and ice-sheet basal temperature negative feedbacks; (9) Ś with weak climate temperature positive and ice-sheet basal temperature negative feedbacks; and (10) $\stackrel{\mathrm{´}}{\mathit{\omega }}$ with weak climate temperature positive and ice-sheet basal temperature negative feedbacks.

3.3 Amplitude of the basal temperature, $\stackrel{\mathrm{´}}{\mathit{\theta }}$

The amplitude spectrum of the θ variable cannot be derived unambiguously from the same simple considerations as we have employed for P and Ś because of the following reasons: (a) we cannot constrain ourselves with only parameters V, a, ε, and P, since the basal temperature θ is measured in degrees Celsius, but none of ε, a, or P contains degrees Celsius and (b) as soon as we disassemble the V number, i.e., use all individual model parameters instead of V, the advantage of using the π theorem is lost. Nevertheless, if we disassemble the V number wisely, we can minimize the number of dimensional parameters, and, as a result, we may be rewarded by discovering the identities of critical groups that define the scaling properties of θ. Accordingly, we will disassemble the V number using not the individual parameters involved but instead using dimensionless groups that are present in the V number: $\mathit{\alpha },\frac{\mathit{\kappa }}{c},\frac{{\mathit{\gamma }}_{\mathrm{1}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}},S_{\mathrm{0}}},$ and$\frac{{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}}}$. If we consider that the group $\frac{{\mathit{\gamma }}_{\mathrm{1}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}}{S}_{\mathrm{0}}}$ is a dimensionless representation of the parameter γ₁ and the group $\frac{{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}}}$ is a dimensionless representation of the parameter β, then the remaining parameters ${\mathit{\gamma }}_{\mathrm{2}},{\mathit{\gamma }}_{\mathrm{3}},$ and S₀ need to be represented individually in the dimensional form. Taking this together, this yields the following hypothesis:

$$\begin{array}{}\text{(10)}& \stackrel{\mathrm{´}}{\mathit{\theta }}=\mathit{\chi }\left(\mathit{\alpha },{\displaystyle \frac{\mathit{\kappa }}{c}},{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{1}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}}{S}_{\mathrm{0}}}},{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}}}},{\mathit{\gamma }}_{\mathrm{2}},{\mathit{\gamma }}_{\mathrm{3}},S_{\mathrm{0}},a,\mathit{\epsilon },P\right).\end{array}$$

Taking γ₂,S₀, and P as independent dimensions, the π theorem implies:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\stackrel{\mathrm{´}}{\mathit{\theta }}& ={\mathit{\gamma }}_{\mathrm{2}}{S}_{\mathrm{0}}P\\ & X\left(\mathit{\alpha },{\displaystyle \frac{\mathit{\kappa }}{c}},{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{1}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}}{S}_{\mathrm{0}}}},{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}}}},{\mathit{\gamma }}_{\mathrm{3}}P,aP{S}_{\mathrm{0}}^{-\mathrm{1}/\mathrm{2}},\mathit{\epsilon }P{S}_{\mathrm{0}}^{-\mathrm{1}/\mathrm{2}}\right),\end{array}\end{array}$$

or, combining groups $\mathit{\alpha },\frac{\mathit{\kappa }}{c},\frac{{\mathit{\gamma }}_{\mathrm{1}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}}{S}_{\mathrm{0}}},$ and$\frac{{\mathit{\gamma }}_{\mathrm{2}}}{\mathit{\beta }{\mathit{\gamma }}_{\mathrm{3}}}$ back into V number as follows:

$$\begin{array}{}\text{(12)}& \stackrel{\mathrm{´}}{\mathit{\theta }}={\mathit{\gamma }}_{\mathrm{2}}{S}_{\mathrm{0}}PX\left(V,{\mathit{\gamma }}_{\mathrm{3}}P,aP{S}_{\mathrm{0}}^{-\mathrm{1}/\mathrm{2}},\mathit{\epsilon }P{S}_{\mathrm{0}}^{-\mathrm{1}/\mathrm{2}}\right).\end{array}$$

Since ${\mathrm{\Pi }}_{\mathrm{2}}=\mathit{\epsilon }/a$,

$$\begin{array}{}\text{(13)}& \stackrel{\mathrm{´}}{\mathit{\theta }}={\mathit{\gamma }}_{\mathrm{2}}{S}_{\mathrm{0}}PX\left({\mathrm{\Pi }}_{\mathrm{1}},{\mathrm{\Pi }}_{\mathrm{2}},{\mathrm{\Pi }}_{\mathrm{3}},{\mathrm{\Pi }}_{\mathrm{4}}\right),\end{array}$$

where Π₃=γ₃P, and ${\mathrm{\Pi }}_{\mathrm{4}}=\mathit{\epsilon }P{S}_{\mathrm{0}}^{-\mathrm{1}/\mathrm{2}}$.

As Π₃ and Π₄ include P, then, generally speaking, the amplitude of basal temperature variations is not expected to be scale invariant. Under some circumstances though, the function $X\left({\mathrm{\Pi }}_{\mathrm{1}},{\mathrm{\Pi }}_{\mathrm{2}},{\mathrm{\Pi }}_{\mathrm{3}},{\mathrm{\Pi }}_{\mathrm{4}}\right)$ may become P-independent and the amplitude of the basal temperature variations may develop the property of scale invariance. For example, we observed experimentally that when Π₄→0 (e.g., the amplitude of the external forcing, ε, is reduced), the Eq. (13) becomes scale invariant with a frequency slope equal to 1. In this case $\stackrel{\mathrm{´}}{\mathit{\theta }}={\mathit{\gamma }}_{\mathrm{2}}{S}_{\mathrm{0}}PX\left({\mathrm{\Pi }}_{\mathrm{1}},{\mathrm{\Pi }}_{\mathrm{2}},{\mathrm{\Pi }}_{\mathrm{3}}/{\mathrm{\Pi }}_{\mathrm{4}}\right)$.

3.4 Amplitude of the climate temperature, $\stackrel{\mathrm{´}}{\mathit{\omega }}$

Since Eq. (3) for ω is linear, it may provide us with a hint about the response scaling characteristics of this variable. In the orbital domain ${\mathrm{\Pi }}_{\mathrm{3}}={\mathit{\gamma }}_{\mathrm{3}}P\gg \mathrm{1}$ so that Eq. (3) may be approximated as ${\mathit{\gamma }}_{\mathrm{3}}\mathit{\omega }\approx {\mathit{\gamma }}_{\mathrm{1}}-{\mathit{\gamma }}_{\mathrm{2}}\left(S-S_{\mathrm{0}}\right)$. Hence, $\stackrel{\mathrm{´}}{\mathit{\omega }}=\frac{{\mathit{\gamma }}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{3}}}\mathit{Ś}$. We may hypothesize therefore that in the orbital domain and possibly even beyond:

$$\begin{array}{}\text{(14)}& \stackrel{\mathrm{´}}{\mathit{\omega }}=\mathit{\nu }\left(V,{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{3}}}},a,\mathit{\epsilon },P\right).\end{array}$$

Taking the dimensions of $\frac{{\mathit{\gamma }}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{3}}},\mathit{\epsilon }$, and P as independent and applying again π-theorem reasoning, we should expect that

$$\begin{array}{}\text{(15)}& \stackrel{\mathrm{´}}{\mathit{\omega }}={\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{3}}}}{\mathit{\epsilon }}^{\mathrm{2}}{P}^{\mathrm{2}}N\left({\mathrm{\Pi }}_{\mathrm{1}},{\mathrm{\Pi }}_{\mathrm{2}}\right).\end{array}$$

At constant amplitude of the external forcing ε, Eq. (15) implies that the amplitude of climate temperature variations $\stackrel{\mathrm{´}}{\mathit{\omega }}$ grows with the square of the response period. The results presented in Fig. 2 ($\stackrel{\mathrm{´}}{\mathit{\omega }}$ with reference parameter values) support the hypothesis (14) and its implication (15): the ω variable amplitude frequency slope is close to 2 (i.e., β_a=1.8) for periods between 5 and 100 kyr. It means that in the orbital and millennial domains, the amplitude of the climate temperature is scale invariant.

4 Discussion

4.1 Scale invariance and a physical nature of the climate system feedbacks

So far, we have based our implications of scaling relationships on the significance of a dimensionless number (in our case, the V number) quantifying a mean ratio between positive and negative feedbacks. That is, the scaling relationships found should be robust across changes in the composition of V, provided that the value of V is unchanged. To illustrate this implication, we conducted four numerical experiments. In the first experiment, we increase coefficients α and κ 2-fold and reduce γ₂ by a half relative to their reference values. This does not change the reference value of the V number (see the Eq. 4, and note that the reference value of γ₁=0) that is V=0.75 but transforms system Eqs. (1)–(3) into a system where the positive feedback is dominated by the climate temperature affecting ice-sheet mass balance and its temperature regime. We then measure the system response to the single-sinusoid forcing of the same amplitude and periods T=5–50 kyr. (Note that periods T=40 and 50 kyr produce system response of periods P=80 and 100 kyr, correspondingly). In the second experiment, we decrease coefficients α and κ by 50 % and increase γ₂ 2-fold relative to their reference values. Again, this does not change the reference value of the V number, V=0.75, but transforms system Eqs. (1)–(3) to a system where the positive feedback is dominated by the albedo feedback. In the third experiment, we increase coefficients α and κ by 50 % as well as the coefficient β, thus creating the system with intensive climate–temperature positive feedback and intensive ice-sheet basal temperature negative feedback, the V number still being equal to 0.75. And finally, we decrease coefficients α, κ, and β by 50 %, making a system with weak climate–temperature positive and ice-sheet basal temperature negative feedbacks. The response of all four systems to the external forcing is shown in Fig. 2. Despite different underlying physics, all four systems demonstrate the same outcomes: in the orbital domain, their amplitudes of glacial area variations are scale invariant with 1.8 frequency slope, and the amplitudes of the climate temperature are scale invariant in the orbital and millennial domains with the same slope.

This robustness is comforting. As we know, the physical interpretation of a low-order dynamical model can be partly ambiguous. For example, the mechanisms responsible for the changes in the effective climate temperature and how it impacts the ice mass balance are not fully described in this model. It is therefore reassuring to have been able to identify what seems to be the key ingredient for the scaling relationship, in this case, that a single quantity (the V number) grossly determines the dynamics of the system response. In other words, it relies on the fact that the number of effective parameters is smaller than is apparent from a more detailed description of the system.

This, incidentally, shows how difficult it is to disambiguate the physical mechanisms responsible for a given behavior. Different assemblages yielding the same V number will, indeed, produce slightly different solutions but less different than one could have perhaps expected. The dimensionless functions like, for example, function $\mathrm{\Phi }\left(V,\mathit{\epsilon }/a\right)$ in the Eq. (9),

$$\mathit{Ś}={\mathit{\epsilon }}^{\mathrm{2}}{P}^{\mathrm{2}}\mathrm{\Phi }\left(V,\mathit{\epsilon }/a\right),$$

and function ${\mathrm{\Phi }}^{\prime }\left(V,\mathit{\epsilon }/a\right)$, corresponding to the same value of the V number but formed by the different physics (different set of parameters),

$$\mathit{Ś}={\mathit{\epsilon }}^{\mathrm{2}}{P}^{\mathrm{2}}{\mathrm{\Phi }}^{\prime }\left(V,\mathit{\epsilon }/a\right),$$

though not identical, yield the same scaling behavior. If the amplitude of the external forcing ε is constant, the period P shows up only as a power-law monomial ∼Pⁿ, and its power n makes the same scale-invariant amplitude-spectrum slope regardless of the specific physics defining the V number. In other words, though the functions $\mathit{\psi }\left(V,a,\mathit{\epsilon },T\right)$, $\mathit{\phi }\left(V,a,\mathit{\epsilon },P\right)$, $\mathit{\chi }\left(V,a,\mathit{\epsilon },P\right)$, and $\mathit{\nu }(V,\frac{{\mathit{\gamma }}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{3}}},a,\mathit{\epsilon },P)$ may change depending on the specific physics forming the V number, their governing parameters always remain the same because they are determined by the structure of the system Eqs. (1)–(3). Accordingly, the functions Ψ(Π₁,Π₂), Φ(Π₁,Π₂), $X\left({\mathrm{\Pi }}_{\mathrm{1}},{\mathrm{\Pi }}_{\mathrm{2}},{\mathrm{\Pi }}_{\mathrm{3}},{\mathrm{\Pi }}_{\mathrm{4}}\right)$, and N(Π₁,Π₂) may also change, but their dimensionless arguments (Π groups) remain unaffected. As long as their groups, like, for example, Π₁ and Π₂, do not contain P, we have a possibility of scale-invariance. This observation makes the scale invariance a very general and expected property of the climate system.

The physical interpretation of the dynamical model we employ in this study (Verbitsky et al., 2018) is very straightforward as far as Eqs. (1) and (2) are concerned; these are scaled equations of mass and energy conservation of viscous ice flow. We must admit, though, that Eq. (3) of the climate temperature is, indeed, ambiguous. In other words, we are uncertain about some key mechanisms that we have chosen to describe using the rest-of-the-climate linear equation. Among others, these may be nonlinear effects related to the carbon cycle, nonlinear effects of sea-level destabilization of ice sheets and related synchronization, nonlinear effects associated with atmospheric circulation, or nonlinear effects related to biogenic calcifiers and their action on alkalinity, etc. A challenger might thus claim that these effects are so important that they should be considered more explicitly. Indeed, we have the hope that even after accounting for these processes, we might end up with a model that still has grossly the same mathematical structure as the Verbitsky et al. (2018) model, even though the meaning of some of the variables will have changed. Specifically, since Eq. (3) is linear, it can be split into several equations:

$$\begin{array}{c}{\displaystyle }\mathit{\omega }={\mathit{\omega }}_{\mathrm{1}}+{\mathit{\omega }}_{\mathrm{2}}+\mathrm{\dots }+{\mathit{\omega }}_n\\ {\displaystyle \frac{\mathrm{d}{\mathit{\omega }}_{\mathrm{1}}}{\mathrm{d}t}}={\mathit{\gamma }}_{\mathrm{11}}-{\mathit{\gamma }}_{\mathrm{21}}\left(S-S_{\mathrm{0}}\right)-{\mathit{\gamma }}_{\mathrm{3}}{\mathit{\omega }}_{\mathrm{1}}\\ {\displaystyle \frac{\mathrm{d}{\mathit{\omega }}_{\mathrm{2}}}{\mathrm{d}t}}={\mathit{\gamma }}_{\mathrm{12}}-{\mathit{\gamma }}_{\mathrm{22}}\left(S-S_{\mathrm{0}}\right)-{\mathit{\gamma }}_{\mathrm{3}}{\mathit{\omega }}_{\mathrm{2}}\\ {\displaystyle }\mathrm{\dots }\\ {\displaystyle \frac{\mathrm{d}{\mathit{\omega }}_n}{\mathrm{d}t}}={\mathit{\gamma }}_{\mathrm{1}n}-{\mathit{\gamma }}_{\mathrm{2}n}\left(S-S_{\mathrm{0}}\right)-{\mathit{\gamma }}_{\mathrm{3}}{\mathit{\omega }}_n.\end{array}$$

Each of the above equations may represent different feedback mechanisms. Therefore our experiments with increased (or reduced) γ₂ may be also understood as experiments with additional feedbacks of different nature (${\mathit{\gamma }}_{\mathrm{2}}={\mathit{\gamma }}_{\mathrm{21}}+{\mathit{\gamma }}_{\mathrm{22}}+\mathrm{\dots }+{\mathit{\gamma }}_{\mathrm{2}n}$), though of the same timescale 1∕γ₃.

4.2 Multi-sinusoid forcing

Thus far we have assumed a single-sinusoid external forcing with an amplitude ε and a period T. When we force our system with normalized mid-July insolation at 65^∘ N (Berger and Loutre, 1991), this assumption is not valid any longer because both the amplitudes and the periods of precession and obliquity are different. Therefore, the hypothesis Eq. (5) must be rewritten as:

$$\begin{array}{}\text{(16)}& P=\mathit{\psi }\left[V,a,{\mathit{\epsilon }}_{\mathrm{1}},T_{\mathrm{1}},{\mathit{\epsilon }}_{\mathrm{2}},T_{\mathrm{2}}\right].\end{array}$$

Here P is a period of the system response to a specific forcing component (a peak of the response spectrum); index “1” corresponds to obliquity, and index “2” corresponds to precession. Taking dimensions of ε₁ and T₁ as independent dimensions and using the π theorem, we obtain:

$$\begin{array}{}\text{(17)}& P_{\mathrm{1}}=T_{\mathrm{1}}{\mathrm{\Psi }}_{\mathrm{1}}\left[V,{\mathit{\epsilon }}_{\mathrm{1}}/a,{\mathit{\epsilon }}_{\mathrm{1}}/{\mathit{\epsilon }}_{\mathrm{2}},T_{\mathrm{1}}/T_{\mathrm{2}}\right].\end{array}$$

Here P₁ is a period of the system response to the obliquity forcing. Similarly, taking dimensions of ε₂ and T₂ as independent dimensions, and using the π theorem, we have

$$\begin{array}{}\text{(18)}& P_{\mathrm{2}}=T_{\mathrm{2}}{\mathrm{\Psi }}_{\mathrm{2}}\left[V,{\mathit{\epsilon }}_{\mathrm{2}}/a,{\mathit{\epsilon }}_{\mathrm{1}}/{\mathit{\epsilon }}_{\mathrm{2}},T_{\mathrm{1}}/T_{\mathrm{2}}\right].\end{array}$$

Here P₂ is a period of the system response to the precession forcing. Since in the case of the orbital forcing ε_1∕ε₂ and T₁∕T₂ are invariant, we can apply generalized π theorem (Sonin, 2004) and to rewrite Eqs. (17) and (18) as

$$\begin{array}{}\text{(19)}& {\displaystyle }{P}_{\mathrm{1}}=T_{\mathrm{1}}{\mathrm{\Psi }}_{\mathrm{1}}\left[V,{\mathit{\epsilon }}_{\mathrm{1}}/a\right]\text{(20)}& {\displaystyle }{P}_{\mathrm{2}}=T_{\mathrm{2}}{\mathrm{\Psi }}_{\mathrm{2}}\left[V,{\mathit{\epsilon }}_{\mathrm{2}}/a\right].\end{array}$$

It can be seen that Eqs. (19) and (20) are identical to Eq. (7), and the response periods to obliquity and to precession do not depend on each other. This result is not by any means intuitive.

We now repeat the same reasoning for the corresponding amplitudes of the system response:

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathit{Ś}}_{\mathrm{1}}={\mathit{\phi }}_{\mathrm{1}}\left(V,a,{\mathit{\epsilon }}_{\mathrm{1}},P_{\mathrm{1}},{\mathit{\epsilon }}_{\mathrm{2}},P_{\mathrm{2}}\right)\text{(22)}& {\displaystyle }{\mathit{Ś}}_{\mathrm{2}}={\mathit{\phi }}_{\mathrm{2}}\left(V,a,{\mathit{\epsilon }}_{\mathrm{1}},P_{\mathrm{1}},{\mathit{\epsilon }}_{\mathrm{2}},P_{\mathrm{2}}\right)\text{(23)}& {\displaystyle }{\mathit{Ś}}_{\mathrm{1}}={\mathit{\epsilon }}_{\mathrm{1}}^{\mathrm{2}}{P}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{1}}\left(V,{\mathit{\epsilon }}_{\mathrm{1}}/a,{\mathit{\epsilon }}_{\mathrm{1}}/{\mathit{\epsilon }}_{\mathrm{2}},P_{\mathrm{1}}/P_{\mathrm{2}}\right)\text{(24)}& {\displaystyle }{\mathit{Ś}}_{\mathrm{2}}={\mathit{\epsilon }}_{\mathrm{2}}^{\mathrm{2}}{P}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{2}}\left(V,{\mathit{\epsilon }}_{\mathrm{2}}/a,{\mathit{\epsilon }}_{\mathrm{1}}/{\mathit{\epsilon }}_{\mathrm{2}},P_{\mathrm{1}}/P_{\mathrm{2}}\right).\end{array}$$

Though in the case of the orbital forcing ε₁∕ε₂ and T₁∕T₂ are invariant, P₁∕P₂ is not an invariant (see Fig. 1); therefore

$$\begin{array}{}\text{(25)}& {\displaystyle }{\mathit{Ś}}_{\mathrm{1}}={\mathit{\epsilon }}_{\mathrm{1}}^{\mathrm{2}}{P}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{1}}\left(V,{\mathit{\epsilon }}_{\mathrm{1}}/a,P_{\mathrm{1}}/P_{\mathrm{2}}\right)\text{(26)}& {\displaystyle }{\mathit{Ś}}_{\mathrm{2}}={\mathit{\epsilon }}_{\mathrm{2}}^{\mathrm{2}}{P}_{\mathrm{2}}^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{2}}\left(V,{\mathit{\epsilon }}_{\mathrm{2}}/a,P_{\mathrm{1}}/P_{\mathrm{2}}\right).\end{array}$$

We can see that although periods of the system response to the precession and obliquity forcings are independent, the amplitudes of the corresponding variations are interdependent and thus may deviate from a pure square-period law. This observation may have an important implication for our understanding of the paleodata. As we demonstrated before (Verbitsky et al., 2018), P₁∕P₂ evolves over time, specifically $P_{\mathrm{1}}/P_{\mathrm{2}}=\mathrm{1}$ for the early Pleistocene due to precession period doubling and $P_{\mathrm{1}}/P_{\mathrm{2}}=\mathrm{4}$ for the late Pleistocene due to obliquity period doubling. It means that the slope of the spectrum of the system response may also evolve.

Introduction of more sinusoids (for example, accounting for the millennial forcing) makes the situation even more complex. In such a case, a period of the system response to a specific forcing component depends on the amplitudes and the periods of all sinusoids:

$$\begin{array}{}\text{(27)}& P=\mathit{\psi }\left[V,a,{\mathit{\epsilon }}_{\mathrm{1}},T_{\mathrm{1}},{\mathit{\epsilon }}_{\mathrm{2}},T_{\mathrm{2}},\mathrm{\dots }{\mathit{\epsilon }}_i,T_i\mathrm{\dots }\right].\end{array}$$

Then, for example, P₁, the period of the system response to obliquity forcing, can be presented as

$$\begin{array}{}\text{(28)}& P_{\mathrm{1}}=T_{\mathrm{1}}{\mathrm{\Psi }}_{\mathrm{1}}\left[V,{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{1}}}{a}},\mathrm{\dots },{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{1}}}{{\mathit{\epsilon }}_i}},{\displaystyle \frac{T_{\mathrm{1}}}{T_i}},\mathrm{\dots }\right],\end{array}$$

and corresponding amplitude of the glaciation area response

$$\begin{array}{}\text{(29)}& {\displaystyle }{\mathit{Ś}}_{\mathrm{1}}={\mathit{\phi }}_{\mathrm{1}}\left[V,a,{\mathit{\epsilon }}_{\mathrm{1}},P_{\mathrm{1}},{\mathit{\epsilon }}_{\mathrm{2}},P_{\mathrm{2}},\mathrm{\dots }{\mathit{\epsilon }}_i,P_i\mathrm{\dots }\right]\text{(30)}& {\displaystyle }{\mathit{Ś}}_{\mathrm{1}}={\mathit{\epsilon }}_{\mathrm{1}}^{\mathrm{2}}{P}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{\Phi }}_{\mathrm{1}}\left[V,{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{1}}}{a}},\mathrm{\dots },{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{1}}}{{\mathit{\epsilon }}_i}},{\displaystyle \frac{P_{\mathrm{1}}}{P_i}},\mathrm{\dots }\right].\end{array}$$

Equations (28) and (30) show that, generally speaking, every peak P and corresponding amplitude Ś of the system response depend on each forcing sinusoid. Such a dependence may break the scale invariance we discussed earlier. For example, we have demonstrated in our previous study (Verbitsky et al., 2019a) that introduction of the millennial variability of significant amplitude (i.e., ${\mathit{\epsilon }}_{\mathrm{1}}/{\mathit{\epsilon }}_i\to \mathrm{0}$) may disrupt the response of the system to the orbital forcing and essentially reduce the slope β_a. The empirical energy density spectrum of Huybers and Curry (2006) has a slope of B≈2 in the orbital domain. Since the energy density slope B relates to the fluctuation amplitude slope β_a as $B=\mathrm{2}{\mathit{\beta }}_{\mathrm{a}}+\mathrm{1}$, B≈2 corresponds to ${\mathit{\beta }}_{\mathrm{a}}=\mathrm{0.5}<\mathrm{2}$. We may therefore speculate that the observed spectrum of the climate variability could be significantly influenced by the millennial forcing propagated into the orbital domain.

4.3 How general is the property of scale invariance?

It is apparent that not every dynamical model has the property of scale invariance, which is encoded in its dynamical equations. As an illustration, let us consider the van der Pol oscillator. It was previously suggested as a minimal model capturing ice-age dynamics (Crucifix, 2012).

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}={\displaystyle \frac{-y+\mathit{\beta }+\mathit{\gamma }F}{\mathit{\tau }}}\text{(32)}& {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}}={\displaystyle \frac{-\mathit{\alpha }\left(\frac{y^{\mathrm{3}}}{\mathrm{3}}-y-x\right)}{\mathit{\tau }}}\end{array}$$

Here all variables and parameters, except τ, are dimensionless; τ is measured in units of time. Variable x is thought to represent the global ice volume, and variable y makes the “rest-of-the climate” response. Using the same π-theorem technique, let us determine the period P and the amplitude x^′ of the system response to the external forcing F of the period T.

$$\begin{array}{}\text{(33)}& {\displaystyle }P=\mathit{\psi }\left(\mathit{\alpha },\mathit{\beta },\mathit{\gamma },\mathit{\tau },T\right)\text{(34)}& {\displaystyle }P=T\mathrm{\Psi }\left(\mathit{\alpha },\mathit{\beta },\mathit{\gamma },\mathit{\tau }/T\right)\end{array}$$

Since α, β, and γare constants,

$$\begin{array}{}\text{(35)}& P=T\mathrm{\Psi }\left(\mathit{\tau }/T\right).\end{array}$$

Similarly,

$$\begin{array}{}\text{(36)}& {\displaystyle }{x}^{\prime }=\mathit{\phi }\left(\mathit{\alpha },\mathit{\beta },\mathit{\gamma },\mathit{\tau },P\right)\text{(37)}& {\displaystyle }{x}^{\prime }=\mathrm{\Phi }\left(\mathit{\alpha },\mathit{\beta },\mathit{\gamma },\mathit{\tau }/P\right)=\mathrm{\Phi }\left(\mathit{\tau }/P\right).\end{array}$$

It means that the amplitudes of forced fluctuations in the van der Pol model are not necessarily scale invariant. We have tested this conclusion experimentally for τ=36.2 kyr and a forcing period T ranging from 5 to 100 kyr. The response shows slope breaks near approximately 90 and 50 kyr, which are clearly related to the auto-oscillation of the 100 kyr dominant period and its 50 kyr overtone.

Therefore, in a search for the most adequate ice-age physics, it would indeed be useful to see whether more sophisticated ice sheet–ocean–atmosphere models have the property of scale invariance. We suspect that potential universality of this property may stem from the universality of the Eq. (1). Equation (1) represents the global ice volume balance and simply says that changes in the ice volume are equal to the mass influx to the ice-sheet surface. This statement is valid for each and every climate model of any complexity. Therefore, if a model can be diagnosed with a single dimensionless number similar to the V number that would effectively capture most of the climate dynamics, then the scale invariance of the glaciation area variations (in square meters) can be reduced from the simple observation that it depends on the mass influx to its surface (in meters per second) and the periodicity of the mass influx variation (in seconds). This might not be too difficult to verify with an adequate set of experiments, but we must obviously leave this task to the scientists who know and develop these models.

5 Conclusions

Twenty-seven years ago, Saltzman and Verbitsky (1993) discussed their model of 18 parameters, nine of which were physically unconstrained (i.e., free parameters) and formulated a challenge “to account for as much of the variance with fewer free parameters. A challenge is thus posed to ourselves and other theoretical paleoclimatologists to construct a more parsimonious model in this regard that can supersede our present effort.” We may now conclude that this challenge has been met. All parameters in our model (Verbitsky et al., 2018) are physically constrained. Moreover, dimensional analysis reveals that there are only two factors that define most of the ice-age dynamics: (a) a balance between intensities of climate positive and ice sheet negative feedbacks, Π₁=V; (b) the period, T, and the amplitude of the external forcing, ε, (specifically, a particular proportion between the external, e.g., orbital, and terrestrial ice sheet mass balance components, ${\mathrm{\Pi }}_{\mathrm{2}}=\mathit{\epsilon }/a)$.

The analysis indicates that the amplitudes of glacial area variations and of climate temperature are scale invariant with a frequency slope of approximately 2. The property of scale invariance does not depend on the physical nature of the underlying positive and negative feedbacks incorporated by the system. It thus turns out to be one of the most fundamental properties of the Pleistocene climate.

Retrospectively, we could have inferred scale invariance from the mere assumption that the behavior of the continental glacial area (measured in square meters) depends on the mass influx to its surface (in meters per second) and the periodicity of the mass influx variation (in seconds), but perhaps these assumptions are too simple to be convincing. In our study, we have chosen a bit more sophisticated but more credible approach. We derived a dynamical model from the scaled conservation equations of viscous non-Newtonian ice combined with an equation describing the evolution of the climate temperature. We observed that most of the dynamical system behavior can be explained by a balance between positive and negative feedbacks. This observation, finally, illuminated the crucial role of the mass influx and its periodicity, making application of the π theorem effective and definitive.

Certainly, we cannot claim to have a full picture of the mechanisms of ice ages, but if ice age physics are well captured by the mathematical structure that we have obtained, then this scale invariance linking response amplitudes and periods applies. We further suggest that a model that would indeed be a bit different than the Verbitsky et al. (2018) model (because it includes some other important, may be nonlinear, mechanisms) might still retain an important property that we have discovered: there is a connection between the sensitivity of the fixed point (since the V number is indeed constructed by consideration of the sensitivity of the fixed point) and a scale invariance linking period and amplitude of response. This seems to be the fundamental proposal, for which we welcome challengers equipped with bigger models.