# Convolutionless Nakajima–Zwanzig equations for stochastic analysis in nonlinear dynamical systems

## Abstract

Determining the statistical properties of stochastic nonlinear systems is of major interest across many disciplines. Currently, there are no general efficient methods to deal with this challenging problem that involves high dimensionality, low regularity and random frequencies. We propose a framework for stochastic analysis in nonlinear dynamical systems based on goal-oriented probability density function (PDF) methods. The key idea stems from techniques of irreversible statistical mechanics, and it relies on deriving evolution equations for the PDF of quantities of interest, e.g. functionals of the solution to systems of stochastic ordinary and partial differential equations. Such quantities could be low-dimensional objects in infinite dimensional phase spaces. We develop the goal-oriented PDF method in the context of the time-convolutionless Nakajima–Zwanzig–Mori formalism. We address the question of approximation of reduced-order density equations by multi-level coarse graining, perturbation series and operator cumulant resummation. Numerical examples are presented for stochastic resonance and stochastic advection–reaction problems.

## 1. Introduction

Computing the statistical properties of a stochastic system relies on representing the functional relation between the state variables and the random input processes. Well-known approaches are polynomial chaos [13], multi-element and sparse adaptive probabilistic collocation [4,5], high-dimensional model representations [6], stochastic biorthogonal expansions [7] and separated representations [8,9]. These techniques can provide considerable speed-up in computational time when compared with Monte Carlo (MC) or quasi-MC methods, for low to moderate number of stochastic dimensions.

In this paper, we propose a different approach based on modeling, via deterministic equations, the probability density function (PDF) of low-dimensional quantities of interest, i.e. phase space functions of high-dimensional stochastic systems. The key idea stems from techniques of irreversible statistical mechanics, in particular the Nakajima–Zwanzig–Mori formalism (e.g. [1013]). To introduce the method, let us consider a general nonlinear system in the form dxˆ(t)dt=G(xˆ,ξ,t)andxˆ(0)=x0(ω),1.1where xˆ(t)Rn is a multi-dimensional stochastic process, ξRm is a time-independent random vector (excitation vector), G:Rn+m+1Rn is a Lipschitz continuous (deterministic) function and xˆ0Rn is a random initial state. We shall denote the components of xˆ and G as xˆi and Gi, respectively. The existence and uniqueness of the solution to (1.1) for each realization of x0 and ξ allows us to consider xˆ(t) as a deterministic function of x0 and ξ, i.e. we have the flow map X(t;x0,ξ). In this hypothesis, it is straightforward to obtain the following exact hyperbolic conservation law for the joint (response–excitation) PDF [14,15] of the random vectors xˆ(t) and ξ (see the electronic supplementary material, S1) p(t,a,b)t=L(t,a,b)p(t,a,b),L(t,a,b)=defG(a,b,t)G(a,b,t),1.2where aRn are the phase space coordinates corresponding to xˆ(t), bRm those corresponding to ξ and ∇=(∂/∂a1,…,∂/∂an). Equation (1.2) is equivalent to the Liouville equation of classical statistical mechanics, with the remarkable difference that the phase variables we consider here can be rather general quantities and not simply the positions and the momenta of random particles. For example, they could be the Galerkin coefficients of the solution to stochastic partial differential equations (SPDEs). Early formulations in this direction are due to Edwards [16], Herring [17] and Montgomery [18].

Nonlinear systems in the form (1.1) can lead to complex dynamics, including bifurcations, fractal attractors, multiple stable steady states and transition scenarios. Consequently, the solution to the joint PDF equation (1.2) could be very complex as well, because it relates directly to the geometry of the phase space. For example, it is well known that the PDF of the solution to the Lorenz three-mode problem lies on a fractal attractor whose Hausdorff dimension has been estimated to be about 2.06 [19]. Chaotic states and existence of strange attractors has been well documented for many other systems, such as the Lorenz-84 [20] and the Lorenz-96 [21] models. Even in the much simpler case of the Duffing equation dxˆ1dt=xˆ2anddxˆ2dt=xˆ15xˆ13xˆ250+8cost2,1.3we can have attractors with fractal structure and chaotic phase similarities [22]. This is clearly illustrated in figure 1 where we plot the Poincaré sections of the two-dimensional phase space at different times. Such sections are obtained by sampling 106 initial states from a zero-mean jointly Gaussian distribution, and then evolving them by using (1.3). As the joint PDF of the phase variables is, in general, a high-dimensional compactly supported distribution with a possibly fractal structure, its numerical approximation is a very challenging task [15], especially in long-time integration.

Figure 1.

Complex structure of Poincaré sections of the phase space defined by equation (1.3) at different times. We sample 106 initial states from a zero-mean jointly Gaussian distribution with covariance C11=C22=14, C12=0. The joint PDF of xˆ1 and xˆ2 closely follows the particles and its value is higher where the particle density is higher. Simple statistical properties such as mean and variance are no longer sufficient to describe the stochastic dynamics in the long time.

The statistical description of system (1.1) via the joint PDF equation (1.2), however, is often far beyond practical needs. For instance, we may be interested only in the PDF of one component of the system, e.g. xˆ1(t), or in the PDF of a phase space function g(xˆ). These quantities can be obtained either by integrating out several phase variables from the solution to equation (1.2) (e.g. equation (1.4)) or, more directly, by applying the projection operator method discussed in subsequent sections. This may yield a low-dimensional PDF equation whose solution is more regular than the one obtained by solving directly equation (1.2), and therefore more amenable to computation. The regularization of the reduced-order PDF is due to multi-dimensional integration (marginalization) of the joint response–excitation PDF.

