Generalized Langevin equations for systems with local interactions

TL;DR

This paper introduces a novel method to approximate the memory integral in generalized Langevin equations for high-dimensional nonlinear systems with local interactions, enabling accurate computation of statistical properties.

Contribution

It develops a new approximation technique for the Mori-Zwanzig memory integral using Faber operator series and an exact combinatorial algorithm, advancing the analysis of complex systems.

Findings

Effective computation of auto-correlation functions demonstrated

Accurate modeling of intermediate scattering functions achieved

New stochastic process representation for equilibrium systems

Abstract

We present a new method to approximate the Mori-Zwanzig (MZ) memory integral in generalized Langevin equations (GLEs) describing the evolution of smooth observables in high-dimensional nonlinear systems with local interactions. Building upon the Faber operator series we recently developed for the orthogonal dynamics propagator, and an exact combinatorial algorithm that allows us to compute memory kernels from first principles, we demonstrate that the proposed method is effective in computing auto-correlation functions, intermediate scattering functions and other important statistical properties of the observable. We also develop a new stochastic process representation of the MZ fluctuation term for systems in statistical equilibrium. Numerical applications are presented for the Fermi-Pasta-Ulam model, and for random wave propagation in homogeneous media.

Figures

Tables

Table 1. Table 1: Series expansions of the Mori-Zwanzig memory kernel ( 14c ). In this Table, J q subscript 𝐽 𝑞 J_{q} is the q 𝑞 q th Bessel function of the first kind, c 0 subscript 𝑐 0 c_{0} and c 1 subscript 𝑐 1 c_{1} are real numbers defining the recurrence relation of the Faber polynomials { F 0 , … ​ F n } subscript 𝐹 0 … subscript 𝐹 𝑛 \{F_{0},\dots F_{n}\} , M 𝑀 M is the number of phase space functions u i ​ ( 𝒙 ) subscript 𝑢 𝑖 𝒙 u_{i}(\bm{x}) (quantities of interest), and G i ​ j − 1 subscript superscript 𝐺 1 𝑖 𝑗 G^{-1}_{ij} is the inverse of the Gram matrix ( 14a ). See [ 68 ] for additional details.

Type	Temporal modes $g_q\left(t\right)$	Coefficients $M_{ijq}$
MZ-Dyson	$\frac{t^q}{q!}$	$\sum _{k=1}^M{G}_{jk}^{-1}{⟨u_k\left(0\right),ℒ{\left(𝒬ℒ\right)}^q𝒬ℒu_i\left(0\right)⟩}_{\rho }$
MZ-Faber	$e^{t{c}_0}\frac{J_q\left(2t\sqrt{-c_1}\right)}{{\left(\sqrt{-c_1}\right)}^q}$	$\sum _{k=1}^M{G}_{jk}^{-1}{⟨u_k\left(0\right),ℒF_q\left(𝒬ℒ\right)𝒬ℒu_i\left(0\right)⟩}_{\rho }$

Equations

$$\frac{d\mathbf{x}}{dt}=\mathbf{F}(\mathbf{x}),\mathbf{x}(0)={\mathbf{x}}_0,$$

$$\mathbf{u}:\mathcal{M}$$

$$\mathbf{x}$$

$$\mathbf{u}(\mathbf{x}(t,{\mathbf{x}}_0))=e^{t\mathcal{L}({\mathbf{x}}_0)}\mathbf{u}({\mathbf{x}}_0),\mathcal{L}({\mathbf{x}}_0)=\sum _{k=1}^N{F}_k({\mathbf{x}}_0)\frac{\partial }{\partial x_{0k}}.$$

$$e^{t\mathcal{L}}=e^{t\mathcal{Q}\mathcal{L}}+{\int }_0^t{e}^{s\mathcal{L}}\mathcal{P}\mathcal{L}{e}^{(t-s)\mathcal{Q}\mathcal{L}}ds$$

$$\frac{d}{dt}{e}^{t\mathcal{L}}=e^{t\mathcal{L}}\mathcal{P}\mathcal{L}+e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}+{\int }_0^t{e}^{s\mathcal{L}}\mathcal{P}\mathcal{L}{e}^{(t-s)\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}ds.$$

$$\frac{\partial }{\partial t}{e}^{t\mathcal{L}}\mathbf{u}(0)$$

$$\frac{\partial }{\partial t}\mathcal{P}{e}^{t\mathcal{L}}\mathbf{u}(0)=\mathcal{P}{e}^{t\mathcal{L}}\mathcal{P}\mathcal{L}\mathbf{u}(0)+{\int }_0^t\mathcal{P}{e}^{s\mathcal{L}}\mathcal{P}\mathcal{L}{e}^{(t-s)\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}\mathbf{u}(0)ds.$$

$$⟨f,g{⟩}_{\rho }={\int }_{\mathcal{M}}f(\mathbf{x})g(\mathbf{x})\rho (\mathbf{x})d\mathbf{x}f,g\in H$$

$$G_{ij}={\int }_{\mathcal{M}}{u}_i(\mathbf{x})u_j(\mathbf{x})\rho (\mathbf{x})d\mathbf{x}.$$

$$\mathcal{P}f=\sum _{i,j=1}^M{G}_{ij}^{-1}⟨u_i(0),f{⟩}_{\rho }{u}_j(0),f\in H.$$

$${\rho }_{eq}(\mathbf{x})=\frac{1}{Z}{e}^{-\beta \mathcal{H}(\mathbf{x})},$$

$$\frac{d\mathbf{u}(t)}{dt}=\mathbf{\Omega }\mathbf{u}(t)+{\int }_0^t\mathbf{K}(t-s)\mathbf{u}(s)ds+\mathbf{f}(t),$$

$$\frac{d{u}_i(t)}{dt}=\sum _{j=1}^M{\Omega }_{ij}{u}_j(t)+\sum _{j=1}^M{\int }_0^t{K}_{ij}(t-s)u_j(s)ds+f_i(t).$$

$$G_{ij}$$

$${\Omega }_{ij}$$

$$K_{ij}(t)$$

$$\mathbf{f}(t)$$

$$\frac{d}{dt}\mathcal{P}\mathbf{u}(t)=\mathbf{\Omega }\mathcal{P}\mathbf{u}(t)+{\int }_0^t\mathbf{K}(t-s)\mathcal{P}\mathbf{u}(s)ds.$$

$$\frac{d}{dt}{C}_{ij}(t)=\sum _{k=1}^M{\Omega }_{ik}{C}_{kj}(t)+\sum _{k=1}^M{\int }_0^t{K}_{ik}(t-s)C_{kj}(s)ds.$$

$$⟨f,\mathcal{L}g{⟩}_{eq}=-⟨\mathcal{L}f,g{⟩}_{eq}f,g\in L^2(\mathcal{M},{\rho }_{eq})\cap \mathcal{D}(\mathcal{L}).$$

$$K_{ij}(t)=$$

$$=$$

$$e^{t\mathcal{L}}=\underset{q\to \infty }{\lim }{\left(1-\frac{t\mathcal{L}}{q}\right)}^{-q}.$$

$$e^{t\mathcal{Q}\mathcal{L}}\simeq \sum _{q=0}^n{g}_q(t){\Phi }_q\left(\mathcal{Q}\mathcal{L}\right).$$

$$K_{ij}(t)=$$

$$\simeq$$

$$=$$

$$M_{ijq}=\sum _{k=1}^M{G}_{jk}^{-1}⟨u_k(0),\mathcal{L}{\Phi }_q(\mathcal{Q}\mathcal{L})\mathcal{Q}\mathcal{L}{u}_i(0){⟩}_{\rho }.$$

$$K(t)\simeq \sum _{q=0}^n{g}_q(t)M_q,\text{where}{M}_q=\frac{⟨u(0),\mathcal{L}{\Phi }_q(\mathcal{Q}\mathcal{L})\mathcal{Q}\mathcal{L}u(0){⟩}_{\rho }}{⟨u(0),u(0){⟩}_{\rho }}.$$

$$\mathcal{P}f=\frac{⟨f,u(0){⟩}_{\rho }}{⟨u(0),u(0){⟩}_{\rho }}u(0).$$

$${\mu }_i=\frac{⟨\mathcal{L}(\mathcal{Q}\mathcal{L}{)}^{i-1}u(0),u(0){⟩}_{\rho }}{⟨u(0),u(0){⟩}_{\rho }},{\gamma }_i=\frac{⟨{\mathcal{L}}^{i}u(0),u(0){⟩}_{\rho }}{⟨u(0),u(0){⟩}_{\rho }}.$$

