## Abstract

A linear response theory of systems of interest in atmospheric and climate dynamics taking fully into account the nonlinearities of the underlying processes is developed. Under the assumption that the source of intrinsic variability can be modeled as a white-noise process, a Fokker–Planck equation approach leads to fluctuation–dissipation-type expressions in the form of time cross-correlation functions, linking the perturbation-induced shift of an observable to the statistical and dynamical properties of the reference system. These expressions feature a generalized potential function and enable one to go systematically beyond the Gaussian approximation usually adopted in the literature. Full solutions and explicit expressions are derived for different subclasses of systems, including a global climate model giving rise to oscillatory behavior.

## 1. Introduction

The response of the atmosphere and climate to external influences—including forcings of anthropogenic origin—is a matter of considerable concern, in connection with the need to produce reliable predictions on their medium- and long-term evolution. Typically, owing to the complexity of the numerical prediction models and the diversity of the possible forcings, the response often appears to be highly system dependent. This makes it difficult to identify basic mechanisms at its origin and, in particular, to relate sensitivity to key properties such as natural variability of the reference system in absence of the forcing.

In a pioneering work, Leith (1975, 1978) suggested that this task could be achieved to the extent that a fluctuation–dissipation-type theorem (FDT) as derived in statistical physics could hold for the atmosphere and climate system. In its best known form (Kubo 1966), this theorem stipulates that in the weak forcing limit the response of a typical set of observables is determined by a linear combination of the forcings present where the proportionality coefficients are integrals over time of autocorrelation or cross-correlation functions (essentially, time-lagged covariances) of these observables. Subsequent investigations focused on the conditions under which such a property could be expected to apply to the atmosphere and climate. Using a truncated version of the barotropic vorticity equation augmented by dissipative terms, Bell (1980) showed that although the FDT in its classic version does not strictly hold, it nevertheless provides a good estimate of the sensitivity of the model. In all these cases, the structure stipulated for the response coefficients involves the cross correlation of the observables of interest and of the state vector itself, which is justified only in a Gaussian approximation. More recently, Majda et al. (2005) presented a general formulation of the fluctuation–dissipation theorem for general classes of complex nonlinear systems along with results testing the validity of the predictions based on the Gaussian approximation for prototypical chaotic systems of interest in atmospheric physics. Further theoretical developments aiming to justify this assumption, complemented with computational algorithms, have been reported (Gritsun and Branstator 2007; Gritsun et al. 2008; Abramov and Majda 2007, 2008; Majda et al. 2010) and validated on models of varying sophistication. Cooper and Haynes (2011) address the non-Gaussian case and show how the FDT can be used to achieve a nonparametric evaluation of the invariant probability density of the reference system. Finally, other investigators (Achatz et al. 2013) assume FDT to be valid and utilize it to carry out parameterizations of subgrid-scale phenomena.

In the present work, FDT is revisited from a different angle. Reference systems in the form of nonlinear dynamical systems subjected to Gaussian Markov noise source terms are considered, which, as argued by many authors (e.g., Imkeller and von Storch 2001), capture a number of aspects of the intrinsic variability of observables of interest in atmospheric and climate dynamics by accounting, in particular, for the effect of subgrid-scale phenomena. We derive systematically general expressions for the response of such systems to weak external forcings in a way that fully accounts for the nonlinearities present in the intrinsic dynamics. We show that the response can be cast into the cross correlation of the vector of observables of interest and the gradient of a generalized potential function, which also determines the structure of the (generally non-Gaussian) probability density of the reference system. We derive a nonlinear differential equation satisfied by this quantity and sketch how, for any given system, this equation can be solved to any desired order by power-series expansion. This formulation leads us to consider a class of dynamical systems referred to as potential (or alternatively integrable) systems, for which the generalized potential generates not only the invariant probability density but also the time evolution of the state variables in absence of noise. Explicit forms of the response functions are obtained for generic subclasses of such systems, in which the role of the nonlinearities and of the strength of the noise can be assessed. In particular, the limitations of the Gaussian approximation in the presence of slowly decaying correlations are brought out. More general classes of systems are subsequently considered, including systems giving rise to instabilities leading to periodic behavior. Linearizing the system’s dynamics in the absence of the forcing around a reference state leads to a Gaussian form of the invariant probability density, whose parameters can be related to those of its intrinsic dynamics. A systematic way to account for nonlinear effects and for non-Gaussian probability densities is also indicated.