Perhaps, the simplest method to determine reduced-order PDF equations relies on integrating (1.2) with respect to different phase variables and then discarding surface integrals at infinity in phase space. This yields a Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY)-type hierarchy [18] of equations for reduced-order distributions. The lowest order ones are (i=1,…,n) pxˆi(t,ai)=p(t)daidb,1.4where dai=da1⋯dai−1 dai+1⋯dan and db=db1⋯dbm. These densities differ from those used in classical BBGKY theory [23], mainly in that they are not necessarily symmetric under interchanges of different phase space coordinates. For instance, pxˆi(t,ai) is not the same function of ai that pxˆj(t,aj) is of aj, if i and j are different.1 Most of the added complexity to the classical BBGKY theory stems from this lack of symmetry. A related approach, owing to Lundgren [24] and Monin [25], yields a hierarchy of PDF equations involving suitable limits of reduced density functions (see also [14,26,27]). The effective computability of both BBGKY-type and Lundgren–Monin hierarchies associated with equation (1.2) relies on appropriate closure schemes, e.g. a truncation of the hierarchy2 based on a suitable decoupling approximation of the joint PDF. For example, the mean field approximation p(t)=pξ(b)i=1npxˆi(t,ai)t0,1.5where pξ(b) is the joint PDF of the random vector ξ, yields the system of conservation laws pxˆi(t,ai)t=ai[pxˆi(t,ai)hi(t,ai)],i=1,,n1.6and hi(t,ai)=defGi(a,b,t)pξ(b)j=1jinpxˆj(t,aj)dajdb.1.7Equation (1.6) is coupled through the coefficients hi, which have to be computed on-the-fly by using the available PDFs at time t. In the following, we consider a different approach to obtain reduced-order PDF equations for quantities of interest. The approach is based on the Nakajima–Zwanzig–Mori representation of irreversible processes [28,29].

This paper is organized as follows. In §2, we introduce the Nakajima–Zwanzig (NZ) projection operator framework for quantities of interest and derive reduced-order PDF equations based on multi-level coarse graining and operator cumulant resummation. In §3, we provide numerical results obtained by applying the Nakajima–Zwanzig PDF (NZ-PDF) method to simple stochastic systems. In particular, we study stochastic resonance driven by coloured random noise and stochastic advection–reaction phenomena subject to high-dimensional random initial conditions and random reaction rates. The main findings are summarized in §4. We also include Appendix A where we discuss the approximation of the NZ-PDF equation in terms of non-commuting operator cumulants.

## 2. Nakajima–Zwanzig projection operator framework

Introducing a time-independent pair of orthogonal projections P and Q such that P+Q=I (identity operator), we obtain from equation (1.2) Pp(t)t=PL(t)[Pp(t)+Qp(t)]2.1and Qp(t)t=QL(t)[Pp(t)+Qp(t)],2.2where we used the short-hand notation p(t) and L(t) for p(t,a,b) and L(t,a,b), respectively. We shall refer to Pp(t) and Qp(t) as the ‘relevant’ and ‘irrelevant’ parts of the joint PDF p(t). We formally solve equation (2.2) as Qp(t)=G(t,0)Qp(0)+0tG(t,s)QL(s)Pp(s)ds,2.3where the operator G(t,s) (forward propagator3 of the orthogonal dynamics) is formally defined by G(t,s)=defTexpstQL(τ)dτ.2.4In the last equation, T denotes the chronological time-ordering operator (latest times to the left). A substitution of (2.3) into (2.4) yields the following equation, first derived by Nakajima [32], Zwanzig [10,33] and Mori [34], Pp(t)t=PL(t)Pp(t)+PL(t)G(t,0)Qp(0)+PL(t)0tG(t,s)QL(s)Pp(s)ds.2.5In the following, we will refer to equation (2.5) as the NZ-PDF equation. Convergence of Dyson or Magnus series representations of the propagator (2.4) is granted only for short time [30]. For large integration times, we can use the semigroup property of G and split the exponential operator into a product of short-term propagators. An interesting feature of equation (2.5) is its ‘irreversibility’. Roughly speaking, the projected distribution function Pp(t), initially in a certain subspace, leaks out of this subspace so that information is lost. This yields the memory term (time convolution) in equation (2.5).

Tokuyama & Mori [35] have shown how the NZ-PDF equation (2.5) can be transformed into a time-convolutionless form, thus avoiding, at least formally, the memory term. This matter was further discussed by several authors [3638]. To derive the convolutionless form of the NZ-PDF equation, let us consider the formal solution (2.3) and replace p(s) with the solution to (1.2), propagated backward from time t to time st, i.e. p(s)=Z(t,s)p(t)andZ(t,s)=defTexpstL(τ)dτ.2.6In the latter definition, T denotes the anti-chronological ordering operator (latest times to the right). Substituting (2.6) into (2.3) yields Qp(t)=[IΣ(t)]1G(t,0)Qp(0)+[IΣ(t)]1Σ(t)Pp(t),2.7where Σ(t)=def0tG(t,s)QL(s)PZ(t,s)ds.2.8

