Homogenization for Generalized Langevin Equations with Applications to Anomalous Diffusion

TL;DR

This paper develops a homogenization theory for generalized Langevin equations with multiple time scales, revealing complex drift corrections and non-Markovian noise in the limits relevant to anomalous diffusion modeling.

Contribution

It introduces a general homogenization theorem for GLEs with multiple scales, including small mass and memory time scales, and applies it to models of anomalous diffusion.

Findings

Limiting SDEs have non-trivial drift corrections.

Limits involve non-Markovian noise processes.

Results applicable to particle diffusion models.

Abstract

We study homogenization for a class of generalized Langevin equations (GLEs) with state-dependent coefficients and exhibiting multiple time scales. In addition to the small mass limit, we focus on homogenization limits, which involve taking to zero the inertial time scale and, possibly, some of the memory time scales and noise correlation time scales. The latter are meaningful limits for a class of GLEs modeling anomalous diffusion. We find that, in general, the limiting stochastic differential equations (SDEs) for the slow degrees of freedom contain non-trivial drift correction terms and are driven by non-Markov noise processes. These results follow from a general homogenization theorem stated and proven here. We illustrate them using stochastic models of particle diffusion.

Equations

$$md{\mathbf{v}}_t={\mathbf{F}}_0\left(t,{\mathbf{x}}_t,{\mathbf{v}}_t,{\mathbf{\eta }}_t\right)dt+{\mathbf{F}}_1\left(t,\{{\mathbf{x}}_s,{\mathbf{v}}_s{\}}_{s\in [0,t]},{\mathbf{\xi }}_t\right)dt+{\mathbf{F}}_e(t,{\mathbf{x}}_t)dt.$$

$${\mathbf{F}}_0\left(t,{\mathbf{x}}_t,{\mathbf{v}}_t,{\mathbf{\eta }}_t\right)dt=-{\mathbf{\gamma }}_0(t,{\mathbf{x}}_t){\mathbf{v}}_{t}dt+{\mathbf{\sigma }}_0(t,{\mathbf{x}}_t)d{\mathbf{W}}_t^{(k)},$$

$${\mathbf{F}}_1\left(t,\{{\mathbf{x}}_s,{\mathbf{v}}_s{\}}_{s\in [0,t]},{\mathbf{\xi }}_t\right)=-\mathbf{g}(t,{\mathbf{x}}_t)\left({\int }_0^t\mathbf{\kappa }(t-s)\mathbf{h}(s,{\mathbf{x}}_s){\mathbf{v}}_{s}ds\right)+\mathbf{\sigma }(t,{\mathbf{x}}_t){\mathbf{\xi }}_t,$$

$$\mathbf{R}(t):=E[{\mathbf{\xi }}_t{\mathbf{\xi }}_0^{\ast }]={\mathbf{R}}^{\ast }(-t)\in {\mathbb{R}}^{r\times r},$$

$$\mathcal{S}(\omega )={\int }_{-\infty }^{\infty }\mathbf{R}(t)e^{-i\omega t}dt.$$

$$d{x}_t$$

$$md{v}_t$$

$$\kappa (t)=\frac{{\Gamma }_2^2({\Gamma }_2{e}^{-{\Gamma }_2|t|}-{\Gamma }_1{e}^{-{\Gamma }_1|t|})}{2({\Gamma }_2^2-{\Gamma }_1^2)},$$

$$\mathcal{S}(\omega )=\frac{{\Gamma }_2^2{\omega }^2}{({\omega }^2+{\Gamma }_1^2)({\omega }^2+{\Gamma }_2^2)}.$$

$$\kappa (t)=\frac{1}{2}(\delta (t)-{\Gamma }_1{e}^{-{\Gamma }_1|t|}),$$

$$\mathcal{S}(\omega )=\frac{{\omega }^2}{{\omega }^2+{\Gamma }_1^2}.$$

$$d{x}_t$$

$$md{v}_t$$

$$\kappa (t)=\frac{{\beta }^2}{2}(\delta (t)-{\Gamma }_1{e}^{-{\Gamma }_1|t|}),$$

$$\mathcal{S}(\omega )=\frac{{\beta }^2{\omega }^2}{{\omega }^2+{\Gamma }_1^2},$$

$$d{X}_t$$

$$d{Y}_t$$

$$d{Z}_t$$

$$S_1(X)$$

$$d{X}_t$$

$$d{U}_t^{\varphi }$$

$$d{x}_t$$

$$md{v}_t$$

$$\kappa (t)=\frac{{\beta }^2{\Gamma }_2^2({\Gamma }_2{e}^{-{\Gamma }_2|t|}-{\Gamma }_1{e}^{-{\Gamma }_1|t|})}{2({\Gamma }_2^2-{\Gamma }_1^2)},$$

$$\mathcal{S}(\omega )=\frac{{\beta }^2{\Gamma }_2^2{\omega }^2}{({\omega }^2+{\Gamma }_1^2)({\omega }^2+{\Gamma }_2^2)}.$$

$$d{x}_t^{ϵ}$$

$$m_0ϵd{v}_t^{ϵ}$$

$${\kappa }^{ϵ}(t)=\frac{{\beta }^2{\gamma }_2^2(\frac{{\gamma }_2}{ϵ}{e}^{-\frac{{\gamma }_2}{ϵ}|t|}-{\Gamma }_1{e}^{-{\Gamma }_1|t|})}{2({\gamma }_2^2-{ϵ}^2{\Gamma }_1^2)},$$

$${\kappa }^{ϵ}(t)\to \kappa (t):=\frac{{\beta }^2}{2}(\delta (t)-{\Gamma }_1{e}^{-{\Gamma }_1|t|})$$

$$d{X}_t$$

$$d{Y}_t$$

$$d{Z}_t$$

$$S_1$$

$$+\frac{\partial }{\partial X}\left(\frac{\sigma }{gh}\right)\frac{4\sigma }{{\beta }^{2}gh+4m_0{\gamma }_2},$$

$$S_2$$

$$d{X}_t$$

Full text

1. 1 Introduction

2. 1.1 Motivation

3. 1.2 Definitions and Models

4. 1.3 Goals, Organization and Summary of Results of the Paper

5. 2 Application to One-Dimensional GLE Models

6. 3 GLEs in Finite Dimensions

7. 3.1 Generalized Langevin Systems

8. 3.2 Generalized Langevin Systems as Input-Output Stochastic Dynamical Systems with Multiple Time Scales

9. 4 On the Homogenization of Generalized Langevin Dynamics

10. 5 Small Mass Limit of Generalized Langevin Dynamics

11. 6 Homogenization for the Case of Vanishing Effective Damping Constant and Effective Diffusion Constant

12. 7 Conclusions and Final Remarks

13. A Homogenization for a Class of SDEs with State-Dependent Coefficients

14. B Proof of Theorem A.6

15. C An Implementation of Algorithm 6.3 under Assumption 6.4

16. D Another Example for Section 2

1. Introduction

1.1. Motivation

Most of the mathematical models of diffusion phenomena use noise which is white (i.e. uncorrelated), or Markovian [53]. The present paper is a step towards removing this limitation. The diffusion models studied here are driven by noises, belonging to a wide class of non-Markov processes. A standard example of Markovian noise is a multidimensional Ornstein-Uhlenbeck process. An important class of Gaussian stochastic processes is obtained by linear transformations of multidimensional Ornstein-Uhlenbeck processes. The covariance (equal to correlation in the case of zero mean) of such a process is a linear combination of exponentials decaying and possibly oscillating on different time scales, and its spectral density (power spectrum) is a ratio of two semi-positive defined polynomials [16]. In cases when the polynomial in the denominator has degenerate zeros, the covariance contains products of exponentials and polynomials in time. This is a very general class of processes—every stationary Gaussian process whose covariance is a Bohl function (see Section 2) can be obtained as a linear transformation of an Ornstein-Uhlenbeck process in some (finite) dimension. In general, these processes are not Markov.

Let us mention here the seminal result by L.A. Khalfin from 1957 [32], who showed, quite generally that in any system with energy spectrum bounded from below (which is a necessary condition for the physical stability), correlations must decay no faster than according to a power law. To this day this result provides inspirations and motivations for further studies in the context of thermalization [70], cooling of atoms in photon reservoirs [40], decay of metastable states as monitored by luminescence [63], or quantum anti-Zeno effect (c.f. [59, 41]), to name a few examples. Khalfin’s result further motivates studying systems with non-Markovian noise, as most natural examples of strongly correlated processes do not satisfy Markov property.

While the noise processes studied here have exponentially decaying covariances, their class is very rich and they may be useful in approximating strongly correlated noises on time intervals, relevant for studied phenomena [67]. In addition, as discussed in more detail later, generalization of the method applied here may lead to a representation of a class of noises whose covariances decay as powers (see Remark 3.7). Also, the representation of spectral density of the noise processes as ratio of two polynomials is convenient in applications, in particular for solving the problem of predicting (in the least mean square sense) a colored noise process given observations on a finite segment of the past or on the full past [16].

1.2. Definitions and Models

We consider the following stochastic model for a particle (for instance, Brownian particle or a tagged tracer particle) interacting with the environment (for instance, a heat bath or a viscous fluid). Let $\boldsymbol{x}_{t}\in\mathbb{R}^{d}$ denote the particle’s position, where $t\geq 0$ denotes time and $d$ is a positive integer. The evolution of the particle’s velocity, $\boldsymbol{v}_{t}:=\dot{\boldsymbol{x}}_{t}\in\mathbb{R}^{d}$, is described by the following generalized Langevin equation (GLE):

[TABLE]

In the above, $m>0$ is the particle’s mass, $\boldsymbol{\eta}_{t}$ is a $k$-dimensional Gaussian white noise satisfying $E[\boldsymbol{\eta}_{t}]=\boldsymbol{0}$ and $E[\boldsymbol{\eta}_{t}\boldsymbol{\eta}_{s}^{*}]=\delta(t-s)\boldsymbol{I}$, and $\boldsymbol{\xi}_{t}$ is a colored noise process independent of $\boldsymbol{\eta}_{t}$. Here and throughout the paper, the superscript ∗ denotes transposition of matrices or vectors, $\boldsymbol{I}$ denotes identity matrix of appropriate dimension, $E$ denotes expectation, and $\mathbb{R}^{+}:=[0,\infty)$. The initial data are random variables, $\boldsymbol{x}_{0}=\boldsymbol{x}$, $\boldsymbol{v}_{0}=\boldsymbol{v}$, independent of $\{\boldsymbol{\xi}_{t},t\in\mathbb{R}^{+}\}$ and $\{\boldsymbol{\eta}_{t},t\in\mathbb{R}^{+}\}$.

The three terms on the right hand side of (1.1) model forces of different physical natures acting on the particle.

• (i)

$\boldsymbol{F}_{e}$ is an external force field, which may be conservative (potential) or not.

• (ii)

$\boldsymbol{F}_{0}$ is a Markovian force of the form

[TABLE]

containing an instantaneous damping term and a multiplicative white noise term. The damping and noise coefficients, $\boldsymbol{\gamma}_{0}:\mathbb{R}^{+}\times\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ and $\boldsymbol{\sigma}_{0}:\mathbb{R}^{+}\times\mathbb{R}^{d}\to\mathbb{R}^{d\times k}$, may depend on the particle’s position and on time. $\boldsymbol{W}^{(k)}_{t}$ denotes a $k$-dimensional Wiener process—the time integral of the white noise $\boldsymbol{\eta}_{t}$.

• (iii)

$\boldsymbol{F}_{1}$ is a non-Markovian force of the form

[TABLE]

containing a non-instantaneous damping term, describing the delayed drag effects by the environment on the particle, and a multiplicative colored noise term. The coefficients, $\boldsymbol{g}:\mathbb{R}^{+}\times\mathbb{R}^{d}\to\mathbb{R}^{d\times q}$, $\boldsymbol{h}:\mathbb{R}^{+}\times\mathbb{R}^{d}\to\mathbb{R}^{q\times d}$ and $\boldsymbol{\sigma}:\mathbb{R}^{+}\times\mathbb{R}^{d}\to\mathbb{R}^{d\times r}$, depend in general on the particle’s position and on time. In the above, $q$ and $r$ are positive integers, and the memory function $\boldsymbol{\kappa}:\mathbb{R}\to\mathbb{R}^{q\times q}$ is a real-valued function that decays sufficiently fast at infinities. $\boldsymbol{\xi}_{t}\in\mathbb{R}^{r}$ is a mean-zero stationary Gaussian vector process, to be defined in detail later. The statistical properties of the process $\boldsymbol{\xi}_{t}$ are completely determined by its (matrix-valued) covariance function,

[TABLE]

or equivalently, by its spectral density, $\boldsymbol{\mathcal{S}}(\omega)$, i.e. the Fourier transform of $\boldsymbol{R}(t)$ defined as:

[TABLE]

For simplicity, we have omitted other forces such as the Basset force [24] from Eqn. (1.1). Note that $\boldsymbol{F}_{0}$ and $\boldsymbol{F}_{1}$ describe two types of forces associated with different physical mechanisms. Of particular interest is when the noise term in $\boldsymbol{F}_{0}$ and $\boldsymbol{F}_{1}$ models environments of different nature (passive bath and active bath respectively [14]) that the particle interacts with.

As the name itself suggests, GLEs are generalized versions of the Markovian Langevin equations, frequently employed to model physical systems. A basic form of the GLEs was first introduced by Mori in [52] and subsequently used in numerous statistical physics models [35, 71, 75]. The studies of GLEs have attracted increasing interest in recent years. We refer to, for instance, [49, 45, 68, 26, 23, 46, 38, 74, 69] for various applications of GLEs and [55, 48, 21, 39] for their asymptotic analysis. The main merit of GLEs from modeling point of view is that they take into account the effects of memory and the colored nature of noise on the dynamics of the system.

Remark 1.1*.*

In general, there need not be any relation between $\boldsymbol{\kappa}(t)$ and $\boldsymbol{R}(t)$, or any relation between the damping coefficients and the noise coefficients appearing in the formula for $\boldsymbol{F}_{0}$ and $\boldsymbol{F}_{1}$. A particular but important case that we will revisit often in this paper is the case when a fluctuation-dissipation relation holds. In this case, $\boldsymbol{\gamma}_{0}$ is proportional to $\boldsymbol{\sigma}_{0}\boldsymbol{\sigma}_{0}^{*}$, $\boldsymbol{h}=\boldsymbol{g}^{*}$, $\boldsymbol{g}$ is proportional to $\boldsymbol{\sigma}$ and (without loss of generality111The factor $k_{B}T$, where $T$ is the absolute temperature and $k_{B}$ denotes the Boltzmann constant, is here set to $1$. In general, it can be absorbed into either one of the coefficients $\boldsymbol{g}$, $\boldsymbol{h}$ or $\boldsymbol{\sigma}$.) $\boldsymbol{R}(t)=\boldsymbol{\kappa}(t)$. Studies of microscopic Hamiltonian models for open classical systems lead to GLEs of the form (1.1) satisfying the above fluctuation-dissipation relation (see, for instance, Appendix A of [42] or [11]). On another note, GLEs of the form (1.1) are extended versions of the ones studied in our previous work [42] – here the GLEs are generalized to include a Markovian force, in addition to the non-Markovian one, as well as explicit time dependence in the coefficients.

As a motivation, we now provide and elaborate on examples of systems that can be modeled by our GLEs.

An important type of diffusion, which has been observed in many physical systems, from charge transport in amorphous materials to intracellular particle motion in cytoplasm of living cells [62], is ballistic diffusion. It is a subclass of anomalous diffusions and is characterized by the property that the particle’s long-time mean-squared displacement grows quadratically in time – in contrast to linear growth in usual diffusion. There are many different theoretical models of anomalous diffusion with diverse properties, coming from different physical assumptions; see [50] for a comprehensive survey. In the following, we provide two GLE models that are employed to study such phenomena. Their properties will be studied in Section 2, as an application of the results proven here.

Example 1*.*

Two GLE models for anomalous diffusion of a free Brownian particle in a heat bath. A large class of models for diffusive systems is described by the system of equations (for simplicity, we restrict to one dimension):

[TABLE]

where $x_{t},\ v_{t}\in\mathbb{R}$ are the position and velocity of the particle, $\kappa(t)$ is called the memory function, and $\xi_{t}$ is a mean-zero stationary Gaussian process.

Two particular GLE models are described by (1.6)-(1.7), with:

• (M1)

memory function of the bi-exponential form:

[TABLE]

where the parameters satisfy $\Gamma_{2}>\Gamma_{1}>0$, and $\xi_{t}$ has the covariance function $R(t)=\kappa(t)$ and thus the spectral density,

[TABLE]

This model is similar to the one first introduced and studied in [3]. The noise with the above covariance function can be realized by the difference between two Ornstein-Uhlenbeck processes, with different damping rates, driven by the same white noise. Various properties as well as applications of GLEs of the form (1.6)-(1.7) were studied in [3, 2, 68].

• (M2)

memory function of the form:

[TABLE]

where $\Gamma_{1}>0$, and $\xi_{t}$ has the covariance function $R(t)=\kappa(t)$ and thus the spectral density,

[TABLE]

This model can be obtained from the one in (M1) by sending $\Gamma_{2}\to\infty$ in the formula for $\kappa(t)$ in (1.8).

Observe that the spectral densities in both models share the same asymptotic behavior near $\omega=0$, i.e. $\mathcal{S}(\omega)\sim\omega^{2}$ as $\omega\to 0$, contributing to the enhanced diffusion (super-diffusion) of the particle with mean-squared displacement growing as $t^{2}$ as $t\to\infty$ [68]. See Proposition 3.5 for a precise argument.

Other examples of systems that can be modeled by our GLEs are multiparticle systems with hydrodynamic interaction [17], active matter systems [65], among others. Although our main results are applicable to these systems, we will not pursue the study of these systems here.

1.3. Goals, Organization and Summary of Results of the Paper

Goals of the Paper. We aim to derive homogenized models for a general class of GLEs (see Section 3), containing the examples (M1) and (M2) as special cases (see Corollary 2.1 and Corollary 2.2). This will allow us to gain insights into the stochastic dynamics of such systems, including many systems that exhibit anomalous diffusion (see discussion in the paragraph before Example 1) – this is, in fact, the main motivation of the present paper. To the best of our knowledge, this is the first work that studies homogenization for GLE models describing anomalous diffusion.