The general formulation is laid down in section 2. Sections 3 and 4 are devoted to linear response for potential systems and for general systems, respectively. The main conclusions are summarized in section 5.

## 2. Formulation

Let be the variables involved in a description afforded by an atmospheric or climatic forecast model. As is well known, their evolution displays a high variability, reflected by irregular deviations around some long-term systematic trend. In what follows we account for this variability by decomposing the associated evolution laws into a deterministic and a random part (Hasselmann 1976; Nicolis and Nicolis 1979; Imkeller and von Storch 2001; Arnold et al. 2003; Nicolis 2005):

Here stands for the long-range, systematic effects. It is, typically, a nonlinear function of the state variables and depends in addition on a number of parameters associated to the different constraints acting in the system. The random part stands for the variability around the evolution generated by **v**. This variability may arise from such effects as the influence of the unresolved scales, energy or momentum transfer between different system components, or, if appropriate, for the gross features of deterministic chaos that would be generated in a more comprehensive description at the level of an enlarged set of variables of which constitute a coarse-grained version. It will be decomposed into a random vector and a coupling (generally state dependent) operator in the form of an matrix. In what follows, the components of will be assimilated to uncorrelated Gaussian white noise,

where is the Kronecker symbol, it being understood that the associated variances are small compared to the magnitude of the deterministic contributions generated by **v**.

Our objective in this work as stated in the introduction is to study the response of a reference system to a weak disturbance. To formulate this problem quantitatively, we decompose **v** into an unperturbed part and a (generally state-dependent) term associated to the action of weak forcings switched on at some stage of the evolution owing, for instance, to the variation of some parameter or to the presence of an external field:

with

As is well known, under the above conditions the instantaneous probability density satisfies a Fokker–Planck equation. According to the theory of stochastic differential equations, when white noise is coupled multiplicatively to the system’s state variables [i.e., when in Eq. (1) depends nontrivially on **x**], the Fokker–Planck equation can be written in two different ways according to whether the Ito or the Stratonovich interpretation is adopted. This relates in turn to the way a stochastic integral—that is, an integral over time of a smooth function multiplied by a white-noise process—is defined as a limit of a sum of terms evaluated at discrete times up to the time of interest. In many physical systems, the Stratonovich interpretation is considered more appropriate, as it obeys to the rules of ordinary calculus when changes of variables are carried out. But in what follows, the Ito interpretation is chosen, as it turns out to be the most appropriate when a continuous-time stochastic process is obtained as the limit of a discrete time process. One can then write the Fokker–Planck equation in the form (Gardiner 1983)

or equivalently

where stands for the unperturbed Fokker–Planck operator in the absence of the forcing and accounts for the presence of the forcing,

To evaluate the response to the disturbance, we decompose in a similar way into an unperturbed part and a perturbation ,

with

taking to be the invariant distribution generated by Eq. (5a) in the absence of the forcing :

Substituting Eqs. (6) and (7) into Eq. (5) and limiting ourselves to linear contributions we obtain the following equation for the perturbation:

The formal solution of this equation satisfying the initial condition is

The average response of the system as measured by a representative observable is thus

To arrive at a more explicit and workable expression, we observe that within the framework of a linear response theory and a weak noise limit adopted here the forcing and the coupling operator can be evaluated at the reference state in the absence of both the noise and the forcing, which satisfies the set of algebraic equations

Stated differently, the system is here taken to operate around a steady-state solution of the unperturbed, noise-free equation. In this limit the forcing and the noise become thus additive and Eq. (7) takes the form

where is given by

Equation (12a) has been studied extensively in the literature. One of the main results of these investigations is that in the weak noise limit the solution takes the form (Kubo et al. 1973; Nicolis and Nicolis 2012)

where *Z* is a normalization factor and the extrema of coincide with the solutions of Eq. (11). Furthermore, is an analytic function of its variables as long as the system is not in the immediate vicinity of a bifurcation point in which several branches of solutions of Eq. (11) coalesce. In analogy with equilibrium thermodynamics we will refer to as a “generalized potential.” It should, however, be kept in mind that differs from traditional thermodynamic potentials as the free energy in the sense that it is generated by the full-scale stochastic dynamics.