Equation (2.7) states that the ‘irrelevant’ part of the PDF Qp(t) can, in principle, be determined from the knowledge of the ‘relevant’ part Pp(t) at time t, and from the initial condition Qp(0). Thus, the dependence on the history of the relevant part which occurs in the classical NZ-PDF equation (2.5) has been formally removed by the introduction of the backward propagator (2.6). By using the orthogonal dynamics equation (2.7), we finally obtain the Markovian (time-convolutionless) NZ-PDF equation Pp(t)t=K(t)Pp(t)+H(t)Qp(0),2.9where K(t)=defPL(t)[IΣ(t)]1andH(t)=defPL(t)[IΣ(t)]1G(t,0).2.10Other equivalent forms of the NZ-PDF equation, different than (2.5) and (2.9), can be derived by using suitable ansatzes for the effective propagator of the system (see the electronic supplementary material, S2).

So far, everything that has been said is exact and has led us to the equation of motion (2.9), which is linear. Unfortunately, such equation is still of little practical use, because the exact determination of the operators K(t) and H(t) is as complicated as the solution of equation (1.2). However, the time-convolutionless form (2.9) is a convenient starting point to construct systematic approximation schemes, e.g. by expanding K(t) and H(t) in terms of operator cumulants relatively to suitable coupling constants [11,3941] (see appendix A). We remark that the NZ equations (2.5) and (2.9) and their discrete versions [29] have been extensively used in many fields, e.g. quantum theory of dissipation [11], optimal prediction [28], nonlinear dispersive waves in thermal equilibrium [42] and mesoscopic/atomistic particle methods [43,44].

### (a) Multi-level coarse graining

The orthogonal projection operator Q=IP can be decomposed further by introducing a new pair of orthogonal projections P1 and Q1 such that P1+Q1=I. This allows us to split the governing equation of the orthogonal dynamics (2.2) into the coupled system P1Qp(t)t=P1QL(t)[Pp(t)+P1Qp(t)+Q1Qp(t)]andQ1Qp(t)t=Q1QL(t)[Pp(t)+P1Qp(t)+Q1Qp(t)].Proceeding similarly, we can split the equation for Q1Qp(t) by using a new pair of orthogonal projections P2 and Q2 satisfying P2+Q2=I. This yields P2Q1Qp(t)t=P2Q1QL(t)[Pp(t)+P1Qp(t)+P2Q1Qp(t)+Q2Q1Qp(t)]andQ2Q1Qp(t)t=Q2Q1QL(t)[Pp(t)+P1Qp(t)+P2Q1Qp(t)+Q2Q1Qp(t)].Obviously, we can repeat this process indefinitely to obtain a hierarchy of PDF equations, which generalizes both the NZ as well as the BBGKY frameworks. The advantage of this formulation relatively to the classical NZ formulation relies on the flexibility we have in selecting the set of projections P1, P2, etc. In particular, if the range of Pi is a low-dimensional manifold then the average of an operator product where Pi appears to the left, e.g. P2Q1Q0p(t), represents a low-dimensional quantity. Therefore, in this generalized framework, the joint PDF equation (1.2) is not simply split into the ‘relevant’ and the ‘irrelevant’ parts by using P and Q, but the dynamics of the irrelevant part is decomposed further in terms of an ensemble of projections. This allows us to coarse-grain the orthogonal dynamics further in terms of low-dimensional quantities. Multi-level coarse graining based on sequences of projections is closely related to the method of subdynamics originally introduced by Prigogine et al. [45] (see also [4648]). Such method has proved to be very powerful for deriving kinetic equations in non-equilibrium statistical mechanics [49].

### (b) Approximations to the Nakajima–Zwanzig probability density function equation

Most approximation schemes for NZ-PDF equations or BBGKY-type hierarchies rely on the identification of some small quantity that serves as basis for perturbation expansion, e.g. the particle density for Boltzmann equations [23], the coupling constant or correlation time for effective Fokker–Planck equations [5052], the Kraichnan absolute equilibrium distribution for turbulent inviscid flows [18,53] (see also [54]) or high-order Fourier coefficients in semi-discrete forms of PDEs [28]. One of the most stubborn impediments for the development of a general theory of reduced-order PDF equations has been the lack of such readily identifiable small parameters. Most of the work that has been done so far refers to the situation where the operator L(t) in equation (1.2) can be decomposed as L(t)=L0+σL1(t),2.11where L0 depends only on the relevant variables of the system, σ is a positive real number (coupling constant in time-dependent quantum perturbation theory) and the norm σL1(t)∥ is somehow small. Note that this does not necessarily mean that σ is itself a small quantity. In fact, the norm of L1(t) could be small as well, because of fast dynamics (small correlation time processes) or other properties of the unresolved phase variables. By using the interaction representation of quantum mechanics [30,31], it is quite straightforward to obtain from (2.9) and (2.11) an effective closure approximation that is valid, e.g. for small noise, small correlation or fast–slow dynamics. One way to do this is to expand the operators K(t) and H(t) in (2.10) in a cumulant series, e.g. in terms of Kubo–Van Kampen operator cumulants [3941,55] (see appendix A). Any finite-order truncation of such series represents an approximation to the exact NZ-PDF equation. In particular, the approximation obtained by retaining only the first two operator cumulants is known as Born approximation in quantum field theory. We remark that from the point of view of perturbation theory, the convolutionless form (2.9) has distinctive advantages over the usual convolution form (2.5). In particular, in the latter case a certain amount of rearrangement is necessary to obtain an expression which is correct up to a certain order in the coupling parameter [56]. More recently, Chorin et al. [28,57] proposed various closure models to deal with situations where there is no clear separation of scales between the resolved and unresolved dynamics.

