Homogenization for a Class of Generalized Langevin Equations with an Application to Thermophoresis

TL;DR

This paper derives a homogenized model for generalized Langevin equations with state-dependent coefficients, revealing additional drift effects, and applies it to analyze thermophoresis in non-equilibrium heat baths.

Contribution

It extends existing homogenization results to a broader class of Langevin equations with weaker spectral assumptions and applies the theory to thermophoresis phenomena.

Findings

Homogenized process includes noise-induced drift terms.

Convergence proven under weaker spectral assumptions.

Application to thermophoresis in non-equilibrium heat baths.

Abstract

We study a class of systems whose dynamics are described by generalized Langevin equations with state-dependent coefficients. We find that in the limit, in which all the characteristic time scales vanish at the same rate, the position variable of the system converges to a homogenized process, described by an equation containing additional drift terms induced by the noise. The convergence results are obtained using the main result in \cite{hottovy2015smoluchowski}, whose version is proven here under a weaker spectral assumption on the damping matrix. We apply our results to study thermophoresis of a Brownian particle in a non-equilibrium heat bath.

Equations

$$m{\ddot{\mathbf{x}}}_t=\mathbf{F}({\mathbf{x}}_t)-\mathbf{g}({\mathbf{x}}_t){\int }_0^t\mathbf{\kappa }(t-s)\mathbf{h}({\mathbf{x}}_s){\dot{\mathbf{x}}}_{s}ds+\mathbf{\sigma }({\mathbf{x}}_t){\mathbf{\xi }}_t,$$

$${\mathbf{x}}_0=\mathbf{x},{\dot{\mathbf{x}}}_0=\mathbf{v}.$$

$$\mathbf{\kappa }(t)={\mathbf{C}}_1{e}^{-{\mathbf{\Gamma }}_1|t|}{\mathbf{M}}_1{\mathbf{C}}_1^{\ast }.$$

$${\mathbf{\xi }}_t={\mathbf{C}}_2{\mathbf{\beta }}_t,$$

$$d{\mathbf{\beta }}_t=-{\mathbf{\Gamma }}_2{\mathbf{\beta }}_{t}dt+{\mathbf{\Sigma }}_{2}d{\mathbf{W}}_t^{(q_2)},$$

$${\mathbf{\Gamma }}_i{\mathbf{M}}_i+{\mathbf{M}}_i{\mathbf{\Gamma }}_i^{\ast }={\mathbf{\Sigma }}_i{\mathbf{\Sigma }}_i^{\ast }.$$

$$\mathbf{R}(t)={\mathbf{C}}_2{e}^{-{\mathbf{\Gamma }}_2|t|}{\mathbf{M}}_2{\mathbf{C}}_2^{\ast },$$

$$({\mathbf{\Gamma }}_i^{\prime }={\mathbf{T}}_i{\mathbf{\Gamma }}_i{\mathbf{T}}_i^{-1},{\mathbf{M}}_i^{\prime }={\mathbf{T}}_i{\mathbf{M}}_i{\mathbf{T}}_i^{\ast },{\mathbf{C}}_i^{\prime }={\mathbf{C}}_i{\mathbf{T}}_i^{-1}),$$

$${\mathbf{K}}_1:={\int }_0^{\infty }\mathbf{\kappa }(t)dt={\mathbf{C}}_1{\mathbf{\Gamma }}_1^{-1}{\mathbf{M}}_1{\mathbf{C}}_1^{\ast }\in {\mathbb{R}}^{q\times q},$$

$${\mathbf{K}}_2:={\int }_0^{\infty }\mathbf{R}(t)dt={\mathbf{C}}_2{\mathbf{\Gamma }}_2^{-1}{\mathbf{M}}_2{\mathbf{C}}_2^{\ast }\in {\mathbb{R}}^{r\times r},$$

$$m{\ddot{\mathbf{x}}}_t=\mathbf{F}({\mathbf{x}}_t)-\mathbf{g}({\mathbf{x}}_t){\int }_0^t{\mathbf{C}}_1{e}^{-{\mathbf{\Gamma }}_1(t-s)}{\mathbf{M}}_1{\mathbf{C}}_1^{\ast }\mathbf{h}({\mathbf{x}}_s){\dot{\mathbf{x}}}_{s}ds+\mathbf{\sigma }({\mathbf{x}}_t){\mathbf{C}}_2{\mathbf{\beta }}_t,$$

$$\mathbf{\kappa }(t-s)={\mathbf{\kappa }}_1(t-s):=\mathbf{A}{e}^{-\mathbf{A}|t-s|},$$

$$d{\mathbf{\eta }}_t=-\mathbf{A}{\mathbf{\eta }}_{t}dt+\mathbf{A}d{\mathbf{W}}_t^{(d)}.$$

$$E[{\mathbf{\eta }}_t]=0,E[{\mathbf{\eta }}_t{\mathbf{\eta }}_s^{\ast }]=\frac{1}{2}{\mathbf{\kappa }}_1(t-s),s,t\ge 0.$$

$$m{\ddot{\mathbf{x}}}_t=\mathbf{F}({\mathbf{x}}_t)-\mathbf{g}({\mathbf{x}}_t){\int }_0^t{\mathbf{\kappa }}_1(t-s)\mathbf{h}({\mathbf{x}}_s){\dot{\mathbf{x}}}_{s}ds+\mathbf{\sigma }({\mathbf{x}}_t){\mathbf{\eta }}_t.$$

