A Statistical Mechanical Approach for the Computation of the Climatic Response to General Forcings

The climate belongs to the class of non-equilibrium forced and dissipative systems, for which most results of quasi-equilibrium statistical mechanics, including the fluctuation-dissipation theorem, do not apply. We show for the first time how the Ruelle linear response theory, developed for studying rigorously the impact of perturbations on general observables of non-equilibrium statistical mechanical systems, can be applied to analyze the climatic response. We choose as test bed the Lorenz 96 model, which has a well-recognized prototypical value. We recapitulate the main aspects of the response theory and propose some new results. We then analyze the frequency dependence of the response of both local and global observables to perturbations with localized as well as global spatial patterns. We derive analytically the asymptotic behaviour, validity of Kramers-Kronig relations, and sum rules for the susceptibilities, and related them to parameters describing the unperturbed properties of the system. We verify the theoretical predictions from the outputs of the simulations with great precision. The theory is used to explain differences in the response of local and global observables, in defining the intensive properties of the system and in generalizing the concept of climate sensitivity to all time scales. We also show how to reconstruct the linear Green function, which maps perturbations of general time patterns into changes in the expectation value of the considered observable. Finally, we propose a general methodology to study Climate Change problems by resorting to few, well selected simulations and discuss the specific case of surface temperature response to changes of the $CO_2$ concentration. This approach may provide a radically new perspective to study rigorously the problem of climate sensitivity and climate change.


Contents 1 Introduction 4 1 Introduction
A crucial goal in the study of general dynamical and statistical mechanical systems is to understand how their statistical properties are altered when we introduce a perturbation related to changes in the external forcing or in the value of some internal parameters. The ability to compute the response of the system is of great relevance for purely mathematical reasons as well as in many fields of science and technology.
The climate system is an outstanding example of a non-equilibrium, forced and dissipative complex system, forced in first instance by spatial differences and temporal variability in the net energy flux at the top of the atmosphere. On a macroscopic level, as a result of being far from equilibrium, the climate system behaves as an engine, driven by the temperature difference between a warm and a cold thermal pool, so that the atmospheric and oceanic motions are at the same time the result of the mechanical work (then dissipated in a turbulent cascade) produced by the engine, and are processes which re-equilibrate the energy balance of the climate system [1,2,3,4].
A primary goal of climate science is to understand how the statistical propertiesmean values, fluctuations, and higher order moments -of the climate system change as a result of modulations to some crucial external (e.g. solar irradiance) or internal (e.g. atmospheric composition) parameters of the system occurring on various time scales. A large class of problems -those involving climate sensitivity, climate variability, climate change, climate tipping points -fall into this category. In a system as complex and as extended as the climate, where lots of feedbacks are active on a variety of spatial and temporal scales, this is in general a very difficult task. The need for scientific advance in this direction is outstanding as one considers that even after several decades of intense scientific efforts, the accurate evaluation of the climate sensitivity par excellence, i.e., the change of the globally averaged surface temperature for doubling of CO 2 concentration with respect to pre-industrial levels (280 ppm to 560 ppm circa), is a tantalizing endeavor, and large uncertainties are still present [5].
Such efforts have significant relevance also in the context of the ever-increasing attention paid by the scientific community to the quest for reliable metrics to be used for the validation of climate models of various degrees of complexity and for the definition of strategies aimed at the radical improvement of their performance [6,7]. The pursuit of a quantum leap in climate modelling -which definitely requires new scientific ideas case.
Ruelle's analysis applies for non-equilibrium steady state systems [19] possessing a Sinai-Ruelle-Bowen (SRB) invariant measure, often referred to as Axiom A system [20,21]. This class of systems, even if mathematically non-generic, includes on the other hand excellent models for general physical systems, as made clear by the chaotic hypothesis [22,23], which can be interpreted as an extension of the ergodic hypothesis to non-Hamiltonian systems [19]. See [24] for an original geophysical perspective.
The Ruelle response theory, with the support of chaotic hypothesis, has interesting conceptual implications for climate studies. In fact, the possibility of defining a response function basically poses the problem of climate response to forcings and of climate change in a well-defined context, and, when considering the procedures aimed at improving climate models, justifies rigorously the procedures of tuning and adjusting the free parameters. Moreover, the response theory allows to compute the climate sensitivity, in the special case when static perturbations to the system parameters are considered. Previously, a response formula was proposed by Cacuci for evaluating the linearized change of the solution of a time-independent generic system of nonlinear equations as a result of a change in the system's parameters [26,27]. This can be interpreted as a special case of Ruelle's theory, where the unperturbed attractor is constituted by a fixed point and a static perturbation to the system evolution equation is considered.
Cacuci proposed to study this problem using the adjoint operator to the original system, which provided an efficient way to determine the impact of small perturbations. Interestingly, early prominent applications of the so-called adjoint method and its extension to time-dependent problems, which allowed for evaluating all possible linear sensitivities of an evolving model in just one simulation, were been proposed for climate related problems. In particular, it was used to evaluate the sensitivities of a simple radiative-convective model possessing an attractor constituted by just one fixed point [28,29], and later, in an inherently heuristic way, for studying the response of a (chaotic) simplified general circulation model to doubling of the CO 2 concentration [30]. Whereas the adjoint method did not find much space in further climatic studies, mostly due to early discouragement for the computational burden of constructing the suitable operators for evaluating the sensitivities, it subsequently reached great success in data assimilation problems for geophysical fluid dynamics [31,32], to the point that a tangent and adjoint model compiler able to automatically generate adjoint model code was has been introduced [33]. More recently, a link between advanced adjoint techniques and the Ruelle theory has been proposed [34].
In the last decade on one side a great effort has been directed at extending the Ruelle response theory for more general classes of dynamical systems (see, e.g., [35,36]), and recent studies [17] have shown that, thanks only to the causal nature of the response, it is possible to apply all the machinery of the theory of Kramers-Kronig (KK) relations [37,38,39,40] for linear and nonlinear processes to study accurately and rigorously the susceptibilities describing in the frequency domain the response of a general observable to a general perturbation.
Moreover, the actual applicability of the theory has been successfully tested in a number of simple dynamical systems case for the linear [41,42] and nonlinear [43] response. Such numerical investigations have clarified that even in systems which are not Axiom A, like the Lorenz 63 system [44], it is possible to successfully use the response theory to construct linear [41] and nonlinear susceptibilities [43] which obey all of the constraints imposed by the KK theory up to a high degree of precision.
These investigations definitely motivate further studies aimed at understanding to what extent the response theory is an efficient tool for analyzing complex and complicated systems. In this paper, we take up such a challenge and consider the Lorenz 96 (L96) system [45,46,47], which provides an excellent and celebrated prototypical model of a one dimensional atmosphere. The variables of the L96 model can be thought as generic meteorological quantities extending around a latitudinal circle and sampled at a regular interval. In spite of not being realistic in the usual sense, the L96 model presents the basic ingredients, such as dissipation, advection and the presence of an external forcing, of the actual atmosphere. For this reason, L96 has quickly become the standard model to be used for predictability studies [48,49,50], when testing data assimilation techniques [14,15,51], and new parameterizations [52].
Although we are unable to prove that the unperturbed L96 is an Axiom A system, in general and for the specific choice of parameters used in our numerical simulations in particular, we adopt the chaotic hypothesis and present the first thorough investigation of a spatially extended system by using the rigorous statistical mechanical methodologies presented in [10,16,17,43]. Moreover, since L96 is a spatially extended system, we also explore the applicability of the response theory in all possible combinations of global/local observables and global/local perturbations. We compute rigorously the corresponding linear susceptibilities, verify the KK relations and the related sum rules, and find an empirical power law, which, as in [53], supports the validity of the chaotic hypothesis, allowing to extend the results obtained for our specific choice of model's parameters to a rather general class of L96 systems. We also show how to go from the frequency back to the time domain, thus deriving from the susceptibility the Green function, which acts as time propagator of the considered perturbation for the considered observable. The Green function allows to predict, in an ensemble mean sense, the change in the observable at any time horizon as a result of a perturbation with the same spatial patter as that considered in the calculation of the susceptibility but featuring a general time modulation. Finally, building upon the results presented here, we propose a simple yet general methodology to study general Climate Change problems on virtually any time scale by resorting to only few, well selected simulations, and by taking full advantage of ensemble methods. The specific case of globally averaged surface temperature response to a general pattern of change of the CO 2 concentration is discussed.
Whereas the paper aims at proposing new methods for tackling classical problems of climate science, most of the results and of the methodologies proposed are of more general interest. In this paper we limit our attention to the linear response. We refer to [17,43] for a theoretical and numerical studies of higher-order effects of perturbations.
The paper is organized as follows. In Sec. 2 we briefly analyze the general theoretical background of the linear response theory and of the properties of the frequency dependent susceptibility and present some new useful results. In Sec. 3 we present the main features of the L96 system, introduce the considered perturbations to the forcing, derive some basic properties of the response of various observables, and present the theoretical predictions. In Sec. 4 we present the results of our numerical investigations and describe how they can be generalized to the entire family of L96 models. In Sec. 5 we provide a relevant example to illustrate how the results presented in this paper can be used to devise simple yet rigorous methods to study the climate response at all time scales on models of any degree of complexity. In Sec. 6 we discuss the conclusions and present perspectives for future work. We consider an Axiom A dynamical system described by the evolution equationẋ = F (x), so that the invariant probability measure ρ 0 of the associated flow is of the SRB type [10]. Let Φ 0 be the expectation value of the general observable Φ defined as ρ 0 (dx) Φ(x). We perturb the flow of the system by adding a on the right hand side of the evolution equation a vector field X(x)f (t), where X(x) defines the pattern of the perturbation, and f (t) is its time modulation. The resulting evolution equation results to beẋ = F (x) + X(x)f (t). Following Ruelle [10], we express the expectation value of Φ(x) in the perturbed system using a perturbative expansions as: Each term of the perturbative series can be expressed as an n-convolution integral of the n th order causal Green function with n delayed time perturbation functions [25,17]. Limiting our attention at the linear case we have: The first order Green function can be expressed as follows: where Λ(•) = X(x)∇(•) takes into account the effects of the perturbative vector field, Θ is the usual Heaviside distribution and Π 0 the unperturbed time evolution operator so that Π 0 F (x) = F (x(t)) for any function F , with x(t) following the unperturbed flow.
Note that it is possible to express the Green function as the expectation value of a nontrivial but computable observable over the unperturbed SRB measure ρ 0 . Therefore the knowledge of the unperturbed features of the flow is sufficient to define the effects of any external perturbation over any observable of our system. In the frequency domain we find that the first term of the perturbative series can be written as: where the Dirac delta implies that we are analyzing the impact of perturbations in the frequency-domain at the frequency ω. The linear susceptibility is defined as: It is important to underline with a thought experiment the computational relevance of the last equations and the importance of the susceptibility function.
Let's suppose we introduce a time dependent perturbation f α (t) to a given pattern of forcing X(x), simulate the system and observe the time response of an arbitrary observable Φ α (1) (t). We now compute the Fourier transform of the observed signal and of the forcing modulation. Inverting Eq.(4), we can find the linear susceptibility . Let's now consider a different time-modulating function of the forcing f β (t) and its corresponding Fourier transform f β (ω). Taking into account Eq. (4), if we multiply f β (ω) times the previously computed function χ the frequency-dependent response of the observable Φ to the forcing X(x) modulated by the new function. By applying the inverse Fourier transform we obtain the timedependent response Φ β (1) (t) without needing any additional simulation.
Moreover, the knowledge of the susceptibility function allows us to reconstruct the G (1) Φ (t) by inverting Eq. (5). Otherwise, the Green function can be obtained directly from observing the response signal by performing a simulation where f (t) = δ(t): in this case (see Eq. (2)) we simply have Φ (1) (t) = G