### (c) Quantities of interest and projection operators

Suppose that we are interested in the PDF of a phase space function g:Rn+1Rl in the form u(t)=gxˆ(t),t,(Quantity of Interest)2.12where xˆ(t) is the stochastic process solving equation (1.1). For example, g could be the hyper-plane defined by the series expansion of the solution to a SPDE, i.e. u(t,x)=j=1nxˆj(t)lj(x),2.13where lj(x) are spatial basis functions (x here are spatial coordinates). By using the NZ projection operator framework, we can obtain the exact evolution equation for the PDF of any phase space function in the form (2.12). The starting point is the Liouville equation for joint PDF of the state variables u, xˆ and ξ, which is similar to equation (1.2), and can be derived by augmenting system (1.1) with additional equations defining the time evolution of (2.12) dudt=gt+gxˆGNow, let aRl and bRn+m be the sets of phase space coordinates corresponding, respectively, to the relevant and the irrelevant variables of the augmented system. In this case, a corresponds to u, while b corresponds to the pair (xˆ,ξ). Let p(t) the joint PDF of u, xˆ and ξ. We define the orthogonal projection operator Pp(t)=defφ(b)p(t)db,2.14where φ is a multi-dimensional PDF with appropriate support4 in Rn+m. Clearly, P is an orthogonal projection since, by the assumptions on φ(b), we have P2=P. Note also that integrating Pp(t) with respect to b yields the PDF of the quantity of interest (2.12) pu(t)=Pp(t)db.The corresponding evolution equation can be easily obtained by substituting the projection operator (2.14) into equation (2.5) or (2.9), and then integrate with respect to b. Clearly, we have a considerable freedom in selecting the projection operator P, i.e. the PDF φ(b). If the relevant variables of the initial condition are separable from the irrelevant ones, i.e. if p(0)=pa(0,a)pb(0,b), then an appropriate choice for φ(b) in (2.14) is the joint PDF of the irrelevant variables, i.e. we can set φ(b)=pb(0,b). In this case, the joint PDF p(0) is in the range of P, i.e. Pp(0)=p(0). This implies that the initial condition term in the NZ-PDF equation is identically zero since Qp(0)=(IP)p(0)=0. In general, the joint PDF of the initial state is not separable. This happens, for example, when we augment system (1.1) with additional equations for the quantities of interest (2.12). In these cases, the projection operator (2.14) yields an initial condition term H(t)Qq(0) in equation (2.9) that does not vanish, and therefore we have to take it into account, e.g. by using operator cumulant series expansions. As an example, by applying the NZ-PDF framework to equation (1.3), we obtain the following exact law governing the PDF of xˆ1(t)5 (figure 2) pxˆ1(t)t=L(t)[IΣ(t)]1pxˆ2(0,a2)da2pxˆ1(t),2.15where the integral is with respect to the PDF of xˆ2(0) (assumed statistically independent of xˆ1(0)), Σ(t) was defined in (2.12) and the operators L(t) and P are, respectively L(t)=a2a1+a2a1+5a13+a2508cost2()andP=pxˆ2(0,a2)()da2.

Figure 2.

Reduced-order PDF dynamics of the Duffing system (1.3). Shown are non-parametric kernel density estimates [58] of the PDF of xˆ1(t) at the same times as in figure 1. The exact evolution equation of pxˆ1(t) is (2.15). (Online version in colour.)

We emphasize that it is rather difficult to determine a computable approximation to the exact NZ-PDF equation (2.15). The main reason is that xˆ1 and xˆ2 in (1.3) have similar dynamical properties and the same order of magnitude. This implies that it is not easy to integrate out the dynamics associated with the phase variable xˆ2, unless we have additional useful information. In general, the problem of integrating out a one-phase variable from a two-dimensional nonlinear system evolving from a random initial state is known as ‘Mori's problem’ [59,60]. Unfortunately, such problem does not have an easily computable solution [61].

## 3. Results and discussion

In this section, we provide numerical results obtained by applying the NZ-PDF approach to simple stochastic systems. In particular, we study stochastic resonance driven by coloured random noise and stochastic advection–reaction phenomena subject to high-dimensional random initial conditions and random reaction rates.

### (a) Stochastic resonance driven by coloured noise