$${\mu }_{i}u(0)=\mathcal{P}\mathcal{L}(\mathcal{Q}\mathcal{L}{)}^{i-1}u(0).$$

$${\mu }_1={\gamma }_1,{\mu }_2={\gamma }_2-{\mu }_1{\gamma }_1,\cdots ,{\mu }_n$$

Full text

††journal: Arxiv

1 Introduction

The Mori-Zwanzig (MZ) formulation is a technique developed in statistical mechanics to formally integrate out phase variables in nonlinear dynamical systems by means of a projection operator. One of the main features of such formulation is that it allows us to systematically derive formally exact generalized Lagevin equations (GLEs) for quantities of interest (macroscopic observables), based on microscopic equations of motion. Such GLEs can be found in a variety of applications, including particle dynamics [52, 66, 37, 38, 35], partial differential equations (PDEs) [11, 54, 8, 10, 59], fluid dynamics [46, 47], and solid-state physics [65, 36, 42]. As an example, consider the Brownian motion of a colloidal particle subject to collision interactions with a large number of fluid particles. By using the MZ formulation it is possible to derive a low-dimensional system of equations characterizing the dynamics (position and momentum) of the colloidal particle alone [27, 35].

Computing the solution to the MZ equation is a challenging task. One of the main difficulties is the approximation of the memory integral (convolution term), and the fluctuation term (noise term). These terms encode the interaction between the so-called orthogonal dynamics and the dynamics of the quantity of interest. The orthogonal dynamics is essentially a high-dimensional nonlinear flow that satisfies a hard-to-solve integro-differential equation. Such flow has, in general, the same order of magnitude and dynamical properties as the quantity of interest, i.e., there is no general scale separation between the so-called resolved and the unresolved variables [14, 53]. As a consequence, approximating the MZ memory integral and the fluctuation term in these cases is often a daunting task, because of the strong coupling between the orthogonal dynamics and the macroscopic observables. Over the years, many techniques have been proposed to address this problem. These techniques can be grouped in two categories: i) data-driven methods; ii) methods based on first-principles. Data-driven methods aim at recovering the MZ memory integral/fluctuation term based on data, usually in the form of sample trajectories of the full system. Typical examples are the NARMAX technique developed by Lu et al. [39], the rational function approximation recently proposed by Lei et al. [35] (see also [15]), and the conditional expectation technique developed by Brennan and Venturi [7]. On the other hand, methods based on first principles aim at approximating the MZ memory integral and fluctuation term based on the structure of the nonlinear system (microscopic equations of motion), without using any simulation data. The first effective method developed within this class is the continued fraction expansion of Mori [43], which can be conveniently formulated in terms of recurrence relations [34, 20]. Other methods based on first-principles include perturbation methods [63, 60], mode coupling techniques, [51, 22], optimal prediction methods [11, 14, 53, 8], and various series expansion [55, 47, 46, 67]. First-principle calculation methods can effectively capture non-Markovian memory effects, e.g., in coarse-grained particle simulations [66, 24]. However, they are often quite involved and they do not generalize well to systems with no scale separation [21]. At the same time, data-driven methods can yield accurate results, but they often require a large number of sample trajectories to faithfully capture memory effects [7, 15, 35, 37, 38]

In this paper, we present a new method to compute the MZ memory integral and the flucuation term from first principles in nonlinear systems with local polynomial interactions. To this end, we build upon the Faber operator series expansion we recently developed in [68], and a new combinatorial algorithm that allows us to compute the MZ memory kernel by using only the structure of the microscopic equations of motion. We also develop a new data-driven stochastic process representation method based on the MZ memory kernel and Karhunen-Loève (KL) series expansions, which allows us to build simple models of the MZ fluctuation term in systems with invariant measures, e.g., Hamiltonian systems or more general systems [6, 19].

This paper is organized as follows. In Section 2 we briefly review the MZ formulation for nonlinear dynamical systems evolving from random initial states. In Section 2.1 we specialize such formulation by introducing Mori’s projection operators. In Section 2.2 and Section 2.3 we develop series expansion of the MZ memory kernel based on the Faber operator series we recently proposed in [68]. Section 2.4 is devoted to the description of an exact combinatorial algorithm to compute the recurrence coefficients of the MZ memory kernel expansion. In Section 3, we develop a new stochastic process representation method to compute the MZ fluctuation term for systems in statistical equilibrium. In Section 4 we demonstrate the accuracy of the MZ memory calculation and the reduced-order stochastic modeling technique in applications to nonlinear random wave propagation described by Hamiltonian partial differential equations. The main findings of the paper are summarized in Section 5. We also include a brief Appendix where prove convergence of KL expansions in representing auto-correlation functions of polynomial observables.

2 The Mori-Zwanzig formulation

Consider the following nonlinear dynamical system evolving on a smooth manifold $\mathcal{M}\subseteq\mathbb{R}^{N}$

[TABLE]

where $\bm{x}_{0}\in\mathcal{M}$ is a random initial state with probability density function $\rho_{0}(\bm{x})$. The dynamics of any vector-valued phase space function

[TABLE]

can be expressed in terms of a semi-group of linear operators acting on $\bm{u}(\bm{x}_{0})$, i.e.,

[TABLE]

In this equation, $\bm{x}(t,\bm{x}_{0})$ represents the flow [64, 25] generated by the system (1), while $e^{t\mathcal{L}}$ is the composition (Koopman) operator of the system [31, 17]. We are interested in deriving the exact evolution equation for the phase space function $\bm{u}(t)=\bm{u}(\bm{x}(t,\bm{x}_{0}))$. To this end, we employ the Mori-Zwanzig formulation [69, 11, 67]. The first step is to introduce an orthogonal projection operator $\mathcal{P}$, and the complementary projection $\mathcal{Q}=\mathcal{I}-\mathcal{P}$, where $\mathcal{I}$ is the identity operator. The mathematical properties of such projections are discussed in detail in [67, 17]. By differentiating the well-known Dyson’s identity

[TABLE]

with respect to time, we obtain the following evolution equation for the Koopman operator $e^{t\mathcal{L}}$

[TABLE]

Applying this equation to any phase space function $\bm{u}(0)=\bm{u}(\bm{x}_{0})$ yields the Mori-Zwanzig (MZ) equation

[TABLE]

The three terms at the right hand side are called, respectively, streaming term, fluctuation (or noise) term, and memory term. It is often more convenient (and tractable) to compute the evolution of the observable $\bm{u}(t)$ within a closed linear space, e.g., the image of the projection operator $\mathcal{P}$. To this end, we apply such projection to both sides of equation (6). This yields the following exact evolution equation111Note that the projected fluctuation term $\mathcal{P}e^{t\mathcal{QL}}\mathcal{QL}\bm{u}(0)$ is identically zero since $\mathcal{P}\mathcal{Q}=0$.

[TABLE]

Depending on the choice of the projection operator, the MZ equation (7) can yield evolution equations for different quantities. For example, if we use Chorin’s projection [11, 12, 67, 59], then (7) is an evolution equation for the conditional mean of $\bm{u}(t)$. Similarly, if we use Mori’s projection [68, 52], then (7) is an evolution equation for the temporal auto-correlation function of $\bm{u}(t)$.

2.1 Mori’s projection operator

Suppose that the phase space function (2) belongs to the weighted Hilbert space $H=L^{2}(\mathcal{M},\rho)$, where $\rho$ is a positive weight function in $\mathcal{M}$. For instance, $\rho$ can be the probability density function of the random initial state $\bm{x}_{0}$ (i.e., $\rho_{0}$, see Eq. (1)), or the equilibrium distribution of the system $\rho_{eq}$ (assuming it exists). Let

[TABLE]

be the inner product in $H$. For any closed linear subspace $V\subset H$ the Mori projection operator $\mathcal{P}$ is defined to be the orthogonal projection onto $V$, relative to the inner product (8). If $V$ is finite-dimensional with dimension $M$, then $\mathcal{P}$ can be effectively constructed if we are given $M$ linearly independent functions $u_{i}(0)=u_{i}(\bm{x})\in V$ ($i=1,...,M$). Clearly, if $\{u_{1}(0),\dots,u_{M}(0)\}$ are linearly independent then $V=\text{span}\{u_{1}(0),\dots,u_{M}(0)\}$. To construct Mori’s projection, we first compute the positive-definite Gram matrix $G_{ij}=\langle u_{i}(0),u_{j}(0)\rangle_{\rho}$, i.e.,

[TABLE]