Kramers-Kronig relations and sum rules
As we see from Eq. (3), for an arbitrary choice of the observable and of the perturbation the corresponding linear Green function is causal. Assuming, on heuristic physical basis, that G Φ (t) ∈ L 2 , we can apply the Titchmarsh theorem [37,38,39,40] and deduce that the linear susceptibility χ According to a general property of Fourier transform we know that the short term behavior of G (1) Φ (t) determines the asymptotic properties of χ (1) Φ (ω). We shall obtain a more quantitative result by exploiting that: where in the second equality we have neglected the fact that the solution is a distribution and considered ω = 0. Therefore, if the Taylor expansion of the Green function in the limit t → 0 + is of the form: the high frequency behavior of the linear susceptibility, i.e. the limit ω → ∞, is: where α =ᾱi (β+1) β!. The parameters β (which is an integer number) andᾱ depend on the observable Φ, on the specific features of the unperturbed system, and on the forcing under consideration. Taking into account (5) and assuming that ω is real, we obtain that χ Φ (−ω)] * , so that Re[χ] is an even function while Im[χ] is odd function of ω. Thus α = α R is real if β is odd, whereas α = iα I is imaginary if β is even.
Taking into account the Titchmarsch theorem, using that χ Φ (−ω)] * , and considering the asymptotic behavior of the susceptibility, it is possible to show that the real and imaginary part of the linear susceptibility obey the following set of general KK dispersion relations [40]: with P indicating integration in principal part and p = 0, . . . , (β − 1)/2 if β is odd and p = 0, . . . , β/2 if β is even. Note that the faster the asymptotic decrease of the susceptibility, the higher the number of independent constraints due to KK relations it has to unavoidably obeys. As thoroughly discussed in [9], in the case of quasi-equilibrium system, the fluctuation-dissipation theorem ensures that the imaginary part of the susceptibility describing the response of a given observable to a perturbation is proportional to a suitably defined power spectrum in the unperturbed system. Therefore, observing the unperturbed system and using Eq.(10) it is possible to reconstruct the entire linear susceptibility, and so know everything about the response properties of the system. In the case of a non-equilibrium system, as discussed in the Introduction, this procedure is not possible.
It is possible to use the KK relations to define specific self-consistency properties of the real and imaginary part of the susceptibility. We first consider the following application: we set p = 0 in Eq. (10) and take the limit ω → 0. We obtain that for any observable: which says that the static susceptibility (i.e., in a more common language, the linear sensitivity of the system) is related to the out-of-phase response of the system at all frequencies. In other terms, Eq. (11) is an exact formula for the linear susceptibility of the system. Note that the static susceptibility is a real number because, thanks to the symmetry properties discussed above, Im[χ (1) (0)] = 0. Moreover, we know that Re[χ (1) (0)] is finite because the susceptibility function is analytic (and so in particular non singular). This is consistent with the fact that as ω → 0, the imaginary part of the susceptibility goes to zero at least as fast as a linear function, as only odd positive integer exponents can appear in its Taylor expansion around ω = 0), so that the integrand in Eq. (11) is not singular. Similarly, we obtain that for ω ∼ 0 the real part of the susceptibility is in general of the form c 1 + c 2 ω. 2 + o(ω 3 ), where c 1 and c 2 are two constants and c 1 is exactly given by Eq. (11). By exploring the ω → ∞ limit in Eqs. (9)-(10) we obtain further integral constraints. By applying the superconvergence theorem [54], we obtain the following set of vanishing sum rules (see [39]): Note that if β = 0 no vanishing sum rules can be written for the susceptibility, whereas if β = 1 only Eq. (13) provides a zero-sum constraint. For each set of KK relations, an additional, non-vanishing sum rule can be obtained. If β is odd, the non-vanishing sum rule is: while if β is even, we have: where the α constants are defined in Eq. (8). These sum rules provide additional general constraints that must be obeyed by any system and can be used to test the quality of the output of any model wishing to describe it. If we are not in the perfect model scenario (e.g., we use a simplified representation of some degrees of freedom) the sum rules can in principle be used to provide a fit for the parametrization. We underline that it is possible to generalize the KK theory for specific classes of nonlinear susceptibilities for both quasi-equilibrium and non-equilibrium systems. Such results, which are particularly suited for studying the fundamental properties of harmonic generation processes, are thoroughly discussed in [17] and will not be reported here.
2.3 A practical formula for the linear susceptibility and consistency relations between susceptibilities of different observables As discussed above, the definition of the linear susceptibility does not depend on the function f (t) modulating the additional forcing, so that it is possible to draw general conclusions on its properties even by choosing a specific function f (t).
Let's consider f (t) = 2ǫ cos(ωt). The impact of the perturbation on the evolution of a general observable Φ(x) is defined as: where x 0 and t 0 are the initial condition and the initial time, and we associate the lower index ǫ to the strength of the forcing. The Ruelle's response theory ensures that δΦ ǫ (t, t 0 , x 0 ) = O(ǫ). Following [41,43] the linear susceptibility results to be: where: is the total susceptibility, affected by the finite time and finite size response of the system. This quantity depends on the initial condition and in principle contains information about the response of the system at all order of nonlinearity. Since d(δΦ ǫ (t, x 0 ))/dt = δ(dΦ ǫ (t, x 0 ))/dt, thanks to the linearity of the time derivative, by considering Eqs. (17)- (18) and performing an integration by parts, we obtain that [41]: Let's now find a different expression for χ where Γ = F · ∇Φ. Similarly, the time derivative for Φ in the case of the perturbed where Ξ = X · ∇Φ. From Eqs. (20)- (21) we obtain that where all terms are of O(ǫ). Furthermore, we integrate each term in Eq. (22) as in Eq. (18), take the limits ǫ → 0 and T → ∞, and obtain: Using the definition in Eq. (17) and the identity given in Eq. (19), we derive: The first limit in Eq. (24) gives, by definition, Ξ 0 , whereas the second limit vanishes as the expression under integral is O(ǫ 2 ), since it is related to second order harmonic generation nonlinear process [43]. Concluding, we obtain the following general consistency relation for the linear susceptibility: Such an identity related the susceptibility of an observable Φ to the susceptibility of the projection of its gradient along the unperturbed flow Γ and to the average value in the unperturbed state of the projection of its gradient along the perturbation flow.
Note that the two terms on the right hand side are radically different. Whereas the first term is related to the projection of the dynamics along the unstable manifold, the second term depends on the structure of the forcing X(x), which may be entirely unrelated to that of the unstable manifold. This is the fundamental reason why the fluctuation-dissipation theorem does not apply in the non-equilibrium case. Moreover, since, as shown in Eq. (7)-(8), the susceptibility of a generic observable decreases to zero at least as fast as ω −1 , for large values of ω we have that χ If ω Ξ 0 = 0, we also have that the leading order of the short-time expansion of the Green function is of the form: in agreement with what can be found by direct inspection of Eq. (3).
3 Application of the Response Theory to the Lorenz 96 Model