$$({\mathbf{\kappa }}_2{)}_{ii}(t-s):=\{\begin{array}{ll}\frac{1}{{\tau }_{ii}}{e}^{-{\omega }_{ii}^2\frac{|t-s|}{2{\tau }_{ii}}}\left[\cos \left(\frac{{\omega }_{ii}^0}{{\tau }_{ii}}(t-s)\right)+\frac{{\omega }_{ii}^1}{2}\sin \left(\frac{{\omega }_{ii}^0}{{\tau }_{ii}}|t-s|\right)\right],& \text{if }|{\omega }_{ii}|<2\\ \frac{1}{{\tau }_{ii}}{e}^{-{\omega }_{ii}^2\frac{|t-s|}{2{\tau }_{ii}}}\left[\cosh \left(\frac{{\tilde{\omega }}_{ii}^0}{{\tau }_{ii}}(t-s)\right)+\frac{{\tilde{\omega }}_{ii}^1}{2}\sinh \left(\frac{{\tilde{\omega }}_{ii}^0}{{\tau }_{ii}}|t-s|\right)\right],& \text{if }|{\omega }_{ii}|>2,\end{array}$$

$$\mathbf{\tau }d{\mathbf{h}}_t$$

$$\mathbf{\tau }d{\mathbf{u}}_t$$

$$E[{\mathbf{h}}_t]=0,E[{\mathbf{h}}_t{\mathbf{h}}_s^{\ast }]=\frac{1}{2}{\mathbf{\kappa }}_2(t-s),s,t\ge 0.$$

$$m{\ddot{\mathbf{x}}}_t=\mathbf{F}({\mathbf{x}}_t)-\mathbf{g}({\mathbf{x}}_t){\int }_0^t{\mathbf{\kappa }}_2(t-s)\mathbf{h}({\mathbf{x}}_s){\dot{\mathbf{x}}}_{s}ds+\mathbf{\sigma }({\mathbf{x}}_t){\mathbf{h}}_t.$$

$$E[{\mathbf{\eta }}_t^0({\mathbf{\eta }}_s^0{)}^{\ast }]=\frac{1}{2}{\mathbf{\kappa }}_1(t-s)$$

$$E[{\mathbf{h}}_t^0({\mathbf{h}}_s^0{)}^{\ast }]=\frac{1}{2}{\mathbf{\kappa }}_2(t-s)$$

$$m_0{ϵ}^{\mu }{\ddot{\mathbf{x}}}_t^{ϵ}=\mathbf{F}({\mathbf{x}}_t^{ϵ})-\frac{\mathbf{g}({\mathbf{x}}_t^{ϵ})}{{\tau }_{\kappa }{ϵ}^{\theta }}{\int }_0^t{\mathbf{C}}_1{e}^{-\frac{{\mathbf{\Gamma }}_1}{{\tau }_{\kappa }{ϵ}^{\theta }}(t-s)}{\mathbf{M}}_1{\mathbf{C}}_1^{\ast }\mathbf{h}({\mathbf{x}}_s^{ϵ}){\dot{\mathbf{x}}}_s^{ϵ}ds+\mathbf{\sigma }({\mathbf{x}}_t^{ϵ}){\mathbf{C}}_2{\mathbf{\beta }}_t^{ϵ},$$

$$d{\mathbf{\beta }}_t^{ϵ}=-\frac{{\mathbf{\Gamma }}_2}{{\tau }_{\xi }{ϵ}^{\nu }}{\mathbf{\beta }}_t^{ϵ}dt+\frac{{\mathbf{\Sigma }}_2}{{\tau }_{\xi }{ϵ}^{\nu }}d{\mathbf{W}}_t^{(q_2)},$$

$$m_0{ϵ}^{\mu }{\ddot{\mathbf{x}}}_t^{ϵ}=\mathbf{F}({\mathbf{x}}_t^{ϵ})-\frac{\mathbf{g}({\mathbf{x}}_t^{ϵ})}{{\tau }_{\kappa }{ϵ}^{\theta }}{\int }_0^t\mathbf{A}{e}^{-\frac{\mathbf{A}}{{\tau }_{\kappa }{ϵ}^{\theta }}(t-s)}\mathbf{h}({\mathbf{x}}_s^{ϵ}){\dot{\mathbf{x}}}_s^{ϵ}ds+\mathbf{\sigma }({\mathbf{x}}_t^{ϵ}){\mathbf{\eta }}_t^{ϵ},$$

$$d{\mathbf{\eta }}_t^{ϵ}=-\frac{\mathbf{A}}{{\tau }_{\eta }{ϵ}^{\nu }}{\mathbf{\eta }}_t^{ϵ}dt+\frac{\mathbf{A}}{{\tau }_{\eta }{ϵ}^{\nu }}d{\mathbf{W}}_t^{(d)}.$$

$$m_0{ϵ}^{\mu }{\ddot{\mathbf{x}}}_t^{ϵ}$$

$$=\mathbf{F}({\mathbf{x}}_t^{ϵ})-\frac{\mathbf{g}({\mathbf{x}}_t^{ϵ})}{{\tau }_{\kappa }{ϵ}^{\theta }}{\int }_0^t{e}^{-{\mathbf{\Omega }}^2\frac{(t-s)}{2{\tau }_{\kappa }{ϵ}^{\theta }}}\left[\cos \left(\frac{{\mathbf{\Omega }}_0}{{\tau }_{\kappa }{ϵ}^{\theta }}(t-s)\right)+\frac{{\mathbf{\Omega }}_1}{2}\sin \left(\frac{{\mathbf{\Omega }}_0}{{\tau }_{\kappa }{ϵ}^{\theta }}(t-s)\right)\right]\mathbf{h}({\mathbf{x}}_s^{ϵ}){\dot{\mathbf{x}}}_s^{ϵ}ds+\mathbf{\sigma }({\mathbf{x}}_t^{ϵ}){\mathbf{h}}_t^{ϵ},$$

$$d{\mathbf{h}}_t^{ϵ}$$

$$d{\mathbf{u}}_t^{ϵ}$$

$$d{\mathbf{x}}_t^m$$

$$md{\mathbf{v}}_t^m$$

Full text

1 Introduction

From physical sciences to social sciences, many phenomena are modeled by noisy dynamical systems. In many such systems, several widely separated time scales are present. The system obtained in the homogenization limit, in which the fast time scales go to zero, is simpler than the original one, while often retaining the essential features of its dynamics majda2001mathematical ; givon2004extracting ; Pavliotis-TwoFast ; Pavliotis . On the other hand, the different fast time scales may compete and this competition is reflected in the homogenized equations.

Of particular interest is the model of a Brownian particle interacting with the environment nelson1967dynamical . The usual model for such system neglects the memory effects, representing the interaction of the particle with the environment as a sum of an instantaneous damping force and a white noise. Although such an idealized model generally gives a good approximate description of the dynamics of the particle, there are situations where the memory effects play an important role, for instance when the particle is subject to a hydrodynamic interaction franosch2011resonances , or when the particle is an atom embedded in a condensed-matter heat bath Groblacher2015 . In addition, the stochastic forcing introduced by the environment is often more accurately modeled by a colored noise than by white noise.

In this paper, we study a class of generalized Langevin equations (GLEs), with state-dependent drift and diffusion coefficients, driven by colored noise. They provide a realistic description of the dynamics of a classical Brownian particle in an inhomogeneous environment; their solutions are not Markov processes. We are interested in the limiting behavior of the particle when the characteristic time scales become small and in how competition of the time scales, as well as inhomogeneity of the environment, impact its limiting dynamics. The main mathematical result of this paper is Theorem 4.1, in which we derive the homogenization limit for a general class of non-Markovian systems. Special cases are studied in some details to obtain more explicit results. Their physical relevance is illustrated by an application to thermophoresis models.

The paper is organized as follows. In Section 2, we introduce and discuss a class of GLEs, as well as its two sub-classes, to be studied in this paper. In Section 3, we revisit the Smoluchowski-Kramers limit for a class of SDEs with state-dependent drift and diffusion coefficients, under a weaker assumption on the spectrum of the damping matrix than that used in earlier work hottovy2015smoluchowski . Using this result of Section 3, we study homogenization for the GLEs in Section 4. We specialize the study to the two sub-classes of models in Section 5. In Section 6, we apply the results obtained in the previous sections to study the thermophoresis of a Brownian particle in a non-equilibrium heat bath. We end the paper by giving the conclusions and final remarks in Section 7. The appendices provide some technical results used in the main paper, as well as physical motivation for the form of the GLEs studied here. In Appendix A we provide a variant of a (heuristic) derivation of the equations studied in this paper from Hamiltonian model of a particle interacting with a system of harmonic oscillators. Appendix B contains a sketch of the proof of Theorem 3.1.

2 Generalized Langevin Equations (GLEs)

2.1 GLEs as Non-Markovian Models

We consider a class of non-Markovian Langevin equations, with state-dependent coefficients, that describe the dynamics of a particle moving in a force field and interacting with the environment. Let $\boldsymbol{x}_{t}\in\mathbb{R}^{d}$, $t\geq 0$, be the position of the particle. The evolution of position, $\boldsymbol{x}_{t}$, is given by the solution to the following stochastic integro-differential equation (SIDE):

[TABLE]

with the initial conditions (here the initial time is chosen to be $t=0$):

[TABLE]

The initial conditions $\boldsymbol{x}$ and $\boldsymbol{v}$ are random variables independent of the process $\{\boldsymbol{\xi}_{t}:\ t\geq 0\}$. Our motivation to study the SIDE (1) is that study of microscopic dynamics leads naturally to equations of this form (see Appendix A).

Here and throughout the paper, overdot denotes derivative with respect to time $t$, the superscript ∗ denotes conjugate transposition of matrices or vectors and $E$ denotes expectation. In the SIDE (1), $m>0$ is the mass of the particle, the matrix-valued functions $\boldsymbol{g}:\mathbb{R}^{d}\to\mathbb{R}^{d\times q}$, $\boldsymbol{h}:\mathbb{R}^{d}\to\mathbb{R}^{q\times d}$ and $\boldsymbol{\sigma}:\mathbb{R}^{r}\to\mathbb{R}^{d\times r}$ are the state-dependent coefficients of the equation, and $\boldsymbol{F}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is a force field acting on the particle. Here $d$, $q$ and $r$ are, possibly distinct, positive integers. The second term on the right hand side of (1) represents the drag experienced by the particle and the last term models the noise.

The matrix-valued function $\boldsymbol{\kappa}:\mathbb{R}\to\mathbb{R}^{q\times q}$ is a memory function which is Bohl, i.e. the matrix elements of $\boldsymbol{\kappa}(t)$ are finite linear combinations of the functions of the form $t^{k}e^{\alpha t}\cos(\omega t)$ and $t^{k}e^{\alpha t}\sin(\omega t)$, where $k$ is an integer and $\alpha$ and $\omega$ are real numbers. For properties of Bohl functions, we refer to Chapter 2 of trentelman2002control . The noise process $\boldsymbol{\xi}_{t}$ is a $r$-dimensional mean zero stationary real-valued Gaussian vector process having a Bohl covariance function, $\boldsymbol{R}(t):=E\boldsymbol{\xi}_{t}\boldsymbol{\xi}_{0}^{*}=\boldsymbol{R}^{*}(-t)$, and, therefore, its spectral density, $\boldsymbol{S}(\omega):=\int_{-\infty}^{\infty}\boldsymbol{R}(t)e^{-i\omega t}dt$, is a rational function willems1980stochastic .

The SIDE (1) is a non-Markovian Langevin equation, since its solution at time $t$ depends on the entire past. Two of its terms are different than those in the usual Langevin equations. One of them is the drag term, which here involves an integral over the particle’s past velocities with a memory kernel $\boldsymbol{\kappa}(t-s)$. It describes the state-dependent dissipation which comprises the back-action effects of the environment up to current time. The other term, involving a Gaussian colored noise $\boldsymbol{\xi}_{t}$, is a multiplicative noise term, also arising from interaction of the particle with the environment. Therefore, (1) is a generalized Langevin equation (GLE), which in its most basic form was first introduced by Mori in mori1965transport and subsequently used to model many systems in statistical physics Kubo_fd ; toda2012statistical ; goychuk2012viscoelastic .

As remarked by van Kampen in van1998remarks , “Non-Markov is the rule, Markov is the exception”. Therefore, it is not surprising that non-Markovian equations (including those of form (1)) find numerous applications and thus have been studied widely in the mathematical, physical and engineering literature. See, for instance, luczka2005non ; samorodnitsky1994stable for surveys of non-Markovian processes, PhysRevB.89.134303 ; mckinley2009transient ; adelman1976generalized for physical applications and Ottobre for asymptotic analysis.

Note that the Gaussian process $\boldsymbol{\xi}_{t}$ which drives the SIDE (1) is not assumed to be Markov. The assumptions we made on its covariance will allow us to present it as a projection of a Markov process in a (typically higher-dimensional) space. This approach, which originated in stochastic control theory kalman1960new , is called stochastic realization. We describe it in detail below.

Let $\boldsymbol{\Gamma}_{1}\in\mathbb{R}^{d_{1}\times d_{1}}$, $\boldsymbol{M}_{1}\in\mathbb{R}^{d_{1}\times d_{1}}$, $\boldsymbol{C}_{1}\in\mathbb{R}^{q\times d_{1}}$, $\boldsymbol{\Sigma}_{1}\in\mathbb{R}^{d_{1}\times q_{1}}$, $\boldsymbol{\Gamma}_{2}\in\mathbb{R}^{d_{2}\times d_{2}}$, $\boldsymbol{M}_{2}\in\mathbb{R}^{d_{2}\times d_{2}}$, $\boldsymbol{C}_{2}\in\mathbb{R}^{r\times d_{2}}$, $\boldsymbol{\Sigma}_{2}\in\mathbb{R}^{d_{2}\times q_{2}}$ be constant matrices, where $d_{1},d_{2},q_{1},q_{2}$, $q$ and $r$ are positive integers. In this paper, we study the class of SIDE (1), with the memory function defined in terms of the triple $(\boldsymbol{\Gamma}_{1},\boldsymbol{M}_{1},\boldsymbol{C}_{1})$ of matrices as follows:

[TABLE]

The noise process is the mean zero, stationary Gaussian vector process, whose covariance will be expressed in terms of the triple $(\boldsymbol{\Gamma}_{2},\boldsymbol{M}_{2},\boldsymbol{C}_{2})$. More precisely, we define it as:

[TABLE]

where $\boldsymbol{\beta}_{t}$ is the solution to the Itô SDE:

[TABLE]

with the initial condition, $\boldsymbol{\beta}_{0}$, normally distributed with zero mean and covariance $\boldsymbol{M}_{2}$. Here, $\boldsymbol{W}^{(q_{2})}_{t}$ denotes a $q_{2}$-dimensional Wiener process and is independent of $\boldsymbol{\beta}$. Throughout the paper the dimension of the Wiener process will be specified by the superscript.

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

[TABLE]

It follows from positive stability of $\boldsymbol{\Gamma}_{i}$ that this equation indeed has a unique solution bellman1997introduction . The covariance matrix, $\boldsymbol{R}(t)\in\mathbb{R}^{r\times r}$, of the noise process is therefore expressed in terms of the matrices $(\boldsymbol{\Gamma}_{2},\boldsymbol{M}_{2},\boldsymbol{C}_{2})$ as follows:

[TABLE]

and therefore the triple $(\boldsymbol{\Gamma}_{2},\boldsymbol{M}_{2},\boldsymbol{C}_{2})$ completely specifies the probability distribution of $\boldsymbol{\xi}_{t}$. It is worth mentioning that the triples that specify the memory function in (3) and the noise process in (4) are only unique up to the following transformations:

[TABLE]

where $i=1,2$ and $\boldsymbol{T}_{i}$ is any invertible matrices of appropriate dimensions.

The triple $(\boldsymbol{\Gamma}_{2},\boldsymbol{M}_{2},\boldsymbol{C}_{2})$ above is called a (weak) stochastic realization of the covariance matrix $\boldsymbol{R}(t)$ in the well established theory of stochastic realization, which is concerned with solving the inverse problem of stationary covariance generation (see lindquist1985realization ; lindquist2015linear ). Any zero mean stationary Gaussian process, $\boldsymbol{\xi}^{\prime}_{t}$, having a Bohl covariance function, can be realized as a projection of a Gaussian Markov process in the above way. Let us remark that Gaussian processes with Bohl covariance functions are precisely those with rational spectral density willems1980stochastic .

Our approach allows to consider the most general Gaussian noises that can be realized in a finite-dimensional state space in the above way (i.e. as a linear transformation of a Gaussian Markov process). In fact, the condition on the covariance function to have entries in the Bohl class is necessary and sufficient for solvability of the problem of stochastic realization of stationary Gaussian processes. We refer to the propositions and theorems on page 303-308 of willems1980stochastic for a brief exposition of stochastic realization problems.

Remark 1.

Physically, the choice of the matrices $\boldsymbol{\Gamma}_{2},\boldsymbol{M}_{2},\boldsymbol{C}_{2}$ specifies the characteristic time scales (eigenvalues of $\boldsymbol{\Gamma}_{2}^{-1}$) present in the environment, introduces the initial state of a stationary Markovian Gaussian noise and selects the parts of the prepared Markovian noise that are (partially) observed, respectively. In other words, we have assumed that the noise in the SIDE (1) is realized or “experimentally prepared” by the above triple of matrices.

For our homogenization study of the equation (1) we need the effective damping constant,

[TABLE]

and the effective diffusion constant,

[TABLE]

to be invertible (see Section 2.2). This is equivalent to the matrices $\boldsymbol{C}_{i}$ having full rank. Homogenization for a class of systems with vanishing effective damping and/or diffusion constant bao2005non will be explored in our future work.

With the above definitions of memory kernel and noise process, the SIDE (1) becomes:

[TABLE]

where $\boldsymbol{\beta}_{t}$ is the solution to the SDE (5). To illustrate the results of the general study in important special cases (which will also be used later in applications), we consider two representative classes of SIDE (1). The driving Gaussian colored noise is Markovian in the first class and non-Markovian in the second. We set $d=d_{1}=d_{2}=q_{1}=q_{2}=q=r$ in the following examples.

• (i)

Example of a SIDE driven by a Markovian colored noise. The memory kernel is given by an exponential function, i.e.

[TABLE]

where $\boldsymbol{A}\in\mathbb{R}^{d\times d}$ is a constant diagonal matrix with positive eigenvalues. The driving noise is the Ornstein-Uhlenbeck (OU) process, $\boldsymbol{\xi}_{t}=\boldsymbol{\eta}_{t}\in\mathbb{R}^{d}$, i.e. a mean zero stationary Gaussian process which is the solution to the SDE:

[TABLE]

In order for the process $\boldsymbol{\eta}_{t}$ to be stationary, the initial condition has to be distributed according to the (unique) stationary measure of the Markov process defined by the above equation, i.e. $\boldsymbol{\eta}_{0}=\boldsymbol{\eta}$ is normally distributed with zero mean and covariance $\boldsymbol{A}/2$. The mean and the covariance of $\boldsymbol{\eta}_{t}$ are given by:

[TABLE]

The resulting SIDE reads:

[TABLE]

Let us note that Ornstein-Uhlenbeck processes are the only stationary, ergodic, Gaussian, Markov processes with continuous covariance functions pavliotis2014stochastic . When all diagonal entries of $\boldsymbol{A}$ go to infinity, the OU process approaches the white noise. For details on OU processes, see for instance pavliotis2014stochastic and Section 2 of hottovy2015small .

• (ii)

Example of a SIDE driven by a non-Markovian colored noise. The memory kernel is given by an oscillatory function whose amplitude is exponentially decaying, i.e. $\boldsymbol{\kappa}(t-s):=\boldsymbol{\kappa}_{2}(t-s)$, a diagonal matrix with the diagonal entries:

[TABLE]

where, for $i=1,\dots,d$, $\tau_{ii}$ is a positive constant, $\omega_{ii}$ is a real constant, $\omega^{0}_{ii}:=\omega_{ii}\sqrt{1-\omega_{ii}^{2}/4}$, $\tilde{\omega}^{0}_{ii}:=\omega_{ii}\sqrt{\omega_{ii}^{2}/4-1}$, $\omega_{ii}^{1}:=\omega_{ii}/\sqrt{1-\omega_{ii}^{2}/4}$, and $\tilde{\omega}_{ii}^{1}:=\omega_{ii}/\sqrt{\omega_{ii}^{2}/4-1}$.

Let $\boldsymbol{\tau}$ be constant diagonal matrix with the positive eigenvalues $(\tau_{jj})_{j=1}^{d}$, $\boldsymbol{\Omega}$ be constant diagonal matrix with the real eigenvalues $(\omega_{jj})_{j=1}^{d}$, $\boldsymbol{\Omega}_{0}$ be constant $d\times d$ diagonal matrix with the eigenvalues $\omega_{jj}\sqrt{1-\omega_{jj}^{2}/4}$ (if $|\omega_{jj}|<2$) and $i\omega_{jj}\sqrt{\omega_{jj}^{2}/4-1}$ (if $|\omega_{jj}|>2$), and $\boldsymbol{\Omega}_{1}$ be constant $d\times d$ diagonal matrix with the eigenvalues $\omega_{jj}/\sqrt{1-\omega_{jj}^{2}/4}$ (if $|\omega_{jj}|<2$) and $-i\omega_{jj}/\sqrt{\omega_{jj}^{2}/4-1}$ (if $|\omega_{jj}|>2$), where $i$ is the imaginary unit.

The driving noise is given by the harmonic noise process, $\boldsymbol{\xi}_{t}=\boldsymbol{h}_{t}\in\mathbb{R}^{d}$, i.e. a mean zero stationary Gaussian process which is the solution to the SDE system:

[TABLE]

with the initial conditions, $\boldsymbol{h}_{0}$ and $\boldsymbol{u}_{0}$, distributed according to the (unique) stationary measure of the above SDE system. The mean and the covariance of $\boldsymbol{h}_{t}$ are given by:

[TABLE]

Note that $\boldsymbol{h}_{t}$ is not a Markov process (but the process $(\boldsymbol{h}_{t},\boldsymbol{u}_{t})$ is).

The resulting SIDE reads:

[TABLE]

The harmonic noise is an approximation of the white noise, smoother than the Ornstein-Uhlenbeck process. It can be shown that in the limit $\omega_{ii}\to\infty$ (for all $i$) the process $\boldsymbol{h}_{t}$ converges to the Ornstein-Uhlenbeck process whose $i$th component process has correlation time $\tau_{ii}$, whereas in the limit $\tau_{ii}\to 0$ (for all $i$) the process $\boldsymbol{h}_{t}$ converges to the white noise. For detailed properties of harmonic noise process, see for instance schimansky1990harmonic or the Appendix in McDaniel14 . We remark that the harmonic noise is one of the simplest examples of a non-Markovian process and its use as the driving noise in the SIDE (1) is a natural choice that models the environment as a bath of damped harmonic oscillators hanggi1993can .

Remark 2.

Note that in the SIDEs for the above two sub-classes, the dimension of the driving Wiener process, $\boldsymbol{W}_{t}^{(d)}$, is the same as that of the colored noise processes $\boldsymbol{\eta}_{t}$ and $\boldsymbol{h}_{t}$, as well as the processes, $\boldsymbol{x}_{t}$ and $\boldsymbol{v}_{t}$. One could as well consider realizing the noise processes using a driving Wiener process of different dimension. Our choice of working with the same dimensions is for the sake of convenience as it will help to simplify the exposition.

Remark 3.

Without loss of generality (due to (8)), we have taken $\boldsymbol{A}$ and $\boldsymbol{\Omega}$ to be diagonal.

Remark 4.

In cases of particular interest in statistical physics, the triples $(\boldsymbol{\Gamma}_{i},\boldsymbol{M}_{i},\boldsymbol{C}_{i})$ coincide, up to the transformations in (8), for $i=1,2$; $\boldsymbol{h}=\boldsymbol{g}^{*}$ and $\boldsymbol{g}$ is proportional to $\boldsymbol{\sigma}$, with the proportionality factor equals $k_{B}T$, where $k_{B}$ denotes the Boltzmann constant and $T>0$ is temperature of the environment (see Appendix A). In this case, we have $d_{1}=d_{2}$ and $q=r$. In particular, for the two sub-classes above we have

[TABLE]

for the first sub-class and

[TABLE]

for the second sub-class. In such cases, the SIDEs describe a particle interacting with an equilibrium heat bath at a temperature $T$, whose dynamics satisfy the fluctuation-dissipation relation toda2012statistical ; zwanzig2001nonequilibrium .

2.2 Homogenization of SIDEs: Discussion and Statement of the Problem

There are three characteristic time scales defining the non-Markovian dynamics described by the SIDE (1):

• (i)

the inertial time scale, $\lambda_{m}$, proportional to $m$, whose physical significance is the relaxation time of the particle velocity process $\boldsymbol{v}_{t}:=\dot{\boldsymbol{x}}_{t}$. The limit $\lambda_{m}\to 0$ is equivalent to the limit $m\to 0$;

• (ii)

the memory time scale, $\lambda_{\kappa}$, defined as the inverse rate of exponential decay of the memory kernel $\boldsymbol{\kappa}(t-s)$;

• (iii)

the noise correlation time scale, $\lambda_{\xi}$.

For the purpose of general multiscale analysis, we set $m=m_{0}\epsilon^{\mu}$, $\lambda_{\kappa}=\tau_{\kappa}\epsilon^{\theta}$ and $\lambda_{\xi}=\tau_{\xi}\epsilon^{\nu}$, where $\epsilon>0$ is a parameter which will be taken to zero, $m_{0}$, $\tau_{\kappa}$, $\tau_{\xi}$ are (fixed) proportionality constants, and $\mu,\theta,\nu$ are positive constants (exponents), specifying the orders at which the time scales $\lambda_{m},\lambda_{\kappa},\lambda_{\xi}$ vanish respectively. We consider a family of SIDEs, parametrized by $\epsilon$, with the inertial time scale $\lambda_{m}$ proportional to $m_{0}\epsilon^{\mu}$, memory time scale $\lambda_{\kappa}=\tau_{\kappa}\epsilon^{\theta}$ and noise correlation time scale $\lambda_{\xi}=\tau_{\xi}\epsilon^{\nu}$, to be defined in the following.

We replace $m$ with $m_{0}\epsilon^{\mu}$, $\boldsymbol{\Gamma}_{1}$ with $\boldsymbol{\Gamma}_{1}/(\tau_{\kappa}\epsilon^{\theta})$, $\boldsymbol{M}_{1}$ with $\boldsymbol{M}_{1}/(\tau_{\kappa}\epsilon^{\theta})$, and $\boldsymbol{x}_{t}$ with $\boldsymbol{x}_{t}^{\epsilon}$ in (11). Also, we substitute $\boldsymbol{\Gamma}_{2}$ with $\boldsymbol{\Gamma}_{2}/(\tau_{\xi}\epsilon^{\nu})$, $\boldsymbol{\Sigma}_{2}$ with $\boldsymbol{\Sigma}_{2}/(\tau_{\xi}\epsilon^{\nu})$, and $\boldsymbol{\beta}_{t}$ with $\boldsymbol{\beta}_{t}^{\epsilon}$ in (5). The SIDE (11) then becomes:

[TABLE]

with the initial conditions, $\boldsymbol{x}^{\epsilon}_{0}=\boldsymbol{x}$ and $\boldsymbol{v}^{\epsilon}_{0}=\boldsymbol{v}$, where $\boldsymbol{\beta}^{\epsilon}_{t}$ is a process, with correlation time $\tau_{\xi}\epsilon^{\nu}$, satisfying the SDE:

[TABLE]

with the initial condition, $\boldsymbol{\beta}^{\epsilon}_{0}$, normally distributed with zero mean and covariance of $\boldsymbol{M}_{2}/(\tau_{\xi}\epsilon^{\nu})$.

We will also perform similar analysis on the two sub-classes of SIDE, in which case:

• (i)

the SIDE (15) becomes (with $m:=m_{0}\epsilon^{\mu}$, $\boldsymbol{A}$ in the formula for $\boldsymbol{\kappa}_{1}$ replaced by $\boldsymbol{A}/(\tau_{\kappa}\epsilon^{\theta})$, $\boldsymbol{A}$ in $\eqref{ou}$ replaced by $\boldsymbol{A}/(\tau_{\eta}\epsilon^{\nu})$, $\boldsymbol{x}_{t}$ replaced by $\boldsymbol{x}_{t}^{\epsilon}$ and $\boldsymbol{\eta}_{t}$ replaced by $\boldsymbol{\eta}_{t}^{\epsilon}$):

[TABLE]

where $\boldsymbol{\eta}^{\epsilon}_{t}$ is the Ornstein-Uhlenbeck process with the correlation time $\tau_{\eta}\epsilon^{\nu}$, i.e. it is a process satisfying the SDE:

[TABLE]

• (ii)

the SIDE (20) becomes (with $m:=m_{0}\epsilon^{\mu}$, $\tau_{ii}:=\tau_{\kappa}\epsilon^{\theta}$ in the formula for ($\boldsymbol{\kappa}_{2})_{ii}$ in (16), $\boldsymbol{x}_{t}$ replaced by $\boldsymbol{x}^{\epsilon}_{t}$, $\boldsymbol{h}_{t}$, $\boldsymbol{u}_{t}$ replaced by $\boldsymbol{h}^{\epsilon}_{t}$, $\boldsymbol{u}^{\epsilon}_{t}$ respectively and $\boldsymbol{\tau}:=\tau_{h}\epsilon^{\nu}\boldsymbol{I}$ in $\eqref{unscaled_har1}$-$\eqref{unscaled_har2}$):

[TABLE]

where $\boldsymbol{h}^{\epsilon}_{t}$ is the harmonic noise process with the correlation time $\tau_{h}\epsilon^{\nu}$, i.e. it is a process satisfying the SDE system:

[TABLE]

Both SIDEs have the initial conditions $\boldsymbol{x}^{\epsilon}_{0}=\boldsymbol{x},\ \dot{\boldsymbol{x}}^{\epsilon}_{0}=\boldsymbol{v}$. The initial conditions $\boldsymbol{\eta}^{\epsilon}_{0}$ (respectively, $\boldsymbol{h}^{\epsilon}_{0}$ and $\boldsymbol{u}^{\epsilon}_{0}$) are distributed according to the stationary measure of the SDE that the process $\boldsymbol{\eta}^{\epsilon}_{t}$ (respectively, $\boldsymbol{h}^{\epsilon}_{t}$ and $\boldsymbol{u}^{\epsilon}_{t}$) satisfies.

In this paper we set $\mu=\theta=\nu$, which is the case when all the characteristic time scales are of comparable magnitude in the limit as $\epsilon\to 0$. Our main goal is to derive a limiting equation for the (slow) $\boldsymbol{x}^{\epsilon}$-component of the process solving the equations (23)-(24), including the special cases $\eqref{goal2}$-$\eqref{rescaled-ou}$ and $\eqref{goal3}$-$\eqref{rescaled_h2}$, in the limit as $\epsilon\to 0$, in a strong pathwise sense.

We explain the motivation behind the above rescalings. The rescaling of the memory kernels $\boldsymbol{\kappa}(t-s)$, $\boldsymbol{\kappa}_{1}(t-s)$, $\boldsymbol{\kappa}_{2}(t-s)$ is such that in the limit $\epsilon\to 0$ the rescaled memory kernels converge to $\boldsymbol{K}_{1}\delta(t)$ formally, where $\delta(t)$ is the Dirac-delta function and $\boldsymbol{K}_{1}$ is the effective damping constant defined in (9). On the other hand, the noise processes $\boldsymbol{\xi}^{\epsilon}_{t}=\boldsymbol{C}_{2}\boldsymbol{\beta}^{\epsilon}_{t}$, $\boldsymbol{\eta}^{\epsilon}_{t}$ and $\boldsymbol{h}^{\epsilon}_{t}$ converge to white noise processes in the limit $\epsilon\to 0$.

3 Smoluchowski-Kramers Limit of SDE’s Revisited

Let $(\boldsymbol{x}_{t}^{m},\boldsymbol{v}^{m}_{t})\in\mathbb{R}^{n}\times\mathbb{R}^{n}$, where $t\in[0,T]$, be a family of solutions (parametrized by a positive constant $m$) to the following SDEs:

[TABLE]

In the SDEs above, $m>0$ is the mass of the particle, $\boldsymbol{F}:\mathbb{R}^{n}\to\mathbb{R}^{n}$, $\boldsymbol{\gamma}:\mathbb{R}^{n}\to\mathbb{R}^{n\times n}$, $\boldsymbol{\sigma}:\mathbb{R}^{n}\to\mathbb{R}^{n\times k}$, and $\boldsymbol{W}^{(k)}$ is a $k$-dimensional Wiener process on the filtered probability space $(\Omega,\mathcal{F},\mathcal{F}_{t},\mathbb{P})$ satisfying the usual conditions (karatzas2012Brownian, ). The initial conditions are given by $\boldsymbol{x}^{m}_{0}=\boldsymbol{x},\ \boldsymbol{v}^{m}_{0}=\boldsymbol{v}^{m}$. The above SDE system models diffusive phenomena in cases where the damping coefficient $\boldsymbol{\gamma}$ and diffusion coefficient $\boldsymbol{\sigma}$ are state-dependent.

The Smoluchowski-Kramers limit (or the small mass limit) of the system (30)-(31) has been studied in (hottovy2015smoluchowski, ; Hottovy12, ; Herzog2016, ; birrell2017small, ; birrell2017homogenization, ). The main result in (birrell2017homogenization, ) says that, under certain assumptions, the $\boldsymbol{x}^{m}$-component of the solution to (30)-(31) converges (in a strong pathwise sense), as $m\to 0$, to the solution of a homogenized SDE that contains in particular the so-called noise-induced drift, that was not present in the pre-limit SDEs (see Theorem 3.1 for a precise statement). The presence of such noise-induced drift in the homogenized equation is a consequence of the state-dependence of the damping coefficient (and therefore also the diffusion coefficient if the system satisfies a fluctuation-dissipation relation). For an overview of the noise-induced drift phenomena, we refer to the review article (volpe2016effective, ).

In all the works mentioned previously, the spectral assumption made on the matrix $\boldsymbol{\gamma}$ was that the symmetrized damping matrix $\frac{1}{2}(\boldsymbol{\gamma}+\boldsymbol{\gamma}^{*})$ is uniformly positive definite (i.e. its smallest eigevalue is positive uniformly in $\boldsymbol{x}$). The same results can be obtained under a weaker assumption that matrix $\boldsymbol{\gamma}$ is uniformly positive stable, i.e. all real parts of the eigenvalues of $\boldsymbol{\gamma}$ are positive uniformly in $\boldsymbol{x}$ (horn1994topics, ).

Notation. Here and in the following, we use Einstein summation convention on repeated indices. The Euclidean norm of a 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})$ denotes the gradient vector $(\frac{\partial[f]_{jk}(\boldsymbol{y})}{\partial[y]_{1}},\dots,\frac{\partial[f]_{jk}(\boldsymbol{y})}{\partial[y]_{n_{1}}})\in\mathbb{R}^{n_{1}}$ for every $j,k$. The symbol $\mathbb{E}$ denotes expectation with respect to $\mathbb{P}$.