With $G_{ij}$ available, we define

[TABLE]

In classical statistical dynamics of Hamiltonian systems, a common choice for the density $\rho$ is the Boltzmann-Gibbs distribution

[TABLE]

where $\mathcal{H}(\bm{x})=\mathcal{H}(\bm{q},\bm{p})$ is the Hamiltonian of the system, $\bm{x}=(\bm{q},\bm{p})$ are generalized coordinates/momenta, and $Z$ is the partition function. For other systems, $\rho$ can be, e.g., the probability density function of the random initial state (see Eq. (1)). Next, suppose that each observable $u_{i}(\bm{x})$ ($i=1,\dots,M$) belongs to the linear space $\mathcal{P}H\cap\mathcal{D}(\mathcal{L})$, where $\mathcal{P}H=V$ and $\mathcal{D}(\mathcal{L})$ denotes the domain of the Liouville operator $\mathcal{L}$ defined in (3). The MZ equation (6), with $\mathcal{P}$ defined in (10), reduces to

[TABLE]

where222Note that the $i$th component of the system (12) can be explicitly written as

$\frac{du_{i}(t)}{dt}=\sum_{j=1}^{M}\Omega_{ij}u_{j}(t)+\sum_{j=1}^{M}\int_{0}^{t}K_{ij}(t-s)u_{j}(s)ds+f_{i}(t).$

(13)

[TABLE]

Equation (12) is often referred to as generalized Langevin equation (GLE) in classical statistical physics and other disciplines [52]. By applying Mori’s projection to (12) we obtain the following linear (and closed) evolution equation for the projected phase space function

[TABLE]

Acting with the inner product $\displaystyle\langle u_{j}(0),\cdot\rangle_{\rho}$ on both sides of equation (15), yields the following exact equation for the temporal auto-correlation matrix $C_{ij}(t)=\langle u_{j}(0),u_{i}(t)\rangle_{\rho}$

[TABLE]

Suppose that the system (1) is Hamiltonian, and that the random initial state $\bm{x}_{0}$ is distributed according to the Boltzmann-Gibbs distribution (11), i.e., $\rho_{0}=\rho_{eq}$. In these assumptions, the Liouville operator $\mathcal{L}$ is skew-adjoint relative to the inner product (8), i.e., we have

[TABLE]

This allows us to simplify the expression of the memory kernel (14c) as

[TABLE]

where $f_{k}(t)$ is the $k$-th component of the fluctuation term (14d). The identity (18) is known as Kubo’s second fluctuation-dissipation theorem [33]. We emphasize there are several advantages in using Mori’s projection (10) over other projection operators, e.g., Chorin’s projection [13]. For example, both MZ equations (12) and (15) are linear and closed, which allows us perform rigorous convergence analysis [68, 67]. Secondly, the streaming matrix (14b) and the memory kernel (14c) are exactly the same for both the projected and the unprojected equations ,i.e., (12) and (15)). Thirdly, we have that the second-fluctuation dissipation theorem (18) holds true, which allows us to express the MZ memory kernel in a relatively simple form in terms of averages of random forces.

2.2 Series expansion of the MZ memory kernel

To compute the solution of the Mori-Zwanzig equation (15) we need to evaluate the memory kernel (14c). This is often a daunting task due to the presence of terms such as $e^{t\mathcal{Q}\mathcal{L}}u_{i}(0)$, i.e., terms involving operator exponentials. A straightforward method to evaluate such terms is to expand the orthogonal dynamics propagator $e^{t\mathcal{Q}\mathcal{L}}$ in a truncated operator polynomial series as333 The series expansion (19) needs to be handled with care if $\mathcal{L}$ is an unbounded operator, e.g., the generator of the Koopman semigroup (3) (see [28], p. 481). In this case, $e^{t\mathcal{L}}$ should be properly defined as

$e^{t\mathcal{L}}=\lim_{q\rightarrow\infty}\left(1-\frac{t\mathcal{L}}{q}\right)^{-q}.$

In fact, $\left(1-t\mathcal{L}/q\right)^{-1}$ is the resolvent of $\mathcal{L}$ (modulus a constant factor), which can be rigorously defined for both bounded and unbounded operators.

[TABLE]

The temporal modes $g_{q}(t)$ can be explicitly computed for any specific choice of the polynomial set $\{\Phi_{0},\dots,\Phi_{n}\}$. For example, if $\Phi_{q}$ ($q=0,\dots,n$) are Faber polynomials [68, 44], then $g_{q}(t)$ are products of exponential functions and Bessel functions of the first kind. Similarly, if $\Phi_{q}\left(\mathcal{Q}\mathcal{L}\right)=(\mathcal{Q}\mathcal{L})^{q}$, then $g_{q}(t)=t^{q}/q!$ (see Table 1). A substitution of (19) into (14c) yields the following series expansion of the MZ memory kernel

[TABLE]

where

[TABLE]

Note that the coefficients $M_{ijq}$ are generalized operator cumulants in the sense of Hegerfeldt and Kubo [23, 32, 52]. Such Coefficients can be computed in a data-driven setting [35, 4], or based on first-principles as we describe in Section 2.3. We also emphasize that, in general, series expansions of the orthogonal dynamics propagator (19) with different polynomial bases $\Phi_{q}$ can yield different approximation errors in the MZ memory kernel (20) for the same polynomial order $n$ (see [68]).

Regarding convergence of (19), our recent analysis [68, 67] demonstrates that it can be rigorously established for linear systems and arbitrary integration times (see §5 in [68]). If the system is nonlinear, then the series expansion of the memory kernel (20) is granted to converge only on a time interval that depends on the system and on the observable $\bm{u}(t)$. In particular, Corollary 3.4.3 in [67] establishes short-time convergence of the MZ memory approximation (20) for a broad class of nonlinear systems of the form (1). This implies that such memory approximation might be accurate only for a short integration time that depends on the system and the observable.

2.3 Calculation of the MZ memory kernel from first principles

In this Section we develop a new algorithm to calculate the MZ memory kernel (14c) based on first-principles, i.e., based on the microscopic equations of motion of the system (1). The algorithm we propose is built upon the combinatorial approach originally proposed by Amati, Meyer and Schilling in [1]. To illustrate the main idea in a simple way, hereafter we study the case where the observable $\bm{u}(t)$ is one-dimensional, i.e., we have only one phase space function $u(t)=u(\bm{x}(t,\bm{x}_{0}))$. In this setting, the series expansion (20) reduces to

[TABLE]

Note that $K(t)$ depends only on the set of parameters $\{M_{0},\dots,M_{n}\}$, since the temporal modes $g_{q}(t)$ are explicitly available given the polynomial set $\{\Phi_{0},\dots,\Phi_{n}\}$ (see Table 1). We aim at determining $\{M_{0},\dots,M_{n}\}$ from first principles. For one-dimensional phase space functions, Mori’s projection (10) reduces to

[TABLE]

At this point, it is convenient to introduce the following notation

[TABLE]

Each coefficient $\mu_{i}$ represents the rescaling of $u(0)$ under the action of the operator $\mathcal{P}\mathcal{L}(\mathcal{Q}\mathcal{L})^{i-1}$, i.e. we have

[TABLE]

Clearly, if we are given $\{\mu_{1},\dots,\mu_{n+2}\}$ then we can easily compute $M_{q}$ in (22), and therefore the memory kernel $K(t)$ for any given polynomial function $\Phi_{q}$. For example, if $\Phi_{q}(\mathcal{Q}\mathcal{L})=(\mathcal{Q}\mathcal{L})^{q}$ then $M_{q}=\mu_{q+2}$ ($q=0,\dots,n$). On the other hand, if $\{\Phi_{0},\dots,\Phi_{n}\}$ are Faber polynomials [68], then we can write each $\Phi_{q}$ as a linear combination of monomials $(\mathcal{Q}\mathcal{L})^{j}$ $(j=0,\dots,q)$ and therefore represent $M_{q}$ as a linear combination of $\{\mu_{1},\dots,\mu_{q+2}\}$. Computing $\mu_{i}$ using the definition (24) involves taking operator powers and averaging, which may be computationally expensive. An alternative effective algorithm relies on the following recursive formula [15, 52, 5]

[TABLE]

In practice, (26) shifts the problem of computing $\{\mu_{1},\dots,\mu_{n}\}$ to the problem of evaluating the coefficients $\{\gamma_{1},\dots,\gamma_{n}\}$ defined in (24). This will be discussed extensively in the subsequent Section 2.4. If the Liouville operator $\mathcal{L}$ is skew-adjoint relative to the inner product (8), then all $\mu_{j}$ and $\gamma_{j}$ corresponding to odd indices are identically zero. This allows us to simplify the recursion (26) as