Statistical properties of the unperturbed Lorenz 96 Model
The Lorenz 96 model [45,46,47] describes the evolution of a generic atmospheric variable defined in N equally spaced grid points along a latitudinal circle and provides a simple, unrealistic but conceptually satisfying representation of some basic atmospheric processes, even if such one-dimensional model it cannot be derived ab-initio from any dynamic equation via subsequent approximations. The evolution equations can be written in a scaled form as follows: where i = 1, 2, ....., N, and the index i is cyclic so that x i+N = x i−N = x i . The quadratic term in the equations simulates advection, the linear one represents thermal or mechanical damping and the constant one is an external forcing. The evolution equations are invariant under i → i + 1, so that the dynamics is the same for all variable. The time scale of the system is given by the damping time, which corresponds to five days. The L96 system shows different features, as different choices of F and N may strongly alter the topology of the attractor, alternating periodic, quasi-periodic and chaotic behavior in a non trivial way. However with a suitable choice of the parameters N and F , the system is markedly chaotic. In particular, as F controls the energy input into the system, we expect that for relatively high values of this parameter the system should simulate a turbulent behavior and live on a strange attractor. As an example, setting N = 40 and F = 8 the system possesses 13 positive Lyapunov exponents, the largest corresponding to a doubling time of 2.1 days, while the fractal dimension of the attractor [21] is about 27.1 [46]. When computing the time derivative of the total energy of the system, defined as E = 1/2 i x 2 i , the advection terms cancel. The evolution equation for E results to be: As the dynamics takes place inside a compact set, Ψ(x i , . . . , x N ) is bounded for any choice of the function Ψ. Therefore the ensemble mean with respect the ρ 0 (or time average) of the temporal derivativeΨ vanishes. Therefore, defining M = i x i as the total momentum of the system, we obtain the following identity: Similarly, we can deduce an additional consistency relation by investigating the expression of the time derivative of M: where where the choice of i is arbitrary and the N and F -dependence is dropped for shortness.
Defining c i = C i 0 /N, c 0 = 2e, we can rewrite Eqs. (29)-(30) as follows: The F -dependence of e and m can be closely approximated in terms of simple power laws. We obtain, within a precision of about 1% in the considered domain, that and, consistently with Eq. (33), with λ ≈ 1.15 and γ ≈ 0.35. Such a smooth dependence of the intensive energy and momentum with respect to the forcing parameter F is indeed in agreement with the hypothesis that the invariant measure is deformed in a very regular fashion not only locally, but over a large range of the parameter's space. Note that, at the fixed point of the system corresponding to a purely zonally symmetric dynamics (x i = F , ∀i) we have m = F and e = F 2 /2. These formulas give much higher values for both m and e than what found with our empirical power laws for the attractor in the chaotic regime. In fact, at such an equilibrium, which is unstable in the parametric range explored here, the energy dissipation is much weaker than in the co-existing chaotic attractor, which corresponds to the case where breaking nonlinear waves and turbulent motions are present. Interestingly, the presence of well-defined scaling laws with respect to the forcing parameters for the energy and momentum of the system with different characteristic exponents in the chaotic regime and in the coexisting unstable equilibrium is in agreement with previous finding recently obtained in a simple baroclinic quasi-geostrophic model [53].