Given a GLE system, it is often desirable to work with simpler, reduced models that capture the essential features of its dynamics. To obtain satisfactory and optimal models, one needs to take into account the trade-off between the simplicity and accuracy of the reduced models sought after. Indeed, one may find that a reduced model, while simplified, fails to give a physically correct model for describing a system of interest [64]. Two successful reductions were carried out in [28] for the case $\boldsymbol{F}_{1}=\boldsymbol{0}$ and in [42] for the case $\boldsymbol{F}_{0}=\boldsymbol{0}$.

One of our main goals in this paper is to devise and study new homogenization procedures that yield reduced models retaining essential features of a more general class of models. This program is of importance for identification, parameter inference and uncertainty quantification of stochastic systems [61, 25, 45, 38] arising in the studies of anomalous diffusion [49, 51], climate modeling [22, 47] and molecular systems [10], among others. There is increasing amount of effort striving to implement this or related programs, starting from microscopic models [60], using various techniques [20, 57, 7, 19, 26], for different systems of interest in the literature. The derived effective SDE models will be of particular interest for modelers of anomalous diffusion.

Organization of the Paper. The paper is organized as follows. We first present the application of the results obtained in the later sections (Section 5 and Section 6) to study homogenization of generalized versions of the one-dimensional models (M1) and (M2) from Example 1 in Section 2. Since these results are easier to state and require minimal notation to understand, we have chosen to present them as early as possible to demonstrate the value of our study to application-oriented readers. The later sections study an extended, multi-dimensional version of the GLEs in Section 2. In Section 3 we introduce the GLEs to be studied and revisit them from the perspective of input-output stochastic dynamical systems exhibiting multiple time scales. In Section 4, we discuss various ways of homogenizing GLEs. Following this discussion, we study the small mass limit of the GLEs in Section 5. We introduce and study novel homogenization procedures for a class of GLEs in Section 6. We state conclusions and make final remarks in Section 7. Relevant technical details and supplementary materials are provided in the appendix. In particular, we state a homogenization theorem for a general class of SDEs with state-dependent coefficients in Appendix A. The proof of this theorem is given in Appendix B.

Summary of the Main Results. For reader’s convenience, below we list (not in exactly the same order as the results appear in the paper) and summarize the main results obtained in the paper.

• •

The first main result is Theorem 5.4. It studies the small mass limit of the GLE described by (5.1)-(5.2). It states that the position process converges, in a strong pathwise sense, to a component of a higher dimensional process satisfying an Itô SDE. The SDE contains non-trivial drift correction terms. We stress that, while being a component of a Markov process, the limiting position process itself is not Markov. This is in constrast to the nature of limiting processes obtained in earlier works, the difference which holds interesting implications from a physical point of view (recall the discussion after Eqn. (1.5)). Therefore, Theorem 5.4 constitutes a novel result, both mathematically and physically.

• •

The second main result is Theorem 6.7. It describes the homogenized behavior of a family of GLEs (Eqns. (6.18)-(6.19)), parametrized by $\epsilon>0$, in the limit as $\epsilon\to 0$. This limit is equivalent to the limit in which the inertial time scale, some of the memory time scales and some of the noise correlation time scales in the pre-limit system, tend to zero at the same rate. As in Theorem 5.4, the result here states that the position process converges, in a strong pathwise sense, to a component of a higher dimensional process satisfying an Itô SDE which contains non-trivial drift correction terms. Again, the limiting position process is non-Markov. However, the structure of the SDE is rather different from the one obtained in Theorem 5.4. As discussed later, this result holds interesting consequences for systems exhibiting anomalous diffusion.

• •

The third and forth main result are Corollary 2.1 and Corollary 2.2. These results specialize the earlier ones to one-dimensional GLE models, which are generalizations of (M1) and (M2), and follow from the earlier theorems. They give explicit expressions for the drift correction terms present in the limiting SDEs and therefore may be used directly for modeling and simulation purposes. Furthermore, we show that, in the important case where the fluctuation-dissipation relation (see Remark 1.1) holds, the two corollaries are intimately connected. Recall that these results are going to be presented first in Section 2.

• •

The last main result is Theorem A.6, on homogenization of a family of parametrized SDEs whose coefficients are state-dependent. These SDEs are variants of the ones studied in earlier works [28, 4, 6]. In comparison with all the earlier studies, the state-dependent coefficients of the pre-limit SDEs (A.3)-(A.4) may depend on the parameter $\epsilon>0$ (to be taken to zero) explicitly. Therefore, this result is new and not simply a minor generalization of earlier results. Moreover, it is important in the context of present paper and is needed here to study various homogenization limits of GLEs, the importance of which is evident in the discussions above, in the main paper.

2. Application to One-Dimensional GLE Models

We first study the small mass limit of a one-dimensional GLE, which is a generalized version of the GLE in model (M2) of Example 1, modeling super-diffusion of a particle in a heat bath. Our models are generalized in that the coefficients of the GLEs are state-dependent. For simplicity, we are going to omit the explicit time dependence in the damping and noise coefficients—but not in the external force.

For $t\in\mathbb{R}^{+}$, $m>0$, let $x_{t},v_{t}\in\mathbb{R}$ be the solutions to the equations:

[TABLE]

where

[TABLE]

where $\Gamma_{1}>0$, and $\xi_{t}$ is the mean-zero stationary Gaussian process with the covariance function $R(t)=\kappa(t)$ and spectral density,

[TABLE]

The initial data $(x,v)$ are random variables independent of $\epsilon$ and have finite moments of all orders.

The following corollary describes the limiting SDE for the particle’s position obtained in the small mass limit of (2.1)-(2.2).

Corollary 2.1**.**

Assume that for every $y\in\mathbb{R}$, $g(y),g^{\prime}(y),h(y),h^{\prime}(y)$, $\sigma(y)$ are bounded continuous functions in $y$, $F_{e}(t,y)$ is bounded and continuous in $t$ and $y$, and all the listed functions have bounded $y$-derivatives. Then in the limit $m\to 0$, the particle’s position, $x_{t}\in\mathbb{R}$, satisfying (2.1)-(2.2), converges to $X_{t}$, where $X_{t}$ solves the following Itô SDE:

[TABLE]

where

[TABLE]

Moreover, if in addition $g:=\phi\sigma$, where $\phi>0$, then the number of limiting SDEs reduces from three to two:

[TABLE]

where $U_{t}^{\phi}=\phi Y_{t}-Z_{t}$.

The convergence is in the sense that for every $T>0$, $\sup_{t\in[0,T]}|x_{t}-X_{t}|\to 0$ in probability as $m\to 0$.

Proof.

We apply Theorem 5.4 by setting $d=1,d_{2}=d_{4}=2$, $\alpha_{1}=\alpha_{3}=0$, $\alpha_{2}=\alpha_{4}=1$, $\boldsymbol{\gamma}_{0}=\beta^{2}gh/2$, $\boldsymbol{\sigma}_{0}=\beta\sigma$, $\boldsymbol{h}=h$, $\boldsymbol{g}=g$, $\boldsymbol{\sigma}=\sigma$, $\boldsymbol{C}_{2}=\boldsymbol{C}_{4}=\beta$, $\boldsymbol{\Gamma}_{2}=\Gamma_{1}$, $\boldsymbol{M}_{2}\boldsymbol{C}_{2}^{*}=-\Gamma_{1}\beta/2$, $\boldsymbol{\Gamma}_{4}=\Gamma_{1}$, $\boldsymbol{\Sigma}_{4}=-\Gamma_{1}$, and $\boldsymbol{F}_{e}=F_{e}$. The assumptions of Theorem 5.4 can be verified in a straightforward way and so the results of the corollary follow. ∎

We next specialize the result of Theorem 6.7 to study homogenization of one-dimensional GLEs which are generalizations of the model (M1) in Example 1: for $t\in\mathbb{R}^{+}$, $m>0$, let $x_{t},v_{t}\in\mathbb{R}$ be the solutions to the equations:

[TABLE]

where

[TABLE]

with $\Gamma_{2}>\Gamma_{1}>0$, and $\xi_{t}$ is the mean-zero stationary Gaussian process with the covariance function $R(t)=\kappa(t)$ and spectral density,

[TABLE]

The initial data $(x,v)$ are random variables independent of $\epsilon$ and have finite moments of all orders.

For $\epsilon>0$, we set $m=m_{0}\epsilon$ and $\Gamma_{2}=\gamma_{2}/\epsilon$ in (2.11)-(2.12), where $m_{0}$ and $\gamma_{2}$ are positive constants. This gives the family of equations:

[TABLE]

where

[TABLE]

and $\xi^{\epsilon}_{t}$ is the family of mean-zero stationary Gaussian processes with the covariance functions, $R^{\epsilon}(t)=\kappa^{\epsilon}(t)$.

Discussion. We discuss the physical meaning behind the above rescaling of parameters. Recall that in the first case of Example 1 (i.e. the model (M1)), the mean-square displacement of the particle grows as $t^{2}$ as $t\to\infty$ and therefore the above model describes a particle exhibiting super-diffusion. As $\epsilon\to 0$, the environment allows for more and more negative correlation and in the limit the covariance function consists of a delta-type peak at $t=0$ and a negative long tail compensating for the positive peak when integrated (see Figure 1 and also page 105 of [71]). Indeed,

[TABLE]

as $\epsilon\to 0$. This is the so-called vanishing effective friction case in [1]. The noise with the covariance function $\kappa^{\epsilon}(t)$ is called harmonic velocity noise, whereas the noise with the covariance function $\kappa(t)$ is the derivative of an Ornstein-Uhlenbeck process.

The following corollary provides the homogenized model in the limit $\epsilon\to 0$ of (2.15)-(2.16).

Corollary 2.2**.**

Assume that for every $y\in\mathbb{R}$, $g(y),g^{\prime}(y),h(y),h^{\prime}(y)$, $\sigma(y)$ are bounded continuous functions in $y$, $F_{e}(t,y)$ is bounded and continuous in $t$ and $y$, and all the listed functions have bounded derivatives in $y$. Then in the limit $\epsilon\to 0$, the particle’s position, $x^{\epsilon}_{t}\in\mathbb{R}$, satisfying (2.15)-(2.16), converges to $X_{t}$, where $X_{t}$ solves the following Itô SDE:

[TABLE]

where $g=g(X_{t})$, $h=h(X_{t})$, $\sigma=\sigma(X_{t})$, $W_{t}$ is a one-dimensional Wiener process, and

[TABLE]

Moreover, if in addition $g:=\phi\sigma$, where $\phi>0$, then the number of limiting SDEs reduces from three to two:

[TABLE]

where $U_{t}^{\phi}=\phi Y_{t}-Z_{t}$.

The convergence is in the sense that for every $T>0$, $\sup_{t\in[0,T]}|x^{\epsilon}_{t}-X_{t}|\to 0$ in probability as $\epsilon\to 0$.

Proof.

Let $d=1$, $d_{2}=d_{4}=2$ and denote the one-dimensional version of the variables, coefficients and parameters in Theorem 6.7 by non-bold letters (for instance, $x_{t}$, $B_{2}$, $\Gamma_{2,2}$ etc.). Furthermore, set $B_{2}=B_{4}=\beta>0$, $\gamma_{2,2}=\gamma_{4,2}=\gamma_{2}>0$ and $\Gamma_{2,1}=\Gamma_{4,1}=\Gamma_{1}$. Then it can be verified that the assumptions of Theorem 6.7 hold and the results follow upon solving a Lyapunov equation. ∎

Remark 2.3*.*

A few remarks on the contents of Corollary 2.2 follow.

• (i)

the homogenized position process is non-Markov, driven by a colored noise process which is the derivative of the Ornstein-Uhlenbeck process. This behavior is expected in view of the asymptotic behavior of the rescaled memory function and spectral density as $\epsilon\to 0$.

• (ii)

similarly to the small mass limit case considered earlier, the limiting equation for the particle’s position not only contains noise-induced drift terms but is also coupled to equations for other slow variables. Moreover, the limiting equations for these other slow variables also contain non-trivial correction terms – the memory induced drift.

Relation between Corollary 2.1 and Corollary 2.2. The limiting SDE systems in Corollaries 2.1 & 2.2 are generally different because of the different correction drift terms $S_{1}$ and $S_{2}$. In other words, sending $\Gamma_{2}\to\infty$ first in (2.11)-(2.12) and then taking $m\to 0$ of the resulting GLE does not, in general, give the same limiting SDE as taking the joint limit of $m\to 0$ and $\Gamma_{2}\to\infty$. However, if one further assumes that $g$ is proportional to $\sigma$, then the limiting SDE systems coincide. An important particular case is when $g=h=\sigma$, in which case a fluctuation-dissipation relation holds and the GLE can be derived from a microscopic Hamiltonian model (see Remark 1.1). In this case, the homogenized model described in both corollaries reduces to:

[TABLE]

To end this section, we remark that one could in principle repeat the above analysis for the case where the spectral density varies as $\omega^{2l}$, for $l=2,4,\dots$ (i.e. the highly nonlinear case) as well as extending the studies done so far in various other directions. To illustrate how non-trivial the calculations and results could become, we work out another example in Appendix D.

3. GLEs in Finite Dimensions

We call a system modeled by GLE of the form (1.1) a generalized Langevin system. Its dynamics will be referred to as generalized Langevin dynamics.

We assume that the memory function $\boldsymbol{\kappa}(t)$ in the GLE (1.1) is a Bohl function, i.e. that each matrix element of $\boldsymbol{\kappa}(t)$ is a finite, real-valued linear combination of exponentials, possibly multiplied by polynomials and/or by trigonometric functions. The noise process, $\{\boldsymbol{\xi}(t),t\in\mathbb{R}^{+}\}$, is a mean-zero, mean-square continuous stationary Gaussian process with Bohl covariance function and, therefore, its spectral density $\boldsymbol{\mathcal{S}}(\omega)$ is a rational function—(see Theorem 2.20 in [72]). In this case, the generalized Langevin dynamics can be realized by an SDE system in a finite-dimensional space (see next subsection for details). The case in which an infinite-dimensional space is required is deferrred to a future work (see also Remark 3.7 and Section 7).

We recall a useful fact: given a rational spectral density $\boldsymbol{\mathcal{S}}(\omega)\in\mathbb{R}^{r\times r}$, there exists a rational function $\boldsymbol{G}(z)\in\mathbb{C}^{r\times l}$, called a spectral factor, such that $\boldsymbol{\mathcal{S}}(\omega)=\boldsymbol{G}(i\omega)\boldsymbol{G}^{*}(-i\omega)$. We emphasize that such factorization is not unique [43].

3.1. Generalized Langevin Systems

Below we define the memory function and the noise process in the GLE (1.1) (see Eqn. (1.3)), and along the way introduce our notation. They are defined in a manner ensuring simplicity as well as providing sufficient parameters for matching the memory function and the correlation function of the noise, thereby preserving the essential statistical properties of the GLE. This provides a systematic framework for our homogenization studies (see the discussion in Section 4).

For $i=1,2,3,4$, let $\boldsymbol{\Gamma}_{i}\in\mathbb{R}^{d_{i}\times d_{i}}$, $\boldsymbol{M}_{i}\in\mathbb{R}^{d_{i}\times d_{i}}$, $\boldsymbol{\Sigma}_{i}\in\mathbb{R}^{d_{i}\times q_{i}}$ be constant matrices. Also, let $\boldsymbol{C}_{i}\in\mathbb{R}^{q\times d_{i}}$ (for $i=1,2$) and $\boldsymbol{C}_{i}\in\mathbb{R}^{r\times d_{i}}$ (for $i=3,4$) be constant matrices. Here, the $d_{i}$ and $q_{i}$ ($i=1,2,3,4$) are positive integers. Let $\alpha_{i}\in\{0,1\}$ be a “switch on or off” parameter. We define the memory function in terms of the sextuple $(\boldsymbol{\Gamma}_{1},\boldsymbol{M}_{1},\boldsymbol{C}_{1};\boldsymbol{\Gamma}_{2},\boldsymbol{M}_{2},\boldsymbol{C}_{2})$ of matrices:

[TABLE]

The noise process is defined as:

[TABLE]

where the $\boldsymbol{\beta}^{j}_{t}\in\mathbb{R}^{d_{j}}$ ($j=3,4$) are independent Ornstein-Uhlenbeck type processes, i.e. solutions of the SDEs:

[TABLE]

with the initial conditions, $\boldsymbol{\beta}^{j}_{0}$, normally distributed with mean-zero and covariance $\boldsymbol{M}_{j}$. Here, $\boldsymbol{W}^{(q_{j})}_{t}$ denotes a $q_{j}$-dimensional Wiener process, independent of $\boldsymbol{\beta}^{j}_{0}$. Also, the Wiener processes $\boldsymbol{W}_{t}^{(q_{3})}$ and $\boldsymbol{W}_{t}^{(q_{4})}$ are independent.

For $i=1,2,3,4$, $\boldsymbol{\Gamma}_{i}$ is positive stable, i.e. all eigenvalues of $\boldsymbol{\Gamma}_{i}$ have positive real parts and $\boldsymbol{M}_{i}=\boldsymbol{M}_{i}^{*}>0$ satisfies the following Lyapunov equation:

[TABLE]

The $\boldsymbol{M}_{i}$ are therefore the steady-state covariances of the systems, i.e. the resulting Ornstein-Uhlenbeck processes are stationary. In control theory, $\boldsymbol{M}_{i}$ is also known as the controllability Gramian for the pair $(\boldsymbol{\Gamma}_{i},\boldsymbol{\Sigma}_{i})$ [72].

The covariance matrix, $\boldsymbol{R}(t)$, of the mean-zero Gaussian noise process is expressed by the sextuple $(\boldsymbol{\Gamma}_{3},\boldsymbol{M}_{3},\boldsymbol{C}_{3};\boldsymbol{\Gamma}_{4},\boldsymbol{M}_{4},\boldsymbol{C}_{4})$ of matrices as follows:

[TABLE]

and so the sextuple $(\boldsymbol{\Gamma}_{3},\boldsymbol{M}_{3},\boldsymbol{C}_{3};\boldsymbol{\Gamma}_{4},\boldsymbol{M}_{4},\boldsymbol{C}_{4})$, together with the parameters $\alpha_{3},\alpha_{4}$, completely determine the probability distributions of $\boldsymbol{\xi}_{t}$. We denote the spectral density of the noise process by $\boldsymbol{\mathcal{S}}(\omega)=\sum_{i=3,4}\alpha_{i}\boldsymbol{\mathcal{S}}_{i}(\omega)$, where $\boldsymbol{\mathcal{S}}_{i}(\omega)$ is the Fourier transform of $\boldsymbol{R}_{i}(t)$ for $i=3,4$.