Taking expressions (12) and (13) into account and recalling that in the setting adopted becomes state independent, one may rewrite Eq. (10) in the form

where the response function is given by the expression

being the adjoint of the unperturbed Fokker–Planck operator .

To give a physical meaning to Eq. (15), we introduce the conditional probability density of the unperturbed system

By definition, for any two observables *A* and *B* the expression

represents then the time cross-correlation function of *A* and *B*,

where the angle brackets represent an ensemble average. Utilizing these expressions along with the definition of the adjoint operator, integrating over , and setting to be the invariant distribution, one can rewrite Eq. (15) as

In summary, from Eqs. (14) and (16), one sees that the linear response along the component *i* of the perturbation can be cast in the form of a cross-correlation function of the observable of interest *h* and a function of the state variables given by the *i*th component of the gradient of . This result, which is reminiscent of the fluctuation–dissipation theorem of statistical mechanics (Kubo 1966), has been rederived in various forms in the context of atmospheric physics (Gritsun and Branstator 2007; Gritsun et al. 2008; Abramov and Majda 2007, 2008; Majda et al. 2010). It should be pointed out that even if the observable *h* is chosen to be one of the state variables , the response function will in general not reduce to a simple autocorrelation function as advanced sometimes in the literature, owing to the presence of nonlinearities and of several coupled variables. As stated in the introduction, the principal difference between our work and previous investigations is that the stochastic dynamics of the reference system is accounted for in its full generality, including nonlinear aspects, and there is no need to impose a Gaussian assumption from the outset. Furthermore, as shown in the sequel, our formulation constitutes a starting point to classify and analyze different situations according to whether the deterministic dynamics derives or not from a potential or satisfies symmetry relationships.

## 3. Potential systems

We shall first apply the formulation developed in the preceding section to the class of systems, where, in the absence of the perturbation and noise, the variables are coupled in a way to satisfy the symmetry property (also referred to as “integrability condition”)

These relations imply in turn that there exists a function generating the unperturbed dynamics through the relation (Thom 1972; Nicolis 1995)

Insofar as is the “flux” associated to the variable , the right-hand side of Eq. (17) can be viewed as a “force” generating this flux. This is reminiscent of irreversible thermodynamics, where the force is the derivative of a thermodynamic potential such as entropy. However, in contrast to thermodynamics, *U* is generated here by the nonlinear dynamical laws responsible for the evolution of . To mark this difference, we shall refer to *U* as the “kinetic potential.” Explicit examples of systems satisfying property (17) will be given later on in this section.

In the above setting, the effective covariance matrix of the noise [Eq. (12b)] will be expressed in the simplified form

The rationale behind this relation is that since “dissipation” (here the deterministic dynamics) is generated by a single mechanism, the “fluctuations” [here the noises and their coupling to the deterministic dynamics] will likewise be governed by a single source. This kind of property linking random forces to systematic effects is referred to as “fluctuation–dissipation theorem of the second kind” (Nicolis 1984a).

which can be integrated straightforwardly to yield

This expression has the form of Eq. (13) in which the generalized potential is replaced by the kinetic potential *U*. The stochastic and the deterministic dynamics are here generated by the same function of the state variables. As a corollary, in Eq. (16) the response function will likewise feature the derivative of the kinetic potential.

In a somewhat more general setting, allowing for couplings between fluxes and forces associated to the different variables of , Eq. (17) is replaced by

Here is a positive-definite symmetric matrix analogous to Onsager’s matrix of phenomenological coefficients governing the coupling between irreversible processes in nonequilibrium thermodynamics (De Groot and Mazur 1962). Equation (12a) admits then still a solution of the form of Eq. (19b), as long as the noise covariance matrix satisfies the relationship

which can again be viewed as a fluctuation–dissipation theorem of the second kind.

### a. Systems involving one variable

Systems involving one variable satisfy automatically condition (17). Here *U* is just the integral over *x*, up to an additive constant, of the right-hand side of the equation of evolution of the variable in absence of both forcing and noise,