Asymptotic properties of the linear susceptibility
We perturb the L96 model by adding a small perturbation modulated by f (t) = 2ǫ cos(ωt). The resulting evolution equation is: is a generic function of variables x i . We adopt the chaotic hypothesis [23,22] and we follow the theory proposed by Ruelle [10] and discussed in Sect. 2 in order to study the linear response of suitably defined observables to the perturbation. We first propose to study the high-frequency, response by analyzing in detail the asymptotic properties of the resulting susceptibilities. As discussed in section 2, this constitutes a crucial step for constructing the set of applicable KK relations and for computing the value of the sum rules. We consider two different forcing patters X i . In the first case, we apply the perturbation over all the grid points and we choose X i = 1 ∀ i. Given an observable Φ, we refer to the linear Green function and linear susceptibility resulting from this choice of X i as G (1) Φ,N and χ (1) Φ,N , respectively, where the lower index N indicates that the perturbation acts over all the variables x i . We refer to this pattern of forcing as global perturbation. In the second case, we apply the perturbation only on the variable x j of the L96 model, and we choose X 1 = 0∀i = j, X j = 1, . Since all the points are equivalent in the unperturbed case, the choice of j is arbitrary. In this case, when referring to the linear Green function and the linear susceptibility, the lower index 1 substitutes the N, indicating that the perturbation is localized to one point. We refer to this pattern of forcing as local perturbation.

Global perturbation
We consider perturbations with spatial pattern given by X i = 1 ∀ i and analyze the response of the observable E. Following Eq. (3), the linear Green function G (1) E,a (t) can be explicitly written as: where x(t) satisfies the unperturbed evolution Equation (27). Taking into account Eq.
(7) and Dq. (8), in order to obtain the asymptotic behavior of the susceptibility, we need to study the short time behavior of the Green function. Therefore, we express E( x(t)) as a Taylor series about t = 0 considering the unperturbed flow, compute the integral of each coefficient of the t-expansion over ρ 0 , and seek the lowest order non-vanishing term [43]. The first two terms of the Taylor expansion of E in Eq. (38) give: Using Eqs. (7)- (8), the leading terms of the asymptotic behavior of linear susceptibility can be written as: Since the symmetry with respect the index i is valid also in the perturbed case, given our choice of the forcing pattern, the linear susceptibility of the total energy is given the sum of N identical contributions, each corresponding to the susceptibility of the observable ε = 1/2x 2 i for each of the N variables x i of the system. Therefore, it is possible to define an intensive linear susceptibility χ (1) describes the response of the local energy to the external perturbation. In particular in the limit ω → ∞ we have: Equations (40) and (41) imply that the imaginary part dominates the asymptotic behavior of the susceptibility, so that at high frequency the response is shifted by about π/2 with respect the forcing. Observing that the leading term of asymptotic χ is of order ω −1 just one sum rules apply for either susceptibilities. Limiting our attention to the intensive quantity e, by applying Eq. (15) we obtain: Along similar lines, if we select as observable the total momentum M, we derive that the asymptotic behavior of its linear susceptibility is: where we have defined the intensive susceptibility χ (1) µ,N (ω), where µ is the intensive momentum of the system. As in Eq. (41), the asymptotic behavior is determined by the imaginary part of χ, and the real part of the susceptibility provides the following sum rule: We now wish to go back to the general consistency equation for linear susceptibilities given in Eq. (46). Considering that in the perturbed system the time derivative of the total energy of the system can be written as: the general result given in Eq. (25) can be written as follows: since in this case Γ = −2E + F M and Ξ = M. It is easy to check that the asymptotic behavior for the susceptibilities given in Eqs. (40)- (43) is in agreement with Eq. (46), which is valid at all frequencies.