We will view the system (3.2)-(3.3) (which is in a statistical steady state) as a representation of the noise process $\boldsymbol{\xi}_{t}$ and call such a representation a (finite-dimensional) stochastic realization of $\boldsymbol{\xi}_{t}$. Similarly, we view (3.1) as a representation of the memory function $\boldsymbol{\kappa}(t)$ and call such a representation a (finite-dimensional, deterministic) memory realization of $\boldsymbol{\kappa}(t)$. We call the Fourier transform of $\boldsymbol{\kappa}(t)$ and $\boldsymbol{R}(t)$ the spectral density of the memory function and spectral density of the noise process respectively.

An important message from the stochastic realization theory is that the system (3.2)-(3.3) is more than a representation of $\boldsymbol{\xi}_{t}$ in terms of a white noise, in that it also contains state variables $\boldsymbol{\beta}^{j}$ ($j=3,4$) which serve as a “dynamical memory”. In contrast to standard treatments, this dynamical memory comes not from one, but from two independent systems of type (3.3). This will be used to include two distinct types of dynamical memory that can be switched on or off using the parameters $\alpha_{i}$ – see Proposition 3.5. This consideration motivates us to define the memory function (and noise) explicitly using two independent systems, with different constraints on their parameters easier to state than if a single higher-dimensional system were used.

The sextuples that define the memory function in (3.1) and the noise process in (3.2) are only unique up to the following transformations:

[TABLE]

where $i=1,2,3,4$ and $\boldsymbol{T}_{i}$ are any invertible matrices of appropriate dimensions [43]. Different choices of $\boldsymbol{T}_{i}$ correspond to different coordinate systems.

Remark 3.1*.*

Realization of the memory function and noise process in terms of the matrix sextuples, as defined above, covers all GLEs driven by Gaussian processes that can be realized in a finite dimension (see the propositions and theorems on page 303-308 of [73]). See also the remarks on the subject in [42].

A summary of the above discussion is included in the following:

Assumption 3.2**.**

The memory function $\boldsymbol{\kappa}(t)$ in the GLE (1.1) is a real-valued Bohl function defined by (3.1) and the noise process, $\{\boldsymbol{\xi}_{t},t\in\mathbb{R}^{+}\}$, is a mean-zero, mean-square continuous, stationary Gaussian process with Bohl covariance function (hence, with rational spectral density), admitting a stochastic realization given by (3.2)-(3.3). Furthermore, we assume that any spectral factors $\boldsymbol{\Phi}_{i}(z)$ ($i=1,2,3,4$) of the spectral densities $\boldsymbol{\mathcal{S}}_{i}(\omega)$ are minimal (see Chapter 10 in [43]).

We introduce a generalized version of the effective damping constant and effective diffusion constant used in [42], which will be useful to study the asymptotic behavior of spectral densities.

Definition 3.3**.**

For $n\in\mathbb{Z}$, the $n$th order effective damping constant is defined as the constant matrix, parametrized by $\alpha_{1},\alpha_{2}\in\{0,1\}$:

[TABLE]

where $\boldsymbol{K}_{i}^{(n)}=\boldsymbol{C}_{i}\boldsymbol{\Gamma}_{i}^{-n}\boldsymbol{M}_{i}\boldsymbol{C}_{i}^{*}$ (for $i=1,2$). Likewise, the $n$th order effective diffusion constant,

[TABLE]

where $\boldsymbol{L}_{j}^{(n)}=\boldsymbol{C}_{j}\boldsymbol{\Gamma}_{j}^{-n}\boldsymbol{M}_{j}\boldsymbol{C}_{j}^{*}$ (for $j=3,4$).

Note that the first order effective damping constant $\boldsymbol{K}^{(1)}(\alpha_{1},\alpha_{2})=\int_{0}^{\infty}\boldsymbol{\kappa}(t)dt$ and the first order effective diffusion constant $\boldsymbol{L}^{(1)}(\alpha_{3},\alpha_{4})=\int_{0}^{\infty}\boldsymbol{R}(t)dt$ are simply the effective damping constant and effective diffusion constant introduced in [42]. The memory function and the covariance function of the noise process can be expressed in terms of these constants:

[TABLE]

Assumption 3.4**.**

The matrix $\boldsymbol{K}_{1}^{(1)}$ in the expression for first order effective damping constant is invertible and the matrix $\boldsymbol{K}_{2}^{(1)}$ equals zero. Similarly, in the expression for the first order effective diffusion constant $\boldsymbol{L}_{3}^{(1)}$, which is invertible, $\boldsymbol{L}_{4}^{(1)}=\boldsymbol{0}$.

In order to develop intuition about general GLEs, it will be helpful to study the following exactly solvable special case.

Example 2*.*

(An exactly solvable case) In the GLE (1.1), set $\boldsymbol{F}_{e}=\boldsymbol{0}$. Let $\boldsymbol{\gamma}_{0}(t,\boldsymbol{x})=\boldsymbol{\gamma}_{0}$, $\boldsymbol{\sigma}_{0}(t,\boldsymbol{x})=\boldsymbol{\sigma}_{0}$, $\boldsymbol{h}(t,\boldsymbol{x})=\boldsymbol{h}$, $\boldsymbol{g}(t,\boldsymbol{x})=\boldsymbol{g}$ and $\boldsymbol{\sigma}(t,\boldsymbol{x})=\boldsymbol{\sigma}$ be constant matrices. The initial data are the random variables, $\boldsymbol{x}(0)=\boldsymbol{x}$, $\boldsymbol{v}(0)=\boldsymbol{v}$, independent of $\{\boldsymbol{\xi}(t),t\in\mathbb{R}^{+}\}$ and of $\{\boldsymbol{W}^{(k)}(t),t\in\mathbb{R}^{+}\}$. The resulting GLE is:

[TABLE]

Of particular interest is the GLE (3.10) with $\boldsymbol{\gamma}_{0}=\boldsymbol{\sigma}_{0}\boldsymbol{\sigma}_{0}^{*}/2\geq 0$, $\boldsymbol{g}=\boldsymbol{h}^{*}=\boldsymbol{\sigma}>0$, and $\boldsymbol{R}(t)=\boldsymbol{\kappa}(t)=\boldsymbol{\kappa}^{*}(t)$, so that the fluctuation-dissipation relations hold (see Remark 1.1 and also Remark 3.6). The resulting GLE gives a simple model describing the motion of a free particle, interacting with a heat bath. Note that generally the process $\boldsymbol{v}(t)$ is not assumed to be stationary, in particular $\boldsymbol{v}(0)$ could be an arbitrarily distributed random variable.

The following proposition gives the asymptotic behavior of the spectral densities (equivalently, covariance functions or memory functions), the regularity222Sample path continuity does not in general imply mean-square continuity. (in the mean-square sense) of the noise process, and, in the exactly solvable case of Example 2, the long-time mean-squared displacement of the particle.

Proposition 3.5**.**

Suppose that the Assumptions 3.2 and 3.4 are satisfied. Let $\boldsymbol{x}(t)=\int_{0}^{t}\boldsymbol{v}(s)ds\in\mathbb{R}^{d}$, where $\boldsymbol{v}(t)$ solves the GLE (3.10).

• (i)

We have $\boldsymbol{\mathcal{S}}_{3}(\omega)=O(1)$ as $\omega\to 0$. Also, let $k\geq 3$ be a positive odd integer and assume that $\boldsymbol{L}_{4}^{(n)}=0$ for $0<n<k$, where $n$ is odd, and $\boldsymbol{L}_{4}^{(k)}\neq 0$. Then $\boldsymbol{\mathcal{S}}_{4}(\omega)=O(\omega^{k-1})$ as $\omega\to 0$. If there exists $h>0$ such that the noise spectral density, $\boldsymbol{\mathcal{S}}(\omega)=O\left(\frac{1}{\omega^{2h+1}}\right)$ as $\omega\to\infty$, then $\boldsymbol{\xi}_{t}$ is $n$-times mean-square differentiable333A process $X(t)$ is mean-square differentiable on a time interval $\mathcal{\tau}$ if for every $t\in\mathcal{\tau}$,

$\left\|\frac{X(t+h)-X(t)}{h}-\frac{dX}{dt}\right\|\to 0,$

as $h\to 0$. * for $n<h$.*

• (ii)

Let $\hat{\boldsymbol{\kappa}}(z)$ denote the Laplace transform of $\boldsymbol{\kappa}(t)$, i.e. $\hat{\boldsymbol{\kappa}}(z):=\int_{0}^{\infty}\boldsymbol{\kappa}(t)e^{-zt}dt$, and $\mathcal{E}=\frac{1}{2}mE[\boldsymbol{v}\boldsymbol{v}^{*}]$ be the particle’s initial average kinetic energy. Assume for simplicity that $\boldsymbol{R}(t)=\boldsymbol{\kappa}(t)$ and $\boldsymbol{\sigma}\boldsymbol{\kappa}(t)\boldsymbol{\sigma}^{*}=\boldsymbol{h}^{*}\boldsymbol{\kappa}^{*}(t)\boldsymbol{g}^{*}$. Then we have the following formula for the particle’s mean-squared displacement (MSD):

[TABLE]

where the Laplace transform of $\boldsymbol{H}(t)$ is given by $\hat{\boldsymbol{H}}(z)=z\hat{\boldsymbol{F}}(z)$, with

[TABLE]

For (iii) and (iv) below, we consider the process $\boldsymbol{x}_{t}$ solving the GLE (3.10) with $\boldsymbol{\gamma}_{0}=\boldsymbol{\sigma}_{0}\boldsymbol{\sigma}_{0}^{*}/2\geq 0$, $\boldsymbol{g}=\boldsymbol{h}^{*}=\boldsymbol{\sigma}>0$, and $\boldsymbol{R}(t)=\boldsymbol{\kappa}(t)=\boldsymbol{\kappa}^{*}(t)$.

• (iii)

Let $\alpha_{1}=\alpha_{3}=1$ ($\alpha_{i}$, for $i=2,4$, can be 0 or 1 and $\boldsymbol{F}_{0}$ can be zero or nonzero). Then $E[\boldsymbol{x}(t)\boldsymbol{x}^{*}(t)]=O(t)$ as $t\to\infty$, in which case we say that the particle diffuses normally.

• (iv)

Let $\alpha_{1}=0$, $\alpha_{2}=1$ and $\boldsymbol{F}_{0}=\boldsymbol{0}$ (the vanishing effective damping constant case). Then $E[\boldsymbol{x}(t)\boldsymbol{x}^{*}(t)]=O(t^{2})$ as $t\to\infty$, in which case we say that the particle exhibits a ballistic (super-diffusive) behavior.

Proof.

• (i)

For $i=3,4$, it is easy to compute that

[TABLE]

and so one has:

[TABLE]

as $\omega\to 0$. The first two statements in (i) then follow by Assumption 3.4. The last statement follows from Lemma 6.11 in [44].

• (ii)

Note that $\dot{\boldsymbol{x}}(t)=\boldsymbol{v}(t)$, with $\boldsymbol{x}(0)=\boldsymbol{0}$ and $\boldsymbol{v}(t)$ solving the GLE (3.10), rewritten as:

[TABLE]

where $\boldsymbol{\eta}(t)dt=d\boldsymbol{W}^{(k)}(t)$, and $\boldsymbol{v}_{0}=\boldsymbol{v}$ is a random variable that is independent of $\{\boldsymbol{\xi}(t),t\in\mathbb{R}^{+}\}$ and of $\{\boldsymbol{\eta}(t),t\in\mathbb{R}^{+}\}$. These equations can be solved analytically by means of Laplace transform. Applying Laplace transform on the equations for $\boldsymbol{x}_{t}$ and $\boldsymbol{v}_{t}$ gives:

[TABLE]

and thus

[TABLE]

where $\hat{\boldsymbol{H}}(z)=(mz^{2}\boldsymbol{I}+z\boldsymbol{\gamma}_{0}+z\boldsymbol{g}\hat{\boldsymbol{\kappa}}(z)\boldsymbol{h})^{-1}$. Taking the inverse transform gives the following formula for $\boldsymbol{x}(t)$:

[TABLE]

where $\boldsymbol{H}(0)=\boldsymbol{0}$.

Therefore, using the mutual independence of $\boldsymbol{v}$, $\{\boldsymbol{\xi}(t),t\in\mathbb{R}^{+}\}$ and $\{\boldsymbol{\eta}(t),t\in\mathbb{R}^{+}\}$, the Itô isometry, and the assumption that $\boldsymbol{R}(t)=\boldsymbol{\kappa}(t)$, we obtain:

[TABLE]

where

[TABLE]

To compute the double integral $\boldsymbol{L}(t)$, we first rewrite it as $\boldsymbol{L}(t)=\boldsymbol{L}_{1}(t)+\boldsymbol{L}_{2}(t)$, with

[TABLE]

We then compute:

[TABLE]

where $\star$ denotes convolution. Now note that, by the convolution theorem, $(\boldsymbol{\sigma}\boldsymbol{\kappa}\boldsymbol{\sigma}^{*}\star\boldsymbol{H}^{*})(u)$ is the inverse Laplace transform of $\boldsymbol{\sigma}\hat{\boldsymbol{\kappa}}(z)\boldsymbol{\sigma}^{*}\hat{\boldsymbol{H}^{*}}(z)$, which can be written as $\boldsymbol{I}/z-(mz\boldsymbol{I}+\boldsymbol{\gamma}_{0}^{*})\hat{\boldsymbol{H}^{*}}(z)$ by using the assumption that $\boldsymbol{\sigma}\boldsymbol{\kappa}(t)\boldsymbol{\sigma}^{*}=\boldsymbol{h}^{*}\boldsymbol{\kappa}^{*}(t)\boldsymbol{g}^{*}$. Computing the inverse transform gives us:

[TABLE]

Similarly, we obtain $\boldsymbol{L}_{2}(t)=\boldsymbol{L}_{1}(t)$, and so $\boldsymbol{L}(t)=2\boldsymbol{L}_{1}(t)$. Therefore, combining (3.22) and (3.30) gives us the desired formula for MSD.

• (iii)

$\&$ (iv) The assumptions that $\boldsymbol{g}=\boldsymbol{h}^{*}=\boldsymbol{\sigma}$ and $\boldsymbol{R}(t)=\boldsymbol{\kappa}(t)=\boldsymbol{\kappa}^{*}(t)$ ensure that we can apply the MSD formula in (ii). The additional assumption that $\boldsymbol{\gamma}_{0}=\boldsymbol{\sigma}_{0}\boldsymbol{\sigma}_{0}^{*}/2$ (fluctuation-dissipation relation of the first kind) implies that $\hat{\boldsymbol{H}}(z)=\hat{\boldsymbol{H}}^{*}(z)$ and simplifies the formula to:

[TABLE]

To determine the behavior of $E[\boldsymbol{x}(t)\boldsymbol{x}^{*}(t)]$ as $t\to\infty$, it suffices to investigate the asymptotic behavior of $\hat{\boldsymbol{H}}(z)$, whose formula is given in (ii), as $z\to 0$. Noting that

[TABLE]

and using Assumption 3.4, we find that, as $z\to 0$,

[TABLE]

Therefore, if $\boldsymbol{\gamma}_{0}=\boldsymbol{\sigma}_{0}\boldsymbol{\sigma}_{0}^{*}/2$ is non-zero, then $\hat{\boldsymbol{H}}(z)\sim 1/z$ as $z\to 0$. Otherwise, if in addition $\alpha_{1}=1$, then $\hat{\boldsymbol{H}}(z)\sim 1/z$ as $z\to 0$, whereas if in addition $\alpha_{1}=0$, $\alpha_{2}=1$, then $\hat{\boldsymbol{H}}(z)\sim 1/z^{2}$ as $z\to 0$. The results in (iii) and (iv) then follow by applying the Tauberian theorems [18], which say, in particular, that if $\hat{\boldsymbol{H}}(z)\sim 1/z^{\beta}$ as $z\to 0$, then $\boldsymbol{H}(t)\sim t^{\beta-1}$ as $t\to\infty$, for $\beta=1,2$ here.

∎

Remark 3.6*.*

We emphasize that superdiffusion with $E[\boldsymbol{x}(t)\boldsymbol{x}^{*}(t)]$ behaving as $t^{\alpha}$ as $t\to\infty$, where $\alpha>2$, cannot take place when the velocity process converges to a stationary state. For a system to behave this way, the velocity itself has to grow with time. Moreover, we remark that one could obtain a richer class of asymptotic behaviors for the MSD by relaxing the assumption of fluctuation-dissipation relations.

To summarize, (i) says that in the case where $\boldsymbol{F}_{0}=\boldsymbol{0}$, $\alpha_{1}=\alpha_{3}=0$, the $n$th order effective constants characterize the asymptotic behavior of the spectral densities at low frequencies; (ii) provides a formula for the particle’s mean-squared displacement, and (iii)-(iv) classify the types of diffusive behavior of the GLE model, in the exactly solvable case of Example 2, satisfying the fluctuation-dissipation relations. We emphasize that in the sequel we go beyond the above exactly solvable case; in particular the coefficients $\boldsymbol{g}$, $\boldsymbol{h}$, $\boldsymbol{\sigma}$, $\boldsymbol{\gamma}_{0}$, $\boldsymbol{\sigma}_{0}$ will depend in general on the particle’s position. However, the GLE in the exactly solvable case can be viewed as linear approximation to the general GLE (1.1) (by expanding these coefficients in a Taylor series about a fixed position $\boldsymbol{x}^{\prime}\in\mathbb{R}^{d}$).

In view of Proposition 3.5, the parameters $\alpha_{i}\in\{0,1\}$ allow us to control diffusive behavior of the generalized Langevin dynamics. Our GLE models are very general and need not satisfy a fluctuation-dissipation relation. As we will see, these different behaviors motivate our introduction and study of various homogenization schemes for the GLE. Depending on the physical systems under consideration, one scheme might be more realistic than the others. It is one of the goals of this paper to explore homogenization schemes for different GLE classes.

Remark 3.7*.*