[TABLE]

As a consequence, the streaming term (14b) in the MZ equation vanishes identically since $\Omega=\mu_{1}=\gamma_{1}=0$. We recall that skew-adjoint Liouville operators arise naturally, e.g., in Hamiltonian dynamical systems at statistical equilibrium.

2.4 Systems with polynomial nonlinearities

In this Section, we address the problem of calculating the coefficients $\{\gamma_{1},\dots,\gamma_{n}\}$ defined in (24) and appearing in the recursion relation (26). With such coefficients available, we can compute $\{\mu_{1},\dots,\mu_{n}\}$ and therefore the MZ memory kernel (22). The calculation we propose is based on first principles, meaning that we do not rely on any assumption or model to evaluate the averages $\gamma_{i}=\langle\mathcal{L}^{i}u(0),u(0)\rangle_{\rho}/\langle u(0),u(0)\rangle_{\rho}$. Instead, we develop a combinatorial algorithm that allows us to track all terms in $\mathcal{L}^{i}u(0)$, hence representing $\gamma_{i}$ exactly as a superimposition of a finite, although possibly large, number of terms. The algorithm we develop is built upon the combinatorial algorithm recently proposed by Amati, Meyer and Schiling in [1]. To describe the algorithm, consider the nonlinear dynamical system (1) and assume that $\bm{F}(\bm{x})$ is a multivariate polynomial in the phase variables $\bm{x}$. A simple example of such system is the Kraichnan-Orszag three-mode problem [45, 62, 7]

[TABLE]

Other examples are the semi-discrete form of the Navier-Stokes equations, or the semi-discrete form of the nonlinear wave equation discussed in Section 4. The key observation to compute $\gamma_{j}$ for systems with polynomial nonlinearities is that the action of the operator power $\mathcal{L}^{i}$ on a polynomial observable $u(\bm{x})$ yields a polynomial function. For instance, consider $u(\bm{x})=x_{1}^{3}$, and the Liouville operator associated with the system (28)

[TABLE]

We have

[TABLE]

Clearly, the number of terms in $\mathcal{L}^{i}x_{1}^{3}$ can rapidly increase, if high-order powers of $\mathcal{L}$ are considered. For higher-dimensional systems with non-local interactions, i.e., for systems where each $F_{i}(\bm{x})$ ($i=1,\dots,N$) depends on all components of $\bm{x}$, this problem is serious, and requires multi-core computer-based combinatorics to systematically track all terms in the expansion of $\mathcal{L}^{i}x_{j}^{q}$. Let us introduce the following notation

[TABLE]

where $\{a_{\bm{b}_{i}}^{(n)}\}$ are polynomial coefficients, and $\{m^{(i)}_{k_{j}}\}$ are polynomial exponents. The set of indexes representing the relevant phase phase variables appearing in $\mathcal{L}^{n}x_{j}^{q}$, i.e., $\{k_{1},\dots,k_{r}\}$, is collected in the index set $K(n,j)=\{k_{1},\dots,k_{r}\}$, which depends on $n$ and $j$. For example, in (30)-(32) we have

[TABLE]

Of course, for low-dimensional dynamical systems, the simplest choice for the relevant variables would be the complete set of variables $\{x_{1},\cdots,x_{N}\}$. However, for high-dimensional systems with local interactions this choice could lead to unnecessary computations. In fact, it can be shown that the variables appearing in the polynomial $\mathcal{L}^{n}x_{j}^{q}$ are usually a (possibly small) subset of the phase variables if the system has local interactions. The vector $\bm{b}_{i}=[m^{(i)}_{k_{1}},\cdots,m^{(i)}_{k_{r}}]$, collects the exponents of the $i$-th monomial appearing in the expansion (33). Similarly, $a^{(n)}_{\bm{b}_{i}}$ is the coefficient multiplying $i$-th monomial in (33). For example, in (30) and (31) we have, respectively,

[TABLE]

[TABLE]

At this point, it is convenient to define the set of polynomial exponents $B^{(n)}=\{\bm{b}_{1},\bm{b}_{2},\cdots\}$, the set polynomial coefficients $A^{(n)}=\{a^{(2)}_{\bm{b}_{1}},a^{(2)}_{\bm{b}_{2}},\cdots\}$, and the combined index set

[TABLE]

Clearly, $\mathcal{I}^{(n)}$ identifies uniquely the polynomial (33), i.e., there is a one-to-one correspondence between $\mathcal{I}^{(n)}$ and $\mathcal{L}^{n}x_{j}^{q}$. For example, in the case of (30)-(32) we have

[TABLE]

If we apply $\mathcal{L}$ to (33) we obtain

[TABLE]

Clearly, if we can compute the mapping $\mathcal{I}^{(n)}\xrightarrow{\mathcal{L}}\mathcal{I}^{(n+1)}$, induced by the action of the Liouville operator $\mathcal{L}$ to the polynomial (33) (represented by $\mathcal{I}^{(n)}$), then we can compute the exact series expansion of $\mathcal{L}^{n}x_{j}^{q}$ for arbitrary $n$. With such expansion available, we can immediately determine the coefficients $\gamma_{j}$ in (24) by averaging over the probability density $\rho$ as

[TABLE]

If the operator $\mathcal{L}$ is skew-adjoint in $L^{2}(\mathcal{M},\rho)$, i.e., if $\langle\mathcal{L}f,g\rangle_{\rho}=-\langle f,\mathcal{L}g\rangle_{\rho}$, then we have

[TABLE]

As pointed out by Maiocchi et al. in [40], the value of the first few coefficients $\{\gamma_{n}\}$ in (40) or (41) can be used to identify non-exponential relaxation patterns to equilibrium.

Remark

To enhance numerical stability when computing the Faber expansion of the MZ memory kernel we scale the Liouville operator $\mathcal{L}$ by a parameter $\delta<1$ (see [68, 26]), i.e., we compute the Faber operator polynomial series relative to $\delta\mathcal{L}$. Correspondingly, the generalized Langevin equation (15) is solved on a time scale $t/\delta$. In this setting, the coefficients (40) are also calculated relative to the rescaled Liouville operator $\delta\mathcal{L}$.

Remark

Computing $\gamma_{j}$ for linear systems reduces to a classical numerical linear algebra problem, i.e., computing the Rayleigh quotient of a matrix power. To show this, consider the $N$-dimensional linear system $\dot{\bm{x}}=\bm{A}\bm{x}$, evolving from the random initial state $\bm{x}_{0}\sim\rho_{0}$ ($\bm{x}_{0}$ column vector). Suppose we are interested in the first component of the system, i.e., set the observable as $u(t)=x_{1}(t,\bm{x}_{0})$. Define the linear subspace $V=\text{span}\{x_{01},x_{02},\cdots,x_{0N}\}\subset L^{2}(\mathcal{M},\rho_{0})$. Clearly $u(t)\in V$ for all $t\geq 0$ [67, 68]. This allows us to calculate $\gamma_{j}$ as

[TABLE]

where $\bm{e}_{1}=[1,0,\dots,0]^{T}$.

2.5 Mapping the index set $\mathcal{I}^{(n)}$

Here we describe the algorithm that allows us to compute the polynomial $\mathcal{L}^{n+1}x_{k}^{q}$ given the polynomial $\mathcal{L}^{n}x_{k}^{q}$ and the Liouville operator $\mathcal{L}$, i.e., the mapping defined by equation (39). This is equivalent to develop a set of algebraic rules to transform the combined index set $\mathcal{I}^{(n)}$ defined in (35) into $\mathcal{I}^{(n+1)}$, for arbitrary $n$. Once such rules are known, we can apply them recursively to compute the polynomial sequence

[TABLE]

to any desired order. In this way, we can determine $\gamma_{n}$ through (40) (or (41)), $\mu_{n}$ through (26) (or (27)), and therefore the MZ memory kernel (22). Before formulating the algorithm in full generality, it is useful to examine how it operates in a concrete example. To this end, consider again the Kraichnan-Orszag system (28), and the transformation between the polynomials (31) and (32) defined by the action of the Liouville operator (29). We are interested in formulating such transformation in terms of a set of algebraic operations mapping the index set $\mathcal{I}^{(2)}$ into $\mathcal{I}^{(3)}$ (Eqs. (37)-(38)). We begin by decomposing the three-dimensional Liouville operator (29) as