Local perturbation
We now perturb the system in a single grid point. The symmetry of the unperturbed system implies the equivalence of every point of the latitude circle. Indicating with x j the grid point where forcing is exerted, the pattern of the perturbation vector field is X 1 = 0∀i = j, X j = 1. We consider the same modulating monochromatic function f (t) = 2ǫ cos(ωt) as in the previous case. We analyze the asymptotic behavior of the linear susceptibility of the total energy of the system E, of the total momentum of the system M, which are global variables, and of the local variables constituted by the energy E j = 1/2x 2 j and momentum M j = x j of the perturbed grid point and of its immediate neighbors.
Since we are looking at the linear response and the global perturbation is given by N spatially shifted copies of the local perturbation, for any observable Φ of the form Φ(x 1 , . . . , φ,N . In the case of the observable E, it is straightforward to verify the previous identity at least in the asymptotic regimes. In fact, the short time behavior of the Green function describing the response of the E to the local perturbation results to be: so that the asymptotic behavior of the corresponding linear susceptibility susceptibility is: which agrees with what found for the intensive energy response when the global perturbation is applied (see Eq. (41)). The sum rule for the real part of the susceptibility is exactly the same as in what given in (42): Analogously, we obtain that the asymptotic behavior of χ M,ǫ,1 can be written as: with the corresponding sum rule: in perfect agreement with Eqs. (43) and (44), respectively. It is rather interesting to look into local energy observables. Considering the energy E j of the perturbed grid point x j we have that its short term Green function can be written as: Since because the ρ 0 -average of the temporal derivative of any observable vanishes, thanks to the compactness of the attractor, we obtain the following asymptotic behavior for the linear susceptibility Since, by linearity, χ E k ,1 , comparing this result with what obtained in Eq. (48), we note that the susceptibility of the energy at the position of the forcing E j provides the leading asymptotic term to the susceptibility of the total energy E. Consequently, in the high-frequency range χ (1) E,1 , and the two susceptibilities obey the same non-vanishing sum rule, so that: Nevertheless, by comparing Eqs. (48)-(54), we discover that contributions to the second leading order (∝ ω −2 ) in the high frequency range of the susceptibility of the total energy do not come just from the response of the energy at the perturbed grid point.
perturbed grid point but some other point give a contribution of order ω −2 . Therefore, the asymptotic behavior of the real part of χ E j ,1 . The locality of the interaction suggests to look at the energy of the closest neighbors of x j . Because of the asymmetry of the nonlinear terms in the L96 evolution equations, we consider the observable ψ = 1/2(E j+1 + E j+2 + E j−1 ). It is possible to prove that: It is easy to observe that the sum of χ E loc ,1 , where E loc = E j + ψ = 1/2x 2 j + 1/2x 2 j+1 + 1/2x 2 j+2 + 1/2x 2 j−1 is the energy of the cluster of points centered in x j .
The analysis of the asymptotic behavior of the susceptibilities related to the local momentum of the system provides additional insights. It is possible to prove that for large frequency the linear susceptibility of the momentum of the perturbed grid point is: which suggests that the response of the local momentum captures the correct asymptotic behavior of both the real and the imaginary part of the total momentum M.
Concluding, we obtain that the following sum rule can be stablished: Therefore, such a constraint is exactly the same whether we analyze the the response of the momentum of a single variable when the perturbation acts over all the grid points, or of the total momentum in the case of a local perturbation, or, in this latter case, of the momentum of the grid point where the local perturbation is applied.
As we have seen in this section, the coefficients of the leading asymptotic terms and the sum rules are given by simple linear functions of m (or equivalently, thanks to Eq. (33), by e) and by F . As we have proposed an efficient parameterization of m and e as functions of F alone in subsection 3.1, our results can be easily applied and numerically verified for a very large class of L96 models.

Simulations and Data Processing
The accurate calculation of the linear susceptibility of the general observable Φ is not as easy task, since the definition given in Eq. (17) requires the evaluation of two limits, whereas we can actually compute only the quantity given in Eq. (18). Averaging the response over a long time T allows for improving the signal-to-noise ratio. Noise is present because, due to the chaotic nature of the flow, we have a continuous spectral background. Instead, considering small values for the perturbation strength ǫ degrades the signal-to-noise ration, but, on the other hand, it is crucial to select a small ǫ in order to keep the perturbations as close as possible to the linear regime. As discussed in [43], we can improve the signal-to-noise ratio without needing to perform very long integrations and to consider large values for ǫ by performing an ergodic averaging of the quantity averaging the quantity χ (1) Φ (ω, x i , ǫ, T ). Therefore, we choose the best estimator of the true susceptibility χ (1) φ as: where the x i are randomly selected initial conditions chosen on the attractor of the unperturbed system. The numerical integrations of the Lorenz 96 system have been performed using the standard configuration where N, the number of degrees of freedom, is set to 40, and F , the intensity of the unperturbed forcing, is set to 8 [47,45]. Equations (27) and (37) are solved using the standard fourth order Runge-Kutta numerical scheme.
For a given observable Φ, the susceptibility at angular frequency ω is computed by applying Eq. (59) to K outputs of Eq. (37), each starting with a different initial condition, where the perturbation has the same angular frequency ω. The angular frequency ω ranges from ω l = 0.2π to ω h = 20π with steps of 0.01π. Each simulation performed with a perturbation of angular frequency ω runs from t = 0 up to t = T = 400π/ω, which corresponds to 200 full periods of the forcing. The length of the simulations depends on the corresponding period of the forcing because we are interested in obtaining a frequency-independent quality for the signal. We have observed that the linear response approximation is obeyed to a good degree of approximation for up to ǫ ≈ 1, which implies that the third order nonlinear effects are relatively small. See [17,43] for further clarifications on this latter point.
When considering the susceptibilities describing the response to the global perturbation, we present results obtained using ǫ = 0.25 and averaging over K = 100 random initial conditions. When assessing the linear response to the local perturbation, a reasonably clear signal is obtained using ǫ = 1 and averaging over K = 300 initial conditions.
Note that, since we are interested in the linear response, it is could have been possible to compute the susceptibility using a generic modulating function f (t) (see Eq. (4)) rather than having to resort to multiple monochromatic perturbations. Nevertheless, for reasons of clarity, and for emphasizing that chaotic dynamical systems can be analyzed using tools typical of spectroscopy, we have used a more cumbersome but probably more convincing approach.
We underline that the numerical results have been obtained using a commercial laptop rather than resorting to HPC. This comes from the motivation of showing that the methodology presented is robust enough that relatively low-end means allow us to see the physical and mathematical properties of our interest. We emphasize that, using HPC, it is rather easy to greatly increase the quality of the signal by increasing K and/or T by a one or two orders of magnitude.

Global Perturbation
We first consider χ (1) ε,N = 1/Nχ (1) E,N , where ε = E/N, and follow up from the discussion in subsection 3.2.1. The measured real and imaginary parts of the susceptibility are depicted with the black lines in Fig. (1)a,b. The imaginary part has a broad spectral feature (with two distinct peaks) spanning from ω ≈ 2 to ω ≈ 4, which corresponds to about twice the time scale (=1) of the system and to four times (see Eq. (28)) the relaxation time of the energy. This hints at the fact that it is not obvious to constrain the spectral features of the response an observable just by performing a scale analysis of its evolution equation. For higher values of ω, the imaginary part decreases in a very regular way, so that in the upper range a very good agreement with the asymptotic behavior ∼ m/ω presented in Eq. (8) is obtained. For low frequencies, the imaginary part appears to decrease towards zero, as expected from symmetry reasons.
Instead, the real part presents a dispersive structure in correspondence with the broad maximum of the imaginary part, and changes sign for ω ≈ 6, so that it is negative for high values of the frequency range. The asymptotic decrease to zero in this range is also in excellent agreement with the estimate ∼ −(F − 2m)/ω 2 given in Eq. (8), whereas for low frequencies the real susceptibility tends to a very high value, this suggesting that the strongest response is obtained for static perturbations.
The measured real and imaginary parts of χ (1) µ,N = 1/Nχ (1) M,N , where µ = M/N, are depicted in black in Fig. (2)a,b. Interestingly, the spectral feature of the imaginary part is shifted to higher frequencies than in the case of the energy susceptibility, so that a well-distinct peak centered on value of ω ≈ 6, which approximately corresponds to the natural time scale of the system. For low frequencies, the susceptibility has almost exclusively a real component. As opposed to the previous case, the largest value for the in-phase response is not obtained for ultralow frequencies, but rather for ω ≈ 4. The asymptotic behavior of both the real and imaginary parts is in perfect agreement with the theoretical result given in Eq. (32), so that they are found to decrease asymptotically for high frequencies as 1/ω 2 and 1/ω, respectively.
We apply the truncated KK relations to the measured data to test the quality of the data inversion process. The estimates of the imaginary part (starting from the measured data of the real part) and of the real part (starting from the measured data of the imaginary part) obtained by applying Eqs. (9)-(10) are shown for χ (1) ε,N in blue in Fig. (1)a,b and for χ (1) µ,N in Fig. (2)a,b. We observe that whereas agreement is very good for the real part for both susceptibilities, only a qualitative match is obtained for the imaginary part, with large discrepancies for ω 2. In this latter case, moreover, the well-known problem of KK divergence at the boundaries of integration [39,40] is very serious for ω = ω l .
It is crucial to test whether the discrepancies are due to the finiteness of the spectral range or are, instead, due to basic problems in the applicability of the Ruelle response theory, related to the fact that the invariant probability measure of the unperturbed system actually features large deviations from an SRB measure.
We proceed testing the first case. In order to widen the spectral range over which the susceptibility is defined, we will exploit the asymptotic properties obtained in section 3.2 as well as the low frequency behavior of the susceptibility discussed in section 2.
We redefine the the imaginary part of the susceptibility of χ (1) ε,N as follows where the measured data are sandwiched between the low and high frequency limit. Whereas we have a rigorous result for the high frequency limit, the low frequency limit is computed by making the reasonable assumption that the leading order of the ω → 0 limit is linear (see discussion after Eq. (11)). Similarly, the real part of the susceptibility χ (1) ε,N can be redefined as follows: where we have used the fact that at low frequencies the real part of the susceptibility is constant in ω up to a quadratic term. A corresponding procedure is used to extend the spectral range of the χ (1) µ,N (ω), where the suitable asymptotic behaviors described in subsection 3.2.1 are adopted. The red lines in Figs. (1)a,b-(2a,b) present the results of such extrapolations, and the magenta lines show the outcome of the data inversion of these functions performed via KK relations. We observe that the agreement is outstanding, with almost perfect overlap inside the region where measurement is performed and remarkable agreement also in the low and high frequency range. This is a very convincing evidence that the Ruelle response theory can be successfully applied for this system. Since the KK relations provide, first and foremost, consistency tests, the agreement the original and the KK-transformed susceptibility automatically confirms that the extrapolation procedure we have adopted is correct. A still better agreement would be found had we taken into account value of ω larger than what considered in the extrapolation used here (up to 100 π).
Furthermore, let's consider the results presented in subsection 3.1. The slopes of the functions e(F ) and m(F ) are given by First, we test the sum rules given in Eq. (42) and (44) for the real part of the extrapolated susceptibilities χ (1) ε,N (ω) and χ (1) µ,N (ω), respectively. Our findings are presented in Fig. (3, where it is shown that an excellent agreement (within 1%) is found between the theoretical values and the numerical results. Since Re[χ (1) e,ǫ,a ](ω) is negative in the high-frequency range, the convergence of the integral to the theoretical value of the sum rule is from above, whereas the opposite occurs for Re[χ (1) m,ǫ,a (ω)]. Extending the integral for even larger values of ω would bring the numerical results to an almost perfect agreement with the theory.
Following the definition given in Eq. (3), the Green function G can be used to compute the time-dependent linearized impact of all perturbations with the same spatial pattern X i (x) but with arbitrary time modulation. Whereas the direct estimate of the Green function from the time dependent dynamics can be obtained by performing an ensemble of simulations where the time modulation of the perturbation is given by a δ(t) pattern (see discussion in section 2, we take the indirect route by considering Eq. (5). By applying the inverse Fourier Transform, we derive the Green functions corresponding to χ (1) ε,a (ω) and χ  41) and (43)).