In finite dimension, it is not possible to realize generalized Langevin dynamics with a noise and/or memory function whose spectral density varies as $1/\omega^{p}$, $p\in(0,1)$, near $\omega=0$ (i.e. the so-called $1/f$-type noise [36]), and, consequently, the noise covariance function and/or memory function decay as a power $1/t^{\alpha}$, $\alpha\in(0,1)$, as $t\to\infty$. In this case one can use the formula in (ii) of Proposition 3.5 to show, at least for the exactly solvable case in Example 2 where the fluctuation-dissipation relations hold, that the asymptotic behavior of the particle is sub-diffusive, i.e. $E[\boldsymbol{x}(t)\boldsymbol{x}^{*}(t)]=O(t^{\beta})$, where $\beta\in(0,1)$, as $t\to\infty$ (see also the related works [48, 15]). Sub-diffusive behavior has been discovered in a wide range of statistical and biological systems [34], and, therefore, making the study in this case relevant. One could, following the ideas in [21, 54], extend the state space of the GLEs to an infinite-dimensional one, in order to study the sub-diffusive case. Homogenization studies, where more technicalities are expected to be encountered due to the infinite-dimensional nature of the systems, for this case will be explored in a future work.

3.2. Generalized Langevin Systems as Input-Output Stochastic Dynamical Systems with Multiple Time Scales

In this subsection, we discuss GLEs of the form (1.1), under Assumptions 3.2-3.4, from the input-output system-theoretic and multiple time scale points of view.

First, we introduce the notion of stochastic dynamical systems.

Definition 3.8**.**

A stochastic dynamical system is a pair $(\boldsymbol{Z},\boldsymbol{\mathcal{F}})$ of vector-valued stochastic processes satisfying equations of the form:

[TABLE]

where $\boldsymbol{A}$, $\boldsymbol{B}$, $\boldsymbol{C}$ are measurable (jointly in $t$ and $\boldsymbol{Z}$) mappings, $\boldsymbol{\eta}(t)$ is a random process (the input). $\boldsymbol{Z}(t)$ is called the state process and $\boldsymbol{\mathcal{F}}(t)$ the output process (observation process). The system is linear if all the mappings are at most linear in $\boldsymbol{Z}$; otherwise the system is nonlinear. The system is time-invariant if all the mappings are independent of $t$.

The equation for the particle’s position, together with the GLE (1.1), can be cast as the system of SDEs for the Markov process

$\boldsymbol{z}_{t}:=(\boldsymbol{x}_{t},\boldsymbol{v}_{t},\boldsymbol{y}^{1}_{t},\boldsymbol{y}^{2}_{t},\boldsymbol{\beta}^{3}_{t},\boldsymbol{\beta}^{4}_{t})\in\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d_{1}}\times\mathbb{R}^{d_{2}}\times\mathbb{R}^{d_{3}}\times\mathbb{R}^{d_{4}}$:

[TABLE]

where we have defined the auxiliary memory processes:

[TABLE]

It is easy to see that the pairs $(\boldsymbol{\beta}^{i},\boldsymbol{\xi}^{i})$, $i=3,4$, defined in the previous subsection, are linear time-invariant Gaussian stochastic dynamical systems with a white noise input (and therefore the state processes $\boldsymbol{\beta}^{j}(t)$ are Markov) in the sense of Definition 3.8. Also, the pairs $(\boldsymbol{y}_{t}^{i},\boldsymbol{C}_{i}\boldsymbol{y}_{t}^{i})$ $(i=1,2$) are linear stochastic dynamical systems driven by the random processes $\boldsymbol{M}_{i}\boldsymbol{C}_{i}^{*}\boldsymbol{h}(t,\boldsymbol{x}_{t})\boldsymbol{v}_{t}$, which depend on the particle’s position and velocity variables. The generalized Langevin system can be viewed as a nonlinear stochastic dynamical system $(\boldsymbol{z},\boldsymbol{\mathcal{F}})$, where the components of $\boldsymbol{z}$ satisfy the SDEs (3.36)-(3.39) and $\boldsymbol{\mathcal{F}}$ is a measurable mapping describing an output process or a quantity of interest, for instance,

[TABLE]

for $p>0$ and $T>0$. In the exactly solvable case of Example 2, the generalized Langevin system reduces to a linear time-invariant stochastic dynamical system and can be viewed as a network of input-output systems consisting of components modeling the memory and noise. One of the goals of homogenization of GLEs is to reduce the number of the components needed to describe the effective dynamics in the considered limit.

It is natural question, what class of GLEs should be taken as the starting point for homogenization. For feasible treatment, the GLEs should be in some sense minimal. In the network interpretation, the original system should be completely described by a minimal number of components, with no redundancies. We will discuss this based on a time scale analysis in the following.

The (discrete) spectrum of the $\boldsymbol{\Gamma}_{i}$ ($i=1,2$) and of the $\boldsymbol{\Gamma}_{j}$ ($j=3,4$) (or equivalently, the spectrum of the Bohl memory function $\boldsymbol{\kappa}(t)$ and that of the covariance function $\boldsymbol{R}(t)$— see Definition 2.5 in [72]) encode information about the memory time scales and noise correlation time scales present in the generalized Langevin system respectively. In realistic experiments, there may be many, possibly infinitely many, time scales (each corresponding to a mode of the environment), but typically they cannot be all observed and/or controlled. When modeling a system, it is important to focus on those time scales that are controllable and observable. This motivates the following definition, closely related to the notions of controllable and observable eigenvalues from the systems theory [72].

Definition 3.9**.**

Consider a linear stochastic dynamical system $(\boldsymbol{Z},\boldsymbol{\mathcal{F}})$, as in Definition 3.8, where $\boldsymbol{A}\in\mathbb{R}^{n\times n}$, $\boldsymbol{B}\in\mathbb{R}^{n\times k}$, $\boldsymbol{C}\in\mathbb{R}^{m\times n}$ are constant matrices. The time scale, $\tau:=1/\lambda$, where $\lambda$ is an eigenvalue of $\boldsymbol{A}$, in the system, is called $(\boldsymbol{A},\boldsymbol{B})$-controllable (or simply controllable) if $rank[\boldsymbol{A}-\lambda\boldsymbol{I}\ \ \boldsymbol{B}]=n$ and $(\boldsymbol{C},\boldsymbol{A})$-observable (or simply observable) if $rank[\boldsymbol{A}-\lambda\boldsymbol{I}\ \ \boldsymbol{C}]^{*}=n$.

The following proposition, which follows from Theorem 3.13 in [72], states well-known results regarding the above notions.

Proposition 3.10**.**

Consider the linear dynamical system defined in Definition 3.9. Then

• (i)

the system is controllable (more precisely, $(\boldsymbol{A},\boldsymbol{B})$-controllable, i.e.

$[\boldsymbol{B}\ \ \boldsymbol{A}\boldsymbol{B}\ \ \cdots\ \ \boldsymbol{A}^{n-1}\boldsymbol{B}]$* is full rank) if and only if every time scale of the system is controllable.*

• (ii)

the system is observable (more precisely, $(\boldsymbol{C},\boldsymbol{A})$-observable, i.e.

$[\boldsymbol{C}\ \ \boldsymbol{C}\boldsymbol{A}\ \ \cdots\ \ \boldsymbol{C}\boldsymbol{A}^{n-1}]^{*}$* is full rank) if and only if every time scale of the system is observable.*

For $i=1,2,3,4$ we define the time scales, $\tau_{i,k_{i}}:=1/\lambda_{i,k_{i}}$, where $\lambda_{i,k_{i}}$ ($k_{i}=1,\dots,d_{i}$) are eigenvalues of $\boldsymbol{\Gamma}_{i}$. We refer to the $\tau_{1,k_{1}},\tau_{2,k_{2}}$ as memory time scales and the $\tau_{3,k_{3}},\tau_{4,k_{4}}$ as noise correlation time scales.

Our consideration of GLEs will be based on the following assumption.

Assumption 3.11**.**

All the memory time scales and the noise correlation time scales in the generalized Langevin systems described by (1.1) are controllable and observable.

From the mathematical point of view, our consideration minimizes the dimension of the state space on which the GLE is realized and therefore minimizes the complexity of the model which will be taken as the starting point for our homogenization studies. Indeed, recall that a stochastic realization is minimal if the realized process has no other stochastic realization of smaller dimension. It follows from our assumptions that all the realizations of the memory function and noise process are minimal, since a sufficient condition for a linear stochastic dynamical system to be minimal is that it is controllable (or reachable in the language of [43]), observable and the spectral factor of its spectral density is minimal [43].

4. On the Homogenization of Generalized Langevin Dynamics

In this section, we discuss some new directions for homogenization of GLEs.

In the case of non-vanishing (first order) effective damping constant and effective diffusion constant, homogenization of a version of the GLE (1.1) was studied in [42], where a limiting SDE for the position process was obtained in the limit, in which all the characteristic time scales of the system (i.e. the inertial time scale, the memory time scale and the noise correlation time scale) tend to zero at the same rate. Extending this result, we are going to focus on the following two cases.

• (A)

The case where an instantaneous damping term is present in the GLE, i.e. $\boldsymbol{F}_{0}\neq\boldsymbol{0}$, or the non-vanishing effective damping constant case, i.e. $\alpha_{1}=1$. Together with the conditions in Example 2, this gives a model for normally diffusing systems; see Proposition 3.5 (iii). One can study the limit in which the inertial time scale and a subset (possibly all or none of) of other characteristic time scales of the system tend to zero; in particular the small mass limit in the case $\boldsymbol{F}_{0}\neq\boldsymbol{0}$ of the generalized Langevin dynamics. We remark that the small mass limit is not well-defined in the case $\boldsymbol{F}_{0}=\boldsymbol{0}$ and $\alpha_{1}=\alpha_{3}=1$ – this was first observed in [49], where it was pointed out that the limit leads to the phenomenon of anomalous gap of the particle’s mean-squared displacement (see also [10, 29]).

• (B)

The vanishing effective damping constant and effective diffusion constant case, i.e. $\boldsymbol{F}_{0}=\boldsymbol{0}$, $\alpha_{1}=\alpha_{3}=0$, $\alpha_{2}=\alpha_{4}=1$. Together with the conditions in Example 2, this gives a model for systems with super-diffusive behavior; see Proposition 3.5 (iv). One can study the limit in which the inertial time scale, a subset of the memory time scales and a subset of the noise correlation time scales tend to zero at the same rate. Such effective models are physically relevant when they preserve the asymptotic behavior of the spectral densities at low and/or high frequencies in the limit. Situations are also possible, where some of the eigenmodes of the memory and noise spectrum are damped much stronger than other, for example due to an injection of monochromatic light from a laser into the system, which is originally in thermal equilibrium. This justifies studying homogenization limits that selectively target a part of frequencies of memory and noise.

We will study homogenization of the GLE (1.1) in the limits described in the above scenarios. In all cases, the inertial time scale is taken to zero – this gives rise to the singular nature of the limit problems. We remark that one could also consider the more interesting scenarios in which the time scales tend to zero at different rates, but we choose not to pursue this in this already lengthy paper.

Notation. Throughout the paper, we denote the variables in the pre-limit equations by small letters (for instance, $\boldsymbol{x}^{\epsilon}(t)$), and those of the limiting equations by capital letters (for instance, $\boldsymbol{X}(t)$). We use Einstein’s summation convention on repeated indices. The Euclidean norm of an arbitrary vector $\boldsymbol{w}$ is denoted by $|\boldsymbol{w}|$ and the (induced operator) norm of a matrix $\boldsymbol{A}$ by $\|\boldsymbol{A}\|$. For an $\mathbb{R}^{n_{2}\times n_{3}}$-valued function $\boldsymbol{f}(\boldsymbol{y}):=([f]_{jk}(\boldsymbol{y}))_{j=1,\dots,n_{2};k=1,\dots,n_{3}}$, $\boldsymbol{y}:=([y]_{1},\dots,[y]_{n_{1}})\in\mathbb{R}^{n_{1}}$, we denote by $(\boldsymbol{f})_{\boldsymbol{y}}(\boldsymbol{y})$ the $n_{1}n_{2}\times n_{3}$ matrix:

[TABLE]

where $\boldsymbol{\nabla}_{\boldsymbol{y}}[f]_{jk}(\boldsymbol{y})$ stands for the gradient vector $\left(\frac{\partial[f]_{jk}(\boldsymbol{y})}{\partial[y]_{1}},\dots,\frac{\partial[f]_{jk}(\boldsymbol{y})}{\partial[y]_{n_{1}}}\right)\in\mathbb{R}^{n_{1}}$ for every $j,k$. We denote by $\boldsymbol{\nabla}\cdot$ the divergence operator which contracts a matrix-valued function to a vector-valued function, i.e. for the matrix-valued function $\boldsymbol{A}(\boldsymbol{X})$, the $i$th component of its divergence is given by $(\boldsymbol{\nabla}\cdot\boldsymbol{A})^{i}=\sum_{j}\frac{\partial A^{ij}}{\partial X^{j}}$. Lastly, the symbol $\mathbb{E}$ denotes expectation with respect to the probability measure $\mathbb{P}$.

5. Small Mass Limit of Generalized Langevin Dynamics

Consider the following family of equations for the processes $(\boldsymbol{x}_{t}^{m},\boldsymbol{v}_{t}^{m})\in\mathbb{R}^{d\times d}$, $t\in[0,T]$, $m>0$:

[TABLE]

where $\boldsymbol{\kappa}(t)$ and $\boldsymbol{\xi}_{t}$ are the memory function and noise process defined in (3.1) and (3.2) respectively, with each of the $\alpha_{i}$ ($i=1,2,3,4$) equal to zero or to one. The equations (5.1)-(5.2) are equivalent to the following system of SDEs for the Markov process $\boldsymbol{z}^{m}_{t}:=(\boldsymbol{x}^{m}_{t},\boldsymbol{v}^{m}_{t},\boldsymbol{y}^{1,m}_{t},\boldsymbol{y}^{2,m}_{t},\boldsymbol{\beta}^{3,m}_{t},\boldsymbol{\beta}^{4,m}_{t})\in\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d_{1}}\times\mathbb{R}^{d_{2}}\times\mathbb{R}^{d_{3}}\times\mathbb{R}^{d_{4}}$:

[TABLE]

where we have defined the auxiliary memory processes:

[TABLE]

Note that the processes $\boldsymbol{\beta}_{t}^{3,m}$ and $\boldsymbol{\beta}_{t}^{4,m}$ do not actually depend on $m$, but we are adding the superscript $m$ for a more homogeneous notation.

We make the following simplifying assumptions concerning (5.3)-(5.6). Let $\boldsymbol{W}^{(q_{j})}$ ($j=3,4$) be independent Wiener processes on a filtered probability space $(\Omega,\mathcal{F},\mathcal{F}_{t},\mathbb{P})$ satisfying the usual conditions and let $\mathbb{E}$ denote expectation with respect to $\mathbb{P}$.

Assumption 5.1**.**

There are no explosions, i.e. almost surely, for every $m>0$ there exists global unique solution to the pre-limit SDE (5.3)-(5.6) and also to the limiting SDEs (5.8)-(5.10) on the time interval $[0,T]$.

Assumption 5.2**.**

For $t\in\mathbb{R}^{+}$, $\boldsymbol{y}\in\mathbb{R}^{d}$, the functions $\boldsymbol{F}_{e}(t,\boldsymbol{y})$, $\boldsymbol{\sigma}_{0}(t,\boldsymbol{y})$ and $\boldsymbol{\sigma}(t,\boldsymbol{y})$ are continuous and bounded (in $t$ and $\boldsymbol{y}$) as well as Lipschitz in $\boldsymbol{y}$, whereas the functions $\boldsymbol{\gamma}_{0}(t,\boldsymbol{y})$, $\boldsymbol{g}(t,\boldsymbol{y})$, $\boldsymbol{h}(t,\boldsymbol{y})$, $(\boldsymbol{\gamma}_{0})_{\boldsymbol{y}}(t,\boldsymbol{y})$, $(\boldsymbol{g})_{\boldsymbol{y}}(t,\boldsymbol{y})$ and $(\boldsymbol{h})_{\boldsymbol{y}}(t,\boldsymbol{y})$ are continuously differentiable and Lipschitz in $\boldsymbol{y}$ as well as bounded (in $t$ and $\boldsymbol{y}$). Moreover, the functions $(\boldsymbol{\gamma}_{0})_{\boldsymbol{y}\boldsymbol{y}}(t,\boldsymbol{y})$, $(\boldsymbol{g})_{\boldsymbol{y}\boldsymbol{y}}(t,\boldsymbol{y})$ and $(\boldsymbol{h})_{\boldsymbol{y}\boldsymbol{y}}(t,\boldsymbol{y})$ are bounded for every $t\in\mathbb{R}^{+}$, $\boldsymbol{y}\in\mathbb{R}^{d}$.

Assumption 5.3**.**

The initial data $\boldsymbol{x},\boldsymbol{v}\in\mathbb{R}^{d}$ are $\mathcal{F}_{0}$-measurable random variables independent of the $\sigma$-algebra generated by the Wiener processes $\boldsymbol{W}^{(q_{j})}$ ($j=3,4$). They are independent of $m$ and have finite moments of all orders.

The following theorem describes the homogenized behavior of the particle’s position modeled by the family of the equations (5.1)-(5.2)—or, equivalently, by the SDE systems (5.3)-(5.6)—in the limit as the particle’s mass tends to zero.

Theorem 5.4**.**

Let $\boldsymbol{z}_{t}^{m}:=(\boldsymbol{x}_{t}^{m},\boldsymbol{v}_{t}^{m},\boldsymbol{y}_{t}^{1,m},\boldsymbol{y}_{t}^{2,m},\boldsymbol{\beta}_{t}^{3,m},\boldsymbol{\beta}_{t}^{4,m})$ be a family of processes solving the SDE system (5.3)-(5.6). Suppose that Assumptions 3.2-3.11 and Assumptions 5.1-5.3 hold. In addition, suppose that for every $m>0$, $\boldsymbol{x}\in\mathbb{R}^{d}$, the family of matrices $\boldsymbol{\gamma}_{0}(t,\boldsymbol{x})$ is positive stable, uniformly in $t$ and $\boldsymbol{x}$. Then as $m\to 0$, the position process $\boldsymbol{x}^{m}_{t}$ converges to $\boldsymbol{X}_{t}$, where $\boldsymbol{X}_{t}$ is the first component of the process $(\boldsymbol{X}_{t},\boldsymbol{Y}_{t}^{1},\boldsymbol{Y}_{t}^{2},\boldsymbol{\beta}_{t}^{3},\boldsymbol{\beta}_{t}^{4})$ satisfying the Itô SDE system:

[TABLE]

where the $i$th component of the $\boldsymbol{S}^{(k)}$ ($k=0,1,2$) is given by:

[TABLE]

and for $k=1,2$,

[TABLE]

with $\boldsymbol{J}\in\mathbb{R}^{d\times d}$ solving the Lyapunov equation, $\boldsymbol{\gamma}_{0}\boldsymbol{J}+\boldsymbol{J}\boldsymbol{\gamma}_{0}^{*}=\boldsymbol{\sigma}_{0}\boldsymbol{\sigma}_{0}^{*}$. The convergence is obtained in the following sense: for all finite $T>0$, $\sup_{t\in[0,T]}|\boldsymbol{x}^{m}_{t}-\boldsymbol{X}_{t}|\to 0$ in probability, as $m\to 0$.

Proof.

We prove the theorem by applying Theorem A.6. Using the notation in the statement of Theorem A.6, let $\epsilon=m$, $n_{1}=d+d_{1}+d_{2}+d_{3}+d_{4}$, $n_{2}=d$, $k_{1}=q_{3}+q_{4}$, $k_{2}=k$, $\boldsymbol{x}^{\epsilon}(t)=(\boldsymbol{x}_{t}^{m},\boldsymbol{y}_{t}^{1,m},\boldsymbol{y}_{t}^{2,m},\boldsymbol{\beta}_{t}^{3,m},\boldsymbol{\beta}_{t}^{4,m})$, $\boldsymbol{v}^{\epsilon}(t)=\boldsymbol{v}_{t}^{m}$,

[TABLE]

$\boldsymbol{W}^{(k_{1})}(t)=(\boldsymbol{W}_{t}^{(q_{3})},\boldsymbol{W}_{t}^{(q_{4})})$ and $\boldsymbol{W}^{(k_{2})}(t)=\boldsymbol{W}_{t}^{(k)}$. The initial conditions are $\boldsymbol{x}(0)=(\boldsymbol{x},\boldsymbol{0},\boldsymbol{0},\boldsymbol{\beta}_{0}^{3},\boldsymbol{\beta}_{0}^{4})$ and $\boldsymbol{v}(0)=\boldsymbol{v}$, where $\boldsymbol{\beta}_{0}^{j}$ $(j=3,4$) are normally distributed with mean-zero and covariance $\boldsymbol{M}_{j}$. They are independent of $m$.

Observe that in the above formula, $\boldsymbol{a}_{i}$, $\boldsymbol{b}_{i}$, $\boldsymbol{\sigma}_{i}$ ($i=1,2$) do not depend explicitly on $\epsilon=m$, so by the convention adopted earlier, we denote them $\boldsymbol{A}_{i}$, $\boldsymbol{B}_{i}$, $\boldsymbol{\Sigma}_{i}$ respectively, and we put $a_{i}=b_{i}=c_{i}=d_{i}=\infty$, where $a_{i},b_{i},c_{i},d_{i}$ are the rates in Assumption A.5.

Next, we verify the assumptions of Theorem A.6. Assumption A.1 clearly follows from the Assumption 5.1. Since the family of matrices $\boldsymbol{\gamma}_{0}(t,\boldsymbol{x})$ is positive stable (uniformly in $t$ and $\boldsymbol{x}$), Assumption A.2 is satisfied. It is straightforward to see that our assumptions on the coefficients of the GLE imply Assumption A.3. As $\boldsymbol{x}(0)$ and $\boldsymbol{v}(0)$ are random variables independent of $m$, Assumption A.4 holds by our assumptions on the initial conditions $\boldsymbol{x}_{0}$, $\boldsymbol{v}_{0}$ and $\boldsymbol{\beta}^{j}_{0}$ ($j=3,4$). Finally, as noted earlier, Assumption A.5 holds with $a_{i}=b_{i}=c_{i}=d_{i}=\infty$. The assumptions of the Theorem A.6 are thus satisfied. Applying it, we obtain the limiting SDE system (5.8)-(5.10).

∎

We remark that the limiting SDE is unique up to transformation in (3.6), as pointed out already in [42].

Remark 5.5*.*

In the special case when $\alpha_{i}=0$ for $i=1,2,3,4$ and the coefficients do not depend on $t$ explicitly, Theorem 5.4 reduces to the result obtained in [28]. In general, by comparing the result with the one obtained in [28], we see that perturbing the original Markovian system by adding a memory and colored noise changes the behavior of the homogenized system obtained in the small mass limit. In particular,

• (i)

the limiting equation for the particle’s position not only contains a correction drift term ($\boldsymbol{S}^{(0)}$) – the noise-induced drift, but is also coupled to equations for other slow variables;

• (ii)

in the case when $\alpha_{1}$ and/or $\alpha_{2}$ equal $1$, the limiting equation for the (slow) auxiliary memory variables contains correction drift terms ($\boldsymbol{S}^{(1)}$ and/or $\boldsymbol{S}^{(2)}$) – which could be called the memory-induced drifts. Interestingly, the memory-induced drifts disappear when $\boldsymbol{h}$ is proportional to $\boldsymbol{\gamma}_{0}$, a phenomenon that can be attributed to the interaction between the forces $\boldsymbol{F}_{0}$ and $\boldsymbol{F}_{1}$.

Note that the highly coupled structure of the limiting SDEs is due to the fact that only one time scale (inertial time scale) was taken to zero in the limit. We expect the structure to simplify when all time scales present in the problem are taken to zero at the same rate.

6. Homogenization for the Case of Vanishing Effective Damping Constant and Effective Diffusion Constant

In this section we consider the GLE (1.1), with $\boldsymbol{F}_{0}=\boldsymbol{0}$, $\alpha_{1}=\alpha_{3}=0$, and $\alpha_{2}=\alpha_{4}=1$. We explore a class of homogenization schemes, aiming to:

• (P1)

reduce the complexity of the generalized Langevin dynamics in a way that the homogenized dynamics can be realized on a state space with minimal dimension and are described by minimal number of effective parameters;

• (P2)

retain non-trivial effects of the memory and the colored noise in the homogenized dynamics by matching the asymptotic behavior of the spectral density of the noise process and memory function in the original and the effective model.

Remark 6.1*.*

Generally, the larger the number of time scales (the eigenvalues of the $\boldsymbol{\Gamma}_{i}$) present in the system, the higher the dimension of the state space needed to realize the generalized Langevin system. On the other hand, in addition to $\boldsymbol{\Gamma}_{i}$, information on $\boldsymbol{C}_{i}$ and $\boldsymbol{M}_{i}$ is needed to determine the asymptotic behavior of the spectral densities (see Proposition 3.5(i)). In other words, although analysis based solely on time scales consideration may reduce the dimension of the model, it does not in general allow one to achieve the model matching in (P2). It is desirable to have homogenization schemes that achieve both goals of dimension reduction (P1) and matching of models (P2). Such a scheme is considered below.

The idea is to consider the limit when the inertial time scale, a proper subset of the memory time scales and a proper subset of the noise correlation time scales tend to zero at the same rate. The case of sending all the characteristic time scales to zero is excluded here as it is uninteresting when the effective damping and diffusion vanish in the limit.

Recall that the notions of controllability and observability are invariant under the trivial equivalence relation of type (3.6). Therefore, one can, without loss of generality, assume that the $\boldsymbol{\Gamma}_{i}$ $(i=1,2,3,4)$ are already in the Jordan normal form and work in Jordan basis. Such form will reveal the slow-fast time scale structure of the system and so give us a rubric to develop homogenization schemes.

Assumption 6.2**.**

Let $i=2,4$. All the $\boldsymbol{\Gamma}_{i}$ are of the following Jordan normal form:

[TABLE]

where $N_{i}<d_{i}$, $\boldsymbol{\Gamma}_{i,k}\in\mathbb{R}^{\nu(\lambda_{i,k})\times\nu(\lambda_{i,k})}$ ($k=1,\dots,N_{i}$) is the Jordan block associated with the (controllable and observable) eigenvalue $\lambda_{i,k}$ (or time scale $\tau_{i,k}=1/\lambda_{i,k}$) and corresponds to the invariant subspace $\mathcal{X}_{i,k}=Ker(\lambda_{i,k}\boldsymbol{I}-\boldsymbol{\Gamma}_{i,k})^{\nu(\lambda_{i,k})}$, where $\nu(\lambda_{i,k})$ is the index of $\lambda_{i,k}$, i.e. the size of the largest Jordan block corresponding to the eigenvalue $\lambda_{i,k}$. Let $1\leq M_{i}<N_{i}$ and the eigenvalues be ordered as $0<\lambda_{i,1}\leq\dots\leq\lambda_{i,M_{i}}<\lambda_{i,M_{i}+1}\leq\dots\leq\lambda_{i,N_{i}}$, so that we have the invariant subspace decomposition, $\mathbb{R}^{d_{i}}=\bigoplus_{j=1}^{N_{i}}\mathcal{X}_{i,j}$, with $d_{i}=\sum_{k=1}^{N_{i}}\nu(\lambda_{i,k})$.

Let $0<l_{i}<d_{i}$. The following procedure studies generalized Langevin dynamics whose spectral densities of the memory and the noise process have the asymptotic behavior, $\boldsymbol{\mathcal{S}}_{i}(\omega)\sim\omega^{2l_{i}}$ for small $\omega$, and $\boldsymbol{\mathcal{S}}_{i}(\omega)\sim 1/\omega^{2d_{i}}$ for large $\omega$, for $i=2,4$. We construct a homogenized version of the model in such a way that its memory and noise processes have spectral densities whose asymptotic behavior at low $\omega$ matches that of the original model (to achieve (P2)), while that at high $\omega$ it varies as $1/\omega^{2l_{i}}$ (to achieve (P1)).

Algorithm 6.3**.**

Procedure to study a class of homogenization problems.

• (1)

Let $\alpha_{1}=\alpha_{3}=0$, $\alpha_{2}=\alpha_{4}=1$ and $\boldsymbol{F}_{0}=\boldsymbol{0}$ in the GLE (1.1). Suppose that Assumption 6.2 holds and there exists $M_{i}$ such that $l_{i}=\sum_{k=1}^{M_{i}}\nu(\lambda_{i,k})$. Take this $M_{i}$.

• (2)

For $i=2,4$, set $m=m^{\prime}\epsilon$ and $\lambda_{i,k}=\lambda^{\prime}_{i,k}/\epsilon$, for $k=M_{i}+1,\dots,N_{i}$ (i.e. we scale the $(d_{2}-l_{2})$ smallest memory time scales and the $(d_{4}-l_{4})$ smallest noise correlation time scales with $\epsilon$), where $m^{\prime}$ and the $\lambda^{\prime}_{i,k}$ are positive constants.

• (3)

Select the $\boldsymbol{C}_{i}$, $\boldsymbol{M}_{i}$, $\boldsymbol{\Sigma}_{i}$ such that the $\boldsymbol{C}_{i}$ are constant matrices independent of the $\lambda_{i,k}$ ($k=1,\dots,N_{i}$), $\boldsymbol{C}_{i}\boldsymbol{\Gamma}_{i}^{-n_{i}}\boldsymbol{M}_{i}\boldsymbol{C}_{i}^{*}=\boldsymbol{0}$ for $0<n_{i}<2l_{i}$, $\boldsymbol{C}_{i}\boldsymbol{\Gamma}_{i}^{-(2l_{i}+1)}\boldsymbol{M}_{i}\boldsymbol{C}_{i}^{*}\neq\boldsymbol{0}$, and upon a suitable rescaling involving the mass, memory time scales and noise correlation time scales the resulting family of GLEs can be cast in the form of the SDEs (A.3)-(A.4). Note that the matrix entries of the $\boldsymbol{M}_{i}$ and/or $\boldsymbol{\Sigma}_{i}$ necessarily depend on the $\lambda_{i,k}$ due to the Lyapunov equations that relate them to the $\boldsymbol{\Gamma}_{i}$.

• (4)

Apply Theorem A.6 to study the limit $\epsilon\to 0$ and obtain the homogenized model, under appropriate assumptions on the coefficients and parameters in the GLEs.

We remark that while one has the above procedure to study homogenization schemes that achieve (P1) and (P2), the derivations and formulas for the limiting equations could become tedious and complicated as the $l_{i}$ and $d_{i}$ become large. To illustrate this, we consider a simple yet still sufficiently general instance of Algorithm 6.3 in the following.

Assumption 6.4**.**

The spectral densities, $\boldsymbol{\mathcal{S}}_{i}(\omega)=\boldsymbol{\Phi}_{i}(i\omega)\boldsymbol{\Phi}_{i}^{*}(-i\omega)$ ($i=2,4$), with the (minimal) spectral factor:

[TABLE]

where the $\boldsymbol{P}_{i}(z)\in\mathbb{R}^{p_{i}\times m_{i}}$ are matrix-valued monomials with degree $l_{i}$ :

[TABLE]

and the $\boldsymbol{Q}_{i}(z)\in\mathbb{R}^{p_{i}\times p_{i}}$ are matrix-valued polynomials with degree $d_{i}$, i.e.

[TABLE]

*Here $p_{2}=q$, $p_{4}=r$, the $m_{i}$ $(i=2,4$) are positive integers, the $\boldsymbol{B}_{l_{i}}\in\mathbb{R}^{p_{i}\times m_{i}}$ are constant matrices, $\boldsymbol{\Gamma}_{i,k}\in\mathbb{R}^{p_{i}\times p_{i}}$ are diagonal matrices with positive entries, and $\boldsymbol{I}$ denotes identity matrix of appropriate dimension. *

Under Assumption 6.4, the spectral densities have the following asymptotic behavior: $\boldsymbol{\mathcal{S}}_{i}(\omega)\sim\omega^{2l_{i}}$ for small $\omega$, and $\boldsymbol{\mathcal{S}}_{i}(\omega)\sim 1/\omega^{2d_{i}}$ for large $\omega$. One can then implement Algorithm 6.3 explicitly to study homogenization for a sufficiently large class of GLEs, where the rescaled spectral densities tend to the ones with the asymptotic behavior mentioned in the paragraph just before Algorithm 6.3 in the limit. We discuss one such implementation in Appendix C. Since the calculations become more complicated as $l_{i}$ and $d_{i}$ become large, we will only study simpler cases and illustrate how things could get complicated in the following.

We assume $d_{2}$ and $d_{4}$ are even integers and consider in detail the case when $l_{2}=l_{4}=l=1$, $d_{2}=d_{4}=h=2$,

[TABLE]

with $\lambda_{2,d_{2}}\geq\dots\geq\lambda_{2,d_{2}/2+1}>\lambda_{2,d_{2}/2}\geq\dots\geq\lambda_{2,1}>0$ and $\lambda_{4,d_{4}}\geq\dots\geq\lambda_{4,d_{4}/2+1}>\lambda_{4,d_{4}/2}\geq\dots\geq\lambda_{4,1}>0$ in Assumption 6.4, so that for $i=2,4$,

[TABLE]

We consider:

[TABLE]

where

[TABLE]

$p_{2}=q$ and $p_{4}=r$ as in Assumption 6.4. One can verify that this is indeed the vanishing effective damping constant and effective diffusion constant case (i.e. $\boldsymbol{C}_{i}\boldsymbol{\Gamma}_{i}^{-1}\boldsymbol{M}_{i}\boldsymbol{C}_{i}^{*}=\boldsymbol{0}$ for $i=2,4$). Also, for $i=2,4$, the memory kernel, $\boldsymbol{\kappa}_{2}(t)$ and covariance function, $\boldsymbol{R}_{4}(t)$, are of the following bi-exponential form:

[TABLE]

and their Fourier transforms are:

[TABLE]

which vary as $\omega^{2}$ near $\omega=0$. Note that in the above the $\boldsymbol{B}_{i}$ do not necessarily commute with the $\boldsymbol{\Gamma}_{i,j}$.

Following step (2) of Algorithm 6.3, we set $m=m_{0}\epsilon$, $\boldsymbol{\Gamma}_{i,2}=\boldsymbol{\gamma}_{i,2}/\epsilon$ for $i=2,4$, where $m_{0}>0$ is a constant and the $\boldsymbol{\gamma}_{i,2}$ are diagonal matrices with positive eigenvalues, in (6.16)-(6.17). We consider the family of GLEs (parametrized by $\epsilon>0$):

[TABLE]

where

[TABLE]

and the covariance function of the noise process $\boldsymbol{\xi}_{t}^{\epsilon}=\boldsymbol{C}_{4}\boldsymbol{\beta}_{t}^{4,\epsilon}$ is given by

[TABLE]

Note that $\boldsymbol{\kappa}_{2}^{\epsilon}(t)$ and $\boldsymbol{R}_{4}^{\epsilon}(t)$ converge (in the sense of distribution), as $\epsilon\to 0$, to

[TABLE]

with $i=2$ and $i=4$ respectively. The corresponding spectral densities are

[TABLE]

with $i=2$ and $i=4$ respectively.

Together with the equation for the particle’s position, the equations (6.18)-(6.19) form the SDE system:

[TABLE]

where

[TABLE]

In the following, we take $\epsilon\in\mathcal{E}$ to be small. We make the following assumptions, similar to those made in Theorem 5.4.

Assumption 6.5**.**

There are no explosions, i.e. almost surely, for every $\epsilon\in\mathcal{E}$, there exist unique solutions on the time interval $[0,T]$ to the pre-limit SDEs (6.24)-(6.29) and to the limiting SDEs (6.36).

Assumption 6.6**.**

The initial data $\boldsymbol{x},\boldsymbol{v}\in\mathbb{R}^{d}$ are $\mathcal{F}_{0}$-measurable random variables independent of the $\sigma$-algebra generated by the Wiener processes $\boldsymbol{W}^{(q_{j})}$ ($j=3,4$). They are independent of $\epsilon$ and have finite moments of all orders.

The following theorem describes the homogenized dynamics of the family of the GLEs (6.18)-(6.19) (or equivalently, of the SDEs (6.24)-(6.29)) in the limit $\epsilon\to 0$, i.e. when the inertial time scale, one half of the memory time scales and one half of the noise correlation time scales in the original generalized Langevin system tend to zero at the same rate.

Theorem 6.7**.**

Consider the family of the GLEs (6.18)-(6.19) (or equivalently, of the SDEs (6.24)-(6.29)). Suppose that Assumption 5.2 and Assumptions 6.4-6.6 hold, with the $\boldsymbol{C}_{i}$, $\boldsymbol{\Sigma}_{i}$, $\boldsymbol{M}_{i}$ and $\boldsymbol{\Gamma}_{i}$ ($i=2,4$) given in (6.7)-(6.9).