[TABLE]

The action of $\mathcal{L}_{i}$ on any monomial generates a polynomial with $S_{i}$ terms. In the present example, we have $S_{1}=S_{2}=1$ and $S_{3}=2$. Let us now consider the first monomial in (31), i.e., $9x_{1}^{3}x_{3}^{2}$. Such monomial is represented by the first element of $A^{(2)}$ and $B^{(2)}$ in (38). The corresponding combined set is $\{9,[3,0,2]\}$. At this point, we apply the operators $\mathcal{L}_{1}$, $\mathcal{L}_{2}$ and $\mathcal{L}_{3}$ to the polynomial $\{9,[3,0,2]\}$. This yields

[TABLE]

The transformation associated with $\mathcal{L}_{3}$ generates the sum of two monomials, namely $18x_{1}^{3}x_{2}^{2}x_{3}-18x_{1}^{5}x_{3}$, which we represent as a union between two index sets. Such union, here denoted as $\biguplus$, is an ordered union that pushes to the left polynomial coefficients and to the right polynomial exponents. Proceeding in a similar manner for all other monomials in (31) and taking ordered unions of all sets, yields the desired mapping $\mathcal{I}^{(2)}\rightarrow\mathcal{I}^{(3)}$. Let us now examine the action of a more general Liouville operator

[TABLE]

on the monomial $ax_{1}^{m_{1}}\cdots x_{N}^{m_{N}}$ represented by the index set $\{a,[m_{1},\dots,m_{N}]\}$. We have

[TABLE]

This defines two linear transformations: a scaling transformation in the space of coefficients, and an addition in the space of exponents

[TABLE]

In a vector notation, upon definition of $\bm{b}=[m_{1},\dots,m_{N}]$, $\bm{\theta}_{j}=[c_{1},\dots,c_{j}-1,\dots,c_{N}]$ and $\alpha_{j}=zm_{j}$, we can write (49) in compact form as

[TABLE]

Let us know consider the general case where the Liouville operator is defined as

[TABLE]

and $F_{k}(\bm{x})$ is a polynomial involving $S_{k}$ monomials in either all variables $\{x_{1},\dots,x_{N}\}$ or a subset of them. The action of $\mathcal{L}$ on each monomial in (39) can be written as

[TABLE]

where $K(n,j)=\{k_{1},\dots,k_{r}\}$ is the set of relevant variables at iteration $n$. The polynomial (52) involves $S_{k_{1}}+\cdots+S_{k_{r}}$ terms, each one of which can be explicitly constructed by applying the linear transformation rules (50). In summary, we have

[TABLE]

where $\#B^{(n)}$ denotes the number of elements in $B^{(n)}$. Note that both $\alpha_{s}^{s}$ and $\bm{\theta}_{s}^{q}$ depend on $q\in K(n,j)$ (index set of relevant variables).

Remark

The recursive algorithm summarized by formula (53) is a modified version of the algorithm originally proposed by Amati, Meyer and Schiling in [1]. The key idea is the same, i.e., to compute the expansion coefficients $\gamma_{n}$ in (40) using polynomial differentiation. However, there are a few differences between our algorithm and the algorithm proposed in [1] which we emphasize hereafter. In [1], the index set $B^{(n)}$ is pre-computed using the so-called spreading operators. Essentially, for each $n$, the iterative scheme generates a new set of polynomial coefficients $A^{(n)}$, which is subsequently matched with the corresponding indexes in $B^{(n)}$. In our algorithm, the sets $B^{(n)}$ and $A^{(n)}$ are computed on-the-fly at each step of the recursion. By doing so, we avoid calculating the spreading operators. This, in turn, allows us to avoid using numerical tensors to store index sets, since in our formulation there is no matching procedure between the polynomial exponents and the polynomial coefficients. Another difference between the two algorithms is that we utilized a rescaled Liouville operator $\delta\mathcal{L}$ ($\delta\in\mathbb{R}$) to enhance numerical stability when computing the operator polynomials $\Phi_{q}(\mathcal{Q}\mathcal{L})$ in (21). The algorithm in [1], on the other hand, is based on a Taylor series expansion of the operator exponential $e^{t\mathcal{L}}$, with unscaled Liouville operator444 In our recent work [68] (§3.1) we proved that a Taylor series of the orthogonal dynamical propagation $e^{t\mathcal{Q}\mathcal{L}}$ yields an expansion of the MZ memory integral that resembles the classical Dyson series in scattering theory..

2.6 An example: the Fermi-Pasta-Ulam (FPU) model

Consider a one-dimensional chain of $N$ anharmonic oscillators with Hamiltonian

[TABLE]

In (54) $\{q_{j},p_{j}\}$ are, respectively, the generalized coordinate and momentum of the $j$-th oscillator, while $V(q_{i}-q_{i-1})$ is the potential energy between two adjacent oscillators. Suppose that the oscillator chain is closed (periodic), i.e., that $q_{0}=q_{N}$ and $p_{0}=p_{N}$. Define the distance between two oscillators as $r_{j}=q_{j}-q_{j-1}$. This allows us to write the Hamilton’s equations of motion as

[TABLE]

The Liouville operator corresponding to this system is

[TABLE]

Setting $V(x)=\alpha x^{2}/2+\beta x^{4}/4$ yields the well-known Fermi-Pasta-Ulam $\beta$-model [42], which we study hereafter. To this end, suppose we are interested in the distance between the oscillators $j$ and $j-1$, i.e., in the polynomial observable $u(\bm{p},\bm{r})=r_{j}$. The action of $\mathcal{L}^{n}$ on $r_{j}$ can be explicitly written as

[TABLE]

where $\{k_{1},\dots,k_{u}\}$ and $\{l_{1},\dots,l_{v}\}$ are the relevant degrees of freedom for the polynomials of $\bm{r}$ and $\bm{p}$, respectively, at iteration $n$. We can explicitly compute the sets of such relevant degrees of freedom as

[TABLE]

The action of the Liouville operator on each monomial appearing in (55) can be written as

[TABLE]

where

[TABLE]

By computing the action of $\mathcal{L}_{r_{v}}$ and $\mathcal{L}_{p_{h}}$ on the monomial $r_{k_{1}}^{m^{(i)}_{k_{1}}}\cdots r_{k_{u}}^{m^{(i)}_{k_{u}}}p_{l_{1}}^{s^{(i)}_{l_{1}}}\cdots_{l_{v}}^{s^{(i)}_{l_{v}}}$ we obtain explicit linear maps of the form (50), involving the polynomial exponents

[TABLE]

and the polynomial coefficients $a_{\bm{b}_{i}}^{(n)}$. With such maps available, we can transform the combined index set $\mathcal{I}^{(n)}$ (representing $\mathcal{L}^{n}r_{j}$) to $\mathcal{I}^{(n+1)}$ (representing $\mathcal{L}^{n+1}r_{j}$) using (53). Specifically, we obtain

[TABLE]

where

[TABLE]

3 Modeling the MZ fluctuation term

In previous Sections we discussed an algorithm to approximate the memory kernel in the MZ equation (12) or (15) based on the microscopic equations of motion (first-principle calculation). In this Section we construct a statistical model for fluctuation term $\bm{f}(t)$ appearing in (12), which will allow us to compute statistical properties of the quantity of interest $\bm{u}(t)$ beyond two-point correlations. A possible way to build such model is to expand (14d) in a finite-dimensional series555Note that $\bm{f}(t)$ is a random process obtained by mapping the random initial state $\bm{u}(0)=\bm{u}(\bm{x}_{0})$ forward in time using the orthogonal dynamics propagator $e^{t\mathcal{Q}\mathcal{L}(\bm{x}_{0})}$. (see Eq. (19)) as

[TABLE]

and evaluate the coefficients $\Phi_{q}(\mathcal{Q}\mathcal{L})\mathcal{Q}\mathcal{L}\bm{u}(0)$ using the combinatorial approach discussed in Section 2.4. However, this may not be straightforward since $\Phi_{q}(\mathcal{Q}\mathcal{L})\mathcal{Q}\mathcal{L}\bm{u}(0)$ is a high-dimensional random field. An alternative approach is to ignore the mathematical structure of $\bm{f}(t)$, i.e., equation (14d) or the series expansions (60), and simply model $\bm{f}(t)$ as a stochastic process. In doing so, we need to make sure that the statistical properties of the reduced-order model, e.g., the equilibrium distribution, are consistent with the full model. Such consistency conditions carry over a certain number of constraints on $\bm{f}(t)$, which allow for its partial identification. As an example, consider the following MZ model recently proposed by Lei et al. in [35] to study the dynamics of a tagged particle in a large particle system