Local Perturbation
The data obtained for the numerical simulations of the response to the local perturbation are, given the much weaker overall strength of the forcing, much noisier that those presented in the previous section. Nevertheless, we shall see that all the theoretical predictions are verified to a surprisingly good degree of approximation.
The global observables E and M are of the form Φ(x 1 , . . . , where φ(x) = x 2 /2 and Φ(x) = x, respectively. We have consistently verified that the identity χ  (1) φ,N discussed in the previous section applies in the whole spectral range explored by our simulations, compatibly with the (slightly) different signal-tonoise ratios in the two sets of simulations. See, e.g., Fig. 5 for the comparison between the two susceptibilities χ (1) E,1 and χ (1) ε,N . We then proceed to analyze more in detail the linear susceptibilities related to local observables. In Fig. (6) we present our results concerning the real and imaginary part of χ (1) E j ,1 (ω). Analogously to what observed in the previous subsection, we have that once the measured susceptibility is extrapolated using the theoretical results obtained via response theory and KK relations, we have an excellent agreement between the original real and imaginary parts and those obtained using the KK inversion. The KK algorithm, instead, provides only a partially satisfying outcome when only data from the measured range are considered. Relatively discrepancies are found near the boundaries of the data range, with an especially serious bias near ω l for the imaginary part of the susceptibility.
When comparing χ (1) E j ,1 (ω) and χ (1) E,1 (ω) (see Fig. 5), we observe that for low frequency the response of the energy at the grid point where the perturbation is applied accounts for about half of the response of the total energy, thus implying that the remaining half is redistributed among the remaining N − 1 grid points. The relevance of the response of grid points other than the directly perturbed one also explains why the peak of Imχ (1) E,1 (ω) (and so of Imχ (1) ε,N (ω)) than that of Imχ (1) E j ,1 (ω) -see the frequency range 2 ω 4. Slower perturbations allow other grid points x k = x j to respond effectively.
Instead, since the leading asymptotic order of χ (1) E j ,1 (ω) and χ (1) E,1 (ω) is the same, at high frequencies the local energy response accounts for most of the energy response of the whole system. In this case, the incoming perturbation is so fast that the internal time scales of the system as bypassed, and mainly a local effect is observed. Nevertheless, the second leading order of the asymptotic expansion of the two susceptibilities has opposite sign (see Eqs. (41) and (54)), which suggests that at any large but finite frequency the local energy response is only a good approximation to the response of the total energy. The changeover between the two regimes occurs around the frequency of the peak of Imχ (1) E j ,1 (ω), which corresponds to a perturbation with period close to 1. Thanks to the asymptotic equivalence between χ (1) E j ,1 (ω) and χ (1) E,1 (ω) (and χ (1) ε,N (ω)), they must obey the same sum for the real part of the susceptibility (see Eqs. (42)-(55)), even if the two real parts, as discussed above, are rather different in value in the low frequency range and even in sign in the high frequency range. Figure 8 confirms that this rather counter-intuitive behavior is actually observed. Note also that sum rules, resulting from an integration, are less sensitive to noise in the data, but this occurs if and only if the underlying signal is correct. Therefore, we understand that in χ (1) E,1 (ω) and χ (1) ε,N (ω) the strong static and quasi-static response and the (rather odd) negative sign for high frequencies of the real part of the linear susceptibility, which are crucially related to the behavior for the grid points different from the perturbed one) compensate each other to guarantee agreement with the sum rule obtained from the real part of χ (1) E j ,1 (ω), which instead has a smaller range and more regular (monotonic) behavior with frequency.
A formally similar -and analogously spectacular -spectral compensation has been observed in a physical process as different from what we are analyzing here as the electromagnetically induced transparency [55]. The result obtained here supports pre-vious findings obtained on quasi-equilibrium systems suggesting that sum rules do not depend on many-particle interactions [39,40].
The investigation of the linear susceptibility of the variable x j is not as insightful as that of E j . We find that linear susceptibility χ (1) M,1 (ω) (and χ (1) µ,N (ω), see Fig. 2) in both the real and imaginary parts at all frequency. The only notable differences are that the static response x j is slightly larger than than of M, and that the imaginary features a secondary peak at slightly larger frequencies than the main spectral feature. We have verified,as in the previous cases, the results of the numerical simulations accurately agree with the theoretical results regarding the asymptotic behavior of both the real and imaginary part and that KK relations map to high degree of precision the real and the imaginary parts into each other. See Fig.  7 for details.
We present as main finding of the analysis of the observable x j that, as predicted by the theory, the real part of χ (1) x j ,1 (ω) obeys the same sum rule as the real part of χ (1) M,1 (ω) or of χ (1) µ,N (ω), because the corresponding imaginary parts feature the same asymptotic behavior. Figure 8 shows that in the case of the momentum variables the cumulative integral is rather similar for the susceptibility of the local and of the global variable, with small discrepancies in the region around the peak of the response.