The importance and genericity of this class of systems is related to the fact that, under appropriate conditions, multivariate systems can be reduced to a description amenable to a single variable. Two well-known instances are time-scale separation allowing for the adiabatic elimination of fast variables and systems operating around a bifurcation point at a simple real eigenvalue of the Jacobian matrix associated to the evolution vector **v**. The first class includes global energy balance models (Nicolis and Nicolis 1981; North and Cahalan 1981). In the second class, a suitable combination of the original variables referred to as the order parameter is shown to satisfy a universal equation referred to as normal form (Nicolis 1995). In the sequel we will be interested in the normal form associated to the pitchfork bifurcation in which two stable branches of solutions are generated beyond the instability of a reference state, as it happens in thermal convection and other related problems involving symmetry breaking:

with

We will focus on the precritical behavior in which is the unique steady-state solution of the reference system (22).

We now derive the expression for the linear response for this class of systems. Utilizing Eqs. (14)–(16) and (19b) and choosing the state variable *x* to be the observable of interest, we obtain

with

where the integral over *x* runs here from to .

Figure 1a depicts the numerically computed response function [Eq. (25b)] as a function of time for , , and (solid line) and (dashed line). After an essentially exponential decay on a time scale of , the function decreases to zero. Subsequent numerical integration over time and multiplication by *F*^{e} taken to be time independent and equal to 0.01 yields the upper curve of Fig. 1b, depicting the time-dependent response for . As can be seen, after a transient stage, the response reaches a plateau value of about 0.058 in excellent agreement with analytic estimates (dotted line; see below). The agreement subsists for higher values of *F*^{e}, at least up to 0.05. For reference the result in which is identified to the time autocorrelation of *x* multiplied by is provided in the dashed line in the same figure. The difference is substantial, showing the importance of nonlinear effects. However, increasing the distance from criticality (or equivalently reducing the correlation time) by, for example, taking leads to a practical merging of the two curves, indicating that nonlinear effects become negligible. In this range a “naive” version of the fluctuation–dissipation theorem in which features the autocorrelation function of *x*, which amounts to keeping only linear contributions in Eq. (22), provides an adequate description.

Actually, as long as one is limited to the asymptotic response an exact analytic expression of the plateau value in Fig. 1 can be derived. Indeed, the noise-free part of the full dynamics [Eq. (22)], including the forcing, derives also from a kinetic potential *W*:

The invariant distribution in the presence of the forcing is then obtained from a straightforward generalization of Eq. (19),

and the mean response of *x* to the disturbance is given to all orders in *F*^{e} by

Expanding the exponential in powers of *F*^{e} and keeping linear terms, noticing that the zeroth-order term vanishes by symmetry, one is left with

The variance is evaluated in the reference (absence of forcing) system and is given by the analytic expression

where is the Bessel function of the second kind. For the parameter values used in Fig. 1, this gives a value of , indicated by the dotted line in the upper part of the figure, in full agreement with the independently carried out numerical evaluation based on the fluctuation–dissipation formalism.

To disentangle the relative roles of the noise and of the linear and nonlinear parts in , we first observe that in the noise-free limit the steady-state solution of Eq. (23) satisfies the relation

This reduces to in the linearized limit and shows that nonlinearities in the precritical regime tend to depress the level of the response. To incorporate the effect of noise we perform an average of Eq. (23) over the realizations of the noise

and decompose into the sum of the deterministic part (see above) and a fluctuating part . Taking symmetries into account we obtain, to the dominant order in ,

showing that noise also tends to depress the level of the response. This action is manifested via the nonlinear part of . In the linearized limit, the deterministic and the full stochastic responses are thus identical.

### b. Transport processes in multivariate systems

Multivariate systems of interest in atmospheric and climate dynamics belonging to the class of potential systems need necessarily be purely dissipative in the absence of forcing and noise—that is, not to include explicitly inertial effects such as advection or wave propagation. From this point of view, they might seem rather atypical, but they have the merit to allow for extensive analytic insights and for exploring situations beyond the Gaussian approximation. Usually, they describe transport processes such as energy transfer or the diffusion of vorticity to which the general vorticity equation reduces for flows with low Reynolds number. In this subsection, we analyze the linear response of a model system proposed by Paltridge (Paltridge 2001) that can be cast in a potential form. In this model, the earth–atmosphere system is described as a two-box model, each box being associated with the adjoining equatorial and polar regions. Boxes 1 and 2 have temperatures and (), emit an infrared radiant energy and (where *σ* is the Stefan–Boltzmann constant), and exchange a heat flux *J* () directed from region 1 to region 2. Additionally, box 1 receives an incoming flux of solar energy *F*:

where are the random forces. In what follows, *J* will be taken to be either a fixed parameter or proportional to the thermodynamic driving force associated to energy transfer,

In both cases the noise-free parts of Eqs. (31) are generated by a potential in the sense of Eq. (17), with

Adopting form (18) for the covariance matrix of the noise leads then to an invariant distribution of the unperturbed system of the form (19b), where the potential *U* is identified to or depending on the form adopted for *J*. To proceed further, we consider case 1 and decompose *F* into a reference part and a perturbation ,

Since this perturbation appears additively in Eq. (31) the perturbed Fokker–Planck operator is simply

with

or, finally,

We see, as also pointed out earlier, that the response cannot be cast in the form of a correlation function of the variables or , even if *h* is chosen to be one of these variables. On the other hand, in Eq. (33b) is just the systematic part of the overall flux contributing to the change of in time [first formula in Eq. (31)]. Choosing the observable *h* to be this flux, one may then cast Eq. (33b) in the form of an autocorrelation function,

The analysis can be extended straightforwardly to case 2, with the associated potential [Eq. (32b)].

Figure 2a depicts the numerically evaluated autocorrelation function [Eq. (34)] for parameter values , , and . After an essentially exponential decay with a characteristic time of about one time unit, the function eventually oscillates irregularly around zero. Subsequent numerical integration over time and multiplication by yields then the response represented by the full line of Fig. 2b. The lower (dashed) line shows the result of a direct numerical evaluation of using the stochastic differential Eqs. (31) with and averaging over the realizations of the noise. The agreement with the result obtained using the fluctuation–dissipation formalism is quite good, considering the inevitable accumulation of small errors for long integration times.

We next proceed to analytic estimates of the plateau value in Fig. 2b. We start with the noise-free limit and case 1, corresponding to potential [Eq. (32a)]. Setting as before and decomposing into the sum of a reference value corresponding to the unforced system () and of a perturbation induced by , we obtain straightforwardly from the linearized form of Eqs. (31)

with

To estimate the effect of noise, we average Eqs. (31) over the different realizations of the random force. This yields

Decomposing *T* into the sum of a mean value and a fluctuation , we rewrite this relation as

The term in accounts for the effect of fluctuations. It can be written as the product of a negative factor and the variance of the random force . On the other hand, may in turn be decomposed into the sum of the reference value and a correction *m* owing to both the forcing and the fluctuations:

One obtains then, to the dominant order,

Comparing with Eq. (35a), we conclude that the plateau value reached by the response tends to be reduced by the fluctuations. This is similar to the conclusion reached in the one-variable case studied in the previous subsection. Using the parameter values adopted in Fig. 2b, one arrives at , in very good agreement with the plateau value found from the linear response formalism and from the numerical integration of the stochastic differential equations.

## 4. General systems

In the majority of the cases encountered in atmospheric and climate physics, the evolution of the variables of interest does not derive from a potential and conditions (18) or (21) do not hold. In this section, we develop a linear response theory for such systems and analyze the structure of the associated response functions in a fluctuation–dissipation theorem perspective.

Our starting point is Eqs. (14) and (15) where, at variance with section 3, the generalized potential appearing in the exponent of the invariant probability of the reference system does not satisfy the integrability condition (17). To obtain a more explicit form of Eq. (15) we need to evaluate this quantity and, to this end, we substitute expression (13) into the Fokker–Planck equation [Eq. (12a)] for . This yields

For any given set , assumed to be of a polynomial form with respect to , this equation allows to determine to any desired order upon expanding in powers and evaluating the expansion coefficients by equating coefficients of equal powers in . These coefficients are expected to depend in a complicated way on the covariance matrix of the noise and on the various indicators (Jacobian matrix of , etc.) of the dynamics of the reference system. As a result, as pointed out already in section 2, expression (15) will not reduce to a simple linear correlation function of variables, even if the observable is chosen to be itself.

### a. A case study: Saltzman’s sea ice–ocean surface temperature model