We make the following assumptions.

Assumption 1**.**

For every $\boldsymbol{x}\in\mathbb{R}^{n}$, the functions $\boldsymbol{F}(\boldsymbol{x})$ and $\boldsymbol{\sigma}(\boldsymbol{x})$ are continuous, bounded and Lipschitz in $\boldsymbol{x}$, whereas the functions $\boldsymbol{\gamma}(\boldsymbol{x})$ and $(\boldsymbol{\gamma})_{\boldsymbol{x}}(\boldsymbol{x})$ are continuously differentiable, bounded and Lipschitz in $\boldsymbol{x}$. Moreover, $(\boldsymbol{\gamma})_{\boldsymbol{x}\boldsymbol{x}}(\boldsymbol{x})$ is bounded for every $\boldsymbol{x}\in\mathbb{R}^{n}$.

Assumption 2**.**

The matrix $\boldsymbol{\gamma}$ is uniformly positive stable, i.e. all real parts of the eigenvalues of $\boldsymbol{\gamma}(\boldsymbol{x})$ are bounded below by some constant $2\kappa>0$, uniformly in $\boldsymbol{x}\in\mathbb{R}^{n}$.

Assumption 3**.**

The initial condition $\boldsymbol{x}^{m}_{0}=\boldsymbol{x}_{0}$ is a random variable independent of $m$ and has finite moments of all orders, i.e. $\mathbb{E}|\boldsymbol{x}|^{p}<\infty$ for all $p>0$. The initial condition $\boldsymbol{v}^{m}_{0}$ is a random variable that possibly depends on $m$ and we assume that for every $p>0$, $\mathbb{E}|m\boldsymbol{v}^{m}|^{p}=O(m^{\alpha})$ as $m\to 0$, for some $\alpha\geq p/2$.

Assumption 4**.**

The global solutions, defined on $[0,T]$, to the pre-limit SDEs (30)-(31) and to the limiting SDE (33) a.s. exist and are unique for all $m>0$ (i.e. there are no explosions).

We now state the result.

Theorem 3.1.