Further implications of Kramers-Kronig relations and sum rules
We now show how the knowledge of the asymptotic behavior of the real and imaginary part and the knowledge of the validity of the KK relations and related sum rules allow to draw general conclusions on the similarities and differences between two given linear susceptibility functions. Let's consider the case that these two susceptibilities feature the same first order asymptotic expansion in the high frequency limit. Let's assume that it is an odd power of ω, so that the real part is negligible for high frequencies.
Therefore, the two susceptibilities will obey the same sum rule for, e.g. the 0 th moment of the real part.
If they agree also in the asymptotic behavior of the real part, they cannot feature large discrepancies in the low frequency range of the real part of the susceptibility either, or otherwise the agreement of the sum rules would be broken. Therefore, the real part of the two susceptibilities are similar, and, as a consequence of the KK relations, the two imaginary parts will also be similar. If, instead, there is a discrepancy in the asymptotic behavior of the real part of the two susceptibilities, the two real parts will necessarily be rather different in the low frequency range, again in order to comply with the sum rule constraint. As the two real parts are different, the imaginary part of the two susceptibility will also be rather different, except, from hypothesis, in the high-frequency range. The first scenario envisioned here pertains to the pair of linear susceptibilities x j ,1 (ω) and χ E,1 (ω). Note that, taking into account the asymptotic properties of the susceptibility of the observable E loc = 1/2x 2 j + 1/2x 2 j+1 + 1/2x 2 j+2 + 1/2x 2 j−1 , discussed in the previous section we conclude that χ Obviously, a similar argument applies if the leading order is real. This discussion further clarifies that the higher the number of independent KK relations and related sum rules verified by a susceptibility functions, the more stringent are the constraints on its properties.

Additional Properties of the Linear Susceptibility
The special mathematical properties of the linear susceptibilities allow to investigate further properties of the response. In particular, we note that for m ≥ 1 the function [χ (1) Φ ] m is analytic in the upper complex ω-plane just as as χ (1) Φ . This allows, as discussed in [40] to derive, in principle, an infinite set of integral relations (KK and sum rules) deriving just from the holomorphic proprieties of the susceptibility. As an example, we have considered the square of the linear susceptibility [χ (1) ε,N (ω)] 2 . From Eq. (41), it is easy to prove that the following asymptotic expansion holds for large values of ω: As shown in the first panel of Fig. 9, KK relations are found to connect up to a high degree of approximation the real and imaginary part of [χ (1) ε,N (ω)] 2 . Moreover, thanks to the asymptotic behavior given in Eq. (64), it is possible to establish, thanks to Eqs.
(13)- (14), the following sum rules: The second panel of Fig. 9 shows that the obtained numerical results are in excellent agreement with the theoretical predictions. Note that these results do not have an obvious physical interpretation, as the inverse Fourier Transform of [χ (1) ε,N (ω)] 2 is given by the convolution product of the Green function G (1) ε,N (t) with itself, while they depend only on the formal properties of the linear susceptibility.

Practical Implications for Climate Change Studies
In this paper we have constructed and verified to a high degree of accuracy the linear response theory for a simple yet prototypical climate model by computing the frequency-dependent susceptibilities of several relevant observables related to localized and global patterns of forcings. These results pave the way for devising a rigorous methodology to be used by climate models of any degree of complexity for studying climate change at, in principle, all time scales using only a very limited set of experiments, and for exploiting effectively the currently adopted ensemble runs methods.
Let's consider, for sake of simplicity, that the observable Φ is the time-dependent globally averaged surface temperature of the planet T S , that F (x) represents the whole set of climate equations in a baseline scenario (e.g., with pre-industrial CO 2 concentration), and that the perturbation field X(x) is nothing but a constant field of CO 2 concentration, which directly impacts only the radiative part of the code. The perturbation is modulated by a time-dependent function f (t) to be specified below. We assume, for simplicity, that the model does not feature daily or seasonal variations in the radiative input at the top of the atmosphere.
From Eq. (2), we have that In practical terms, the left hand side of this equation is nothing but the ensemble average of the time series of the change between the globally averaged surface temperature of the planet at a time t after the perturbation has started. Note that the direct estimate of G T S (σ) is likely to be overwhelmingly difficult. Using Eq. (4), we have that: which implies that once we compute the Fourier Transform of the time series mentioned above and we know the modulating function f (t) (and so its Fourier Transform f (ω)), we can reconstruct χ T S (ω). Let's select a particularly simple example of modulating function f (t) = ǫ(Θ(t) − Θ(t − τ )). This is just a rectangular function of width τ , of height ǫ, and shifted from the origin by a forward time translation τ /2. In practical terms, this corresponds to changing abruptly the field CO 2 concentration by ǫ at time t = 0 and taking it back to its original value at t = τ . we then obtain: ǫ(sin(ωτ ) + i(1 − cos(ωτ ))) .
Once we know χ T S (ω), as widely discussed in this paper, we can compute G T S (t), and we know everything about the response of the system at all time scales, including the static response. Note that any choice of f (t) is equally valid to set up this procedure as long as f (t) is square integrable. This implies that, in a very profound way, the kind of forcing scenarios used in the various assessment Reports of IPCC, where the CO 2 concentration typically stabilizes at a different value from the preindustrial one (so that f (t) does not tend to 0 as t → ∞) are not necessarily the only nor the best ones, in spite of what could be intuitively guessed, to study even the steady state response of the system.
Obviously, a similar set of experiments could be devised for studying rather thoroughly the response of the climate system to a variety of forcings, such as changes in the O 3 concentration, aerosols, solar radiance, as well as to changes in the parameteri-zations. In the case of uncoupled models of one subdomain of the climate system (e.g. atmospheric and oceanic GCMs, land-surface models), this strategy could be used to study the impact of perturbations to the boundary conditions provided by the other subdomains of the climate system.