As an illustration of the general linear response formalism for nonpotential systems, we revisit in this subsection the classic global climate model proposed by Saltzman (Saltzman et al. 1982) to describe the coupled sea ice extent–ocean surface temperature dynamics. Let be the deviation of the sine of the latitude of sea ice extent from a reference state and be the excess mean ocean temperature. In the absence of forcing and noise, the Saltzman equations expressed in dimensionless form for both the variables and the parameters read

These equations have been analyzed extensively in the literature. We briefly summarize the main results (Nicolis 1984b).

#### 1) Steady states

These solutions arise from the trivial one through a pitchfork bifurcation occurring at . They exist as long as ; otherwise, the trivial state is the unique steady-state solution.

#### 2) Stability

Setting

and linearizing Eqs. (38) with respect to yields

or, more explicitly,

The Jacobian matrix of around is given by

Computation of the eigenvalues of shows that as long as and , the system undergoes damped oscillations around the trivial state for *b* close to but less than 1. A Hopf bifurcation takes place at and, in the range , the system undergoes sustained (limit cycle) oscillations around the trivial state. In the range , the trivial state becomes unstable and the system performs damped or sustained oscillations around one of the nontrivial steady states according to whether *a* is less or greater than 1, being a value at which a Hopf bifurcation is taking place.

We next consider the stochastic dynamics in the presence of an external forcing. Equations (38) are replaced by

### b. Gaussian limit

Our next step is to study Eq. (37) and its specific form [Eq. (43)] in the first nontrivial order, where are limited to their linear parts with respect to the deviations from a reference state [Eqs. (40)]. Expanding in these excess variables will then yield a power series of the form with . In this subsection we will evaluate explicitly in the limit where only the quadratic terms are retained or, alternatively, where the invariant density reduces to a Gaussian.

Setting

Introducing the inverse matrix of one obtains the more compact form

where ^{+} is the adjoint of .

Starting from Eq. (46), we now derive explicit values of for Saltzman’s model, focusing on the case of damped oscillations around the trivial state ( and *b* close to but less than 1). Using Eqs. (41a) and (43) and taking in addition , , and , we arrive at

where

Using Eq. (44), we may thus write the response function in the Gaussian limit in the form

where satisfy (46) and the angle brackets denote as before average over expressed in terms of the excess variable . Again, as pointed out in the beginning of this section, this expression does not reduce to a simple linear correlation of variables even if is chosen to be itself. In this latter case one obtains, rather, a linear combination of correlation functions of variables whose coefficients depend on the parameters in a complicated and not a priori prescribed way. Nevertheless, the following points are worth noticing. First, to the extent that can be regarded as a generalized potential and by analogy with thermodynamics (De Groot and Mazur 1962), the sum over *k* in Eq. (48)—which is minus the derivative of with respect to —can be regarded as a generalized force acting as long as is not zero. Second, the linear part of the evolution of the unforced noise-free reference system [Eq. (40b)] can be regarded as a flux that is again active as soon as is not zero and can in this respect be regarded as driven by the aforementioned generalized force. If is now chosen to be one of these two quantities, then the response function can be cast in the form of a generalized force–force autocorrelation function or of a generalized flux–force cross-correlation function.

As an example, for Saltzman’s model choosing , constant, and , we obtain

with

where and are given by Eqs. (47), or

displaying the autocorrelation as and the cross correlation of and as . These quantities can be evaluated explicitly by Laplace transform methods as seen in detail in the appendix. This results in the following expression for the response in the long time limit:

Here the quantities in angle brackets are the second moments of the invariant distribution in the Gaussian limit, which under the additional assumption are given by

We notice that all moments diverge as *b* approaches the Hopf bifurcation point from below. This holds true for the correlation functions and as well. On the other hand, the asymptotic response remains finite, owing to the cancellation of the factor in Eq. (50); see Eqs. (47) and (51). This can be checked independently by averaging Eq. (42) over the realizations of the noise. Linearizing in and (consistent with the Gaussian limit considered here), taking , and solving for as a function of yields

which remains, indeed, finite when *b* crosses the value .

### c. Beyond the Gaussian: Nonlinear effects