In bistable systems, the cooperation between random noise and deterministic periodic signals of small amplitude can yield a phenomenon known as stochastic resonance [62,63]. Roughly speaking, random noise can enhance significantly the transmission of the deterministic signals, provided we tune system's parameters appropriately. The mechanism that makes this possible is explained in figure 3, with reference to the system dxˆ(t)dt=dU(xˆ(t))dxˆ+σf(t,ξ)+ϵcos(Ωt),xˆ(0)=x0(ω).3.1Specifically, the potential U(xˆ) is assumed to be in the form U(xˆ)=def11+xˆ2μxˆ22+νxˆ44μ,ν>0.3.2This function is a modification of the classical double-well potential [52,64] because it has two symmetric minima at xˆ1,2=±1+2μ/ν1, a local maximum at xˆ=0 and it grows quadratically as xˆ±. A similar system has been studied in [64,65] by using the unified coloured noise approximation [6669]. The small time-periodic (deterministic) signal ϵcos(Ωt) in (3.1) can be included in the potential U(xˆ), yielding the time-dependent modulation U(xˆ,t)=defU(xˆ)ϵxˆcos(Ωt).3.3Correspondingly, system (3.1) can be rewritten as dxˆ(t)/dt=U(xˆ,t)/xˆ+σf(t,ξ). Modulation (3.3) is at the basis of the stochastic resonance phenomenon, which is classically studied by using white-noise excitation and Fokker–Planck theory [62]. Here, we employ the convolutionless NZ approach to determine the evolution of the PDF of xˆ(t) in the case where f(t,ξ) in (3.1) is a zero-mean Gaussian process with correlation function C(t,s). Specifically, we consider

• (i) Exponentially correlated noise6 C(t,s)=12τe|ts|/τ(τ is the noise correlation time)3.4

• (ii) Fractional Brownian noise [70] C(t,s)=12(|t|2h+|s|2h|ts|2h)0<h<1(Hurst index).3.5

Figure 3.

Evidence of stochastic resonance. We study system (3.1) with parameters μ=10, ν=3, Ω=2 and ϵ=0.2 subject to weakly coloured random noise (τ=0.01) of different amplitudes: (a) σ=0.0; (b) σ=0.2; (c) σ=0.4 and (d) σ=0.8. Each plot shows only one solution sample. At the low values of the noise, the average residence time in the two states is much longer than the driving period. However, if we increase the noise level to σ=0.8 (c), then we observe almost periodic transitions between the two meta-stable states. In most cases, we have a jump from one state to the other and back again approximately once per modulation period.

In both cases, a Karhunen–Loève series expansion of f(t;ξ) can be determined (e.g. [71]). In particular, the Karhunen–Loève eigenvalue problem for the exponential correlation function admit an analytical solution [72] (see the electronic supplementary material, S3.1).

The exact evolution equation for the PDF of xˆ(t), denoted as pxˆ(t), can be obtained by applying the convolutionless projection operator method of §2, with projection operator defined as P=defpξ(b)()db,3.6pξ(b) being the joint PDF of the random vector ξ. Such NZ-PDF equation is a linear partial differential equation involving derivatives of infinite-order. If we consider the Born approximation, i.e. if we expand the propagator of pxˆ in terms of Kubo–Van Kampen operator cumulants and truncate the expansion at the second order (see appendix A), we obtain pxˆ(t)t=L0pxˆϵcos(Ωt)pxˆ(t)a+σ20tC(t,s)ae(ts)L0ae(st)L0dspxˆ(t),3.7where L0=defa2μa2νa3νa52(1+a2)2().3.8

It can be shown that equation (3.7) is in the form of an advection–diffusion equation, with diffusion coefficient that depends on the phase variable a as well as on time t (see the electronic supplementary material, S3.3). The specific form of the diffusion coefficient has been the objective of extensive investigation [13,52,73,74]. The rationale behind the Born approximation (3.7) is that higher order operator cumulants can be neglected. This happens, in particular, if both ϵ and σ are small. However when the noise is coloured, it is in general quite difficult to estimate the magnitude of the contribution of higher order operator cumulants. For example, Faetti et al. [13] have shown that for ϵ=0 the correction to (3.7) due to the fourth-order cumulant is O(σ4τ2) for exponentially correlated Gaussian noise (see eqn (4.18) in [13]). Note that if the noise correlation time τ goes to zero (white-noise limit), then equation (3.7), with C(t,s) defined in (3.4), consistently reduces to the classical Fokker–Planck equation.7

Next, we study the transient dynamics of pxˆ(t) within the time interval [0,3]. To this end, we consider the following set of parameters μ=1, ν=1, Ω=10 and ϵ=0.5, leading to a slow relaxation to statistical equilibrium. This allows us to study the transient dynamics more carefully and compare results with MC simulation. This is done in figure 4 for exponentially correlated noise. There it is seen that random noise with very small amplitude does not influence significantly the response of the system. In fact, the Gaussian ensemble of initial states is mainly advected by the operator L0, yielding accumulation nearby the meta-stable states ±31. For larger noise amplitudes, the probability of switching between the meta-stable states increases and therefore the strong bi-modality observed in figure 4a is attenuated. In figure 5, we study the effects of the noise correlation time τ on the transient dynamics and compare results with MC. It is seen that random noise with larger correlation time inhibits the switching rate between the two meta-stable states. This means that stochastic resonance is suppressed with increasing noise colour, in agreement with the findings of Hänggi et al. [65].

Figure 4.

Stochastic resonance simulation. Time snapshots of the PDF of xˆ(t) as predicted by equation (3.7) (continuous lines) and MC simulation (105 samples) (dashed lines). The Gaussian noise here is exponentially correlated with correlation time τ=0.1 and amplitude σ=0.1 (a) and σ=0.5 (b). The Karhunen–Loève series representation of such noise requires 280 Gaussian random variables to achieve 99% of the correlation energy in the time interval [0,3]. (Online version in colour.)

Figure 5.