[TABLE]

It was shown in [35] that if $\bm{f}(t)$ is Gaussian white noise (in time) with auto-correlation function

[TABLE]

then the equilibrium distribution of the particle system has the Boltzmann-Gibbs form

[TABLE]

$V(\bm{q})$ being the inter-particle potential. However, modeling $\bm{f}(t)$ as a Gaussian process does not provide satisfactory statistics in MZ equations is we use Mori’s projection. In fact, equation (12) is linear and therefore the equilibrium distribution of $\bm{u}(t)$ (assuming it exists) under Gaussian noise $\bm{f}(t)$ will be necessarily Gaussian. In most applications, however, the equilibrium distribution of $\bm{u}(t)$ is strongly non-Gaussian. To overcome this difficulty Chu and Li [15] recently developed a multiplicative Gaussian noise model that generalizes (12) in the sense that it is not based on additive noise, and it allows for non-Gaussian responses.

In this Section we propose a different stochastic modeling approach for $\bm{f}(t)$ based on bi-orthogonal representations random processes [57, 61, 56, 3, 2]. To describe the method, we study the case where the observable $\bm{u}(t)$ is real valued (one-dimensional) and square integrable. This allows us to develop the theory in a clear and concise way. We also assume that the system is in statistical equilibrium, i.e., that there exists an equilibrium distribution $\rho_{eq}(\bm{x})$ (or more generally an invariant measure) for the phase variables $\bm{x}(t,\bm{x}_{0})$, and that $\bm{x}_{0}$ is sampled from such distribution. The MZ equation (12) for one-dimensional phase space functions $u(t)=u(\bm{x}(t,\bm{x}_{0}))$ reduces to

[TABLE]

where

[TABLE]

Since $u(t)$ is assumed to be a second-order random process in the time interval $[0,T]$, we can expand it in a truncated Karhunen-Loéve series

[TABLE]

where $\overline{u}$ denotes the mean of $u(t)$ relative to the equilibrium distribution666The mean of $u(t)=u(\bm{x}(t,\bm{x}_{0}))$ is necessarily time-independent at statistical equilibrium. In fact, at equilibrium we have that $\bm{x}_{0}\sim\rho_{eq}$ implies that $\bm{x}(t)\sim\rho_{eq}$ for all $t\geq 0$. A statistically stationary process however, may not be stationary in phase space. Indeed, $\bm{x}(t)$ evolves in time, eventually in a chaotic way as it happens for systems with strange attractors and invariant measures., $\{\xi_{1},\dots,\xi_{K}\}$ are uncorrelated random variables ($\langle\xi_{i}\xi_{j}\rangle_{eq}=\delta_{ij}$), and $\{\lambda_{k},e_{k}(t)\}$ ($k=1,\dots,K$) are, respectively, eigenvalues and eigenfunctions of the homogeneous Fredholm integral equation of the second kind

[TABLE]

We recall that for ergodic systems in statistical equilibrium the auto-correlation function $\langle u(t)u(s)\rangle_{eq}$ decays to zero as $|t-s|\rightarrow\infty$. Also, the integral operator at the left hand side of (67) is positive-definite and compact [2]. The orthogonal random variables $\xi_{k}$ and the temporal modes $e_{k}(t)$ are related to each other by the following dispersion relations [3, 57]

[TABLE]

Equation (68) suggests that if $u(t)$ is a Gaussian random process (e.g., a Wiener process) then $\{\xi_{1},\dots,\xi_{K}\}$ are necessarily independent Gaussian random variables. On the other hand, if $u(t)$ is non-Gaussian then the joint distribution of $\{\xi_{1},\dots,\xi_{K}\}$ is unknown, although it can be in principle computed by using the transformation $u(t)\rightarrow\xi_{k}$ ($k=1,..,K$) defined in (68), given $\lambda_{k}$ and $e_{k}(t)$.

An alternative approach to identify the PDF of $\{\xi_{1},\dots,\xi_{K}\}$ relies on sampling. In particular, as recently shown by Phoon et al. [49, 50], it is possible to develop effective sampling algorithms for the KL expansion (66). Such algorithms allow to sample the uncorrelated variables $\{\xi_{1},\dots,\xi_{K}\}$ in a way that makes the PDF of $u(t)$ consistent with the equilibrium distribution, which can be calculated by mapping $\bm{x}_{0}\sim\rho_{eq}(\bm{x}_{0})$ to $u(\bm{x}_{0})$. At this point, we have available a consistent bi-orthogonal representation of the random process $u(t)$ defined by the series expansion (66). It is straightforward to see that such representation yields a corresponding series expansion of the fluctuation term $f(t)$ in (64). In fact we have the following

Proposition 1**.**

For any bi-orthogonal series expansion (66) of the solution to the MZ-equation (64), there exists a unique series expansion of the fluctuation term $f(t)$ of the form

[TABLE]

Proof.

It is sufficient to prove the theorem for zero-mean processes. To this end, we set $\overline{u}=0$ and $\overline{f}=0$ in (66) and (69). A substitution of (66) into (64) yields, for all $t\in[0,T]$

[TABLE]

Define,

[TABLE]

This equation does not allow us to compute $h_{k}$ explicitly quite yet. In fact, the MZ memory kernel $K(t-s)$ depends on $f(t)$ (see Eq. (65)). However, a substitution of (69) (with $\overline{f}=0$) into the analytical expression of $K(t)$ yields

[TABLE]

To evaluate $\langle\xi_{i},\mathcal{L}\xi_{j}\rangle_{eq}$ we need to express $\{\xi_{1},\dots,\xi_{K}\}$ as a function of $\bm{x}_{0}$ (recall that $\mathcal{L}$ operates on functions of $\bm{x}_{0}$, see Eq. (3)), and then integrate over $\rho_{eq}(\bm{x}_{0})$. This is easily achieved by using the dispersion relation (68). Specifically, we have

[TABLE]

At this point, we substitute (72) into (71) to obtain

[TABLE]

Given $\{e_{1}(t),\dots,e_{K}(t)\}$, $\Omega$ and $v_{ij}$, this equation can be solved uniquely for $\{h_{1}(t),\dots,h_{K}(t)\}$ by using Laplace transforms. Note that $\{h_{1}(t),\dots,h_{K}(t)\}$ are not necessarily orthogonal in $L^{2}([0,T])$.

∎

Remark

If the dynamical system (1) is Hamiltonian then the MZ steaming term vanishes, and the MZ memory kernel can be written in terms the fluctuation term as (see Eq. (18))

[TABLE]

A substitution of this expression into (64) yields, after projection onto $\xi_{k}$

[TABLE]

This equation establishes a one-to-one correspondence between the temporal modes of the KL expansion (66) and the temporal modes of the fluctuation term (70). In particular, given $\{e_{1}(t),\dots,e_{K}(t)\}$, we can determine $\{h_{1}(t),\dots,h_{K}(t)\}$ directly by using Laplace transforms, without building the MZ memory kernel (72).

3.1 Building MZ-KL stochastic models from first principles

Proposition 1 establishes a one-to-one correspondence between the noise process in the MZ equation (64) and the biorthogonal series expansion of the solution. This new paradigm allows us to build stochastic models for the observable $u(t)$ at statistical equilibrium from first principles. To this end,

Compute the solution to the MZ equation for the temporal correlation function of $u(t)$ (see Eq. (16))

[TABLE]

The memory kernel $K(t-s)$ can be expanded as in (22), and computed from first-principles using the combinatorial approach we discussed in Section 2.4. 2. 2.

Build the Karhunen-Loève expansion (66) by spectrally decomposing the correlation function $C(t)=\langle u(0)u(t)\rangle_{eq}$ obtained at point 1. Recall that at statistical equilibrium we have $C(t-s)=\langle u(0)u(t-s)\rangle_{eq}=\langle u(s)u(t)\rangle_{eq}$. This yields eigenvalues $\{\lambda_{j}\}$ and the eigenfunctions $e_{j}(t)$. The uncorrelated random variables $\xi_{k}$ appearing in (66) can be sampled consistently with the equilibrium distribution $\rho_{eq}$ by using, e.g., the iterative algorithm recently proposed by Phoon et al. [50, 49]. 3. 3.