Figure 3a depicts the normalized correlation functions and of Saltzman’s model as functions of time for , , and . The figure has been obtained by solving numerically for and the full nonlinear stochastic differential Eqs. (42) in the absence of the forcing and averaging subsequently the products and over (or equivalently over the realizations of the noise, as the system is ergodic). Both functions tend to zero by performing damped oscillations, as expected. The average trajectory of as computed numerically from the full Eqs. (42), now including the forcing, is shown in Fig. 3b for and . As can be seen, the response fluctuates around a plateau value of about −0.0083.

These behaviors are in qualitative agreement with the theoretical predictions derived in the previous subsection in the Gaussian limit. In particular, the shape of the correlation functions in Fig. 3a follows the analytic expressions (Eqs. (A5)] derived in the appendix and the plateau value in Fig. 3b is close to the theoretical estimate of −0.01 for the same parameter values provided by Eqs. (50)–(52). Still, the comparison suggests that nonlinear and noise effects tend to induce deviations from the Gaussian behavior and, in particular, to reduce the magnitude of the response. The origin of this trend can be seen analytically. Upon averaging Eqs. (42) over the realizations of the noise and expressing as a function of , one obtains (setting again )

We notice that the first term in the right-hand side is identical to the result of Eq. (52) obtained in the Gaussian limit. The second term, which represents the nonlinearity and fluctuation-induced correction manifested beyond the Gaussian limit, can be estimated by decomposing and into the sum of their averages plus the fluctuations. One obtains, keeping terms up to second-order variances,

Substituting as a zeroth-order estimate , from Eqs. (52), we may further write this expression as

or, using the Gaussian limit expressions (51) as a zeroth-order estimate,

For *a* and *b* values compatible with the occurrence of damped oscillations around the trivial state (, , ), both terms in Eq. (54) are negative, entailing that the nonlinearity and noise-induced correction in Eq. (53) is positive and hence that the level of the response tends to be reduced. This is in accord with the results in Fig. 3b. We notice that for given value of the strength of noise the magnitude of the nonlinear effects increases as *b* approaches the bifurcation point and/or the value of *a*. In this range the behavior is markedly non-Gaussian and a full-scale nonlinear analysis using Eq. (43) becomes necessary.

## 5. Conclusions

In this work we applied linear response theory to nonlinear dynamical systems of interest in atmospheric and climate dynamics, under the assumption that the source of their intrinsic variability can be modeled as a Gaussian Markov noise process. Using a Fokker–Planck equation approach, we derived formal expressions linking the perturbation-induced shift of an observable to the dynamical properties of the reference, perturbation-free system. These expressions have the general structure of time cross-correlation functions featuring a generalized potential–like function, which determines the structure of the invariant probability density and satisfies a nonlinear partial differential equation whose structure depends on the underlying dynamics. As a corollary, it turned out that, typically, the response cannot be cast in the form of autocorrelation and cross-correlation functions of the state variables as in classic versions of the fluctuation–dissipation theorem.

Our approach was carried out along the same lines as in recent reformulations and extensions of linear response theory (Majda et al. 2005; Baiesi and Maes 2013). A linear response approach to study the effect of perturbations in atmospheric physics–related problems was also developed by Lucarini (2012) and Lucarini and Sarno (2011), with emphasis on the structure of the power spectrum. While sharing a similar perspective, the present work has focused on different aspects such as the connection between the statistical properties of the response and the nature of the underlying dynamics beyond the Gaussian approximation usually adopted in the literature. Full solutions and explicit expressions were obtained for the class of potential systems, for which the roles of the parameters and of the nonlinearities in the deviations from Gaussian behavior and the sensitivity of the response have been determined in detail. Systematic approaches for more general systems were subsequently outlined, based on power-series solutions of the equation for the generalized potential. Explicit relations between the quantities of interest in linear response and the underlying dynamics have been obtained in the Gaussian limit and illustrated on a representative example. Corrections associated with deviations from Gaussian behavior have also been identified.

Broadly speaking, there are two qualitatively different situations in which deviations from Gaussian behavior can be expected. In the first and simplest one, the system operates around a single stable steady state. Nonlinearities induce then a slight (in the weak noise limit) shift of its probability density from the Gaussian in, for instance, the form of skewness around the most probable value entailing in turn a slight difference between the mean and the most probable value. These corrections are amplified if the system is brought close to criticality, as seen in section 3a when parameter *μ* tends to zero starting from negative values.

