A theory of average response to large jump perturbations
Rafail V. Abramov
TL;DR
This paper extends the Fluctuation Dissipation theorem to large, instantaneous jumps in nonlinear dynamical systems, showing that average responses can still be expressed via time correlation functions of unperturbed dynamics, including stochastic cases.
Contribution
It develops a theoretical framework for understanding the average response of nonlinear systems to large, instantaneous perturbations, both deterministic and stochastic, using correlation functions.
Findings
Average response formulas derived for large jumps
Applicable to stochastic and deterministic jumps
Efficient approximations for practical use
Abstract
A key feature of the classical Fluctuation Dissipation theorem is its ability to approximate the average response of a dynamical system to a sufficiently small external perturbation from an appropriate time correlation function of the unperturbed dynamics of this system. In the present work, we examine the situation where the state of a nonlinear dynamical system is perturbed by a finitely large, instantaneous external perturbation (jump) -- for example, the Earth climate perturbed by an extinction level event. Such jump can be either deterministic or stochastic, and in the case of a stochastic jump its randomness can be spatial, or temporal, or both. We show that, even for large instantaneous jumps, the average response of the system can be expressed in the form of a suitable time correlation function of the corresponding unperturbed dynamics. For stochastic jumps, we consider two…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
A theory of average response to large jump perturbations
Rafail V. Abramov
Department of Mathematics, Statistics and Computer Science, University of Illinois at Chicago, 851 S. Morgan st., Chicago, IL 60607
Abstract.
A key feature of the classical Fluctuation Dissipation theorem is its ability to approximate the average response of a dynamical system to a sufficiently small external perturbation from an appropriate time correlation function of the unperturbed dynamics of this system. In the present work, we examine the situation where the state of a nonlinear dynamical system is perturbed by a finitely large, instantaneous external perturbation (jump) – for example, the Earth climate perturbed by an extinction level event. Such jump can be either deterministic or stochastic, and in the case of a stochastic jump its randomness can be spatial, or temporal, or both. We show that, even for large instantaneous jumps, the average response of the system can be expressed in the form of a suitable time correlation function of the corresponding unperturbed dynamics. For stochastic jumps, we consider two situations: one where a single spatially random jump of a system state occurs at a predetermined time, and another where jumps occur randomly in time with small space-time dependent statistical intensity. For all studied configurations, we compute the corresponding average response formulas in the form of suitable time correlation functions of the unperturbed dynamics. Some efficiently computable approximations are derived for practical modeling scenarios.
1. Introduction
The classical Fluctuation Dissipation theorem (FDT) [54, 55, 76, 67] provides a leading order approximation to the statistical response of a dynamical system to a small deterministic external perturbation via statistical correlations of the unperturbed dynamics. The FDT offered more insight into statistical properties of dynamical processes near equilibrium in various scientific applications, such as statistical mechanics of identical particles [25, 34, 56], Ornstein–Uhlenbeck Brownian motion [54, 55, 68, 79], motion of electric charges [67, 80], turbulence [52, 53], quantum field theory [26, 36, 38, 81], chemical physics [77, 78], and physical chemistry [65]. In geophysical science, the FDT was proposed as a sensible approximation for appropriate variables in a complex climate system [46, 57, 58] despite the absence of a classical Gaussian equilibrium state of the traditional statistical mechanics. This observation spurred a series of works [22, 30, 29, 57, 27, 28, 41, 40, 42, 43, 44, 49, 61, 62, 66, 70], where various applications of the FDT in the weather and climate modeling were proposed. In the author’s past works [14, 15, 16, 1, 2, 3, 6, 8], a computational framework predicting the average response of both deterministic and stochastic dynamical systems to a small deterministic or stochastic external perturbation was developed, studied, and used in a new method for the parameterization of unresolved processes in reduced models of multiscale dynamics [4, 5, 7, 10].
Thus far, the extensively studied types of external perturbations were largely limited to small bounded forcing perturbations – either deterministic, or in the form of a Brownian motion. However, in many practical applications the external perturbations are, first, not necessarily small, and, second, not necessarily in the form of a bounded external forcing. For example, a physical system may experience an external “impulse forcing” – that is, the external forcing in the form of a delta-function, which instantaneously changes the state of the system. Moreover, this change can be finitely large – that is, the delta-function of the impulse forcing may not be necessarily scaled by a small parameter.
An ubiquitous example of such a global event is known as the Permian–Triassic extinction [33, 50, 82]. Due to an unknown cause, catastrophic global climate changes occurred on Earth some 250 million years ago between Permian and Triassic periods. As a result, about 90% of all living species has become extinct [45, 72, 74, 75]. One of the accepted scientific hypotheses is that this event was triggered by a large meteorite impact [20, 21]. It is strikingly obvious that any mathematically and physically adequate explanation or model of this event cannot possibly assume that this impact was “small” in any reasonable sense – although one can likely reasonably assume that this event was “instantaneous”, at least relative to the time scale of the subsequent global climate change. In addition, another reasonable assumption is that such events tend to occur randomly in time, although they are statistically unlikely to occur frequently.
In the current work, we focus on the scenario where the state of a nonlinear dynamical system is perturbed by a finitely large, instantaneous external perturbation, or “jump”. Such a scenario was examined previously [24] for the simple case of a fixed (that is, independent of the system state) perturbation at a prescribed time along a single component of the system state vector. For generality, here we assume that this jump may depend on the state of the system which immediately precedes the time of the jump. Furthermore, we assume that the jump can be either deterministic or stochastic, and for a stochastic jump we further assume that its randomness can manifest in the spatial configuration of the jump, or its temporal frequency, or both. We find that the average response of the system can be expressed in the form of a suitable time correlation function of the corresponding unperturbed dynamics even for large instantaneous jumps. In the case of jumps which incorporate a random component, we examine two scenarios. The first one is where a single spatially random jump of a system state occurs at a predetermined time, while the second scenario is the one where jumps occur randomly in time with small space-time dependent statistical intensity. We also find that, for all studied configurations, the corresponding average response formulas are computable in the form of suitable time correlation functions of the unperturbed dynamics, just like in the classical FDT formulation. For practical modeling scenarios, we also derive suitable approximations of the general response formulas which can be efficiently computed in the context of a numerical simulation.
The work is organized as follows. In Section 2 we state the general form of the underlying, unperturbed nonlinear dynamical system, together with its corresponding forward Kolmogorov (or Fokker–Planck) equation [19, 69, 76]. In Section 3 we derive the explicit formula of the average response of the system to a deterministic instantaneous jump perturbation. In Section 4 we extend the previously derived average response formula onto the case of a random instantaneous perturbation, which nonetheless occurs at a prescribed time. In Section 5 we derive the average response formula for the most general perturbation scenario, where spatially random jump perturbations also occur at random times. In Section 7 we discuss the results and suggest future research directions.
2. Unperturbed dynamics
In this section, we need to specify the general form of an unperturbed dynamical system for a typical application in natural sciences. From what is to be presented below, the basic requirements on the form of the unperturbed dynamical system can be summarized as follows. First, the unperturbed dynamical system may incorporate a deterministic vector field, as well as stochastic effects. We will, however, restrict the form of the unperturbed dynamical system so that its solutions are continuous, as it is implied that the discontinuities will be introduced via the jump perturbations. Additionally, we require that the unperturbed system has an explicitly formulated forward Kolmogorov equation [69, 76] with a differentiable stationary solution. The corresponding perturbed system must also have a forward Kolmogorov equation, although we are not going to consider its stationary solutions, if any.
Other than the points listed above, there do not appear to be any other fundamental constraints which would impose further unmitigable issues. In particular, the presence of the stochastic component of the dynamics is not necessarily required for the stationary solution of the forward Kolmogorov equation to be differentiable – for example, if the deterministic vector field preserves a quadratic energy, the uniform distribution on a constant energy surface is usually both the stationary statistical state for the system and a smooth distribution [12, 13, 11].
Thus, throughout the work, we will assume that the unperturbed dynamical system is described, in general, via the following Itô stochastic differential equation:
[TABLE]
Here, is the state vector of the system, is a -dimensional Wiener process, while and are smooth vector fields. As motivated above, the form of (2.1) is chosen so that it is the most general form of the Lévy-type Feller process [35] whose solutions are almost surely continuous, and which, under a random jump perturbation, retains the form of the infinitesimal generator compatible with Courrège’s theorem [31], so that the forward Kolmogorov equation of the jump-perturbed process can be obtained in an explicit manner [19].
Let be a twice differentiable function with bounded second derivatives. Let denote the conditional expectation of at time , provided that the initial state of the system at time is . Then, the infinitesimal generator [19, 39, 69] of is given via:
[TABLE]
where “” denotes the Frobenius matrix product.
Let denote the probability distribution of solutions of (2.1). Then, we can relate the -average of at time to the probability density via
[TABLE]
The forward partial differential equation for is known either as the forward Kolmogorov equation [19, 39, 69] or as the Fokker-Planck equation [76], and is given via
[TABLE]
In what follows, we will assume that (2.4) has a stationary solution ,
[TABLE]
If the union of all solutions to (2.5) consists of multiple ergodic components, we will assume below that , whose support we denote as , is the indecomposable ergodic component which is “physically relevant” – that is, a generic initial condition to (2.1) almost certainly falls into , and the subsequent jump perturbations are such that the state of the system in (2.1) never leaves . This assumption is necessitated by the need to use the Birkhoff–Khinchin theorem [23, 51] for the practical computation of the average response further below.
As an example of a suitable invariant state, can in its entirety be ergodic and supported on the whole phase space, which typically happens when and have bounded derivatives of all orders in , and, in addition to that, the matrix product is uniformly positive definite in [73]. An elementary example of such a system is the Ornstein–Uhlenbeck process [79]. Besides, the assumption of ergodicity of is a standard assumption in many works on this topic [14, 15, 16, 1, 2, 3, 6, 8]. In such a case, the Birkhoff–Khinchin theorem applies for any finite jump perturbation of such a system.
Another example is a system whose (possibly deterministic) vector field preserves the quadratic kinetic energy, in which case belongs to a family of the corresponding microcanonical Gibbs states for each constant energy surface [11, 12, 13]. One can then introduce the jump perturbation of the state of the system in the form of a “molecular collision” [9, 47], given via an explicit formula which preserves the quadratic energy, which leaves the state of the system on its assigned constant energy surface after the jump. More precisely, let us assume that (2.1) is deterministic, possesses the Liouville property
[TABLE]
and preserves the quadratic energy
[TABLE]
along any trajectory . An example of such a system is the truncated Burgers–Hopf model [11, 62]. It can be shown [11, 62] that the microcanonical Gibbs state, which is a uniform probability distribution on the surface of constant energy, is automatically the invariant statistical state of (2.1).
Let for some positive integer , with , and let and be computed via the transformation
[TABLE]
For , the transformation above in (2.8) describes the change of velocities during a collision of two hard spheres [9, 47]. It can be shown trivially that the sum of squared norms is preserved by this transformation. Indeed, observe that
[TABLE]
Clearly, if and are themselves two distinct subsets of components of , then their replacement with and does not change the norm , and, therefore, leaves the perturbed vector on its constant energy surface. Thus, the transformation in (2.8) is the example of an energy-preserving jump perturbation of a subset of components of . Note that, aside from , is otherwise arbitrary, and can even be a random variable.
2.1. Conditional probability density
In what follows, it is convenient to introduce the conditional probability density , defined via
[TABLE]
From (2.3), and observing that is arbitrary, it follows that
[TABLE]
Substituting the above equation into (2.4) and stripping the integral over , we obtain the equation for :
[TABLE]
Observing that the above equation is autonomous with respect to , it is clear that
[TABLE]
that is, the conditional probability density is a function of difference between the starting and ending times. Therefore, without loss of generality, we can assume that , which yields
[TABLE]
Consequently, according to (2.10), the conditional expectation is also a function of the time difference: .
2.2. Special case: the Ornstein–Uhlenbeck process
Here we consider a special case of the unperturbed dynamical system in (2.1), which has the convenience of being exactly solvable. Consider the following centered Ornstein–Uhlenbeck process [79]:
[TABLE]
Above, is a constant positive definite matrix, while is a constant matrix, with being positive definite. The solution to the Ornstein–Uhlenbeck process above is given via Duhamel’s principle:
[TABLE]
Via the properties of the Itô integral [48, 39, 69], it can be readily seen that the expectation of the state variable satisfies
[TABLE]
The stationary probability density of the Ornstein–Uhlenbeck process in (2.15) is given via
[TABLE]
where is the covariance matrix of the process, given via the equation
[TABLE]
More details on this topic can be found in [76].
3. A large deterministic jump perturbation
In this work, the definition of the average response is identical to that in our earlier works on this topic [15, 14, 16, 62, 1, 2, 3, 6, 8]. Namely, assume that we have a large statistical ensemble of solutions of (2.1). This ensemble is initially distributed according to . At the prescribed time , the following instantaneous jump perturbation is applied to the states of all ensemble members:
[TABLE]
Above in (3.1), is a continuous function, and denotes the left-limit at . Observe that the action of can be interpreted as the instantaneous jump perturbation of a given trajectory, which, generally, depends on the pre-jump state (but not on the time of the jump). As noted above, here we assume that for any , – that is, the state of the system does not leave the support of the physically relevant ergodic component of (2.1) as a result of the jump.
The resulting statistical discrepancy between the perturbed ensemble and the unperturbed ensemble, as a function of time past the moment of the perturbation, is what we refer to as the “average response”. The key difference between the previous studies and the current work is that in our previous works a time-dependent forcing was applied at the initial time and past that, whereas here the perturbed trajectories continue according to (2.1) once the jump perturbation has been applied.
The average response of the statistical ensemble is measured via the difference between the perturbed and unperturbed conditional expectations of a test function , averaged over the equilibrium probability density :
[TABLE]
where is the time when the external perturbation has occurred.
3.1. An exactly solvable example
As an example, we can easily compute the mean state response of the Ornstein–Uhlenbeck process in (2.15), with help of (2.17):
[TABLE]
where we observe that is indecomposable and supported on the whole space, and thus is the same as . Above, we only need to be able to evaluate the integral in terms of elementary functions, which, given the form of in (2.18), leaves a broad choice for (it can be, for example, a polynomial in ).
3.2. A practical average response formula for computational
modeling
It is also a possibility that the dynamical system in (2.1) is not exactly solvable, although can be numerically simulated or modeled. In this case, we need a suitable adaptation of the response formula in (3.2) to the typical practical restrictions of numerical modeling. Here, we will assume that, first, there exists an approximation to the stationary probability distribution in terms of elementary functions, and, second, there does not exist a similar approximation to . These assumptions are reasonable for a range of prototype nonlinear dynamical systems, such as the Lorenz 96 model [62, 59, 60], the barotropic model of Earth atmosphere [16, 18, 37], or the quasigeostrophic 1.5-layer wind driven double gyre ocean circulation model [17, 63, 64].
In such a situation, let be the inverse of . Clearly, satisfies
[TABLE]
Then, in the perturbed expectation integral, we can change the variables as follows:
[TABLE]
where is the Jacobian of . The average response formula in (3.2) can thus be written in the form
[TABLE]
where we can divide by since it is nonzero in . Finally, for the purposes of practical computation, we use the ergodicity property of in and replace, with help of the Birkhoff–Khinchin theorem [23, 51], the measure average with the following time correlation function over the long-term trajectory of the unperturbed system in (2.1):
[TABLE]
where the initial condition for the time series is presumed to be taken in . The practical computational constraint here (aside from the necessity for an accurate approximation for ) is the ability to invert the function . While, obviously, a large variety of possible forms of the jump function is available in general, here we point out a simple form of which is explicitly invertible and is likely broad enough for the majority of practical applications:
[TABLE]
where is a constant -vector, and is a constant matrix. For such a form of , we obtain, explicitly,
[TABLE]
A special case of this scenario was studied previously [24] with a fixed perturbation of the form
[TABLE]
that is, the fixed jump perturbation was applied to a single component of the system. It is easy to see that, in such a case, , with its Jacobian being equal to 1.
Finally, we would like to emphasize that none of the average response formulas above rely on a “small parameter” of any kind – any imprecisions of the average response computation above will manifest due to, for example, inaccuracy of the approximation for , or the statistical undersampling of the long-term trajectory of (2.1).
3.3. The quasi-Gaussian approximation
The quasi-Gaussian approximation [15, 14, 16, 62] for the average response formula in (3.7) emerges when the stationary distribution is replaced by the corresponding Gaussian distribution with the same mean state and covariance matrix. More precisely, let be the mean state of , and let be its covariance matrix. Then, the quasi-Gaussian approximation for (3.7) is given via
[TABLE]
where is given explicitly via
[TABLE]
The above expression, together with (3.9), renders the integrand of the time correlation function in (3.11) explicitly computable for a given (computed or observed) time series , and thus the average response in (3.11) can be computed numerically using standard methods [15, 14, 16, 62].
It may seem that, with the replacement of with in (3.11), one may lose the information of the relevant ergodic component of (that is, if the latter is not ergodic on the whole space ). However, observe that it is not the case – due to the fact that the time average in (3.11) is computed over the exact time series of the original dynamical system in (2.1), the statistical sampling of the integrand of (3.11) still occurs on with the distribution . Moreover, the most remarkable property of the time correlation function in (3.11) is that an observed or historically recorded time series can be used for its evaluation, without any need to numerically simulate the dynamical process. This suggests that the average response could be computed with sufficient degree of accuracy even for complex geophysical phenomena, as long as the recorded time series are sufficiently detailed.
4. A large random jump perturbation
An interesting generalization of the previous scenario is the randomization of the jump (which still occurs at the time ). Namely, let be a random variable with the distribution measure , and let be a continuous function of , and bounded of . Then, it is natural to define the average response of not only by averaging over the states of the system (that is, over ), but also over the possible jumps (that is, over ):
[TABLE]
Observe that the only difference between (3.2) and (4.1) is the additional average over , which allows to easily extend the average response formulas in (3.3) and (3.7) onto the random perturbation in a straightforward fashion.
4.1. An exactly solvable example
The corresponding exactly solvable average response of the mean state of the Ornstein–Uhlenbeck process is given via
[TABLE]
which can be easily seen by taking and recalling the formula for the expectation in (2.17). The practical computability of the mean state response formula in (4.2) now depends, in addition to the choice of , also on the choice of the intensity measure . In the case of being of the form
[TABLE]
we find that the double integral becomes the sum of the product of single integrals:
[TABLE]
Here, observe that, for each term in the sum, we obtained the formula for the mean state response of an Ornstein–Uhlenbeck process to a deterministic perturbation in (3.3), additionally multiplied by the integral over . The computability of this integral depends entirely on the choice of and , and is independent of the rest of the set-up.
4.2. A practical average response formula for computational
modeling
In order to obtain an analog of the time-averaged response formula in (3.7) for the response to the random perturbation in (4.1), we follow the same steps as in the previous section. We let be the -inverse of for a given :
[TABLE]
Then, we can again replace in the perturbed integral and write
[TABLE]
where in the unperturbed integral we used the fact that is normalized to 1, and, as before, is the Jacobian of in the -variable. The average response formula in (4.1) can thus be written in the form
[TABLE]
where
[TABLE]
Upon the application of the Birkhoff–Khinchin theorem [23, 51], the average response formula becomes
[TABLE]
Here, again observe that the average response formula above does not rely on a small parameter of any kind. Any imprecisions of the average response computation above will manifest due to, for example, inaccuracy of the approximation for , or a possible approximation for , or the statistical undersampling of the long-term trajectory of (2.1).
However, note that for an efficient numerical computation of the time correlation function in (4.9), the term in (4.8) must be expressed in terms of elementary functions. This requirement, obviously, places restrictions on the choice of , , and the approximation for . Below, we consider a few special cases where the -integral above is explicitly computable.
4.3. Special case 1: a set of spatially pre-determined jumps
As a simplest example where the -integral is explicitly computable, we consider the following special case. Assume that, instead of a single deterministic jump, several different, yet spatially pre-determined jumps may occur randomly with prescribed probabilities. In other words, the random variable may return only a finite set of values , each with probability , . In this case, we can define
[TABLE]
and, therefore,
[TABLE]
Observe that, as long as and are available in the form of explicit formulas, is also an explicit function of , which allows to compute the time correlation function in (4.9) efficiently. For a practically computable example of , we can generalize (3.8) to include the -dependence:
[TABLE]
where the -vector and matrix are known, explicit functions of . This leads to
[TABLE]
and, subsequently,
[TABLE]
4.4. Special case 2: a Gaussian jump distribution
As a practical example, let us consider the scenario where has a Gaussian density, with its own mean state vector and covariance matrix :
[TABLE]
For practical computation, we restrict the -dependent jump function in (4.12) to
[TABLE]
with , and being the constant -vector, matrix and matrix, respectively. This yields in the form
[TABLE]
The form of in (4.17) leads to the following expression for in (4.8):
[TABLE]
To be able to compute the -integral above in the form of an explicit formula, we will assume that is also given by its quasi-Gaussian approximation in (3.12). This leads to
[TABLE]
where the terms , and are given via
[TABLE]
The details of the computation are given in the Appendix A.
Observe that the whole expression in (4.19) is an explicit function of ( is constant, is linear in , and is quadratic in ), which means that it can be evaluated efficiently for the computation of the time correlation function in (4.9).
4.5. Special case 3: a linear combination of Gaussian densities
The natural generalization of the above formula can be easily derived for a scenario where both and are linear combinations of Gaussian densities:
[TABLE]
for some integers , . It is easy to verify that the corresponding integral in (4.8) becomes
[TABLE]
with
[TABLE]
where and are the mean state and covariance matrix of , respectively. Again, observe that are constants, while and are explicit functions of . This means that the measure integral can be evaluated efficiently in the numerical computation of the time correlation function in (4.9).
5. A sequence of large random jump perturbations at random times
An important generalization of the previously studied random jump process is the extension of the randomness of the perturbation onto time. Recall that the traditional fluctuation dissipation setting [44, 42, 1, 2, 3, 6, 10, 15, 14, 16, 17, 62] typically consists of a statistical ensemble of solutions of the unperturbed system in (2.1), which are initially distributed according to the invariant probability measure , given via (2.5). Then, at the initial time , a small deterministic perturbation is added to each member of this statistical ensemble. This perturbation may be an explicit function of both time and the state of that ensemble member. Then, the average response of a function is the difference of the ensemble average values between the perturbed and unperturbed statistical ensembles, starting at time :
[TABLE]
Above, is the solution of the forward Kolmogorov equation of the corresponding perturbed system, with the initial condition at given via .
In the current work, we retain the definition of the average response as in (5.1), but replace the small deterministic perturbation of the traditional setting with a new, different perturbation process, which consists of random jumps of finitely large magnitude. As follows below, we formulate the random perturbation process in such a way that its jumps satisfy the following general properties:
- (1)
The spacing between the times of jumps is random. However, the statistical frequency of the jumps (the so-called “intensity”) is explicitly prescribed, and may conditionally depend, in a specified way, on both the time and the pre-jump state of the system. We will assume that the intensity of the jump times is small. 2. (2)
The jumps themselves are also random, whenever they occur. However, their spatial statistical distribution can be explicitly prescribed, and may conditionally depend, in a specified way, on the pre-jump state of the system. It may not, however, depend on the time at the moment of the jump. There will be no restriction on how large the jumps can be (except that they must be finite).
As we can see, the properties of the jumps above are general enough to allow for a wide variety of prototype scenarios in practical applications. Next, we are going to formulate the requisite jump perturbation process so that it indeed satisfies the properties delineated above.
5.1. The formulation of the random perturbation process
Recall that the unperturbed system in (2.1) may include a Wiener process , which in itself is a random variable. Thus, we introduce the corresponding probability space equipped with a filtration , such that for , . The Wiener process in (2.1), as well as the solution , are chosen to be adapted to .
At this point, we proceed with the definition of the random process, which will perturb the trajectory of the dynamical system in (2.1) via instantaneous random jumps. Our goal here is to formulate the process in such a way so as to have explicit control over the statistical distribution of jump times, as well as the spatial distribution of the jumps themselves when they occur. We start with defining the properties of the statistical distribution of jump times. Here, we choose the statistical intensity of the times, at which jumps occur, to be of the form
[TABLE]
Above, is a specified bounded strictly positive function on a real line, is a specified bounded strictly positive function on , and is the left-limit at of the perturbed solution. The small constant scaling parameter signifies that the temporal intensity of jumps must be small.
As we can see, the intensity of jump times in (5.2) is conditional – that is, is by itself a random variable, adapted to . This definition of the intensity in (5.2) is specifically chosen to be very broad, since it allows to specify explicitly both the deterministic dependence on the time via , as well as the conditional dependence on the perturbed solution via . Of course, one can set if only the deterministic dependence of the intensity on time is required.
To perform the actual jump perturbations of the state of the dynamical system in (2.1), we introduce a Poisson point process with the corresponding jump process :
[TABLE]
We denote the corresponding Poisson random measure of by :
[TABLE]
As this Poisson process is a part of the external perturbation to (2.1) and thus its properties can be specified however needed, we choose the intensity measure of , given via
[TABLE]
to be normalized to 1, so that the average intensity of jumps of is one jump per unit of time.
In order to adjust the intensity of jumps of the homogeneous Poisson point process to the conditional intensity in (5.2), we introduce the following time-inhomogeneous point process [32, 71] via
[TABLE]
where the compensator is given via
[TABLE]
From the perspective of the theory for point processes with conditional intensities [32, 71], what is presented above should be understood in reverse order – that is, for a point process , adapted to and with the conditional intensity of jump times in (5.2), the corresponding process , defined via (5.6), is a time-homogeneous Poisson point process, whose jump times have statistical intensity 1. However, above we chose a more “constructive” order of presentation for better clarity.
Thus far, we have arranged for the conditional intensity of jump times in (5.2) to be specified as a function of both time and pre-jump state. However, observe that the actual jumps of are identical to those of – and, in particular, their spatial statistical distribution is given via the intensity measure , which is independent of the system state.
To define the conditional dependence of a perturbation on the pre-jump state of the system, we use the same approach as in the previous section – namely, we introduce the jump function , which depends both on the state of the system , as well as the random jump, generated by . We now define the perturbed dynamical system via the following Lévy-type Feller process [19, 35, 31]:
[TABLE]
where is the corresponding Poisson random measure of . Observe that the properties of random jumps are explicitly specified via the following quantities:
- (1)
The statistical intensity of jump times is determined through the choice of the functions and in (5.2), where the former defines the explicit dependence of the intensity on time, while the letter defines the conditional dependence of the intensity on the pre-jump system state. 2. (2)
The spatial distribution of the perturbation jumps is defined via the intensity measure of the Poisson point process (5.6) in conjunction with the jump function . Here, the choice of the measure affects the random spread of the jumps, which can be further adjusted by the choice of , while the latter also allows to specify the conditional dependence of the jumps on the pre-jump state of the system.
Additionally, observe that there is no provision for the perturbation jumps to be small – as mentioned before, it is assumed that the jumps can be finitely large in magnitude. Instead, we require that the statistical frequency of jump times is small, for which a small constant scaling parameter is provided in (5.2).
5.2. The infinitesimal generator
To compute the perturbed probability density in the average response formula (5.1), one first needs to obtain the infinitesimal generator of the perturbed process in (5.8). When the infinitesimal generator is known, its integration by parts against the probability density of the system leads to the forward Kolmogorov equation for [62, 8, 76].
In order to obtain the infinitesimal generator of (5.8), we transform the stochastic jump integral in (5.8) into an integral of a time-homogeneous Poisson process, so that the standard Itô formula [19] for Lévy-type Feller processes [35] could be used. With (5.6), we can write the stochastic -integral in the perturbed process (5.8) as
[TABLE]
Now that the stochastic jump integral in the right-hand side of (5.8) has been expressed via a time-homogeneous Poisson point process, we can obtain the infinitesimal generator of (5.8) using the Itô formula [19]. For a differentiable test function we have the following Itô formula for (5.8):
[TABLE]
Above, denotes the usual Itô integral from to for the drift-diffusion process in (2.1), while denotes the following expression:
[TABLE]
In order to obtain the infinitesimal generator, we need to compute the expectation of the above integral, conditioned on . The expectation of is, of course, given via the usual Itô isometry for drift-diffusion processes. For the stochastic integral over the Poisson random measure , we observe that is time-stationary with independent increments, which leads to the following expression (keep in mind that is -adapted):
[TABLE]
Recalling (5.2), we obtain the following infinitesimal generator for (5.8):
[TABLE]
Above, , and denotes the infinitesimal generator of the unperturbed system in (2.2).
5.3. The forward Kolmogorov equation
Observe that the only difference between the unperturbed (2.2) and perturbed (5.13) infinitesimal generators is that the latter has an additional term corresponding to the random jump perturbation. Therefore, it is easy to extend the unperturbed Kolmogorov equation in (2.4) onto the perturbed dynamics, as long as the intensity integral in (5.13) can be integrated by parts against . Indeed, omitting for convenience (as this factor is independent of ), we have
[TABLE]
where in the second term we use the fact that is normalized to 1.
In the second term above, the integral in can be readily stripped, however, the first term requires an appropriate change of variables. With help of , we change the variables in the first integral in the right-hand side above:
[TABLE]
Combining the terms, stripping the integrals in , and multiplying by , we obtain the perturbed Kolmogorov equation in the form
[TABLE]
Despite the fact that the type of the external perturbation here is completely different from what was studied previously [22, 30, 29, 57, 27, 28, 41, 40, 42, 43, 44, 49, 61, 62, 66, 70], the forward Kolmogorov equation (5.16) of the perturbed process in (5.8) has a very similar structure to that of a deterministic [62, 76] or Brownian motion [8] perturbation. Namely, observe that there is the unperturbed part (for which we know the steady solution ), and the perturbed part, scaled by a small parameter . Thus, to recover an approximate solution of (5.16), we will proceed in a standard way, which is to expand the solution in -power series around , and recover the first-order approximation via Duhamel’s principle.
5.4. The leading order response formula
We now can look for a solution of (5.16) around the stationary solution in the power series form
[TABLE]
which, upon substitution, yields for
[TABLE]
where by we denote the expression
[TABLE]
Assuming that the initial condition is given via , it is clear that can be found via Duhamel’s principle:
[TABLE]
where is the conditional distribution of the unperturbed Kolmogorov equation in (2.14). Here, the terms under the -integral are known. On the other hand, the transition probability is usually unavailable in the form of an explicit formula.
However, we can circumvent the direct computation of if the average response of a test function is needed. Indeed, observe that, according to (5.1), we can express
[TABLE]
In turn, the integral of against can be expressed via
[TABLE]
where we denote the average response operator via
[TABLE]
5.5. An example of an exactly computable average response operator
In order to obtain the exact formula for the average response operator (5.23) of the Ornstein–Uhlenbeck process, we, first, observe that the part of the integral in (5.23) with can be written as
[TABLE]
which leads to in the form
[TABLE]
Now, for and the Ornstein–Uhlenbeck process (2.15), we have
[TABLE]
For the jump function of the form in (4.3), we have
[TABLE]
Observe that the only difference between (4.4) and the formula above is the presence of in the integral over . As the invariant measure in (2.18) is Gaussian, certain forms of and allow to compute the integral explicitly – for example, if is a polynomial, and is a Gaussian function itself (remember that , since it is part of the statistical intensity of random jumps).
5.6. A practical computational formula for the average response
To compute (5.23) in practice, one can employ the Birkhoff–Khinchin theorem and replace the integral over with the long-term average over time series of the unperturbed system in (2.1) in the same way it was done above in (4.6)–(4.9). This leads to the average response operator in the form of the following time correlation function:
[TABLE]
Observe that, although the external perturbation here is of a completely different nature than what was considered traditionally [22, 30, 29, 57, 27, 28, 41, 40, 42, 43, 44, 49, 61, 62, 66, 70], the form of is similar to that of the classical linear response operator [44, 62, 76]. Namely, the response operator above in (5.28) is also a time correlation function, computed over a long-time series of a solution of the unperturbed system (2.1), except that the response function is multiplied by a different term under the correlation integral in (5.28).
For an efficient computation of the time correlation function above in (5.28), it is necessary for the term to be in the form of an explicit formula. Below we consider a few special cases where the corresponding integral in (5.19) is explicitly computable.
5.7. Special case 1: a set of spatially pre-determined jumps
As a simplest example where the -integral in (5.19) is explicitly computable, we consider the following special case. Assume that the time-homogeneous compound Poisson process (and, subsequently, the time-inhomogeneous compound Poisson process ) produce a finite set of jumps – that is, the jump process may assume only a finite set of values , each with probability , . In this case, the intensity measure of is the same as in (4.10), and, therefore,
[TABLE]
Observe that, as long as and are available in the form of explicit formulas, is also an explicit function of , which allows to compute the time correlation function in (4.9) efficiently. For a practically computable example, we can take of the form (4.12), which leads to
[TABLE]
5.8. Special case 2: a Gaussian jump distribution
As another practical example, let us consider the scenario where has a Gaussian density, with its own mean state vector and covariance matrix as above in (4.15). In addition, we choose the same form for jump function as above in (4.16), which yields the same as in (4.17). This form of and leads to the following expression for in (5.19):
[TABLE]
Here we will further assume that is also given by the Gaussian density in (3.12) (that is, we use the quasi-Gaussian approximation for ).
At this point, the main difference between the current set-up and the previously examined scenario in (4.19) is the additional presence of the function , whose range must lie above zero. Clearly, the most simple choice of leads back to (4.19). Here, we are going to assume that is given via the Gaussian function of the form
[TABLE]
with the -vector and symmetric positive definite matrix being constant parameters. Observe that such form of maximizes the statistical intensity of jumps in the vicinity of the state , with the “width” of the intensity bump described via . This leads to of the form
[TABLE]
with , and given via
[TABLE]
The details of the computation are given in the Appendix A.
As before, observe that the whole expression above is an explicit function of ( is constant, is linear in , and is quadratic in ), which means that it can be evaluated efficiently for the computation of the time correlation function in (5.28).
5.9. Special case 3: a linear combination of Gaussian densities
The natural generalization of the above formula can be easily derived for a scenario where , and are linear combinations of Gaussian densities:
[TABLE]
for some integers , , . It is easy to verify that the corresponding integral in (5.19) becomes
[TABLE]
with
[TABLE]
where and are the mean state and covariance matrix of , respectively. Again, observe that are constants, while and are explicit functions of . This means that the measure integral can be evaluated efficiently in the numerical computation of the time correlation function in (5.28).
6. Accuracy estimates of the leading order average response
Recall that, while the average response formulas in in Sections 3 and 4 are exact, the leading order average response formula for random-time perturbations, described in Section 5, is approximate. The accuracy of the formula in question, specified in (5.21)–(5.23), clearly depends on the magnitude of the scaling parameter .
Recall that in the conventional setting with deterministic external perturbations [44, 62] it is often intuitively clear how large the external perturbation can be; for example, if an external forcing is already present in the dynamical system (say, the constant forcing in the Lorenz 96 system [59, 60, 15], or the vorticity forcing in the barotropic model of the atmosphere [16, 37]), then one can reason that the external perturbation should be much smaller than the already present forcing [42].
On the other hand, observe that, in the present situation, does not determine the magnitude of the external perturbation – in fact, the latter is presumed to be finitely large. Instead, regulates the statistical frequency of (potentially large) external perturbations, and it is thus not immediately clear how small the scaling parameter should be chosen to ensure the accuracy of the leading order average response. In what follows, we address this question to the extent allowed by the generality of problem setting.
6.1. Accuracy of the average response operator of the
Ornstein–Uhlenbeck process
First, we examine a special case where some explicit accuracy estimates can be made. Here, the unperturbed dynamics are given via the Ornstein–Uhlenbeck process (2.15), and the response function is itself.
The perturbed Ornstein–Uhlenbeck process is given via
[TABLE]
Applying the expectation to both sides, we arrive at
[TABLE]
Here, observe that the expectation in the last term is not conditioned on , but rather on . Therefore, in order to split this expectation into the product of independent expectations of and separately, we must set , which renders independent of . In this case, we obtain
[TABLE]
Next, we assume that is of the form (4.16), which further yields
[TABLE]
and, subsequently,
[TABLE]
where we denote
[TABLE]
We thus arrive at the system of ordinary differential equations
[TABLE]
To obtain the exact solution of this system, we further assume, for convenience, that , which, together with the earlier imposed condition , sets the temporal intensity of jumps to jumps per unit of time, on average. Subsequently, the equation above is transformed into
[TABLE]
and its solution is now given via Duhamel’s principle:
[TABLE]
The expectation of the corresponding unperturbed solution is, of course, obtained by setting above:
[TABLE]
with the difference between the two given via
[TABLE]
The exact average response of the mean state is, therefore, given via
[TABLE]
where we use the fact that
[TABLE]
The leading order response, obtained via the average response operator for the same , and , is, however, given via
[TABLE]
The immediate qualitative difference here is that the leading order response, computed via the average response operator, is always bounded regardless of what and are. On the other hand, the exact average response can become unbounded as , as long as and are such that has negative eigenvalues.
Next, assume that is positive definite, so that the exact average response is also bounded as . Then, it is convenient to regard the difference between the infinite-time responses as a metric of accuracy, since it no longer involves the matrix exponential. At the infinite time, with the help of Neumann’s series we have
[TABLE]
Assuming that the series converge, and that , it is clear that the error remains small as long as .
To connect the estimate above with the statistical properties of the dynamics, observe that
[TABLE]
At the same time, observe that, for the Ornstein–Uhlenbeck process in (2.15), the time autocorrelation function of the solution is given via the regression theorem [76]
[TABLE]
where is the covariance matrix (2.19) of the Ornstein–Uhlenbeck process. Therefore, for we have
[TABLE]
that is, the largest eigenvalue of is the statistical autocorrelation time of the solution of the Ornstein–Uhlenbeck process. Thus, the error of the leading order average response operator remains small as long as
[TABLE]
Recalling that, for and , the average time between the jumps is given via , it follows that the average time between the jumps should be much larger than the autodecorrelation time of the solution of the unperturbed system.
6.2. A crude accuracy estimate of the general leading order response
Here we consider a more general setting with the unspecified response function , and the unperturbed dynamics given via (2.1). For convenience, let us, as above for the Ornstein–Uhlenbeck process, assume that , so that the temporal intensity of jumps is time-independent. Just as in the case of the Ornstein–Uhlenbeck process, here it is convenient to look at the response of a function at infinite time. Assuming that the power series in (5.17) converge as , for the exact response and for the leading order approximation , given via (5.23), we have, respectively,
[TABLE]
The condition for the accuracy of the leading order response is, obviously, given via
[TABLE]
which means that the remainder of the infinite series, starting with the second order correction, must be much smaller than the leading order term itself:
[TABLE]
In order to obtain a somewhat simpler requirement for , let us presume that there exists a constant , such that, for all , the following condition holds:
[TABLE]
Then, we can estimate
[TABLE]
Therefore, choosing so as to impose the condition
[TABLE]
is sufficient to ensure that (6.21) holds automatically.
It remains to determine the physical meaning of the parameter , for which we need to formally solve the perturbed Kolmogorov equation in (5.16). For convenience, we rewrite the latter as
[TABLE]
where is the infinitesimal generator of the unperturbed dynamics in (2.2), while is the perturbation component of the perturbed infinitesimal generator in (5.13), given, for , via
[TABLE]
Since the perturbed Kolmogorov equation in (6.26) no longer depends explicitly on time, we can set without loss of generality. Then, in the notations of the introduced operators, we have, via Duhamel’s principle,
[TABLE]
Now, observe that, for a solution of the unperturbed Kolmogorov equation in (2.4), we can write
[TABLE]
which, in turn, leads to
[TABLE]
Applying the formula recursively, we arrive at
[TABLE]
where the integration occurs over the -dimensional simplex with the edge length . Sending and assuming that all are finite, we separate the integrals with the help of Fubini’s theorem:
[TABLE]
where we make use of the identity
[TABLE]
Recalling the definition of in (6.23), we find that, for all ,
[TABLE]
Now, observe that the expressions above can be written as time correlation functions
[TABLE]
where, with the help of the Birkhoff–Khinchin theorem, we can compute the correlation functions via the long-term time averages along a trajectory of the unperturbed system in (2.1):
[TABLE]
Thus, it is clear that is the supremum of all decorrelation times for the correlation functions in (6.36) across all ,
[TABLE]
where is given via (6.32).
Clearly, it is impossible to compute above exactly for a general dynamical system of the form (2.1). However, if the given system is known to have a “decorrelation time scale”, in the sense that all “typical” time correlation functions decay roughly on the same scale, then one can use this time scale as a guidance when estimating suitable values of the parameter for practical computations.
7. Discussion
In the current work, we develop a theory for the average response of a general nonlinear, possibly stochastic, dynamical system to instantaneous external jump perturbations of its state. In real-world physical processes, such perturbations could be a result of an “impulse forcing” – that is, an external forcing in the form of a delta-function. We consider three distinct scenarios of the jump perturbation; the first is where a deterministic jump perturbation is applied at a prescribed time, the second is where a spatially random jump perturbation is applied at a prescribed time, and the third is where the jump perturbation is spatially random, and also occurs at random times (although statistically infrequently).
In all scenarios, we derive the formulas for the average response of a chosen general test function through suitable time correlation functions of the unperturbed dynamics. Throughout the study, we never assume that these jump perturbations are “small” in any reasonable or perceived sense; on the contrary, these perturbations are presumed to be finitely large in all studied scenarios. Below we compare and contrast the developed formulas for the average response to jump perturbations with the response formulas of the classical Fluctuation Dissipation theorem.
7.1. Similarities to the classical FDT
We find that, in the same manner as with the classical Fluctuation Dissipation theorem, the average response of a chosen test function to a finitely large jump perturbation can be computed in terms of a suitable correlation function of the time series of the unperturbed dynamics – in particular, no explicit knowledge of the unperturbed dynamical system is required as long as a sufficiently detailed long-term historical record of its time series is available. Thus, in both frameworks, one could theoretically use the observed time series of the corresponding unperturbed real-world process directly to compute the average response, potentially reducing the need in direct numerical simulations of the (likely complex, large scale and computationally expensive) real-world dynamics.
Also, just like with the classical FDT, a suitable approximation of the probability density of the invariant statistical state of the unperturbed dynamics is necessary to compute the time correlation function. In particular, this results directly in the analog of the quasi-Gaussian FDT formula [15, 14, 16, 62] for the jump perturbations.
7.2. Differences from the classical FDT
First, recall that the classical Fluctuation Dissipation theorem furnishes a leading order response formula [8, 62, 76] to an external perturbation, which, in turn, is accurate only if the magnitude of the perturbation is small enough. In the response framework for jump perturbations, developed here, the accuracy of the average response formulas does not depend on the magnitude of the perturbation jumps. The only leading order approximation used above is with respect to the statistical intensity of randomly triggered jumps – in order for the response formula to be accurate, the jumps should occur statistically infrequently.
Second, in context of the classical Fluctuation Dissipation theorem, the average response is computable for a largely arbitrary additive external perturbation. In the jump-response framework, the average response formula is not necessarily available for an arbitrary jump function – instead, the jump function must satisfy the inverse condition in (3.4) or (4.5), and the inverse needs to be known explicitly. Also, in the case of random jumps, the invariant probability density of the unperturbed dynamics, the statistical intensity of external jump perturbations, and the spatial distribution of the random jumps must all be such that corresponding intensity integral in the average response formula is computable in terms of elementary functions. The latter requirement is necessary for the efficient numerical computation of the time correlation function for the average response.
At the same time, recall that the average response formula for the classical FDT involves the gradient of the invariant probability density of the unperturbed dynamics, which means that any approximation to the latter must also approximate the gradient sufficiently well. In contrast, no differentiation of the invariant probability density is present in any of the jump-response formulas developed here. Potentially, this may mean that less accurate approximations for the invariant probability density could be used without compromising the precision of the predicted response to a considerable extent.
One of the “nice” properties of the classical FDT response formula is that it is linear with respect to the external forcing – usually it is a matrix of certain time correlation functions multiplying the vector of external perturbations. This makes the classical FDT particularly suitable for the “inverse problem”, that is, the computation of the perturbation vector from the known or measured response. Observe, however, that all of the average response formulas, developed above for the external jump perturbations, are nonlinear with respect to the jump function. As a result, the computation of the jump perturbation which caused a known response – for example, the computation of the approximate location and magnitude of the meteorite impact from the measured climate response – is naturally a more challenging and interesting problem.
7.3. Future research directions
Given the introductory nature of the present work, here we feel compelled to point out a few relevant directions of future research.
The basic research directions are those driven by the applicability of the present framework to real-world situations. Arguably, the most limiting restriction in the presented theory of the average response is the need of the jump function to have an explicit inverse (whereas there is no similar requirement in the classical FDT framework). While the linear jump function used in the present work is likely a viable option in a variety of scenarios, there are situations where such a jump function does not work. As an example, one can consider a scenario where the jumps must occur between the states of equal energy in an otherwise energy-conserving system – clearly, the jump function must be nonlinear to ensure such a property. Thus, the study and categorization of possible nonlinear (for example, quadratic) jump functions with explicit inverses should likely be of high priority.
Another question of importance is the already mentioned “inverse problem”, that is, the reconstruction of the perturbation from a known response. While in the classical FDT framework the inverse problem is technically trivial due to the linearity of the response with respect to the external forcing, here we can see that the average response even to a simplest form of the jump perturbation – a constant deterministic jump which occurred at a prescribed time – is inherently nonlinear. Thus, recovering the jump function from the measured response is clearly an important problem in the present framework.
Unlike the classical FDT framework, where the response is obtained by taking the time-convolution of the linear response operator with the external forcing, here observe that, in the scenario with a single jump perturbation at a deterministic time the average response is the time correlation function itself. Naturally, one expects the response to such a perturbation to decay in time, so that the dynamical system eventually returns to is statistical equilibrium state. However, recall that, unlike the classical linear response, here the time correlation function is the exact average response to the external perturbation, provided that the invariant probability state is known exactly. Thus, if the jump-response time correlation function decays to zero, it means that the system relaxes back toward its equilibrium state, and vice versa.
If a dynamical system is strongly mixing (in addition to being ergodic), then the decay of its time correlation functions is guaranteed, and so is the response to a single jump perturbation. However, for an ergodic (and even chaotic) nonlinear dynamical system, being additionally strongly mixing is not a necessary requirement. Moreover, there are empirical examples of chaotic systems which do not appear to be strongly mixing – for example, the Lorenz 96 system [59, 60, 62] in weakly chaotic regimes does not exhibit noticeable decay of time correlation functions at linearly unstable wavenumbers. In such a case, it could be possible that the corresponding dynamical system never “settles down” after a single jump perturbation.
Thus, the question is: are there real-world scenarios in which the geophysical dynamics do not rapidly return back to statistical equilibrium after being subjected to an external jump perturbation? In particular, could the event, which caused the Permian–Triassic extinction, be an example of such a perturbation, such that the occasional planetary climate shifts due to that event continue to occur to this day? The long-time response of the planetary climate dynamics to jump perturbations is clearly a very interesting problem in multiple aspects.
Appendix A Computation of the response for the Gaussian distribution
of jumps
For (4.17), we have
[TABLE]
To obtain (4.19), we first multiply the above expression by (4.15). The resulting form of is given via
[TABLE]
Multiplying the Gaussian densities under the integral, in the argument of the resulting exponential we have
[TABLE]
where , and are given via (4.20). Observing that the integral
[TABLE]
we arrive at
[TABLE]
and, subsequently, to (4.19).
To obtain (5.33), we first multiply in (5.32) by (4.15) and (A.1). The resulting form of is given via
[TABLE]
Multiplying the Gaussian densities under the integral, in the argument of the resulting exponential we have
[TABLE]
where , and are given via (5.34). Observing that the integral
[TABLE]
we arrive at
[TABLE]
and, subsequently, to (5.33).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] R.V. Abramov. Short-time linear response with reduced-rank tangent map. Chin. Ann. Math. , 30B(5):447–462, 2009.
- 2[2] R.V. Abramov. Approximate linear response for slow variables of deterministic or stochastic dynamics with time scale separation. J. Comput. Phys. , 229(20):7739–7746, 2010.
- 3[3] R.V. Abramov. Improved linear response for stochastically driven systems. Front. Math. China , 7(2):199–216, 2012.
- 4[4] R.V. Abramov. A simple linear response closure approximation for slow dynamics of a multiscale system with linear coupling. Multiscale Model. Simul. , 10(1):28–47, 2012.
- 5[5] R.V. Abramov. A simple closure approximation for slow dynamics of a multiscale system: Nonlinear and multiplicative coupling. Multiscale Model. Simul. , 11(1):134–151, 2013.
- 6[6] R.V. Abramov. Linear response of the Lyapunov exponent to a small constant perturbation. Comm. Math. Sci , 14(4):1155–1167, 2016.
- 7[7] R.V. Abramov. A simple stochastic parameterization for reduced models of multiscale dynamics. Fluids , 1(1), 2016.
- 8[8] R.V. Abramov. Leading order response of statistical averages of a dynamical system to small stochastic perturbations. J. Stat. Phys. , 166(6):1483–1508, 2017.