With $\{\xi_{1},\dots,\xi_{K}\}$, $\{e_{1}(t),\dots,e_{K}(t)\}$ and $\{\lambda_{1},\dots,\lambda_{K}\}$ available, we can uniquely identify the noise process $f(t)$ in the MZ equation (64). To this end, we simply use Proposition 1, with the temporal modes $h_{k}(t)$ obtained by solving equation (74) or (76) with the Laplace transform. 4. 4.

With $K(t)$ computed from first principles, and $f(t)$ modeled based on the auto-correlation function $C(t)$, we can generate samples of the observable $u(t)$ by solving equation (64).

Remark

We emphasize that the correlation function $C(t)$ can be also computed directly from data, e.g., by using a Monte-Carlo or a quasi Monte Carlo method [16]. With $C(t)$ available it is possible to determine the fluctuation term $f(t)$ with equation (69) and the MZ memory kernel $K(t)$ using equation (72).

The results of this Section can be generalized to vector-valued phase space functions $\bm{u}(t)$ at statistical equilibrium. The starting point is the KL expansion for multi-correlated stochastic processes we recently proposed in [9]. Such expansion is constructed based on cross-correlation information777At statistical equilibrium the cross correlation functions are invariant under temporal shifts. This means that $\langle u_{i}(s),u_{j}(t)\rangle_{eq}=\langle u_{i}(0),u_{j}(t-s)\rangle_{eq}$ for all $t\geq s$. Hence, the solution to the projected MZ equation (16) is sufficient to compute the KL expansion of the multi-correlated process $\bm{u}(t)$, e.g., using the series expansion method proposed in [9]., and can be made consistent with the equilibrium distribution of $\bm{u}(t)$, e.g., by using the sampling strategy proposed in [50, 49]. The correspondence between the KL expansions of $\bm{u}(t)$ and the vector-valued fluctuation term $\bm{f}(t)$ can be established by following the same arguments we used in the proof of Proposition 1.

4 Applications to nonlinear systems with local interactions

In this Section, we demonstrate the accuracy of the MZ memory calculation method and the reduced-order stochastic modeling technique we discussed in Section 2.3 and Section 3, respectively. To this end, we study nonlinear random wave propagation described by Hamiltonian partial differential equations (PDEs). To derive such PDEs consider the nonlinear functional

[TABLE]

where $u=u(x,t)$ represents the wave displacement, $p=p(x,t)$ is the canonical momentum (field variable conjugate to $u(x,t)$), $u_{x}=\partial u/\partial x$, and $G(p,u_{x},u)$ is the nonlinear interaction term. By taking functional derivatives of (78) with respect to $p$ and $u$ (see, e.g., [58]) we obtain the Hamilton’s equations of motion

[TABLE]

The corresponding nonlinear wave equation is

[TABLE]

This equation has been studied extensively in mathematical physics [30, 18, 29, 48], in particular in general relativity, statistical mechanics, and in the theory of viscoelastic fluids. In Figure 1 and Figure 2, we plot a few sample numerical solutions to (80) corresponding to different initial conditions and different nonlinear interaction term $G(p,u_{x},u)$. These solutions are computed by an accurate Fourier spectral method with $N=512$ modes (periodic boundary conditions in $x\in[0,2\pi]$). Throughout this Section, we assume that the initial state $\{u(x,0),p(x,0)\}$ is random and distributed according to the functional Boltzmann-Gibbs equilibrium distribution888The partition function $Z(\alpha,\gamma)$ is defined as a functional integral over $u(x)$ and $p(x)$ (see, e.g., [58]).

[TABLE]

We emphasize that $\rho_{eq}([p],[u])$ is invariant under the infinite-dimensional flow generated by (80) with periodic boundary conditions, since the Hamiltonian (78) is a constant of motion (conserved quantity) in this case.

4.1 Linear waves

Setting the interaction term $G(p,u_{x},u)$ in (78) and (80) equal to zero yields the well-known linear wave equation

[TABLE]

We discretize (82) in space using second-order finite differences on the (periodic) grid $x_{j}=2\pi j/N$ ($j=0,\dots,N$). This yields the following linear dynamical system

[TABLE]

where $u_{j}(t)=u(x_{j},t)$, $p_{j}(t)=p(x_{j},t)$, and $h=2\pi/N$ is the mesh size. The Hamilton’s function corresponding to the finite-difference scheme (83) is obtained by discretizing the integral (78), e.g., with the rectangle rule. This yields

[TABLE]

where we defined $\alpha_{1}=\alpha/h^{2}$. The corresponding finite-dimensional Gibbs distribution can be written as

[TABLE]

$Z_{1}(\alpha_{1},\gamma)$ being the partition function (normalization constant). Note that we absorbed the scaling factor $h$ in the parameter $\gamma>0$. It is straightforward to verify that the lattice Hamiltonian (84) is preserved if $u_{0}=u_{N}$ and $p_{0}=p_{N}$ (periodic boundary conditions). This implies that the PDF (85) is invariant under the flow generated by the linear ODE (83). Note that the lattice Hamiltonian (84) coincides with the Hamiltonian of a one-dimensional chain of harmonic oscillators with uniform mass $m=1$ and spring constants $k=\alpha_{1}$. We set $N=100$ and $\alpha=(2\pi/100)^{2}$ in equation (83). In this way, the system (83) is $200$-dimensional and the modeling parameter $\alpha_{1}$ in (84)-(85) is equal to $1$.

MZ memory kernel and auto-correlation functions

The Hamiltonian system (83) with periodic boundary conditions has many symmetries. In particular, the statistical properties of wave displacement $u(x,t)$ at any point $x_{j}$ are the same, if the initial state is distributed according to (85). In addition, the PDF of the wave momentum999Note that for linear waves the wave momentum $p(x,t)$ is equal to $\partial u(x,t)/\partial t$ (see Eq. (79)). $p(x_{j},t)$ and the wave displacement $r(x_{j},t)=u(x_{j+1},t)-u(x_{j},t)$ are both Gaussian (see Eq. (85)). Suppose we are interested in the temporal auto-correlation function of the wave momentum $p(x_{j},t)=p_{j}$, at an arbitrary location $x_{j}$, i.e.,

[TABLE]

where $\langle,\rangle_{eq}$ is an integral over the equilibrium distribution (85). Such correlation function admits the analytical expression (see [20])

[TABLE]

where $J_{0}$ is the zero-order Bessel function of the first kind. With $C_{p_{j}}(t)$ available, we can solve the MZ equation

[TABLE]

for the memory kernel $K(t)$ by using Laplace transforms. This yields the exact MZ kernel

[TABLE]

where $J_{1}$ is the first-order Bessel function of the first kind. In Figure 3, we compare the exact memory kernel (89) and the correlation function (88) with the results we obtained using the iterative algorithm discussed in Section 2.5. Note that the system (83) is linear. Therefore, we can use the formula (42) to compute the coefficients $\{\gamma_{1},\dots,\gamma_{n+2}\}$. With such coefficients available, we then compute $\{\mu_{1},\dots,\mu_{n+2}\}$ using the recurrence relation (27), and the MZ memory kernel (22).

In Figure (3) we demonstrate that the MZ-Faber expansion rapidly converges to the exact auto-correlation function (86) of the wave momentum as we increase the Faber expansion order $n$. This is not surprising since the linear wave equation is a well-known integrable system for which convergence of the MZ-Faber series can be rigorously established (§5 in [68]).

Reduced-order stochastic modeling

Suppose we are interested in building a consistent reduced-order stochastic model for the wave momentum $p(x_{j},t)=\partial u(x_{j},t)/\partial t$ at statistical equilibrium. To this end, we employ the spectral expansion technique we discussed in Section 3. The auto-correlation function of the process $p(t)=p(x_{j},t)$ (at any location $x_{j}$), i.e., (86), is obtained by solving the MZ equation (88) with the kernel computed using the combinatorial algorithm described in Section 2.5. Following the stochastic modeling paradigm we developed in Section 3, we expand $p(t)$ as

[TABLE]

where $(\lambda_{k},e_{k}(t))$ are eigenvalues and eigenfunctions of (86). By enforcing consistency of (90) with the equilibrium distribution (85) at each fixed time we obtain that the random variables $p(t_{j})$ are normally distributed with zero mean and variance $1/\gamma$, for all $t_{j}\in[0,10]$. In other words $p(t)$ is a centered, stationary Gaussian random process with correlation function (86). In Figure 4, we plot the auto-correlation functions

[TABLE]

we obtained with an MZ-Faber expansion of degree $n=6$. Convergence of KL expansions representing high-order correlation functions such as (91) is established in Appendix.

4.2 Nonlinear waves