Assume that for every $t\in\mathbb{R}^{+}$, $\epsilon>0$, $\boldsymbol{x}\in\mathbb{R}^{d}$,

[TABLE]

are invertible for all $\lambda$ in the right half plane $\{\lambda\in\mathbb{C}:Re(\lambda)>0\}$, where

[TABLE]

Also, assume that $\boldsymbol{\nu}(t,\boldsymbol{x}):=\frac{1}{2}\boldsymbol{g}(t,\boldsymbol{x})\boldsymbol{B}_{2}\boldsymbol{B}_{2}^{*}\boldsymbol{h}(t,\boldsymbol{x})$ is invertible for every $t\in\mathbb{R}^{+}$, $\boldsymbol{x}\in\mathbb{R}^{d}$.

Then the particle’s position, $\boldsymbol{x}^{\epsilon}_{t}\in\mathbb{R}^{d}$, solving the family of GLEs, converges as $\epsilon\to 0$, to $\boldsymbol{X}_{t}\in\mathbb{R}^{d}$, where $\boldsymbol{X}_{t}$ is the first component of the process $\boldsymbol{\theta}_{t}:=(\boldsymbol{X}_{t},\boldsymbol{Y}_{t},\boldsymbol{Z}_{t})\in\mathbb{R}^{d+d_{2}/2+d_{4}/2}$, satisfying the Itô SDE:

[TABLE]

where

[TABLE]

[TABLE]

and the $i$th component of $\boldsymbol{Q}$, $i=1,\dots,d+d_{2}/2+d_{4}/2$, is given by:

[TABLE]

with $\boldsymbol{H}(t,\boldsymbol{X})=\boldsymbol{T}(t,\boldsymbol{X})\boldsymbol{U}^{-1}(t,\boldsymbol{X})\in\mathbb{R}^{(d+d_{2}/2+d_{4}/2)\times(d+d_{2}/2+d_{4}/2)}$ and

$\boldsymbol{J}\in\mathbb{R}^{(d+d_{2}/2+d_{4}/2)\times(d+d_{2}/2+d_{4}/2)}$* is the solution to the Lyapunov equation $\boldsymbol{U}\boldsymbol{J}+\boldsymbol{J}\boldsymbol{U}^{*}=diag(\boldsymbol{0},\boldsymbol{0},\boldsymbol{\gamma}_{4,2}^{2})$, where*

[TABLE]

The convergence holds in the same sense as in Theorem 5.4, i.e. for all finite $T>0$, $\sup_{t\in[0,T]}|\boldsymbol{x}^{\epsilon}_{t}-\boldsymbol{X}_{t}|\to 0$ in probability, as $\epsilon\to 0$.

Proof.

We apply Theorem A.6 to the SDEs (6.24)-(6.29). To this end, we set, in Theorem A.6, $n_{1}=n_{2}=d+d_{2}/2+d_{4}/2$, $k_{1}=k_{2}=d_{4}/2$ and

[TABLE]

[TABLE]

The initial conditions are $\boldsymbol{x}^{\epsilon}(0)=(\boldsymbol{x},\boldsymbol{0},\boldsymbol{\beta}_{0}^{4,1,\epsilon})$ and $\boldsymbol{v}^{\epsilon}(0)=(\boldsymbol{v},\boldsymbol{0},\boldsymbol{\beta}_{0}^{4,2,\epsilon})$; both depend on $\epsilon$.

We now verify each of the assumptions of Theorem A.6. Assumption A.1 clearly holds by our assumptions on the GLE. The assumptions on the coefficients in the SDEs follow easily from the Assumptions 5.2-5.3 and therefore Assumption A.3 holds.

Next, note that $\boldsymbol{\beta}_{0}^{4,\epsilon}=(\boldsymbol{\beta}_{0}^{4,1,\epsilon},\boldsymbol{\beta}_{0}^{4,2,\epsilon})$ is a random variable normally distributed with mean-zero and covariance:

[TABLE]

where

[TABLE]

as $\epsilon\to 0$. Using the bound $\mathbb{E}[|\boldsymbol{z}|^{p}]\leq C_{p}(\mathbb{E}[|\boldsymbol{z}|^{2}])^{p/2}$, where $\boldsymbol{z}$ is a mean-zero Gaussian random variable, $C_{p}>0$ is a constant and $p>0$, it is straightforward to see that Assumption A.4 is satisfied.

Note that $\boldsymbol{B}_{i}=\boldsymbol{b}_{i}$ (for $i=1,2$) by our convention, as the $\boldsymbol{b}_{i}$ do not depend explicitly on $\epsilon$. The uniform convergence of $\boldsymbol{a}_{i}(t,\boldsymbol{x},\epsilon)$, $(\boldsymbol{a}_{i})_{\boldsymbol{x}}(t,\boldsymbol{x},\epsilon)$ and $\boldsymbol{\sigma}_{i}(t,\boldsymbol{x},\epsilon)$ (in $\boldsymbol{x}$) to $\boldsymbol{A}_{i}(t,\boldsymbol{x})$, $(\boldsymbol{A}_{i})_{\boldsymbol{x}}(t,\boldsymbol{x})$ and $\boldsymbol{\Sigma}_{i}(t,\boldsymbol{x})$ respectively in the limit $\epsilon\to 0$ can be shown easily and, in fact, we see that $\boldsymbol{A}_{1}=\boldsymbol{T}$, $\boldsymbol{A}_{2}=-\boldsymbol{U}$, where $\boldsymbol{T}$ and $\boldsymbol{U}$ are given in the theorem,

[TABLE]

and $a_{1}=a_{2}=c_{1}=c_{2}=d_{1}=d_{2}=1$, $b_{1}=b_{2}=\infty$, where the $a_{i}$, $b_{i}$, $c_{i}$ and $d_{i}$ are from Assumption A.5 of Theorem A.6. Therefore, the first part of Assumption A.5 is satisfied.

It remains to verify the (uniform) Hurwitz stability of $\boldsymbol{a}_{2}$ and $\boldsymbol{A}_{2}$ (i.e. Assumption A.2 and the last part of Assumption A.5). This can be done using the methods of the proof of Theorem 2 in [42] and we omit the details here. The results then follow by applying Theorem A.6 and (6.36)-(6.40) follow from matrix algebraic calculations. ∎

It is clear from Theorem 6.7 that the homogenized position process is a component of the (slow) Markov process $\boldsymbol{\theta}_{t}$. In general, it is not a Markov process itself. Also, the components of $\boldsymbol{\theta}_{t}$ are coupled in a non-trivial way. We emphasize that one could use Theorem A.6 to study cases in which the different time scales are taken to zero in a different manner.

The limiting SDE for the position process may simplify under additional assumptions. In particular, in the one-dimensional case, i.e. with $d=1$ (or when all the matrix-valued coefficients and the parameters are diagonal in the multi-dimensional case), the formula for the limiting SDEs becomes more explicit. This special case has been studied in an earlier section in the context of the models (M1) and (M2) from Example 1.

7. Conclusions and Final Remarks

We have explored various homogenization schemes for a wide class of generalized Langevin equations. The relevance of the studied limit problems in the context of usual and anomalous diffusion of a particle in a heat bath. Our explorations here open up a wide range of possibilities and provide insights in the model reduction of and effective drifts in generalized Langevin systems.

The following summarizes the main conclusions of the paper:

• (i)

(stochastic modeling point of view) Homogenization schemes producing effective SDEs, driven by white noise, should be the exception rather than the rule. This is particularly important if one seeks to reduce the original model, retain its non-trivial features;

• (ii)

(complexity reduction point of view) There is a trade-off in simplifying GLE models with state-dependent coefficients: the greater the level of model reduction, the more complicated the correction drift terms, entering the homogenized model;

• (iii)

(statistical physics point of view) Homogenized equation obtained could be further simplified, i.e. number of effective equations could be reduced and the drift terms become simplified, when certain special conditions such as a fluctuation-dissipation theorem holds.

We conclude this paper by mentioning a very interesting future direction. As mentioned in Remark 3.7, one could extend the current GLE studies to the infinite-dimensional setting so that a larger class of memory functions and covariance functions can be covered. To this end, one can define the noise process as an appropriate linear functional of a Hilbert space valued process solving a stochastic evolution equation [12, 13]. This way, one can approach a class of GLEs, driven by noises having a completely monotone covariance function. This large class of functions contains covariances with power decay and thus the method outlined above can be viewed as an extension of those considered in [21, 54], where the memory function and covariance of the driving noise are represented as suitable infinite series with a power-law tail (these works are, to our knowledge, among the few works that study rigorously GLEs with a power-law memory). This approach to systems driven by strongly correlated noise, which is our future project, is expected to involve substantial technical difficulties. More importantly, one can expect that power decay of correlations leads to new phenomena, altering the nature of noise-induced drift.

Appendix A Homogenization for a Class of SDEs with State-Dependent Coefficients

In this section, we study homogenization for a general class of perturbed SDEs with state-dependent coefficients. Homogenization of differential equations has been extensively studied, from the seminal works of Kurtz [37], Papanicolaou [56] and Khasminksy [33] to the more recent works [58, 57, 28, 27, 5, 4, 9]. Here we are going to present yet another variant of homogenization result that will be needed for studying homogenization for our GLEs (see the last paragraph in Section 1.3 for comments on novelty of this result).

Let $n_{1}$, $n_{2}$, $k_{1}$, $k_{2}$ be positive integers. Let $\epsilon\in(0,\epsilon_{0}]=:\mathcal{E}$ be a small parameter and $\boldsymbol{x}^{\epsilon}(t)\in\mathbb{R}^{n_{1}}$, $\boldsymbol{v}^{\epsilon}(t)\in\mathbb{R}^{n_{2}}$ for $t\in[0,T]$, where $\epsilon_{0}>0$ and $T>0$ are finite constants. Let $\boldsymbol{W}^{(k_{1})}$ and $\boldsymbol{W}^{(k_{2})}$ denote independent Wiener processes, which are $\mathbb{R}^{k_{1}}$-valued and $\mathbb{R}^{k_{2}}$-valued respectively, on a filtered probability space $(\Omega,\mathcal{F},\mathcal{F}_{t},\mathbb{P})$ satisfying the usual conditions [31].

With respect to the standard bases of $\mathbb{R}^{n_{1}}$ and $\mathbb{R}^{n_{2}}$ respectively, we write:

[TABLE]

We consider the following family of perturbed SDE systems444Note that here the variables $\boldsymbol{x}^{\epsilon}(t)$ and $\boldsymbol{v}^{\epsilon}(t)$ are general and they do not necessarily represent position and velocity variables of a physical system. for

$(\boldsymbol{x}^{\epsilon}(t),\boldsymbol{v}^{\epsilon}(t))\in\mathbb{R}^{n_{1}+n_{2}}$:

[TABLE]

with the initial conditions, $\boldsymbol{x}^{\epsilon}(0)=\boldsymbol{x}^{\epsilon}$ and $\boldsymbol{v}^{\epsilon}(0)=\boldsymbol{v}^{\epsilon}$, where $\boldsymbol{x}^{\epsilon}$ and $\boldsymbol{v}^{\epsilon}$ are random variables that possibly depend on $\epsilon$. In the SDEs (A.3)-(A.4), the coefficients $\boldsymbol{a}_{1}:\mathbb{R}^{+}\times\mathbb{R}^{n_{1}}\times\mathcal{E}\to\mathbb{R}^{n_{1}\times n_{2}}$, $\boldsymbol{a}_{2}:\mathbb{R}^{+}\times\mathbb{R}^{n_{1}}\times\mathcal{E}\to\mathbb{R}^{n_{2}\times n_{2}}$, $\boldsymbol{\sigma}_{2}:\mathbb{R}^{+}\times\mathbb{R}^{n_{1}}\times\mathcal{E}\to\mathbb{R}^{n_{2}\times k_{2}}$ are non-zero matrix-valued functions, whereas $\boldsymbol{b}_{1}:\mathbb{R}^{+}\times\mathbb{R}^{n_{1}}\times\mathcal{E}\to\mathbb{R}^{n_{1}}$, $\boldsymbol{b}_{2}:\mathbb{R}^{+}\times\mathbb{R}^{n_{1}}\times\mathcal{E}\to\mathbb{R}^{n_{2}}$, $\boldsymbol{\sigma}_{1}:\mathbb{R}^{+}\times\mathbb{R}^{n_{1}}\times\mathcal{E}\to\mathbb{R}^{n_{1}\times k_{1}}$ are matrix-valued or vector-valued functions, which may depend on $\boldsymbol{x}^{\epsilon}$, as well as on $t$ and $\epsilon$ explicitly, as indicated by the parenthesis $(t,\boldsymbol{x}^{\epsilon}(t),\epsilon)$. In the case where the coefficients do not depend on $\epsilon$ explicitly, we will denote them by the corresponding capital letters (for instance, if $\boldsymbol{a}_{i}(t,\boldsymbol{x},\epsilon)=\boldsymbol{a}_{i}(t,\boldsymbol{x})$, then $\boldsymbol{a}_{i}(t,\boldsymbol{x}):=\boldsymbol{A}_{i}(t,\boldsymbol{x})$ etc.).

We are interested in the limit as $\epsilon\to 0$ of the SDEs (A.3)-(A.4), in particular the limiting behavior of the process $\boldsymbol{x}^{\epsilon}(t)$, under appropriate assumptions555We forewarn the readers that our assumptions can be relaxed in various directions (see later remarks) but we will not pursue these generalizations here. on the coefficients. In this section, we present a homogenization theorem that studies this limit and delay its proof and applications to later sections.

We make the following assumptions concerning the SDEs (A.3)-(A.4) and (A.10).

Assumption A.1**.**

The global solutions, defined on $[0,T]$, to the pre-limit SDEs (A.3)-(A.4) and to the limiting SDE (A.10) a.s. exist and are unique for all $\epsilon\in\mathcal{E}$ (i.e. there are no explosions).

Assumption A.2**.**

The matrix-valued functions

[TABLE]

are uniformly positive stable, i.e. all real parts of the eigenvalues of $-\boldsymbol{a}_{2}(t,\boldsymbol{y},\epsilon)$ are bounded from below, uniformly in $t$, $\boldsymbol{y}$ and $\epsilon$, by a positive constant (or, equivalently, the matrix-valued functions $\{\boldsymbol{a}_{2}(t,\boldsymbol{y},\epsilon);t\in[0,T],\boldsymbol{y}\in\mathbb{R}^{n_{1}},\epsilon\in\mathcal{E}\}$ are uniformly Hurwitz stable). They are $O(1)$ as $\epsilon\to 0$ (see Assumption A.5).

Assumption A.3**.**

For $t\in[0,T]$, $\boldsymbol{y}\in\mathbb{R}^{n_{1}}$, $\epsilon\in\mathcal{E}$, and $i=1,2$, the functions $\boldsymbol{b}_{i}(t,\boldsymbol{y},\epsilon)$ and $\boldsymbol{\sigma}_{i}(t,\boldsymbol{y},\epsilon)$ are continuous and bounded in $t$ and $\boldsymbol{y}$, and Lipschitz in $\boldsymbol{y}$, whereas the functions $\boldsymbol{a}_{i}(t,\boldsymbol{y},\epsilon)$ and $(\boldsymbol{a}_{i})_{\boldsymbol{y}}(t,\boldsymbol{y},\epsilon)$ are continuous in $t$, continuously differentiable in $\boldsymbol{y}$, bounded in $t$ and $\boldsymbol{y}$, and Lipschitz in $\boldsymbol{y}$. Moreover, the functions $(\boldsymbol{a}_{i})_{\boldsymbol{y}\boldsymbol{y}}(t,\boldsymbol{y},\epsilon)$ ($i=1,2$) are bounded for every $t\in[0,T]$, $\boldsymbol{y}\in\mathbb{R}^{n_{1}}$ and $\epsilon\in\mathcal{E}$.

We assume that the (global) Lipschitz constants are bounded by $L(\epsilon)$, where $L(\epsilon)=O(1)$ as $\epsilon\to 0$, i.e. for every $t\in[0,T]$, $\boldsymbol{x}$, $\boldsymbol{y}\in\mathbb{R}^{n_{1}}$,

[TABLE]

Assumption A.4**.**

The initial condition $\boldsymbol{x}^{\epsilon}_{0}=\boldsymbol{x}^{\epsilon}\in\mathbb{R}^{n_{1}}$ is an $\mathcal{F}_{0}$-measurable random variable that may depend on $\epsilon$, and we assume that $\mathbb{E}[|\boldsymbol{x}^{\epsilon}|^{p}]=O(1)$ as $\epsilon\to 0$ for all $p>0$. Also, $\boldsymbol{x}^{\epsilon}$ converges, in the limit as $\epsilon\to 0$, to a random variable $\boldsymbol{x}$ as follows: $\mathbb{E}\left[|\boldsymbol{x}^{\epsilon}-\boldsymbol{x}|^{p}\right]=O(\epsilon^{pr_{0}})$, where $r_{0}>1/2$ is a constant, as $\epsilon\to 0$. The initial condition $\boldsymbol{v}^{\epsilon}_{0}=\boldsymbol{v}^{\epsilon}\in\mathbb{R}^{n_{2}}$ is an $\mathcal{F}_{0}$-measurable random variable that may depend on $\epsilon$, and we assume that for every $p>0$, $\mathbb{E}[|\epsilon\boldsymbol{v}^{\epsilon}|^{p}]=O(\epsilon^{\alpha})$ as $\epsilon\to 0$, for some $\alpha\geq p/2$.

Assumption A.5**.**

For $i=1,2$, $t\in[0,T]$, and every $\boldsymbol{x}\in\mathbb{R}^{n_{1}}$, each of the matrix or vector entries of the (non-zero) functions $\boldsymbol{a}_{i}(t,\boldsymbol{x},\epsilon)$, $(\boldsymbol{a}_{i})_{\boldsymbol{x}}(t,\boldsymbol{x},\epsilon)$, $\boldsymbol{b}_{i}(t,\boldsymbol{x},\epsilon)$ and $\boldsymbol{\sigma}_{i}(t,\boldsymbol{x},\epsilon)$, converges, uniformly in $\boldsymbol{x}$, to a unique non-zero element, in the limit as $\epsilon\to 0$. Their limits are denoted by $\boldsymbol{A}_{i}(t,\boldsymbol{x})$, $(\boldsymbol{A}_{i})_{\boldsymbol{x}}(t,\boldsymbol{x})$, $\boldsymbol{B}_{i}(t,\boldsymbol{x})$ and $\boldsymbol{\Sigma}_{i}(t,\boldsymbol{x})$ respectively. Their rate of convergence is assumed to satisfy the following power law bounds: for every $t\in[0,T]$, $\boldsymbol{x}\in\mathbb{R}^{n_{1}}$ and $i=1,2$,