The second, more intricate situation pertains to the postcritical behavior of a system beyond a bifurcation point. Depending on the case, the system will then possess multiple steady states or periodic or aperiodic (e.g., chaotic) attractors. Accordingly, the corresponding probability densities will be multihamped or complex crater-like surfaces peaked in state space around more than one fixed point or along a fractal manifold. They will thus be manifestly non-Gaussian and their mean values will generally correspond to unrepresentative (e.g., unstable) states. To determine their analytic form, solutions of equations like Eq. (37) or (43) will have to be carried out. This can be achieved systematically in the vicinity of transcritical, pitchfork, and Hopf bifurcations, where the center manifold theorem (Guckenheimer and Holmes 1983) allows one to limit the expression of the generalized potential to quartic terms (Nicolis and Nicolis 1987). Linear response theory will then be expressed in terms of correlation functions involving nonlinear combinations of observables up to order 3.

In most of the analyses reported in the context of linear response, the random forces accounting for internal variability are assumed to satisfy certain prescribed relations. In reality, as alluded in section 3, it can be legitimately expected that, to the extent that these forces have an intrinsic origin, they should satisfy certain constraints imposed by the dynamics in the noise-free limit. We have provided some examples of such links, to which we referred as fluctuation–dissipation theorem of a second kind (Nicolis 1984a; see also Nicolis et al. 1985; Nicolis 2002), but it would certainly be worth undertaking more extensive work along these lines. Last but not least, fundamentally, the laws of nature are deterministic. Noise source terms should therefore be related in one way or the other to dynamical processes present on different scales in the atmosphere and climate. In some cases, notably in the presence of wide time-scale separation, it has been shown (Arnold et al. 2003) that the full-scale deterministic (and, typically, chaotic) dynamics can indeed be mapped into the stochastic dynamics of a subset of the original variables. Still, the general problem remains open. It is expected that low-dimensional chaotic dynamical systems should behave radically differently from high-dimensional ones. Indeed, in the first class the invariant probability density displays a markedly singular structure (Ruelle 1998) while for the second class a central limit theorem and a regularization of singularities are expected to hold (Kaneko 1989). It would be very interesting to use multidimensional chaos as the link to establish a connection between the stochastic approach followed in this work and a workable, purely deterministic formulation.

## Acknowledgments

This work is supported, in part, by the European Space Agency and the Belgian Federal Science Policy Office.

### APPENDIX

#### Analytic Evaluation of the Correlation Functions and [Eq. (49b)] for Saltzman’s Model

We start with the linearized Saltzman equations of the reference system () around the trivial state in the presence of noise,

Taking a Laplace transform on both sides, solving for the transforms and of and , multiplying by and , and averaging over the realizations of the noise, we obtain

where the angle brackets on the right-hand sides denote averages over the invariant distribution in the Gaussian limit. Taking the inverse Laplace transform to deduce and involves integration in the complex plane on a closed contour surrounding the singularities of and —that is, the solutions of the equation

The solutions of this equation are just the eigenvalues of the linear stability matrix in Eq. (41a) and are given by

For *b* close to but less than one, , and , they are complex conjugate numbers with negative real parts.

Performing the integration in the complex plane using the calculus of residues, we obtain straightforwardly

Integrating over *t* from zero to infinity and utilizing Eq. (A3) leads to the following expressions for the integral correlation functions

and ultimately to the expression in Eq. (50) for the linear response .

## REFERENCES

*Nonequilibrium Thermodynamics.*North Holland, 510 pp.

*Handbook of Stochastic Methods.*Springer, 442 pp.

*J. Atmos. Sci.,*

*Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields.*Springer, 459 pp.

*Information Theory and Stochastics for Multiscale Systems.*CRM Monograph Series, Vol. 25, American Mathematical Society, 141 pp.

*J. Atmos. Sci.,*

*Irreversible Phenomena and Dynamical Systems Analysis in the Geosciences*. Reidel, 578 pp.

*Introduction to Nonlinear Science*. Cambridge University Press, 254 pp.

*Foundations of Complex Systems*. World Scientific, 367 pp.

*Stabilité Structurelle est Mophpogenése.*Benjamin, 362 pp.