Suppose that the SDE system $\eqref{gsk1}$-$\eqref{gsk2}$ satisfies Assumption 1-4. Let $(\boldsymbol{x}^{m}_{t},\boldsymbol{v}^{m}_{t})\in\mathbb{R}^{n}\times\mathbb{R}^{n}$ be its solution, with the initial condition $(\boldsymbol{x},\boldsymbol{v}^{m})$. Let $\boldsymbol{X}_{t}\in\mathbb{R}^{n}$ be the solution to the following Itô SDE with initial position $\boldsymbol{X}_{0}=\boldsymbol{x}$:

[TABLE]

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

[TABLE]

and $\boldsymbol{J}$ is the unique matrix solving the Lyapunov equation

[TABLE]

Then the process $\boldsymbol{x}^{m}_{t}$ converges, as $m\to 0$, to the solution $\boldsymbol{X}_{t}$, of the Itô SDE (33), in the following sense: for all finite $T>0$,

[TABLE]

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

We end this section with a few remarks concerning the statements in Theorem 3.1.

Remark 5.

• (i)

Because of the relaxed spectral assumption on $\boldsymbol{\gamma}$, a new idea has to be used to prove decay estimates for solutions of the velocity equation. Once this is done, Theorem 1 can be proven using the technique of birrell2017homogenization (note that Assumption 1 is essentially the same as the assumptions in Appendix A of birrell2017homogenization , specialized to the present case). In Appendix B we give a sketch of the proof of Theorem 3.1, pointing out the necessary modifications. The reader is referred to birrell2017homogenization for more details.

• (ii)

Our assumption on the initial variable $\boldsymbol{v}_{0}^{m}$ implies that the initial average kinetic energy, $K(\boldsymbol{v}^{m}):=\mathbb{E}\frac{1}{2}m|\boldsymbol{v}^{m}|^{2}$, does not blow up (but can possibly vanish) as $m\to 0$. This is analogous to the Assumption 2.4 in birrell2017homogenization and it is more general than the one in hottovy2015smoluchowski .

• (iii)

With slightly more work and additional assumptions, one could prove the statement in Assumption 4 from Assumptions 1-3. However, such existence and uniqueness result is not the focus of this paper and, therefore, we choose to take the existence and uniqueness result for granted in Assumption 4.

• (iv)

We make no claim that Assumptions 2-4 are as weak or as general as possible. In particular, the boundedness assumption on the coefficients of the SDEs could be relaxed (for instance, using the techniques in Herzog2016 ) and the initial condition $\boldsymbol{x}$ could have some dependence on $m$ (see, for instance, birrell2017homogenization ) at the cost of more technicalities.

The main focus of our revisit here is to point out that the result in hottovy2015smoluchowski still holds with a relaxed spectral assumption on the matrix $\boldsymbol{\gamma}$ and with the initial condition $\boldsymbol{v}_{0}^{m}$ possibly dependent on $m$ – this will be important for applications in later sections (see also Remark 10).

4 Homogenization for Generalized Langevin Dynamics

In this section, we study homogenization for the system of equations (23)-(24) (with $\mu=\theta=\nu$) by taking the limit as $\epsilon\to 0$, under appropriate assumptions.

Without loss of generality, we set $\mu=\theta=\nu=1$. We cast (23)-(24) as the system of SDEs for the Markov process $(\boldsymbol{x}^{\epsilon}_{t},\boldsymbol{v}^{\epsilon}_{t},\boldsymbol{z}^{\epsilon}_{t},\boldsymbol{y}^{\epsilon}_{t},\boldsymbol{\zeta}^{\epsilon}_{t},\boldsymbol{\beta}^{\epsilon}_{t})$ on the state space $\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d_{1}}\times\mathbb{R}^{d_{1}}\times\mathbb{R}^{d_{2}}\times\mathbb{R}^{d_{2}}$:

[TABLE]

where we have defined the auxiliary process

[TABLE]

Here, the initial conditions $\boldsymbol{x}^{\epsilon}_{0}=\boldsymbol{x}$, $\boldsymbol{v}^{\epsilon}_{0}=\boldsymbol{v}$, $\boldsymbol{z}^{\epsilon}_{0}=\boldsymbol{z}$ and $\boldsymbol{\zeta}^{\epsilon}_{0}=\boldsymbol{\zeta}$ are random variables. Note that $\boldsymbol{y}^{\epsilon}_{0}=\boldsymbol{0}$, and $\boldsymbol{\beta}^{\epsilon}_{0}$ is a zero mean Gaussian random variable with covariance $\boldsymbol{M}_{2}/\tau_{\xi}\epsilon$.

Let $\boldsymbol{W}^{(q_{2})}_{t}$ be an $\mathbb{R}^{q_{2}}$-valued Wiener process on the filtered probability space $(\Omega,\mathcal{F},\mathcal{F}_{t},\mathbb{P})$ satisfying the usual conditions karatzas2012Brownian and $\mathbb{E}$ denotes expectation with respect to $\mathbb{P}$.

We use the notation introduced in Section 3 and make the following assumptions.

Assumption 5**.**

For every $\boldsymbol{x}\in\mathbb{R}^{d}$, the vector-valued function $\boldsymbol{F}(\boldsymbol{x})$ is continuous, bounded and Lipschitz in $\boldsymbol{x}$, whereas the matrix-valued functions $\boldsymbol{g}(\boldsymbol{x})$, $\boldsymbol{h}(\boldsymbol{x})$, $\boldsymbol{\sigma}(\boldsymbol{x})$, $(\boldsymbol{g})_{\boldsymbol{x}}(\boldsymbol{x})$, $(\boldsymbol{h})_{\boldsymbol{x}}(\boldsymbol{x})$ and $(\boldsymbol{\sigma})_{\boldsymbol{x}}(\boldsymbol{x})$ are continuously differentiable, bounded and Lipschitz in $\boldsymbol{x}$. Moreover, $(\boldsymbol{g})_{\boldsymbol{x}\boldsymbol{x}}(\boldsymbol{x})$, $(\boldsymbol{h})_{\boldsymbol{x}\boldsymbol{x}}(\boldsymbol{x})$ and $(\boldsymbol{\sigma})_{\boldsymbol{x}\boldsymbol{x}}(\boldsymbol{x})$ are bounded for every $\boldsymbol{x}\in\mathbb{R}^{d}$.

Assumption 6**.**

The initial conditions $\boldsymbol{x}$, $\boldsymbol{v}$, $\boldsymbol{z}$, $\boldsymbol{\zeta}$ are random variables independent of $\epsilon$. We assume that they have finite moments of all orders, i.e. $\mathbb{E}|\boldsymbol{x}|^{p},\ \mathbb{E}|\boldsymbol{v}|^{p},\ \mathbb{E}|\boldsymbol{z}|^{p},\ \mathbb{E}|\boldsymbol{\zeta}|^{p}<\infty$ for all $p>0$.

Assumption 7**.**

There are no explosions, i.e. almost surely, for any value of the parameter $\epsilon$ there exists a unique solution on the compact time interval $[0,T]$ to the pre-limit equations (23)-(24), and also to the limiting equation (45).

The following convergence theorem is the main result of this paper. It provides a homogenized SDE for the particle’s position in the limit as the inertial time scale, the memory time scale and the noise correlation time scale go to zero at the same rate in the case when the pre-limit dynamics are described by the family of equations (23)-(24) (with $\mu=\theta=\nu=1$), or equivalently by the SDEs (37)-(42). In the following, $(\boldsymbol{D})_{ij}$ denotes the $(i,j)$-entry of the matrix $\boldsymbol{D}$.

Theorem 4.1.

Let $\boldsymbol{x}^{\epsilon}_{t}\in\mathbb{R}^{d}$ be the solution to the SDEs (37)-(42). Suppose that Assumptions 5-7 are satisfied and the effective damping and diffusion (constant) matrices, $\boldsymbol{K}_{1}$, $\boldsymbol{K}_{2}$, defined in (9) and (10) respectively, are invertible. Moreover, we assume that for every $\boldsymbol{x}\in\mathbb{R}^{d}$,

[TABLE]

is invertible for all $\lambda$ in the right half plane $\{\lambda\in\mathbb{C}:Re(\lambda)>0\}$, where $\tilde{\boldsymbol{\kappa}}(z):=\boldsymbol{C}_{1}(z\boldsymbol{I}+\boldsymbol{\Gamma}_{1})^{-1}\boldsymbol{M}_{1}\boldsymbol{C}_{1}^{*}$, for $z\in\mathbb{C}$, is the Laplace transform of the memory function.

Denote $\boldsymbol{\theta}(\boldsymbol{X}):=\boldsymbol{g}(\boldsymbol{X})\boldsymbol{K}_{1}\boldsymbol{h}(\boldsymbol{X})\in\mathbb{R}^{d\times d}$ for $\boldsymbol{X}\in\mathbb{R}^{d}$. Then as $\epsilon\to 0$, the process $\boldsymbol{x}^{\epsilon}_{t}$ converges to the solution, $\boldsymbol{X}_{t}$, of the following Itô SDE:

[TABLE]

with $\boldsymbol{S}(\boldsymbol{X}_{t})=\boldsymbol{S}^{(1)}(\boldsymbol{X}_{t})+\boldsymbol{S}^{(2)}(\boldsymbol{X}_{t})+\boldsymbol{S}^{(3)}(\boldsymbol{X}_{t}),$ where the $\boldsymbol{S}^{(k)}$ are the noise-induced drifts whose $i$th components are given by

[TABLE]

Here $\boldsymbol{J}_{11}=\boldsymbol{J}_{11}^{*}\in\mathbb{R}^{d\times d}$, $\boldsymbol{J}_{21}=\boldsymbol{J}_{12}^{*}\in\mathbb{R}^{d_{1}\times d}$ and $\boldsymbol{J}_{31}=\boldsymbol{J}_{13}^{*}\in\mathbb{R}^{d_{2}\times d}$ satisfy the following system of five matrix equations:

[TABLE]

The convergence is obtained in the following sense: for all finite $T>0$,

[TABLE]

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

Remark 6.

Invertibility of the matrices $\boldsymbol{B}_{\lambda}(\boldsymbol{x})$ (the assumption (44)) is a technical condition which will be used in the proof of the theorem. We are going to verify it in the special cases and applications discussed later (see Corollary 3 and Corollary 5). In particular, it follows from the stronger spectral condition, namely that $\boldsymbol{g}(\boldsymbol{x})\tilde{\boldsymbol{\kappa}}(\mu)\boldsymbol{h}(\boldsymbol{x})$ has no spectrum in the right half plane for any $\mu$ with $Re(\mu)>0$. See also Remark 7.

Proof.

We denote $\boldsymbol{\hat{x}}^{\epsilon}_{t}:=(\boldsymbol{x}^{\epsilon}_{t},\boldsymbol{z}^{\epsilon}_{t},\boldsymbol{\zeta}^{\epsilon}_{t})\in\mathbb{R}^{d+d_{1}+d_{2}}$ and $\boldsymbol{\hat{v}}^{\epsilon}_{t}:=(\boldsymbol{v}^{\epsilon}_{t},\boldsymbol{y}^{\epsilon}_{t},\boldsymbol{\eta}^{\epsilon}_{t})\in\mathbb{R}^{d+d_{1}+d_{2}}$ and rewrite the above SDE system in the form $\eqref{gsk1}$-$\eqref{gsk2}$:

[TABLE]

with

[TABLE]

where $\boldsymbol{\hat{\gamma}}\in\mathbb{R}^{(d+d_{1}+d_{2})\times(d+d_{1}+d_{2})}$ is a 3 by 3 block matrix with each block a matrix of appropriate dimension; $\boldsymbol{\hat{F}}\in\mathbb{R}^{d+d_{1}+d_{2}}$, $\boldsymbol{\hat{\sigma}}\in\mathbb{R}^{(d+d_{1}+d_{2})\times q_{2}}$ and the $\boldsymbol{0}$ appearing in $\boldsymbol{\hat{\gamma}}$, $\boldsymbol{\hat{F}}$ and $\boldsymbol{\hat{\sigma}}$ is a zero vector or matrix of appropriate dimension.

We now want to apply Theorem 3.1 (with $m:=\epsilon$, $n:=d+d_{1}+d_{2}$, $\boldsymbol{\gamma}$ replaced by $\boldsymbol{\hat{\gamma}}$, $\boldsymbol{F}$ replaced by $\boldsymbol{\hat{F}}$, $\boldsymbol{\sigma}$ replaced by $\boldsymbol{\hat{\sigma}}$, etc.) to $\eqref{gen_nsk1}$-$\eqref{gen_nsk2}$.

It is straightforward to see that Assumption 5 implies Assumption 1 and Assumption 7 implies Assumption 4.

To verify Assumption 3, note again that $\boldsymbol{y}^{\epsilon}_{0}=\boldsymbol{0}$ and so by Assumption 6, we only need to show that for every $p>0$, $\mathbb{E}|\epsilon\boldsymbol{\beta}^{\epsilon}_{0}|^{p}=O(\epsilon^{\alpha})$ as $\epsilon\to 0$, for some constant $\alpha\geq p/2$. To show this, we use the fact that for a mean zero Gaussian random variable, $X\in\mathbb{R}$, with variance $\sigma^{2}$,

[TABLE]

for every $p>0$, where $\Gamma$ denotes the gamma function winkelbauer2012moments . Applying this to $\boldsymbol{\beta}^{\epsilon}_{0}$, we obtain, for every $p>0$, $\mathbb{E}|\boldsymbol{\beta}_{0}^{\epsilon}|^{p}=O(1/\epsilon^{p/2})$ as $\epsilon\to 0$ and so $\mathbb{E}|\epsilon\boldsymbol{\beta}_{0}^{\epsilon}|^{p}=O(\epsilon^{p/2})$ as $\epsilon\to 0$. Therefore, Assumption 3 is verified.

It remains to verify Assumption 2, i.e. that $\boldsymbol{\hat{\gamma}}$ is positive stable. Note that $\boldsymbol{\Gamma}_{2}$ is positive stable by assumption and the triangular-block structure of $\boldsymbol{\hat{\gamma}}$ implies that one only needs to verify that the upper left 2 by 2 block matrix of $\boldsymbol{\hat{\gamma}}$:

[TABLE]

is positive stable, where $\boldsymbol{x}\in\mathbb{R}^{d}$.