Here we study the nonlinear wave equation (80) with interaction term $G(p,u_{x},u)=\beta u_{x}^{4}/4$, i.e.,

[TABLE]

In Figure 1 and Figure 2 we plot sample solutions of (92) corresponding to different initial conditions. It is clearly seen that the nonlinearity $u_{x}^{2}u_{xx}$ breaks the periodicity of traveling wave. This effect is more pronounced if the initial condition is rougher in $x$, as $u^{2}_{x}$ and $u_{xx}$ are larger in this case, thereby increasing magnitude of the nonlinear term in (92). As before, we discretize (92) and the Hamiltonian (78) with finite differences on a periodic spatial grid ($N$ points in $[0,2\pi]$). This yields

[TABLE]

where $u_{j}(t)=u(x_{j},t)$ and $p_{j}(t)=\partial u(x_{j},t)/\partial t$ represent the wave amplitude and momentum at location $x_{j}=hj$ ($j=0,\dots,N$, $h=2\pi/N$), $\alpha_{1}=\alpha/h^{2}$ and $\beta_{1}=\beta/h^{4}$. The discretized equilibrium distribution (81) then becomes

[TABLE]

As before, we absorbed the factor $h$ into the parameter $\gamma$. Note that the lattice Hamiltonian (93) coincides with the Hamiltonian of the Fermi-Pasta-Ulam $\beta$-model (54), with $m_{j}=1$. We emphasize that if a different scheme is used to discretize the wave equation (92), then the lattice Hamiltonian (93) may not be a conserved quantity.

MZ memory term and auto-correlation functions

We choose the wave momentum $p_{j}(t)$ and the wave displacement $r_{j}(t)=u_{j+1}(t)-u_{j}(t)$ as quantities of interest. Moreover, we set $N=100$ and $\alpha=(2\pi/100)^{2}$. To study the effects of the nonlinear interaction term, we consider different values of $\beta=\beta_{1}\alpha^{2}$, with $\beta_{1}$ ranging from $0.01$ to $1$. This corresponds to the FPU models with mild and strong nonlinearities, respectively. Based on the structure of the Hamiltonian (93) and the equilibrium distribution (94), we expect that the dynamics of $p_{j}(t)$ and $r_{j}(t)$ will be different for different parameters $\beta$. To calculate the temporal auto-correlation function of $p_{j}(t)$ and $r_{j}(t)$ at an arbitrary spatial point $x_{j}$, we solve the corresponding MZ equations. Such equations are of the form (88), where the memory kernel $K(t-s)$ is computed from first-principles (i.e., from the microscopic equations of motion) using the algorithm we presented in Section 2.5. In Figure 5, we compare the temporal auto-correlation function we obtained for the wave displacement $r_{j}(t)$ with results of Markov-Chain-Monte-Carlo (MCMC) ($10^{6}$ sample paths) for FPU systems with mild nonlinearities ($\beta_{1}=0.01$ and $\beta_{1}=0.1$ ) at different temperatures ($\gamma=1$ and $\gamma=40$). It is seen that the MZ-Faber approximation of the MZ memory kernel yields relatively accurate results for FPU systems with mild nonlinearties at both low ($\gamma=40$) and high temperature ($\gamma=1$) as we increase the polynomial order $n$.

Reduced-order stochastic modeling

We employ the spectral approach of Section 3 to build stochastic low-dimensional models of the wave momentum $p_{j}(t)$ and wave displacement $r_{j}(t)=u_{j+1}(t)-u_{j}(t)$ at statistical equilibrium. Since we assumed that we are at statistical equilibrium, the statistical properties of the random processes representing $p_{j}(t)$ and $r_{j}(t)$ are time-independent. For instance, by integrating (94) we obtain the following expression for the one-time PDF of $r_{j}(t)$

[TABLE]

Clearly, $r_{j}(t)$ is a stationary non-Gaussian process. To sample the KL expansion of $r_{j}(t)$ in a way that is consistent with the PDF (95) we used the algorithm discussed in [49, 50]. For the FPU system with $\alpha_{1}=\beta_{1}=1$, it is straightforward to show that for all $m\in\mathbb{N}$

[TABLE]

where $\Gamma(x)$ is the Gamma function, $K_{n}(z)$ is the modified Bessel function of the second kind and $U(x,y,z)$ is Tricomi’s confluent hypergeometric function. Therefore, for all positive $\gamma$ and finite $m$ we have that $\mathbb{E}\{r_{j}^{2m}(t)\}<\infty$, i.e., $r_{j}(t)$ is $L^{2m}$ process. This condition guarantees convergence of the KL expansion to temporal correlation functions of order greater than two (see Appendix A). In Figure 6 we plot the temporal auto-correlation function of various polynomial observables of the nonlinear wave momentum and displacement at an arbitrary spatial point $x_{j}$. We compare results we obtained from Markov Chain Monte Carlo simulation (dashed line), with the MZ-KL expansion method based the first-principle memory calculation (continuous line). We also provide results we obtained by using KL expansions with covariance kernel estimated from data (dotted line).

5 Summary

We developed a new method to approximate the Mori-Zwanzig (MZ) memory integral in generalized Langevin equations (GLEs) describing the evolution of smooth observables in high-dimensional nonlinear systems with local interactions. The new method is based on Faber operator series expansions [68], and a formally exact combinatorial algorithm that allows us to compute the expansion coefficients of the MZ memory from first principles, i.e., based on the microscopic equations of motion. We also developed a new stochastic modeling technique that employs Karhunen-Loève expansions to represent the MZ fluctuation term (random noise) for systems in statistical equilibrium. We demonstrated the MZ memory calculation method and the MZ-KL stochastic modeling technique in applications to random wave propagation and prototype problems in classical statistical mechanics such as the Fermi-Pasta-Ulam $\beta$-model. We found that the proposed algorithms can accurately capture relaxation to statistical equilibrium in systems with mild nonlinearities, and in strongly nonlinear systems at high-temperature. At low temperature the Faber expansion of the MZ memory kernel is granted to converge only on a time interval that depends on the system and on the observable. In particular, Corollary 3.4.3 in [67] establishes short-time convergence of the MZ-Faber memory approximation for a broad class of nonlinear systems of the form (1). This implies that the MZ-Faber cumulant expansion can exhibit short-time convergence, meaning that it produces first-principle results that are accurate only for relatively short integration times.

We conclude by emphasizing that the mathematical techniques we presented can be readily applied to more general systems with local interactions such as particle systems modeling the microscopic dynamics of solids and liquids [66, 37, 38]. This opens the possibility to build new approximation schemes for MZ equations and derive new types of coarse-grained models where the MZ memory is constructed from first-principles and the fluctuation term is modeled stochastically.

Acknowledgements This research was supported by the Air Force Office of Scientific Research (AFOSR) grant FA9550-16-586-1-0092.

Appendix A Auto-correlation function of polynomial observables

In this Appendix we prove that the temporal auto-correlation function of phase space functions of the form $u^{n}(t)=u^{n}(\bm{x}(t,\bm{x}_{0}))$, i.e.,

[TABLE]

can be represented by replacing $u(t)$ with the KL expansion (66), and then sending $K$ to infinity. This result allows us to compute the auto-correlation function of $u^{n}(t)$ based on the KL expansions of $u(t)$.

Theorem A.1 Consider a zero-mean stationary stochastic process $u(t)$, $t\in[0,T]$, and assume that it has finite joint moments up to any desired order. Let

[TABLE]

be the truncated Karhunen-Loève expansion of $u(t)$. Then

[TABLE]

*i.e., $\langle u_{K}^{n}(0),u_{K}^{n}(t)\rangle_{\rho}$ converges uniformly to $\langle u^{n}(0),u^{n}(t)\rangle_{\rho}$ as $K\rightarrow\infty$. *

Proof.

Let us define

[TABLE]

The first term at the right hand side is of the form

[TABLE]

By using the Cauchy-Schwarz inequality, we obtain

[TABLE]

where we defined $\epsilon_{K}(t)=\|u_{K}(t)-u(t)\|_{L^{2}}$. It is well-known that $\epsilon_{K}(t)\rightarrow 0$ as $K\rightarrow\infty$ (see, e.g., [41]). By using the generalized Hölder’s inequality $\|fg\|_{L^{p}}\leq\|f\|_{L^{q}}\|g\|_{L^{q}}$, where $2p=q$ and the Minkowski inequality, we obtain

[TABLE]

where

[TABLE]

Similarly, we have

[TABLE]

By combining (99), (A) and (102), we obtain

[TABLE]

which proves the theorem.

∎