Summary, Discussion and Conclusions
The climate can be seen as a complex, non-equilibrium system, which generates entropy by irreversible processes, transforms moist static energy into mechanical energy [1,2] as if it were a heat engine [3,4], and, when the external and internal parameters have fixed values, achieves a steady state by balancing the input and output of energy and entropy with the surrounding environment [56,4]. For such basic reasons, the tool of equilibrium and quasi-equilibrium statistical mechanics cannot provide suitable tools for studying the fundamental properties of the climate system. In particular, the fluctuation-dissipation theorem, which allows for deriving the properties of the response of the system to external perturbations from the observations of its internal variability cannot be applied.
It is reasonable to ask whether is possible to evaluate how far from equilibrium the climate system actually is. It is possible to evaluate such distance in a mathematically sound way by assessing the ratio of the dimensionality of the attractor of the system over the total number of degrees of freedom. Whereas a ratio close to one indicates that only small deviations from equilibrium are present, a small ratio suggests that strongly non-equilibrium conditions are established. See [57] for a detailed treatment of this problem in the classical case of heat conduction. In the case of a quasi-geostrophic atmospheric model forced by Earth-like boundary conditions, the dimensionality of the attractor of the model is about one order of magnitude smaller than the total number of degrees of freedom [58]. While not conclusive, this seems to suggest that the best framework to interpret the climate is that of a far from equilibrium system. Following either explicitly or implicitly the programme of the Catastrophe theory [59], many authors have approached the problem of understanding the fundamental properties of the climate system by looking at the detailed structure of the bifurcations of the deterministic dynamical system constructed heuristically in order to represent the dynamics of the main climate modes using as few degrees of freedom as possible. Such an approach often hardly allows to efficiently represent the fluctuations and the statistical properties of the system. The introduction of stochastic forcing provides a relatively simple but conceptually rich partial solution to some of these draw-backs, even if the Hasselmann programme [60] suffers from the need for a -usually beyond reachclosure theory for the properties of noise. Therefore, the stochastic component is usually introduced ad hoc, with the ensuing lack of universality and/or robustness when various levels of truncations are considered. These strategies have anyway brought to outstanding scientific results and has been suggested the existence of generic mathematical structures present in hierarchies of CMs [61]. Recently, the unified treatment of chaotic and stochastic dynamics using the results of the mathematical theory of random dynamical systems is emerging as new, promising paradigm for the investigation of the structural properties of the climate system [62].
We have proposed a different perspective. In agreement with the view given above, we have taken as mathematical framework for the analysis of the climate system that of non-equilibrium statistical mechanics, and have focused on the steady state properties of ergodic dynamical systems [20] possessing the special property of having an invariant measure of the SRB type [21]. As proposed by the chaotic hypothesis [23,22], this mathematical framework is well suited for analyzing general non-equilibrium physical systems.
In this context, the impact on the system of general perturbation can be treated using the response theory recently introduced by Ruelle [10,25,16], which allows to compute the change in the expectation value of a generic observable as a perturbative series where each term is given by the average over the unperturbed invariant measure of a function of the phase space which depends on the considered observable and on the applied perturbation. In other terms, even if the internal dynamics of the system is nonlinear and chaotic, the leading order of the response is in general linear with the strength of the added perturbation. This approach overcomes the difficulties related to the singularity of the invariant measure discussed in [63].
At each order, the propagator of the perturbation, i.e. the Green function, is causal. This allows for applying dispersion theory and establish general integral constraints -KK relations -connecting the real and imaginary parts of the susceptibility, i.e. the Fourier Transform of the Green function [17,43] In this paper we have first recapitulated the main aspects of the general response theory and have propose some new general results, which boil down to consistency relations between the linear susceptibilities of different observables. The obtained equation provides the basic idea why the fluctuation-dissipation theorem does not apply in nonequilibrium cases. We have showed for the first time that the Ruelle linear response theory can be applied with great success to analyze the climatic response to general forcings. We have chosen as test bed the Lorenz 96 model [45,46,47], which, in spite of its simplicity, has a well-recognized prototypical value as it is a spatially extended one-dimensional model and presents the basic ingredients, such as dissipation, advection and the presence of an external forcing, of the actual atmosphere. Such a model features a different level of complexity with respect to those adopted in previous numerical investigations of Ruelle's theory [41,42,43] We have analyzed the frequency dependence of the response of the local and global energy and momentum of the system to perturbations having a global spatial pattern and to perturbations acting only on one grid point. We have derived analytically several properties of the corresponding susceptibilities, such as asymptotic behavior, validity of KK relations, and sum rules. We have shown that all the coefficients of the leading asymptotic expansions as well as the integral constraints can be written as linear functions of parameters that describe unperturbed properties of the system, and in particular its average energy and average momentum. The theory has been used to explain differences in the response of local and global observables, in defining the intensive properties of the system and in generalizing the concept of climate sensitivity to all time scales.
We have then verified the theoretical predictions from the outputs of the simulations up to a high degree of precision, even if we have used rather modest computational resources (a total of about 30 cpu days of a mid-range commercial laptop). We have verified that the linear response theory holds for perturbations of intensity accounting to up to about 10% of the unperturbed forcing terms. Even when local perturbation and local observables are considered it is possible to achieve a signal-to-noise ratio which permits rather satisfactory comparisons with the theory. We have proved that the combined use of KK relations and the knowledge of the asymptotic behavior of the susceptibilities allows for extrapolating in a rigorous way the observed data. We also have shown how to reconstruct the linear Green function, which can be used to map perturbations of general time modulation into changes in the expectation value of the considered observable for finite as well as infinite time.
Our numerical experiments have been performed using one of the standard settings of the Lorenz 96 model, namely the version identified by having N = 40 degrees of freedom and forcing F = 8. Nevertheless, some newly obtained empirical closure equations expressing the average energy and the average momentum of the unperturbed system as simple power laws of F (with no dependence on N) have allowed to extend our results to the entire class of chaotic Lorenz 96 models.
In this paper we have only used the KK relations in the most simplistic framework, i.e., computing the KK transforms and evaluating their agreement with the original data. Actually, several more sophisticated analysis techniques are available, such as recursive self-consistent algorithms, where the measured data are taken as first guess, exploiting the fact that multiple applications of KK relations, combined with sum rules, automatically filter our the noise and remove most of the spurious signal [40].
We believe that the proposed approach, which we may dub as spectroscopy of the climate system, may constitute a mathematically rigorous and practically very effective way to approach the problem of evaluating climate sensitivity and climate change from a radically new perspective. In this regard, we have proposed a rigorous way to compute the surface temperature response to changes in the the CO 2 concentration at all time scales using only a specific set of simulations, and taking advantage of the theoretical results presented here. We underline that our approach takes into account all the (linear and nonlinear) feedbacks of the system, as they are included in the definition of the Green function. This, at a very practical level, is the great advantage of using Ruelle's formulas. At a more basic level, whereas considering more complex models requires heavier computational resources, the modest cost of the present set of simulations suggests that, at least for global or regional climatic observables, it is feasible to test the theory discussed here for simplified yet Earth-like climate models without resorting to top-notch computing facilities. Moreover, while in this paper we have computed the susceptibilities using, on purpose, a very cumbersome method, more efficient strategies can be devised, at least when the linear regime of the response is considered. Apart from the practical example given for the case of the impact of the CO 2 concentration, these include studying the response of the system to δ(t) like perturbations, which gives directly the Green function of the system, and including in the forcing various monochromatic signals. Of course, in all cases, a Monte Carlo approach is needed in order to sample effectively the attractor of the unperturbed system in terms of the initial conditions of the simulations.
These results pave the way for future investigations aimed at improving and extending the theoretical framework presented here, at finding results of general applicability in the context of the modelling of geophysical fluid dynamics, and, finally at answering specific questions of relevance for climate dynamics. In this paper we have analyzed the simple case of the linear response, but, as discussed in [25,17,43], we have the algorithm to compute higher order terms, so that the treatment of the nonlinear response in entirely feasible.
In the first direction, we foresee the possibility of writing out explicitly the linear susceptibility of a general observable by projecting the perturbation onto the unstable, neutral and stable manifolds and analyzing separately the contributions to the total response. This will probably require the adoption of adjoint techniques. Moreover, we will be testing the radius of convergence of the Ruelle response theory is some specific examples.
Along the second direction, we propose to study the impact of stochastic forcing to deterministic chaotic models by treating the (additive or multiplicative) noise as a perturbation to be analyzed using the linear and nonlinear Ruelle response theory and related spectral methods. Moreover, we shall look into the spectral peaks of the susceptibilities and try to understand how the amplification of the response is related to resonances of the system and to the activation of positive feedbacks.
Along the third direction, we envision the analysis of the impact of topography on the statistical properties of the circulation in a quasi-geostrophic setting, thus extending in a climatic perspective what presented in [64]. Moreover, we will tackle in an idealized setting the problem of computing the response of the storm track to changes in the surface temperature [65]. Moreover, we will try to compute along the way discussed here the climate response to changes in CO 2 and solar irradiance using simplified but rather valuable climate models like PLASIM [66].
Finally, we would like to remark that the theory and the practical recipes proposed here could be of direct interest for all projects aimed at auditing climate models' performances and at studying practical problems connected to climate change, such as