[TABLE]

*where $\alpha_{i}(\epsilon)=O(\epsilon^{a_{i}})$, $\beta_{i}(\epsilon)=O(\epsilon^{b_{i}})$, $\gamma_{i}(\epsilon)=O(\epsilon^{c_{i}})$ and $\theta_{i}(\epsilon)=O(\epsilon^{d_{i}})$, as $\epsilon\to 0$, for some positive exponents $a_{i}$, $b_{i}$, $c_{i}$ and $d_{i}$. Moreover, we assume that $\boldsymbol{A}_{2}(t,\boldsymbol{x})$ is Hurwitz stable for every $t$ and $\boldsymbol{x}$.

Convention.* In the case where the coefficients do not show explicit dependence on $\epsilon$ or the case when any of the coefficients $\boldsymbol{b}_{1}$, $\boldsymbol{b}_{2}$ and $\boldsymbol{\sigma}_{1}$ is zero, we set the exponent, describing the corresponding rate of convergence, to infinity. For instance, if $\boldsymbol{a}_{i}(t,\boldsymbol{x},\epsilon)=\boldsymbol{A}_{i}(t,\boldsymbol{x})$, we set $a_{i}=\infty$. Meanwhile, if $\boldsymbol{\sigma}_{1}=\boldsymbol{0}$, we set $c_{1}=\infty$, etc..*

We now state our homogenization theorem.

Theorem A.6**.**

Suppose that the family of SDE systems $\eqref{sde1}$-$\eqref{sde2}$ satisfies Assumption A.1-A.5. Let $(\boldsymbol{x}^{\epsilon}(t),\boldsymbol{v}^{\epsilon}(t))\in\mathbb{R}^{n_{1}}\times\mathbb{R}^{n_{2}}$ be their solutions, with the initial conditions $(\boldsymbol{x}^{\epsilon},\boldsymbol{v}^{\epsilon})$. Let $\boldsymbol{X}(t)\in\mathbb{R}^{n_{1}}$ be the solution to the following Itô SDE with the initial position $\boldsymbol{X}(0)=\boldsymbol{x}$:

[TABLE]

where $\boldsymbol{S}(t,\boldsymbol{X}(t))$ is the noise-induced drift vector whose $i$th component is given by

[TABLE]

where $i,l=1,\dots,n_{1},\ j,k=1,\dots,n_{2}$, or in index-free notation,

[TABLE]

and $\boldsymbol{J}\in\mathbb{R}^{n_{2}\times n_{2}}$ is the unique solution to the Lyapunov equation:

[TABLE]

Then the process $\boldsymbol{x}^{\epsilon}(t)$ converges, as $\epsilon\to 0$, to the solution $\boldsymbol{X}(t)$, of the Itô SDE (A.10), in the following sense: for all finite $T>0$, $p>0$, there exists a positive random variable $\epsilon_{1}$ such that

[TABLE]

in the limit as $\epsilon\to 0$, with $r>0$ is the rate determined to be:

[TABLE]

where the $a_{i}$, $b_{i}$, $c_{i}$, $d_{i}$ ($i=1,2$) are the positive constants from Assumption A.5. In particular, for all finite $T>0$,

[TABLE]

in probability, in the limit as $\epsilon\to 0$.

Remark A.7*.*

With more work and additional assumptions, one could prove the statements in Assumption A.1 from Assumption A.2-A.5. However, we choose to incorporate such existence and uniquess results into our assumptions and work with the assumptions as stated above. Moreover, as we have forewarned the readers, our assumptions can be relaxed in various directions at the cost of more technicalities. For instance, the boundedness assumption on the coefficients of the SDEs may be removed to obtain still a pathwise convergence result by adapting the techniques in [27] – see also analogous remarks in Remark 5 in [42]. However, we choose not to pursue the above technical details in this already lengthy paper.

Appendix B Proof of Theorem A.6

Proof of Theorem A.6 uses techniques developed in earlier works [28, 6, 42], but here one needs to additionally take into account the $\epsilon$-dependence of the coefficients in the SDEs (A.3)-(A.4). As a preparation for the proof, we need a few lemmas and propositions.

We start from an elementary calculus result.

Lemma B.1**.**

For $i=1,\dots,N$, let $\boldsymbol{f}_{i}(\boldsymbol{y},\epsilon):\mathbb{R}^{n}\times(0,\infty)\to\mathbb{R}^{m_{i}\times n}$ be bounded and globally Lipschitz in $\boldsymbol{y}$ for every $\epsilon>0$, with a Lipschitz constant that is bounded as $\epsilon\to 0$, i.e. for every $\boldsymbol{y},\boldsymbol{z}\in\mathbb{R}^{n}$, there exists a constant $M_{i}(\epsilon)>0$ such that

[TABLE]

where $M_{i}(\epsilon)=O(1)$ as $\epsilon\to 0$.

• (i)

Suppose that for each $i$ and $\boldsymbol{y}\in\mathbb{R}^{n}$, there exists a unique bounded $\boldsymbol{F}_{i}(\boldsymbol{y}):\mathbb{R}^{n}\to\mathbb{R}^{m_{i}\times n}$ and a constant $C_{i}>0$ such that $\|\boldsymbol{f}_{i}(\boldsymbol{y},\epsilon)-\boldsymbol{F}_{i}(\boldsymbol{y})\|\leq C_{i}\epsilon^{r_{i}}$, for some positive constant $r_{i}$, as $\epsilon\to 0$ (i.e. the left-hand side is of order $O(\epsilon^{r_{i}})$ as $\epsilon\to 0$). Then there exists constants $D$, $K_{1},\dots,K_{N}>0$, such that

[TABLE]

as $\epsilon\to 0$. If, in addition, $n=m_{1}$, $\boldsymbol{f}_{1}(\boldsymbol{y},\epsilon)$ and $\boldsymbol{F}_{1}(\boldsymbol{y})$ are invertible for every $\boldsymbol{y}\in\mathbb{R}^{n}$ and $\epsilon>0$, then $\|\boldsymbol{f}_{1}^{-1}(\boldsymbol{y},\epsilon)-\boldsymbol{F}_{1}^{-1}(\boldsymbol{y})\|=O(\epsilon^{r_{1}})$ as $\epsilon\to 0$.

• (ii)

Let $c_{i}\in\mathbb{R}$, $i=1,\dots,N$. For every $\epsilon>0$ and $\boldsymbol{y}\in\mathbb{R}^{n}$, $\sum_{i=1}^{N}c_{i}\boldsymbol{f}_{i}(\boldsymbol{y},\epsilon)$ and $\prod_{i=1}^{N}c_{i}\boldsymbol{f}_{i}(\boldsymbol{y},\epsilon)$ are globally Lipschitz with a Lipschitz constant that is $O(1)$ as $\epsilon\to 0$. Moreover, if $m_{1}=n$ and for every $\epsilon>0$, $\boldsymbol{y}\in\mathbb{R}^{n}$, $\boldsymbol{f}_{1}(\boldsymbol{y},\epsilon)$ is invertible, then for every $\epsilon>0$, $\boldsymbol{y}\in\mathbb{R}^{n}$, $\boldsymbol{f}^{-1}_{1}(\boldsymbol{y},\epsilon)$ is globally Lipschitz in $\boldsymbol{y}$ with a Lipschitz constant that is $O(1)$ as $\epsilon\to 0$.

Proof.

• (i)

We prove this inductively. The base case of $N=1$ clearly holds with $D=C_{1}$. Let $k\in\{1,\dots,N-1\}$. Assume that (B.2) holds with $N:=k$ and $D:=D_{k}$. Then

[TABLE]

as $\epsilon\to 0$, where $C$, $D_{k+1}$ are positive constants and we have used the inductive hypothesis and assumptions of the lemma in the last two lines above. The last statement follows from:

[TABLE]

as $\epsilon\to 0$, where $C$ is a positive constant.

• (ii)

The statements can be proven using the same techniques used for (i) and so we omit the proof.

∎

Let $\boldsymbol{x}^{\epsilon}(t)\in\mathbb{R}^{n_{1}}$, $\boldsymbol{v}^{\epsilon}(t)\in\mathbb{R}^{n_{2}}$ and $T>0$. For $t\in[0,T]$, let $\boldsymbol{p}^{\epsilon}(t):=\epsilon\boldsymbol{v}^{\epsilon}(t)$ denote a solution of the SDE:

[TABLE]

We provide estimates for the moments concerning the process $\boldsymbol{p}^{\epsilon}(t)$, under appropriate assumptions on the coefficients and the initial conditions, in the limit as $\epsilon\to 0$.

We need the following lemma, adapted from Proposition A.2.3 of [30], to obtain an exponential bound on certain fundamental matrix solution.

Lemma B.2**.**

Fix a filtered probability space $(\Omega,\mathcal{F},\mathcal{F}_{t},\mathbb{P})$. For each $\epsilon>0$, let $\boldsymbol{B}^{\epsilon}:[0,T]\times\Omega\to\mathbb{R}^{n\times n}$ be a bounded (uniformly in $\epsilon$, $\omega\in\Omega$ and $t\in[0,T]$), pathwise continuous process. Assume that the real parts of all eigenvalues of $\boldsymbol{B}$ are bounded from above by $-2\kappa$, uniformly in $\epsilon$, $\omega\in\Omega$ and $t\in[0,T]$, where $\kappa$ is a positive constant. Let $\boldsymbol{\Phi}^{\epsilon}(t,s,\omega)$ be the fundamental matrix that solves the initial value problem (IVP):

[TABLE]

Then there exists a constant $C>0$ and an (in general random666See also Remark 14 in [42].) $\epsilon_{1}=\epsilon_{1}(\omega)$ such that

[TABLE]

for all $\epsilon\leq\epsilon_{1}$ and for all $s,t\in[0,T]$.

Proof.

Let $u\in[s,t]$. We rewrite for $\omega\in\Omega$, $s,t\in[0,T]$:

[TABLE]

and represent the solution to the IVP as:

[TABLE]

Denote $\boldsymbol{W}^{\epsilon}(t,s,\omega):=e^{\kappa(t-s)/\epsilon}\boldsymbol{\Phi}^{\epsilon}(t,s,\omega)$. Setting $u=t$ in the above representation and multiplying both sides by $e^{\kappa(t-s)/\epsilon}$, we obtain:

[TABLE]

Since $\boldsymbol{B}^{\epsilon}$ is bounded (uniformly in $\omega$, $t$ and $\epsilon$), by assumption on the spectrum of $\boldsymbol{B}^{\epsilon}$, there exists a constant $C>0$, such that for all $s,t\in[0,T]$ we have

[TABLE]

Using this, we obtain:

[TABLE]

This leads to the estimate:

[TABLE]

where

[TABLE]

For a fixed $\omega\in\Omega$, $A_{\epsilon}(\omega)$ can be made arbitrary small as $\epsilon\to 0$. Therefore, there exists an $\epsilon_{1}=\epsilon_{1}(\omega)>0$ (generally dependent on $\omega$) such that

[TABLE]

for all $\epsilon\leq\epsilon_{1}$. This implies that $\sup_{s,t\in[0,T]}\|\boldsymbol{W}^{\epsilon}(t,s,\omega)\|\leq 2C$, which is the claimed bound. ∎

We now prove a lemma that gives a bound on a class of stochastic integrals. It is modification of Lemma 5.1 in [4]. In both cases, the main idea is to rewrite some of the stochastic integrals in terms of ordinary ones.

Lemma B.3**.**

Let $\boldsymbol{H}_{t}:=\boldsymbol{H}_{0}+\boldsymbol{M}_{t}+\boldsymbol{A}_{t}$ be the Doob-Meyer decomposition of a continuous $\mathbb{R}^{k}$-valued semimartingale on $(\Omega,\mathcal{F},\mathcal{F}_{t},P)$ with a local martingale $\boldsymbol{M}_{t}$ and a process of locally bounded variation $\boldsymbol{A}_{t}$. Let $\boldsymbol{V}\in L_{loc}^{1}(A)\cap L_{loc}^{2}(M)$ be $\mathbb{R}^{n\times k}$-valued and let $\boldsymbol{B}^{\epsilon}(t)$ be an adapted process whose values are $n\times n$ matrices, satisfying the assumptions of Lemma B.2. Let $\boldsymbol{\Phi}^{\epsilon}(t):=\boldsymbol{\Phi}^{\epsilon}(t,0)$ be the adapted $C^{1}$ process that pathwise solves the IVP (B.12). Then for every $T\geq\delta>0$ and for every $\epsilon\leq\epsilon_{1}$, we have the $\mathbb{P}$-a.s bound:

[TABLE]

where $N=\max\{k\in\mathbb{Z}:k\delta<T\}$, $\epsilon_{1}$, $\kappa$ and $C$ are from Lemma B.2, and $l_{2}$-norm is used on every $\mathbb{R}^{k}$.

Proof.

The proof is identical to that of Lemma 5.1 in [4] up to line (5.10), with the constant $\alpha$ there replaced by $\kappa$, etc. We let $\epsilon\leq\epsilon_{1}$ and replace the bound in line (5.11) there by the following bound, which follows from the semigroup property of the fundamental matrix process and Lemma B.2:

[TABLE]

Then we proceed as in the proof of Lemma 5.1 in [4] to get the desired bound. ∎

In particular, (B.13) and (B.23) hold for $\boldsymbol{B}^{\epsilon}=\boldsymbol{a}_{2}(t,\boldsymbol{x}^{\epsilon}(t),\epsilon)$.

Proposition B.4**.**

Suppose that Assumptions A.1-A.5 hold. For all $p\geq 1$, $T>0$, $0<\beta<p/2$, there exists a positive random variable $\epsilon_{1}$ such that:

[TABLE]

as $\epsilon\to 0$, where $\boldsymbol{p}^{\epsilon}(t)$ solves the SDE (B.11). Therefore, for any $p\geq 1$, $T>0$, $\beta>0$, we have

[TABLE]

as $\epsilon\to 0$, where $\|\cdot\|_{F}$ denotes the Frobenius norm.

Proof.

Let $\boldsymbol{\Phi}_{\epsilon}(t)$ be the matrix-valued process solving the IVP:

[TABLE]

Then,

[TABLE]

Therefore, for $T>0$ and $p\geq 1$, using the bound

[TABLE]

for $p\geq 1$ (here the $a_{i}\in\mathbb{R}$ and $N$ is a positive integer), taking supremum on both sides, and applying Lemma B.2 (with $\boldsymbol{B}^{\epsilon}=\boldsymbol{a}_{2}(t,\boldsymbol{x}^{\epsilon}(t),\epsilon)$), we estimate:

[TABLE]

for $\epsilon\leq\epsilon_{1}$, where $C>0$, $\kappa>0$, and $\epsilon_{1}>0$ is the random variable whose existence was proven in Lemma B.2.

Note that $\sup_{s\in[0,T]}|\boldsymbol{B}_{2}(s,\boldsymbol{x}^{\epsilon}(s))|^{p}<\infty$ and Assumption A.5 implies that

[TABLE]

where $\beta_{2}(\epsilon)\leq K\epsilon^{b_{2}}$.

Denote $\mathbb{E}_{1}[\cdot]=\mathbb{E}[\cdot;\epsilon\leq\epsilon_{1}]$, i.e. the expectation is taken on $\{\omega:\epsilon\leq\epsilon_{1}(\omega)\}$. We are going to estimate $\mathbb{E}_{1}\left[\sup_{t\in[0,T]}|\boldsymbol{p}^{\epsilon}(t)|^{p}\right]$.

By Assumption A.4, we have $\mathbb{E}_{1}[\sup_{t\in[0,T]}|\epsilon\boldsymbol{v}^{\epsilon}|^{p}]=O(\epsilon^{\alpha})$ as $\epsilon\to 0$, for some $\alpha\geq p/2$. Therefore, combining the above estimates, we obtain:

[TABLE]

where $C_{1}(p),C_{2}(p)>0$ are constants.

Next, the idea is to use Lemma B.3 and the Burkholder-Davis-Gundy inequality (see Theorem 3.28 in [31]) to estimate the last term on the right hand side above. This is analogous to the technique used in the proof of Proposition 5.1 in [4].

Let $\delta$ be a constant such that $0<\delta<T$. Applying Lemma B.3, we estimate, using (B.30):

[TABLE]

where $\|\boldsymbol{a}_{2}(t,\boldsymbol{x}^{\epsilon}(t),\epsilon)\|_{\infty}:=\sup_{t\in[0,T],\boldsymbol{y}\in\mathbb{R}^{n_{1}},\epsilon\in\mathcal{E}}\|\boldsymbol{a}_{2}(t,\boldsymbol{y},\epsilon)\|$ and

[TABLE]

We estimate:

[TABLE]

with $N:=\max\{k\in\mathbb{Z}:k\delta<T\}$, where we have used the fact that the $l^{\infty}$-norm on $\mathbb{R}^{N}$ is bounded by the $l^{q}$ norm for every $q\geq 1$ and then applied Hölder’s inequality to get the last two lines above.

Now, letting $\delta=\epsilon^{1-h}$ for $0<h<1$, and using the Burkholder-Davis-Gundy inequality,

[TABLE]

where $C_{p,q}$ is some constant and

[TABLE]

Since $N\delta<T$, we have $N\delta^{pq/2}<T\delta^{pq/2-1}=T\epsilon^{(1-h)(pq/2-1)}$. Therefore, $\mathbb{E}_{1}[\Pi]=O(\epsilon^{(1-h)(p/2-1/q)})$. For all $0<\beta<p/2$, one can choose $0<h<1$ and $q>1$ such that $(1-h)(p/2-1/q)=\beta$.

Therefore, we have

[TABLE]

as $\epsilon\to 0$, for all $0<\beta<p/2$.

Combining all the estimates obtained, one has:

[TABLE]

where the $C_{i}$ are positive constants, $\alpha\geq p/2$ is some constant, and $b_{2}>0$ is the constant from Assumption A.5. The statement of the proposition follows.