We thus need to show that the resolvent set of $-\boldsymbol{L(\boldsymbol{x})}$, $\rho(-\boldsymbol{L}(\boldsymbol{x})):=\{\lambda\in\mathbb{C}:(\lambda\boldsymbol{I}+\boldsymbol{L}(\boldsymbol{x}))^{-1}\text{ exists}\}$, contains the right half plane $\{\lambda\in\mathbb{C}:Re(\lambda)>0\}$ for every $\boldsymbol{x}\in\mathbb{R}^{d}$.

Let $\lambda\in\mathbb{C}$ such that $Re(\lambda)>0$. We will use the following formula for blockwise inversion of a block matrix: provided that $\boldsymbol{S}$ and $\boldsymbol{P}-\boldsymbol{Q}\boldsymbol{S}^{-1}\boldsymbol{R}$ are nonsingular, we have

[TABLE]

where $\boldsymbol{P}$, $\boldsymbol{Q}$, $\boldsymbol{R}$, $\boldsymbol{S}$ are matrix sub-blocks of arbitrary dimension.

Since the matrices $\boldsymbol{A}_{\lambda}:=\boldsymbol{\Gamma}_{1}/\tau_{\kappa}+\lambda\boldsymbol{I}$ and $\boldsymbol{B}_{\lambda}(\boldsymbol{x}):=\boldsymbol{I}+\boldsymbol{g}(\boldsymbol{x})\tilde{\boldsymbol{\kappa}}(\lambda\tau_{\kappa})\boldsymbol{h}(\boldsymbol{x})/\lambda m_{0}$ are invertible for all $\lambda$ in the right half plane by assumption, $\lambda\boldsymbol{I}+\boldsymbol{L}(\boldsymbol{x})$ is indeed invertible for every $\boldsymbol{x}$ and in fact using the above formula for the inverse of a block matrix, we have:

[TABLE]

Therefore, $\boldsymbol{\hat{\gamma}}$ is invertible and one can compute:

[TABLE]

where $\boldsymbol{\theta}=\boldsymbol{g}\boldsymbol{K}_{1}\boldsymbol{h}$. The result follows by applying Theorem 3.1 to the SDE systems (55)-(56). In particular, a rewriting of the resulting Lyapunov equation (35) gives the system of matrix equations (49)-(53). ∎

Remark 7.

In the above proof, the condition of invertibility of $\boldsymbol{B}_{\lambda}(\boldsymbol{x})$ is only used to guarantee the positive stability of the matrix $\hat{\boldsymbol{\gamma}}$. Therefore, the conclusion of the theorem holds also when the latter can be established in another way. This can indeed be done in a number of concrete examples.

Remark 8.

Our SIDEs belong to a special class of non-Markovian equations, the so-called quasi-Markovian Langevin equations eckmann1999non . For these equations one can introduce a finite number of auxiliary variables in such a way that the evolution of particle’s position and velocity, together with these auxiliary variables, is described by a usual SDE system and is thus Markovian. We remark that such “Markovianization” procedure works here because the colored noise can be generated by a linear system of SDEs and the memory kernel satisfies a linear system of ordinary differential equations—both with constant coefficients. If, on the other hand, the memory kernel decays as a power, then there is no finite dimensional extension of the space which would make the solution process Markovian luczka2005non .

The following corollary uses a linear change of variables in a given SIDE, to arrive at an alternative form of the corresponding homogenized SDE of the form (45).

Corollary 1.

For $i=1,2$, let $\boldsymbol{T}_{i}$ be arbitrary $d_{i}\times d_{i}$ constant invertible matrix, where $d_{1},d_{2}$ are positive integers. For $t\geq 0$, denote $\boldsymbol{\Gamma}^{\prime}_{i}=\boldsymbol{T}_{i}\boldsymbol{\Gamma}_{i}\boldsymbol{T}^{-1}_{i}$, $\boldsymbol{M}_{i}^{\prime}=\boldsymbol{T}_{i}\boldsymbol{M}_{i}\boldsymbol{T}_{i}^{*}$, $\boldsymbol{C}^{\prime}_{i}=\boldsymbol{C}_{i}\boldsymbol{T}_{i}^{-1}$, $(\boldsymbol{\beta}^{\epsilon}_{t})^{\prime}=\boldsymbol{T}_{2}\boldsymbol{\beta}^{\epsilon}_{t}$, $\boldsymbol{\Sigma}_{i}^{\prime}=\boldsymbol{T}_{i}\boldsymbol{\Sigma}_{i}$ and consider the equations:

[TABLE]

where $\boldsymbol{W}^{\prime}_{t}$ is a $q_{2}$-dimensional Wiener process and the initial condition, $(\boldsymbol{\beta}^{\epsilon}_{0})^{\prime}$, is normally distributed with zero mean and covariance of $\boldsymbol{M}^{\prime}_{2}/(\tau_{\xi}\epsilon^{\nu})$.

Suppose that Assumptions 5-7 are satisfied and the effective damping and diffusion constants, $\boldsymbol{K}^{\prime}_{i}=\boldsymbol{C}^{\prime}_{i}(\boldsymbol{\Gamma}_{i}^{\prime})^{-1}\boldsymbol{M}_{i}^{\prime}(\boldsymbol{C}_{i}^{\prime})^{*}=\boldsymbol{K}_{i}$ ($i=1,2$), are invertible. Moreover, we assume that $\boldsymbol{I}+\boldsymbol{g}(\boldsymbol{x})\tilde{\boldsymbol{\kappa}^{\prime}}(\lambda\tau_{\kappa})\boldsymbol{h}(\boldsymbol{x})/\lambda m_{0}$ is invertible for all $\lambda$ in the right half plane $\{\lambda\in\mathbb{C}:Re(\lambda)>0\}$ and $\boldsymbol{x}\in\mathbb{R}^{d}$, where $\tilde{\boldsymbol{\kappa}^{\prime}}(z):=\boldsymbol{C}^{\prime}_{1}(z\boldsymbol{I}+\boldsymbol{\Gamma}^{\prime}_{1})^{-1}\boldsymbol{M}^{\prime}_{1}(\boldsymbol{C}^{\prime}_{1})^{*}=\boldsymbol{\tilde{\kappa}}(z)$ for $z\in\mathbb{C}$.

Let $\mu=\theta=\nu$ in (63)-(64). Then as $\epsilon\to 0$, the process $\boldsymbol{x}^{\epsilon}_{t}$ converges, in the similar sense as in Theorem 4.1, to $\boldsymbol{X}_{t}$, where $\boldsymbol{X}_{t}$ is the solution of the SDE (45) with the $\boldsymbol{C}_{i}$, $\boldsymbol{\Gamma}_{i}$, $\boldsymbol{M}_{i}$, $\boldsymbol{\Sigma}_{i}$ replaced by $\boldsymbol{C}^{\prime}_{i}$, $\boldsymbol{\Gamma}^{\prime}_{i}$, $\boldsymbol{M}^{\prime}_{i}$, $\boldsymbol{\Sigma}^{\prime}_{i}$ respectively, and the driving Wiener process $\boldsymbol{W}^{(q_{2})}_{t}$ replaced by $\boldsymbol{W}_{t}^{\prime}$.

Corollary 1 is an easy consequence of Theorem 4.1.

Next, we discuss a particular, but very important, case when a fluctuation-dissipation relation holds. This is, for instance, the case when the pre-limit dynamics are (heuristically) derived from Hamiltonian dynamics (see Appendix A). We will further explore similar cases of fluctuation-dissipation relations for the two sub-classes.

Corollary 2.

Let $\boldsymbol{x}^{\epsilon}_{t}\in\mathbb{R}^{d}$ be the solution to the SDEs (37)-(42). Suppose that the assumptions of Theorem 4.1 holds. Moreover, we assume that:

[TABLE]

where $\tau$ is a positive constant, and

[TABLE]

(so that $q=r$ and $d_{1}=d_{2}$). Denote $\boldsymbol{K}:=\boldsymbol{C}\boldsymbol{\Gamma}^{-1}\boldsymbol{M}\boldsymbol{C}^{*}$. Then as $\epsilon\to 0$, the process $\boldsymbol{x}^{\epsilon}_{t}$ converges to the solution, $\boldsymbol{X}_{t}$, of the following Itô SDE:

[TABLE]

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

[TABLE]

where $\boldsymbol{J}_{11}$ solves the following system of three matrix equations:

[TABLE]

The convergence is obtained in the same sense as in Theorem 4.1.

Eqns. (65)-(66) are a form of fluctuation-dissipation relation familiar from non-equilibrium statistical mechanics toda2012statistical . As stationary measures of systems satisfying fluctuation-dissipation relations are in equilibrium with respect to the underlying dynamics, this result is relevant for describing equilibrium properties of such systems in the small mass limit.

Remark 9.

Therefore, if the fluctuation-dissipation relation holds, the noise-induced drift in the limiting SDE reduces to a single term (later we will see how this term simplifies in some special cases). This result may have interesting implications for nanoscale systems in equilibrium. We remark that the conditions for the fluctuation-dissipation relation in Corollary 2 can be written in other equivalent forms, up to the transformations in (8) and multiplication by a constant.

Proof.

The above corollary follows from applying Theorem 4.1. Indeed, by assumptions of the corollary, (49) simplifies to:

[TABLE]

This implies that $\boldsymbol{J}_{12}=\boldsymbol{J}_{13}$ and, therefore, $\boldsymbol{S}^{(2)}$ and $\boldsymbol{S}^{(3)}$ cancel. Rewriting the resulting system of matrix equations in (49)-(53) give (69)-(71). ∎

5 Homogenization for Models of the Two Sub-Classes

We now return to the two sub-classes of SIDEs (11) introduced in Section 2. In this section, we study the effective dynamics described by SIDEs (25) and (27) in the limit as $\epsilon\to 0$. By specializing to these two sub-classes, the convergence result of Theorem 4.1, in particular the expressions in (45)-(53), can be made more explicit under certain assumptions on the matrix-valued coefficients and therefore the limiting equation obtained may be useful for modeling purposes.

5.1 SIDEs Driven by a Markovian Colored Noise

The following convergence result provides a homogenized SDE for the particle’s position in the limit as the inertial time scale, the memory time scale and the noise correlation time scale vanish at the same rate in the case when the pre-limit dynamics are driven by an Ornstein-Uhlenbeck noise.

Corollary 3.

Let $d=d_{1}=d_{2}=q_{1}=q_{2}=q=r$. We set, in the SDEs (37)-(42): $\boldsymbol{\beta}^{\epsilon}_{t}=\boldsymbol{\eta}^{\epsilon}_{t}$, $\tau_{\xi}=\tau_{\eta}$, $\boldsymbol{W}_{t}^{(q_{2})}=\boldsymbol{W}^{(d)}_{t}:=\boldsymbol{W}_{t}$ and

[TABLE]

to obtain SDEs equivalent to equations (25)-(26) with $\mu=\theta=\nu=1$. Let $\boldsymbol{x}^{\epsilon}_{t}\in\mathbb{R}^{d}$ be the solution to these equations, with the matrices $\boldsymbol{g}(\boldsymbol{x})$ and $\boldsymbol{h}(\boldsymbol{x})$ positive definite for every $\boldsymbol{x}\in\mathbb{R}^{d}$. Suppose that Assumptions 5-7 are satisfied and, moreover, that $\boldsymbol{g}(\boldsymbol{x})$, $\boldsymbol{h}(\boldsymbol{x})$ and the diagonal matrix $\boldsymbol{A}$ are commuting.

Then as $\epsilon\to 0$, the process $\boldsymbol{x}^{\epsilon}_{t}$ converges to the solution, $\boldsymbol{X}_{t}$, of the following Itô SDE:

[TABLE]

with $\boldsymbol{S}(\boldsymbol{X}_{t})=\boldsymbol{S}^{(1)}(\boldsymbol{X}_{t})+\boldsymbol{S}^{(2)}(\boldsymbol{X}_{t})+\boldsymbol{S}^{(3)}(\boldsymbol{X}_{t}),$ where the $\boldsymbol{S}^{(k)}$ are the noise-induced drifts whose $i$th components are given by

[TABLE]

Here $\boldsymbol{J}_{11}=\boldsymbol{J}_{11}^{*}$, $\boldsymbol{J}_{21}=\boldsymbol{J}_{12}^{*}$ and $\boldsymbol{J}_{31}=\boldsymbol{J}_{13}^{*}$ are $d$ by $d$ block matrices satisfying the following system of matrix equations:

[TABLE]

The convergence is obtained in the same sense as in Theorem 4.1.

Proof.

We will apply Theorem 4.1. As $\boldsymbol{K}_{1}=2\boldsymbol{K}_{2}=\boldsymbol{I}$, clearly they are invertible. Also, being positive definite, $\boldsymbol{g}$ and $\boldsymbol{h}$ are invertible.

Since $\boldsymbol{g}$, $\boldsymbol{h}$ and $\boldsymbol{A}$ are positive definite and commuting matrices, the matrix $\boldsymbol{B}_{\lambda}(\boldsymbol{x})$, defined in (44), is invertible for all $\lambda$ such that $Re(\lambda)>0$. Indeed, in this case $\boldsymbol{B}_{\lambda}(\boldsymbol{x})=\boldsymbol{I}+\boldsymbol{g}(\boldsymbol{x})(\lambda\tau_{\kappa}\boldsymbol{I}+\boldsymbol{A})^{-1}\boldsymbol{A}\boldsymbol{h}(\boldsymbol{x})/\lambda m_{0}$. Since $\boldsymbol{g}$, $\boldsymbol{h}$ and $\boldsymbol{A}$ are positive definite and commuting, they have positive eigenvalues and can be simultaneously diagonalized. Therefore, all the eigenvalues of $\boldsymbol{B}_{\lambda}(\boldsymbol{x})$ are nonzero for every $\lambda$ with $Re(\lambda)>0$ and $\boldsymbol{x}\in\mathbb{R}^{d}$, so the invertibility condition is verified. Therefore, the block matrix:

[TABLE]

is positive stable (see Remark 7). The result then follows by applying Theorem 4.1. ∎

For special one-dimensional systems, the form of the limiting equation can be made even more explicit.

Corollary 4.