Suppression of stochastic resonance with increasing noise colour. It is seen that as we increase the correlation time τ from 0.01 (a) to 0.4 (b), the probability of switching between the meta-stable states ±31 decreases. This is even more emphasized in figure (c), where we study stochastic resonance driven by fractional Brownian motion with Hurst index h=0.1 and 0.6. (Online version in colour.)

The response of system (3.1) to fractional Brownian motion of small amplitude can be computed similarly, by substituting the non-stationary covariance function (3.5) into equation (3.7). The numerical results are shown in figure 5c, where it is seen that fractional Brownian motion inhibits stochastic resonance in a more pronounced way than exponentially correlated processes. Also, larger Hurst indices induce larger diffusion effects in the PDF at later times.

The connection between the statistical properties of random noise and the structure of the NZ-PDF equation can be further generalized to finite-amplitude perturbations (e.g. [75,76]).

Let us consider the advection–reaction equation ut+V(x)u=[κ0(x)+σκ1(x;ξ)]R(u),3.12where V(x) is a deterministic velocity field, R(u) is a nonlinear reaction term and κ1(x;ξ) is a zero-mean random perturbation of the deterministic reaction rate κ0(x). In [77,78], we studied equation (3.9) by using the joint response–excitation PDF method as well as the large-eddy-diffusivity closure [79]. Here, we consider a different approach based on the NZ-PDF equation. To this end, let us assume that σ in (3.9) is reasonably small and that the concentration field u is statistically independent of κ1 at initial time. Under these hypotheses, by using the projection operator P=defpξ(b)()db3.13and the perturbation analysis presented in appendix A, we obtain the following Born approximation to the exact NZ-PDF equation for the concentration field (see the electronic supplementary material, S4) pu(t)t=L0pu(t)+σ20tκ1esL0κ1esL0dsF2pu(t),3.14where L0=defκ0(x)FV(x)andF=defR(a)a+R(a)a,3.15and the average 〈⋅〉 is defined to be a multi-dimensional integral with respect to the joint PDF of ξ. In figure 6, we plot a few realizations of the concentration field solving equation (3.9) in one-space dimension, where we have set V(x)=32+12esin(x)+cos(x)/2cos(x),3.16 κ0(x)=3225(esin(x)/2+cos(x)),3.17 R(u)=120(2e-ecos(2πu)2-ecos(2πu)),3.18 κ1(x;ξ)=1mΣj=1mξj(1)(ω)sin(jx)+Σj=1mξj(2)(ω)cos(jx)3.19and u(x,0;ζ)=120sin(x)+120m1j=1m1ζj(1)(ω)sin(jx)+j=1m1ζj(2)(ω)cos(jx),3.20where {ξj(i)} and {ζj(i)} being two sets of uncorrelated normal random variables. In figure 7, we show the time snapshots of the PDF of the concentration field as predicted by the NZ-PDF equation (3.11) and compare the results with MC simulation. We note that, in this case, the Born approximation provides accurate results for a quite large degree of perturbation.

Figure 6.

Stochastic advection–reaction. Samples of the concentration field for different number of dimensions in the initial condition (3.17). Specifically, we consider m1=4 (a), 10 (b), 20 (c), 40 (d) 80 (e). The perturbation amplitude in the reaction rate is set to σ=0.5, and we consider m=5 (10 random variables) in (3.16). (Online version in colour.)

Figure 7.

Stochastic advection–reaction simulation. Time snapshots of the concentration PDF predicted by the NZ-PDF equation (3.11) (a) and PDF dynamics at x=1 (b). Shown is the comparison between NZ-PDF solution and a non-parametric kernel density estimation [58] of the PDF based on 103 MC samples. The zero-order approximation is obtained by neglecting the term multiplying σ2 in equation (3.11). We set m=5 in the random reaction rate (3.14) (10 random variables) and m1=80 in the initial concentration field (3.17) (160 random variables). (Online version in colour.)

## 4. Summary

We presented a general methodology to derive exact PDF equations for quantities of interest in nonlinear stochastic dynamical system with parametric uncertainty. The key idea stems from techniques of irreversible statistical mechanics, and it relies on defining suitable phase space functions and corresponding projection operators, yielding formally exact reduced-order PDF equations. The effective numerical simulation of such reduced-order PDF equations relies on appropriate approximations. Most of the schemes proposed so far are based on the identification of a small quantity that serves as a basis for perturbation expansion in terms of generalized operator cumulants. Approximation of reduced-order PDF equations beyond perturbation is still an open question. Any advancement in this direction would open the possibility to tackle high-dimensional stochastic systems where the relevant and irrelevant phase variables have the same order of magnitude and similar dynamical properties. There have been recent work in this direction, e.g. by Chorin et al. [28,57] and Darve et al. [29], in the context of model reduction of deterministic autonomous dynamical systems. In particular, various models such as the t-model [80] and its renormalized version [81] have been proposed to deal with situations where there is no clear separation of scales between the resolved and the unresolved dynamics.

## Funding statement

This work was supported by OSD-MURI grant no. FA9550-09-1-0613, by DOE grant no. DE-SC0009247 though the Collaboratory on Mathematics for Mesoscopic Modelling of Materials (CM4) and NSF/DMS-1216437.

## Appendix A. Operator cumulant approximation