∎

We also need the following estimate on a class of integrals with respect to products of the coordinates of the process $\boldsymbol{p}^{\epsilon}(t)$.

Proposition B.5**.**

Suppose that Assumptions A.1-A.5 hold and $\epsilon\in\mathcal{E}$. Let $h^{\epsilon}:\mathbb{R}^{+}\times\mathbb{R}^{n_{1}}\to\mathbb{R}$ be a family of functions, continuously differentiable in $\boldsymbol{y}\in\mathbb{R}^{n_{1}}$ and bounded (in $s\in\mathbb{R}^{+}$ and $\boldsymbol{y}\in\mathbb{R}^{n_{1}}$), with bounded first derivatives $\boldsymbol{\nabla}_{\boldsymbol{y}}h^{\epsilon}(\boldsymbol{y})$ for $\boldsymbol{y}\in\mathbb{R}^{n_{1}}$. Assume that $h^{\epsilon}$ and $\boldsymbol{\nabla}_{\boldsymbol{y}}h^{\epsilon}(\boldsymbol{y})$ are $O(1)$ as $\epsilon\to 0$. Moreover, assume that $\frac{\partial}{\partial s}h^{\epsilon}$ is bounded (in all variables) and is $O(1)$ as $\epsilon\to 0$.

Then for any $p\geq 1$, $T>0$, $0<\beta<p/2$, $i,j=1,\dots,n_{2}$, in the limit as $\epsilon\to 0$ we have

[TABLE]

where $\boldsymbol{x}^{\epsilon}(t)$ and $\boldsymbol{p}^{\epsilon}(t)$ solve the SDEs (A.3)-(A.4) and the SDE (B.11) respectively, and $\epsilon_{1}$ is from Proposition B.4.

Proof.

Let $\epsilon\in\mathcal{E}$, $t\in[0,T]$, and $i,j=1,\dots,n_{2}$. An integration by parts gives:

[TABLE]

Using the notation $\mathbb{E}_{1}[\cdot]=\mathbb{E}[\cdot;\epsilon\leq\epsilon_{1}]$, we estimate, for $p\geq 1$,

[TABLE]

where $C(p,T)>0$ is a constant, $\|g^{\epsilon}\|_{\infty}:=\sup_{s\in[0,T],\boldsymbol{y}\in\mathbb{R}^{n_{1}}}|g^{\epsilon}(s,\boldsymbol{y})|$, and we have used Einstein’s summation over repeated indices convention.

Now, estimating as before, we obtain:

[TABLE]

where $D(p,T)>0$ is a constant.

By our assumptions, all the quantities in the form $\|\cdot\|_{\infty}$ are bounded and are $O(1)$ as $\epsilon\to 0$. Therefore, collecting the above estimates, using Assumption A.4, and applying Proposition B.4, we have, for $p\geq 1$, $T>0$, $i,j=1,\dots,n_{2}$,

[TABLE]

for every $0<\beta<p/2$.

∎

Now we proceed to prove Theorem A.6. Using the above moment estimates and the proof techniques in [4, 6], we are going to first obtain the convergence of $\boldsymbol{x}^{\epsilon}_{t}$ to $\boldsymbol{X}_{t}$ in the limit as $\epsilon\to 0$ in the following sense: for all finite $T>0$, $p\geq 1$,

[TABLE]

as $\epsilon\to 0$, where the $\epsilon_{1}$ is from Proposition B.4. The main tools are well known ordinary and stochastic integral inequalities, as well as a Gronwall type argument. This result would then imply that for all finite $T>0$, $\sup_{t\in[0,T]}|\boldsymbol{x}_{t}^{\epsilon}-\boldsymbol{X}_{t}|\to 0$ in probability, in the limit as $\epsilon\to 0$ (see Lemma 1 in [42]).

Proof.

(Proof of Theorem A.6) Let $T>0$ and recall that $[\boldsymbol{B}]_{i,j}$ denotes the $(i,j)$-entry of a matrix $\boldsymbol{B}$. First, we assume that $p>2$.

From (A.4), we have, for every $\epsilon>0$, $t\in[0,T]$,

[TABLE]

Substituting this into (A.3), we obtain:

[TABLE]

In integral form, we have:

[TABLE]

Its $i$th component, $[\boldsymbol{x}^{\epsilon}]_{i}(t)$ ($i=1,2,\dots,n_{1}$) is (recall that we are employing Einstein’s summation convention):

[TABLE]

Next, we perform integration by parts in the second term on the right hand side above:

[TABLE]

Substituting the following expression for $d[\boldsymbol{x}^{\epsilon}]_{l}(s)$:

[TABLE]

into (B.58), we obtain:

[TABLE]

Next, we apply Itô formula to $\epsilon\boldsymbol{v}^{\epsilon}(t)(\epsilon\boldsymbol{v}^{\epsilon}(t))^{*}\in\mathbb{R}^{n_{2}\times n_{2}}$:

[TABLE]

Denoting $\boldsymbol{J}^{\epsilon}(t):=\epsilon\boldsymbol{v}^{\epsilon}(t)(\boldsymbol{v}^{\epsilon}(t))^{*}$, we can rewrite the above as:

[TABLE]

where

[TABLE]

Since $-\boldsymbol{a}_{2}(t,\boldsymbol{x}^{\epsilon}(t),\epsilon)$ is positive stable uniformly (in $t$, $\boldsymbol{x}^{\epsilon}$ and $\epsilon$) by Assumption A.2, the solution of the Lyapunov equation (B.63) can be represented as:

[TABLE]

where

[TABLE]

for $n=1,2,3$.

Therefore, for $s\in[0,T]$,

[TABLE]

On the other hand, by (A.10),

[TABLE]

We use again the notation $\mathbb{E}_{1}[\cdot]:=\mathbb{E}[\cdot;\epsilon\leq\epsilon_{1}]$, where $\epsilon_{1}>0$ is the random variable from Proposition B.4.

For any $p>2$, $T>0$, $i=1,\dots,n_{1}$ (recall that $[\boldsymbol{b}]_{i}$ denotes the $i$th component of vector $\boldsymbol{b}$), we estimate:

[TABLE]

By Assumption A.4, $R_{0}=\mathbb{E}_{1}\left[|\boldsymbol{x}^{\epsilon}-\boldsymbol{x}|^{p}\right]\leq\mathbb{E}\left[|\boldsymbol{x}^{\epsilon}-\boldsymbol{x}|^{p}\right]=O(\epsilon^{pr_{0}})$ as $\epsilon\to 0$, where $r_{0}>1/2$ is a constant. We now estimate each of the $R_{k}$, $k=1,\dots,5$.

We have:

[TABLE]

on the set $S_{1}:=\{\epsilon:\epsilon\leq\epsilon_{1}\}$, where $\mathbb{1}_{A}$ denotes the indicator function of a set $A$, $L_{3}(\epsilon,p,T)=O(1)$ as $\epsilon\to 0$ and $C_{3}(p,T)$ is a constant dependent on $p$ and $T$. In the last two lines of the above estimate, we have used Assumption A.3, Assumption A.5, and the inequality:

[TABLE]

where $\boldsymbol{u}(s)\in\mathbb{R}^{n_{1}}$ for $s\in[0,T]$ (recall that $L(\epsilon)=O(1)$ as $\epsilon\to 0$ by the Assumption A.3).

Using again the above techniques, together with Lemma B.1, one obtains:

[TABLE]

on $S_{1}$, where $\alpha_{1}(\epsilon)$, $\alpha_{2}(\epsilon)$, $\beta_{2}(\epsilon)$ are from Assumption A.3, $L_{2}(\epsilon,p,T)=O(1)$ as $\epsilon\to 0$ and $C_{2}(p,T)$ is a constant.

To estimate $R_{5}$, we use the Burkholder-Davis-Gundy inequality:

[TABLE]

where $C^{\prime}_{p}$ is a positive constant and $\|\cdot\|_{F}$ denotes the Frobenius norm. Using Hölder’s inequality, Assumption A.3, Assumption A.5, and the above techniques, we obtain:

[TABLE]

on the set $S_{1}$, where $C_{p}^{\prime\prime}$ and $C_{p}^{\prime\prime\prime}$ are constants, $\gamma_{1}(\epsilon)$ is from Assumption A.3, $L_{5}(\epsilon,p,T)=O(1)$ as $\epsilon\to 0$ and $C_{5}(p,T)$ is a constant.

Similarly, using the above techniques and Lemma B.1, one can show:

[TABLE]

on $S_{1}$, where $\gamma_{2}(\epsilon)$ is from Assumption A.3, $L_{4}(\epsilon,p,T)=O(1)$ as $\epsilon\to 0$ and $C_{4}(p,T)$ is a constant.

To obtain a bound for $R_{1}$, first we estimate:

[TABLE]

and so $R_{1}\leq 6^{p-1}\sum_{k=0}^{6}\left(\mathbb{E}_{1}\sup_{t\in[0,T]}|\Pi_{k}|^{p}\right)=:6^{p-1}\sum_{k=0}^{6}M_{k}.$

It is straightforward to show, using the boundedness assumptions of the theorem, that for $k=0,1,2,3,5$:

[TABLE]

where the $C_{k}$ are positive constants.

Applying Proposition B.5, we obtain:

[TABLE]

on $S_{1}$, for all $0<\beta<p/2$, as $\epsilon\to 0$, where $C_{4}(p,T)$ is a positive constant.

We now estimate $M_{6}$:

[TABLE]

where the $C(p)$ are constants may vary from one expression to another.

Note that in the above, $\boldsymbol{J}_{3}^{\epsilon}(s)$ and $\boldsymbol{J}(s)$ are solutions to the Lyapunov equation

[TABLE]

and

[TABLE]

respectively.

Let $\boldsymbol{H}^{\epsilon}(s):=\boldsymbol{J}_{3}^{\epsilon}(s)-\boldsymbol{J}(s)$ and $\boldsymbol{G}^{\epsilon}(s):=\boldsymbol{a}_{2}(s,\boldsymbol{x}^{\epsilon}(s),\epsilon)-\boldsymbol{A}_{2}(s,\boldsymbol{X}(s))$. After some algebraic manipulations with the above pair of Lyapunov equations, we obtain another Lyapunov equation:

[TABLE]

By the last statement in Assumption A.5, $\boldsymbol{A}_{2}$ is positive stable uniformly (in $\boldsymbol{X}$ and $s$), therefore the above Lyapunov equation has a unique solution:

[TABLE]

Using (B.94), the assumptions of the theorem, and estimating as before, we obtain:

[TABLE]

on the set $S_{1}$, where $C(\epsilon,p,T)=O(1)$ as $\epsilon\to 0$ and $D(p,T)$ is a positive constant, $\alpha_{2}(\epsilon)$ and $\gamma_{2}(\epsilon)$ are from Assumption A.5.

Applying the above estimates, Lemma B.1 and techniques used earlier, one obtains from (B.90):

[TABLE]

on $S_{1}$, where $L_{6}(\epsilon,p,T)=O(1)$ as $\epsilon\to 0$, $C_{6}(p,T)$ is a positive constant, and $\alpha_{i}(\epsilon)$, $\theta_{i}(\epsilon)$ ($i=1,2$) and $\gamma_{2}(\epsilon)$ are from Assumption A.5.

Collecting the above estimates for the $M_{k}$, we obtain:

[TABLE]

on $S_{1}$, where $C_{1}(p,T)$ and $C_{3}(p,T)$ are constants, $C_{2}(\epsilon,p,T)=O(1)$ as $\epsilon\to 0$, and $M_{4}$ satisfies the bound in (B.87).

Using all the estimates for the $R_{i}$, we have:

[TABLE]

on $S_{1}$, where $L(\epsilon,p,T,n_{1})=O(1)$ as $\epsilon\to 0$, $r$ is the rate of convergence (A.15) in the statement of the theorem, $C(p,T,n_{1})$ is a constant that changes from line to line, and we have applied Proposition B.4, Lemma B.1 and Assumption A.5 to get the last expression in the above estimate.

Finally, applying the Gronwall lemma gives:

[TABLE]

on $S_{1}$.

(A.14) then follows for the case $p>2$. The result for $0<p\leq 2$ follows by an application of the Hölder’s inequality: for $0<p\leq 2$, taking $q>2$ so that $p/q<1$, we have

[TABLE]

for all $0<\beta<p^{\prime}$, as $\epsilon\to 0$. The last statement on convergence in probabiity in the theorem follows from Lemma 1 in [42]. ∎

Appendix C An Implementation of Algorithm 6.3 under Assumption 6.4

We describe how Algorithm 6.3 can be applied to a large class of GLEs, satisfying Assumption 6.4. For $i=2,4$, one can write

[TABLE]

where the $\boldsymbol{a}_{i,k}$ are related to the $\boldsymbol{\Gamma}_{i,k}$ as follows:

[TABLE]

Then it can be shown that $\boldsymbol{\Phi}_{i}(z)$ admits the following (controllable) realization [8]: $\boldsymbol{\Phi}_{i}(z)=\boldsymbol{H}_{i}(z\boldsymbol{I}+\boldsymbol{F}_{i})^{-1}\boldsymbol{G}_{i}$, with

[TABLE]

where $\boldsymbol{B}_{l_{i}}$ is in the $l_{i}$th slot,

[TABLE]

[TABLE]

Then the realization of the memory function (for the case $i=2$) and noise process (for the case $i=4$) can be obtained by taking $\boldsymbol{\Gamma}_{i}=\boldsymbol{F}_{i}$, $\boldsymbol{C}_{i}=\boldsymbol{H}_{i}$ and solving the following linear matrix inequality:

[TABLE]

for $\boldsymbol{M}_{i}=\boldsymbol{M}_{i}^{*}$ [73].

The above realization gives us the desired spectral densities. Indeed, let us use the transformation of type (3.6) to diagonalize the $\boldsymbol{M}_{i}$, i.e. $\boldsymbol{M}_{i}^{\prime}=\boldsymbol{T}_{i}\boldsymbol{M}_{i}\boldsymbol{T}_{i}^{*}=\boldsymbol{I}$, $\boldsymbol{\Gamma}_{i}^{\prime}=\boldsymbol{T}_{i}\boldsymbol{\Gamma}_{i}\boldsymbol{T}_{i}^{-1}$, $\boldsymbol{\Sigma}_{i}=\boldsymbol{T}_{i}\boldsymbol{\Sigma}_{i}$, $\boldsymbol{C}^{\prime}_{i}=\boldsymbol{C}_{i}\boldsymbol{T}_{i}^{-1}$. In this case, for $i=4$ we have: $(\boldsymbol{\xi}^{i})^{\prime}_{t}=\boldsymbol{C}_{i}^{\prime}(\boldsymbol{\beta}^{i})^{\prime}_{t}=\boldsymbol{C}_{i}\boldsymbol{\beta}^{i}_{t}=\boldsymbol{\xi}^{i}_{t}$, where $(\boldsymbol{\beta}^{i})^{\prime}_{t}$ solves the SDE:

[TABLE]

and one can compute the spectral density to be:

[TABLE]

A similar discussion applies to the realization of the memory function.

For $i=2,4$, set $m=\epsilon m_{0}$, $\boldsymbol{\Gamma}_{i,k}=\boldsymbol{\gamma}_{i,k}/\epsilon$ for $k=l_{i}+1,\dots,d_{i}$ and rescale the $\boldsymbol{B}_{l_{i}}$ with $\epsilon$ accordingly, so that the limit as $\epsilon\to 0$ of the rescaled spectral densities gives us the desired asymptotic behavior. The choice of which and how many of the $\boldsymbol{\Gamma}_{i,k}$ to rescale as well as the smallness of $\epsilon$ (i.e. what determines the wide separation of time scales and their magnitude) depends on the physical system under study. The resulting family of GLEs can then be cast in a form suitable for application of Theorem A.6 and the homogenized SDE for the particle’s position can be obtained, under appropriate assumptions on the coefficients of the GLE.

Appendix D Another Example for Section 2

Consider the case of $l_{2}=l_{4}=l=2$, $d_{2}=d_{4}=d=3$ in Assumption 6.4 and specialize to one-dimensional models as before. In this case the covariance function has a stronger singularity near $t=0$ than in cases studied previously. The spectral density of the driving noise in the GLE is taken to be:

[TABLE]

in which case the memory kernel (and covariance function) are

[TABLE]

where $0<\Gamma_{1}<\Gamma_{2}<\Gamma_{3}$ (see Figure 1 for a plot of $\kappa(t)$). This gives a model for hyper-diffusion of a particle in a heat bath [66].

Rescale the parameters by setting $m=m_{0}\epsilon$, $\Gamma_{3}=\gamma_{3}/\epsilon$, where $m_{0}$ and $\gamma_{3}$ are positive constants, and study the limit $\epsilon\to 0$ of the resulting family of GLEs as before. The resulting rescaled versions of $\kappa(t)$ and $\mathcal{S}(\omega)$ have the following asymptotic behavior as $\epsilon\to 0$:

[TABLE]

We outline, omitting details, a convergence result similar to Corollary 2.2, focusing on the particular case $g=h$. In this case, the particle’s position, $x^{\epsilon}_{t}\in\mathbb{R}$, converges, as $\epsilon\to 0$, to $X_{t}$, satisfying the following Itô SDE system:

[TABLE]

Inspecting the above SDEs in detail, one can make similar remarks to those made in Section 2. In particular, taking $g$ to be proportional to $\sigma$ (so that a fluctuation-dissipation relation holds) again allows us to reduce the number of effective SDEs. Physically, this means that homogenized GLEs for models of hyper-diffusion of particles in a non-equilibrium bath may be highly non-trivial but they simplify when the fluctuation-dissipation relation is satisfied.

Acknowledgment

S.H.Lim and J.Wehr were partially supported by the NSF grant DMS 1615045. S.H.Lim is grateful for the support provided by the Michael Tabor Fellowship from the Program in Applied Mathematics at the University of Arizona during the academic year 2017-2018. M.L. acknowledges the Spanish Ministry MINECO (National Plan 15 Grant: FISICATEAMO No. FIS2016-79508-P, SEVERO OCHOA No. SEV-2015-0522, FPI), European Social Fund, Fundació Cellex, Generalitat de Catalunya (AGAUR Grant No. 2017 SGR 1341 and CERCA/Program), ERC AdG OSYRIS, EU FETPRO QUIC, and the National Science Centre, Poland-Symfonia Grant No. 2016/20/W/ST4/00314.