In the one-dimensional case, we drop the boldface and write $\boldsymbol{X}_{t}:=X_{t}\in\mathbb{R},\ \boldsymbol{g}(\boldsymbol{X}):=g(X),$ with $g:\mathbb{R}\to\mathbb{R}$, etc.. We assume that $h=g$ and $\boldsymbol{A}:=\alpha>0$ is a constant. The homogenized equation is given by:

[TABLE]

with $S(X_{t})=S^{(1)}(X_{t})+S^{(2)}(X_{t})+S^{(3)}(X_{t}),$ where the noise-induced drift terms $S^{(k)}(X_{t})$ have the following explicit expressions that depend on the parameters $m_{0},\tau_{\kappa}$ and $\tau_{\eta}$:

[TABLE]

where the prime ′ denotes derivative with respect to $X_{t}$.

Proof.

With $\boldsymbol{x}^{\epsilon}_{t}:=(x^{\epsilon}_{t},z^{\epsilon}_{t},\zeta^{\epsilon}_{t})\in\mathbb{R}^{3}$ and $\boldsymbol{v}^{\epsilon}_{t}:=(v^{\epsilon}_{t},y^{\epsilon}_{t},\eta^{\epsilon}_{t})\in\mathbb{R}^{3}$, SDEs (55)-(56) become:

[TABLE]

where

[TABLE]

It follows from Corollary 3 that the matrix $\boldsymbol{\gamma}$ is positive stable; one can also calculate its eigenvalues explicitly and see that their real parts are positive. The eigenvalues of $\boldsymbol{\gamma}$ are

[TABLE]

and so their real parts are indeed positive.

On the other hand, the solution, $\boldsymbol{J}\in\mathbb{R}^{3\times 3}$, to the Lyapunov equation,

[TABLE]

can be computed (using Mathematica®) to be:

[TABLE]

The result then follows from Corollary 3. ∎

Remark 10.

Note that here the matrix $\boldsymbol{\gamma}$ in (90) is not symmetric and the smallest eigenvalue of its symmetric part can be negative. Moreover, the initial condition $\boldsymbol{v}^{\epsilon}_{0}$ depends on $\epsilon$ through the component $\eta^{\epsilon}_{0}$ (which is a zero mean Gaussian random variable with variance $\alpha/2\epsilon$). Thus, we cannot apply the main results in hottovy2015smoluchowski to obtain the convergence result. This is our main motivation to revisit the Smoluchowski-Kramers limit of SDEs in Section III under a weakened spectral assumption on the matrix $\boldsymbol{\gamma}$ (or $\boldsymbol{\hat{\gamma}}$ in the multidimensional case) and a relaxed assumption concerning the $\epsilon$ dependence of $\boldsymbol{v}^{\epsilon}_{0}$ (or $\hat{\boldsymbol{v}}^{\epsilon}_{0}$ in the multidimensional case).

Remark 11.

In the important case when the fluctuation-dissipation relation (i.e. $\tau_{\kappa}=\tau_{\eta}$, $h=g$ and $g$ is proportional to $\sigma$) holds for the one-dimensional models of the first sub-class, the correction drift terms $S^{(2)}$ and $S^{(3)}$ cancel each other and the resulting (single) noise-induced drift term coincides with that obtained in the limit as $m\to 0$ of the systems with no memory, driven by a white noise to which Theorem 3.1 applies directly! However, when the relation fails, we obtain three different drift corrections induced by vanishing of all time scales. Again, the presence of these correction terms may have significant consequences for the dynamics of the systems (see Section 6).

5.2 SIDEs Driven by a Non-Markovian Colored Noise

The following corollary provides a homogenized SDE for the particle position in the limit, in which the inertial time scale, the memory time scale and the noise correlation time scale vanish at the same rate in the case when the pre-limit dynamics are driven by the harmonic noise. We emphasize that in this case the original system is driven by a noise which is not a Markov process.

Corollary 5.

Let $d=d_{1}=d_{2}=q_{1}=q_{2}=q=r$. We set, in the SDEs (37)-(42): $\tau_{\xi}=\tau_{h}$, $\boldsymbol{W}_{t}^{(q_{2})}=\boldsymbol{W}^{(d)}_{t}=\boldsymbol{W}_{t}$ and

[TABLE]

[TABLE]

[TABLE]

to obtain SDEs equivalent to equations (27)-(29) with $\mu=\theta=\nu=1$. Let $\boldsymbol{x}^{\epsilon}_{t}\in\mathbb{R}^{d}$ be the solution to the SDEs (37)-(42), with the matrices $\boldsymbol{g}(\boldsymbol{x})$ and $\boldsymbol{h}(\boldsymbol{x})$ positive definite for every $\boldsymbol{x}\in\mathbb{R}^{d}$. Moreover, $\boldsymbol{g}(\boldsymbol{x})$, $\boldsymbol{h}(\boldsymbol{x})$ and the diagonal matrix $\boldsymbol{\Omega}^{2}$ are commuting. Suppose that Assumptions 5-7 are satisfied.

Then as $\epsilon\to 0$, the process $\boldsymbol{x}^{\epsilon}_{t}$ converges to the solution, $\boldsymbol{X}_{t}$, of the following Itô SDE

[TABLE]

with $\boldsymbol{S}(\boldsymbol{X}_{t})=\boldsymbol{S}^{(1)}(\boldsymbol{X}_{t})+\boldsymbol{S}^{(2)}(\boldsymbol{X}_{t})+\boldsymbol{S}^{(3)}(\boldsymbol{X}_{t}),$ where the $\boldsymbol{S}^{(k)}$ are the noise-induced drift terms whose $i$th components are given by the expressions

[TABLE]

where $i,j,l=1,\dots,d$. In the above,

[TABLE]

is the block matrix solving the Lyapunov equation

[TABLE]

where

[TABLE]

In the above, $\boldsymbol{\hat{\gamma}}\in\mathbb{R}^{5d\times 5d}$ is a 5 by 5 block matrix with each block an $\mathbb{R}^{d\times d}$-valued matrix, $\boldsymbol{\hat{\sigma}}\in\mathbb{R}^{5d\times d}$ is a 5 by 1 block matrix with each block a $\mathbb{R}^{d\times d}$-valued matrix, $\boldsymbol{I}$ is a $d\times d$ identity matrix, $\boldsymbol{0}$ in $\boldsymbol{\hat{\gamma}}$ and $\boldsymbol{\hat{\sigma}}$ is a $d\times d$ zero matrix, and $\boldsymbol{W}$ is a $d$-dimensional Wiener process. The convergence is obtained in the same sense as in Theorem 4.1.

Note that the oscillatory nature of covariance function of the harmonic noise in the pre-limit SIDE makes the noise-induced drift in the resulting limiting SDE more complicated (there are more terms) compared to the case of OU process in the first sub-class. Therefore, we write the system of matrix equations that the $\boldsymbol{J}_{kl}$ satisfy in the form of a matrix Lyapunov equation in Corollary 5, without breaking it up into equations for individual blocks. This could of course be done, leading to a (more complicated) analog of (78)-(82). The proof of Corollary 5 is essentially identical to the proof of Corollary 3, so we omit it.

Again, for special one-dimensional systems, we are going to make the result more explicit.

Corollary 6.

In the one-dimensional case, we drop the boldface and write $\boldsymbol{X}_{t}:=X_{t}\in\mathbb{R},\ \boldsymbol{g}(\boldsymbol{x}):=g(x),$ with $g:\mathbb{R}\to\mathbb{R}$, etc.. We assume that $h=g$ and $\boldsymbol{\Omega}:=\Omega$ is a real constant. The homogenized equation is given by:

[TABLE]

with $S(X_{t})=S^{(1)}(X_{t})+S^{(2)}(X_{t})+S^{(3)}(X_{t}),$ where the noise-induced drift terms $S^{(k)}(X_{t})$ have the following explicit expressions (computed using Mathematica®) that depend on the parameters $m_{0},\tau_{\kappa}$ and $\tau_{h}$:

[TABLE]

where the prime ′ denotes derivative with respect to $X$ and the $J_{kl}(X)$ are given by:

[TABLE]

where $g=g(X)$, $\sigma=\sigma(X)$ and

[TABLE]

Note that if we send $\Omega\to\infty$ in the expressions for the $S^{(i)}(X)$ $(i=1,2,3)$ above, we recover the corresponding expressions given in Corollary 4 (with $\alpha=1$). This is not surprising, since in this limit the harmonic noise becomes an OU process (with $\alpha=1)$.

Moreover, when $\tau_{\kappa}=\tau_{h}=\tau$, the noise-induced drift becomes $S(X)=S^{(1)}(X)+S^{(2)}(X)+S^{(3)}(X),$ where

[TABLE]

Again, in the case when the fluctuation-dissipation relation holds we see that the noise-induced drift coincides with that obtained in the limit as $m\to 0$ of the Markovian model in Section III.

6 Application to the Study of Thermophoresis

6.1 Introduction

We revisit the dynamics of a free Brownian particle immersed in a heat bath where a temperature gradient is present. This was previously studied in Hottovy2012a . It was found there that the particle experiences a drift in response to the temperature gradient, due to the interplay between the inertial time scale and the noise correlation time scale. Such phenomenon is called thermophoresis. We refer to Hottovy2012a ; piazza2008thermophoresis and the references therein for further descriptions of this phenomenon, including references to experiments.

Here, we will study the dynamics of the particle in a non-equilibrium heat bath, where a generalized fluctuation-dissipation relation holds, in which both the diffusion coefficient and the temperature of the heat bath vary with the position. In contrast to Hottovy2012a , we take into account also the memory time scale (in addition to the inertial time scale and the noise correlation time scale) and model the position of the particle as the solution to a SIDE of the form (11). Unlike the model used in Hottovy2012a , the model can be derived heuristically from microscopic dynamics by an argument very similar to that of Appendix A.

For a spherical particle of radius $R$ immersed in a fluid of viscosity $\mu$, which in general is a function of the temperature $T=T(x)$ (and thus depends on $x$ as well), the friction (or damping) coefficient $\gamma$ satisfies the Stokes law toda2012statistical :

[TABLE]

On the other hand, the damping coefficient $\gamma(x)$ and the noise coefficient $\sigma(x)$ are expressed in terms of the diffusion coefficient $D(x)$ and the temperature $T(x)$ as follows:

[TABLE]

In the following, we study two one-dimensional non-Markovian models of thermophoresis. The first model is driven by a Markovian colored noise and the second model by a non-Markovian one.

6.2 A Thermophoresis Model with Ornstein-Uhlenbeck Noise

In this section we model evolution of the position, $x_{t}\in\mathbb{R}$, of a particle by the following SIDE:

[TABLE]

where $\eta_{t}$ is a stationary process, satisfying the SDE:

[TABLE]

The above equations are obtained by setting $d=1$, $\boldsymbol{F}=0$, $\boldsymbol{h}=\boldsymbol{g}=g:=\sqrt{\gamma}$, $\boldsymbol{\sigma}=\sigma$ in $\eqref{side2}$ and $\boldsymbol{A}=\alpha$ in $\eqref{ou}$, where $\gamma$ and $\sigma$ are given by (118). Note that the noise correlation function is proportional to the memory kernel in the SIDE (119), i.e.

[TABLE]

as in (21). Together with (118), this implies that (119) satisfies the generalized fluctuation-dissipation relation (see the statement of Corollary 2 and Remark 11). Note also that $g$ is a constant multiple of $\sigma$ if and only if $T$ is position-independent.

We now consider the effective dynamics of the particle in the limit when all the three characteristic time scales vanish at the same rate. In the following, the prime ′ denotes derivative with respect to the argument of the function.

Corollary 7.

Let $\epsilon>0$ be a small parameter and let the particle’s position, $x^{\epsilon}_{t}\in\mathbb{R}$ ($t\geq 0$), satisfy the following rescaled version of (119)-(120):

[TABLE]

where $m_{0}$, $\alpha$, $\tau$ are positive constants, and $(W_{t})$ is a one-dimensional Wiener process. The initial conditions are random variables $x^{\epsilon}_{0}=x$, $v^{\epsilon}_{0}=v$, independent of $\epsilon$ and (statistically) independent of $(W_{t})$, and $\eta^{\epsilon}_{0}$ is distributed according to the invariant distribution of the SDE (124).

Assume that the assumptions of Corollary 3 are satisfied (in particular, $\gamma(x)>0$ for every $x\in\mathbb{R}$). Then, in the limit as $\epsilon\to 0$, $x^{\epsilon}_{t}$ converges (in the same sense as in Corollary 3) to the process $X_{t}\in\mathbb{R}$, satisfying the SDE:

[TABLE]

with the noise-induced drift, $b_{1}(X)=S^{(1)}(X)+S^{(2)}(X)+S^{(3)}(X)$, where

[TABLE]

Proof.

The corollary follows from Corollary 4. In particular, the expressions for $S^{(1)}$, $S^{(2)}$ and $S^{(3)}$ follow from applying Corollary 4 to the present system (see (85)-(87)). ∎

We give some remarks and discussions of the contents of Corollary 7 before we end this subsection.

Remark 12.

We see that in this case a part of $S^{(2)}$ cancels $S^{(3)}$ and therefore the noise-induced drift simplifies to:

[TABLE]

Using the Stokes law (117) which gives

[TABLE]

where $\mu(T)=\mu(T(X))$, we have

[TABLE]

Equation (129) gives the thermophoretic drift in the limit when the three characteristic time scales vanish. Since it arises in the absence of an external force acting on the particle, it is a “spurious drift” caused by the presence of the temperature gradient and the state-dependence of the diffusion coefficient. Compared to eqn. (101) in hottovy2015smoluchowski , the drift term derived here contains a correction term due to the temperature profile.

Discussion. We discuss some physical implications of the thermophoretic drift given in $\eqref{thermodrift}$. As discussed in Hottovy2012a , the sign of $b_{1}(X)$ determines the direction in which the particle is expected to travel. The particle will eventually reach some boundaries, which can be either absorbing or reflecting. We are going to consider the reflecting boundaries case. The position of the particle reaches a steady-state distribution $\rho_{\infty}(X)$ in the limit $t\to\infty$. Assuming that the particle is confined to the interval $(a,b)$, $a<b$, one can compute the stationary density:

[TABLE]