Let us assume that the operator L(t) in equation (1.2) can be decomposed as L(t)=L0+σL1(t),A1where L0 is a time-independent operator depending only on the relevant phase variables of the system, while σ is a positive coupling constant, for example the amplitude of an external random noise. In this case, we can integrate out exactly the dynamics associated with L0 first so as to circumscribe the approximation problem to L1(t). This can be carried out by means of a preliminary time-dependent transformation. In quantum mechanics, such transformation corresponds to a new evolution picture, e.g. the interaction or the adiabatic picture [30,31,82]. In particular, in the interaction picture, we first define the auxiliary field q(t)=defetL0p(t)A2and then rewrite equation (1.2) as qt=σN(t)q(t)q(0)=p(0),A3where N(t)=defetL0L1(t)etL0.A4By using the assumption that L0 depends only the relevant variables of the system (i.e. it commutes with the projection P defined in (2.14), we have from (A2) Pp(t)t=L0Pp(t)+etL0Pq(t)t.A5At this point, by following the same steps that led us to equation (2.9), we represent the evolution of Pq(t) in terms of an exact convolutionless NZ equation Pq(t)t=Kˆ(t)Pq(t)+Hˆ(t)Qp(0),A6where Kˆ(t)=defσPN(t)[IσΣˆ(t)]1A7and Hˆ(t)=defσPN(t)[IσΣˆ(t)]1Gˆ(t,0),A8and Σˆ(t)=def0tGˆ(t,s)QN(s)PZˆ(t,s)ds,A9 Gˆ(t,s)=defTexpσstQN(τ)dτA10and Zˆ(t,s)=defTexpσstN(τ)dτ.A11Equations (A6)–(A11) can be used as a starting point of a perturbation theory.

(a) Perturbation series

Let us assume that the norm of the operator σN(t) appearing in (A3) is small and expand the generator (A7) and the inhomogeneity (A8) in the formal power series Kˆ(t)=n=1σnKˆn(t)andHˆ(t)=n=1σnHˆn(t).A12To determine the operators Kˆn and Hˆn, we first expand [IσΣˆ(t)]1 in equations (A7)–(A8) in a Neumann series8 as [IσΣˆ(t)]1=k=0σkΣˆ(t)k.A13We also expand in power series both Σˆ(t) and Gˆ(t,0) in equations (A9) and (A10) Σˆ(t)=n=0σnΣˆn(t)andGˆ(t,0)=n=0σnGˆn(t).A14A substitution of (A14) into (A13) yields [IσΣˆ]1=I+σΣˆ0+σ2[Σˆ1+Σˆ02]+A15and [IσΣˆ]1Gˆ=I+σ[Σˆ0+Gˆ1]+σ2[Σˆ02+Σˆ1+Gˆ2+Σˆ0Gˆ1]+A16where the coefficients are obtained as Σˆ0(t)=0tdt1QN(t1)P,A17 Σˆ1(t)=0tdt10t1dt2[QN(t1)QN(t2)PQN(t2)PN(t1)],A18 Gˆ0(t)=I,A19 Gˆ1(t)=0tdt1QN(t1)A20 Gˆ2(t)=0tdt10t1dt2QN(t1)QN(t2),A21Finally, by substituting (A14)–(A21) into (A7) and (A8), we obtain the following coefficients of the power series expansions (A12): Kˆ1(t)=PN(t),A22 Kˆ2(t)=0tdt1[PN(t)N(t1)PN(t)PN(t1)],A23 Kˆ3(t)=0tdt10t1dt2[PN(t)QN(t1)QN(t2)PN(t)QN(t2)PN(t1)],A24 Hˆ1(t)=PN(t),A25 Hˆ2(t)=0tdt1[PN(t)N(t1)PN(t)PN(t1)+PN(t)QN(t1)],A26 Hˆ3(t)=0tdt10t1dt2[PN(t)QN(t1)QN(t2)PN(t)QN(t2)PN(t1)]+0tdt10t1dt2PN(t)QN(t1)QN(t2),A27Any finite-order truncation of the series in (A12) represents a closure approximation of the NZ equation. In particular, the approximation obtained by retaining only the first two terms in the expansion of Kˆ and Hˆ is known as Born approximation in quantum field theory.

(b) Convergence

The series expansions we have considered so far are based on the hypothesis that the norm of the operator σN(t) is somehow small. Although this is true in some instances, we cannot obviously expect that such assumption is satisfied in general. In addition, the smallness of the norm σN(t)∥ is not the only requirement for the convergence of the perturbation theory we have discussed so far. In fact, we also have a finite radius of convergence (in time) in the Dyson and Magnus series representations of the propagators (A10) and (A11).9 This issue can be overcome by using the semigroup property of Gˆ and Zˆ. For example, Gˆ(tn,t0) can be represented as a superimposition of short-term propagators, for which convergence is granted, i.e. Gˆ(tn,t0)=Gˆ(tn,tn1)Gˆ(t2,t1)Gˆ(t1,t0).A28

(c) Kubo–Van Kampen, Waldenfels and Speicher operator cumulants

Let U(t,s) be the forward propagator of the solution to equation (A3), i.e. q(t)=U(t,s)q(s)ts0.A29Such propagator satisfies the operator equation dU(t,s)dt=σN(t)U(t,s),U(s,s)=I.A30In many applications, the initial condition of the joint PDF p(0) is also separable, i.e. we have p(0,a,b)=pa(0,a)pb(0,b),A31where a and b are the relevant and irrelevant phase variables of the system (pa and pb are, respectively, the joint PDFs of the relevant and irrelevant variables of the system). This leads us to a formal reduced-order solution to equation (A5) in the form pa(t,a)=etL0U(t,0)pa(0,a),U(t,0)=defU(t,0)pb(0,b)db.A32This equation defines the exact time dynamics of the relevant part of the joint PDF p(t,a,b) in terms of the averaged (effective) propagatorU(t,0)〉. Only rarely this operator can be calculated explicitly. The tempting idea to average the terms of the Dyson series, U(t,s)=I+n=1σnstdt1st1dt2stn1dtnN(t1)N(tn),A33is, in general, not appropriate if one needs 〈U(t,s)〉 for large (ts). In fact, there appear the so-called secular terms growing like (ts)n if the nth-order correlation 〈N(t1)⋯N(tn)〉 does not decay rapidly enough. As pointed out by Kubo [55,85], Hegerfeldt & Schulze [39] and Neu & Speicher [41], this problem sometimes can be overcome by imposing an ansatz for 〈U(t,s)〉, in particular a differential equation characterizing the evolution of the averaged propagator. This yields a resummation of the various correlations appearing in equation (A33), i.e. the introduction of different types of operator cumulants. Known ansatzes are U(t,s)dt=F(t,s)U(t,s),A34 U(t,s)dt=σN(t)U(t,s)+stA(t,τ)U(τ,s)dτA35 and U(t,s)dt=n=0stdt1st1dt2stn1dtnSn(t,t1,,tn)U(t,t1)U(t1,t2)U(tn,s).A36It can be shown that (A34) is equivalent to (A33) if we choose F(t,s)=σN(t)+n=1σn+1stdt1st1dt2stn1dtnN(t)N(t1)N(t2)N(tn)K,A37where 〈N(t)N(t1)N(t2)⋯N(tn)〉K are Kubo–Van Kampen (K-cumulants). By using the shorthand notation 〈01⋯k〉=〈N(t)N(t1)⋯N(tk)〉, the first three K-cumulants can be written as N(t)K=0,A38 N(t)N(t1)K=0101A39 and N(t)N(t1)N(t2)K=012012021012+012+021.A40Note that these cumulants coincide with those appearing in equations (A22)–(A24), if we define the average as in (A32) and we select the projection operator as P()=defpb(0,b)()db.A41The second choice (A35), introduced by Terwiel [56] (see also [39,55]), yields a solution in the form (A33) if we select A(t,τ)=σ2N(t)N(s)W+n=1σn+2τtdt1τt1dt2τtn1dtnN(t)N(t1)N(tn)N(τ)WA42where 〈N(t)N(t1)N(t2)⋯N(tn)〉W are Waldenfels cumulants (W-cumulants) N(t)W=0,A43 N(t)N(t1)W=0101A44and N(t)N(t1)N(t2)W=012012012+012.A45Note that differently from the K-cumulants, the W-cumulants are time-ordered (compare (A38)–(A40) with (A43)–(A45)). The last ansatz (A36), introduced by Neu & Speicher [41], yields another type of cumulants, i.e. a specific class of non-crossing cumulants. All these representations can be obtained by setting a specific ordering prescription in the definition of the propagator U(t,s) [55,85,86]. This also justifies the existence of relations between different types of cumulants. It is important to emphasize that all expansions we have discussed so far are completely equivalent if we consider an infinite number of terms in the series. However, the convergence properties of truncated series could be significantly different and could also lead to positivity problems [39], p. 704.

## Footnotes

• 1 In the classical BBGKY framework, the phase coordinates of the systems are positions and momenta of identical particles. Therefore, the reduced-order multi-point densities are invariant under interchanges of phase space coordinates of the same type, e.g. positions or momenta.

• 2 A low-order truncation of the BBGKY hierarchy for particle systems can be used to derive the Boltzmann equation [23].

• 3 The propagator G(t,s) is a linear operator (semi-group) that admits many different representations, for instance Magnus series or Dyson series [30,31].

• 4 The distribution function φ(b) could be, for example, a uniform distribution in a domain B, i.e. φ(b)=1/|B|, where |B| denotes the volume of the set. It could also be a multivariate Gaussian distribution or the equilibrium distribution of the irrelevant phase variables of the system.

• 5 In this case, the mapping (2.12) is trivial and it reduces to g(xˆ1,xˆ2,t)=xˆ1.

• 6 In the limit τ0, f(t;ξ) becomes Gaussian white noise with correlation function δ(ts).

• 7 The proof is simple, and it relies on the limits limτ00t12τes/τds=12andlimτ00t12τes/τskds=0kN. In fact, these equations allow us to conclude that limτ00t12τes/τaesL0aesL0ds=limτ00t12τes/τds2a2=122a2, i.e. for τ0, equation (3.7) reduces to the Fokker–Planck equation pxˆt=L0pxˆϵcos(Ωt)pxˆa+σ222pxˆa2.

• 8 The convergence radius of the Neumann series (A13) is related to magnitude of σΣˆ(t), i.e. to the magnitude of σN(t)∥ (see [83]).

• 9 If L1(t) in equation (A1) is time-independent then the propagators (A10) and (A11) are standard exponential operators for which uniform convergence is granted on the whole real axis [83,84].