where in terms of the original parameters of the model, $r:=\tau/m_{0}>0$, and $C$ is a normalizing constant. In particular, in absence of temperature gradient ($T^{\prime}(y)=0$), the particle is equally likely to be found anywhere in $(a,b)$, whereas when a temperature gradient is present, the distribution of the particle’s position is not uniform. In the limit $r\to\infty$, the particle’s position is again distributed uniformly on $(a,b)$. On the other hand, in the limit $r\to 0$ the stationary density is inversely proportional to the temperature, i.e. $\rho_{\infty}(X)=\tilde{C}T(X)^{-1},$ where $\tilde{C}$ is a normalizing constant. Thus, the particle is more likely to be found in the colder region. In the special case when $D(X)$ is proportional to $T(X)$, so that $\gamma$ is independent of $X$, we have

[TABLE]

where $\tilde{C}$ is a normalizing constant, so the particle is more likely to be found in the colder region, with the likelihood decreasing as $r$ increases.

Next, we are going to study the sign of the thermophoretic drift directly using (129) (this is in contrast to the approach in Hottovy2012a , where $\mu(T)$ is expanded around a fixed temperature). We find that $b_{1}(X)>0$ if and only if $r>r_{c}$ and $r_{c}$ is the critical ratio of $\tau/m_{0}$, given by:

[TABLE]

where $\mu(T)=\mu(T(X))$ is obtained from the Stokes law. For $r=r_{c}$, the stationary density (132) reduces to:

[TABLE]

where $C$ is a normalizing constant. Importantly, note that the drift does not change sign if $T$ is independent of $X$.

Finally, we discuss a special case. When $\mu(T)=\mu_{0}>0$ is a constant (so that $\gamma(X)$ is a constant), the thermophoretic drift is given by:

[TABLE]

In agreement with the result in Hottovy2012a , $b_{1}(X)$ has the same sign as $T^{\prime}(X)$, leading to a flow towards the hotter region. The steady-state density is

[TABLE]

where $C$ is a normalizing constant, and the particle is more likely to be found in the colder region for all $r>0$, even though the thermophoretic drift actually directs the particle towards the hotter regions. This effect is in agreement with experiments, and is explained by the presence of reflecting boundary conditions.

6.3 A Thermophoresis Model with Non-Markovian (Harmonic) Noise

We repeat the analysis of the previous subsection in the case when the colored noise is a harmonic noise. We set $d=1$, $\boldsymbol{F}=0$, $\boldsymbol{h}=\boldsymbol{g}=g:=\sqrt{\gamma}$, $\boldsymbol{\sigma}=\sigma$, $\boldsymbol{\Omega}=\Omega$, $\boldsymbol{\Omega}_{0}=\Omega_{0}:=\Omega\sqrt{1-\Omega^{2}/4}$, $\boldsymbol{\Omega}_{1}=\Omega_{1}:=\Omega/\sqrt{1-\Omega^{2}/4}$ (where $|\Omega|<2$) in the SIDE $\eqref{side3}$ and study the effective dynamics of the resulting system as before. The case where $|\Omega|>2$ can be studied analogously. The following result follows from Corollary 6.

Corollary 8.

Let $\epsilon>0$ be a small parameter and the particle’s position, $x^{\epsilon}_{t}\in\mathbb{R}\ (t\geq 0)$, satisfy the following rescaled SDEs:

[TABLE]

where $m_{0}$ and $\tau$ are positive constants, $\Omega$, $\Omega_{0}$ and $\Omega_{1}$ are constants defined as before, and $(W_{t})$ is a one-dimensional Wiener process. The initial conditions are given by the random variables $x^{\epsilon}_{0}=x$, $v^{\epsilon}_{0}=v$, independent of $\epsilon$, and $(h^{\epsilon}_{0},u^{\epsilon}_{0})$ are distributed according to the invariant measure of the SDEs (140)-(141).

Assume that the assumptions in Corollary 5 are satisfied. Then, in the limit as $\epsilon\to 0$, the process $x^{\epsilon}_{t}$ converges (in the same sense as Corollary 5) to the process $X_{t}\in\mathbb{R}$, satisfying the SDE:

[TABLE]

where the noise-induced drift term is given by:

[TABLE]

We next discuss the contents of Corollary 8.

Remark 13.

Note that $b_{2}(X)$ differs from $b_{1}(X)$ obtained previously and $b_{2}(X)\to b_{1}(X)$, with $\alpha=1$ in the expression for $b_{1}(X)$, in the limit $\Omega\to\infty$.

Discussion. In the reflecting boundaries case, the stationary distribution of the particle’s position is

[TABLE]

where $r:=\tau/m_{0}>0$ and $C$ is a normalizing constant. Similarly to the previous model, in the absence of temperature gradient (i.e. when $T$ is a constant), the particle is equally likely to be found anywhere in $(a,b)$. When a temperature gradient is present, distribution of the particle’s position is not uniform. However, in contrast to the previous model, in the limit $r\to\infty$ the particle is not distributed uniformly on $(a,b)$ and in the limit $r\to 0$ the stationary density is no longer inversely proportional to the temperature. Both distributions depend on the diffusion coefficient $D(X)$ as well as on the temperature profile $T(X)$.

We can also study the sign of the thermophoretic drift. In this case there can be up to two critical ratios, $r_{c}$, at which $b_{2}(X)$ changes sign, as the equation $b_{2}(X)=0$ is a quadratic equation in $r$. In the special case when $\mu(T)=\mu_{0}>0$ is a constant (and thus so is $\gamma(X)$), the thermophoretic drift is given by:

[TABLE]

In contrast to the result in previous model, $b_{2}(X)$ has the same sign as $T^{\prime}(X)$ provided that

[TABLE]

Thus, $b_{2}(X)$ and $T^{\prime}(X)$ do not share the same sign for all $r>0$, unless $|\Omega|\geq 1.$ Thus, according to this model, presence of a temperature gradient allows us to tune the parameters $(m_{0},\tau,\Omega)$ to control the direction which the particle travels. The steady-state density in this case is

[TABLE]

where $C$ is a normalizing constant. The particle will be more likely found in the colder region for all $r>0$ if $|\Omega|\geq 1$, whereas this might not be true for all $r>0$ if $|\Omega|<1$.

7 Conclusions and Final Remarks

We have studied homogenization of a class of GLEs in the limit when three characteristic time scales, i.e. the inertial time, the characteristic memory time in the damping term, and the correlation time of colored noise driving the equations, vanish at the same rate. We have derived effective equations, which are simpler in three respects:

The velocity variables have been homogenized. As a result, the number of degrees of freedom is reduced and there are no fast variables left.

• 2.

The equations become regular SDEs, since the memory time has been taken to zero.

• 3.

The system is driven by a white noise.

Importantly, noise-induced drifts are present in the limiting equations, resulting from the dependence of the coefficients of the original model on the state of the system. We have applied the general results to a study of thermophoretic drift, correcting the formulae obtained in an earlier work Hottovy2012a . In systems, satisfying a fluctuation-dissipation relation, the noise-induced drifts in the limiting SDEs for the particle’s position reduce to a single term, and for special cases the limiting SDEs coincide with that of hottovy2015smoluchowski . However, in the more general case, new terms appear, absent in the case without memory. To prove the main theorem, we have employed the main result of hottovy2015smoluchowski , proven here in a different version under a relaxed assumption on the damping matrix and the initial conditions.

Homogenization of other specific non-Markovian models can also be studied using the methods of this paper. An example is a system with exponentially decaying memory kernel, driven by white noise in the limit as the inertial and memory time scales vanish at the same rate. In this case the noise-induced drift in the limiting equation will consist of two terms, not three, as in the case studied here.

The colored noises considered in this paper have correlations decaying exponentially (short-range memory). It would be interesting to study cases where the GLE is driven by other colored noises such as fractional Gaussian noises, with covariances decaying as a power, relevant for modeling anomalous diffusion phenomena in fields ranging from biology to finance Kou (2011). As mentioned in Section 2, we will explore homogenization for GLEs with vanishing effective damping and diffusion constant in a future work.

Acknowledgements.

The authors were partially supported by NSF grant DMS-1615045. S. 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. The authors learned the method of introducing additional variables to eliminate the memory term from E. Vanden-Eijnden. They would like to thank Maciej Lewenstein for insightful discussion on the GLEs.

Appendices

Appendix A Derivation of SIDEs From a Hamiltonian Model

For completeness, we provide a derivation for a special case of SIDEs (11), (15) and (20) from a Hamiltonian model of a small system (Brownian particle) in contact with a heat bath in thermal equilibrium. The particle is moving in a potential $U$. The heat bath is modeled as a system of non-interacting harmonic oscillators whose initial energy is distributed according to the Gibbs distribution at temperature $T$. The Brownian particle is coupled to each oscillator in the bath. This model is used widely to study many systems in statistical physics (ford1965statistical, ; mori1965transport, ). Our goal is to derive, heuristically, a stochastic integro-differential equation (SIDE) for the position and momentum variables of the particle from the Hamiltonian dynamics. This derivation serves to motivate the class of SIDEs that we are studying in this paper. We emphasize that our derivation here is certainly not original and follows closely that in (hanggi1997generalized, ) (see also an abstract approach in (Zwanzig1973, )).

One approach to derive the equations is to assume first that there are finitely many harmonic oscillators in the bath (Kac-Zwanzig model (zwanzig2001nonequilibrium, ; ariel2008strong, )). We then take the thermodynamic limit by sending the number of oscillators to infinity in the resulting equations (replacing finite sum over oscillator frequencies by an integral), arguing that the set of frequencies must be dense to allow dissipation of energy from the system to the bath and to eliminate Poincaré recurrence. Another approach, which is more technically involved, is to replace the finite system of oscillator equations by a system modeled by a wave equation (rey2006open, ; pavliotis2014stochastic, ). We will derive the SIDEs by adopting the former approach in the multi-dimensional case.

We consider the situation where the coupling is nonlinear in the particle’s position and linear in the bath variables. Let $\boldsymbol{\hat{x}}=(\boldsymbol{x},\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N})\in\mathbb{R}^{d+Np}$ and $\boldsymbol{\hat{p}}=(\boldsymbol{p},\boldsymbol{p}_{1},\dots,\boldsymbol{p}_{N})\in\mathbb{R}^{d+Np}$ (here, $\boldsymbol{x},\boldsymbol{p}\in\mathbb{R}^{d}$ and $\boldsymbol{x}_{i},\boldsymbol{p}_{i}\in\mathbb{R}^{p}$ for $i=1,\dots,N$). Hereafter, the superscript ∗ denotes transposition and $|\boldsymbol{b}|^{2}:=\boldsymbol{b}^{*}\boldsymbol{b}=\sum_{k=1}^{n}b_{k}^{2}$ denotes square of the norm of vector $\boldsymbol{b}:=(b_{1},\dots,b_{n})\in\mathbb{R}^{n}$.

The Hamiltonian of the system plus bath is:

[TABLE]

where $m$ is the particle’s mass, $\boldsymbol{x}\in\mathbb{R}^{d}$ and $\boldsymbol{p}\in\mathbb{R}^{d}$ are respectively the particle’s position and momentum, $\boldsymbol{x}_{k}\in\mathbb{R}^{p}$, $\boldsymbol{p}_{k}\in\mathbb{R}^{p}$ and $\omega_{k}\in\mathbb{R}^{+}$ $(k=1,\dots,N)$ are respectively the position, momentum and frequency of the $k$th bath oscillator (with unit mass), $\boldsymbol{f}(\boldsymbol{x}):=(f_{1}(\boldsymbol{x}),\dots,f_{r}(\boldsymbol{x}))\in\mathbb{R}^{r}$ is a vector function of $\boldsymbol{x}:=(x^{(1)},\dots,x^{(d)})$ and the $\boldsymbol{c}_{k}\in\mathbb{R}^{r\times p}$ (so $\boldsymbol{c}_{k}^{*}\in\mathbb{R}^{p\times r}$) are coupling matrices that specify the coupling strength between the system and the $k$th bath oscillator.

To derive an equation for the particle’s position and momentum, we start by substituting the expression for $H(\hat{\boldsymbol{\boldsymbol{x}}},\hat{\boldsymbol{\boldsymbol{p}}})$ into the Hamilton’s equations to obtain:

[TABLE]

where $\boldsymbol{g}(\boldsymbol{x})\in\mathbb{R}^{d\times r}$ denotes the Jacobian matrix $\left(\frac{\partial f_{i}}{\partial x^{(j)}}\right)_{ij}.$

Next, we eliminate the bath variables $\boldsymbol{x}_{k},\boldsymbol{p}_{k}$, $k=1,\dots,N$, from the system’s dynamics. Solving for $\boldsymbol{x}_{k}(t)$ in terms of $\boldsymbol{x}(t)$:

[TABLE]

Substituting this into $\eqref{momeq}$, we obtain:

[TABLE]

where

[TABLE]

In the integral term above, we integrate by parts to obtain:

[TABLE]

Using this, the equation for $\boldsymbol{p}(t)$ becomes the generalized Langevin equation (GLE):

[TABLE]

where

[TABLE]

and

[TABLE]

Note that $\boldsymbol{\xi}(t)\in\mathbb{R}^{r}$ is expressed in terms of the initial values of the variables $\boldsymbol{x}^{\prime}_{k}(0):=\boldsymbol{x}_{k}(0)-\frac{\boldsymbol{c}^{*}_{k}}{\omega_{k}^{2}}\boldsymbol{f}(\boldsymbol{x}(0))\in\mathbb{R}^{p}$ and $\boldsymbol{p}_{k}(0)\in\mathbb{R}^{p}$. If all these initial values are known, then $\boldsymbol{\xi}(t)$ is a deterministic force. However, one rarely has a complete information about these initial values and this is where the introduction of randomness can help to simplify the model. In view of this, we assume that the variables $\boldsymbol{x}^{\prime}_{k}(0)$ and $\boldsymbol{p}_{k}(0)$ are random and are distributed according to a Gibbs measure, with the density:

[TABLE]

where $\beta=1/(k_{B}T)$ and $Z$ is the partition function. Taking the averages of the bath variables with respect to the above density:

[TABLE]

where $E_{\rho}$ denotes mathematical expectation with respect to $\rho$ and $\boldsymbol{I}\in\mathbb{R}^{p\times p}$ is identity matrix.

Note that $\boldsymbol{\xi}(t)$ is a stationary, Gaussian process, if it is conditionally averaged with respect to $\rho$ zwanzig2001nonequilibrium . It follows from this distribution of the bath variables that we have the fluctuation-dissipation relation:

[TABLE]

where $\boldsymbol{\kappa}(t-s)$ is the memory kernel whose formula is given in (159). Later, we will generalize the resulting covariance of the process $\boldsymbol{\xi}(t)$ to an integral expression. We remark that the memory function $\boldsymbol{\kappa}(t)$ and the “color” of the noise $\boldsymbol{\xi}(t)$ are determined by the bath spectrum and the system-bath coupling.

Now we pass to the continuum limit by replacing the sum over $k$ in $\boldsymbol{\kappa}(t)$ by an integral $\int_{\mathbb{R}^{+}}d\omega n(\omega)$, where $n(\omega)$ is a density of states. Then, if the $\boldsymbol{c}_{k}$ are replaced by $\boldsymbol{c}(\omega)\in\mathbb{R}^{r\times p}$, the memory function $\boldsymbol{\kappa}(t)$ becomes the function:

[TABLE]

where $\hat{\boldsymbol{\kappa}}_{c}(\omega):=n(\omega)\boldsymbol{c}(\omega)\boldsymbol{c}(\omega)^{*}/\omega^{2}\in L^{1}(\mathbb{R}^{+})$. Repeating the same procedure for the noise process and also replacing the $\boldsymbol{x}^{\prime}_{k}(0)$ and $\boldsymbol{p}_{k}(0)$ by $\boldsymbol{x}^{\prime}(\omega)$ and $\boldsymbol{p}(\omega)$ respectively, $\boldsymbol{\xi}(t)$ becomes:

[TABLE]

The choice of the $n(\omega)$ and $\boldsymbol{c}(\omega)$ specifies the memory function and therefore (by the fluctuation-dissipation relation) the statistical properties of the noise process. We write $\boldsymbol{\kappa}(t)$ as an inverse Fourier transform of a measure:

[TABLE]

where the measure is absolutely continuous with respect to the Lebesgue measure, with the density $\boldsymbol{S}(\omega)=\pi\hat{\boldsymbol{\kappa}}_{c}(\omega)\geq 0$. The density $\boldsymbol{S}(\omega)$ is known as the spectral density of the bath.

In the following examples, we take $n(\omega)=2\omega^{2}/\pi$ (Debye-type spectrum for phonon bath).

Example 1.

If we choose $\boldsymbol{c}(\omega)\in\mathbb{R}^{d\times d}$ to be the identity matrix $\boldsymbol{I}$ multiplied by a scalar constant that is independent of $\omega$, then $\boldsymbol{\kappa}(t)$ is proportional to $\delta(t)\boldsymbol{I}$. This leads to a Langevin equation driven by white noise, in which the damping term is instantaneous. In this case, we have the SDE system for $(\boldsymbol{x}_{t},\boldsymbol{v}_{t})\in\mathbb{R}^{d\times d}$:

[TABLE]

where $\boldsymbol{\xi}_{t}$ is a white noise.

Example 2.

If we choose $\boldsymbol{c}(\omega)\in\mathbb{R}^{d\times d}$ to be the diagonal matrix with the $k$th entry

[TABLE]

where the $\alpha_{k}>0$, then we have:

[TABLE]

where $\boldsymbol{A}$ is the constant diagonal matrix with the $k$th entry equal $\alpha_{k}$. This gives SIDE (15). On the other hand, choosing $\boldsymbol{c}(\omega)$ to be the diagonal matrix with the $k$th entry

[TABLE]

allows us to obtain the covariance function of a harmonic noise process, where the $\omega_{kk}$ and $\tau_{kk}$ are the diagonal entries of the matrices $\boldsymbol{\Omega}$ and $\boldsymbol{\tau}$ respectively. In the general case where $\boldsymbol{\kappa}(t)$ is of the form (3), one may take $\boldsymbol{M}_{1}=\boldsymbol{I}$, $\boldsymbol{\Gamma}_{1}$ to be positive definite (so that the Lyapunov equation gives $\boldsymbol{\Gamma}_{1}=\boldsymbol{\Sigma}_{1}\boldsymbol{\Sigma}_{1}^{*}/2$) and choose

[TABLE]

Conclusion. The goal of this paper is to study effective dynamics of systems which can be modeled by equations of the form (11). As argued above, Hamiltonian systems describing particles interacting with heat baths can be modeled by equations of the form (11).

Appendix B Outline of the Proof of Theorem 3.1

We provide minimal outline of the proof of Theorem 3.1 in the following.

B.1 Derivation of Limiting SDE

We start by rewriting SDE (31) as

[TABLE]

The integral form of the above is given by

[TABLE]

We are interested in the limit as $m\to 0$ of the process $\boldsymbol{x}^{m}_{t}$. As $m\to 0$, we expect the sum of the second and third integral terms in the right hand side above to converge to $\int_{0}^{t}\boldsymbol{\gamma}^{-1}(\boldsymbol{X}_{s})\boldsymbol{F}(\boldsymbol{X}_{s})ds+\int_{0}^{t}\boldsymbol{\gamma}^{-1}(\boldsymbol{X}_{s})\boldsymbol{\sigma}(\boldsymbol{X}_{s})d\boldsymbol{W}_{s}$.

To examine the asymptotics of the first integral term when $m$ becomes small, we integrate by parts to write its $i$th component as:

[TABLE]

Note that the product $m(v^{m}_{s})_{j}(v^{m}_{s})_{l}$ is the $(j,l)$-entry of the matrix $m\boldsymbol{v}^{m}_{s}(\boldsymbol{v}^{m}_{s})^{*}$.

We now examine the asymptotic behavior of the above expression in the limit as $m\to 0$. Following hottovy2015smoluchowski , we express the matrix $m\boldsymbol{v}^{m}_{s}(\boldsymbol{v}^{m}_{s})^{*}$, $s\in[0,t]$, as a solution to an equation by applying Itô’s formula to the matrix $(m\boldsymbol{v}^{m}_{s})(m(\boldsymbol{v}^{m}_{s})^{*})$. This leads to:

[TABLE]

where

[TABLE]

Denoting $m\boldsymbol{v}^{m}_{s}(\boldsymbol{v}^{m}_{s})^{*}ds$ by $\boldsymbol{V}$, $-\boldsymbol{\gamma}(\boldsymbol{x}^{m}_{s})$ by $\boldsymbol{Q}$ and letting

[TABLE]

we can write the above equation as the following Lyapunov equation

[TABLE]

By Assumption 2, all real parts of the eigenvalues of $\boldsymbol{Q}$ are negative, thus the Lyapunov equation has a unique solution given by:

[TABLE]

Writing this out explicitly, we obtain:

[TABLE]

Based on the prior result in hottovy2015smoluchowski , we expect that only the second term on the right hand side has a nonzero limit as $m\to 0$. The other terms are expected to vanish as $m\to 0$. Thus, in the limit $m\to 0$, we expect that $m\boldsymbol{v}^{m}_{t}(\boldsymbol{v}^{m}_{t})^{*}$ converges to the solution, $\boldsymbol{J}$, of the Lyapunov equation:

[TABLE]

given in the statement of Theorem 3.1 (see eqn. (35)).

B.2 Moment Estimates

To justify the above convergence arguments, we provide estimate on the $p$th moment of the momentum process $\boldsymbol{p}^{m}_{t}:=m\boldsymbol{v}^{m}_{t}$, in the limit $m\to 0$.

Proposition 1.

Suppose that Assumption 1-4 hold. For all $p\geq 1$, $T>0$, there exists a positive random variable $m_{1}$ such that:

[TABLE]

as $m\to 0$.

Proof.

For $t\in[0,T]$, $m>0$, the process $\boldsymbol{p}_{t}^{m}$ satisfies the SDE:

[TABLE]

Let $\tau\in[0,t]$ and rewrite the above equation as:

[TABLE]

which admits the following solution representation:

[TABLE]

We set $\tau=t$ in the above representation to obtain:

[TABLE]

Therefore,

[TABLE]

By assumption on the boundedness and spectrum of $\boldsymbol{\gamma}(\boldsymbol{x}_{\tau})\in\mathbb{R}^{n\times n}$ (see Assumption 1-2), there exists positive constants $\kappa>0$ and $C>0$ such that

[TABLE]

for all $s,\tau\in[0,T]$, $m>0$. Indeed, applying the formula (A.2.4) in kabanov2013two , one has:

[TABLE]

where $\Lambda:=\min_{k}Re(\lambda_{k})>0$ and the $\lambda_{k}$ are the eigenvalues of $\boldsymbol{\gamma}$. Therefore, there exists a constant $C>0$ such that $\left\|e^{-\frac{\boldsymbol{\gamma}(\boldsymbol{x}^{m}_{\tau})}{m}s}\right\|\leq Ce^{-\frac{\Lambda}{2m}s}$.

Using this, we obtain the following $\mathbb{P}$-a.s. estimate:

[TABLE]

where $\tilde{C}>0$ is a constant.

Next, the key observation on (193) is that the term in parenthesis in the second line above can be made small by choosing a sufficiently small $m$. More precisely, using an adapted version of Lemma A.2.4 in kabanov2013two , there exists a (generally random) $m_{1}>0$ such that for $m\leq m_{1}$,

[TABLE]

and so:

[TABLE]

Therefore, for $p\geq 1$,

[TABLE]

where $C_{1}(p)$, $C_{2}(p)$ and $C_{3}(p)$ are some positive constants.

We estimate the last term in the above, using in particular the Burkholder-Davis-Gundy inequality karatzas2012Brownian :

[TABLE]

where $C(p,n)$, $\tilde{C}(p,n)$ are positive constants dependent on $p$ and $n$, and $\|\cdot\|_{F}$ denotes Frobenius norm. Using this estimate, (197) and Assumption 3, we see that $\mathbb{E}\left[\sup_{t\in[0,T]}|\boldsymbol{p}_{t}^{m}|^{p};m\leq m_{1}\right]\to 0$ as $m\to 0$. ∎

We also need the following estimate on a class of integrals with respect to products of the components of the momentum process $\boldsymbol{p}^{m}_{t}=m\boldsymbol{v}^{m}_{t}$. The estimate is a straightforward modification of the one given in Proposition 2.3 in birrell2017homogenization .

Proposition 2.

Suppose that Assumption 1-4 hold. Let $h:\mathbb{R}^{n}\to\mathbb{R}$ be a $C_{b}^{1}$ function (i.e. continuously differentiable and bounded function) on $[0,T]$, with bounded first derivative $\boldsymbol{\nabla}_{\boldsymbol{x}}h(\boldsymbol{x})$ for every $\boldsymbol{x}\in\mathbb{R}^{n}$. Then for any $p\geq 1$, $T>0$, $i,j=1,\dots,n$,

[TABLE]

as $m\to 0$, where the $m_{1}$ is from Proposition 1. Here $i,j$ denote the components of the momentum process $\boldsymbol{p}^{m}_{t}$ in the standard basis for $\mathbb{R}^{n}$.

Using the above moment estimates and the proof techniques (the main tools are well known ordinary and stochastic integral inequalities as well as a Gronwall type argument) in birrell2017small ; birrell2017homogenization , one obtain the convergence of $\boldsymbol{x}^{m}_{t}$ to $\boldsymbol{X}_{t}$ in the limit as $m\to 0$ in the following sense: for all finite $T>0$, $p\geq 1$,

[TABLE]

as $m\to 0$, where the $m_{1}$ is from Proposition 1. This implies that for all finite $T>0$, $\sup_{t\in[0,T]}|\boldsymbol{x}_{t}^{m}-\boldsymbol{X}_{t}|\to 0$ in probability, in the limit as $m\to 0$ (see Lemma 1 below).

Lemma 1.

Let $Y_{m}\geq 0$ ($m\geq 0$) be a family of random variables. Suppose that there exists a strictly positive random variable $m^{*}(\omega)$ such that

[TABLE]

for some $p\geq 1$. Then $Y_{m}\to 0$ in probability as $m\to 0$.

Proof.

Let $g$ be a bounded continuous function on $\mathbb{R}_{+}$. We have

[TABLE]

The first term is bounded from above by $(\sup|g|)\mathbb{P}\left[m^{*}\leq m\right]$ which goes to zero as $m\to 0$, since $\cap_{m>0}\{m^{*}\leq m\}=\emptyset$. To estimate the second term, let $\epsilon>0$. Choose $\delta>0$ such that $|g(y)|<\epsilon$ whenever $y<\delta$. We have

[TABLE]

The first term is bounded from above by $\epsilon$. The function ${g(y)\over y}$ on $\{y:y\geq\delta\}$ is bounded by some constant $M$, so the second term on the right-hand side of the above equation is bounded from above by $M\mathbb{E}\left[Y_{m};m^{*}\geq m\right]$ which goes to zero as $m\to 0$ by assumption. The lemma is proven. ∎

We end this appendix with a remark on the mode of convergence stated in Theorem 1.

Remark 14.

Provided that for $p\geq 1$, there exists a constant $C>0$ such that $\mathbb{E}\sup_{t\in[0,T]}|\boldsymbol{p}_{t}^{m}|^{p}<C$ for all $m>0$ (so that the family of random variables $(\sup_{t\in[0,T]}|\boldsymbol{p}_{t}^{m}|^{p})_{m>0}$ is uniformly integrable), one could, using Proposition 1, Lemma 1 and Theorem 13.7 in williams1991probability , obtain $L^{p}$-convergence of $\sup_{t\in[0,T]}|\boldsymbol{p}_{t}^{m}|$ to zero and hence strengthen the convergence result stated in Theorem 1 to $L^{p}$-convergence. The uniform integrability condition is satisfied if $\mathbb{E}\sup_{t\in[0,T]}\|\boldsymbol{\Phi}^{m}(t)\|^{p}$ is bounded uniformly in $m$, where $\boldsymbol{\Phi}^{m}(t)$ is the fundamental matrix that solves the random initial value problem:

[TABLE]

However, it is not obvious how one could verify the latter condition from our assumptions on $\boldsymbol{\gamma}$. Roughly speaking, one does not have a good control of $\sup_{t\in[0,T]}|\boldsymbol{p}_{t}^{m}|^{p}$ outside of the set $\{m\leq m_{1}\}$. If $\boldsymbol{\gamma}$ was, in addition, symmetric (and so all the (real) eigenvalues of $\boldsymbol{\gamma}$ are bounded from below by a positive constant – c.f. birrell2017homogenization ), the condition can be easily verified.